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 
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 
59 void vField(double& v1, double& v2, double& v3,
60  const double& x, const double& y, const double& z);
61 
62 using namespace std;
63 using namespace Intrepid;
64 
65 int 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
200  CellTools::mapToReferenceSubcell(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
240  CellTools::getPhysicalFaceTangents(worksetFaceTu,
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 
371 void 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 
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
vField
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.
Definition: example_04.cpp:371
Intrepid::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:111
Generated by   doxygen 1.8.5