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