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"
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){
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){
97 template<
typename Scalar>
100 const int num_real_cells,
101 const int num_virtual_cells,
102 const shards::CellTopology & cell_topology,
103 const SubcellConnectivity & face_connectivity)
105 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::correctVirtualNormals()",corr_virt_norms);
109 const int space_dim = cell_topology.getDimension();
110 const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1);
111 const int points = normals.extent_int(1);
112 const int points_per_face = points / faces_per_cell;
114 Kokkos::parallel_for(num_virtual_cells, KOKKOS_LAMBDA(
const int & virtual_cell_ordinal){
116 const int virtual_cell = virtual_cell_ordinal+num_real_cells;
120 int face, virtual_lidx;
121 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
123 face = face_connectivity.subcellForCell(virtual_cell, local_face_id);
125 virtual_lidx = local_face_id;
131 const int real_side = (face_connectivity.cellForSubcell(face, 0) == virtual_cell) ? 1 : 0;
132 const int real_cell = face_connectivity.cellForSubcell(face,real_side);
133 const int real_lidx = face_connectivity.localSubcellForSubcell(face,real_side);
136 KOKKOS_ASSERT(real_cell < num_real_cells);
138 for(
int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
141 const int virtual_cell_point = points_per_face * virtual_lidx + point_ordinal;
142 const int real_cell_point = points_per_face * real_lidx + point_ordinal;
144 for (
int d=0; d<space_dim; d++)
145 normals(virtual_cell,virtual_cell_point,d) = -normals(real_cell,real_cell_point,d);
150 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
152 if (local_face_id == virtual_lidx)
continue;
154 for (
int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
155 const int point = local_face_id * points_per_face + point_ordinal;
156 for (
int dim=0; dim<space_dim; dim++)
157 normals(virtual_cell,point,dim) = 0.0;
161 PHX::Device::execution_space().fence();
165 template<
typename Scalar>
168 const int num_real_cells,
169 const int num_virtual_cells,
170 const shards::CellTopology & cell_topology,
171 const SubcellConnectivity & face_connectivity)
173 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::correctVirtualRotationMatrices()",corr_virt_rotmat);
177 const int space_dim = cell_topology.getDimension();
178 const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1);
179 const int points = rotation_matrices.extent_int(1);
180 const int points_per_face = points / faces_per_cell;
182 Kokkos::parallel_for(num_virtual_cells, KOKKOS_LAMBDA(
const int & virtual_cell_ordinal){
184 const int virtual_cell = virtual_cell_ordinal+num_real_cells;
188 int face, virtual_lidx;
189 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
191 face = face_connectivity.subcellForCell(virtual_cell, local_face_id);
193 virtual_lidx = local_face_id;
201 for (
int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
203 if (local_face_id == virtual_lidx)
continue;
205 for (
int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
206 const int point = local_face_id * points_per_face + point_ordinal;
207 for (
int dim=0; dim<3; dim++)
208 for (
int dim2=0; dim2<3; dim2++)
209 rotation_matrices(virtual_cell,point,dim,dim2) = 0.0;
213 PHX::Device::execution_space().fence();
216 template<
typename Scalar>
221 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::applyBasePermutation(rank 1)",app_base_perm_r1);
222 MDFieldArrayFactory af(
"",
true);
224 const int num_ip = field.extent(0);
226 auto scratch = af.template buildStaticArray<Scalar,IP>(
"scratch", num_ip);
227 scratch.deep_copy(field);
228 Kokkos::parallel_for(1, KOKKOS_LAMBDA(
const int & ){
229 for(
int ip=0; ip<num_ip; ++ip)
230 if (ip != permutations(0,ip))
231 field(ip) = scratch(permutations(0,ip));
233 PHX::Device::execution_space().fence();
236 template<
typename Scalar>
241 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::applyBasePermutation(rank 2)",app_base_perm_r2);
242 MDFieldArrayFactory af(
"",
true);
244 const int num_ip = field.extent(0);
245 const int num_dim = field.extent(1);
247 auto scratch = af.template buildStaticArray<Scalar,IP,Dim>(
"scratch", num_ip,num_dim);
248 scratch.deep_copy(field);
249 Kokkos::parallel_for(1, KOKKOS_LAMBDA(
const int & ){
250 for(
int ip=0; ip<num_ip; ++ip)
251 if (ip != permutations(0,ip))
252 for(
int dim=0; dim<num_dim; ++dim)
253 field(ip,dim) = scratch(permutations(0,ip),dim);
255 PHX::Device::execution_space().fence();
258 template<
typename Scalar>
263 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::applyPermutation(rank 2)",app_perm_r2);
264 MDFieldArrayFactory af(
"",
true);
266 const int num_cells = field.extent(0);
267 const int num_ip = field.extent(1);
269 auto scratch = af.template buildStaticArray<Scalar,Cell,IP>(
"scratch", num_cells, num_ip);
270 scratch.deep_copy(field);
271 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
272 for(
int ip=0; ip<num_ip; ++ip)
273 if (ip != permutations(cell,ip))
274 field(cell,ip) = scratch(cell,permutations(cell,ip));
276 PHX::Device::execution_space().fence();
279 template<
typename Scalar>
284 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::applyPermutation(rank 3)",app_perm_r3);
285 MDFieldArrayFactory af(
"",
true);
287 const int num_cells = field.extent(0);
288 const int num_ip = field.extent(1);
289 const int num_dim = field.extent(2);
291 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"scratch", num_cells, num_ip, num_dim);
292 scratch.deep_copy(field);
293 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
294 for(
int ip=0; ip<num_ip; ++ip)
295 if (ip != permutations(cell,ip))
296 for(
int dim=0; dim<num_dim; ++dim)
297 field(cell,ip,dim) = scratch(cell,permutations(cell,ip),dim);
299 PHX::Device::execution_space().fence();
302 template<
typename Scalar>
307 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::applyPermutation(rank 4)",app_perm_r4);
308 MDFieldArrayFactory af(
"",
true);
310 const int num_cells = field.extent(0);
311 const int num_ip = field.extent(1);
312 const int num_dim = field.extent(2);
313 const int num_dim2 = field.extent(3);
315 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"scratch", num_cells, num_ip, num_dim, num_dim2);
316 scratch.deep_copy(field);
317 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
318 for(
int ip=0; ip<num_ip; ++ip)
319 if (ip != permutations(cell,ip))
320 for(
int dim=0; dim<num_dim; ++dim)
321 for(
int dim2=0; dim2<num_dim2; ++dim2)
322 field(cell,ip,dim,dim2) = scratch(cell,permutations(cell,ip),dim,dim2);
324 PHX::Device::execution_space().fence();
331 template <
typename Scalar>
333 generatePermutations(
const int num_cells,
337 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::generatePermutations()",gen_perms);
339 const int num_ip = coords.extent(1);
340 const int num_dim = coords.extent(2);
342 MDFieldArrayFactory af(
"",
true);
345 auto permutation = af.template buildStaticArray<int,Cell,IP>(
"permutation", num_cells, num_ip);
346 permutation.deep_copy(0);
349 auto taken = af.template buildStaticArray<int,Cell,IP>(
"taken", num_cells, num_ip);
352 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(
const int & cell){
354 for (
int ip = 0; ip < num_ip; ++ip) {
358 for (
int other_ip = 0; other_ip < num_ip; ++other_ip) {
360 if(taken(cell,other_ip))
continue;
363 for (
int dim = 0; dim < num_dim; ++dim) {
364 const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim);
367 if (d_min < 0 || d < d_min) {
373 permutation(cell,ip) = i_min;
375 taken(cell,i_min) = 1;
378 PHX::Device::execution_space().fence();
384 template <
typename Scalar>
386 generateSurfacePermutations(
const int num_cells,
387 const SubcellConnectivity face_connectivity,
392 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::generateSurfacePermutations()",gen_surf_perms);
399 const int num_points = surface_points.extent_int(1);
400 const int num_faces_per_cell = face_connectivity.numSubcellsOnCellHost(0);
401 const int num_points_per_face = num_points / num_faces_per_cell;
402 const int cell_dim = surface_points.extent(2);
404 MDFieldArrayFactory af(
"",
true);
407 auto permutation = af.template buildStaticArray<int,Cell,IP>(
"permutation", num_cells, num_points);
408 permutation.deep_copy(0);
411 Kokkos::parallel_for(num_cells,KOKKOS_LAMBDA (
const int & cell) {
412 for(
int point=0; point<num_points; ++point)
413 permutation(cell,point) = point;
418 if(num_points_per_face != 1) {
424 #define PANZER_DOT(a,b) (a[0]*b[0] + a[1]*b[1] + a[2]*b[2])
425 #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];}
428 Kokkos::parallel_for(
"face iteration",face_connectivity.numSubcells(),KOKKOS_LAMBDA (
const int face) {
430 const int cell_0 = face_connectivity.cellForSubcell(face,0);
431 const int cell_1 = face_connectivity.cellForSubcell(face,1);
438 const int lidx_0 = face_connectivity.localSubcellForSubcell(face,0);
439 const int lidx_1 = face_connectivity.localSubcellForSubcell(face,1);
442 KOKKOS_ASSERT(lidx_1 >= 0);
447 Scalar xc0[3] = {0}, xc1[3] = {0}, r2 = 0;
448 for(
int face_point=0; face_point<num_points_per_face; ++face_point){
450 for(
int dim=0; dim<cell_dim; ++dim){
451 xc0[dim] += surface_points(cell_0,lidx_0*num_points_per_face + face_point,dim);
452 xc1[dim] += surface_points(cell_1,lidx_1*num_points_per_face + face_point,dim);
453 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);
458 r2 = (r2 < dx2) ? dx2 : r2;
460 for(
int dim=0; dim<cell_dim; ++dim){
461 xc0[dim] /= double(num_points_per_face);
462 xc1[dim] /= double(num_points_per_face);
471 const int example_point_0 = lidx_0*num_points_per_face;
472 const int example_point_1 = lidx_1*num_points_per_face;
475 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)};
476 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)};
482 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)};
483 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)};
491 if(Kokkos::fabs(n0_dot_n1 - 1.) < 1.e-8)
495 if(Kokkos::fabs(n0_dot_n1 + 1.) > 1.e-2){
502 const Scalar mag_t = Kokkos::sqrt(
PANZER_DOT(t,t));
508 b[0] = n0[0] + n1[0];
509 b[1] = n0[1] + n1[1];
510 b[2] = n0[2] + n1[2];
513 const Scalar mag_b = Kokkos::sqrt(
PANZER_DOT(b,b));
531 for(
int face_point_1=0; face_point_1<num_points_per_face; ++face_point_1){
534 const int point_1 = lidx_1*num_points_per_face + face_point_1;
537 for(
int dim=0; dim<cell_dim; ++dim)
538 x1[dim] = surface_points(cell_1,point_1,dim) - xc1[dim];
545 for(
int face_point_0=0; face_point_0<num_points_per_face; ++face_point_0){
548 const int point_0 = lidx_0*num_points_per_face + face_point_0;
551 for(
int dim=0; dim<cell_dim; ++dim)
552 x0[dim] = surface_points(cell_0,point_0,dim) - xc0[dim];
559 const Scalar p012 = (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]);
563 if(p012 / r2 < 1.e-12){
564 permutation(cell_1,lidx_1*num_points_per_face+face_point_0) = point_1;
571 PHX::Device::execution_space().fence();
585 template <
typename Scalar>
588 const bool allocArrays):
590 num_evaluate_cells_(0),
591 num_virtual_cells_(-1),
592 requires_permutation_(false),
593 alloc_arrays_(allocArrays),
600 template <
typename Scalar>
605 cub_points_evaluated_ =
false;
606 side_cub_points_evaluated_ =
false;
607 cub_weights_evaluated_ =
false;
608 node_coordinates_evaluated_ =
false;
609 jac_evaluated_ =
false;
610 jac_inv_evaluated_ =
false;
611 jac_det_evaluated_ =
false;
612 weighted_measure_evaluated_ =
false;
613 weighted_normals_evaluated_ =
false;
614 surface_normals_evaluated_ =
false;
615 surface_rotation_matrices_evaluated_ =
false;
616 covarient_evaluated_ =
false;
617 contravarient_evaluated_ =
false;
618 norm_contravarient_evaluated_ =
false;
619 ip_coordinates_evaluated_ =
false;
620 ref_ip_coordinates_evaluated_ =
false;
625 template <
typename Scalar>
629 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::setupArrays()",setup_arrays);
637 const int num_nodes = ir->
topology->getNodeCount();
639 const int num_space_dim = ir->
topology->getDimension();
642 const bool is_node_rule = (num_space_dim==1) and ir->
isSide();
643 if(not is_node_rule) {
645 intrepid_cubature = getIntrepidCubature(*ir);
648 const int num_ip = is_node_rule ? 1 : ir->
num_points;
650 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
653 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip,ir->
side_topology->getDimension());
655 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
657 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
659 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells, num_ip, num_space_dim,num_space_dim);
661 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
663 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells, num_ip);
665 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells, num_ip);
667 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells, num_ip, num_space_dim,num_space_dim);
669 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
671 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells, num_ip);
673 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells, num_ip,num_space_dim);
675 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells, num_ip,num_space_dim);
677 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted_normal",num_cells,num_ip,num_space_dim);
679 surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells, num_ip,num_space_dim);
681 surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells, num_ip,3,3);
688 template <
typename Scalar>
691 const int in_num_cells,
693 const int num_virtual_cells)
695 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::evaluateValues(with virtual cells)",eval_vals_with_virts);
697 setup(int_rule, in_node_coordinates, in_num_cells);
701 if((int_rule->getType() == ID::SURFACE) and (face_connectivity != Teuchos::null))
702 setupPermutations(face_connectivity, num_virtual_cells);
705 evaluateEverything();
708 template <
typename Scalar>
713 const int in_num_cells)
715 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::evaluateValues(no virtual cells)",eval_vals_no_virts);
717 setup(int_rule, in_node_coordinates, in_num_cells);
720 setupPermutations(other_ip_coordinates);
723 evaluateEverything();
726 template <
typename Scalar>
730 const int num_virtual_cells)
732 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::setupPermutations(connectivity)",setup_perms_conn);
737 "IntegrationValues2::setupPermutations : Face connectivity based permutations are only required for surface integration schemes");
739 "IntegrationValues2::setupPermutations : Number of virtual cells for surface permuations must be >= 0");
741 side_connectivity_ = face_connectivity;
742 num_virtual_cells_ = num_virtual_cells;
743 requires_permutation_ =
false;
744 permutations_ = generateSurfacePermutations<Scalar>(num_evaluate_cells_,*face_connectivity, getCubaturePoints(
false,
true), getSurfaceRotationMatrices(
false,
true));
745 requires_permutation_ =
true;
748 template <
typename Scalar>
753 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::setupPermutations(other_coords)",setup_perms_coords);
755 requires_permutation_ =
false;
756 permutations_ = generatePermutations<Scalar>(num_evaluate_cells_, getCubaturePoints(
false,
true), other_ip_coordinates);
757 requires_permutation_ =
true;
761 template <
typename Scalar>
768 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::setup()",setup);
773 num_cells_ = cell_node_coordinates.extent(0);
774 num_evaluate_cells_ = num_cells < 0 ? cell_node_coordinates.extent(0) : num_cells;
778 intrepid_cubature = getIntrepidCubature(*ir);
784 const int num_space_dim = int_rule->topology->getDimension();
785 const int num_nodes = int_rule->topology->getNodeCount();
786 TEUCHOS_ASSERT(static_cast<int>(cell_node_coordinates.extent(1)) == num_nodes);
788 auto aux = af.template buildStaticArray<Scalar,Cell,NODE,Dim>(
"node_coordinates",num_cells_,num_nodes, num_space_dim);
789 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_nodes,num_space_dim});
790 Kokkos::parallel_for(policy,KOKKOS_LAMBDA(
const int & cell,
const int & point,
const int & dim){
791 aux(cell,point,dim) = cell_node_coordinates(cell,point,dim);
793 PHX::Device::execution_space().fence();
794 node_coordinates = aux;
799 template <
typename Scalar>
804 const bool apply_permutation)
const
806 if(cub_points_evaluated_ and not force)
810 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getUniformCubaturePointsRef()",get_uniform_cub_pts_ref);
812 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
815 int num_space_dim = int_rule->topology->getDimension();
816 int num_ip = int_rule->num_points;
818 auto aux = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
820 if (int_rule->isSide() && num_space_dim==1) {
821 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
822 <<
"non-natural integration rules.";
827 "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme.");
830 "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for surface integration scheme.");
832 auto weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
834 if (!int_rule->isSide())
835 intrepid_cubature->getCubature(aux.get_view(), weights.get_view());
837 auto s_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip, num_space_dim-1);
839 intrepid_cubature->getCubature(s_cub_points.get_view(), weights.get_view());
840 cell_tools.mapToReferenceSubcell(aux.get_view(), s_cub_points.get_view(), num_space_dim-1, int_rule->getSide(), *(int_rule->topology));
843 PHX::Device::execution_space().fence();
845 if(apply_permutation and requires_permutation_)
846 applyBasePermutation(aux, permutations_);
850 cub_points_evaluated_ =
true;
857 template <
typename Scalar>
862 const bool apply_permutation)
const
864 if(side_cub_points_evaluated_ and not force)
865 return side_cub_points;
868 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getUniformSideCubaturePointsRef()",get_uniform_side_cub_pts_ref);
872 int num_space_dim = int_rule->topology->getDimension();
873 int num_ip = int_rule->num_points;
875 auto aux = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip, num_space_dim-1);
877 if (int_rule->isSide() && num_space_dim==1) {
878 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
879 <<
"non-natural integration rules.";
884 "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme.");
887 "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for surface integration scheme.");
890 "IntegrationValues2::getUniformSideCubaturePointsRef : Requested side points, which is not supported by integration rule.");
892 auto weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
894 intrepid_cubature->getCubature(aux.get_view(), weights.get_view());
896 PHX::Device::execution_space().fence();
898 if(apply_permutation and requires_permutation_)
899 applyBasePermutation(aux, permutations_);
902 side_cub_points = aux;
903 side_cub_points_evaluated_ =
true;
909 template <
typename Scalar>
914 const bool apply_permutation)
const
916 if(cub_weights_evaluated_ and not force)
920 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getUniformCubatureWeightRef()",get_uniform_cub_weights_ref);
924 int num_space_dim = int_rule->topology->getDimension();
925 int num_ip = int_rule->num_points;
927 auto aux = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
929 if (int_rule->isSide() && num_space_dim==1) {
930 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
931 <<
"non-natural integration rules.";
936 "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for control volume integration scheme.");
939 "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for surface integration scheme.");
941 auto points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip,num_space_dim);
943 intrepid_cubature->getCubature(points.get_view(), aux.get_view());
945 PHX::Device::execution_space().fence();
947 if(apply_permutation and requires_permutation_)
948 applyBasePermutation(aux, permutations_);
952 cub_weights_evaluated_ =
true;
959 template <
typename Scalar>
964 return node_coordinates;
967 template <
typename Scalar>
971 const bool force)
const
973 if(jac_evaluated_ and not force)
977 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getJacobian()",get_jacobian);
979 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
982 int num_space_dim = int_rule->topology->getDimension();
983 int num_ip = int_rule->num_points;
986 auto const_ref_coord = getCubaturePointsRef(
false,force);
987 auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord);
988 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
989 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells_, num_ip, num_space_dim,num_space_dim);
991 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
992 auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
993 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
994 auto s_jac = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
996 cell_tools.setJacobian(s_jac, s_ref_coord, s_node_coord,*(int_rule->topology));
998 PHX::Device::execution_space().fence();
1002 jac_evaluated_ =
true;
1008 template <
typename Scalar>
1012 const bool force)
const
1014 if(jac_inv_evaluated_ and not force)
1018 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getJacobianInverse()",get_jacobian_inv);
1020 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1023 const int num_space_dim = int_rule->topology->getDimension();
1024 const int num_ip = int_rule->num_points;
1026 auto jacobian = getJacobian(
false,force);
1027 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells_, num_ip, num_space_dim,num_space_dim);
1029 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1030 auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1031 auto s_jac_inv = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1033 cell_tools.setJacobianInv(s_jac_inv, s_jac);
1035 PHX::Device::execution_space().fence();
1039 jac_inv_evaluated_ =
true;
1045 template <
typename Scalar>
1049 const bool force)
const
1051 if(jac_det_evaluated_ and not force)
1055 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getJacobianDeterminant()",get_jacobian_det);
1057 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1060 const int num_ip = int_rule->num_points;
1062 auto jacobian = getJacobian(
false,force);
1063 auto aux = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells_, num_ip);
1065 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1066 auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1067 auto s_jac_det = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL());
1069 cell_tools.setJacobianDet(s_jac_det, s_jac);
1071 PHX::Device::execution_space().fence();
1075 jac_det_evaluated_ =
true;
1081 template <
typename Scalar>
1085 const bool force)
const
1087 if(weighted_measure_evaluated_ and not force)
1088 return weighted_measure;
1091 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getWeightedMeasure()",get_wt_meas);
1095 const int num_space_dim = int_rule->topology->getDimension();
1096 const int num_ip = int_rule->num_points;
1098 auto aux = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells_, num_ip);
1100 if(int_rule->cv_type !=
"none"){
1103 "IntegrationValues2::getWeightedMeasure : Weighted measure not available for side control volume integrators. Use getWeightedNormals instead.");
1107 auto s_cub_points = af.template buildStaticArray<Scalar, Cell, IP, Dim>(
"cub_points",num_evaluate_cells_,num_ip,num_space_dim);
1109 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1111 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1112 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1113 auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL());
1115 intrepid_cubature->getCubature(s_cub_points.get_view(),s_weighted_measure,s_node_coord);
1119 const auto & cell_topology = *int_rule->topology;
1120 const int cell_dim = cell_topology.getDimension();
1121 const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount();
1122 const int cubature_order = int_rule->order();
1123 const int num_points_on_side = num_ip / num_sides;
1125 Intrepid2::DefaultCubatureFactory cubature_factory;
1126 auto jacobian = getJacobian(
false,force);
1128 for(
int side=0; side<num_sides; ++side) {
1130 const int point_offset=side*num_points_on_side;
1133 Kokkos::DynRankView<double,PHX::Device> side_cub_weights;
1135 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"side_cub_weights",num_points_on_side);
1136 auto tmp_side_cub_weights_host = Kokkos::create_mirror_view(side_cub_weights);
1137 tmp_side_cub_weights_host(0)=1.;
1138 Kokkos::deep_copy(side_cub_weights,tmp_side_cub_weights_host);
1142 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(cell_dim-1,side));
1144 auto ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(face_topology,cubature_order);
1146 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"side_cub_weights",num_points_on_side);
1147 auto subcell_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"subcell_cub_points",num_points_on_side,cell_dim-1);
1150 ic->getCubature(subcell_cub_points, side_cub_weights);
1153 PHX::Device::execution_space().fence();
1156 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_side});
1159 auto side_weighted_measure = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_weighted_measure",num_evaluate_cells_,num_points_on_side);
1161 auto side_cub_weights_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),side_cub_weights);
1162 Kokkos::deep_copy(side_weighted_measure, side_cub_weights_host(0));
1166 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_jac",num_evaluate_cells_,num_points_on_side,cell_dim,cell_dim);
1168 Kokkos::parallel_for(
"Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1169 for(
int dim=0;dim<cell_dim;++dim)
1170 for(
int dim1=0;dim1<cell_dim;++dim1)
1171 side_jacobian(cell,point,dim,dim1) = jacobian(cell,point_offset+point,dim,dim1);
1174 auto scratch = af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_measure", num_evaluate_cells_*num_points_on_side*num_space_dim*num_space_dim);
1177 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1178 computeEdgeMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
1180 scratch.get_view());
1181 PHX::Device::execution_space().fence();
1182 }
else if(cell_dim == 3){
1183 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1184 computeFaceMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
1186 scratch.get_view());
1187 PHX::Device::execution_space().fence();
1193 Kokkos::parallel_for(
"copy surface weighted measure values",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1194 aux(cell,point_offset + point) = side_weighted_measure(cell,point);
1196 PHX::Device::execution_space().fence();
1201 auto cell_range = std::make_pair(0,num_evaluate_cells_);
1202 auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL());
1203 auto cubature_weights = getUniformCubatureWeightsRef(
false,force,
false);
1205 if (!int_rule->isSide()) {
1207 auto s_jac_det = Kokkos::subview(getJacobianDeterminant(
false,force).get_view(),cell_range,Kokkos::ALL());
1208 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1209 computeCellMeasure(s_weighted_measure, s_jac_det, cubature_weights.get_view());
1211 }
else if(int_rule->spatial_dimension==3) {
1213 auto s_jac = Kokkos::subview(getJacobian(
false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1214 auto scratch = af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim);
1215 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1216 computeFaceMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(),
1217 int_rule->side, *int_rule->topology,
1218 scratch.get_view());
1220 }
else if(int_rule->spatial_dimension==2) {
1222 auto s_jac = Kokkos::subview(getJacobian(
false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1223 auto scratch = af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim);
1224 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1225 computeEdgeMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(),
1226 int_rule->side,*int_rule->topology,
1227 scratch.get_view());
1235 PHX::Device::execution_space().fence();
1238 if(requires_permutation_)
1239 applyPermutation(aux, permutations_);
1242 weighted_measure = aux;
1243 weighted_measure_evaluated_ =
true;
1249 template <
typename Scalar>
1253 const bool force)
const
1255 if(weighted_normals_evaluated_ and not force)
1256 return weighted_normals;
1259 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getWeightedNormals()",get_wt_normals);
1263 const int num_space_dim = int_rule->topology->getDimension();
1264 const int num_ip = int_rule->num_points;
1266 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted_normals",num_cells_,num_ip,num_space_dim);
1269 "IntegrationValues2::getWeightedNormals : Cannot build reference weighted normals for surface integration scheme.");
1271 auto points = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"cub_points",num_cells_,num_ip,num_space_dim);
1273 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1275 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1276 auto s_cub_points = Kokkos::subview(points.get_view(),cell_range, Kokkos::ALL(), Kokkos::ALL());
1277 auto s_weighted_normals = Kokkos::subview(aux.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1278 auto s_node_coord = Kokkos::subview(node_coord, cell_range, Kokkos::ALL(), Kokkos::ALL());
1280 intrepid_cubature->getCubature(s_cub_points,s_weighted_normals,s_node_coord);
1282 PHX::Device::execution_space().fence();
1285 if(requires_permutation_)
1286 applyPermutation(aux, permutations_);
1289 weighted_normals = aux;
1290 weighted_normals_evaluated_ =
true;
1296 template <
typename Scalar>
1300 const bool force)
const
1302 if(surface_normals_evaluated_ and not force)
1303 return surface_normals;
1306 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getSurfaceNormals()",get_surf_normals);
1309 "IntegrationValues2::getSurfaceNormals : This call doesn't work with sides (only surfaces).");
1312 "IntegrationValues2::getSurfaceNormals : This call does not support control volume integration schemes.");
1315 "IntegrationValues2::getSurfaceNormals : Can only build for surface integrators.");
1317 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1320 const shards::CellTopology & cell_topology = *(int_rule->topology);
1321 const int cell_dim = cell_topology.getDimension();
1322 const int subcell_dim = cell_topology.getDimension()-1;
1323 const int num_subcells = cell_topology.getSubcellCount(subcell_dim);
1324 const int num_space_dim = int_rule->topology->getDimension();
1325 const int num_ip = int_rule->num_points;
1326 const int num_points_on_face = num_ip / num_subcells;
1328 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells_,num_ip,num_space_dim);
1333 jacobian = getJacobian(
false,force);
1336 for(
int subcell_index=0; subcell_index<num_subcells; ++subcell_index) {
1338 const int point_offset = subcell_index * num_points_on_face;;
1341 auto side_normals = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_normals",num_evaluate_cells_,num_points_on_face,cell_dim);
1344 const int other_subcell_index = (subcell_index==0) ? 1 : 0;
1345 auto in_node_coordinates_k = getNodeCoordinates().get_view();
1346 const auto min = std::numeric_limits<typename Sacado::ScalarType<Scalar>::type>::min();
1348 Kokkos::parallel_for(
"compute normals 1D",num_evaluate_cells_,KOKKOS_LAMBDA (
const int cell) {
1349 Scalar norm = (in_node_coordinates_k(cell,subcell_index,0) - in_node_coordinates_k(cell,other_subcell_index,0));
1350 side_normals(cell,0,0) = norm / fabs(norm + min);
1356 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_face});
1358 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_jac",num_evaluate_cells_,num_points_on_face,cell_dim,cell_dim);
1359 Kokkos::parallel_for(
"Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1360 for(
int dim=0;dim<cell_dim;++dim)
1361 for(
int dim1=0;dim1<cell_dim;++dim1)
1362 side_jacobian(cell,point,dim,dim1) = jacobian(cell,point_offset+point,dim,dim1);
1366 cell_tools.getPhysicalSideNormals(side_normals,side_jacobian,subcell_index,cell_topology);
1369 Kokkos::parallel_for(
"Normalize the normals",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1371 for(
int dim=0;dim<cell_dim;++dim){
1372 n += side_normals(cell,point,dim)*side_normals(cell,point,dim);
1376 n = Kokkos::sqrt(n);
1377 for(
int dim=0;dim<cell_dim;++dim)
1378 side_normals(cell,point,dim) /= n;
1383 PHX::Device::execution_space().fence();
1385 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_face});
1386 Kokkos::parallel_for(
"copy surface normals", policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1387 for(
int dim=0;dim<cell_dim;++dim)
1388 aux(cell,point_offset + point,dim) = side_normals(cell,point,dim);
1390 PHX::Device::execution_space().fence();
1395 TEUCHOS_TEST_FOR_EXCEPT_MSG(side_connectivity_ == Teuchos::null,
1396 "IntegrationValues2::getSurfaceNormals : Surface normals require a SubcellConnectivity object pass in through the setupPermutations call");
1397 TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ < 0,
1398 "IntegrationValues2::getSurfaceNormals : Number of virtual cells (see setupPermutations or evaluateValues) must be greater or equal to zero");
1401 if(num_virtual_cells_ > 0)
1402 correctVirtualNormals(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_);
1406 surface_normals = aux;
1407 surface_normals_evaluated_ =
true;
1414 template <
typename Scalar>
1418 const bool force)
const
1420 if(surface_rotation_matrices_evaluated_ and not force)
1421 return surface_rotation_matrices;
1424 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getSurfaceRotationMatrices()",get_surf_rot_mat);
1428 const int num_ip = int_rule->num_points;
1429 const int cell_dim = int_rule->topology->getDimension();
1433 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells_, num_ip, 3, 3);
1435 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1436 Kokkos::parallel_for(
"create surface rotation matrices",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
1438 for(
int i=0;i<3;i++)
1440 for(
int dim=0; dim<cell_dim; ++dim)
1441 normal[dim] = normals(cell,point,dim);
1443 Scalar transverse[3];
1447 for(
int dim=0; dim<3; ++dim){
1448 aux(cell,point,0,dim) = normal[dim];
1449 aux(cell,point,1,dim) = transverse[dim];
1450 aux(cell,point,2,dim) = binormal[dim];
1453 PHX::Device::execution_space().fence();
1458 "IntegrationValues2::getSurfaceRotationMatrices : Surface normals require a SubcellConnectivity object pass in through the setupPermutations call");
1460 "IntegrationValues2::getSurfaceRotationMatrices : Number of virtual cells (see setupPermutations or evaluateValues) must be greater or equal to zero");
1463 if(num_virtual_cells_ > 0)
1464 correctVirtualRotationMatrices(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_);
1468 surface_rotation_matrices = aux;
1469 surface_rotation_matrices_evaluated_ =
true;
1475 template <
typename Scalar>
1479 const bool force)
const
1481 if(covarient_evaluated_ and not force)
1485 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getCovariantMatrix()",get_cov_mat);
1489 const int num_space_dim = int_rule->topology->getDimension();
1490 const int num_ip = int_rule->num_points;
1492 auto jacobian = getJacobian(
false,force).get_static_view();
1493 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells_, num_ip, num_space_dim,num_space_dim);
1495 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1496 Kokkos::parallel_for(
"evalaute covarient metric tensor",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
1498 for (
int i = 0; i < num_space_dim; ++i) {
1499 for (
int j = 0; j < num_space_dim; ++j) {
1501 for (
int alpha = 0; alpha < num_space_dim; ++alpha)
1502 value += jacobian(cell,ip,i,alpha) * jacobian(cell,ip,j,alpha);
1503 aux(cell,ip,i,j) = value;
1507 PHX::Device::execution_space().fence();
1511 covarient_evaluated_ =
true;
1517 template <
typename Scalar>
1521 const bool force)
const
1523 if(contravarient_evaluated_ and not force)
1524 return contravarient;
1527 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getContravarientMatrix()",get_contra_mat);
1531 const int num_space_dim = int_rule->topology->getDimension();
1532 const int num_ip = int_rule->num_points;
1534 auto cov = getCovarientMatrix(
false,force);
1535 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells_, num_ip, num_space_dim,num_space_dim);
1537 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1538 auto s_contravarient = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
1539 auto s_covarient = Kokkos::subview(cov.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
1541 Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
1542 PHX::Device::execution_space().fence();
1545 contravarient = aux;
1546 contravarient_evaluated_ =
true;
1552 template <
typename Scalar>
1556 const bool force)
const
1558 if(norm_contravarient_evaluated_ and not force)
1559 return norm_contravarient;
1562 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getNormContravarientMatrix()",get_norm_contr_mat);
1566 const int num_space_dim = int_rule->topology->getDimension();
1567 const int num_ip = int_rule->num_points;
1569 auto con = getContravarientMatrix(
false,force).get_static_view();
1570 auto aux = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells_, num_ip);
1573 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1574 Kokkos::parallel_for(
"evaluate norm_contravarient",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
1576 for (
int i = 0; i < num_space_dim; ++i) {
1577 for (
int j = 0; j < num_space_dim; ++j) {
1578 aux(cell,ip) += con(cell,ip,i,j) * con(cell,ip,i,j);
1581 aux(cell,ip) = Kokkos::sqrt(aux(cell,ip));
1583 PHX::Device::execution_space().fence();
1586 norm_contravarient = aux;
1587 norm_contravarient_evaluated_ =
true;
1593 template <
typename Scalar>
1597 const bool force)
const
1599 if(ip_coordinates_evaluated_ and not force)
1600 return ip_coordinates;
1603 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getCubaturePoints()",get_cub_pts);
1607 const int num_space_dim = int_rule->topology->getDimension();
1608 const int num_ip = int_rule->num_points;
1610 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordinates",num_cells_, num_ip, num_space_dim);
1613 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1615 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1620 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1621 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1622 auto s_cub_points = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1625 if(int_rule->cv_type ==
"side"){
1626 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"scratch",num_evaluate_cells_,num_ip,num_space_dim);
1627 intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord);
1630 TEUCHOS_ASSERT((int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_BOUNDARY));
1631 auto scratch = af.template buildStaticArray<Scalar,Cell,IP>(
"scratch",num_evaluate_cells_,num_ip);
1632 intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord);
1638 auto const_ref_coord = getCubaturePointsRef(
false,force);
1639 auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord);
1641 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1642 auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1643 auto s_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1644 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1646 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1647 cell_tools.mapToPhysicalFrame(s_coord, s_ref_coord, s_node_coord, *(int_rule->topology));
1651 PHX::Device::execution_space().fence();
1654 ip_coordinates = aux;
1655 ip_coordinates_evaluated_ =
true;
1662 template <
typename Scalar>
1666 const bool force)
const
1668 if(ref_ip_coordinates_evaluated_)
1669 return ref_ip_coordinates;
1672 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::getCubaturePointsRef()",get_cub_pts_ref);
1675 const bool is_surface = int_rule->
getType() == ID::SURFACE;
1676 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1678 const int num_space_dim = int_rule->topology->getDimension();
1679 const int num_ip = int_rule->num_points;
1683 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1685 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells_, num_ip, num_space_dim);
1691 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1694 auto const_coord = getCubaturePoints(
false,force);
1695 auto coord = PHX::getNonConstDynRankViewFromConstMDField(const_coord);
1697 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1698 auto s_ref_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1699 auto s_coord = Kokkos::subview(coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1700 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1702 cell_tools.mapToReferenceFrame(s_ref_coord, s_coord, s_node_coord, *(int_rule->topology));
1704 }
else if(is_surface){
1706 const auto & cell_topology = *int_rule->topology;
1707 const int cell_dim = cell_topology.getDimension();
1708 const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount();
1709 const int subcell_dim = cell_dim-1;
1710 const int num_points_on_face = num_ip / num_sides;
1711 const int order = int_rule->getOrder();
1714 auto side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_points_on_face,cell_dim);
1716 Intrepid2::DefaultCubatureFactory cubature_factory;
1719 for(
int side=0; side<num_sides; ++side) {
1721 const int point_offset = side*num_points_on_face;
1726 auto side_cub_points_host = Kokkos::create_mirror_view(side_cub_points.get_view());
1727 side_cub_points_host(0,0) = (side==0)? -1. : 1.;
1728 Kokkos::deep_copy(side_cub_points.get_view(),side_cub_points_host);
1732 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,side));
1735 auto ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(face_topology,order);
1736 auto tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"tmp_side_cub_weights",num_points_on_face);
1737 auto tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"tmp_side_cub_points",num_points_on_face,subcell_dim);
1740 ic->getCubature(tmp_side_cub_points, tmp_side_cub_weights);
1743 cell_tools.mapToReferenceSubcell(side_cub_points.get_view(), tmp_side_cub_points, subcell_dim, side, cell_topology);
1746 PHX::Device::execution_space().fence();
1749 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_points_on_face, num_space_dim});
1750 Kokkos::parallel_for(
"copy values",policy,KOKKOS_LAMBDA (
const int cell,
const int point,
const int dim) {
1751 aux(cell,point_offset + point,dim) = side_cub_points(point,dim);
1753 PHX::Device::execution_space().fence();
1760 "ERROR: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
1761 <<
"non-natural integration rules.");
1763 auto cub_points = getUniformCubaturePointsRef(
false,force,
false);
1765 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_ip,num_space_dim});
1766 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(
const int & cell,
const int & ip,
const int & dim){
1767 aux(cell,ip,dim) = cub_points(ip,dim);
1771 PHX::Device::execution_space().fence();
1773 if(requires_permutation_)
1774 applyPermutation(aux, permutations_);
1777 ref_ip_coordinates = aux;
1778 ref_ip_coordinates_evaluated_ =
true;
1784 template <
typename Scalar>
1789 PANZER_FUNC_TIME_MONITOR_DIFF(
"panzer::integrationValues2::evaluateEverything()",eval_everything);
1792 const bool is_surface = int_rule->
getType() == ID::SURFACE;
1793 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1794 const bool is_side = int_rule->isSide();
1800 getCubaturePoints(
true,
true);
1801 getCubaturePointsRef(
true,
true);
1804 getUniformCubaturePointsRef(
true,
true);
1805 getUniformCubatureWeightsRef(
true,
true);
1807 getUniformSideCubaturePointsRef(
true,
true);
1809 getCubaturePointsRef(
true,
true);
1810 getCubaturePoints(
true,
true);
1814 getJacobian(
true,
true);
1815 getJacobianDeterminant(
true,
true);
1816 getJacobianInverse(
true,
true);
1817 if(int_rule->cv_type ==
"side")
1818 getWeightedNormals(
true,
true);
1820 getWeightedMeasure(
true,
true);
1824 getSurfaceNormals(
true,
true);
1825 getSurfaceRotationMatrices(
true,
true);
1829 if(not (is_surface or is_cv)){
1830 getContravarientMatrix(
true,
true);
1831 getCovarientMatrix(
true,
true);
1832 getNormContravarientMatrix(
true,
true);
1836 #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \
1837 template class IntegrationValues2<SCALAR>;
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.
void setupPermutations(const Teuchos::RCP< const SubcellConnectivity > &face_connectivity, const int num_virtual_cells)
Initialize the permutation arrays given a face connectivity.
ConstArray_CellIPDimDim getContravarientMatrix(const bool cache=true, const bool force=false) const
Get the contravarient matrix.
void setup(const Teuchos::RCP< const panzer::IntegrationRule > &ir, const PHX::MDField< Scalar, Cell, NODE, Dim > &cell_node_coordinates, const int num_cells=-1)
Main setup call for the lazy evaluation interface.
ConstArray_CellIP getWeightedMeasure(const bool cache=true, const bool force=false) const
Get the weighted measure (integration weights)
const int & getType() const
Get type of integrator.
void evaluateEverything()
ConstArray_CellIPDimDim getJacobianInverse(const bool cache=true, const bool force=false) const
Get the inverse of the Jacobian matrix evaluated at the cubature points.
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
ConstArray_CellIPDimDim getCovarientMatrix(const bool cache=true, const bool force=false) const
Get the covarient matrix.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
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 getSurfaceRotationMatrices(const bool cache=true, const bool force=false) const
Get the surface rotation matrices.
#define PANZER_CROSS(a, b, c)
Teuchos::RCP< const shards::CellTopology > topology
KOKKOS_INLINE_FUNCTION void convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])
ConstArray_CellIPDim getCubaturePoints(const bool cache=true, const bool force=false) const
Get the cubature points in physical space.
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
ConstArray_CellIPDim getCubaturePointsRef(const bool cache=true, const bool force=false) const
Get the cubature points in the reference space.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
ConstArray_IPDim getUniformCubaturePointsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature points.
ConstArray_CellIPDimDim getJacobian(const bool cache=true, const bool force=false) const
Get 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.
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...
KOKKOS_FORCEINLINE_FUNCTION array_type get_static_view()
IntegrationValues2(const std::string &pre="", const bool allocArrays=false)
Base constructor.
Teuchos::RCP< shards::CellTopology > side_topology
#define TEUCHOS_ASSERT(assertion_test)
ConstArray_CellBASISDim getNodeCoordinates() const
Get the node coordinates describing the geometry of the mesh.
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &cell_node_coordinates, const int num_cells=-1, const Teuchos::RCP< const SubcellConnectivity > &face_connectivity=Teuchos::null, const int num_virtual_cells=-1)
Evaluate basis values.
ConstArray_IP getUniformCubatureWeightsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature weights.
ConstArray_CellIPDim getSurfaceNormals(const bool cache=true, const bool force=false) const
Get the surface normals.
const int & getOrder() const
Get order of integrator.