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;
120 tangentFields_.resize(numFields);
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;
371 tangentFields_.resize(numFields);
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;
504 string blockId(this->wda(workset).block_id);
505 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
506 int numFields(gatherFields_.size()), numCells(localCellIds.size());
511 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
513 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
514 int indexerId(indexerIds_[fieldInd]),
515 subFieldNum(subFieldIds_[fieldInd]);
518 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
519 auto subRowIndexer = indexers_[indexerId];
520 const vector<int>& elmtOffset =
521 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
522 int numBases(elmtOffset.size());
525 for (
int cell(0); cell < numCells; ++cell)
527 LO cellLocalId = localCellIds[cell];
528 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
531 for (
int basis(0); basis < numBases; ++basis)
533 int offset(elmtOffset[basis]), lid(LIDs[offset]);
534 field(cell, basis) = (*xEvRoGed)[lid];
542 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
544 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
545 int indexerId(indexerIds_[fieldInd]),
546 subFieldNum(subFieldIds_[fieldInd]);
550 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
551 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
552 auto subRowIndexer = indexers_[indexerId];
553 const vector<int>& elmtOffset =
554 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
555 int numBases(elmtOffset.size());
558 for (
int cell(0); cell < numCells; ++cell)
560 LO cellLocalId = localCellIds[cell];
561 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
564 for (
int basis(0); basis < numBases; ++basis)
566 int offset(elmtOffset[basis]), lid(LIDs[offset]);
567 field(cell, basis) = x[lid];
574 if (hasTangentFields_)
577 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
579 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
580 int indexerId(indexerIds_[fieldInd]),
581 subFieldNum(subFieldIds_[fieldInd]);
582 auto subRowIndexer = indexers_[indexerId];
583 const vector<int>& elmtOffset =
584 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
585 int numBases(elmtOffset.size());
588 for (
int cell(0); cell < numCells; ++cell)
591 for (
int basis(0); basis < numBases; ++basis)
593 int numTangentFields(tangentFields_[fieldInd].size());
594 for (
int i(0); i < numTangentFields; ++i)
595 field(cell, basis).fastAccessDx(i) =
596 tangentFields_[fieldInd][i](cell, basis).val();
608 template<
typename TRAITS,
typename LO,
typename GO>
620 using PHX::typeAsString;
641 MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
642 gatherFields_[fd] = f;
643 this->addEvaluatedField(gatherFields_[fd]);
647 string firstName(
"<none>"), n(
"GatherSolution (BlockedEpetra");
649 firstName = names[0];
650 if (disableSensitivities_)
651 n +=
", No Sensitivities";
652 n +=
"): " + firstName +
" (" + typeAsString<EvalT>() +
")";
661 template<
typename TRAITS,
typename LO,
typename GO>
666 typename TRAITS::SetupData ,
678 const string& fieldName(indexerNames_[fd]);
680 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
683 indexerNames_.clear();
691 template<
typename TRAITS,
typename LO,
typename GO>
696 typename TRAITS::PreEvalData d)
699 using std::logic_error;
702 using Teuchos::rcp_dynamic_cast;
710 applySensitivities_ =
false;
711 if ((not disableSensitivities_ ) and
712 (d.first_sensitivities_name == sensitivitiesName_))
713 applySensitivities_ =
true;
717 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
718 if (d.gedc->containsDataObject(globalDataKey_ + post))
720 ged = d.gedc->getDataObject(globalDataKey_ + post);
721 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
727 ged = d.gedc->getDataObject(globalDataKey_);
730 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
731 auto beLoc = rcp_dynamic_cast<
const BELOC>(ged);
732 if (not roGed.is_null())
733 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
734 else if (not beLoc.is_null())
736 if (useTimeDerivativeSolutionVector_)
737 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
739 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
741 "Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ <<
742 "\" (" << post <<
"). A cast failed for " << ged <<
". Type is " <<
753 template<
typename TRAITS,
typename LO,
typename GO>
758 typename TRAITS::EvalData workset)
765 using Teuchos::ptrFromRef;
767 using Teuchos::rcp_dynamic_cast;
769 using Thyra::SpmdVectorBase;
773 string blockId(this->wda(workset).block_id);
774 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
775 int numFields(gatherFields_.size()), numCells(localCellIds.size());
777 if (applySensitivities_)
779 if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
780 seedValue = workset.alpha;
781 else if (gatherSeedIndex_ < 0)
782 seedValue = workset.beta;
783 else if (not useTimeDerivativeSolutionVector_)
784 seedValue = workset.gather_seeds[gatherSeedIndex_];
792 if (not applySensitivities_)
795 vector<int> blockOffsets;
797 int numDerivs(blockOffsets[blockOffsets.size() - 1]);
805 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
807 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
808 int indexerId(indexerIds_[fieldInd]),
809 subFieldNum(subFieldIds_[fieldInd]);
812 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
813 auto subRowIndexer = indexers_[indexerId];
814 const vector<int>& elmtOffset =
815 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
816 int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
819 for (
int cell(0); cell < numCells; ++cell)
821 LO cellLocalId = localCellIds[cell];
822 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
825 for (
int basis(0); basis < numBases; ++basis)
828 int offset(elmtOffset[basis]), lid(LIDs[offset]);
829 field(cell, basis) =
ScalarT(numDerivs, (*xEvRoGed)[lid]);
830 field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
838 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
840 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
841 int indexerId(indexerIds_[fieldInd]),
842 subFieldNum(subFieldIds_[fieldInd]);
846 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
847 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
848 auto subRowIndexer = indexers_[indexerId];
849 const vector<int>& elmtOffset =
850 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
851 int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
854 for (
int cell(0); cell < numCells; ++cell)
856 LO cellLocalId = localCellIds[cell];
857 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
860 for (
int basis(0); basis < numBases; ++basis)
863 int offset(elmtOffset[basis]), lid(LIDs[offset]);
865 field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
872 #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)