|
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 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00030 00035 #include "Intrepid_CellTools.hpp" 00036 #include "Intrepid_FieldContainer.hpp" 00037 #include "Intrepid_DefaultCubatureFactory.hpp" 00038 #include "Shards_CellTopology.hpp" 00039 #include "Teuchos_GlobalMPISession.hpp" 00040 00045 void vField(double& v1, double& v2, double& v3, 00046 const double& x, const double& y, const double& z); 00047 00048 using namespace std; 00049 using namespace Intrepid; 00050 00051 int main(int argc, char *argv[]) { 00052 00053 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00054 00055 typedef CellTools<double> CellTools; 00056 typedef RealSpaceTools<double> RealSpaceTools; 00057 typedef shards::CellTopology CellTopology; 00058 00059 std::cout \ 00060 << "===============================================================================\n" \ 00061 << "| |\n" \ 00062 << "| Example use of the CellTools class |\n" \ 00063 << "| |\n" \ 00064 << "| 1) Computation of face flux, for a given vector field, on a face workset |\n" \ 00065 << "| 2) Computation of edge circulation, for a given vector field, on a face |\n" \ 00066 << "| workset. |\n" \ 00067 << "| |\n" \ 00068 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \ 00069 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \ 00070 << "| Kara Peterson (kjpeter@sandia.gov) |\n" \ 00071 << "| |\n" \ 00072 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00073 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00074 << "| |\n" \ 00075 << "===============================================================================\n"\ 00076 << "| EXAMPLE 1: Computation of face flux on a face workset |\n"\ 00077 << "===============================================================================\n"; 00078 00079 00094 /************************************************************************************************* 00095 * 00096 * Step 0. Face workset comprising of 1 face of a Hexahedron<8> cell 00097 * 00098 ************************************************************************************************/ 00099 00100 // Step 0.a: Specify cell topology of the parent cell 00101 CellTopology hexahedron_8( shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00102 00103 // Step 0.b: Specify the vertices of the parent Hexahedron<8> cell 00104 int worksetSize = 2; 00105 int pCellNodeCount = hexahedron_8.getVertexCount(); 00106 int pCellDim = hexahedron_8.getDimension(); 00107 00108 FieldContainer<double> hexNodes(worksetSize, pCellNodeCount, pCellDim); 00109 // cell 0 bottom face vertices: 00110 hexNodes(0, 0, 0) = 0.00; hexNodes(0, 0, 1) = 0.00, hexNodes(0, 0, 2) = 0.00; 00111 hexNodes(0, 1, 0) = 1.00; hexNodes(0, 1, 1) = 0.00, hexNodes(0, 1, 2) = 0.00; 00112 hexNodes(0, 2, 0) = 1.00; hexNodes(0, 2, 1) = 1.00, hexNodes(0, 2, 2) = 0.00; 00113 hexNodes(0, 3, 0) = 0.00; hexNodes(0, 3, 1) = 1.00, hexNodes(0, 3, 2) = 0.00; 00114 // cell 0 top face vertices 00115 hexNodes(0, 4, 0) = 0.00; hexNodes(0, 4, 1) = 0.00, hexNodes(0, 4, 2) = 1.00; 00116 hexNodes(0, 5, 0) = 1.00; hexNodes(0, 5, 1) = 0.00, hexNodes(0, 5, 2) = 1.00; 00117 hexNodes(0, 6, 0) = 1.00; hexNodes(0, 6, 1) = 1.00, hexNodes(0, 6, 2) = 1.00; 00118 hexNodes(0, 7, 0) = 0.00; hexNodes(0, 7, 1) = 1.00, hexNodes(0, 7, 2) = 1.00; 00119 00120 // cell 1 bottom face vertices: 00121 hexNodes(1, 0, 0) = 0.00; hexNodes(1, 0, 1) = 0.00, hexNodes(1, 0, 2) = 0.00; 00122 hexNodes(1, 1, 0) = 1.00; hexNodes(1, 1, 1) = 0.00, hexNodes(1, 1, 2) = 0.00; 00123 hexNodes(1, 2, 0) = 1.00; hexNodes(1, 2, 1) = 1.00, hexNodes(1, 2, 2) = 0.00; 00124 hexNodes(1, 3, 0) = 0.00; hexNodes(1, 3, 1) = 1.00, hexNodes(1, 3, 2) = 0.00; 00125 // cell 1 top face vertices 00126 hexNodes(1, 4, 0) = 0.00; hexNodes(1, 4, 1) = 0.00, hexNodes(1, 4, 2) = 1.00; 00127 hexNodes(1, 5, 0) = 1.00; hexNodes(1, 5, 1) = 0.00, hexNodes(1, 5, 2) = 1.00; 00128 hexNodes(1, 6, 0) = 1.00; hexNodes(1, 6, 1) = 1.00, hexNodes(1, 6, 2) = 0.75; 00129 hexNodes(1, 7, 0) = 0.00; hexNodes(1, 7, 1) = 1.00, hexNodes(1, 7, 2) = 1.00; 00130 00131 00132 // Step 0.c: Specify the face ordinal, relative to the reference cell, of the face workset 00133 int subcellDim = 2; 00134 int subcellOrd = 1; 00135 00136 00137 00138 /************************************************************************************************* 00139 * 00140 * Step 1: Obtain Gauss points on workset faces for Hexahedron<8> topology 00141 * 1.1 Define cubature factory, face parametrization domain and arrays for cubature points 00142 * 1.2 Select Gauss rule on D = [-1,1]x[-1,1] 00143 * 1.3 Map Gauss points from D to reference face and workset faces 00144 * 00145 ************************************************************************************************/ 00146 00147 // Step 1.1.a: Define CubatureFactory 00148 DefaultCubatureFactory<double> cubFactory; 00149 00150 // Step 1.1.b: Define topology of the face parametrization domain as [-1,1]x[-1,1] 00151 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00152 00153 // Step 1.1.c: Define storage for cubature points and weights on [-1,1]x[-1,1] 00154 FieldContainer<double> paramGaussWeights; 00155 FieldContainer<double> paramGaussPoints; 00156 00157 // Step 1.1.d: Define storage for cubature points on a reference face 00158 FieldContainer<double> refGaussPoints; 00159 00160 // Step 1.1.f: Define storage for cubature points on workset faces 00161 FieldContainer<double> worksetGaussPoints; 00162 00163 //---------------- 00164 00165 // Step 1.2.a: selects Gauss rule of order 3 on [-1,1]x[-1,1] 00166 Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactory.create(paramQuadFace, 3); 00167 00168 // Step 1.2.b allocate storage for cubature points on [-1,1]x[-1,1] 00169 int cubDim = hexFaceCubature -> getDimension(); 00170 int numCubPoints = hexFaceCubature -> getNumPoints(); 00171 00172 // Arrays must be properly sized for the specified set of Gauss points 00173 paramGaussPoints.resize(numCubPoints, cubDim); 00174 paramGaussWeights.resize(numCubPoints); 00175 hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights); 00176 00177 //---------------- 00178 00179 // Step 1.3.a: Allocate storage for Gauss points on the reference face 00180 refGaussPoints.resize(numCubPoints, pCellDim); 00181 00182 // Step 1.3.b: Allocate storage for Gauss points on the face in the workset 00183 worksetGaussPoints.resize(worksetSize, numCubPoints, pCellDim); 00184 00185 // Step 1.3.c: Map Gauss points to reference face: paramGaussPoints -> refGaussPoints 00186 CellTools::mapToReferenceSubcell(refGaussPoints, 00187 paramGaussPoints, 00188 subcellDim, 00189 subcellOrd, 00190 hexahedron_8); 00191 00192 // Step 1.3.d: Map Gauss points from ref. face to face workset: refGaussPoints -> worksetGaussPoints 00193 CellTools::mapToPhysicalFrame(worksetGaussPoints, 00194 refGaussPoints, 00195 hexNodes, 00196 hexahedron_8); 00197 00198 00199 00200 /************************************************************************************************* 00201 * 00202 * Step 2. Obtain (non-normalized) face normals at cubature points on workset faces 00203 * 2.1 Compute parent cell Jacobians at Gauss points on workset faces 00204 * 2.2 Compute face tangents on workset faces and their vector product 00205 * 00206 ************************************************************************************************/ 00207 00208 // Step 2.1.a: Define and allocate storage for workset Jacobians 00209 FieldContainer<double> worksetJacobians(worksetSize, numCubPoints, pCellDim, pCellDim); 00210 00211 // Step 2.1.b: Compute Jacobians at Gauss pts. on reference face for all parent cells: 00212 CellTools::setJacobian(worksetJacobians, 00213 refGaussPoints, 00214 hexNodes, 00215 hexahedron_8); 00216 00217 //---------------- 00218 00219 // Step 2.2.a: Allocate storage for face tangents and face normals 00220 FieldContainer<double> worksetFaceTu(worksetSize, numCubPoints, pCellDim); 00221 FieldContainer<double> worksetFaceTv(worksetSize, numCubPoints, pCellDim); 00222 FieldContainer<double> worksetFaceN(worksetSize, numCubPoints, pCellDim); 00223 00224 // Step 2.2.b: Compute face tangents 00225 CellTools::getPhysicalFaceTangents(worksetFaceTu, 00226 worksetFaceTv, 00227 worksetJacobians, 00228 subcellOrd, 00229 hexahedron_8); 00230 00231 // Step 2.2.c: Face outer normals (relative to parent cell) are uTan x vTan: 00232 RealSpaceTools::vecprod(worksetFaceN, worksetFaceTu, worksetFaceTv); 00233 00234 00235 00236 /************************************************************************************************* 00237 * 00238 * Step 3. Evaluate the vector field u(x,y,z) at cubature points on workset faces 00239 * 00240 ************************************************************************************************/ 00241 00242 // Step 3.a: Allocate storage for vector field values at Gauss points on workset faces 00243 FieldContainer<double> worksetVFieldVals(worksetSize, numCubPoints, pCellDim); 00244 00245 // Step 3.b: Compute vector field at Gauss points: here we take u(x,y,z) = (x,y,z) 00246 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){ 00247 for(int ptOrd = 0; ptOrd < numCubPoints; ptOrd++){ 00248 00249 double x = worksetGaussPoints(pCellOrd, ptOrd, 0); 00250 double y = worksetGaussPoints(pCellOrd, ptOrd, 1); 00251 double z = worksetGaussPoints(pCellOrd, ptOrd, 2); 00252 00253 vField(worksetVFieldVals(pCellOrd, ptOrd, 0), 00254 worksetVFieldVals(pCellOrd, ptOrd, 1), 00255 worksetVFieldVals(pCellOrd, ptOrd, 2), x, y, z); 00256 00257 }// pt 00258 }//cell 00259 00260 00261 /************************************************************************************************* 00262 * 00263 * Step 4. Compute dot product of u(x,y,z) and the face normals, times the cubature weights 00264 * 00265 ************************************************************************************************/ 00266 00267 // Allocate storage for dot product of vector field and face normals at Gauss points 00268 FieldContainer<double> worksetFieldDotNormal(worksetSize, numCubPoints); 00269 00270 // Compute the dot product 00271 RealSpaceTools::dot(worksetFieldDotNormal, worksetVFieldVals, worksetFaceN); 00272 00273 // Allocate storage for face fluxes on the workset 00274 FieldContainer<double> worksetFluxes(worksetSize); 00275 00276 //---------------- 00277 00278 // Integration loop (temporary) 00279 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){ 00280 worksetFluxes(pCellOrd) = 0.0; 00281 for(int pt = 0; pt < numCubPoints; pt++ ){ 00282 worksetFluxes(pCellOrd) += worksetFieldDotNormal(pCellOrd, pt)*paramGaussWeights(pt); 00283 }// pt 00284 }//cell 00285 00286 std::cout << "Face fluxes on workset faces : \n\n"; 00287 for(int pCellOrd = 0; pCellOrd < worksetSize; pCellOrd++){ 00288 00289 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCellOrd, subcellDim, subcellOrd); 00290 std::cout << " Flux = " << worksetFluxes(pCellOrd) << "\n\n"; 00291 00292 } 00293 00294 00295 00296 /************************************************************************************************* 00297 * 00298 * Optional: print Gauss points and face normals at Gauss points 00299 * 00300 ************************************************************************************************/ 00301 00302 // Print Gauss points on [-1,1]x[-1,1] and their images on workset faces 00303 std::cout \ 00304 << "===============================================================================\n" \ 00305 << "| Gauss points on workset faces: |\n" \ 00306 << "===============================================================================\n"; 00307 for(int pCell = 0; pCell < worksetSize; pCell++){ 00308 00309 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd); 00310 00311 for(int pt = 0; pt < numCubPoints; pt++){ 00312 std::cout << "\t 2D Gauss point (" 00313 << std::setw(8) << std::right << paramGaussPoints(pt, 0) << ", " 00314 << std::setw(8) << std::right << paramGaussPoints(pt, 1) << ") " 00315 << std::setw(8) << " --> " << "(" 00316 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", " 00317 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", " 00318 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")\n"; 00319 } 00320 std::cout << "\n\n"; 00321 }//pCell 00322 00323 00324 // Print face normals at Gauss points on workset faces 00325 std::cout \ 00326 << "===============================================================================\n" \ 00327 << "| Face normals (non-unit) at Gauss points on workset faces: |\n" \ 00328 << "===============================================================================\n"; 00329 for(int pCell = 0; pCell < worksetSize; pCell++){ 00330 00331 CellTools::printWorksetSubcell(hexNodes, hexahedron_8, pCell, subcellDim, subcellOrd); 00332 00333 for(int pt = 0; pt < numCubPoints; pt++){ 00334 std::cout << "\t 3D Gauss point: (" 00335 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 0) << ", " 00336 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 1) << ", " 00337 << std::setw(8) << std::right << worksetGaussPoints(pCell, pt, 2) << ")" 00338 << std::setw(8) << " out. normal: " << "(" 00339 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 0) << ", " 00340 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 1) << ", " 00341 << std::setw(8) << std::right << worksetFaceN(pCell, pt, 2) << ")\n"; 00342 } 00343 std::cout << "\n"; 00344 }//pCell 00345 00346 return 0; 00347 } 00348 00349 /************************************************************************************************* 00350 * 00351 * Definition of the vector field function 00352 * 00353 ************************************************************************************************/ 00354 00355 00356 void vField(double& v1, double& v2, double& v3, const double& x, const double& y, const double& z) 00357 { 00358 v1 = x; 00359 v2 = y; 00360 v3 = z; 00361 } 00362 00363
1.7.4