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"
63 #include "Phalanx_DataLayout_MDALayout.hpp"
65 #include "Teuchos_FancyOStream.hpp"
71 template<
typename TRAITS,
typename LO,
typename GO>
76 bool useDiscreteAdjoint)
77 : globalIndexer_(indexer)
78 , globalDataKey_(
"Residual Scatter Container")
79 , useDiscreteAdjoint_(useDiscreteAdjoint)
81 std::string scatterName = p.
get<std::string>(
"Scatter Name");
86 const std::vector<std::string>& names =
96 scatterFields_.resize(names.size());
97 for (std::size_t eq = 0; eq < names.size(); ++eq) {
98 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
101 this->addDependentField(scatterFields_[eq]);
105 this->addEvaluatedField(*scatterHolder_);
107 if (p.
isType<std::string>(
"Global Data Key"))
108 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
110 this->setName(scatterName+
" Scatter Residual");
114 template<
typename TRAITS,
typename LO,
typename GO>
119 fieldIds_.resize(scatterFields_.size());
121 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
123 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
124 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
129 template<
typename TRAITS,
typename LO,
typename GO>
136 if(epetraContainer_==Teuchos::null) {
144 template<
typename TRAITS,
typename LO,
typename GO>
149 std::string blockId = this->wda(workset).block_id;
150 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
160 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
161 std::size_t cellLocalId = localCellIds[worksetCellIndex];
163 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
166 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
167 int fieldNum = fieldIds_[fieldIndex];
168 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
171 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
172 int offset = elmtOffset[basis];
173 int lid = LIDs[offset];
174 (*r)[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis);
184 template<
typename TRAITS,
typename LO,
typename GO>
190 : globalIndexer_(indexer)
192 std::string scatterName = p.
get<std::string>(
"Scatter Name");
197 const std::vector<std::string>& names =
207 scatterFields_.resize(names.size());
208 for (std::size_t eq = 0; eq < names.size(); ++eq) {
209 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
212 this->addDependentField(scatterFields_[eq]);
216 this->addEvaluatedField(*scatterHolder_);
218 this->setName(scatterName+
" Scatter Tangent");
222 template<
typename TRAITS,
typename LO,
typename GO>
227 fieldIds_.resize(scatterFields_.size());
229 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
231 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
232 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
237 template<
typename TRAITS,
typename LO,
typename GO>
242 using Teuchos::rcp_dynamic_cast;
245 std::vector<std::string> activeParameters =
248 dfdp_vectors_.clear();
249 for(std::size_t i=0;i<activeParameters.size();i++) {
251 dfdp_vectors_.push_back(vec);
256 template<
typename TRAITS,
typename LO,
typename GO>
261 std::string blockId = this->wda(workset).block_id;
262 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
270 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
271 std::size_t cellLocalId = localCellIds[worksetCellIndex];
273 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
276 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
277 int fieldNum = fieldIds_[fieldIndex];
278 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
281 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
282 int offset = elmtOffset[basis];
283 int lid = LIDs[offset];
286 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
287 for(
int d=0;d<value.size();d++)
288 (*dfdp_vectors_[d])[lid] += value.fastAccessDx(d);
298 template<
typename TRAITS,
typename LO,
typename GO>
303 bool useDiscreteAdjoint)
304 : globalIndexer_(indexer)
305 , colGlobalIndexer_(cIndexer)
306 , globalDataKey_(
"Residual Scatter Container")
307 , useDiscreteAdjoint_(useDiscreteAdjoint)
309 std::string scatterName = p.
get<std::string>(
"Scatter Name");
314 const std::vector<std::string>& names =
324 scatterFields_.resize(names.size());
325 for (std::size_t eq = 0; eq < names.size(); ++eq) {
326 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
329 this->addDependentField(scatterFields_[eq]);
333 this->addEvaluatedField(*scatterHolder_);
335 if (p.
isType<std::string>(
"Global Data Key"))
336 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
337 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
338 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
341 if(useDiscreteAdjoint)
344 this->setName(scatterName+
" Scatter Residual Epetra (Jacobian)");
348 template<
typename TRAITS,
typename LO,
typename GO>
355 fieldIds_.resize(scatterFields_.size());
357 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
359 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
360 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
365 template<
typename TRAITS,
typename LO,
typename GO>
372 if(epetraContainer_==Teuchos::null) {
380 template<
typename TRAITS,
typename LO,
typename GO>
384 std::vector<double> jacRow;
386 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
389 std::string blockId = this->wda(workset).block_id;
390 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
396 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
404 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
405 std::size_t cellLocalId = localCellIds[worksetCellIndex];
407 auto rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
408 auto initial_cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
409 std::vector<int> cLIDs;
410 for (
int i(0); i < static_cast<int>(initial_cLIDs.extent(0)); ++i)
411 cLIDs.push_back(initial_cLIDs(i));
413 const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
414 auto other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
415 for (
int i(0); i < static_cast<int>(other_cLIDs.extent(0)); ++i)
416 cLIDs.push_back(other_cLIDs(i));
420 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
421 int fieldNum = fieldIds_[fieldIndex];
422 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
425 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
426 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
427 int rowOffset = elmtOffset[rowBasisNum];
428 int row = rLIDs[rowOffset];
432 r->SumIntoMyValue(row,0,scatterField.val());
435 if(useDiscreteAdjoint_) {
437 jacRow.resize(scatterField.size());
439 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
440 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
442 for(std::size_t c=0;c<cLIDs.size();c++) {
450 std::min(cLIDs.size(),
static_cast<size_t>(scatterField.size())),
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.
panzer::Traits::Jacobian::ScalarT ScalarT
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)