43 #ifndef PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP 44 #define PANZER_SCATTER_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" 62 #include "Phalanx_DataLayout_MDALayout.hpp" 64 #include "Teuchos_FancyOStream.hpp" 70 template<
typename TRAITS,
typename LO,
typename GO>
74 const Teuchos::ParameterList& p,
75 bool useDiscreteAdjoint)
76 : globalIndexer_(indexer)
77 , globalDataKey_(
"Residual Scatter Container")
78 , useDiscreteAdjoint_(useDiscreteAdjoint)
80 std::string scatterName = p.get<std::string>(
"Scatter Name");
82 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
85 const std::vector<std::string>& names =
86 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
89 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
91 Teuchos::RCP<PHX::DataLayout> dl =
92 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
96 for (std::size_t eq = 0; eq < names.size(); ++eq) {
97 scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
104 this->addEvaluatedField(*scatterHolder_);
106 if (p.isType<std::string>(
"Global Data Key"))
107 globalDataKey_ = p.get<std::string>(
"Global Data Key");
109 this->setName(scatterName+
" Scatter Residual");
113 template<
typename TRAITS,
typename LO,
typename GO>
122 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
123 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
131 template<
typename TRAITS,
typename LO,
typename GO>
138 if(epetraContainer_==Teuchos::null) {
140 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
146 template<
typename TRAITS,
typename LO,
typename GO>
150 std::vector<int> LIDs;
153 std::string blockId = this->wda(workset).block_id;
154 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
156 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
164 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
165 std::size_t cellLocalId = localCellIds[worksetCellIndex];
167 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
170 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
171 int fieldNum = fieldIds_[fieldIndex];
172 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
188 template<
typename TRAITS,
typename LO,
typename GO>
192 const Teuchos::ParameterList& p,
193 bool useDiscreteAdjoint)
194 : globalIndexer_(indexer)
196 std::string scatterName = p.get<std::string>(
"Scatter Name");
198 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
201 const std::vector<std::string>& names =
202 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
205 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
207 Teuchos::RCP<PHX::DataLayout> dl =
208 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
212 for (std::size_t eq = 0; eq < names.size(); ++eq) {
213 scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
220 this->addEvaluatedField(*scatterHolder_);
222 this->setName(scatterName+
" Scatter Tangent");
226 template<
typename TRAITS,
typename LO,
typename GO>
235 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
236 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
244 template<
typename TRAITS,
typename LO,
typename GO>
249 using Teuchos::rcp_dynamic_cast;
252 std::vector<std::string> activeParameters =
255 dfdp_vectors_.clear();
256 for(std::size_t i=0;i<activeParameters.size();i++) {
257 RCP<Epetra_Vector> vec = rcp_dynamic_cast<
EpetraLinearObjContainer>(d.gedc.getDataObject(activeParameters[i]),
true)->get_f();
258 dfdp_vectors_.push_back(vec);
263 template<
typename TRAITS,
typename LO,
typename GO>
267 std::vector<int> LIDs;
270 std::string blockId = this->wda(workset).block_id;
271 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
279 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
280 std::size_t cellLocalId = localCellIds[worksetCellIndex];
282 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
285 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
286 int fieldNum = fieldIds_[fieldIndex];
287 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
296 for(
int d=0;d<
value.size();d++)
297 (*dfdp_vectors_[d])[lid] +=
value.fastAccessDx(d);
307 template<
typename TRAITS,
typename LO,
typename GO>
311 const Teuchos::ParameterList& p,
312 bool useDiscreteAdjoint)
313 : globalIndexer_(indexer)
314 , colGlobalIndexer_(cIndexer)
315 , globalDataKey_(
"Residual Scatter Container")
316 , useDiscreteAdjoint_(useDiscreteAdjoint)
318 std::string scatterName = p.get<std::string>(
"Scatter Name");
320 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
323 const std::vector<std::string>& names =
324 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
327 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
329 Teuchos::RCP<PHX::DataLayout> dl =
330 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
334 for (std::size_t eq = 0; eq < names.size(); ++eq) {
335 scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
342 this->addEvaluatedField(*scatterHolder_);
344 if (p.isType<std::string>(
"Global Data Key"))
345 globalDataKey_ = p.get<std::string>(
"Global Data Key");
346 if (p.isType<
bool>(
"Use Discrete Adjoint"))
347 useDiscreteAdjoint = p.get<
bool>(
"Use Discrete Adjoint");
350 if(useDiscreteAdjoint)
351 { TEUCHOS_ASSERT(colGlobalIndexer_==globalIndexer_); }
353 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
357 template<
typename TRAITS,
typename LO,
typename GO>
368 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
369 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
377 template<
typename TRAITS,
typename LO,
typename GO>
384 if(epetraContainer_==Teuchos::null) {
386 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
392 template<
typename TRAITS,
typename LO,
typename GO>
396 std::vector<int> cLIDs, rLIDs;
397 std::vector<double> jacRow;
399 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
402 std::string blockId = this->wda(workset).block_id;
403 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
405 Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
406 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
408 const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> >&
409 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
417 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
418 std::size_t cellLocalId = localCellIds[worksetCellIndex];
420 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
421 cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
422 if (Teuchos::nonnull(workset.other)) {
423 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
424 const std::vector<int> other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
425 cLIDs.insert(cLIDs.end(), other_cLIDs.begin(), other_cLIDs.end());
429 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
430 int fieldNum = fieldIds_[fieldIndex];
431 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
434 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
436 int rowOffset = elmtOffset[rowBasisNum];
437 int row = rLIDs[rowOffset];
441 r->SumIntoMyValue(row,0,scatterField.val());
444 if(useDiscreteAdjoint_) {
446 jacRow.resize(scatterField.size());
448 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
449 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
451 for(std::size_t c=0;c<cLIDs.size();c++) {
452 int err = Jac->SumIntoMyValues(cLIDs[c], 1, &jacRow[c],&row);
453 TEUCHOS_ASSERT_EQUALITY(err,0);
457 int err = Jac->SumIntoMyValues(
459 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
462 TEUCHOS_ASSERT_EQUALITY(err,0);
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
panzer::Traits::Tangent::ScalarT ScalarT
T * ptrFromStlVector(std::vector< T > &v)
Pushes residual values into the residual vector for a Newton-based solve.
panzer::Traits::Jacobian::ScalarT ScalarT
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.