43 #ifndef __Panzer_ScatterResidual_Epetra_Hessian_impl_hpp__ 44 #define __Panzer_ScatterResidual_Epetra_Hessian_impl_hpp__ 47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT 50 #include "Teuchos_Assert.hpp" 52 #include "Phalanx_DataLayout.hpp" 54 #include "Epetra_Map.h" 55 #include "Epetra_Vector.h" 56 #include "Epetra_CrsMatrix.h" 66 #include "Phalanx_DataLayout_MDALayout.hpp" 68 #include "Teuchos_FancyOStream.hpp" 74 template<
typename TRAITS,
typename LO,
typename GO>
79 bool useDiscreteAdjoint)
80 : globalIndexer_(indexer)
81 , colGlobalIndexer_(cIndexer)
82 , globalDataKey_(
"Residual Scatter Container")
83 , useDiscreteAdjoint_(useDiscreteAdjoint)
85 std::string scatterName = p.
get<std::string>(
"Scatter Name");
90 const std::vector<std::string>& names =
101 for (std::size_t eq = 0; eq < names.size(); ++eq) {
102 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
109 this->addEvaluatedField(*scatterHolder_);
111 if (p.
isType<std::string>(
"Global Data Key"))
112 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
113 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
114 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
117 if(useDiscreteAdjoint)
120 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
124 template<
typename TRAITS,
typename LO,
typename GO>
133 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
134 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
142 template<
typename TRAITS,
typename LO,
typename GO>
147 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc.getDataObject(globalDataKey_));
149 if(epetraContainer_==Teuchos::null) {
152 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
157 template<
typename TRAITS,
typename LO,
typename GO>
161 std::vector<int> cLIDs, rLIDs;
162 std::vector<double> jacRow;
164 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
167 std::string blockId = this->wda(workset).block_id;
168 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
174 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
182 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
183 std::size_t cellLocalId = localCellIds[worksetCellIndex];
185 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
186 cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
188 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
189 const std::vector<int> other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
190 cLIDs.insert(cLIDs.end(), other_cLIDs.begin(), other_cLIDs.end());
194 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
195 int fieldNum = fieldIds_[fieldIndex];
196 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
199 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
200 const ScalarT scatterField = (
scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
201 int rowOffset = elmtOffset[rowBasisNum];
202 int row = rLIDs[rowOffset];
205 jacRow.resize(scatterField.size());
207 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
208 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
211 int err = Jac->SumIntoMyValues(
213 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T & get(const std::string &name, T def_value)
bool nonnull(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)
Pushes residual values into the residual vector for a Newton-based solve.
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)