Intrepid
Intrepid_CubaturePolygonDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
52 #include "Intrepid_CellTools.hpp"
53 #include <vector>
54 #include <iostream>
55 
56 namespace Intrepid{
57 
58  template<class Scalar, class ArrayPoint, class ArrayWeight>
59  CubaturePolygon<Scalar,ArrayPoint,ArrayWeight>::CubaturePolygon(const shards::CellTopology& cellTopology,
60  const ArrayPoint& cellVertices,
61  int degree)
62  : degree_(degree), cubDimension_(2), cellTopology_(cellTopology), cellVertices_(cellVertices){
63 
64  TEUCHOS_TEST_FOR_EXCEPTION( (degree < 0) || degree > INTREPID_CUBATURE_TRI_DEFAULT_MAX, std::out_of_range,
65  ">>> ERROR (CubaturePolygon): No direct cubature rule implemented for the desired polynomial degree.");
66  // compute area and centroid of polygon
67  Scalar area;
68  std::vector<Scalar> centroid(2,0);
69  int numNodes = cellTopology_.getNodeCount();
70 
71  for (int i=0;i<numNodes;i++){
72  int first = cellTopology_.getNodeMap(1,i,0);
73  int second = cellTopology_.getNodeMap(1,i,1);
74  area += cellVertices_(first,0)*cellVertices_(second,1) - cellVertices_(second,0)*cellVertices_(first,1);
75  centroid[0] += (cellVertices_(first,0) + cellVertices_(second,0))*(cellVertices_(first,0)*cellVertices_(second,1) - cellVertices_(second,0)*cellVertices_(first,1));
76  centroid[1] += (cellVertices_(first,1) + cellVertices_(second,1))*(cellVertices_(first,0)*cellVertices_(second,1) - cellVertices_(second,0)*cellVertices_(first,1));
77  }
78  area /= 2;
79  centroid[0] /= (6*area);
80  centroid[1] /= (6*area);
81 
82  // get cubature for reference triangle
84  int numCubPointsPerTri = cubatureTri.getNumPoints();
85  int cubDim = cubatureTri.getDimension();
86  cubDimension_ = cubDim;
87  FieldContainer<Scalar> cubatureTriPoints(numCubPointsPerTri,cubDim);
88 
89  FieldContainer<Scalar> cubatureTriWeights(numCubPointsPerTri);
90  cubatureTri.getCubature(cubatureTriPoints,cubatureTriWeights);
91 
92  // copy into (C,P,D) sized field container where C is the number of triangles in polygon
93  int numCells = cellTopology_.getEdgeCount();
94  FieldContainer<Scalar> cubatureCellPoints(numCells,numCubPointsPerTri,cubDim);
95  for (int k=0;k<numCells;k++){
96  for (int i=0;i<numCubPointsPerTri;i++){
97  for (int j=0;j<cubDim;j++){
98  cubatureCellPoints(k,i,j) = cubatureTriPoints(i,j);
99  }
100  }
101  }
102 
103 
104  // now map cubature to each triangle cell
105  shards::CellTopology triangleTopology(shards::getCellTopologyData<shards::Triangle<3> >());
106  int totalCubPoints = numCubPointsPerTri*cellTopology_.getEdgeCount();
107  numPoints_ = totalCubPoints;
108  cubaturePoints_.resize(totalCubPoints,cubDim);
109  cubatureWeights_.resize(totalCubPoints);
110 
111  FieldContainer<Scalar> physicalPoints(numCells,numCubPointsPerTri,cubDim);
112  FieldContainer<Scalar> trianglePoints(numCells,3,cubDim);
113  int currPoint = 0;
114  for (int i=0;i<numCells;i++){
115  for (int j=0;j<cubDim;j++){
116  trianglePoints(i,0,j) = cellVertices_(cellTopology_.getNodeMap(1,i,0),j);
117  trianglePoints(i,1,j) = cellVertices_(cellTopology_.getNodeMap(1,i,1),j);
118  trianglePoints(i,2,j) = centroid[j];
119  }
120  }
121 
122  CellTools<Scalar>::mapToPhysicalFrame(physicalPoints,cubatureTriPoints,trianglePoints,triangleTopology);
123 
124  // compute area of each triangle cell -- need when computing new weights
125  FieldContainer<Scalar> jacobians(numCells,numCubPointsPerTri,cubDim,cubDim);
126  FieldContainer<Scalar> detJacobians(numCells, numCubPointsPerTri);
127 
128  CellTools<Scalar>::setJacobian(jacobians,physicalPoints,trianglePoints,triangleTopology);
129  CellTools<Scalar>::setJacobianDet(detJacobians,jacobians);
130 
131  for (int i=0;i<numCells;i++){
132  for (int j=0;j<numCubPointsPerTri;j++){
133  for (int k=0;k<cubDim;k++){
134  cubaturePoints_(currPoint,k) = physicalPoints(i,j,k);
135  }
136  cubatureWeights_(currPoint++) = cubatureTriWeights(j)*detJacobians(i,j);
137  }
138  }
139  } // end Constructor
140 
141 
142 template<class Scalar, class ArrayPoint, class ArrayWeight>
144  ArrayWeight& cubWeights)const{
145  int numCubPoints = numPoints_;
146  int cellDim = cubDimension_;
147 
148  TEUCHOS_TEST_FOR_EXCEPTION ( ( cubPoints.size() < numCubPoints*cellDim || cubWeights.size() < numCubPoints ),
149  std::out_of_range,
150  ">>> ERROR (CubaturePolygon): Insufficient space allocated for cubature points or weights.");
151 
152  for (int pointId = 0; pointId < numCubPoints; pointId++){
153  for (int dim = 0; dim < cellDim; dim++){
154  cubPoints(pointId,dim) = cubaturePoints_(pointId,dim);
155  }
156  cubWeights(pointId) = cubatureWeights_(pointId);
157  }
158 } // end getCubature
159 
160 template<class Scalar, class ArrayPoint, class ArrayWeight>
162  ArrayWeight& cubWeights,
163  ArrayPoint& cellCoords) const
164 {
165  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
166  ">>> ERROR (CubaturePolygon): Cubature defined in reference space calling method for physical space cubature.");
167 }
168 
169 
170 
171 template<class Scalar, class ArrayPoint, class ArrayWeight>
173  return numPoints_;
174 } // end getNumPoints
175 
176 template<class Scalar, class ArrayPoint, class ArrayWeight>
178  return cubDimension_;
179 } // end getDimension
180 
181 template<class Scalar, class ArrayPoint, class ArrayWeight>
183  accuracy.assign(1,degree_);
184 } // end getAccuracy
185 
186 } // namespace Intrepid
int cubDimension_
Dimension of integration domain.
Header file for the Intrepid::CellTools class.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly on each triangle. The return vector ha...
virtual int getDimension() const
Returns dimension of integration domain.
Header file for utility class to provide multidimensional containers.
#define INTREPID_CUBATURE_TRI_DEFAULT_MAX
The maximum degree of the polynomial that can be integrated exactly by a direct triangle rule of the ...
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
int getDimension() const
Returns dimension of integration domain.
int getNumPoints() const
Returns the number of cubature points.
shards::CellTopology cellTopology_
The topology of the polygon.
int degree_
The degree of the polynomials that are integrated exactly on each triangle.
FieldContainer< Scalar > cubatureWeights_
Local copy of cubature weights.
FieldContainer< Scalar > cubaturePoints_
Local copy of cubature points.
Header file for the Intrepid::CubatureDirectTriDefault class.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
Defines direct integration rules on a triangle.
virtual int getNumPoints() const
Returns the number of cubature points.
CubaturePolygon(const shards::CellTopology &cellTopology, const ArrayPoint &cellVertices, int degree)
int numPoints_
The number of cubature points.
A stateless class for operations on cell data. Provides methods for:
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F...
ArrayPoint cellVertices_
The vertices of the polygon.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).