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 
52 namespace Intrepid{
53 
54 template<class Scalar, class ArrayPoint, class ArrayWeight>
55 CubatureControlVolume<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 
75 template<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 
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  // 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 
149 template<class Scalar, class ArrayPoint, class ArrayWeight>
151  return numPoints_;
152 } // end getNumPoints
153 
154 template<class Scalar, class ArrayPoint, class ArrayWeight>
156  return cubDimension_;
157 } // end getNumPoints
158 
159 template<class Scalar, class ArrayPoint, class ArrayWeight>
161  accuracy.assign(1,degree_);
162 } // end getAccuracy
163 
164 } // end namespace Intrepid
165 
166 #endif
167 
168 
169 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
170 #ifdef __GNUC__
171 #warning "The Intrepid package is deprecated"
172 #endif
173 #endif
174 
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)
int getDimension() const
Returns dimension of integration domain.
int getNumPoints() const
Returns the number of cubature points.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1...
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.
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...