43 #ifndef __Panzer_GatherSolution_BlockedEpetra_Hessian_impl_hpp__ 44 #define __Panzer_GatherSolution_BlockedEpetra_Hessian_impl_hpp__ 47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT 49 #include "Teuchos_Assert.hpp" 50 #include "Phalanx_DataLayout.hpp" 60 #include "Teuchos_FancyOStream.hpp" 62 #include "Thyra_SpmdVectorBase.hpp" 63 #include "Thyra_ProductVectorBase.hpp" 71 template<
typename TRAITS,
typename LO,
typename GO>
77 typedef std::vector< std::vector<std::string> > vvstring;
79 GatherSolution_Input
input;
80 input.setParameterList(p);
82 const std::vector<std::string> & names =
input.getDofNames();
85 indexerNames_ =
input.getIndexerNames();
86 useTimeDerivativeSolutionVector_ =
input.useTimeDerivativeSolutionVector();
87 globalDataKey_ =
input.getGlobalDataKey();
89 gatherSeedIndex_ =
input.getGatherSeedIndex();
90 sensitivitiesName_ =
input.getSensitivitiesName();
91 firstSensitivitiesAvailable_ =
input.firstSensitivitiesAvailable();
93 secondSensitivitiesAvailable_ =
input.secondSensitivitiesAvailable();
94 sensitivities2ndPrefix_ =
input.getSecondSensitivityDataKeyPrefix();
97 gatherFields_.resize(names.size());
98 for (std::size_t fd = 0; fd < names.size(); ++fd) {
100 gatherFields_[fd] = f;
101 this->addEvaluatedField(gatherFields_[fd]);
105 std::string firstName =
"<none>";
107 firstName = names[0];
110 if(!firstSensitivitiesAvailable_) {
111 std::string n =
"GatherSolution (BlockedEpetra, No First Sensitivities): "+firstName+
" (Hessian)";
115 std::string n =
"GatherSolution (BlockedEpetra): "+firstName+
" ("+PHX::typeAsString<EvalT>()+
") ";
121 template<
typename TRAITS,
typename LO,
typename GO>
128 indexerIds_.resize(gatherFields_.size());
129 subFieldIds_.resize(gatherFields_.size());
131 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
133 const std::string& fieldName = indexerNames_[fd];
136 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
141 this->utils.setFieldData(gatherFields_[fd],fm);
144 indexerNames_.clear();
148 template<
typename TRAITS,
typename LO,
typename GO>
152 typedef BlockedEpetraLinearObjContainer BLOC;
155 using Teuchos::rcp_dynamic_cast;
159 if(firstSensitivitiesAvailable_) {
160 if(d.first_sensitivities_name==sensitivitiesName_)
161 firstApplySensitivities_ =
true;
163 firstApplySensitivities_ =
false;
166 firstApplySensitivities_ =
false;
168 if(secondSensitivitiesAvailable_) {
169 if(d.second_sensitivities_name==sensitivitiesName_)
170 secondApplySensitivities_ =
true;
172 secondApplySensitivities_ =
false;
175 secondApplySensitivities_ =
false;
179 RCP<GlobalEvaluationData> ged;
182 std::string post = useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X";
183 if(d.gedc.containsDataObject(globalDataKey_+post)) {
184 ged = d.gedc.getDataObject(globalDataKey_+post);
186 RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,
true);
192 "Gather Hessian: Can't find x vector in GEDC \"" << globalDataKey_ << post <<
"\"" 193 "A cast failed for " << ged <<
". Type is " <<
Teuchos::typeName(*ged) << std::endl);
195 else if(d.gedc.containsDataObject(globalDataKey_)) {
196 ged = d.gedc.getDataObject(globalDataKey_);
199 RCP<const BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<
const BlockedVector_ReadOnly_GlobalEvaluationData>(ged);
200 RCP<const BlockedEpetraLinearObjContainer> blockedContainer = rcp_dynamic_cast<
const BLOC>(ged);
202 if(ro_ged!=Teuchos::null) {
203 RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,
true);
207 else if(blockedContainer!=Teuchos::null) {
208 if (useTimeDerivativeSolutionVector_)
216 "Gather Hessian: Can't find x vector in GEDC \"" << globalDataKey_ <<
"\" (" << post <<
"). " 224 if(!secondApplySensitivities_) {
233 if(d.gedc.containsDataObject(sensitivities2ndPrefix_+globalDataKey_)) {
234 ged = d.gedc.getDataObject(sensitivities2ndPrefix_+globalDataKey_);
236 RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,
true);
243 "Cannot find sensitivity vector associated with \""+sensitivities2ndPrefix_+globalDataKey_+
"\" and \""+post+
"\"");
247 template<
typename TRAITS,
typename LO,
typename GO>
252 using Teuchos::rcp_dynamic_cast;
253 using Teuchos::ptrFromRef;
257 using Thyra::SpmdVectorBase;
261 std::string blockId = this->wda(workset).block_id;
262 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
264 double seed_value = 0.0;
265 if(firstApplySensitivities_) {
266 if (useTimeDerivativeSolutionVector_ && gatherSeedIndex_<0) {
267 seed_value = workset.alpha;
269 else if (gatherSeedIndex_<0) {
270 seed_value = workset.beta;
272 else if(!useTimeDerivativeSolutionVector_) {
273 seed_value = workset.gather_seeds[gatherSeedIndex_];
283 if(!firstApplySensitivities_)
286 std::vector<int> blockOffsets;
290 for(std::size_t fieldIndex=0;
291 fieldIndex<gatherFields_.size();fieldIndex++) {
293 PHX::MDField<ScalarT,Cell,NODE> &
field = gatherFields_[fieldIndex];
295 int indexerId = indexerIds_[fieldIndex];
296 int subFieldNum = subFieldIds_[fieldIndex];
300 rcp_dynamic_cast<SpmdVectorBase<double> >(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(local_x));
301 if(secondApplySensitivities_)
302 rcp_dynamic_cast<SpmdVectorBase<double> >(dx_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(local_dx));
304 auto subRowIndexer = indexers_[indexerId];
305 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
307 int startBlkOffset = blockOffsets[indexerId];
310 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
311 std::size_t cellLocalId = localCellIds[worksetCellIndex];
313 const std::vector<int> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
315 if(!firstApplySensitivities_) {
332 field(worksetCellIndex,
basis).val() = local_x[lid];
333 field(worksetCellIndex,
basis).fastAccessDx(startBlkOffset+
offset) = seed_value;
338 if(secondApplySensitivities_) {
342 field(worksetCellIndex,
basis).val().fastAccessDx(0) = local_dx[lid];
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)
PHX::MDField< const ScalarT > input
PHX::MDField< ScalarT > vector
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
PHX::MDField< const ScalarT, Cell, IP > field
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
void evaluateFields(typename TRAITS::EvalData d)
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>
std::string typeName(const T &t)