BasisFunction.cc 3.29 KB
 Peter Gottschling committed Feb 15, 2008 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 ``````#include #include "FixVec.h" #include "DOFVector.h" #include "BasisFunction.h" #include "Lagrange.h" namespace AMDiS { /****************************************************************************/ /* Lagrangian basisfunctions of order 0-4; these */ /* functions are evaluated in barycentric coordinates; the derivatives */ /* are those corresponding to these barycentric coordinates. */ /****************************************************************************/ BasisFunction::~BasisFunction() { DELETE nDOF; } BasisFunction::BasisFunction(const ::std::string& name_, int dim_, int degree_) : name(name_), degree(degree_), dim(dim_) { FUNCNAME("BasisFunction::BasisFunction()"); nDOF = NEW DimVec(dim, DEFAULT_VALUE, -1); }; /****************************************************************************/ /* some routines for evaluation of a finite element function, its gradient */ /* and second derivatives; all those with _fast use the preevaluated */ /* basis functions at that point. */ /****************************************************************************/ double BasisFunction::evalUh(const DimVec& lambda, const double *uh_loc) const { double val = 0.0; for (int i = 0; i < nBasFcts; i++) val += uh_loc[i] * (*(*phi)[i])(lambda); return(val); } const WorldVector& BasisFunction::evalGrdUh(const DimVec& lambda, const DimVec >& grd_lambda, const double *uh_loc, WorldVector* grd_uh) const { static WorldVector grd; DimVec grd_b(dim, DEFAULT_VALUE, 0.); WorldVector *val; DimVec grd1(dim, DEFAULT_VALUE, 0.); val = grd_uh ? grd_uh : &grd; for (int j = 0; j < dim + 1; j++) grd1[j] = 0.0; for (int i = 0; i < nBasFcts; i++) { grd_b = (*(*grdPhi)[i])(lambda); for (int j = 0; j < dim + 1; j++) grd1[j] += uh_loc[i] * grd_b[j]; } int dow = Global::getGeo(WORLD); for (int i = 0; i < dow; i++) { (*val)[i] = 0; for (int j = 0; j < dim + 1; j++) (*val)[i] += grd_lambda[j][i] * grd1[j]; } return((*val)); } const WorldMatrix& BasisFunction::evalD2Uh(const DimVec& lambda, const DimVec >& grd_lambda, const double *uh_loc, WorldMatrix* D2_uh) const { static WorldMatrix D2(DEFAULT_VALUE, 0.); DimMat D2_b(dim, NO_INIT); WorldMatrix *val; DimMat D2_tmp(dim, DEFAULT_VALUE, 0.0); val = D2_uh ? D2_uh : &D2; for (int i = 0; i < nBasFcts; i++) { D2_b = (*(*d2Phi)[i])(lambda); for (int k = 0; k < dim + 1; k++) for (int l = 0; l < dim + 1; l++) D2_tmp[k][l] += uh_loc[i] * D2_b[k][l]; } int dow = Global::getGeo(WORLD); for (int i = 0; i < dow; i++) for (int j = 0; j < dow; j++) { (*val)[i][j] = 0.0; for (int k = 0; k < dim+1; k++) for (int l = 0; l < dim+1; l++) (*val)[i][j] += grd_lambda[k][i] * grd_lambda[l][j] * D2_tmp[k][l]; } return((*val)); } int BasisFunction::getNumberOfDOFs(int i) const { return (*nDOF)[i]; } }``````