43 #ifndef __Panzer_ScatterResidual_Epetra_Hessian_impl_hpp__ 44 #define __Panzer_ScatterResidual_Epetra_Hessian_impl_hpp__ 47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT 49 #include "Teuchos_RCP.hpp" 50 #include "Teuchos_Assert.hpp" 52 #include "Phalanx_DataLayout.hpp" 54 #include "Epetra_Map.h" 55 #include "Epetra_Vector.h" 56 #include "Epetra_CrsMatrix.h" 66 #include "Phalanx_DataLayout_MDALayout.hpp" 68 #include "Teuchos_FancyOStream.hpp" 74 template<
typename TRAITS,
typename LO,
typename GO>
78 const Teuchos::ParameterList& p,
79 bool useDiscreteAdjoint)
80 : globalIndexer_(indexer)
81 , colGlobalIndexer_(cIndexer)
82 , globalDataKey_(
"Residual Scatter Container")
83 , useDiscreteAdjoint_(useDiscreteAdjoint)
85 std::string scatterName = p.get<std::string>(
"Scatter Name");
87 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
90 const std::vector<std::string>& names =
91 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
94 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
96 Teuchos::RCP<PHX::DataLayout> dl =
97 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
101 for (std::size_t eq = 0; eq < names.size(); ++eq) {
102 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
109 this->addEvaluatedField(*scatterHolder_);
111 if (p.isType<std::string>(
"Global Data Key"))
112 globalDataKey_ = p.get<std::string>(
"Global Data Key");
113 if (p.isType<
bool>(
"Use Discrete Adjoint"))
114 useDiscreteAdjoint = p.get<
bool>(
"Use Discrete Adjoint");
117 if(useDiscreteAdjoint)
118 { TEUCHOS_ASSERT(colGlobalIndexer_==globalIndexer_); }
120 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
124 template<
typename TRAITS,
typename LO,
typename GO>
133 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
134 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
142 template<
typename TRAITS,
typename LO,
typename GO>
147 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc.getDataObject(globalDataKey_));
149 if(epetraContainer_==Teuchos::null) {
151 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
152 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
157 template<
typename TRAITS,
typename LO,
typename GO>
161 std::vector<int> cLIDs, rLIDs;
162 std::vector<double> jacRow;
164 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
167 std::string blockId = this->wda(workset).block_id;
168 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
170 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
171 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
173 const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> >&
174 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
182 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
183 std::size_t cellLocalId = localCellIds[worksetCellIndex];
185 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
186 cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
187 if (Teuchos::nonnull(workset.other)) {
188 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
189 const std::vector<int> other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
190 cLIDs.insert(cLIDs.end(), other_cLIDs.begin(), other_cLIDs.end());
194 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
195 int fieldNum = fieldIds_[fieldIndex];
196 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
199 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
200 const ScalarT scatterField = (
scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
201 int rowOffset = elmtOffset[rowBasisNum];
202 int row = rLIDs[rowOffset];
205 jacRow.resize(scatterField.size());
207 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
208 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
211 int err = Jac->SumIntoMyValues(
213 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
216 TEUCHOS_ASSERT_EQUALITY(err,0);
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T * ptrFromStlVector(std::vector< T > &v)
Pushes residual values into the residual vector for a Newton-based solve.