43 #ifndef PANZER_EVALUATOR_GRADBASISCROSSVECTOR_IMPL_HPP 44 #define PANZER_EVALUATOR_GRADBASISCROSSVECTOR_IMPL_HPP 53 #include "Intrepid2_FunctionSpaceTools.hpp" 56 #include "Kokkos_ViewFactory.hpp" 70 template<
typename EvalT,
typename Traits>
74 const std::vector<std::string>& resNames,
75 const std::string& vecName,
79 const std::vector<std::string>& fmNames,
81 const Teuchos::RCP<PHX::DataLayout>& vecDL )
85 numDims_(resNames.size()),
86 numGradDims_(ir.dl_vector->extent(2)),
87 basisName_(basis.name())
94 using PHX::DataLayout;
98 using std::invalid_argument;
99 using std::logic_error;
105 TEUCHOS_TEST_FOR_EXCEPTION(
numDims_ != 3, invalid_argument,
"Error: " \
106 "Integrator_GradBasisCrossVector called with the number of residual " \
107 "names not equal to three.")
108 for (
const auto& name : resNames)
109 TEUCHOS_TEST_FOR_EXCEPTION(name ==
"", invalid_argument,
"Error: " \
110 "Integrator_GradBasisCrossVector called with an empty residual name.")
111 TEUCHOS_TEST_FOR_EXCEPTION(vecName ==
"", invalid_argument,
"Error: " \
112 "Integrator_GradBasisCrossVector called with an empty vector name.")
113 RCP<const PureBasis> tmpBasis = basis.
getBasis();
114 TEUCHOS_TEST_FOR_EXCEPTION(not tmpBasis->supportsGrad(), logic_error,
115 "Error: Integrator_GradBasisCrossVector: Basis of type \"" 116 << tmpBasis->name() <<
"\" does not support the gradient operator.")
118 if (not vecDL.is_null())
121 TEUCHOS_TEST_FOR_EXCEPTION(
122 tmpVecDL->extent(2) < ir.
dl_vector->extent(2), logic_error,
123 "Error: Integrator_GradBasisCrossVector: Dimension of space " \
124 "exceeds dimension of Vector Data Layout.")
125 TEUCHOS_TEST_FOR_EXCEPTION(
numDims_ !=
126 static_cast<int>(vecDL->extent(2)), logic_error,
"Error: " \
127 "Integrator_GradBasisCrossVector: The vector must be the same " \
128 "length as the number of residuals.")
131 "Error: Integrator_GradBasisCrossVector: The vector must have at " \
132 "least as many components as there are dimensions in the mesh.")
135 vector_ = MDField<const ScalarT, Cell, IP, Dim>(vecName, tmpVecDL);
136 this->addDependentField(
vector_);
140 fields_ = PHX::View<PHX::MDField<ScalarT, Cell, BASIS>*>(
"Integrator_GradBasisCrossVector::fields_", resNames.size());
144 for (
const auto& name : resNames)
147 for (std::size_t i=0; i<
fields_.extent(0); ++i) {
150 this->addContributedField(
field);
152 this->addEvaluatedField(
field);
157 kokkosFieldMults_ = View<View<const ScalarT**>*>(
"GradBasisCrossVector::KokkosFieldMultipliers", fmNames.size());
158 for (
const auto& name : fmNames)
165 string n(
"Integrator_GradBasisCrossVector (");
171 for (
size_t j=0; j <
fields_.extent(0) - 1; ++j)
172 n +=
fields_[j].fieldTag().name() +
", ";
182 template<
typename EvalT,
typename Traits>
185 const Teuchos::ParameterList& p)
189 p.get<const
std::vector<
std::string>>(
"Residual Names"),
190 p.get<
std::string>(
"Vector Name"),
193 p.get<double>(
"Multiplier"),
195 (
"Field Multipliers") ?
197 (
"Field Multipliers")) :
std::vector<
std::string>(),
198 p.isType<
Teuchos::RCP<
PHX::DataLayout>>(
"Data Layout Vector") ?
199 p.get<
Teuchos::RCP<
PHX::DataLayout>>(
"Data Layout Vector") :
202 using Teuchos::ParameterList;
207 p.validateParameters(*validParams);
215 template<
typename EvalT,
typename Traits>
226 for (
size_t i(0); i < fieldMults_.size(); ++i)
227 kokkosFieldMults_(i) = fieldMults_[i].get_static_view();
238 template<
typename EvalT,
typename Traits>
239 template<
int NUM_FIELD_MULT>
240 KOKKOS_INLINE_FUNCTION
245 const size_t& cell)
const 250 const int numBases(fields_[0].extent(1)), numQP(vector_.extent(1));
252 for (
int dim(0); dim < numDims_; ++dim)
253 for (
int basis(0); basis < numBases; ++basis)
254 fields_[dim](cell, basis) = 0.0;
259 const int X(0), Y(1), Z(2);
260 if (NUM_FIELD_MULT == 0)
262 if (numGradDims_ == 1)
264 for (
int qp(0); qp < numQP; ++qp)
266 tmp[Y] = multiplier_ * vector_(cell, qp, Y);
267 tmp[Z] = multiplier_ * vector_(cell, qp, Z);
268 for (
int basis(0); basis < numBases; ++basis)
270 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X);
271 fields_[Z](cell, basis) += -tmp[Y] * basis_(cell, basis, qp, X);
275 else if (numGradDims_ == 2)
277 for (
int qp(0); qp < numQP; ++qp)
279 tmp[X] = multiplier_ * vector_(cell, qp, X);
280 tmp[Y] = multiplier_ * vector_(cell, qp, Y);
281 tmp[Z] = multiplier_ * vector_(cell, qp, Z);
282 for (
int basis(0); basis < numBases; ++basis)
284 fields_[X](cell, basis) += -tmp[Z] * basis_(cell, basis, qp, Y);
285 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X);
286 fields_[Z](cell, basis) += tmp[X] * basis_(cell, basis, qp, Y) -
287 tmp[Y] * basis_(cell, basis, qp, X);
291 else if (numGradDims_ == 3)
293 for (
int qp(0); qp < numQP; ++qp)
295 tmp[X] = multiplier_ * vector_(cell, qp, X);
296 tmp[Y] = multiplier_ * vector_(cell, qp, Y);
297 tmp[Z] = multiplier_ * vector_(cell, qp, Z);
298 for (
int basis(0); basis < numBases; ++basis)
300 fields_[X](cell, basis) += tmp[Y] * basis_(cell, basis, qp, Z) -
301 tmp[Z] * basis_(cell, basis, qp, Y);
302 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X) -
303 tmp[X] * basis_(cell, basis, qp, Z);
304 fields_[Z](cell, basis) += tmp[X] * basis_(cell, basis, qp, Y) -
305 tmp[Y] * basis_(cell, basis, qp, X);
310 else if (NUM_FIELD_MULT == 1)
312 if (numGradDims_ == 1)
314 for (
int qp(0); qp < numQP; ++qp)
316 tmp[Y] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
317 vector_(cell, qp, Y);
318 tmp[Z] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
319 vector_(cell, qp, Z);
320 for (
int basis(0); basis < numBases; ++basis)
322 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X);
323 fields_[Z](cell, basis) += -tmp[Y] * basis_(cell, basis, qp, X);
327 else if (numGradDims_ == 2)
329 for (
int qp(0); qp < numQP; ++qp)
331 tmp[X] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
332 vector_(cell, qp, X);
333 tmp[Y] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
334 vector_(cell, qp, Y);
335 tmp[Z] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
336 vector_(cell, qp, Z);
337 for (
int basis(0); basis < numBases; ++basis)
339 fields_[X](cell, basis) += -tmp[Z] * basis_(cell, basis, qp, Y);
340 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X);
341 fields_[Z](cell, basis) += tmp[X] * basis_(cell, basis, qp, Y) -
342 tmp[Y] * basis_(cell, basis, qp, X);
346 else if (numGradDims_ == 3)
348 for (
int qp(0); qp < numQP; ++qp)
350 tmp[X] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
351 vector_(cell, qp, X);
352 tmp[Y] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
353 vector_(cell, qp, Y);
354 tmp[Z] = multiplier_ * kokkosFieldMults_(0)(cell, qp) *
355 vector_(cell, qp, Z);
356 for (
int basis(0); basis < numBases; ++basis)
358 fields_[X](cell, basis) += tmp[Y] * basis_(cell, basis, qp, Z) -
359 tmp[Z] * basis_(cell, basis, qp, Y);
360 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X) -
361 tmp[X] * basis_(cell, basis, qp, Z);
362 fields_[Z](cell, basis) += tmp[X] * basis_(cell, basis, qp, Y) -
363 tmp[Y] * basis_(cell, basis, qp, X);
370 const int numFieldMults(kokkosFieldMults_.extent(0));
371 if (numGradDims_ == 1)
373 for (
int qp(0); qp < numQP; ++qp)
375 tmp[Y] = multiplier_ * vector_(cell, qp, Y);
376 tmp[Z] = multiplier_ * vector_(cell, qp, Z);
377 for (
int fm(0); fm < numFieldMults; ++fm)
379 tmp[Y] *= kokkosFieldMults_(fm)(cell, qp);
380 tmp[Z] *= kokkosFieldMults_(fm)(cell, qp);
382 for (
int basis(0); basis < numBases; ++basis)
384 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X);
385 fields_[Z](cell, basis) += -tmp[Y] * basis_(cell, basis, qp, X);
389 else if (numGradDims_ == 2)
391 for (
int qp(0); qp < numQP; ++qp)
393 tmp[X] = multiplier_ * vector_(cell, qp, X);
394 tmp[Y] = multiplier_ * vector_(cell, qp, Y);
395 tmp[Z] = multiplier_ * vector_(cell, qp, Z);
396 for (
int fm(0); fm < numFieldMults; ++fm)
398 tmp[X] *= kokkosFieldMults_(fm)(cell, qp);
399 tmp[Y] *= kokkosFieldMults_(fm)(cell, qp);
400 tmp[Z] *= kokkosFieldMults_(fm)(cell, qp);
402 for (
int basis(0); basis < numBases; ++basis)
404 fields_[X](cell, basis) += -tmp[Z] * basis_(cell, basis, qp, Y);
405 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X);
406 fields_[Z](cell, basis) += tmp[X] * basis_(cell, basis, qp, Y) -
407 tmp[Y] * basis_(cell, basis, qp, X);
411 else if (numGradDims_ == 3)
413 for (
int qp(0); qp < numQP; ++qp)
415 tmp[X] = multiplier_ * vector_(cell, qp, X);
416 tmp[Y] = multiplier_ * vector_(cell, qp, Y);
417 tmp[Z] = multiplier_ * vector_(cell, qp, Z);
418 for (
int fm(0); fm < numFieldMults; ++fm)
420 tmp[X] *= kokkosFieldMults_(fm)(cell, qp);
421 tmp[Y] *= kokkosFieldMults_(fm)(cell, qp);
422 tmp[Z] *= kokkosFieldMults_(fm)(cell, qp);
424 for (
int basis(0); basis < numBases; ++basis)
426 fields_[X](cell, basis) += tmp[Y] * basis_(cell, basis, qp, Z) -
427 tmp[Z] * basis_(cell, basis, qp, Y);
428 fields_[Y](cell, basis) += tmp[Z] * basis_(cell, basis, qp, X) -
429 tmp[X] * basis_(cell, basis, qp, Z);
430 fields_[Z](cell, basis) += tmp[X] * basis_(cell, basis, qp, Y) -
431 tmp[Y] * basis_(cell, basis, qp, X);
443 template<
typename EvalT,
typename Traits>
449 using Kokkos::parallel_for;
450 using Kokkos::RangePolicy;
453 basis_ = this->wda(workset).bases[basisIndex_]->weighted_grad_basis;
458 if (fieldMults_.size() == 0)
460 else if (fieldMults_.size() == 1)
471 template<
typename EvalT,
typename TRAITS>
472 Teuchos::RCP<Teuchos::ParameterList>
478 using PHX::DataLayout;
481 using Teuchos::ParameterList;
486 RCP<ParameterList> p = rcp(
new ParameterList);
488 RCP<const vector<string>> resNames;
489 p->set(
"Residual Names", resNames);
490 p->set<
string>(
"Vector Name",
"?");
491 RCP<BasisIRLayout> basis;
492 p->set(
"Basis", basis);
493 RCP<IntegrationRule> ir;
495 p->set<
double>(
"Multiplier", 1.0);
496 RCP<const vector<string>> fms;
497 p->set(
"Field Multipliers", fms);
498 RCP<DataLayout> vecDL;
499 p->set(
"Data Layout Vector", vecDL);
506 #endif // PANZER_EVALUATOR_GRADBASISCROSSVECTOR_IMPL_HPP int num_cells
DEPRECATED - use: numCells()
std::vector< PHX::MDField< const ScalarT, Cell, IP > > fieldMults_
The (possibly empty) list of fields that are multipliers out in front of the integral ( ...
Integrator_GradBasisCrossVector(const panzer::EvaluatorStyle &evalStyle, const std::vector< std::string > &resNames, const std::string &vecName, const panzer::BasisIRLayout &basis, const panzer::IntegrationRule &ir, const double &multiplier=1, const std::vector< std::string > &fmNames=std::vector< std::string >(), const Teuchos::RCP< PHX::DataLayout > &vecDL=Teuchos::null)
Main Constructor.
This empty struct allows us to optimize operator()() depending on the number of field multipliers...
Teuchos::RCP< Teuchos::ParameterList > getValidParameters() const
Get Valid Parameters.
Teuchos::RCP< const PureBasis > getBasis() const
EvaluatorStyle
An indication of how an Evaluator will behave.
const panzer::EvaluatorStyle evalStyle_
An enum determining the behavior of this Evaluator.
void evaluateFields(typename Traits::EvalData d)
Evaluate Fields.
PHX::MDField< const ScalarT, Cell, IP, Dim > vector_
A field representing the vector-valued function we're integrating ( ).
double multiplier
The scalar multiplier out in front of the integral ( ).
Kokkos::DynRankView< typename InputArray::value_type, PHX::Device > createDynRankView(const InputArray &a, const std::string &name, const DimensionPack... dims)
Wrapper to simplify Panzer use of Sacado ViewFactory.
PHX::View< PHX::View< const ScalarT ** > * > kokkosFieldMults_
The PHX::View representation of the (possibly empty) list of fields that are multipliers out in front...
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
panzer::EvaluatorStyle evalStyle
The EvaluatorStyle of the parent Integrator_CurlBasisDotVector object.
typename EvalT::ScalarT ScalarT
The scalar type.
std::vector< std::string >::size_type getBasisIndex(std::string basis_name, const panzer::Workset &workset, WorksetDetailsAccessor &wda)
Returns the index in the workset bases for a particular BasisIRLayout name.
PHX::View< PHX::MDField< ScalarT, Cell, BASIS > * > fields_
The fields to which we'll contribute, or in which we'll store, the result of computing this integral...
KOKKOS_INLINE_FUNCTION void operator()(const FieldMultTag< NUM_FIELD_MULT > &tag, const std::size_t &cell) const
Perform the integration.
Teuchos::RCP< PHX::DataLayout > dl_vector
Data layout for vector fields.
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral...
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
Post-Registration Setup.
int numDims_
The number of dimensions associated with the vector.
int numGradDims_
The number of dimensions associated with the gradient.
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis>
Teuchos::RCP< const std::vector< panzer::Workset > > worksets_