Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Integration/Intrepid_CubaturePolygonDef.hpp
Go to the documentation of this file.
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 
00050 #include "Intrepid_CubatureDirectTriDefault.hpp"
00051 #include "Intrepid_FieldContainer.hpp"
00052 #include "Intrepid_CellTools.hpp"
00053 #include <vector>
00054 #include <iostream>
00055 
00056 namespace Intrepid{
00057 
00058   template<class Scalar, class ArrayPoint, class ArrayWeight>
00059   CubaturePolygon<Scalar,ArrayPoint,ArrayWeight>::CubaturePolygon(const shards::CellTopology& cellTopology,
00060                                                                   const ArrayPoint& cellVertices,
00061                                                                   int degree)
00062     : degree_(degree), cubDimension_(2), cellTopology_(cellTopology), cellVertices_(cellVertices){
00063     
00064     TEUCHOS_TEST_FOR_EXCEPTION( (degree < 0) || degree > INTREPID_CUBATURE_TRI_DEFAULT_MAX, std::out_of_range,
00065                         ">>> ERROR (CubaturePolygon): No direct cubature rule implemented for the desired polynomial degree.");
00066     // compute area and centroid of polygon
00067     Scalar area;
00068     std::vector<Scalar> centroid(2,0);
00069     int numNodes = cellTopology_.getNodeCount();
00070     
00071     for (int i=0;i<numNodes;i++){
00072       int first = cellTopology_.getNodeMap(1,i,0);
00073       int second = cellTopology_.getNodeMap(1,i,1);
00074       area += cellVertices_(first,0)*cellVertices_(second,1) - cellVertices_(second,0)*cellVertices_(first,1);
00075       centroid[0] += (cellVertices_(first,0) + cellVertices_(second,0))*(cellVertices_(first,0)*cellVertices_(second,1) - cellVertices_(second,0)*cellVertices_(first,1));
00076       centroid[1] += (cellVertices_(first,1) + cellVertices_(second,1))*(cellVertices_(first,0)*cellVertices_(second,1) - cellVertices_(second,0)*cellVertices_(first,1));
00077     }
00078     area /= 2;
00079     centroid[0] /= (6*area);
00080     centroid[1] /= (6*area);
00081         
00082     // get cubature for reference triangle
00083     CubatureDirectTriDefault<Scalar,ArrayPoint,ArrayWeight> cubatureTri(degree_);
00084     int numCubPointsPerTri = cubatureTri.getNumPoints();
00085     int cubDim = cubatureTri.getDimension();
00086     cubDimension_ = cubDim;
00087     FieldContainer<Scalar> cubatureTriPoints(numCubPointsPerTri,cubDim);
00088     
00089     FieldContainer<Scalar> cubatureTriWeights(numCubPointsPerTri);
00090     cubatureTri.getCubature(cubatureTriPoints,cubatureTriWeights);
00091     
00092     // copy into (C,P,D) sized field container where C is the number of triangles in polygon
00093     int numCells = cellTopology_.getEdgeCount();
00094     FieldContainer<Scalar> cubatureCellPoints(numCells,numCubPointsPerTri,cubDim);
00095     for (int k=0;k<numCells;k++){
00096       for (int i=0;i<numCubPointsPerTri;i++){
00097         for (int j=0;j<cubDim;j++){
00098           cubatureCellPoints(k,i,j) = cubatureTriPoints(i,j);
00099         }
00100       }
00101     }
00102     
00103     
00104     // now map cubature to each triangle cell
00105     shards::CellTopology triangleTopology(shards::getCellTopologyData<shards::Triangle<3> >());
00106     int totalCubPoints = numCubPointsPerTri*cellTopology_.getEdgeCount();
00107     numPoints_ = totalCubPoints;
00108     cubaturePoints_.resize(totalCubPoints,cubDim);
00109     cubatureWeights_.resize(totalCubPoints);
00110     
00111     FieldContainer<Scalar> physicalPoints(numCells,numCubPointsPerTri,cubDim);
00112     FieldContainer<Scalar> trianglePoints(numCells,3,cubDim);
00113     int currPoint = 0;
00114     for (int i=0;i<numCells;i++){
00115       for (int j=0;j<cubDim;j++){
00116         trianglePoints(i,0,j) = cellVertices_(cellTopology_.getNodeMap(1,i,0),j);
00117         trianglePoints(i,1,j) = cellVertices_(cellTopology_.getNodeMap(1,i,1),j);
00118         trianglePoints(i,2,j) = centroid[j];
00119       }
00120     }
00121     
00122     CellTools<Scalar>::mapToPhysicalFrame(physicalPoints,cubatureTriPoints,trianglePoints,triangleTopology);
00123     
00124     // compute area of each triangle cell -- need when computing new weights
00125     FieldContainer<Scalar> jacobians(numCells,numCubPointsPerTri,cubDim,cubDim);
00126     FieldContainer<Scalar> detJacobians(numCells, numCubPointsPerTri);
00127 
00128     CellTools<Scalar>::setJacobian(jacobians,physicalPoints,trianglePoints,triangleTopology);
00129     CellTools<Scalar>::setJacobianDet(detJacobians,jacobians);
00130     
00131     for (int i=0;i<numCells;i++){
00132       for (int j=0;j<numCubPointsPerTri;j++){
00133         for (int k=0;k<cubDim;k++){
00134           cubaturePoints_(currPoint,k) = physicalPoints(i,j,k);
00135         }
00136         cubatureWeights_(currPoint++) = cubatureTriWeights(j)*detJacobians(i,j);
00137       }
00138     }
00139   }  // end Constructor
00140 
00141 
00142 template<class Scalar, class ArrayPoint, class ArrayWeight>
00143 void CubaturePolygon<Scalar,ArrayPoint,ArrayWeight>::getCubature(ArrayPoint& cubPoints,
00144                                                                  ArrayWeight& cubWeights)const{
00145   int numCubPoints = numPoints_;
00146   int cellDim = cubDimension_;
00147     
00148   TEUCHOS_TEST_FOR_EXCEPTION ( ( cubPoints.size() < numCubPoints*cellDim || cubWeights.size() < numCubPoints ),
00149                        std::out_of_range,
00150                        ">>> ERROR (CubaturePolygon): Insufficient space allocated for cubature points or weights.");
00151 
00152   for (int pointId = 0; pointId < numCubPoints; pointId++){
00153     for (int dim = 0; dim < cellDim; dim++){
00154       cubPoints(pointId,dim) = cubaturePoints_(pointId,dim);
00155     }
00156     cubWeights(pointId) = cubatureWeights_(pointId);
00157   }
00158 } // end getCubature
00159   
00160 
00161 template<class Scalar, class ArrayPoint, class ArrayWeight>
00162 int CubaturePolygon<Scalar,ArrayPoint,ArrayWeight>::getNumPoints()const{
00163   return numPoints_;
00164 } // end getNumPoints
00165 
00166 template<class Scalar, class ArrayPoint, class ArrayWeight>
00167 int CubaturePolygon<Scalar,ArrayPoint,ArrayWeight>::getDimension()const{
00168   return cubDimension_;
00169 } // end getDimension
00170 
00171 template<class Scalar, class ArrayPoint, class ArrayWeight>
00172 void CubaturePolygon<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(std::vector<int>& accuracy)const{
00173   accuracy.assign(1,degree_);
00174 } // end getAccuracy
00175   
00176 } // namespace Intrepid