Intrepid
Intrepid_HGRAD_TRI_C2_FEMDef.hpp
1 #ifndef INTREPID_HGRAD_TRI_C2_FEMDEF_HPP
2 #define INTREPID_HGRAD_TRI_C2_FEMDEF_HPP
3 // @HEADER
4 // ************************************************************************
5 //
6 // Intrepid Package
7 // Copyright (2007) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40 // Denis Ridzal (dridzal@sandia.gov), or
41 // Kara Peterson (kjpeter@sandia.gov)
42 //
43 // ************************************************************************
44 // @HEADER
45 
51 namespace Intrepid {
52 
53 
54 template<class Scalar, class ArrayScalar>
56  {
57  this -> basisCardinality_ = 6;
58  this -> basisDegree_ = 2;
59  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
60  this -> basisType_ = BASIS_FEM_DEFAULT;
61  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
62  this -> basisTagsAreSet_ = false;
63  }
64 
65 
66 
67 template<class Scalar, class ArrayScalar>
69 
70  // Basis-dependent initializations
71  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
72  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
73  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
74  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
75 
76  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
77  int tags[] = { 0, 0, 0, 1,
78  0, 1, 0, 1,
79  0, 2, 0, 1,
80  1, 0, 0, 1,
81  1, 1, 0, 1,
82  1, 2, 0, 1};
83 
84  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
85  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
86  this -> ordinalToTag_,
87  tags,
88  this -> basisCardinality_,
89  tagSize,
90  posScDim,
91  posScOrd,
92  posDfOrd);
93 }
94 
95 
96 
97 template<class Scalar, class ArrayScalar>
99  const ArrayScalar & inputPoints,
100  const EOperator operatorType) const {
101 
102  // Verify arguments
103 #ifdef HAVE_INTREPID_DEBUG
104  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
105  inputPoints,
106  operatorType,
107  this -> getBaseCellTopology(),
108  this -> getCardinality() );
109 #endif
110 
111  // Number of evaluation points = dim 0 of inputPoints
112  int dim0 = inputPoints.dimension(0);
113 
114  // Temporaries: (x,y) coordinates of the evaluation point
115  Scalar x = 0.0;
116  Scalar y = 0.0;
117 
118  switch (operatorType) {
119 
120  case OPERATOR_VALUE:
121  for (int i0 = 0; i0 < dim0; i0++) {
122  x = inputPoints(i0, 0);
123  y = inputPoints(i0, 1);
124 
125  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
126  outputValues(0, i0) = (x + y - 1.0)*(2.0*x + 2.0*y - 1.0);
127  outputValues(1, i0) = x*(2.0*x - 1.0);
128  outputValues(2, i0) = y*(2.0*y - 1.0);
129  outputValues(3, i0) = -4.0*x*(x + y - 1.0);
130  outputValues(4, i0) = 4.0*x*y;
131  outputValues(5, i0) = -4.0*y*(x + y - 1.0);
132 
133  }
134  break;
135 
136  case OPERATOR_GRAD:
137  case OPERATOR_D1:
138  for (int i0 = 0; i0 < dim0; i0++) {
139  x = inputPoints(i0, 0);
140  y = inputPoints(i0, 1);
141 
142  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
143  outputValues(0, i0, 0) = 4.0*x + 4.0*y - 3.0;
144  outputValues(0, i0, 1) = 4.0*x + 4.0*y - 3.0;
145 
146  outputValues(1, i0, 0) = 4.0*x - 1.0;
147  outputValues(1, i0, 1) = 0.0;
148 
149  outputValues(2, i0, 0) = 0.0;
150  outputValues(2, i0, 1) = 4.0*y - 1.0;
151 
152  outputValues(3, i0, 0) = -4.0*(2.0*x + y - 1.0);
153  outputValues(3, i0, 1) = -4.0*x;
154 
155  outputValues(4, i0, 0) = 4.0*y;
156  outputValues(4, i0, 1) = 4.0*x;
157 
158  outputValues(5, i0, 0) = -4.0*y;
159  outputValues(5, i0, 1) = -4.0*(x + 2.0*y - 1.0);
160  }
161  break;
162 
163  case OPERATOR_CURL:
164  for (int i0 = 0; i0 < dim0; i0++) {
165  x = inputPoints(i0, 0);
166  y = inputPoints(i0, 1);
167 
168  // CURL(u) = (u_y, -u_x), is rotated GRAD
169  outputValues(0, i0, 1) =-(4.0*x + 4.0*y - 3.0);
170  outputValues(0, i0, 0) = 4.0*x + 4.0*y - 3.0;
171 
172  outputValues(1, i0, 1) =-(4.0*x - 1.0);
173  outputValues(1, i0, 0) = 0.0;
174 
175  outputValues(2, i0, 1) = 0.0;
176  outputValues(2, i0, 0) = 4.0*y - 1.0;
177 
178  outputValues(3, i0, 1) = 4.0*(2.0*x + y - 1.0);
179  outputValues(3, i0, 0) = -4.0*x;
180 
181  outputValues(4, i0, 1) = -4.0*y;
182  outputValues(4, i0, 0) = 4.0*x;
183 
184  outputValues(5, i0, 1) = 4.0*y;
185  outputValues(5, i0, 0) = -4.0*(x + 2.0*y - 1.0);
186  }
187  break;
188 
189  case OPERATOR_DIV:
190  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
191  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): DIV is invalid operator for rank-0 (scalar) fields in 2D.");
192  break;
193 
194  case OPERATOR_D2:
195  for (int i0 = 0; i0 < dim0; i0++) {
196  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
197  // D2 -> (2,0) -> dx^2.
198  outputValues(0, i0, 0) = 4.0;
199  outputValues(1, i0, 0) = 4.0;
200  outputValues(2, i0, 0) = 0.0;
201  outputValues(3, i0, 0) =-8.0;
202  outputValues(4, i0, 0) = 0.0;
203  outputValues(5, i0, 0) = 0.0;
204 
205  // D2 -> (1,1) -> dx dy
206  outputValues(0, i0, 1) = 4.0;
207  outputValues(1, i0, 1) = 0.0;
208  outputValues(2, i0, 1) = 0.0;
209  outputValues(3, i0, 1) =-4.0;
210  outputValues(4, i0, 1) = 4.0;
211  outputValues(5, i0, 1) =-4.0;
212 
213  // D2 -> (0,2) -> dy^2
214  outputValues(0, i0, 2) = 4.0;
215  outputValues(1, i0, 2) = 0.0;
216  outputValues(2, i0, 2) = 4.0;
217  outputValues(3, i0, 2) = 0.0;
218  outputValues(4, i0, 2) = 0.0;
219  outputValues(5, i0, 2) =-8.0;
220  }// for i0
221  break;
222 
223  case OPERATOR_D3:
224  case OPERATOR_D4:
225  case OPERATOR_D5:
226  case OPERATOR_D6:
227  case OPERATOR_D7:
228  case OPERATOR_D8:
229  case OPERATOR_D9:
230  case OPERATOR_D10:
231  {
232  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
233  int DkCardinality = Intrepid::getDkCardinality(operatorType,
234  this -> basisCellTopology_.getDimension() );
235  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
236  for (int i0 = 0; i0 < dim0; i0++) {
237  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
238  outputValues(dofOrd, i0, dkOrd) = 0.0;
239  }
240  }
241  }
242  }
243  break;
244 
245  default:
246  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
247  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): Invalid operator type");
248  }
249 }
250 
251 
252 
253 template<class Scalar, class ArrayScalar>
255  const ArrayScalar & inputPoints,
256  const ArrayScalar & cellVertices,
257  const EOperator operatorType) const {
258  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
259  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): FEM Basis calling an FVD member function");
260 }
261 
262 template<class Scalar, class ArrayScalar>
264 #ifdef HAVE_INTREPID_DEBUG
265  // Verify rank of output array.
266  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
267  ">>> ERROR: (Intrepid::Basis_HGRAD_TRI_C2_FEM::getDofCoords) rank = 2 required for DofCoords array");
268  // Verify 0th dimension of output array.
269  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
270  ">>> ERROR: (Intrepid::Basis_HGRAD_TRI_C2_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
271  // Verify 1st dimension of output array.
272  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
273  ">>> ERROR: (Intrepid::Basis_HGRAD_TRI_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
274 #endif
275 
276  DofCoords(0,0) = 0.0; DofCoords(0,1) = 0.0;
277  DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0;
278  DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0;
279  DofCoords(3,0) = 0.5; DofCoords(3,1) = 0.0;
280  DofCoords(4,0) = 0.5; DofCoords(4,1) = 0.5;
281  DofCoords(5,0) = 0.0; DofCoords(5,1) = 0.5;
282 }
283 
284 }// namespace Intrepid
285 #endif
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
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.
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.