43 #ifndef PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP 44 #define PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 51 #include "Epetra_Map.h" 52 #include "Epetra_Vector.h" 53 #include "Epetra_CrsMatrix.h" 62 #include "Phalanx_DataLayout_MDALayout.hpp" 64 #include "Teuchos_FancyOStream.hpp" 70 template<
typename TRAITS,
typename LO,
typename GO>
75 bool useDiscreteAdjoint)
76 : globalIndexer_(indexer)
77 , globalDataKey_(
"Residual Scatter Container")
78 , useDiscreteAdjoint_(useDiscreteAdjoint)
80 std::string scatterName = p.
get<std::string>(
"Scatter Name");
85 const std::vector<std::string>& names =
96 for (std::size_t eq = 0; eq < names.size(); ++eq) {
97 scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
104 this->addEvaluatedField(*scatterHolder_);
106 if (p.
isType<std::string>(
"Global Data Key"))
107 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
109 this->setName(scatterName+
" Scatter Residual");
113 template<
typename TRAITS,
typename LO,
typename GO>
122 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
123 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
131 template<
typename TRAITS,
typename LO,
typename GO>
138 if(epetraContainer_==Teuchos::null) {
146 template<
typename TRAITS,
typename LO,
typename GO>
150 std::vector<int> LIDs;
153 std::string blockId = this->wda(workset).block_id;
154 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
164 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
165 std::size_t cellLocalId = localCellIds[worksetCellIndex];
167 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
170 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
171 int fieldNum = fieldIds_[fieldIndex];
172 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
188 template<
typename TRAITS,
typename LO,
typename GO>
193 bool useDiscreteAdjoint)
194 : globalIndexer_(indexer)
196 std::string scatterName = p.
get<std::string>(
"Scatter Name");
201 const std::vector<std::string>& names =
212 for (std::size_t eq = 0; eq < names.size(); ++eq) {
213 scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
220 this->addEvaluatedField(*scatterHolder_);
222 this->setName(scatterName+
" Scatter Tangent");
226 template<
typename TRAITS,
typename LO,
typename GO>
235 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
236 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
244 template<
typename TRAITS,
typename LO,
typename GO>
249 using Teuchos::rcp_dynamic_cast;
252 std::vector<std::string> activeParameters =
255 dfdp_vectors_.clear();
256 for(std::size_t i=0;i<activeParameters.size();i++) {
258 dfdp_vectors_.push_back(vec);
263 template<
typename TRAITS,
typename LO,
typename GO>
267 std::vector<int> LIDs;
270 std::string blockId = this->wda(workset).block_id;
271 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
279 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
280 std::size_t cellLocalId = localCellIds[worksetCellIndex];
282 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
285 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
286 int fieldNum = fieldIds_[fieldIndex];
287 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
296 for(
int d=0;d<
value.size();d++)
297 (*dfdp_vectors_[d])[lid] +=
value.fastAccessDx(d);
307 template<
typename TRAITS,
typename LO,
typename GO>
312 bool useDiscreteAdjoint)
313 : globalIndexer_(indexer)
314 , colGlobalIndexer_(cIndexer)
315 , globalDataKey_(
"Residual Scatter Container")
316 , useDiscreteAdjoint_(useDiscreteAdjoint)
318 std::string scatterName = p.
get<std::string>(
"Scatter Name");
323 const std::vector<std::string>& names =
334 for (std::size_t eq = 0; eq < names.size(); ++eq) {
335 scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
342 this->addEvaluatedField(*scatterHolder_);
344 if (p.
isType<std::string>(
"Global Data Key"))
345 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
346 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
347 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
350 if(useDiscreteAdjoint)
353 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
357 template<
typename TRAITS,
typename LO,
typename GO>
368 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
369 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
377 template<
typename TRAITS,
typename LO,
typename GO>
384 if(epetraContainer_==Teuchos::null) {
392 template<
typename TRAITS,
typename LO,
typename GO>
396 std::vector<int> cLIDs, rLIDs;
397 std::vector<double> jacRow;
399 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
402 std::string blockId = this->wda(workset).block_id;
403 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
409 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
417 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
418 std::size_t cellLocalId = localCellIds[worksetCellIndex];
420 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
421 cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
423 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
424 const std::vector<int> other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
425 cLIDs.insert(cLIDs.end(), other_cLIDs.begin(), other_cLIDs.end());
429 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
430 int fieldNum = fieldIds_[fieldIndex];
431 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
434 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
436 int rowOffset = elmtOffset[rowBasisNum];
437 int row = rLIDs[rowOffset];
441 r->SumIntoMyValue(row,0,scatterField.val());
444 if(useDiscreteAdjoint_) {
446 jacRow.resize(scatterField.size());
448 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
449 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
451 for(std::size_t c=0;c<cLIDs.size();c++) {
452 int err = Jac->SumIntoMyValues(cLIDs[c], 1, &jacRow[c],&row);
457 int err = Jac->SumIntoMyValues(
459 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)
panzer::Traits::Tangent::ScalarT ScalarT
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
panzer::Traits::Jacobian::ScalarT ScalarT
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)