Intrepid
Intrepid_CellToolsDef.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 
44 
49 #ifndef INTREPID_CELLTOOLSDEF_HPP
50 #define INTREPID_CELLTOOLSDEF_HPP
51 
52 // disable clang warnings
53 #if defined (__clang__) && !defined (__INTEL_COMPILER)
54 #pragma clang system_header
55 #endif
56 
57 namespace Intrepid {
58 
59  template<class Scalar>
60  const FieldContainer<double>& CellTools<Scalar>::getSubcellParametrization(const int subcellDim,
61  const shards::CellTopology& parentCell){
62 
63 #ifdef HAVE_INTREPID_DEBUG
64  TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument,
65  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): the specified cell topology does not have a reference cell.");
66 
67  TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= subcellDim) && (subcellDim <= 2 ) ), std::invalid_argument,
68  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): parametrization defined only for 1 and 2-dimensional subcells.");
69 #endif
70 
71  // Coefficients of the coordinate functions defining the parametrization maps are stored in
72  // rank-3 arrays with dimensions (SC, PCD, COEF) where:
73  // - SC is the subcell count of subcells with the specified dimension in the parent cell
74  // - PCD is Parent Cell Dimension, which gives the number of coordinate functions in the map:
75  // PCD = 2 for standard 2D cells and non-standard 2D cells: shell line and beam
76  // PCD = 3 for standard 3D cells and non-standard 3D cells: shell Tri and Quad
77  // - COEF is number of coefficients needed to specify a coordinate function:
78  // COEFF = 2 for edge parametrizations
79  // COEFF = 3 for both Quad and Tri face parametrizations. Because all Quad reference faces
80  // are affine, the coefficient of the bilinear term u*v is zero and is not stored, i.e.,
81  // 3 coefficients are sufficient to store Quad face parameterization maps.
82  //
83  // Arrays are sized and filled only when parametrization of a particular subcell is requested
84  // by setSubcellParametrization.
85 
86  // Edge maps for 2D non-standard cells: ShellLine and Beam
87  static FieldContainer<double> lineEdges; static int lineEdgesSet = 0;
88 
89  // Edge maps for 2D standard cells: Triangle and Quadrilateral
90  static FieldContainer<double> triEdges; static int triEdgesSet = 0;
91  static FieldContainer<double> quadEdges; static int quadEdgesSet = 0;
92 
93  // Edge maps for 3D non-standard cells: Shell Tri and Quad
94  static FieldContainer<double> shellTriEdges; static int shellTriEdgesSet = 0;
95  static FieldContainer<double> shellQuadEdges; static int shellQuadEdgesSet = 0;
96 
97  // Edge maps for 3D standard cells:
98  static FieldContainer<double> tetEdges; static int tetEdgesSet = 0;
99  static FieldContainer<double> hexEdges; static int hexEdgesSet = 0;
100  static FieldContainer<double> pyrEdges; static int pyrEdgesSet = 0;
101  static FieldContainer<double> wedgeEdges; static int wedgeEdgesSet = 0;
102 
103 
104  // Face maps for 3D non-standard cells: Shell Triangle and Quadrilateral
105  static FieldContainer<double> shellTriFaces; static int shellTriFacesSet = 0;
106  static FieldContainer<double> shellQuadFaces; static int shellQuadFacesSet = 0;
107 
108  // Face maps for 3D standard cells:
109  static FieldContainer<double> tetFaces; static int tetFacesSet = 0;
110  static FieldContainer<double> hexFaces; static int hexFacesSet = 0;
111  static FieldContainer<double> pyrFaces; static int pyrFacesSet = 0;
112  static FieldContainer<double> wedgeFaces; static int wedgeFacesSet = 0;
113 
114  // Select subcell parametrization according to its parent cell type
115  switch(parentCell.getKey() ) {
116 
117  // Tet cells
118  case shards::Tetrahedron<4>::key:
119  case shards::Tetrahedron<8>::key:
120  case shards::Tetrahedron<10>::key:
121  case shards::Tetrahedron<11>::key:
122  if(subcellDim == 2) {
123  if(!tetFacesSet){
124  setSubcellParametrization(tetFaces, subcellDim, parentCell);
125  tetFacesSet = 1;
126  }
127  return tetFaces;
128  }
129  else if(subcellDim == 1) {
130  if(!tetEdgesSet){
131  setSubcellParametrization(tetEdges, subcellDim, parentCell);
132  tetEdgesSet = 1;
133  }
134  return tetEdges;
135  }
136  else{
137  TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument,
138  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Tet parametrizations defined for 1 and 2-subcells only");
139  }
140  break;
141 
142  // Hex cells
143  case shards::Hexahedron<8>::key:
144  case shards::Hexahedron<20>::key:
145  case shards::Hexahedron<27>::key:
146  if(subcellDim == 2) {
147  if(!hexFacesSet){
148  setSubcellParametrization(hexFaces, subcellDim, parentCell);
149  hexFacesSet = 1;
150  }
151  return hexFaces;
152  }
153  else if(subcellDim == 1) {
154  if(!hexEdgesSet){
155  setSubcellParametrization(hexEdges, subcellDim, parentCell);
156  hexEdgesSet = 1;
157  }
158  return hexEdges;
159  }
160  else{
161  TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument,
162  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Hex parametrizations defined for 1 and 2-subcells only");
163  }
164  break;
165 
166  // Pyramid cells
167  case shards::Pyramid<5>::key:
168  case shards::Pyramid<13>::key:
169  case shards::Pyramid<14>::key:
170  if(subcellDim == 2) {
171  if(!pyrFacesSet){
172  setSubcellParametrization(pyrFaces, subcellDim, parentCell);
173  pyrFacesSet = 1;
174  }
175  return pyrFaces;
176  }
177  else if(subcellDim == 1) {
178  if(!pyrEdgesSet){
179  setSubcellParametrization(pyrEdges, subcellDim, parentCell);
180  pyrEdgesSet = 1;
181  }
182  return pyrEdges;
183  }
184  else {
185  TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument,
186  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Pyramid parametrizations defined for 1 and 2-subcells only");
187  }
188  break;
189 
190  // Wedge cells
191  case shards::Wedge<6>::key:
192  case shards::Wedge<15>::key:
193  case shards::Wedge<18>::key:
194  if(subcellDim == 2) {
195  if(!wedgeFacesSet){
196  setSubcellParametrization(wedgeFaces, subcellDim, parentCell);
197  wedgeFacesSet = 1;
198  }
199  return wedgeFaces;
200  }
201  else if(subcellDim == 1) {
202  if(!wedgeEdgesSet){
203  setSubcellParametrization(wedgeEdges, subcellDim, parentCell);
204  wedgeEdgesSet = 1;
205  }
206  return wedgeEdges;
207  }
208  else {
209  TEUCHOS_TEST_FOR_EXCEPTION( (subcellDim != 1 || subcellDim != 2), std::invalid_argument,
210  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Wedge parametrization defined for 1 and 2-subcells only");
211  }
212  break;
213  //
214  // Standard 2D cells have only 1-subcells
215  //
216  case shards::Triangle<3>::key:
217  case shards::Triangle<4>::key:
218  case shards::Triangle<6>::key:
219  if(subcellDim == 1) {
220  if(!triEdgesSet){
221  setSubcellParametrization(triEdges, subcellDim, parentCell);
222  triEdgesSet = 1;
223  }
224  return triEdges;
225  }
226  else{
227  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
228  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Triangle parametrizations defined for 1-subcells only");
229  }
230  break;
231 
232  case shards::Quadrilateral<4>::key:
233  case shards::Quadrilateral<8>::key:
234  case shards::Quadrilateral<9>::key:
235  if(subcellDim == 1) {
236  if(!quadEdgesSet){
237  setSubcellParametrization(quadEdges, subcellDim, parentCell);
238  quadEdgesSet = 1;
239  }
240  return quadEdges;
241  }
242  else{
243  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
244  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Quad parametrizations defined for 1-subcells only");
245  }
246  break;
247  //
248  // Non-standard 3D Shell Tri and Quad cells have 1 and 2-subcells. Because they are 3D cells
249  // can't reuse edge parametrization arrays for 2D Triangle and Quadrilateral.
250  //
251  case shards::ShellTriangle<3>::key:
252  case shards::ShellTriangle<6>::key:
253  if(subcellDim == 2) {
254  if(!shellTriFacesSet){
255  setSubcellParametrization(shellTriFaces, subcellDim, parentCell);
256  shellTriFacesSet = 1;
257  }
258  return shellTriFaces;
259  }
260  else if(subcellDim == 1) {
261  if(!shellTriEdgesSet){
262  setSubcellParametrization(shellTriEdges, subcellDim, parentCell);
263  shellTriEdgesSet = 1;
264  }
265  return shellTriEdges;
266  }
267  else if( subcellDim != 1 || subcellDim != 2){
268  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
269  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Shell Triangle parametrizations defined for 1 and 2-subcells only");
270  }
271  break;
272 
273  case shards::ShellQuadrilateral<4>::key:
274  case shards::ShellQuadrilateral<8>::key:
275  case shards::ShellQuadrilateral<9>::key:
276  if(subcellDim == 2) {
277  if(!shellQuadFacesSet){
278  setSubcellParametrization(shellQuadFaces, subcellDim, parentCell);
279  shellQuadFacesSet = 1;
280  }
281  return shellQuadFaces;
282  }
283  else if(subcellDim == 1) {
284  if(!shellQuadEdgesSet){
285  setSubcellParametrization(shellQuadEdges, subcellDim, parentCell);
286  shellQuadEdgesSet = 1;
287  }
288  return shellQuadEdges;
289  }
290  else if( subcellDim != 1 || subcellDim != 2){
291  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
292  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): Shell Quad parametrizations defined for 1 and 2-subcells only");
293  }
294  break;
295 
296  //
297  // Non-standard 2D cells: Shell Lines and Beams have 1-subcells
298  //
299  case shards::ShellLine<2>::key:
300  case shards::ShellLine<3>::key:
301  case shards::Beam<2>::key:
302  case shards::Beam<3>::key:
303  if(subcellDim == 1) {
304  if(!lineEdgesSet){
305  setSubcellParametrization(lineEdges, subcellDim, parentCell);
306  lineEdgesSet = 1;
307  }
308  return lineEdges;
309  }
310  else{
311  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
312  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): shell line/beam parametrizations defined for 1-subcells only");
313  }
314  break;
315 
316  default:
317  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
318  ">>> ERROR (Intrepid::CellTools::getSubcellParametrization): invalid cell topology.");
319  }//cell key
320  // To disable compiler warning, should never be reached
321  return lineEdges;
322  }
323 
324 
325 
326  template<class Scalar>
327  void CellTools<Scalar>::setSubcellParametrization(FieldContainer<double>& subcellParametrization,
328  const int subcellDim,
329  const shards::CellTopology& parentCell)
330  {
331 #ifdef HAVE_INTREPID_DEBUG
332  TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument,
333  ">>> ERROR (Intrepid::CellTools::setSubcellParametrization): the specified cell topology does not have a reference cell.");
334 
335  TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= subcellDim) && (subcellDim <= 2 ) ), std::invalid_argument,
336  ">>> ERROR (Intrepid::CellTools::setSubcellParametrization): parametrization defined only for 1 and 2-dimensional subcells.");
337 #endif
338  // subcellParametrization is rank-3 FieldContainer with dimensions (SC, PCD, COEF) where:
339  // - SC is the subcell count of subcells with the specified dimension in the parent cell
340  // - PCD is Parent Cell Dimension, which gives the number of coordinate functions in the map
341  // PCD = 2 for standard 2D cells and non-standard 2D cells: shell line and beam
342  // PCD = 3 for standard 3D cells and non-standard 3D cells: shell Tri and Quad
343  // - COEF is number of coefficients needed to specify a coordinate function:
344  // COEFF = 2 for edge parametrizations
345  // COEFF = 3 for both Quad and Tri face parametrizations. Because all Quad reference faces
346  // are affine, the coefficient of the bilinear term u*v is zero and is not stored, i.e.,
347  // 3 coefficients are sufficient to store Quad face parameterization maps.
348  //
349  // Edge parametrization maps [-1,1] to edge defined by (v0, v1)
350  // Face parametrization maps [-1,1]^2 to quadrilateral face (v0, v1, v2, v3), or
351  // standard 2-simplex {(0,0),(1,0),(0,1)} to traingle face (v0, v1, v2).
352  // This defines orientation-preserving parametrizations with respect to reference edge and
353  // face orientations induced by their vertex order.
354 
355  // get subcellParametrization dimensions: (sc, pcd, coeff)
356  unsigned sc = parentCell.getSubcellCount(subcellDim);
357  unsigned pcd = parentCell.getDimension();
358  unsigned coeff = (subcellDim == 1) ? 2 : 3;
359 
360 
361  // Resize container
362  subcellParametrization.resize(sc, pcd, coeff);
363 
364  // Edge parametrizations of 2D and 3D cells (shell lines and beams are 2D cells with edges)
365  if(subcellDim == 1){
366  for(unsigned subcellOrd = 0; subcellOrd < sc; subcellOrd++){
367 
368  int v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
369  int v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
370 
371  // vertexK[0] = x_k; vertexK[1] = y_k; vertexK[2] = z_k; z_k = 0 for 2D cells
372  // Note that ShellLine and Beam are 2D cells!
373  const double* v0 = getReferenceVertex(parentCell, v0ord);
374  const double* v1 = getReferenceVertex(parentCell, v1ord);
375 
376  // x(t) = (x0 + x1)/2 + t*(x1 - x0)/2
377  subcellParametrization(subcellOrd, 0, 0) = (v0[0] + v1[0])/2.0;
378  subcellParametrization(subcellOrd, 0, 1) = (v1[0] - v0[0])/2.0;
379 
380  // y(t) = (y0 + y1)/2 + t*(y1 - y0)/2
381  subcellParametrization(subcellOrd, 1, 0) = (v0[1] + v1[1])/2.0;
382  subcellParametrization(subcellOrd, 1, 1) = (v1[1] - v0[1])/2.0;
383 
384  if( pcd == 3 ) {
385  // z(t) = (z0 + z1)/2 + t*(z1 - z0)/2
386  subcellParametrization(subcellOrd, 2, 0) = (v0[2] + v1[2])/2.0;
387  subcellParametrization(subcellOrd, 2, 1) = (v1[2] - v0[2])/2.0;
388  }
389  }// for loop over 1-subcells
390  }
391 
392  // Face parametrizations of 3D cells: (shell Tri and Quad are 3D cells with faces)
393  // A 3D cell can have both Tri and Quad faces, but because they are affine images of the
394  // parametrization domain, 3 coefficients are enough to store them in both cases.
395  else if(subcellDim == 2) {
396  for(unsigned subcellOrd = 0; subcellOrd < sc; subcellOrd++){
397 
398  switch(parentCell.getKey(subcellDim,subcellOrd)){
399 
400  // Admissible triangular faces for 3D cells in Shards:
401  case shards::Triangle<3>::key:
402  case shards::Triangle<4>::key:
403  case shards::Triangle<6>::key:
404  {
405  int v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
406  int v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
407  int v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
408  const double* v0 = getReferenceVertex(parentCell, v0ord);
409  const double* v1 = getReferenceVertex(parentCell, v1ord);
410  const double* v2 = getReferenceVertex(parentCell, v2ord);
411 
412  // x(u,v) = x0 + (x1 - x0)*u + (x2 - x0)*v
413  subcellParametrization(subcellOrd, 0, 0) = v0[0];
414  subcellParametrization(subcellOrd, 0, 1) = v1[0] - v0[0];
415  subcellParametrization(subcellOrd, 0, 2) = v2[0] - v0[0];
416 
417  // y(u,v) = y0 + (y1 - y0)*u + (y2 - y0)*v
418  subcellParametrization(subcellOrd, 1, 0) = v0[1];
419  subcellParametrization(subcellOrd, 1, 1) = v1[1] - v0[1];
420  subcellParametrization(subcellOrd, 1, 2) = v2[1] - v0[1];
421 
422  // z(u,v) = z0 + (z1 - z0)*u + (z2 - z0)*v
423  subcellParametrization(subcellOrd, 2, 0) = v0[2];
424  subcellParametrization(subcellOrd, 2, 1) = v1[2] - v0[2];
425  subcellParametrization(subcellOrd, 2, 2) = v2[2] - v0[2];
426  }
427  break;
428 
429  // Admissible quadrilateral faces for 3D cells in Shards:
430  case shards::Quadrilateral<4>::key:
431  case shards::Quadrilateral<8>::key:
432  case shards::Quadrilateral<9>::key:
433  {
434  int v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
435  int v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
436  int v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
437  int v3ord = parentCell.getNodeMap(subcellDim, subcellOrd, 3);
438  const double* v0 = getReferenceVertex(parentCell, v0ord);
439  const double* v1 = getReferenceVertex(parentCell, v1ord);
440  const double* v2 = getReferenceVertex(parentCell, v2ord);
441  const double* v3 = getReferenceVertex(parentCell, v3ord);
442 
443  // x(u,v) = (x0+x1+x2+x3)/4+u*(-x0+x1+x2-x3)/4+v*(-x0-x1+x2+x3)/4+uv*(0=x0-x1+x2-x3)/4
444  subcellParametrization(subcellOrd, 0, 0) = ( v0[0] + v1[0] + v2[0] + v3[0])/4.0;
445  subcellParametrization(subcellOrd, 0, 1) = (-v0[0] + v1[0] + v2[0] - v3[0])/4.0;
446  subcellParametrization(subcellOrd, 0, 2) = (-v0[0] - v1[0] + v2[0] + v3[0])/4.0;
447 
448  // y(u,v) = (y0+y1+y2+y3)/4+u*(-y0+y1+y2-y3)/4+v*(-y0-y1+y2+y3)/4+uv*(0=y0-y1+y2-y3)/4
449  subcellParametrization(subcellOrd, 1, 0) = ( v0[1] + v1[1] + v2[1] + v3[1])/4.0;
450  subcellParametrization(subcellOrd, 1, 1) = (-v0[1] + v1[1] + v2[1] - v3[1])/4.0;
451  subcellParametrization(subcellOrd, 1, 2) = (-v0[1] - v1[1] + v2[1] + v3[1])/4.0;
452 
453  // z(u,v) = (z0+z1+z2+z3)/4+u*(-z0+z1+z2-z3)/4+v*(-z0-z1+z2+z3)/4+uv*(0=z0-z1+z2-z3)/4
454  subcellParametrization(subcellOrd, 2, 0) = ( v0[2] + v1[2] + v2[2] + v3[2])/4.0;
455  subcellParametrization(subcellOrd, 2, 1) = (-v0[2] + v1[2] + v2[2] - v3[2])/4.0;
456  subcellParametrization(subcellOrd, 2, 2) = (-v0[2] - v1[2] + v2[2] + v3[2])/4.0;
457  }
458  break;
459  default:
460  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
461  ">>> ERROR (Intrepid::CellTools::setSubcellParametrization): parametrization not defined for the specified face topology.");
462 
463  }// switch face topology key
464  }// for subcellOrd
465  }
466  }
467 
468 
469 
470  template<class Scalar>
471  const double* CellTools<Scalar>::getReferenceVertex(const shards::CellTopology& cell,
472  const int vertexOrd){
473 
474 #ifdef HAVE_INTREPID_DEBUG
475  TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(cell) ), std::invalid_argument,
476  ">>> ERROR (Intrepid::CellTools::getReferenceVertex): the specified cell topology does not have a reference cell.");
477 
478  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= vertexOrd) && (vertexOrd < (int)cell.getVertexCount() ) ), std::invalid_argument,
479  ">>> ERROR (Intrepid::CellTools::getReferenceVertex): invalid node ordinal for the specified cell topology. ");
480 #endif
481 
482  // Simply call getReferenceNode with the base topology of the cell
483  return getReferenceNode(cell.getBaseCellTopologyData(), vertexOrd);
484  }
485 
486 
487 
488  template<class Scalar>
489  template<class ArraySubcellVert>
490  void CellTools<Scalar>::getReferenceSubcellVertices(ArraySubcellVert & subcellVertices,
491  const int subcellDim,
492  const int subcellOrd,
493  const shards::CellTopology& parentCell){
494 #ifdef HAVE_INTREPID_DEBUG
495  TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument,
496  ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices): the specified cell topology does not have a reference cell.");
497 
498  // subcellDim can equal the cell dimension because the cell itself is a valid subcell! In this case
499  // the method will return all cell cellWorkset.
500  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellDim) && (subcellDim <= (int)parentCell.getDimension()) ), std::invalid_argument,
501  ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices): subcell dimension out of range.");
502 
503  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellOrd) && (subcellOrd < (int)parentCell.getSubcellCount(subcellDim) ) ), std::invalid_argument,
504  ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices): subcell ordinal out of range.");
505 
506  // Verify subcellVertices rank and dimensions
507  {
508  std::string errmsg = ">>> ERROR (Intrepid::CellTools::getReferenceSubcellVertices):";
509  TEUCHOS_TEST_FOR_EXCEPTION( !( requireRankRange(errmsg, subcellVertices, 2, 2) ), std::invalid_argument, errmsg);
510 
511  int subcVertexCount = parentCell.getVertexCount(subcellDim, subcellOrd);
512  int spaceDim = parentCell.getDimension();
513 
514  TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellVertices, 0, subcVertexCount, subcVertexCount) ),
515  std::invalid_argument, errmsg);
516 
517  TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellVertices, 1, spaceDim, spaceDim) ),
518  std::invalid_argument, errmsg);
519  }
520 #endif
521 
522  // Simply call getReferenceNodes with the base topology
523  getReferenceSubcellNodes(subcellVertices, subcellDim, subcellOrd, parentCell.getBaseCellTopologyData() );
524  }
525 
526 
527 
528  template<class Scalar>
529  const double* CellTools<Scalar>::getReferenceNode(const shards::CellTopology& cell,
530  const int nodeOrd){
531 
532 #ifdef HAVE_INTREPID_DEBUG
533  TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(cell) ), std::invalid_argument,
534  ">>> ERROR (Intrepid::CellTools::getReferenceNode): the specified cell topology does not have a reference cell.");
535 
536  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= nodeOrd) && (nodeOrd < (int)cell.getNodeCount() ) ), std::invalid_argument,
537  ">>> ERROR (Intrepid::CellTools::getReferenceNode): invalid node ordinal for the specified cell topology. ");
538 #endif
539 
540  // Cartesian coordinates of supported reference cell cellWorkset, padded to three-dimensions.
541  // Node order follows cell topology definition in Shards
542  static const double line[2][3] ={
543  {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}
544  };
545  static const double line_3[3][3] = {
546  {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0},
547  // Extension node: edge midpoint
548  { 0.0, 0.0, 0.0}
549  };
550 
551 
552  // Triangle topologies
553  static const double triangle[3][3] = {
554  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}
555  };
556  static const double triangle_4[4][3] = {
557  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
558  // Extension node: cell center
559  { 1/3, 1/3, 0.0}
560  };
561  static const double triangle_6[6][3] = {
562  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
563  // Extension cellWorkset: 3 edge midpoints
564  { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}
565  };
566 
567 
568  // Quadrilateral topologies
569  static const double quadrilateral[4][3] = {
570  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}
571  };
572  static const double quadrilateral_8[8][3] = {
573  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
574  // Extension cellWorkset: 4 edge midpoints
575  { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}
576  };
577  static const double quadrilateral_9[9][3] = {
578  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
579  // Extension cellWorkset: 4 edge midpoints + 1 cell center
580  { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}, { 0.0, 0.0, 0.0}
581  };
582 
583 
584  // Tetrahedron topologies
585  static const double tetrahedron[4][3] = {
586  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0}
587  };
588  static const double tetrahedron_8[8][3] = {
589  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
590  // Extension cellWorkset: 4 face centers (do not follow natural face order - see the cell topology!)
591  { 1/3, 0.0, 1/3}, { 1/3, 1/3, 1/3}, { 1/3, 1/3, 0.0}, { 0.0, 1/3, 1/3}
592  };
593  static const double tetrahedron_10[10][3] = {
594  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
595  // Extension cellWorkset: 6 edge midpoints
596  { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}, { 0.0, 0.0, 0.5}, { 0.5, 0.0, 0.5}, { 0.0, 0.5, 0.5}
597  };
598 
599  static const double tetrahedron_11[10][3] = {
600  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
601  // Extension cellWorkset: 6 edge midpoints
602  { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}, { 0.0, 0.0, 0.5}, { 0.5, 0.0, 0.5}, { 0.0, 0.5, 0.5}
603  };
604 
605 
606  // Hexahedron topologies
607  static const double hexahedron[8][3] = {
608  {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
609  {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0}
610  };
611  static const double hexahedron_20[20][3] = {
612  {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
613  {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0},
614  // Extension cellWorkset: 12 edge midpoints (do not follow natural edge order - see cell topology!)
615  { 0.0,-1.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, {-1.0, 0.0,-1.0},
616  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
617  { 0.0,-1.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}, {-1.0, 0.0, 1.0}
618  };
619  static const double hexahedron_27[27][3] = {
620  {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
621  {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0},
622  // Extension cellWorkset: 12 edge midpoints + 1 cell center + 6 face centers (do not follow natural subcell order!)
623  { 0.0,-1.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, {-1.0, 0.0,-1.0},
624  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
625  { 0.0,-1.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}, {-1.0, 0.0, 1.0},
626  { 0.0, 0.0, 0.0},
627  { 0.0, 0.0,-1.0}, { 0.0, 0.0, 1.0}, {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, {0.0,-1.0, 0.0}, {0.0, 1.0, 0.0}
628  };
629 
630 
631  // Pyramid topologies
632  static const double pyramid[5][3] = {
633  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0}
634  };
635  static const double pyramid_13[13][3] = {
636  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
637  // Extension cellWorkset: 8 edge midpoints
638  { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0},
639  {-0.5,-0.5, 0.5}, { 0.5,-0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5}
640  };
641  static const double pyramid_14[14][3] = {
642  {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
643  // Extension cellWorkset: 8 edge midpoints + quadrilateral face midpoint
644  { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0},
645  {-0.5,-0.5, 0.5}, { 0.5,-0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5}, { 0.0, 0.0, 0.0}
646  };
647 
648 
649  // Wedge topologies
650  static const double wedge[6][3] = {
651  { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}
652  };
653  static const double wedge_15[15][3] = {
654  { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0},
655  // Extension cellWorkset: 9 edge midpoints (do not follow natural edge order - see cell topology!)
656  { 0.5, 0.0,-1.0}, { 0.5, 0.5,-1.0}, { 0.0, 0.5,-1.0}, { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
657  { 0.5, 0.0, 1.0}, { 0.5, 0.5, 1.0}, { 0.0, 0.5, 1.0}
658  };
659  static const double wedge_18[18][3] = {
660  { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0},
661  // Extension cellWorkset: 9 edge midpoints + 3 quad face centers (do not follow natural subcell order - see cell topology!)
662  { 0.5, 0.0,-1.0}, { 0.5, 0.5,-1.0}, { 0.0, 0.5,-1.0}, { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
663  { 0.5, 0.0, 1.0}, { 0.5, 0.5, 1.0}, { 0.0, 0.5, 1.0},
664  { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}
665  };
666 
667 
668  switch(cell.getKey() ) {
669 
670  // Base line topologies
671  case shards::Line<2>::key:
672  case shards::ShellLine<2>::key:
673  case shards::Beam<2>::key:
674  return line[nodeOrd];
675  break;
676 
677  // Extended line topologies
678  case shards::Line<3>::key:
679  case shards::ShellLine<3>::key:
680  case shards::Beam<3>::key:
681  return line_3[nodeOrd];
682  break;
683 
684 
685  // Base triangle topologies
686  case shards::Triangle<3>::key:
687  case shards::ShellTriangle<3>::key:
688  return triangle[nodeOrd];
689  break;
690 
691  // Extened Triangle topologies
692  case shards::Triangle<4>::key:
693  return triangle_4[nodeOrd];
694  break;
695  case shards::Triangle<6>::key:
696  case shards::ShellTriangle<6>::key:
697  return triangle_6[nodeOrd];
698  break;
699 
700 
701  // Base Quadrilateral topologies
702  case shards::Quadrilateral<4>::key:
703  case shards::ShellQuadrilateral<4>::key:
704  return quadrilateral[nodeOrd];
705  break;
706 
707  // Extended Quadrilateral topologies
708  case shards::Quadrilateral<8>::key:
709  case shards::ShellQuadrilateral<8>::key:
710  return quadrilateral_8[nodeOrd];
711  break;
712  case shards::Quadrilateral<9>::key:
713  case shards::ShellQuadrilateral<9>::key:
714  return quadrilateral_9[nodeOrd];
715  break;
716 
717 
718  // Base Tetrahedron topology
719  case shards::Tetrahedron<4>::key:
720  return tetrahedron[nodeOrd];
721  break;
722 
723  // Extended Tetrahedron topologies
724  case shards::Tetrahedron<8>::key:
725  return tetrahedron_8[nodeOrd];
726  break;
727  case shards::Tetrahedron<10>::key:
728  return tetrahedron_10[nodeOrd];
729  break;
730  case shards::Tetrahedron<11>::key:
731  return tetrahedron_11[nodeOrd];
732  break;
733 
734 
735  // Base Hexahedron topology
736  case shards::Hexahedron<8>::key:
737  return hexahedron[nodeOrd];
738  break;
739 
740  // Extended Hexahedron topologies
741  case shards::Hexahedron<20>::key:
742  return hexahedron_20[nodeOrd];
743  break;
744  case shards::Hexahedron<27>::key:
745  return hexahedron_27[nodeOrd];
746  break;
747 
748 
749  // Base Pyramid topology
750  case shards::Pyramid<5>::key:
751  return pyramid[nodeOrd];
752  break;
753 
754  // Extended pyramid topologies
755  case shards::Pyramid<13>::key:
756  return pyramid_13[nodeOrd];
757  break;
758  case shards::Pyramid<14>::key:
759  return pyramid_14[nodeOrd];
760  break;
761 
762 
763  // Base Wedge topology
764  case shards::Wedge<6>::key:
765  return wedge[nodeOrd];
766  break;
767 
768  // Extended Wedge topologies
769  case shards::Wedge<15>::key:
770  return wedge_15[nodeOrd];
771  break;
772  case shards::Wedge<18>::key:
773  return wedge_18[nodeOrd];
774  break;
775 
776  default:
777  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
778  ">>> ERROR (Intrepid::CellTools::getReferenceNode): invalid cell topology.");
779  }
780  // To disable compiler warning, should never be reached
781  return line[0];
782  }
783 
784 
785 
786  template<class Scalar>
787  template<class ArraySubcellNode>
788  void CellTools<Scalar>::getReferenceSubcellNodes(ArraySubcellNode & subcellNodes,
789  const int subcellDim,
790  const int subcellOrd,
791  const shards::CellTopology& parentCell){
792 #ifdef HAVE_INTREPID_DEBUG
793  TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument,
794  ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes): the specified cell topology does not have a reference cell.");
795 
796  // subcellDim can equal the cell dimension because the cell itself is a valid subcell! In this case
797  // the method will return all cell cellWorkset.
798  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellDim) && (subcellDim <= (int)parentCell.getDimension()) ), std::invalid_argument,
799  ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes): subcell dimension out of range.");
800 
801  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellOrd) && (subcellOrd < (int)parentCell.getSubcellCount(subcellDim) ) ), std::invalid_argument,
802  ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes): subcell ordinal out of range.");
803 
804  // Verify subcellNodes rank and dimensions
805  {
806  std::string errmsg = ">>> ERROR (Intrepid::CellTools::getReferenceSubcellNodes):";
807  TEUCHOS_TEST_FOR_EXCEPTION( !( requireRankRange(errmsg, subcellNodes, 2, 2) ), std::invalid_argument, errmsg);
808 
809  int subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
810  int spaceDim = parentCell.getDimension();
811 
812  TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellNodes, 0, subcNodeCount, subcNodeCount) ),
813  std::invalid_argument, errmsg);
814 
815  TEUCHOS_TEST_FOR_EXCEPTION( !( requireDimensionRange(errmsg, subcellNodes, 1, spaceDim, spaceDim) ),
816  std::invalid_argument, errmsg);
817  }
818 #endif
819 
820  // Find how many cellWorkset does the specified subcell have.
821  int subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
822 
823  // Loop over subcell cellWorkset
824  for(int subcNodeOrd = 0; subcNodeOrd < subcNodeCount; subcNodeOrd++){
825 
826  // Get the node number relative to the parent reference cell
827  int cellNodeOrd = parentCell.getNodeMap(subcellDim, subcellOrd, subcNodeOrd);
828 
829  // Loop over node's Cartesian coordinates
830  for(int dim = 0; dim < (int)parentCell.getDimension(); dim++){
831  subcellNodes(subcNodeOrd, dim) = CellTools::getReferenceNode(parentCell, cellNodeOrd)[dim];
832  }
833  }
834  }
835 
836 
837 
838  template<class Scalar>
839  int CellTools<Scalar>::hasReferenceCell(const shards::CellTopology& cell) {
840 
841  switch(cell.getKey() ) {
842  case shards::Line<2>::key:
843  case shards::Line<3>::key:
844  case shards::ShellLine<2>::key:
845  case shards::ShellLine<3>::key:
846  case shards::Beam<2>::key:
847  case shards::Beam<3>::key:
848 
849  case shards::Triangle<3>::key:
850  case shards::Triangle<4>::key:
851  case shards::Triangle<6>::key:
852  case shards::ShellTriangle<3>::key:
853  case shards::ShellTriangle<6>::key:
854 
855  case shards::Quadrilateral<4>::key:
856  case shards::Quadrilateral<8>::key:
857  case shards::Quadrilateral<9>::key:
858  case shards::ShellQuadrilateral<4>::key:
859  case shards::ShellQuadrilateral<8>::key:
860  case shards::ShellQuadrilateral<9>::key:
861 
862  case shards::Tetrahedron<4>::key:
863  case shards::Tetrahedron<8>::key:
864  case shards::Tetrahedron<10>::key:
865  case shards::Tetrahedron<11>::key:
866 
867  case shards::Hexahedron<8>::key:
868  case shards::Hexahedron<20>::key:
869  case shards::Hexahedron<27>::key:
870 
871  case shards::Pyramid<5>::key:
872  case shards::Pyramid<13>::key:
873  case shards::Pyramid<14>::key:
874 
875  case shards::Wedge<6>::key:
876  case shards::Wedge<15>::key:
877  case shards::Wedge<18>::key:
878  return 1;
879  break;
880 
881  default:
882  return 0;
883  }
884  return 0;
885  }
886 
887  //============================================================================================//
888  // //
889  // Jacobian, inverse Jacobian and Jacobian determinant //
890  // //
891  //============================================================================================//
892 
893 
894 
895 
896  template<class Scalar>
897  template<class ArrayJac, class ArrayPoint, class ArrayCell>
898  void CellTools<Scalar>::setJacobian(ArrayJac & jacobian,
899  const ArrayPoint & points,
900  const ArrayCell & cellWorkset,
901  const shards::CellTopology & cellTopo,
902  const int & whichCell)
903  {
904  INTREPID_VALIDATE( validateArguments_setJacobian(jacobian, points, cellWorkset, whichCell, cellTopo) );
905 
906  ArrayWrapper<Scalar,ArrayJac, Rank<ArrayJac >::value, false>jacobianWrap(jacobian);
907  ArrayWrapper<Scalar,ArrayPoint, Rank<ArrayPoint >::value, true>pointsWrap(points);
908  ArrayWrapper<Scalar,ArrayCell, Rank<ArrayCell >::value, true>cellWorksetWrap(cellWorkset);
909  int spaceDim = (size_t)cellTopo.getDimension();
910  size_t numCells = static_cast<size_t>(cellWorkset.dimension(0));
911  //points can be rank-2 (P,D), or rank-3 (C,P,D)
912  size_t numPoints = (getrank(points) == 2) ? static_cast<size_t>(points.dimension(0)) : static_cast<size_t>(points.dimension(1));
913 
914  // Jacobian is computed using gradients of an appropriate H(grad) basis function: define RCP to the base class
915  Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
916 
917  // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
918  switch( cellTopo.getKey() ){
919 
920  // Standard Base topologies (number of cellWorkset = number of vertices)
921  case shards::Line<2>::key:
922  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
923  break;
924 
925  case shards::Triangle<3>::key:
926  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
927  break;
928 
929  case shards::Quadrilateral<4>::key:
930  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
931  break;
932 
933  case shards::Tetrahedron<4>::key:
934  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
935  break;
936 
937  case shards::Hexahedron<8>::key:
938  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
939  break;
940 
941  case shards::Wedge<6>::key:
942  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
943  break;
944 
945  case shards::Pyramid<5>::key:
946  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
947  break;
948 
949  // Standard Extended topologies
950  case shards::Triangle<6>::key:
951  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
952  break;
953  case shards::Quadrilateral<9>::key:
954  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
955  break;
956 
957  case shards::Tetrahedron<10>::key:
958  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
959  break;
960 
961  case shards::Tetrahedron<11>::key:
962  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
963  break;
964 
965  case shards::Hexahedron<20>::key:
966  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<Scalar, FieldContainer<Scalar> >() );
967  break;
968 
969  case shards::Hexahedron<27>::key:
970  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
971  break;
972 
973  case shards::Wedge<15>::key:
974  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<Scalar, FieldContainer<Scalar> >() );
975  break;
976 
977  case shards::Wedge<18>::key:
978  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
979  break;
980 
981  case shards::Pyramid<13>::key:
982  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<Scalar, FieldContainer<Scalar> >() );
983  break;
984 
985  // These extended topologies are not used for mapping purposes
986  case shards::Quadrilateral<8>::key:
987  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
988  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
989  break;
990 
991  // Base and Extended Line, Beam and Shell topologies
992  case shards::Line<3>::key:
993  case shards::Beam<2>::key:
994  case shards::Beam<3>::key:
995  case shards::ShellLine<2>::key:
996  case shards::ShellLine<3>::key:
997  case shards::ShellTriangle<3>::key:
998  case shards::ShellTriangle<6>::key:
999  case shards::ShellQuadrilateral<4>::key:
1000  case shards::ShellQuadrilateral<8>::key:
1001  case shards::ShellQuadrilateral<9>::key:
1002  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1003  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
1004  break;
1005  default:
1006  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1007  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported.");
1008  }// switch
1009 
1010  // Temp (F,P,D) array for the values of basis functions gradients at the reference points
1011  int basisCardinality = HGRAD_Basis -> getCardinality();
1012  FieldContainer<Scalar> basisGrads(basisCardinality, numPoints, spaceDim);
1013 
1014 
1015 if(getrank(jacobian)==4){
1016  for (size_t i=0; i< static_cast<size_t>(jacobian.dimension(0)); i++){
1017  for (size_t j=0; j< static_cast<size_t>(jacobian.dimension(1)); j++){
1018  for (size_t k=0; k< static_cast<size_t>(jacobian.dimension(2)); k++){
1019  for (size_t l=0; l< static_cast<size_t>(jacobian.dimension(3)); l++){
1020  jacobianWrap(i,j,k,l)=0.0;
1021  }
1022  }
1023  }
1024  }
1025 }
1026 
1027 if(getrank(jacobian)==3){
1028  for (size_t i=0; i< static_cast<size_t>(jacobian.dimension(0)); i++){
1029  for (size_t j=0; j< static_cast<size_t>(jacobian.dimension(1)); j++){
1030  for (size_t k=0; k< static_cast<size_t>(jacobian.dimension(2)); k++){
1031  jacobianWrap(i,j,k)=0.0;
1032  }
1033  }
1034  }
1035 }
1036  // Handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of points arrays.
1037  switch(getrank(points)) {
1038 
1039  // refPoints is (P,D): a single or multiple cell jacobians computed for a single set of ref. points
1040  case 2:
1041  {
1042  // getValues requires rank-2 (P,D) input array, but points cannot be passed directly as argument because they are a user type
1043  FieldContainer<Scalar> tempPoints( static_cast<size_t>(points.dimension(0)), static_cast<size_t>(points.dimension(1)) );
1044  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1045  for(size_t pt = 0; pt < static_cast<size_t>(points.dimension(0)); pt++){
1046  for(size_t dm = 0; dm < static_cast<size_t>(points.dimension(1)) ; dm++){
1047  tempPoints(pt, dm) = pointsWrap(pt, dm);
1048  }//dm
1049  }//pt
1050 
1051  HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
1052 
1053  // The outer loops select the multi-index of the Jacobian entry: cell, point, row, col
1054  // If whichCell = -1, all jacobians are computed, otherwise a single cell jacobian is computed
1055  size_t cellLoop = (whichCell == -1) ? numCells : 1 ;
1056 
1057  if(whichCell == -1) {
1058  for(size_t cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
1059  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1060  for(int row = 0; row < spaceDim; row++){
1061  for(int col = 0; col < spaceDim; col++){
1062 
1063  // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
1064  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1065  jacobianWrap(cellOrd, pointOrd, row, col) += cellWorksetWrap(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
1066  } // bfOrd
1067  } // col
1068  } // row
1069  } // pointOrd
1070  } // cellOrd
1071 
1072  //}
1073 
1074  }
1075  else {
1076  for(size_t cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
1077  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1078  for(int row = 0; row < spaceDim; row++){
1079  for(int col = 0; col < spaceDim; col++){
1080 
1081  // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
1082  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1083  jacobianWrap(pointOrd, row, col) += cellWorksetWrap(whichCell, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
1084  } // bfOrd
1085  } // col
1086  } // row
1087  } // pointOrd
1088  } // cellOrd
1089 // }
1090  } // if whichcell
1091  }// case 2
1092  break;
1093 
1094  // points is (C,P,D): multiple jacobians computed at multiple point sets, one jacobian per cell
1095  case 3:
1096  {
1097  // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
1098  FieldContainer<Scalar> tempPoints( static_cast<size_t>(points.dimension(1)), static_cast<size_t>(points.dimension(2)) );
1099  for(size_t cellOrd = 0; cellOrd < numCells; cellOrd++) {
1100 
1101  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1102  for(size_t pt = 0; pt < static_cast<size_t>(points.dimension(1)); pt++){
1103  for(size_t dm = 0; dm < static_cast<size_t>(points.dimension(2)) ; dm++){
1104  tempPoints(pt, dm) = pointsWrap(cellOrd, pt, dm);
1105  }//dm
1106  }//pt
1107 
1108  // Compute gradients of basis functions at this set of ref. points
1109  HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
1110 
1111  // Compute jacobians for the point set corresponding to the current cellordinal
1112  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1113  for(int row = 0; row < spaceDim; row++){
1114  for(int col = 0; col < spaceDim; col++){
1115 
1116  // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same
1117  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1118  jacobianWrap(cellOrd, pointOrd, row, col) += cellWorksetWrap(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
1119  } // bfOrd
1120  } // col
1121  } // row
1122  } // pointOrd
1123  }//cellOrd
1124 // }
1125  }// case 3
1126 
1127  break;
1128 
1129  default:
1130  TEUCHOS_TEST_FOR_EXCEPTION( !( (getrank(points) == 2) && (getrank(points) == 3) ), std::invalid_argument,
1131  ">>> ERROR (Intrepid::CellTools::setJacobian): rank 2 or 3 required for points array. ");
1132  }//switch
1133  }
1134 
1135 
1136  template<class Scalar>
1137  template<class ArrayJac, class ArrayPoint, class ArrayCell>
1138  void CellTools<Scalar>::setJacobian(ArrayJac & jacobian,
1139  const ArrayPoint & points,
1140  const ArrayCell & cellWorkset,
1141  const Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis,
1142  const int & whichCell)
1143  {
1144  //IKT, 10/7/15: OK to not validate arguments for this implementation?
1145  //INTREPID_VALIDATE( validateArguments_setJacobian(jacobian, points, cellWorkset, whichCell, cellTopo) );
1146 
1147  ArrayWrapper<Scalar,ArrayJac, Rank<ArrayJac >::value, false>jacobianWrap(jacobian);
1148  ArrayWrapper<Scalar,ArrayPoint, Rank<ArrayPoint >::value, true>pointsWrap(points);
1149  ArrayWrapper<Scalar,ArrayCell, Rank<ArrayCell >::value, true>cellWorksetWrap(cellWorkset);
1150  int spaceDim = (size_t)HGRAD_Basis->getBaseCellTopology().getDimension();
1151  size_t numCells = static_cast<size_t>(cellWorkset.dimension(0));
1152  //points can be rank-2 (P,D), or rank-3 (C,P,D)
1153  size_t numPoints = (getrank(points) == 2) ? static_cast<size_t>(points.dimension(0)) : static_cast<size_t>(points.dimension(1));
1154 
1155  // Temp (F,P,D) array for the values of basis functions gradients at the reference points
1156  int basisCardinality = HGRAD_Basis -> getCardinality();
1157  FieldContainer<Scalar> basisGrads(basisCardinality, numPoints, spaceDim);
1158 
1159 
1160 if(getrank(jacobian)==4){
1161  for (size_t i=0; i< static_cast<size_t>(jacobian.dimension(0)); i++){
1162  for (size_t j=0; j< static_cast<size_t>(jacobian.dimension(1)); j++){
1163  for (size_t k=0; k< static_cast<size_t>(jacobian.dimension(2)); k++){
1164  for (size_t l=0; l< static_cast<size_t>(jacobian.dimension(3)); l++){
1165  jacobianWrap(i,j,k,l)=0.0;
1166  }
1167  }
1168  }
1169  }
1170 }
1171 
1172 if(getrank(jacobian)==3){
1173  for (size_t i=0; i< static_cast<size_t>(jacobian.dimension(0)); i++){
1174  for (size_t j=0; j< static_cast<size_t>(jacobian.dimension(1)); j++){
1175  for (size_t k=0; k< static_cast<size_t>(jacobian.dimension(2)); k++){
1176  jacobianWrap(i,j,k)=0.0;
1177  }
1178  }
1179  }
1180 }
1181  // Handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of points arrays.
1182  switch(getrank(points)) {
1183 
1184  // refPoints is (P,D): a single or multiple cell jacobians computed for a single set of ref. points
1185  case 2:
1186  {
1187  // getValues requires rank-2 (P,D) input array, but points cannot be passed directly as argument because they are a user type
1188  FieldContainer<Scalar> tempPoints( static_cast<size_t>(points.dimension(0)), static_cast<size_t>(points.dimension(1)) );
1189  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1190  for(size_t pt = 0; pt < static_cast<size_t>(points.dimension(0)); pt++){
1191  for(size_t dm = 0; dm < static_cast<size_t>(points.dimension(1)) ; dm++){
1192  tempPoints(pt, dm) = pointsWrap(pt, dm);
1193  }//dm
1194  }//pt
1195 
1196  HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
1197 
1198  // The outer loops select the multi-index of the Jacobian entry: cell, point, row, col
1199  // If whichCell = -1, all jacobians are computed, otherwise a single cell jacobian is computed
1200  size_t cellLoop = (whichCell == -1) ? numCells : 1 ;
1201 
1202  if(whichCell == -1) {
1203  for(size_t cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
1204  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1205  for(int row = 0; row < spaceDim; row++){
1206  for(int col = 0; col < spaceDim; col++){
1207 
1208  // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
1209  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1210  jacobianWrap(cellOrd, pointOrd, row, col) += cellWorksetWrap(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
1211  } // bfOrd
1212  } // col
1213  } // row
1214  } // pointOrd
1215  } // cellOrd
1216 
1217  //}
1218 
1219  }
1220  else {
1221  for(size_t cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
1222  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1223  for(int row = 0; row < spaceDim; row++){
1224  for(int col = 0; col < spaceDim; col++){
1225 
1226  // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same.
1227  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1228  jacobianWrap(pointOrd, row, col) += cellWorksetWrap(whichCell, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
1229  } // bfOrd
1230  } // col
1231  } // row
1232  } // pointOrd
1233  } // cellOrd
1234 // }
1235  } // if whichcell
1236  }// case 2
1237  break;
1238 
1239  // points is (C,P,D): multiple jacobians computed at multiple point sets, one jacobian per cell
1240  case 3:
1241  {
1242  // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
1243  FieldContainer<Scalar> tempPoints( static_cast<size_t>(points.dimension(1)), static_cast<size_t>(points.dimension(2)) );
1244  for(size_t cellOrd = 0; cellOrd < numCells; cellOrd++) {
1245 
1246  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1247  for(size_t pt = 0; pt < static_cast<size_t>(points.dimension(1)); pt++){
1248  for(size_t dm = 0; dm < static_cast<size_t>(points.dimension(2)) ; dm++){
1249  tempPoints(pt, dm) = pointsWrap(cellOrd, pt, dm);
1250  }//dm
1251  }//pt
1252 
1253  // Compute gradients of basis functions at this set of ref. points
1254  HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
1255 
1256  // Compute jacobians for the point set corresponding to the current cellordinal
1257  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1258  for(int row = 0; row < spaceDim; row++){
1259  for(int col = 0; col < spaceDim; col++){
1260 
1261  // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same
1262  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1263  jacobianWrap(cellOrd, pointOrd, row, col) += cellWorksetWrap(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
1264  } // bfOrd
1265  } // col
1266  } // row
1267  } // pointOrd
1268  }//cellOrd
1269 // }
1270  }// case 3
1271 
1272  break;
1273 
1274  default:
1275  TEUCHOS_TEST_FOR_EXCEPTION( !( (getrank(points) == 2) && (getrank(points) == 3) ), std::invalid_argument,
1276  ">>> ERROR (Intrepid::CellTools::setJacobian): rank 2 or 3 required for points array. ");
1277  }//switch
1278 
1279  }
1280 
1281 template<class Scalar>
1282 template<class ArrayJacInv, class ArrayJac>
1283 void CellTools<Scalar>::setJacobianInv(ArrayJacInv & jacobianInv,
1284  const ArrayJac & jacobian)
1285 {
1286  INTREPID_VALIDATE( validateArguments_setJacobianInv(jacobianInv, jacobian) );
1287 
1288  RealSpaceTools<Scalar>::inverse(jacobianInv, jacobian);
1289 }
1290 /*
1291 template<class Scalar>
1292 template<class ArrayJacInv, class ArrayJac>
1293 void CellTools<Scalar>::setJacobianInvTemp(ArrayJacInv & jacobianInv,
1294  const ArrayJac & jacobian)
1295 {
1296 
1297 
1298  RealSpaceTools<Scalar>::inverseTemp(jacobianInv, jacobian);
1299 }
1300 */
1301 /*
1302 template<class Scalar>
1303 template<class ArrayJacDet, class ArrayJac>
1304 void CellTools<Scalar>::setJacobianDet(ArrayJacDet & jacobianDet,
1305  const ArrayJac & jacobian)
1306 {
1307  INTREPID_VALIDATE( validateArguments_setJacobianDetArgs(jacobianDet, jacobian) );
1308 
1309  RealSpaceTools<Scalar>::det(jacobianDet, jacobian);
1310 }
1311 */
1312 template<class Scalar>
1313 template<class ArrayJacDet, class ArrayJac>
1314 void CellTools<Scalar>::setJacobianDet(ArrayJacDet & jacobianDet,
1315  const ArrayJac & jacobian)
1316 {
1317  INTREPID_VALIDATE( validateArguments_setJacobianDetArgs(jacobianDet, jacobian) );
1318  RealSpaceTools<Scalar>::det(jacobianDet, jacobian);
1319 }
1320 
1321 //============================================================================================//
1322 // //
1323 // Reference-to-physical frame mapping and its inverse //
1324 // //
1325 //============================================================================================//
1326 
1327 template<class Scalar>
1328 template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
1329 void CellTools<Scalar>::mapToPhysicalFrame(ArrayPhysPoint & physPoints,
1330  const ArrayRefPoint & refPoints,
1331  const ArrayCell & cellWorkset,
1332  const Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis,
1333  const int & whichCell)
1334 {
1335  //INTREPID_VALIDATE(validateArguments_mapToPhysicalFrame( physPoints, refPoints, cellWorkset, cellTopo, whichCell) );
1336 
1337  ArrayWrapper<Scalar,ArrayPhysPoint, Rank<ArrayPhysPoint >::value, false>physPointsWrap(physPoints);
1338  ArrayWrapper<Scalar,ArrayRefPoint, Rank<ArrayRefPoint >::value, true>refPointsWrap(refPoints);
1339  ArrayWrapper<Scalar,ArrayCell, Rank<ArrayCell >::value,true>cellWorksetWrap(cellWorkset);
1340 
1341  size_t spaceDim = (size_t)HGRAD_Basis->getBaseCellTopology().getDimension();
1342 
1343  size_t numCells = static_cast<size_t>(cellWorkset.dimension(0));
1344  //points can be rank-2 (P,D), or rank-3 (C,P,D)
1345  size_t numPoints = (getrank(refPoints) == 2) ? static_cast<size_t>(refPoints.dimension(0)) : static_cast<size_t>(refPoints.dimension(1));
1346 
1347  // Temp (F,P) array for the values of nodal basis functions at the reference points
1348  int basisCardinality = HGRAD_Basis -> getCardinality();
1349  FieldContainer<Scalar> basisVals(basisCardinality, numPoints);
1350 
1351  // Initialize physPoints
1352  if(getrank(physPoints)==3){
1353 for(size_t i = 0; i < static_cast<size_t>(physPoints.dimension(0)); i++) {
1354  for(size_t j = 0; j < static_cast<size_t>(physPoints.dimension(1)); j++){
1355  for(size_t k = 0; k < static_cast<size_t>(physPoints.dimension(2)); k++){
1356  physPointsWrap(i,j,k) = 0.0;
1357  }
1358  }
1359 }
1360  }else if(getrank(physPoints)==2){
1361  for(size_t i = 0; i < static_cast<size_t>(physPoints.dimension(0)); i++){
1362  for(size_t j = 0; j < static_cast<size_t>(physPoints.dimension(1)); j++){
1363  physPointsWrap(i,j) = 0.0;
1364  }
1365  }
1366 
1367  }
1368 
1369  // handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of refPoints
1370  switch(getrank(refPoints)) {
1371 
1372  // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
1373  case 2:
1374  {
1375 
1376  // getValues requires rank-2 (P,D) input array, but refPoints cannot be passed directly as argument because they are a user type
1377  FieldContainer<Scalar> tempPoints( static_cast<size_t>(refPoints.dimension(0)), static_cast<size_t>(refPoints.dimension(1)) );
1378  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1379  for(size_t pt = 0; pt < static_cast<size_t>(refPoints.dimension(0)); pt++){
1380  for(size_t dm = 0; dm < static_cast<size_t>(refPoints.dimension(1)) ; dm++){
1381  tempPoints(pt, dm) = refPointsWrap(pt, dm);
1382  }//dm
1383  }//pt
1384  HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
1385 
1386  // If whichCell = -1, ref pt. set is mapped to all cells, otherwise, the set is mapped to one cell only
1387  size_t cellLoop = (whichCell == -1) ? numCells : 1 ;
1388 
1389  // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
1390  for(size_t cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
1391  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1392  for(size_t dim = 0; dim < spaceDim; dim++){
1393  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1394 
1395  if(whichCell == -1){
1396  physPointsWrap(cellOrd, pointOrd, dim) += cellWorksetWrap(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
1397  }
1398  else{
1399  physPointsWrap(pointOrd, dim) += cellWorksetWrap(whichCell, bfOrd, dim)*basisVals(bfOrd, pointOrd);
1400  }
1401  } // bfOrd
1402  }// dim
1403  }// pointOrd
1404  }//cellOrd
1405  }// case 2
1406 
1407  break;
1408 
1409  // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.
1410  case 3:
1411  {
1412 
1413  // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
1414  FieldContainer<Scalar> tempPoints( static_cast<size_t>(refPoints.dimension(1)), static_cast<size_t>(refPoints.dimension(2)) );
1415 
1416  // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
1417  for(size_t cellOrd = 0; cellOrd < numCells; cellOrd++) {
1418 
1419  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1420  for(size_t pt = 0; pt < static_cast<size_t>(refPoints.dimension(1)); pt++){
1421  for(size_t dm = 0; dm < static_cast<size_t>(refPoints.dimension(2)) ; dm++){
1422  tempPoints(pt, dm) = refPointsWrap(cellOrd, pt, dm);
1423  }//dm
1424  }//pt
1425 
1426  // Compute basis values for this set of ref. points
1427  HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
1428 
1429  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1430  for(size_t dim = 0; dim < spaceDim; dim++){
1431  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1432 
1433  physPointsWrap(cellOrd, pointOrd, dim) += cellWorksetWrap(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
1434 
1435  } // bfOrd
1436  }// dim
1437  }// pointOrd
1438  }//cellOrd
1439  }// case 3
1440  break;
1441 
1442  }
1443 }
1444 
1445 template<class Scalar>
1446 template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
1447 void CellTools<Scalar>::mapToPhysicalFrame(ArrayPhysPoint & physPoints,
1448  const ArrayRefPoint & refPoints,
1449  const ArrayCell & cellWorkset,
1450  const shards::CellTopology & cellTopo,
1451  const int & whichCell)
1452 {
1453  INTREPID_VALIDATE(validateArguments_mapToPhysicalFrame( physPoints, refPoints, cellWorkset, cellTopo, whichCell) );
1454 
1455  ArrayWrapper<Scalar,ArrayPhysPoint, Rank<ArrayPhysPoint >::value, false>physPointsWrap(physPoints);
1456  ArrayWrapper<Scalar,ArrayRefPoint, Rank<ArrayRefPoint >::value, true>refPointsWrap(refPoints);
1457  ArrayWrapper<Scalar,ArrayCell, Rank<ArrayCell >::value,true>cellWorksetWrap(cellWorkset);
1458 
1459  size_t spaceDim = (size_t)cellTopo.getDimension();
1460  size_t numCells = static_cast<size_t>(cellWorkset.dimension(0));
1461  //points can be rank-2 (P,D), or rank-3 (C,P,D)
1462  size_t numPoints = (getrank(refPoints) == 2) ? static_cast<size_t>(refPoints.dimension(0)) : static_cast<size_t>(refPoints.dimension(1));
1463 
1464  // Mapping is computed using an appropriate H(grad) basis function: define RCP to the base class
1465  Teuchos::RCP<Basis<Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
1466 
1467  // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
1468  switch( cellTopo.getKey() ){
1469 
1470  // Standard Base topologies (number of cellWorkset = number of vertices)
1471  case shards::Line<2>::key:
1472  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
1473  break;
1474 
1475  case shards::Triangle<3>::key:
1476  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
1477  break;
1478 
1479  case shards::Quadrilateral<4>::key:
1480  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
1481  break;
1482 
1483  case shards::Tetrahedron<4>::key:
1484  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
1485  break;
1486 
1487  case shards::Hexahedron<8>::key:
1488  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
1489  break;
1490 
1491  case shards::Wedge<6>::key:
1492  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
1493  break;
1494 
1495  case shards::Pyramid<5>::key:
1496  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
1497  break;
1498 
1499  // Standard Extended topologies
1500  case shards::Triangle<6>::key:
1501  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
1502  break;
1503 
1504  case shards::Quadrilateral<9>::key:
1505  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
1506  break;
1507 
1508  case shards::Tetrahedron<10>::key:
1509  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
1510  break;
1511 
1512  case shards::Tetrahedron<11>::key:
1513  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
1514  break;
1515 
1516  case shards::Hexahedron<20>::key:
1517  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<Scalar, FieldContainer<Scalar> >() );
1518  break;
1519 
1520  case shards::Hexahedron<27>::key:
1521  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
1522  break;
1523 
1524  case shards::Wedge<15>::key:
1525  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<Scalar, FieldContainer<Scalar> >() );
1526  break;
1527 
1528  case shards::Wedge<18>::key:
1529  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
1530  break;
1531 
1532  case shards::Pyramid<13>::key:
1533  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<Scalar, FieldContainer<Scalar> >() );
1534  break;
1535 
1536  // These extended topologies are not used for mapping purposes
1537  case shards::Quadrilateral<8>::key:
1538  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1539  ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
1540  break;
1541 
1542  // Base and Extended Line, Beam and Shell topologies
1543  case shards::Line<3>::key:
1544  case shards::Beam<2>::key:
1545  case shards::Beam<3>::key:
1546  case shards::ShellLine<2>::key:
1547  case shards::ShellLine<3>::key:
1548  case shards::ShellTriangle<3>::key:
1549  case shards::ShellTriangle<6>::key:
1550  case shards::ShellQuadrilateral<4>::key:
1551  case shards::ShellQuadrilateral<8>::key:
1552  case shards::ShellQuadrilateral<9>::key:
1553  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1554  ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
1555  break;
1556  default:
1557  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1558  ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported.");
1559  }// switch
1560 
1561  // Temp (F,P) array for the values of nodal basis functions at the reference points
1562  int basisCardinality = HGRAD_Basis -> getCardinality();
1563  FieldContainer<Scalar> basisVals(basisCardinality, numPoints);
1564 
1565  // Initialize physPoints
1566  if(getrank(physPoints)==3){
1567 for(size_t i = 0; i < static_cast<size_t>(physPoints.dimension(0)); i++) {
1568  for(size_t j = 0; j < static_cast<size_t>(physPoints.dimension(1)); j++){
1569  for(size_t k = 0; k < static_cast<size_t>(physPoints.dimension(2)); k++){
1570  physPointsWrap(i,j,k) = 0.0;
1571  }
1572  }
1573 }
1574  }else if(getrank(physPoints)==2){
1575  for(size_t i = 0; i < static_cast<size_t>(physPoints.dimension(0)); i++){
1576  for(size_t j = 0; j < static_cast<size_t>(physPoints.dimension(1)); j++){
1577  physPointsWrap(i,j) = 0.0;
1578  }
1579  }
1580 
1581  }
1582 
1583  // handle separately rank-2 (P,D) and rank-3 (C,P,D) cases of refPoints
1584  switch(getrank(refPoints)) {
1585 
1586  // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
1587  case 2:
1588  {
1589 
1590  // getValues requires rank-2 (P,D) input array, but refPoints cannot be passed directly as argument because they are a user type
1591  FieldContainer<Scalar> tempPoints( static_cast<size_t>(refPoints.dimension(0)), static_cast<size_t>(refPoints.dimension(1)) );
1592  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1593  for(size_t pt = 0; pt < static_cast<size_t>(refPoints.dimension(0)); pt++){
1594  for(size_t dm = 0; dm < static_cast<size_t>(refPoints.dimension(1)) ; dm++){
1595  tempPoints(pt, dm) = refPointsWrap(pt, dm);
1596  }//dm
1597  }//pt
1598  HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
1599 
1600  // If whichCell = -1, ref pt. set is mapped to all cells, otherwise, the set is mapped to one cell only
1601  size_t cellLoop = (whichCell == -1) ? numCells : 1 ;
1602 
1603  // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
1604  for(size_t cellOrd = 0; cellOrd < cellLoop; cellOrd++) {
1605  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1606  for(size_t dim = 0; dim < spaceDim; dim++){
1607  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1608 
1609  if(whichCell == -1){
1610  physPointsWrap(cellOrd, pointOrd, dim) += cellWorksetWrap(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
1611  }
1612  else{
1613  physPointsWrap(pointOrd, dim) += cellWorksetWrap(whichCell, bfOrd, dim)*basisVals(bfOrd, pointOrd);
1614  }
1615  } // bfOrd
1616  }// dim
1617  }// pointOrd
1618  }//cellOrd
1619  }// case 2
1620 
1621  break;
1622 
1623  // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.
1624  case 3:
1625  {
1626 
1627  // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
1628  FieldContainer<Scalar> tempPoints( static_cast<size_t>(refPoints.dimension(1)), static_cast<size_t>(refPoints.dimension(2)) );
1629 
1630  // Compute the map F(refPoints) = sum node_coordinate*basis(refPoints)
1631  for(size_t cellOrd = 0; cellOrd < numCells; cellOrd++) {
1632 
1633  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
1634  for(size_t pt = 0; pt < static_cast<size_t>(refPoints.dimension(1)); pt++){
1635  for(size_t dm = 0; dm < static_cast<size_t>(refPoints.dimension(2)) ; dm++){
1636  tempPoints(pt, dm) = refPointsWrap(cellOrd, pt, dm);
1637  }//dm
1638  }//pt
1639 
1640  // Compute basis values for this set of ref. points
1641  HGRAD_Basis -> getValues(basisVals, tempPoints, OPERATOR_VALUE);
1642 
1643  for(size_t pointOrd = 0; pointOrd < numPoints; pointOrd++) {
1644  for(size_t dim = 0; dim < spaceDim; dim++){
1645  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
1646 
1647  physPointsWrap(cellOrd, pointOrd, dim) += cellWorksetWrap(cellOrd, bfOrd, dim)*basisVals(bfOrd, pointOrd);
1648 
1649  } // bfOrd
1650  }// dim
1651  }// pointOrd
1652  }//cellOrd
1653  }// case 3
1654  break;
1655 
1656 
1657  }
1658 }
1659 
1660 template<class Scalar>
1661 template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
1662 void CellTools<Scalar>::mapToReferenceFrame(ArrayRefPoint & refPoints,
1663  const ArrayPhysPoint & physPoints,
1664  const ArrayCell & cellWorkset,
1665  const shards::CellTopology & cellTopo,
1666  const int & whichCell)
1667 {
1668  INTREPID_VALIDATE( validateArguments_mapToReferenceFrame(refPoints, physPoints, cellWorkset, cellTopo, whichCell) );
1669 
1670  size_t spaceDim = (size_t)cellTopo.getDimension();
1671  size_t numPoints;
1672  size_t numCells;
1673 
1674  // Define initial guesses to be the Cell centers of the reference cell topology
1675  FieldContainer<Scalar> cellCenter(spaceDim);
1676  switch( cellTopo.getKey() ){
1677  // Standard Base topologies (number of cellWorkset = number of vertices)
1678  case shards::Line<2>::key:
1679  cellCenter(0) = 0.0; break;
1680 
1681  case shards::Triangle<3>::key:
1682  case shards::Triangle<6>::key:
1683  cellCenter(0) = 1./3.; cellCenter(1) = 1./3.; break;
1684 
1685  case shards::Quadrilateral<4>::key:
1686  case shards::Quadrilateral<9>::key:
1687  cellCenter(0) = 0.0; cellCenter(1) = 0.0; break;
1688 
1689  case shards::Tetrahedron<4>::key:
1690  case shards::Tetrahedron<10>::key:
1691  case shards::Tetrahedron<11>::key:
1692  cellCenter(0) = 1./6.; cellCenter(1) = 1./6.; cellCenter(2) = 1./6.; break;
1693 
1694  case shards::Hexahedron<8>::key:
1695  case shards::Hexahedron<20>::key:
1696  case shards::Hexahedron<27>::key:
1697  cellCenter(0) = 0.0; cellCenter(1) = 0.0; cellCenter(2) = 0.0; break;
1698 
1699  case shards::Wedge<6>::key:
1700  case shards::Wedge<15>::key:
1701  case shards::Wedge<18>::key:
1702  cellCenter(0) = 1./3.; cellCenter(1) = 1./3.; cellCenter(2) = 0.0; break;
1703 
1704  case shards::Pyramid<5>::key:
1705  case shards::Pyramid<13>::key:
1706  cellCenter(0) = 0.; cellCenter(1) = 0.; cellCenter(2) = 0.25; break;
1707 
1708  // These extended topologies are not used for mapping purposes
1709  case shards::Quadrilateral<8>::key:
1710  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1711  ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported. ");
1712  break;
1713 
1714  // Base and Extended Line, Beam and Shell topologies
1715  case shards::Line<3>::key:
1716  case shards::Beam<2>::key:
1717  case shards::Beam<3>::key:
1718  case shards::ShellLine<2>::key:
1719  case shards::ShellLine<3>::key:
1720  case shards::ShellTriangle<3>::key:
1721  case shards::ShellTriangle<6>::key:
1722  case shards::ShellQuadrilateral<4>::key:
1723  case shards::ShellQuadrilateral<8>::key:
1724  case shards::ShellQuadrilateral<9>::key:
1725  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1726  ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported. ");
1727  break;
1728  default:
1729  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
1730  ">>> ERROR (Intrepid::CellTools::mapToReferenceFrame): Cell topology not supported.");
1731  }// switch key
1732 
1733  // Resize initial guess depending on the rank of the physical points array
1734  FieldContainer<Scalar> initGuess;
1735 
1736  // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array. Requires (C,P,D) initial guess.
1737  if(whichCell == -1){
1738  numPoints = static_cast<size_t>(physPoints.dimension(1));
1739  numCells = static_cast<size_t>(cellWorkset.dimension(0));
1740  initGuess.resize(numCells, numPoints, spaceDim);
1741  // Set initial guess:
1742  for(size_t c = 0; c < numCells; c++){
1743  for(size_t p = 0; p < numPoints; p++){
1744  for(size_t d = 0; d < spaceDim; d++){
1745  initGuess(c, p, d) = cellCenter(d);
1746  }// d
1747  }// p
1748  }// c
1749  }
1750  // Custom: map (P,D) array of physical pts. to (P,D) array. Requires (P,D) initial guess.
1751  else {
1752  numPoints = static_cast<size_t>(physPoints.dimension(0));
1753  initGuess.resize(numPoints, spaceDim);
1754  // Set initial guess:
1755  for(size_t p = 0; p < numPoints; p++){
1756  for(size_t d = 0; d < spaceDim; d++){
1757  initGuess(p, d) = cellCenter(d);
1758  }// d
1759  }// p
1760  }
1761  // Call method with initial guess
1762  mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, cellWorkset, cellTopo, whichCell);
1763 
1764 }
1765 
1766 
1767 template<class Scalar>
1768 template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
1769 void CellTools<Scalar>::mapToReferenceFrameInitGuess(ArrayRefPoint & refPoints,
1770  const ArrayInitGuess & initGuess,
1771  const ArrayPhysPoint & physPoints,
1772  const ArrayCell & cellWorkset,
1773  const Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis,
1774  const int & whichCell)
1775 {
1776 ArrayWrapper<Scalar,ArrayInitGuess, Rank<ArrayInitGuess >::value, true>initGuessWrap(initGuess);
1777 ArrayWrapper<Scalar,ArrayRefPoint, Rank<ArrayRefPoint >::value, false>refPointsWrap(refPoints);
1778 // INTREPID_VALIDATE( validateArguments_mapToReferenceFrame(refPoints, initGuess, physPoints, cellWorkset, cellTopo, whichCell) );
1779  size_t spaceDim = (size_t)HGRAD_Basis->getBaseCellTopology().getDimension();
1780  size_t numPoints;
1781  size_t numCells=0;
1782 
1783  // Temp arrays for Newton iterates and Jacobians. Resize according to rank of ref. point array
1784  FieldContainer<Scalar> xOld;
1785  FieldContainer<Scalar> xTem;
1786  FieldContainer<Scalar> jacobian;
1787  FieldContainer<Scalar> jacobInv;
1788  FieldContainer<Scalar> error;
1789  FieldContainer<Scalar> cellCenter(spaceDim);
1790 
1791  // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array. Requires (C,P,D) temp arrays and (C,P,D,D) Jacobians.
1792  if(whichCell == -1){
1793  numPoints = static_cast<size_t>(physPoints.dimension(1));
1794  numCells = static_cast<size_t>(cellWorkset.dimension(0));
1795  xOld.resize(numCells, numPoints, spaceDim);
1796  xTem.resize(numCells, numPoints, spaceDim);
1797  jacobian.resize(numCells,numPoints, spaceDim, spaceDim);
1798  jacobInv.resize(numCells,numPoints, spaceDim, spaceDim);
1799  error.resize(numCells,numPoints);
1800  // Set initial guess to xOld
1801  for(size_t c = 0; c < numCells; c++){
1802  for(size_t p = 0; p < numPoints; p++){
1803  for(size_t d = 0; d < spaceDim; d++){
1804  xOld(c, p, d) = initGuessWrap(c, p, d);
1805  }// d
1806  }// p
1807  }// c
1808  }
1809  // Custom: map (P,D) array of physical pts. to (P,D) array. Requires (P,D) temp arrays and (P,D,D) Jacobians.
1810  else {
1811  numPoints = static_cast<size_t>(physPoints.dimension(0));
1812  xOld.resize(numPoints, spaceDim);
1813  xTem.resize(numPoints, spaceDim);
1814  jacobian.resize(numPoints, spaceDim, spaceDim);
1815  jacobInv.resize(numPoints, spaceDim, spaceDim);
1816  error.resize(numPoints);
1817  // Set initial guess to xOld
1818  for(size_t c = 0; c < numCells; c++){
1819  for(size_t p = 0; p < numPoints; p++){
1820  for(size_t d = 0; d < spaceDim; d++){
1821  xOld(c, p, d) = initGuessWrap(c, p, d);
1822  }// d
1823  }// p
1824  }// c
1825  }
1826 
1827  // Newton method to solve the equation F(refPoints) - physPoints = 0:
1828  // refPoints = xOld - DF^{-1}(xOld)*(F(xOld) - physPoints) = xOld + DF^{-1}(xOld)*(physPoints - F(xOld))
1829  for(int iter = 0; iter < INTREPID_MAX_NEWTON; ++iter) {
1830 
1831  // Jacobians at the old iterates and their inverses.
1832  setJacobian(jacobian, xOld, cellWorkset, HGRAD_Basis, whichCell);
1833  setJacobianInv(jacobInv, jacobian);
1834  // The Newton step.
1835  mapToPhysicalFrame( xTem, xOld, cellWorkset, HGRAD_Basis->getBaseCellTopology(), whichCell ); // xTem <- F(xOld)
1836  RealSpaceTools<Scalar>::subtract( xTem, physPoints, xTem ); // xTem <- physPoints - F(xOld)
1837  RealSpaceTools<Scalar>::matvec( refPoints, jacobInv, xTem); // refPoints <- DF^{-1}( physPoints - F(xOld) )
1838  RealSpaceTools<Scalar>::add( refPoints, xOld ); // refPoints <- DF^{-1}( physPoints - F(xOld) ) + xOld
1839 
1840  // l2 error (Euclidean distance) between old and new iterates: |xOld - xNew|
1841  RealSpaceTools<Scalar>::subtract( xTem, xOld, refPoints );
1842  RealSpaceTools<Scalar>::vectorNorm( error, xTem, NORM_TWO );
1843 
1844  // Average L2 error for a multiple sets of physical points: error is rank-2 (C,P) array
1845  Scalar totalError;
1846  if(whichCell == -1) {
1847  FieldContainer<Scalar> cellWiseError(numCells);
1848  // error(C,P) -> cellWiseError(P)
1849 
1850  RealSpaceTools<Scalar>::vectorNorm( cellWiseError, error, NORM_ONE );
1851  totalError = RealSpaceTools<Scalar>::vectorNorm( cellWiseError, NORM_ONE );
1852  }
1853  //Average L2 error for a single set of physical points: error is rank-1 (P) array
1854  else{
1855 
1856  totalError = RealSpaceTools<Scalar>::vectorNorm( error, NORM_ONE );
1857  totalError = totalError;
1858  }
1859 
1860  // Stopping criterion:
1861  if (totalError < INTREPID_TOL) {
1862  break;
1863  }
1864  else if ( iter > INTREPID_MAX_NEWTON) {
1865  INTREPID_VALIDATE(std::cout << " Intrepid::CellTools::mapToReferenceFrameInitGuess failed to converge to desired tolerance within "
1866  << INTREPID_MAX_NEWTON << " iterations\n" );
1867  break;
1868  }
1869 
1870  // initialize next Newton step
1871 // xOld = refPoints;
1872 int refPointsRank=getrank(refPoints);
1873 if (refPointsRank==3){
1874  for(size_t i=0;i<static_cast<size_t>(refPoints.dimension(0));i++){
1875  for(size_t j=0;j<static_cast<size_t>(refPoints.dimension(1));j++){
1876  for(size_t k=0;k<static_cast<size_t>(refPoints.dimension(2));k++){
1877  xOld(i,j,k) = refPointsWrap(i,j,k);
1878  }
1879  }
1880  }
1881 }else if(refPointsRank==2){
1882  for(size_t i=0;i<static_cast<size_t>(refPoints.dimension(0));i++){
1883  for(size_t j=0;j<static_cast<size_t>(refPoints.dimension(1));j++){
1884  xOld(i,j) = refPointsWrap(i,j);
1885  }
1886  }
1887 
1888 }
1889 
1890 
1891 
1892  } // for(iter)
1893 }
1894 
1895 
1896 
1897 template<class Scalar>
1898 template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
1899 void CellTools<Scalar>::mapToReferenceFrameInitGuess(ArrayRefPoint & refPoints,
1900  const ArrayInitGuess & initGuess,
1901  const ArrayPhysPoint & physPoints,
1902  const ArrayCell & cellWorkset,
1903  const shards::CellTopology & cellTopo,
1904  const int & whichCell)
1905 {
1906 ArrayWrapper<Scalar,ArrayInitGuess, Rank<ArrayInitGuess >::value, true>initGuessWrap(initGuess);
1907 ArrayWrapper<Scalar,ArrayRefPoint, Rank<ArrayRefPoint >::value, false>refPointsWrap(refPoints);
1908  INTREPID_VALIDATE( validateArguments_mapToReferenceFrame(refPoints, initGuess, physPoints, cellWorkset, cellTopo, whichCell) );
1909  size_t spaceDim = (size_t)cellTopo.getDimension();
1910  size_t numPoints;
1911  size_t numCells=0;
1912 
1913  // Temp arrays for Newton iterates and Jacobians. Resize according to rank of ref. point array
1914  FieldContainer<Scalar> xOld;
1915  FieldContainer<Scalar> xTem;
1916  FieldContainer<Scalar> jacobian;
1917  FieldContainer<Scalar> jacobInv;
1918  FieldContainer<Scalar> error;
1919  FieldContainer<Scalar> cellCenter(spaceDim);
1920 
1921  // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array. Requires (C,P,D) temp arrays and (C,P,D,D) Jacobians.
1922  if(whichCell == -1){
1923  numPoints = static_cast<size_t>(physPoints.dimension(1));
1924  numCells = static_cast<size_t>(cellWorkset.dimension(0));
1925  xOld.resize(numCells, numPoints, spaceDim);
1926  xTem.resize(numCells, numPoints, spaceDim);
1927  jacobian.resize(numCells,numPoints, spaceDim, spaceDim);
1928  jacobInv.resize(numCells,numPoints, spaceDim, spaceDim);
1929  error.resize(numCells,numPoints);
1930  // Set initial guess to xOld
1931  for(size_t c = 0; c < numCells; c++){
1932  for(size_t p = 0; p < numPoints; p++){
1933  for(size_t d = 0; d < spaceDim; d++){
1934  xOld(c, p, d) = initGuessWrap(c, p, d);
1935  }// d
1936  }// p
1937  }// c
1938  }
1939  // Custom: map (P,D) array of physical pts. to (P,D) array. Requires (P,D) temp arrays and (P,D,D) Jacobians.
1940  else {
1941  numPoints = static_cast<size_t>(physPoints.dimension(0));
1942  xOld.resize(numPoints, spaceDim);
1943  xTem.resize(numPoints, spaceDim);
1944  jacobian.resize(numPoints, spaceDim, spaceDim);
1945  jacobInv.resize(numPoints, spaceDim, spaceDim);
1946  error.resize(numPoints);
1947  // Set initial guess to xOld
1948  for(size_t p = 0; p < numPoints; p++){
1949  for(size_t d = 0; d < spaceDim; d++){
1950  xOld(p, d) = initGuessWrap(p, d);
1951  }// d
1952  }// p
1953  }
1954 
1955  // Newton method to solve the equation F(refPoints) - physPoints = 0:
1956  // refPoints = xOld - DF^{-1}(xOld)*(F(xOld) - physPoints) = xOld + DF^{-1}(xOld)*(physPoints - F(xOld))
1957  for(int iter = 0; iter < INTREPID_MAX_NEWTON; ++iter) {
1958 
1959  // Jacobians at the old iterates and their inverses.
1960  setJacobian(jacobian, xOld, cellWorkset, cellTopo, whichCell);
1961  setJacobianInv(jacobInv, jacobian);
1962  // The Newton step.
1963  mapToPhysicalFrame( xTem, xOld, cellWorkset, cellTopo, whichCell ); // xTem <- F(xOld)
1964  RealSpaceTools<Scalar>::subtract( xTem, physPoints, xTem ); // xTem <- physPoints - F(xOld)
1965  RealSpaceTools<Scalar>::matvec( refPoints, jacobInv, xTem); // refPoints <- DF^{-1}( physPoints - F(xOld) )
1966  RealSpaceTools<Scalar>::add( refPoints, xOld ); // refPoints <- DF^{-1}( physPoints - F(xOld) ) + xOld
1967 
1968  // l2 error (Euclidean distance) between old and new iterates: |xOld - xNew|
1969  RealSpaceTools<Scalar>::subtract( xTem, xOld, refPoints );
1970  RealSpaceTools<Scalar>::vectorNorm( error, xTem, NORM_TWO );
1971 
1972  // Average L2 error for a multiple sets of physical points: error is rank-2 (C,P) array
1973  Scalar totalError;
1974  if(whichCell == -1) {
1975  FieldContainer<Scalar> cellWiseError(numCells);
1976  // error(C,P) -> cellWiseError(P)
1977 
1978  RealSpaceTools<Scalar>::vectorNorm( cellWiseError, error, NORM_ONE );
1979  totalError = RealSpaceTools<Scalar>::vectorNorm( cellWiseError, NORM_ONE );
1980  }
1981  //Average L2 error for a single set of physical points: error is rank-1 (P) array
1982  else{
1983 
1984  totalError = RealSpaceTools<Scalar>::vectorNorm( error, NORM_ONE );
1985  totalError = totalError;
1986  }
1987 
1988  // Stopping criterion:
1989  if (totalError < INTREPID_TOL) {
1990  break;
1991  }
1992  else if ( iter > INTREPID_MAX_NEWTON) {
1993  INTREPID_VALIDATE(std::cout << " Intrepid::CellTools::mapToReferenceFrameInitGuess failed to converge to desired tolerance within "
1994  << INTREPID_MAX_NEWTON << " iterations\n" );
1995  break;
1996  }
1997 
1998  // initialize next Newton step
1999 // xOld = refPoints;
2000 int refPointsRank=getrank(refPoints);
2001 if (refPointsRank==3){
2002  for(size_t i=0;i<static_cast<size_t>(refPoints.dimension(0));i++){
2003  for(size_t j=0;j<static_cast<size_t>(refPoints.dimension(1));j++){
2004  for(size_t k=0;k<static_cast<size_t>(refPoints.dimension(2));k++){
2005  xOld(i,j,k) = refPointsWrap(i,j,k);
2006  }
2007  }
2008  }
2009 }else if(refPointsRank==2){
2010  for(size_t i=0;i<static_cast<size_t>(refPoints.dimension(0));i++){
2011  for(size_t j=0;j<static_cast<size_t>(refPoints.dimension(1));j++){
2012  xOld(i,j) = refPointsWrap(i,j);
2013  }
2014  }
2015 
2016 }
2017 
2018 
2019 
2020  } // for(iter)
2021 }
2022 
2023 
2024 
2025 template<class Scalar>
2026 template<class ArraySubcellPoint, class ArrayParamPoint>
2027 void CellTools<Scalar>::mapToReferenceSubcell(ArraySubcellPoint & refSubcellPoints,
2028  const ArrayParamPoint & paramPoints,
2029  const int subcellDim,
2030  const int subcellOrd,
2031  const shards::CellTopology & parentCell){
2032 
2033  int cellDim = parentCell.getDimension();
2034  size_t numPts = static_cast<size_t>(paramPoints.dimension(0));
2035 
2036 #ifdef HAVE_INTREPID_DEBUG
2037  TEUCHOS_TEST_FOR_EXCEPTION( !(hasReferenceCell(parentCell) ), std::invalid_argument,
2038  ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
2039 
2040  TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= subcellDim) && (subcellDim <= 2 ) ), std::invalid_argument,
2041  ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-dimensional subcells.");
2042 
2043  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= subcellOrd) && (subcellOrd < (int)parentCell.getSubcellCount(subcellDim) ) ), std::invalid_argument,
2044  ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
2045 
2046  // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
2047  std::string errmsg = ">>> ERROR (Intrepid::mapToReferenceSubcell):";
2048  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, refSubcellPoints, 2,2), std::invalid_argument, errmsg);
2049  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, refSubcellPoints, 1, cellDim, cellDim), std::invalid_argument, errmsg);
2050 
2051  // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
2052  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, paramPoints, 2,2), std::invalid_argument, errmsg);
2053  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, paramPoints, 1, subcellDim, subcellDim), std::invalid_argument, errmsg);
2054 
2055  // cross check: refSubcellPoints and paramPoints: dimension 0 must match
2056  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, refSubcellPoints, 0, paramPoints, 0), std::invalid_argument, errmsg);
2057 #endif
2058 
2059 
2060  // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
2061  const FieldContainer<double>& subcellMap = getSubcellParametrization(subcellDim, parentCell);
2062 
2063  // Apply the parametrization map to every point in parameter domain
2064  if(subcellDim == 2) {
2065  for(size_t pt = 0; pt < numPts; pt++){
2066  double u = paramPoints(pt,0);
2067  double v = paramPoints(pt,1);
2068 
2069  // map_dim(u,v) = c_0(dim) + c_1(dim)*u + c_2(dim)*v because both Quad and Tri ref faces are affine!
2070  for(int dim = 0; dim < cellDim; dim++){
2071  refSubcellPoints(pt, dim) = subcellMap(subcellOrd, dim, 0) + \
2072  subcellMap(subcellOrd, dim, 1)*u + \
2073  subcellMap(subcellOrd, dim, 2)*v;
2074  }
2075  }
2076  }
2077  else if(subcellDim == 1) {
2078  for(size_t pt = 0; pt < numPts; pt++){
2079  for(int dim = 0; dim < cellDim; dim++) {
2080  refSubcellPoints(pt, dim) = subcellMap(subcellOrd, dim, 0) + subcellMap(subcellOrd, dim, 1)*paramPoints(pt,0);
2081  }
2082  }
2083  }
2084  else{
2085  TEUCHOS_TEST_FOR_EXCEPTION( !( (subcellDim == 1) || (subcellDim == 2) ), std::invalid_argument,
2086  ">>> ERROR (Intrepid::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-subcells");
2087  }
2088 }
2089 
2090 
2091 
2092 template<class Scalar>
2093 template<class ArrayEdgeTangent>
2094 void CellTools<Scalar>::getReferenceEdgeTangent(ArrayEdgeTangent & refEdgeTangent,
2095  const int & edgeOrd,
2096  const shards::CellTopology & parentCell){
2097 
2098  int spaceDim = parentCell.getDimension();
2099 
2100 #ifdef HAVE_INTREPID_DEBUG
2101 
2102  TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument,
2103  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): two or three-dimensional parent cell required");
2104 
2105  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= edgeOrd) && (edgeOrd < (int)parentCell.getSubcellCount(1) ) ), std::invalid_argument,
2106  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): edge ordinal out of bounds");
2107 
2108  TEUCHOS_TEST_FOR_EXCEPTION( !( refEdgeTangent.size() == spaceDim ), std::invalid_argument,
2109  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): output array size is required to match space dimension");
2110 #endif
2111  // Edge parametrizations are computed in setSubcellParametrization and stored in rank-3 array
2112  // (subcOrd, coordinate, coefficient)
2113  const FieldContainer<double>& edgeMap = getSubcellParametrization(1, parentCell);
2114 
2115  // All ref. edge maps have affine coordinate functions: f_dim(u) = C_0(dim) + C_1(dim)*u,
2116  // => edge Tangent: -> C_1(*)
2117  refEdgeTangent(0) = edgeMap(edgeOrd, 0, 1);
2118  refEdgeTangent(1) = edgeMap(edgeOrd, 1, 1);
2119 
2120  // Skip last coordinate for 2D parent cells
2121  if(spaceDim == 3) {
2122  refEdgeTangent(2) = edgeMap(edgeOrd, 2, 1);
2123  }
2124 }
2125 
2126 
2127 
2128 template<class Scalar>
2129 template<class ArrayFaceTangentU, class ArrayFaceTangentV>
2130 void CellTools<Scalar>::getReferenceFaceTangents(ArrayFaceTangentU & uTan,
2131  ArrayFaceTangentV & vTan,
2132  const int & faceOrd,
2133  const shards::CellTopology & parentCell){
2134 
2135 #ifdef HAVE_INTREPID_DEBUG
2136  int spaceDim = parentCell.getDimension();
2137  TEUCHOS_TEST_FOR_EXCEPTION( !(spaceDim == 3), std::invalid_argument,
2138  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): three-dimensional parent cell required");
2139 
2140  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= faceOrd) && (faceOrd < (int)parentCell.getSubcellCount(2) ) ), std::invalid_argument,
2141  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): face ordinal out of bounds");
2142 
2143  TEUCHOS_TEST_FOR_EXCEPTION( !( (getrank(uTan) == 1) && (getrank(vTan) == 1) ), std::invalid_argument,
2144  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): rank = 1 required for output arrays");
2145 
2146  TEUCHOS_TEST_FOR_EXCEPTION( !( uTan.dimension(0) == spaceDim ), std::invalid_argument,
2147  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");
2148 
2149  TEUCHOS_TEST_FOR_EXCEPTION( !( vTan.dimension(0) == spaceDim ), std::invalid_argument,
2150  ">>> ERROR (Intrepid::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");
2151 #endif
2152 
2153  // Face parametrizations are computed in setSubcellParametrization and stored in rank-3 array
2154  // (subcOrd, coordinate, coefficient): retrieve this array
2155  const FieldContainer<double>& faceMap = getSubcellParametrization(2, parentCell);
2156 
2157  /* All ref. face maps have affine coordinate functions: f_dim(u,v) = C_0(dim) + C_1(dim)*u + C_2(dim)*v
2158  * ` => Tangent vectors are: uTan -> C_1(*); vTan -> C_2(*)
2159  */
2160  // set uTan -> C_1(*)
2161  uTan(0) = faceMap(faceOrd, 0, 1);
2162  uTan(1) = faceMap(faceOrd, 1, 1);
2163  uTan(2) = faceMap(faceOrd, 2, 1);
2164 
2165  // set vTan -> C_2(*)
2166  vTan(0) = faceMap(faceOrd, 0, 2);
2167  vTan(1) = faceMap(faceOrd, 1, 2);
2168  vTan(2) = faceMap(faceOrd, 2, 2);
2169 }
2170 
2171 
2172 
2173 template<class Scalar>
2174 template<class ArraySideNormal>
2175 void CellTools<Scalar>::getReferenceSideNormal(ArraySideNormal & refSideNormal,
2176  const int & sideOrd,
2177  const shards::CellTopology & parentCell){
2178  int spaceDim = parentCell.getDimension();
2179 
2180  #ifdef HAVE_INTREPID_DEBUG
2181 
2182  TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument,
2183  ">>> ERROR (Intrepid::CellTools::getReferenceSideNormal): two or three-dimensional parent cell required");
2184 
2185  // Check side ordinal: by definition side is subcell whose dimension = spaceDim-1
2186  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= sideOrd) && (sideOrd < (int)parentCell.getSubcellCount(spaceDim - 1) ) ), std::invalid_argument,
2187  ">>> ERROR (Intrepid::CellTools::getReferenceSideNormal): side ordinal out of bounds");
2188 #endif
2189  if(spaceDim == 2){
2190 
2191  // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
2192  getReferenceEdgeTangent(refSideNormal, sideOrd, parentCell);
2193 
2194  // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
2195  Scalar temp = refSideNormal(0);
2196  refSideNormal(0) = refSideNormal(1);
2197  refSideNormal(1) = -temp;
2198  }
2199  else{
2200  // 3D parent cell: side = 2D subcell (face), call the face normal method.
2201  getReferenceFaceNormal(refSideNormal, sideOrd, parentCell);
2202  }
2203 }
2204 
2205 
2206 
2207 template<class Scalar>
2208 template<class ArrayFaceNormal>
2209 void CellTools<Scalar>::getReferenceFaceNormal(ArrayFaceNormal & refFaceNormal,
2210  const int & faceOrd,
2211  const shards::CellTopology & parentCell){
2212  int spaceDim = parentCell.getDimension();
2213  #ifdef HAVE_INTREPID_DEBUG
2214 
2215  TEUCHOS_TEST_FOR_EXCEPTION( !(spaceDim == 3), std::invalid_argument,
2216  ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): three-dimensional parent cell required");
2217 
2218  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= faceOrd) && (faceOrd < (int)parentCell.getSubcellCount(2) ) ), std::invalid_argument,
2219  ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): face ordinal out of bounds");
2220 
2221  TEUCHOS_TEST_FOR_EXCEPTION( !( getrank(refFaceNormal) == 1 ), std::invalid_argument,
2222  ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): rank = 1 required for output array");
2223 
2224  TEUCHOS_TEST_FOR_EXCEPTION( !( static_cast<size_t>(refFaceNormal.dimension(0)) == static_cast<size_t>(spaceDim) ), std::invalid_argument,
2225  ">>> ERROR (Intrepid::CellTools::getReferenceFaceNormal): dim0 (spatial dim) must match that of parent cell");
2226 #endif
2227 
2228  // Reference face normal = vector product of reference face tangents. Allocate temp FC storage:
2229  FieldContainer<Scalar> uTan(spaceDim);
2230  FieldContainer<Scalar> vTan(spaceDim);
2231  getReferenceFaceTangents(uTan, vTan, faceOrd, parentCell);
2232 
2233  // Compute the vector product of the reference face tangents:
2234  RealSpaceTools<Scalar>::vecprod(refFaceNormal, uTan, vTan);
2235 }
2236 
2237 template<class Scalar>
2238 template<class ArrayEdgeTangent, class ArrayJac>
2239 void CellTools<Scalar>::getPhysicalEdgeTangents(ArrayEdgeTangent & edgeTangents,
2240  const ArrayJac & worksetJacobians,
2241  const int & worksetEdgeOrd,
2242  const shards::CellTopology & parentCell){
2243  size_t worksetSize = static_cast<size_t>(worksetJacobians.dimension(0));
2244  size_t edgePtCount = static_cast<size_t>(worksetJacobians.dimension(1));
2245  int pCellDim = parentCell.getDimension();
2246  #ifdef HAVE_INTREPID_DEBUG
2247  std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalEdgeTangents):";
2248 
2249  TEUCHOS_TEST_FOR_EXCEPTION( !( (pCellDim == 3) || (pCellDim == 2) ), std::invalid_argument,
2250  ">>> ERROR (Intrepid::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required");
2251 
2252  // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
2253  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, edgeTangents, 3,3), std::invalid_argument, errmsg);
2254  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, edgeTangents, 2, 2,3), std::invalid_argument, errmsg);
2255 
2256  // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
2257  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
2258  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 2,3), std::invalid_argument, errmsg);
2259  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 2,3), std::invalid_argument, errmsg);
2260 
2261  // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
2262  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, edgeTangents, 0,1,2,2, worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);
2263 
2264 #endif
2265 
2266  // Temp storage for constant reference edge tangent: rank-1 (D) arrays
2267  FieldContainer<double> refEdgeTan(pCellDim);
2268  getReferenceEdgeTangent(refEdgeTan, worksetEdgeOrd, parentCell);
2269 
2270  // Loop over workset faces and edge points
2271  for(size_t pCell = 0; pCell < worksetSize; pCell++){
2272  for(size_t pt = 0; pt < edgePtCount; pt++){
2273 
2274  // Apply parent cell Jacobian to ref. edge tangent
2275  for(int i = 0; i < pCellDim; i++){
2276  edgeTangents(pCell, pt, i) = 0.0;
2277  for(int j = 0; j < pCellDim; j++){
2278  edgeTangents(pCell, pt, i) += worksetJacobians(pCell, pt, i, j)*refEdgeTan(j);
2279  }// for j
2280  }// for i
2281  }// for pt
2282  }// for pCell
2283 }
2284 template<class Scalar>
2285 template<class ArrayFaceTangentU, class ArrayFaceTangentV, class ArrayJac>
2286 void CellTools<Scalar>::getPhysicalFaceTangents(ArrayFaceTangentU & faceTanU,
2287  ArrayFaceTangentV & faceTanV,
2288  const ArrayJac & worksetJacobians,
2289  const int & worksetFaceOrd,
2290  const shards::CellTopology & parentCell){
2291  size_t worksetSize = static_cast<size_t>(worksetJacobians.dimension(0));
2292  size_t facePtCount = static_cast<size_t>(worksetJacobians.dimension(1));
2293  int pCellDim = parentCell.getDimension();
2294  #ifdef HAVE_INTREPID_DEBUG
2295  std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalFaceTangents):";
2296 
2297  TEUCHOS_TEST_FOR_EXCEPTION( !(pCellDim == 3), std::invalid_argument,
2298  ">>> ERROR (Intrepid::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");
2299 
2300  // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
2301  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceTanU, 3,3), std::invalid_argument, errmsg);
2302  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceTanU, 2, 3,3), std::invalid_argument, errmsg);
2303 
2304  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceTanV, 3,3), std::invalid_argument, errmsg);
2305  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceTanV, 2, 3,3), std::invalid_argument, errmsg);
2306 
2307  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceTanU, faceTanV), std::invalid_argument, errmsg);
2308 
2309  // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
2310  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
2311  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 3,3), std::invalid_argument, errmsg);
2312  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 3,3), std::invalid_argument, errmsg);
2313 
2314  // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
2315  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceTanU, 0,1,2,2, worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);
2316 
2317 #endif
2318 
2319  // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
2320  FieldContainer<double> refFaceTanU(pCellDim);
2321  FieldContainer<double> refFaceTanV(pCellDim);
2322  getReferenceFaceTangents(refFaceTanU, refFaceTanV, worksetFaceOrd, parentCell);
2323 
2324  // Loop over workset faces and face points
2325  for(size_t pCell = 0; pCell < worksetSize; pCell++){
2326  for(size_t pt = 0; pt < facePtCount; pt++){
2327 
2328  // Apply parent cell Jacobian to ref. face tangents
2329  for(int dim = 0; dim < pCellDim; dim++){
2330  faceTanU(pCell, pt, dim) = 0.0;
2331  faceTanV(pCell, pt, dim) = 0.0;
2332 
2333  // Unroll loops: parent cell dimension can only be 3
2334  faceTanU(pCell, pt, dim) = \
2335  worksetJacobians(pCell, pt, dim, 0)*refFaceTanU(0) + \
2336  worksetJacobians(pCell, pt, dim, 1)*refFaceTanU(1) + \
2337  worksetJacobians(pCell, pt, dim, 2)*refFaceTanU(2);
2338  faceTanV(pCell, pt, dim) = \
2339  worksetJacobians(pCell, pt, dim, 0)*refFaceTanV(0) + \
2340  worksetJacobians(pCell, pt, dim, 1)*refFaceTanV(1) + \
2341  worksetJacobians(pCell, pt, dim, 2)*refFaceTanV(2);
2342  }// for dim
2343  }// for pt
2344  }// for pCell
2345 }
2346 
2347 template<class Scalar>
2348 template<class ArraySideNormal, class ArrayJac>
2349 void CellTools<Scalar>::getPhysicalSideNormals(ArraySideNormal & sideNormals,
2350  const ArrayJac & worksetJacobians,
2351  const int & worksetSideOrd,
2352  const shards::CellTopology & parentCell){
2353  size_t worksetSize = static_cast<size_t>(worksetJacobians.dimension(0));
2354  size_t sidePtCount = static_cast<size_t>(worksetJacobians.dimension(1));
2355  int spaceDim = parentCell.getDimension();
2356  #ifdef HAVE_INTREPID_DEBUG
2357  TEUCHOS_TEST_FOR_EXCEPTION( !( (spaceDim == 2) || (spaceDim == 3) ), std::invalid_argument,
2358  ">>> ERROR (Intrepid::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
2359 
2360  // Check side ordinal: by definition side is subcell whose dimension = spaceDim-1
2361  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 <= worksetSideOrd) && (worksetSideOrd < (int)parentCell.getSubcellCount(spaceDim - 1) ) ), std::invalid_argument,
2362  ">>> ERROR (Intrepid::CellTools::getPhysicalSideNormals): side ordinal out of bounds");
2363 #endif
2364 
2365  if(spaceDim == 2){
2366 
2367  // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
2368  getPhysicalEdgeTangents(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
2369 
2370  // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
2371  for(size_t cell = 0; cell < worksetSize; cell++){
2372  for(size_t pt = 0; pt < sidePtCount; pt++){
2373  Scalar temp = sideNormals(cell, pt, 0);
2374  sideNormals(cell, pt, 0) = sideNormals(cell, pt, 1);
2375  sideNormals(cell, pt, 1) = -temp;
2376  }// for pt
2377  }// for cell
2378  }
2379  else{
2380  // 3D parent cell: side = 2D subcell (face), call the face normal method.
2381  getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
2382  }
2383 }
2384 
2385 
2386 template<class Scalar>
2387 template<class ArrayFaceNormal, class ArrayJac>
2388 void CellTools<Scalar>::getPhysicalFaceNormals(ArrayFaceNormal & faceNormals,
2389  const ArrayJac & worksetJacobians,
2390  const int & worksetFaceOrd,
2391  const shards::CellTopology & parentCell){
2392  size_t worksetSize = static_cast<size_t>(worksetJacobians.dimension(0));
2393  size_t facePtCount = static_cast<size_t>(worksetJacobians.dimension(1));
2394  int pCellDim = parentCell.getDimension();
2395  #ifdef HAVE_INTREPID_DEBUG
2396  std::string errmsg = ">>> ERROR (Intrepid::CellTools::getPhysicalFaceNormals):";
2397 
2398  TEUCHOS_TEST_FOR_EXCEPTION( !(pCellDim == 3), std::invalid_argument,
2399  ">>> ERROR (Intrepid::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required");
2400 
2401  // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
2402  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, faceNormals, 3,3), std::invalid_argument, errmsg);
2403  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, faceNormals, 2, 3,3), std::invalid_argument, errmsg);
2404 
2405  // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
2406  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, worksetJacobians, 4,4), std::invalid_argument, errmsg);
2407  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 2, 3,3), std::invalid_argument, errmsg);
2408  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, worksetJacobians, 3, 3,3), std::invalid_argument, errmsg);
2409 
2410  // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
2411  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, faceNormals, 0,1,2,2, worksetJacobians, 0,1,2,3), std::invalid_argument, errmsg);
2412 #endif
2413 
2414  // Temp storage for physical face tangents: rank-3 (C,P,D) arrays
2415  FieldContainer<Scalar> faceTanU(worksetSize, facePtCount, pCellDim);
2416  FieldContainer<Scalar> faceTanV(worksetSize, facePtCount, pCellDim);
2417  getPhysicalFaceTangents(faceTanU, faceTanV, worksetJacobians, worksetFaceOrd, parentCell);
2418 
2419  // Compute the vector product of the physical face tangents:
2420  RealSpaceTools<Scalar>::vecprod(faceNormals, faceTanU, faceTanV);
2421 
2422 
2423 }
2424 //============================================================================================//
2425 // //
2426 // Inclusion tests //
2427 // //
2428 //============================================================================================//
2429 
2430 
2431 template<class Scalar>
2432 int CellTools<Scalar>::checkPointInclusion(const Scalar* point,
2433  const int pointDim,
2434  const shards::CellTopology & cellTopo,
2435  const double & threshold) {
2436 #ifdef HAVE_INTREPID_DEBUG
2437  TEUCHOS_TEST_FOR_EXCEPTION( !(pointDim == (int)cellTopo.getDimension() ), std::invalid_argument,
2438  ">>> ERROR (Intrepid::CellTools::checkPointInclusion): Point and cell dimensions do not match. ");
2439 #endif
2440  int testResult = 1;
2441 
2442  // Using these values in the tests effectievly inflates the reference element to a larger one
2443  double minus_one = -1.0 - threshold;
2444  double plus_one = 1.0 + threshold;
2445  double minus_zero = - threshold;
2446 
2447  // A cell with extended topology has the same reference cell as a cell with base topology.
2448  // => testing for inclusion in a reference Triangle<> and a reference Triangle<6> relies on
2449  // on the same set of inequalities. To eliminate unnecessary cases we switch on the base topology
2450  unsigned key = cellTopo.getBaseCellTopologyData() -> key ;
2451  switch( key ) {
2452 
2453  case shards::Line<>::key :
2454  if( !(minus_one <= point[0] && point[0] <= plus_one)) testResult = 0;
2455  break;
2456 
2457  case shards::Triangle<>::key : {
2458  Scalar distance = std::max( std::max( -point[0], -point[1] ), point[0] + point[1] - 1.0 );
2459  if( distance > threshold ) testResult = 0;
2460  break;
2461  }
2462 
2463  case shards::Quadrilateral<>::key :
2464  if(!( (minus_one <= point[0] && point[0] <= plus_one) && \
2465  (minus_one <= point[1] && point[1] <= plus_one) ) ) testResult = 0;
2466  break;
2467 
2468  case shards::Tetrahedron<>::key : {
2469  Scalar distance = std::max( std::max(-point[0],-point[1]), \
2470  std::max(-point[2], point[0] + point[1] + point[2] - 1) );
2471  if( distance > threshold ) testResult = 0;
2472  break;
2473  }
2474 
2475  case shards::Hexahedron<>::key :
2476  if(!((minus_one <= point[0] && point[0] <= plus_one) && \
2477  (minus_one <= point[1] && point[1] <= plus_one) && \
2478  (minus_one <= point[2] && point[2] <= plus_one))) \
2479  testResult = 0;
2480  break;
2481 
2482  // The base of the reference prism is the same as the reference triangle => apply triangle test
2483  // to X and Y coordinates and test whether Z is in [-1,1]
2484  case shards::Wedge<>::key : {
2485  Scalar distance = std::max( std::max( -point[0], -point[1] ), point[0] + point[1] - 1 );
2486  if( distance > threshold || \
2487  point[2] < minus_one || point[2] > plus_one) \
2488  testResult = 0;
2489  break;
2490  }
2491 
2492  // The base of the reference pyramid is the same as the reference quad cell => a horizontal plane
2493  // through a point P(x,y,z) intersects the pyramid at a quadrilateral that equals the base quad
2494  // scaled by (1-z) => P(x,y,z) is inside the pyramid <=> (x,y) is in [-1+z,1-z]^2 && 0 <= Z <= 1
2495  case shards::Pyramid<>::key : {
2496  Scalar left = minus_one + point[2];
2497  Scalar right = plus_one - point[2];
2498  if(!( (left <= point[0] && point[0] <= right) && \
2499  (left <= point[1] && point[1] <= right) &&
2500  (minus_zero <= point[2] && point[2] <= plus_one) ) ) \
2501  testResult = 0;
2502  break;
2503  }
2504 
2505  default:
2506  TEUCHOS_TEST_FOR_EXCEPTION( !( (key == shards::Line<>::key ) ||
2507  (key == shards::Triangle<>::key) ||
2508  (key == shards::Quadrilateral<>::key) ||
2509  (key == shards::Tetrahedron<>::key) ||
2510  (key == shards::Hexahedron<>::key) ||
2511  (key == shards::Wedge<>::key) ||
2512  (key == shards::Pyramid<>::key) ),
2513  std::invalid_argument,
2514  ">>> ERROR (Intrepid::CellTools::checkPointInclusion): Invalid cell topology. ");
2515  }
2516  return testResult;
2517 }
2518 
2519 
2520 
2521 template<class Scalar>
2522 template<class ArrayPoint>
2523 int CellTools<Scalar>::checkPointsetInclusion(const ArrayPoint& points,
2524  const shards::CellTopology & cellTopo,
2525  const double & threshold) {
2526 
2527  int rank = points.rank();
2528 
2529 #ifdef HAVE_INTREPID_DEBUG
2530  TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <=getrank(points) ) && (getrank(points) <= 3) ), std::invalid_argument,
2531  ">>> ERROR (Intrepid::CellTools::checkPointsetInclusion): rank-1, 2 or 3 required for input points array. ");
2532 
2533  // The last dimension of points array at (rank - 1) is the spatial dimension. Must equal the cell dimension.
2534  TEUCHOS_TEST_FOR_EXCEPTION( !((size_t) points.dimension(rank - 1) == (size_t)cellTopo.getDimension() ), std::invalid_argument,
2535  ">>> ERROR (Intrepid::CellTools::checkPointsetInclusion): Point and cell dimensions do not match. ");
2536 #endif
2537 
2538  // create temp output array depending on the rank of the input array
2539  FieldContainer<int> inRefCell;
2540  switch(rank) {
2541  case 1: inRefCell.resize(1); break;
2542  case 2: inRefCell.resize( static_cast<size_t>(points.dimension(0)) ); break;
2543  case 3: inRefCell.resize( static_cast<size_t>(points.dimension(0)), static_cast<size_t>(points.dimension(1)) ); break;
2544  }
2545 
2546  // Call the inclusion method which returns inclusion results for all points
2547  checkPointwiseInclusion(inRefCell, points, cellTopo, threshold);
2548 
2549  // Check if any points were outside, break when finding the first one
2550  int allInside = 1;
2551  for(int i = 0; i < inRefCell.size(); i++ ){
2552  if (inRefCell[i] == 0) {
2553  allInside = 0;
2554  break;
2555  }
2556  }
2557  return allInside;
2558 }
2559 
2560 
2561 
2562 template<class Scalar>
2563 template<class ArrayIncl, class ArrayPoint>
2564 void CellTools<Scalar>::checkPointwiseInclusion(ArrayIncl & inRefCell,
2565  const ArrayPoint & points,
2566  const shards::CellTopology & cellTopo,
2567  const double & threshold) {
2568  int apRank = points.rank();
2569 
2570 #ifdef HAVE_INTREPID_DEBUG
2571 
2572  // Verify that points and inRefCell have correct ranks and dimensions
2573  std::string errmsg = ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion):";
2574  if(getrank(points) == 1) {
2575  TEUCHOS_TEST_FOR_EXCEPTION( !(getrank(inRefCell) == 1 ), std::invalid_argument,
2576  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1 input array requires rank-1 output array.");
2577  TEUCHOS_TEST_FOR_EXCEPTION( !(static_cast<size_t>(inRefCell.dimension(0)) == 1), std::invalid_argument,
2578  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1 input array requires dim0 = 1 for output array.");
2579  }
2580  else if(getrank(points) == 2){
2581  TEUCHOS_TEST_FOR_EXCEPTION( !(getrank(inRefCell) == 1 ), std::invalid_argument,
2582  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-2 input array requires rank-1 output array.");
2583  // dimension 0 of the arrays must match
2584  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch( errmsg, inRefCell, 0, points, 0), std::invalid_argument, errmsg);
2585  }
2586  else if (getrank(points) == 3) {
2587  TEUCHOS_TEST_FOR_EXCEPTION( !(getrank(inRefCell) == 2 ), std::invalid_argument,
2588  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-3 input array requires rank-2 output array.");
2589  // dimensions 0 and 1 of the arrays must match
2590  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch( errmsg, inRefCell, 0,1, points, 0,1), std::invalid_argument, errmsg);
2591  }
2592  else{
2593  TEUCHOS_TEST_FOR_EXCEPTION( !( (getrank(points) == 1) || (getrank(points) == 2) || (getrank(points) == 3) ), std::invalid_argument,
2594  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1, 2 or 3 required for input points array. ");
2595  }
2596 
2597  // The last dimension of points array at (rank - 1) is the spatial dimension. Must equal the cell dimension.
2598  TEUCHOS_TEST_FOR_EXCEPTION( !((size_t)points.dimension(apRank - 1) == (size_t)cellTopo.getDimension() ), std::invalid_argument,
2599  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): Point and cell dimensions do not match. ");
2600 
2601 #endif
2602 
2603  // Initializations
2604  int dim0 = 1;
2605  int dim1 = 1;
2606  int pointDim = 0;
2607  switch(apRank) {
2608  case 1:
2609  pointDim = static_cast<size_t>(points.dimension(0));
2610  break;
2611  case 2:
2612  dim1 = static_cast<size_t>(points.dimension(0));
2613  pointDim = static_cast<size_t>(points.dimension(1));
2614  break;
2615  case 3:
2616  dim0 = static_cast<size_t>(points.dimension(0));
2617  dim1 = static_cast<size_t>(points.dimension(1));
2618  pointDim = static_cast<size_t>(points.dimension(2));
2619  break;
2620  default:
2621  TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= getrank(points) ) && (getrank(points) <= 3) ), std::invalid_argument,
2622  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): rank-1, 2 or 3 required for input points array. ");
2623  }// switch
2624 
2625 
2626  // This method can handle up to rank-3 input arrays. The spatial dim must be the last dimension.
2627  // The method uses [] accessor because array rank is determined at runtime and the appropriate
2628  // (i,j,..,k) accessor is not known. Use of [] requires the following offsets:
2629  // for input array = i0*dim1*pointDim + i1*dim1 (computed in 2 pieces: inPtr0 and inPtr1, resp)
2630  // for output array = i0*dim1 (computed in one piece: outPtr0)
2631  int inPtr0 = 0;
2632  int inPtr1 = 0;
2633  int outPtr0 = 0;
2634  Scalar point[3] = {0.0, 0.0, 0.0};
2635 
2636  for(int i0 = 0; i0 < dim0; i0++){
2637  outPtr0 = i0*dim1;
2638  inPtr0 = outPtr0*pointDim;
2639 
2640  for(int i1 = 0; i1 < dim1; i1++) {
2641  inPtr1 = inPtr0 + i1*pointDim;
2642  point[0] = points[inPtr1];
2643  if(pointDim > 1) {
2644  point[1] = points[inPtr1 + 1];
2645  if(pointDim > 2) {
2646  point[2] = points[inPtr1 + 2];
2647  if(pointDim > 3) {
2648  TEUCHOS_TEST_FOR_EXCEPTION( !( (1 <= pointDim) && (pointDim <= 3)), std::invalid_argument,
2649  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): Input array specifies invalid point dimension ");
2650  }
2651  }
2652  } //if(pointDim > 1)
2653  inRefCell[outPtr0 + i1] = checkPointInclusion(point, pointDim, cellTopo, threshold);
2654  } // for (i1)
2655  } // for(i2)
2656 
2657 }
2658 
2659 
2660 template<class Scalar>
2661 template<class ArrayIncl, class ArrayPoint, class ArrayCell>
2662 void CellTools<Scalar>::checkPointwiseInclusion(ArrayIncl & inCell,
2663  const ArrayPoint & points,
2664  const ArrayCell & cellWorkset,
2665  const shards::CellTopology & cell,
2666  const int & whichCell,
2667  const double & threshold)
2668 {
2669  INTREPID_VALIDATE( validateArguments_checkPointwiseInclusion(inCell, points, cellWorkset, whichCell, cell) );
2670 
2671  // For cell topologies with reference cells this test maps the points back to the reference cell
2672  // and uses the method for reference cells
2673  unsigned baseKey = cell.getBaseCellTopologyData() -> key;
2674 
2675  switch(baseKey){
2676 
2677  case shards::Line<>::key :
2678  case shards::Triangle<>::key:
2679  case shards::Quadrilateral<>::key :
2680  case shards::Tetrahedron<>::key :
2681  case shards::Hexahedron<>::key :
2682  case shards::Wedge<>::key :
2683  case shards::Pyramid<>::key :
2684  {
2685  FieldContainer<Scalar> refPoints;
2686 
2687  if(getrank(points) == 2){
2688  refPoints.resize(static_cast<size_t>(points.dimension(0)), static_cast<size_t>(points.dimension(1)) );
2689  mapToReferenceFrame(refPoints, points, cellWorkset, cell, whichCell);
2690  checkPointwiseInclusion(inCell, refPoints, cell, threshold );
2691  }
2692  else if(getrank(points) == 3){
2693  refPoints.resize(static_cast<size_t>(points.dimension(0)), static_cast<size_t>(points.dimension(1)), static_cast<size_t>(points.dimension(2)) );
2694  mapToReferenceFrame(refPoints, points, cellWorkset, cell, whichCell);
2695  checkPointwiseInclusion(inCell, refPoints, cell, threshold );
2696  }
2697  break;
2698  }
2699  default:
2700  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
2701  ">>> ERROR (Intrepid::CellTools::checkPointwiseInclusion): cell topology not supported");
2702  }// switch
2703 
2704 }
2705 
2706 
2707 //============================================================================================//
2708 // //
2709 // Validation of input/output arguments for CellTools methods //
2710 // //
2711 //============================================================================================//
2712 
2713 template<class Scalar>
2714 template<class ArrayJac, class ArrayPoint, class ArrayCell>
2715 void CellTools<Scalar>::validateArguments_setJacobian(const ArrayJac & jacobian,
2716  const ArrayPoint & points,
2717  const ArrayCell & cellWorkset,
2718  const int & whichCell,
2719  const shards::CellTopology & cellTopo){
2720 
2721  // Validate cellWorkset array
2722  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(cellWorkset) != 3), std::invalid_argument,
2723  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 3 required for cellWorkset array");
2724 
2725  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(0)) <= 0), std::invalid_argument,
2726  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) >= 1 required for cellWorkset array");
2727 
2728  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(1)) != (size_t)cellTopo.getSubcellCount(0) ), std::invalid_argument,
2729  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
2730 
2731  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(2)) != (size_t)cellTopo.getDimension() ), std::invalid_argument,
2732  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of cellWorkset array does not match cell dimension");
2733 
2734  // validate whichCell. It can be either -1 (default value) or a valid cell ordinal.
2735  TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (static_cast<size_t>(whichCell) < static_cast<size_t>(cellWorkset.dimension(0)) ) ) || (whichCell == -1) ), std::invalid_argument,
2736  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): whichCell = -1 or a valid cell ordinal is required.");
2737 
2738 
2739  // Validate points array: can be rank-2 (P,D) or rank-3 (C,P,D)
2740  // If rank-2: admissible jacobians: rank-3 (P,D,D) or rank-4 (C,P,D,D); admissible whichCell: -1 (default) or cell ordinal.
2741  if(getrank(points) == 2) {
2742  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(points.dimension(0)) <= 0), std::invalid_argument,
2743  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of points) >= 1 required for points array ");
2744 
2745  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(points.dimension(1)) != (size_t)cellTopo.getDimension() ), std::invalid_argument,
2746  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (spatial dimension) of points array does not match cell dimension");
2747 
2748  // Validate the output array for the Jacobian: if whichCell == -1 all Jacobians are computed, rank-4 (C,P,D,D) required
2749  if(whichCell == -1) {
2750  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(jacobian) != 4), std::invalid_argument,
2751  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 4 required for jacobian array");
2752 
2753  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(0)) != static_cast<size_t>(cellWorkset.dimension(0))), std::invalid_argument,
2754  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of jacobian array must equal dim 0 of cellWorkset array");
2755 
2756  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(1)) != static_cast<size_t>(points.dimension(0))), std::invalid_argument,
2757  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) of jacobian array must equal dim 0 of points array");
2758 
2759  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(2)) != static_cast<size_t>(points.dimension(1))), std::invalid_argument,
2760  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of jacobian array must equal dim 1 of points array");
2761 
2762  TEUCHOS_TEST_FOR_EXCEPTION( !(static_cast<size_t>(jacobian.dimension(2)) == static_cast<size_t>(jacobian.dimension(3)) ), std::invalid_argument,
2763  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 = dim 3 (same spatial dimensions) required for jacobian array. ");
2764 
2765  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < static_cast<size_t>(jacobian.dimension(3)) ) && (static_cast<size_t>(jacobian.dimension(3)) < 4) ), std::invalid_argument,
2766  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 and dim 3 (spatial dimensions) must be between 1 and 3. ");
2767  }
2768  // A single cell Jacobian is computed when whichCell != -1 (whichCell has been already validated), rank-3 (P,D,D) required
2769  else {
2770  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(jacobian) != 3), std::invalid_argument,
2771  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 3 required for jacobian array");
2772 
2773  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(0)) != static_cast<size_t>(points.dimension(0))), std::invalid_argument,
2774  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of points) of jacobian array must equal dim 0 of points array");
2775 
2776  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(1)) != static_cast<size_t>(points.dimension(1))), std::invalid_argument,
2777  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (spatial dimension) of jacobian array must equal dim 1 of points array");
2778 
2779  TEUCHOS_TEST_FOR_EXCEPTION( !(static_cast<size_t>(jacobian.dimension(1)) == static_cast<size_t>(jacobian.dimension(2)) ), std::invalid_argument,
2780  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 = dim 2 (same spatial dimensions) required for jacobian array. ");
2781 
2782  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < static_cast<size_t>(jacobian.dimension(1)) ) && (static_cast<size_t>(jacobian.dimension(1)) < 4) ), std::invalid_argument,
2783  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 and dim 2 (spatial dimensions) must be between 1 and 3. ");
2784  }
2785  }
2786  // Point array is rank-3 (C,P,D): requires whichCell = -1 and rank-4 (C,P,D,D) jacobians
2787  else if(getrank(points) ==3){
2788  std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian):";
2789  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(points.dimension(0)) != static_cast<size_t>(cellWorkset.dimension(0)) ), std::invalid_argument,
2790  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of points array must equal dim 0 of cellWorkset array");
2791 
2792  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(points.dimension(1)) <= 0), std::invalid_argument,
2793  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) >= 1 required for points array ");
2794 
2795  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(points.dimension(2)) != (size_t)cellTopo.getDimension() ), std::invalid_argument,
2796  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of points array does not match cell dimension");
2797 
2798  TEUCHOS_TEST_FOR_EXCEPTION( (whichCell != -1), std::invalid_argument,
2799  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): default value whichCell=-1 required for rank-3 input points");
2800 
2801  // rank-4 (C,P,D,D) jacobian required for rank-3 (C,P,D) input points
2802  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, jacobian, 4, 4), std::invalid_argument,errmsg);
2803 
2804  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(0)) != static_cast<size_t>(points.dimension(0))), std::invalid_argument,
2805  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 0 (number of cells) of jacobian array must equal dim 0 of points array");
2806 
2807  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(1)) != static_cast<size_t>(points.dimension(1))), std::invalid_argument,
2808  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 1 (number of points) of jacobian array must equal dim 1 of points array");
2809 
2810  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(jacobian.dimension(2)) != static_cast<size_t>(points.dimension(2))), std::invalid_argument,
2811  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 (spatial dimension) of jacobian array must equal dim 2 of points array");
2812 
2813  TEUCHOS_TEST_FOR_EXCEPTION( !(static_cast<size_t>(jacobian.dimension(2)) == static_cast<size_t>(jacobian.dimension(3)) ), std::invalid_argument,
2814  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 = dim 3 (same spatial dimensions) required for jacobian array. ");
2815 
2816  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < static_cast<size_t>(jacobian.dimension(3)) ) && (static_cast<size_t>(jacobian.dimension(3)) < 4) ), std::invalid_argument,
2817  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): dim 2 and dim 3 (spatial dimensions) must be between 1 and 3. ");
2818  }
2819  else {
2820  TEUCHOS_TEST_FOR_EXCEPTION( !( (getrank(points) == 2) && (getrank(points) ==3) ), std::invalid_argument,
2821  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobian): rank = 2 or 3 required for points array");
2822  }
2823 }
2824 
2825 
2826 
2827 template<class Scalar>
2828 template<class ArrayJacInv, class ArrayJac>
2829 void CellTools<Scalar>::validateArguments_setJacobianInv(const ArrayJacInv & jacobianInv,
2830  const ArrayJac & jacobian)
2831 {
2832  // Validate input jacobian array: admissible ranks & dimensions are:
2833  // - rank-4 with dimensions (C,P,D,D), or rank-3 with dimensions (P,D,D).
2834  int jacobRank = getrank(jacobian);
2835  TEUCHOS_TEST_FOR_EXCEPTION( !( (jacobRank == 4) || (jacobRank == 3) ), std::invalid_argument,
2836  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): rank = 4 or 3 required for jacobian array. ");
2837 
2838  // Verify correctness of spatial dimensions - they are the last two dimensions of the array: rank-2 and rank-1
2839  TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(jacobRank - 1) == jacobian.dimension(jacobRank - 2) ), std::invalid_argument,
2840  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-2) = dim(rank-2) (same spatial dimensions) required for jacobian array. ");
2841 
2842  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(jacobRank - 1) ) && (jacobian.dimension(jacobRank - 1) < 4) ), std::invalid_argument,
2843  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-1) and dim(rank-2) (spatial dimensions) must be between 1 and 3. ");
2844 
2845  // Validate output jacobianInv array: must have the same rank and dimensions as the input array.
2846  std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv):";
2847 
2848  TEUCHOS_TEST_FOR_EXCEPTION( !(requireRankMatch(errmsg, jacobianInv, jacobian) ), std::invalid_argument, errmsg);
2849 
2850  TEUCHOS_TEST_FOR_EXCEPTION( !(requireDimensionMatch(errmsg, jacobianInv, jacobian) ), std::invalid_argument, errmsg);
2851 }
2852 
2853 
2854 
2855 template<class Scalar>
2856 template<class ArrayJacDet, class ArrayJac>
2857 void CellTools<Scalar>::validateArguments_setJacobianDetArgs(const ArrayJacDet & jacobianDet,
2858  const ArrayJac & jacobian)
2859 {
2860  // Validate input jacobian array: admissible ranks & dimensions are:
2861  // - rank-4 with dimensions (C,P,D,D), or rank-3 with dimensions (P,D,D).
2862  int jacobRank = getrank(jacobian);
2863  TEUCHOS_TEST_FOR_EXCEPTION( !( (jacobRank == 4) || (jacobRank == 3) ), std::invalid_argument,
2864  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): rank = 4 or 3 required for jacobian array. ");
2865 
2866  // Verify correctness of spatial dimensions - they are the last two dimensions of the array: rank-2 and rank-1
2867  TEUCHOS_TEST_FOR_EXCEPTION( !(jacobian.dimension(jacobRank - 1) == jacobian.dimension(jacobRank - 2) ), std::invalid_argument,
2868  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-2) = dim(rank-2) (same spatial dimensions) required for jacobian array. ");
2869 
2870  TEUCHOS_TEST_FOR_EXCEPTION( !( (0 < jacobian.dimension(jacobRank - 1) ) && (jacobian.dimension(jacobRank - 1) < 4) ), std::invalid_argument,
2871  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianInv): dim(rank-1) and dim(rank-2) (spatial dimensions) must be between 1 and 3. ");
2872 
2873 
2874  // Validate output jacobianDet array: must be rank-2 with dimensions (C,P) if jacobian was rank-4:
2875  if(jacobRank == 4){
2876  TEUCHOS_TEST_FOR_EXCEPTION( !(getrank(jacobianDet) == 2), std::invalid_argument,
2877  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): rank = 2 required for jacobianDet if jacobian is rank-4. ");
2878 
2879  TEUCHOS_TEST_FOR_EXCEPTION( !(static_cast<size_t>(jacobianDet.dimension(0)) == static_cast<size_t>(jacobian.dimension(0)) ), std::invalid_argument,
2880  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 0 (number of cells) of jacobianDet array must equal dim 0 of jacobian array. ");
2881 
2882  TEUCHOS_TEST_FOR_EXCEPTION( !(static_cast<size_t>(jacobianDet.dimension(1)) == static_cast<size_t>(jacobian.dimension(1)) ), std::invalid_argument,
2883  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 1 (number of points) of jacobianDet array must equal dim 1 of jacobian array.");
2884  }
2885 
2886  // must be rank-1 with dimension (P) if jacobian was rank-3
2887  else {
2888  TEUCHOS_TEST_FOR_EXCEPTION( !(getrank(jacobianDet) == 1), std::invalid_argument,
2889  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): rank = 1 required for jacobianDet if jacobian is rank-3. ");
2890 
2891  TEUCHOS_TEST_FOR_EXCEPTION( !(static_cast<size_t>(jacobianDet.dimension(0)) == static_cast<size_t>(jacobian.dimension(0)) ), std::invalid_argument,
2892  ">>> ERROR (Intrepid::CellTools::validateArguments_setJacobianDetArgs): dim 0 (number of points) of jacobianDet array must equal dim 0 of jacobian array.");
2893  }
2894 }
2895 
2896 
2897 
2898 template<class Scalar>
2899 template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
2900 void CellTools<Scalar>::validateArguments_mapToPhysicalFrame(const ArrayPhysPoint & physPoints,
2901  const ArrayRefPoint & refPoints,
2902  const ArrayCell & cellWorkset,
2903  const shards::CellTopology & cellTopo,
2904  const int& whichCell)
2905 {
2906  std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame):";
2907 
2908  // Validate cellWorkset array
2909  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(cellWorkset) != 3), std::invalid_argument,
2910  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 3 required for cellWorkset array");
2911 
2912  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(0)) <= 0), std::invalid_argument,
2913  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) >= 1 required for cellWorkset array");
2914 
2915  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(1)) != (size_t)cellTopo.getSubcellCount(0) ), std::invalid_argument,
2916  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
2917 
2918  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(2)) != (size_t)cellTopo.getDimension() ), std::invalid_argument,
2919  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) of cellWorkset array does not match cell dimension");
2920 
2921 
2922 
2923 
2924 TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && ((size_t)whichCell < (size_t)cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
2925  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): whichCell = -1 or a valid cell ordinal is required.");
2926 
2927  // Validate refPoints array: can be rank-2 (P,D) or rank-3 (C,P,D) array
2928  // If rank-2: admissible output array is (P,D) or (C,P,D); admissible whichCell: -1 (default) or cell ordinal
2929  if(getrank(refPoints) == 2) {
2930  TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(0) <= 0), std::invalid_argument,
2931  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of points) >= 1 required for refPoints array ");
2932 
2933  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)refPoints.dimension(1) != (size_t)cellTopo.getDimension() ), std::invalid_argument,
2934  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (spatial dimension) of refPoints array does not match cell dimension");
2935 
2936  // Validate output array: whichCell = -1 requires rank-3 array with dimensions (C,P,D)
2937  if(whichCell == -1) {
2938  TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(physPoints) != 3) && (whichCell == -1) ), std::invalid_argument,
2939  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 3 required for physPoints array for the default whichCell value");
2940 
2941  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)physPoints.dimension(0) != (size_t)cellWorkset.dimension(0)), std::invalid_argument,
2942  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) of physPoints array must equal dim 0 of cellWorkset array");
2943 
2944  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)physPoints.dimension(1) != (size_t)refPoints.dimension(0)), std::invalid_argument,
2945  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of points) of physPoints array must equal dim 0 of refPoints array");
2946 
2947  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)physPoints.dimension(2) != (size_t)cellTopo.getDimension()), std::invalid_argument,
2948  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) does not match cell dimension ");
2949  }
2950  // 0 <= whichCell < num cells requires rank-2 (P,D) arrays for both refPoints and physPoints
2951  else{
2952  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(physPoints) != 2), std::invalid_argument,
2953  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 2 required for physPoints array");
2954 
2955  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)physPoints.dimension(0) != (size_t)refPoints.dimension(0)), std::invalid_argument,
2956  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of points) of physPoints array must equal dim 0 of refPoints array");
2957 
2958  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)physPoints.dimension(1) != (size_t)cellTopo.getDimension()), std::invalid_argument,
2959  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (spatial dimension) does not match cell dimension ");
2960  }
2961  }
2962  // refPoints is (C,P,D): requires physPoints to be (C,P,D) and whichCell=-1 (because all cell mappings are applied)
2963  else if(getrank(refPoints) == 3) {
2964 
2965  // 1. validate refPoints dimensions and rank
2966  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)refPoints.dimension(0) !=(size_t) cellWorkset.dimension(0) ), std::invalid_argument,
2967  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 0 (number of cells) of refPoints and cellWorkset arraya are required to match ");
2968 
2969  TEUCHOS_TEST_FOR_EXCEPTION( (refPoints.dimension(1) <= 0), std::invalid_argument,
2970  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 1 (number of points) >= 1 required for refPoints array ");
2971 
2972  TEUCHOS_TEST_FOR_EXCEPTION( ((size_t)refPoints.dimension(2) != (size_t)cellTopo.getDimension() ), std::invalid_argument,
2973  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): dim 2 (spatial dimension) of refPoints array does not match cell dimension");
2974 
2975  // 2. whichCell must be -1
2976  TEUCHOS_TEST_FOR_EXCEPTION( (whichCell != -1), std::invalid_argument,
2977  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): default value is required for rank-3 refPoints array");
2978 
2979  // 3. physPoints must match rank and dimensions of refPoints
2980  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankMatch(errmsg, refPoints, physPoints), std::invalid_argument, errmsg );
2981  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, refPoints, physPoints), std::invalid_argument, errmsg);
2982  }
2983  // if rank is not 2 or 3 throw exception
2984  else {
2985  TEUCHOS_TEST_FOR_EXCEPTION( !( (getrank(refPoints) == 2) || (getrank(refPoints) == 3) ), std::invalid_argument,
2986  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToPhysicalFrame): rank = 2 or 3 required for refPoints array");
2987  }
2988 }
2989 template<class Scalar>
2990 template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
2991 void CellTools<Scalar>::validateArguments_mapToReferenceFrame(const ArrayRefPoint & refPoints,
2992  const ArrayPhysPoint & physPoints,
2993  const ArrayCell & cellWorkset,
2994  const shards::CellTopology & cellTopo,
2995  const int& whichCell)
2996 {
2997  std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
2998  std::string errmsg1 = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
2999 
3000  // Validate cellWorkset array
3001  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(cellWorkset) != 3), std::invalid_argument,
3002  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): rank = 3 required for cellWorkset array");
3003 
3004  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(0)) <= 0), std::invalid_argument,
3005  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 0 (number of cells) >= 1 required for cellWorkset array");
3006 
3007  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(1)) != (size_t)cellTopo.getSubcellCount(0) ), std::invalid_argument,
3008  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
3009 
3010  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(2)) != (size_t)cellTopo.getDimension() ), std::invalid_argument,
3011  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): dim 2 (spatial dimension) of cellWorkset array does not match cell dimension");
3012 
3013  // Validate whichCell. It can be either -1 (default value) or a valid celli ordinal.
3014  TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && ((size_t)whichCell <(size_t) cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
3015  ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame): whichCell = -1 or a valid cell ordinal is required.");
3016  // Admissible ranks and dimensions of refPoints and physPoints depend on whichCell value:
3017  // default is to map multiple sets of points to multiple sets of points. (C,P,D) arrays required
3018  int validRank;
3019  if(whichCell == -1) {
3020  validRank = 3;
3021  errmsg1 += " default value of whichCell requires rank-3 arrays:";
3022  }
3023  // whichCell is valid cell ordinal => we map single set of pts to a single set of pts. (P,D) arrays required
3024  else{
3025  errmsg1 += " rank-2 arrays required when whichCell is valid cell ordinal";
3026  validRank = 2;
3027  }
3028  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg1, refPoints, validRank,validRank), std::invalid_argument, errmsg1);
3029  TEUCHOS_TEST_FOR_EXCEPTION( !requireRankMatch(errmsg1, physPoints, refPoints), std::invalid_argument, errmsg1);
3030  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg1, refPoints, physPoints), std::invalid_argument, errmsg1);
3031 }
3032 
3033 
3034 
3035 template<class Scalar>
3036 template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
3037 void CellTools<Scalar>::validateArguments_mapToReferenceFrame(const ArrayRefPoint & refPoints,
3038  const ArrayInitGuess & initGuess,
3039  const ArrayPhysPoint & physPoints,
3040  const ArrayCell & cellWorkset,
3041  const shards::CellTopology & cellTopo,
3042  const int& whichCell)
3043 {
3044  // Call the method that validates arguments with the default initial guess selection
3045  validateArguments_mapToReferenceFrame(refPoints, physPoints, cellWorkset, cellTopo, whichCell);
3046 
3047  // Then check initGuess: its rank and dimensions must match those of physPoints.
3048  std::string errmsg = ">>> ERROR (Intrepid::CellTools::validateArguments_mapToReferenceFrame):";
3049  TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, initGuess, physPoints), std::invalid_argument, errmsg);
3050 }
3051 
3052 
3053 template<class Scalar>
3054 template<class ArrayIncl, class ArrayPoint, class ArrayCell>
3055 void CellTools<Scalar>::validateArguments_checkPointwiseInclusion(ArrayIncl & inCell,
3056  const ArrayPoint & physPoints,
3057  const ArrayCell & cellWorkset,
3058  const int & whichCell,
3059  const shards::CellTopology & cell)
3060 {
3061  // Validate cellWorkset array
3062  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(cellWorkset) != 3), std::invalid_argument,
3063  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 3 required for cellWorkset array");
3064 
3065  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(0)) <= 0), std::invalid_argument,
3066  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells) >= 1 required for cellWorkset array");
3067 
3068  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(1)) != (size_t)cell.getSubcellCount(0) ), std::invalid_argument,
3069  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of cell nodes) of cellWorkset array does not match cell topology");
3070 
3071  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(cellWorkset.dimension(2)) != (size_t)cell.getDimension() ), std::invalid_argument,
3072  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 2 (spatial dimension) of cellWorkset array does not match cell dimension");
3073 
3074 
3075  // Validate whichCell It can be either -1 (default value) or a valid cell ordinal.
3076  TEUCHOS_TEST_FOR_EXCEPTION( !( ( (0 <= whichCell ) && (whichCell < cellWorkset.dimension(0) ) ) || (whichCell == -1) ), std::invalid_argument,
3077  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = -1 or a valid cell ordinal is required.");
3078 
3079  // Validate points array: can be rank-2 (P,D) or rank-3 (C,P,D)
3080  // If rank-2: admissible inCell is rank-1 (P); admissible whichCell is valid cell ordinal but not -1.
3081  if(getrank(physPoints) == 2) {
3082 
3083  TEUCHOS_TEST_FOR_EXCEPTION( (whichCell == -1), std::invalid_argument,
3084  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = a valid cell ordinal is required with rank-2 input array.");
3085 
3086  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(physPoints.dimension(0)) <= 0), std::invalid_argument,
3087  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of points) >= 1 required for physPoints array ");
3088 
3089  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(physPoints.dimension(1)) != (size_t)cell.getDimension() ), std::invalid_argument,
3090  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (spatial dimension) of physPoints array does not match cell dimension");
3091 
3092  // Validate inCell
3093  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inCell) != 1), std::invalid_argument,
3094  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 1 required for inCell array");
3095 
3096  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(inCell.dimension(0)) != static_cast<size_t>(physPoints.dimension(0))), std::invalid_argument,
3097  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of points) of inCell array must equal dim 0 of physPoints array");
3098  }
3099  // If rank-3: admissible inCell is rank-2 (C,P); admissible whichCell = -1.
3100  else if (getrank(physPoints) == 3){
3101 
3102  TEUCHOS_TEST_FOR_EXCEPTION( !(whichCell == -1), std::invalid_argument,
3103  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): whichCell = -1 is required with rank-3 input array.");
3104 
3105  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(physPoints.dimension(0)) != static_cast<size_t>(cellWorkset.dimension(0)) ), std::invalid_argument,
3106  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells) of physPoints array must equal dim 0 of cellWorkset array ");
3107 
3108  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(physPoints.dimension(1)) <= 0), std::invalid_argument,
3109  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of points) >= 1 required for physPoints array ");
3110 
3111  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(physPoints.dimension(2)) != (size_t)cell.getDimension() ), std::invalid_argument,
3112  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 2 (spatial dimension) of physPoints array does not match cell dimension");
3113 
3114  // Validate inCell
3115  TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inCell) != 2), std::invalid_argument,
3116  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 2 required for inCell array");
3117 
3118  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(inCell.dimension(0)) != static_cast<size_t>(physPoints.dimension(0))), std::invalid_argument,
3119  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 0 (number of cells) of inCell array must equal dim 0 of physPoints array");
3120 
3121  TEUCHOS_TEST_FOR_EXCEPTION( (static_cast<size_t>(inCell.dimension(1)) != static_cast<size_t>(physPoints.dimension(1))), std::invalid_argument,
3122  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): dim 1 (number of points) of inCell array must equal dim 1 of physPoints array");
3123  }
3124  else {
3125  TEUCHOS_TEST_FOR_EXCEPTION( !( (getrank(physPoints) == 2) && (getrank(physPoints) ==3) ), std::invalid_argument,
3126  ">>> ERROR (Intrepid::CellTools::validateArguments_checkPointwiseInclusion): rank = 2 or 3 required for points array");
3127  }
3128 }
3129 
3130 
3131 
3132 //============================================================================================//
3133 // //
3134 // Debug //
3135 // //
3136 //============================================================================================//
3137 
3138 
3139 template<class Scalar>
3140 void CellTools<Scalar>::printSubcellVertices(const int subcellDim,
3141  const int subcellOrd,
3142  const shards::CellTopology & parentCell){
3143 
3144  // Get number of vertices for the specified subcell and parent cell dimension
3145  int subcVertexCount = parentCell.getVertexCount(subcellDim, subcellOrd);
3146  int cellDim = parentCell.getDimension();
3147 
3148  // Allocate space for the subcell vertex coordinates
3149  FieldContainer<double> subcellVertices(subcVertexCount, cellDim);
3150 
3151  // Retrieve the vertex coordinates
3152  getReferenceSubcellVertices(subcellVertices,
3153  subcellDim,
3154  subcellOrd,
3155  parentCell);
3156 
3157  // Print the vertices
3158  std::cout
3159  << " Subcell " << std::setw(2) << subcellOrd
3160  << " is " << parentCell.getName(subcellDim, subcellOrd) << " with vertices = {";
3161 
3162  // Loop over subcell vertices
3163  for(int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){
3164  std::cout<< "(";
3165 
3166  // Loop over vertex Cartesian coordinates
3167  for(int dim = 0; dim < (int)parentCell.getDimension(); dim++){
3168  std::cout << subcellVertices(subcVertOrd, dim);
3169  if(dim < (int)parentCell.getDimension()-1 ) { std::cout << ","; }
3170  }
3171  std::cout<< ")";
3172  if(subcVertOrd < subcVertexCount - 1) { std::cout << ", "; }
3173  }
3174  std::cout << "}\n";
3175 }
3176 
3177 
3178 template<class Scalar>
3179 template<class ArrayCell>
3180 void CellTools<Scalar>::printWorksetSubcell(const ArrayCell & cellWorkset,
3181  const shards::CellTopology & parentCell,
3182  const int& pCellOrd,
3183  const int& subcellDim,
3184  const int& subcellOrd,
3185  const int& fieldWidth){
3186 
3187  // Get the ordinals, relative to reference cell, of subcell cellWorkset
3188  int subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
3189  int pCellDim = parentCell.getDimension();
3190  std::vector<int> subcNodeOrdinals(subcNodeCount);
3191 
3192  for(int i = 0; i < subcNodeCount; i++){
3193  subcNodeOrdinals[i] = parentCell.getNodeMap(subcellDim, subcellOrd, i);
3194  }
3195 
3196  // Loop over parent cells and print subcell cellWorkset
3197 
3198  std::cout
3199  << " Subcell " << subcellOrd << " on parent cell " << pCellOrd << " is "
3200  << parentCell.getName(subcellDim, subcellOrd) << " with node(s) \n ({";
3201 
3202  for(int i = 0; i < subcNodeCount; i++){
3203 
3204  // print Cartesian coordinates of the node
3205  for(int dim = 0; dim < pCellDim; dim++){
3206  std::cout
3207  << std::setw(fieldWidth) << std::right << cellWorkset(pCellOrd, subcNodeOrdinals[i], dim);
3208  if(dim < pCellDim - 1){ std::cout << ","; }
3209  }
3210  std::cout << "}";
3211  if(i < subcNodeCount - 1){ std::cout <<", {"; }
3212  }
3213  std::cout << ")\n\n";
3214 }
3215 //============================================================================================//
3216 // //
3217 // Control Volume Coordinates //
3218 // //
3219 //============================================================================================//
3220 
3221  template<class Scalar>
3222  template<class ArrayCVCoord, class ArrayCellCoord>
3223  void CellTools<Scalar>::getSubCVCoords(ArrayCVCoord & subCVCoords,
3224  const ArrayCellCoord & cellCoords,
3225  const shards::CellTopology& primaryCell)
3226  {
3227 
3228  // get array dimensions
3229  int numCells = cellCoords.dimension(0);
3230  int numNodesPerCell = cellCoords.dimension(1);
3231  int spaceDim = cellCoords.dimension(2);
3232 
3233  // num edges per primary cell
3234  int numEdgesPerCell = primaryCell.getEdgeCount();
3235 
3236  // num faces per primary cell
3237  int numFacesPerCell = 0;
3238  if (spaceDim > 2){
3239  numFacesPerCell = primaryCell.getFaceCount();
3240  }
3241 
3242  // get cell centroids
3243  Intrepid::FieldContainer<Scalar> barycenter(numCells,spaceDim);
3244  getBarycenter(barycenter,cellCoords);
3245 
3246  // loop over cells
3247  for (int icell = 0; icell < numCells; icell++){
3248 
3249  // get primary edge midpoints
3250  Intrepid::FieldContainer<Scalar> edgeMidpts(numEdgesPerCell,spaceDim);
3251  for (int iedge = 0; iedge < numEdgesPerCell; iedge++){
3252  for (int idim = 0; idim < spaceDim; idim++){
3253 
3254  int node0 = primaryCell.getNodeMap(1,iedge,0);
3255  int node1 = primaryCell.getNodeMap(1,iedge,1);
3256  edgeMidpts(iedge,idim) = (cellCoords(icell,node0,idim) +
3257  cellCoords(icell,node1,idim))/2.0;
3258 
3259  } // end loop over dimensions
3260  } // end loop over cell edges
3261 
3262  // get primary face midpoints in 3-D
3263  int numNodesPerFace;
3264  Intrepid::FieldContainer<Scalar> faceMidpts(numFacesPerCell,spaceDim);
3265  if (spaceDim > 2) {
3266  for (int iface = 0; iface < numFacesPerCell; iface++){
3267  numNodesPerFace = primaryCell.getNodeCount(2,iface);
3268 
3269  for (int idim = 0; idim < spaceDim; idim++){
3270 
3271  for (int inode0 = 0; inode0 < numNodesPerFace; inode0++) {
3272  int node1 = primaryCell.getNodeMap(2,iface,inode0);
3273  faceMidpts(iface,idim) += cellCoords(icell,node1,idim)/numNodesPerFace;
3274  }
3275 
3276  } // end loop over dimensions
3277  } // end loop over cell faces
3278  }
3279 
3280  // define coordinates for subcontrol volumes
3281  switch(primaryCell.getKey() ) {
3282 
3283  // 2-d parent cells
3284  case shards::Triangle<3>::key:
3285  case shards::Quadrilateral<4>::key:
3286 
3287  for (int inode = 0; inode < numNodesPerCell; inode++){
3288  for (int idim = 0; idim < spaceDim; idim++){
3289 
3290  // set first node to primary cell node
3291  subCVCoords(icell,inode,0,idim) = cellCoords(icell,inode,idim);
3292 
3293  // set second node to adjacent edge midpoint
3294  subCVCoords(icell,inode,1,idim) = edgeMidpts(inode,idim);
3295 
3296  // set third node to cell barycenter
3297  subCVCoords(icell,inode,2,idim) = barycenter(icell,idim);
3298 
3299  // set fourth node to other adjacent edge midpoint
3300  int jnode = numNodesPerCell-1;
3301  if (inode > 0) jnode = inode - 1;
3302  subCVCoords(icell,inode,3,idim) = edgeMidpts(jnode,idim);
3303 
3304  } // dim loop
3305  } // node loop
3306 
3307  break;
3308 
3309  case shards::Hexahedron<8>::key:
3310 
3311  for (int idim = 0; idim < spaceDim; idim++){
3312 
3313  // loop over the horizontal quads that define the subcontrol volume coords
3314  for (int icount = 0; icount < 4; icount++){
3315 
3316  // set first node of bottom hex to primary cell node
3317  // and fifth node of upper hex
3318  subCVCoords(icell,icount,0,idim) = cellCoords(icell,icount,idim);
3319  subCVCoords(icell,icount+4,4,idim) = cellCoords(icell,icount+4,idim);
3320 
3321  // set second node of bottom hex to adjacent edge midpoint
3322  // and sixth node of upper hex
3323  subCVCoords(icell,icount,1,idim) = edgeMidpts(icount,idim);
3324  subCVCoords(icell,icount+4,5,idim) = edgeMidpts(icount+4,idim);
3325 
3326  // set third node of bottom hex to bottom face midpoint (number 4)
3327  // and seventh node of upper hex to top face midpoint
3328  subCVCoords(icell,icount,2,idim) = faceMidpts(4,idim);
3329  subCVCoords(icell,icount+4,6,idim) = faceMidpts(5,idim);
3330 
3331  // set fourth node of bottom hex to other adjacent edge midpoint
3332  // and eight node of upper hex to other adjacent edge midpoint
3333  int jcount = 3;
3334  if (icount > 0) jcount = icount - 1;
3335  subCVCoords(icell,icount,3,idim) = edgeMidpts(jcount,idim);
3336  subCVCoords(icell,icount+4,7,idim) = edgeMidpts(jcount+4,idim);
3337 
3338  // set fifth node to vertical edge
3339  // same as first node of upper hex
3340  subCVCoords(icell,icount,4,idim) = edgeMidpts(icount+numNodesPerCell,idim);
3341  subCVCoords(icell,icount+4,0,idim) = edgeMidpts(icount+numNodesPerCell,idim);
3342 
3343  // set sixth node to adjacent face midpoint
3344  // same as second node of upper hex
3345  subCVCoords(icell,icount,5,idim) = faceMidpts(icount,idim);
3346  subCVCoords(icell,icount+4,1,idim) = faceMidpts(icount,idim);
3347 
3348  // set seventh node to barycenter
3349  // same as third node of upper hex
3350  subCVCoords(icell,icount,6,idim) = barycenter(icell,idim);
3351  subCVCoords(icell,icount+4,2,idim) = barycenter(icell,idim);
3352 
3353  // set eighth node to other adjacent face midpoint
3354  // same as fourth node of upper hex
3355  jcount = 3;
3356  if (icount > 0) jcount = icount - 1;
3357  subCVCoords(icell,icount,7,idim) = faceMidpts(jcount,idim);
3358  subCVCoords(icell,icount+4,3,idim) = faceMidpts(jcount,idim);
3359 
3360  } // count loop
3361 
3362  } // dim loop
3363 
3364  break;
3365 
3366  case shards::Tetrahedron<4>::key:
3367 
3368  for (int idim = 0; idim < spaceDim; idim++){
3369 
3370  // loop over the three bottom nodes
3371  for (int icount = 0; icount < 3; icount++){
3372 
3373  // set first node of bottom hex to primary cell node
3374  subCVCoords(icell,icount,0,idim) = cellCoords(icell,icount,idim);
3375 
3376  // set second node of bottom hex to adjacent edge midpoint
3377  subCVCoords(icell,icount,1,idim) = edgeMidpts(icount,idim);
3378 
3379  // set third node of bottom hex to bottom face midpoint (number 3)
3380  subCVCoords(icell,icount,2,idim) = faceMidpts(3,idim);
3381 
3382  // set fourth node of bottom hex to other adjacent edge midpoint
3383  int jcount = 2;
3384  if (icount > 0) jcount = icount - 1;
3385  subCVCoords(icell,icount,3,idim) = edgeMidpts(jcount,idim);
3386 
3387  // set fifth node to vertical edge
3388  subCVCoords(icell,icount,4,idim) = edgeMidpts(icount+3,idim);
3389 
3390  // set sixth node to adjacent face midpoint
3391  subCVCoords(icell,icount,5,idim) = faceMidpts(icount,idim);
3392 
3393  // set seventh node to barycenter
3394  subCVCoords(icell,icount,6,idim) = barycenter(icell,idim);
3395 
3396  // set eighth node to other adjacent face midpoint
3397  jcount = 2;
3398  if (icount > 0) jcount = icount - 1;
3399  subCVCoords(icell,icount,7,idim) = faceMidpts(jcount,idim);
3400 
3401  } //count loop
3402 
3403  // Control volume attached to fourth node
3404  // set first node of bottom hex to primary cell node
3405  subCVCoords(icell,3,0,idim) = cellCoords(icell,3,idim);
3406 
3407  // set second node of bottom hex to adjacent edge midpoint
3408  subCVCoords(icell,3,1,idim) = edgeMidpts(3,idim);
3409 
3410  // set third node of bottom hex to bottom face midpoint (number 3)
3411  subCVCoords(icell,3,2,idim) = faceMidpts(2,idim);
3412 
3413  // set fourth node of bottom hex to other adjacent edge midpoint
3414  subCVCoords(icell,3,3,idim) = edgeMidpts(5,idim);
3415 
3416  // set fifth node to vertical edge
3417  subCVCoords(icell,3,4,idim) = edgeMidpts(4,idim);
3418 
3419  // set sixth node to adjacent face midpoint
3420  subCVCoords(icell,3,5,idim) = faceMidpts(0,idim);
3421 
3422  // set seventh node to barycenter
3423  subCVCoords(icell,3,6,idim) = barycenter(icell,idim);
3424 
3425  // set eighth node to other adjacent face midpoint
3426  subCVCoords(icell,3,7,idim) = faceMidpts(1,idim);
3427 
3428  } // dim loop
3429 
3430  break;
3431 
3432  default:
3433  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
3434  ">>> ERROR (getSubCVCoords: invalid cell topology.");
3435  } // cell key
3436 
3437  } // cell loop
3438 
3439 } // getSubCVCoords
3440 
3441  template<class Scalar>
3442  template<class ArrayCent, class ArrayCellCoord>
3443  void CellTools<Scalar>::getBarycenter(ArrayCent & barycenter, const ArrayCellCoord & cellCoords)
3444 {
3445  // get array dimensions
3446  int numCells = cellCoords.dimension(0);
3447  int numVertsPerCell = cellCoords.dimension(1);
3448  int spaceDim = cellCoords.dimension(2);
3449 
3450  if (spaceDim == 2)
3451  {
3452  // Method for general polygons
3453  for (int icell = 0; icell < numCells; icell++){
3454 
3455  Intrepid::FieldContainer<Scalar> cell_centroid(spaceDim);
3456  Scalar area = 0;
3457 
3458  for (int inode = 0; inode < numVertsPerCell; inode++){
3459 
3460  int jnode = inode + 1;
3461  if (jnode >= numVertsPerCell) {
3462  jnode = 0;
3463  }
3464 
3465  Scalar area_mult = cellCoords(icell,inode,0)*cellCoords(icell,jnode,1)
3466  - cellCoords(icell,jnode,0)*cellCoords(icell,inode,1);
3467  cell_centroid(0) += (cellCoords(icell,inode,0) + cellCoords(icell,jnode,0))*area_mult;
3468  cell_centroid(1) += (cellCoords(icell,inode,1) + cellCoords(icell,jnode,1))*area_mult;
3469 
3470  area += 0.5*area_mult;
3471  }
3472 
3473  barycenter(icell,0) = cell_centroid(0)/(6.0*area);
3474  barycenter(icell,1) = cell_centroid(1)/(6.0*area);
3475  }
3476 
3477  }
3478  else
3479  {
3480  // This method works fine for simplices, but for other 3-d shapes
3481  // is not precisely accurate. Could replace with approximate integration
3482  // perhaps.
3483  for (int icell = 0; icell < numCells; icell++){
3484 
3485  Intrepid::FieldContainer<Scalar> cell_centroid(spaceDim);
3486 
3487  for (int inode = 0; inode < numVertsPerCell; inode++){
3488  for (int idim = 0; idim < spaceDim; idim++){
3489  cell_centroid(idim) += cellCoords(icell,inode,idim)/numVertsPerCell;
3490  }
3491  }
3492  for (int idim = 0; idim < spaceDim; idim++){
3493  barycenter(icell,idim) = cell_centroid(idim);
3494  }
3495  }
3496  }
3497 
3498  } // get Barycenter
3499 } // namespace Intrepid
3500 #endif
Intrepid::FieldContainer::size
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
Definition: Intrepid_FieldContainerDef.hpp:601
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:65
Intrepid::CellTools::getReferenceNode
static const double * getReferenceNode(const shards::CellTopology &cell, const int nodeOrd)
Retrieves the Cartesian coordinates of a reference cell node.
Definition: Intrepid_CellToolsDef.hpp:529
Intrepid::CellTools::validateArguments_mapToReferenceFrame
static void validateArguments_mapToReferenceFrame(const ArrayRefPoint &refPoints, const ArrayPhysPoint &physPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell)
Validates arguments to Intrepid::CellTools::mapToReferenceFrame with default initial guess...
Definition: Intrepid_CellToolsDef.hpp:2991
Intrepid::Basis_HGRAD_HEX_I2_FEM
Implementation of the serendipity-family H(grad)-compatible FEM basis of degree 2 on a Hexahedron cel...
Definition: Intrepid_HGRAD_HEX_I2_FEM.hpp:123
Intrepid::CellTools::getReferenceEdgeTangent
static void getReferenceEdgeTangent(ArrayEdgeTangent &refEdgeTangent, const int &edgeOrd, const shards::CellTopology &parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
Definition: Intrepid_CellToolsDef.hpp:2094
Intrepid::CellTools::mapToReferenceSubcell
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
Definition: Intrepid_CellToolsDef.hpp:2027
Intrepid::Basis_HGRAD_WEDGE_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Wedge cell.
Definition: Intrepid_HGRAD_WEDGE_C1_FEM.hpp:85
Intrepid::CellTools::getPhysicalFaceNormals
static void getPhysicalFaceNormals(ArrayFaceNormal &faceNormals, const ArrayJac &worksetJacobians, const int &worksetFaceOrd, const shards::CellTopology &parentCell)
Computes non-normalized normal vectors to physical faces in a face workset ; (see Subcell worksets fo...
Definition: Intrepid_CellToolsDef.hpp:2388
Intrepid::CellTools::setJacobianInv
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F...
Definition: Intrepid_CellToolsDef.hpp:1283
Intrepid::Basis_HGRAD_WEDGE_C2_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Wedge cell.
Definition: Intrepid_HGRAD_WEDGE_C2_FEM.hpp:119
Intrepid::CellTools::getPhysicalEdgeTangents
static void getPhysicalEdgeTangents(ArrayEdgeTangent &edgeTangents, const ArrayJac &worksetJacobians, const int &worksetEdgeOrd, const shards::CellTopology &parentCell)
Computes non-normalized tangent vectors to physical edges in an edge workset ; (see Subcell worksets ...
Definition: Intrepid_CellToolsDef.hpp:2239
Intrepid::CellTools::mapToReferenceFrame
static void mapToReferenceFrame(ArrayRefPoint &refPoints, const ArrayPhysPoint &physPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes , the inverse of the reference-to-physical frame map using a default initial guess...
Definition: Intrepid_CellToolsDef.hpp:1662
INTREPID_MAX_NEWTON
#define INTREPID_MAX_NEWTON
Maximum number of Newton iterations used internally in methods such as computing the action of the in...
Definition: Intrepid_Types.hpp:103
Intrepid::CellTools::getReferenceSideNormal
static void getReferenceSideNormal(ArraySideNormal &refSideNormal, const int &sideOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
Definition: Intrepid_CellToolsDef.hpp:2175
Intrepid::CellTools::getReferenceFaceTangents
static void getReferenceFaceTangents(ArrayFaceTangentU &refFaceTanU, ArrayFaceTangentV &refFaceTanV, const int &faceOrd, const shards::CellTopology &parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
Definition: Intrepid_CellToolsDef.hpp:2130
Intrepid::Basis_HGRAD_HEX_C2_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Definition: Intrepid_HGRAD_HEX_C2_FEM.hpp:138
Intrepid::CellTools::printSubcellVertices
static void printSubcellVertices(const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Prints the reference space coordinates of the vertices of the specified subcell.
Definition: Intrepid_CellToolsDef.hpp:3140
Intrepid::RealSpaceTools::add
static void add(Scalar *sumArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Adds contiguous data inArray1 and inArray2 of size size: sumArray = inArray1 + inArray2.
Definition: Intrepid_RealSpaceToolsDef.hpp:1506
Intrepid::CellTools::getReferenceFaceNormal
static void getReferenceFaceNormal(ArrayFaceNormal &refFaceNormal, const int &faceOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to faces of 3D reference cell.
Definition: Intrepid_CellToolsDef.hpp:2209
Intrepid::Basis_HGRAD_LINE_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Line cell.
Definition: Intrepid_HGRAD_LINE_C1_FEM.hpp:79
Intrepid::CellTools::mapToPhysicalFrame
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
Definition: Intrepid_CellToolsDef.hpp:1447
Intrepid::RealSpaceTools::vectorNorm
static Scalar vectorNorm(const Scalar *inVec, const size_t dim, const ENorm normType)
Computes norm (1, 2, infinity) of the vector inVec of size dim.
Definition: Intrepid_RealSpaceToolsDef.hpp:139
Intrepid::Basis_HGRAD_TRI_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Triangle cell...
Definition: Intrepid_HGRAD_TRI_C1_FEM.hpp:81
Intrepid::CellTools::getReferenceSubcellVertices
static void getReferenceSubcellVertices(ArraySubcellVert &subcellVertices, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Retrieves the Cartesian coordinates of all vertices of a reference subcell.
Definition: Intrepid_CellToolsDef.hpp:490
Intrepid::CellTools::getPhysicalSideNormals
static void getPhysicalSideNormals(ArraySideNormal &sideNormals, const ArrayJac &worksetJacobians, const int &worksetSideOrd, const shards::CellTopology &parentCell)
Computes non-normalized normal vectors to physical sides in a side workset .
Definition: Intrepid_CellToolsDef.hpp:2349
Intrepid::CellTools::checkPointsetInclusion
static int checkPointsetInclusion(const ArrayPoint &points, const shards::CellTopology &cellTopo, const double &threshold=INTREPID_THRESHOLD)
Checks if a set of points belongs to a reference cell.
Definition: Intrepid_CellToolsDef.hpp:2523
Intrepid::CellTools::hasReferenceCell
static int hasReferenceCell(const shards::CellTopology &cellTopo)
Checks if a cell topology has reference cell.
Definition: Intrepid_CellToolsDef.hpp:839
Intrepid::Basis_HGRAD_WEDGE_I2_FEM
Implementation of an H(grad)-compatible FEM basis of degree 2 on Wedge cell.
Definition: Intrepid_HGRAD_WEDGE_I2_FEM.hpp:110
Intrepid::Basis_HGRAD_QUAD_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Quadrilateral cell...
Definition: Intrepid_HGRAD_QUAD_C1_FEM.hpp:81
Intrepid::RealSpaceTools::inverse
static void inverse(Scalar *inverseMat, const Scalar *inMat, const size_t dim)
Computes inverse of the square matrix inMat of size dim by dim.
Definition: Intrepid_RealSpaceToolsDef.hpp:653
Intrepid::Basis_HGRAD_TET_C2_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Tetrahedron cell...
Definition: Intrepid_HGRAD_TET_C2_FEM.hpp:99
Intrepid::CellTools::mapToReferenceFrameInitGuess
static void mapToReferenceFrameInitGuess(ArrayRefPoint &refPoints, const ArrayInitGuess &initGuess, const ArrayPhysPoint &physPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computation of , the inverse of the reference-to-physical frame map using user-supplied initial guess...
Definition: Intrepid_CellToolsDef.hpp:1899
Intrepid::Basis_HGRAD_TET_COMP12_FEM
Definition: Intrepid_HGRAD_TET_COMP12_FEM.hpp:99
Intrepid::CellTools::setSubcellParametrization
static void setSubcellParametrization(FieldContainer< double > &subcellParametrization, const int subcellDim, const shards::CellTopology &parentCell)
Defines orientation-preserving parametrizations of reference edges and faces of cell topologies with ...
Definition: Intrepid_CellToolsDef.hpp:327
Intrepid::CellTools::getReferenceSubcellNodes
static void getReferenceSubcellNodes(ArraySubcellNode &subcellNodes, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Retrieves the Cartesian coordinates of all nodes of a reference subcell.
Definition: Intrepid_CellToolsDef.hpp:788
Intrepid::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition: Intrepid_Basis.hpp:89
Intrepid::CellTools::checkPointInclusion
static int checkPointInclusion(const Scalar *point, const int pointDim, const shards::CellTopology &cellTopo, const double &threshold=INTREPID_THRESHOLD)
Checks if a point belongs to a reference cell.
Definition: Intrepid_CellToolsDef.hpp:2432
Intrepid::CellTools::getReferenceVertex
static const double * getReferenceVertex(const shards::CellTopology &cell, const int vertexOrd)
Retrieves the Cartesian coordinates of a reference cell vertex.
Definition: Intrepid_CellToolsDef.hpp:471
Intrepid::CellTools::printWorksetSubcell
static void printWorksetSubcell(const ArrayCell &cellWorkset, const shards::CellTopology &parentCell, const int &pCellOrd, const int &subcellDim, const int &subcellOrd, const int &fieldWidth=3)
Prints the nodes of a subcell from a cell workset.
Definition: Intrepid_CellToolsDef.hpp:3180
Intrepid::CellTools::validateArguments_checkPointwiseInclusion
static void validateArguments_checkPointwiseInclusion(ArrayIncl &inCell, const ArrayPoint &physPoints, const ArrayCell &cellWorkset, const int &whichCell, const shards::CellTopology &cell)
Validates arguments to Intrepid::CellTools::checkPointwiseInclusion.
Definition: Intrepid_CellToolsDef.hpp:3055
Intrepid::Basis_HGRAD_QUAD_C2_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Quadrilateral cell...
Definition: Intrepid_HGRAD_QUAD_C2_FEM.hpp:91
Intrepid::FieldContainer::resize
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Definition: Intrepid_FieldContainerDef.hpp:1137
Intrepid::CellTools::validateArguments_setJacobianDetArgs
static void validateArguments_setJacobianDetArgs(const ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Validates arguments to Intrepid::CellTools::setJacobianDet.
Definition: Intrepid_CellToolsDef.hpp:2857
Intrepid::Basis_HGRAD_PYR_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Pyramid cell...
Definition: Intrepid_HGRAD_PYR_C1_FEM.hpp:83
Intrepid::RealSpaceTools::matvec
static void matvec(Scalar *matVec, const Scalar *inMat, const Scalar *inVec, const size_t dim)
Matrix-vector left multiply using contiguous data: matVec = inMat * inVec.
Definition: Intrepid_RealSpaceToolsDef.hpp:1954
Intrepid::RealSpaceTools::vecprod
static void vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight
Definition: Intrepid_RealSpaceToolsDef.hpp:2039
Intrepid::Basis_HGRAD_HEX_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell...
Definition: Intrepid_HGRAD_HEX_C1_FEM.hpp:91
Intrepid::Basis_HGRAD_PYR_I2_FEM
Implementation of an H(grad)-compatible FEM basis of degree 2 on a Pyramid cell.
Definition: Intrepid_HGRAD_PYR_I2_FEM.hpp:105
Intrepid::CellTools::validateArguments_setJacobianInv
static void validateArguments_setJacobianInv(const ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Validates arguments to Intrepid::CellTools::setJacobianInv.
Definition: Intrepid_CellToolsDef.hpp:2829
Intrepid::CellTools::validateArguments_mapToPhysicalFrame
static void validateArguments_mapToPhysicalFrame(const ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell)
Validates arguments to Intrepid::CellTools::mapToPhysicalFrame.
Definition: Intrepid_CellToolsDef.hpp:2900
Intrepid::CellTools::checkPointwiseInclusion
static void checkPointwiseInclusion(ArrayIncl &inRefCell, const ArrayPoint &points, const shards::CellTopology &cellTopo, const double &threshold=INTREPID_THRESHOLD)
Checks every point in a set for inclusion in a reference cell.
Definition: Intrepid_CellToolsDef.hpp:2564
Intrepid::CellTools::validateArguments_setJacobian
static void validateArguments_setJacobian(const ArrayJac &jacobian, const ArrayPoint &points, const ArrayCell &cellWorkset, const int &whichCell, const shards::CellTopology &cellTopo)
Validates arguments to Intrepid::CellTools::setJacobian.
Definition: Intrepid_CellToolsDef.hpp:2715
Intrepid::RealSpaceTools::subtract
static void subtract(Scalar *diffArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Subtracts contiguous data inArray2 from inArray1 of size size: diffArray = inArray1 - inArray2...
Definition: Intrepid_RealSpaceToolsDef.hpp:1638
Intrepid::CellTools::getSubcellParametrization
static const FieldContainer< double > & getSubcellParametrization(const int subcellDim, const shards::CellTopology &parentCell)
Returns array with the coefficients of the parametrization maps for the edges or faces of a reference...
Definition: Intrepid_CellToolsDef.hpp:60
Intrepid::CellTools::getBarycenter
static void getBarycenter(ArrayCent &barycenter, const ArrayCellCoord &cellCoords)
Compute cell barycenters.
Definition: Intrepid_CellToolsDef.hpp:3443
Intrepid::CellTools::getSubCVCoords
static void getSubCVCoords(ArrayCVCoord &subCVcoords, const ArrayCellCoord &cellCoords, const shards::CellTopology &primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
Definition: Intrepid_CellToolsDef.hpp:3223
Intrepid::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:109
Intrepid::CellTools::setJacobianDet
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F...
Definition: Intrepid_CellToolsDef.hpp:1314
Intrepid::CellTools::getPhysicalFaceTangents
static void getPhysicalFaceTangents(ArrayFaceTangentU &faceTanU, ArrayFaceTangentV &faceTanV, const ArrayJac &worksetJacobians, const int &worksetFaceOrd, const shards::CellTopology &parentCell)
Computes non-normalized tangent vector pairs to physical faces in a face workset ; (see Subcell works...
Definition: Intrepid_CellToolsDef.hpp:2286
Intrepid::Basis_HGRAD_TRI_C2_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Triangle cell...
Definition: Intrepid_HGRAD_TRI_C2_FEM.hpp:87
Intrepid::Basis_HGRAD_TET_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell...
Definition: Intrepid_HGRAD_TET_C1_FEM.hpp:81
Generated by   doxygen 1.8.5