Intrepid
example_02.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
53#include "Teuchos_GlobalMPISession.hpp"
54
55#include "Shards_CellTopology.hpp"
56
57#include "Teuchos_RCP.hpp"
58
59using namespace std;
60using namespace Intrepid;
61using namespace shards;
62
63int main(int argc, char *argv[]) {
64
65 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
66
67 typedef CellTools<double> CellTools;
68 typedef shards::CellTopology CellTopology;
69
70 cout \
71 << "===============================================================================\n" \
72 << "| |\n" \
73 << "| Example use of the CellTools class |\n" \
74 << "| |\n" \
75 << "| 1) Reference edge parametrizations |\n" \
76 << "| 2) Reference face parametrizations |\n" \
77 << "| |\n" \
78 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
79 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \
80 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
81 << "| |\n" \
82 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
83 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
84 << "| |\n" \
85 << "===============================================================================\n"\
86 << "| Summary: |\n"\
87 << "| Reference edge parametrizations map [-1,1] to the edges of reference cells. |\n"\
88 << "| They are used to define, e.g., integration points on the edges of 2D and 3D |\n"\
89 << "| reference cells. Edge parametrizations for special 2D cells such as Beam |\n"\
90 << "| and ShellLine, are also supported. |\n"\
91 << "===============================================================================\n";
92
93 /*
94 Specification of integration points on 1-subcells (edges) of reference cells. Edges are
95 parametrized by [-1,1] and integration points on an edge are defined by mapping integration
96 points from the parametrization domain [-1,1] to a specific edge on the reference cell.
97
98 1. Common tasks: definition of integration points in the edge parametrization domain [-1,1]
99 These steps are independent of parent cell topology:
100
101 a. Instantiate a CubatureFactory object to create cubatures (needed for face maps too)
102 b. Define parametrization domain for the edges as having Line<2> cell topology. This is
103 required by the CubatureFactory in order to select cubature points and weights from
104 the reference line [-1,1]
105 c. Use CubatureFactory to select cubature of the desired degree for the Line<2> topology
106 d. Allocate containers for the cubature points and weights.
107
108 2. Parent cell topology specific tasks
109
110 a. Select the parent cell topology
111 b. Allocate containers for the images of the integration points on [-1,1] on the edges
112 c. Apply the edge parametrization map to the pointss in [-1,1]
113 */
114
115 // Step 1.a (Define CubatureFactory)
116 DefaultCubatureFactory<double> cubFactory;
117
118
119 // Step 1.b (Define the topology of the parametrization domain)
120 CellTopology edgeParam(shards::getCellTopologyData<shards::Line<2> >() );
121
122
123 // Step 1.c (selects Gauss rule of order 6 on [-1,1])
124 Teuchos::RCP<Cubature<double> > edgeParamCubature = cubFactory.create(edgeParam, 6);
125
126
127 // Step 1.d (allocate storage for cubature points)
128 int cubDim = edgeParamCubature -> getDimension();
129 int numCubPoints = edgeParamCubature -> getNumPoints();
130
131 FieldContainer<double> edgeParamCubPts(numCubPoints, cubDim);
132 FieldContainer<double> edgeParamCubWts(numCubPoints);
133 edgeParamCubature -> getCubature(edgeParamCubPts, edgeParamCubWts);
134
135
136
137 std::cout \
138 << "===============================================================================\n"\
139 << "| EXAMPLE 1.1 |\n"
140 << "| Edge parametrizations for standard 2D cells: Triangle |\n"\
141 << "===============================================================================\n";
142
143 // Step 2.a (select reference cell topology)
144 CellTopology triangle_3(getCellTopologyData<Triangle<3> >() );
145
146
147 // Step 2.b (allocate storage for points on an edge of the reference cell)
148 FieldContainer<double> triEdgePoints(numCubPoints, triangle_3.getDimension() );
149
150
151 // Step 2.c (same points are mapped to all edges, can also map different set to each edge)
152 for(int edgeOrd = 0; edgeOrd < (int)triangle_3.getEdgeCount(); edgeOrd++){
153
155 edgeParamCubPts,
156 1,
157 edgeOrd,
158 triangle_3);
159
160 // Optional: print the vertices of the reference edge
161 CellTools::printSubcellVertices(1, edgeOrd, triangle_3);
162
163 for(int pt = 0; pt < numCubPoints; pt++){
164 std::cout << "\t Parameter point "
165 << std::setw(12) << std::right << edgeParamCubPts(pt, 0) << std::setw(10) << " --> " << "("
166 << std::setw(10) << std::right << triEdgePoints(pt, 0) << " , "
167 << std::setw(10) << std::right << triEdgePoints(pt, 1) << ")\n";
168 }
169 std::cout << "\n";
170 }
171
172
173
174 std::cout \
175 << "===============================================================================\n"\
176 << "| EXAMPLE 1.2 |\n"
177 << "| Edge parametrizations for standard 2D cells: Quadrilateral |\n"\
178 << "===============================================================================\n";
179
180 // Step 2.a (select reference cell topology)
181 CellTopology quad_4(getCellTopologyData<Quadrilateral<4> >() );
182
183
184 // Step 2.b (allocate storage for points on an edge of the reference cell)
185 FieldContainer<double> quadEdgePoints(numCubPoints, quad_4.getDimension() );
186
187
188 // Step 2.c (same points are mapped to all edges, can also map different set to each edge)
189 for(int edgeOrd = 0; edgeOrd < (int)quad_4.getEdgeCount(); edgeOrd++){
190
191 CellTools::mapToReferenceSubcell(quadEdgePoints,
192 edgeParamCubPts,
193 1,
194 edgeOrd,
195 quad_4);
196
197 // Optional: print the vertices of the reference edge
198 CellTools::printSubcellVertices(1, edgeOrd, quad_4);
199
200 for(int pt = 0; pt < numCubPoints; pt++){
201 std::cout << "\t Parameter point "
202 << std::setw(12) << std::right << edgeParamCubPts(pt, 0) << std::setw(10) << " --> " << "("
203 << std::setw(10) << std::right << quadEdgePoints(pt, 0) << " , "
204 << std::setw(10) << std::right << quadEdgePoints(pt, 1) << ")\n";
205 }
206 std::cout << "\n";
207 }
208
209
210 /*
211 Specification of integration points on 2-subcells (faces) of reference cells. Reference cells
212 can have triangular, quadrilateral or a mixture of triangular and quadrilateral faces. Thus,
213 parametrization domain of a face depends on that face's topology and is either the standard
214 2-simplex {(0,0), (1,0), (0,1)} for triangular faces or the standard 2-cube [-1,1]^2 for
215 quadrilateral faces.
216
217 1. Common tasks: definition of integration points in the standard 2-simplex and the standard
218 2-cube. These steps are independent of parent cell topology:
219
220 a. Instantiate a CubatureFactory object to create cubatures (already done!)
221 b. Define parametrization domain for traingular faces as having Triangle<3> topology and for
222 quadrilateral faces - as having Quadrilateral<4> topology. This is required by the
223 CubatureFactory in order to select cubature points and weights from the appropriate
224 face parametrization domain.
225 c. Use CubatureFactory to select cubature of the desired degree for Triangle<3> and
226 Quadrilateral<4> topologies
227 d. Allocate containers for the cubature points and weights on the parametrization domains.
228
229 2. Parent cell topology specific tasks
230
231 a. Select the parent cell topology
232 b. Allocate containers for the images of the integration points from the parametrization
233 domains on the reference faces
234 c. Apply the face parametrization map to the points in the parametrization domain
235 */
236
237 // Step 1.b (Define the topology of the parametrization domain)
238 CellTopology triFaceParam(shards::getCellTopologyData<shards::Triangle<3> >() );
239 CellTopology quadFaceParam(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
240
241
242 // Step 1.c (selects Gauss rule of order 3 on [-1,1]^2 and a rule of order 3 on Triangle)
243 Teuchos::RCP<Cubature<double> > triFaceParamCubature = cubFactory.create(triFaceParam, 3);
244 Teuchos::RCP<Cubature<double> > quadFaceParamCubature = cubFactory.create(quadFaceParam, 3);
245
246
247 // Step 1.d - Triangle faces (allocate storage for cubature points)
248 int triFaceCubDim = triFaceParamCubature -> getDimension();
249 int triFaceNumCubPts = triFaceParamCubature -> getNumPoints();
250
251 FieldContainer<double> triFaceParamCubPts(triFaceNumCubPts, triFaceCubDim);
252 FieldContainer<double> triFaceParamCubWts(triFaceNumCubPts);
253 triFaceParamCubature -> getCubature(triFaceParamCubPts, triFaceParamCubWts);
254
255
256 // Step 1.d - Quadrilateral faces (allocate storage for cubature points)
257 int quadFaceCubDim = quadFaceParamCubature -> getDimension();
258 int quadFaceNumCubPts = quadFaceParamCubature -> getNumPoints();
259
260 FieldContainer<double> quadFaceParamCubPts(quadFaceNumCubPts, quadFaceCubDim);
261 FieldContainer<double> quadFaceParamCubWts(quadFaceNumCubPts);
262 quadFaceParamCubature -> getCubature(quadFaceParamCubPts, quadFaceParamCubWts);
263
264
265
266 std::cout \
267 << "===============================================================================\n"\
268 << "| EXAMPLE 2.1 |\n"
269 << "| Face parametrizations for standard 3D cells: Tetrahedron |\n"\
270 << "===============================================================================\n";
271
272 // Step 2.a (select reference cell topology)
273 CellTopology tet_4(getCellTopologyData<Tetrahedron<4> >() );
274
275
276 // Step 2.b (allocate storage for points on a face of the reference cell)
277 FieldContainer<double> tetFacePoints(triFaceNumCubPts, tet_4.getDimension() );
278
279
280 // Step 2.c (same points are mapped to all faces, can also map different set to each face)
281 for(int faceOrd = 0; faceOrd < (int)tet_4.getSideCount(); faceOrd++){
282
284 triFaceParamCubPts,
285 2,
286 faceOrd,
287 tet_4);
288
289 // Optional: print the vertices of the reference face
290 CellTools::printSubcellVertices(2, faceOrd, tet_4);
291
292 for(int pt = 0; pt < triFaceNumCubPts; pt++){
293 std::cout << "\t Parameter point ("
294 << std::setw(10) << std::right << triFaceParamCubPts(pt, 0) << " , "
295 << std::setw(10) << std::right << triFaceParamCubPts(pt, 1) << ") "
296 << std::setw(10) << " --> " << "("
297 << std::setw(10) << std::right << tetFacePoints(pt, 0) << " , "
298 << std::setw(10) << std::right << tetFacePoints(pt, 1) << " , "
299 << std::setw(10) << std::right << tetFacePoints(pt, 2) << ")\n";
300 }
301 std::cout << "\n";
302 }
303
304
305
306 std::cout \
307 << "===============================================================================\n"\
308 << "| EXAMPLE 2.2 |\n"
309 << "| Face parametrizations for standard 3D cells: Wedge |\n"\
310 << "| Example of a reference cell that has two different kinds of faces |\n"\
311 << "===============================================================================\n";
312
313
314 // Step 2.a (select reference cell topology)
315 CellTopology wedge_6(getCellTopologyData<Wedge<6> >() );
316
317
318 // Step 2.b (allocate storage for points on a face of the reference cell)
319 // Wedge<6> has Triangle<3> and Quadrilateral<4> faces. Two different arrays are needed
320 // to store the points on each face because different types integration rules are used
321 // on their respective parametrization domains and numbers of points defined by these
322 // rules do not necessarily match.
323 FieldContainer<double> wedgeTriFacePoints(triFaceNumCubPts, wedge_6.getDimension() );
324 FieldContainer<double> wedgeQuadFacePoints(quadFaceNumCubPts, wedge_6.getDimension() );
325
326
327 // Step 2.c (for Wedge<6> need to distinguish Triangle<3> and Quadrilateral<4> faces)
328 for(int faceOrd = 0; faceOrd < (int)wedge_6.getSideCount(); faceOrd++){
329
330 // Optional: print the vertices of the reference face
331 CellTools::printSubcellVertices(2, faceOrd, wedge_6);
332
333 if( wedge_6.getKey(2, faceOrd) == shards::Triangle<3>::key ){
334 CellTools::mapToReferenceSubcell(wedgeTriFacePoints,
335 triFaceParamCubPts,
336 2,
337 faceOrd,
338 wedge_6);
339
340 for(int pt = 0; pt < triFaceNumCubPts; pt++){
341 std::cout << "\t Parameter point ("
342 << std::setw(10) << std::right << triFaceParamCubPts(pt, 0) << " , "
343 << std::setw(10) << std::right << triFaceParamCubPts(pt, 1) << ") "
344 << std::setw(10) << " --> " << "("
345 << std::setw(10) << std::right << wedgeTriFacePoints(pt, 0) << " , "
346 << std::setw(10) << std::right << wedgeTriFacePoints(pt, 1) << " , "
347 << std::setw(10) << std::right << wedgeTriFacePoints(pt, 2) << ")\n";
348 }
349 std::cout << "\n";
350 }
351 else if(wedge_6.getKey(2, faceOrd) == shards::Quadrilateral<4>::key) {
352 CellTools::mapToReferenceSubcell(wedgeQuadFacePoints,
353 quadFaceParamCubPts,
354 2,
355 faceOrd,
356 wedge_6);
357
358 for(int pt = 0; pt < quadFaceNumCubPts; pt++){
359 std::cout << "\t Parameter point ("
360 << std::setw(10) << std::right << quadFaceParamCubPts(pt, 0) << " , "
361 << std::setw(10) << std::right << quadFaceParamCubPts(pt, 1) << ") "
362 << std::setw(10) << " --> " << "("
363 << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 0) << " , "
364 << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 1) << " , "
365 << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 2) << ")\n";
366 }
367 std::cout << "\n";
368 }
369 else {
370 std::cout << " Invalid face encountered \n";
371 }
372 }
373
374 return 0;
375}
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 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 printSubcellVertices(const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Prints the reference space coordinates of the vertices of the specified subcell.