43 #ifndef __Panzer_ScatterDirichletResidual_BlockedEpetra_Hessian_impl_hpp__ 44 #define __Panzer_ScatterDirichletResidual_BlockedEpetra_Hessian_impl_hpp__ 47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT 58 template<
typename TRAITS,
typename LO,
typename GO>
63 bool useDiscreteAdjoint)
64 : rowIndexers_(rIndexers)
65 , colIndexers_(cIndexers)
66 , globalDataKey_(
"Residual Scatter Container")
68 std::string scatterName = p.
get<std::string>(
"Scatter Name");
73 const std::vector<std::string>& names =
82 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
83 local_side_id_ = p.
get<
int>(
"Local Side ID");
87 for (std::size_t eq = 0; eq < names.size(); ++eq) {
88 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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)");
115 template<
typename TRAITS,
typename LO,
typename GO>
127 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
130 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
136 this->utils.setFieldData(applyBC_[fd],fm);
144 template<
typename TRAITS,
typename LO,
typename GO>
146 ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
147 preEvaluate(
typename TRAITS::PreEvalData d)
149 typedef BlockedEpetraLinearObjContainer BLOC;
151 using Teuchos::rcp_dynamic_cast;
155 = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
161 blockContainer = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_),
true);
167 template<
typename TRAITS,
typename LO,
typename GO>
174 using Teuchos::ptrFromRef;
175 using Teuchos::rcp_dynamic_cast;
177 using Thyra::SpmdVectorBase;
180 std::string blockId = this->wda(workset).block_id;
181 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
183 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
185 std::vector<int> blockOffsets;
195 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
196 int rowIndexer = indexerIds_[fieldIndex];
197 int subFieldNum = subFieldIds_[fieldIndex];
201 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
202 ->getNonconstLocalData(ptrFromRef(local_dc));
204 auto subRowIndexer = rowIndexers_[rowIndexer];
205 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
206 const std::vector<int> & subElmtOffset = subIndicePair.first;
207 const std::vector<int> & subBasisIdMap = subIndicePair.second;
210 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
211 std::size_t cellLocalId = localCellIds[worksetCellIndex];
213 const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
222 int basisId = subBasisIdMap[
basis];
225 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
229 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
230 int start = blockOffsets[colIndexer];
231 int end = blockOffsets[colIndexer+1];
237 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
240 if(subJac==Teuchos::null) {
249 jacEpetraBlocks[blockIndex] = subJac;
253 int * rowIndices = 0;
254 double * rowValues = 0;
258 for(
int i=0;i<numEntries;i++)
262 const ScalarT scatterField = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
267 std::vector<double> jacRow(scatterField.size(),0.0);
269 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
270 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
272 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
273 int start = blockOffsets[colIndexer];
274 int end = blockOffsets[colIndexer+1];
279 auto subColIndexer = colIndexers_[colIndexer];
280 const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
285 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
289 if(subJac==Teuchos::null) {
298 jacEpetraBlocks[blockIndex] = subJac;
302 int err = subJac->
ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
304 std::stringstream ss;
305 ss <<
"Failed inserting row: " <<
" (" << lid <<
"): ";
306 for(
int i=0;i<end-start;i++)
307 ss << cLIDs[i] <<
" ";
309 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
311 ss <<
"scatter field = ";
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)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
PHX::MDField< ScalarT > vector
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evaluateFields(typename TRAITS::EvalData d)
ScatterDirichletResidual_BlockedEpetra(const Teuchos::ParameterList &p)
bool isType(const std::string &name) const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
#define TEUCHOS_ASSERT(assertion_test)
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)