45 #include "Shards_CellTopology.hpp" 47 #include "Kokkos_DynRankView.hpp" 48 #include "Intrepid2_FunctionSpaceTools.hpp" 49 #include "Intrepid2_RealSpaceTools.hpp" 50 #include "Intrepid2_CellTools.hpp" 51 #include "Intrepid2_ArrayTools.hpp" 52 #include "Intrepid2_CubatureControlVolume.hpp" 53 #include "Intrepid2_CubatureControlVolumeSide.hpp" 54 #include "Intrepid2_CubatureControlVolumeBoundary.hpp" 66 template <
typename Scalar>
74 int num_space_dim = ir->
topology->getDimension();
78 dyn_cub_points = af.template buildArray<double,IP,Dim>(
"cub_points",
num_ip, num_space_dim);
79 dyn_cub_weights = af.template buildArray<double,IP>(
"cub_weights",
num_ip);
81 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",
num_ip, num_space_dim);
84 dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
"side_cub_points",
num_ip, ir->
side_topology->getDimension());
85 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",
num_ip,ir->
side_topology->getDimension());
89 dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_points",num_cells,
num_ip, num_space_dim);
90 dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>(
"phys_cub_weights",num_cells,
num_ip);
92 dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_norms",num_cells,
num_ip, num_space_dim);
96 dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>(
"node_coordinates",num_cells,
num_nodes, num_space_dim);
98 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",
num_ip);
100 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells,
num_nodes, num_space_dim);
102 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells,
num_ip, num_space_dim,num_space_dim);
104 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells,
num_ip, num_space_dim,num_space_dim);
106 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells,
num_ip);
108 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells,
num_ip);
110 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells,
num_ip, num_space_dim,num_space_dim);
112 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells,
num_ip, num_space_dim,num_space_dim);
114 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells,
num_ip);
116 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells,
num_ip,num_space_dim);
118 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells,
num_ip,num_space_dim);
120 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted normal",num_cells,
num_ip,num_space_dim);
124 template <
typename Scalar>
134 int num_space_dim = ir->
topology->getDimension();
137 if(num_space_dim==1 && ir->
isSide()) {
138 setupArraysForNodeRule(ir);
142 Intrepid2::DefaultCubatureFactory<double,DblArrayDynamic>
146 intrepid_cubature =
Teuchos::rcp(
new Intrepid2::CubatureControlVolumeSide<double,DblArrayDynamic,DblArrayDynamic>(ir->
topology));
148 else if (ir->
cv_type ==
"volume")
149 intrepid_cubature =
Teuchos::rcp(
new Intrepid2::CubatureControlVolume<double,DblArrayDynamic,DblArrayDynamic>(ir->
topology));
152 intrepid_cubature =
Teuchos::rcp(
new Intrepid2::CubatureControlVolumeBoundary<double,DblArrayDynamic,DblArrayDynamic>(ir->
topology,ir->
side));
155 intrepid_cubature = cubature_factory.create(*(ir->
side_topology),
158 intrepid_cubature = cubature_factory.create(*(ir->
topology),
162 int num_ip = intrepid_cubature->getNumPoints();
164 dyn_cub_points = af.template buildArray<double,IP,Dim>(
"cub_points",
num_ip, num_space_dim);
165 dyn_cub_weights = af.template buildArray<double,IP>(
"cub_weights",
num_ip);
167 cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"cub_points",
num_ip, num_space_dim);
170 dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
"side_cub_points",
num_ip, ir->
side_topology->getDimension());
171 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>(
"side_cub_points",
num_ip,ir->
side_topology->getDimension());
175 dyn_phys_cub_points = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_points",num_cells,
num_ip, num_space_dim);
176 dyn_phys_cub_weights = af.template buildArray<double,Cell,IP>(
"phys_cub_weights",num_cells,
num_ip);
178 dyn_phys_cub_norms = af.template buildArray<double,Cell,IP,Dim>(
"phys_cub_norms",num_cells,
num_ip, num_space_dim);
182 dyn_node_coordinates = af.template buildArray<double,Cell,BASIS,Dim>(
"node_coordinates",num_cells,
num_nodes, num_space_dim);
184 cub_weights = af.template buildStaticArray<Scalar,IP>(
"cub_weights",
num_ip);
186 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>(
"node_coordinates",num_cells,
num_nodes, num_space_dim);
188 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac",num_cells,
num_ip, num_space_dim,num_space_dim);
190 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"jac_inv",num_cells,
num_ip, num_space_dim,num_space_dim);
192 jac_det = af.template buildStaticArray<Scalar,Cell,IP>(
"jac_det",num_cells,
num_ip);
194 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
"weighted_measure",num_cells,
num_ip);
196 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"covarient",num_cells,
num_ip, num_space_dim,num_space_dim);
198 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>(
"contravarient",num_cells,
num_ip, num_space_dim,num_space_dim);
200 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>(
"norm_contravarient",num_cells,
num_ip);
202 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ip_coordiantes",num_cells,
num_ip,num_space_dim);
204 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"ref_ip_coordinates",num_cells,
num_ip,num_space_dim);
206 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
"weighted_normal",num_cells,
num_ip,num_space_dim);
213 template <
typename Scalar>
217 if (int_rule->cv_type !=
"none") {
218 getCubatureCV(in_node_coordinates);
219 evaluateValuesCV(in_node_coordinates);
222 getCubature(in_node_coordinates);
223 evaluateRemainingValues(in_node_coordinates);
227 template <
typename Scalar>
229 getCubature(
const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates)
231 int num_space_dim = int_rule->topology->getDimension();
232 if (int_rule->isSide() && num_space_dim==1) {
233 std::cout <<
"WARNING: 0-D quadrature rule ifrastructure does not exist!!! Will not be able to do " 234 <<
"non-natural integration rules.";
238 Intrepid2::CellTools<Scalar> cell_tools;
240 if (!int_rule->isSide())
241 intrepid_cubature->getCubature(dyn_cub_points, dyn_cub_weights);
243 intrepid_cubature->getCubature(dyn_side_cub_points, dyn_cub_weights);
245 cell_tools.mapToReferenceSubcell(dyn_cub_points,
247 int_rule->spatial_dimension-1,
249 *(int_rule->topology));
253 cell_tools.mapToPhysicalFrame(ip_coordinates, dyn_cub_points, in_node_coordinates, *(int_rule->topology));
257 template <
typename Scalar>
261 Intrepid2::CellTools<Scalar> cell_tools;
266 size_type num_dims = dyn_cub_points.dimension(1);
269 cub_weights(ip) = dyn_cub_weights(ip);
270 for (
size_type dim = 0; dim < num_dims; ++dim)
271 cub_points(ip,dim) = dyn_cub_points(ip,dim);
275 if (int_rule->isSide()) {
276 const size_type num_ip = dyn_cub_points.dimension(0), num_side_dims = dyn_side_cub_points.dimension(1);
278 for (
size_type dim = 0; dim < num_side_dims; ++dim)
279 side_cub_points(ip,dim) = dyn_side_cub_points(ip,dim);
283 size_type num_cells = in_node_coordinates.dimension(0);
285 size_type num_dims = in_node_coordinates.dimension(2);
287 for (
size_type cell = 0; cell < num_cells; ++cell) {
289 for (
size_type dim = 0; dim < num_dims; ++dim) {
290 node_coordinates(cell,node,dim) =
291 in_node_coordinates(cell,node,dim);
297 cell_tools.setJacobian(
jac, cub_points, node_coordinates,
298 *(int_rule->topology));
300 cell_tools.setJacobianInv(jac_inv,
jac);
302 cell_tools.setJacobianDet(jac_det,
jac);
304 if (!int_rule->isSide()) {
305 Intrepid2::FunctionSpaceTools::
306 computeCellMeasure<Scalar>(weighted_measure, jac_det, cub_weights);
308 else if(int_rule->spatial_dimension==3) {
309 Intrepid2::FunctionSpaceTools::
310 computeFaceMeasure<Scalar>(weighted_measure,
jac, cub_weights,int_rule->side,*int_rule->topology);
312 else if(int_rule->spatial_dimension==2) {
313 Intrepid2::FunctionSpaceTools::
314 computeEdgeMeasure<Scalar>(weighted_measure,
jac, cub_weights,int_rule->side,*int_rule->topology);
319 for (
size_type cell = 0; cell < contravarient.dimension(0); ++cell) {
320 for (
size_type ip = 0; ip < contravarient.dimension(1); ++ip) {
323 for (
size_type i = 0; i < contravarient.dimension(2); ++i)
324 for (
size_type j = 0; j < contravarient.dimension(3); ++j)
325 covarient(cell,ip,i,j) = 0.0;
328 for (
size_type i = 0; i < contravarient.dimension(2); ++i) {
329 for (
size_type j = 0; j < contravarient.dimension(3); ++j) {
330 for (
size_type alpha = 0; alpha < contravarient.dimension(2); ++alpha) {
331 covarient(cell,ip,i,j) +=
jac(cell,ip,i,alpha) *
jac(cell,ip,j,alpha);
339 Intrepid2::RealSpaceTools<Scalar>::inverse(contravarient, covarient);
342 for (
size_type cell = 0; cell < contravarient.dimension(0); ++cell) {
343 for (
size_type ip = 0; ip < contravarient.dimension(1); ++ip) {
344 norm_contravarient(cell,ip) = 0.0;
345 for (
size_type i = 0; i < contravarient.dimension(2); ++i) {
346 for (
size_type j = 0; j < contravarient.dimension(3); ++j) {
347 norm_contravarient(cell,ip) += contravarient(cell,ip,i,j) * contravarient(cell,ip,i,j);
350 norm_contravarient(cell,ip) = std::sqrt(norm_contravarient(cell,ip));
358 template <
typename Scalar>
361 const PHX::MDField<Scalar,Cell,IP,Dim>& other_coords,
368 const size_type cell = 0;
369 const size_type
num_ip = coords.dimension(1),
num_dim = coords.dimension(2);
370 permutation.resize(
num_ip);
371 std::vector<char> taken(
num_ip, 0);
372 for (size_type ip = 0; ip <
num_ip; ++ip) {
376 for (size_type other_ip = 0; other_ip <
num_ip; ++other_ip) {
378 if (taken[other_ip])
continue;
381 for (size_type dim = 0; dim <
num_dim; ++dim) {
382 const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim);
385 if (d_min < 0 || d < d_min) {
391 permutation[ip] = i_min;
397 template <
typename Scalar>
400 const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates)
402 if (int_rule->cv_type ==
"none") {
404 getCubature(in_node_coordinates);
408 std::vector<size_type> permutation(other_ip_coordinates.dimension(1));
409 permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
415 DblArrayDynamic old_dyn_side_cub_points = af.template buildArray<double,IP,Dim>(
417 old_dyn_side_cub_points.deep_copy(dyn_side_cub_points);
419 if (ip != permutation[ip])
421 dyn_side_cub_points(ip, dim) = old_dyn_side_cub_points(permutation[ip], dim);
425 DblArrayDynamic old_dyn_cub_points = af.template buildArray<double,IP,Dim>(
427 old_dyn_cub_points.deep_copy(dyn_cub_points);
429 if (ip != permutation[ip])
431 dyn_cub_points(ip, dim) = old_dyn_cub_points(permutation[ip], dim);
434 DblArrayDynamic old_dyn_cub_weights = af.template buildArray<double,IP>(
435 "old_dyn_cub_weights",
num_ip);
436 old_dyn_cub_weights.deep_copy(dyn_cub_weights);
437 for (
size_type ip = 0; ip < dyn_cub_weights.dimension(0); ++ip)
438 if (ip != permutation[ip])
439 dyn_cub_weights(ip) = old_dyn_cub_weights(permutation[ip]);
442 const size_type num_cells = ip_coordinates.dimension(0),
num_ip = ip_coordinates.dimension(1),
443 num_dim = ip_coordinates.dimension(2);
444 Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
446 Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
447 for (
size_type cell = 0; cell < num_cells; ++cell)
449 if (ip != permutation[ip])
451 ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
456 evaluateRemainingValues(in_node_coordinates);
461 getCubatureCV(in_node_coordinates);
464 std::vector<size_type> permutation(other_ip_coordinates.dimension(1));
465 permuteToOther(ip_coordinates, other_ip_coordinates, permutation);
471 const size_type num_cells = ip_coordinates.dimension(0),
num_ip = ip_coordinates.dimension(1),
472 num_dim = ip_coordinates.dimension(2);
473 Array_CellIPDim old_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
475 Kokkos::deep_copy(old_ip_coordinates.get_static_view(), ip_coordinates.get_static_view());
476 Array_CellIPDim old_weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>(
478 Array_CellIP old_weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>(
479 "old_weighted_measure", num_cells,
num_ip);
480 if (int_rule->cv_type ==
"side")
481 Kokkos::deep_copy(old_weighted_normals.get_static_view(), weighted_normals.get_static_view());
483 Kokkos::deep_copy(old_weighted_measure.get_static_view(), weighted_measure.get_static_view());
484 for (
size_type cell = 0; cell < num_cells; ++cell)
488 if (ip != permutation[ip]) {
489 if (int_rule->cv_type ==
"boundary" || int_rule->cv_type ==
"volume")
490 weighted_measure(cell, ip) = old_weighted_measure(cell, permutation[ip]);
493 ip_coordinates(cell, ip, dim) = old_ip_coordinates(cell, permutation[ip], dim);
494 if (int_rule->cv_type ==
"side")
495 weighted_normals(cell, ip, dim) = old_weighted_normals(cell, permutation[ip], dim);
503 evaluateValuesCV(in_node_coordinates);
507 template <
typename Scalar>
509 getCubatureCV(
const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates)
511 int num_space_dim = int_rule->topology->getDimension();
512 if (int_rule->isSide() && num_space_dim==1) {
513 std::cout <<
"WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do " 514 <<
"non-natural integration rules.";
518 size_type num_cells = in_node_coordinates.dimension(0);
520 size_type num_dims = in_node_coordinates.dimension(2);
522 for (
size_type cell = 0; cell < num_cells; ++cell) {
524 for (
size_type dim = 0; dim < num_dims; ++dim) {
525 node_coordinates(cell,node,dim) =
526 in_node_coordinates(cell,node,dim);
527 dyn_node_coordinates(cell,node,dim) =
528 Sacado::ScalarValue<Scalar>::eval(in_node_coordinates(cell,node,dim));
534 if (int_rule->cv_type ==
"side")
535 intrepid_cubature->getCubature(dyn_phys_cub_points,dyn_phys_cub_norms,dyn_node_coordinates);
537 intrepid_cubature->getCubature(dyn_phys_cub_points,dyn_phys_cub_weights,dyn_node_coordinates);
539 size_type num_cells = dyn_phys_cub_points.dimension(0);
541 size_type num_dims = dyn_phys_cub_points.dimension(2);
543 for (
size_type cell = 0; cell < num_cells; ++cell) {
545 if (int_rule->cv_type !=
"side")
546 weighted_measure(cell,ip) = dyn_phys_cub_weights(cell,ip);
547 for (
size_type dim = 0; dim < num_dims; ++dim) {
548 ip_coordinates(cell,ip,dim) = dyn_phys_cub_points(cell,ip,dim);
549 if (int_rule->cv_type ==
"side")
550 weighted_normals(cell,ip,dim) = dyn_phys_cub_norms(cell,ip,dim);
557 template <
typename Scalar>
562 Intrepid2::CellTools<Scalar> cell_tools;
564 cell_tools.mapToReferenceFrame(ref_ip_coordinates, ip_coordinates, node_coordinates,
565 *(int_rule->topology),-1);
567 cell_tools.setJacobian(
jac, ref_ip_coordinates, node_coordinates,
568 *(int_rule->topology));
570 cell_tools.setJacobianInv(jac_inv,
jac);
572 cell_tools.setJacobianDet(jac_det,
jac);
576 #define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \ 577 template class IntegrationValues2<SCALAR>;
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates)
Cell vertex coordinates, not basis coordinates.
PHX::MDField< Scalar, Cell, IP > Array_CellIP
int side
Defaults to -1 if this is volume and not sideset.
void evaluateRemainingValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &in_node_coordinates)
PHX::MDField< ScalarT > vector
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
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)