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
57namespace Intrepid {
58
59 template<class Scalar>
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>
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);
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
1015if(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
1027if(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);
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
1160if(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
1172if(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
1281template<class Scalar>
1282template<class ArrayJacInv, class ArrayJac>
1283void 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/*
1291template<class Scalar>
1292template<class ArrayJacInv, class ArrayJac>
1293void CellTools<Scalar>::setJacobianInvTemp(ArrayJacInv & jacobianInv,
1294 const ArrayJac & jacobian)
1295{
1296
1297
1298 RealSpaceTools<Scalar>::inverseTemp(jacobianInv, jacobian);
1299}
1300*/
1301/*
1302template<class Scalar>
1303template<class ArrayJacDet, class ArrayJac>
1304void 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*/
1312template<class Scalar>
1313template<class ArrayJacDet, class ArrayJac>
1314void 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
1327template<class Scalar>
1328template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
1329void 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
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){
1353for(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
1445template<class Scalar>
1446template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
1447void 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
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){
1567for(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
1660template<class Scalar>
1661template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
1662void 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
1767template<class Scalar>
1768template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
1769void 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{
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
1786 FieldContainer<Scalar> jacobian;
1787 FieldContainer<Scalar> jacobInv;
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;
1872int refPointsRank=getrank(refPoints);
1873if (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
1897template<class Scalar>
1898template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
1900 const ArrayInitGuess & initGuess,
1901 const ArrayPhysPoint & physPoints,
1902 const ArrayCell & cellWorkset,
1903 const shards::CellTopology & cellTopo,
1904 const int & whichCell)
1905{
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
1916 FieldContainer<Scalar> jacobian;
1917 FieldContainer<Scalar> jacobInv;
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;
2000int refPointsRank=getrank(refPoints);
2001if (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
2025template<class Scalar>
2026template<class ArraySubcellPoint, class ArrayParamPoint>
2027void 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
2092template<class Scalar>
2093template<class ArrayEdgeTangent>
2094void 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
2128template<class Scalar>
2129template<class ArrayFaceTangentU, class ArrayFaceTangentV>
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
2173template<class Scalar>
2174template<class ArraySideNormal>
2175void 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
2207template<class Scalar>
2208template<class ArrayFaceNormal>
2209void 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
2237template<class Scalar>
2238template<class ArrayEdgeTangent, class ArrayJac>
2239void 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}
2284template<class Scalar>
2285template<class ArrayFaceTangentU, class ArrayFaceTangentV, class ArrayJac>
2286void 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
2347template<class Scalar>
2348template<class ArraySideNormal, class ArrayJac>
2349void 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
2386template<class Scalar>
2387template<class ArrayFaceNormal, class ArrayJac>
2388void 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
2431template<class Scalar>
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
2521template<class Scalar>
2522template<class ArrayPoint>
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
2562template<class Scalar>
2563template<class ArrayIncl, class ArrayPoint>
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
2660template<class Scalar>
2661template<class ArrayIncl, class ArrayPoint, class ArrayCell>
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
2713template<class Scalar>
2714template<class ArrayJac, class ArrayPoint, class ArrayCell>
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
2827template<class Scalar>
2828template<class ArrayJacInv, class ArrayJac>
2829void 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
2855template<class Scalar>
2856template<class ArrayJacDet, class ArrayJac>
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
2898template<class Scalar>
2899template<class ArrayPhysPoint, class ArrayRefPoint, class ArrayCell>
2900void 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
2924TEUCHOS_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}
2989template<class Scalar>
2990template<class ArrayRefPoint, class ArrayPhysPoint, class ArrayCell>
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
3035template<class Scalar>
3036template<class ArrayRefPoint, class ArrayInitGuess, class ArrayPhysPoint, class ArrayCell>
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
3053template<class Scalar>
3054template<class ArrayIncl, class ArrayPoint, class ArrayCell>
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
3139template<class Scalar>
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
3178template<class Scalar>
3179template<class ArrayCell>
3180void 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
#define INTREPID_MAX_NEWTON
Maximum number of Newton iterations used internally in methods such as computing the action of the in...
bool requireRankRange(std::string &errmsg, const Array &array, const int lowerBound, const int upperBound)
Checks if the rank of the array argument is in the required range.
bool requireDimensionRange(std::string &errmsg, const Array &array, const int dim, const int lowerBound, const int upperBound)
Checks if the specified array dimension is in the required range.
bool requireRankMatch(std::string &errmsg, const Array1 &array1, const Array2 &array2)
Checks if two arrays have matching ranks.
bool requireDimensionMatch(std::string &errmsg, const Array1 &array1, const int a1_dim0, const Array2 &array2, const int a2_dim0)
Checks arrays for a single matching dimension.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
Implementation of the serendipity-family H(grad)-compatible FEM basis of degree 2 on a Hexahedron cel...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Line cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Pyramid cell.
Implementation of an H(grad)-compatible FEM basis of degree 2 on a Pyramid cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Quadrilateral cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Quadrilateral cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Tetrahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Triangle cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Triangle cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Wedge cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Wedge cell.
Implementation of an H(grad)-compatible FEM basis of degree 2 on Wedge cell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
A stateless class for operations on cell data. Provides methods for:
static void validateArguments_setJacobianInv(const ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Validates arguments to Intrepid::CellTools::setJacobianInv.
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...
static const double * getReferenceVertex(const shards::CellTopology &cell, const int vertexOrd)
Retrieves the Cartesian coordinates of a reference cell vertex.
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.
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.
static void getReferenceEdgeTangent(ArrayEdgeTangent &refEdgeTangent, const int &edgeOrd, const shards::CellTopology &parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
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.
static void getReferenceSideNormal(ArraySideNormal &refSideNormal, const int &sideOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
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 ...
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.
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.
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.
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.
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.
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 .
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...
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void getSubCVCoords(ArrayCVCoord &subCVcoords, const ArrayCellCoord &cellCoords, const shards::CellTopology &primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
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...
static const double * getReferenceNode(const shards::CellTopology &cell, const int nodeOrd)
Retrieves the Cartesian coordinates of a reference cell node.
static void validateArguments_setJacobianDetArgs(const ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Validates arguments to Intrepid::CellTools::setJacobianDet.
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.
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.
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.
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...
static void getReferenceFaceNormal(ArrayFaceNormal &refFaceNormal, const int &faceOrd, const shards::CellTopology &parentCell)
Computes constant normal vectors to faces of 3D reference cell.
static int hasReferenceCell(const shards::CellTopology &cellTopo)
Checks if a cell topology has reference cell.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
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.
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 ...
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.
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.
static void getBarycenter(ArrayCent &barycenter, const ArrayCellCoord &cellCoords)
Compute cell barycenters.
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.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
Implementation of basic linear algebra functionality in Euclidean space.
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.
static void inverse(Scalar *inverseMat, const Scalar *inMat, const size_t dim)
Computes inverse of the square matrix inMat of size dim by dim.
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.
static void vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight
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.
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.