Intrepid
Intrepid_HGRAD_POLY_C1_FEMDef.hpp
2 #include "Sacado_Fad_DFad.hpp"
3 
4 
5 namespace 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 
214 
215 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
216 #ifdef __GNUC__
217 #warning "The Intrepid package is deprecated"
218 #endif
219 #endif
220 
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void evaluateWeightFunctions(ArrayScalar1 &outputValues, const ArrayScalar1 &inputValues, const ArrayScalar1 &cellVertices) const
Evaluation of the Wachspress weight functions.
Header file for utility class to provide multidimensional containers.
void shapeFunctions(ArrayScalar1 &outputValues, const ArrayScalar1 &inputValues, const ArrayScalar1 &cellVertices) const
Evaluation of Wachspress shape functions.
Scalar1 computeArea(const ArrayScalar1 &p1, const ArrayScalar1 &p2, const ArrayScalar1 &p3) const
Helper function to compute area of triangle formed by 3 points.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM reference basis evaluation: invocation of this method throws an exception.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
Basis_HGRAD_POLY_C1_FEM(const shards::CellTopology &cellTopology)
Constructor.