BasisFunction.cc 3.29 KB
Newer Older
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 <algorithm>

#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<int>(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<double>& 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<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda, 
						      const DimVec<WorldVector<double> >& grd_lambda,
						      const double *uh_loc,  WorldVector<double>* grd_uh) const
  {
    static WorldVector<double>  grd;

    DimVec<double> grd_b(dim, DEFAULT_VALUE, 0.);
    WorldVector<double>     *val;
    DimVec<double> 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<double>& BasisFunction::evalD2Uh(const DimVec<double>& lambda,
						     const DimVec<WorldVector<double> >& grd_lambda,
						     const double *uh_loc, WorldMatrix<double>* D2_uh) const
  {
    static WorldMatrix<double> D2(DEFAULT_VALUE, 0.);
    DimMat<double> D2_b(dim, NO_INIT);
    WorldMatrix<double>    *val;
    DimMat<double> 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];
  }

}