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 
52 namespace Intrepid{
53 
54 template<class Scalar, class ArrayPoint, class ArrayWeight>
55 CubatureControlVolumeSide<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 
75 template<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 
83 template<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 
262 template<class Scalar, class ArrayPoint, class ArrayWeight>
264  return numPoints_;
265 } // end getNumPoints
266 
267 template<class Scalar, class ArrayPoint, class ArrayWeight>
269  return cubDimension_;
270 } // end getNumPoints
271 
272 template<class Scalar, class ArrayPoint, class ArrayWeight>
274  accuracy.assign(1,degree_);
275 } // end getAccuracy
276 
277 } // end namespace Intrepid
278 
279 #endif
280 
281 
282 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
283 #ifdef __GNUC__
284 #warning "The Intrepid package is deprecated"
285 #endif
286 #endif
287 
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.
CubatureControlVolumeSide(const Teuchos::RCP< const shards::CellTopology > &cellTopology)
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 .
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...
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights Method for reference space cubature - throws an exception...
static void getSubCVCoords(ArrayCVCoord &subCVcoords, const ArrayCellCoord &cellCoords, const shards::CellTopology &primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
A stateless class for operations on cell data. Provides methods for:
int getNumPoints() const
Returns the number of cubature points.