46#include "Shards_CellTopology.hpp"
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"
63#include "Phalanx_GetNonConstDynRankViewFromConstMDField.hpp"
69Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double>>
73 Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double> > ic;
75 Intrepid2::DefaultCubatureFactory cubature_factory;
77 if(ir.
getType() == ID::CV_SIDE){
78 ic = Teuchos::rcp(
new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir.
topology));
79 }
else if(ir.
getType() == ID::CV_VOLUME){
80 ic = Teuchos::rcp(
new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir.
topology));
81 }
else if(ir.
getType() == ID::CV_BOUNDARY){
82 TEUCHOS_ASSERT(ir.
isSide());
83 ic = Teuchos::rcp(
new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir.
topology,ir.
getSide()));
84 }
else if(ir.
getType() == ID::VOLUME){
85 ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(*(ir.
topology),ir.
getOrder());
86 }
else if(ir.
getType() == ID::SIDE){
87 ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(*(ir.
side_topology),ir.
getOrder());
88 }
else if(ir.
getType() == ID::SURFACE){
91 TEUCHOS_ASSERT(
false);
97template<
typename Scalar>
99correctVirtualNormals(PHX::MDField<Scalar,Cell,IP,Dim> normals,
100 const int num_real_cells,
101 const int num_virtual_cells,
102 const shards::CellTopology & cell_topology,
103 const SubcellConnectivity & face_connectivity)
108 const int space_dim = cell_topology.getDimension();
109 const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1);
110 const int points = normals.extent_int(1);
111 const int points_per_face = points / faces_per_cell;
113 Kokkos::parallel_for(num_virtual_cells, KOKKOS_LAMBDA(
const int & virtual_cell_ordinal){
115 const int virtual_cell = virtual_cell_ordinal+num_real_cells;
119 int face, virtual_lidx;
120 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
122 face = face_connectivity.subcellForCell(virtual_cell, local_face_id);
124 virtual_lidx = local_face_id;
130 const int real_side = (face_connectivity.cellForSubcell(face, 0) == virtual_cell) ? 1 : 0;
131 const int real_cell = face_connectivity.cellForSubcell(face,real_side);
132 const int real_lidx = face_connectivity.localSubcellForSubcell(face,real_side);
135 KOKKOS_ASSERT(real_cell < num_real_cells);
137 for(
int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
140 const int virtual_cell_point = points_per_face * virtual_lidx + point_ordinal;
141 const int real_cell_point = points_per_face * real_lidx + point_ordinal;
143 for (
int d=0; d<space_dim; d++)
144 normals(virtual_cell,virtual_cell_point,d) = -normals(real_cell,real_cell_point,d);
149 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
151 if (local_face_id == virtual_lidx)
continue;
153 for (
int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
154 const int point = local_face_id * points_per_face + point_ordinal;
155 for (
int dim=0; dim<space_dim; dim++)
156 normals(virtual_cell,point,dim) = 0.0;
160 PHX::Device::execution_space().fence();
164template<
typename Scalar>
166correctVirtualRotationMatrices(PHX::MDField<Scalar,Cell,IP,Dim,Dim> rotation_matrices,
167 const int num_real_cells,
168 const int num_virtual_cells,
169 const shards::CellTopology & cell_topology,
170 const SubcellConnectivity & face_connectivity)
175 const int space_dim = cell_topology.getDimension();
176 const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1);
177 const int points = rotation_matrices.extent_int(1);
178 const int points_per_face = points / faces_per_cell;
180 Kokkos::parallel_for(num_virtual_cells, KOKKOS_LAMBDA(
const int & virtual_cell_ordinal){
182 const int virtual_cell = virtual_cell_ordinal+num_real_cells;
186 int face, virtual_lidx;
187 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
189 face = face_connectivity.subcellForCell(virtual_cell, local_face_id);
191 virtual_lidx = local_face_id;
199 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
201 if (local_face_id == virtual_lidx)
continue;
203 for (
int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
204 const int point = local_face_id * points_per_face + point_ordinal;
205 for (
int dim=0; dim<3; dim++)
206 for (
int dim2=0; dim2<3; dim2++)
207 rotation_matrices(virtual_cell,point,dim,dim2) = 0.0;
211 PHX::Device::execution_space().fence();
214template<
typename Scalar>
216applyBasePermutation(PHX::MDField<Scalar,IP>
field,
217 PHX::MDField<const int,Cell,IP> permutations)
219 MDFieldArrayFactory af(
"",
true);
221 const int num_ip =
field.extent(0);
223 auto scratch = af.template buildStaticArray<Scalar,IP>(
"scratch", num_ip);
224 scratch.deep_copy(
field);
225 Kokkos::parallel_for(1, KOKKOS_LAMBDA(
const int & ){
226 for(
int ip=0; ip<num_ip; ++ip)
227 if (ip != permutations(0,ip))
228 field(ip) = scratch(permutations(0,ip));
230 PHX::Device::execution_space().fence();
233template<
typename Scalar>
235applyBasePermutation(PHX::MDField<Scalar,IP,Dim>
field,
236 PHX::MDField<const int,Cell,IP> permutations)
238 MDFieldArrayFactory af(
"",
true);
240 const int num_ip =
field.extent(0);
241 const int num_dim =
field.extent(1);
243 auto scratch = af.template buildStaticArray<Scalar,IP,Dim>(
"scratch", num_ip,num_dim);
244 scratch.deep_copy(
field);
245 Kokkos::parallel_for(1, KOKKOS_LAMBDA(
const int & ){
246 for(
int ip=0; ip<num_ip; ++ip)
247 if (ip != permutations(0,ip))
248 for(
int dim=0; dim<num_dim; ++dim)
249 field(ip,dim) = scratch(permutations(0,ip),dim);
251 PHX::Device::execution_space().fence();
254template<
typename Scalar>
256applyPermutation(PHX::MDField<Scalar,Cell,IP>
field,
257 PHX::MDField<const int,Cell,IP> permutations)
259 MDFieldArrayFactory af(
"",
true);
261 const int num_cells =
field.extent(0);
262 const int num_ip =
field.extent(1);
264 auto scratch = af.template buildStaticArray<Scalar,Cell,IP>(
"scratch", num_cells, num_ip);
265 scratch.deep_copy(
field);
266 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
267 for(
int ip=0; ip<num_ip; ++ip)
268 if (ip != permutations(cell,ip))
269 field(cell,ip) = scratch(cell,permutations(cell,ip));
271 PHX::Device::execution_space().fence();
274template<
typename Scalar>
276applyPermutation(PHX::MDField<Scalar,Cell,IP,Dim>
field,
277 PHX::MDField<const int,Cell,IP> permutations)
279 MDFieldArrayFactory af(
"",
true);
281 const int num_cells =
field.extent(0);
282 const int num_ip =
field.extent(1);
283 const int num_dim =
field.extent(2);
285 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"scratch", num_cells, num_ip, num_dim);
286 scratch.deep_copy(
field);
287 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
288 for(
int ip=0; ip<num_ip; ++ip)
289 if (ip != permutations(cell,ip))
290 for(
int dim=0; dim<num_dim; ++dim)
291 field(cell,ip,dim) = scratch(cell,permutations(cell,ip),dim);
293 PHX::Device::execution_space().fence();
296template<
typename Scalar>
298applyPermutation(PHX::MDField<Scalar,Cell,IP,Dim,Dim>
field,
299 PHX::MDField<const int,Cell,IP> permutations)
301 MDFieldArrayFactory af(
"",
true);
303 const int num_cells =
field.extent(0);
304 const int num_ip =
field.extent(1);
305 const int num_dim =
field.extent(2);
306 const int num_dim2 =
field.extent(3);
308 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"scratch", num_cells, num_ip, num_dim, num_dim2);
309 scratch.deep_copy(
field);
310 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
311 for(
int ip=0; ip<num_ip; ++ip)
312 if (ip != permutations(cell,ip))
313 for(
int dim=0; dim<num_dim; ++dim)
314 for(
int dim2=0; dim2<num_dim2; ++dim2)
315 field(cell,ip,dim,dim2) = scratch(cell,permutations(cell,ip),dim,dim2);
317 PHX::Device::execution_space().fence();
324template <
typename Scalar>
325PHX::MDField<const int,Cell,IP>
326generatePermutations(
const int num_cells,
327 PHX::MDField<const Scalar,Cell,IP,Dim> coords,
328 PHX::MDField<const Scalar,Cell,IP,Dim> other_coords)
330 const int num_ip = coords.extent(1);
331 const int num_dim = coords.extent(2);
333 MDFieldArrayFactory af(
"",
true);
336 auto permutation = af.template buildStaticArray<int,Cell,IP>(
"permutation", num_cells, num_ip);
337 permutation.deep_copy(0);
340 auto taken = af.template buildStaticArray<int,Cell,IP>(
"taken", num_cells, num_ip);
343 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
345 for (
int ip = 0; ip < num_ip; ++ip) {
349 for (
int other_ip = 0; other_ip < num_ip; ++other_ip) {
351 if(taken(cell,other_ip))
continue;
354 for (
int dim = 0; dim < num_dim; ++dim) {
355 const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim);
358 if (d_min < 0 || d < d_min) {
364 permutation(cell,ip) = i_min;
366 taken(cell,i_min) = 1;
369 PHX::Device::execution_space().fence();
375template <
typename Scalar>
376PHX::MDField<const int,Cell,IP>
377generateSurfacePermutations(
const int num_cells,
378 const SubcellConnectivity face_connectivity,
379 PHX::MDField<const Scalar,Cell,IP,Dim> surface_points,
380 PHX::MDField<const Scalar,Cell,IP,Dim,Dim> surface_rotation_matrices)
389 const int num_points = surface_points.extent_int(1);
390 const int num_faces_per_cell = face_connectivity.numSubcellsOnCellHost(0);
391 const int num_points_per_face = num_points / num_faces_per_cell;
392 const int cell_dim = surface_points.extent(2);
394 MDFieldArrayFactory af(
"",
true);
397 auto permutation = af.template buildStaticArray<int,Cell,IP>(
"permutation", num_cells, num_points);
398 permutation.deep_copy(0);
401 Kokkos::parallel_for(num_cells,KOKKOS_LAMBDA (
const int & cell) {
402 for(
int point=0; point<num_points; ++point)
403 permutation(cell,point) = point;
408 if(num_points_per_face != 1) {
414#define PANZER_DOT(a,b) (a[0]*b[0] + a[1]*b[1] + a[2]*b[2])
415#define PANZER_CROSS(a,b,c) {c[0] = a[1]*b[2] - a[2]*b[1]; c[1] = a[2]*b[0] - a[0]*b[2]; c[2] = a[0]*b[1] - a[1]*b[0];}
418 Kokkos::parallel_for(
"face iteration",face_connectivity.numSubcells(),KOKKOS_LAMBDA (
const int face) {
420 const int cell_0 = face_connectivity.cellForSubcell(face,0);
421 const int cell_1 = face_connectivity.cellForSubcell(face,1);
428 const int lidx_0 = face_connectivity.localSubcellForSubcell(face,0);
429 const int lidx_1 = face_connectivity.localSubcellForSubcell(face,1);
432 KOKKOS_ASSERT(lidx_1 >= 0);
437 Scalar xc0[3] = {0}, xc1[3] = {0}, r2 = 0;
438 for(
int face_point=0; face_point<num_points_per_face; ++face_point){
440 for(
int dim=0; dim<cell_dim; ++dim){
441 xc0[dim] += surface_points(cell_0,lidx_0*num_points_per_face + face_point,dim);
442 xc1[dim] += surface_points(cell_1,lidx_1*num_points_per_face + face_point,dim);
443 const Scalar dx = surface_points(cell_0,lidx_0*num_points_per_face + face_point,dim) - surface_points(cell_0,lidx_0*num_points_per_face,dim);
448 r2 = (r2 < dx2) ? dx2 : r2;
450 for(
int dim=0; dim<cell_dim; ++dim){
451 xc0[dim] /= double(num_points_per_face);
452 xc1[dim] /= double(num_points_per_face);
461 const int example_point_0 = lidx_0*num_points_per_face;
462 const int example_point_1 = lidx_1*num_points_per_face;
465 Scalar t[3] = {surface_rotation_matrices(cell_0,example_point_0,1,0), surface_rotation_matrices(cell_0,example_point_0,1,1), surface_rotation_matrices(cell_0,example_point_0,1,2)};
466 Scalar b[3] = {surface_rotation_matrices(cell_0,example_point_0,2,0), surface_rotation_matrices(cell_0,example_point_0,2,1), surface_rotation_matrices(cell_0,example_point_0,2,2)};
472 const Scalar n0[3] = {surface_rotation_matrices(cell_0,example_point_0,0,0), surface_rotation_matrices(cell_0,example_point_0,0,1), surface_rotation_matrices(cell_0,example_point_0,0,2)};
473 const Scalar n1[3] = {surface_rotation_matrices(cell_1,example_point_1,0,0), surface_rotation_matrices(cell_1,example_point_1,0,1), surface_rotation_matrices(cell_1,example_point_1,0,2)};
481 if(Kokkos::fabs(n0_dot_n1 - 1.) < 1.e-8)
485 if(Kokkos::fabs(n0_dot_n1 + 1.) > 1.e-2){
492 const Scalar mag_t = Kokkos::sqrt(
PANZER_DOT(t,t));
498 b[0] = n0[0] + n1[0];
499 b[1] = n0[1] + n1[1];
500 b[2] = n0[2] + n1[2];
503 const Scalar mag_b = Kokkos::sqrt(
PANZER_DOT(b,b));
521 for(
int face_point_1=0; face_point_1<num_points_per_face; ++face_point_1){
524 const int point_1 = lidx_1*num_points_per_face + face_point_1;
527 for(
int dim=0; dim<cell_dim; ++dim)
528 x1[dim] = surface_points(cell_1,point_1,dim) - xc1[dim];
535 for(
int face_point_0=0; face_point_0<num_points_per_face; ++face_point_0){
538 const int point_0 = lidx_0*num_points_per_face + face_point_0;
541 for(
int dim=0; dim<cell_dim; ++dim)
542 x0[dim] = surface_points(cell_0,point_0,dim) - xc0[dim];
549 const Scalar p012 = (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]);
553 if(p012 / r2 < 1.e-12){
554 permutation(cell_1,lidx_1*num_points_per_face+face_point_0) = point_1;
561 PHX::Device::execution_space().fence();
576template <
typename Scalar>
579 const bool allocArrays):
581 num_evaluate_cells_(0),
582 num_virtual_cells_(-1),
583 requires_permutation_(false),
584 alloc_arrays_(allocArrays),
591template <
typename Scalar>
596 cub_points_evaluated_ =
false;
597 side_cub_points_evaluated_ =
false;
598 cub_weights_evaluated_ =
false;
599 node_coordinates_evaluated_ =
false;
600 jac_evaluated_ =
false;
601 jac_inv_evaluated_ =
false;
602 jac_det_evaluated_ =
false;
603 weighted_measure_evaluated_ =
false;
604 weighted_normals_evaluated_ =
false;
605 surface_normals_evaluated_ =
false;
606 surface_rotation_matrices_evaluated_ =
false;
607 covarient_evaluated_ =
false;
608 contravarient_evaluated_ =
false;
609 norm_contravarient_evaluated_ =
false;
610 ip_coordinates_evaluated_ =
false;
611 ref_ip_coordinates_evaluated_ =
false;
616template <
typename Scalar>
618setupArrays(
const Teuchos::RCP<const panzer::IntegrationRule>& ir)
626 const int num_nodes = ir->topology->getNodeCount();
627 const int num_cells = ir->workset_size;
628 const int num_space_dim = ir->topology->getDimension();
631 const bool is_node_rule = (num_space_dim==1) and ir->isSide();
632 if(not is_node_rule) {
633 TEUCHOS_ASSERT(ir->getType() != ID::NONE);
634 intrepid_cubature = getIntrepidCubature(*ir);
637 const int num_ip = is_node_rule ? 1 : ir->num_points;
639 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
641 if (ir->isSide() && ir->cv_type ==
"none")
642 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip,ir->side_topology->getDimension());
644 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
646 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
648 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells, num_ip, num_space_dim,num_space_dim);
650 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
652 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells, num_ip);
654 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells, num_ip);
656 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells, num_ip, num_space_dim,num_space_dim);
658 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
660 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells, num_ip);
662 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells, num_ip,num_space_dim);
664 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells, num_ip,num_space_dim);
666 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted_normal",num_cells,num_ip,num_space_dim);
668 surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells, num_ip,num_space_dim);
670 surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells, num_ip,3,3);
678template <
typename Scalar>
680evaluateValues(
const PHX::MDField<Scalar,Cell,NODE,Dim> & in_node_coordinates,
681 const int in_num_cells,
682 const Teuchos::RCP<const SubcellConnectivity> & face_connectivity,
683 const int num_virtual_cells)
686 setup(int_rule, in_node_coordinates, in_num_cells);
690 if((int_rule->getType() == ID::SURFACE) and (face_connectivity != Teuchos::null))
691 setupPermutations(face_connectivity, num_virtual_cells);
694 evaluateEverything();
698template <
typename Scalar>
701evaluateValues(
const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
702 const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates,
703 const int in_num_cells)
706 setup(int_rule, in_node_coordinates, in_num_cells);
709 setupPermutations(other_ip_coordinates);
712 evaluateEverything();
716template <
typename Scalar>
719setupPermutations(
const Teuchos::RCP<const SubcellConnectivity> & face_connectivity,
720 const int num_virtual_cells)
722 TEUCHOS_ASSERT(not int_rule->isSide());
723 TEUCHOS_ASSERT(face_connectivity != Teuchos::null);
725 "IntegrationValues2::setupPermutations : Face connectivity based permutations are only required for surface integration schemes");
726 TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ >= 0,
727 "IntegrationValues2::setupPermutations : Number of virtual cells for surface permuations must be >= 0");
729 side_connectivity_ = face_connectivity;
730 num_virtual_cells_ = num_virtual_cells;
731 requires_permutation_ =
false;
732 permutations_ = generateSurfacePermutations<Scalar>(num_evaluate_cells_,*face_connectivity, getCubaturePoints(
false,
true), getSurfaceRotationMatrices(
false,
true));
733 requires_permutation_ =
true;
736template <
typename Scalar>
742 requires_permutation_ =
false;
743 permutations_ = generatePermutations<Scalar>(num_evaluate_cells_, getCubaturePoints(
false,
true), other_ip_coordinates);
744 requires_permutation_ =
true;
748template <
typename Scalar>
751setup(
const Teuchos::RCP<const panzer::IntegrationRule>& ir,
752 const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
759 num_cells_ = vertex_coordinates.extent(0);
760 num_evaluate_cells_ = num_cells < 0 ? vertex_coordinates.extent(0) : num_cells;
764 intrepid_cubature = getIntrepidCubature(*ir);
770 const int num_space_dim = int_rule->topology->getDimension();
771 const int num_vertexes = int_rule->topology->getNodeCount();
772 TEUCHOS_ASSERT(
static_cast<int>(vertex_coordinates.extent(1)) == num_vertexes);
774 auto aux = af.template buildStaticArray<Scalar,Cell,NODE,Dim>(
"node_coordinates",num_cells_,num_vertexes, num_space_dim);
775 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_vertexes,num_space_dim});
776 Kokkos::parallel_for(policy,KOKKOS_LAMBDA(
const int & cell,
const int & point,
const int & dim){
777 aux(cell,point,dim) = vertex_coordinates(cell,point,dim);
779 PHX::Device::execution_space().fence();
780 node_coordinates = aux;
785template <
typename Scalar>
790 const bool apply_permutation)
const
793 if(cub_points_evaluated_ and not force)
796 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
799 int num_space_dim = int_rule->topology->getDimension();
800 int num_ip = int_rule->num_points;
802 auto aux = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
804 if (int_rule->isSide() && num_space_dim==1) {
805 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
806 <<
"non-natural integration rules.";
810 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type !=
"none",
811 "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme.");
814 "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for surface integration scheme.");
816 auto weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
818 if (!int_rule->isSide())
819 intrepid_cubature->getCubature(aux.get_view(), weights.get_view());
821 auto s_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip, num_space_dim-1);
823 intrepid_cubature->getCubature(s_cub_points.get_view(), weights.get_view());
824 cell_tools.mapToReferenceSubcell(aux.get_view(), s_cub_points.get_view(), num_space_dim-1, int_rule->getSide(), *(int_rule->topology));
827 PHX::Device::execution_space().fence();
829 if(apply_permutation and requires_permutation_)
830 applyBasePermutation(aux, permutations_);
834 cub_points_evaluated_ =
true;
841template <
typename Scalar>
846 const bool apply_permutation)
const
849 if(side_cub_points_evaluated_ and not force)
850 return side_cub_points;
852 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
855 int num_space_dim = int_rule->topology->getDimension();
856 int num_ip = int_rule->num_points;
858 auto aux = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip, num_space_dim-1);
860 if (int_rule->isSide() && num_space_dim==1) {
861 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
862 <<
"non-natural integration rules.";
866 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type !=
"none",
867 "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme.");
870 "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for surface integration scheme.");
872 TEUCHOS_TEST_FOR_EXCEPT_MSG(!int_rule->isSide(),
873 "IntegrationValues2::getUniformSideCubaturePointsRef : Requested side points, which is not supported by integration rule.");
875 auto weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
877 intrepid_cubature->getCubature(aux.get_view(), weights.get_view());
879 PHX::Device::execution_space().fence();
881 if(apply_permutation and requires_permutation_)
882 applyBasePermutation(aux, permutations_);
885 side_cub_points = aux;
886 side_cub_points_evaluated_ =
true;
892template <
typename Scalar>
897 const bool apply_permutation)
const
900 if(cub_weights_evaluated_ and not force)
903 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
906 int num_space_dim = int_rule->topology->getDimension();
907 int num_ip = int_rule->num_points;
909 auto aux = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
911 if (int_rule->isSide() && num_space_dim==1) {
912 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
913 <<
"non-natural integration rules.";
917 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type !=
"none",
918 "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for control volume integration scheme.");
921 "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for surface integration scheme.");
923 auto points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip,num_space_dim);
925 intrepid_cubature->getCubature(points.get_view(), aux.get_view());
927 PHX::Device::execution_space().fence();
929 if(apply_permutation and requires_permutation_)
930 applyBasePermutation(aux, permutations_);
934 cub_weights_evaluated_ =
true;
941template <
typename Scalar>
946 return node_coordinates;
949template <
typename Scalar>
953 const bool force)
const
955 if(jac_evaluated_ and not force)
958 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
961 int num_space_dim = int_rule->topology->getDimension();
962 int num_ip = int_rule->num_points;
965 auto const_ref_coord = getCubaturePointsRef(
false,force);
966 auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord);
967 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
968 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells_, num_ip, num_space_dim,num_space_dim);
970 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
971 auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
972 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
973 auto s_jac = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
975 cell_tools.setJacobian(s_jac, s_ref_coord, s_node_coord,*(int_rule->topology));
977 PHX::Device::execution_space().fence();
981 jac_evaluated_ =
true;
988template <
typename Scalar>
992 const bool force)
const
995 if(jac_inv_evaluated_ and not force)
998 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1001 const int num_space_dim = int_rule->topology->getDimension();
1002 const int num_ip = int_rule->num_points;
1004 auto jacobian = getJacobian(
false,force);
1005 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells_, num_ip, num_space_dim,num_space_dim);
1007 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1008 auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1009 auto s_jac_inv = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1011 cell_tools.setJacobianInv(s_jac_inv, s_jac);
1013 PHX::Device::execution_space().fence();
1017 jac_inv_evaluated_ =
true;
1024template <
typename Scalar>
1028 const bool force)
const
1031 if(jac_det_evaluated_ and not force)
1034 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1037 const int num_ip = int_rule->num_points;
1039 auto jacobian = getJacobian(
false,force);
1040 auto aux = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells_, num_ip);
1042 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1043 auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1044 auto s_jac_det = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL());
1046 cell_tools.setJacobianDet(s_jac_det, s_jac);
1048 PHX::Device::execution_space().fence();
1052 jac_det_evaluated_ =
true;
1059template <
typename Scalar>
1063 const bool force)
const
1066 if(weighted_measure_evaluated_ and not force)
1067 return weighted_measure;
1069 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1072 const int num_space_dim = int_rule->topology->getDimension();
1073 const int num_ip = int_rule->num_points;
1075 auto aux = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells_, num_ip);
1077 if(int_rule->cv_type !=
"none"){
1079 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type ==
"side",
1080 "IntegrationValues2::getWeightedMeasure : Weighted measure not available for side control volume integrators. Use getWeightedNormals instead.");
1084 auto s_cub_points = af.template buildStaticArray<Scalar, Cell, IP, Dim>(
"cub_points",num_evaluate_cells_,num_ip,num_space_dim);
1086 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1088 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1089 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1090 auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL());
1092 intrepid_cubature->getCubature(s_cub_points.get_view(),s_weighted_measure,s_node_coord);
1096 const auto & cell_topology = *int_rule->topology;
1097 const int cell_dim = cell_topology.getDimension();
1098 const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount();
1099 const int cubature_order = int_rule->order();
1100 const int num_points_on_side = num_ip / num_sides;
1102 Intrepid2::DefaultCubatureFactory cubature_factory;
1103 auto jacobian = getJacobian(
false,force);
1105 for(
int side=0; side<num_sides; ++side) {
1107 const int point_offset=side*num_points_on_side;
1110 Kokkos::DynRankView<double,PHX::Device> side_cub_weights;
1112 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"side_cub_weights",num_points_on_side);
1113 auto tmp_side_cub_weights_host = Kokkos::create_mirror_view(side_cub_weights);
1114 tmp_side_cub_weights_host(0)=1.;
1115 Kokkos::deep_copy(side_cub_weights,tmp_side_cub_weights_host);
1119 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(cell_dim-1,side));
1121 auto ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(face_topology,cubature_order);
1123 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"side_cub_weights",num_points_on_side);
1124 auto subcell_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"subcell_cub_points",num_points_on_side,cell_dim-1);
1127 ic->getCubature(subcell_cub_points, side_cub_weights);
1130 PHX::Device::execution_space().fence();
1133 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_side});
1136 auto side_weighted_measure = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_weighted_measure",num_evaluate_cells_,num_points_on_side);
1138 auto side_cub_weights_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),side_cub_weights);
1139 Kokkos::deep_copy(side_weighted_measure, side_cub_weights_host(0));
1143 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_jac",num_evaluate_cells_,num_points_on_side,cell_dim,cell_dim);
1145 Kokkos::parallel_for(
"Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1146 for(
int dim=0;dim<cell_dim;++dim)
1147 for(
int dim1=0;dim1<cell_dim;++dim1)
1148 side_jacobian(cell,point,dim,dim1) = jacobian(cell,point_offset+point,dim,dim1);
1151 auto scratch = af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_measure", num_evaluate_cells_*num_points_on_side*num_space_dim*num_space_dim);
1154 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1155 computeEdgeMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
1157 scratch.get_view());
1158 PHX::Device::execution_space().fence();
1159 }
else if(cell_dim == 3){
1160 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1161 computeFaceMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
1163 scratch.get_view());
1164 PHX::Device::execution_space().fence();
1170 Kokkos::parallel_for(
"copy surface weighted measure values",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1171 aux(cell,point_offset + point) = side_weighted_measure(cell,point);
1173 PHX::Device::execution_space().fence();
1178 auto cell_range = std::make_pair(0,num_evaluate_cells_);
1179 auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL());
1180 auto cubature_weights = getUniformCubatureWeightsRef(
false,force,
false);
1182 if (!int_rule->isSide()) {
1184 auto s_jac_det = Kokkos::subview(getJacobianDeterminant(
false,force).get_view(),cell_range,Kokkos::ALL());
1185 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1186 computeCellMeasure(s_weighted_measure, s_jac_det, cubature_weights.get_view());
1188 }
else if(int_rule->spatial_dimension==3) {
1190 auto s_jac = Kokkos::subview(getJacobian(
false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1191 auto scratch = af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim);
1192 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1193 computeFaceMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(),
1194 int_rule->side, *int_rule->topology,
1195 scratch.get_view());
1197 }
else if(int_rule->spatial_dimension==2) {
1199 auto s_jac = Kokkos::subview(getJacobian(
false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1200 auto scratch = af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim);
1201 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1202 computeEdgeMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(),
1203 int_rule->side,*int_rule->topology,
1204 scratch.get_view());
1207 TEUCHOS_ASSERT(
false);
1212 PHX::Device::execution_space().fence();
1215 if(requires_permutation_)
1216 applyPermutation(aux, permutations_);
1219 weighted_measure = aux;
1220 weighted_measure_evaluated_ =
true;
1227template <
typename Scalar>
1231 const bool force)
const
1234 if(weighted_normals_evaluated_ and not force)
1235 return weighted_normals;
1237 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1240 const int num_space_dim = int_rule->topology->getDimension();
1241 const int num_ip = int_rule->num_points;
1243 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted_normals",num_cells_,num_ip,num_space_dim);
1246 "IntegrationValues2::getWeightedNormals : Cannot build reference weighted normals for surface integration scheme.");
1248 auto points = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"cub_points",num_cells_,num_ip,num_space_dim);
1250 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1252 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1253 auto s_cub_points = Kokkos::subview(points.get_view(),cell_range, Kokkos::ALL(), Kokkos::ALL());
1254 auto s_weighted_normals = Kokkos::subview(aux.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1255 auto s_node_coord = Kokkos::subview(node_coord, cell_range, Kokkos::ALL(), Kokkos::ALL());
1257 intrepid_cubature->getCubature(s_cub_points,s_weighted_normals,s_node_coord);
1259 PHX::Device::execution_space().fence();
1262 if(requires_permutation_)
1263 applyPermutation(aux, permutations_);
1266 weighted_normals = aux;
1267 weighted_normals_evaluated_ =
true;
1274template <
typename Scalar>
1278 const bool force)
const
1281 if(surface_normals_evaluated_ and not force)
1282 return surface_normals;
1284 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->isSide(),
1285 "IntegrationValues2::getSurfaceNormals : This call doesn't work with sides (only surfaces).");
1287 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type !=
"none",
1288 "IntegrationValues2::getSurfaceNormals : This call does not support control volume integration schemes.");
1291 "IntegrationValues2::getSurfaceNormals : Can only build for surface integrators.");
1293 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1296 const shards::CellTopology & cell_topology = *(int_rule->topology);
1297 const int cell_dim = cell_topology.getDimension();
1298 const int subcell_dim = cell_topology.getDimension()-1;
1299 const int num_subcells = cell_topology.getSubcellCount(subcell_dim);
1300 const int num_space_dim = int_rule->topology->getDimension();
1301 const int num_ip = int_rule->num_points;
1302 const int num_points_on_face = num_ip / num_subcells;
1304 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells_,num_ip,num_space_dim);
1309 jacobian = getJacobian(
false,force);
1312 for(
int subcell_index=0; subcell_index<num_subcells; ++subcell_index) {
1314 const int point_offset = subcell_index * num_points_on_face;;
1317 auto side_normals = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_normals",num_evaluate_cells_,num_points_on_face,cell_dim);
1320 const int other_subcell_index = (subcell_index==0) ? 1 : 0;
1321 auto in_node_coordinates_k = getNodeCoordinates().get_view();
1322 const auto min = std::numeric_limits<typename Sacado::ScalarType<Scalar>::type>::min();
1324 Kokkos::parallel_for(
"compute normals 1D",num_evaluate_cells_,KOKKOS_LAMBDA (
const int cell) {
1325 Scalar norm = (in_node_coordinates_k(cell,subcell_index,0) - in_node_coordinates_k(cell,other_subcell_index,0));
1326 side_normals(cell,0,0) = norm / fabs(norm + min);
1332 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_face});
1334 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_jac",num_evaluate_cells_,num_points_on_face,cell_dim,cell_dim);
1335 Kokkos::parallel_for(
"Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1336 for(
int dim=0;dim<cell_dim;++dim)
1337 for(
int dim1=0;dim1<cell_dim;++dim1)
1338 side_jacobian(cell,point,dim,dim1) = jacobian(cell,point_offset+point,dim,dim1);
1342 cell_tools.getPhysicalSideNormals(side_normals,side_jacobian,subcell_index,cell_topology);
1345 Kokkos::parallel_for(
"Normalize the normals",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1347 for(
int dim=0;dim<cell_dim;++dim){
1348 n += side_normals(cell,point,dim)*side_normals(cell,point,dim);
1352 n = Kokkos::sqrt(n);
1353 for(
int dim=0;dim<cell_dim;++dim)
1354 side_normals(cell,point,dim) /= n;
1359 PHX::Device::execution_space().fence();
1361 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_face});
1362 Kokkos::parallel_for(
"copy surface normals", policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1363 for(
int dim=0;dim<cell_dim;++dim)
1364 aux(cell,point_offset + point,dim) = side_normals(cell,point,dim);
1366 PHX::Device::execution_space().fence();
1371 TEUCHOS_TEST_FOR_EXCEPT_MSG(side_connectivity_ == Teuchos::null,
1372 "IntegrationValues2::getSurfaceNormals : Surface normals require a SubcellConnectivity object pass in through the setupPermutations call");
1373 TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ < 0,
1374 "IntegrationValues2::getSurfaceNormals : Number of virtual cells (see setupPermutations or evaluateValues) must be greater or equal to zero");
1377 if(num_virtual_cells_ > 0)
1378 correctVirtualNormals(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_);
1382 surface_normals = aux;
1383 surface_normals_evaluated_ =
true;
1390template <
typename Scalar>
1394 const bool force)
const
1397 if(surface_rotation_matrices_evaluated_ and not force)
1398 return surface_rotation_matrices;
1402 const int num_ip = int_rule->num_points;
1403 const int cell_dim = int_rule->topology->getDimension();
1406 auto normals = getSurfaceNormals(
false,force).get_static_view();
1407 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells_, num_ip, 3, 3);
1409 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1410 Kokkos::parallel_for(
"create surface rotation matrices",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1412 for(
int i=0;i<3;i++)
1414 for(
int dim=0; dim<cell_dim; ++dim)
1415 normal[dim] = normals(cell,point,dim);
1417 Scalar transverse[3];
1421 for(
int dim=0; dim<3; ++dim){
1422 aux(cell,point,0,dim) = normal[dim];
1423 aux(cell,point,1,dim) = transverse[dim];
1424 aux(cell,point,2,dim) = binormal[dim];
1427 PHX::Device::execution_space().fence();
1431 TEUCHOS_TEST_FOR_EXCEPT_MSG(side_connectivity_ == Teuchos::null,
1432 "IntegrationValues2::getSurfaceRotationMatrices : Surface normals require a SubcellConnectivity object pass in through the setupPermutations call");
1433 TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ < 0,
1434 "IntegrationValues2::getSurfaceRotationMatrices : Number of virtual cells (see setupPermutations or evaluateValues) must be greater or equal to zero");
1437 if(num_virtual_cells_ > 0)
1438 correctVirtualRotationMatrices(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_);
1442 surface_rotation_matrices = aux;
1443 surface_rotation_matrices_evaluated_ =
true;
1449template <
typename Scalar>
1453 const bool force)
const
1456 if(covarient_evaluated_ and not force)
1461 const int num_space_dim = int_rule->topology->getDimension();
1462 const int num_ip = int_rule->num_points;
1464 auto jacobian = getJacobian(
false,force).get_static_view();
1465 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells_, num_ip, num_space_dim,num_space_dim);
1467 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1468 Kokkos::parallel_for(
"evalaute covarient metric tensor",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
1470 for (
int i = 0; i < num_space_dim; ++i) {
1471 for (
int j = 0; j < num_space_dim; ++j) {
1473 for (
int alpha = 0; alpha < num_space_dim; ++alpha)
1474 value += jacobian(cell,ip,i,alpha) * jacobian(cell,ip,j,alpha);
1475 aux(cell,ip,i,j) = value;
1479 PHX::Device::execution_space().fence();
1483 covarient_evaluated_ =
true;
1490template <
typename Scalar>
1494 const bool force)
const
1497 if(contravarient_evaluated_ and not force)
1498 return contravarient;
1502 const int num_space_dim = int_rule->topology->getDimension();
1503 const int num_ip = int_rule->num_points;
1505 auto cov = getCovarientMatrix(
false,force);
1506 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells_, num_ip, num_space_dim,num_space_dim);
1508 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1509 auto s_contravarient = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
1510 auto s_covarient = Kokkos::subview(cov.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
1512 Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
1513 PHX::Device::execution_space().fence();
1516 contravarient = aux;
1517 contravarient_evaluated_ =
true;
1524template <
typename Scalar>
1528 const bool force)
const
1531 if(norm_contravarient_evaluated_ and not force)
1532 return norm_contravarient;
1536 const int num_space_dim = int_rule->topology->getDimension();
1537 const int num_ip = int_rule->num_points;
1539 auto con = getContravarientMatrix(
false,force).get_static_view();
1540 auto aux = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells_, num_ip);
1543 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1544 Kokkos::parallel_for(
"evaluate norm_contravarient",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
1546 for (
int i = 0; i < num_space_dim; ++i) {
1547 for (
int j = 0; j < num_space_dim; ++j) {
1548 aux(cell,ip) += con(cell,ip,i,j) * con(cell,ip,i,j);
1551 aux(cell,ip) = Kokkos::sqrt(aux(cell,ip));
1553 PHX::Device::execution_space().fence();
1556 norm_contravarient = aux;
1557 norm_contravarient_evaluated_ =
true;
1564template <
typename Scalar>
1568 const bool force)
const
1571 if(ip_coordinates_evaluated_ and not force)
1572 return ip_coordinates;
1576 const int num_space_dim = int_rule->topology->getDimension();
1577 const int num_ip = int_rule->num_points;
1579 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordinates",num_cells_, num_ip, num_space_dim);
1582 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1584 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1589 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1590 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1591 auto s_cub_points = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1594 if(int_rule->cv_type ==
"side"){
1595 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"scratch",num_evaluate_cells_,num_ip,num_space_dim);
1596 intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord);
1599 TEUCHOS_ASSERT((int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_BOUNDARY));
1600 auto scratch = af.template buildStaticArray<Scalar,Cell,IP>(
"scratch",num_evaluate_cells_,num_ip);
1601 intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord);
1607 auto const_ref_coord = getCubaturePointsRef(
false,force);
1608 auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord);
1610 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1611 auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1612 auto s_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1613 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1615 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1616 cell_tools.mapToPhysicalFrame(s_coord, s_ref_coord, s_node_coord, *(int_rule->topology));
1620 PHX::Device::execution_space().fence();
1623 ip_coordinates = aux;
1624 ip_coordinates_evaluated_ =
true;
1632template <
typename Scalar>
1636 const bool force)
const
1639 if(ref_ip_coordinates_evaluated_)
1640 return ref_ip_coordinates;
1643 const bool is_surface = int_rule->
getType() == ID::SURFACE;
1644 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1646 const int num_space_dim = int_rule->topology->getDimension();
1647 const int num_ip = int_rule->num_points;
1651 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1653 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells_, num_ip, num_space_dim);
1659 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1662 auto const_coord = getCubaturePoints(
false,force);
1663 auto coord = PHX::getNonConstDynRankViewFromConstMDField(const_coord);
1665 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1666 auto s_ref_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1667 auto s_coord = Kokkos::subview(coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1668 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1670 cell_tools.mapToReferenceFrame(s_ref_coord, s_coord, s_node_coord, *(int_rule->topology));
1672 }
else if(is_surface){
1674 const auto & cell_topology = *int_rule->topology;
1675 const int cell_dim = cell_topology.getDimension();
1676 const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount();
1677 const int subcell_dim = cell_dim-1;
1678 const int num_points_on_face = num_ip / num_sides;
1679 const int order = int_rule->getOrder();
1682 auto side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_points_on_face,cell_dim);
1684 Intrepid2::DefaultCubatureFactory cubature_factory;
1687 for(
int side=0; side<num_sides; ++side) {
1689 const int point_offset = side*num_points_on_face;
1694 auto side_cub_points_host = Kokkos::create_mirror_view(side_cub_points.get_view());
1695 side_cub_points_host(0,0) = (side==0)? -1. : 1.;
1696 Kokkos::deep_copy(side_cub_points.get_view(),side_cub_points_host);
1700 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,side));
1703 auto ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(face_topology,order);
1704 auto tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"tmp_side_cub_weights",num_points_on_face);
1705 auto tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"tmp_side_cub_points",num_points_on_face,subcell_dim);
1708 ic->getCubature(tmp_side_cub_points, tmp_side_cub_weights);
1711 cell_tools.mapToReferenceSubcell(side_cub_points.get_view(), tmp_side_cub_points, subcell_dim, side, cell_topology);
1714 PHX::Device::execution_space().fence();
1717 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_points_on_face, num_space_dim});
1718 Kokkos::parallel_for(
"copy values",policy,KOKKOS_LAMBDA (
const int cell,
const int point,
const int dim) {
1719 aux(cell,point_offset + point,dim) = side_cub_points(point,dim);
1721 PHX::Device::execution_space().fence();
1727 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->isSide() && num_space_dim==1,
1728 "ERROR: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
1729 <<
"non-natural integration rules.");
1731 auto cub_points = getUniformCubaturePointsRef(
false,force,
false);
1733 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_ip,num_space_dim});
1734 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(
const int & cell,
const int & ip,
const int & dim){
1735 aux(cell,ip,dim) = cub_points(ip,dim);
1739 PHX::Device::execution_space().fence();
1741 if(requires_permutation_)
1742 applyPermutation(aux, permutations_);
1745 ref_ip_coordinates = aux;
1746 ref_ip_coordinates_evaluated_ =
true;
1753template <
typename Scalar>
1760 const bool is_surface = int_rule->
getType() == ID::SURFACE;
1761 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1762 const bool is_side = int_rule->isSide();
1768 getCubaturePoints(
true,
true);
1769 getCubaturePointsRef(
true,
true);
1772 getUniformCubaturePointsRef(
true,
true);
1773 getUniformCubatureWeightsRef(
true,
true);
1775 getUniformSideCubaturePointsRef(
true,
true);
1777 getCubaturePointsRef(
true,
true);
1778 getCubaturePoints(
true,
true);
1782 getJacobian(
true,
true);
1783 getJacobianDeterminant(
true,
true);
1784 getJacobianInverse(
true,
true);
1785 if(int_rule->cv_type ==
"side")
1786 getWeightedNormals(
true,
true);
1788 getWeightedMeasure(
true,
true);
1792 getSurfaceNormals(
true,
true);
1793 getSurfaceRotationMatrices(
true,
true);
1797 if(not (is_surface or is_cv)){
1798 getContravarientMatrix(
true,
true);
1799 getCovarientMatrix(
true,
true);
1800 getNormContravarientMatrix(
true,
true);
1805#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \
1806 template class IntegrationValues2<SCALAR>;
#define PANZER_CROSS(a, b, c)
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
@ SURFACE
Integral over volume.
const int & getOrder() const
Get order of integrator.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
const int & getType() const
Get type of integrator.
ConstArray_IPDim getUniformCubaturePointsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature points.
void setupPermutations(const Teuchos::RCP< const SubcellConnectivity > &face_connectivity, const int num_virtual_cells)
Initialize the permutation arrays given a face connectivity.
ConstArray_CellIP getWeightedMeasure(const bool cache=true, const bool force=false) const
Get the weighted measure (integration weights)
ConstArray_CellIPDimDim getJacobian(const bool cache=true, const bool force=false) const
Get the Jacobian matrix evaluated at the cubature points.
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int num_cells=-1, const Teuchos::RCP< const SubcellConnectivity > &face_connectivity=Teuchos::null, const int num_virtual_cells=-1)
Evaluate basis values.
PHX::MDField< const Scalar, Cell, IP > ConstArray_CellIP
ConstArray_CellIPDim getCubaturePointsRef(const bool cache=true, const bool force=false) const
Get the cubature points in the reference space.
ConstArray_CellIPDimDim getSurfaceRotationMatrices(const bool cache=true, const bool force=false) const
Get the surface rotation matrices.
PHX::MDField< const Scalar, IP, Dim > ConstArray_IPDim
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
PHX::MDField< const Scalar, Cell, IP, Dim, Dim > ConstArray_CellIPDimDim
PHX::MDField< const Scalar, IP > ConstArray_IP
IntegrationValues2(const std::string &pre="", const bool allocArrays=false)
Base constructor.
ConstArray_CellIPDim getWeightedNormals(const bool cache=true, const bool force=false) const
Get the weighted normals.
ConstArray_CellIP getNormContravarientMatrix(const bool cache=true, const bool force=false) const
Get the contravarient matrix.
ConstArray_CellIPDimDim getContravarientMatrix(const bool cache=true, const bool force=false) const
Get the contravarient matrix.
ConstArray_CellIPDim getCubaturePoints(const bool cache=true, const bool force=false) const
Get the cubature points in physical space.
PHX::MDField< const Scalar, Cell, IP, Dim > ConstArray_CellIPDim
void evaluateEverything()
void setup(const Teuchos::RCP< const panzer::IntegrationRule > &ir, const PHX::MDField< Scalar, Cell, NODE, Dim > &node_coordinates, const int num_cells=-1)
Main setup call for the lazy evaluation interface.
ConstArray_CellIPDimDim getCovarientMatrix(const bool cache=true, const bool force=false) const
Get the covarient matrix.
ConstArray_CellBASISDim getNodeCoordinates() const
Get the node coordinates describing the geometry of the mesh.
ConstArray_IPDim getUniformSideCubaturePointsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature points for a side.
ConstArray_CellIPDim getSurfaceNormals(const bool cache=true, const bool force=false) const
Get the surface normals.
PHX::MDField< const Scalar, Cell, BASIS, Dim > ConstArray_CellBASISDim
ConstArray_IP getUniformCubatureWeightsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature weights.
ConstArray_CellIPDimDim getJacobianInverse(const bool cache=true, const bool force=false) const
Get the inverse of the Jacobian matrix evaluated at the cubature points.
ConstArray_CellIP getJacobianDeterminant(const bool cache=true, const bool force=false) const
Get the determinant of the Jacobian matrix evaluated at the cubature points.
Teuchos::RCP< shards::CellTopology > side_topology
Teuchos::RCP< const shards::CellTopology > topology
KOKKOS_INLINE_FUNCTION void convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])