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 
49 #include "Intrepid_CellTools.hpp"
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Intrepid_DefaultCubatureFactory.hpp"
52 #include "Shards_CellTopology.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 
55 using namespace std;
56 using namespace Intrepid;
57 
58 int 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
172  CellTools::mapToReferenceSubcell(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
265  CellTools::mapToReferenceSubcell(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:
331  CellTools::getPhysicalEdgeTangents(edgeTangents,
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:
372  CellTools::getPhysicalSideNormals(sideNormals,
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
474  CellTools::mapToReferenceSubcell(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:
542  CellTools::getPhysicalEdgeTangents(edgeTangents,
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
629  CellTools::mapToReferenceSubcell(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
699  CellTools::getPhysicalFaceNormals(faceNormals,
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
735  CellTools::getPhysicalFaceTangents(uFaceTan,
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:
781  CellTools::getPhysicalSideNormals(sideNormals,
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 }
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:68
Intrepid_CellTools.hpp
Header file for the Intrepid::CellTools class.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Intrepid::FieldContainer::resize
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Definition: Intrepid_FieldContainerDef.hpp:1137
Intrepid::DefaultCubatureFactory
A factory class that generates specific instances of cubatures.
Definition: Intrepid_DefaultCubatureFactory.hpp:77
Intrepid::DefaultCubatureFactory::create
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Definition: Intrepid_DefaultCubatureFactoryDef.hpp:53
Intrepid::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:111
Generated by   doxygen 1.8.5