43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 46 #include "Teuchos_RCP.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>
69 const Teuchos::ParameterList& p)
70 : globalIndexer_(indexer)
71 , globalDataKey_(
"Residual Scatter Container")
73 std::string scatterName = p.get<std::string>(
"Scatter Name");
75 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
78 const std::vector<std::string>& names =
79 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
82 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
85 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
87 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
88 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
89 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional ;
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) {
157 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
158 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
160 dirichletCounter_ = Teuchos::null;
164 Teuchos::RCP<LOC> tpetraContainer
165 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
167 dirichletCounter_ = tpetraContainer->
get_f();
168 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
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;
184 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
185 tpetraContainer_->get_f() :
186 tpetraContainer_->get_x();
188 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
189 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
198 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
199 std::size_t cellLocalId = localCellIds[worksetCellIndex];
201 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
204 LIDs.resize(GIDs.size());
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>
266 const Teuchos::ParameterList& p)
267 : globalIndexer_(indexer)
268 , globalDataKey_(
"Residual Scatter Container")
270 std::string scatterName = p.get<std::string>(
"Scatter Name");
272 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
275 const std::vector<std::string>& names =
276 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
279 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
282 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
284 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
285 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
286 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional ;
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) {
354 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
355 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
357 dirichletCounter_ = Teuchos::null;
361 Teuchos::RCP<LOC> tpetraContainer
362 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
364 dirichletCounter_ = tpetraContainer->
get_f();
365 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
369 using Teuchos::rcp_dynamic_cast;
372 std::vector<std::string> activeParameters =
376 TEUCHOS_ASSERT(!scatterIC_);
377 dfdp_vectors_.clear();
378 for(std::size_t i=0;i<activeParameters.size();i++) {
379 RCP<typename LOC::VectorType> vec =
380 rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(activeParameters[i]),
true)->get_f();
381 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
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;
398 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
399 tpetraContainer_->get_f() :
400 tpetraContainer_->get_x();
402 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
403 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
412 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
413 std::size_t cellLocalId = localCellIds[worksetCellIndex];
415 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
418 LIDs.resize(GIDs.size());
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>
499 const Teuchos::ParameterList& p)
500 : globalIndexer_(indexer)
501 , globalDataKey_(
"Residual Scatter Container")
503 std::string scatterName = p.get<std::string>(
"Scatter Name");
505 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
508 const std::vector<std::string>& names =
509 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
512 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
514 Teuchos::RCP<PHX::DataLayout> dl =
515 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
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) {
582 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
583 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
585 dirichletCounter_ = Teuchos::null;
589 Teuchos::RCP<LOC> tpetraContainer
590 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
592 dirichletCounter_ = tpetraContainer->
get_f();
593 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
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;
608 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
609 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
611 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
612 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
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;
653 Teuchos::Array<LO> rowIndices(sz);
654 Teuchos::Array<double> rowValues(sz);
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);
676 TEUCHOS_ASSERT(jacRow.size()==GIDs.size());
678 Jac->replaceGlobalValues(gid, GIDs, jacRow);
panzer::Traits::Tangent::ScalarT ScalarT
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
Pushes residual values into the residual vector for a Newton-based solve.
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
const Teuchos::RCP< VectorType > get_f() const
panzer::Traits::Jacobian::ScalarT ScalarT