Intrepid
test_03.cpp
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 #include "Intrepid_CellTools.hpp"
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_GlobalMPISession.hpp"
55 #include "Teuchos_ScalarTraits.hpp"
56 
57 using namespace std;
58 using namespace Intrepid;
59 
60 
79 void testSubcellParametrizations(int& errorFlag,
80  const shards::CellTopology& parentCell,
81  const FieldContainer<double>& subcParamVert_A,
82  const FieldContainer<double>& subcParamVert_B,
83  const int subcDim,
84  const Teuchos::RCP<std::ostream>& outStream);
85 
86 
87 //IKT, 10/8/15: The following is a helper function that takes in a shards::CellTopology
88 //and returns a pointer to the Intrepid::Basis associated with this topology.
89 Teuchos::RCP<Intrepid::Basis<double, Intrepid::FieldContainer<double> > >
90 getIntrepidBasis(const shards::CellTopology &cellTopo)
91 {
92  Teuchos::RCP< Basis< double, FieldContainer<double> > > HGRAD_Basis;
93  // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
94  switch( cellTopo.getKey() ){
95 
96  // Standard Base topologies (number of cellWorkset = number of vertices)
97  case shards::Line<2>::key:
98  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<double, FieldContainer<double> >() );
99  break;
100 
101  case shards::Triangle<3>::key:
102  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<double, FieldContainer<double> >() );
103  break;
104 
105  case shards::Quadrilateral<4>::key:
106  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<double, FieldContainer<double> >() );
107  break;
108 
109  case shards::Tetrahedron<4>::key:
110  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<double, FieldContainer<double> >() );
111  break;
112 
113  case shards::Hexahedron<8>::key:
114  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> >() );
115  break;
116 
117  case shards::Wedge<6>::key:
118  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<double, FieldContainer<double> >() );
119  break;
120 
121  case shards::Pyramid<5>::key:
122  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<double, FieldContainer<double> >() );
123  break;
124 
125  // Standard Extended topologies
126  case shards::Triangle<6>::key:
127  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<double, FieldContainer<double> >() );
128  break;
129 
130  case shards::Quadrilateral<9>::key:
131  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<double, FieldContainer<double> >() );
132  break;
133 
134  case shards::Tetrahedron<10>::key:
135  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<double, FieldContainer<double> >() );
136  break;
137 
138  case shards::Tetrahedron<11>::key:
139  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<double, FieldContainer<double> >() );
140  break;
141 
142  case shards::Hexahedron<20>::key:
143  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<double, FieldContainer<double> >() );
144  break;
145 
146  case shards::Hexahedron<27>::key:
147  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<double, FieldContainer<double> >() );
148  break;
149 
150  case shards::Wedge<15>::key:
151  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<double, FieldContainer<double> >() );
152  break;
153 
154  case shards::Wedge<18>::key:
155  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<double, FieldContainer<double> >() );
156  break;
157 
158  case shards::Pyramid<13>::key:
159  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<double, FieldContainer<double> >() );
160  break;
161 
162  // These extended topologies are not used for mapping purposes
163  case shards::Quadrilateral<8>::key:
164  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
165  ">>> ERROR (getIntrepidBasis): Cell topology not supported. ");
166  break;
167 
168  // Base and Extended Line, Beam and Shell topologies
169  case shards::Line<3>::key:
170  case shards::Beam<2>::key:
171  case shards::Beam<3>::key:
172  case shards::ShellLine<2>::key:
173  case shards::ShellLine<3>::key:
174  case shards::ShellTriangle<3>::key:
175  case shards::ShellTriangle<6>::key:
176  case shards::ShellQuadrilateral<4>::key:
177  case shards::ShellQuadrilateral<8>::key:
178  case shards::ShellQuadrilateral<9>::key:
179  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
180  ">>> ERROR (getIntrepidBasis): Cell topology not supported. ");
181  break;
182  default:
183  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
184  ">>> ERROR (getIntrepidBasis): Cell topology not supported.");
185  }// switch
186  return HGRAD_Basis;
187 }
188 
189 int main(int argc, char *argv[]) {
190 
191  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
192 
194  typedef shards::CellTopology CellTopology;
195 
196  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
197  int iprint = argc - 1;
198 
199  Teuchos::RCP<std::ostream> outStream;
200  Teuchos::oblackholestream bhs; // outputs nothing
201 
202  if (iprint > 0)
203  outStream = Teuchos::rcp(&std::cout, false);
204  else
205  outStream = Teuchos::rcp(&bhs, false);
206 
207  // Save the format state of the original std::cout.
208  Teuchos::oblackholestream oldFormatState;
209  oldFormatState.copyfmt(std::cout);
210 
211  *outStream \
212  << "===============================================================================\n" \
213  << "| |\n" \
214  << "| Unit Test CellTools |\n" \
215  << "| |\n" \
216  << "| 1) Edge parametrizations |\n" \
217  << "| 2) Face parametrizations |\n" \
218  << "| 3) Edge tangents |\n" \
219  << "| 4) Face tangents and normals |\n" \
220  << "| |\n" \
221  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
222  << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \
223  << "| Kara Peterson (kjpeter@sandia.gov) |\n" \
224  << "| |\n" \
225  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
226  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
227  << "| |\n" \
228  << "===============================================================================\n";
229 
230  int errorFlag = 0;
231 
232 
233  // Vertices of the parametrization domain for 1-subcells: standard 1-cube [-1,1]
234  FieldContainer<double> cube_1(2, 1);
235  cube_1(0,0) = -1.0;
236  cube_1(1,0) = 1.0;
237 
238 
239  // Vertices of the parametrization domain for triangular faces: the standard 2-simplex
240  FieldContainer<double> simplex_2(3, 2);
241  simplex_2(0, 0) = 0.0; simplex_2(0, 1) = 0.0;
242  simplex_2(1, 0) = 1.0; simplex_2(1, 1) = 0.0;
243  simplex_2(2, 0) = 0.0; simplex_2(2, 1) = 1.0;
244 
245 
246  // Vertices of the parametrization domain for quadrilateral faces: the standard 2-cube
247  FieldContainer<double> cube_2(4, 2);
248  cube_2(0, 0) = -1.0; cube_2(0, 1) = -1.0;
249  cube_2(1, 0) = 1.0; cube_2(1, 1) = -1.0;
250  cube_2(2, 0) = 1.0; cube_2(2, 1) = 1.0;
251  cube_2(3, 0) = -1.0; cube_2(3, 1) = 1.0;
252 
253 
254  // Pull all available topologies from Shards
255  std::vector<shards::CellTopology> allTopologies;
256  shards::getTopologies(allTopologies);
257 
258 
259  // Set to 1 for edge and 2 for face tests
260  int subcDim;
261 
262  try{
263 
264  *outStream \
265  << "\n"
266  << "===============================================================================\n"\
267  << "| Test 1: edge parametrizations: |\n"\
268  << "===============================================================================\n\n";
269 
270  subcDim = 1;
271 
272  // Loop over the cell topologies
273  for(int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){
274 
275  // Test only 2D and 3D topologies that have reference cells, e.g., exclude Line, Pentagon, etc.
276  if(allTopologies[topoOrd].getDimension() > 1 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){
277  *outStream << " Testing edge parametrization for " << allTopologies[topoOrd].getName() <<"\n";
278  testSubcellParametrizations(errorFlag,
279  allTopologies[topoOrd],
280  cube_1,
281  cube_1,
282  subcDim,
283  outStream);
284  }
285  }
286 
287 
288  *outStream \
289  << "\n"
290  << "===============================================================================\n"\
291  << "| Test 2: face parametrizations: |\n"\
292  << "===============================================================================\n\n";
293 
294  subcDim = 2;
295 
296  // Loop over the cell topologies
297  for(int topoOrd = 0; topoOrd < (int)allTopologies.size(); topoOrd++){
298 
299  // Test only 3D topologies that have reference cells
300  if(allTopologies[topoOrd].getDimension() > 2 && CellTools::hasReferenceCell(allTopologies[topoOrd]) ){
301  *outStream << " Testing face parametrization for cell topology " << allTopologies[topoOrd].getName() <<"\n";
302  testSubcellParametrizations(errorFlag,
303  allTopologies[topoOrd],
304  simplex_2,
305  cube_2,
306  subcDim,
307  outStream);
308  }
309  }
310 
311 
312  /***********************************************************************************************
313  *
314  * Common for test 3 and 4: edge tangents and face normals for standard cells with base topo
315  *
316  **********************************************************************************************/
317 
318  // Allocate storage and extract all standard cells with base topologies
319  std::vector<shards::CellTopology> standardBaseTopologies;
320  shards::getTopologies(standardBaseTopologies, 4, shards::STANDARD_CELL, shards::BASE_TOPOLOGY);
321 
322  // Define topologies for the edge and face parametrization domains. (faces are Tri or Quad)
323  CellTopology paramEdge (shards::getCellTopologyData<shards::Line<2> >() );
324  CellTopology paramTriFace (shards::getCellTopologyData<shards::Triangle<3> >() );
325  CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
326 
327  // Define CubatureFactory:
328  DefaultCubatureFactory<double> cubFactory;
329 
330  *outStream \
331  << "\n"
332  << "===============================================================================\n"\
333  << "| Test 3: edge tangents/normals for stand. cells with base topologies: |\n"\
334  << "===============================================================================\n\n";
335  // This test loops over standard cells with base topologies, creates a set of nodes and tests tangents/normals
336  std::vector<shards::CellTopology>::iterator cti;
337 
338  // Define cubature on the edge parametrization domain:
339  Teuchos::RCP<Cubature<double> > edgeCubature = cubFactory.create(paramEdge, 6);
340  int cubDim = edgeCubature -> getDimension();
341  int numCubPoints = edgeCubature -> getNumPoints();
342 
343  // Allocate storage for cubature points and weights on edge parameter domain and fill with points:
344  FieldContainer<double> paramEdgePoints(numCubPoints, cubDim);
345  FieldContainer<double> paramEdgeWeights(numCubPoints);
346  edgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
347 
348  // Loop over admissible topologies
349  for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){
350 
351  // Exclude 0D (node), 1D (Line) and Pyramid<5> cells
352  if( ( (*cti).getDimension() >= 2) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){
353 
354  int cellDim = (*cti).getDimension();
355  int vCount = (*cti).getVertexCount();
356  FieldContainer<double> refCellVertices(vCount, cellDim);
357  CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) );
358 
359  *outStream << " Testing edge tangents";
360  if(cellDim == 2) { *outStream << " and normals"; }
361  *outStream <<" for cell topology " << (*cti).getName() <<"\n";
362 
363 
364  // Array for physical cell vertices ( must have rank 3 for setJacobians)
365  FieldContainer<double> physCellVertices(1, vCount, cellDim);
366 
367  // Randomize reference cell vertices by moving them up to +/- (1/8) units along their
368  // coordinate axis. Guaranteed to be non-degenerate for standard cells with base topology
369  for(int v = 0; v < vCount; v++){
370  for(int d = 0; d < cellDim; d++){
371  double delta = Teuchos::ScalarTraits<double>::random()/8.0;
372  physCellVertices(0, v, d) = refCellVertices(v, d) + delta;
373  } //for d
374  }// for v
375 
376  // Allocate storage for cub. points on a ref. edge; Jacobians, phys. edge tangents/normals
377  FieldContainer<double> refEdgePoints(numCubPoints, cellDim);
378  FieldContainer<double> edgePointsJacobians(1, numCubPoints, cellDim, cellDim);
379  FieldContainer<double> edgePointTangents(1, numCubPoints, cellDim);
380  FieldContainer<double> edgePointNormals(1, numCubPoints, cellDim);
381 
382  // Loop over edges:
383  for(int edgeOrd = 0; edgeOrd < (int)(*cti).getEdgeCount(); edgeOrd++){
384  /*
385  * Compute tangents on the specified physical edge using CellTools:
386  * 1. Map points from edge parametrization domain to ref. edge with specified ordinal
387  * 2. Compute parent cell Jacobians at ref. edge points
388  * 3. Compute physical edge tangents
389  */
390  CellTools::mapToReferenceSubcell(refEdgePoints, paramEdgePoints, 1, edgeOrd, (*cti) );
391  CellTools::setJacobian(edgePointsJacobians, refEdgePoints, physCellVertices, getIntrepidBasis(*cti) );
392  CellTools::getPhysicalEdgeTangents(edgePointTangents, edgePointsJacobians, edgeOrd, (*cti));
393  /*
394  * Compute tangents directly using parametrization of phys. edge and compare with CellTools tangents.
395  * 1. Get edge vertices
396  * 2. For affine edges tangent coordinates are given by F'(t) = (V1-V0)/2
397  * (for now we only test affine edges, but later we will test edges for cells
398  * with extended topologies.)
399  */
400  int v0ord = (*cti).getNodeMap(1, edgeOrd, 0);
401  int v1ord = (*cti).getNodeMap(1, edgeOrd, 1);
402 
403  for(int pt = 0; pt < numCubPoints; pt++){
404 
405  // Temp storage for directly computed edge tangents
406  FieldContainer<double> edgeBenchmarkTangents(3);
407 
408  for(int d = 0; d < cellDim; d++){
409  edgeBenchmarkTangents(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d))/2.0;
410 
411  // Compare with d-component of edge tangent by CellTools
412  if( abs(edgeBenchmarkTangents(d) - edgePointTangents(0, pt, d)) > INTREPID_THRESHOLD ){
413  errorFlag++;
414  *outStream
415  << std::setw(70) << "^^^^----FAILURE!" << "\n"
416  << " Edge tangent computation by CellTools failed for: \n"
417  << " Cell Topology = " << (*cti).getName() << "\n"
418  << " Edge ordinal = " << edgeOrd << "\n"
419  << " Edge point number = " << pt << "\n"
420  << " Tangent coordinate = " << d << "\n"
421  << " CellTools value = " << edgePointTangents(0, pt, d) << "\n"
422  << " Benchmark value = " << edgeBenchmarkTangents(d) << "\n\n";
423  }
424  } // for d
425 
426  // Test side normals for 2D cells only: edge normal has coordinates (t1, -t0)
427  if(cellDim == 2) {
428  CellTools::getPhysicalSideNormals(edgePointNormals, edgePointsJacobians, edgeOrd, (*cti));
429  if( abs(edgeBenchmarkTangents(1) - edgePointNormals(0, pt, 0)) > INTREPID_THRESHOLD ){
430  errorFlag++;
431  *outStream
432  << std::setw(70) << "^^^^----FAILURE!" << "\n"
433  << " Edge Normal computation by CellTools failed for: \n"
434  << " Cell Topology = " << (*cti).getName() << "\n"
435  << " Edge ordinal = " << edgeOrd << "\n"
436  << " Edge point number = " << pt << "\n"
437  << " Normal coordinate = " << 0 << "\n"
438  << " CellTools value = " << edgePointNormals(0, pt, 0) << "\n"
439  << " Benchmark value = " << edgeBenchmarkTangents(1) << "\n\n";
440  }
441  if( abs(edgeBenchmarkTangents(0) + edgePointNormals(0, pt, 1)) > INTREPID_THRESHOLD ){
442  errorFlag++;
443  *outStream
444  << std::setw(70) << "^^^^----FAILURE!" << "\n"
445  << " Edge Normal computation by CellTools failed for: \n"
446  << " Cell Topology = " << (*cti).getName() << "\n"
447  << " Edge ordinal = " << edgeOrd << "\n"
448  << " Edge point number = " << pt << "\n"
449  << " Normal coordinate = " << 1 << "\n"
450  << " CellTools value = " << edgePointNormals(0, pt, 1) << "\n"
451  << " Benchmark value = " << -edgeBenchmarkTangents(0) << "\n\n";
452  }
453  } // edge normals
454  } // for pt
455  }// for edgeOrd
456  }// if admissible cell
457  }// for cti
458 
459 
460 
461  *outStream \
462  << "\n"
463  << "===============================================================================\n"\
464  << "| Test 4: face/side normals for stand. 3D cells with base topologies: | |\n"\
465  << "===============================================================================\n\n";
466  // This test loops over standard 3D cells with base topologies, creates a set of nodes and tests normals
467 
468  // Define cubature on the edge parametrization domain:
469  Teuchos::RCP<Cubature<double> > triFaceCubature = cubFactory.create(paramTriFace, 6);
470  Teuchos::RCP<Cubature<double> > quadFaceCubature = cubFactory.create(paramQuadFace, 6);
471 
472  int faceCubDim = triFaceCubature -> getDimension();
473  int numTriFaceCubPoints = triFaceCubature -> getNumPoints();
474  int numQuadFaceCubPoints = quadFaceCubature -> getNumPoints();
475 
476  // Allocate storage for cubature points and weights on face parameter domain and fill with points:
477  FieldContainer<double> paramTriFacePoints(numTriFaceCubPoints, faceCubDim);
478  FieldContainer<double> paramTriFaceWeights(numTriFaceCubPoints);
479  FieldContainer<double> paramQuadFacePoints(numQuadFaceCubPoints, faceCubDim);
480  FieldContainer<double> paramQuadFaceWeights(numQuadFaceCubPoints);
481 
482  triFaceCubature -> getCubature(paramTriFacePoints, paramTriFaceWeights);
483  quadFaceCubature -> getCubature(paramQuadFacePoints, paramQuadFaceWeights);
484 
485  // Loop over admissible topologies
486  for(cti = standardBaseTopologies.begin(); cti !=standardBaseTopologies.end(); ++cti){
487 
488  // Exclude 2D and Pyramid<5> cells
489  if( ( (*cti).getDimension() == 3) && ( (*cti).getKey() != shards::Pyramid<5>::key) ){
490 
491  int cellDim = (*cti).getDimension();
492  int vCount = (*cti).getVertexCount();
493  FieldContainer<double> refCellVertices(vCount, cellDim);
494  CellTools::getReferenceSubcellVertices(refCellVertices, cellDim, 0, (*cti) );
495 
496  *outStream << " Testing face/side normals for cell topology " << (*cti).getName() <<"\n";
497 
498  // Array for physical cell vertices ( must have rank 3 for setJacobians)
499  FieldContainer<double> physCellVertices(1, vCount, cellDim);
500 
501  // Randomize reference cell vertices by moving them up to +/- (1/8) units along their
502  // coordinate axis. Guaranteed to be non-degenerate for standard cells with base topology
503  for(int v = 0; v < vCount; v++){
504  for(int d = 0; d < cellDim; d++){
505  double delta = Teuchos::ScalarTraits<double>::random()/8.0;
506  physCellVertices(0, v, d) = refCellVertices(v, d) + delta;
507  } //for d
508  }// for v
509 
510  // Allocate storage for cub. points on a ref. face; Jacobians, phys. face normals and
511  // benchmark normals.
512  FieldContainer<double> refTriFacePoints(numTriFaceCubPoints, cellDim);
513  FieldContainer<double> refQuadFacePoints(numQuadFaceCubPoints, cellDim);
514  FieldContainer<double> triFacePointsJacobians(1, numTriFaceCubPoints, cellDim, cellDim);
515  FieldContainer<double> quadFacePointsJacobians(1, numQuadFaceCubPoints, cellDim, cellDim);
516  FieldContainer<double> triFacePointNormals(1, numTriFaceCubPoints, cellDim);
517  FieldContainer<double> triSidePointNormals(1, numTriFaceCubPoints, cellDim);
518  FieldContainer<double> quadFacePointNormals(1, numQuadFaceCubPoints, cellDim);
519  FieldContainer<double> quadSidePointNormals(1, numQuadFaceCubPoints, cellDim);
520 
521 
522  // Loop over faces:
523  for(int faceOrd = 0; faceOrd < (int)(*cti).getSideCount(); faceOrd++){
524 
525  // This test presently includes only Triangle<3> and Quadrilateral<4> faces. Once we support
526  // cells with extended topologies we will add their faces as well.
527  switch( (*cti).getCellTopologyData(2, faceOrd) -> key ) {
528 
529  case shards::Triangle<3>::key:
530  {
531  // Compute face normals using CellTools
532  CellTools::mapToReferenceSubcell(refTriFacePoints, paramTriFacePoints, 2, faceOrd, (*cti) );
533  CellTools::setJacobian(triFacePointsJacobians, refTriFacePoints, physCellVertices, getIntrepidBasis(*cti) );
534  CellTools::getPhysicalFaceNormals(triFacePointNormals, triFacePointsJacobians, faceOrd, (*cti));
535  CellTools::getPhysicalSideNormals(triSidePointNormals, triFacePointsJacobians, faceOrd, (*cti));
536  /*
537  * Compute face normals using direct linear parametrization of the face: the map from
538  * standard 2-simplex to physical Triangle<3> face in 3D is
539  * F(x,y) = V0 + (V1-V0)x + (V2-V0)*y
540  * Face normal is vector product Tx X Ty where Tx = (V1-V0); Ty = (V2-V0)
541  */
542  int v0ord = (*cti).getNodeMap(2, faceOrd, 0);
543  int v1ord = (*cti).getNodeMap(2, faceOrd, 1);
544  int v2ord = (*cti).getNodeMap(2, faceOrd, 2);
545 
546  // Loop over face points: redundant for affine faces, but CellTools gives one vector
547  // per point so need to check all points anyways.
548  for(int pt = 0; pt < numTriFaceCubPoints; pt++){
549  FieldContainer<double> tanX(3), tanY(3), faceNormal(3);
550  for(int d = 0; d < cellDim; d++){
551  tanX(d) = (physCellVertices(0, v1ord, d) - physCellVertices(0, v0ord, d));
552  tanY(d) = (physCellVertices(0, v2ord, d) - physCellVertices(0, v0ord, d));
553  }// for d
554 
555  RealSpaceTools<double>::vecprod(faceNormal, tanX, tanY);
556 
557  // Compare direct normal with d-component of the face/side normal by CellTools
558  for(int d = 0; d < cellDim; d++){
559 
560  // face normal method
561  if( abs(faceNormal(d) - triFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
562  errorFlag++;
563  *outStream
564  << std::setw(70) << "^^^^----FAILURE!" << "\n"
565  << " Face normal computation by CellTools failed for: \n"
566  << " Cell Topology = " << (*cti).getName() << "\n"
567  << " Face Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name << "\n"
568  << " Face ordinal = " << faceOrd << "\n"
569  << " Face point number = " << pt << "\n"
570  << " Normal coordinate = " << d << "\n"
571  << " CellTools value = " << triFacePointNormals(0, pt, d)
572  << " Benchmark value = " << faceNormal(d) << "\n\n";
573  }
574  //side normal method
575  if( abs(faceNormal(d) - triSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
576  errorFlag++;
577  *outStream
578  << std::setw(70) << "^^^^----FAILURE!" << "\n"
579  << " Side normal computation by CellTools failed for: \n"
580  << " Cell Topology = " << (*cti).getName() << "\n"
581  << " Side Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name << "\n"
582  << " Side ordinal = " << faceOrd << "\n"
583  << " Side point number = " << pt << "\n"
584  << " Normal coordinate = " << d << "\n"
585  << " CellTools value = " << triSidePointNormals(0, pt, d)
586  << " Benchmark value = " << faceNormal(d) << "\n\n";
587  }
588  } // for d
589  } // for pt
590  }
591  break;
592 
593  case shards::Quadrilateral<4>::key:
594  {
595  // Compute face normals using CellTools
596  CellTools::mapToReferenceSubcell(refQuadFacePoints, paramQuadFacePoints, 2, faceOrd, (*cti) );
597  CellTools::setJacobian(quadFacePointsJacobians, refQuadFacePoints, physCellVertices, getIntrepidBasis(*cti) );
598  CellTools::getPhysicalFaceNormals(quadFacePointNormals, quadFacePointsJacobians, faceOrd, (*cti));
599  CellTools::getPhysicalSideNormals(quadSidePointNormals, quadFacePointsJacobians, faceOrd, (*cti));
600  /*
601  * Compute face normals using direct bilinear parametrization of the face: the map from
602  * [-1,1]^2 to physical Quadrilateral<4> face in 3D is
603  * F(x,y) = ((V0+V1+V2+V3) + (-V0+V1+V2-V3)*X + (-V0-V1+V2+V3)*Y + (V0-V1+V2-V3)*X*Y)/4
604  * Face normal is vector product Tx X Ty where
605  * Tx = ((-V0+V1+V2-V3) + (V0-V1+V2-V3)*Y)/4
606  * Ty = ((-V0-V1+V2+V3) + (V0-V1+V2-V3)*X)/4
607  */
608  int v0ord = (*cti).getNodeMap(2, faceOrd, 0);
609  int v1ord = (*cti).getNodeMap(2, faceOrd, 1);
610  int v2ord = (*cti).getNodeMap(2, faceOrd, 2);
611  int v3ord = (*cti).getNodeMap(2, faceOrd, 3);
612 
613  // Loop over face points (redundant for affine faces, but needed for later when we handle non-affine ones)
614  for(int pt = 0; pt < numTriFaceCubPoints; pt++){
615  FieldContainer<double> tanX(3), tanY(3), faceNormal(3);
616  for(int d = 0; d < cellDim; d++){
617  tanX(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,1) ) +
618  physCellVertices(0, v1ord, d)*( 1.0 - paramQuadFacePoints(pt,1) ) +
619  physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,1) ) +
620  physCellVertices(0, v3ord, d)*(-1.0 - paramQuadFacePoints(pt,1) ) )/4.0;
621 
622  tanY(d) = (physCellVertices(0, v0ord, d)*(-1.0 + paramQuadFacePoints(pt,0) ) +
623  physCellVertices(0, v1ord, d)*(-1.0 - paramQuadFacePoints(pt,0) ) +
624  physCellVertices(0, v2ord, d)*( 1.0 + paramQuadFacePoints(pt,0) ) +
625  physCellVertices(0, v3ord, d)*( 1.0 - paramQuadFacePoints(pt,0) ) )/4.0;
626  }// for d
627 
628  RealSpaceTools<double>::vecprod(faceNormal, tanX, tanY);
629  // Compare direct normal with d-component of the face/side normal by CellTools
630  for(int d = 0; d < cellDim; d++){
631 
632  // face normal method
633  if( abs(faceNormal(d) - quadFacePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
634  errorFlag++;
635  *outStream
636  << std::setw(70) << "^^^^----FAILURE!" << "\n"
637  << " Face normal computation by CellTools failed for: \n"
638  << " Cell Topology = " << (*cti).getName() << "\n"
639  << " Face Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name << "\n"
640  << " Face ordinal = " << faceOrd << "\n"
641  << " Face point number = " << pt << "\n"
642  << " Normal coordinate = " << d << "\n"
643  << " CellTools value = " << quadFacePointNormals(0, pt, d)
644  << " Benchmark value = " << faceNormal(d) << "\n\n";
645  }
646  //side normal method
647  if( abs(faceNormal(d) - quadSidePointNormals(0, pt, d)) > INTREPID_THRESHOLD ){
648  errorFlag++;
649  *outStream
650  << std::setw(70) << "^^^^----FAILURE!" << "\n"
651  << " Side normal computation by CellTools failed for: \n"
652  << " Cell Topology = " << (*cti).getName() << "\n"
653  << " Side Topology = " << (*cti).getCellTopologyData(2, faceOrd) -> name << "\n"
654  << " Side ordinal = " << faceOrd << "\n"
655  << " Side point number = " << pt << "\n"
656  << " Normal coordinate = " << d << "\n"
657  << " CellTools value = " << quadSidePointNormals(0, pt, d)
658  << " Benchmark value = " << faceNormal(d) << "\n\n";
659  }
660  } // for d
661  }// for pt
662  }// case Quad
663  break;
664  default:
665  errorFlag++;
666  *outStream << " Face normals test failure: face topology not supported \n\n";
667  } // switch
668  }// for faceOrd
669  }// if admissible
670  }// for cti
671  }// try
672 
673  //============================================================================================//
674  // Wrap up test: check if the test broke down unexpectedly due to an exception //
675  //============================================================================================//
676  catch (std::logic_error err) {
677  *outStream << err.what() << "\n";
678  errorFlag = -1000;
679  };
680 
681 
682  if (errorFlag != 0)
683  std::cout << "End Result: TEST FAILED\n";
684  else
685  std::cout << "End Result: TEST PASSED\n";
686 
687  // reset format state of std::cout
688  std::cout.copyfmt(oldFormatState);
689 
690  return errorFlag;
691 }
692 
693 
694 
695 void testSubcellParametrizations(int& errorFlag,
696  const shards::CellTopology& parentCell,
697  const FieldContainer<double>& subcParamVert_A,
698  const FieldContainer<double>& subcParamVert_B,
699  const int subcDim,
700  const Teuchos::RCP<std::ostream>& outStream){
701 
702  // Get cell dimension and subcell count
703  int cellDim = parentCell.getDimension();
704  int subcCount = parentCell.getSubcellCount(subcDim);
705 
706 
707  // Loop over subcells of the specified dimension
708  for(int subcOrd = 0; subcOrd < subcCount; subcOrd++){
709  int subcVertexCount = parentCell.getVertexCount(subcDim, subcOrd);
710 
711 
712  // Storage for correct reference subcell vertices and for the images of the parametrization domain points
713  FieldContainer<double> refSubcellVertices(subcVertexCount, cellDim);
714  FieldContainer<double> mappedParamVertices(subcVertexCount, cellDim);
715 
716 
717  // Retrieve correct reference subcell vertices
718  CellTools<double>::getReferenceSubcellVertices(refSubcellVertices, subcDim, subcOrd, parentCell);
719 
720 
721  // Map vertices of the parametrization domain to 1 or 2-subcell with ordinal subcOrd
722  // For edges parametrization domain is always 1-cube passed as "subcParamVert_A"
723  if(subcDim == 1) {
724  CellTools<double>::mapToReferenceSubcell(mappedParamVertices,
725  subcParamVert_A,
726  subcDim,
727  subcOrd,
728  parentCell);
729  }
730  // For faces need to treat Triangle and Quadrilateral faces separately
731  else if(subcDim == 2) {
732 
733  // domain "subcParamVert_A" is the standard 2-simplex
734  if(subcVertexCount == 3){
735  CellTools<double>::mapToReferenceSubcell(mappedParamVertices,
736  subcParamVert_A,
737  subcDim,
738  subcOrd,
739  parentCell);
740  }
741  // Domain "subcParamVert_B" is the standard 2-cube
742  else if(subcVertexCount == 4){
743  CellTools<double>::mapToReferenceSubcell(mappedParamVertices,
744  subcParamVert_B,
745  subcDim,
746  subcOrd,
747  parentCell);
748  }
749  }
750 
751  // Compare the images of the parametrization domain vertices with the true vertices.
752  for(int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){
753  for(int dim = 0; dim < cellDim; dim++){
754 
755  if(mappedParamVertices(subcVertOrd, dim) != refSubcellVertices(subcVertOrd, dim) ) {
756  errorFlag++;
757  *outStream
758  << std::setw(70) << "^^^^----FAILURE!" << "\n"
759  << " Cell Topology = " << parentCell.getName() << "\n"
760  << " Parametrization of subcell " << subcOrd << " which is "
761  << parentCell.getName(subcDim,subcOrd) << " failed for vertex " << subcVertOrd << ":\n"
762  << " parametrization map fails to map correctly coordinate " << dim << " of that vertex\n\n";
763 
764  }//if
765  }// for dim
766  }// for subcVertOrd
767  }// for subcOrd
768 
769 }
770 
771 
772 
773 
774 
775 
Implementation of basic linear algebra functionality in Euclidean space.
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 Wedge cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Wedge cell.
Header file for the Intrepid::CellTools class.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Header file for utility class to provide multidimensional containers.
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 Triangle cell...
Implementation of an H(grad)-compatible FEM basis of degree 2 on Wedge cell.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
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 Tetrahedron 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 Pyramid cell...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell...
Implementation of an H(grad)-compatible FEM basis of degree 2 on a Pyramid cell.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
A stateless class for operations on cell data. Provides methods for:
void testSubcellParametrizations(int &errorFlag, const shards::CellTopology &parentCell, const FieldContainer< double > &subcParamVert_A, const FieldContainer< double > &subcParamVert_B, const int subcDim, const Teuchos::RCP< std::ostream > &outStream)
Maps the vertices of the subcell parametrization domain to that subcell.
Definition: test_01.cpp:595
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 Tetrahedron cell...