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 
45 #include "Shards_CellTopology.hpp"
46 
47 #include "Kokkos_DynRankView.hpp"
48 #include "Intrepid2_FunctionSpaceTools.hpp"
49 #include "Intrepid2_RealSpaceTools.hpp"
50 #include "Intrepid2_CellTools.hpp"
51 #include "Intrepid2_ArrayTools.hpp"
52 #include "Intrepid2_CubatureControlVolume.hpp"
53 #include "Intrepid2_CubatureControlVolumeSide.hpp"
54 #include "Intrepid2_CubatureControlVolumeBoundary.hpp"
55 
57 #include "Panzer_Traits.hpp"
58 
59 // ***********************************************************
60 // * Specializations of setupArrays() for different array types
61 // ***********************************************************
62 
63 namespace panzer {
64 
65  // * Specialization for Kokkos::DynRankView<double,PHX::Device>
66  template <typename Scalar>
68  setupArraysForNodeRule(const Teuchos::RCP<const panzer::IntegrationRule>& ir)
69  {
70  MDFieldArrayFactory af(prefix,alloc_arrays);
71 
72  int num_nodes = ir->topology->getNodeCount();
73  int num_cells = ir->workset_size;
74  int num_space_dim = ir->topology->getDimension();
75 
76  int num_ip = 1;
77 
78  dyn_cub_points = af.template buildArray<double,IP,Dim>("cub_points",num_ip, num_space_dim);
79  dyn_cub_weights = af.template buildArray<double,IP>("cub_weights",num_ip);
80 
81  cub_points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
82 
83  if (ir->cv_type == "none" && ir->isSide()) {
84  dyn_side_cub_points = af.template buildArray<double,IP,Dim>("side_cub_points",num_ip, ir->side_topology->getDimension());
85  side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip,ir->side_topology->getDimension());
86  }
87 
88  if (ir->cv_type != "none") {
89  dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>("phys_cub_points",num_cells, num_ip, num_space_dim);
90  dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>("phys_cub_weights",num_cells, num_ip);
91  if (ir->cv_type == "side") {
92  dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>("phys_cub_norms",num_cells, num_ip, num_space_dim);
93  }
94  }
95 
96  dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
97 
98  cub_weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
99 
100  node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
101 
102  jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells, num_ip, num_space_dim,num_space_dim);
103 
104  jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
105 
106  jac_det = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells, num_ip);
107 
108  weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells, num_ip);
109 
110  covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells, num_ip, num_space_dim,num_space_dim);
111 
112  contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
113 
114  norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells, num_ip);
115 
116  ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordiantes",num_cells, num_ip,num_space_dim);
117 
118  ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells, num_ip,num_space_dim);
119 
120  weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted normal",num_cells, num_ip,num_space_dim);
121 
122  }
123 
124  template <typename Scalar>
126  setupArrays(const Teuchos::RCP<const panzer::IntegrationRule>& ir)
127  {
128  MDFieldArrayFactory af(prefix,alloc_arrays);
129 
130  int_rule = ir;
131 
132  int num_nodes = ir->topology->getNodeCount();
133  int num_cells = ir->workset_size;
134  int num_space_dim = ir->topology->getDimension();
135 
136  // specialize content if this is quadrature at anode
137  if(num_space_dim==1 && ir->isSide()) {
138  setupArraysForNodeRule(ir);
139  return;
140  }
141 
142  Intrepid2::DefaultCubatureFactory<double,DblArrayDynamic>
143  cubature_factory;
144 
145  if (ir->cv_type == "side")
146  intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<double,DblArrayDynamic,DblArrayDynamic>(ir->topology));
147 
148  else if (ir->cv_type == "volume")
149  intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolume<double,DblArrayDynamic,DblArrayDynamic>(ir->topology));
150 
151  else if (ir->cv_type == "boundary" && ir->isSide())
152  intrepid_cubature = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary<double,DblArrayDynamic,DblArrayDynamic>(ir->topology,ir->side));
153 
154  else if (ir->cv_type == "none" && ir->isSide())
155  intrepid_cubature = cubature_factory.create(*(ir->side_topology),
156  ir->cubature_degree);
157  else
158  intrepid_cubature = cubature_factory.create(*(ir->topology),
159  ir->cubature_degree);
160 
161 
162  int num_ip = intrepid_cubature->getNumPoints();
163 
164  dyn_cub_points = af.template buildArray<double,IP,Dim>("cub_points",num_ip, num_space_dim);
165  dyn_cub_weights = af.template buildArray<double,IP>("cub_weights",num_ip);
166 
167  cub_points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
168 
169  if (ir->isSide() && ir->cv_type == "none") {
170  dyn_side_cub_points = af.template buildArray<double,IP,Dim>("side_cub_points",num_ip, ir->side_topology->getDimension());
171  side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip,ir->side_topology->getDimension());
172  }
173 
174  if (ir->cv_type != "none") {
175  dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>("phys_cub_points",num_cells, num_ip, num_space_dim);
176  dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>("phys_cub_weights",num_cells, num_ip);
177  if (ir->cv_type == "side") {
178  dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>("phys_cub_norms",num_cells, num_ip, num_space_dim);
179  }
180  }
181 
182  dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
183 
184  cub_weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
185 
186  node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
187 
188  jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells, num_ip, num_space_dim,num_space_dim);
189 
190  jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
191 
192  jac_det = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells, num_ip);
193 
194  weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells, num_ip);
195 
196  covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells, num_ip, num_space_dim,num_space_dim);
197 
198  contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
199 
200  norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells, num_ip);
201 
202  ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordiantes",num_cells, num_ip,num_space_dim);
203 
204  ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells, num_ip,num_space_dim);
205 
206  weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted_normal",num_cells,num_ip,num_space_dim);
207 
208  }
209 
210 // ***********************************************************
211 // * Evaluation of values - NOT specialized
212 // ***********************************************************
213  template <typename Scalar>
215  evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim> & in_node_coordinates)
216  {
217  if (int_rule->cv_type != "none") {
218  getCubatureCV(in_node_coordinates);
219  evaluateValuesCV(in_node_coordinates);
220  }
221  else {
222  getCubature(in_node_coordinates);
223  evaluateRemainingValues(in_node_coordinates);
224  }
225  }
226 
227  template <typename Scalar>
229  getCubature(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates)
230  {
231  int num_space_dim = int_rule->topology->getDimension();
232  if (int_rule->isSide() && num_space_dim==1) {
233  std::cout << "WARNING: 0-D quadrature rule ifrastructure does not exist!!! Will not be able to do "
234  << "non-natural integration rules.";
235  return;
236  }
237 
238  Intrepid2::CellTools<Scalar> cell_tools;
239 
240  if (!int_rule->isSide())
241  intrepid_cubature->getCubature(dyn_cub_points, dyn_cub_weights);
242  else {
243  intrepid_cubature->getCubature(dyn_side_cub_points, dyn_cub_weights);
244 
245  cell_tools.mapToReferenceSubcell(dyn_cub_points,
246  dyn_side_cub_points,
247  int_rule->spatial_dimension-1,
248  int_rule->side,
249  *(int_rule->topology));
250  }
251 
252  // IP coordinates
253  cell_tools.mapToPhysicalFrame(ip_coordinates, dyn_cub_points, in_node_coordinates, *(int_rule->topology));
254  }
255 
256 
257  template <typename Scalar>
259  evaluateRemainingValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates)
260  {
261  Intrepid2::CellTools<Scalar> cell_tools;
262 
263  // copy the dynamic data structures into the static data structures
264  {
265  size_type num_ip = dyn_cub_points.dimension(0);
266  size_type num_dims = dyn_cub_points.dimension(1);
267 
268  for (size_type ip = 0; ip < num_ip; ++ip) {
269  cub_weights(ip) = dyn_cub_weights(ip);
270  for (size_type dim = 0; dim < num_dims; ++dim)
271  cub_points(ip,dim) = dyn_cub_points(ip,dim);
272  }
273  }
274 
275  if (int_rule->isSide()) {
276  const size_type num_ip = dyn_cub_points.dimension(0), num_side_dims = dyn_side_cub_points.dimension(1);
277  for (size_type ip = 0; ip < num_ip; ++ip)
278  for (size_type dim = 0; dim < num_side_dims; ++dim)
279  side_cub_points(ip,dim) = dyn_side_cub_points(ip,dim);
280  }
281 
282  {
283  size_type num_cells = in_node_coordinates.dimension(0);
284  size_type num_nodes = in_node_coordinates.dimension(1);
285  size_type num_dims = in_node_coordinates.dimension(2);
286 
287  for (size_type cell = 0; cell < num_cells; ++cell) {
288  for (size_type node = 0; node < num_nodes; ++node) {
289  for (size_type dim = 0; dim < num_dims; ++dim) {
290  node_coordinates(cell,node,dim) =
291  in_node_coordinates(cell,node,dim);
292  }
293  }
294  }
295  }
296 
297  cell_tools.setJacobian(jac, cub_points, node_coordinates,
298  *(int_rule->topology));
299 
300  cell_tools.setJacobianInv(jac_inv, jac);
301 
302  cell_tools.setJacobianDet(jac_det, jac);
303 
304  if (!int_rule->isSide()) {
305  Intrepid2::FunctionSpaceTools::
306  computeCellMeasure<Scalar>(weighted_measure, jac_det, cub_weights);
307  }
308  else if(int_rule->spatial_dimension==3) {
309  Intrepid2::FunctionSpaceTools::
310  computeFaceMeasure<Scalar>(weighted_measure, jac, cub_weights,int_rule->side,*int_rule->topology);
311  }
312  else if(int_rule->spatial_dimension==2) {
313  Intrepid2::FunctionSpaceTools::
314  computeEdgeMeasure<Scalar>(weighted_measure, jac, cub_weights,int_rule->side,*int_rule->topology);
315  }
316  else TEUCHOS_ASSERT(false);
317 
318  // Shakib contravarient metric tensor
319  for (size_type cell = 0; cell < contravarient.dimension(0); ++cell) {
320  for (size_type ip = 0; ip < contravarient.dimension(1); ++ip) {
321 
322  // zero out matrix
323  for (size_type i = 0; i < contravarient.dimension(2); ++i)
324  for (size_type j = 0; j < contravarient.dimension(3); ++j)
325  covarient(cell,ip,i,j) = 0.0;
326 
327  // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha}
328  for (size_type i = 0; i < contravarient.dimension(2); ++i) {
329  for (size_type j = 0; j < contravarient.dimension(3); ++j) {
330  for (size_type alpha = 0; alpha < contravarient.dimension(2); ++alpha) {
331  covarient(cell,ip,i,j) += jac(cell,ip,i,alpha) * jac(cell,ip,j,alpha);
332  }
333  }
334  }
335 
336  }
337  }
338 
339  Intrepid2::RealSpaceTools<Scalar>::inverse(contravarient, covarient);
340 
341  // norm of g_ij
342  for (size_type cell = 0; cell < contravarient.dimension(0); ++cell) {
343  for (size_type ip = 0; ip < contravarient.dimension(1); ++ip) {
344  norm_contravarient(cell,ip) = 0.0;
345  for (size_type i = 0; i < contravarient.dimension(2); ++i) {
346  for (size_type j = 0; j < contravarient.dimension(3); ++j) {
347  norm_contravarient(cell,ip) += contravarient(cell,ip,i,j) * contravarient(cell,ip,i,j);
348  }
349  }
350  norm_contravarient(cell,ip) = std::sqrt(norm_contravarient(cell,ip));
351  }
352  }
353  }
354 
355  // Find the permutation that maps the set of points coords to other_coords. To
356  // avoid possible finite precision issues, == is not used, but rather
357  // min(norm(.)).
358  template <typename Scalar>
359  static void
360  permuteToOther(const PHX::MDField<Scalar,Cell,IP,Dim>& coords,
361  const PHX::MDField<Scalar,Cell,IP,Dim>& other_coords,
362  std::vector<typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type>& permutation)
363  {
364  typedef typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type size_type;
365  // We can safely assume: (1) The permutation is the same for every cell in
366  // the workset. (2) The first workset has valid data. Hence we operate only
367  // on cell 0.
368  const size_type cell = 0;
369  const size_type num_ip = coords.dimension(1), num_dim = coords.dimension(2);
370  permutation.resize(num_ip);
371  std::vector<char> taken(num_ip, 0);
372  for (size_type ip = 0; ip < num_ip; ++ip) {
373  // Find an other point to associate with ip.
374  size_type i_min = 0;
375  Scalar d_min = -1;
376  for (size_type other_ip = 0; other_ip < num_ip; ++other_ip) {
377  // For speed, skip other points that are already associated.
378  if (taken[other_ip]) continue;
379  // Compute the distance between the two points.
380  Scalar d(0);
381  for (size_type dim = 0; dim < num_dim; ++dim) {
382  const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim);
383  d += diff*diff;
384  }
385  if (d_min < 0 || d < d_min) {
386  d_min = d;
387  i_min = other_ip;
388  }
389  }
390  // Record the match.
391  permutation[ip] = i_min;
392  // This point is now taken.
393  taken[i_min] = 1;
394  }
395  }
396 
397  template <typename Scalar>
399  evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
400  const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates)
401  {
402  if (int_rule->cv_type == "none") {
403 
404  getCubature(in_node_coordinates);
405 
406  {
407  // Determine the permutation.
408  std::vector<size_type> permutation(other_ip_coordinates.dimension(1));
409  permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
410  // Apply the permutation to the cubature arrays.
411  MDFieldArrayFactory af(prefix, alloc_arrays);
412  const size_type num_ip = dyn_cub_points.dimension(0);
413  {
414  const size_type num_dim = dyn_side_cub_points.dimension(1);
415  DblArrayDynamic old_dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
416  "old_dyn_side_cub_points", num_ip, num_dim);
417  old_dyn_side_cub_points.deep_copy(dyn_side_cub_points);
418  for (size_type ip = 0; ip < num_ip; ++ip)
419  if (ip != permutation[ip])
420  for (size_type dim = 0; dim < num_dim; ++dim)
421  dyn_side_cub_points(ip, dim) = old_dyn_side_cub_points(permutation[ip], dim);
422  }
423  {
424  const size_type num_dim = dyn_cub_points.dimension(1);
425  DblArrayDynamic old_dyn_cub_points = af.template buildArray<double,IP,Dim>(
426  "old_dyn_cub_points", num_ip, num_dim);
427  old_dyn_cub_points.deep_copy(dyn_cub_points);
428  for (size_type ip = 0; ip < num_ip; ++ip)
429  if (ip != permutation[ip])
430  for (size_type dim = 0; dim < num_dim; ++dim)
431  dyn_cub_points(ip, dim) = old_dyn_cub_points(permutation[ip], dim);
432  }
433  {
434  DblArrayDynamic old_dyn_cub_weights = af.template buildArray<double,IP>(
435  "old_dyn_cub_weights", num_ip);
436  old_dyn_cub_weights.deep_copy(dyn_cub_weights);
437  for (size_type ip = 0; ip < dyn_cub_weights.dimension(0); ++ip)
438  if (ip != permutation[ip])
439  dyn_cub_weights(ip) = old_dyn_cub_weights(permutation[ip]);
440  }
441  {
442  const size_type num_cells = ip_coordinates.dimension(0), num_ip = ip_coordinates.dimension(1),
443  num_dim = ip_coordinates.dimension(2);
444  Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
445  "old_ip_coordinates", num_cells, num_ip, num_dim);
446  Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
447  for (size_type cell = 0; cell < num_cells; ++cell)
448  for (size_type ip = 0; ip < num_ip; ++ip)
449  if (ip != permutation[ip])
450  for (size_type dim = 0; dim < num_dim; ++dim)
451  ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
452  }
453  // All subsequent calculations inherit the permutation.
454  }
455 
456  evaluateRemainingValues(in_node_coordinates);
457  }
458 
459  else {
460 
461  getCubatureCV(in_node_coordinates);
462 
463  // Determine the permutation.
464  std::vector<size_type> permutation(other_ip_coordinates.dimension(1));
465  permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
466 
467  // Apply the permutation to the cubature arrays.
468  MDFieldArrayFactory af(prefix, alloc_arrays);
469  const size_type num_ip = dyn_cub_points.dimension(0);
470  {
471  const size_type num_cells = ip_coordinates.dimension(0), num_ip = ip_coordinates.dimension(1),
472  num_dim = ip_coordinates.dimension(2);
473  Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
474  "old_ip_coordinates", num_cells, num_ip, num_dim);
475  Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
476  Array_CellIPDim old_weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
477  "old_weighted_normals", num_cells, num_ip, num_dim);
478  Array_CellIP old_weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
479  "old_weighted_measure", num_cells, num_ip);
480  if (int_rule->cv_type == "side")
481  Kokkos::deep_copy(old_weighted_normals.get_static_view(), weighted_normals.get_static_view());
482  else
483  Kokkos::deep_copy(old_weighted_measure.get_static_view(), weighted_measure.get_static_view());
484  for (size_type cell = 0; cell < num_cells; ++cell)
485  {
486  for (size_type ip = 0; ip < num_ip; ++ip)
487  {
488  if (ip != permutation[ip]) {
489  if (int_rule->cv_type == "boundary" || int_rule->cv_type == "volume")
490  weighted_measure(cell, ip) = old_weighted_measure(cell, permutation[ip]);
491  for (size_type dim = 0; dim < num_dim; ++dim)
492  {
493  ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
494  if (int_rule->cv_type == "side")
495  weighted_normals(cell, ip, dim) = old_weighted_normals(cell, permutation[ip], dim);
496 
497  }
498  }
499  }
500  }
501  }
502 
503  evaluateValuesCV(in_node_coordinates);
504  }
505  }
506 
507  template <typename Scalar>
509  getCubatureCV(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates)
510  {
511  int num_space_dim = int_rule->topology->getDimension();
512  if (int_rule->isSide() && num_space_dim==1) {
513  std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
514  << "non-natural integration rules.";
515  return;
516  }
517  {
518  size_type num_cells = in_node_coordinates.dimension(0);
519  size_type num_nodes = in_node_coordinates.dimension(1);
520  size_type num_dims = in_node_coordinates.dimension(2);
521 
522  for (size_type cell = 0; cell < num_cells; ++cell) {
523  for (size_type node = 0; node < num_nodes; ++node) {
524  for (size_type dim = 0; dim < num_dims; ++dim) {
525  node_coordinates(cell,node,dim) =
526  in_node_coordinates(cell,node,dim);
527  dyn_node_coordinates(cell,node,dim) =
528  Sacado::ScalarValue<Scalar>::eval(in_node_coordinates(cell,node,dim));
529  }
530  }
531  }
532  }
533 
534  if (int_rule->cv_type == "side")
535  intrepid_cubature->getCubature(dyn_phys_cub_points,dyn_phys_cub_norms,dyn_node_coordinates);
536  else
537  intrepid_cubature->getCubature(dyn_phys_cub_points,dyn_phys_cub_weights,dyn_node_coordinates);
538 
539  size_type num_cells = dyn_phys_cub_points.dimension(0);
540  size_type num_ip =dyn_phys_cub_points.dimension(1);
541  size_type num_dims = dyn_phys_cub_points.dimension(2);
542 
543  for (size_type cell = 0; cell < num_cells; ++cell) {
544  for (size_type ip = 0; ip < num_ip; ++ip) {
545  if (int_rule->cv_type != "side")
546  weighted_measure(cell,ip) = dyn_phys_cub_weights(cell,ip);
547  for (size_type dim = 0; dim < num_dims; ++dim) {
548  ip_coordinates(cell,ip,dim) = dyn_phys_cub_points(cell,ip,dim);
549  if (int_rule->cv_type == "side")
550  weighted_normals(cell,ip,dim) = dyn_phys_cub_norms(cell,ip,dim);
551  }
552  }
553  }
554 
555  }
556 
557  template <typename Scalar>
559  evaluateValuesCV(const PHX::MDField<Scalar,Cell,NODE,Dim> & in_node_coordinates)
560  {
561 
562  Intrepid2::CellTools<Scalar> cell_tools;
563 
564  cell_tools.mapToReferenceFrame(ref_ip_coordinates, ip_coordinates, node_coordinates,
565  *(int_rule->topology),-1);
566 
567  cell_tools.setJacobian(jac, ref_ip_coordinates, node_coordinates,
568  *(int_rule->topology));
569 
570  cell_tools.setJacobianInv(jac_inv, jac);
571 
572  cell_tools.setJacobianDet(jac_det, jac);
573 
574  }
575 
576 #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \
577 template class IntegrationValues2<SCALAR>;
578 
581 
582 }
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates)
Cell vertex coordinates, not basis coordinates.
PHX::MDField< Scalar, Cell, IP > Array_CellIP
void evaluateRemainingValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates)
PHX::MDField< ScalarT > vector
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
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
PHX::MDField< double > DblArrayDynamic
void getCubature(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates)
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates)
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)