43 #ifndef __Panzer_GatherSolution_Epetra_Hessian_impl_hpp__
44 #define __Panzer_GatherSolution_Epetra_Hessian_impl_hpp__
47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
56 #include "Epetra_Map.h"
57 #include "Epetra_Vector.h"
68 #include "Phalanx_DataLayout.hpp"
71 #include "Teuchos_Assert.hpp"
72 #include "Teuchos_FancyOStream.hpp"
75 #include "Thyra_SpmdVectorBase.hpp"
82 template<
typename TRAITS,
typename LO,
typename GO>
88 globalIndexer_(indexer)
92 using PHX::typeAsString;
116 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
117 this->addEvaluatedField(gatherFields_[fd]);
121 string firstName(
"<none>"), n(
"Gather Solution (Epetra");
123 firstName = names[0];
124 if (not firstSensitivitiesAvailable_)
125 n +=
", No First Sensitivities";
126 n +=
"): " + firstName +
" (" + typeAsString<EvalT>() +
") ";
135 template<
typename TRAITS,
typename LO,
typename GO>
139 typename TRAITS::SetupData ,
142 using std::logic_error;
151 const string& fieldName(indexerNames_[fd]);
152 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
156 "GatherSolution_Epetra<Hessian>: Could not find field \"" + fieldName
157 +
"\" in the global indexer. ")
159 indexerNames_.clear();
167 template<
typename TRAITS,
typename LO,
typename GO>
171 typename TRAITS::PreEvalData d)
173 using std::logic_error;
176 using Teuchos::rcp_dynamic_cast;
184 firstApplySensitivities_ =
false;
185 if ((firstSensitivitiesAvailable_ ) and
186 (d.first_sensitivities_name == sensitivitiesName_))
187 firstApplySensitivities_ =
true;
188 secondApplySensitivities_ =
false;
189 if ((secondSensitivitiesAvailable_ ) and
190 (d.second_sensitivities_name == sensitivitiesName_))
191 secondApplySensitivities_ =
true;
195 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
196 if (d.gedc->containsDataObject(globalDataKey_ + post))
198 ged = d.gedc->getDataObject(globalDataKey_ + post);
199 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged,
true);
203 else if (d.gedc->containsDataObject(globalDataKey_))
205 ged = d.gedc->getDataObject(globalDataKey_);
208 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged);
209 auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
213 auto locPair = rcp_dynamic_cast<LPGED>(ged);
214 if (not locPair.is_null())
216 RCP<LOC> loc = locPair->getGhostedLOC();
217 epetraContainer = rcp_dynamic_cast<ELOC>(loc);
220 if ((xEvRoGed_.is_null()) and (not epetraContainer.is_null()))
222 if (useTimeDerivativeSolutionVector_)
223 x_ = epetraContainer->get_dxdt();
225 x_ = epetraContainer->get_x();
231 logic_error,
"GatherSolution_Epetra_Hessian::preEvaluate(): Unable to " \
232 "find solution vector.")
235 if (not secondApplySensitivities_)
239 if (d.gedc->containsDataObject(sensitivities2ndPrefix_ + globalDataKey_))
241 ged = d.gedc->getDataObject(sensitivities2ndPrefix_ + globalDataKey_);
242 dxEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged,
true);
246 "find sensitivity vector associated with \"" + sensitivities2ndPrefix_ +
247 globalDataKey_ +
"\" and \"" + post +
"\".")
256 template<typename TRAITS, typename LO, typename GO>
260 typename TRAITS::EvalData workset)
267 using Teuchos::ptrFromRef;
269 using Teuchos::rcp_dynamic_cast;
270 using Thyra::SpmdVectorBase;
273 string blockId(this->wda(workset).block_id);
274 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
275 int numFields(gatherFields_.size()), numCells(localCellIds.size());
279 if (firstApplySensitivities_)
281 if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
282 seedValue = workset.alpha;
283 else if (gatherSeedIndex_ < 0)
284 seedValue = workset.beta;
285 else if (not useTimeDerivativeSolutionVector_)
286 seedValue = workset.gather_seeds[gatherSeedIndex_];
294 if (not firstApplySensitivities_)
301 if (this->wda.getDetailsIndex() == 1)
304 dos = globalIndexer_->getElementBlockGIDCount(workset.details(0).block_id);
310 for (
int fieldIndex(0); fieldIndex <
numFields; ++fieldIndex)
312 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldIndex];
313 int fieldNum(fieldIds_[fieldIndex]);
314 const vector<int>& elmtOffset =
315 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
316 int numBases(elmtOffset.size());
319 for (
int cell(0); cell < numCells; ++cell)
321 size_t cellLocalId(localCellIds[cell]);
322 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
325 for (
int basis(0); basis < numBases; ++basis)
327 int offset(elmtOffset[basis]), lid(LIDs[offset]);
328 field(cell, basis) = (*xEvRoGed_)[lid];
336 for (
int fieldIndex(0); fieldIndex <
numFields; ++fieldIndex)
338 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldIndex];
339 int fieldNum(fieldIds_[fieldIndex]);
340 const vector<int>& elmtOffset =
341 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
342 int numBases(elmtOffset.size());
345 for (
int cell(0); cell < numCells; ++cell)
347 size_t cellLocalId(localCellIds[cell]);
348 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
351 for (
int basis(0); basis < numBases; ++basis)
353 int offset(elmtOffset[basis]), lid(LIDs[offset]);
354 field(cell, basis) = (*x_)[lid];
361 if (firstApplySensitivities_)
364 for (
int fieldIndex(0); fieldIndex <
numFields; ++fieldIndex)
366 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldIndex];
367 int fieldNum(fieldIds_[fieldIndex]);
368 const vector<int>& elmtOffset =
369 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
370 int numBases(elmtOffset.size());
373 for (
int cell(0); cell < numCells; ++cell)
376 for (
int basis(0); basis < numBases; ++basis)
378 int offset(elmtOffset[basis]);
379 field(cell, basis).fastAccessDx(dos + offset) = seedValue;
386 if (secondApplySensitivities_)
389 for (
int fieldIndex(0); fieldIndex <
numFields; ++fieldIndex)
391 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldIndex];
392 int fieldNum(fieldIds_[fieldIndex]);
393 const vector<int>& elmtOffset =
394 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
395 int numBases(elmtOffset.size());
398 for (
int cell(0); cell < numCells; ++cell)
400 size_t cellLocalId(localCellIds[cell]);
401 auto LIDs = globalIndexer_->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_Epetra_Hessian_impl_hpp__
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
This class provides a boundary exchange communication mechanism for vectors.
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...
GatherSolution_Epetra (Hessian Specialization)
Description and data layouts associated with a particular basis.
#define TEUCHOS_ASSERT(assertion_test)