Panzer  Version of the Day
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_Assert.hpp"
47 #include "Phalanx_DataLayout_MDALayout.hpp"
48 #include "Intrepid2_DefaultCubatureFactory.hpp"
49 #include "Intrepid2_CubatureControlVolume.hpp"
50 #include "Intrepid2_CubatureControlVolumeSide.hpp"
51 #include "Intrepid2_CubatureControlVolumeBoundary.hpp"
52 #include "Panzer_Dimension.hpp"
53 #include "Panzer_CellData.hpp"
54 
56 IntegrationRule(int in_cubature_degree, const panzer::CellData& cell_data)
57  : PointRule()
58 {
59  setup(in_cubature_degree,cell_data);
60 }
61 
63 IntegrationRule(const panzer::CellData& cell_data, std::string in_cv_type)
64  : PointRule()
65 {
66  setup_cv(cell_data,in_cv_type);
67 }
68 
69 void panzer::IntegrationRule::setup(int in_cubature_degree, const panzer::CellData& cell_data)
70 {
71  cubature_degree = in_cubature_degree;
72  cv_type = "none";
73  int spatialDimension = cell_data.baseCellDimension();
74 
75  std::stringstream ss;
76  ss << "CubaturePoints (Degree=" << cubature_degree;
77 
78  // Intrepid2 does not support a quadrature on a 0-dimensional object
79  // (which doesn't make much sense anyway) to work around this we
80  // will adjust the integration rule manually
81  if(cell_data.isSide() && spatialDimension==1) {
82  ss << ",side)";
83  PointRule::setup(ss.str(),1,cell_data);
84 
85  return;
86  }
87 
89  Teuchos::RCP<shards::CellTopology> sideTopo = getSideTopology(cell_data);
90 
91  Intrepid2::DefaultCubatureFactory<double,Kokkos::DynRankView<double,PHX::Device> > cubature_factory;
93 
94  // get side topology
95  if (Teuchos::is_null(sideTopo)) {
96  ss << ",volume)";
97  intrepid_cubature = cubature_factory.create(*topo, cubature_degree);
98  }
99  else {
100  ss << ",side)";
101  intrepid_cubature = cubature_factory.create(*sideTopo, cubature_degree);
102  }
103 
104  PointRule::setup(ss.str(),intrepid_cubature->getNumPoints(),cell_data);
105 }
106 
107 void panzer::IntegrationRule::setup_cv(const panzer::CellData& cell_data, std::string in_cv_type)
108 {
109  // set cubature degree to arbitrary constant for indexing
110  cubature_degree = 1;
111  cv_type = in_cv_type;
112  if (cv_type == "volume") {
113  cubature_degree = 75;
114  }
115  if (cv_type == "side") {
116  cubature_degree = 85;
117  }
118  if (cv_type == "boundary") {
119  cubature_degree = 95;
120  }
121 
122  //int spatialDimension = cell_data.baseCellDimension();
123 
124  std::stringstream ss;
125  ss << "CubaturePoints ControlVol (Index=" << cubature_degree;
126 
128 
130 
131  int num_points(0);
132  if (cv_type == "volume") {
133  ss << ",volume)";
134  intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolume<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device> >(topo));
135  num_points = intrepid_cubature->getNumPoints();
136  }
137  else if (cv_type == "side") {
138  ss << ",side)";
139  intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device> >(topo));
140  num_points = intrepid_cubature->getNumPoints();
141  }
142  else if (cv_type == "boundary") {
143  ss << ",boundary)";
144  intrepid_cubature = Teuchos::rcp(new
145  Intrepid2::CubatureControlVolumeBoundary<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device> >(topo,cell_data.side()));
146  num_points = intrepid_cubature->getNumPoints();
147  }
148 
149  PointRule::setup(ss.str(),num_points,cell_data);
150 }
151 
153 { return cubature_degree; }
154 
155 void panzer::IntegrationRule::print(std::ostream & os)
156 {
157  os << "IntegrationRule ( "
158  << "Name = " << getName()
159  << ", Degree = " << cubature_degree
160  << ", Dimension = " << spatial_dimension
161  << ", Workset Size = " << workset_size
162  << ", Num Points = " << num_points
163  << ", Side = " << side
164  << " )";
165 }
166 
167 void panzer::IntegrationRule::referenceCoordinates(Kokkos::DynRankView<double,PHX::Device> & cub_points)
168 {
169  // build an interpid cubature rule
171  Intrepid2::DefaultCubatureFactory<double,Kokkos::DynRankView<double,PHX::Device> > cubature_factory;
172 
173  if (!isSide())
174  intrepid_cubature = cubature_factory.create(*(topology),cubature_degree);
175  else
176  intrepid_cubature = cubature_factory.create(*(side_topology),cubature_degree);
177 
178  int num_ip = intrepid_cubature->getNumPoints();
179  Kokkos::DynRankView<double,PHX::Device> cub_weights("cub_weights",num_ip);
180 
181  // now compute weights (and throw them out) as well as reference points
182  if (!isSide()) {
183  cub_points = Kokkos::DynRankView<double,PHX::Device>("cub_points", num_ip, topology->getDimension());
184  intrepid_cubature->getCubature(cub_points, cub_weights);
185  }
186  else {
187  cub_points = Kokkos::DynRankView<double,PHX::Device>("cub_points", num_ip, side_topology->getDimension());
188  intrepid_cubature->getCubature(cub_points, cub_weights);
189  }
190 }
void setup(const std::string &ptName, int np, const panzer::CellData &cell_data)
bool is_null(const std::shared_ptr< T > &p)
void referenceCoordinates(Kokkos::DynRankView< double, PHX::Device > &container)
Construct an array containing the reference coordinates.
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.
void setup(int cubature_degree, const panzer::CellData &cell_data)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Data for determining cell topology and dimensionality.
bool isSide() const
void setup_cv(const panzer::CellData &cell_data, std::string cv_type)
virtual void print(std::ostream &os)
print information about the integration rule
int order() const
Returns the order of integration (cubature degree in intrepid lingo)
IntegrationRule(int cubature_degree, const panzer::CellData &cell_data)
if side = -1 then we use the cell volume integration rule.