|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00038 // Denis Ridzal (dridzal@sandia.gov), or 00039 // Kara Peterson (kjpeter@sandia.gov) 00040 // 00041 // ************************************************************************ 00042 // @HEADER 00043 00044 00049 #include "Intrepid_CellTools.hpp" 00050 #include "Intrepid_FieldContainer.hpp" 00051 #include "Intrepid_DefaultCubatureFactory.hpp" 00052 #include "Shards_CellTopology.hpp" 00053 #include "Teuchos_GlobalMPISession.hpp" 00054 00059 void vField(double& v1, double& v2, double& v3, 00060 const double& x, const double& y, const double& z); 00061 00062 using namespace std; 00063 using namespace Intrepid; 00064 00065 int main(int argc, char *argv[]) { 00066 00067 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00068 00069 typedef CellTools<double> CellTools; 00070 typedef RealSpaceTools<double> RealSpaceTools; 00071 typedef shards::CellTopology CellTopology; 00072 00073 std::cout \ 00074 << "===============================================================================\n" \ 00075 << "| |\n" \ 00076 << "| Example use of the CellTools class |\n" \ 00077 << "| |\n" \ 00078 << "| 1) Computation of face flux, for a given vector field, on a face workset |\n" \ 00079 << "| 2) Computation of edge circulation, for a given vector field, on a face |\n" \ 00080 << "| workset. |\n" \ 00081 << "| |\n" \ 00082 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \ 00083 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \ 00084 << "| Kara Peterson (kjpeter@sandia.gov) |\n" \ 00085 << "| |\n" \ 00086 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00087 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00088 << "| |\n" \ 00089 << "===============================================================================\n"\ 00090 << "| EXAMPLE 1: Computation of face flux on a face workset |\n"\ 00091 << "===============================================================================\n"; 00092 00093 00108 /************************************************************************************************* 00109 * 00110 * Step 0. Face workset comprising of 1 face of a Hexahedron<8> cell 00111 * 00112 ************************************************************************************************/ 00113 00114 // Step 0.a: Specify cell topology of the parent cell 00115 CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00116 00117 // Step 0.b: Specify the vertices of the parent Hexahedron<8> cell 00118 int worksetSize = 2; 00119 int pCellNodeCount = hexahedron_8.getVertexCount(); 00120 int pCellDim = hexahedron_8.getDimension(); 00121 00122 FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim); 00123 // cell 0 bottom face vertices: 00124 hexNodes(0, 0, 0) = 0.00; hexNodes(0, 0, 1) = 0.00, hexNodes(0, 0, 2) = 0.00; 00125 hexNodes(0, 1, 0) = 1.00; hexNodes(0, 1, 1) = 0.00, hexNodes(0, 1, 2) = 0.00; 00126 hexNodes(0, 2, 0) = 1.00; hexNodes(0, 2, 1) = 1.00, hexNodes(0, 2, 2) = 0.00; 00127 hexNodes(0, 3, 0) = 0.00; hexNodes(0, 3, 1) = 1.00, hexNodes(0, 3, 2) = 0.00; 00128 // cell 0 top face vertices 00129 hexNodes(0, 4, 0) = 0.00; hexNodes(0, 4, 1) = 0.00, hexNodes(0, 4, 2) = 1.00; 00130 hexNodes(0, 5, 0) = 1.00; hexNodes(0, 5, 1) = 0.00, hexNodes(0, 5, 2) = 1.00; 00131 hexNodes(0, 6, 0) = 1.00; hexNodes(0, 6, 1) = 1.00, hexNodes(0, 6, 2) = 1.00; 00132 hexNodes(0, 7, 0) = 0.00; hexNodes(0, 7, 1) = 1.00, hexNodes(0, 7, 2) = 1.00; 00133 00134 // cell 1 bottom face vertices: 00135 hexNodes(1, 0, 0) = 0.00; hexNodes(1, 0, 1) = 0.00, hexNodes(1, 0, 2) = 0.00; 00136 hexNodes(1, 1, 0) = 1.00; hexNodes(1, 1, 1) = 0.00, hexNodes(1, 1, 2) = 0.00; 00137 hexNodes(1, 2, 0) = 1.00; hexNodes(1, 2, 1) = 1.00, hexNodes(1, 2, 2) = 0.00; 00138 hexNodes(1, 3, 0) = 0.00; hexNodes(1, 3, 1) = 1.00, hexNodes(1, 3, 2) = 0.00; 00139 // cell 1 top face vertices 00140 hexNodes(1, 4, 0) = 0.00; hexNodes(1, 4, 1) = 0.00, hexNodes(1, 4, 2) = 1.00; 00141 hexNodes(1, 5, 0) = 1.00; hexNodes(1, 5, 1) = 0.00, hexNodes(1, 5, 2) = 1.00; 00142 hexNodes(1, 6, 0) = 1.00; hexNodes(1, 6, 1) = 1.00, hexNodes(1, 6, 2) = 0.75; 00143 hexNodes(1, 7, 0) = 0.00; hexNodes(1, 7, 1) = 1.00, hexNodes(1, 7, 2) = 1.00; 00144 00145 00146 // Step 0.c: Specify the face ordinal, relative to the reference cell, of the face workset 00147 int subcellDim = 2; 00148 int subcellOrd = 1; 00149 00150 00151 00152 /************************************************************************************************* 00153 * 00154 * Step 1: Obtain Gauss points on workset faces for Hexahedron<8> topology 00155 * 1.1 Define cubature factory, face parametrization domain and arrays for cubature points 00156 * 1.2 Select Gauss rule on D = [-1,1]x[-1,1] 00157 * 1.3 Map Gauss points from D to reference face and workset faces 00158 * 00159 ************************************************************************************************/ 00160 00161 // Step 1.1.a: Define CubatureFactory 00162 DefaultCubatureFactory<double> cubFactory; 00163 00164 // Step 1.1.b: Define topology of the face parametrization domain as [-1,1]x[-1,1] 00165 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00166 00167 // Step 1.1.c: Define storage for cubature points and weights on [-1,1]x[-1,1] 00168 FieldContainer<double> paramGaussWeights; 00169 FieldContainer<double> paramGaussPoints; 00170 00171 // Step 1.1.d: Define storage for cubature points on a reference face 00172 FieldContainer<double> refGaussPoints; 00173 00174 // Step 1.1.f: Define storage for cubature points on workset faces 00175 FieldContainer<double> worksetGaussPoints; 00176 00177 //---------------- 00178 00179 // Step 1.2.a: selects Gauss rule of order 3 on [-1,1]x[-1,1] 00180 Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3); 00181 00182 // Step 1.2.b allocate storage for cubature points on [-1,1]x[-1,1] 00183 int cubDim = hexFaceCubature -> getDimension(); 00184 int numCubPoints = hexFaceCubature -> getNumPoints(); 00185 00186 // Arrays must be properly sized for the specified set of Gauss points 00187 paramGaussPoints.resize(numCubPoints, cubDim); 00188 paramGaussWeights.resize(numCubPoints); 00189 hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights); 00190 00191 //---------------- 00192 00193 // Step 1.3.a: Allocate storage for Gauss points on the reference face 00194 refGaussPoints.resize(numCubPoints, pCellDim); 00195 00196 // Step 1.3.b: Allocate storage for Gauss points on the face in the workset 00197 worksetGaussPoints.resize(worksetSize, numCubPoints, pCellDim); 00198 00199 // Step 1.3.c: Map Gauss points to reference face: paramGaussPoints -> refGaussPoints 00200 CellTools::mapToReferenceSubcell(refGaussPoints, 00201 paramGaussPoints, 00202 subcellDim, 00203 subcellOrd, 00204 hexahedron_8); 00205 00206 // Step 1.3.d: Map Gauss points from ref. face to face workset: refGaussPoints -> worksetGaussPoints 00207 CellTools::mapToPhysicalFrame(worksetGaussPoints, 00208 refGaussPoints, 00209 hexNodes, 00210 hexahedron_8); 00211 00212 00213 00214 /************************************************************************************************* 00215 * 00216 * Step 2. Obtain (non-normalized) face normals at cubature points on workset faces 00217 * 2.1 Compute parent cell Jacobians at Gauss points on workset faces 00218 * 2.2 Compute face tangents on workset faces and their vector product 00219 * 00220 ************************************************************************************************/ 00221 00222 // Step 2.1.a: Define and allocate storage for workset Jacobians 00223 FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim); 00224 00225 // Step 2.1.b: Compute Jacobians at Gauss pts. on reference face for all parent cells: 00226 CellTools::setJacobian(worksetJacobians, 00227 refGaussPoints, 00228 hexNodes, 00229 hexahedron_8); 00230 00231 //---------------- 00232 00233 // Step 2.2.a: Allocate storage for face tangents and face normals 00234 FieldContainer<double> worksetFaceTu(worksetSize, numCubPoints, pCellDim); 00235 FieldContainer<double> worksetFaceTv(worksetSize, numCubPoints, pCellDim); 00236 FieldContainer<double> worksetFaceN(worksetSize, numCubPoints, pCellDim); 00237 00238 // Step 2.2.b: Compute face tangents 00239 CellTools::getPhysicalFaceTangents(worksetFaceTu, 00240 worksetFaceTv, 00241 worksetJacobians, 00242 subcellOrd, 00243 hexahedron_8); 00244 00245 // Step 2.2.c: Face outer normals (relative to parent cell) are uTan x vTan: 00246 RealSpaceTools::vecprod(worksetFaceN, worksetFaceTu, worksetFaceTv); 00247 00248 00249 00250 /************************************************************************************************* 00251 * 00252 * Step 3. Evaluate the vector field u(x,y,z) at cubature points on workset faces 00253 * 00254 ************************************************************************************************/ 00255 00256 // Step 3.a: Allocate storage for vector field values at Gauss points on workset faces 00257 FieldContainer<double> worksetVFieldVals(worksetSize, numCubPoints, pCellDim); 00258 00259 // Step 3.b: Compute vector field at Gauss points: here we take u(x,y,z) = (x,y,z) 00260 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){ 00261 for(int ptOrd = 0; ptOrd < numCubPoints; ptOrd++){ 00262 00263 double x = worksetGaussPoints(pCellOrd, ptOrd, 0); 00264 double y = worksetGaussPoints(pCellOrd, ptOrd, 1); 00265 double z = worksetGaussPoints(pCellOrd, ptOrd, 2); 00266 00267 vField(worksetVFieldVals(pCellOrd, ptOrd, 0), 00268 worksetVFieldVals(pCellOrd, ptOrd, 1), 00269 worksetVFieldVals(pCellOrd, ptOrd, 2), x, y, z); 00270 00271 }// pt 00272 }//cell 00273 00274 00275 /************************************************************************************************* 00276 * 00277 * Step 4. Compute dot product of u(x,y,z) and the face normals, times the cubature weights 00278 * 00279 ************************************************************************************************/ 00280 00281 // Allocate storage for dot product of vector field and face normals at Gauss points 00282 FieldContainer<double> worksetFieldDotNormal(worksetSize, numCubPoints); 00283 00284 // Compute the dot product 00285 RealSpaceTools::dot(worksetFieldDotNormal, worksetVFieldVals, worksetFaceN); 00286 00287 // Allocate storage for face fluxes on the workset 00288 FieldContainer<double> worksetFluxes(worksetSize); 00289 00290 //---------------- 00291 00292 // Integration loop (temporary) 00293 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){ 00294 worksetFluxes(pCellOrd) = 0.0; 00295 for(int pt = 0; pt < numCubPoints; pt++ ){ 00296 worksetFluxes(pCellOrd) += worksetFieldDotNormal(pCellOrd, pt)*paramGaussWeights(pt); 00297 }// pt 00298 }//cell 00299 00300 std::cout << "Face fluxes on workset faces : \n\n"; 00301 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){ 00302 00303 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCellOrd, subcellDim, subcellOrd); 00304 std::cout << " Flux = " << worksetFluxes(pCellOrd) << "\n\n"; 00305 00306 } 00307 00308 00309 00310 /************************************************************************************************* 00311 * 00312 * Optional: print Gauss points and face normals at Gauss points 00313 * 00314 ************************************************************************************************/ 00315 00316 // Print Gauss points on [-1,1]x[-1,1] and their images on workset faces 00317 std::cout \ 00318 << "===============================================================================\n" \ 00319 << "| Gauss points on workset faces: |\n" \ 00320 << "===============================================================================\n"; 00321 for(int pCell = 0; pCell < worksetSize; pCell++){ 00322 00323 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd); 00324 00325 for(int pt = 0; pt < numCubPoints; pt++){ 00326 std::cout << "\t 2D Gauss point (" 00327 << std::setw(8) << std::right << paramGaussPoints(pt, 0) << ", " 00328 << std::setw(8) << std::right << paramGaussPoints(pt, 1) << ") " 00329 << std::setw(8) << " --> " << "(" 00330 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", " 00331 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", " 00332 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")\n"; 00333 } 00334 std::cout << "\n\n"; 00335 }//pCell 00336 00337 00338 // Print face normals at Gauss points on workset faces 00339 std::cout \ 00340 << "===============================================================================\n" \ 00341 << "| Face normals (non-unit) at Gauss points on workset faces: |\n" \ 00342 << "===============================================================================\n"; 00343 for(int pCell = 0; pCell < worksetSize; pCell++){ 00344 00345 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd); 00346 00347 for(int pt = 0; pt < numCubPoints; pt++){ 00348 std::cout << "\t 3D Gauss point: (" 00349 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", " 00350 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", " 00351 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")" 00352 << std::setw(8) << " out. normal: " << "(" 00353 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 0) << ", " 00354 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 1) << ", " 00355 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 2) << ")\n"; 00356 } 00357 std::cout << "\n"; 00358 }//pCell 00359 00360 return 0; 00361 } 00362 00363 /************************************************************************************************* 00364 * 00365 * Definition of the vector field function 00366 * 00367 ************************************************************************************************/ 00368 00369 00370 void vField(double& v1, double& v2, double& v3, const double& x, const double& y, const double& z) 00371 { 00372 v1 = x; 00373 v2 = y; 00374 v3 = z; 00375 } 00376 00377
1.7.6.1