Panzer  Version of the Day
Panzer_Dirichlet_Residual_EdgeBasis_impl.hpp
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 #ifndef PANZER_DIRICHLET_RESIDUAL_EDGEBASIS_IMPL_HPP
44 #define PANZER_DIRICHLET_RESIDUAL_EDGEBASIS_IMPL_HPP
45 
46 #include <cstddef>
47 #include <string>
48 #include <vector>
49 
50 #include "Intrepid2_Kernels.hpp"
51 #include "Intrepid2_CellTools.hpp"
52 #include "Intrepid2_OrientationTools.hpp"
53 
54 #include "Phalanx_Print.hpp"
55 
57 #include "Kokkos_ViewFactory.hpp"
58 
59 namespace panzer {
60 
61 //**********************************************************************
62 template<typename EvalT, typename Traits>
65  const Teuchos::ParameterList& p)
66 {
67  std::string residual_name = p.get<std::string>("Residual Name");
68 
69  basis = p.get<Teuchos::RCP<const panzer::PureBasis> >("Basis");
70  pointRule = p.get<Teuchos::RCP<const panzer::PointRule> >("Point Rule");
71 
72  std::string field_name = p.get<std::string>("DOF Name");
73  std::string dof_name = field_name+"_"+pointRule->getName();
74  std::string value_name = p.get<std::string>("Value Name");
75 
76  Teuchos::RCP<PHX::DataLayout> basis_layout = basis->functional;
77  Teuchos::RCP<PHX::DataLayout> vector_layout_dof = pointRule->dl_vector;
78  Teuchos::RCP<PHX::DataLayout> vector_layout_vector = basis->functional_grad;
79 
80  // some sanity checks
81  TEUCHOS_ASSERT(basis->isVectorBasis());
82  TEUCHOS_ASSERT(basis_layout->extent(0)==vector_layout_dof->extent(0));
83  TEUCHOS_ASSERT(basis_layout->extent(1)==vector_layout_dof->extent(1));
84  TEUCHOS_ASSERT(Teuchos::as<unsigned>(basis->dimension())==vector_layout_dof->extent(2));
85  TEUCHOS_ASSERT(vector_layout_vector->extent(0)==vector_layout_dof->extent(0));
86  TEUCHOS_ASSERT(vector_layout_vector->extent(1)==vector_layout_dof->extent(1));
87  TEUCHOS_ASSERT(vector_layout_vector->extent(2)==vector_layout_dof->extent(2));
88 
89  residual = PHX::MDField<ScalarT,Cell,BASIS>(residual_name, basis_layout);
90  dof = PHX::MDField<const ScalarT,Cell,Point,Dim>(dof_name, vector_layout_dof);
91  value = PHX::MDField<const ScalarT,Cell,Point,Dim>(value_name, vector_layout_vector);
92 
93  // setup all basis fields that are required
94 
95  // setup all fields to be evaluated and constructed
96  pointValues = PointValues2<double>(pointRule->getName()+"_",false);
97  pointValues.setupArrays(pointRule);
98 
99  // the field manager will allocate all of these field
100  this->addDependentField(pointValues.jac);
101 
102  this->addEvaluatedField(residual);
103  this->addDependentField(dof);
104  this->addDependentField(value);
105 
106  std::string n = "Dirichlet Residual Edge Basis Evaluator";
107  this->setName(n);
108 }
109 
110 //**********************************************************************
111 template<typename EvalT, typename Traits>
112 void
115  typename Traits::SetupData sd,
117 {
118  orientations = sd.orientations_;
119  this->utils.setFieldData(pointValues.jac,fm);
120 
121  const auto cellTopo = *basis->getCellTopology();
122  const int edgeDim = 1;
123  const int faceDim = 2;
124  if(cellTopo.getDimension() > edgeDim)
125  edgeParam = Intrepid2::RefSubcellParametrization<Kokkos::HostSpace>::get(edgeDim, cellTopo.getKey());
126 
127  if(cellTopo.getDimension() > faceDim)
128  faceParam = Intrepid2::RefSubcellParametrization<Kokkos::HostSpace>::get(faceDim, cellTopo.getKey());
129 }
130 
131 //**********************************************************************
132 template<typename EvalT, typename Traits>
133 void
136  typename Traits::EvalData workset)
137 {
138  const int numCells = workset.num_cells;
139  if(numCells <= 0)
140  return;
141  else {
142  residual.deep_copy(ScalarT(0.0));
143 
144  // dofs are already oriented but tangent directions are not oriented
145 
146  const int subcellDim = workset.subcell_dim;
147  const int subcellOrd = this->wda(workset).subcell_index;
148 
149  const auto cellTopo = *basis->getCellTopology();
150  const auto worksetJacobians = pointValues.jac.get_view();
151 
152  const int cellDim = cellTopo.getDimension();
153  const int edgeDim = 1;
154  const int faceDim = 2;
155 
156  auto intrepid_basis = basis->getIntrepid2Basis();
157  const WorksetDetails & details = workset;
158 
159  //const bool is_normalize = true;
160  auto work = Kokkos::create_mirror_view(Kokkos::createDynRankView(residual.get_static_view(),"work", 4, cellDim));
161 
162  // compute residual
163  auto residual_h = Kokkos::create_mirror_view(residual.get_static_view());
164  auto dof_h = Kokkos::create_mirror_view(dof.get_static_view());
165  auto value_h = Kokkos::create_mirror_view(value.get_static_view());
166  Kokkos::deep_copy(dof_h, dof.get_static_view());
167  Kokkos::deep_copy(value_h, value.get_static_view());
168  switch (subcellDim) {
169  case 1: { // 2D element Tri and Quad
170  if (intrepid_basis->getDofCount(1, subcellOrd)) {
171  auto ortEdgeTan = Kokkos::subview(work, 0, Kokkos::ALL());
172  auto phyEdgeTan = Kokkos::subview(work, 1, Kokkos::ALL());
173 
174  const int ndofsEdge = intrepid_basis->getDofCount(1, subcellOrd);
175  const int numEdges = cellTopo.getEdgeCount();
176  /* */ int edgeOrts[4] = {};
177 
178  for(index_t c=0;c<workset.num_cells;c++) {
179  orientations->at(details.cell_local_ids[c]).getEdgeOrientation(edgeOrts, numEdges);
180 
181  Intrepid2::Impl::OrientationTools::getRefSubcellTangents(work,
182  edgeParam,
183  cellTopo.getKey(edgeDim,subcellOrd),
184  subcellOrd,
185  edgeOrts[subcellOrd]);
186 
187  for (int i=0;i<ndofsEdge;++i) {
188  const int b = intrepid_basis->getDofOrdinal(1, subcellOrd, i);
189  auto J = Kokkos::create_mirror_view(Kokkos::subview(worksetJacobians, c, b, Kokkos::ALL(), Kokkos::ALL()));
190  Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan, J, ortEdgeTan);
191 
192  for(int d=0;d<cellDim;d++) {
193  residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyEdgeTan(d);
194  }
195  }
196  }
197  }
198  break;
199  }
200  case 2: { // 3D element Tet and Hex
201  const int numEdges = cellTopo.getEdgeCount();
202  const int numFaces = cellTopo.getFaceCount();
203 
204  {
205 
206  auto ortEdgeTan = Kokkos::subview(work, 0, Kokkos::ALL());
207  auto phyEdgeTan = Kokkos::subview(work, 1, Kokkos::ALL());
208 
209  const int numEdgesOfFace= cellTopo.getEdgeCount(2, subcellOrd);
210 
211  int edgeOrts[12] = {};
212  for(index_t c=0;c<workset.num_cells;c++) {
213  for (int i=0;i<numEdgesOfFace;++i) {
214 
215  const int edgeOrd = Intrepid2::Orientation::getEdgeOrdinalOfFace(i, subcellOrd, cellTopo);
216  const int b = edgeOrd;
217  orientations->at(details.cell_local_ids[c]).getEdgeOrientation(edgeOrts, numEdges);
218 
219  Intrepid2::Impl::OrientationTools::getRefSubcellTangents(work,
220  edgeParam,
221  cellTopo.getKey(edgeDim,edgeOrd),
222  edgeOrd,
223  edgeOrts[edgeOrd]);
224 
225  {
226  auto J = Kokkos::create_mirror_view(Kokkos::subview(worksetJacobians, c, b, Kokkos::ALL(), Kokkos::ALL()));
227  Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan, J, ortEdgeTan);
228 
229  for(int d=0;d<dof.extent_int(2);d++) {
230  residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyEdgeTan(d);
231  }
232  }
233  }
234  }
235  }
236 
237  if (intrepid_basis->getDofCount(2, subcellOrd)) {
238  auto ortFaceTanU = Kokkos::subview(work, 0, Kokkos::ALL());
239  auto ortFaceTanV = Kokkos::subview(work, 1, Kokkos::ALL());
240  auto phyFaceTanU = Kokkos::subview(work, 2, Kokkos::ALL());
241  auto phyFaceTanV = Kokkos::subview(work, 3, Kokkos::ALL());
242 
243  int faceOrts[6] = {};
244  for(index_t c=0;c<workset.num_cells;c++) {
245  orientations->at(details.cell_local_ids[c]).getFaceOrientation(faceOrts, numFaces);
246 
247  Intrepid2::Impl::OrientationTools::getRefSubcellTangents(work,
248  faceParam,
249  cellTopo.getKey(faceDim,subcellOrd),
250  subcellOrd,
251  faceOrts[subcellOrd]);
252 
253  for(int b=0;b<dof.extent_int(1);b++) {
254  auto J = Kokkos::create_mirror_view(Kokkos::subview(worksetJacobians, c, b, Kokkos::ALL(), Kokkos::ALL()));
255  Intrepid2::Kernels::Serial::matvec_product(phyFaceTanU, J, ortFaceTanU);
256  Intrepid2::Kernels::Serial::matvec_product(phyFaceTanV, J, ortFaceTanV);
257 
258  for(int d=0;d<dof.extent_int(2);d++) {
259  residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyFaceTanU(d);
260  residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyFaceTanV(d);
261  }
262  }
263  }
264  }
265 
266  break;
267  }
268  }
269  Kokkos::deep_copy(residual.get_static_view(), residual_h);
270  }
271 
272 }
273 
274 //**********************************************************************
275 
276 }
277 
278 #endif
Teuchos::RCP< const std::vector< Intrepid2::Orientation > > orientations_
int num_cells
DEPRECATED - use: numCells()
int subcell_dim
DEPRECATED - use: getSubcellDimension()
void setupArrays(const Teuchos::RCP< const panzer::PointRule > &pr)
Sizes/allocates memory for arrays.
Kokkos::DynRankView< typename InputArray::value_type, PHX::Device > createDynRankView(const InputArray &a, const std::string &name, const DimensionPack... dims)
Wrapper to simplify Panzer use of Sacado ViewFactory.
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)