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...
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)