43 #include "PanzerDiscFE_config.hpp" 49 #include "Intrepid2_FunctionSpaceTools.hpp" 54 template <
typename Scalar>
56 evaluateValues(
const PHX::MDField<Scalar,IP,Dim,void,void,void,void,void,void> & cub_points,
57 const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> &
jac,
58 const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
59 const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv)
61 PHX::MDField<Scalar,Cell,IP> weighted_measure;
62 PHX::MDField<Scalar,Cell,NODE,Dim> vertex_coordinates;
63 build_weighted =
false;
64 evaluateValues(cub_points,
jac,jac_det,jac_inv,weighted_measure,vertex_coordinates,
false);
67 template <
typename Scalar>
69 evaluateValues(
const PHX::MDField<Scalar,IP,Dim,void,void,void,void,void,void> & cub_points,
70 const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> &
jac,
71 const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
72 const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv,
73 const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
74 const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
75 bool use_vertex_coordinates)
79 int num_dim = basis_layout->dimension();
83 evaluateReferenceValues(cub_points,compute_derivatives,use_vertex_coordinates);
88 Intrepid2::FunctionSpaceTools::
89 HGRADtransformVALUE<Scalar>(basis_scalar,
93 Intrepid2::FunctionSpaceTools::
94 multiplyMeasure<Scalar>(weighted_basis_scalar,
100 Intrepid2::FunctionSpaceTools::
101 HCURLtransformVALUE<Scalar>(basis_vector,
106 Intrepid2::FunctionSpaceTools::
114 Intrepid2::FunctionSpaceTools::
115 HDIVtransformVALUE<Scalar>(basis_vector,
121 Intrepid2::FunctionSpaceTools::
130 Intrepid2::FunctionSpaceTools::
131 HGRADtransformGRAD<Scalar>(grad_basis,
136 Intrepid2::FunctionSpaceTools::
137 multiplyMeasure<Scalar>(weighted_grad_basis,
143 Intrepid2::FunctionSpaceTools::
144 HDIVtransformDIV<Scalar>(curl_basis_scalar,
148 curl_basis_ref_scalar);
151 Intrepid2::FunctionSpaceTools::
152 multiplyMeasure<Scalar>(weighted_curl_basis_scalar,
158 Intrepid2::FunctionSpaceTools::
159 HCURLtransformCURL<Scalar>(curl_basis_vector,
162 curl_basis_ref_vector);
165 Intrepid2::FunctionSpaceTools::
166 multiplyMeasure<Scalar>(weighted_curl_basis_vector,
172 Intrepid2::FunctionSpaceTools::
173 HDIVtransformDIV<Scalar>(div_basis,
178 Intrepid2::FunctionSpaceTools::
179 multiplyMeasure<Scalar>(weighted_div_basis,
187 if(use_vertex_coordinates) {
189 = Teuchos::rcp_dynamic_cast<Intrepid2::DofCoordsInterface<ArrayDynamic> >(intrepid_basis);
201 Intrepid2::CellTools<Scalar> cell_tools;
202 cell_tools.mapToPhysicalFrame(basis_coordinates,
203 basis_coordinates_ref,
205 intrepid_basis->getBaseCellTopology());
210 template <
typename Scalar>
212 evaluateValuesCV(
const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cell_cub_points,
213 const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> &
jac,
214 const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
215 const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv)
219 int num_ip = basis_layout->numPoints();
220 int num_card = basis_layout->cardinality();
221 int num_dim = basis_layout->dimension();
230 for (
size_type icell = 0; icell < num_cells; ++icell)
232 for (
int ip = 0; ip <
num_ip; ++ip)
233 for (
int d = 0; d <
num_dim; ++d)
234 dyn_cub_points(ip,d) = cell_cub_points(icell,ip,d);
239 intrepid_basis->getValues(dyn_basis_ref_scalar, dyn_cub_points,
240 Intrepid2::OPERATOR_VALUE);
243 for (
int b = 0; b < num_card; ++b)
244 for (
int ip = 0; ip <
num_ip; ++ip)
245 basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip);
251 intrepid_basis->getValues(dyn_basis_ref_scalar, dyn_cub_points,
252 Intrepid2::OPERATOR_VALUE);
255 for (
int b = 0; b < num_card; ++b)
256 for (
int ip = 0; ip <
num_ip; ++ip)
257 basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip);
259 if(compute_derivatives) {
266 intrepid_basis->getValues(dyn_grad_basis_ref, dyn_cub_points,
267 Intrepid2::OPERATOR_GRAD);
270 for (
int ip = 0; ip <
num_ip; ++ip)
271 for (
int d1 = 0; d1 <
num_dim; ++d1)
272 for (
int d2 = 0; d2 <
num_dim; ++d2)
273 dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2);
275 Intrepid2::FunctionSpaceTools::HGRADtransformGRAD<Scalar>(dyn_grad_basis,
279 for (
int b = 0; b < num_card; ++b)
280 for (
int ip = 0; ip <
num_ip; ++ip)
281 for (
int d = 0; d <
num_dim; ++d)
282 grad_basis(icell,b,ip,d) = dyn_grad_basis(0,b,ip,d);
289 intrepid_basis->getValues(dyn_basis_ref_vector, dyn_cub_points,
290 Intrepid2::OPERATOR_VALUE);
297 for (
int ip = 0; ip <
num_ip; ++ip)
298 for (
int d1 = 0; d1 <
num_dim; ++d1)
299 for (
int d2 = 0; d2 <
num_dim; ++d2)
300 dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2);
302 Intrepid2::FunctionSpaceTools::HCURLtransformVALUE<Scalar>(dyn_basis_vector,
304 dyn_basis_ref_vector);
306 for (
int b = 0; b < num_card; ++b)
307 for (
int ip = 0; ip <
num_ip; ++ip)
308 for (
int d = 0; d <
num_dim; ++d)
309 basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d);
311 if(compute_derivatives &&
num_dim ==2) {
318 intrepid_basis->getValues(dyn_curl_basis_ref_scalar, dyn_cub_points,
319 Intrepid2::OPERATOR_CURL);
322 for (
int ip = 0; ip <
num_ip; ++ip)
323 dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
325 Intrepid2::FunctionSpaceTools::HDIVtransformDIV<Scalar>(dyn_curl_basis_scalar,
327 dyn_curl_basis_ref_scalar);
329 for (
int b = 0; b < num_card; ++b)
330 for (
int ip = 0; ip <
num_ip; ++ip)
331 curl_basis_scalar(icell,b,ip) = dyn_curl_basis_scalar(0,b,ip);
334 if(compute_derivatives &&
num_dim ==3) {
342 intrepid_basis->getValues(dyn_curl_basis_ref, dyn_cub_points,
343 Intrepid2::OPERATOR_CURL);
346 for (
int ip = 0; ip <
num_ip; ++ip)
348 dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
349 for (
int d1 = 0; d1 <
num_dim; ++d1)
350 for (
int d2 = 0; d2 <
num_dim; ++d2)
351 dyn_jac(cellInd,ip,d1,d2) =
jac(icell,ip,d1,d2);
354 Intrepid2::FunctionSpaceTools::HCURLtransformCURL<Scalar>(dyn_curl_basis,
359 for (
int b = 0; b < num_card; ++b)
360 for (
int ip = 0; ip <
num_ip; ++ip)
361 for (
int d = 0; d <
num_dim; ++d)
362 curl_basis_vector(icell,b,ip,d) = dyn_curl_basis(0,b,ip,d);
371 intrepid_basis->getValues(dyn_basis_ref_vector, dyn_cub_points,
372 Intrepid2::OPERATOR_VALUE);
380 for (
int ip = 0; ip <
num_ip; ++ip)
382 dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
383 for (
int d1 = 0; d1 <
num_dim; ++d1)
384 for (
int d2 = 0; d2 <
num_dim; ++d2)
385 dyn_jac(cellInd,ip,d1,d2) =
jac(icell,ip,d1,d2);
388 Intrepid2::FunctionSpaceTools::HDIVtransformVALUE<Scalar>(dyn_basis_vector,
390 dyn_basis_ref_vector);
392 for (
int b = 0; b < num_card; ++b)
393 for (
int ip = 0; ip <
num_ip; ++ip)
394 for (
int d = 0; d <
num_dim; ++d)
395 basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d);
397 if(compute_derivatives) {
402 intrepid_basis->getValues(dyn_div_basis_ref, dyn_cub_points,
403 Intrepid2::OPERATOR_DIV);
405 Intrepid2::FunctionSpaceTools::HDIVtransformDIV<Scalar>(dyn_div_basis,
409 for (
int b = 0; b < num_card; ++b)
410 for (
int ip = 0; ip <
num_ip; ++ip)
411 div_basis(icell,b,ip) = dyn_div_basis(0,b,ip);
422 template <
typename Scalar>
428 int num_quad = basis_layout->numPoints();
429 int num_dim = basis_layout->dimension();
430 int num_card = basis_layout->cardinality();
434 for (
int ip = 0; ip < num_quad; ++ip)
435 for (
int d = 0; d <
num_dim; ++d)
436 dyn_cub_points(ip,d) = cub_points(ip,d);
440 ArrayDynamic dyn_basis_ref_scalar = af.
buildArray<Scalar,BASIS,IP>(
"dyn_basis_ref_scalar",num_card,num_quad);
442 intrepid_basis->getValues(dyn_basis_ref_scalar, dyn_cub_points,
443 Intrepid2::OPERATOR_VALUE);
445 for (
int b = 0; b < num_card; ++b)
446 for (
int ip = 0; ip < num_quad; ++ip)
447 basis_ref_scalar(b,ip) = dyn_basis_ref_scalar(b,ip);
452 intrepid_basis->getValues(dyn_basis_ref_vector, dyn_cub_points,
453 Intrepid2::OPERATOR_VALUE);
455 for (
int b = 0; b < num_card; ++b)
456 for (
int ip = 0; ip < num_quad; ++ip)
457 for (
int d = 0; d <
num_dim; ++d)
458 basis_ref_vector(b,ip,d) = dyn_basis_ref_vector(b,ip,d);
465 intrepid_basis->getValues(dyn_grad_basis_ref, dyn_cub_points,
466 Intrepid2::OPERATOR_GRAD);
468 for (
int b = 0; b < num_card; ++b)
469 for (
int ip = 0; ip < num_quad; ++ip)
470 for (
int d = 0; d <
num_dim; ++d)
471 grad_basis_ref(b,ip,d) = dyn_grad_basis_ref(b,ip,d);
474 ArrayDynamic dyn_curl_basis_ref = af.
buildArray<Scalar,BASIS,IP>(
"dyn_curl_basis_ref_scalar",num_card,num_quad);
476 intrepid_basis->getValues(dyn_curl_basis_ref, dyn_cub_points,
477 Intrepid2::OPERATOR_CURL);
479 for (
int b = 0; b < num_card; ++b)
480 for (
int ip = 0; ip < num_quad; ++ip)
481 curl_basis_ref_scalar(b,ip) = dyn_curl_basis_ref(b,ip);
486 intrepid_basis->getValues(dyn_curl_basis_ref, dyn_cub_points,
487 Intrepid2::OPERATOR_CURL);
489 for (
int b = 0; b < num_card; ++b)
490 for (
int ip = 0; ip < num_quad; ++ip)
491 for (
int d = 0; d <
num_dim; ++d)
492 curl_basis_ref_vector(b,ip,d) = dyn_curl_basis_ref(b,ip,d);
495 ArrayDynamic dyn_div_basis_ref = af.
buildArray<Scalar,BASIS,IP>(
"dyn_div_basis_ref_scalar",num_card,num_quad);
497 intrepid_basis->getValues(dyn_div_basis_ref, dyn_cub_points,
498 Intrepid2::OPERATOR_DIV);
500 for (
int b = 0; b < num_card; ++b)
501 for (
int ip = 0; ip < num_quad; ++ip)
502 div_basis_ref(b,ip) = dyn_div_basis_ref(b,ip);
506 if(use_vertex_coordinates) {
508 = Teuchos::rcp_dynamic_cast<Intrepid2::DofCoordsInterface<ArrayDynamic> >(intrepid_basis);
512 ArrayDynamic dyn_basis_coordinates_ref = af.
buildArray<Scalar,BASIS,Dim>(
"basis_coordinates_ref",basis_coordinates_ref.dimension(0),basis_coordinates_ref.dimension(1));
513 coords->getDofCoords(dyn_basis_coordinates_ref);
516 for (
int i = 0; i < basis_coordinates_ref.extent_int(0); ++i)
517 for (
int j = 0; j < basis_coordinates_ref.extent_int(1); ++j)
518 basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j);
522 references_evaluated =
true;
526 template <
typename Scalar>
530 int num_cell = orientations.dimension(0);
531 int num_basis = orientations.dimension(1);
532 int num_dim = basis_layout->dimension();
533 int num_ip = basis_layout->numPoints();
541 for (
int c=0; c<num_cell; c++)
542 for (
int b=0; b<num_basis; b++)
543 for (
int p=0; p<
num_ip; p++)
545 basis_vector(c, b, p, d) *= orientations(c, b);
547 if(compute_derivatives) {
549 for (
int c=0; c<num_cell; c++)
550 for (
int b=0; b<num_basis; b++)
551 for (
int p=0; p<
num_ip; p++)
552 curl_basis_scalar(c, b, p) *= orientations(c, b);
559 if(compute_derivatives)
560 Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_curl_basis_scalar,orientations);
568 for (
int c=0; c<num_cell; c++)
569 for (
int b=0; b<num_basis; b++)
570 for (
int p=0; p<
num_ip; p++)
572 basis_vector(c, b, p, d) *= orientations(c, b);
574 if(compute_derivatives) {
576 for (
int c=0; c<num_cell; c++)
577 for (
int b=0; b<num_basis; b++)
578 for (
int p=0; p<
num_ip; p++)
580 curl_basis_vector(c, b, p,d) *= orientations(c, b);
587 if(compute_derivatives)
588 Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_curl_basis_vector,orientations);
595 for (
int c=0; c<num_cell; c++)
596 for (
int b=0; b<num_basis; b++)
597 for (
int p=0; p<
num_ip; p++)
599 basis_vector(c, b, p, d) *= orientations(c, b);
601 if(compute_derivatives) {
604 for (
int c=0; c<num_cell; c++)
605 for (
int b=0; b<num_basis; b++)
606 for (
int p=0; p<
num_ip; p++)
607 div_basis(c, b, p) *= orientations(c, b);
614 if(compute_derivatives)
615 Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_div_basis,orientations);
620 template <
typename Scalar>
622 {
return basis_layout->getBasis()->getElementSpace(); }
624 template <
typename Scalar>
627 bool computeDerivatives)
631 compute_derivatives = computeDerivatives;
632 basis_layout = layout;
639 int numcells = basisDesc->
numCells();
654 basis_ref_scalar = af.
buildStaticArray<Scalar,BASIS,IP>(
"basis_ref",card,num_quad);
655 basis_scalar = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"basis",numcells,card,num_quad);
658 weighted_basis_scalar = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"weighted_basis",numcells,card,num_quad);
663 if(compute_derivatives) {
664 grad_basis_ref = af.
buildStaticArray<Scalar,BASIS,IP,Dim>(
"grad_basis_ref",card,num_quad,dim);
665 grad_basis = af.
buildStaticArray<Scalar,Cell,BASIS,IP,Dim>(
"grad_basis",numcells,card,num_quad,dim);
668 weighted_grad_basis = af.
buildStaticArray<Scalar,Cell,BASIS,IP,Dim>(
"weighted_grad_basis",numcells,card,num_quad,dim);
682 basis_ref_vector = af.
buildStaticArray<Scalar,BASIS,IP,Dim>(
"basis_ref",card,num_quad,dim);
683 basis_vector = af.
buildStaticArray<Scalar,Cell,BASIS,IP,Dim>(
"basis",numcells,card,num_quad,dim);
696 if(compute_derivatives) {
699 curl_basis_ref_scalar = af.
buildStaticArray<Scalar,BASIS,IP>(
"curl_basis_ref",card,num_quad);
700 curl_basis_scalar = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"curl_basis",numcells,card,num_quad);
703 weighted_curl_basis_scalar = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"weighted_curl_basis",numcells,card,num_quad);
706 curl_basis_ref_vector = af.
buildStaticArray<Scalar,BASIS,IP,Dim>(
"curl_basis_ref",card,num_quad,dim);
707 curl_basis_vector = af.
buildStaticArray<Scalar,Cell,BASIS,IP,Dim>(
"curl_basis",numcells,card,num_quad,dim);
710 weighted_curl_basis_vector = af.
buildStaticArray<Scalar,Cell,BASIS,IP,Dim>(
"weighted_curl_basis",numcells,card,num_quad,dim);
721 basis_ref_vector = af.
buildStaticArray<Scalar,BASIS,IP,Dim>(
"basis_ref",card,num_quad,dim);
722 basis_vector = af.
buildStaticArray<Scalar,Cell,BASIS,IP,Dim>(
"basis",numcells,card,num_quad,dim);
740 if(compute_derivatives) {
741 div_basis_ref = af.
buildStaticArray<Scalar,BASIS,IP>(
"div_basis_ref",card,num_quad);
742 div_basis = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"div_basis",numcells,card,num_quad);
745 weighted_div_basis = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"weighted_div_basis",numcells,card,num_quad);
753 basis_ref_scalar = af.
buildStaticArray<Scalar,BASIS,IP>(
"basis_ref",card,num_quad);
754 basis_scalar = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"basis",numcells,card,num_quad);
757 weighted_basis_scalar = af.
buildStaticArray<Scalar,Cell,BASIS,IP>(
"weighted_basis",numcells,card,num_quad);
776 basis_coordinates_ref = af.
buildStaticArray<Scalar,BASIS,Dim>(
"basis_coordinates_ref",card,dim);
777 basis_coordinates = af.
buildStaticArray<Scalar,Cell,BASIS,Dim>(
"basis_coordinates",numcells,card,dim);
782 #define BASIS_VALUES_INSTANTIATION(SCALAR) \ 783 template class BasisValues2<SCALAR>; 787 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
void evaluateValues(const PHX::MDField< Scalar, IP, Dim, void, void, void, void, void, void > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv)
void evaluateReferenceValues(const PHX::MDField< Scalar, IP, Dim > &cub_points, bool compute_derivatives, bool use_vertex_coordinates)
bool is_null(const std::shared_ptr< T > &p)
ArrayTraits< Scalar, PHX::MDField< Scalar, void, void, void, void, void, void, void, void > >::size_type size_type
Teuchos::RCP< const PureBasis > getBasis() const
PHX::MDField< Scalar > buildArray(const std::string &str, int d0) const
PHX::MDField< double, Cell, BASIS, IP, Dim > weighted_basis_vector
EElementSpace getElementSpace() const
#define BASIS_VALUES_INSTANTIATION(SCALAR)
int dimension() const
Returns the dimension of the basis from the topology.
Teuchos::RCP< Intrepid2::Basis< double, Kokkos::DynRankView< double, PHX::Device > > > getIntrepid2Basis() const
int numCells() const
Returns the number of cells in the data layouts.
void setupArrays(const Teuchos::RCP< const panzer::BasisIRLayout > &basis, bool computeDerivatives=true)
Sizes/allocates memory for arrays.
Teuchos::RCP< const shards::CellTopology > getCellTopology() const
void applyOrientations(const PHX::MDField< Scalar, Cell, BASIS > &orientations)
Method to apply orientations to a basis values container.
#define TEUCHOS_ASSERT(assertion_test)
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, IP, Dim, void, void, void, void, void > &cell_cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv)
PureBasis::EElementSpace getElementSpace() const
int cardinality() const
Returns the number of basis coefficients.
PHX::MDField< Scalar > ArrayDynamic