Intrepid
Intrepid_CubatureControlVolumeBoundaryDef.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 
49 #ifndef INTREPID_CUBATURE_CONTROLVOLUMEBOUNDARYDEF_HPP
50 #define INTREPID_CUBATURE_CONTROLVOLUMEBOUNDARYDEF_HPP
51 
52 namespace Intrepid{
53 
54 template<class Scalar, class ArrayPoint, class ArrayWeight>
55 CubatureControlVolumeBoundary<Scalar,ArrayPoint,ArrayWeight>::CubatureControlVolumeBoundary(const Teuchos::RCP<const shards::CellTopology> & cellTopology, int cellSide)
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  cubDimension_ = primaryCellTopo_->getDimension();
70 
71  sideIndex_ = cellSide;
72 
73  // one control volume boundary cubature point per subcell node (for now)
74  numPoints_ = primaryCellTopo_->getNodeCount(primaryCellTopo_->getDimension()-1,cellSide);
75 
76 }
77 
78 template<class Scalar, class ArrayPoint, class ArrayWeight>
80  ArrayWeight& cubWeights) const
81 {
82  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
83  ">>> ERROR (CubatureControlVolumeBoundary): Cubature defined in physical space calling method for reference space cubature.");
84 }
85 
86 template<class Scalar, class ArrayPoint, class ArrayWeight>
88  ArrayWeight& cubWeights,
89  ArrayPoint& cellCoords) const
90 {
91  // get array dimensions
92  int numCells = cellCoords.dimension(0);
93  int numNodesPerCell = cellCoords.dimension(1);
94  int spaceDim = cellCoords.dimension(2);
95  int numNodesPerSubCV = subCVCellTopo_->getNodeCount();
96 
97  // get sub-control volume coordinates (one sub-control volume per node of primary cell)
98  Intrepid::FieldContainer<Scalar> subCVCoords(numCells,numNodesPerCell,numNodesPerSubCV,spaceDim);
99  Intrepid::CellTools<Scalar>::getSubCVCoords(subCVCoords,cellCoords,*(primaryCellTopo_));
100 
101  // define subcontrol volume side index corresponding to primary cell side
102  int numPrimarySideNodes = primaryCellTopo_->getNodeCount(spaceDim-1,sideIndex_);
103  int numCVSideNodes = subCVCellTopo_->getNodeCount(spaceDim-1,0);
104  int numPrimarySides = primaryCellTopo_->getSubcellCount(spaceDim-1);
105  Intrepid::FieldContainer<int> CVSideonBoundary(numPrimarySides,numCVSideNodes);
106 
107  switch(primaryCellTopo_->getKey() ) {
108 
109  case shards::Triangle<3>::key:
110  case shards::Quadrilateral<4>::key:
111 
112  for (int iside=0; iside<numPrimarySides; iside++) {
113  CVSideonBoundary(iside,0) = 0; CVSideonBoundary(iside,1) = 3;
114  }
115 
116  break;
117 
118  case shards::Hexahedron<8>::key:
119 
120  // sides 0-3
121  for (int iside=0; iside<4; iside++) {
122  CVSideonBoundary(iside,0) = 0; CVSideonBoundary(iside,1) = 3;
123  CVSideonBoundary(iside,2) = 3; CVSideonBoundary(iside,3) = 0;
124  }
125  // side 4
126  CVSideonBoundary(4,0) = 4; CVSideonBoundary(4,1) = 4;
127  CVSideonBoundary(4,2) = 4; CVSideonBoundary(4,3) = 4;
128 
129  // side 5
130  CVSideonBoundary(5,0) = 5; CVSideonBoundary(5,1) = 5;
131  CVSideonBoundary(5,2) = 5; CVSideonBoundary(5,3) = 5;
132 
133  break;
134 
135  case shards::Tetrahedron<4>::key:
136 
137  CVSideonBoundary(0,0) = 0; CVSideonBoundary(0,1) = 3; CVSideonBoundary(0,2) = 0;
138  CVSideonBoundary(1,0) = 0; CVSideonBoundary(1,1) = 3; CVSideonBoundary(1,2) = 3;
139  CVSideonBoundary(2,0) = 3; CVSideonBoundary(2,1) = 4; CVSideonBoundary(2,2) = 0;
140  CVSideonBoundary(3,0) = 4; CVSideonBoundary(3,1) = 4; CVSideonBoundary(3,2) = 4;
141 
142  break;
143 
144  default:
145  TEUCHOS_TEST_FOR_EXCEPTION( true, std::invalid_argument,
146  ">>> ERROR (CubatureControlVolumeBoundary: invalid cell topology.");
147  } // cell key
148 
149  Intrepid::FieldContainer<double> sideCenterLocal(1,spaceDim-1);
150  for (int idim = 0; idim < spaceDim-1; idim++){
151  sideCenterLocal(0,idim) = 0.0;
152  }
153 
154  // get side cubature points
155  for (int icell = 0; icell < numCells; icell++)
156  {
157 
158  for (int inode=0; inode<numPrimarySideNodes; inode++){
159 
160  int cvind = primaryCellTopo_->getNodeMap(spaceDim-1,sideIndex_,inode);
161  int cvside = CVSideonBoundary(sideIndex_,inode);
162 
163  Intrepid::FieldContainer<double> cubpoint(spaceDim);
164  for (int idim=0; idim<spaceDim; idim++) {
165  for (int icvnode=0; icvnode<numCVSideNodes; icvnode++) {
166  int cvnode = subCVCellTopo_->getNodeMap(spaceDim-1,cvside,icvnode);
167  cubpoint(idim) = cubpoint(idim) + subCVCoords(icell,cvind,cvnode,idim);
168  }
169  cubPoints(icell,inode,idim) = cubpoint(idim)/numCVSideNodes;
170  }
171 
172  // map side center to reference subcell
173  Intrepid::FieldContainer<Scalar> refSidePoints(1,spaceDim);
175  sideCenterLocal,
176  spaceDim-1, cvside, *(subCVCellTopo_));
177 
178  // array of sub-control volume coordinates
179  Intrepid::FieldContainer<Scalar> cellCVCoords(1, numNodesPerSubCV, spaceDim);
180  for (int icvnode = 0; icvnode < numNodesPerSubCV; icvnode++){
181  for (int idim = 0; idim < spaceDim; idim++){
182  cellCVCoords(0,icvnode,idim) = subCVCoords(icell,cvind,icvnode,idim);
183  }
184  }
185 
186  // calculate Jacobian at side centers
187  Intrepid::FieldContainer<Scalar> subCVsideJacobian(1, 1, spaceDim, spaceDim);
188  Intrepid::FieldContainer<Scalar> subCVsideJacobianDet(1, 1);
189  Intrepid::CellTools<Scalar>::setJacobian(subCVsideJacobian, refSidePoints, cellCVCoords, *(subCVCellTopo_));
190  Intrepid::CellTools<Scalar>::setJacobianDet(subCVsideJacobianDet, subCVsideJacobian);
191 
192  // calculate Jacobian at side centers
193  Intrepid::FieldContainer<Scalar> measure(1, 1);
194  Intrepid::FieldContainer<Scalar> weights(1, 1);
195  if (spaceDim == 3){
196  weights(0,0) = 4.0;
197  Intrepid::FunctionSpaceTools::computeFaceMeasure<Scalar>(measure,subCVsideJacobian,weights,cvside,*(subCVCellTopo_));
198  }
199  else if (spaceDim == 2){
200  weights(0,0) = 2.0;
201  Intrepid::FunctionSpaceTools::computeEdgeMeasure<Scalar>(measure,subCVsideJacobian,weights,cvside,*(subCVCellTopo_));
202  }
203 
204  cubWeights(icell,inode) = measure(0,0);
205 
206  } // end loop over primary side nodes
207 
208  } // end cell loop
209 
210 } // end getCubature
211 
212 
213 template<class Scalar, class ArrayPoint, class ArrayWeight>
215  return numPoints_;
216 } // end getNumPoints
217 
218 template<class Scalar, class ArrayPoint, class ArrayWeight>
220  return cubDimension_;
221 } // end getDimension
222 
223 template<class Scalar, class ArrayPoint, class ArrayWeight>
225  accuracy.assign(1,degree_);
226 } // end getAccuracy
227 
228 } // end namespace Intrepid
229 
230 #endif
231 
232 
233 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
234 #ifdef __GNUC__
235 #warning "The Intrepid package is deprecated"
236 #endif
237 #endif
238 
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.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1...
CubatureControlVolumeBoundary(const Teuchos::RCP< const shards::CellTopology > &cellTopology, int cellSide=0)
int getDimension() const
Returns dimension of integration domain.
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...
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:
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F...