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" 67 template <
typename Scalar>
73 int num_nodes = ir->
topology->getNodeCount();
75 int num_space_dim = ir->
topology->getDimension();
79 dyn_cub_points = af.template buildArray<double,IP,Dim>(
"cub_points",num_ip, num_space_dim);
80 dyn_cub_weights = af.template buildArray<double,IP>(
"cub_weights",num_ip);
82 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
85 dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
"side_cub_points",num_ip, ir->
side_topology->getDimension());
86 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip,ir->
side_topology->getDimension());
90 dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_points",num_cells, num_ip, num_space_dim);
91 dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>(
"phys_cub_weights",num_cells, num_ip);
93 dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_norms",num_cells, num_ip, num_space_dim);
97 dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
99 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
101 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
103 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells, num_ip, num_space_dim,num_space_dim);
105 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
107 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells, num_ip);
109 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells, num_ip);
111 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells, num_ip, num_space_dim,num_space_dim);
113 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
115 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells, num_ip);
117 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells, num_ip,num_space_dim);
119 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells, num_ip,num_space_dim);
121 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted normal",num_cells, num_ip,num_space_dim);
123 surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells, num_ip,num_space_dim);
125 surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells, num_ip,3,3);
129 template <
typename Scalar>
139 int num_nodes = ir->
topology->getNodeCount();
141 int num_space_dim = ir->
topology->getDimension();
144 if(num_space_dim==1 && ir->
isSide()) {
145 setupArraysForNodeRule(ir);
150 intrepid_cubature = getIntrepidCubature(*ir);
173 dyn_cub_points = af.template buildArray<double,IP,Dim>(
"cub_points",num_ip, num_space_dim);
174 dyn_cub_weights = af.template buildArray<double,IP>(
"cub_weights",num_ip);
176 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",num_ip, num_space_dim);
179 dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
"side_cub_points",num_ip, ir->
side_topology->getDimension());
180 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",num_ip,ir->
side_topology->getDimension());
184 dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_points",num_cells, num_ip, num_space_dim);
185 dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>(
"phys_cub_weights",num_cells, num_ip);
187 dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_norms",num_cells, num_ip, num_space_dim);
191 dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
193 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",num_ip);
195 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells, num_nodes, num_space_dim);
197 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells, num_ip, num_space_dim,num_space_dim);
199 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
201 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells, num_ip);
203 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells, num_ip);
205 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells, num_ip, num_space_dim,num_space_dim);
207 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
209 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells, num_ip);
211 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells, num_ip,num_space_dim);
213 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells, num_ip,num_space_dim);
215 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted_normal",num_cells,num_ip,num_space_dim);
217 surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"surface_normals",num_cells, num_ip,num_space_dim);
219 surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"surface_rotation_matrices",num_cells, num_ip,3,3);
221 scratch_for_compute_side_measure =
222 af.template buildStaticArray<Scalar,Point>(
"scratch_for_compute_side_measure",
jac.get_view().span());
226 template <
typename Scalar>
234 Intrepid2::DefaultCubatureFactory cubature_factory;
236 if(ir.
getType() == ID::CV_SIDE){
237 ic =
Teuchos::rcp(
new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir.
topology));
238 }
else if(ir.
getType() == ID::CV_VOLUME){
239 ic =
Teuchos::rcp(
new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir.
topology));
240 }
else if(ir.
getType() == ID::CV_BOUNDARY){
241 ic =
Teuchos::rcp(
new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir.
topology,ir.
getSide()));
243 if(ir.
getType() == ID::VOLUME){
244 ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(*(ir.
topology),ir.
getOrder());
245 }
else if(ir.
getType() == ID::SIDE){
246 ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(*(ir.
side_topology),ir.
getOrder());
247 }
else if(ir.
getType() == ID::SURFACE){
261 template <
typename Scalar>
266 const bool is_surface = int_rule->
getType() == ID::SURFACE;
267 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
272 generateSurfaceCubatureValues(in_node_coordinates);
274 getCubatureCV(in_node_coordinates);
275 evaluateValuesCV(in_node_coordinates);
277 getCubature(in_node_coordinates);
278 evaluateRemainingValues(in_node_coordinates);
282 template <
typename Scalar>
284 getCubature(
const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates)
287 int num_space_dim = int_rule->topology->getDimension();
288 if (int_rule->isSide() && num_space_dim==1) {
289 std::cout <<
"WARNING: 0-D quadrature rule ifrastructure does not exist!!! Will not be able to do " 290 <<
"non-natural integration rules.";
294 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
296 if (!int_rule->isSide())
297 intrepid_cubature->getCubature(dyn_cub_points.get_view(), dyn_cub_weights.get_view());
299 intrepid_cubature->getCubature(dyn_side_cub_points.get_view(), dyn_cub_weights.get_view());
301 cell_tools.mapToReferenceSubcell(dyn_cub_points.get_view(),
302 dyn_side_cub_points.get_view(),
303 int_rule->spatial_dimension-1,
305 *(int_rule->topology));
309 cell_tools.mapToPhysicalFrame(ip_coordinates.get_view(),
310 dyn_cub_points.get_view(),
311 in_node_coordinates.get_view(),
312 *(int_rule->topology));
322 template <
typename array_t,
typename scalar_t>
327 point_sorter_t() =
delete;
328 point_sorter_t(
const array_t & array,
const int cell,
const int offset):
339 bool operator()(
const int & point_a,
const int & point_b)
const 347 const scalar_t rel = std::max(std::fabs(x_a),std::fabs(x_b));
349 return test_less(x_a,x_b,rel);
359 const scalar_t rel_x = std::max(std::fabs(x_a),std::fabs(x_b));
360 const scalar_t rel_y = std::max(std::fabs(y_a),std::fabs(y_b));
361 const scalar_t rel = std::max(rel_x,rel_y);
363 if(test_eq(y_a,y_b,rel)){
364 if(test_less(x_a,x_b,rel)){
368 }
else if(test_less(y_a,y_b,rel)){
387 const scalar_t rel_x = std::max(std::fabs(x_a),std::fabs(x_b));
388 const scalar_t rel_y = std::max(std::fabs(y_a),std::fabs(y_b));
389 const scalar_t rel_z = std::max(std::fabs(z_a),std::fabs(z_b));
390 const scalar_t rel = std::max(rel_x,std::max(rel_y,rel_z));
392 if(test_less(z_a,z_b,rel)){
395 }
else if(test_eq(z_a,z_b,rel)){
396 if(test_eq(y_a,y_b,rel)){
397 if(test_less(x_a,x_b,rel)){
401 }
else if(test_less(y_a,y_b,rel)){
417 test_eq(
const scalar_t & a,
const scalar_t & b,
const scalar_t & rel)
const 422 return std::fabs(a-b) <
_rel_tol * rel;
426 test_less(
const scalar_t & a,
const scalar_t & b,
const scalar_t & rel)
const 444 convertNormalToRotationMatrix(
const T normal[3], T transverse[3], T binormal[3])
447 const T n = sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
454 transverse[0]=0.;transverse[1]=1.;transverse[2]=0.;
455 if(std::fabs(normal[0]*transverse[0]+normal[1]*transverse[1])>0.9){
456 transverse[0]=1.;transverse[1]=0.;
459 const T nt = normal[0]*transverse[0]+normal[1]*transverse[1]+normal[2]*transverse[2];
462 const T mult = nt/(n*n);
465 for(
int dim=0;dim<3;++dim){
466 transverse[dim] = transverse[dim] - mult * normal[dim];
469 const T t = sqrt(transverse[0]*transverse[0]+transverse[1]*transverse[1]+transverse[2]*transverse[2]);
471 for(
int dim=0;dim<3;++dim){
472 transverse[dim] /= t;
476 binormal[0] = (normal[1] * transverse[2] - normal[2] * transverse[1]);
477 binormal[1] = (normal[2] * transverse[0] - normal[0] * transverse[2]);
478 binormal[2] = (normal[0] * transverse[1] - normal[1] * transverse[0]);
481 const T b = sqrt(binormal[0]*binormal[0]+binormal[1]*binormal[1]+binormal[2]*binormal[2]);
482 for(
int dim=0;dim<3;++dim){
500 template <
typename Scalar>
506 const int new_cell_point = a;
507 const int old_cell_point = b;
509 const int cell_dim = ref_ip_coordinates.extent(2);
513 hold = weighted_measure(cell,new_cell_point);
514 weighted_measure(cell,new_cell_point) = weighted_measure(cell,old_cell_point);
515 weighted_measure(cell,old_cell_point) = hold;
517 hold = jac_det(cell,new_cell_point);
518 jac_det(cell,new_cell_point) = jac_det(cell,old_cell_point);
519 jac_det(cell,old_cell_point) = hold;
521 for(
int dim=0;dim<cell_dim;++dim){
523 hold = ref_ip_coordinates(cell,new_cell_point,dim);
524 ref_ip_coordinates(cell,new_cell_point,dim) = ref_ip_coordinates(cell,old_cell_point,dim);
525 ref_ip_coordinates(cell,old_cell_point,dim) = hold;
527 hold = ip_coordinates(cell,new_cell_point,dim);
528 ip_coordinates(cell,new_cell_point,dim) = ip_coordinates(cell,old_cell_point,dim);
529 ip_coordinates(cell,old_cell_point,dim) = hold;
531 hold = surface_normals(cell,new_cell_point,dim);
532 surface_normals(cell,new_cell_point,dim) = surface_normals(cell,old_cell_point,dim);
533 surface_normals(cell,old_cell_point,dim) = hold;
535 for(
int dim2=0;dim2<cell_dim;++dim2){
537 hold =
jac(cell,new_cell_point,dim,dim2);
538 jac(cell,new_cell_point,dim,dim2) =
jac(cell,old_cell_point,dim,dim2);
539 jac(cell,old_cell_point,dim,dim2) = hold;
541 hold = jac_inv(cell,new_cell_point,dim,dim2);
542 jac_inv(cell,new_cell_point,dim,dim2) = jac_inv(cell,old_cell_point,dim,dim2);
543 jac_inv(cell,old_cell_point,dim,dim2) = hold;
548 template <
typename Scalar>
553 std::vector<int> & order)
555 for(
size_t point_index=0;point_index<order.size();++point_index){
556 order[point_index] = point_index;
562 point_sorter_t<Array_CellIPDim,Scalar> sorter(coords,cell,offset);
563 std::sort(order.begin(),order.end(),sorter);
566 template <
typename Scalar>
573 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
575 const shards::CellTopology & cell_topology = *(int_rule->topology);
580 const int num_cells = in_node_coordinates.extent(0);
581 const int num_nodes = in_node_coordinates.extent(1);
582 const int num_dims = in_node_coordinates.extent(2);
584 for(
int cell=0; cell<num_cells; ++cell){
585 for(
int node=0; node<num_nodes; ++node){
586 for(
int dim=0; dim<num_dims; ++dim){
587 node_coordinates(cell,node,dim) = in_node_coordinates(cell,node,dim);
596 const int num_cells = in_node_coordinates.extent(0);
597 const int cell_dim = cell_topology.getDimension();
598 const int subcell_dim = cell_topology.getDimension()-1;
599 const int num_subcells = cell_topology.getSubcellCount(subcell_dim);
601 Intrepid2::DefaultCubatureFactory cubature_factory;
605 for(
int subcell_index=0; subcell_index<num_subcells; ++subcell_index) {
608 int num_points_on_face = 1;
611 Kokkos::DynRankView<double,PHX::Device> side_cub_weights;
612 Kokkos::DynRankView<double,PHX::Device> side_cub_points;
614 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"side_cub_weights",num_points_on_face);
615 side_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"cell_side_cub_points",num_points_on_face,cell_dim);
616 side_cub_weights(0)=1.;
617 side_cub_points(0,0) = (subcell_index==0)? -1. : 1.;
621 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,subcell_index));
623 auto ic = cubature_factory.create<PHX::Device::execution_space,double,
double>(face_topology,ir.
getOrder());
624 num_points_on_face = ic->getNumPoints();
626 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>(
"side_cub_weights",num_points_on_face);
627 side_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"cell_side_cub_points",num_points_on_face,cell_dim);
629 auto subcell_cub_points = Kokkos::DynRankView<double,PHX::Device>(
"side_cub_points",num_points_on_face,subcell_dim);
632 ic->getCubature(subcell_cub_points, side_cub_weights);
635 cell_tools.mapToReferenceSubcell(side_cub_points,
643 for(
int local_point=0;local_point<num_points_on_face;++local_point){
644 const int point = point_offset + local_point;
645 for(
int dim=0;dim<cell_dim;++dim){
646 cub_points(point,dim) = side_cub_points(local_point,dim);
652 auto side_ip_coordinates = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_ip_coordinates",num_cells,num_points_on_face,cell_dim);
653 cell_tools.mapToPhysicalFrame(side_ip_coordinates,
655 node_coordinates.get_view(),
659 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
660 cell_tools.setJacobian(side_jacobian,
662 node_coordinates.get_view(),
665 auto side_inverse_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_inv_jac",num_cells,num_points_on_face,cell_dim,cell_dim);
666 cell_tools.setJacobianInv(side_inverse_jacobian, side_jacobian);
668 auto side_det_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_det_jac",num_cells,num_points_on_face);
669 cell_tools.setJacobianDet(side_det_jacobian, side_jacobian);
672 auto side_weighted_measure = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_weighted_measure",num_cells,num_points_on_face);
674 Kokkos::deep_copy(side_weighted_measure, side_cub_weights(0));
675 }
else if(cell_dim == 2){
676 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
677 computeEdgeMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
678 subcell_index,cell_topology,
679 scratch_for_compute_side_measure.get_view());
680 }
else if(cell_dim == 3){
681 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
682 computeFaceMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
683 subcell_index,cell_topology,
684 scratch_for_compute_side_measure.get_view());
688 auto side_normals = Kokkos::DynRankView<Scalar,PHX::Device>(
"side_normals",num_cells,num_points_on_face,cell_dim);
691 int other_subcell_index = (subcell_index==0) ? 1 : 0;
693 for(
int cell=0;cell<num_cells;++cell){
694 Scalar norm = (in_node_coordinates(cell,subcell_index,0) - in_node_coordinates(cell,other_subcell_index,0));
695 side_normals(cell,0,0) = norm / fabs(norm);
700 cell_tools.getPhysicalSideNormals(side_normals,side_jacobian,subcell_index,cell_topology);
703 for(
int cell=0;cell<num_cells;++cell){
704 for(
int point=0;point<num_points_on_face;++point){
706 for(
int dim=0;dim<cell_dim;++dim){
707 n += side_normals(cell,point,dim)*side_normals(cell,point,dim);
712 for(
int dim=0;dim<cell_dim;++dim){
713 side_normals(cell,point,dim) /= n;
723 for(
int cell=0;cell<num_cells;++cell){
724 for(
int side_point=0; side_point<num_points_on_face;++side_point){
725 const int cell_point = point_offset + side_point;
727 weighted_measure(cell,cell_point) = side_weighted_measure(cell,side_point);
728 jac_det(cell,cell_point) = side_det_jacobian(cell,side_point);
729 for(
int dim=0;dim<cell_dim;++dim){
730 ref_ip_coordinates(cell,cell_point,dim) = cub_points(cell_point,dim);
731 ip_coordinates(cell,cell_point,dim) = side_ip_coordinates(cell,side_point,dim);
732 surface_normals(cell,cell_point,dim) = side_normals(cell,side_point,dim);
734 for(
int dim2=0;dim2<cell_dim;++dim2){
735 jac(cell,cell_point,dim,dim2) = side_jacobian(cell,side_point,dim,dim2);
736 jac_inv(cell,cell_point,dim,dim2) = side_inverse_jacobian(cell,side_point,dim,dim2);
741 point_offset += num_points_on_face;
746 for(
int subcell_index=0; subcell_index<num_subcells;++subcell_index){
749 const int num_points_on_face = ir.
getPointOffset(subcell_index+1) - point_offset;
750 std::vector<int> point_indexes(num_points_on_face,-1);
752 for(
int cell=0; cell<num_cells; ++cell){
755 uniqueCoordOrdering(ip_coordinates,cell,point_offset,point_indexes);
758 reorder(point_indexes,[=](
int a,
int b) { swapQuadraturePoints(cell,point_offset+a,point_offset+b); });
766 Scalar transverse[3];
768 for(
int i=0;i<3;i++){normal[i]=0.;}
769 for(
int i=0;i<3;i++){transverse[i]=0.;}
770 for(
int i=0;i<3;i++){binormal[i]=0.;}
771 for(
int cell=0; cell<num_cells; ++cell){
772 for(
int subcell_index=0; subcell_index<num_subcells; ++subcell_index){
773 for(
int point=0; point<num_points; ++point){
775 for(
int dim=0; dim<3; ++dim)
778 for(
int dim=0; dim<cell_dim; ++dim){
779 normal[dim] = surface_normals(cell,point,dim);
782 convertNormalToRotationMatrix<Scalar>(normal,transverse,binormal);
784 for(
int dim=0; dim<3; ++dim){
785 surface_rotation_matrices(cell,point,0,dim) = normal[dim];
786 surface_rotation_matrices(cell,point,1,dim) = transverse[dim];
787 surface_rotation_matrices(cell,point,2,dim) = binormal[dim];
797 for (
size_type cell = 0; cell < contravarient.extent(0); ++cell) {
798 for (
size_type ip = 0; ip < contravarient.extent(1); ++ip) {
801 for (
size_type i = 0; i < contravarient.extent(2); ++i)
802 for (
size_type j = 0; j < contravarient.extent(3); ++j)
803 covarient(cell,ip,i,j) = 0.0;
806 for (
size_type i = 0; i < contravarient.extent(2); ++i) {
807 for (
size_type j = 0; j < contravarient.extent(3); ++j) {
808 for (
size_type alpha = 0; alpha < contravarient.extent(2); ++alpha) {
809 covarient(cell,ip,i,j) +=
jac(cell,ip,i,alpha) *
jac(cell,ip,j,alpha);
817 Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(contravarient.get_view(), covarient.get_view());
820 for (
size_type cell = 0; cell < contravarient.extent(0); ++cell) {
821 for (
size_type ip = 0; ip < contravarient.extent(1); ++ip) {
822 norm_contravarient(cell,ip) = 0.0;
823 for (
size_type i = 0; i < contravarient.extent(2); ++i) {
824 for (
size_type j = 0; j < contravarient.extent(3); ++j) {
825 norm_contravarient(cell,ip) += contravarient(cell,ip,i,j) * contravarient(cell,ip,i,j);
828 norm_contravarient(cell,ip) = std::sqrt(norm_contravarient(cell,ip));
835 template <
typename Scalar>
839 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
843 size_type num_ip = dyn_cub_points.extent(0);
844 size_type num_dims = dyn_cub_points.extent(1);
846 for (
size_type ip = 0; ip < num_ip; ++ip) {
847 cub_weights(ip) = dyn_cub_weights(ip);
848 for (
size_type dim = 0; dim < num_dims; ++dim)
849 cub_points(ip,dim) = dyn_cub_points(ip,dim);
853 if (int_rule->isSide()) {
854 const size_type num_ip = dyn_cub_points.extent(0), num_side_dims = dyn_side_cub_points.extent(1);
855 for (
size_type ip = 0; ip < num_ip; ++ip)
856 for (
size_type dim = 0; dim < num_side_dims; ++dim)
857 side_cub_points(ip,dim) = dyn_side_cub_points(ip,dim);
861 size_type num_cells = in_node_coordinates.extent(0);
862 size_type num_nodes = in_node_coordinates.extent(1);
863 size_type num_dims = in_node_coordinates.extent(2);
865 for (
size_type cell = 0; cell < num_cells; ++cell) {
866 for (
size_type node = 0; node < num_nodes; ++node) {
867 for (
size_type dim = 0; dim < num_dims; ++dim) {
868 node_coordinates(cell,node,dim) =
869 in_node_coordinates(cell,node,dim);
875 cell_tools.setJacobian(
jac.get_view(),
876 cub_points.get_view(),
877 node_coordinates.get_view(),
878 *(int_rule->topology));
880 cell_tools.setJacobianInv(jac_inv.get_view(),
jac.get_view());
882 cell_tools.setJacobianDet(jac_det.get_view(),
jac.get_view());
884 if (!int_rule->isSide()) {
885 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
886 computeCellMeasure(weighted_measure.get_view(), jac_det.get_view(), cub_weights.get_view());
888 else if(int_rule->spatial_dimension==3) {
889 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
890 computeFaceMeasure(weighted_measure.get_view(),
jac.get_view(), cub_weights.get_view(),
891 int_rule->side, *int_rule->topology,
892 scratch_for_compute_side_measure.get_view());
894 else if(int_rule->spatial_dimension==2) {
895 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
896 computeEdgeMeasure(weighted_measure.get_view(),
jac.get_view(), cub_weights.get_view(),
897 int_rule->side,*int_rule->topology,
898 scratch_for_compute_side_measure.get_view());
903 for (
size_type cell = 0; cell < contravarient.extent(0); ++cell) {
904 for (
size_type ip = 0; ip < contravarient.extent(1); ++ip) {
907 for (
size_type i = 0; i < contravarient.extent(2); ++i)
908 for (
size_type j = 0; j < contravarient.extent(3); ++j)
909 covarient(cell,ip,i,j) = 0.0;
912 for (
size_type i = 0; i < contravarient.extent(2); ++i) {
913 for (
size_type j = 0; j < contravarient.extent(3); ++j) {
914 for (
size_type alpha = 0; alpha < contravarient.extent(2); ++alpha) {
915 covarient(cell,ip,i,j) +=
jac(cell,ip,i,alpha) *
jac(cell,ip,j,alpha);
923 Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(contravarient.get_view(), covarient.get_view());
926 for (
size_type cell = 0; cell < contravarient.extent(0); ++cell) {
927 for (
size_type ip = 0; ip < contravarient.extent(1); ++ip) {
928 norm_contravarient(cell,ip) = 0.0;
929 for (
size_type i = 0; i < contravarient.extent(2); ++i) {
930 for (
size_type j = 0; j < contravarient.extent(3); ++j) {
931 norm_contravarient(cell,ip) += contravarient(cell,ip,i,j) * contravarient(cell,ip,i,j);
934 norm_contravarient(cell,ip) = std::sqrt(norm_contravarient(cell,ip));
942 template <
typename Scalar>
945 const PHX::MDField<Scalar,Cell,IP,Dim>& other_coords,
946 std::vector<
typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type>& permutation)
952 const size_type cell = 0;
953 const size_type num_ip = coords.extent(1), num_dim = coords.extent(2);
954 permutation.resize(num_ip);
955 std::vector<char> taken(num_ip, 0);
956 for (size_type ip = 0; ip < num_ip; ++ip) {
960 for (size_type other_ip = 0; other_ip < num_ip; ++other_ip) {
962 if (taken[other_ip])
continue;
965 for (size_type dim = 0; dim < num_dim; ++dim) {
966 const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim);
969 if (d_min < 0 || d < d_min) {
975 permutation[ip] = i_min;
981 template <
typename Scalar>
984 const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates)
986 if (int_rule->cv_type ==
"none") {
988 getCubature(in_node_coordinates);
992 std::vector<size_type> permutation(other_ip_coordinates.extent(1));
993 permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
996 const size_type num_ip = dyn_cub_points.extent(0);
998 const size_type num_dim = dyn_side_cub_points.extent(1);
999 DblArrayDynamic old_dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
1000 "old_dyn_side_cub_points", num_ip, num_dim);
1001 old_dyn_side_cub_points.deep_copy(dyn_side_cub_points);
1002 for (
size_type ip = 0; ip < num_ip; ++ip)
1003 if (ip != permutation[ip])
1004 for (
size_type dim = 0; dim < num_dim; ++dim)
1005 dyn_side_cub_points(ip, dim) = old_dyn_side_cub_points(permutation[ip], dim);
1008 const size_type num_dim = dyn_cub_points.extent(1);
1009 DblArrayDynamic old_dyn_cub_points = af.template buildArray<double,IP,Dim>(
1010 "old_dyn_cub_points", num_ip, num_dim);
1011 old_dyn_cub_points.deep_copy(dyn_cub_points);
1012 for (
size_type ip = 0; ip < num_ip; ++ip)
1013 if (ip != permutation[ip])
1014 for (
size_type dim = 0; dim < num_dim; ++dim)
1015 dyn_cub_points(ip, dim) = old_dyn_cub_points(permutation[ip], dim);
1018 DblArrayDynamic old_dyn_cub_weights = af.template buildArray<double,IP>(
1019 "old_dyn_cub_weights", num_ip);
1020 old_dyn_cub_weights.deep_copy(dyn_cub_weights);
1021 for (
size_type ip = 0; ip < dyn_cub_weights.extent(0); ++ip)
1022 if (ip != permutation[ip])
1023 dyn_cub_weights(ip) = old_dyn_cub_weights(permutation[ip]);
1026 const size_type num_cells = ip_coordinates.extent(0), num_ip = ip_coordinates.extent(1),
1027 num_dim = ip_coordinates.extent(2);
1028 Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1029 "old_ip_coordinates", num_cells, num_ip, num_dim);
1030 Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
1031 for (
size_type cell = 0; cell < num_cells; ++cell)
1032 for (
size_type ip = 0; ip < num_ip; ++ip)
1033 if (ip != permutation[ip])
1034 for (
size_type dim = 0; dim < num_dim; ++dim)
1035 ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
1040 evaluateRemainingValues(in_node_coordinates);
1045 getCubatureCV(in_node_coordinates);
1048 std::vector<size_type> permutation(other_ip_coordinates.extent(1));
1049 permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
1054 const size_type num_cells = ip_coordinates.extent(0), num_ip = ip_coordinates.extent(1),
1055 num_dim = ip_coordinates.extent(2);
1056 Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1057 "old_ip_coordinates", num_cells, num_ip, num_dim);
1058 Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
1059 Array_CellIPDim old_weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
1060 "old_weighted_normals", num_cells, num_ip, num_dim);
1061 Array_CellIP old_weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
1062 "old_weighted_measure", num_cells, num_ip);
1063 if (int_rule->cv_type ==
"side")
1064 Kokkos::deep_copy(old_weighted_normals.get_static_view(), weighted_normals.get_static_view());
1066 Kokkos::deep_copy(old_weighted_measure.get_static_view(), weighted_measure.get_static_view());
1067 for (
size_type cell = 0; cell < num_cells; ++cell)
1069 for (
size_type ip = 0; ip < num_ip; ++ip)
1071 if (ip != permutation[ip]) {
1072 if (int_rule->cv_type ==
"boundary" || int_rule->cv_type ==
"volume")
1073 weighted_measure(cell, ip) = old_weighted_measure(cell, permutation[ip]);
1074 for (
size_type dim = 0; dim < num_dim; ++dim)
1076 ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
1077 if (int_rule->cv_type ==
"side")
1078 weighted_normals(cell, ip, dim) = old_weighted_normals(cell, permutation[ip], dim);
1086 evaluateValuesCV(in_node_coordinates);
1090 template <
typename Scalar>
1094 int num_space_dim = int_rule->topology->getDimension();
1095 if (int_rule->isSide() && num_space_dim==1) {
1096 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " 1097 <<
"non-natural integration rules.";
1101 size_type num_cells = in_node_coordinates.extent(0);
1102 size_type num_nodes = in_node_coordinates.extent(1);
1103 size_type num_dims = in_node_coordinates.extent(2);
1105 for (
size_type cell = 0; cell < num_cells; ++cell) {
1106 for (
size_type node = 0; node < num_nodes; ++node) {
1107 for (
size_type dim = 0; dim < num_dims; ++dim) {
1108 node_coordinates(cell,node,dim) =
1109 in_node_coordinates(cell,node,dim);
1110 dyn_node_coordinates(cell,node,dim) =
1111 Sacado::ScalarValue<Scalar>::eval(in_node_coordinates(cell,node,dim));
1117 if (int_rule->cv_type ==
"side")
1118 intrepid_cubature->getCubature(dyn_phys_cub_points.get_view(),dyn_phys_cub_norms.get_view(),dyn_node_coordinates.get_view());
1120 intrepid_cubature->getCubature(dyn_phys_cub_points.get_view(),dyn_phys_cub_weights.get_view(),dyn_node_coordinates.get_view());
1122 size_type num_cells = dyn_phys_cub_points.extent(0);
1123 size_type num_ip =dyn_phys_cub_points.extent(1);
1124 size_type num_dims = dyn_phys_cub_points.extent(2);
1126 for (
size_type cell = 0; cell < num_cells; ++cell) {
1127 for (
size_type ip = 0; ip < num_ip; ++ip) {
1128 if (int_rule->cv_type !=
"side")
1129 weighted_measure(cell,ip) = dyn_phys_cub_weights(cell,ip);
1130 for (
size_type dim = 0; dim < num_dims; ++dim) {
1131 ip_coordinates(cell,ip,dim) = dyn_phys_cub_points(cell,ip,dim);
1132 if (int_rule->cv_type ==
"side")
1133 weighted_normals(cell,ip,dim) = dyn_phys_cub_norms(cell,ip,dim);
1140 template <
typename Scalar>
1145 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1147 cell_tools.mapToReferenceFrame(ref_ip_coordinates.get_view(),
1148 ip_coordinates.get_view(),
1149 node_coordinates.get_view(),
1150 *(int_rule->topology));
1152 cell_tools.setJacobian(
jac.get_view(),
1153 ref_ip_coordinates.get_view(),
1154 node_coordinates.get_view(),
1155 *(int_rule->topology));
1157 cell_tools.setJacobianInv(jac_inv.get_view(),
jac.get_view());
1159 cell_tools.setJacobianDet(jac_det.get_view(),
jac.get_view());
1163 #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \ 1164 template class IntegrationValues2<SCALAR>;
void generateSurfaceCubatureValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates)
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates)
Cell vertex coordinates, not basis coordinates.
static void uniqueCoordOrdering(Array_CellIPDim &coords, int cell, int offset, std::vector< int > &order)
Using coordinate build an arrray that specifies a unique ordering.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
PHX::MDField< Scalar, Cell, IP > Array_CellIP
const int & getType() const
Get type of integrator.
void evaluateRemainingValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates)
void swapQuadraturePoints(int cell, int a, int b)
Swap the ordering of quadrature points in a specified cell.
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int getPointOffset(const int subcell_index) const
Returns the integration point offset for a given subcell_index (i.e. local face index) ...
Teuchos::RCP< Intrepid2::Cubature< PHX::Device::execution_space, double, double > > getIntrepidCubature(const panzer::IntegrationRule &ir) const
Teuchos::RCP< const shards::CellTopology > topology
ArrayTraits< Scalar, PHX::MDField< Scalar > >::size_type size_type
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
PHX::MDField< Scalar, Cell, IP, Dim > Array_CellIPDim
void reorder(std::vector< int > &order, std::function< void(int, int)> swapper)
Using a functor, reorder an array using a order vector.
PHX::MDField< double > DblArrayDynamic
void getCubature(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates)
Teuchos::RCP< shards::CellTopology > side_topology
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates)
#define TEUCHOS_ASSERT(assertion_test)
static void permuteToOther(const PHX::MDField< Scalar, Cell, IP, Dim > &coords, const PHX::MDField< Scalar, Cell, IP, Dim > &other_coords, std::vector< typename ArrayTraits< Scalar, PHX::MDField< Scalar > >::size_type > &permutation)
void getCubatureCV(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates)
void setupArraysForNodeRule(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
const int & getOrder() const
Get order of integrator.