Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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_TypeStrings.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  constJac_ = pointValues.jac;
101  this->addDependentField(constJac_);
102 
103  this->addEvaluatedField(residual);
104  this->addDependentField(dof);
105  this->addDependentField(value);
106 
107  std::string n = "Dirichlet Residual Edge Basis Evaluator";
108  this->setName(n);
109 }
110 
111 //**********************************************************************
112 template<typename EvalT, typename Traits>
113 void
116  typename Traits::SetupData sd,
118 {
119  orientations = sd.orientations_;
120  this->utils.setFieldData(pointValues.jac,fm);
121 }
122 
123 //**********************************************************************
124 template<typename EvalT, typename Traits>
125 void
128  typename Traits::EvalData workset)
129 {
130  const int numCells = workset.num_cells;
131  if(numCells <= 0)
132  return;
133  else {
134  residual.deep_copy(ScalarT(0.0));
135 
136  // dofs are already oriented but tangent directions are not oriented
137 
138  const int subcellDim = workset.subcell_dim;
139  const int subcellOrd = this->wda(workset).subcell_index;
140 
141  const auto cellTopo = *basis->getCellTopology();
142  const auto worksetJacobians = pointValues.jac.get_view();
143 
144  const int cellDim = cellTopo.getDimension();
145 
146  auto intrepid_basis = basis->getIntrepid2Basis();
147  const WorksetDetails & details = workset;
148 
149  const bool is_normalize = true;
150  auto work = Kokkos::createDynRankView(residual.get_static_view(),"work", 4, cellDim);
151 
152  // compute residual
153  switch (subcellDim) {
154  case 1: { // 2D element Tri and Quad
155  if (intrepid_basis->getDofCount(1, subcellOrd)) {
156  auto phyEdgeTan = Kokkos::subview(work, 0, Kokkos::ALL());
157  auto ortEdgeTan = Kokkos::subview(work, 1, Kokkos::ALL());
158 
159  const int ndofsEdge = intrepid_basis->getDofCount(1, subcellOrd);
160  const int numEdges = cellTopo.getEdgeCount();
161  /* */ int edgeOrts[4] = {};
162  for(index_t c=0;c<workset.num_cells;c++) {
163  orientations->at(details.cell_local_ids[c]).getEdgeOrientation(edgeOrts, numEdges);
164 
165  Intrepid2::Orientation::getReferenceEdgeTangent(ortEdgeTan,
166  subcellOrd,
167  cellTopo,
168  edgeOrts[subcellOrd],
169  is_normalize);
170 
171  for (int i=0;i<ndofsEdge;++i) {
172  const int b = intrepid_basis->getDofOrdinal(1, subcellOrd, i);
173  auto J = Kokkos::subview(worksetJacobians, c, b, Kokkos::ALL(), Kokkos::ALL());
174  Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan, J, ortEdgeTan);
175 
176  for(int d=0;d<cellDim;d++) {
177  residual(c,b) += (dof(c,b,d)-value(c,b,d))*phyEdgeTan(d);
178  }
179  }
180  }
181  }
182  break;
183  }
184  case 2: { // 3D element Tet and Hex
185  const int numEdges = cellTopo.getEdgeCount();
186  const int numFaces = cellTopo.getFaceCount();
187 
188  {
189  auto phyEdgeTan = Kokkos::subview(work, 0, Kokkos::ALL());
190  auto ortEdgeTan = Kokkos::subview(work, 1, Kokkos::ALL());
191 
192  const int numEdgesOfFace= cellTopo.getEdgeCount(2, subcellOrd);
193 
194  int edgeOrts[12] = {};
195  for(index_t c=0;c<workset.num_cells;c++) {
196  for (int i=0;i<numEdgesOfFace;++i) {
197 
198  const int edgeOrd = Intrepid2::Orientation::getEdgeOrdinalOfFace(i, subcellOrd, cellTopo);
199  const int b = edgeOrd;
200  orientations->at(details.cell_local_ids[c]).getEdgeOrientation(edgeOrts, numEdges);
201 
202  Intrepid2::Orientation::getReferenceEdgeTangent(ortEdgeTan,
203  edgeOrd,
204  cellTopo,
205  edgeOrts[edgeOrd],
206  is_normalize);
207 
208  // for(int b=0;b<dof.extent_int(1);b++)
209  {
210  auto J = Kokkos::subview(worksetJacobians, c, b, Kokkos::ALL(), Kokkos::ALL());
211  Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan, J, ortEdgeTan);
212 
213  for(int d=0;d<dof.extent_int(2);d++) {
214  residual(c,b) += (dof(c,b,d)-value(c,b,d))*phyEdgeTan(d);
215  }
216  }
217  }
218  }
219  }
220 
221  if (intrepid_basis->getDofCount(2, subcellOrd)) {
222  auto phyFaceTanU = Kokkos::subview(work, 0, Kokkos::ALL());
223  auto ortFaceTanU = Kokkos::subview(work, 1, Kokkos::ALL());
224  auto phyFaceTanV = Kokkos::subview(work, 2, Kokkos::ALL());
225  auto ortFaceTanV = Kokkos::subview(work, 3, Kokkos::ALL());
226 
227  int faceOrts[6] = {};
228  for(index_t c=0;c<workset.num_cells;c++) {
229  orientations->at(details.cell_local_ids[c]).getFaceOrientation(faceOrts, numFaces);
230  Intrepid2::Orientation::getReferenceFaceTangents(ortFaceTanU,
231  ortFaceTanV,
232  subcellOrd,
233  cellTopo,
234  faceOrts[subcellOrd],
235  is_normalize);
236 
237  for(int b=0;b<dof.extent_int(1);b++) {
238  auto J = Kokkos::subview(worksetJacobians, c, b, Kokkos::ALL(), Kokkos::ALL());
239  Intrepid2::Kernels::Serial::matvec_product(phyFaceTanU, J, ortFaceTanU);
240  Intrepid2::Kernels::Serial::matvec_product(phyFaceTanV, J, ortFaceTanV);
241 
242  for(int d=0;d<dof.extent_int(2);d++) {
243  residual(c,b) += (dof(c,b,d)-value(c,b,d))*phyFaceTanU(d);
244  residual(c,b) += (dof(c,b,d)-value(c,b,d))*phyFaceTanV(d);
245  }
246  }
247  }
248  }
249 
250  break;
251  }
252  }
253  }
254 
255 }
256 
257 //**********************************************************************
258 
259 }
260 
261 #endif
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.
Teuchos::RCP< const std::vector< Intrepid2::Orientation > > orientations_
T & get(const std::string &name, T def_value)
void setupArrays(const Teuchos::RCP< const panzer::PointRule > &pr)
Sizes/allocates memory for arrays.
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
#define TEUCHOS_ASSERT(assertion_test)
std::vector< GO > cell_local_ids