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