Intrepid
Intrepid_CubatureControlVolumeDef.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_CONTROLVOLUMEDEF_HPP
50#define INTREPID_CUBATURE_CONTROLVOLUMEDEF_HPP
51
52namespace Intrepid{
53
54template<class Scalar, class ArrayPoint, class ArrayWeight>
55CubatureControlVolume<Scalar,ArrayPoint,ArrayWeight>::CubatureControlVolume(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 subCVCellTopo_ = Teuchos::rcp(new shards::CellTopology(&myCellData));
65
66 degree_ = 1;
67
68 // one control volume cubature point per primary cell node
69 numPoints_ = primaryCellTopo_->getNodeCount();
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 (CubatureControlVolume): 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 // Integration points and weights for calculating sub-control volumes
100 int subcvCubDegree = 2;
101 Teuchos::RCP<Intrepid::Cubature<double,Intrepid::FieldContainer<double> > > subCVCubature;
102 subCVCubature = subCVCubFactory.create(*(subCVCellTopo_), subcvCubDegree);
103
104 int subcvCubDim = subCVCubature -> getDimension();
105 int numSubcvCubPoints = subCVCubature -> getNumPoints();
106
107 // Get numerical integration points and weights
108 Intrepid::FieldContainer<double> subcvCubPoints (numSubcvCubPoints, subcvCubDim);
109 Intrepid::FieldContainer<double> subcvCubWeights(numSubcvCubPoints);
110
111 subCVCubature -> getCubature(subcvCubPoints, subcvCubWeights);
112
113 // Loop over cells
114 for (int icell = 0; icell < numCells; icell++){
115
116 // get sub-control volume centers (integration points)
117 Intrepid::FieldContainer<Scalar> subCVCenter(numNodesPerCell,1,spaceDim);
118 Intrepid::FieldContainer<Scalar> cellCVCoords(numNodesPerCell,numNodesPerSubCV,spaceDim);
119 for (int isubcv = 0; isubcv < numNodesPerCell; isubcv++){
120 for (int idim = 0; idim < spaceDim; idim++){
121 for (int inode = 0; inode < numNodesPerSubCV; inode++){
122 subCVCenter(isubcv,0,idim) += subCVCoords(icell,isubcv,inode,idim)/numNodesPerSubCV;
123 cellCVCoords(isubcv,inode,idim) = subCVCoords(icell,isubcv,inode,idim);
124 }
125 cubPoints(icell,isubcv,idim) = subCVCenter(isubcv,0,idim);
126 }
127 }
128
129 // calculate Jacobian and determinant for each subCV quadrature point
130 Intrepid::FieldContainer<Scalar> subCVJacobian(numNodesPerCell, numSubcvCubPoints, spaceDim, spaceDim);
131 Intrepid::FieldContainer<Scalar> subCVJacobDet(numNodesPerCell, numSubcvCubPoints);
132 Intrepid::CellTools<Scalar>::setJacobian(subCVJacobian, subcvCubPoints, cellCVCoords, *(subCVCellTopo_));
133 Intrepid::CellTools<Scalar>::setJacobianDet(subCVJacobDet, subCVJacobian );
134
135 // fill array with sub control volumes (the sub control volume cell measure)
136 for (int inode = 0; inode < numNodesPerCell; inode++){
137 Scalar vol = 0;
138 for (int ipt = 0; ipt < numSubcvCubPoints; ipt++){
139 vol += subcvCubWeights(ipt)*subCVJacobDet(inode,ipt);
140 }
141 cubWeights(icell,inode) = vol;
142 }
143
144 } // end cell loop
145
146} // end getCubature
147
148
149template<class Scalar, class ArrayPoint, class ArrayWeight>
151 return numPoints_;
152} // end getNumPoints
153
154template<class Scalar, class ArrayPoint, class ArrayWeight>
156 return cubDimension_;
157} // end getNumPoints
158
159template<class Scalar, class ArrayPoint, class ArrayWeight>
161 accuracy.assign(1,degree_);
162} // end getAccuracy
163
164} // end namespace Intrepid
165
166#endif
167
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.
static void getSubCVCoords(ArrayCVCoord &subCVcoords, const ArrayCellCoord &cellCoords, const shards::CellTopology &primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
int getNumPoints() const
Returns the number of cubature points.
int getDimension() const
Returns dimension of integration domain.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights Method for reference space cubature - throws an exception.
CubatureControlVolume(const Teuchos::RCP< const shards::CellTopology > &cellTopology)
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....