43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 57 #include "Phalanx_DataLayout_MDALayout.hpp" 59 #include "Teuchos_FancyOStream.hpp" 66 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
70 : globalIndexer_(indexer)
71 , globalDataKey_(
"Residual Scatter Container")
73 std::string scatterName = p.
get<std::string>(
"Scatter Name");
78 const std::vector<std::string>& names =
85 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
91 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
92 local_side_id_ = p.
get<
int>(
"Local Side ID");
97 for (std::size_t eq = 0; eq < names.size(); ++eq) {
98 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
104 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
106 applyBC_.resize(names.size());
107 for (std::size_t eq = 0; eq < names.size(); ++eq) {
108 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
109 this->addDependentField(applyBC_[eq]);
114 this->addEvaluatedField(*scatterHolder_);
116 if (p.
isType<std::string>(
"Global Data Key"))
117 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
119 this->setName(scatterName+
" Scatter Residual");
123 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
133 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
134 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
140 this->utils.setFieldData(applyBC_[fd],fm);
148 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
153 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(globalDataKey_));
155 if(tpetraContainer_==Teuchos::null) {
158 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
160 dirichletCounter_ = Teuchos::null;
165 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
167 dirichletCounter_ = tpetraContainer->get_f();
173 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
177 std::vector<GO> GIDs;
178 std::vector<LO> LIDs;
181 std::string blockId = this->wda(workset).block_id;
182 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
185 tpetraContainer_->get_f() :
186 tpetraContainer_->get_x();
198 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
199 std::size_t cellLocalId = localCellIds[worksetCellIndex];
201 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
205 for(std::size_t i=0;i<GIDs.size();i++)
206 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
209 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
210 int fieldNum = fieldIds_[fieldIndex];
214 const std::pair<std::vector<int>,std::vector<int> > & indicePair
215 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
216 const std::vector<int> & elmtOffset = indicePair.first;
217 const std::vector<int> & basisIdMap = indicePair.second;
226 int basisId = basisIdMap[
basis];
229 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
232 r_array[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
239 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
263 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
267 : globalIndexer_(indexer)
268 , globalDataKey_(
"Residual Scatter Container")
270 std::string scatterName = p.
get<std::string>(
"Scatter Name");
275 const std::vector<std::string>& names =
282 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
288 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
289 local_side_id_ = p.
get<
int>(
"Local Side ID");
294 for (std::size_t eq = 0; eq < names.size(); ++eq) {
295 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
301 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
303 applyBC_.resize(names.size());
304 for (std::size_t eq = 0; eq < names.size(); ++eq) {
305 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
306 this->addDependentField(applyBC_[eq]);
311 this->addEvaluatedField(*scatterHolder_);
313 if (p.
isType<std::string>(
"Global Data Key"))
314 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
316 this->setName(scatterName+
" Scatter Tangent");
320 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
330 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
331 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
337 this->utils.setFieldData(applyBC_[fd],fm);
345 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
350 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(globalDataKey_));
352 if(tpetraContainer_==Teuchos::null) {
355 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
357 dirichletCounter_ = Teuchos::null;
362 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
364 dirichletCounter_ = tpetraContainer->get_f();
369 using Teuchos::rcp_dynamic_cast;
372 std::vector<std::string> activeParameters =
377 dfdp_vectors_.clear();
378 for(std::size_t i=0;i<activeParameters.size();i++) {
380 rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(activeParameters[i]),
true)->get_f();
382 dfdp_vectors_.push_back(vec_array);
387 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
391 std::vector<GO> GIDs;
392 std::vector<LO> LIDs;
395 std::string blockId = this->wda(workset).block_id;
396 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
399 tpetraContainer_->get_f() :
400 tpetraContainer_->get_x();
412 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
413 std::size_t cellLocalId = localCellIds[worksetCellIndex];
415 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
419 for(std::size_t i=0;i<GIDs.size();i++)
420 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
423 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
424 int fieldNum = fieldIds_[fieldIndex];
428 const std::pair<std::vector<int>,std::vector<int> > & indicePair
429 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
430 const std::vector<int> & elmtOffset = indicePair.first;
431 const std::vector<int> & basisIdMap = indicePair.second;
440 int basisId = basisIdMap[
basis];
443 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
451 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
452 dfdp_vectors_[d][lid] = 0.0;
454 for(
int d=0;d<
value.size();d++) {
455 dfdp_vectors_[d][lid] =
value.fastAccessDx(d);
463 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
477 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
478 dfdp_vectors_[d][lid] = 0.0;
480 for(
int d=0;d<
value.size();d++) {
481 dfdp_vectors_[d][lid] =
value.fastAccessDx(d);
496 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
500 : globalIndexer_(indexer)
501 , globalDataKey_(
"Residual Scatter Container")
503 std::string scatterName = p.
get<std::string>(
"Scatter Name");
508 const std::vector<std::string>& names =
517 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
518 local_side_id_ = p.
get<
int>(
"Local Side ID");
522 for (std::size_t eq = 0; eq < names.size(); ++eq) {
523 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
529 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
531 applyBC_.resize(names.size());
532 for (std::size_t eq = 0; eq < names.size(); ++eq) {
533 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
534 this->addDependentField(applyBC_[eq]);
539 this->addEvaluatedField(*scatterHolder_);
541 if (p.
isType<std::string>(
"Global Data Key"))
542 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
544 this->setName(scatterName+
" Scatter Residual (Jacobian)");
548 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
557 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
558 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
564 this->utils.setFieldData(applyBC_[fd],fm);
573 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
578 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(globalDataKey_));
580 if(tpetraContainer_==Teuchos::null) {
583 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
585 dirichletCounter_ = Teuchos::null;
590 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
592 dirichletCounter_ = tpetraContainer->get_f();
598 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
602 std::vector<GO> GIDs;
605 std::string blockId = this->wda(workset).block_id;
606 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
620 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
621 std::size_t cellLocalId = localCellIds[worksetCellIndex];
623 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
624 const std::vector<LO> & LIDs = globalIndexer_->getElementLIDs(cellLocalId);
627 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
628 int fieldNum = fieldIds_[fieldIndex];
631 const std::pair<std::vector<int>,std::vector<int> > & indicePair
632 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
633 const std::vector<int> & elmtOffset = indicePair.first;
634 const std::vector<int> & basisIdMap = indicePair.second;
643 int basisId = basisIdMap[
basis];
646 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
651 std::size_t sz = Jac->getNumEntriesInLocalRow(lid);
652 std::size_t numEntries = 0;
657 Jac->getLocalRowCopy(lid,rowIndices,rowValues,numEntries);
659 for(std::size_t i=0;i<numEntries;i++)
662 Jac->replaceLocalValues(lid,rowIndices,rowValues);
668 r_array[lid] = scatterField.val();
672 std::vector<double> jacRow(scatterField.size(),0.0);
674 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
675 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
678 Jac->replaceGlobalValues(gid, GIDs, jacRow);
panzer::Traits::Tangent::ScalarT ScalarT
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void resize(const size_type n, const T &val=T())
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pushes residual values into the residual vector for a Newton-based solve.
bool isType(const std::string &name) const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
bool isParameter(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>
panzer::Traits::Jacobian::ScalarT ScalarT