|
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), 00025 // Denis Ridzal (dridzal@sandia.gov), or 00026 // Kara Peterson (kjpeter@sandia.gov) 00027 // 00028 // ************************************************************************ 00029 // @HEADER 00030 00031 00037 #include "Intrepid_CellTools.hpp" 00038 #include "Intrepid_FieldContainer.hpp" 00039 #include "Intrepid_DefaultCubatureFactory.hpp" 00040 #include "Teuchos_GlobalMPISession.hpp" 00041 00042 #include "Shards_CellTopology.hpp" 00043 00044 #include "Teuchos_RCP.hpp" 00045 00046 using namespace std; 00047 using namespace Intrepid; 00048 using namespace shards; 00049 00050 int main(int argc, char *argv[]) { 00051 00052 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00053 00054 typedef CellTools<double> CellTools; 00055 typedef shards::CellTopology CellTopology; 00056 00057 cout \ 00058 << "===============================================================================\n" \ 00059 << "| |\n" \ 00060 << "| Example use of the CellTools class |\n" \ 00061 << "| |\n" \ 00062 << "| 1) Reference edge parametrizations |\n" \ 00063 << "| 2) Reference face parametrizations |\n" \ 00064 << "| |\n" \ 00065 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) |\n" \ 00066 << "| Denis Ridzal (dridzal@sandia.gov), or |\n" \ 00067 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00068 << "| |\n" \ 00069 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00070 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00071 << "| |\n" \ 00072 << "===============================================================================\n"\ 00073 << "| Summary: |\n"\ 00074 << "| Reference edge parametrizations map [-1,1] to the edges of reference cells. |\n"\ 00075 << "| They are used to define, e.g., integration points on the edges of 2D and 3D |\n"\ 00076 << "| reference cells. Edge parametrizations for special 2D cells such as Beam |\n"\ 00077 << "| and ShellLine, are also supported. |\n"\ 00078 << "===============================================================================\n"; 00079 00080 /* 00081 Specification of integration points on 1-subcells (edges) of reference cells. Edges are 00082 parametrized by [-1,1] and integration points on an edge are defined by mapping integration 00083 points from the parametrization domain [-1,1] to a specific edge on the reference cell. 00084 00085 1. Common tasks: definition of integration points in the edge parametrization domain [-1,1] 00086 These steps are independent of parent cell topology: 00087 00088 a. Instantiate a CubatureFactory object to create cubatures (needed for face maps too) 00089 b. Define parametrization domain for the edges as having Line<2> cell topology. This is 00090 required by the CubatureFactory in order to select cubature points and weights from 00091 the reference line [-1,1] 00092 c. Use CubatureFactory to select cubature of the desired degree for the Line<2> topology 00093 d. Allocate containers for the cubature points and weights. 00094 00095 2. Parent cell topology specific tasks 00096 00097 a. Select the parent cell topology 00098 b. Allocate containers for the images of the integration points on [-1,1] on the edges 00099 c. Apply the edge parametrization map to the pointss in [-1,1] 00100 */ 00101 00102 // Step 1.a (Define CubatureFactory) 00103 DefaultCubatureFactory<double> cubFactory; 00104 00105 00106 // Step 1.b (Define the topology of the parametrization domain) 00107 CellTopology edgeParam(shards::getCellTopologyData<shards::Line<2> >() ); 00108 00109 00110 // Step 1.c (selects Gauss rule of order 6 on [-1,1]) 00111 Teuchos::RCP<Cubature<double> > edgeParamCubature = cubFactory.create(edgeParam, 6); 00112 00113 00114 // Step 1.d (allocate storage for cubature points) 00115 int cubDim = edgeParamCubature -> getDimension(); 00116 int numCubPoints = edgeParamCubature -> getNumPoints(); 00117 00118 FieldContainer<double> edgeParamCubPts(numCubPoints, cubDim); 00119 FieldContainer<double> edgeParamCubWts(numCubPoints); 00120 edgeParamCubature -> getCubature(edgeParamCubPts, edgeParamCubWts); 00121 00122 00123 00124 std::cout \ 00125 << "===============================================================================\n"\ 00126 << "| EXAMPLE 1.1 |\n" 00127 << "| Edge parametrizations for standard 2D cells: Triangle |\n"\ 00128 << "===============================================================================\n"; 00129 00130 // Step 2.a (select reference cell topology) 00131 CellTopology triangle_3(getCellTopologyData<Triangle<3> >() ); 00132 00133 00134 // Step 2.b (allocate storage for points on an edge of the reference cell) 00135 FieldContainer<double> triEdgePoints(numCubPoints, triangle_3.getDimension() ); 00136 00137 00138 // Step 2.c (same points are mapped to all edges, can also map different set to each edge) 00139 for(int edgeOrd = 0; edgeOrd < (int)triangle_3.getEdgeCount(); edgeOrd++){ 00140 00141 CellTools::mapToReferenceSubcell(triEdgePoints, 00142 edgeParamCubPts, 00143 1, 00144 edgeOrd, 00145 triangle_3); 00146 00147 // Optional: print the vertices of the reference edge 00148 CellTools::printSubcellVertices(1, edgeOrd, triangle_3); 00149 00150 for(int pt = 0; pt < numCubPoints; pt++){ 00151 std::cout << "\t Parameter point " 00152 << std::setw(12) << std::right << edgeParamCubPts(pt, 0) << std::setw(10) << " --> " << "(" 00153 << std::setw(10) << std::right << triEdgePoints(pt, 0) << " , " 00154 << std::setw(10) << std::right << triEdgePoints(pt, 1) << ")\n"; 00155 } 00156 std::cout << "\n"; 00157 } 00158 00159 00160 00161 std::cout \ 00162 << "===============================================================================\n"\ 00163 << "| EXAMPLE 1.2 |\n" 00164 << "| Edge parametrizations for standard 2D cells: Quadrilateral |\n"\ 00165 << "===============================================================================\n"; 00166 00167 // Step 2.a (select reference cell topology) 00168 CellTopology quad_4(getCellTopologyData<Quadrilateral<4> >() ); 00169 00170 00171 // Step 2.b (allocate storage for points on an edge of the reference cell) 00172 FieldContainer<double> quadEdgePoints(numCubPoints, quad_4.getDimension() ); 00173 00174 00175 // Step 2.c (same points are mapped to all edges, can also map different set to each edge) 00176 for(int edgeOrd = 0; edgeOrd < (int)quad_4.getEdgeCount(); edgeOrd++){ 00177 00178 CellTools::mapToReferenceSubcell(quadEdgePoints, 00179 edgeParamCubPts, 00180 1, 00181 edgeOrd, 00182 quad_4); 00183 00184 // Optional: print the vertices of the reference edge 00185 CellTools::printSubcellVertices(1, edgeOrd, quad_4); 00186 00187 for(int pt = 0; pt < numCubPoints; pt++){ 00188 std::cout << "\t Parameter point " 00189 << std::setw(12) << std::right << edgeParamCubPts(pt, 0) << std::setw(10) << " --> " << "(" 00190 << std::setw(10) << std::right << quadEdgePoints(pt, 0) << " , " 00191 << std::setw(10) << std::right << quadEdgePoints(pt, 1) << ")\n"; 00192 } 00193 std::cout << "\n"; 00194 } 00195 00196 00197 /* 00198 Specification of integration points on 2-subcells (faces) of reference cells. Reference cells 00199 can have triangular, quadrilateral or a mixture of triangular and quadrilateral faces. Thus, 00200 parametrization domain of a face depends on that face's topology and is either the standard 00201 2-simplex {(0,0), (1,0), (0,1)} for triangular faces or the standard 2-cube [-1,1]^2 for 00202 quadrilateral faces. 00203 00204 1. Common tasks: definition of integration points in the standard 2-simplex and the standard 00205 2-cube. These steps are independent of parent cell topology: 00206 00207 a. Instantiate a CubatureFactory object to create cubatures (already done!) 00208 b. Define parametrization domain for traingular faces as having Triangle<3> topology and for 00209 quadrilateral faces - as having Quadrilateral<4> topology. This is required by the 00210 CubatureFactory in order to select cubature points and weights from the appropriate 00211 face parametrization domain. 00212 c. Use CubatureFactory to select cubature of the desired degree for Triangle<3> and 00213 Quadrilateral<4> topologies 00214 d. Allocate containers for the cubature points and weights on the parametrization domains. 00215 00216 2. Parent cell topology specific tasks 00217 00218 a. Select the parent cell topology 00219 b. Allocate containers for the images of the integration points from the parametrization 00220 domains on the reference faces 00221 c. Apply the face parametrization map to the points in the parametrization domain 00222 */ 00223 00224 // Step 1.b (Define the topology of the parametrization domain) 00225 CellTopology triFaceParam(shards::getCellTopologyData<shards::Triangle<3> >() ); 00226 CellTopology quadFaceParam(shards::getCellTopologyData<shards::Quadrilateral<4> >() ); 00227 00228 00229 // Step 1.c (selects Gauss rule of order 3 on [-1,1]^2 and a rule of order 3 on Triangle) 00230 Teuchos::RCP<Cubature<double> > triFaceParamCubature = cubFactory.create(triFaceParam, 3); 00231 Teuchos::RCP<Cubature<double> > quadFaceParamCubature = cubFactory.create(quadFaceParam, 3); 00232 00233 00234 // Step 1.d - Triangle faces (allocate storage for cubature points) 00235 int triFaceCubDim = triFaceParamCubature -> getDimension(); 00236 int triFaceNumCubPts = triFaceParamCubature -> getNumPoints(); 00237 00238 FieldContainer<double> triFaceParamCubPts(triFaceNumCubPts, triFaceCubDim); 00239 FieldContainer<double> triFaceParamCubWts(triFaceNumCubPts); 00240 triFaceParamCubature -> getCubature(triFaceParamCubPts, triFaceParamCubWts); 00241 00242 00243 // Step 1.d - Quadrilateral faces (allocate storage for cubature points) 00244 int quadFaceCubDim = quadFaceParamCubature -> getDimension(); 00245 int quadFaceNumCubPts = quadFaceParamCubature -> getNumPoints(); 00246 00247 FieldContainer<double> quadFaceParamCubPts(quadFaceNumCubPts, quadFaceCubDim); 00248 FieldContainer<double> quadFaceParamCubWts(quadFaceNumCubPts); 00249 quadFaceParamCubature -> getCubature(quadFaceParamCubPts, quadFaceParamCubWts); 00250 00251 00252 00253 std::cout \ 00254 << "===============================================================================\n"\ 00255 << "| EXAMPLE 2.1 |\n" 00256 << "| Face parametrizations for standard 3D cells: Tetrahedron |\n"\ 00257 << "===============================================================================\n"; 00258 00259 // Step 2.a (select reference cell topology) 00260 CellTopology tet_4(getCellTopologyData<Tetrahedron<4> >() ); 00261 00262 00263 // Step 2.b (allocate storage for points on a face of the reference cell) 00264 FieldContainer<double> tetFacePoints(triFaceNumCubPts, tet_4.getDimension() ); 00265 00266 00267 // Step 2.c (same points are mapped to all faces, can also map different set to each face) 00268 for(int faceOrd = 0; faceOrd < (int)tet_4.getSideCount(); faceOrd++){ 00269 00270 CellTools::mapToReferenceSubcell(tetFacePoints, 00271 triFaceParamCubPts, 00272 2, 00273 faceOrd, 00274 tet_4); 00275 00276 // Optional: print the vertices of the reference face 00277 CellTools::printSubcellVertices(2, faceOrd, tet_4); 00278 00279 for(int pt = 0; pt < triFaceNumCubPts; pt++){ 00280 std::cout << "\t Parameter point (" 00281 << std::setw(10) << std::right << triFaceParamCubPts(pt, 0) << " , " 00282 << std::setw(10) << std::right << triFaceParamCubPts(pt, 1) << ") " 00283 << std::setw(10) << " --> " << "(" 00284 << std::setw(10) << std::right << tetFacePoints(pt, 0) << " , " 00285 << std::setw(10) << std::right << tetFacePoints(pt, 1) << " , " 00286 << std::setw(10) << std::right << tetFacePoints(pt, 2) << ")\n"; 00287 } 00288 std::cout << "\n"; 00289 } 00290 00291 00292 00293 std::cout \ 00294 << "===============================================================================\n"\ 00295 << "| EXAMPLE 2.2 |\n" 00296 << "| Face parametrizations for standard 3D cells: Wedge |\n"\ 00297 << "| Example of a reference cell that has two different kinds of faces |\n"\ 00298 << "===============================================================================\n"; 00299 00300 00301 // Step 2.a (select reference cell topology) 00302 CellTopology wedge_6(getCellTopologyData<Wedge<6> >() ); 00303 00304 00305 // Step 2.b (allocate storage for points on a face of the reference cell) 00306 // Wedge<6> has Triangle<3> and Quadrilateral<4> faces. Two different arrays are needed 00307 // to store the points on each face because different types integration rules are used 00308 // on their respective parametrization domains and numbers of points defined by these 00309 // rules do not necessarily match. 00310 FieldContainer<double> wedgeTriFacePoints(triFaceNumCubPts, wedge_6.getDimension() ); 00311 FieldContainer<double> wedgeQuadFacePoints(quadFaceNumCubPts, wedge_6.getDimension() ); 00312 00313 00314 // Step 2.c (for Wedge<6> need to distinguish Triangle<3> and Quadrilateral<4> faces) 00315 for(int faceOrd = 0; faceOrd < (int)wedge_6.getSideCount(); faceOrd++){ 00316 00317 // Optional: print the vertices of the reference face 00318 CellTools::printSubcellVertices(2, faceOrd, wedge_6); 00319 00320 if( wedge_6.getKey(2, faceOrd) == shards::Triangle<3>::key ){ 00321 CellTools::mapToReferenceSubcell(wedgeTriFacePoints, 00322 triFaceParamCubPts, 00323 2, 00324 faceOrd, 00325 wedge_6); 00326 00327 for(int pt = 0; pt < triFaceNumCubPts; pt++){ 00328 std::cout << "\t Parameter point (" 00329 << std::setw(10) << std::right << triFaceParamCubPts(pt, 0) << " , " 00330 << std::setw(10) << std::right << triFaceParamCubPts(pt, 1) << ") " 00331 << std::setw(10) << " --> " << "(" 00332 << std::setw(10) << std::right << wedgeTriFacePoints(pt, 0) << " , " 00333 << std::setw(10) << std::right << wedgeTriFacePoints(pt, 1) << " , " 00334 << std::setw(10) << std::right << wedgeTriFacePoints(pt, 2) << ")\n"; 00335 } 00336 std::cout << "\n"; 00337 } 00338 else if(wedge_6.getKey(2, faceOrd) == shards::Quadrilateral<4>::key) { 00339 CellTools::mapToReferenceSubcell(wedgeQuadFacePoints, 00340 quadFaceParamCubPts, 00341 2, 00342 faceOrd, 00343 wedge_6); 00344 00345 for(int pt = 0; pt < quadFaceNumCubPts; pt++){ 00346 std::cout << "\t Parameter point (" 00347 << std::setw(10) << std::right << quadFaceParamCubPts(pt, 0) << " , " 00348 << std::setw(10) << std::right << quadFaceParamCubPts(pt, 1) << ") " 00349 << std::setw(10) << " --> " << "(" 00350 << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 0) << " , " 00351 << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 1) << " , " 00352 << std::setw(10) << std::right << wedgeQuadFacePoints(pt, 2) << ")\n"; 00353 } 00354 std::cout << "\n"; 00355 } 00356 else { 00357 std::cout << " Invalid face encountered \n"; 00358 } 00359 } 00360 00361 return 0; 00362 }
1.7.4