Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
59 
60 // ***********************************************************
61 // * Specializations of setupArrays() for different array types
62 // ***********************************************************
63 
64 namespace panzer {
65 
66 // * Specialization for Kokkos::DynRankView<double,PHX::Device>
67 template <typename Scalar>
70 {
71  MDFieldArrayFactory af(prefix,alloc_arrays);
72 
73  int num_nodes = ir->topology->getNodeCount();
74  int num_cells = ir->workset_size;
75  int num_space_dim = ir->topology->getDimension();
76 
77  int num_ip = 1;
78 
79  dyn_cub_points = af.template buildArray<double,IP,Dim>("cub_points",num_ip, num_space_dim);
80  dyn_cub_weights = af.template buildArray<double,IP>("cub_weights",num_ip);
81 
82  cub_points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
83 
84  if (ir->cv_type == "none" && ir->isSide()) {
85  dyn_side_cub_points = af.template buildArray<double,IP,Dim>("side_cub_points",num_ip, ir->side_topology->getDimension());
86  side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip,ir->side_topology->getDimension());
87  }
88 
89  if (ir->cv_type != "none") {
90  dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>("phys_cub_points",num_cells, num_ip, num_space_dim);
91  dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>("phys_cub_weights",num_cells, num_ip);
92  if (ir->cv_type == "side") {
93  dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>("phys_cub_norms",num_cells, num_ip, num_space_dim);
94  }
95  }
96 
97  dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
98 
99  cub_weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
100 
101  node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
102 
103  jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells, num_ip, num_space_dim,num_space_dim);
104 
105  jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
106 
107  jac_det = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells, num_ip);
108 
109  weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells, num_ip);
110 
111  covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells, num_ip, num_space_dim,num_space_dim);
112 
113  contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
114 
115  norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells, num_ip);
116 
117  ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordiantes",num_cells, num_ip,num_space_dim);
118 
119  ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells, num_ip,num_space_dim);
120 
121  weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted normal",num_cells, num_ip,num_space_dim);
122 
123  surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("surface_normals",num_cells, num_ip,num_space_dim);
124 
125  surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("surface_rotation_matrices",num_cells, num_ip,3,3);
126 
127 }
128 
129 template <typename Scalar>
132 {
133  MDFieldArrayFactory af(prefix,alloc_arrays);
134 
136 
137  int_rule = ir;
138 
139  int num_nodes = ir->topology->getNodeCount();
140  int num_cells = ir->workset_size;
141  int num_space_dim = ir->topology->getDimension();
142 
143  // specialize content if this is quadrature at anode
144  if(num_space_dim==1 && ir->isSide()) {
145  setupArraysForNodeRule(ir);
146  return;
147  }
148 
149  TEUCHOS_ASSERT(ir->getType() != ID::NONE);
150  intrepid_cubature = getIntrepidCubature(*ir);
151 
152  int num_ip = ir->num_points;
153 
154 // Intrepid2::DefaultCubatureFactory cubature_factory;
155 //
156 // if (ir->cv_type == "side")
157 // intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir->topology));
158 //
159 // else if (ir->cv_type == "volume")
160 // intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir->topology));
161 //
162 // else if (ir->cv_type == "boundary" && ir->isSide())
163 // intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir->topology,ir->side));
164 //
165 // else if (ir->cv_type == "none" && ir->isSide())
166 // intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir->side_topology),
167 // ir->cubature_degree);
168 // else
169 // intrepid_cubature = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir->topology),
170 // ir->cubature_degree);
171 // int num_ip = intrepid_cubature->getNumPoints();
172 
173  dyn_cub_points = af.template buildArray<double,IP,Dim>("cub_points",num_ip, num_space_dim);
174  dyn_cub_weights = af.template buildArray<double,IP>("cub_weights",num_ip);
175 
176  cub_points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
177 
178  if (ir->isSide() && ir->cv_type == "none") {
179  dyn_side_cub_points = af.template buildArray<double,IP,Dim>("side_cub_points",num_ip, ir->side_topology->getDimension());
180  side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip,ir->side_topology->getDimension());
181  }
182 
183  if (ir->cv_type != "none") {
184  dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>("phys_cub_points",num_cells, num_ip, num_space_dim);
185  dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>("phys_cub_weights",num_cells, num_ip);
186  if (ir->cv_type == "side") {
187  dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>("phys_cub_norms",num_cells, num_ip, num_space_dim);
188  }
189  }
190 
191  dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
192 
193  cub_weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
194 
195  node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
196 
197  jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells, num_ip, num_space_dim,num_space_dim);
198 
199  jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
200 
201  jac_det = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells, num_ip);
202 
203  weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells, num_ip);
204 
205  covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells, num_ip, num_space_dim,num_space_dim);
206 
207  contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
208 
209  norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells, num_ip);
210 
211  ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordiantes",num_cells, num_ip,num_space_dim);
212 
213  ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells, num_ip,num_space_dim);
214 
215  weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted_normal",num_cells,num_ip,num_space_dim);
216 
217  surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("surface_normals",num_cells, num_ip,num_space_dim);
218 
219  surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("surface_rotation_matrices",num_cells, num_ip,3,3);
220 
221  scratch_for_compute_side_measure =
222  af.template buildStaticArray<Scalar,Point>("scratch_for_compute_side_measure", jac.get_view().span());
223 
224 }
225 
226 template <typename Scalar>
230 {
233 
234  Intrepid2::DefaultCubatureFactory cubature_factory;
235 
236  if(ir.getType() == ID::CV_SIDE){
237  ic = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir.topology));
238  } else if(ir.getType() == ID::CV_VOLUME){
239  ic = Teuchos::rcp(new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir.topology));
240  } else if(ir.getType() == ID::CV_BOUNDARY){
241  ic = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir.topology,ir.getSide()));
242  } else {
243  if(ir.getType() == ID::VOLUME){
244  ic = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir.topology),ir.getOrder());
245  } else if(ir.getType() == ID::SIDE){
246  ic = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir.side_topology),ir.getOrder());
247  } else if(ir.getType() == ID::SURFACE){
248  // closed surface integrals don't exist in intrepid.
249  } else {
250  TEUCHOS_ASSERT(false);
251  }
252  }
253 
254  return ic;
255 }
256 
257 
258 // ***********************************************************
259 // * Evaluation of values - NOT specialized
260 // ***********************************************************
261 template <typename Scalar>
263 evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim> & in_node_coordinates,
264  const int in_num_cells)
265 {
267  const bool is_surface = int_rule->getType() == ID::SURFACE;
268  const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
269 
270  TEUCHOS_ASSERT(not (is_surface and is_cv));
271 
272  if(is_surface){
273  generateSurfaceCubatureValues(in_node_coordinates,in_num_cells);
274  } else if (is_cv) {
275  getCubatureCV(in_node_coordinates, in_num_cells);
276  evaluateValuesCV(in_node_coordinates, in_num_cells);
277  } else {
278  getCubature(in_node_coordinates, in_num_cells);
279  evaluateRemainingValues(in_node_coordinates, in_num_cells);
280  }
281 }
282 
283 template <typename Scalar>
285 getCubature(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
286  const int in_num_cells)
287 {
288 
289  int num_space_dim = int_rule->topology->getDimension();
290  if (int_rule->isSide() && num_space_dim==1) {
291  std::cout << "WARNING: 0-D quadrature rule ifrastructure does not exist!!! Will not be able to do "
292  << "non-natural integration rules.";
293  return;
294  }
295 
296  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
297 
298  if (!int_rule->isSide())
299  intrepid_cubature->getCubature(dyn_cub_points.get_view(), dyn_cub_weights.get_view());
300  else {
301  intrepid_cubature->getCubature(dyn_side_cub_points.get_view(), dyn_cub_weights.get_view());
302 
303  cell_tools.mapToReferenceSubcell(dyn_cub_points.get_view(),
304  dyn_side_cub_points.get_view(),
305  int_rule->spatial_dimension-1,
306  int_rule->side,
307  *(int_rule->topology));
308  }
309 
310  // IP coordinates
311  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
312  auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
313  auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
314  cell_tools.mapToPhysicalFrame(s_ip_coordinates,
315  dyn_cub_points.get_view(),
316  s_in_node_coordinates,
317  *(int_rule->topology));
318 }
319 
320 
321 
322 
323 
324 namespace
325 {
326 
327 template <typename array_t, typename scalar_t>
328 class point_sorter_t
329 {
330 public:
331 
332  point_sorter_t() = delete;
333  point_sorter_t(const array_t & array, const int cell, const int offset):
334  _array(array),
335  _cell(cell),
336  _offset(offset),
337  _rel_tol(1.e-12)
338  {
339  _num_dims=_array.extent(2);
340  }
341 
342 
343  // This needs to be optimized
344  bool operator()(const int & point_a, const int & point_b) const
345  {
346 
347  if(_num_dims==1){
348 
349  const scalar_t & x_a = _array(_cell,_offset+point_a,0);
350  const scalar_t & x_b = _array(_cell,_offset+point_b,0);
351 
352  const scalar_t rel = std::max(std::fabs(x_a),std::fabs(x_b));
353 
354  return test_less(x_a,x_b,rel);
355 
356  } else if(_num_dims==2){
357 
358  const scalar_t & x_a = _array(_cell,_offset+point_a,0);
359  const scalar_t & x_b = _array(_cell,_offset+point_b,0);
360 
361  const scalar_t & y_a = _array(_cell,_offset+point_a,1);
362  const scalar_t & y_b = _array(_cell,_offset+point_b,1);
363 
364  const scalar_t rel_x = std::max(std::fabs(x_a),std::fabs(x_b));
365  const scalar_t rel_y = std::max(std::fabs(y_a),std::fabs(y_b));
366  const scalar_t rel = std::max(rel_x,rel_y);
367 
368  if(test_eq(y_a,y_b,rel)){
369  if(test_less(x_a,x_b,rel)){
370  // Sort by x
371  return true;
372  }
373  } else if(test_less(y_a,y_b,rel)){
374  // Sort by y
375  return true;
376  }
377 
378  // Otherwise b < a
379  return false;
380 
381  } else if(_num_dims==3){
382 
383  const scalar_t & x_a = _array(_cell,_offset+point_a,0);
384  const scalar_t & x_b = _array(_cell,_offset+point_b,0);
385 
386  const scalar_t & y_a = _array(_cell,_offset+point_a,1);
387  const scalar_t & y_b = _array(_cell,_offset+point_b,1);
388 
389  const scalar_t & z_a = _array(_cell,_offset+point_a,2);
390  const scalar_t & z_b = _array(_cell,_offset+point_b,2);
391 
392  const scalar_t rel_x = std::max(std::fabs(x_a),std::fabs(x_b));
393  const scalar_t rel_y = std::max(std::fabs(y_a),std::fabs(y_b));
394  const scalar_t rel_z = std::max(std::fabs(z_a),std::fabs(z_b));
395  const scalar_t rel = std::max(rel_x,std::max(rel_y,rel_z));
396 
397  if(test_less(z_a,z_b,rel)){
398  // Sort by z
399  return true;
400  } else if(test_eq(z_a,z_b,rel)){
401  if(test_eq(y_a,y_b,rel)){
402  if(test_less(x_a,x_b,rel)){
403  // Sort by x
404  return true;
405  }
406  } else if(test_less(y_a,y_b,rel)){
407  // Sort by y
408  return true;
409  }
410  }
411  // Otherwise b < a
412  return false;
413 
414  } else {
415  TEUCHOS_ASSERT(false);
416  }
417  }
418 
419 protected:
420 
421  bool
422  test_eq(const scalar_t & a, const scalar_t & b, const scalar_t & rel) const
423  {
424  if(rel==0){
425  return true;
426  }
427  return std::fabs(a-b) < _rel_tol * rel;
428  }
429 
430  bool
431  test_less(const scalar_t & a, const scalar_t & b, const scalar_t & rel) const
432  {
433  if(rel==0){
434  return false;
435  }
436  return (a-b) < -_rel_tol * rel;
437  }
438 
439  const array_t & _array;
440  int _cell;
441  int _offset;
443  scalar_t _rel_tol;
444 
445 };
446 
447 template<typename T>
448 void
449 convertNormalToRotationMatrix(const T normal[3], T transverse[3], T binormal[3])
450 {
451 
452  const T n = sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
453 
454  // If this fails then the geometry for this cell is probably undefined
455  if(n > 0.){
456 
457 
458  // Make sure transverse is not parallel to normal within some margin of error
459  transverse[0]=0.;transverse[1]=1.;transverse[2]=0.;
460  if(std::fabs(normal[0]*transverse[0]+normal[1]*transverse[1])>0.9){
461  transverse[0]=1.;transverse[1]=0.;
462  }
463 
464  const T nt = normal[0]*transverse[0]+normal[1]*transverse[1]+normal[2]*transverse[2];
465 
466  // Note normal has unit length
467  const T mult = nt/(n*n); // = nt
468 
469  // Remove normal projection from transverse
470  for(int dim=0;dim<3;++dim){
471  transverse[dim] = transverse[dim] - mult * normal[dim];
472  }
473 
474  const T t = sqrt(transverse[0]*transverse[0]+transverse[1]*transverse[1]+transverse[2]*transverse[2]);
475  TEUCHOS_ASSERT(t != 0.);
476  for(int dim=0;dim<3;++dim){
477  transverse[dim] /= t;
478  }
479 
480  // We assume a right handed system such that b = n \times t
481  binormal[0] = (normal[1] * transverse[2] - normal[2] * transverse[1]);
482  binormal[1] = (normal[2] * transverse[0] - normal[0] * transverse[2]);
483  binormal[2] = (normal[0] * transverse[1] - normal[1] * transverse[0]);
484 
485  // Normalize binormal
486  const T b = sqrt(binormal[0]*binormal[0]+binormal[1]*binormal[1]+binormal[2]*binormal[2]);
487  for(int dim=0;dim<3;++dim){
488  binormal[dim] /= b;
489  }
490  } else {
491  transverse[0] = 0.;
492  transverse[1] = 0.;
493  transverse[2] = 0.;
494  binormal[0] = 0.;
495  binormal[1] = 0.;
496  binormal[2] = 0.;
497 
498  // TEUCHOS_ASSERT(false);
499  }
500 
501 }
502 
503 }
504 
505 template <typename Scalar>
508  int a,
509  int b)
510 {
511  const int new_cell_point = a;
512  const int old_cell_point = b;
513 
514  const int cell_dim = ref_ip_coordinates.extent(2);
515 
516  Scalar hold;
517 
518  hold = weighted_measure(cell,new_cell_point);
519  weighted_measure(cell,new_cell_point) = weighted_measure(cell,old_cell_point);
520  weighted_measure(cell,old_cell_point) = hold;
521 
522  hold = jac_det(cell,new_cell_point);
523  jac_det(cell,new_cell_point) = jac_det(cell,old_cell_point);
524  jac_det(cell,old_cell_point) = hold;
525 
526  for(int dim=0;dim<cell_dim;++dim){
527 
528  hold = ref_ip_coordinates(cell,new_cell_point,dim);
529  ref_ip_coordinates(cell,new_cell_point,dim) = ref_ip_coordinates(cell,old_cell_point,dim);
530  ref_ip_coordinates(cell,old_cell_point,dim) = hold;
531 
532  hold = ip_coordinates(cell,new_cell_point,dim);
533  ip_coordinates(cell,new_cell_point,dim) = ip_coordinates(cell,old_cell_point,dim);
534  ip_coordinates(cell,old_cell_point,dim) = hold;
535 
536  hold = surface_normals(cell,new_cell_point,dim);
537  surface_normals(cell,new_cell_point,dim) = surface_normals(cell,old_cell_point,dim);
538  surface_normals(cell,old_cell_point,dim) = hold;
539 
540  for(int dim2=0;dim2<cell_dim;++dim2){
541 
542  hold = jac(cell,new_cell_point,dim,dim2);
543  jac(cell,new_cell_point,dim,dim2) = jac(cell,old_cell_point,dim,dim2);
544  jac(cell,old_cell_point,dim,dim2) = hold;
545 
546  hold = jac_inv(cell,new_cell_point,dim,dim2);
547  jac_inv(cell,new_cell_point,dim,dim2) = jac_inv(cell,old_cell_point,dim,dim2);
548  jac_inv(cell,old_cell_point,dim,dim2) = hold;
549  }
550  }
551 }
552 
553 template <typename Scalar>
556  int cell,
557  int offset,
558  std::vector<int> & order)
559 {
560  for(size_t point_index=0;point_index<order.size();++point_index){
561  order[point_index] = point_index;
562  }
563 
564  // We need to sort the indexes in point_indexes by their ip_coordinate's position in space.
565  // We will then use that to sort all of our arrays.
566 
567  point_sorter_t<Array_CellIPDim,Scalar> sorter(coords,cell,offset);
568  std::sort(order.begin(),order.end(),sorter);
569 }
570 
571 template <typename Scalar>
573 generateSurfaceCubatureValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
574  const int in_num_cells)
575 {
576 
577  TEUCHOS_ASSERT(int_rule->getType() == IntegrationDescriptor::SURFACE);
578 
579  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
580 
581  const shards::CellTopology & cell_topology = *(int_rule->topology);
582  const panzer::IntegrationRule & ir = *int_rule;
583 
584  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
585 
586  // Copy over coordinates
587  {
588  const int num_nodes = in_node_coordinates.extent(1);
589  const int num_dims = in_node_coordinates.extent(2);
590 
591  for(int cell=0; cell<num_cells; ++cell){
592  for(int node=0; node<num_nodes; ++node){
593  for(int dim=0; dim<num_dims; ++dim){
594  node_coordinates(cell,node,dim) = in_node_coordinates(cell,node,dim);
595  }
596  }
597  }
598  }
599 
600  // NOTE: We are assuming that each face can have a different number of points.
601  // Not sure if this is necessary, but it requires a lot of additional allocations
602 
603  const int cell_dim = cell_topology.getDimension();
604  const int subcell_dim = cell_topology.getDimension()-1;
605  const int num_subcells = cell_topology.getSubcellCount(subcell_dim);
606 
607  Intrepid2::DefaultCubatureFactory cubature_factory;
608 
609  // We get to build up our cubature one face at a time
610  {
611  int point_offset=0;
612  for(int subcell_index=0; subcell_index<num_subcells; ++subcell_index) {
613 
614  // Default for 1D
615  int num_points_on_face = 1;
616 
617  // Get the cubature for the side
618  Kokkos::DynRankView<double,PHX::Device> tmp_side_cub_weights;
619  Kokkos::DynRankView<double,PHX::Device> tmp_side_cub_points;
620  if(cell_dim==1){
621  tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>("tmp_side_cub_weights",num_points_on_face);
622  tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>("cell_tmp_side_cub_points",num_points_on_face,cell_dim);
623  tmp_side_cub_weights(0)=1.;
624  tmp_side_cub_points(0,0) = (subcell_index==0)? -1. : 1.;
625  } else {
626 
627  // Get the face topology from the cell topology
628  const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,subcell_index));
629 
630  auto ic = cubature_factory.create<PHX::Device::execution_space,double,double>(face_topology,ir.getOrder());
631  num_points_on_face = ic->getNumPoints();
632 
633  tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>("tmp_side_cub_weights",num_points_on_face);
634  tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>("cell_tmp_side_cub_points",num_points_on_face,cell_dim);
635 
636  auto subcell_cub_points = Kokkos::DynRankView<double,PHX::Device>("subcell_cub_points",num_points_on_face,subcell_dim);
637 
638  // Get the reference face points
639  ic->getCubature(subcell_cub_points, tmp_side_cub_weights);
640 
641  // Convert from reference face points to reference cell points
642  cell_tools.mapToReferenceSubcell(tmp_side_cub_points,
643  subcell_cub_points,
644  subcell_dim,
645  subcell_index,
646  cell_topology);
647  }
648 
649 
650  for(int local_point=0;local_point<num_points_on_face;++local_point){
651  const int point = point_offset + local_point;
652  for(int dim=0;dim<cell_dim;++dim){
653  cub_points(point,dim) = tmp_side_cub_points(local_point,dim);
654  }
655  }
656 
657 
658  // Map from side points to physical points
659  auto side_ip_coordinates = Kokkos::DynRankView<Scalar,PHX::Device>("side_ip_coordinates",num_cells,num_points_on_face,cell_dim);
660  auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL);
661  cell_tools.mapToPhysicalFrame(side_ip_coordinates,
662  tmp_side_cub_points,
663  s_node_coordinates,
664  cell_topology);
665 
666  // Create a jacobian and his friends for this side
667  auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
668  cell_tools.setJacobian(side_jacobian,
669  tmp_side_cub_points,
670  s_node_coordinates,
671  cell_topology);
672 
673  auto side_inverse_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_inv_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
674  cell_tools.setJacobianInv(side_inverse_jacobian, side_jacobian);
675 
676  auto side_det_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_det_jac",num_cells,num_points_on_face);
677  cell_tools.setJacobianDet(side_det_jacobian, side_jacobian);
678 
679  // Calculate measures (quadrature weights in physical space) for this side
680  auto side_weighted_measure = Kokkos::DynRankView<Scalar,PHX::Device>("side_weighted_measure",num_cells,num_points_on_face);
681  if(cell_dim == 1){
682  Kokkos::deep_copy(side_weighted_measure, tmp_side_cub_weights(0));
683  } else if(cell_dim == 2){
684  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
685  computeEdgeMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights,
686  subcell_index,cell_topology,
687  scratch_for_compute_side_measure.get_view());
688  } else if(cell_dim == 3){
689  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
690  computeFaceMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights,
691  subcell_index,cell_topology,
692  scratch_for_compute_side_measure.get_view());
693  }
694 
695  // Calculate normals
696  auto side_normals = Kokkos::DynRankView<Scalar,PHX::Device>("side_normals",num_cells,num_points_on_face,cell_dim);
697  if(cell_dim == 1){
698 
699  int other_subcell_index = (subcell_index==0) ? 1 : 0;
700 
701  for(int cell=0;cell<num_cells;++cell){
702  Scalar norm = (in_node_coordinates(cell,subcell_index,0) - in_node_coordinates(cell,other_subcell_index,0));
703  side_normals(cell,0,0) = norm / fabs(norm);
704  }
705 
706  } else {
707 
708  cell_tools.getPhysicalSideNormals(side_normals,side_jacobian,subcell_index,cell_topology);
709 
710  // Normalize each normal
711  for(int cell=0;cell<num_cells;++cell){
712  for(int point=0;point<num_points_on_face;++point){
713  Scalar n = 0.;
714  for(int dim=0;dim<cell_dim;++dim){
715  n += side_normals(cell,point,dim)*side_normals(cell,point,dim);
716  }
717  // If n is zero then this is - hopefully - a virtual cell
718  if(n > 0.){
719  n = std::sqrt(n);
720  for(int dim=0;dim<cell_dim;++dim){
721  side_normals(cell,point,dim) /= n;
722  }
723  }
724  }
725  }
726 
727  }
728 
729 
730  // Now that we have all these wonderful values, lets copy them to the actual arrays
731  for(int cell=0;cell<num_cells;++cell){
732  for(int side_point=0; side_point<num_points_on_face;++side_point){
733  const int cell_point = point_offset + side_point;
734 
735  weighted_measure(cell,cell_point) = side_weighted_measure(cell,side_point);
736  jac_det(cell,cell_point) = side_det_jacobian(cell,side_point);
737  for(int dim=0;dim<cell_dim;++dim){
738  ref_ip_coordinates(cell,cell_point,dim) = cub_points(cell_point,dim);
739  ip_coordinates(cell,cell_point,dim) = side_ip_coordinates(cell,side_point,dim);
740  surface_normals(cell,cell_point,dim) = side_normals(cell,side_point,dim);
741 
742  for(int dim2=0;dim2<cell_dim;++dim2){
743  jac(cell,cell_point,dim,dim2) = side_jacobian(cell,side_point,dim,dim2);
744  jac_inv(cell,cell_point,dim,dim2) = side_inverse_jacobian(cell,side_point,dim,dim2);
745  }
746  }
747  }
748  }
749  point_offset += num_points_on_face;
750  }
751  }
752 
753  // Now we need to sort the cubature points for each face so that they will line up between cells
754  {
755  for(int subcell_index=0; subcell_index<num_subcells;++subcell_index){
756 
757  const int point_offset = ir.getPointOffset(subcell_index);
758  const int num_points_on_face = ir.getPointOffset(subcell_index+1) - point_offset;
759  std::vector<int> point_indexes(num_points_on_face,-1);
760 
761  for(int cell=0; cell<num_cells; ++cell){
762 
763  // build a point index array based on point coordinates
764  uniqueCoordOrdering(ip_coordinates,cell,point_offset,point_indexes);
765 
766  // Indexes are now sorted, now we swap everything around
767  reorder(point_indexes,[=](int a,int b) { swapQuadraturePoints(cell,point_offset+a,point_offset+b); });
768  }
769  }
770  }
771 
772  // We also need surface rotation matrices
773  const int num_points = ir.getPointOffset(num_subcells);
774  Scalar normal[3];
775  Scalar transverse[3];
776  Scalar binormal[3];
777  for(int i=0;i<3;i++){normal[i]=0.;}
778  for(int i=0;i<3;i++){transverse[i]=0.;}
779  for(int i=0;i<3;i++){binormal[i]=0.;}
780  for(int cell=0; cell<num_cells; ++cell){
781  for(int subcell_index=0; subcell_index<num_subcells; ++subcell_index){
782  for(int point=0; point<num_points; ++point){
783 
784  for(int dim=0; dim<3; ++dim)
785  normal[dim] = 0.0;
786 
787  for(int dim=0; dim<cell_dim; ++dim){
788  normal[dim] = surface_normals(cell,point,dim);
789  }
790 
791  convertNormalToRotationMatrix<Scalar>(normal,transverse,binormal);
792 
793  for(int dim=0; dim<3; ++dim){
794  surface_rotation_matrices(cell,point,0,dim) = normal[dim];
795  surface_rotation_matrices(cell,point,1,dim) = transverse[dim];
796  surface_rotation_matrices(cell,point,2,dim) = binormal[dim];
797  }
798  }
799  }
800  }
801 
802 
803  // I'm not sure if these should exist for surface integrals, but here we go!
804 
805  // Shakib contravarient metric tensor
806  for (int cell = 0; cell < num_cells; ++cell) {
807  for (size_type ip = 0; ip < contravarient.extent(1); ++ip) {
808 
809  // zero out matrix
810  for (size_type i = 0; i < contravarient.extent(2); ++i)
811  for (size_type j = 0; j < contravarient.extent(3); ++j)
812  covarient(cell,ip,i,j) = 0.0;
813 
814  // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha}
815  for (size_type i = 0; i < contravarient.extent(2); ++i) {
816  for (size_type j = 0; j < contravarient.extent(3); ++j) {
817  for (size_type alpha = 0; alpha < contravarient.extent(2); ++alpha) {
818  covarient(cell,ip,i,j) += jac(cell,ip,i,alpha) * jac(cell,ip,j,alpha);
819  }
820  }
821  }
822 
823  }
824  }
825 
826  auto s_contravarient = Kokkos::subview(contravarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
827  auto s_covarient = Kokkos::subview(covarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
828  Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
829 
830  // norm of g_ij
831  for (int cell = 0; cell < num_cells; ++cell) {
832  for (size_type ip = 0; ip < contravarient.extent(1); ++ip) {
833  norm_contravarient(cell,ip) = 0.0;
834  for (size_type i = 0; i < contravarient.extent(2); ++i) {
835  for (size_type j = 0; j < contravarient.extent(3); ++j) {
836  norm_contravarient(cell,ip) += contravarient(cell,ip,i,j) * contravarient(cell,ip,i,j);
837  }
838  }
839  norm_contravarient(cell,ip) = std::sqrt(norm_contravarient(cell,ip));
840  }
841  }
842 
843 }
844 
845 
846 template <typename Scalar>
848 evaluateRemainingValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
849  const int in_num_cells)
850 {
851  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
852 
853  // copy the dynamic data structures into the static data structures
854  {
855  size_type num_ip = dyn_cub_points.extent(0);
856  size_type num_dims = dyn_cub_points.extent(1);
857 
858  for (size_type ip = 0; ip < num_ip; ++ip) {
859  cub_weights(ip) = dyn_cub_weights(ip);
860  for (size_type dim = 0; dim < num_dims; ++dim)
861  cub_points(ip,dim) = dyn_cub_points(ip,dim);
862  }
863  }
864 
865  if (int_rule->isSide()) {
866  const size_type num_ip = dyn_cub_points.extent(0), num_side_dims = dyn_side_cub_points.extent(1);
867  for (size_type ip = 0; ip < num_ip; ++ip)
868  for (size_type dim = 0; dim < num_side_dims; ++dim)
869  side_cub_points(ip,dim) = dyn_side_cub_points(ip,dim);
870  }
871 
872  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
873 
874  {
875  size_type num_nodes = in_node_coordinates.extent(1);
876  size_type num_dims = in_node_coordinates.extent(2);
877 
878  for (int cell = 0; cell < num_cells; ++cell) {
879  for (size_type node = 0; node < num_nodes; ++node) {
880  for (size_type dim = 0; dim < num_dims; ++dim) {
881  node_coordinates(cell,node,dim) =
882  in_node_coordinates(cell,node,dim);
883  }
884  }
885  }
886  }
887 
888  auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
889  auto s_jac = Kokkos::subview(jac.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
890  cell_tools.setJacobian(jac.get_view(),
891  cub_points.get_view(),
892  node_coordinates.get_view(),
893  *(int_rule->topology));
894 
895  auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
896  cell_tools.setJacobianInv(s_jac_inv, s_jac);
897 
898  auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair(0,num_cells),Kokkos::ALL());
899  cell_tools.setJacobianDet(s_jac_det, s_jac);
900 
901  auto s_weighted_measure = Kokkos::subview(weighted_measure.get_view(),std::make_pair(0,num_cells),Kokkos::ALL());
902  if (!int_rule->isSide()) {
903  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
904  computeCellMeasure(s_weighted_measure, s_jac_det, cub_weights.get_view());
905  }
906  else if(int_rule->spatial_dimension==3) {
907  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
908  computeFaceMeasure(s_weighted_measure, s_jac, cub_weights.get_view(),
909  int_rule->side, *int_rule->topology,
910  scratch_for_compute_side_measure.get_view());
911  }
912  else if(int_rule->spatial_dimension==2) {
913  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
914  computeEdgeMeasure(s_weighted_measure, s_jac, cub_weights.get_view(),
915  int_rule->side,*int_rule->topology,
916  scratch_for_compute_side_measure.get_view());
917  }
918  else TEUCHOS_ASSERT(false);
919 
920  // Shakib contravarient metric tensor
921  for (int cell = 0; cell < num_cells; ++cell) {
922  for (size_type ip = 0; ip < contravarient.extent(1); ++ip) {
923 
924  // zero out matrix
925  for (size_type i = 0; i < contravarient.extent(2); ++i)
926  for (size_type j = 0; j < contravarient.extent(3); ++j)
927  covarient(cell,ip,i,j) = 0.0;
928 
929  // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha}
930  for (size_type i = 0; i < contravarient.extent(2); ++i) {
931  for (size_type j = 0; j < contravarient.extent(3); ++j) {
932  for (size_type alpha = 0; alpha < contravarient.extent(2); ++alpha) {
933  covarient(cell,ip,i,j) += jac(cell,ip,i,alpha) * jac(cell,ip,j,alpha);
934  }
935  }
936  }
937 
938  }
939  }
940 
941  auto s_covarient = Kokkos::subview(covarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
942  auto s_contravarient = Kokkos::subview(contravarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
943  Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
944 
945  // norm of g_ij
946  for (int cell = 0; cell < num_cells; ++cell) {
947  for (size_type ip = 0; ip < contravarient.extent(1); ++ip) {
948  norm_contravarient(cell,ip) = 0.0;
949  for (size_type i = 0; i < contravarient.extent(2); ++i) {
950  for (size_type j = 0; j < contravarient.extent(3); ++j) {
951  norm_contravarient(cell,ip) += contravarient(cell,ip,i,j) * contravarient(cell,ip,i,j);
952  }
953  }
954  norm_contravarient(cell,ip) = std::sqrt(norm_contravarient(cell,ip));
955  }
956  }
957 }
958 
959 // Find the permutation that maps the set of points coords to other_coords. To
960 // avoid possible finite precision issues, == is not used, but rather
961 // min(norm(.)).
962 template <typename Scalar>
963 static void
964 permuteToOther(const PHX::MDField<Scalar,Cell,IP,Dim>& coords,
965  const PHX::MDField<Scalar,Cell,IP,Dim>& other_coords,
966  std::vector<typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type>& permutation)
967 {
968  typedef typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type size_type;
969  // We can safely assume: (1) The permutation is the same for every cell in
970  // the workset. (2) The first workset has valid data. Hence we operate only
971  // on cell 0.
972  const size_type cell = 0;
973  const size_type num_ip = coords.extent(1), num_dim = coords.extent(2);
974  permutation.resize(num_ip);
975  std::vector<char> taken(num_ip, 0);
976  for (size_type ip = 0; ip < num_ip; ++ip) {
977  // Find an other point to associate with ip.
978  size_type i_min = 0;
979  Scalar d_min = -1;
980  for (size_type other_ip = 0; other_ip < num_ip; ++other_ip) {
981  // For speed, skip other points that are already associated.
982  if (taken[other_ip]) continue;
983  // Compute the distance between the two points.
984  Scalar d(0);
985  for (size_type dim = 0; dim < num_dim; ++dim) {
986  const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim);
987  d += diff*diff;
988  }
989  if (d_min < 0 || d < d_min) {
990  d_min = d;
991  i_min = other_ip;
992  }
993  }
994  // Record the match.
995  permutation[ip] = i_min;
996  // This point is now taken.
997  taken[i_min] = 1;
998  }
999 }
1000 
1001 template <typename Scalar>
1003 evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
1004  const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates,
1005  const int in_num_cells)
1006 {
1007  const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
1008 
1009  if (int_rule->cv_type == "none") {
1010 
1011  getCubature(in_node_coordinates, in_num_cells);
1012 
1013  {
1014  // Determine the permutation.
1015  std::vector<size_type> permutation(other_ip_coordinates.extent(1));
1016  permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
1017  // Apply the permutation to the cubature arrays.
1018  MDFieldArrayFactory af(prefix, alloc_arrays);
1019  {
1020  const size_type num_ip = dyn_cub_points.extent(0);
1021  {
1022  const size_type num_dim = dyn_side_cub_points.extent(1);
1023  DblArrayDynamic old_dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
1024  "old_dyn_side_cub_points", num_ip, num_dim);
1025  old_dyn_side_cub_points.deep_copy(dyn_side_cub_points);
1026  for (size_type ip = 0; ip < num_ip; ++ip)
1027  if (ip != permutation[ip])
1028  for (size_type dim = 0; dim < num_dim; ++dim)
1029  dyn_side_cub_points(ip, dim) = old_dyn_side_cub_points(permutation[ip], dim);
1030  }
1031  {
1032  const size_type num_dim = dyn_cub_points.extent(1);
1033  DblArrayDynamic old_dyn_cub_points = af.template buildArray<double,IP,Dim>(
1034  "old_dyn_cub_points", num_ip, num_dim);
1035  old_dyn_cub_points.deep_copy(dyn_cub_points);
1036  for (size_type ip = 0; ip < num_ip; ++ip)
1037  if (ip != permutation[ip])
1038  for (size_type dim = 0; dim < num_dim; ++dim)
1039  dyn_cub_points(ip, dim) = old_dyn_cub_points(permutation[ip], dim);
1040  }
1041  {
1042  DblArrayDynamic old_dyn_cub_weights = af.template buildArray<double,IP>(
1043  "old_dyn_cub_weights", num_ip);
1044  old_dyn_cub_weights.deep_copy(dyn_cub_weights);
1045  for (size_type ip = 0; ip < dyn_cub_weights.extent(0); ++ip)
1046  if (ip != permutation[ip])
1047  dyn_cub_weights(ip) = old_dyn_cub_weights(permutation[ip]);
1048  }
1049  }
1050  {
1051  const size_type num_ip = ip_coordinates.extent(1);
1052  const size_type num_dim = ip_coordinates.extent(2);
1053  Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1054  "old_ip_coordinates", ip_coordinates.extent(0), num_ip, num_dim);
1055  Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
1056  for (int cell = 0; cell < num_cells; ++cell)
1057  for (size_type ip = 0; ip < num_ip; ++ip)
1058  if (ip != permutation[ip])
1059  for (size_type dim = 0; dim < num_dim; ++dim)
1060  ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
1061  }
1062  // All subsequent calculations inherit the permutation.
1063  }
1064 
1065  evaluateRemainingValues(in_node_coordinates, in_num_cells);
1066  }
1067 
1068  else {
1069 
1070  getCubatureCV(in_node_coordinates, in_num_cells);
1071 
1072  // Determine the permutation.
1073  std::vector<size_type> permutation(other_ip_coordinates.extent(1));
1074  permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
1075 
1076  // Apply the permutation to the cubature arrays.
1077  MDFieldArrayFactory af(prefix, alloc_arrays);
1078  {
1079  const size_type workset_size = ip_coordinates.extent(0), num_ip = ip_coordinates.extent(1),
1080  num_dim = ip_coordinates.extent(2);
1081  Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1082  "old_ip_coordinates", workset_size, num_ip, num_dim);
1083  Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
1084  Array_CellIPDim old_weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1085  "old_weighted_normals", workset_size, num_ip, num_dim);
1086  Array_CellIP old_weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
1087  "old_weighted_measure", workset_size, num_ip);
1088  if (int_rule->cv_type == "side")
1089  Kokkos::deep_copy(old_weighted_normals.get_static_view(), weighted_normals.get_static_view());
1090  else
1091  Kokkos::deep_copy(old_weighted_measure.get_static_view(), weighted_measure.get_static_view());
1092  for (int cell = 0; cell < num_cells; ++cell)
1093  {
1094  for (size_type ip = 0; ip < num_ip; ++ip)
1095  {
1096  if (ip != permutation[ip]) {
1097  if (int_rule->cv_type == "boundary" || int_rule->cv_type == "volume")
1098  weighted_measure(cell, ip) = old_weighted_measure(cell, permutation[ip]);
1099  for (size_type dim = 0; dim < num_dim; ++dim)
1100  {
1101  ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
1102  if (int_rule->cv_type == "side")
1103  weighted_normals(cell, ip, dim) = old_weighted_normals(cell, permutation[ip], dim);
1104 
1105  }
1106  }
1107  }
1108  }
1109  }
1110 
1111  evaluateValuesCV(in_node_coordinates, in_num_cells);
1112  }
1113 }
1114 
1115 template <typename Scalar>
1117 getCubatureCV(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
1118  const int in_num_cells)
1119 {
1120  int num_space_dim = int_rule->topology->getDimension();
1121  if (int_rule->isSide() && num_space_dim==1) {
1122  std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
1123  << "non-natural integration rules.";
1124  return;
1125  }
1126 
1127  size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (size_type) in_num_cells;
1128  std::pair<int,int> cell_range(0,num_cells);
1129  {
1130  size_type num_nodes = in_node_coordinates.extent(1);
1131  size_type num_dims = in_node_coordinates.extent(2);
1132 
1133  for (size_type cell = 0; cell < num_cells; ++cell) {
1134  for (size_type node = 0; node < num_nodes; ++node) {
1135  for (size_type dim = 0; dim < num_dims; ++dim) {
1136  node_coordinates(cell,node,dim) =
1137  in_node_coordinates(cell,node,dim);
1138  dyn_node_coordinates(cell,node,dim) =
1139  Sacado::ScalarValue<Scalar>::eval(in_node_coordinates(cell,node,dim));
1140  }
1141  }
1142  }
1143  }
1144 
1145  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());
1146  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());
1147  if (int_rule->cv_type == "side") {
1148  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());
1149  intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_norms,s_dyn_node_coordinates);
1150  }
1151  else {
1152  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());
1153  intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_weights,s_dyn_node_coordinates);
1154  }
1155 
1156  // size_type num_cells = dyn_phys_cub_points.extent(0);
1157  size_type num_ip =dyn_phys_cub_points.extent(1);
1158  size_type num_dims = dyn_phys_cub_points.extent(2);
1159 
1160  for (size_type cell = 0; cell < num_cells; ++cell) {
1161  for (size_type ip = 0; ip < num_ip; ++ip) {
1162  if (int_rule->cv_type != "side")
1163  weighted_measure(cell,ip) = dyn_phys_cub_weights(cell,ip);
1164  for (size_type dim = 0; dim < num_dims; ++dim) {
1165  ip_coordinates(cell,ip,dim) = dyn_phys_cub_points(cell,ip,dim);
1166  if (int_rule->cv_type == "side")
1167  weighted_normals(cell,ip,dim) = dyn_phys_cub_norms(cell,ip,dim);
1168  }
1169  }
1170  }
1171 
1172 }
1173 
1174 template <typename Scalar>
1176 evaluateValuesCV(const PHX::MDField<Scalar, Cell, NODE, Dim>& in_node_coordinates,
1177  const int in_num_cells)
1178 {
1179 
1180  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1181 
1182  size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (size_type) in_num_cells;
1183 
1184  auto s_ref_ip_coordinates = Kokkos::subview(ref_ip_coordinates.get_view(),std::make_pair(0,(int)num_cells),Kokkos::ALL(),Kokkos::ALL());
1185  auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
1186  auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
1187 
1188  cell_tools.mapToReferenceFrame(s_ref_ip_coordinates,
1189  s_ip_coordinates,
1190  s_node_coordinates,
1191  *(int_rule->topology));
1192 
1193  auto s_jac = Kokkos::subview(jac.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1194 
1195  cell_tools.setJacobian(s_jac,
1196  s_ref_ip_coordinates,
1197  s_node_coordinates,
1198  *(int_rule->topology));
1199 
1200  auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1201 
1202  cell_tools.setJacobianInv(s_jac_inv, s_jac);
1203 
1204  auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL());
1205 
1206  cell_tools.setJacobianDet(s_jac_det, s_jac);
1207 }
1208 
1209 #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \
1210  template class IntegrationValues2<SCALAR>;
1211 
1213 INTEGRATION_VALUES2_INSTANTIATION(panzer::Traits::FadType)
1214 
1215 }
static void uniqueCoordOrdering(Array_CellIPDim &coords, int cell, int offset, std::vector< int > &order)
Using coordinate build an arrray that specifies a unique ordering.
PHX::MDField< Scalar, Cell, IP > Array_CellIP
const array_t & _array
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int in_num_cells)
const int & getType() const
Get type of integrator.
void generateSurfaceCubatureValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
void swapQuadraturePoints(int cell, int a, int b)
Swap the ordering of quadrature points in a specified cell.
void getCubatureCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
int getPointOffset(const int subcell_index) const
Returns the integration point offset for a given subcell_index (i.e. local face index) ...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int num_cells=-1)
Evaluate basis values.
Teuchos::RCP< const shards::CellTopology > topology
void evaluateRemainingValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
ArrayTraits< Scalar, PHX::MDField< Scalar > >::size_type size_type
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
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
void reorder(std::vector< int > &order, std::function< void(int, int)> swapper)
Using a functor, reorder an array using a order vector.
Teuchos::RCP< Intrepid2::Cubature< PHX::Device::execution_space, double, double > > getIntrepidCubature(const panzer::IntegrationRule &ir) const
PHX::MDField< double > DblArrayDynamic
Teuchos::RCP< shards::CellTopology > side_topology
#define TEUCHOS_ASSERT(assertion_test)
scalar_t _rel_tol
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)
void setupArraysForNodeRule(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
const int & getOrder() const
Get order of integrator.