Intrepid
Intrepid_CubatureControlVolumeSideDef.hpp
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
49#ifndef INTREPID_CUBATURE_CONTROLVOLUMESIDEDEF_HPP
50#define INTREPID_CUBATURE_CONTROLVOLUMESIDEDEF_HPP
51
52namespace Intrepid{
53
54template<class Scalar, class ArrayPoint, class ArrayWeight>
55CubatureControlVolumeSide<Scalar,ArrayPoint,ArrayWeight>::CubatureControlVolumeSide(const Teuchos::RCP<const shards::CellTopology>& cellTopology)
56{
57 // topology of primary cell
58 primaryCellTopo_ = cellTopology;
59
60 // get topology of sub-control volume (will be quad or hex depending on dimension)
61 const CellTopologyData &myCellData =
62 (primaryCellTopo_->getDimension() > 2) ? *shards::getCellTopologyData<shards::Hexahedron<8> >() :
63 *shards::getCellTopologyData<shards::Quadrilateral<4> >();
64
65 subCVCellTopo_ = Teuchos::rcp(new shards::CellTopology(&myCellData));
66
67 degree_ = 1;
68
69 numPoints_ = primaryCellTopo_->getEdgeCount();
70
71 cubDimension_ = primaryCellTopo_->getDimension();
72
73}
74
75template<class Scalar, class ArrayPoint, class ArrayWeight>
77 ArrayWeight& cubWeights) const
78{
79 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
80 ">>> ERROR (CubatureControlVolumeSide): Cubature defined in physical space calling method for reference space cubature.");
81}
82
83template<class Scalar, class ArrayPoint, class ArrayWeight>
85 ArrayWeight& cubWeights,
86 ArrayPoint& cellCoords) const
87{
88 // get array dimensions
89 int numCells = cellCoords.dimension(0);
90 int numNodesPerCell = cellCoords.dimension(1);
91 int spaceDim = cellCoords.dimension(2);
92 int numNodesPerSubCV = subCVCellTopo_->getNodeCount();
93
94 // get sub-control volume coordinates (one sub-control volume per node of primary cell)
95 Intrepid::FieldContainer<Scalar> subCVCoords(numCells,numNodesPerCell,numNodesPerSubCV,spaceDim);
96 Intrepid::CellTools<Scalar>::getSubCVCoords(subCVCoords,cellCoords,*(primaryCellTopo_));
97
98 // num edges per primary cell
99 int numEdgesPerCell = primaryCellTopo_->getEdgeCount();
100
101 // Loop over cells
102 for (int icell = 0; icell < numCells; icell++){
103
104 // Get subcontrol volume side midpoints and normals
105 int iside = 1;
106 int numNodesPerSide = subCVCellTopo_->getNodeCount(spaceDim-1,iside);
107 Intrepid::FieldContainer<int> sideNodes(numNodesPerSide);
108 for (int i=0; i<numNodesPerSide; i++){
109 sideNodes(i) = subCVCellTopo_->getNodeMap(spaceDim-1,iside,i);
110 }
111
112 // Loop over primary cell nodes and get side midpoints
113 // In each primary cell the number of control volume side integration
114 // points is equal to the number of primary cell edges. In 2d the
115 // number of edges = number of nodes and this loop defines all side
116 // points. In 3d this loop computes the side points for all
117 // subcontrol volume sides for iside = 1. Additional code below
118 // computes the remaining points for particular 3d topologies.
119 for (int inode=0; inode < numNodesPerCell; inode++){
120 for(int idim=0; idim < spaceDim; idim++){
121 Scalar midpt = 0.0;
122 for (int i=0; i<numNodesPerSide; i++){
123 midpt += subCVCoords(icell,inode,sideNodes(i),idim);
124 }
125 cubPoints(icell,inode,idim) = midpt/numNodesPerSide;
126 }
127 }
128
129 // Map side center to reference subcell
130 //Intrepid::FieldContainer<Scalar> sideCenterLocal(1,spaceDim-1);
131 Intrepid::FieldContainer<double> sideCenterLocal(1,spaceDim-1);
132 for (int idim = 0; idim < spaceDim-1; idim++){
133 sideCenterLocal(0,idim) = 0.0;
134 }
135
136 Intrepid::FieldContainer<Scalar> refSidePoints(1,spaceDim);
137 iside = 1;
139 sideCenterLocal,
140 spaceDim-1, iside, *(subCVCellTopo_));
141
142 // Array of cell control volume coordinates
143 Intrepid::FieldContainer<Scalar> cellCVCoords(numNodesPerCell, numNodesPerSubCV, spaceDim);
144 for (int isubcv = 0; isubcv < numNodesPerCell; isubcv++) {
145 for (int inode = 0; inode < numNodesPerSubCV; inode++){
146 for (int idim = 0; idim < spaceDim; idim++){
147 cellCVCoords(isubcv,inode,idim) = subCVCoords(icell,isubcv,inode,idim);
148 }
149 }
150 }
151
152 // calculate Jacobian at side centers
153 Intrepid::FieldContainer<Scalar> subCVsideJacobian(numNodesPerCell, 1, spaceDim, spaceDim);
154 Intrepid::CellTools<Scalar>::setJacobian(subCVsideJacobian, refSidePoints, cellCVCoords, *(subCVCellTopo_));
155
156 // Get subcontrol volume side normals
157 Intrepid::FieldContainer<Scalar> normals(numNodesPerCell, 1, spaceDim);
158 Intrepid::CellTools<Scalar>::getPhysicalSideNormals(normals,subCVsideJacobian,iside,*(subCVCellTopo_));
159
160 for (int inode = 0; inode < numNodesPerCell; inode++) {
161 for (int idim = 0; idim < spaceDim; idim++){
162 cubWeights(icell,inode,idim) = normals(inode,0,idim)*pow(2,spaceDim-1);
163 }
164 }
165
166 if (primaryCellTopo_->getKey()==shards::Hexahedron<8>::key)
167 {
168 // first set of side midpoints and normals (above) associated with
169 // primary cell edges 0-7 are obtained from side 1 of the
170 // eight control volumes
171
172 // second set of side midpoints and normals associated with
173 // primary cell edges 8-11 are obtained from side 5 of the
174 // first four control volumes.
175 iside = 5;
176 for (int i=0; i<numNodesPerSide; i++){
177 sideNodes(i) = subCVCellTopo_->getNodeMap(spaceDim-1,iside,i);
178 }
179 int numExtraSides = numEdgesPerCell - numNodesPerCell;
180 for (int icount=0; icount < numExtraSides; icount++){
181 int iedge = icount + numNodesPerCell;
182 for(int idim=0; idim < spaceDim; idim++){
183 Scalar midpt = 0.0;
184 for (int i=0; i<numNodesPerSide; i++){
185 midpt += subCVCoords(icell,icount,sideNodes(i),idim)/numNodesPerSide;
186 }
187 cubPoints(icell,iedge,idim) = midpt;
188 }
189 }
190
191 // Map side center to reference subcell
192 iside = 5;
194 sideCenterLocal,
195 spaceDim-1, iside, *(subCVCellTopo_));
196
197 // calculate Jacobian at side centers
198 Intrepid::CellTools<Scalar>::setJacobian(subCVsideJacobian, refSidePoints, cellCVCoords, *(subCVCellTopo_));
199
200 // Get subcontrol volume side normals
201 Intrepid::CellTools<Scalar>::getPhysicalSideNormals(normals,subCVsideJacobian,iside,*(subCVCellTopo_));
202
203 for (int icount = 0; icount < numExtraSides; icount++) {
204 int iedge = icount + numNodesPerCell;
205 for (int idim = 0; idim < spaceDim; idim++){
206 cubWeights(icell,iedge,idim) = normals(icount,0,idim)*pow(2,spaceDim-1);
207 }
208 }
209
210 } // end if Hex
211
212 if (primaryCellTopo_->getKey()==shards::Tetrahedron<4>::key)
213 {
214 // first set of side midpoints and normals associated with
215 // primary cell edges 0-2 are obtained from side 1 of the
216 // eight control volumes (above)
217
218 // second set of side midpoints and normals associated with
219 // primary cell edges 3-5 are obtained from side 5 of the
220 // first three control volumes.
221 iside = 5;
222 for (int i=0; i<numNodesPerSide; i++){
223 sideNodes(i) = subCVCellTopo_->getNodeMap(spaceDim-1,iside,i);
224 }
225 for (int icount=0; icount < 3; icount++){
226 int iedge = icount + 3;
227 for(int idim=0; idim < spaceDim; idim++){
228 Scalar midpt = 0.0;
229 for (int i=0; i<numNodesPerSide; i++){
230 midpt += subCVCoords(icell,icount,sideNodes(i),idim)/numNodesPerSide;
231 }
232 cubPoints(icell,iedge,idim) = midpt;
233 }
234 }
235
236 // Map side center to reference subcell
237 iside = 5;
239 sideCenterLocal,
240 spaceDim-1, iside, *(subCVCellTopo_));
241
242 // calculate Jacobian at side centers
243 Intrepid::CellTools<Scalar>::setJacobian(subCVsideJacobian, refSidePoints, cellCVCoords, *(subCVCellTopo_));
244
245 // Get subcontrol volume side normals
246 Intrepid::CellTools<Scalar>::getPhysicalSideNormals(normals,subCVsideJacobian,iside,*(subCVCellTopo_));
247
248 for (int icount = 0; icount < 3; icount++) {
249 int iedge = icount + 3;
250 for (int idim = 0; idim < spaceDim; idim++){
251 cubWeights(icell,iedge,idim) = normals(icount,0,idim)*pow(2,spaceDim-1);
252 }
253 }
254
255 }// if tetrahedron
256
257 } // end loop over cells
258
259} // end getCubature
260
261
262template<class Scalar, class ArrayPoint, class ArrayWeight>
264 return numPoints_;
265} // end getNumPoints
266
267template<class Scalar, class ArrayPoint, class ArrayWeight>
269 return cubDimension_;
270} // end getNumPoints
271
272template<class Scalar, class ArrayPoint, class ArrayWeight>
274 accuracy.assign(1,degree_);
275} // end getAccuracy
276
277} // end namespace Intrepid
278
279#endif
280
A stateless class for operations on cell data. Provides methods for:
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
static void getPhysicalSideNormals(ArraySideNormal &sideNormals, const ArrayJac &worksetJacobians, const int &worksetSideOrd, const shards::CellTopology &parentCell)
Computes non-normalized normal vectors to physical sides in a side workset .
static void getSubCVCoords(ArrayCVCoord &subCVcoords, const ArrayCellCoord &cellCoords, const shards::CellTopology &primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
CubatureControlVolumeSide(const Teuchos::RCP< const shards::CellTopology > &cellTopology)
int getDimension() const
Returns dimension of integration domain.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly on each triangle. The return vector ha...
int getNumPoints() const
Returns the number of cubature points.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights Method for reference space cubature - throws an exception.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....