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