Intrepid
Intrepid_HGRAD_POLY_C1_FEMDef.hpp
2#include "Sacado_Fad_DFad.hpp"
3
4
5namespace Intrepid{
6
7 template<class Scalar, class ArrayScalar>
9 this -> basisCardinality_ = cellTopology.getNodeCount();
10 this -> basisDegree_ = 1;
11 this -> basisCellTopology_ = cellTopology;
12 this -> basisType_ = BASIS_FEM_DEFAULT;
13 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
14 this -> basisTagsAreSet_ = false;
15 }
16
17
18 template<class Scalar, class ArrayScalar>
20 // Basis-dependent initializations
21 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
22 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
23 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
24 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
25
26 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
27
28 int *tags = new int[tagSize * this->getCardinality()];
29 for (int i=0;i<this->getCardinality();i++){
30 tags[4*i] = 0;
31 tags[4*i+1] = i;
32 tags[4*i+2] = 0;
33 tags[4*i+3] = 1;
34 }
35
36
37 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
38 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
39 this -> ordinalToTag_,
40 tags,
41 this -> basisCardinality_,
42 tagSize,
43 posScDim,
44 posScOrd,
45 posDfOrd);
46
47 delete [] tags;
48 }
49
50 // this is the FEM reference element version, this should not be called for polygonal basis
51 // since polygonal basis is defined on physical element.
52 template<class Scalar, class ArrayScalar>
54 const ArrayScalar& inputPoints,
55 const EOperator operatorType) const{
56 TEUCHOS_TEST_FOR_EXCEPTION ( true, std::logic_error,
57 ">>>ERROR (Basis_HGRAD_POLY_C1_FEM): Polygonal basis calling FEM member function");
58 }
59
60
61 template<class Scalar, class ArrayScalar>
63 const ArrayScalar& inputPoints,
64 const ArrayScalar& cellVertices,
65 const EOperator operatorType) const{
66
67
68 // implement wachspress basis functions
69 switch (operatorType) {
70 case OPERATOR_VALUE:
71 {
72 shapeFunctions<Scalar,ArrayScalar>(outputValues,inputPoints,cellVertices);
73 }
74 break;
75 case OPERATOR_GRAD:
76 {
77 FieldContainer<Sacado::Fad::DFad<Scalar> > dInput(inputPoints.dimension(0),inputPoints.dimension(1));
78 for (int i=0;i<inputPoints.dimension(0);i++){
79 for (int j=0;j<2;j++){
80 dInput(i,j) = Sacado::Fad::DFad<Scalar>( inputPoints(i,j));
81 dInput(i,j).diff(j,2);
82 }
83 }
84 FieldContainer<Sacado::Fad::DFad<Scalar> > cellVerticesFad(cellVertices.dimension(0),cellVertices.dimension(1));
85 for (int i=0;i<cellVertices.dimension(0);i++){
86 for (int j=0;j<cellVertices.dimension(1);j++){
87 cellVerticesFad(i,j) = Sacado::Fad::DFad<Scalar>( cellVertices(i,j) );
88 }
89 }
90
91 FieldContainer<Sacado::Fad::DFad<Scalar> > dOutput(outputValues.dimension(0),outputValues.dimension(1));
92
93 shapeFunctions<Sacado::Fad::DFad<Scalar>, FieldContainer<Sacado::Fad::DFad<Scalar> > >(dOutput,dInput,cellVerticesFad);
94
95 for (int i=0;i<outputValues.dimension(0);i++){
96 for (int j=0;j<outputValues.dimension(1);j++){
97 for (int k=0;k<outputValues.dimension(2);k++){
98 outputValues(i,j,k) = dOutput(i,j).dx(k);
99 }
100 }
101 }
102 }
103 break;
104 case OPERATOR_D1:
105 case OPERATOR_D2:
106 case OPERATOR_D3:
107 case OPERATOR_D4:
108 case OPERATOR_D5:
109 case OPERATOR_D6:
110 case OPERATOR_D7:
111 case OPERATOR_D8:
112 case OPERATOR_D9:
113 case OPERATOR_D10:
114 {
115 TEUCHOS_TEST_FOR_EXCEPTION ( true, std::invalid_argument,
116 ">>> ERROR (Basis_HGRAD_POLY_C1_FEM): operator not implemented yet");
117 }
118 break;
119 default:
120 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType)), std::invalid_argument,
121 ">>> ERROR (Basis_HGRAD_POLY_C1_FEM): Invalid operator type");
122 break;
123 }
124
125 }
126
127
128 template<class Scalar, class ArrayScalar>
129 template<class Scalar1, class ArrayScalar1>
131 const ArrayScalar1& inputPoints,
132 const ArrayScalar1& cellVertices)const{
133 int numPoints = inputPoints.dimension(0);
134 FieldContainer<Scalar1> weightFunctions( this->basisCardinality_,numPoints );
135 evaluateWeightFunctions<Scalar1, ArrayScalar1>(weightFunctions,inputPoints,cellVertices);
136 for (int pointIndex = 0;pointIndex<outputValues.dimension(1);pointIndex++){
137 Scalar1 denominator=0;
138
139 for (int k=0;k<weightFunctions.dimension(0);k++){
140 denominator += weightFunctions(k,pointIndex);
141 }
142 for (int k=0;k<outputValues.dimension(0);k++){
143 outputValues(k,pointIndex) = weightFunctions(k,pointIndex)/denominator;
144 }
145 }
146 }
147
148
149
150 template<class Scalar, class ArrayScalar>
151 template<class Scalar1, class ArrayScalar1>
153 const ArrayScalar1& p2,
154 const ArrayScalar1& p3) const{
155 Scalar1 area = 0.5*(p1(0)*(p2(1)-p3(1))
156 -p2(0)*(p1(1)-p3(1))
157 +p3(0)*(p1(1)-p2(1)));
158 return area;
159 }
160
161
162 template<class Scalar, class ArrayScalar>
163 template<class Scalar1, class ArrayScalar1>
165 const ArrayScalar1& inputValues,
166 const ArrayScalar1& cellVertices)const{
167
168
169 int spaceDim = this->basisCellTopology_.getDimension();
170 for (int k=0;k<outputValues.dimension(0);k++){
171
172 // compute a_k for each weight function
173 // need a_k = A(p_i,p_j,p_k) where nodes i and j are adjacent to node k
174 int adjIndex1 = -1, adjIndex2 = -1;
175 for (int i = 0;i < this->basisCellTopology_.getEdgeCount();i++){
176 if ( this->basisCellTopology_.getNodeMap(1,i,0) == k )
177 adjIndex1 = this->basisCellTopology_.getNodeMap(1,i,1);
178 else if ( this->basisCellTopology_.getNodeMap(1,i,1) == k )
179 adjIndex2 = this->basisCellTopology_.getNodeMap(1,i,0);
180 }
181 TEUCHOS_TEST_FOR_EXCEPTION( (adjIndex1 == -1 || adjIndex2 == -1), std::invalid_argument,
182 ">>> ERROR (Intrepid_HGRAD_POLY_C1_FEM): cannot find adjacent nodes when evaluating Wachspress weight function.");
183 FieldContainer<Scalar1> p1(spaceDim);
184 FieldContainer<Scalar1> p2(spaceDim);
185 FieldContainer<Scalar1> p3(spaceDim);
186 for (int i=0;i<spaceDim;i++){
187 p1(i) = cellVertices(adjIndex1,i);
188 p2(i) = cellVertices(k,i);
189 p3(i) = cellVertices(adjIndex2,i);
190 }
191 Scalar1 a_k = computeArea<Scalar1, ArrayScalar1>(p1,p2,p3);
192 // now compute prod_{ij!=k} r_ij
193 for (int pointIndex = 0;pointIndex < inputValues.dimension(0);pointIndex++){
194 Scalar1 product = a_k;
195 for (int edgeIndex = 0;edgeIndex < this->basisCellTopology_.getEdgeCount();edgeIndex++){
196 int edgeNode1 = this->basisCellTopology_.getNodeMap(1,edgeIndex,0);
197 int edgeNode2 = this->basisCellTopology_.getNodeMap(1,edgeIndex,1);
198 if ( edgeNode1 != k && edgeNode2 != k ){
199 for (int i=0;i<spaceDim;i++){
200 p1(i) = inputValues(pointIndex,i);
201 p2(i) = cellVertices(edgeNode1,i);
202 p3(i) = cellVertices(edgeNode2,i);
203 }
204 product *= computeArea<Scalar1, ArrayScalar1>(p1,p2,p3);
205 }
206 }
207 outputValues(k,pointIndex) = product;
208 }
209 }
210 }
211} // namespace Intrepid
212
213
Header file for utility class to provide multidimensional containers.
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
Scalar1 computeArea(const ArrayScalar1 &p1, const ArrayScalar1 &p2, const ArrayScalar1 &p3) const
Helper function to compute area of triangle formed by 3 points.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM reference basis evaluation: invocation of this method throws an exception.
void evaluateWeightFunctions(ArrayScalar1 &outputValues, const ArrayScalar1 &inputValues, const ArrayScalar1 &cellVertices) const
Evaluation of the Wachspress weight functions.
void shapeFunctions(ArrayScalar1 &outputValues, const ArrayScalar1 &inputValues, const ArrayScalar1 &cellVertices) const
Evaluation of Wachspress shape functions.
Basis_HGRAD_POLY_C1_FEM(const shards::CellTopology &cellTopology)
Constructor.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int dimension(const int whichDim) const
Returns the specified dimension.