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"
62 #include "Phalanx_DataLayout_MDALayout.hpp"
64 #include "Thyra_SpmdVectorBase.hpp"
65 #include "Thyra_ProductVectorBase.hpp"
66 #include "Thyra_DefaultProductVector.hpp"
67 #include "Thyra_BlockedLinearOpBase.hpp"
68 #include "Thyra_get_Epetra_Operator.hpp"
70 #include "Teuchos_FancyOStream.hpp"
72 #include <unordered_map>
79 template<
typename TRAITS,
typename LO,
typename GO>
85 : rowIndexers_(rIndexers)
86 , colIndexers_(cIndexers)
87 , globalDataKey_(
"Residual Scatter Container")
89 std::string scatterName = p.
get<std::string>(
"Scatter Name");
94 const std::vector<std::string>& names =
101 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
107 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
108 local_side_id_ = p.
get<
int>(
"Local Side ID");
112 scatterFields_.resize(names.size());
113 for (std::size_t eq = 0; eq < names.size(); ++eq) {
114 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
117 this->addDependentField(scatterFields_[eq]);
120 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
122 applyBC_.resize(names.size());
123 for (std::size_t eq = 0; eq < names.size(); ++eq) {
124 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
125 this->addDependentField(applyBC_[eq]);
130 this->addEvaluatedField(*scatterHolder_);
132 if (p.
isType<std::string>(
"Global Data Key"))
133 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
135 this->setName(scatterName+
" Scatter Residual");
139 template<
typename TRAITS,
typename LO,
typename GO>
144 indexerIds_.resize(scatterFields_.size());
145 subFieldIds_.resize(scatterFields_.size());
148 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
150 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
153 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
157 num_nodes = scatterFields_[0].extent(1);
161 template<
typename TRAITS,
typename LO,
typename GO>
168 using Teuchos::rcp_dynamic_cast;
173 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
179 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_));
180 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc->getDataObject(globalDataKey_));
183 if(blockedContainer!=Teuchos::null)
185 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
186 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
187 else if(epetraContainer!=Teuchos::null)
189 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
190 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
196 template<
typename TRAITS,
typename LO,
typename GO>
202 using Teuchos::ptrFromRef;
203 using Teuchos::rcp_dynamic_cast;
206 using Thyra::SpmdVectorBase;
210 std::string blockId = this->wda(workset).block_id;
211 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
221 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
222 int rowIndexer = indexerIds_[fieldIndex];
223 int subFieldNum = subFieldIds_[fieldIndex];
225 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
226 ->getNonconstLocalData(ptrFromRef(local_dc));
229 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
230 ->getNonconstLocalData(ptrFromRef(local_r));
232 auto subRowIndexer = rowIndexers_[rowIndexer];
235 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
236 std::size_t cellLocalId = localCellIds[worksetCellIndex];
238 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
242 const std::pair<std::vector<int>,std::vector<int> > & indicePair
243 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
244 const std::vector<int> & elmtOffset = indicePair.first;
245 const std::vector<int> & basisIdMap = indicePair.second;
248 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
249 int offset = elmtOffset[basis];
250 int lid = LIDs[offset];
254 int basisId = basisIdMap[basis];
257 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
260 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
267 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
270 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
271 int offset = elmtOffset[basis];
272 int lid = LIDs[offset];
276 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
291 template<
typename TRAITS,
typename LO,
typename GO>
297 : rowIndexers_(rIndexers)
298 , colIndexers_(cIndexers)
299 , globalDataKey_(
"Residual Scatter Container")
301 std::string scatterName = p.
get<std::string>(
"Scatter Name");
306 const std::vector<std::string>& names =
313 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
319 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
320 local_side_id_ = p.
get<
int>(
"Local Side ID");
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]);
332 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
334 applyBC_.resize(names.size());
335 for (std::size_t eq = 0; eq < names.size(); ++eq) {
336 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
337 this->addDependentField(applyBC_[eq]);
342 this->addEvaluatedField(*scatterHolder_);
344 if (p.
isType<std::string>(
"Global Data Key"))
345 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
347 this->setName(scatterName+
" Scatter Tangent");
351 template<
typename TRAITS,
typename LO,
typename GO>
356 indexerIds_.resize(scatterFields_.size());
357 subFieldIds_.resize(scatterFields_.size());
360 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
362 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
365 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
369 num_nodes = scatterFields_[0].extent(1);
373 template<
typename TRAITS,
typename LO,
typename GO>
380 using Teuchos::rcp_dynamic_cast;
385 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
391 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_));
392 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc->getDataObject(globalDataKey_));
395 if(blockedContainer!=Teuchos::null)
397 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
398 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
399 else if(epetraContainer!=Teuchos::null)
401 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
402 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
408 template<
typename TRAITS,
typename LO,
typename GO>
416 using Teuchos::ptrFromRef;
417 using Teuchos::rcp_dynamic_cast;
420 using Thyra::SpmdVectorBase;
424 std::string blockId = this->wda(workset).block_id;
425 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
435 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
436 int rowIndexer = indexerIds_[fieldIndex];
437 int subFieldNum = subFieldIds_[fieldIndex];
439 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
440 ->getNonconstLocalData(ptrFromRef(local_dc));
443 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
444 ->getNonconstLocalData(ptrFromRef(local_r));
446 auto subRowIndexer = rowIndexers_[rowIndexer];
449 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
450 std::size_t cellLocalId = localCellIds[worksetCellIndex];
452 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
456 const std::pair<std::vector<int>,std::vector<int> > & indicePair
457 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
458 const std::vector<int> & elmtOffset = indicePair.first;
459 const std::vector<int> & basisIdMap = indicePair.second;
462 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
463 int offset = elmtOffset[basis];
464 int lid = LIDs[offset];
468 int basisId = basisIdMap[basis];
471 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
474 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
481 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
484 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
485 int offset = elmtOffset[basis];
486 int lid = LIDs[offset];
490 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
504 template<
typename TRAITS,
typename LO,
typename GO>
510 : rowIndexers_(rIndexers)
511 , colIndexers_(cIndexers)
512 , globalDataKey_(
"Residual Scatter Container")
514 std::string scatterName = p.
get<std::string>(
"Scatter Name");
519 const std::vector<std::string>& names =
528 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
529 local_side_id_ = p.
get<
int>(
"Local Side ID");
532 scatterFields_.resize(names.size());
533 for (std::size_t eq = 0; eq < names.size(); ++eq) {
534 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
537 this->addDependentField(scatterFields_[eq]);
540 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
542 applyBC_.resize(names.size());
543 for (std::size_t eq = 0; eq < names.size(); ++eq) {
544 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
545 this->addDependentField(applyBC_[eq]);
550 this->addEvaluatedField(*scatterHolder_);
552 if (p.
isType<std::string>(
"Global Data Key"))
553 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
555 if(colIndexers_.size()==0)
556 colIndexers_ = rowIndexers_;
558 this->setName(scatterName+
" Scatter Residual (Jacobian)");
562 template<
typename TRAITS,
typename LO,
typename GO>
567 indexerIds_.resize(scatterFields_.size());
568 subFieldIds_.resize(scatterFields_.size());
571 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
573 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
576 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
580 num_nodes = scatterFields_[0].extent(1);
581 num_eq = scatterFields_.size();
585 template<
typename TRAITS,
typename LO,
typename GO>
591 using Teuchos::rcp_dynamic_cast;
595 = rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
601 blockContainer = rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_),
true);
609 template<
typename TRAITS,
typename LO,
typename GO>
615 using Teuchos::ptrFromRef;
616 using Teuchos::rcp_dynamic_cast;
618 using Thyra::SpmdVectorBase;
621 std::string blockId = this->wda(workset).block_id;
622 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
624 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
626 std::vector<int> blockOffsets;
636 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
637 int rowIndexer = indexerIds_[fieldIndex];
638 int subFieldNum = subFieldIds_[fieldIndex];
642 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
643 ->getNonconstLocalData(ptrFromRef(local_dc));
647 if(r_!=Teuchos::null)
648 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
649 ->getNonconstLocalData(ptrFromRef(local_r));
651 auto subRowIndexer = rowIndexers_[rowIndexer];
652 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
653 const std::vector<int> & subElmtOffset = subIndicePair.first;
654 const std::vector<int> & subBasisIdMap = subIndicePair.second;
657 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
658 std::size_t cellLocalId = localCellIds[worksetCellIndex];
660 auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
663 for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
664 int offset = subElmtOffset[basis];
665 int lid = rLIDs[offset];
669 int basisId = subBasisIdMap[basis];
672 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
676 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
677 int start = blockOffsets[colIndexer];
678 int end = blockOffsets[colIndexer+1];
684 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
685 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
687 if(subJac==Teuchos::null) {
696 jacEpetraBlocks[blockIndex] = subJac;
700 int * rowIndices = 0;
701 double * rowValues = 0;
705 for(
int i=0;i<numEntries;i++)
709 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
711 if(r_!=Teuchos::null)
712 local_r[lid] = scatterField.val();
716 std::vector<double> jacRow(scatterField.size(),0.0);
718 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
719 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
721 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
722 int start = blockOffsets[colIndexer];
723 int end = blockOffsets[colIndexer+1];
728 auto subColIndexer = colIndexers_[colIndexer];
729 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
734 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
735 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
738 if(subJac==Teuchos::null) {
747 jacEpetraBlocks[blockIndex] = subJac;
751 int err = subJac->
ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
753 std::stringstream ss;
754 ss <<
"Failed inserting row: " <<
" (" << lid <<
"): ";
755 for(
int i=0;i<end-start;i++)
756 ss << cLIDs[i] <<
" ";
758 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
760 ss <<
"scatter field = ";
761 scatterFields_[fieldIndex].print(ss);
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis, std::vector< int > &blockOffsets)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Pushes residual values into the residual vector for a Newton-based solve.
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evaluateFields(typename TRAITS::EvalData d)
ScatterDirichletResidual_BlockedEpetra(const Teuchos::ParameterList &p)
panzer::Traits::Jacobian::ScalarT ScalarT
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)