Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_IntegrationRule.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
43
45
46#include "Teuchos_ArrayRCP.hpp"
47#include "Teuchos_Assert.hpp"
48#include "Phalanx_DataLayout_MDALayout.hpp"
49#include "Intrepid2_DefaultCubatureFactory.hpp"
50#include "Intrepid2_CubatureControlVolume.hpp"
51#include "Intrepid2_CubatureControlVolumeSide.hpp"
52#include "Intrepid2_CubatureControlVolumeBoundary.hpp"
53#include "Panzer_Dimension.hpp"
54#include "Panzer_CellData.hpp"
55
56
58IntegrationRule(int in_cubature_degree, const panzer::CellData& cell_data)
60{
61 if(cell_data.isSide()){
62 IntegrationDescriptor::setup(in_cubature_degree, IntegrationDescriptor::SIDE, cell_data.side());
63 } else {
65 }
66 setup(in_cubature_degree,cell_data);
67}
68
70IntegrationRule(const panzer::CellData& cell_data, const std::string & in_cv_type)
72{
73 // Cubature orders are only used for indexing so we make them large enough not to interfere with other rules.
74 // TODO: This requirement (on arbitrary cubature order) will be dropped with the new Workset design (using descriptors to find integration rules)
75 if(in_cv_type == "volume"){
76 TEUCHOS_TEST_FOR_EXCEPT_MSG(cell_data.isSide(),
77 "IntegrationRule::IntegrationRule : Control Volume 'volume' type requested, but CellData is setup for sides.");
79 } else if(in_cv_type == "side"){
81 } else if(in_cv_type == "boundary"){
82 TEUCHOS_TEST_FOR_EXCEPT_MSG(not cell_data.isSide(),
83 "IntegrationRule::IntegrationRule : Control Volume 'boundary' type requested, but CellData is not setup for sides.");
85 } else {
86 TEUCHOS_ASSERT(false);
87 }
88 setup_cv(cell_data,in_cv_type);
89}
90
93 const Teuchos::RCP<const shards::CellTopology> & cell_topology,
94 const int num_cells,
95 const int num_faces)
96 : PointRule(), IntegrationDescriptor(description)
97{
98
99 TEUCHOS_ASSERT(description.getType() != panzer::IntegrationDescriptor::NONE);
100
101 cubature_degree = description.getOrder();
102 cv_type = "none";
103
104 panzer::CellData cell_data;
105 if(isSide()){
106 cell_data = panzer::CellData(num_cells, getSide(), cell_topology);
107 } else {
108 cell_data = panzer::CellData(num_cells, cell_topology);
109 }
110
112 setup(getOrder(), cell_data);
113 } else if(description.getType() == panzer::IntegrationDescriptor::SIDE){
114 setup(getOrder(), cell_data);
115 } else if(description.getType() == panzer::IntegrationDescriptor::SURFACE){
116 TEUCHOS_ASSERT(num_faces!=-1);
117 setup_surface(cell_topology, num_cells, num_faces);
118 } else if(description.getType() == panzer::IntegrationDescriptor::CV_VOLUME){
119 setup_cv(cell_data, "volume");
120 } else if(description.getType() == panzer::IntegrationDescriptor::CV_SIDE){
121 setup_cv(cell_data, "side");
122 } else if(description.getType() == panzer::IntegrationDescriptor::CV_BOUNDARY){
123 setup_cv(cell_data, "boundary");
124 } else {
125 TEUCHOS_ASSERT(false);
126 }
127
128}
129
130void panzer::IntegrationRule::setup(int in_cubature_degree, const panzer::CellData& cell_data)
131{
132 cubature_degree = in_cubature_degree;
133 cv_type = "none";
134 int spatialDimension = cell_data.baseCellDimension();
135
136 std::stringstream ss;
137 ss << "CubaturePoints (Degree=" << cubature_degree;
138
139 // Intrepid2 does not support a quadrature on a 0-dimensional object
140 // (which doesn't make much sense anyway) to work around this we
141 // will adjust the integration rule manually
142 if(cell_data.isSide() && spatialDimension==1) {
143 ss << ",side)";
144 PointRule::setup(ss.str(),1,cell_data);
145
146 return;
147 }
148
149 const shards::CellTopology & topo = *cell_data.getCellTopology();
150 Teuchos::RCP<shards::CellTopology> sideTopo = getSideTopology(cell_data);
151
152 Intrepid2::DefaultCubatureFactory cubature_factory;
153 Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double>> intrepid_cubature;
154
155 // get side topology
156 if (Teuchos::is_null(sideTopo)) {
157 ss << ",volume)";
158 intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(topo, cubature_degree);
159 }
160 else {
161 ss << ",side)";
162 intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*sideTopo, cubature_degree);
163 }
164
165 PointRule::setup(ss.str(),intrepid_cubature->getNumPoints(),cell_data);
166}
167
168void panzer::IntegrationRule::setup_surface(const Teuchos::RCP<const shards::CellTopology> & cell_topology, const int num_cells, const int num_faces)
169{
170
171 const int cell_dim = cell_topology->getDimension();
172 const int subcell_dim = cell_dim-1;
173 const int num_faces_per_cell = cell_topology->getSubcellCount(subcell_dim);
174
175 panzer::CellData cell_data(num_cells, cell_topology);
176
177 std::string point_rule_name;
178 {
179 std::stringstream ss;
180 ss << "CubaturePoints (Degree=" << getOrder() << ",surface)";
181 point_rule_name = ss.str();
182 }
183
184 // We can skip some steps for 1D
185 if(cell_dim == 1){
186 const int num_points_per_cell = num_faces_per_cell;
187 const int num_points_per_face = 1;
188 PointRule::setup(point_rule_name, num_cells, num_points_per_cell, num_faces, num_points_per_face, cell_topology);
189 _point_offsets.resize(3,0);
190 _point_offsets[0] = 0;
191 _point_offsets[1] = num_points_per_face;
192 _point_offsets[2] = _point_offsets[1]+num_points_per_face;
193 return;
194 }
195
196 Intrepid2::DefaultCubatureFactory cubature_factory;
197
198 _point_offsets.resize(num_faces_per_cell+1,0);
199 int test_face_size = -1;
200 for(int subcell_index=0; subcell_index<num_faces_per_cell; ++subcell_index){
201 Teuchos::RCP<shards::CellTopology> face_topology = Teuchos::rcp(new shards::CellTopology(cell_topology->getCellTopologyData(subcell_dim,subcell_index)));
202 const auto & intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*face_topology, getOrder());
203 const int num_face_points = intrepid_cubature->getNumPoints();
204 _point_offsets[subcell_index+1] = _point_offsets[subcell_index] + num_face_points;
205
206 // Right now we only support each face having the same number of points
207 if(test_face_size==-1){
208 test_face_size = num_face_points;
209 } else {
210 TEUCHOS_ASSERT(num_face_points == test_face_size);
211 }
212 }
213
214 const int num_points_per_cell = _point_offsets.back();
215 const int num_points_per_face = _point_offsets[1];
216
217 PointRule::setup(point_rule_name, num_cells, num_points_per_cell, num_faces, num_points_per_face, cell_topology);
218
219}
220
221void panzer::IntegrationRule::setup_cv(const panzer::CellData& cell_data, std::string in_cv_type)
222{
223 // set cubature degree to arbitrary constant for indexing
224 cubature_degree = 1;
225 cv_type = in_cv_type;
226 if (cv_type == "volume") {
227 cubature_degree = 75;
228 }
229 if (cv_type == "side") {
230 cubature_degree = 85;
231 }
232 if (cv_type == "boundary") {
233 cubature_degree = 95;
234 }
235
236 //int spatialDimension = cell_data.baseCellDimension();
237
238 std::stringstream ss;
239 ss << "CubaturePoints ControlVol (Index=" << cubature_degree;
240
241 const shards::CellTopology & topo = *cell_data.getCellTopology();
242
243 Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double> > intrepid_cubature;
244
245 int tmp_num_points = 0;
246 if (cv_type == "volume") {
247 ss << ",volume)";
248 intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(topo));
249 tmp_num_points = intrepid_cubature->getNumPoints();
250 }
251 else if (cv_type == "side") {
252 ss << ",side)";
253 intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(topo));
254 tmp_num_points = intrepid_cubature->getNumPoints();
255 }
256 else if (cv_type == "boundary") {
257 ss << ",boundary)";
258 intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(topo,cell_data.side()));
259 tmp_num_points = intrepid_cubature->getNumPoints();
260 }
261
262 PointRule::setup(ss.str(),tmp_num_points,cell_data);
263}
264
266{ return cubature_degree; }
267
268
269int panzer::IntegrationRule::getPointOffset(const int subcell_index) const
270{
271 // Need to make sure this is a surface integrator
272 TEUCHOS_ASSERT(getType() == SURFACE);
273 return _point_offsets[subcell_index];
274}
275
276
277void panzer::IntegrationRule::print(std::ostream & os)
278{
279 os << "IntegrationRule ( "
280 << "Name = " << getName()
281 << ", Degree = " << cubature_degree
282 << ", Dimension = " << spatial_dimension
283 << ", Workset Size = " << workset_size
284 << ", Num Points = " << num_points
285 << ", Side = " << side
286 << " )";
287}
288
289void panzer::IntegrationRule::referenceCoordinates(Kokkos::DynRankView<double,PHX::Device> & cub_points)
290{
291 // build an interpid cubature rule
292 Teuchos::RCP< Intrepid2::Cubature<PHX::Device::execution_space,double,double> > intrepid_cubature;
293 Intrepid2::DefaultCubatureFactory cubature_factory;
294
295 if (!isSide())
296 intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*(topology),cubature_degree);
297 else
298 intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*(side_topology),cubature_degree);
299
300 int num_ip = intrepid_cubature->getNumPoints();
301 Kokkos::DynRankView<double,PHX::Device> cub_weights("cub_weights",num_ip);
302
303 // now compute weights (and throw them out) as well as reference points
304 if (!isSide()) {
305 cub_points = Kokkos::DynRankView<double,PHX::Device>("cub_points", num_ip, topology->getDimension());
306 intrepid_cubature->getCubature(cub_points, cub_weights);
307 }
308 else {
309 cub_points = Kokkos::DynRankView<double,PHX::Device>("cub_points", num_ip, side_topology->getDimension());
310 intrepid_cubature->getCubature(cub_points, cub_weights);
311 }
312}
Data for determining cell topology and dimensionality.
int baseCellDimension() const
Dimension of the base cell. NOT the dimension of the local side, even if the side() method returns tr...
bool isSide() const
Teuchos::RCP< const shards::CellTopology > getCellTopology() const
Get CellTopology for the base cell.
@ VOLUME
No integral specified - default state.
@ CV_VOLUME
Integral over a specific side of cells (side must be set)
@ CV_BOUNDARY
Control volume side integral.
@ SIDE
Integral over all sides of cells (closed surface integral)
const int & getOrder() const
Get order of integrator.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
const int & getType() const
Get type of integrator.
void setup(const int cubature_order, const int integration_type, const int side=-1)
Setup function.
IntegrationRule(int cubature_degree, const panzer::CellData &cell_data)
if side = -1 then we use the cell volume integration rule.
void setup_cv(const panzer::CellData &cell_data, std::string cv_type)
int getPointOffset(const int subcell_index) const
Returns the integration point offset for a given subcell_index (i.e. local face index)
void setup_surface(const Teuchos::RCP< const shards::CellTopology > &cell_topology, const int num_cells, const int num_faces)
Setup a surface integration.
virtual void print(std::ostream &os)
print information about the integration rule
void setup(int cubature_degree, const panzer::CellData &cell_data)
void referenceCoordinates(Kokkos::DynRankView< double, PHX::Device > &container)
Construct an array containing the reference coordinates.
int order() const
Returns the order of integration (cubature degree in intrepid lingo)
void setup(const std::string &ptName, int np, const panzer::CellData &cell_data)