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>
51 int 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 
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>
100 inline int Basis<Scalar, ArrayScalar>::getCardinality() const {
101  return basisCardinality_;
102 }
103 
104 
105 template<class Scalar, class ArrayScalar>
106 inline EBasis Basis<Scalar, ArrayScalar>::getBasisType() const {
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>
118 inline int Basis<Scalar,ArrayScalar>::getDegree() const {
119  return basisDegree_;
120 }
121 
122 
123 template<class Scalar, class ArrayScalar>
124 inline ECoordinates Basis<Scalar, ArrayScalar>::getCoordinateSystem() const {
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>
403 DofCoordsInterface<ArrayScalar>::~DofCoordsInterface() {}
Intrepid::ECoordinates
ECoordinates
Enumeration of coordinate systems for geometrical entities (cells, points).
Definition: Intrepid_Types.hpp:127
Intrepid::EBasis
EBasis
Enumeration of basis types for discrete spaces in Intrepid.
Definition: Intrepid_Types.hpp:417
Intrepid::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition: Intrepid_Basis.hpp:89
Intrepid::getValues_HDIV_Args
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....
Intrepid::getDkCardinality
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Definition: Intrepid_Utils.cpp:400
Intrepid::getValues_HCURL_Args
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....
Intrepid::EOperator
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
Definition: Intrepid_Types.hpp:206
Intrepid::getValues_HGRAD_Args
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....
Intrepid::DofCoordsInterface
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
Definition: Intrepid_Basis.hpp:353
Generated by   doxygen 1.8.16