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>
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>
101 return basisCardinality_;
102}
103
104
105template<class Scalar, class ArrayScalar>
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>
119 return basisDegree_;
120}
121
122
123template<class Scalar, class ArrayScalar>
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.");
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}
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>
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....
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
This is an interface class for bases whose degrees of freedom can be associated with spatial location...