62 const Teuchos::ParameterList& p,
64 : rowIndexers_(rIndexers)
65 , colIndexers_(cIndexers)
66 , globalDataKey_(
"Residual Scatter Container")
68 std::string scatterName = p.get<std::string>(
"Scatter Name");
70 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
73 const std::vector<std::string>& names =
74 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
77 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
79 Teuchos::RCP<PHX::DataLayout> dl =
80 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
82 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
83 local_side_id_ = p.get<
int>(
"Local Side ID");
86 scatterFields_.resize(names.size());
87 for (std::size_t eq = 0; eq < names.size(); ++eq) {
88 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
91 this->addDependentField(scatterFields_[eq]);
94 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
96 applyBC_.resize(names.size());
97 for (std::size_t eq = 0; eq < names.size(); ++eq) {
98 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
99 this->addDependentField(applyBC_[eq]);
104 this->addEvaluatedField(*scatterHolder_);
106 if (p.isType<std::string>(
"Global Data Key"))
107 globalDataKey_ = p.get<std::string>(
"Global Data Key");
109 if(colIndexers_.size()==0)
110 colIndexers_ = rowIndexers_;
112 this->setName(scatterName+
" Scatter Dirichlet Residual : BlockedEpetra (Hessian)");
167 using Teuchos::ArrayRCP;
168 using Teuchos::ptrFromRef;
169 using Teuchos::rcp_dynamic_cast;
171 using Thyra::SpmdVectorBase;
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;
182 std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,
panzer::pair_hash> jacEpetraBlocks;
189 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
190 int rowIndexer = indexerIds_[fieldIndex];
191 int subFieldNum = subFieldIds_[fieldIndex];
194 Teuchos::ArrayRCP<double> local_dc;
195 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
196 ->getNonconstLocalData(ptrFromRef(local_dc));
198 auto subRowIndexer = rowIndexers_[rowIndexer];
199 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
200 const std::vector<int> & subElmtOffset = subIndicePair.first;
201 const std::vector<int> & subBasisIdMap = subIndicePair.second;
204 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
205 std::size_t cellLocalId = localCellIds[worksetCellIndex];
207 auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
210 for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
211 int offset = subElmtOffset[basis];
212 int lid = rLIDs[offset];
216 int basisId = subBasisIdMap[basis];
219 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
223 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
224 int start = blockOffsets[colIndexer];
225 int end = blockOffsets[colIndexer+1];
231 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
232 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
234 if(subJac==Teuchos::null) {
235 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
238 if(Teuchos::is_null(tOp))
241 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
242 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
243 jacEpetraBlocks[blockIndex] = subJac;
247 int * rowIndices = 0;
248 double * rowValues = 0;
250 subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
252 for(
int i=0;i<numEntries;i++)
256 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
261 std::vector<double> jacRow(scatterField.size(),0.0);
263 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
264 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
266 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
267 int start = blockOffsets[colIndexer];
268 int end = blockOffsets[colIndexer+1];
273 auto subColIndexer = colIndexers_[colIndexer];
274 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
276 TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
279 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
280 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
283 if(subJac==Teuchos::null) {
284 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
287 if(Teuchos::is_null(tOp))
290 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
291 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
292 jacEpetraBlocks[blockIndex] = subJac;
296 int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
298 std::stringstream ss;
299 ss <<
"Failed inserting row: " <<
" (" << lid <<
"): ";
300 for(
int i=0;i<end-start;i++)
301 ss << cLIDs[i] <<
" ";
303 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
305 ss <<
"scatter field = ";
306 scatterFields_[fieldIndex].print(ss);
309 TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());