43 #ifndef __Panzer_GatherSolution_BlockedEpetra_impl_hpp__ 44 #define __Panzer_GatherSolution_BlockedEpetra_impl_hpp__ 62 #include "Phalanx_DataLayout.hpp" 65 #include "Teuchos_Assert.hpp" 66 #include "Teuchos_FancyOStream.hpp" 69 #include "Thyra_ProductVectorBase.hpp" 70 #include "Thyra_SpmdVectorBase.hpp" 77 template<
typename TRAITS,
typename LO,
typename GO>
86 hasTangentFields_(false)
90 using PHX::typeAsString;
95 using vvstring = std::vector<std::vector<std::string>>;
111 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
112 this->addEvaluatedField(gatherFields_[fd]);
116 if (tangentFieldNames.size() > 0)
119 hasTangentFields_ =
true;
123 int numTangentFields(tangentFieldNames[fd].size());
124 tangentFields_[fd].resize(numTangentFields);
125 for (
int i(0); i < numTangentFields; ++i)
127 tangentFields_[fd][i] =
128 MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
130 this->addDependentField(tangentFields_[fd][i]);
136 string firstName(
"<none>");
138 firstName = names[0];
139 string n(
"GatherSolution (BlockedEpetra): " + firstName +
" (" +
140 typeAsString<EvalT>() +
")");
149 template<
typename TRAITS,
typename LO,
typename GO>
154 typename TRAITS::SetupData ,
166 const string& fieldName(indexerNames_[fd]);
168 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
171 indexerNames_.clear();
179 template<
typename TRAITS,
typename LO,
typename GO>
184 typename TRAITS::PreEvalData d)
186 using std::logic_error;
189 using Teuchos::rcp_dynamic_cast;
198 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
199 if (d.gedc->containsDataObject(globalDataKey_ + post))
201 ged = d.gedc->getDataObject(globalDataKey_ + post);
202 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
208 ged = d.gedc->getDataObject(globalDataKey_);
211 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
212 auto beLoc = rcp_dynamic_cast<
const BELOC>(ged);
213 if (not roGed.is_null())
214 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
215 else if (not beLoc.is_null())
217 if (useTimeDerivativeSolutionVector_)
218 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
220 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
222 "Residual: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
223 "\" (" << post <<
"). A cast failed for " << ged <<
". Type is " <<
234 template<
typename TRAITS,
typename LO,
typename GO>
239 typename TRAITS::EvalData workset)
246 using Teuchos::ptrFromRef;
248 using Teuchos::rcp_dynamic_cast;
250 using Thyra::SpmdVectorBase;
254 string blockId(this->wda(workset).block_id);
255 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
256 int numFields(gatherFields_.size()), numCells(localCellIds.size());
261 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
263 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
264 int indexerId(indexerIds_[fieldInd]),
265 subFieldNum(subFieldIds_[fieldInd]);
268 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
269 auto subRowIndexer = indexers_[indexerId];
270 const vector<int>& elmtOffset =
271 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
272 int numBases(elmtOffset.size());
275 for (
int cell(0); cell < numCells; ++cell)
277 LO cellLocalId = localCellIds[cell];
278 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
281 for (
int basis(0); basis < numBases; ++basis)
283 int offset(elmtOffset[basis]), lid(LIDs[offset]);
284 field(cell, basis) = (*xEvRoGed)[lid];
292 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
294 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
295 int indexerId(indexerIds_[fieldInd]),
296 subFieldNum(subFieldIds_[fieldInd]);
300 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
301 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
302 auto subRowIndexer = indexers_[indexerId];
303 const vector<int>& elmtOffset =
304 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
305 int numBases(elmtOffset.size());
308 for (
int cell(0); cell < numCells; ++cell)
310 LO cellLocalId = localCellIds[cell];
311 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
314 for (
int basis(0); basis < numBases; ++basis)
316 int offset(elmtOffset[basis]), lid(LIDs[offset]);
317 field(cell, basis) = x[lid];
329 template<
typename TRAITS,
typename LO,
typename GO>
337 hasTangentFields_(false)
341 using PHX::typeAsString;
346 using vvstring = std::vector<std::vector<std::string>>;
362 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
363 this->addEvaluatedField(gatherFields_[fd]);
367 if (tangentFieldNames.size() > 0)
370 hasTangentFields_ =
true;
374 int numTangentFields(tangentFieldNames[fd].size());
375 tangentFields_[fd].resize(numTangentFields);
376 for (
int i(0); i < numTangentFields; ++i)
378 tangentFields_[fd][i] =
379 MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
381 this->addDependentField(tangentFields_[fd][i]);
387 string firstName(
"<none>");
389 firstName = names[0];
390 string n(
"GatherSolution Tangent (BlockedEpetra): " + firstName +
" (" +
391 typeAsString<EvalT>() +
")");
400 template<
typename TRAITS,
typename LO,
typename GO>
404 typename TRAITS::SetupData ,
416 const string& fieldName(indexerNames_[fd]);
418 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
421 indexerNames_.clear();
429 template<
typename TRAITS,
typename LO,
typename GO>
433 typename TRAITS::PreEvalData d)
435 using std::logic_error;
438 using Teuchos::rcp_dynamic_cast;
447 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
448 if (d.gedc->containsDataObject(globalDataKey_ + post))
450 ged = d.gedc->getDataObject(globalDataKey_ + post);
451 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
457 ged = d.gedc->getDataObject(globalDataKey_);
460 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
461 auto beLoc = rcp_dynamic_cast<
const BELOC>(ged);
462 if (not roGed.is_null())
463 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
464 else if (not beLoc.is_null())
466 if (useTimeDerivativeSolutionVector_)
467 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
469 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
471 "Tangent: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
472 "\" (" << post <<
"). A cast failed for " << ged <<
". Type is " <<
483 template<
typename TRAITS,
typename LO,
typename GO>
487 typename TRAITS::EvalData workset)
495 using Teuchos::ptrFromRef;
497 using Teuchos::rcp_dynamic_cast;
499 using Thyra::SpmdVectorBase;
503 vector<pair<int, GO>> GIDs;
505 string blockId(this->wda(workset).block_id);
506 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
507 int numFields(gatherFields_.size()), numCells(localCellIds.size());
512 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
514 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
515 int indexerId(indexerIds_[fieldInd]),
516 subFieldNum(subFieldIds_[fieldInd]);
519 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
520 auto subRowIndexer = indexers_[indexerId];
521 const vector<int>& elmtOffset =
522 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
523 int numBases(elmtOffset.size());
526 for (
int cell(0); cell < numCells; ++cell)
528 LO cellLocalId = localCellIds[cell];
529 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
532 for (
int basis(0); basis < numBases; ++basis)
534 int offset(elmtOffset[basis]), lid(LIDs[offset]);
535 field(cell, basis) = (*xEvRoGed)[lid];
543 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
545 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
546 int indexerId(indexerIds_[fieldInd]),
547 subFieldNum(subFieldIds_[fieldInd]);
551 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
552 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
553 auto subRowIndexer = indexers_[indexerId];
554 const vector<int>& elmtOffset =
555 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
556 int numBases(elmtOffset.size());
559 for (
int cell(0); cell < numCells; ++cell)
561 LO cellLocalId = localCellIds[cell];
562 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
565 for (
int basis(0); basis < numBases; ++basis)
567 int offset(elmtOffset[basis]), lid(LIDs[offset]);
568 field(cell, basis) = x[lid];
575 if (hasTangentFields_)
578 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
580 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
581 int indexerId(indexerIds_[fieldInd]),
582 subFieldNum(subFieldIds_[fieldInd]);
583 auto subRowIndexer = indexers_[indexerId];
584 const vector<int>& elmtOffset =
585 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
586 int numBases(elmtOffset.size());
589 for (
int cell(0); cell < numCells; ++cell)
592 for (
int basis(0); basis < numBases; ++basis)
594 int numTangentFields(tangentFields_[fieldInd].size());
595 for (
int i(0); i < numTangentFields; ++i)
596 field(cell, basis).fastAccessDx(i) =
597 tangentFields_[fieldInd][i](cell, basis).val();
609 template<
typename TRAITS,
typename LO,
typename GO>
621 using PHX::typeAsString;
642 MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
643 gatherFields_[fd] = f;
644 this->addEvaluatedField(gatherFields_[fd]);
648 string firstName(
"<none>"), n(
"GatherSolution (BlockedEpetra");
650 firstName = names[0];
651 if (disableSensitivities_)
652 n +=
", No Sensitivities";
653 n +=
"): " + firstName +
" (" + typeAsString<EvalT>() +
")";
662 template<
typename TRAITS,
typename LO,
typename GO>
667 typename TRAITS::SetupData ,
679 const string& fieldName(indexerNames_[fd]);
681 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
684 indexerNames_.clear();
692 template<
typename TRAITS,
typename LO,
typename GO>
697 typename TRAITS::PreEvalData d)
700 using std::logic_error;
703 using Teuchos::rcp_dynamic_cast;
711 applySensitivities_ =
false;
712 if ((not disableSensitivities_ ) and
713 (d.first_sensitivities_name == sensitivitiesName_))
714 applySensitivities_ =
true;
718 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
719 if (d.gedc->containsDataObject(globalDataKey_ + post))
721 ged = d.gedc->getDataObject(globalDataKey_ + post);
722 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
728 ged = d.gedc->getDataObject(globalDataKey_);
731 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
732 auto beLoc = rcp_dynamic_cast<
const BELOC>(ged);
733 if (not roGed.is_null())
734 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
735 else if (not beLoc.is_null())
737 if (useTimeDerivativeSolutionVector_)
738 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
740 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
742 "Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ <<
743 "\" (" << post <<
"). A cast failed for " << ged <<
". Type is " <<
754 template<
typename TRAITS,
typename LO,
typename GO>
759 typename TRAITS::EvalData workset)
766 using Teuchos::ptrFromRef;
768 using Teuchos::rcp_dynamic_cast;
770 using Thyra::SpmdVectorBase;
774 string blockId(this->wda(workset).block_id);
775 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
776 int numFields(gatherFields_.size()), numCells(localCellIds.size());
778 if (applySensitivities_)
780 if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
781 seedValue = workset.alpha;
782 else if (gatherSeedIndex_ < 0)
783 seedValue = workset.beta;
784 else if (not useTimeDerivativeSolutionVector_)
785 seedValue = workset.gather_seeds[gatherSeedIndex_];
793 if (not applySensitivities_)
796 vector<int> blockOffsets;
798 int numDerivs(blockOffsets[blockOffsets.size() - 1]);
806 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
808 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
809 int indexerId(indexerIds_[fieldInd]),
810 subFieldNum(subFieldIds_[fieldInd]);
813 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
814 auto subRowIndexer = indexers_[indexerId];
815 const vector<int>& elmtOffset =
816 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
817 int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
820 for (
int cell(0); cell < numCells; ++cell)
822 LO cellLocalId = localCellIds[cell];
823 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
826 for (
int basis(0); basis < numBases; ++basis)
829 int offset(elmtOffset[basis]), lid(LIDs[offset]);
830 field(cell, basis) =
ScalarT(numDerivs, (*xEvRoGed)[lid]);
831 field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
839 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
841 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
842 int indexerId(indexerIds_[fieldInd]),
843 subFieldNum(subFieldIds_[fieldInd]);
847 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
848 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
849 auto subRowIndexer = indexers_[indexerId];
850 const vector<int>& elmtOffset =
851 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
852 int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
855 for (
int cell(0); cell < numCells; ++cell)
857 LO cellLocalId = localCellIds[cell];
858 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
861 for (
int basis(0); basis < numBases; ++basis)
864 int offset(elmtOffset[basis]), lid(LIDs[offset]);
866 field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
873 #endif // __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
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)
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Post-Registration Setup.
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
Constructor.
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral...
void evaluateFields(typename TRAITS::EvalData d)
Evaluate Fields.
Description and data layouts associated with a particular basis.
#define TEUCHOS_ASSERT(assertion_test)
panzer::Traits::Jacobian::ScalarT ScalarT
The scalar type.
std::string typeName(const T &t)