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
52namespace Intrepid{
53
54template<class Scalar, class ArrayPoint, class ArrayWeight>
55CubatureControlVolumeBoundary<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
78template<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
86template<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
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
213template<class Scalar, class ArrayPoint, class ArrayWeight>
215 return numPoints_;
216} // end getNumPoints
217
218template<class Scalar, class ArrayPoint, class ArrayWeight>
220 return cubDimension_;
221} // end getDimension
222
223template<class Scalar, class ArrayPoint, class ArrayWeight>
225 accuracy.assign(1,degree_);
226} // end getAccuracy
227
228} // end namespace Intrepid
229
230#endif
231
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 setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void getSubCVCoords(ArrayCVCoord &subCVcoords, const ArrayCellCoord &cellCoords, const shards::CellTopology &primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
CubatureControlVolumeBoundary(const Teuchos::RCP< const shards::CellTopology > &cellTopology, int cellSide=0)
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1.
int getDimension() const
Returns dimension of integration domain.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights Method for reference space cubature - throws an exception.
int getNumPoints() const
Returns the number of cubature points.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....