BasisFunction.cc 4.02 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#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.               */
  /****************************************************************************/

Thomas Witkowski's avatar
Thomas Witkowski committed
29
  BasisFunction::BasisFunction(std::string name_, int dim_, int degree_)
Thomas Witkowski's avatar
Thomas Witkowski committed
30
31
32
    : name(name_), 
      degree(degree_), 
      dim(dim_)
33
34
  {
    FUNCNAME("BasisFunction::BasisFunction()");
Thomas Witkowski's avatar
Thomas Witkowski committed
35

Thomas Witkowski's avatar
Thomas Witkowski committed
36
    nDOF = new DimVec<int>(dim, DEFAULT_VALUE, -1);
Thomas Witkowski's avatar
Thomas Witkowski committed
37
    dow = Global::getGeo(WORLD);
38
  }
39

Thomas Witkowski's avatar
Thomas Witkowski committed
40
41
  BasisFunction::~BasisFunction()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
42
    delete nDOF;
Thomas Witkowski's avatar
Thomas Witkowski committed
43
  }
44
45
46
47
48
49
50
51

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

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

    return(val);
  }


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


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

90
91
    mtl::dense_vector<double> grdTmp1(dim + 1);
    mtl::dense_vector<double> grdTmp2(dim + 1, 0.0);
92
93

    for (int i = 0; i < nBasFcts; i++) {
94
      (*(*grdPhi)[i])(lambda, grdTmp1);
95

96
      for (int j = 0; j < dim + 1; j++)
97
	grdTmp2[j] += uh_loc[i] * grdTmp1[j];
98
99
100
    }

    for (int i = 0; i < dow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
101
      (*val)[i] = 0.0;
102
      for (int j = 0; j < dim + 1; j++)
103
	(*val)[i] += grd_lambda[j][i] * grdTmp2[j];
104
105
    }
  
106
    return ((*val));
107
108
  }

109

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

138

139
  int BasisFunction::getNumberOfDofs(int i) const
140
141
142
143
144
  { 
    return (*nDOF)[i];
  }

}