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 =
84 for (std::size_t eq = 0; eq < names.size(); ++eq) {
85 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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>
121 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
124 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
131 template<
typename TRAITS,
typename LO,
typename GO>
133 ScatterResidual_BlockedEpetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
134 preEvaluate(
typename TRAITS::PreEvalData d)
137 using Teuchos::rcp_dynamic_cast;
139 typedef BlockedEpetraLinearObjContainer BLOC;
140 typedef BlockedEpetraLinearObjContainer ELOC;
143 RCP<const BLOC> blockedContainer = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
144 RCP<const ELOC> epetraContainer = rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
147 if(blockedContainer!=Teuchos::null) {
150 else if(epetraContainer!=Teuchos::null) {
152 RCP<Thyra::LinearOpBase<double> > J = blockedContainer->get_A_th();
159 template<
typename TRAITS,
typename LO,
typename GO>
166 using Teuchos::ptrFromRef;
167 using Teuchos::rcp_dynamic_cast;
170 using Thyra::SpmdVectorBase;
174 std::vector<double> jacRow;
177 std::string blockId = this->wda(workset).block_id;
178 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
180 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
182 std::vector<int> blockOffsets;
188 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
189 int rowIndexer = indexerIds_[fieldIndex];
190 int subFieldNum = subFieldIds_[fieldIndex];
192 auto subRowIndexer = rowIndexers_[rowIndexer];
193 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
196 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
197 std::size_t cellLocalId = localCellIds[worksetCellIndex];
199 const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
202 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
203 const ScalarT scatterField = (
scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
204 int rowOffset = elmtOffset[rowBasisNum];
205 int r_lid = rLIDs[rowOffset];
208 jacRow.resize(scatterField.size());
211 if(scatterField.size() == 0)
214 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
215 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
218 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
219 int start = blockOffsets[colIndexer];
220 int end = blockOffsets[colIndexer+1];
225 auto subColIndexer = colIndexers_[colIndexer];
226 const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
231 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
235 if(subJac==Teuchos::null) {
244 jacEpetraBlocks[blockIndex] = subJac;
249 int err = subJac->
SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs[0]);
252 std::stringstream ss;
253 ss <<
"Failed inserting row: " <<
"LID = " << r_lid <<
": ";
254 for(
int i=0;i<end-start;i++)
255 ss << cLIDs[i] <<
" ";
257 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
259 ss <<
"scatter field = ";
264 for(
int i=start;i<end;i++)
265 ss << jacRow[i] <<
" ";
268 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)
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
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)
PHX::MDField< ScalarT > vector
void evaluateFields(typename TRAITS::EvalData d)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isType(const std::string &name) const
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
#define TEUCHOS_ASSERT(assertion_test)