43 #ifndef __Panzer_ScatterDirichletResidual_Epetra_Hessian_impl_hpp__ 44 #define __Panzer_ScatterDirichletResidual_Epetra_Hessian_impl_hpp__ 47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT 59 template<
typename TRAITS,
typename LO,
typename GO>
63 const Teuchos::ParameterList& p)
64 : globalIndexer_(indexer)
65 , colGlobalIndexer_(cIndexer)
66 , globalDataKey_(
"Residual Scatter Container")
67 , preserveDiagonal_(false)
69 std::string scatterName = p.get<std::string>(
"Scatter Name");
71 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
74 const std::vector<std::string>& names =
75 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
78 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
80 Teuchos::RCP<PHX::DataLayout> dl =
81 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
83 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
84 local_side_id_ = p.get<
int>(
"Local Side ID");
87 scatterFields_.resize(names.size());
88 for (std::size_t eq = 0; eq < names.size(); ++eq) {
89 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
92 this->addDependentField(scatterFields_[eq]);
95 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
97 applyBC_.resize(names.size());
98 for (std::size_t eq = 0; eq < names.size(); ++eq) {
99 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
100 this->addDependentField(applyBC_[eq]);
105 this->addEvaluatedField(*scatterHolder_);
107 if (p.isType<std::string>(
"Global Data Key"))
108 globalDataKey_ = p.get<std::string>(
"Global Data Key");
110 if (p.isType<
bool>(
"Preserve Diagonal"))
111 preserveDiagonal_ = p.get<
bool>(
"Preserve Diagonal");
113 this->setName(scatterName+
" Scatter Residual (Hessian)");
116 template<
typename TRAITS,
typename LO,
typename GO>
122 fieldIds_.resize(scatterFields_.size());
125 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
127 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
128 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
132 num_nodes = scatterFields_[0].extent(1);
133 num_eq = scatterFields_.size();
136 template<
typename TRAITS,
typename LO,
typename GO>
144 if(epetraContainer_==Teuchos::null) {
146 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
149 dirichletCounter_ = Teuchos::null;
153 Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc->getDataObject(
"Dirichlet Counter");
154 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<
EpetraLinearObjContainer>(dataContainer,
true);
156 dirichletCounter_ = epetraContainer->
get_f();
157 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
161 template<
typename TRAITS,
typename LO,
typename GO>
172 Kokkos::View<const int*, Kokkos::LayoutRight, PHX::Device> cLIDs, rLIDs;
173 std::vector<double> jacRow;
175 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
178 std::string blockId = this->wda(workset).block_id;
179 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
181 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
182 TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
183 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
190 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
191 std::size_t cellLocalId = localCellIds[worksetCellIndex];
193 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
195 cLIDs = colGlobalIndexer_->getElementLIDs(cellLocalId);
200 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
201 int fieldNum = fieldIds_[fieldIndex];
204 const std::pair<std::vector<int>,std::vector<int> > & indicePair
205 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
206 const std::vector<int> & elmtOffset = indicePair.first;
207 const std::vector<int> & basisIdMap = indicePair.second;
210 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
211 int offset = elmtOffset[basis];
212 int row = rLIDs[offset];
216 int basisId = basisIdMap[basis];
219 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
225 int * rowIndices = 0;
226 double * rowValues = 0;
228 Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
230 for(
int i=0;i<numEntries;i++) {
231 if(preserveDiagonal_) {
232 if(row!=rowIndices[i])
241 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
243 if(dirichletCounter_!=Teuchos::null) {
245 (*dirichletCounter_)[row] = 1.0;
249 jacRow.resize(scatterField.size());
250 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
251 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
253 if(!preserveDiagonal_) {
254 int err = Jac->ReplaceMyValues(row, cLIDs.size(),
256 TEUCHOS_ASSERT(err==0);
panzer::Traits::Hessian::ScalarT ScalarT
T * ptrFromStlVector(std::vector< T > &v)
ScatterDirichletResidual_Epetra()
const Teuchos::RCP< Epetra_Vector > get_f() const
Pushes residual values into the residual vector for a Newton-based solve.