Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

BasisFunction.cc 4.26 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
18

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(const ::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

25
    nDOF = NEW DimVec<int>(dim, DEFAULT_VALUE, -1);
Thomas Witkowski's avatar
Thomas Witkowski committed
26
27
28
29
30
31
32
33
34
    dow = Global::getGeo(WORLD);

    grdTmpVec1.resize(omp_get_max_threads());
    grdTmpVec2.resize(omp_get_max_threads());

    for (int i = 0; i < omp_get_max_threads(); i++) {
      grdTmpVec1[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      grdTmpVec2[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
    }
35
36
  };

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

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
44
45
      DELETE grdTmpVec1[i];
      DELETE grdTmpVec2[i];
    }
  }
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

  /****************************************************************************/
  /*  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);
  }


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

  const WorldMatrix<double>& BasisFunction::evalD2Uh(const DimVec<double>& lambda,
						     const DimVec<WorldVector<double> >& grd_lambda,
						     const double *uh_loc, WorldMatrix<double>* D2_uh) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
    static WorldMatrix<double> D2(DEFAULT_VALUE, 0.0);
121
    DimMat<double> D2_b(dim, DEFAULT_VALUE, 0.0);
122
    DimMat<double> D2_tmp(dim, DEFAULT_VALUE, 0.0);
123
    WorldMatrix<double> *val = D2_uh ? D2_uh : &D2;
124
125
  
    for (int i = 0; i < nBasFcts; i++) {
126
      (*(*d2Phi)[i])(lambda, D2_b);
127
128
129
130
131
132
133
134
      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;
135
136
	for (int k = 0; k < dim + 1; k++)
	  for (int l = 0; l < dim + 1; l++)
137
138
139
	    (*val)[i][j] += grd_lambda[k][i] * grd_lambda[l][j] * D2_tmp[k][l];
      }
    
140
    return ((*val));
141
142
143
144
145
146
147
148
149
  }


  int BasisFunction::getNumberOfDOFs(int i) const
  { 
    return (*nDOF)[i];
  }

}