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>
60 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
61 ScatterDirichletResidual_Epetra(
const Teuchos::RCP<
const UniqueGlobalIndexer<LO,GO> > & indexer,
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");
88 for (std::size_t eq = 0; eq < names.size(); ++eq) {
89 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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>
118 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
119 postRegistrationSetup(
typename TRAITS::SetupData d,
127 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
128 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
133 this->utils.setFieldData(applyBC_[fd],fm);
141 template<
typename TRAITS,
typename LO,
typename GO>
143 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
144 preEvaluate(
typename TRAITS::PreEvalData d)
147 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc.getDataObject(globalDataKey_));
149 if(epetraContainer_==Teuchos::null) {
151 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
152 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc,
true);
154 dirichletCounter_ = Teuchos::null;
158 Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc.getDataObject(
"Dirichlet Counter");
159 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(dataContainer,
true);
161 dirichletCounter_ = epetraContainer->get_f();
162 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
166 template<
typename TRAITS,
typename LO,
typename GO>
168 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
169 evaluateFields(
typename TRAITS::EvalData workset)
174 std::vector<int> cLIDs, rLIDs;
175 std::vector<double> jacRow;
177 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
180 std::string blockId = this->wda(workset).block_id;
181 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
183 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
184 TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
185 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
192 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
193 std::size_t cellLocalId = localCellIds[worksetCellIndex];
195 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
197 cLIDs = colGlobalIndexer_->getElementLIDs(cellLocalId);
202 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
203 int fieldNum = fieldIds_[fieldIndex];
206 const std::pair<std::vector<int>,std::vector<int> > & indicePair
207 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
208 const std::vector<int> & elmtOffset = indicePair.first;
209 const std::vector<int> & basisIdMap = indicePair.second;
218 int basisId = basisIdMap[
basis];
221 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
227 int * rowIndices = 0;
228 double * rowValues = 0;
230 Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
232 for(
int i=0;i<numEntries;i++) {
233 if(preserveDiagonal_) {
234 if(row!=rowIndices[i])
243 const ScalarT scatterField = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
245 if(dirichletCounter_!=Teuchos::null) {
247 (*dirichletCounter_)[row] = 1.0;
251 jacRow.resize(scatterField.size());
252 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
253 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
255 if(!preserveDiagonal_) {
256 int err = Jac->ReplaceMyValues(row, cLIDs.size(),
259 TEUCHOS_ASSERT(err==0);
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T * ptrFromStlVector(std::vector< T > &v)
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.