Intrepid
example_04.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
59void vField(double& v1, double& v2, double& v3,
60 const double& x, const double& y, const double& z);
61
62using namespace std;
63using namespace Intrepid;
64
65int main(int argc, char *argv[]) {
66
67 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
68
69 typedef CellTools<double> CellTools;
70 typedef RealSpaceTools<double> RealSpaceTools;
71 typedef shards::CellTopology CellTopology;
72
73 std::cout \
74 << "===============================================================================\n" \
75 << "| |\n" \
76 << "| Example use of the CellTools class |\n" \
77 << "| |\n" \
78 << "| 1) Computation of face flux, for a given vector field, on a face workset |\n" \
79 << "| 2) Computation of edge circulation, for a given vector field, on a face |\n" \
80 << "| workset. |\n" \
81 << "| |\n" \
82 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \
83 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \
84 << "| Kara Peterson (kjpeter@sandia.gov) |\n" \
85 << "| |\n" \
86 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
87 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
88 << "| |\n" \
89 << "===============================================================================\n"\
90 << "| EXAMPLE 1: Computation of face flux on a face workset |\n"\
91 << "===============================================================================\n";
92
93
108 /*************************************************************************************************
109 *
110 * Step 0. Face workset comprising of 1 face of a Hexahedron<8> cell
111 *
112 ************************************************************************************************/
113
114 // Step 0.a: Specify cell topology of the parent cell
115 CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() );
116
117 // Step 0.b: Specify the vertices of the parent Hexahedron<8> cell
118 int worksetSize = 2;
119 int pCellNodeCount = hexahedron_8.getVertexCount();
120 int pCellDim = hexahedron_8.getDimension();
121
122 FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim);
123 // cell 0 bottom face vertices:
124 hexNodes(0, 0, 0) = 0.00; hexNodes(0, 0, 1) = 0.00, hexNodes(0, 0, 2) = 0.00;
125 hexNodes(0, 1, 0) = 1.00; hexNodes(0, 1, 1) = 0.00, hexNodes(0, 1, 2) = 0.00;
126 hexNodes(0, 2, 0) = 1.00; hexNodes(0, 2, 1) = 1.00, hexNodes(0, 2, 2) = 0.00;
127 hexNodes(0, 3, 0) = 0.00; hexNodes(0, 3, 1) = 1.00, hexNodes(0, 3, 2) = 0.00;
128 // cell 0 top face vertices
129 hexNodes(0, 4, 0) = 0.00; hexNodes(0, 4, 1) = 0.00, hexNodes(0, 4, 2) = 1.00;
130 hexNodes(0, 5, 0) = 1.00; hexNodes(0, 5, 1) = 0.00, hexNodes(0, 5, 2) = 1.00;
131 hexNodes(0, 6, 0) = 1.00; hexNodes(0, 6, 1) = 1.00, hexNodes(0, 6, 2) = 1.00;
132 hexNodes(0, 7, 0) = 0.00; hexNodes(0, 7, 1) = 1.00, hexNodes(0, 7, 2) = 1.00;
133
134 // cell 1 bottom face vertices:
135 hexNodes(1, 0, 0) = 0.00; hexNodes(1, 0, 1) = 0.00, hexNodes(1, 0, 2) = 0.00;
136 hexNodes(1, 1, 0) = 1.00; hexNodes(1, 1, 1) = 0.00, hexNodes(1, 1, 2) = 0.00;
137 hexNodes(1, 2, 0) = 1.00; hexNodes(1, 2, 1) = 1.00, hexNodes(1, 2, 2) = 0.00;
138 hexNodes(1, 3, 0) = 0.00; hexNodes(1, 3, 1) = 1.00, hexNodes(1, 3, 2) = 0.00;
139 // cell 1 top face vertices
140 hexNodes(1, 4, 0) = 0.00; hexNodes(1, 4, 1) = 0.00, hexNodes(1, 4, 2) = 1.00;
141 hexNodes(1, 5, 0) = 1.00; hexNodes(1, 5, 1) = 0.00, hexNodes(1, 5, 2) = 1.00;
142 hexNodes(1, 6, 0) = 1.00; hexNodes(1, 6, 1) = 1.00, hexNodes(1, 6, 2) = 0.75;
143 hexNodes(1, 7, 0) = 0.00; hexNodes(1, 7, 1) = 1.00, hexNodes(1, 7, 2) = 1.00;
144
145
146 // Step 0.c: Specify the face ordinal, relative to the reference cell, of the face workset
147 int subcellDim = 2;
148 int subcellOrd = 1;
149
150
151
152 /*************************************************************************************************
153 *
154 * Step 1: Obtain Gauss points on workset faces for Hexahedron<8> topology
155 * 1.1 Define cubature factory, face parametrization domain and arrays for cubature points
156 * 1.2 Select Gauss rule on D = [-1,1]x[-1,1]
157 * 1.3 Map Gauss points from D to reference face and workset faces
158 *
159 ************************************************************************************************/
160
161 // Step 1.1.a: Define CubatureFactory
162 DefaultCubatureFactory<double> cubFactory;
163
164 // Step 1.1.b: Define topology of the face parametrization domain as [-1,1]x[-1,1]
165 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
166
167 // Step 1.1.c: Define storage for cubature points and weights on [-1,1]x[-1,1]
168 FieldContainer<double> paramGaussWeights;
169 FieldContainer<double> paramGaussPoints;
170
171 // Step 1.1.d: Define storage for cubature points on a reference face
172 FieldContainer<double> refGaussPoints;
173
174 // Step 1.1.f: Define storage for cubature points on workset faces
175 FieldContainer<double> worksetGaussPoints;
176
177 //----------------
178
179 // Step 1.2.a: selects Gauss rule of order 3 on [-1,1]x[-1,1]
180 Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3);
181
182 // Step 1.2.b allocate storage for cubature points on [-1,1]x[-1,1]
183 int cubDim = hexFaceCubature -> getDimension();
184 int numCubPoints = hexFaceCubature -> getNumPoints();
185
186 // Arrays must be properly sized for the specified set of Gauss points
187 paramGaussPoints.resize(numCubPoints, cubDim);
188 paramGaussWeights.resize(numCubPoints);
189 hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights);
190
191 //----------------
192
193 // Step 1.3.a: Allocate storage for Gauss points on the reference face
194 refGaussPoints.resize(numCubPoints, pCellDim);
195
196 // Step 1.3.b: Allocate storage for Gauss points on the face in the workset
197 worksetGaussPoints.resize(worksetSize, numCubPoints, pCellDim);
198
199 // Step 1.3.c: Map Gauss points to reference face: paramGaussPoints -> refGaussPoints
201 paramGaussPoints,
202 subcellDim,
203 subcellOrd,
204 hexahedron_8);
205
206 // Step 1.3.d: Map Gauss points from ref. face to face workset: refGaussPoints -> worksetGaussPoints
207 CellTools::mapToPhysicalFrame(worksetGaussPoints,
208 refGaussPoints,
209 hexNodes,
210 hexahedron_8);
211
212
213
214
215 /*************************************************************************************************
216 *
217 * Step 2. Obtain (non-normalized) face normals at cubature points on workset faces
218 * 2.1 Compute parent cell Jacobians at Gauss points on workset faces
219 * 2.2 Compute face tangents on workset faces and their vector product
220 *
221 ************************************************************************************************/
222
223 // Step 2.1.a: Define and allocate storage for workset Jacobians
224 FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim);
225
226 // Step 2.1.b: Compute Jacobians at Gauss pts. on reference face for all parent cells:
227 CellTools::setJacobian(worksetJacobians,
228 refGaussPoints,
229 hexNodes,
230 hexahedron_8);
231
232 //----------------
233
234 // Step 2.2.a: Allocate storage for face tangents and face normals
235 FieldContainer<double> worksetFaceTu(worksetSize, numCubPoints, pCellDim);
236 FieldContainer<double> worksetFaceTv(worksetSize, numCubPoints, pCellDim);
237 FieldContainer<double> worksetFaceN(worksetSize, numCubPoints, pCellDim);
238
239 // Step 2.2.b: Compute face tangents
241 worksetFaceTv,
242 worksetJacobians,
243 subcellOrd,
244 hexahedron_8);
245
246 // Step 2.2.c: Face outer normals (relative to parent cell) are uTan x vTan:
247 RealSpaceTools::vecprod(worksetFaceN, worksetFaceTu, worksetFaceTv);
248
249
250
251 /*************************************************************************************************
252 *
253 * Step 3. Evaluate the vector field u(x,y,z) at cubature points on workset faces
254 *
255 ************************************************************************************************/
256
257 // Step 3.a: Allocate storage for vector field values at Gauss points on workset faces
258 FieldContainer<double> worksetVFieldVals(worksetSize, numCubPoints, pCellDim);
259
260 // Step 3.b: Compute vector field at Gauss points: here we take u(x,y,z) = (x,y,z)
261 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
262 for(int ptOrd = 0; ptOrd < numCubPoints; ptOrd++){
263
264 double x = worksetGaussPoints(pCellOrd, ptOrd, 0);
265 double y = worksetGaussPoints(pCellOrd, ptOrd, 1);
266 double z = worksetGaussPoints(pCellOrd, ptOrd, 2);
267
268 vField(worksetVFieldVals(pCellOrd, ptOrd, 0),
269 worksetVFieldVals(pCellOrd, ptOrd, 1),
270 worksetVFieldVals(pCellOrd, ptOrd, 2), x, y, z);
271
272 }// pt
273 }//cell
274
275
276 /*************************************************************************************************
277 *
278 * Step 4. Compute dot product of u(x,y,z) and the face normals, times the cubature weights
279 *
280 ************************************************************************************************/
281
282 // Allocate storage for dot product of vector field and face normals at Gauss points
283 FieldContainer<double> worksetFieldDotNormal(worksetSize, numCubPoints);
284
285 // Compute the dot product
286 RealSpaceTools::dot(worksetFieldDotNormal, worksetVFieldVals, worksetFaceN);
287
288 // Allocate storage for face fluxes on the workset
289 FieldContainer<double> worksetFluxes(worksetSize);
290
291 //----------------
292
293 // Integration loop (temporary)
294 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
295 worksetFluxes(pCellOrd) = 0.0;
296 for(int pt = 0; pt < numCubPoints; pt++ ){
297 worksetFluxes(pCellOrd) += worksetFieldDotNormal(pCellOrd, pt)*paramGaussWeights(pt);
298 }// pt
299 }//cell
300
301 std::cout << "Face fluxes on workset faces : \n\n";
302 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){
303
304 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCellOrd, subcellDim, subcellOrd);
305 std::cout << " Flux = " << worksetFluxes(pCellOrd) << "\n\n";
306
307 }
308
309
310
311 /*************************************************************************************************
312 *
313 * Optional: print Gauss points and face normals at Gauss points
314 *
315 ************************************************************************************************/
316
317 // Print Gauss points on [-1,1]x[-1,1] and their images on workset faces
318 std::cout \
319 << "===============================================================================\n" \
320 << "| Gauss points on workset faces: |\n" \
321 << "===============================================================================\n";
322 for(int pCell = 0; pCell < worksetSize; pCell++){
323
324 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
325
326 for(int pt = 0; pt < numCubPoints; pt++){
327 std::cout << "\t 2D Gauss point ("
328 << std::setw(8) << std::right << paramGaussPoints(pt, 0) << ", "
329 << std::setw(8) << std::right << paramGaussPoints(pt, 1) << ") "
330 << std::setw(8) << " --> " << "("
331 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", "
332 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", "
333 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")\n";
334 }
335 std::cout << "\n\n";
336 }//pCell
337
338
339 // Print face normals at Gauss points on workset faces
340 std::cout \
341 << "===============================================================================\n" \
342 << "| Face normals (non-unit) at Gauss points on workset faces: |\n" \
343 << "===============================================================================\n";
344 for(int pCell = 0; pCell < worksetSize; pCell++){
345
346 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd);
347
348 for(int pt = 0; pt < numCubPoints; pt++){
349 std::cout << "\t 3D Gauss point: ("
350 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", "
351 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", "
352 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")"
353 << std::setw(8) << " out. normal: " << "("
354 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 0) << ", "
355 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 1) << ", "
356 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 2) << ")\n";
357 }
358 std::cout << "\n";
359 }//pCell
360
361 return 0;
362}
363
364/*************************************************************************************************
365 *
366 * Definition of the vector field function
367 *
368 ************************************************************************************************/
369
370
371void vField(double& v1, double& v2, double& v3, const double& x, const double& y, const double& z)
372{
373 v1 = x;
374 v2 = y;
375 v3 = z;
376}
377
378
int main(int argc, char *argv[])
void vField(double &v1, double &v2, double &v3, const double &x, const double &y, const double &z)
Evaluation of a 3D vector field in physical coordinate frame.
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 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 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.
Implementation of basic linear algebra functionality in Euclidean space.
static Scalar dot(const Scalar *inArray1, const Scalar *inArray2, const int size)
Computes dot product of contiguous data inArray1 and inArray2 of size size.
static void vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight