Skip to content
Snippets Groups Projects
Select Git revision
  • 07070da8458b307f08234cdc4843ef7a0c4aaa74
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

TetNodeBasis.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    TetNodeBasis.cpp 4.86 KiB
    #include "TetNodeBasis.h"
    #include "TetReferenceSpace.h"
    #include "Legendre.h"
    
    using namespace std;
    
    TetNodeBasis::TetNodeBasis(size_t order){
      // Reference Space //
      refSpace  = new TetReferenceSpace;
      nRefSpace = getReferenceSpace().getNReferenceSpace();
    
      const vector<vector<vector<size_t> > >&
        edgeIdx = refSpace->getEdgeNodeIndex();
    
      const vector<vector<vector<size_t> > >&
        faceIdx = refSpace->getFaceNodeIndex();
    
      // Set Basis Type //
      this->order = order;
    
      type = 0;
      dim  = 3;
    
      nVertex   = 4;
      nEdge     = 6 * (order - 1);
      nFace     = 2 * (order - 1) * (order - 2);
      nCell     =     (order - 1) * (order - 2) * (order - 3) / 6;
      nFunction = nVertex + nEdge + nFace + nCell;
    
      // Alloc Temporary Space //
      const int orderMinus      = order - 1;
      const int orderMinusTwo   = order - 2;
      const int orderMinusThree = order - 3;
    
      Polynomial* legendre     = new Polynomial[order];
      Polynomial* sclLegendre  = new Polynomial[order];
      Polynomial* intLegendre  = new Polynomial[order];
    
      // Legendre Polynomial //
      Legendre::legendre(legendre, orderMinus);
      Legendre::scaled(sclLegendre, orderMinus);
      Legendre::intScaled(intLegendre, order);
    
      // Lagrange Polynomial //
      const Polynomial lagrange[4] =
        {
          Polynomial(Polynomial(1, 0, 0, 0) -
                     Polynomial(1, 1, 0, 0) -
                     Polynomial(1, 0, 1, 0) -
                     Polynomial(1, 0, 0, 1)),
    
          Polynomial(Polynomial(1, 1, 0, 0)),
    
          Polynomial(Polynomial(1, 0, 1, 0)),
    
          Polynomial(Polynomial(1, 0, 0, 1))
        };
    
    
      // Basis //
      basis = new Polynomial**[nRefSpace];
    
      for(size_t s = 0; s < nRefSpace; s++)
        basis[s] = new Polynomial*[nFunction];
    
      // Vertex Based //
      for(size_t s = 0; s < nRefSpace; s++){
        basis[s][0] = new Polynomial(lagrange[0]);
        basis[s][1] = new Polynomial(lagrange[1]);
        basis[s][2] = new Polynomial(lagrange[2]);
        basis[s][3] = new Polynomial(lagrange[3]);
      }
    
      // Edge Based //
      for(size_t s = 0; s < nRefSpace; s++){
        size_t i = nVertex;
    
        for(int e = 0; e < 6; e++){
          for(size_t l = 1; l < order; l++){
            basis[s][i] =
              new Polynomial(intLegendre[l].compose
                             (lagrange[edgeIdx[s][e][0]] -
                              lagrange[edgeIdx[s][e][1]]
                              ,
                              lagrange[edgeIdx[s][e][0]] +
                              lagrange[edgeIdx[s][e][1]]));
            i++;
          }
        }
      }
    
      // Face Based //
      // TO CHECK: Are Triangles face matching tets ?
    
      for(size_t s = 0; s < nRefSpace; s++){
        size_t i = nVertex + nEdge;
    
        for(int f = 0; f < 4; f++){
          for(int l1 = 1; l1 < orderMinus; l1++){
            for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
              Polynomial sum =
                lagrange[faceIdx[s][f][0]] +
                lagrange[faceIdx[s][f][1]] +
                lagrange[faceIdx[s][f][2]];
    
              basis[s][i] =
                new Polynomial(intLegendre[l1].compose
                               (lagrange[faceIdx[s][f][0]] -
                                lagrange[faceIdx[s][f][1]]
                                ,
                                lagrange[faceIdx[s][f][0]] +
                                lagrange[faceIdx[s][f][1]])
    
                               *
    
                               lagrange[faceIdx[s][f][2]]
                               *
                               sclLegendre[l2].compose
                               (lagrange[faceIdx[s][f][2]] * 2 - sum, sum));
              i++;
            }
          }
        }
      }
    
      // Cell Based //
      const Polynomial oneMinusFour         = Polynomial(1, 0, 0, 0) - lagrange[3];
      const Polynomial twoThreeOneMinusFour = lagrange[2] * 2 - oneMinusFour;
      const Polynomial twoFourMinusOne      = lagrange[3] * 2 - Polynomial(1, 0, 0, 0);
    
      const Polynomial sub = lagrange[0] - lagrange[1];
      const Polynomial add = lagrange[0] + lagrange[1];
    
      for(size_t s = 0; s < nRefSpace; s++){
        size_t i = nVertex + nEdge + nFace;
    
        for(int l1 = 1; l1 < orderMinusTwo; l1++){
          for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
            for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
              basis[s][i] =
                new Polynomial(intLegendre[l1].compose(sub, add)             *
                               lagrange[2]                                   *
                               sclLegendre[l2].compose(twoThreeOneMinusFour,
                                                       oneMinusFour)         *
                               lagrange[3]                                   *
                               legendre[l3].compose(twoFourMinusOne));
              i++;
            }
          }
        }
      }
    
      // Free Temporary Sapce //
      delete[] legendre;
      delete[] sclLegendre;
      delete[] intLegendre;
    }
    
    TetNodeBasis::~TetNodeBasis(void){
      // ReferenceSpace //
      delete refSpace;
    
      // Basis //
      for(size_t i = 0; i < nRefSpace; i++){
        for(size_t j = 0; j < nFunction; j++)
          delete basis[i][j];
    
        delete[] basis[i];
      }
    
      delete[] basis;
    }