43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_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" 60 #include "Phalanx_DataLayout_MDALayout.hpp" 62 #include "Thyra_SpmdVectorBase.hpp" 63 #include "Thyra_ProductVectorBase.hpp" 64 #include "Thyra_BlockedLinearOpBase.hpp" 65 #include "Thyra_get_Epetra_Operator.hpp" 67 #include "Teuchos_FancyOStream.hpp" 69 #include <unordered_map> 71 template <
typename EvalT,
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
76 std::string scatterName = p.
get<std::string>(
"Scatter Name");
81 const std::vector<std::string>& names =
88 for (std::size_t eq = 0; eq < names.size(); ++eq) {
89 PHX::MDField<const ScalarT,Cell,NODE> scatterField = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
92 this->addDependentField(scatterField.fieldTag());
96 this->addEvaluatedField(*scatterHolder);
98 this->setName(scatterName+
" Scatter Residual");
106 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
110 : globalIndexer_(indexer)
111 , globalDataKey_(
"Residual Scatter Container")
113 std::string scatterName = p.
get<std::string>(
"Scatter Name");
118 const std::vector<std::string>& names =
125 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
131 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
132 local_side_id_ = p.
get<
int>(
"Local Side ID");
137 for (std::size_t eq = 0; eq < names.size(); ++eq) {
138 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
144 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
146 applyBC_.resize(names.size());
147 for (std::size_t eq = 0; eq < names.size(); ++eq) {
148 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
149 this->addDependentField(applyBC_[eq]);
154 this->addEvaluatedField(*scatterHolder_);
156 if (p.
isType<std::string>(
"Global Data Key"))
157 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
159 this->setName(scatterName+
" Scatter Residual");
163 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
172 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
173 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
179 this->utils.setFieldData(applyBC_[fd],fm);
187 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
193 = Teuchos::rcp_dynamic_cast<
ContainerType>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
199 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc.getDataObject(globalDataKey_),
true);
204 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
210 using Teuchos::ptrFromRef;
211 using Teuchos::rcp_dynamic_cast;
214 using Thyra::SpmdVectorBase;
221 std::vector<std::pair<int,GO> > GIDs;
222 std::vector<LO> LIDs;
225 std::string blockId = this->wda(workset).block_id;
226 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
229 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f(),
true) :
230 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x(),
true);
238 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
239 std::size_t cellLocalId = localCellIds[worksetCellIndex];
241 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
244 LIDs.resize(GIDs.size());
245 for(std::size_t i=0;i<GIDs.size();i++) {
249 LIDs[i] = r_map->getLocalElement(GIDs[i].second);
255 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
256 int fieldNum = fieldIds_[fieldIndex];
257 int indexerId = globalIndexer_->getFieldBlock(fieldNum);
259 RCP<SpmdVectorBase<double> > dc = rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(indexerId));
260 dc->getNonconstLocalData(ptrFromRef(local_dc));
264 block_r->getNonconstLocalData(ptrFromRef(local_r));
268 const std::pair<std::vector<int>,std::vector<int> > & indicePair
269 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
270 const std::vector<int> & elmtOffset = indicePair.first;
271 const std::vector<int> & basisIdMap = indicePair.second;
280 int basisId = basisIdMap[
basis];
283 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
286 local_r[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
293 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
316 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
320 : globalIndexer_(indexer)
321 , globalDataKey_(
"Residual Scatter Container")
323 std::string scatterName = p.
get<std::string>(
"Scatter Name");
328 const std::vector<std::string>& names =
337 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
338 local_side_id_ = p.
get<
int>(
"Local Side ID");
342 for (std::size_t eq = 0; eq < names.size(); ++eq) {
343 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
349 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
351 applyBC_.resize(names.size());
352 for (std::size_t eq = 0; eq < names.size(); ++eq) {
353 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
354 this->addDependentField(applyBC_[eq]);
359 this->addEvaluatedField(*scatterHolder_);
361 if (p.
isType<std::string>(
"Global Data Key"))
362 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
364 this->setName(scatterName+
" Scatter Residual (Jacobian)");
368 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
377 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
378 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
384 this->utils.setFieldData(applyBC_[fd],fm);
393 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
399 = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
405 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc.getDataObject(globalDataKey_),
true);
410 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
416 using Teuchos::ptrFromRef;
417 using Teuchos::rcp_dynamic_cast;
420 using Thyra::SpmdVectorBase;
424 std::vector<std::pair<int,GO> > GIDs;
425 std::vector<LO> LIDs;
428 std::string blockId = this->wda(workset).block_id;
429 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
434 int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
435 std::vector<int> blockOffsets(numFieldBlocks+1);
436 for(
int blk=0;blk<numFieldBlocks;blk++) {
437 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
438 blockOffsets[blk] = blockOffset;
449 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
450 std::size_t cellLocalId = localCellIds[worksetCellIndex];
452 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
453 blockOffsets[numFieldBlocks] = GIDs.size();
456 LIDs.resize(GIDs.size());
457 for(std::size_t i=0;i<GIDs.size();i++) {
461 LIDs[i] = r_map->getLocalElement(GIDs[i].second);
466 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
467 int fieldNum = fieldIds_[fieldIndex];
468 int blockRowIndex = globalIndexer_->getFieldBlock(fieldNum);
470 RCP<SpmdVectorBase<double> > dc = rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(blockRowIndex));
471 dc->getNonconstLocalData(ptrFromRef(local_dc));
475 block_r->getNonconstLocalData(ptrFromRef(local_r));
478 const std::pair<std::vector<int>,std::vector<int> > & indicePair
479 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
480 const std::vector<int> & elmtOffset = indicePair.first;
481 const std::vector<int> & basisIdMap = indicePair.second;
490 int basisId = basisIdMap[
basis];
493 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
497 for(
int blockColIndex=0;blockColIndex<numFieldBlocks;blockColIndex++) {
498 int start = blockOffsets[blockColIndex];
499 int end = blockOffsets[blockColIndex+1];
505 std::pair<int,int> blockIndex = std::make_pair(blockRowIndex,blockColIndex);
509 if(subJac==Teuchos::null) {
518 jacTpetraBlocks[blockIndex] = subJac;
521 std::size_t sz = subJac->getNumEntriesInLocalRow(lid);
522 std::size_t numEntries = 0;
526 subJac->getLocalRowCopy(lid,rowIndices,rowValues,numEntries);
528 for(std::size_t i=0;i<numEntries;i++)
531 subJac->replaceLocalValues(lid,rowIndices,rowValues);
536 local_r[lid] = scatterField.val();
540 std::vector<double> jacRow(scatterField.size(),0.0);
542 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
543 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
546 for(
int blockColIndex=0;blockColIndex<numFieldBlocks;blockColIndex++) {
547 int start = blockOffsets[blockColIndex];
548 int end = blockOffsets[blockColIndex+1];
554 std::pair<int,int> blockIndex = std::make_pair(blockRowIndex,blockColIndex);
558 if(subJac==Teuchos::null) {
567 jacTpetraBlocks[blockIndex] = subJac;
571 subJac->replaceLocalValues(lid, Teuchos::arrayViewFromVector(LIDs).view(start,end-start),
572 Teuchos::arrayViewFromVector(jacRow).view(start,end-start));
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
Tpetra::CrsMatrix< RealType, LO, GO, NodeT > CrsMatrixType
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evaluateFields(typename TRAITS::EvalData d)
Thyra::TpetraLinearOp< RealType, LO, GO, NodeT > ThyraLinearOp
bool isType(const std::string &name) const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
bool isParameter(const std::string &name) const
Pushes residual values into the residual vector for a Newton-based solve.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>
panzer::Traits::Jacobian::ScalarT ScalarT
ScatterDirichletResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &indexer)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)