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"
74template<
typename TRAITS,
typename LO,
typename GO>
77 const Teuchos::RCP<const panzer::GlobalIndexer> & cIndexer,
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;
100 scatterFields_.resize(names.size());
101 for (std::size_t eq = 0; eq < names.size(); ++eq) {
102 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
105 this->addDependentField(scatterFields_[eq]);
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)");
124template<
typename TRAITS,
typename LO,
typename GO>
129 fieldIds_.resize(scatterFields_.size());
131 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
133 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
134 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
139template<
typename TRAITS,
typename LO,
typename GO>
144 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
146 if(epetraContainer_==Teuchos::null) {
148 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
149 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
154template<
typename TRAITS,
typename LO,
typename GO>
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::GlobalIndexer>&
174 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
182 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
183 std::size_t cellLocalId = localCellIds[worksetCellIndex];
185 auto rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
186 auto initial_cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
188 for (
int i(0); i < static_cast<int>(initial_cLIDs.extent(0)); ++i)
189 cLIDs.push_back(initial_cLIDs(i));
190 if (Teuchos::nonnull(workset.other)) {
191 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
192 auto other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
193 for (
int i(0); i < static_cast<int>(other_cLIDs.extent(0)); ++i)
194 cLIDs.push_back(other_cLIDs(i));
198 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
199 int fieldNum = fieldIds_[fieldIndex];
200 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
203 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
204 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
205 int rowOffset = elmtOffset[rowBasisNum];
206 int row = rLIDs[rowOffset];
209 jacRow.resize(scatterField.size());
211 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
212 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
215 int err = Jac->SumIntoMyValues(
217 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
220 TEUCHOS_ASSERT_EQUALITY(err,0);
panzer::Traits::Hessian::ScalarT ScalarT
Pushes residual values into the residual vector for a Newton-based solve.
T * ptrFromStlVector(std::vector< T > &v)