43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_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" 61 #include "Phalanx_DataLayout_MDALayout.hpp" 63 #include "Teuchos_FancyOStream.hpp" 70 template<
typename TRAITS,
typename LO,
typename GO>
75 : globalIndexer_(indexer)
76 , globalDataKey_(
"Residual Scatter Container")
78 std::string scatterName = p.
get<std::string>(
"Scatter Name");
83 const std::vector<std::string>& names =
90 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
96 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
97 local_side_id_ = p.
get<
int>(
"Local Side ID");
102 for (std::size_t eq = 0; eq < names.size(); ++eq) {
103 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
109 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
111 applyBC_.resize(names.size());
112 for (std::size_t eq = 0; eq < names.size(); ++eq) {
113 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
114 this->addDependentField(applyBC_[eq]);
119 this->addEvaluatedField(*scatterHolder_);
121 if (p.
isType<std::string>(
"Global Data Key"))
122 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
124 this->setName(scatterName+
" Scatter Residual");
128 template<
typename TRAITS,
typename LO,
typename GO>
138 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
139 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
144 this->utils.setFieldData(applyBC_[fd],fm);
152 template<
typename TRAITS,
typename LO,
typename GO>
159 if(epetraContainer_==Teuchos::null) {
164 dirichletCounter_ = Teuchos::null;
171 dirichletCounter_ = epetraContainer->
get_f();
177 template<
typename TRAITS,
typename LO,
typename GO>
181 std::vector<int> LIDs;
184 std::string blockId = this->wda(workset).block_id;
185 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
189 epetraContainer->
get_f() :
190 epetraContainer->
get_x();
197 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
198 std::size_t cellLocalId = localCellIds[worksetCellIndex];
200 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
203 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
204 int fieldNum = fieldIds_[fieldIndex];
208 const std::pair<std::vector<int>,std::vector<int> > & indicePair
209 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
210 const std::vector<int> & elmtOffset = indicePair.first;
211 const std::vector<int> & basisIdMap = indicePair.second;
220 int basisId = basisIdMap[
basis];
223 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
226 (*r)[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
229 if(dirichletCounter_!=Teuchos::null)
230 (*dirichletCounter_)[lid] = 1.0;
234 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
246 if(dirichletCounter_!=Teuchos::null)
247 (*dirichletCounter_)[lid] = 1.0;
259 template<
typename TRAITS,
typename LO,
typename GO>
264 : globalIndexer_(indexer)
265 , globalDataKey_(
"Residual Scatter Container")
267 std::string scatterName = p.
get<std::string>(
"Scatter Name");
272 const std::vector<std::string>& names =
279 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
285 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
286 local_side_id_ = p.
get<
int>(
"Local Side ID");
291 for (std::size_t eq = 0; eq < names.size(); ++eq) {
292 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
298 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
300 applyBC_.resize(names.size());
301 for (std::size_t eq = 0; eq < names.size(); ++eq) {
302 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
303 this->addDependentField(applyBC_[eq]);
308 this->addEvaluatedField(*scatterHolder_);
310 if (p.
isType<std::string>(
"Global Data Key"))
311 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
313 this->setName(scatterName+
" Scatter Tangent");
317 template<
typename TRAITS,
typename LO,
typename GO>
327 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
328 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
333 this->utils.setFieldData(applyBC_[fd],fm);
341 template<
typename TRAITS,
typename LO,
typename GO>
348 if(epetraContainer_==Teuchos::null) {
353 dirichletCounter_ = Teuchos::null;
360 dirichletCounter_ = epetraContainer->
get_f();
365 using Teuchos::rcp_dynamic_cast;
368 std::vector<std::string> activeParameters =
373 dfdp_vectors_.clear();
374 for(std::size_t i=0;i<activeParameters.size();i++) {
376 dfdp_vectors_.push_back(vec);
381 template<
typename TRAITS,
typename LO,
typename GO>
385 std::vector<int> LIDs;
388 std::string blockId = this->wda(workset).block_id;
389 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
393 epetraContainer->
get_f() :
394 epetraContainer->
get_x();
401 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
402 std::size_t cellLocalId = localCellIds[worksetCellIndex];
405 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
408 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
409 int fieldNum = fieldIds_[fieldIndex];
413 const std::pair<std::vector<int>,std::vector<int> > & indicePair
414 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
415 const std::vector<int> & elmtOffset = indicePair.first;
416 const std::vector<int> & basisIdMap = indicePair.second;
425 int basisId = basisIdMap[
basis];
428 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
436 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
437 (*dfdp_vectors_[d])[lid] = 0.0;
439 for(
int d=0;d<
value.size();d++) {
440 (*dfdp_vectors_[d])[lid] =
value.fastAccessDx(d);
444 if(dirichletCounter_!=Teuchos::null)
445 (*dirichletCounter_)[lid] = 1.0;
449 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
463 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
464 (*dfdp_vectors_[d])[lid] = 0.0;
466 for(
int d=0;d<
value.size();d++) {
467 (*dfdp_vectors_[d])[lid] =
value.fastAccessDx(d);
471 if(dirichletCounter_!=Teuchos::null)
472 (*dirichletCounter_)[lid] = 1.0;
483 template<
typename TRAITS,
typename LO,
typename GO>
488 : globalIndexer_(indexer)
489 , colGlobalIndexer_(cIndexer)
490 , globalDataKey_(
"Residual Scatter Container")
491 , preserveDiagonal_(false)
493 std::string scatterName = p.
get<std::string>(
"Scatter Name");
498 const std::vector<std::string>& names =
507 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
508 local_side_id_ = p.
get<
int>(
"Local Side ID");
512 for (std::size_t eq = 0; eq < names.size(); ++eq) {
513 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
519 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
521 applyBC_.resize(names.size());
522 for (std::size_t eq = 0; eq < names.size(); ++eq) {
523 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
524 this->addDependentField(applyBC_[eq]);
529 this->addEvaluatedField(*scatterHolder_);
531 if (p.
isType<std::string>(
"Global Data Key"))
532 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
534 if (p.
isType<
bool>(
"Preserve Diagonal"))
535 preserveDiagonal_ = p.
get<
bool>(
"Preserve Diagonal");
537 this->setName(scatterName+
" Scatter Residual (Jacobian)");
541 template<
typename TRAITS,
typename LO,
typename GO>
551 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
552 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
557 this->utils.setFieldData(applyBC_[fd],fm);
566 template<
typename TRAITS,
typename LO,
typename GO>
573 if(epetraContainer_==Teuchos::null) {
578 dirichletCounter_ = Teuchos::null;
585 dirichletCounter_ = epetraContainer->
get_f();
591 template<
typename TRAITS,
typename LO,
typename GO>
595 std::vector<int> cLIDs, rLIDs;
596 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
599 std::string blockId = this->wda(workset).block_id;
600 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
612 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
613 std::size_t cellLocalId = localCellIds[worksetCellIndex];
615 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
617 cLIDs = colGlobalIndexer_->getElementLIDs(cellLocalId);
622 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
623 int fieldNum = fieldIds_[fieldIndex];
626 const std::pair<std::vector<int>,std::vector<int> > & indicePair
627 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
628 const std::vector<int> & elmtOffset = indicePair.first;
629 const std::vector<int> & basisIdMap = indicePair.second;
638 int basisId = basisIdMap[
basis];
641 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
647 int * rowIndices = 0;
648 double * rowValues = 0;
652 for(
int i=0;i<numEntries;i++) {
653 if(preserveDiagonal_) {
654 if(row!=rowIndices[i])
666 (*r)[row] = scatterField.val();
667 if(dirichletCounter_!=Teuchos::null) {
669 (*dirichletCounter_)[row] = 1.0;
673 std::vector<double> jacRow(scatterField.size(),0.0);
675 if(!preserveDiagonal_) {
676 int err = Jac->
ReplaceMyValues(row, cLIDs.size(), scatterField.dx(), &cLIDs[0]);
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
panzer::Traits::Tangent::ScalarT ScalarT
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const Teuchos::RCP< Epetra_Vector > get_x() const
panzer::Traits::Jacobian::ScalarT ScalarT
bool isType(const std::string &name) const
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>
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
const Teuchos::RCP< Epetra_Vector > get_f() const
Pushes residual values into the residual vector for a Newton-based solve.