43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP 46 #include "Teuchos_RCP.hpp" 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 51 #include "Epetra_Map.h" 52 #include "Epetra_Vector.h" 53 #include "Epetra_CrsMatrix.h" 61 #include "Phalanx_DataLayout_MDALayout.hpp" 63 #include "Teuchos_FancyOStream.hpp" 70 template<
typename TRAITS,
typename LO,
typename GO>
74 const Teuchos::ParameterList& p)
75 : globalIndexer_(indexer)
76 , globalDataKey_(
"Residual Scatter Container")
78 std::string scatterName = p.get<std::string>(
"Scatter Name");
80 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
83 const std::vector<std::string>& names =
84 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
87 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
90 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
92 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
93 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
94 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
96 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
97 local_side_id_ = p.get<
int>(
"Local Side ID");
102 for (std::size_t eq = 0; eq < names.size(); ++eq) {
103 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
109 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
111 applyBC_.resize(names.size());
112 for (std::size_t eq = 0; eq < names.size(); ++eq) {
113 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
114 this->addDependentField(applyBC_[eq]);
119 this->addEvaluatedField(*scatterHolder_);
121 if (p.isType<std::string>(
"Global Data Key"))
122 globalDataKey_ = p.get<std::string>(
"Global Data Key");
124 this->setName(scatterName+
" Scatter Residual");
128 template<
typename TRAITS,
typename LO,
typename GO>
138 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
139 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
144 this->utils.setFieldData(applyBC_[fd],fm);
152 template<
typename TRAITS,
typename LO,
typename GO>
159 if(epetraContainer_==Teuchos::null) {
161 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
164 dirichletCounter_ = Teuchos::null;
168 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer
171 dirichletCounter_ = epetraContainer->
get_f();
172 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
177 template<
typename TRAITS,
typename LO,
typename GO>
181 std::vector<int> LIDs;
184 std::string blockId = this->wda(workset).block_id;
185 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
187 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
188 Teuchos::RCP<Epetra_Vector> r = (!scatterIC_) ?
189 epetraContainer->get_f() :
190 epetraContainer->get_x();
197 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
198 std::size_t cellLocalId = localCellIds[worksetCellIndex];
200 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
203 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
204 int fieldNum = fieldIds_[fieldIndex];
208 const std::pair<std::vector<int>,std::vector<int> > & indicePair
209 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
210 const std::vector<int> & elmtOffset = indicePair.first;
211 const std::vector<int> & basisIdMap = indicePair.second;
220 int basisId = basisIdMap[
basis];
223 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
226 (*r)[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
229 if(dirichletCounter_!=Teuchos::null)
230 (*dirichletCounter_)[lid] = 1.0;
234 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
246 if(dirichletCounter_!=Teuchos::null)
247 (*dirichletCounter_)[lid] = 1.0;
259 template<
typename TRAITS,
typename LO,
typename GO>
263 const Teuchos::ParameterList& p)
264 : globalIndexer_(indexer)
265 , globalDataKey_(
"Residual Scatter Container")
267 std::string scatterName = p.get<std::string>(
"Scatter Name");
269 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
272 const std::vector<std::string>& names =
273 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
276 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
279 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
281 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
282 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
283 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
285 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
286 local_side_id_ = p.get<
int>(
"Local Side ID");
291 for (std::size_t eq = 0; eq < names.size(); ++eq) {
292 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
298 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
300 applyBC_.resize(names.size());
301 for (std::size_t eq = 0; eq < names.size(); ++eq) {
302 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
303 this->addDependentField(applyBC_[eq]);
308 this->addEvaluatedField(*scatterHolder_);
310 if (p.isType<std::string>(
"Global Data Key"))
311 globalDataKey_ = p.get<std::string>(
"Global Data Key");
313 this->setName(scatterName+
" Scatter Tangent");
317 template<
typename TRAITS,
typename LO,
typename GO>
327 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
328 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
333 this->utils.setFieldData(applyBC_[fd],fm);
341 template<
typename TRAITS,
typename LO,
typename GO>
348 if(epetraContainer_==Teuchos::null) {
350 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
353 dirichletCounter_ = Teuchos::null;
357 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer
360 dirichletCounter_ = epetraContainer->
get_f();
361 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
365 using Teuchos::rcp_dynamic_cast;
368 std::vector<std::string> activeParameters =
372 TEUCHOS_ASSERT(!scatterIC_);
373 dfdp_vectors_.clear();
374 for(std::size_t i=0;i<activeParameters.size();i++) {
375 RCP<Epetra_Vector> vec = rcp_dynamic_cast<
EpetraLinearObjContainer>(d.gedc.getDataObject(activeParameters[i]),
true)->get_f();
376 dfdp_vectors_.push_back(vec);
381 template<
typename TRAITS,
typename LO,
typename GO>
385 std::vector<int> LIDs;
388 std::string blockId = this->wda(workset).block_id;
389 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
391 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
392 Teuchos::RCP<Epetra_Vector> r = (!scatterIC_) ?
393 epetraContainer->get_f() :
394 epetraContainer->get_x();
401 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
402 std::size_t cellLocalId = localCellIds[worksetCellIndex];
405 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
408 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
409 int fieldNum = fieldIds_[fieldIndex];
413 const std::pair<std::vector<int>,std::vector<int> > & indicePair
414 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
415 const std::vector<int> & elmtOffset = indicePair.first;
416 const std::vector<int> & basisIdMap = indicePair.second;
425 int basisId = basisIdMap[
basis];
428 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
436 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
437 (*dfdp_vectors_[d])[lid] = 0.0;
439 for(
int d=0;d<
value.size();d++) {
440 (*dfdp_vectors_[d])[lid] =
value.fastAccessDx(d);
444 if(dirichletCounter_!=Teuchos::null)
445 (*dirichletCounter_)[lid] = 1.0;
449 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
463 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
464 (*dfdp_vectors_[d])[lid] = 0.0;
466 for(
int d=0;d<
value.size();d++) {
467 (*dfdp_vectors_[d])[lid] =
value.fastAccessDx(d);
471 if(dirichletCounter_!=Teuchos::null)
472 (*dirichletCounter_)[lid] = 1.0;
483 template<
typename TRAITS,
typename LO,
typename GO>
487 const Teuchos::ParameterList& p)
488 : globalIndexer_(indexer)
489 , colGlobalIndexer_(cIndexer)
490 , globalDataKey_(
"Residual Scatter Container")
491 , preserveDiagonal_(false)
493 std::string scatterName = p.get<std::string>(
"Scatter Name");
495 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
498 const std::vector<std::string>& names =
499 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
502 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
504 Teuchos::RCP<PHX::DataLayout> dl =
505 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
507 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
508 local_side_id_ = p.get<
int>(
"Local Side ID");
512 for (std::size_t eq = 0; eq < names.size(); ++eq) {
513 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
519 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
521 applyBC_.resize(names.size());
522 for (std::size_t eq = 0; eq < names.size(); ++eq) {
523 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
524 this->addDependentField(applyBC_[eq]);
529 this->addEvaluatedField(*scatterHolder_);
531 if (p.isType<std::string>(
"Global Data Key"))
532 globalDataKey_ = p.get<std::string>(
"Global Data Key");
534 if (p.isType<
bool>(
"Preserve Diagonal"))
535 preserveDiagonal_ = p.get<
bool>(
"Preserve Diagonal");
537 this->setName(scatterName+
" Scatter Residual (Jacobian)");
541 template<
typename TRAITS,
typename LO,
typename GO>
551 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
552 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
557 this->utils.setFieldData(applyBC_[fd],fm);
566 template<
typename TRAITS,
typename LO,
typename GO>
573 if(epetraContainer_==Teuchos::null) {
575 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
578 dirichletCounter_ = Teuchos::null;
582 Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc.getDataObject(
"Dirichlet Counter");
583 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<
EpetraLinearObjContainer>(dataContainer,
true);
585 dirichletCounter_ = epetraContainer->
get_f();
586 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
591 template<
typename TRAITS,
typename LO,
typename GO>
595 std::vector<int> cLIDs, rLIDs;
596 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
599 std::string blockId = this->wda(workset).block_id;
600 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
602 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
603 TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
604 Teuchos::RCP<Epetra_Vector> r = epetraContainer->get_f();
605 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
612 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
613 std::size_t cellLocalId = localCellIds[worksetCellIndex];
615 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
617 cLIDs = colGlobalIndexer_->getElementLIDs(cellLocalId);
622 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
623 int fieldNum = fieldIds_[fieldIndex];
626 const std::pair<std::vector<int>,std::vector<int> > & indicePair
627 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
628 const std::vector<int> & elmtOffset = indicePair.first;
629 const std::vector<int> & basisIdMap = indicePair.second;
638 int basisId = basisIdMap[
basis];
641 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
647 int * rowIndices = 0;
648 double * rowValues = 0;
650 Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
652 for(
int i=0;i<numEntries;i++) {
653 if(preserveDiagonal_) {
654 if(row!=rowIndices[i])
666 (*r)[row] = scatterField.val();
667 if(dirichletCounter_!=Teuchos::null) {
669 (*dirichletCounter_)[row] = 1.0;
673 std::vector<double> jacRow(scatterField.size(),0.0);
675 if(!preserveDiagonal_) {
676 int err = Jac->ReplaceMyValues(row, cLIDs.size(), scatterField.dx(), &cLIDs[0]);
677 TEUCHOS_ASSERT(err==0);
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
panzer::Traits::Tangent::ScalarT ScalarT
panzer::Traits::Jacobian::ScalarT ScalarT
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
const Teuchos::RCP< Epetra_Vector > get_f() const
Pushes residual values into the residual vector for a Newton-based solve.