43 #ifndef __Panzer_GatherSolution_BlockedEpetra_Hessian_impl_hpp__
44 #define __Panzer_GatherSolution_BlockedEpetra_Hessian_impl_hpp__
47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
65 #include "Phalanx_DataLayout.hpp"
68 #include "Teuchos_Assert.hpp"
69 #include "Teuchos_FancyOStream.hpp"
72 #include "Thyra_ProductVectorBase.hpp"
73 #include "Thyra_SpmdVectorBase.hpp"
80 template<
typename TRAITS,
typename LO,
typename GO>
91 using PHX::typeAsString;
115 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
116 this->addEvaluatedField(gatherFields_[fd]);
120 string firstName(
"<none>"), n(
"GatherSolution (BlockedEpetra");
122 firstName = names[0];
123 if (not firstSensitivitiesAvailable_)
124 n +=
", No First Sensitivities";
125 n +=
"): " + firstName +
" (" + typeAsString<EvalT>() +
")";
134 template<
typename TRAITS,
typename LO,
typename GO>
138 typename TRAITS::SetupData ,
150 const string& fieldName(indexerNames_[fd]);
152 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
155 indexerNames_.clear();
163 template<
typename TRAITS,
typename LO,
typename GO>
167 typename TRAITS::PreEvalData d)
170 using std::logic_error;
173 using Teuchos::rcp_dynamic_cast;
180 firstApplySensitivities_ =
false;
181 if ((firstSensitivitiesAvailable_ ) and
182 (d.first_sensitivities_name == sensitivitiesName_))
183 firstApplySensitivities_ =
true;
184 secondApplySensitivities_ =
false;
185 if ((secondSensitivitiesAvailable_ ) and
186 (d.second_sensitivities_name == sensitivitiesName_))
187 secondApplySensitivities_ =
true;
191 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
192 if (d.gedc->containsDataObject(globalDataKey_ + post))
194 ged = d.gedc->getDataObject(globalDataKey_ + post);
195 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
199 else if (d.gedc->containsDataObject(globalDataKey_))
201 ged = d.gedc->getDataObject(globalDataKey_);
204 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
205 auto blockedContainer = rcp_dynamic_cast<
const BELOC>(ged);
206 if (not roGed.is_null())
207 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
208 else if (not blockedContainer.is_null())
210 if (useTimeDerivativeSolutionVector_)
211 x_ = rcp_dynamic_cast<ProductVectorBase<double>>
212 (blockedContainer->get_dxdt());
214 x_ = rcp_dynamic_cast<ProductVectorBase<double>>
215 (blockedContainer->get_x());
223 if (not secondApplySensitivities_)
227 if (d.gedc->containsDataObject(sensitivities2ndPrefix_ + globalDataKey_))
229 ged = d.gedc->getDataObject(sensitivities2ndPrefix_ + globalDataKey_);
230 dxBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
235 "Cannot find sensitivity vector associated with \"" +
236 sensitivities2ndPrefix_ + globalDataKey_ +
"\" and \"" + post +
"\".");
244 template<
typename TRAITS,
typename LO,
typename GO>
248 typename TRAITS::EvalData workset)
255 using Teuchos::ptrFromRef;
257 using Teuchos::rcp_dynamic_cast;
259 using Thyra::SpmdVectorBase;
263 string blockId(this->wda(workset).block_id);
264 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
265 int numFields(gatherFields_.size()), numCells(localCellIds.size());
267 if (firstApplySensitivities_)
269 if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
270 seedValue = workset.alpha;
271 else if (gatherSeedIndex_ < 0)
272 seedValue = workset.beta;
273 else if (not useTimeDerivativeSolutionVector_)
274 seedValue = workset.gather_seeds[gatherSeedIndex_];
282 if (not firstApplySensitivities_)
285 vector<int> blockOffsets;
290 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
292 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
293 int indexerId(indexerIds_[fieldInd]),
294 subFieldNum(subFieldIds_[fieldInd]);
297 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
298 auto subRowIndexer = indexers_[indexerId];
299 const vector<int>& elmtOffset =
300 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
301 int numBases(elmtOffset.size());
304 for (
int cell(0); cell < numCells; ++cell)
306 size_t cellLocalId(localCellIds[cell]);
307 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
310 for (
int basis(0); basis < numBases; ++basis)
312 int offset(elmtOffset[basis]), lid(LIDs[offset]);
313 field(cell, basis) = (*xEvRoGed)[lid];
321 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
323 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
324 int indexerId(indexerIds_[fieldInd]),
325 subFieldNum(subFieldIds_[fieldInd]);
329 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
330 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
331 auto subRowIndexer = indexers_[indexerId];
332 const vector<int>& elmtOffset =
333 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
334 int numBases(elmtOffset.size());
337 for (
int cell(0); cell < numCells; ++cell)
339 size_t cellLocalId(localCellIds[cell]);
340 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
343 for (
int basis(0); basis < numBases; ++basis)
345 int offset(elmtOffset[basis]), lid(LIDs[offset]);
346 field(cell, basis) = x[lid];
353 if (firstApplySensitivities_)
356 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
358 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
359 int indexerId(indexerIds_[fieldInd]),
360 subFieldNum(subFieldIds_[fieldInd]);
361 auto subRowIndexer = indexers_[indexerId];
362 const vector<int>& elmtOffset =
363 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
364 int startBlkOffset(blockOffsets[indexerId]),
365 numBases(elmtOffset.size());
368 for (
int cell(0); cell < numCells; ++cell)
371 for (
int basis(0); basis < numBases; ++basis)
373 int offset(elmtOffset[basis]);
374 field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
381 if (secondApplySensitivities_)
384 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
386 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
387 int indexerId(indexerIds_[fieldInd]),
388 subFieldNum(subFieldIds_[fieldInd]);
391 auto dxEvRoGed = dxBvRoGed_->getGEDBlock(indexerId);
392 auto subRowIndexer = indexers_[indexerId];
393 const vector<int>& elmtOffset =
394 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
395 int numBases(elmtOffset.size());
398 for (
int cell(0); cell < numCells; ++cell)
400 size_t cellLocalId(localCellIds[cell]);
401 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
404 for (
int basis(0); basis < numBases; ++basis)
406 int offset(elmtOffset[basis]), lid(LIDs[offset]);
407 field(cell, basis).val().fastAccessDx(0) = (*dxEvRoGed)[lid];
414 #endif // Panzer_BUILD_HESSIAN_SUPPORT
416 #endif // __Panzer_GatherSolution_BlockedEpetra_Hessian_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)
std::string typeName(const T &t)