Intrepid
Intrepid_BasisDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
50template<class Scalar, class ArrayScalar>
51int Basis<Scalar, ArrayScalar>::getDofOrdinal(const int subcDim,
52 const int subcOrd,
53 const int subcDofOrd) {
54 if (!basisTagsAreSet_) {
55 initializeTags();
56 basisTagsAreSet_ = true;
57 }
58 // Use .at() for bounds checking
59 int dofOrdinal = tagToOrdinal_.at(subcDim).at(subcOrd).at(subcDofOrd);
60 TEUCHOS_TEST_FOR_EXCEPTION( (dofOrdinal == -1), std::invalid_argument,
61 ">>> ERROR (Basis): Invalid DoF tag");
62 return dofOrdinal;
63}
64
65
66template<class Scalar,class ArrayScalar>
67const std::vector<std::vector<std::vector<int> > > & Basis<Scalar, ArrayScalar>::getDofOrdinalData( )
68{
69 if (!basisTagsAreSet_) {
70 initializeTags();
71 basisTagsAreSet_ = true;
72 }
73 return tagToOrdinal_;
74}
75
76
77template<class Scalar, class ArrayScalar>
78const std::vector<int>& Basis<Scalar, ArrayScalar>::getDofTag(int dofOrd) {
79 if (!basisTagsAreSet_) {
80 initializeTags();
81 basisTagsAreSet_ = true;
82 }
83 // Use .at() for bounds checking
84 return ordinalToTag_.at(dofOrd);
85}
86
87
88template<class Scalar, class ArrayScalar>
89const std::vector<std::vector<int> > & Basis<Scalar, ArrayScalar>::getAllDofTags() {
90 if (!basisTagsAreSet_) {
91 initializeTags();
92 basisTagsAreSet_ = true;
93 }
94 return ordinalToTag_;
95}
96
97
98
99template<class Scalar, class ArrayScalar>
100inline int Basis<Scalar, ArrayScalar>::getCardinality() const {
101 return basisCardinality_;
102}
103
104
105template<class Scalar, class ArrayScalar>
106inline EBasis Basis<Scalar, ArrayScalar>::getBasisType() const {
107 return basisType_;
108}
109
110
111template<class Scalar, class ArrayScalar>
112inline const shards::CellTopology Basis<Scalar, ArrayScalar>::getBaseCellTopology() const {
113 return basisCellTopology_;
114}
115
116
117template<class Scalar, class ArrayScalar>
118inline int Basis<Scalar,ArrayScalar>::getDegree() const {
119 return basisDegree_;
120}
121
122
123template<class Scalar, class ArrayScalar>
124inline ECoordinates Basis<Scalar, ArrayScalar>::getCoordinateSystem() const {
125 return basisCoordinates_;
126}
127
128
129//--------------------------------------------------------------------------------------------//
130// //
131// Helper functions of the Basis class //
132// //
133//--------------------------------------------------------------------------------------------//
134
135template<class Scalar, class ArrayScalar>
136void getValues_HGRAD_Args(ArrayScalar & outputValues,
137 const ArrayScalar & inputPoints,
138 const EOperator operatorType,
139 const shards::CellTopology& cellTopo,
140 const int basisCard){
141
142 int spaceDim = cellTopo.getDimension();
143
144 // Verify inputPoints array
145 TEUCHOS_TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument,
146 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 2 required for inputPoints array");
147
148 TEUCHOS_TEST_FOR_EXCEPTION( (inputPoints.dimension(0) <= 0), std::invalid_argument,
149 ">>> ERROR (Intrepid::getValues_HGRAD_Args): dim 0 (number of points) > 0 required for inputPoints array");
150
151 TEUCHOS_TEST_FOR_EXCEPTION( !(inputPoints.dimension(1) == spaceDim), std::invalid_argument,
152 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
153
154
155 // Verify that all inputPoints are in the reference cell
156 /*
157 TEUCHOS_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
158 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) One or more points are outside the "
159 << cellTopo <<" reference cell");
160 */
161
162
163 // Verify that operatorType is admissible for HGRAD fields
164 TEUCHOS_TEST_FOR_EXCEPTION( ( (spaceDim == 2) && (operatorType == OPERATOR_DIV) ), std::invalid_argument,
165 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) DIV is invalid operator for rank-0 (scalar) fields in 2D.");
166
167 TEUCHOS_TEST_FOR_EXCEPTION( ( (spaceDim == 3) && ( (operatorType == OPERATOR_DIV) ||
168 (operatorType == OPERATOR_CURL) ) ), std::invalid_argument,
169 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) DIV and CURL are invalid operators for rank-0 (scalar) fields in 3D.");
170
171
172 // Check rank of outputValues (all operators are admissible in 1D) and its dim 2 when operator is
173 // GRAD, CURL (only in 2D), or Dk.
174
175 if(spaceDim == 1) {
176 switch(operatorType){
177 case OPERATOR_VALUE:
178 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
179 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
180 break;
181 case OPERATOR_GRAD:
182 case OPERATOR_CURL:
183 case OPERATOR_DIV:
184 case OPERATOR_D1:
185 case OPERATOR_D2:
186 case OPERATOR_D3:
187 case OPERATOR_D4:
188 case OPERATOR_D5:
189 case OPERATOR_D6:
190 case OPERATOR_D7:
191 case OPERATOR_D8:
192 case OPERATOR_D9:
193 case OPERATOR_D10:
194 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
195 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 3 required for outputValues in 1D when operator = GRAD, CURL, DIV, or Dk.");
196
197 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == 1 ),
198 std::invalid_argument,
199 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 2 of outputValues must equal 1 when operator = GRAD, CURL, DIV, or Dk.");
200 break;
201 default:
202 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HGRAD_Args) Invalid operator");
203 }
204 }
205 else if(spaceDim > 1) {
206 switch(operatorType){
207 case OPERATOR_VALUE:
208 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
209 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 2 required for outputValues when operator = VALUE.");
210 break;
211 case OPERATOR_GRAD:
212 case OPERATOR_CURL:
213 case OPERATOR_D1:
214 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
215 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
216
217 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim ),
218 std::invalid_argument,
219 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 2 of outputValues must equal cell dimension when operator = GRAD, CURL (in 2D), or D1.");
220 break;
221 case OPERATOR_D2:
222 case OPERATOR_D3:
223 case OPERATOR_D4:
224 case OPERATOR_D5:
225 case OPERATOR_D6:
226 case OPERATOR_D7:
227 case OPERATOR_D8:
228 case OPERATOR_D9:
229 case OPERATOR_D10:
230 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
231 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) rank = 3 required for outputValues in 2D and 3D when operator = GRAD, CURL (in 2D), or Dk.");
232
233 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == Intrepid::getDkCardinality(operatorType, spaceDim) ),
234 std::invalid_argument,
235 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 2 of outputValues must equal cardinality of the Dk multiset.");
236 break;
237 default:
238 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HGRAD_Args) Invalid operator");
239 }
240 }
241
242
243 // Verify dim 0 and dim 1 of outputValues
244 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(1) == inputPoints.dimension(0) ),
245 std::invalid_argument,
246 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
247
248 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(0) == basisCard ),
249 std::invalid_argument,
250 ">>> ERROR: (Intrepid::getValues_HGRAD_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
251}
252
253
254
255template<class Scalar, class ArrayScalar>
256void getValues_HCURL_Args(ArrayScalar & outputValues,
257 const ArrayScalar & inputPoints,
258 const EOperator operatorType,
259 const shards::CellTopology& cellTopo,
260 const int basisCard){
261
262 int spaceDim = cellTopo.getDimension();
263
264 // Verify that cell is 2D or 3D (this is redundant for default bases where we use correct cells,
265 // but will force any user-defined basis for HCURL spaces to use 2D or 3D cells
266 TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument,
267 ">>> ERROR: (Intrepid::getValues_HCURL_Args) cell dimension = 2 or 3 required for HCURL basis");
268
269
270 // Verify inputPoints array
271 TEUCHOS_TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument,
272 ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 2 required for inputPoints array");
273 TEUCHOS_TEST_FOR_EXCEPTION( (inputPoints.dimension(0) <= 0), std::invalid_argument,
274 ">>> ERROR (Intrepid::getValues_HCURL_Args): dim 0 (number of points) > 0 required for inputPoints array");
275
276 TEUCHOS_TEST_FOR_EXCEPTION( !(inputPoints.dimension(1) == spaceDim), std::invalid_argument,
277 ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
278
279 // Verify that all inputPoints are in the reference cell
280 /*
281 TEUCHOS_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
282 ">>> ERROR: (Intrepid::getValues_HCURL_Args) One or more points are outside the "
283 << cellTopo <<" reference cell");
284 */
285
286 // Verify that operatorType is admissible for HCURL fields
287 TEUCHOS_TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_CURL) ), std::invalid_argument,
288 ">>> ERROR: (Intrepid::getValues_HCURL_Args) operator = VALUE or CURL required for HCURL fields.");
289
290
291 // Check rank of outputValues: for VALUE should always be rank-3 array with (F,P,D) layout
292 switch(operatorType) {
293
294 case OPERATOR_VALUE:
295 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
296 ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 3 required for outputValues when operator is VALUE");
297 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim ),
298 std::invalid_argument,
299 ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension when operator is VALUE.");
300 break;
301
302 case OPERATOR_CURL:
303
304 // in 3D we need an (F,P,D) container because CURL gives a vector field:
305 if(spaceDim == 3) {
306 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3 ) ,
307 std::invalid_argument,
308 ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 3 required for outputValues in 3D when operator is CURL");
309 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim),
310 std::invalid_argument,
311 ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 2 of outputValues must equal cell dimension in 3D when operator is CURL.");
312 }
313 // In 2D we need an (F,P) container because CURL gives a scalar field
314 else if(spaceDim == 2) {
315 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2 ) ,
316 std::invalid_argument,
317 ">>> ERROR: (Intrepid::getValues_HCURL_Args) rank = 2 required for outputValues in 2D when operator is CURL");
318 }
319 break;
320
321 default:
322 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HCURL_Args) Invalid operator");
323 }
324
325
326 // Verify dim 0 and dim 1 of outputValues
327 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(1) == inputPoints.dimension(0) ),
328 std::invalid_argument,
329 ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
330
331 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(0) == basisCard ),
332 std::invalid_argument,
333 ">>> ERROR: (Intrepid::getValues_HCURL_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
334
335}
336
337
338
339template<class Scalar, class ArrayScalar>
340void getValues_HDIV_Args(ArrayScalar & outputValues,
341 const ArrayScalar & inputPoints,
342 const EOperator operatorType,
343 const shards::CellTopology& cellTopo,
344 const int basisCard){
345
346 int spaceDim = cellTopo.getDimension();
347
348 // Verify inputPoints array
349 TEUCHOS_TEST_FOR_EXCEPTION( !(inputPoints.rank() == 2), std::invalid_argument,
350 ">>> ERROR: (Intrepid::getValues_HDIV_Args) rank = 2 required for inputPoints array");
351 TEUCHOS_TEST_FOR_EXCEPTION( (inputPoints.dimension(0) <= 0), std::invalid_argument,
352 ">>> ERROR (Intrepid::getValues_HDIV_Args): dim 0 (number of points) > 0 required for inputPoints array");
353
354 TEUCHOS_TEST_FOR_EXCEPTION( !(inputPoints.dimension(1) == spaceDim), std::invalid_argument,
355 ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 1 (spatial dimension) of inputPoints array does not match cell dimension");
356
357 // Verify that all inputPoints are in the reference cell
358 /*
359 TEUCHOS_TEST_FOR_EXCEPTION( !CellTools<Scalar>::checkPointSetInclusion(inputPoints, cellTopo), std::invalid_argument,
360 ">>> ERROR: (Intrepid::getValues_HDIV_Args) One or more points are outside the "
361 << cellTopo <<" reference cell");
362 */
363
364 // Verify that operatorType is admissible for HDIV fields
365 TEUCHOS_TEST_FOR_EXCEPTION( !( (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_DIV) ), std::invalid_argument,
366 ">>> ERROR: (Intrepid::getValues_HDIV_Args) operator = VALUE or DIV required for HDIV fields.");
367
368
369 // Check rank of outputValues
370 switch(operatorType) {
371 case OPERATOR_VALUE:
372 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 3), std::invalid_argument,
373 ">>> ERROR: (Intrepid::getValues_HDIV_Args) rank = 3 required for outputValues when operator is VALUE.");
374
375 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(2) == spaceDim ),
376 std::invalid_argument,
377 ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 2 of outputValues must equal cell dimension for operator VALUE.");
378 break;
379 case OPERATOR_DIV:
380 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.rank() == 2), std::invalid_argument,
381 ">>> ERROR: (Intrepid::getValues_HDIV_Args) rank = 2 required for outputValues when operator is DIV.");
382 break;
383
384 default:
385 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument, ">>> ERROR: (Intrepid::getValues_HDIV_Args) Invalid operator");
386 }
387
388
389 // Verify dim 0 and dim 1 of outputValues
390 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(1) == inputPoints.dimension(0) ),
391 std::invalid_argument,
392 ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 1 (number of points) of outputValues must equal dim 0 of inputPoints.");
393
394 TEUCHOS_TEST_FOR_EXCEPTION( !(outputValues.dimension(0) == basisCard ),
395 std::invalid_argument,
396 ">>> ERROR: (Intrepid::getValues_HDIV_Args) dim 0 (number of basis functions) of outputValues must equal basis cardinality.");
397}
398
399// Pure virtual destructor (gives warnings if not included).
400// Following "Effective C++: 3rd Ed." item 7 the implementation
401// is included in the definition file.
402template<class ArrayScalar>
403DofCoordsInterface<ArrayScalar>::~DofCoordsInterface() {}
void getValues_HDIV_Args(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType, const shards::CellTopology &cellTopo, const int basisCard)
Runtime check of the arguments for the getValues method in an HDIV-conforming FEM basis....
void getValues_HCURL_Args(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType, const shards::CellTopology &cellTopo, const int basisCard)
Runtime check of the arguments for the getValues method in an HCURL-conforming FEM basis....
void getValues_HGRAD_Args(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType, const shards::CellTopology &cellTopo, const int basisCard)
Runtime check of the arguments for the getValues method in an HGRAD-conforming FEM basis....