Intrepid
example_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 "Shards_CellTopology.hpp"
53#include "Teuchos_GlobalMPISession.hpp"
54
55using namespace std;
56using namespace Intrepid;
57
58int main(int argc, char *argv[]) {
59
60 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
61
62 typedef CellTools<double> CellTools;
63 typedef shards::CellTopology CellTopology;
64
65 std::cout \
66 << "===============================================================================\n" \
67 << "| |\n" \
68 << "| Example use of the CellTools class |\n" \
69 << "| |\n" \
70 << "| 1) Definition of integration points on edge and face worksets |\n" \
71 << "| 2) Computation of face normals and edge tangents on face and edge worksets |\n" \
72 << "| 2) Computation of side normals on face and edge worksets |\n" \
73 << "| |\n" \
74 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
75 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \
76 << "| Kara Peterson (kjpeter@sandia.gov) |\n" \
77 << "| |\n" \
78 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
79 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
80 << "| |\n" \
81 << "===============================================================================\n"\
82 << "| EXAMPLE 1: Definition of integration points on a Triangle edge workset |\n"\
83 << "===============================================================================\n";
84
85 /*
86 * 1. Common tasks for getting integration points on edge worksets
87 */
88
89 // Step 1.a: Define CubatureFactory
90 DefaultCubatureFactory<double> cubFactory;
91
92 // Step 1.b: Define topology of the edge parametrization domain [-1,1]
93 CellTopology paramEdge(shards::getCellTopologyData<shards::Line<2> >() );
94
95 // Step 1.c: Define storage for cubature points and weights on [-1,1]
96 FieldContainer<double> paramEdgeWeights;
97 FieldContainer<double> paramEdgePoints;
98
99 // Step 1.d: Define storage for cubature points on a reference edge
100 FieldContainer<double> refEdgePoints;
101
102 // Step 1.f: Define storage for cubature points on workset edges
103 FieldContainer<double> worksetEdgePoints;
104
105
106
107 /*
108 * 2. Define edge workset comprising of 4 edges corresponding to reference edge 2 on Triangle<3>
109 */
110
111 // Step 2.a: Specify cell topology of the parent cell
112 CellTopology triangle_3( shards::getCellTopologyData<shards::Triangle<3> >() );
113
114 // Step 2.b: Specify the vertices of 4 parent Triangle<3> cells
115 int worksetSize = 4;
116 int pCellNodeCount = triangle_3.getVertexCount();
117 int pCellDim = triangle_3.getDimension();
118
119 FieldContainer<double> triNodes(worksetSize, pCellNodeCount, pCellDim);
120 // Triangle 0
121 triNodes(0, 0, 0) = 0.5; triNodes(0, 0, 1) = 2.0;
122 triNodes(0, 1, 0) = 0.0; triNodes(0, 1, 1) = 1.0;
123 triNodes(0, 2, 0) = 1.0; triNodes(0, 2, 1) = 1.0;
124 // Triangle 1
125 triNodes(1, 0, 0) = 1.0; triNodes(1, 0, 1) = 1.0;
126 triNodes(1, 1, 0) = 1.0; triNodes(1, 1, 1) = 0.0;
127 triNodes(1, 2, 0) = 0.0; triNodes(1, 2, 1) = 0.0;
128 // Triangle 2
129 triNodes(2, 0, 0) = 0.0; triNodes(2, 0, 1) = 0.0;
130 triNodes(2, 1, 0) = 1.0; triNodes(2, 1, 1) = 1.0;
131 triNodes(2, 2, 0) = 0.0; triNodes(2, 2, 1) = 1.0;
132 // Triangle 3
133 triNodes(3, 0, 0) = 1.0; triNodes(3, 0, 1) = 1.0;
134 triNodes(3, 1, 0) = 2.5; triNodes(3, 1, 1) = 1.5;
135 triNodes(3, 2, 0) = 0.5; triNodes(3, 2, 1) = 2.0;
136
137 // Step 2.c: Specify dimension and ordinal (relative to ref. cell) of the subcells in the workset
138 int subcellDim = 1;
139 int subcellOrd = 2;
140
141
142
143 /*
144 * 3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1]
145 */
146
147 // Step 3.a: selects Gauss rule of order 6 on [-1,1]
148 Teuchos::RCP<Cubature<double> > triEdgeCubature = cubFactory.create(paramEdge, 6);
149
150 // Step 3.b allocate storage for cubature points on [-1,1]
151 int cubDim = triEdgeCubature -> getDimension();
152 int numCubPoints = triEdgeCubature -> getNumPoints();
153
154 // Arrays must be properly sized for the specified set of Gauss points
155 paramEdgePoints.resize(numCubPoints, cubDim);
156 paramEdgeWeights.resize(numCubPoints);
157 triEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
158
159
160
161 /*
162 * 4. Map Gauss points from [-1,1] to every edge in the workset
163 */
164
165 // Step 4.a: Allocate storage for integration points on the reference edge
166 refEdgePoints.resize(numCubPoints, pCellDim);
167
168 // Step 4.b: Allocate storage for integration points on the edge workset
169 worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
170
171 // Step 4.c: Map Gauss points to reference edge: paramEdgePoints -> refEdgePoints
173 paramEdgePoints,
174 subcellDim,
175 subcellOrd,
176 triangle_3);
177
178 // Step 4.d: Map Gauss points from ref. edge to edge workset: refEdgePoints -> worksetEdgePoints
179 CellTools::mapToPhysicalFrame(worksetEdgePoints,
180 refEdgePoints,
181 triNodes,
182 triangle_3);
183
184
185
186 /*
187 * 5. Print Gauss points on the edges of the workset
188 */
189 for(int pCell = 0; pCell < worksetSize; pCell++){
190
191 CellTools::printWorksetSubcell(triNodes, triangle_3, pCell, subcellDim, subcellOrd);
192
193 for(int pt = 0; pt < numCubPoints; pt++){
194 std::cout << "\t 1D Gauss point "
195 << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) << " --> " << "("
196 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) << ", "
197 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) << ")\n";
198 }
199 std::cout << "\n";
200 }//pCell
201
202
203
204 std::cout << "\n" \
205 << "===============================================================================\n"\
206 << "| EXAMPLE 2: Definition of integration points on a Quadrilateral edge workset|\n"\
207 << "===============================================================================\n";
208
209 /*
210 * 2. Define edge workset comprising of 2 edges corresponding to reference edge 2 on Quadrilateral<4>
211 * (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
212 */
213
214 // Step 2.a: Specify cell topology of the parent cell
215 CellTopology quadrilateral_4( shards::getCellTopologyData<shards::Quadrilateral<4> >() );
216
217 // Step 2.b: Specify the vertices of 2 parent Quadrilateral<4> cells
218 worksetSize = 2;
219 pCellNodeCount = quadrilateral_4.getVertexCount();
220 pCellDim = quadrilateral_4.getDimension();
221
222 FieldContainer<double> quadNodes(worksetSize, pCellNodeCount, pCellDim);
223 // Quadrilateral 0
224 quadNodes(0, 0, 0) = 1.00; quadNodes(0, 0, 1) = 1.00;
225 quadNodes(0, 1, 0) = 2.00; quadNodes(0, 1, 1) = 0.75;
226 quadNodes(0, 2, 0) = 1.75; quadNodes(0, 2, 1) = 2.00;
227 quadNodes(0, 3, 0) = 1.25; quadNodes(0, 3, 1) = 2.00;
228 // Quadrilateral 1
229 quadNodes(1, 0, 0) = 2.00; quadNodes(1, 0, 1) = 0.75;
230 quadNodes(1, 1, 0) = 3.00; quadNodes(1, 1, 1) = 1.25;
231 quadNodes(1, 2, 0) = 2.75; quadNodes(1, 2, 1) = 2.25;
232 quadNodes(1, 3, 0) = 1.75; quadNodes(1, 3, 1) = 2.00;
233
234 // Step 2.c: Specify dimension and ordinal (relative to ref. cell) of the subcells in the workset
235 subcellDim = 1;
236 subcellOrd = 2;
237
238 /*
239 * 3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1]
240 */
241
242 // Step 3.a: selects Gauss rule of order 4 on [-1,1]
243 Teuchos::RCP<Cubature<double> > quadEdgeCubature = cubFactory.create(paramEdge, 6);
244
245 // Step 3.b: allocate storage for cubature points
246 cubDim = quadEdgeCubature -> getDimension();
247 numCubPoints = quadEdgeCubature -> getNumPoints();
248
249 // Arrays must be properly sized for the specified set of Gauss points
250 paramEdgePoints.resize(numCubPoints, cubDim);
251 paramEdgeWeights.resize(numCubPoints);
252 quadEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
253
254 /*
255 * 4. Map Gauss points from [-1,1] to every edge in the workset
256 */
257
258 // Step 4.a: Allocate storage for integration points on the reference edge
259 refEdgePoints.resize(numCubPoints, pCellDim);
260
261 // Step 4.b: Allocate storage for integration points on the edge workset
262 worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
263
264 // Step 4.c: Map Gauss points to reference edge: paramEdgePoints -> refEdgePoints
266 paramEdgePoints,
267 subcellDim,
268 subcellOrd,
269 quadrilateral_4);
270
271 // Step 4.d: Map Gauss points from ref. edge to edge workset: refEdgePoints -> worksetEdgePoints
272 CellTools::mapToPhysicalFrame(worksetEdgePoints,
273 refEdgePoints,
274 quadNodes,
275 quadrilateral_4);
276
277 /*
278 * 5. Print Gauss points on the edges of the workset
279 */
280 for(int pCell = 0; pCell < worksetSize; pCell++){
281
282 CellTools::printWorksetSubcell(quadNodes, quadrilateral_4, pCell, subcellDim, subcellOrd);
283
284 for(int pt = 0; pt < numCubPoints; pt++){
285 std::cout << "\t 1D Gauss point "
286 << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) << " --> " << "("
287 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) << ", "
288 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) << ")\n";
289 }
290 std::cout << "\n";
291 }//pCell
292
293
294
295 std::cout << "\n" \
296 << "===============================================================================\n"\
297 << "| EXAMPLE 3: Edge tangents at Gauss points on a Quadrilateral edge workset |\n"\
298 << "===============================================================================\n";
299
300 /* This task requires Gauss points on edge parametrization domain [-1,1] and on the
301 * reference edge whose ordinal matches the edge workset ordinal. This repeats the first few
302 * steps from Example 2:
303 *
304 * 1. Define cubature factory and topology for edge parametrization domain [-1,1];
305 * (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
306 * 2. Define an edge workset;
307 * 3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1];
308 * 4. Map Gauss points from [-1,1] to reference edge
309 * NOTE: this example only demonstrates computation of the edge tangents and so, Gauss
310 * points on the edge workset are not needed. Thus we skip mapping Gauss points from
311 * reference edge to edge workset
312 *
313 * 5. Compute Jacobians at Gauss points on reference edge for all cells in the workset:
314 */
315
316 // Step 5.a: Define and allocate storage for workset Jacobians
317 FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim);
318
319 // Step 5.b: Compute Jacobians at Gauss pts. on reference edge for all parent cells:
320 CellTools::setJacobian(worksetJacobians,
321 refEdgePoints,
322 quadNodes,
323 quadrilateral_4);
324 /*
325 * 6. Get the (non-normalized) edge tangents for the edge workset:
326 */
327 // Step 6.a: Allocate storage for edge tangents
328 FieldContainer<double> edgeTangents(worksetSize, numCubPoints, pCellDim);
329
330 // Step 6.b: Compute the edge tangents:
332 worksetJacobians,
333 subcellOrd,
334 quadrilateral_4);
335
336 // Step 6.c: Print edge tangents at Gauss points on workset edges (these Gauss points were computed in Example 2)
337 std::cout
338 << "Edge tangents computed by CellTools::getPhysicalEdgeTangents.\n"
339 << "Local edge ordinal = " << subcellOrd <<"\n";
340
341 for(int pCell = 0; pCell < worksetSize; pCell++){
342
343 CellTools::printWorksetSubcell(quadNodes, quadrilateral_4, pCell, subcellDim, subcellOrd);
344
345 for(int pt = 0; pt < numCubPoints; pt++){
346 std::cout << "\t 2D Gauss point: ("
347 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) << ", "
348 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) << ") "
349 << std::setw(8) << " edge tangent: " << "("
350 << std::setw(8) << std::right << edgeTangents(pCell, pt, 0) << ", "
351 << std::setw(8) << std::right << edgeTangents(pCell, pt, 1) << ")\n";
352 }
353 std::cout << "\n";
354 }//pCell
355
356 std::cout << "\n" \
357 << "===============================================================================\n"\
358 << "| EXAMPLE 4: Side normals at Gauss points on a Quadrilateral side workset |\n"\
359 << "===============================================================================\n";
360
361 /* For this task we reuse the edge workset from Example 3 as a side workset and compute
362 * normals to the sides in that set. This task repeats the first 5 steps from Example 3
363 * The only difference is that in Step 6 we call CellTools::getPhysicalSideNormals
364 */
365 /*
366 * 6. Get the (non-normalized) side normals for the side (edge) workset:
367 */
368 // Step 6.a: Allocate storage for side normals
369 FieldContainer<double> sideNormals(worksetSize, numCubPoints, pCellDim);
370
371 // Step 6.b: Compute the side normals:
373 worksetJacobians,
374 subcellOrd,
375 quadrilateral_4);
376
377 // Step 6.c: Print side normals at Gauss points on workset sides (these Gauss points were computed in Example 2)
378 std::cout
379 << "Side normals computed by CellTools::getPhysicalSideNormals.\n"
380 << "Local side (edge) ordinal = " << subcellOrd <<"\n";
381
382 for(int pCell = 0; pCell < worksetSize; pCell++){
383
384 CellTools::printWorksetSubcell(quadNodes, quadrilateral_4, pCell, subcellDim, subcellOrd);
385
386 for(int pt = 0; pt < numCubPoints; pt++){
387 std::cout << "\t 2D Gauss point: ("
388 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) << ", "
389 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) << ") "
390 << std::setw(8) << " side normal: " << "("
391 << std::setw(8) << std::right << sideNormals(pCell, pt, 0) << ", "
392 << std::setw(8) << std::right << sideNormals(pCell, pt, 1) << ")\n";
393 }
394 std::cout << "\n";
395 }//pCell
396
397
398 std::cout << "\n" \
399 << "===============================================================================\n"\
400 << "| EXAMPLE 5: Definition of integration points on a Hexahedron edge workset |\n"\
401 << "===============================================================================\n";
402
403 /*
404 * 2. Define edge workset comprising of 1 edge corresponding to reference edge 10 on Hexahedron<8>
405 * (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
406 */
407
408 // Step 2.a: Specify cell topology of the parent cell
409 CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() );
410
411 // Step 2.b: Specify the vertices of the parent Hexahedron<8> cell
412 worksetSize = 1;
413 pCellNodeCount = hexahedron_8.getVertexCount();
414 pCellDim = hexahedron_8.getDimension();
415
416 FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim);
417 // bottom face vertices
418 hexNodes(0, 0, 0) = 0.00; hexNodes(0, 0, 1) = 0.00, hexNodes(0, 0, 2) = 0.00;
419 hexNodes(0, 1, 0) = 1.00; hexNodes(0, 1, 1) = 0.00, hexNodes(0, 1, 2) = 0.00;
420 hexNodes(0, 2, 0) = 1.00; hexNodes(0, 2, 1) = 1.00, hexNodes(0, 2, 2) = 0.00;
421 hexNodes(0, 3, 0) = 0.00; hexNodes(0, 3, 1) = 1.00, hexNodes(0, 3, 2) = 0.00;
422 // top face vertices
423 hexNodes(0, 4, 0) = 0.00; hexNodes(0, 4, 1) = 0.00, hexNodes(0, 4, 2) = 1.00;
424 hexNodes(0, 5, 0) = 1.00; hexNodes(0, 5, 1) = 0.00, hexNodes(0, 5, 2) = 1.00;
425 hexNodes(0, 6, 0) = 1.00; hexNodes(0, 6, 1) = 1.00, hexNodes(0, 6, 2) = 0.75;
426 hexNodes(0, 7, 0) = 0.00; hexNodes(0, 7, 1) = 1.00, hexNodes(0, 7, 2) = 1.00;
427
428 // An alternative hex obtained by intersection of the unit cube [0,1]^3 and the plane
429 // z = 1 - 1/4x - 1/4y. The top face (local ordinal 5) of the resulting hex lies in this plane
430 // and has the same normal vector parallel to (1/4,1/4,1). This workset allows to test face normals
431 FieldContainer<double> hexNodesAlt(worksetSize, pCellNodeCount, pCellDim);
432 // bottom face vertices
433 hexNodesAlt(0, 0, 0) = 0.00; hexNodesAlt(0, 0, 1) = 0.00, hexNodesAlt(0, 0, 2) = 0.00;
434 hexNodesAlt(0, 1, 0) = 1.00; hexNodesAlt(0, 1, 1) = 0.00, hexNodesAlt(0, 1, 2) = 0.00;
435 hexNodesAlt(0, 2, 0) = 1.00; hexNodesAlt(0, 2, 1) = 1.00, hexNodesAlt(0, 2, 2) = 0.00;
436 hexNodesAlt(0, 3, 0) = 0.00; hexNodesAlt(0, 3, 1) = 1.00, hexNodesAlt(0, 3, 2) = 0.00;
437 // top face vertices
438 hexNodesAlt(0, 4, 0) = 0.00; hexNodesAlt(0, 4, 1) = 0.00, hexNodesAlt(0, 4, 2) = 1.00;
439 hexNodesAlt(0, 5, 0) = 1.00; hexNodesAlt(0, 5, 1) = 0.00, hexNodesAlt(0, 5, 2) = 0.75;
440 hexNodesAlt(0, 6, 0) = 1.00; hexNodesAlt(0, 6, 1) = 1.00, hexNodesAlt(0, 6, 2) = 0.50;
441 hexNodesAlt(0, 7, 0) = 0.00; hexNodesAlt(0, 7, 1) = 1.00, hexNodesAlt(0, 7, 2) = 0.75;
442
443 // Step 2.c: Specify the edge ordinal, relative to the reference cell, of the edge workset
444 subcellDim = 1;
445 subcellOrd = 5;
446
447 /*
448 * 3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1]
449 */
450
451 // Step 3.a: selects Gauss rule of order 4 on [-1,1]
452 Teuchos::RCP<Cubature<double> > hexEdgeCubature = cubFactory.create(paramEdge, 4);
453
454 // Step 3.b allocate storage for cubature points
455 cubDim = hexEdgeCubature -> getDimension();
456 numCubPoints = hexEdgeCubature -> getNumPoints();
457
458 // Arrays must be properly sized for the specified set of Gauss points
459 paramEdgePoints.resize(numCubPoints, cubDim);
460 paramEdgeWeights.resize(numCubPoints);
461 hexEdgeCubature -> getCubature(paramEdgePoints, paramEdgeWeights);
462
463 /*
464 * 4. Map Gauss points from [-1,1] to every edge in the workset
465 */
466
467 // Step 4.a: Allocate storage for integration points on the reference edge
468 refEdgePoints.resize(numCubPoints, pCellDim);
469
470 // Step 4.b: Allocate storage for integration points on the edge workset
471 worksetEdgePoints.resize(worksetSize, numCubPoints, pCellDim);
472
473 // Step 4.c: Map Gauss points to reference edge: paramEdgePoints -> refEdgePoints
475 paramEdgePoints,
476 subcellDim,
477 subcellOrd,
478 hexahedron_8);
479
480 // Step 4.d: Map Gauss points from ref. edge to edge workset: refEdgePoints -> worksetEdgePoints
481 CellTools::mapToPhysicalFrame(worksetEdgePoints,
482 refEdgePoints,
483 hexNodes,
484 hexahedron_8);
485
486 /*
487 * 5. Print Gauss points on the edges of the workset
488 */
489 for(int pCell = 0; pCell < worksetSize; pCell++){
490
491 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
492
493 for(int pt = 0; pt < numCubPoints; pt++){
494 std::cout << "\t 1D Gauss point "
495 << std::setw(12) << std::right << paramEdgePoints(pt, 0) << std::setw(10) << " --> " << "("
496 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 0) << ", "
497 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 1) << ", "
498 << std::setw(10) << std::right << worksetEdgePoints(pCell, pt, 2) << ")\n";
499 }
500 std::cout << "\n";
501 }//pCell
502
503
504
505 std::cout << "\n" \
506 << "===============================================================================\n"\
507 << "| EXAMPLE 6: Edge tangents at Gauss points on a Hexahedron edge workset |\n"\
508 << "===============================================================================\n";
509
510 /* This task requires Gauss points on edge parametrization domain [-1,1] and on the
511 * reference edge whose ordinal matches the edge workset ordinal. This repeats the first few
512 * steps from Example 5:
513 *
514 * 1. Define cubature factory and topology for edge parametrization domain [-1,1];
515 * (Can reuse Step 1, Example 1 because edge parametrization domain is always [-1,1])
516 * 2. Define an edge workset;
517 * 3. Obtain the desired Gauss rule on the edge parametrization domain [-1,1];
518 * 4. Map Gauss points from [-1,1] to reference edge
519 * NOTE: this example only demonstrates computation of the edge tangents and so, Gauss
520 * points on the edge workset are not needed. Thus we skip mapping Gauss points from
521 * reference edge to edge workset
522 *
523 * 5. Compute Jacobians at Gauss points on reference edge for all cells in the workset:
524 */
525
526 // Step 5.a: Define and allocate storage for workset Jacobians
527 worksetJacobians.resize(worksetSize, numCubPoints, pCellDim, pCellDim);
528
529 // Step 5.b: Compute Jacobians at Gauss pts. on reference edge for all parent cells:
530 CellTools::setJacobian(worksetJacobians,
531 refEdgePoints,
532 hexNodes,
533 hexahedron_8);
534
535 /*
536 * 6. Get the (non-normalized) edge tangents for the edge workset:
537 */
538 // Step 6.a: Allocate storage for edge tangents
539 edgeTangents.resize(worksetSize, numCubPoints, pCellDim);
540
541 // Step 6.b: Compute the edge tangents:
543 worksetJacobians,
544 subcellOrd,
545 hexahedron_8);
546
547 // Step 6.c: Print edge tangents at Gauss points on workset edges (these Gauss points were computed in Example 5)
548 std::cout
549 << "Edge tangents computed by CellTools::getPhysicalEdgeTangents.\n"
550 << "Local edge ordinal = " << subcellOrd <<"\n";
551
552 for(int pCell = 0; pCell < worksetSize; pCell++){
553
554 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
555
556 for(int pt = 0; pt < numCubPoints; pt++){
557 std::cout << "\t 3D Gauss point: ("
558 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 0) << ", "
559 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 1) << ", "
560 << std::setw(8) << std::right << worksetEdgePoints(pCell, pt, 2) << ") "
561 << std::setw(8) << " edge tangent: " << "("
562 << std::setw(8) << std::right << edgeTangents(pCell, pt, 0) << ", "
563 << std::setw(8) << std::right << edgeTangents(pCell, pt, 1) << ", "
564 << std::setw(8) << std::right << edgeTangents(pCell, pt, 2) << ")\n";
565 }
566 std::cout << "\n";
567 }//pCell
568
569
570 std::cout << "\n" \
571 << "===============================================================================\n"\
572 << "| EXAMPLE 7: Definition of integration points on a Hexahedron face workset |\n"\
573 << "===============================================================================\n";
574 /*
575 * 1. Common tasks for getting integration points on face worksets:
576 * 1.a: can reuse the cubature factory, but need parametrization domain for the faces
577 */
578
579 // Step 1.b: Define topology of the face parametrization domain as [-1,1]x[-1,1]
580 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
581
582 // Step 1.c: Define storage for cubature points and weights on [-1,1]x[-1,1]
583 FieldContainer<double> paramFaceWeights;
584 FieldContainer<double> paramFacePoints;
585
586 // Step 1.d: Define storage for cubature points on a reference face
587 FieldContainer<double> refFacePoints;
588
589 // Step 1.f: Define storage for cubature points on workset faces
590 FieldContainer<double> worksetFacePoints;
591
592 /*
593 * 2. Define face workset comprising of 1 face corresponding to reference face 5 on Hexahedron<8>
594 * 2.a: Reuse the parent cell topology from Example 3
595 * 2.b: Reuse the vertices from Example 3
596 */
597
598 // Step 2.c: Specify dimension and ordinal (relative to ref. cell) of the subcells in the workset
599 subcellDim = 2;
600 subcellOrd = 5;
601
602 /*
603 * 3. Obtain the desired Gauss rule on the face parametrization domain [-1,1]x[-1,1]
604 */
605
606 // Step 3.a: selects Gauss rule of order 3 on [-1,1]x[-1,1]
607 Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3);
608
609 // Step 3.b allocate storage for cubature points on [-1,1]x[-1,1]
610 cubDim = hexFaceCubature -> getDimension();
611 numCubPoints = hexFaceCubature -> getNumPoints();
612
613 // Arrays must be properly sized for the specified set of Gauss points
614 paramFacePoints.resize(numCubPoints, cubDim);
615 paramFaceWeights.resize(numCubPoints);
616 hexFaceCubature -> getCubature(paramFacePoints, paramFaceWeights);
617
618 /*
619 * 4. Map Gauss points from [-1,1]x[-1,1] to every face in the workset
620 */
621
622 // Step 4.a: Allocate storage for integration points on the reference face
623 refFacePoints.resize(numCubPoints, pCellDim);
624
625 // Step 4.b: Allocate storage for integration points on the face in the workset
626 worksetFacePoints.resize(worksetSize, numCubPoints, pCellDim);
627
628 // Step 4.c: Map Gauss points to reference face: paramFacePoints -> refFacePoints
630 paramFacePoints,
631 subcellDim,
632 subcellOrd,
633 hexahedron_8);
634
635 // Step 4.d: Map Gauss points from ref. face to face workset: refFacePoints -> worksetFacePoints
636 CellTools::mapToPhysicalFrame(worksetFacePoints,
637 refFacePoints,
638 hexNodesAlt,
639 hexahedron_8);
640
641
642 /*
643 * 5. Print Gauss points on the faces of the workset
644 */
645 for(int pCell = 0; pCell < worksetSize; pCell++){
646
647 CellTools::printWorksetSubcell(hexNodesAlt, hexahedron_8, pCell, subcellDim, subcellOrd);
648
649 for(int pt = 0; pt < numCubPoints; pt++){
650 std::cout << "\t 2D Gauss point ("
651 << std::setw(8) << std::right << paramFacePoints(pt, 0) << ", "
652 << std::setw(8) << std::right << paramFacePoints(pt, 1) << ") "
653 << std::setw(8) << " --> " << "("
654 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", "
655 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", "
656 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ")\n";
657 }
658 std::cout << "\n";
659 }//pCell
660
661
662 std::cout << "\n" \
663 << "===============================================================================\n"\
664 << "| EXAMPLE 8: Face normals at Gauss points on a Hexahedron face workset |\n"\
665 << "===============================================================================\n";
666
667 /* This task requires Gauss points on face parametrization domain [-1,1]x[-1,1] and on the
668 * reference face whose ordinal matches the face workset ordinal. This repeats the first few
669 * steps from Example 5:
670 *
671 * 1. Define cubature factory and topology for face parametrization domain [-1,1]x[-1,1];
672 * 2. Define a face workset;
673 * 3. Obtain the desired Gauss rule on the face parametrization domain [-1,1]x[-1,1];
674 * 4. Map Gauss points from [-1,1]x[-1,1] to reference face
675 * NOTE: this example only demonstrates computation of the face normals and so, Gaus
676 * points on the face workset are not needed. Thus we skip mapping Gauss points from
677 * reference face to face workset
678 *
679 * 5. Compute Jacobians at Gauss points on reference face for all cells in the workset:
680 */
681
682 // Step 5.a: Define and allocate storage for workset Jacobians (reuse FC from example 5)
683 worksetJacobians.resize(worksetSize, numCubPoints, pCellDim, pCellDim);
684
685 // Step 5.b: Compute Jacobians at Gauss pts. on reference face for all parent cells:
686 CellTools::setJacobian(worksetJacobians,
687 refFacePoints,
688 hexNodesAlt,
689 hexahedron_8);
690
691
692 /*
693 * 6. Get the (non-normalized) face normals for the face workset directly
694 */
695 // Step 6.a: Allocate storage for face normals
696 FieldContainer<double> faceNormals(worksetSize, numCubPoints, pCellDim);
697
698 // Step 6.b: Compute the face normals
700 worksetJacobians,
701 subcellOrd,
702 hexahedron_8);
703
704 // Step 6.c: Print face normals at Gauss points on workset faces (these Gauss points were computed in Example 5)
705 std::cout
706 << "Face normals computed by CellTools::getPhysicalFaceNormals\n"
707 << "Local face ordinal = " << subcellOrd <<"\n";
708
709 for(int pCell = 0; pCell < worksetSize; pCell++){
710
711 CellTools::printWorksetSubcell(hexNodesAlt, hexahedron_8, pCell, subcellDim, subcellOrd);
712
713 for(int pt = 0; pt < numCubPoints; pt++){
714 std::cout << "\t 3D Gauss point: ("
715 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", "
716 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", "
717 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ") "
718 << std::setw(8) << " outer normal: " << "("
719 << std::setw(8) << std::right << faceNormals(pCell, pt, 0) << ", "
720 << std::setw(8) << std::right << faceNormals(pCell, pt, 1) << ", "
721 << std::setw(8) << std::right << faceNormals(pCell, pt, 2) << ")\n";
722 }
723 std::cout << "\n";
724 }//pCell
725
726 /*
727 * 7. Get the (non-normalized) face normals for the face workset using the face tangents. This may
728 * be useful if, for whatever reason, face tangents are needed independently
729 */
730 // Step 7.a: Allocate storage for face tangents
731 FieldContainer<double> uFaceTan(worksetSize, numCubPoints, pCellDim);
732 FieldContainer<double> vFaceTan(worksetSize, numCubPoints, pCellDim);
733
734 // Step 7.b: Compute face tangents
736 vFaceTan,
737 worksetJacobians,
738 subcellOrd,
739 hexahedron_8);
740
741 // Step 7.c: Face outer normals (relative to parent cell) are uTan x vTan:
742 RealSpaceTools<double>::vecprod(faceNormals, uFaceTan, vFaceTan);
743
744 // Step 7.d: Print face normals at Gauss points on workset faces (these points were computed in Example 7)
745 std::cout << "Face normals computed by CellTools::getPhysicalFaceTangents followed by vecprod\n";
746 for(int pCell = 0; pCell < worksetSize; pCell++){
747
748 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
749
750 for(int pt = 0; pt < numCubPoints; pt++){
751 std::cout << "\t 3D Gauss point: ("
752 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", "
753 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", "
754 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ") "
755 << std::setw(8) << " outer normal: " << "("
756 << std::setw(8) << std::right << faceNormals(pCell, pt, 0) << ", "
757 << std::setw(8) << std::right << faceNormals(pCell, pt, 1) << ", "
758 << std::setw(8) << std::right << faceNormals(pCell, pt, 2) << ")\n";
759 }
760 std::cout << "\n";
761 }//pCell
762
763
764 std::cout << "\n" \
765 << "===============================================================================\n"\
766 << "| EXAMPLE 8: Side normals at Gauss points on a Hexahedron side workset |\n"\
767 << "===============================================================================\n";
768
769 /* For this task we reuse the edge workset from Example 7 as a side workset and compute
770 * normals to the sides in that set. This task repeats the first 5 steps from Example 3
771 * The only difference is that in Step 6 we call CellTools::getPhysicalSideNormals
772 */
773
774 /*
775 * 6. Get the (non-normalized) side normals for the side (face) workset:
776 */
777 // Step 6.a: Allocate storage for side normals
778 sideNormals.resize(worksetSize, numCubPoints, pCellDim);
779
780 // Step 6.b: Compute the side normals:
782 worksetJacobians,
783 subcellOrd,
784 hexahedron_8);
785
786 // Step 7.d: Print side normals at Gauss points on workset sides (these points were computed in Example 7)
787 std::cout << "Side normals computed by CellTools::getPhysicalSideNormals\n";
788 for(int pCell = 0; pCell < worksetSize; pCell++){
789
790 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
791
792 for(int pt = 0; pt < numCubPoints; pt++){
793 std::cout << "\t 3D Gauss point: ("
794 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 0) << ", "
795 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 1) << ", "
796 << std::setw(8) << std::right << worksetFacePoints(pCell, pt, 2) << ") "
797 << std::setw(8) << " side normal: " << "("
798 << std::setw(8) << std::right << sideNormals(pCell, pt, 0) << ", "
799 << std::setw(8) << std::right << sideNormals(pCell, pt, 1) << ", "
800 << std::setw(8) << std::right << sideNormals(pCell, pt, 2) << ")\n";
801 }
802 std::cout << "\n";
803 }//pCell
804
805 return 0;
806}
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 mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
static void getPhysicalSideNormals(ArraySideNormal &sideNormals, const ArrayJac &worksetJacobians, const int &worksetSideOrd, const shards::CellTopology &parentCell)
Computes non-normalized normal vectors to physical sides in a side workset .
static void getPhysicalFaceTangents(ArrayFaceTangentU &faceTanU, ArrayFaceTangentV &faceTanV, const ArrayJac &worksetJacobians, const int &worksetFaceOrd, const shards::CellTopology &parentCell)
Computes non-normalized tangent vector pairs to physical faces in a face workset ; (see Subcell works...
static void printWorksetSubcell(const ArrayCell &cellWorkset, const shards::CellTopology &parentCell, const int &pCellOrd, const int &subcellDim, const int &subcellOrd, const int &fieldWidth=3)
Prints the nodes of a subcell from a cell workset.
static void 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 vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight