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
52#include "Teuchos_oblackholestream.hpp"
53#include "Teuchos_RCP.hpp"
54#include "Teuchos_GlobalMPISession.hpp"
55#include "Teuchos_ScalarTraits.hpp"
56
57using namespace std;
58using namespace Intrepid;
59
60
79void 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
87int main(int argc, char *argv[]) {
88
89 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
90
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";
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";
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:
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
595void 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) {
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){
636 subcParamVert_A,
637 subcDim,
638 subcOrd,
639 parentCell);
640 }
641 // Domain "subcParamVert_B" is the standard 2-cube
642 else if(subcVertexCount == 4){
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
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
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
A stateless class for operations on cell data. Provides methods for:
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 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 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 int hasReferenceCell(const shards::CellTopology &cellTopo)
Checks if a cell topology has reference cell.
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 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.
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.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
static void vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight