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" 72 template <
typename Scalar>
78 int num_nodes = ir->topology->getNodeCount();
79 int num_cells = ir->workset_size;
80 int num_space_dim = ir->topology->getDimension();
84 dyn_cub_points = af.template buildArray<double,IP,Dim>(
"cub_points",num_ip, num_space_dim);
85 dyn_cub_weights = af.template buildArray<double,IP>(
"cub_weights",num_ip);
87 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
89 if (ir->cv_type ==
"none" && ir->isSide()) {
90 dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
"side_cub_points",num_ip, ir->side_topology->getDimension());
91 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip,ir->side_topology->getDimension());
94 if (ir->cv_type !=
"none") {
95 dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_points",num_cells, num_ip, num_space_dim);
96 dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>(
"phys_cub_weights",num_cells, num_ip);
97 if (ir->cv_type ==
"side") {
98 dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_norms",num_cells, num_ip, num_space_dim);
102 dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
104 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
106 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
108 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells, num_ip, num_space_dim,num_space_dim);
110 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
112 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells, num_ip);
114 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells, num_ip);
116 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells, num_ip, num_space_dim,num_space_dim);
118 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
120 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells, num_ip);
122 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells, num_ip,num_space_dim);
124 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells, num_ip,num_space_dim);
126 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted normal",num_cells, num_ip,num_space_dim);
128 surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells, num_ip,num_space_dim);
130 surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells, num_ip,3,3);
134 template <
typename Scalar>
136 setupArrays(
const Teuchos::RCP<const panzer::IntegrationRule>& ir)
144 int num_nodes = ir->topology->getNodeCount();
145 int num_cells = ir->workset_size;
146 int num_space_dim = ir->topology->getDimension();
149 if(num_space_dim==1 && ir->isSide()) {
150 setupArraysForNodeRule(ir);
154 TEUCHOS_ASSERT(ir->getType() != ID::NONE);
155 intrepid_cubature = getIntrepidCubature(*ir);
157 int num_ip = ir->num_points;
178 dyn_cub_points = af.template buildArray<double,IP,Dim>(
"cub_points",num_ip, num_space_dim);
179 dyn_cub_weights = af.template buildArray<double,IP>(
"cub_weights",num_ip);
181 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
183 if (ir->isSide() && ir->cv_type ==
"none") {
184 dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
"side_cub_points",num_ip, ir->side_topology->getDimension());
185 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip,ir->side_topology->getDimension());
188 if (ir->cv_type !=
"none") {
189 dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_points",num_cells, num_ip, num_space_dim);
190 dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>(
"phys_cub_weights",num_cells, num_ip);
191 if (ir->cv_type ==
"side") {
192 dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_norms",num_cells, num_ip, num_space_dim);
196 dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
198 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
200 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
202 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells, num_ip, num_space_dim,num_space_dim);
204 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
206 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells, num_ip);
208 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells, num_ip);
210 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells, num_ip, num_space_dim,num_space_dim);
212 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
214 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells, num_ip);
216 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells, num_ip,num_space_dim);
218 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells, num_ip,num_space_dim);
220 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted_normal",num_cells,num_ip,num_space_dim);
222 surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells, num_ip,num_space_dim);
224 surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells, num_ip,3,3);
226 scratch_for_compute_side_measure =
227 af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_side_measure",
jac.get_view().span());
231 template <
typename Scalar>
232 Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double>>
237 Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double> > ic;
239 Intrepid2::DefaultCubatureFactory cubature_factory;
241 if(ir.
getType() == ID::CV_SIDE){
242 ic = Teuchos::rcp(
new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir.
topology));
243 }
else if(ir.
getType() == ID::CV_VOLUME){
244 ic = Teuchos::rcp(
new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir.
topology));
245 }
else if(ir.
getType() == ID::CV_BOUNDARY){
246 ic = Teuchos::rcp(
new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir.
topology,ir.
getSide()));
248 if(ir.
getType() == ID::VOLUME){
249 ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(*(ir.
topology),ir.
getOrder());
250 }
else if(ir.
getType() == ID::SIDE){
251 ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(*(ir.
side_topology),ir.
getOrder());
252 }
else if(ir.
getType() == ID::SURFACE){
255 TEUCHOS_ASSERT(
false);
266 template <
typename Scalar>
269 const int in_num_cells,
270 const Teuchos::RCP<const SubcellConnectivity> & face_connectivity)
273 const bool is_surface = int_rule->
getType() == ID::SURFACE;
274 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
276 TEUCHOS_ASSERT(not (is_surface and is_cv));
279 TEUCHOS_TEST_FOR_EXCEPT_MSG(face_connectivity == Teuchos::null,
280 "IntegrationValues2::evaluateValues : Surface integration requires the face connectivity");
281 generateSurfaceCubatureValues(in_node_coordinates,in_num_cells,*face_connectivity);
283 getCubatureCV(in_node_coordinates, in_num_cells);
284 evaluateValuesCV(in_node_coordinates, in_num_cells);
286 getCubature(in_node_coordinates, in_num_cells);
287 evaluateRemainingValues(in_node_coordinates, in_num_cells);
291 template <
typename Scalar>
293 getCubature(
const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
294 const int in_num_cells)
297 int num_space_dim = int_rule->topology->getDimension();
298 if (int_rule->isSide() && num_space_dim==1) {
299 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " 300 <<
"non-natural integration rules.";
304 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
306 if (!int_rule->isSide())
307 intrepid_cubature->getCubature(dyn_cub_points.get_view(), dyn_cub_weights.get_view());
309 intrepid_cubature->getCubature(dyn_side_cub_points.get_view(), dyn_cub_weights.get_view());
311 cell_tools.mapToReferenceSubcell(dyn_cub_points.get_view(),
312 dyn_side_cub_points.get_view(),
313 int_rule->spatial_dimension-1,
315 *(int_rule->topology));
319 const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
320 auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
321 auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
322 cell_tools.mapToPhysicalFrame(s_ip_coordinates,
323 dyn_cub_points.get_view(),
324 s_in_node_coordinates,
325 *(int_rule->topology));
328 template <
typename Scalar>
331 const int in_num_cells,
337 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
339 const shards::CellTopology & cell_topology = *(int_rule->topology);
342 const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
346 const int num_nodes = in_node_coordinates.extent(1);
347 const int num_dims = in_node_coordinates.extent(2);
348 auto node_coordinates_k = node_coordinates.get_view();
349 auto in_node_coordinates_k = in_node_coordinates.get_view();
351 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<3>> policy({0,0,0},{num_cells,num_nodes,num_dims});
352 Kokkos::parallel_for(
"node_coordinates",policy,KOKKOS_LAMBDA (
const int cell,
const int node,
const int dim) {
353 node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim);
355 PHX::Device::execution_space().fence();
361 const int cell_dim = cell_topology.getDimension();
362 const int subcell_dim = cell_topology.getDimension()-1;
363 const int num_subcells = cell_topology.getSubcellCount(subcell_dim);
365 Intrepid2::DefaultCubatureFactory cubature_factory;
370 for(
int subcell_index=0; subcell_index<num_subcells; ++subcell_index) {
373 int num_points_on_face = 1;
376 Kokkos::DynRankView<double,PHX::Device> tmp_side_cub_weights;
377 Kokkos::DynRankView<double,PHX::Device> tmp_side_cub_points;
379 tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"tmp_side_cub_weights",num_points_on_face);
380 tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"cell_tmp_side_cub_points",num_points_on_face,cell_dim);
381 auto tmp_side_cub_weights_host = Kokkos::create_mirror_view(tmp_side_cub_weights);
382 auto tmp_side_cub_points_host = Kokkos::create_mirror_view(tmp_side_cub_points);
383 tmp_side_cub_weights_host(0)=1.;
384 tmp_side_cub_points_host(0,0) = (subcell_index==0)? -1. : 1.;
385 Kokkos::deep_copy(tmp_side_cub_weights,tmp_side_cub_weights_host);
386 Kokkos::deep_copy(tmp_side_cub_points,tmp_side_cub_points_host);
390 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,subcell_index));
392 auto ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(face_topology,ir.
getOrder());
393 num_points_on_face = ic->getNumPoints();
395 tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"tmp_side_cub_weights",num_points_on_face);
396 tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"cell_tmp_side_cub_points",num_points_on_face,cell_dim);
398 auto subcell_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"subcell_cub_points",num_points_on_face,subcell_dim);
401 ic->getCubature(subcell_cub_points, tmp_side_cub_weights);
404 cell_tools.mapToReferenceSubcell(tmp_side_cub_points,
413 auto tmp_side_cub_points_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),tmp_side_cub_points);
414 auto cub_points_host = Kokkos::create_mirror_view(cub_points.get_static_view());
415 for(
int local_point=0;local_point<num_points_on_face;++local_point){
416 const int point = point_offset + local_point;
417 for(
int dim=0;dim<cell_dim;++dim){
418 cub_points_host(point,dim) = tmp_side_cub_points_host(local_point,dim);
421 Kokkos::deep_copy(cub_points.get_static_view(),cub_points_host);
426 auto side_ip_coordinates = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_ip_coordinates",num_cells,num_points_on_face,cell_dim);
427 auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL);
428 cell_tools.mapToPhysicalFrame(side_ip_coordinates,
434 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
435 cell_tools.setJacobian(side_jacobian,
440 PHX::Device::execution_space().fence();
442 auto side_inverse_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_inv_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
443 cell_tools.setJacobianInv(side_inverse_jacobian, side_jacobian);
445 auto side_det_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_det_jac",num_cells,num_points_on_face);
446 cell_tools.setJacobianDet(side_det_jacobian, side_jacobian);
448 PHX::Device::execution_space().fence();
451 auto side_weighted_measure = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_weighted_measure",num_cells,num_points_on_face);
453 Kokkos::deep_copy(side_weighted_measure, tmp_side_cub_weights(0));
454 }
else if(cell_dim == 2){
455 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
456 computeEdgeMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights,
457 subcell_index,cell_topology,
458 scratch_for_compute_side_measure.get_view());
459 PHX::Device::execution_space().fence();
460 }
else if(cell_dim == 3){
461 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
462 computeFaceMeasure(side_weighted_measure, side_jacobian, tmp_side_cub_weights,
463 subcell_index,cell_topology,
464 scratch_for_compute_side_measure.get_view());
465 PHX::Device::execution_space().fence();
469 auto side_normals = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_normals",num_cells,num_points_on_face,cell_dim);
472 const int other_subcell_index = (subcell_index==0) ? 1 : 0;
473 auto in_node_coordinates_k = in_node_coordinates.get_view();
474 const auto min = std::numeric_limits<typename Sacado::ScalarType<Scalar>::type>::min();
476 Kokkos::parallel_for(
"compute normals 1D",num_cells,KOKKOS_LAMBDA (
const int cell) {
477 Scalar norm = (in_node_coordinates_k(cell,subcell_index,0) - in_node_coordinates_k(cell,other_subcell_index,0));
478 side_normals(cell,0,0) = norm / fabs(norm + min);
483 cell_tools.getPhysicalSideNormals(side_normals,side_jacobian,subcell_index,cell_topology);
486 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_cells,num_points_on_face});
487 Kokkos::parallel_for(
"Normalize the normals",policy,KOKKOS_LAMBDA (
const int cell,
const int point) {
489 for(
int dim=0;dim<cell_dim;++dim){
490 n += side_normals(cell,point,dim)*side_normals(cell,point,dim);
494 n = Kokkos::Experimental::sqrt(n);
495 for(
int dim=0;dim<cell_dim;++dim){
496 side_normals(cell,point,dim) /= n;
502 PHX::Device::execution_space().fence();
506 auto weighted_measure_k = weighted_measure.get_view();
507 auto jac_k =
jac.get_view();
508 auto jac_inv_k = jac_inv.get_view();
509 auto jac_det_k = jac_det.get_view();
510 auto ref_ip_coordinates_k = ref_ip_coordinates.get_view();
511 auto ip_coordinates_k = ip_coordinates.get_view();
512 auto surface_normals_k = surface_normals.get_view();
513 auto cub_points_k = cub_points.get_view();
514 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_cells,num_points_on_face});
515 Kokkos::parallel_for(
"copy values",policy,KOKKOS_LAMBDA (
const int cell,
const int side_point) {
516 const int cell_point = point_offset + side_point;
518 weighted_measure_k(cell,cell_point) = side_weighted_measure(cell,side_point);
519 jac_det_k(cell,cell_point) = side_det_jacobian(cell,side_point);
520 for(
int dim=0;dim<cell_dim;++dim){
521 ref_ip_coordinates_k(cell,cell_point,dim) = cub_points_k(cell_point,dim);
522 ip_coordinates_k(cell,cell_point,dim) = side_ip_coordinates(cell,side_point,dim);
523 surface_normals_k(cell,cell_point,dim) = side_normals(cell,side_point,dim);
525 for(
int dim2=0;dim2<cell_dim;++dim2){
526 jac_k(cell,cell_point,dim,dim2) = side_jacobian(cell,side_point,dim,dim2);
527 jac_inv_k(cell,cell_point,dim,dim2) = side_inverse_jacobian(cell,side_point,dim,dim2);
531 PHX::Device::execution_space().fence();
533 point_offset += num_points_on_face;
540 auto surface_normals_k = surface_normals.get_view();
541 auto surface_rotation_matrices_k = surface_rotation_matrices.get_view();
542 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<3>> policy({0,0,0},{num_cells,num_subcells,num_points});
543 Kokkos::parallel_for(
"create surface rotation matrices",policy,KOKKOS_LAMBDA (
const int cell,
const int subcell_index,
const int point) {
546 for(
int i=0;i<3;i++){normal[i]=0.;}
548 for(
int dim=0; dim<cell_dim; ++dim){
549 normal[dim] = surface_normals_k(cell,point,dim);
552 Scalar transverse[3];
556 for(
int dim=0; dim<3; ++dim){
557 surface_rotation_matrices_k(cell,point,0,dim) = normal[dim];
558 surface_rotation_matrices_k(cell,point,1,dim) = transverse[dim];
559 surface_rotation_matrices_k(cell,point,2,dim) = binormal[dim];
562 PHX::Device::execution_space().fence();
568 const int num_points = ip_coordinates.extent_int(1);
570 const int num_points_per_face = num_points / num_faces_per_cell;
574 if(num_points_per_face != 1){
580 #define PANZER_DOT(a,b) (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) 581 #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];} 589 PHX::View<int**> point_order(
"scratch: point_order",face_connectivity.
numSubcells(),num_points_per_face);
592 auto ref_ip_coordinates_k = ref_ip_coordinates.get_view();
593 auto ip_coordinates_k = ip_coordinates.get_view();
594 auto weighted_measure_k = weighted_measure.get_view();
595 auto jac_k =
jac.get_view();
596 auto jac_det_k = jac_det.get_view();
597 auto jac_inv_k = jac_inv.get_view();
598 auto surface_normals_k = surface_normals.get_view();
599 auto surface_rotation_matrices_k = surface_rotation_matrices.get_view();
600 Kokkos::parallel_for(
"face iteration",face_connectivity.
numSubcells(),KOKKOS_LAMBDA (
const int face) {
614 KOKKOS_ASSERT(lidx_1 >= 0);
619 Scalar xc0[3] = {0}, xc1[3] = {0}, r2 = 0;
620 for(
int face_point=0; face_point<num_points_per_face; ++face_point){
622 for(
int dim=0; dim<cell_dim; ++dim){
623 xc0[dim] += ip_coordinates_k(cell_0,lidx_0*num_points_per_face + face_point,dim);
624 xc1[dim] += ip_coordinates_k(cell_1,lidx_1*num_points_per_face + face_point,dim);
625 const Scalar dx = ip_coordinates_k(cell_0,lidx_0*num_points_per_face + face_point,dim) - ip_coordinates_k(cell_0,lidx_0*num_points_per_face,dim);
630 r2 = (r2 < dx2) ? dx2 : r2;
632 for(
int dim=0; dim<cell_dim; ++dim){
633 xc0[dim] /= double(num_points_per_face);
634 xc1[dim] /= double(num_points_per_face);
643 const int example_point_0 = lidx_0*num_points_per_face;
644 const int example_point_1 = lidx_1*num_points_per_face;
647 Scalar t[3] = {surface_rotation_matrices_k(cell_0,example_point_0,1,0), surface_rotation_matrices_k(cell_0,example_point_0,1,1), surface_rotation_matrices_k(cell_0,example_point_0,1,2)};
648 Scalar b[3] = {surface_rotation_matrices_k(cell_0,example_point_0,2,0), surface_rotation_matrices_k(cell_0,example_point_0,2,1), surface_rotation_matrices_k(cell_0,example_point_0,2,2)};
654 const Scalar n0[3] = {surface_rotation_matrices_k(cell_0,example_point_0,0,0), surface_rotation_matrices_k(cell_0,example_point_0,0,1), surface_rotation_matrices_k(cell_0,example_point_0,0,2)};
655 const Scalar n1[3] = {surface_rotation_matrices_k(cell_1,example_point_1,0,0), surface_rotation_matrices_k(cell_1,example_point_1,0,1), surface_rotation_matrices_k(cell_1,example_point_1,0,2)};
663 if(Kokkos::Experimental::fabs(n0_dot_n1 - 1.) < 1.e-8)
667 if(Kokkos::Experimental::fabs(n0_dot_n1 + 1.) > 1.e-2){
674 const Scalar mag_t = Kokkos::Experimental::sqrt(
PANZER_DOT(t,t));
680 b[0] = n0[0] + n1[0];
681 b[1] = n0[1] + n1[1];
682 b[2] = n0[2] + n1[2];
685 const Scalar mag_b = Kokkos::Experimental::sqrt(
PANZER_DOT(b,b));
703 for(
int face_point_1=0; face_point_1<num_points_per_face; ++face_point_1){
706 const int point_1 = lidx_1*num_points_per_face + face_point_1;
709 for(
int dim=0; dim<cell_dim; ++dim)
710 x1[dim] = ip_coordinates_k(cell_1,point_1,dim) - xc1[dim];
717 point_order(face,face_point_1) = face_point_1;
720 for(
int face_point_0=0; face_point_0<num_points_per_face; ++face_point_0){
723 const int point_0 = lidx_0*num_points_per_face + face_point_0;
726 for(
int dim=0; dim<cell_dim; ++dim)
727 x0[dim] = ip_coordinates_k(cell_0,point_0,dim) - xc0[dim];
734 const Scalar p012 = (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]);
737 if(p012 / r2 < 1.e-12){
738 point_order(face,face_point_1) = face_point_0;
743 KOKKOS_ASSERT(face_point_0 != num_points_per_face-1);
749 const int point_offset = lidx_1*num_points_per_face;
750 for(
int face_point_1 = 0; face_point_1 < num_points_per_face - 1; ++face_point_1 ){
752 while( face_point_1 != point_order(face,face_point_1) ){
754 const int face_point_0 = point_order(face,face_point_1);
755 panzer::swapQuadraturePoints<double>(cell_1,point_offset+face_point_1,point_offset+face_point_0,
756 ref_ip_coordinates_k,
763 surface_rotation_matrices_k);
764 Scalar tmp = point_order(face,face_point_1);
765 point_order(face,face_point_1) = point_order(face,face_point_0);
766 point_order(face,face_point_0) = tmp;
771 PHX::Device::execution_space().fence();
783 auto contravarient_k = this->contravarient.get_static_view();
784 auto covarient_k = this->covarient.get_static_view();
785 auto jac_k = this->
jac.get_static_view();
786 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{(int)num_cells,(
int)contravarient.extent(1)});
787 Kokkos::parallel_for(
"covarient metric tensor",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
789 for (
size_type i = 0; i < contravarient_k.extent(2); ++i)
790 for (
size_type j = 0; j < contravarient_k.extent(3); ++j)
791 covarient_k(cell,ip,i,j) = 0.0;
794 for (
size_type i = 0; i < contravarient_k.extent(2); ++i) {
795 for (
size_type j = 0; j < contravarient_k.extent(3); ++j) {
796 for (
size_type alpha = 0; alpha < contravarient_k.extent(2); ++alpha) {
797 covarient_k(cell,ip,i,j) += jac_k(cell,ip,i,alpha) * jac_k(cell,ip,j,alpha);
802 PHX::Device::execution_space().fence();
806 auto s_contravarient = Kokkos::subview(contravarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
807 auto s_covarient = Kokkos::subview(covarient.get_view(), std::make_pair(0,num_cells),Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
808 Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
809 PHX::Device::execution_space().fence();
814 auto contravarient_k = this->contravarient.get_static_view();
815 auto norm_contravarient_k = this->norm_contravarient.get_static_view();
816 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{(int)num_cells,(
int)contravarient.extent(1)});
817 Kokkos::parallel_for(
"covarient metric tensor",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
818 norm_contravarient_k(cell,ip) = 0.0;
819 for (
size_type i = 0; i < contravarient_k.extent(2); ++i) {
820 for (
size_type j = 0; j < contravarient_k.extent(3); ++j) {
821 norm_contravarient_k(cell,ip) += contravarient_k(cell,ip,i,j) * contravarient_k(cell,ip,i,j);
824 norm_contravarient_k(cell,ip) = Kokkos::Experimental::sqrt(norm_contravarient_k(cell,ip));
826 PHX::Device::execution_space().fence();
832 template <
typename Scalar>
835 const int in_num_cells)
837 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
841 Kokkos::deep_copy(cub_weights.get_static_view(),dyn_cub_weights.get_view());
842 Kokkos::deep_copy(cub_points.get_static_view(),dyn_cub_points.get_view());
845 if (int_rule->isSide()) {
846 Kokkos::deep_copy(side_cub_points.get_static_view(),dyn_side_cub_points.get_view());
849 const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
852 size_type num_nodes = in_node_coordinates.extent(1);
853 size_type num_dims = in_node_coordinates.extent(2);
854 auto node_coordinates_k = node_coordinates.get_view();
855 auto in_node_coordinates_k = in_node_coordinates.get_view();
857 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{(int)num_cells,(
int)num_nodes,(int)num_dims});
858 Kokkos::parallel_for(
"node coordinates",policy,KOKKOS_LAMBDA (
const int cell,
const int node,
const int dim) {
859 node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim);
861 PHX::Device::execution_space().fence();
864 auto s_in_node_coordinates = Kokkos::subview(in_node_coordinates.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
865 auto s_jac = Kokkos::subview(
jac.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
866 cell_tools.setJacobian(
jac.get_view(),
867 cub_points.get_view(),
868 node_coordinates.get_view(),
869 *(int_rule->topology));
870 PHX::Device::execution_space().fence();
872 auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
873 cell_tools.setJacobianInv(s_jac_inv, s_jac);
875 auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair(0,num_cells),Kokkos::ALL());
876 cell_tools.setJacobianDet(s_jac_det, s_jac);
878 PHX::Device::execution_space().fence();
880 auto s_weighted_measure = Kokkos::subview(weighted_measure.get_view(),std::make_pair(0,num_cells),Kokkos::ALL());
881 if (!int_rule->isSide()) {
882 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
883 computeCellMeasure(s_weighted_measure, s_jac_det, cub_weights.get_view());
885 else if(int_rule->spatial_dimension==3) {
886 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
887 computeFaceMeasure(s_weighted_measure, s_jac, cub_weights.get_view(),
888 int_rule->side, *int_rule->topology,
889 scratch_for_compute_side_measure.get_view());
891 else if(int_rule->spatial_dimension==2) {
892 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
893 computeEdgeMeasure(s_weighted_measure, s_jac, cub_weights.get_view(),
894 int_rule->side,*int_rule->topology,
895 scratch_for_compute_side_measure.get_view());
897 else TEUCHOS_ASSERT(
false);
899 PHX::Device::execution_space().fence();
903 auto covarient_k = covarient.get_view();
904 auto contravarient_k = contravarient.get_view();
905 auto jac_k =
jac.get_view();
906 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_cells,
static_cast<int>(contravarient.extent(1))});
907 Kokkos::parallel_for(
"evalaute covarient metric tensor",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
909 for (
int i = 0; i < static_cast<int>(contravarient_k.extent(2)); ++i)
910 for (
int j = 0; j < static_cast<int>(contravarient_k.extent(3)); ++j)
911 covarient_k(cell,ip,i,j) = 0.0;
914 for (
int i = 0; i < static_cast<int>(contravarient_k.extent(2)); ++i) {
915 for (
int j = 0; j < static_cast<int>(contravarient_k.extent(3)); ++j) {
916 for (
int alpha = 0; alpha < static_cast<int>(contravarient_k.extent(2)); ++alpha) {
917 covarient_k(cell,ip,i,j) += jac_k(cell,ip,i,alpha) * jac_k(cell,ip,j,alpha);
922 PHX::Device::execution_space().fence();
925 auto s_covarient = Kokkos::subview(covarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
926 auto s_contravarient = Kokkos::subview(contravarient.get_view(),std::make_pair(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
927 Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
928 PHX::Device::execution_space().fence();
932 auto contravarient_k = contravarient.get_view();
933 auto norm_contravarient_k = norm_contravarient.get_view();
934 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_cells,
static_cast<int>(contravarient.extent(1))});
935 Kokkos::parallel_for(
"evaluate norm_contravarient",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
936 norm_contravarient_k(cell,ip) = 0.0;
937 for (
int i = 0; i < static_cast<int>(contravarient_k.extent(2)); ++i) {
938 for (
int j = 0; j < static_cast<int>(contravarient_k.extent(3)); ++j) {
939 norm_contravarient_k(cell,ip) += contravarient_k(cell,ip,i,j) * contravarient_k(cell,ip,i,j);
942 norm_contravarient_k(cell,ip) = Kokkos::Experimental::sqrt(norm_contravarient_k(cell,ip));
944 PHX::Device::execution_space().fence();
951 template <
typename Scalar>
954 const PHX::MDField<Scalar,Cell,IP,Dim>& other_coords,
955 std::vector<
typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type>& permutation)
961 const size_type cell = 0;
962 const size_type num_ip = coords.extent(1), num_dim = coords.extent(2);
963 permutation.resize(num_ip);
964 std::vector<char> taken(num_ip, 0);
966 auto coords_view = coords.get_view();
967 auto coords_h = Kokkos::create_mirror_view(coords_view);
968 Kokkos::deep_copy(coords_h, coords_view);
970 auto other_coords_view = other_coords.get_view();
971 auto other_coords_h = Kokkos::create_mirror_view(other_coords_view);
972 Kokkos::deep_copy(other_coords_h, other_coords_view);
974 for (size_type ip = 0; ip < num_ip; ++ip) {
978 for (size_type other_ip = 0; other_ip < num_ip; ++other_ip) {
980 if (taken[other_ip])
continue;
983 for (size_type dim = 0; dim < num_dim; ++dim) {
984 const Scalar diff = coords_h(cell, ip, dim) - other_coords_h(cell, other_ip, dim);
987 if (d_min < 0 || d < d_min) {
993 permutation[ip] = i_min;
999 template <
typename Scalar>
1002 const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates,
1003 const int in_num_cells)
1005 const int num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : in_num_cells;
1007 if (int_rule->cv_type ==
"none") {
1009 getCubature(in_node_coordinates, in_num_cells);
1013 std::vector<size_type> permutation(other_ip_coordinates.extent(1));
1014 permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
1018 const size_type num_ip = dyn_cub_points.extent(0);
1020 const size_type num_dim = dyn_side_cub_points.extent(1);
1021 DblArrayDynamic old_dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
1022 "old_dyn_side_cub_points", num_ip, num_dim);
1023 old_dyn_side_cub_points.deep_copy(dyn_side_cub_points);
1025 auto dyn_side_cub_points_h = Kokkos::create_mirror_view(as_view(dyn_side_cub_points));
1026 auto old_dyn_side_cub_points_h = Kokkos::create_mirror_view(as_view(old_dyn_side_cub_points));
1027 Kokkos::deep_copy(dyn_side_cub_points_h, as_view(dyn_side_cub_points));
1028 Kokkos::deep_copy(old_dyn_side_cub_points_h, as_view(old_dyn_side_cub_points));
1030 for (
size_type ip = 0; ip < num_ip; ++ip)
1031 if (ip != permutation[ip])
1032 for (
size_type dim = 0; dim < num_dim; ++dim)
1033 dyn_side_cub_points_h(ip, dim) = old_dyn_side_cub_points_h(permutation[ip], dim);
1035 Kokkos::deep_copy(as_view(dyn_side_cub_points), dyn_side_cub_points_h);
1038 const size_type num_dim = dyn_cub_points.extent(1);
1039 DblArrayDynamic old_dyn_cub_points = af.template buildArray<double,IP,Dim>(
1040 "old_dyn_cub_points", num_ip, num_dim);
1041 old_dyn_cub_points.deep_copy(dyn_cub_points);
1043 auto dyn_cub_points_h = Kokkos::create_mirror_view(as_view(dyn_cub_points));
1044 auto old_dyn_cub_points_h = Kokkos::create_mirror_view(as_view(old_dyn_cub_points));
1045 Kokkos::deep_copy(dyn_cub_points_h, as_view(dyn_cub_points));
1046 Kokkos::deep_copy(old_dyn_cub_points_h, as_view(old_dyn_cub_points));
1048 for (
size_type ip = 0; ip < num_ip; ++ip)
1049 if (ip != permutation[ip])
1050 for (
size_type dim = 0; dim < num_dim; ++dim)
1051 dyn_cub_points_h(ip, dim) = old_dyn_cub_points_h(permutation[ip], dim);
1053 Kokkos::deep_copy(as_view(dyn_cub_points), dyn_cub_points_h);
1056 DblArrayDynamic old_dyn_cub_weights = af.template buildArray<double,IP>(
1057 "old_dyn_cub_weights", num_ip);
1058 old_dyn_cub_weights.deep_copy(dyn_cub_weights);
1060 auto dyn_cub_weights_h = Kokkos::create_mirror_view(as_view(dyn_cub_weights));
1061 auto old_dyn_cub_weights_h = Kokkos::create_mirror_view(as_view(old_dyn_cub_weights));
1062 Kokkos::deep_copy(dyn_cub_weights_h, as_view(dyn_cub_weights));
1063 Kokkos::deep_copy(old_dyn_cub_weights_h, as_view(old_dyn_cub_weights));
1065 for (
size_type ip = 0; ip < dyn_cub_weights.extent(0); ++ip)
1066 if (ip != permutation[ip])
1067 dyn_cub_weights_h(ip) = old_dyn_cub_weights_h(permutation[ip]);
1069 Kokkos::deep_copy(as_view(dyn_cub_weights), old_dyn_cub_weights_h);
1073 const size_type num_ip = ip_coordinates.extent(1);
1074 const size_type num_dim = ip_coordinates.extent(2);
1075 Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1076 "old_ip_coordinates", ip_coordinates.extent(0), num_ip, num_dim);
1078 auto ip_coordinates_h = Kokkos::create_mirror_view(as_view(ip_coordinates));
1079 auto old_ip_coordinates_h = Kokkos::create_mirror_view(as_view(old_ip_coordinates));
1080 Kokkos::deep_copy(old_ip_coordinates_h, as_view(ip_coordinates));
1081 Kokkos::deep_copy(ip_coordinates_h, as_view(ip_coordinates));
1083 for (
int cell = 0; cell < num_cells; ++cell)
1084 for (
size_type ip = 0; ip < num_ip; ++ip)
1085 if (ip != permutation[ip])
1086 for (
size_type dim = 0; dim < num_dim; ++dim)
1087 ip_coordinates_h(cell, ip, dim) = old_ip_coordinates_h(cell, permutation[ip], dim);
1089 Kokkos::deep_copy(as_view(ip_coordinates), ip_coordinates_h);
1094 evaluateRemainingValues(in_node_coordinates, in_num_cells);
1099 getCubatureCV(in_node_coordinates, in_num_cells);
1102 std::vector<size_type> permutation(other_ip_coordinates.extent(1));
1103 permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
1108 const size_type workset_size = ip_coordinates.extent(0), num_ip = ip_coordinates.extent(1),
1109 num_dim = ip_coordinates.extent(2);
1110 Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1111 "old_ip_coordinates", workset_size, num_ip, num_dim);
1112 Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
1113 Array_CellIPDim old_weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1114 "old_weighted_normals", workset_size, num_ip, num_dim);
1115 Array_CellIP old_weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
1116 "old_weighted_measure", workset_size, num_ip);
1117 if (int_rule->cv_type ==
"side")
1118 Kokkos::deep_copy(old_weighted_normals.get_static_view(), weighted_normals.get_static_view());
1120 Kokkos::deep_copy(old_weighted_measure.get_static_view(), weighted_measure.get_static_view());
1121 for (
int cell = 0; cell < num_cells; ++cell)
1123 for (
size_type ip = 0; ip < num_ip; ++ip)
1125 if (ip != permutation[ip]) {
1126 if (int_rule->cv_type ==
"boundary" || int_rule->cv_type ==
"volume")
1127 weighted_measure(cell, ip) = old_weighted_measure(cell, permutation[ip]);
1128 for (
size_type dim = 0; dim < num_dim; ++dim)
1130 ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
1131 if (int_rule->cv_type ==
"side")
1132 weighted_normals(cell, ip, dim) = old_weighted_normals(cell, permutation[ip], dim);
1140 evaluateValuesCV(in_node_coordinates, in_num_cells);
1144 template <
typename Scalar>
1147 const int in_num_cells)
1149 int num_space_dim = int_rule->topology->getDimension();
1150 if (int_rule->isSide() && num_space_dim==1) {
1151 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " 1152 <<
"non-natural integration rules.";
1156 size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (
size_type) in_num_cells;
1157 std::pair<int,int> cell_range(0,num_cells);
1159 size_type num_nodes = in_node_coordinates.extent(1);
1160 size_type num_dims = in_node_coordinates.extent(2);
1161 auto node_coordinates_k = node_coordinates.get_view();
1162 auto dyn_node_coordinates_k = dyn_node_coordinates.get_view();
1163 auto in_node_coordinates_k = in_node_coordinates.get_view();
1165 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_cells,num_nodes,num_dims});
1166 Kokkos::parallel_for(
"getCubatureCV: node coordinates",policy,KOKKOS_LAMBDA (
const int cell,
const int node,
const int dim) {
1167 node_coordinates_k(cell,node,dim) = in_node_coordinates_k(cell,node,dim);
1168 dyn_node_coordinates_k(cell,node,dim) = Sacado::ScalarValue<Scalar>::eval(in_node_coordinates_k(cell,node,dim));
1170 PHX::Device::execution_space().fence();
1173 auto s_dyn_phys_cub_points = Kokkos::subdynrankview(dyn_phys_cub_points.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1174 auto s_dyn_node_coordinates = Kokkos::subdynrankview(dyn_node_coordinates.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1175 if (int_rule->cv_type ==
"side") {
1176 auto s_dyn_phys_cub_norms = Kokkos::subdynrankview(dyn_phys_cub_norms.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1177 intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_norms,s_dyn_node_coordinates);
1180 auto s_dyn_phys_cub_weights = Kokkos::subdynrankview(dyn_phys_cub_weights.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1181 intrepid_cubature->getCubature(s_dyn_phys_cub_points,s_dyn_phys_cub_weights,s_dyn_node_coordinates);
1185 size_type num_ip =dyn_phys_cub_points.extent(1);
1186 size_type num_dims = dyn_phys_cub_points.extent(2);
1187 auto weighted_measure_k = weighted_measure.get_view();
1188 auto dyn_phys_cub_weights_k = dyn_phys_cub_weights.get_view();
1189 auto ip_coordinates_k = ip_coordinates.get_view();
1190 auto dyn_phys_cub_points_k = dyn_phys_cub_points.get_view();
1191 auto weighted_normals_k = weighted_normals.get_view();
1192 auto dyn_phys_cub_norms_k = dyn_phys_cub_norms.get_view();
1193 bool is_side =
false;
1194 if (int_rule->cv_type ==
"side")
1197 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_cells,num_ip});
1198 Kokkos::parallel_for(
"getCubatureCV: weighted measure",policy,KOKKOS_LAMBDA (
const int cell,
const int ip) {
1200 weighted_measure_k(cell,ip) = dyn_phys_cub_weights_k(cell,ip);
1201 for (
size_type dim = 0; dim < num_dims; ++dim) {
1202 ip_coordinates_k(cell,ip,dim) = dyn_phys_cub_points_k(cell,ip,dim);
1204 weighted_normals_k(cell,ip,dim) = dyn_phys_cub_norms_k(cell,ip,dim);
1207 PHX::Device::execution_space().fence();
1210 template <
typename Scalar>
1213 const int in_num_cells)
1216 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1218 size_type num_cells = in_num_cells < 0 ? in_node_coordinates.extent(0) : (
size_type) in_num_cells;
1220 auto s_ref_ip_coordinates = Kokkos::subview(ref_ip_coordinates.get_view(),std::make_pair(0,(
int)num_cells),Kokkos::ALL(),Kokkos::ALL());
1221 auto s_ip_coordinates = Kokkos::subview(ip_coordinates.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
1222 auto s_node_coordinates = Kokkos::subview(node_coordinates.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL());
1224 cell_tools.mapToReferenceFrame(s_ref_ip_coordinates,
1227 *(int_rule->topology));
1229 auto s_jac = Kokkos::subview(
jac.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1231 cell_tools.setJacobian(s_jac,
1232 s_ref_ip_coordinates,
1234 *(int_rule->topology));
1236 auto s_jac_inv = Kokkos::subview(jac_inv.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1238 cell_tools.setJacobianInv(s_jac_inv, s_jac);
1240 auto s_jac_det = Kokkos::subview(jac_det.get_view(),std::make_pair<int,int>(0,num_cells),Kokkos::ALL());
1242 cell_tools.setJacobianDet(s_jac_det, s_jac);
1245 #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \ 1246 template class IntegrationValues2<SCALAR>;
KOKKOS_INLINE_FUNCTION int cellForSubcell(const int subcell, const int local_cell_index) const
Get the cell for a given subcell and a local_cell_index.
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)
Evaluate basis values.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
void generateSurfaceCubatureValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells, const SubcellConnectivity &face_connectivity)
This should be a private method, but using lambdas on cuda forces this to be public.
KOKKOS_INLINE_FUNCTION int numSubcells() const
Gives number of subcells (e.g. faces) in connectivity.
PHX::MDField< Scalar, Cell, IP > Array_CellIP
const int & getType() const
Get type of integrator.
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int in_num_cells)
void getCubatureCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
This should be a private method, but using lambdas on cuda forces this to be public.
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
int getPointOffset(const int subcell_index) const
Returns the integration point offset for a given subcell_index (i.e. local face index) ...
#define PANZER_CROSS(a, b, c)
Teuchos::RCP< Intrepid2::Cubature< PHX::Device::execution_space, double, double > > getIntrepidCubature(const panzer::IntegrationRule &ir) const
Teuchos::RCP< const shards::CellTopology > topology
void evaluateRemainingValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
This should be a private method, but using lambdas on cuda forces this to be public.
ArrayTraits< Scalar, PHX::MDField< Scalar > >::size_type size_type
KOKKOS_INLINE_FUNCTION void convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
void getCubature(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates, const int in_num_cells)
PHX::MDField< Scalar, Cell, IP, Dim > Array_CellIPDim
PHX::MDField< double > DblArrayDynamic
Teuchos::RCP< shards::CellTopology > side_topology
static void permuteToOther(const PHX::MDField< Scalar, Cell, IP, Dim > &coords, const PHX::MDField< Scalar, Cell, IP, Dim > &other_coords, std::vector< typename ArrayTraits< Scalar, PHX::MDField< Scalar > >::size_type > &permutation)
int numSubcellsOnCellHost(const int cell) const
void setupArraysForNodeRule(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
const int & getOrder() const
Get order of integrator.
KOKKOS_INLINE_FUNCTION int localSubcellForSubcell(const int subcell, const int local_cell_index) const
Get the local subcell index given a subcell and a local cell index.