BasisFunction.cc 4.33 KB
Newer Older
1 2 3 4 5 6
#include <algorithm>

#include "FixVec.h"
#include "DOFVector.h"
#include "BasisFunction.h"
#include "Lagrange.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
7
#include "OpenMP.h"
8 9 10 11 12 13 14 15 16 17

namespace AMDiS {


  /****************************************************************************/
  /*  Lagrangian basisfunctions of order 0-4; these                           */
  /*  functions are evaluated in barycentric coordinates; the derivatives     */
  /*  are those corresponding to these barycentric coordinates.               */
  /****************************************************************************/

Thomas Witkowski's avatar
Thomas Witkowski committed
18
  BasisFunction::BasisFunction(std::string name_, int dim_, int degree_)
Thomas Witkowski's avatar
Thomas Witkowski committed
19 20 21
    : name(name_), 
      degree(degree_), 
      dim(dim_)
22 23
  {
    FUNCNAME("BasisFunction::BasisFunction()");
Thomas Witkowski's avatar
Thomas Witkowski committed
24

Thomas Witkowski's avatar
Thomas Witkowski committed
25
    nDOF = new DimVec<int>(dim, DEFAULT_VALUE, -1);
Thomas Witkowski's avatar
Thomas Witkowski committed
26 27
    dow = Global::getGeo(WORLD);

28 29
    grdTmpVec1.resize(omp_get_overall_max_threads());
    grdTmpVec2.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
30

31
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
32 33
      grdTmpVec1[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      grdTmpVec2[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
34
    }
35
  }
36

Thomas Witkowski's avatar
Thomas Witkowski committed
37 38
  BasisFunction::~BasisFunction()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
39
    delete nDOF;
Thomas Witkowski's avatar
Thomas Witkowski committed
40

Thomas Witkowski's avatar
Thomas Witkowski committed
41
    for (int i = 0; i < static_cast<int>(grdTmpVec1.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
42 43
      delete grdTmpVec1[i];
      delete grdTmpVec2[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
44 45
    }
  }
46 47 48 49 50 51 52 53

  /****************************************************************************/
  /*  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<double>& lambda,
54
			       const ElementVector& uh_loc) const
55 56 57 58 59 60 61 62 63 64
  {
    double val = 0.0;

    for (int i = 0; i < nBasFcts; i++)
      val += uh_loc[i] * (*(*phi)[i])(lambda);

    return(val);
  }


65
  const WorldVector<double>& BasisFunction::evalUh(const DimVec<double>& lambda, 
66
						   const mtl::dense_vector<WorldVector<double> >& uh_loc, 
67 68
						   WorldVector<double>* values) const 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
70 71 72
    WorldVector<double> *val = (NULL != values) ? values : &Values;   
    
    for (int n = 0; n < dow; n++)
Thomas Witkowski's avatar
Thomas Witkowski committed
73
      (*val)[n] = 0.0;
74 75 76 77 78 79 80 81 82 83 84
    
    for (int i = 0; i < nBasFcts; i++) {
      double phil = (*(*phi)[i])(lambda);
      for (int n = 0; n < dow; n++)
	(*val)[n] += uh_loc[i][n] * phil;
    }
    
    return((*val));
  }


85 86
  const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda, 
						      const DimVec<WorldVector<double> >& grd_lambda,
87
						      const ElementVector& uh_loc,  
Thomas Witkowski's avatar
Thomas Witkowski committed
88
						      WorldVector<double>* val) const
89
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
90 91
    TEST_EXIT_DBG(val)("return value is NULL\n");

92 93 94 95
    int myRank = omp_get_thread_num();

    DimVec<double> *grdTmp1 = grdTmpVec1[myRank];
    DimVec<double> *grdTmp2 = grdTmpVec2[myRank];
96 97

    for (int j = 0; j < dim + 1; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
98
      (*grdTmp2)[j] = 0.0;
99 100

    for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
101
      (*(*grdPhi)[i])(lambda, *grdTmp1);
102

103
      for (int j = 0; j < dim + 1; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
104
	(*grdTmp2)[j] += uh_loc[i] * (*grdTmp1)[j];
105 106 107
    }

    for (int i = 0; i < dow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
108
      (*val)[i] = 0.0;
109
      for (int j = 0; j < dim + 1; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
110
	(*val)[i] += grd_lambda[j][i] * (*grdTmp2)[j];
111 112
    }
  
113
    return ((*val));
114 115
  }

116

117 118
  const WorldMatrix<double>& BasisFunction::evalD2Uh(const DimVec<double>& lambda,
						     const DimVec<WorldVector<double> >& grd_lambda,
Thomas Witkowski's avatar
Thomas Witkowski committed
119 120
						     const ElementVector& uh_loc, 
						     WorldMatrix<double>* D2_uh) const
121
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
122
    static WorldMatrix<double> D2(DEFAULT_VALUE, 0.0);
123
    DimMat<double> D2_b(dim, DEFAULT_VALUE, 0.0);
124
    DimMat<double> D2_tmp(dim, DEFAULT_VALUE, 0.0);
125
    WorldMatrix<double> *val = D2_uh ? D2_uh : &D2;
126 127
  
    for (int i = 0; i < nBasFcts; i++) {
128
      (*(*d2Phi)[i])(lambda, D2_b);
129 130 131 132 133 134 135 136
      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];
    }

    for (int i = 0; i < dow; i++)
      for (int j = 0; j < dow; j++) {
	(*val)[i][j] = 0.0;
137 138
	for (int k = 0; k < dim + 1; k++)
	  for (int l = 0; l < dim + 1; l++)
139 140 141
	    (*val)[i][j] += grd_lambda[k][i] * grd_lambda[l][j] * D2_tmp[k][l];
      }
    
142
    return ((*val));
143 144
  }

145

146
  int BasisFunction::getNumberOfDofs(int i) const
147 148 149 150 151
  { 
    return (*nDOF)[i];
  }

}