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,
64 : globalIndexer_(indexer)
65 , colGlobalIndexer_(cIndexer)
66 , globalDataKey_(
"Residual Scatter Container")
67 , preserveDiagonal_(false)
69 std::string scatterName = p.
get<std::string>(
"Scatter Name");
74 const std::vector<std::string>& names =
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) {
152 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc,
true);
154 dirichletCounter_ = Teuchos::null;
161 dirichletCounter_ = epetraContainer->
get_f();
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;
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(),
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T * ptrFromStlVector(std::vector< T > &v)
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)
const Teuchos::RCP< Epetra_Vector > get_f() const