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 
50 template<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 
66 template<class Scalar,class ArrayScalar>
67 const 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 
77 template<class Scalar, class ArrayScalar>
78 const 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 
88 template<class Scalar, class ArrayScalar>
89 const 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 
99 template<class Scalar, class ArrayScalar>
101  return basisCardinality_;
102 }
103 
104 
105 template<class Scalar, class ArrayScalar>
107  return basisType_;
108 }
109 
110 
111 template<class Scalar, class ArrayScalar>
112 inline const shards::CellTopology Basis<Scalar, ArrayScalar>::getBaseCellTopology() const {
113  return basisCellTopology_;
114 }
115 
116 
117 template<class Scalar, class ArrayScalar>
119  return basisDegree_;
120 }
121 
122 
123 template<class Scalar, class ArrayScalar>
125  return basisCoordinates_;
126 }
127 
128 
129 //--------------------------------------------------------------------------------------------//
130 // //
131 // Helper functions of the Basis class //
132 // //
133 //--------------------------------------------------------------------------------------------//
134 
135 template<class Scalar, class ArrayScalar>
136 void 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 
255 template<class Scalar, class ArrayScalar>
256 void 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 
339 template<class Scalar, class ArrayScalar>
340 void 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.
402 template<class ArrayScalar>
404 
405 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
406 #ifdef __GNUC__
407 #warning "The Intrepid package is deprecated"
408 #endif
409 #endif
410 
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
An abstract base class that defines interface for concrete basis implementations for Finite Element (...