43 #ifndef PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP 44 #define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 58 #include "Thyra_SpmdVectorBase.hpp" 59 #include "Thyra_ProductVectorBase.hpp" 60 #include "Thyra_BlockedLinearOpBase.hpp" 62 #include "Phalanx_DataLayout_MDALayout.hpp" 64 #include "Teuchos_FancyOStream.hpp" 66 template <
typename EvalT,
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
71 std::string scatterName = p.
get<std::string>(
"Scatter Name");
76 const std::vector<std::string>& names =
83 for (std::size_t eq = 0; eq < names.size(); ++eq) {
84 PHX::MDField<const ScalarT,Cell,NODE>
field = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
87 this->addDependentField(
field.fieldTag());
91 this->addEvaluatedField(*scatterHolder);
93 this->setName(scatterName+
" Scatter Residual");
100 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
104 : globalIndexer_(indexer)
105 , globalDataKey_(
"Residual Scatter Container")
107 std::string scatterName = p.
get<std::string>(
"Scatter Name");
112 const std::vector<std::string>& names =
123 for (std::size_t eq = 0; eq < names.size(); ++eq) {
124 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
131 this->addEvaluatedField(*scatterHolder_);
133 if (p.
isType<std::string>(
"Global Data Key"))
134 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
136 this->setName(scatterName+
" Scatter Residual");
140 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
150 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
151 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
159 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
164 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc.getDataObject(globalDataKey_),
true);
168 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
174 using Teuchos::ptrFromRef;
175 using Teuchos::rcp_dynamic_cast;
178 using Thyra::SpmdVectorBase;
181 std::vector<std::pair<int,GO> > GIDs;
182 std::vector<LO> LIDs;
185 std::string blockId = this->wda(workset).block_id;
186 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
198 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
199 std::size_t cellLocalId = localCellIds[worksetCellIndex];
201 globalIndexer_->getElementGIDs(cellLocalId,GIDs,blockId);
204 LIDs.resize(GIDs.size());
205 for(std::size_t i=0;i<GIDs.size();i++) {
209 LIDs[i] = r_map->getLocalElement(GIDs[i].second);
214 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
215 int fieldNum = fieldIds_[fieldIndex];
216 int indexerId = globalIndexer_->getFieldBlock(fieldNum);
220 block_r->getNonconstLocalData(ptrFromRef(local_r));
222 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
238 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
242 : globalIndexer_(indexer)
243 , globalDataKey_(
"Residual Scatter Container")
245 std::string scatterName = p.
get<std::string>(
"Scatter Name");
250 const std::vector<std::string>& names =
261 for (std::size_t eq = 0; eq < names.size(); ++eq) {
262 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
269 this->addEvaluatedField(*scatterHolder_);
271 if (p.
isType<std::string>(
"Global Data Key"))
272 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
274 this->setName(scatterName+
" Scatter Residual (Jacobian)");
278 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
288 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
289 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
297 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
302 using Teuchos::rcp_dynamic_cast;
305 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(d.gedc.getDataObject(globalDataKey_));
307 if(blockedContainer_==Teuchos::null) {
314 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
320 using Teuchos::ptrFromRef;
321 using Teuchos::rcp_dynamic_cast;
324 using Thyra::SpmdVectorBase;
328 std::vector<std::pair<int,GO> > GIDs;
329 std::vector<LO> LIDs;
330 std::vector<double> jacRow;
333 std::string blockId = this->wda(workset).block_id;
334 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
341 int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
342 std::vector<int> blockOffsets(numFieldBlocks+1);
343 for(
int blk=0;blk<numFieldBlocks;blk++) {
344 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
345 blockOffsets[blk] = blockOffset;
356 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
357 std::size_t cellLocalId = localCellIds[worksetCellIndex];
359 globalIndexer_->getElementGIDs(cellLocalId,GIDs,blockId);
362 LIDs.resize(GIDs.size());
363 for(std::size_t i=0;i<GIDs.size();i++) {
367 LIDs[i] = r_map->getLocalElement(GIDs[i].second);
372 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
373 int fieldNum = fieldIds_[fieldIndex];
374 int blockRowIndex = globalIndexer_->getFieldBlock(fieldNum);
377 if(r!=Teuchos::null) {
379 block_r->getNonconstLocalData(ptrFromRef(local_r));
382 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
385 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
387 int rowOffset = elmtOffset[rowBasisNum];
388 int r_lid = LIDs[rowOffset];
391 if(local_r!=Teuchos::null)
392 local_r[r_lid] += (scatterField.val());
394 blockOffsets[numFieldBlocks] = scatterField.
size();
397 jacRow.resize(scatterField.size());
399 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex) {
400 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
403 for(
int blockColIndex=0;blockColIndex<numFieldBlocks;blockColIndex++) {
404 int start = blockOffsets[blockColIndex];
405 int end = blockOffsets[blockColIndex+1];
411 std::pair<int,int> blockIndex = std::make_pair(blockRowIndex,blockColIndex);
415 if(subJac==Teuchos::null) {
424 jacTpetraBlocks[blockIndex] = subJac;
428 subJac->sumIntoLocalValues(r_lid, Teuchos::arrayViewFromVector(LIDs).view(start,end-start),
429 Teuchos::arrayViewFromVector(jacRow).view(start,end-start));
Pushes residual values into the residual vector for a Newton-based solve.
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
ScatterResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &indexer)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isType(const std::string &name) const
PHX::MDField< const ScalarT, Cell, IP > field
void evaluateFields(typename TRAITS::EvalData d)
Teuchos::RCP< LinearObjContainer > getGhostedLOC() const
Thyra::TpetraLinearOp< RealType, LO, GO, NodeT > ThyraLinearOp
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
Tpetra::CrsMatrix< RealType, LO, GO, NodeT > CrsMatrixType
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
panzer::Traits::Jacobian::ScalarT ScalarT