Panzer Version of the Day
Loading...
Searching...
No Matches
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
59namespace panzer {
60
61//**********************************************************************
62template<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->addNonConstDependentField(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//**********************************************************************
111template<typename EvalT, typename Traits>
112void
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//**********************************************************************
132template<typename EvalT, typename Traits>
133void
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 auto worksetJacobians_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),worksetJacobians);
169 switch (subcellDim) {
170 case 1: { // 2D element Tri and Quad
171 if (intrepid_basis->getDofCount(1, subcellOrd)) {
172 auto ortEdgeTan = Kokkos::subview(work, 0, Kokkos::ALL());
173 auto phyEdgeTan = Kokkos::subview(work, 1, Kokkos::ALL());
174
175 const int ndofsEdge = intrepid_basis->getDofCount(1, subcellOrd);
176 const int numEdges = cellTopo.getEdgeCount();
177 /* */ int edgeOrts[4] = {};
178
179 for(index_t c=0;c<workset.num_cells;c++) {
180 orientations->at(details.cell_local_ids[c]).getEdgeOrientation(edgeOrts, numEdges);
181
182 Intrepid2::Impl::OrientationTools::getRefSubcellTangents(work,
183 edgeParam,
184 cellTopo.getKey(edgeDim,subcellOrd),
185 subcellOrd,
186 edgeOrts[subcellOrd]);
187
188 for (int i=0;i<ndofsEdge;++i) {
189 const int b = intrepid_basis->getDofOrdinal(1, subcellOrd, i);
190 auto J = Kokkos::subview(worksetJacobians_h, c, b, Kokkos::ALL(), Kokkos::ALL());
191 Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan, J, ortEdgeTan);
192
193 for(int d=0;d<cellDim;d++) {
194 residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyEdgeTan(d);
195 }
196 }
197 }
198 }
199 break;
200 }
201 case 2: { // 3D element Tet and Hex
202 const int numEdges = cellTopo.getEdgeCount();
203 const int numFaces = cellTopo.getFaceCount();
204
205 {
206
207 auto ortEdgeTan = Kokkos::subview(work, 0, Kokkos::ALL());
208 auto phyEdgeTan = Kokkos::subview(work, 1, Kokkos::ALL());
209
210 const int numEdgesOfFace= cellTopo.getEdgeCount(2, subcellOrd);
211
212 int edgeOrts[12] = {};
213 for(index_t c=0;c<workset.num_cells;c++) {
214 for (int i=0;i<numEdgesOfFace;++i) {
215
216 const int edgeOrd = Intrepid2::Orientation::getEdgeOrdinalOfFace(i, subcellOrd, cellTopo);
217 const int b = edgeOrd;
218 orientations->at(details.cell_local_ids[c]).getEdgeOrientation(edgeOrts, numEdges);
219
220 Intrepid2::Impl::OrientationTools::getRefSubcellTangents(work,
221 edgeParam,
222 cellTopo.getKey(edgeDim,edgeOrd),
223 edgeOrd,
224 edgeOrts[edgeOrd]);
225
226 {
227 auto J = Kokkos::subview(worksetJacobians_h, c, b, Kokkos::ALL(), Kokkos::ALL());
228 Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan, J, ortEdgeTan);
229
230 for(int d=0;d<dof.extent_int(2);d++) {
231 residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyEdgeTan(d);
232 }
233 }
234 }
235 }
236 }
237
238 if (intrepid_basis->getDofCount(2, subcellOrd)) {
239 auto ortFaceTanU = Kokkos::subview(work, 0, Kokkos::ALL());
240 auto ortFaceTanV = Kokkos::subview(work, 1, Kokkos::ALL());
241 auto phyFaceTanU = Kokkos::subview(work, 2, Kokkos::ALL());
242 auto phyFaceTanV = Kokkos::subview(work, 3, Kokkos::ALL());
243
244 int faceOrts[6] = {};
245 for(index_t c=0;c<workset.num_cells;c++) {
246 orientations->at(details.cell_local_ids[c]).getFaceOrientation(faceOrts, numFaces);
247
248 Intrepid2::Impl::OrientationTools::getRefSubcellTangents(work,
249 faceParam,
250 cellTopo.getKey(faceDim,subcellOrd),
251 subcellOrd,
252 faceOrts[subcellOrd]);
253
254 for(int b=0;b<dof.extent_int(1);b++) {
255 auto J = Kokkos::subview(worksetJacobians, c, b, Kokkos::ALL(), Kokkos::ALL());
256 Intrepid2::Kernels::Serial::matvec_product(phyFaceTanU, J, ortFaceTanU);
257 Intrepid2::Kernels::Serial::matvec_product(phyFaceTanV, J, ortFaceTanV);
258
259 for(int d=0;d<dof.extent_int(2);d++) {
260 residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyFaceTanU(d);
261 residual_h(c,b) += (dof_h(c,b,d)-value_h(c,b,d))*phyFaceTanV(d);
262 }
263 }
264 }
265 }
266
267 break;
268 }
269 }
270 Kokkos::deep_copy(residual.get_static_view(), residual_h);
271 }
272
273}
274
275//**********************************************************************
276
277}
278
279#endif
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
int subcell_dim
DEPRECATED - use: getSubcellDimension()
int num_cells
DEPRECATED - use: numCells()
Teuchos::RCP< const std::vector< Intrepid2::Orientation > > orientations_