Skip to content
Snippets Groups Projects
Select Git revision
  • 94d80bb82f523abcbf26b89346cef4d3d9d9a997
  • master default protected
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

BergotBasis.cpp

Blame
  • BergotBasis.cpp 5.74 KiB
    // Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to the public mailing list <gmsh@onelab.info>.
    
    #include <cmath>
    #include "BergotBasis.h"
    #include "MElement.h"
    #include "orthogonalBasis.h"
    
    BergotBasis::BergotBasis(int p, bool incpl) : order(p), incomplete(incpl)
    {
      if(incomplete && order > 2) {
        Msg::Error("Incomplete pyramids of order %i not yet implemented", order);
      }
    }
    
    BergotBasis::~BergotBasis() {}
    
    bool BergotBasis::validIJ(int i, int j) const
    {
      if(!incomplete) return (i <= order) && (j <= order);
      if(i + j <= order) return true;
      if(i + j == order + 1) return i == 1 || j == 1;
      return false;
    }
    
    // Values of Bergot basis functions for coordinates (u,v,w) in the unit pyramid:
    // f = L_i(uhat)*L_j(vhat)*(1-w)^max(i,j)*P_k^{2*max(i,j)+2,0}(what)
    // with i,j = 0...order and k = 0...order-max(i,j)
    // and (uhat,vhat,what) = (u/(1-w),v/(1-w),2*w-1) reduced coordinates in the
    // unit cube [-1,1]^3
    void BergotBasis::f(double u, double v, double w, double *val) const
    {
      // Compute Legendre polynomials at uhat
      double uhat = (w == 1.) ? 0. : u / (1. - w);
      std::vector<double> uFcts(order + 1);
      LegendrePolynomials::f(order, uhat, &(uFcts[0]));
    
      // Compute Legendre polynomials at vhat
      double vhat = (w == 1.) ? 0. : v / (1. - w);
      std::vector<double> vFcts(order + 1);
      LegendrePolynomials::f(order, vhat, &(vFcts[0]));
    
      // Compute Jacobi polynomials at what
      double what = 2. * w - 1.;
      std::vector<std::vector<double> > wFcts(order + 1), wGrads(order + 1);
      for(int mIJ = 0; mIJ <= order; mIJ++) {
        int kMax = order - mIJ;
        std::vector<double> &wf = wFcts[mIJ];
        wf.resize(kMax + 1);
        JacobiPolynomials::f(kMax, 2. * mIJ + 2., 0., what, &(wf[0]));
      }
    
      // Recombine to find shape function values
      int index = 0;
      for(int i = 0; i <= order; i++) {
        for(int j = 0; j <= order; j++) {
          if(validIJ(i, j)) {
            int mIJ = std::max(i, j);
            double fact = pow(1. - w, mIJ);
            std::vector<double> &wf = wFcts[mIJ];
            for(int k = 0; k <= order - mIJ; k++, index++)
              val[index] = uFcts[i] * vFcts[j] * wf[k] * fact;
          }
        }
      }
    }
    
    // Derivatives of Bergot basis functions for coordinates (u,v,w) in the unit