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" 61 #include "Phalanx_DataLayout_MDALayout.hpp" 63 #include "Thyra_SpmdVectorBase.hpp" 64 #include "Thyra_ProductVectorBase.hpp" 65 #include "Thyra_DefaultProductVector.hpp" 66 #include "Thyra_BlockedLinearOpBase.hpp" 67 #include "Thyra_get_Epetra_Operator.hpp" 69 #include "Teuchos_FancyOStream.hpp" 71 #include <unordered_map> 78 template<
typename TRAITS,
typename LO,
typename GO>
83 bool useDiscreteAdjoint)
84 : rowIndexers_(rIndexers)
85 , colIndexers_(cIndexers)
86 , globalDataKey_(
"Residual Scatter Container")
88 std::string scatterName = p.
get<std::string>(
"Scatter Name");
93 const std::vector<std::string>& names =
100 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
106 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
107 local_side_id_ = p.
get<
int>(
"Local Side ID");
112 for (std::size_t eq = 0; eq < names.size(); ++eq) {
113 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
119 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
121 applyBC_.resize(names.size());
122 for (std::size_t eq = 0; eq < names.size(); ++eq) {
123 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
124 this->addDependentField(applyBC_[eq]);
129 this->addEvaluatedField(*scatterHolder_);
131 if (p.
isType<std::string>(
"Global Data Key"))
132 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
134 this->setName(scatterName+
" Scatter Residual");
138 template<
typename TRAITS,
typename LO,
typename GO>
149 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
152 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
158 this->utils.setFieldData(applyBC_[fd],fm);
166 template<
typename TRAITS,
typename LO,
typename GO>
173 using Teuchos::rcp_dynamic_cast;
178 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
184 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
185 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
188 if(blockedContainer!=Teuchos::null)
190 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
191 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
192 else if(epetraContainer!=Teuchos::null)
194 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
195 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
201 template<
typename TRAITS,
typename LO,
typename GO>
207 using Teuchos::ptrFromRef;
208 using Teuchos::rcp_dynamic_cast;
211 using Thyra::SpmdVectorBase;
215 std::string blockId = this->wda(workset).block_id;
216 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
226 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
227 int rowIndexer = indexerIds_[fieldIndex];
228 int subFieldNum = subFieldIds_[fieldIndex];
230 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
231 ->getNonconstLocalData(ptrFromRef(local_dc));
234 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
235 ->getNonconstLocalData(ptrFromRef(local_r));
237 auto subRowIndexer = rowIndexers_[rowIndexer];
240 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
241 std::size_t cellLocalId = localCellIds[worksetCellIndex];
243 const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
247 const std::pair<std::vector<int>,std::vector<int> > & indicePair
248 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
249 const std::vector<int> & elmtOffset = indicePair.first;
250 const std::vector<int> & basisIdMap = indicePair.second;
259 int basisId = basisIdMap[
basis];
262 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
265 local_r[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
272 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
296 template<
typename TRAITS,
typename LO,
typename GO>
301 bool useDiscreteAdjoint)
302 : rowIndexers_(rIndexers)
303 , colIndexers_(cIndexers)
304 , globalDataKey_(
"Residual Scatter Container")
306 std::string scatterName = p.
get<std::string>(
"Scatter Name");
311 const std::vector<std::string>& names =
318 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
324 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
325 local_side_id_ = p.
get<
int>(
"Local Side ID");
330 for (std::size_t eq = 0; eq < names.size(); ++eq) {
331 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
337 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
339 applyBC_.resize(names.size());
340 for (std::size_t eq = 0; eq < names.size(); ++eq) {
341 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
342 this->addDependentField(applyBC_[eq]);
347 this->addEvaluatedField(*scatterHolder_);
349 if (p.
isType<std::string>(
"Global Data Key"))
350 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
352 this->setName(scatterName+
" Scatter Tangent");
356 template<
typename TRAITS,
typename LO,
typename GO>
367 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
370 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
376 this->utils.setFieldData(applyBC_[fd],fm);
384 template<
typename TRAITS,
typename LO,
typename GO>
391 using Teuchos::rcp_dynamic_cast;
396 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
402 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
403 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
406 if(blockedContainer!=Teuchos::null)
408 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
409 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
410 else if(epetraContainer!=Teuchos::null)
412 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
413 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
419 template<
typename TRAITS,
typename LO,
typename GO>
427 using Teuchos::ptrFromRef;
428 using Teuchos::rcp_dynamic_cast;
431 using Thyra::SpmdVectorBase;
435 std::string blockId = this->wda(workset).block_id;
436 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
446 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
447 int rowIndexer = indexerIds_[fieldIndex];
448 int subFieldNum = subFieldIds_[fieldIndex];
450 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
451 ->getNonconstLocalData(ptrFromRef(local_dc));
454 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
455 ->getNonconstLocalData(ptrFromRef(local_r));
457 auto subRowIndexer = rowIndexers_[rowIndexer];
460 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
461 std::size_t cellLocalId = localCellIds[worksetCellIndex];
463 const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
467 const std::pair<std::vector<int>,std::vector<int> > & indicePair
468 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
469 const std::vector<int> & elmtOffset = indicePair.first;
470 const std::vector<int> & basisIdMap = indicePair.second;
479 int basisId = basisIdMap[
basis];
482 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
485 local_r[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
492 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
515 template<
typename TRAITS,
typename LO,
typename GO>
520 bool useDiscreteAdjoint)
521 : rowIndexers_(rIndexers)
522 , colIndexers_(cIndexers)
523 , globalDataKey_(
"Residual Scatter Container")
525 std::string scatterName = p.
get<std::string>(
"Scatter Name");
530 const std::vector<std::string>& names =
539 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
540 local_side_id_ = p.
get<
int>(
"Local Side ID");
544 for (std::size_t eq = 0; eq < names.size(); ++eq) {
545 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
551 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
553 applyBC_.resize(names.size());
554 for (std::size_t eq = 0; eq < names.size(); ++eq) {
555 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
556 this->addDependentField(applyBC_[eq]);
561 this->addEvaluatedField(*scatterHolder_);
563 if (p.
isType<std::string>(
"Global Data Key"))
564 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
566 if(colIndexers_.size()==0)
567 colIndexers_ = rowIndexers_;
569 this->setName(scatterName+
" Scatter Residual (Jacobian)");
573 template<
typename TRAITS,
typename LO,
typename GO>
584 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
587 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
593 this->utils.setFieldData(applyBC_[fd],fm);
602 template<
typename TRAITS,
typename LO,
typename GO>
608 using Teuchos::rcp_dynamic_cast;
612 = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
618 blockContainer = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_),
true);
626 template<
typename TRAITS,
typename LO,
typename GO>
632 using Teuchos::ptrFromRef;
633 using Teuchos::rcp_dynamic_cast;
635 using Thyra::SpmdVectorBase;
638 std::string blockId = this->wda(workset).block_id;
639 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
641 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
643 std::vector<int> blockOffsets;
653 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
654 int rowIndexer = indexerIds_[fieldIndex];
655 int subFieldNum = subFieldIds_[fieldIndex];
659 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
660 ->getNonconstLocalData(ptrFromRef(local_dc));
664 if(r_!=Teuchos::null)
665 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
666 ->getNonconstLocalData(ptrFromRef(local_r));
668 auto subRowIndexer = rowIndexers_[rowIndexer];
669 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
670 const std::vector<int> & subElmtOffset = subIndicePair.first;
671 const std::vector<int> & subBasisIdMap = subIndicePair.second;
674 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
675 std::size_t cellLocalId = localCellIds[worksetCellIndex];
677 const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
686 int basisId = subBasisIdMap[
basis];
689 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
693 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
694 int start = blockOffsets[colIndexer];
695 int end = blockOffsets[colIndexer+1];
701 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
704 if(subJac==Teuchos::null) {
712 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
713 jacEpetraBlocks[blockIndex] = subJac;
717 int * rowIndices = 0;
718 double * rowValues = 0;
720 subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
722 for(
int i=0;i<numEntries;i++)
728 if(r_!=Teuchos::null)
729 local_r[lid] = scatterField.val();
733 std::vector<double> jacRow(scatterField.size(),0.0);
735 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
736 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
738 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
739 int start = blockOffsets[colIndexer];
740 int end = blockOffsets[colIndexer+1];
745 auto subColIndexer = colIndexers_[colIndexer];
746 const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
751 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
755 if(subJac==Teuchos::null) {
763 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
764 jacEpetraBlocks[blockIndex] = subJac;
768 int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
770 std::stringstream ss;
771 ss <<
"Failed inserting row: " <<
" (" << lid <<
"): ";
772 for(
int i=0;i<end-start;i++)
773 ss << cLIDs[i] <<
" ";
775 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
777 ss <<
"scatter field = ";
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
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)
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)
PHX::MDField< ScalarT > vector
Pushes residual values into the residual vector for a Newton-based solve.
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)
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.
bool isParameter(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>