Panzer  Version of the Day
Panzer_IntegrationValues2.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 
44 #include "Panzer_UtilityAlgs.hpp"
45 
46 #include "Shards_CellTopology.hpp"
47 
48 #include "Kokkos_DynRankView.hpp"
49 #include "Intrepid2_FunctionSpaceTools.hpp"
50 #include "Intrepid2_RealSpaceTools.hpp"
51 #include "Intrepid2_CellTools.hpp"
52 #include "Intrepid2_ArrayTools.hpp"
53 #include "Intrepid2_CubatureControlVolume.hpp"
54 #include "Intrepid2_CubatureControlVolumeSide.hpp"
55 #include "Intrepid2_CubatureControlVolumeBoundary.hpp"
56 
58 #include "Panzer_Traits.hpp"
61 
62 // For device member functions to avoid CUDA RDC in unit testing
64 
65 // ***********************************************************
66 // * Specializations of setupArrays() for different array types
67 // ***********************************************************
68 
69 namespace panzer {
70 
71 // * Specialization for Kokkos::DynRankView<double,PHX::Device>
72 template <typename Scalar>
74 setupArraysForNodeRule(const Teuchos::RCP<const panzer::IntegrationRule>& ir)
75 {
76  MDFieldArrayFactory af(prefix,alloc_arrays);
77 
78  int num_nodes = ir->topology->getNodeCount();
79  int num_cells = ir->workset_size;
80  int num_space_dim = ir->topology->getDimension();
81 
82  int num_ip = 1;
83 
84  dyn_cub_points = af.template buildArray<double,IP,Dim>("cub_points",num_ip, num_space_dim);
85  dyn_cub_weights = af.template buildArray<double,IP>("cub_weights",num_ip);
86 
87  cub_points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
88 
89  if (ir->cv_type == "none" && ir->isSide()) {
90  dyn_side_cub_points = af.template buildArray<double,IP,Dim>("side_cub_points",num_ip, ir->side_topology->getDimension());
91  side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip,ir->side_topology->getDimension());
92  }
93 
94  if (ir->cv_type != "none") {
95  dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>("phys_cub_points",num_cells, num_ip, num_space_dim);
96  dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>("phys_cub_weights",num_cells, num_ip);
97  if (ir->cv_type == "side") {
98  dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>("phys_cub_norms",num_cells, num_ip, num_space_dim);
99  }
100  }
101 
102  dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
103 
104  cub_weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
105 
106  node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
107 
108  jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells, num_ip, num_space_dim,num_space_dim);
109 
110  jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
111 
112  jac_det = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells, num_ip);
113 
114  weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells, num_ip);
115 
116  covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells, num_ip, num_space_dim,num_space_dim);
117 
118  contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
119 
120  norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells, num_ip);
121 
122  ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordiantes",num_cells, num_ip,num_space_dim);
123 
124  ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells, num_ip,num_space_dim);
125 
126  weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted normal",num_cells, num_ip,num_space_dim);
127 
128  surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("surface_normals",num_cells, num_ip,num_space_dim);
129 
130  surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("surface_rotation_matrices",num_cells, num_ip,3,3);
131 
132 }
133 
134 template <typename Scalar>
136 setupArrays(const Teuchos::RCP<const panzer::IntegrationRule>& ir)
137 {
138  MDFieldArrayFactory af(prefix,alloc_arrays);
139 
141 
142  int_rule = ir;
143 
144  int num_nodes = ir->topology->getNodeCount();
145  int num_cells = ir->workset_size;
146  int num_space_dim = ir->topology->getDimension();
147 
148  // specialize content if this is quadrature at anode
149  if(num_space_dim==1 && ir->isSide()) {
150  setupArraysForNodeRule(ir);
151  return;
152  }
153 
154  TEUCHOS_ASSERT(ir->getType() != ID::NONE);
155  intrepid_cubature = getIntrepidCubature(*ir);
156 
157  int num_ip = ir->num_points;
158 
159 // Intrepid2::DefaultCubatureFactory cubature_factory;
160 //
161 // if (ir->cv_type == "side")
162 // intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir->topology));
163 //
164 // else if (ir->cv_type == "volume")
165 // intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir->topology));
166 //
167 // else if (ir->cv_type == "boundary" && ir->isSide())
168 // intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir->topology,ir->side));
169 //
170 // else if (ir->cv_type == "none" && ir->isSide())
171 // intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir->side_topology),
172 // ir->cubature_degree);
173 // else
174 // intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir->topology),
175 // ir->cubature_degree);
176 // int num_ip = intrepid_cubature->getNumPoints();
177 
178  dyn_cub_points = af.template buildArray<double,IP,Dim>("cub_points",num_ip, num_space_dim);
179  dyn_cub_weights = af.template buildArray<double,IP>("cub_weights",num_ip);
180 
181  cub_points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
182 
183  if (ir->isSide() && ir->cv_type == "none") {
184  dyn_side_cub_points = af.template buildArray<double,IP,Dim>("side_cub_points",num_ip, ir->side_topology->getDimension());
185  side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip,ir->side_topology->getDimension());
186  }
187 
188  if (ir->cv_type != "none") {
189  dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>("phys_cub_points",num_cells, num_ip, num_space_dim);
190  dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>("phys_cub_weights",num_cells, num_ip);
191  if (ir->cv_type == "side") {
192  dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>("phys_cub_norms",num_cells, num_ip, num_space_dim);
193  }
194  }
195 
196  dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
197 
198  cub_weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
199 
200  node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
201 
202  jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells, num_ip, num_space_dim,num_space_dim);
203 
204  jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
205 
206  jac_det = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells, num_ip);
207 
208  weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells, num_ip);
209 
210  covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells, num_ip, num_space_dim,num_space_dim);
211 
212  contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
213 
214  norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells, num_ip);
215 
216  ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordiantes",num_cells, num_ip,num_space_dim);
217 
218  ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells, num_ip,num_space_dim);
219 
220  weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted_normal",num_cells,num_ip,num_space_dim);
221 
222  surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("surface_normals",num_cells, num_ip,num_space_dim);
223 
224  surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("surface_rotation_matrices",num_cells, num_ip,3,3);
225 
226  scratch_for_compute_side_measure =
227  af.template buildStaticArray<Scalar,Point>("scratch_for_compute_side_measure", jac.get_view().span());
228 
229 }
230 
231 template <typename Scalar>
232 Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double>>
235 {
237  Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double> > ic;
238 
239  Intrepid2::DefaultCubatureFactory cubature_factory;
240 
241  if(ir.getType() == ID::CV_SIDE){
242  ic = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir.topology));
243  } else if(ir.getType() == ID::CV_VOLUME){
244  ic = Teuchos::rcp(new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir.topology));
245  } else if(ir.getType() == ID::CV_BOUNDARY){
246  ic = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir.topology,ir.getSide()));
247  } else {
248  if(ir.getType() == ID::VOLUME){
249  ic = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir.topology),ir.getOrder());
250  } else if(ir.getType() == ID::SIDE){
251  ic = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir.side_topology),ir.getOrder());
252  } else if(ir.getType() == ID::SURFACE){
253  // closed surface integrals don't exist in intrepid.
254  } else {
255  TEUCHOS_ASSERT(false);
256  }
257  }
258 
259  return ic;
260 }
261 
262 
263 // ***********************************************************
264 // * Evaluation of values - NOT specialized
265 // ***********************************************************
266 template <typename Scalar>
268 evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim> & in_node_coordinates,
269  const int in_num_cells,
270  const Teuchos::RCP<const SubcellConnectivity> & face_connectivity)
271 {
273  const bool is_surface = int_rule->getType() == ID::SURFACE;
274  const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
275 
276  TEUCHOS_ASSERT(not (is_surface and is_cv));
277 
278  if(is_surface){
279  TEUCHOS_TEST_FOR_EXCEPT_MSG(face_connectivity == Teuchos::null,
280  "IntegrationValues2::evaluateValues : Surface integration requires the face connectivity");
281  generateSurfaceCubatureValues(in_node_coordinates,in_num_cells,*face_connectivity);
282  } else if (is_cv) {
283  getCubatureCV(in_node_coordinates, in_num_cells);
284  evaluateValuesCV(in_node_coordinates, in_num_cells);
285  } else {
286  getCubature(in_node_coordinates, in_num_cells);
287  evaluateRemainingValues(in_node_coordinates, in_num_cells);
288  }
289 }
290 
291 template <typename Scalar>
293 getCubature(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
294  const int in_num_cells)
295 {
296 
297  int num_space_dim = int_rule->topology->getDimension();
298  if (int_rule->isSide() && num_space_dim==1) {
299  std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
300  << "non-natural integration rules.";
301  return;
302  }
303 
304  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
305 
306  if (!int_rule->isSide())
307  intrepid_cubature->getCubature(dyn_cub_points.get_view(), dyn_cub_weights.get_view());
308  else {
309  intrepid_cubature->getCubature(dyn_side_cub_points.get_view(), dyn_cub_weights.get_view());
310 
311  cell_tools.mapToReferenceSubcell(dyn_cub_points.get_view(),
312  dyn_side_cub_points.get_view(),
313  int_rule->spatial_dimension-1,
314  int_rule->side,
315  *(int_rule->topology));
316  }
317 
318  // IP coordinates
319  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
320  auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
321  auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
322  cell_tools.mapToPhysicalFrame(s_ip_coordinates,
323  dyn_cub_points.get_view(),
324  s_in_node_coordinates,
325  *(int_rule->topology));
326 }
327 
328 template <typename Scalar>
330 generateSurfaceCubatureValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
331  const int in_num_cells,
332  const SubcellConnectivity & face_connectivity)
333 {
334 
335  TEUCHOS_ASSERT(int_rule->getType() == IntegrationDescriptor::SURFACE);
336 
337  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
338 
339  const shards::CellTopology & cell_topology = *(int_rule->topology);
340  const panzer::IntegrationRule & ir = *int_rule;
341 
342  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
343 
344  // Copy over coordinates
345  {
346  const int num_nodes = in_node_coordinates.extent(1);
347  const int num_dims = in_node_coordinates.extent(2);
348  auto node_coordinates_k = node_coordinates.get_view();
349  auto in_node_coordinates_k = in_node_coordinates.get_view();
350 
351  Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<3>> policy({0,0,0},{num_cells,num_nodes,num_dims});
352  Kokkos::parallel_for("node_coordinates",policy,KOKKOS_LAMBDA (const int cell, const int node, const int dim) {
353  node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim);
354  });
355  PHX::Device::execution_space().fence();
356  }
357 
358  // NOTE: We are assuming that each face can have a different number of points.
359  // Not sure if this is necessary, but it requires a lot of additional allocations
360 
361  const int cell_dim = cell_topology.getDimension();
362  const int subcell_dim = cell_topology.getDimension()-1;
363  const int num_subcells = cell_topology.getSubcellCount(subcell_dim);
364 
365  Intrepid2::DefaultCubatureFactory cubature_factory;
366 
367  // We get to build up our cubature one face at a time
368  {
369  int point_offset=0;
370  for(int subcell_index=0; subcell_index<num_subcells; ++subcell_index) {
371 
372  // Default for 1D
373  int num_points_on_face = 1;
374 
375  // Get the cubature for the side
376  Kokkos::DynRankView<double,PHX::Device> tmp_side_cub_weights;
377  Kokkos::DynRankView<double,PHX::Device> tmp_side_cub_points;
378  if(cell_dim==1){
379  tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>("tmp_side_cub_weights",num_points_on_face);
380  tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>("cell_tmp_side_cub_points",num_points_on_face,cell_dim);
381  auto tmp_side_cub_weights_host = Kokkos::create_mirror_view(tmp_side_cub_weights);
382  auto tmp_side_cub_points_host = Kokkos::create_mirror_view(tmp_side_cub_points);
383  tmp_side_cub_weights_host(0)=1.;
384  tmp_side_cub_points_host(0,0) = (subcell_index==0)? -1. : 1.;
385  Kokkos::deep_copy(tmp_side_cub_weights,tmp_side_cub_weights_host);
386  Kokkos::deep_copy(tmp_side_cub_points,tmp_side_cub_points_host);
387  } else {
388 
389  // Get the face topology from the cell topology
390  const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,subcell_index));
391 
392  auto ic = cubature_factory.create<PHX::Device::execution_space,double,double>(face_topology,ir.getOrder());
393  num_points_on_face = ic->getNumPoints();
394 
395  tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>("tmp_side_cub_weights",num_points_on_face);
396  tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>("cell_tmp_side_cub_points",num_points_on_face,cell_dim);
397 
398  auto subcell_cub_points = Kokkos::DynRankView<double,PHX::Device>("subcell_cub_points",num_points_on_face,subcell_dim);
399 
400  // Get the reference face points
401  ic->getCubature(subcell_cub_points, tmp_side_cub_weights);
402 
403  // Convert from reference face points to reference cell points
404  cell_tools.mapToReferenceSubcell(tmp_side_cub_points,
405  subcell_cub_points,
406  subcell_dim,
407  subcell_index,
408  cell_topology);
409  }
410 
411  // Do this on host
412  {
413  auto tmp_side_cub_points_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),tmp_side_cub_points);
414  auto cub_points_host = Kokkos::create_mirror_view(cub_points.get_static_view());
415  for(int local_point=0;local_point<num_points_on_face;++local_point){
416  const int point = point_offset + local_point;
417  for(int dim=0;dim<cell_dim;++dim){
418  cub_points_host(point,dim) = tmp_side_cub_points_host(local_point,dim);
419  }
420  }
421  Kokkos::deep_copy(cub_points.get_static_view(),cub_points_host);
422  }
423 
424 
425  // Map from side points to physical points
426  auto side_ip_coordinates = Kokkos::DynRankView<Scalar,PHX::Device>("side_ip_coordinates",num_cells,num_points_on_face,cell_dim);
427  auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL);
428  cell_tools.mapToPhysicalFrame(side_ip_coordinates,
429  tmp_side_cub_points,
430  s_node_coordinates,
431  cell_topology);
432 
433  // Create a jacobian and his friends for this side
434  auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
435  cell_tools.setJacobian(side_jacobian,
436  tmp_side_cub_points,
437  s_node_coordinates,
438  cell_topology);
439 
440  PHX::Device::execution_space().fence();
441 
442  auto side_inverse_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_inv_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
443  cell_tools.setJacobianInv(side_inverse_jacobian, side_jacobian);
444 
445  auto side_det_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_det_jac",num_cells,num_points_on_face);
446  cell_tools.setJacobianDet(side_det_jacobian, side_jacobian);
447 
448  PHX::Device::execution_space().fence();
449 
450  // Calculate measures (quadrature weights in physical space) for this side
451  auto side_weighted_measure = Kokkos::DynRankView<Scalar,PHX::Device>("side_weighted_measure",num_cells,num_points_on_face);
452  if(cell_dim == 1){
453  Kokkos::deep_copy(side_weighted_measure, tmp_side_cub_weights(0));
454  } else if(cell_dim == 2){
455  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
456  computeEdgeMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights,
457  subcell_index,cell_topology,
458  scratch_for_compute_side_measure.get_view());
459  PHX::Device::execution_space().fence();
460  } else if(cell_dim == 3){
461  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
462  computeFaceMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights,
463  subcell_index,cell_topology,
464  scratch_for_compute_side_measure.get_view());
465  PHX::Device::execution_space().fence();
466  }
467 
468  // Calculate normals
469  auto side_normals = Kokkos::DynRankView<Scalar,PHX::Device>("side_normals",num_cells,num_points_on_face,cell_dim);
470  if(cell_dim == 1){
471 
472  const int other_subcell_index = (subcell_index==0) ? 1 : 0;
473  auto in_node_coordinates_k = in_node_coordinates.get_view();
474  const auto min = std::numeric_limits<typename Sacado::ScalarType<Scalar>::type>::min();
475 
476  Kokkos::parallel_for("compute normals 1D",num_cells,KOKKOS_LAMBDA (const int cell) {
477  Scalar norm = (in_node_coordinates_k(cell,subcell_index,0) - in_node_coordinates_k(cell,other_subcell_index,0));
478  side_normals(cell,0,0) = norm / fabs(norm + min);
479  });
480 
481  } else {
482 
483  cell_tools.getPhysicalSideNormals(side_normals,side_jacobian,subcell_index,cell_topology);
484 
485  // Normalize each normal
486  Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_cells,num_points_on_face});
487  Kokkos::parallel_for("Normalize the normals",policy,KOKKOS_LAMBDA (const int cell,const int point) {
488  Scalar n = 0.;
489  for(int dim=0;dim<cell_dim;++dim){
490  n += side_normals(cell,point,dim)*side_normals(cell,point,dim);
491  }
492  // If n is zero then this is - hopefully - a virtual cell
493  if(n > 0.){
494  n = Kokkos::Experimental::sqrt(n);
495  for(int dim=0;dim<cell_dim;++dim){
496  side_normals(cell,point,dim) /= n;
497  }
498  }
499  });
500 
501  }
502  PHX::Device::execution_space().fence();
503 
504  // Now that we have all these wonderful values, lets copy them to the actual arrays
505  {
506  auto weighted_measure_k = weighted_measure.get_view();
507  auto jac_k = jac.get_view();
508  auto jac_inv_k = jac_inv.get_view();
509  auto jac_det_k = jac_det.get_view();
510  auto ref_ip_coordinates_k = ref_ip_coordinates.get_view();
511  auto ip_coordinates_k = ip_coordinates.get_view();
512  auto surface_normals_k = surface_normals.get_view();
513  auto cub_points_k = cub_points.get_view();
514  Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_cells,num_points_on_face});
515  Kokkos::parallel_for("copy values",policy,KOKKOS_LAMBDA (const int cell,const int side_point) {
516  const int cell_point = point_offset + side_point;
517 
518  weighted_measure_k(cell,cell_point) = side_weighted_measure(cell,side_point);
519  jac_det_k(cell,cell_point) = side_det_jacobian(cell,side_point);
520  for(int dim=0;dim<cell_dim;++dim){
521  ref_ip_coordinates_k(cell,cell_point,dim) = cub_points_k(cell_point,dim);
522  ip_coordinates_k(cell,cell_point,dim) = side_ip_coordinates(cell,side_point,dim);
523  surface_normals_k(cell,cell_point,dim) = side_normals(cell,side_point,dim);
524 
525  for(int dim2=0;dim2<cell_dim;++dim2){
526  jac_k(cell,cell_point,dim,dim2) = side_jacobian(cell,side_point,dim,dim2);
527  jac_inv_k(cell,cell_point,dim,dim2) = side_inverse_jacobian(cell,side_point,dim,dim2);
528  }
529  }
530  });
531  PHX::Device::execution_space().fence();
532  }
533  point_offset += num_points_on_face;
534  }
535  }
536 
537  // We also need surface rotation matrices in order to enforce alignment
538  {
539  const int num_points = ir.getPointOffset(num_subcells);
540  auto surface_normals_k = surface_normals.get_view();
541  auto surface_rotation_matrices_k = surface_rotation_matrices.get_view();
542  Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<3>> policy({0,0,0},{num_cells,num_subcells,num_points});
543  Kokkos::parallel_for("create surface rotation matrices",policy,KOKKOS_LAMBDA (const int cell,const int subcell_index,const int point) {
544  Scalar normal[3];
545 
546  for(int i=0;i<3;i++){normal[i]=0.;}
547 
548  for(int dim=0; dim<cell_dim; ++dim){
549  normal[dim] = surface_normals_k(cell,point,dim);
550  }
551 
552  Scalar transverse[3];
553  Scalar binormal[3];
554  panzer::convertNormalToRotationMatrix(normal,transverse,binormal);
555 
556  for(int dim=0; dim<3; ++dim){
557  surface_rotation_matrices_k(cell,point,0,dim) = normal[dim];
558  surface_rotation_matrices_k(cell,point,1,dim) = transverse[dim];
559  surface_rotation_matrices_k(cell,point,2,dim) = binormal[dim];
560  }
561  });
562  PHX::Device::execution_space().fence();
563  }
564 
565  // =========================================================
566  // Enforce alignment across surface quadrature points
567 
568  const int num_points = ip_coordinates.extent_int(1);
569  const int num_faces_per_cell = face_connectivity.numSubcellsOnCellHost(0);
570  const int num_points_per_face = num_points / num_faces_per_cell;
571 
572  // Now we need to align the cubature points for each face
573  // If there is only one point there is no need to align things
574  if(num_points_per_face != 1){
575  // To enforce that quadrature points on faces are aligned properly we will iterate through faces,
576  // map points to a plane shared by the faces, then re-order quadrature points on the "1" face to
577  // line up with the "0" face
578 
579  // Utility calls
580 #define PANZER_DOT(a,b) (a[0]*b[0] + a[1]*b[1] + a[2]*b[2])
581 #define PANZER_CROSS(a,b,c) {c[0] = a[1]*b[2] - a[2]*b[1]; c[1] = a[2]*b[0] - a[0]*b[2]; c[2] = a[0]*b[1] - a[1]*b[0];}
582 
583  // Reorder index vector. This is scratch space that needs to be
584  // allocated per thread now. Should use team size instead of num
585  // faces, but then we need to convert the lambda below to a
586  // functor to query the policy using the functor. We'll just live
587  // with the extra memory since it is temporary scratch.
588  // std::vector<int> point_order(num_points_per_face);
589  PHX::View<int**> point_order("scratch: point_order",face_connectivity.numSubcells(),num_points_per_face);
590 
591  // Iterate through faces
592  auto ref_ip_coordinates_k = ref_ip_coordinates.get_view();
593  auto ip_coordinates_k = ip_coordinates.get_view();
594  auto weighted_measure_k = weighted_measure.get_view();
595  auto jac_k = jac.get_view();
596  auto jac_det_k = jac_det.get_view();
597  auto jac_inv_k = jac_inv.get_view();
598  auto surface_normals_k = surface_normals.get_view();
599  auto surface_rotation_matrices_k = surface_rotation_matrices.get_view();
600  Kokkos::parallel_for("face iteration",face_connectivity.numSubcells(),KOKKOS_LAMBDA (const int face) {
601  // Cells for sides 0 and 1
602  const int cell_0 = face_connectivity.cellForSubcell(face,0);
603  const int cell_1 = face_connectivity.cellForSubcell(face,1);
604 
605  // If this face doesn't connect to anything we don't need to worry about alignment
606  if(cell_1 < 0)
607  return;
608 
609  // Local face index for sides 0 and 1
610  const int lidx_0 = face_connectivity.localSubcellForSubcell(face,0);
611  const int lidx_1 = face_connectivity.localSubcellForSubcell(face,1);
612 
613  // If the cell exists, then the local face index needs to exist
614  KOKKOS_ASSERT(lidx_1 >= 0);
615 
616  // To compare points on planes, it is good to center the planes around the origin
617  // Find the face center for the left and right cell (may not be the same point - e.g. periodic conditions)
618  // We also define a length scale r2 which gives an order of magnitude approximation to the cell size squared
619  Scalar xc0[3] = {0}, xc1[3] = {0}, r2 = 0;
620  for(int face_point=0; face_point<num_points_per_face; ++face_point){
621  Scalar dx2 = 0.;
622  for(int dim=0; dim<cell_dim; ++dim){
623  xc0[dim] += ip_coordinates_k(cell_0,lidx_0*num_points_per_face + face_point,dim);
624  xc1[dim] += ip_coordinates_k(cell_1,lidx_1*num_points_per_face + face_point,dim);
625  const Scalar dx = ip_coordinates_k(cell_0,lidx_0*num_points_per_face + face_point,dim) - ip_coordinates_k(cell_0,lidx_0*num_points_per_face,dim);
626  dx2 += dx*dx;
627  }
628 
629  // Check if the distance squared between these two points is larger than the others (doesn't need to be perfect)
630  r2 = (r2 < dx2) ? dx2 : r2;
631  }
632  for(int dim=0; dim<cell_dim; ++dim){
633  xc0[dim] /= double(num_points_per_face);
634  xc1[dim] /= double(num_points_per_face);
635  }
636 
637  // TODO: This needs to be adaptable to having curved faces - for now this will work
638 
639  // We need to define a plane with two vectors (transverse and binormal)
640  // These will be used with the face centers to construct a local reference frame for aligning points
641 
642  // We use the first point on the face to define the normal (may break for curved faces)
643  const int example_point_0 = lidx_0*num_points_per_face;
644  const int example_point_1 = lidx_1*num_points_per_face;
645 
646  // Load the transverse and binormal for the 0 cell (default)
647  Scalar t[3] = {surface_rotation_matrices_k(cell_0,example_point_0,1,0), surface_rotation_matrices_k(cell_0,example_point_0,1,1), surface_rotation_matrices_k(cell_0,example_point_0,1,2)};
648  Scalar b[3] = {surface_rotation_matrices_k(cell_0,example_point_0,2,0), surface_rotation_matrices_k(cell_0,example_point_0,2,1), surface_rotation_matrices_k(cell_0,example_point_0,2,2)};
649 
650  // In case the faces are not antiparallel (e.g. periodic wedge), we may need to change the transverse and binormal
651  {
652 
653  // Get the normals for each face for constructing one of the plane vectors (transverse)
654  const Scalar n0[3] = {surface_rotation_matrices_k(cell_0,example_point_0,0,0), surface_rotation_matrices_k(cell_0,example_point_0,0,1), surface_rotation_matrices_k(cell_0,example_point_0,0,2)};
655  const Scalar n1[3] = {surface_rotation_matrices_k(cell_1,example_point_1,0,0), surface_rotation_matrices_k(cell_1,example_point_1,0,1), surface_rotation_matrices_k(cell_1,example_point_1,0,2)};
656 
657  // n_0*n_1 == -1 (antiparallel), n_0*n_1 == 1 (parallel - bad), |n_0*n_1| < 1 (other)
658  const Scalar n0_dot_n1 = PANZER_DOT(n0,n1);
659 
660  // FIXME: Currently virtual cells will set their surface normal along the same direction as the cell they "reflect"
661  // This causes a host of issues (e.g. identifying 180 degree periodic wedges), but we have to support virtual cells as a priority
662  // Therefore, we will just assume that the ordering is fine (not valid for 180 degree periodic wedges)
663  if(Kokkos::Experimental::fabs(n0_dot_n1 - 1.) < 1.e-8)
664  return;
665 
666  // Rotated faces case (e.g. periodic wedge) we need to check for non-antiparallel face normals
667  if(Kokkos::Experimental::fabs(n0_dot_n1 + 1.) > 1.e-2){
668 
669  // Now we need to define an arbitrary transverse and binormal in the plane across which the faces are anti-symmetric
670  // We can do this by setting t = n_0 \times n_1
671  PANZER_CROSS(n0,n1,t);
672 
673  // Normalize the transverse vector
674  const Scalar mag_t = Kokkos::Experimental::sqrt(PANZER_DOT(t,t));
675  t[0] /= mag_t;
676  t[1] /= mag_t;
677  t[2] /= mag_t;
678 
679  // The binormal will be the sum of the normals (does not need to be a right handed system)
680  b[0] = n0[0] + n1[0];
681  b[1] = n0[1] + n1[1];
682  b[2] = n0[2] + n1[2];
683 
684  // Normalize the binormal vector
685  const Scalar mag_b = Kokkos::Experimental::sqrt(PANZER_DOT(b,b));
686  b[0] /= mag_b;
687  b[1] /= mag_b;
688  b[2] /= mag_b;
689 
690  }
691  }
692 
693  // Now that we have a reference coordinate plane in which to align our points
694  // Point on the transverse/binormal plane
695  Scalar p0[2] = {0};
696  Scalar p1[2] = {0};
697 
698  // Differential position in mesh
699  Scalar x0[3] = {0};
700  Scalar x1[3] = {0};
701 
702  // Iterate through points in cell 1 and find which point they align with in cell 0
703  for(int face_point_1=0; face_point_1<num_points_per_face; ++face_point_1){
704 
705  // Get the point index in the 1 cell
706  const int point_1 = lidx_1*num_points_per_face + face_point_1;
707 
708  // Load position shifted by face center
709  for(int dim=0; dim<cell_dim; ++dim)
710  x1[dim] = ip_coordinates_k(cell_1,point_1,dim) - xc1[dim];
711 
712  // Convert position to transverse/binormal plane
713  p1[0] = PANZER_DOT(x1,t);
714  p1[1] = PANZER_DOT(x1,b);
715 
716  // Set the default for the point order
717  point_order(face,face_point_1) = face_point_1;
718 
719  // Compare to points on the other surface
720  for(int face_point_0=0; face_point_0<num_points_per_face; ++face_point_0){
721 
722  // Get the point index in the 0 cell
723  const int point_0 = lidx_0*num_points_per_face + face_point_0;
724 
725  // Load position shifted by face center
726  for(int dim=0; dim<cell_dim; ++dim)
727  x0[dim] = ip_coordinates_k(cell_0,point_0,dim) - xc0[dim];
728 
729  // Convert position to transverse/binormal plane
730  p0[0] = PANZER_DOT(x0,t);
731  p0[1] = PANZER_DOT(x0,b);
732 
733  // Find the distance squared between p0 and p1
734  const Scalar p012 = (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]);
735 
736  // If the distance, compared to the size of the cell, is small, we assume these are the same points
737  if(p012 / r2 < 1.e-12){
738  point_order(face,face_point_1) = face_point_0;
739  break;
740  }
741 
742  // Big problem - there wan't a point that aligned properly
743  KOKKOS_ASSERT(face_point_0 != num_points_per_face-1);
744 
745  }
746  }
747 
748  // Now re-order the points on face 1 to correct the alignment
749  const int point_offset = lidx_1*num_points_per_face;
750  for( int face_point_1 = 0; face_point_1 < num_points_per_face - 1; ++face_point_1 ){
751  // While the point is not yet in place - keep swapping until it is in place (N^2 operations - not great)
752  while( face_point_1 != point_order(face,face_point_1) ){
753  // We need to swap with the component in this position
754  const int face_point_0 = point_order(face,face_point_1);
755  panzer::swapQuadraturePoints<double>(cell_1,point_offset+face_point_1,point_offset+face_point_0,
756  ref_ip_coordinates_k,
757  ip_coordinates_k,
758  weighted_measure_k,
759  jac_k,
760  jac_det_k,
761  jac_inv_k,
762  surface_normals_k,
763  surface_rotation_matrices_k);
764  Scalar tmp = point_order(face,face_point_1);
765  point_order(face,face_point_1) = point_order(face,face_point_0);
766  point_order(face,face_point_0) = tmp;
767  }
768  }
769 
770  });
771  PHX::Device::execution_space().fence();
772  }
773 
774 #undef PANZER_DOT
775 #undef PANZER_CROSS
776 
777  // =========================================================
778 
779  // I'm not sure if these should exist for surface integrals, but here we go!
780 
781  // Shakib contravarient metric tensor
782  {
783  auto contravarient_k = this->contravarient.get_static_view();
784  auto covarient_k = this->covarient.get_static_view();
785  auto jac_k = this->jac.get_static_view();
786  Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{(int)num_cells,(int)contravarient.extent(1)});
787  Kokkos::parallel_for("covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) {
788  // zero out matrix
789  for (size_type i = 0; i < contravarient_k.extent(2); ++i)
790  for (size_type j = 0; j < contravarient_k.extent(3); ++j)
791  covarient_k(cell,ip,i,j) = 0.0;
792 
793  // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha}
794  for (size_type i = 0; i < contravarient_k.extent(2); ++i) {
795  for (size_type j = 0; j < contravarient_k.extent(3); ++j) {
796  for (size_type alpha = 0; alpha < contravarient_k.extent(2); ++alpha) {
797  covarient_k(cell,ip,i,j) += jac_k(cell,ip,i,alpha) * jac_k(cell,ip,j,alpha);
798  }
799  }
800  }
801  });
802  PHX::Device::execution_space().fence();
803  }
804 
805  {
806  auto s_contravarient = Kokkos::subview(contravarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
807  auto s_covarient = Kokkos::subview(covarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
808  Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
809  PHX::Device::execution_space().fence();
810  }
811 
812  // norm of g_ij
813  {
814  auto contravarient_k = this->contravarient.get_static_view();
815  auto norm_contravarient_k = this->norm_contravarient.get_static_view();
816  Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{(int)num_cells,(int)contravarient.extent(1)});
817  Kokkos::parallel_for("covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) {
818  norm_contravarient_k(cell,ip) = 0.0;
819  for (size_type i = 0; i < contravarient_k.extent(2); ++i) {
820  for (size_type j = 0; j < contravarient_k.extent(3); ++j) {
821  norm_contravarient_k(cell,ip) += contravarient_k(cell,ip,i,j) * contravarient_k(cell,ip,i,j);
822  }
823  }
824  norm_contravarient_k(cell,ip) = Kokkos::Experimental::sqrt(norm_contravarient_k(cell,ip));
825  });
826  PHX::Device::execution_space().fence();
827  }
828 
829 }
830 
831 
832 template <typename Scalar>
834 evaluateRemainingValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
835  const int in_num_cells)
836 {
837  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
838 
839  // copy the dynamic data structures into the static data structures
840  {
841  Kokkos::deep_copy(cub_weights.get_static_view(),dyn_cub_weights.get_view());
842  Kokkos::deep_copy(cub_points.get_static_view(),dyn_cub_points.get_view());
843  }
844 
845  if (int_rule->isSide()) {
846  Kokkos::deep_copy(side_cub_points.get_static_view(),dyn_side_cub_points.get_view());
847  }
848 
849  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
850 
851  {
852  size_type num_nodes = in_node_coordinates.extent(1);
853  size_type num_dims = in_node_coordinates.extent(2);
854  auto node_coordinates_k = node_coordinates.get_view();
855  auto in_node_coordinates_k = in_node_coordinates.get_view();
856 
857  Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{(int)num_cells,(int)num_nodes,(int)num_dims});
858  Kokkos::parallel_for("node coordinates",policy,KOKKOS_LAMBDA (const int cell,const int node,const int dim) {
859  node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim);
860  });
861  PHX::Device::execution_space().fence();
862  }
863 
864  auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
865  auto s_jac = Kokkos::subview(jac.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
866  cell_tools.setJacobian(jac.get_view(),
867  cub_points.get_view(),
868  node_coordinates.get_view(),
869  *(int_rule->topology));
870  PHX::Device::execution_space().fence();
871 
872  auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
873  cell_tools.setJacobianInv(s_jac_inv, s_jac);
874 
875  auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair(0,num_cells),Kokkos::ALL());
876  cell_tools.setJacobianDet(s_jac_det, s_jac);
877 
878  PHX::Device::execution_space().fence();
879 
880  auto s_weighted_measure = Kokkos::subview(weighted_measure.get_view(),std::make_pair(0,num_cells),Kokkos::ALL());
881  if (!int_rule->isSide()) {
882  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
883  computeCellMeasure(s_weighted_measure, s_jac_det, cub_weights.get_view());
884  }
885  else if(int_rule->spatial_dimension==3) {
886  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
887  computeFaceMeasure(s_weighted_measure, s_jac, cub_weights.get_view(),
888  int_rule->side, *int_rule->topology,
889  scratch_for_compute_side_measure.get_view());
890  }
891  else if(int_rule->spatial_dimension==2) {
892  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
893  computeEdgeMeasure(s_weighted_measure, s_jac, cub_weights.get_view(),
894  int_rule->side,*int_rule->topology,
895  scratch_for_compute_side_measure.get_view());
896  }
897  else TEUCHOS_ASSERT(false);
898 
899  PHX::Device::execution_space().fence();
900 
901  // Shakib contravarient metric tensor
902  {
903  auto covarient_k = covarient.get_view();
904  auto contravarient_k = contravarient.get_view();
905  auto jac_k = jac.get_view();
906  Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_cells,static_cast<int>(contravarient.extent(1))});
907  Kokkos::parallel_for("evalaute covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) {
908  // zero out matrix
909  for (int i = 0; i < static_cast<int>(contravarient_k.extent(2)); ++i)
910  for (int j = 0; j < static_cast<int>(contravarient_k.extent(3)); ++j)
911  covarient_k(cell,ip,i,j) = 0.0;
912 
913  // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha}
914  for (int i = 0; i < static_cast<int>(contravarient_k.extent(2)); ++i) {
915  for (int j = 0; j < static_cast<int>(contravarient_k.extent(3)); ++j) {
916  for (int alpha = 0; alpha < static_cast<int>(contravarient_k.extent(2)); ++alpha) {
917  covarient_k(cell,ip,i,j) += jac_k(cell,ip,i,alpha) * jac_k(cell,ip,j,alpha);
918  }
919  }
920  }
921  });
922  PHX::Device::execution_space().fence();
923  }
924 
925  auto s_covarient = Kokkos::subview(covarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
926  auto s_contravarient = Kokkos::subview(contravarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
927  Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
928  PHX::Device::execution_space().fence();
929 
930  // norm of g_ij
931  {
932  auto contravarient_k = contravarient.get_view();
933  auto norm_contravarient_k = norm_contravarient.get_view();
934  Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_cells,static_cast<int>(contravarient.extent(1))});
935  Kokkos::parallel_for("evaluate norm_contravarient",policy,KOKKOS_LAMBDA (const int cell,const int ip) {
936  norm_contravarient_k(cell,ip) = 0.0;
937  for (int i = 0; i < static_cast<int>(contravarient_k.extent(2)); ++i) {
938  for (int j = 0; j < static_cast<int>(contravarient_k.extent(3)); ++j) {
939  norm_contravarient_k(cell,ip) += contravarient_k(cell,ip,i,j) * contravarient_k(cell,ip,i,j);
940  }
941  }
942  norm_contravarient_k(cell,ip) = Kokkos::Experimental::sqrt(norm_contravarient_k(cell,ip));
943  });
944  PHX::Device::execution_space().fence();
945  }
946 }
947 
948 // Find the permutation that maps the set of points coords to other_coords. To
949 // avoid possible finite precision issues, == is not used, but rather
950 // min(norm(.)).
951 template <typename Scalar>
952 static void
953 permuteToOther(const PHX::MDField<Scalar,Cell,IP,Dim>& coords,
954  const PHX::MDField<Scalar,Cell,IP,Dim>& other_coords,
955  std::vector<typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type>& permutation)
956 {
957  typedef typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type size_type;
958  // We can safely assume: (1) The permutation is the same for every cell in
959  // the workset. (2) The first workset has valid data. Hence we operate only
960  // on cell 0.
961  const size_type cell = 0;
962  const size_type num_ip = coords.extent(1), num_dim = coords.extent(2);
963  permutation.resize(num_ip);
964  std::vector<char> taken(num_ip, 0);
965 
966  auto coords_view = coords.get_view();
967  auto coords_h = Kokkos::create_mirror_view(coords_view);
968  Kokkos::deep_copy(coords_h, coords_view);
969 
970  auto other_coords_view = other_coords.get_view();
971  auto other_coords_h = Kokkos::create_mirror_view(other_coords_view);
972  Kokkos::deep_copy(other_coords_h, other_coords_view);
973 
974  for (size_type ip = 0; ip < num_ip; ++ip) {
975  // Find an other point to associate with ip.
976  size_type i_min = 0;
977  Scalar d_min = -1;
978  for (size_type other_ip = 0; other_ip < num_ip; ++other_ip) {
979  // For speed, skip other points that are already associated.
980  if (taken[other_ip]) continue;
981  // Compute the distance between the two points.
982  Scalar d(0);
983  for (size_type dim = 0; dim < num_dim; ++dim) {
984  const Scalar diff = coords_h(cell, ip, dim) - other_coords_h(cell, other_ip, dim);
985  d += diff*diff;
986  }
987  if (d_min < 0 || d < d_min) {
988  d_min = d;
989  i_min = other_ip;
990  }
991  }
992  // Record the match.
993  permutation[ip] = i_min;
994  // This point is now taken.
995  taken[i_min] = 1;
996  }
997 }
998 
999 template <typename Scalar>
1001 evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
1002  const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates,
1003  const int in_num_cells)
1004 {
1005  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
1006 
1007  if (int_rule->cv_type == "none") {
1008 
1009  getCubature(in_node_coordinates, in_num_cells);
1010 
1011  {
1012  // Determine the permutation.
1013  std::vector<size_type> permutation(other_ip_coordinates.extent(1));
1014  permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
1015  // Apply the permutation to the cubature arrays.
1016  MDFieldArrayFactory af(prefix, alloc_arrays);
1017  {
1018  const size_type num_ip = dyn_cub_points.extent(0);
1019  {
1020  const size_type num_dim = dyn_side_cub_points.extent(1);
1021  DblArrayDynamic old_dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
1022  "old_dyn_side_cub_points", num_ip, num_dim);
1023  old_dyn_side_cub_points.deep_copy(dyn_side_cub_points);
1024 
1025  auto dyn_side_cub_points_h = Kokkos::create_mirror_view(as_view(dyn_side_cub_points));
1026  auto old_dyn_side_cub_points_h = Kokkos::create_mirror_view(as_view(old_dyn_side_cub_points));
1027  Kokkos::deep_copy(dyn_side_cub_points_h, as_view(dyn_side_cub_points));
1028  Kokkos::deep_copy(old_dyn_side_cub_points_h, as_view(old_dyn_side_cub_points));
1029 
1030  for (size_type ip = 0; ip < num_ip; ++ip)
1031  if (ip != permutation[ip])
1032  for (size_type dim = 0; dim < num_dim; ++dim)
1033  dyn_side_cub_points_h(ip, dim) = old_dyn_side_cub_points_h(permutation[ip], dim);
1034 
1035  Kokkos::deep_copy(as_view(dyn_side_cub_points), dyn_side_cub_points_h);
1036  }
1037  {
1038  const size_type num_dim = dyn_cub_points.extent(1);
1039  DblArrayDynamic old_dyn_cub_points = af.template buildArray<double,IP,Dim>(
1040  "old_dyn_cub_points", num_ip, num_dim);
1041  old_dyn_cub_points.deep_copy(dyn_cub_points);
1042 
1043  auto dyn_cub_points_h = Kokkos::create_mirror_view(as_view(dyn_cub_points));
1044  auto old_dyn_cub_points_h = Kokkos::create_mirror_view(as_view(old_dyn_cub_points));
1045  Kokkos::deep_copy(dyn_cub_points_h, as_view(dyn_cub_points));
1046  Kokkos::deep_copy(old_dyn_cub_points_h, as_view(old_dyn_cub_points));
1047 
1048  for (size_type ip = 0; ip < num_ip; ++ip)
1049  if (ip != permutation[ip])
1050  for (size_type dim = 0; dim < num_dim; ++dim)
1051  dyn_cub_points_h(ip, dim) = old_dyn_cub_points_h(permutation[ip], dim);
1052 
1053  Kokkos::deep_copy(as_view(dyn_cub_points), dyn_cub_points_h);
1054  }
1055  {
1056  DblArrayDynamic old_dyn_cub_weights = af.template buildArray<double,IP>(
1057  "old_dyn_cub_weights", num_ip);
1058  old_dyn_cub_weights.deep_copy(dyn_cub_weights);
1059 
1060  auto dyn_cub_weights_h = Kokkos::create_mirror_view(as_view(dyn_cub_weights));
1061  auto old_dyn_cub_weights_h = Kokkos::create_mirror_view(as_view(old_dyn_cub_weights));
1062  Kokkos::deep_copy(dyn_cub_weights_h, as_view(dyn_cub_weights));
1063  Kokkos::deep_copy(old_dyn_cub_weights_h, as_view(old_dyn_cub_weights));
1064 
1065  for (size_type ip = 0; ip < dyn_cub_weights.extent(0); ++ip)
1066  if (ip != permutation[ip])
1067  dyn_cub_weights_h(ip) = old_dyn_cub_weights_h(permutation[ip]);
1068 
1069  Kokkos::deep_copy(as_view(dyn_cub_weights), old_dyn_cub_weights_h);
1070  }
1071  }
1072  {
1073  const size_type num_ip = ip_coordinates.extent(1);
1074  const size_type num_dim = ip_coordinates.extent(2);
1075  Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1076  "old_ip_coordinates", ip_coordinates.extent(0), num_ip, num_dim);
1077 
1078  auto ip_coordinates_h = Kokkos::create_mirror_view(as_view(ip_coordinates));
1079  auto old_ip_coordinates_h = Kokkos::create_mirror_view(as_view(old_ip_coordinates));
1080  Kokkos::deep_copy(old_ip_coordinates_h, as_view(ip_coordinates));
1081  Kokkos::deep_copy(ip_coordinates_h, as_view(ip_coordinates));
1082 
1083  for (int cell = 0; cell < num_cells; ++cell)
1084  for (size_type ip = 0; ip < num_ip; ++ip)
1085  if (ip != permutation[ip])
1086  for (size_type dim = 0; dim < num_dim; ++dim)
1087  ip_coordinates_h(cell, ip, dim) = old_ip_coordinates_h(cell, permutation[ip], dim);
1088 
1089  Kokkos::deep_copy(as_view(ip_coordinates), ip_coordinates_h);
1090  }
1091  // All subsequent calculations inherit the permutation.
1092  }
1093 
1094  evaluateRemainingValues(in_node_coordinates, in_num_cells);
1095  }
1096 
1097  else {
1098 
1099  getCubatureCV(in_node_coordinates, in_num_cells);
1100 
1101  // Determine the permutation.
1102  std::vector<size_type> permutation(other_ip_coordinates.extent(1));
1103  permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
1104 
1105  // Apply the permutation to the cubature arrays.
1106  MDFieldArrayFactory af(prefix, alloc_arrays);
1107  {
1108  const size_type workset_size = ip_coordinates.extent(0), num_ip = ip_coordinates.extent(1),
1109  num_dim = ip_coordinates.extent(2);
1110  Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1111  "old_ip_coordinates", workset_size, num_ip, num_dim);
1112  Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
1113  Array_CellIPDim old_weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1114  "old_weighted_normals", workset_size, num_ip, num_dim);
1115  Array_CellIP old_weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
1116  "old_weighted_measure", workset_size, num_ip);
1117  if (int_rule->cv_type == "side")
1118  Kokkos::deep_copy(old_weighted_normals.get_static_view(), weighted_normals.get_static_view());
1119  else
1120  Kokkos::deep_copy(old_weighted_measure.get_static_view(), weighted_measure.get_static_view());
1121  for (int cell = 0; cell < num_cells; ++cell)
1122  {
1123  for (size_type ip = 0; ip < num_ip; ++ip)
1124  {
1125  if (ip != permutation[ip]) {
1126  if (int_rule->cv_type == "boundary" || int_rule->cv_type == "volume")
1127  weighted_measure(cell, ip) = old_weighted_measure(cell, permutation[ip]);
1128  for (size_type dim = 0; dim < num_dim; ++dim)
1129  {
1130  ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
1131  if (int_rule->cv_type == "side")
1132  weighted_normals(cell, ip, dim) = old_weighted_normals(cell, permutation[ip], dim);
1133 
1134  }
1135  }
1136  }
1137  }
1138  }
1139 
1140  evaluateValuesCV(in_node_coordinates, in_num_cells);
1141  }
1142 }
1143 
1144 template <typename Scalar>
1146 getCubatureCV(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
1147  const int in_num_cells)
1148 {
1149  int num_space_dim = int_rule->topology->getDimension();
1150  if (int_rule->isSide() && num_space_dim==1) {
1151  std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
1152  << "non-natural integration rules.";
1153  return;
1154  }
1155 
1156  size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (size_type) in_num_cells;
1157  std::pair<int,int> cell_range(0,num_cells);
1158  {
1159  size_type num_nodes = in_node_coordinates.extent(1);
1160  size_type num_dims = in_node_coordinates.extent(2);
1161  auto node_coordinates_k = node_coordinates.get_view();
1162  auto dyn_node_coordinates_k = dyn_node_coordinates.get_view();
1163  auto in_node_coordinates_k = in_node_coordinates.get_view();
1164 
1165  Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_cells,num_nodes,num_dims});
1166  Kokkos::parallel_for("getCubatureCV: node coordinates",policy,KOKKOS_LAMBDA (const int cell,const int node,const int dim) {
1167  node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim);
1168  dyn_node_coordinates_k(cell,node,dim) = Sacado::ScalarValue<Scalar>::eval(in_node_coordinates_k(cell,node,dim));
1169  });
1170  PHX::Device::execution_space().fence();
1171  }
1172 
1173  auto s_dyn_phys_cub_points = Kokkos::subdynrankview(dyn_phys_cub_points.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1174  auto s_dyn_node_coordinates = Kokkos::subdynrankview(dyn_node_coordinates.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1175  if (int_rule->cv_type == "side") {
1176  auto s_dyn_phys_cub_norms = Kokkos::subdynrankview(dyn_phys_cub_norms.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1177  intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_norms,s_dyn_node_coordinates);
1178  }
1179  else {
1180  auto s_dyn_phys_cub_weights = Kokkos::subdynrankview(dyn_phys_cub_weights.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1181  intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_weights,s_dyn_node_coordinates);
1182  }
1183 
1184  // size_type num_cells = dyn_phys_cub_points.extent(0);
1185  size_type num_ip =dyn_phys_cub_points.extent(1);
1186  size_type num_dims = dyn_phys_cub_points.extent(2);
1187  auto weighted_measure_k = weighted_measure.get_view();
1188  auto dyn_phys_cub_weights_k = dyn_phys_cub_weights.get_view();
1189  auto ip_coordinates_k = ip_coordinates.get_view();
1190  auto dyn_phys_cub_points_k = dyn_phys_cub_points.get_view();
1191  auto weighted_normals_k = weighted_normals.get_view();
1192  auto dyn_phys_cub_norms_k = dyn_phys_cub_norms.get_view();
1193  bool is_side = false;
1194  if (int_rule->cv_type == "side")
1195  is_side = true;
1196 
1197  Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_cells,num_ip});
1198  Kokkos::parallel_for("getCubatureCV: weighted measure",policy,KOKKOS_LAMBDA (const int cell,const int ip) {
1199  if (!is_side)
1200  weighted_measure_k(cell,ip) = dyn_phys_cub_weights_k(cell,ip);
1201  for (size_type dim = 0; dim < num_dims; ++dim) {
1202  ip_coordinates_k(cell,ip,dim) = dyn_phys_cub_points_k(cell,ip,dim);
1203  if (is_side)
1204  weighted_normals_k(cell,ip,dim) = dyn_phys_cub_norms_k(cell,ip,dim);
1205  }
1206  });
1207  PHX::Device::execution_space().fence();
1208 }
1209 
1210 template <typename Scalar>
1212 evaluateValuesCV(const PHX::MDField<Scalar, Cell, NODE, Dim>& in_node_coordinates,
1213  const int in_num_cells)
1214 {
1215 
1216  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1217 
1218  size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (size_type) in_num_cells;
1219 
1220  auto s_ref_ip_coordinates = Kokkos::subview(ref_ip_coordinates.get_view(),std::make_pair(0,(int)num_cells),Kokkos::ALL(),Kokkos::ALL());
1221  auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
1222  auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
1223 
1224  cell_tools.mapToReferenceFrame(s_ref_ip_coordinates,
1225  s_ip_coordinates,
1226  s_node_coordinates,
1227  *(int_rule->topology));
1228 
1229  auto s_jac = Kokkos::subview(jac.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1230 
1231  cell_tools.setJacobian(s_jac,
1232  s_ref_ip_coordinates,
1233  s_node_coordinates,
1234  *(int_rule->topology));
1235 
1236  auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1237 
1238  cell_tools.setJacobianInv(s_jac_inv, s_jac);
1239 
1240  auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL());
1241 
1242  cell_tools.setJacobianDet(s_jac_det, s_jac);
1243 }
1244 
1245 #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \
1246  template class IntegrationValues2<SCALAR>;
1247 
1249 
1250 // Disabled FAD support due to long build times on cuda (in debug mode
1251 // it takes multiple hours on some platforms). If we need
1252 // sensitivities wrt coordinates, we can reenable.
1253 
1254 // INTEGRATION_VALUES2_INSTANTIATION(panzer::Traits::FadType)
1255 
1256 }
KOKKOS_INLINE_FUNCTION int cellForSubcell(const int subcell, const int local_cell_index) const
Get the cell for a given subcell and a local_cell_index.
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int num_cells=-1, const Teuchos::RCP< const SubcellConnectivity > &face_connectivity=Teuchos::null)
Evaluate basis values.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
void generateSurfaceCubatureValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells, const SubcellConnectivity &face_connectivity)
This should be a private method, but using lambdas on cuda forces this to be public.
KOKKOS_INLINE_FUNCTION int numSubcells() const
Gives number of subcells (e.g. faces) in connectivity.
PHX::MDField< Scalar, Cell, IP > Array_CellIP
const int & getType() const
Get type of integrator.
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int in_num_cells)
void getCubatureCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
This should be a private method, but using lambdas on cuda forces this to be public.
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
#define PANZER_DOT(a, b)
int getPointOffset(const int subcell_index) const
Returns the integration point offset for a given subcell_index (i.e. local face index) ...
#define PANZER_CROSS(a, b, c)
Teuchos::RCP< Intrepid2::Cubature< PHX::Device::execution_space, double, double > > getIntrepidCubature(const panzer::IntegrationRule &ir) const
Teuchos::RCP< const shards::CellTopology > topology
void evaluateRemainingValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
This should be a private method, but using lambdas on cuda forces this to be public.
ArrayTraits< Scalar, PHX::MDField< Scalar > >::size_type size_type
KOKKOS_INLINE_FUNCTION void convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
void getCubature(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
PHX::MDField< Scalar, Cell, IP, Dim > Array_CellIPDim
PHX::MDField< double > DblArrayDynamic
Teuchos::RCP< shards::CellTopology > side_topology
static void permuteToOther(const PHX::MDField< Scalar, Cell, IP, Dim > &coords, const PHX::MDField< Scalar, Cell, IP, Dim > &other_coords, std::vector< typename ArrayTraits< Scalar, PHX::MDField< Scalar > >::size_type > &permutation)
int numSubcellsOnCellHost(const int cell) const
void setupArraysForNodeRule(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
const int & getOrder() const
Get order of integrator.
KOKKOS_INLINE_FUNCTION int localSubcellForSubcell(const int subcell, const int local_cell_index) const
Get the local subcell index given a subcell and a local cell index.