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
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
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.
89Teuchos::RCP<Intrepid::Basis<double, Intrepid::FieldContainer<double> > >
90getIntrepidBasis(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
189int main(int argc, char *argv[]) {
190
191 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
192
193 typedef CellTools<double> CellTools;
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";
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";
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 (const 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
695void 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) {
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){
736 subcParamVert_A,
737 subcDim,
738 subcOrd,
739 parentCell);
740 }
741 // Domain "subcParamVert_B" is the standard 2-cube
742 else if(subcVertexCount == 4){
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
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_03.cpp:695
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.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
Implementation of the serendipity-family H(grad)-compatible FEM basis of degree 2 on a Hexahedron cel...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Line cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Pyramid cell.
Implementation of an H(grad)-compatible FEM basis of degree 2 on a Pyramid cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Quadrilateral cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Quadrilateral cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Tetrahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Triangle cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Triangle cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Wedge cell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Wedge cell.
Implementation of an H(grad)-compatible FEM basis of degree 2 on Wedge cell.
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.
static void vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight