43 #ifndef __Panzer_ScatterResidual_BlockedEpetra_Hessian_impl_hpp__
44 #define __Panzer_ScatterResidual_BlockedEpetra_Hessian_impl_hpp__
47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
57 template<
typename TRAITS,
typename LO,
typename GO>
62 bool useDiscreteAdjoint)
63 : rowIndexers_(rIndexers)
64 , colIndexers_(cIndexers)
65 , globalDataKey_(
"Residual Scatter Container")
66 , useDiscreteAdjoint_(useDiscreteAdjoint)
68 std::string scatterName = p.
get<std::string>(
"Scatter Name");
73 const std::vector<std::string>& names =
83 scatterFields_.resize(names.size());
84 for (std::size_t eq = 0; eq < names.size(); ++eq) {
85 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
88 this->addDependentField(scatterFields_[eq]);
92 this->addEvaluatedField(*scatterHolder_);
94 if (p.
isType<std::string>(
"Global Data Key"))
95 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
96 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
97 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
100 if(useDiscreteAdjoint)
103 if(colIndexers_.size()==0)
104 colIndexers_ = rowIndexers_;
106 this->setName(scatterName+
" Scatter Residual BlockedEpetra (Hessian)");
109 template<
typename TRAITS,
typename LO,
typename GO>
115 indexerIds_.resize(scatterFields_.size());
116 subFieldIds_.resize(scatterFields_.size());
119 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
121 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
124 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
128 template<
typename TRAITS,
typename LO,
typename GO>
134 using Teuchos::rcp_dynamic_cast;
140 RCP<const BLOC> blockedContainer = rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_));
141 RCP<const ELOC> epetraContainer = rcp_dynamic_cast<
const ELOC>(d.gedc->getDataObject(globalDataKey_));
144 if(blockedContainer!=Teuchos::null) {
147 else if(epetraContainer!=Teuchos::null) {
156 template<
typename TRAITS,
typename LO,
typename GO>
163 using Teuchos::ptrFromRef;
164 using Teuchos::rcp_dynamic_cast;
167 using Thyra::SpmdVectorBase;
171 std::vector<double> jacRow;
174 std::string blockId = this->
wda(workset).block_id;
175 const std::vector<std::size_t> & localCellIds = this->
wda(workset).cell_local_ids;
177 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
179 std::vector<int> blockOffsets;
185 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
186 int rowIndexer = indexerIds_[fieldIndex];
187 int subFieldNum = subFieldIds_[fieldIndex];
189 auto subRowIndexer = rowIndexers_[rowIndexer];
190 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
193 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
194 std::size_t cellLocalId = localCellIds[worksetCellIndex];
196 auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
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 r_lid = rLIDs[rowOffset];
205 jacRow.resize(scatterField.size());
208 if(scatterField.size() == 0)
211 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
212 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
215 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
216 int start = blockOffsets[colIndexer];
217 int end = blockOffsets[colIndexer+1];
222 auto subColIndexer = colIndexers_[colIndexer];
223 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
228 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
229 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
232 if(subJac==Teuchos::null) {
241 jacEpetraBlocks[blockIndex] = subJac;
246 int err = subJac->
SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs[0]);
249 std::stringstream ss;
250 ss <<
"Failed inserting row: " <<
"LID = " << r_lid <<
": ";
251 for(
int i=0;i<end-start;i++)
252 ss << cLIDs[i] <<
" ";
254 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
256 ss <<
"scatter field = ";
257 scatterFields_[fieldIndex].print(ss);
261 for(
int i=start;i<end;i++)
262 ss << jacRow[i] <<
" ";
265 std::cout << ss.str() << std::endl;
ScatterResidual_BlockedEpetra(const Teuchos::ParameterList &p)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
panzer::Traits::Hessian::ScalarT ScalarT
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis, std::vector< int > &blockOffsets)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void evaluateFields(typename TRAITS::EvalData d)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Pushes residual values into the residual vector for a Newton-based solve.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
bool isType(const std::string &name) const
WorksetDetailsAccessor wda
#define TEUCHOS_ASSERT(assertion_test)