Panzer  Version of the Day
Panzer_GatherSolution_BlockedEpetra_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
44 #define PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
45 
46 #include "Teuchos_Assert.hpp"
47 #include "Phalanx_DataLayout.hpp"
48 
51 #include "Panzer_PureBasis.hpp"
52 #include "Panzer_EpetraLinearObjFactory.hpp"
56 
57 #include "Teuchos_FancyOStream.hpp"
58 
59 #include "Thyra_SpmdVectorBase.hpp"
60 #include "Thyra_ProductVectorBase.hpp"
61 
62 // **********************************************************************
63 // Specialization: Residual
64 // **********************************************************************
65 
66 template<typename TRAITS,typename LO,typename GO>
69  const Teuchos::ParameterList& p)
70  : indexers_(indexers)
71  , has_tangent_fields_(false)
72 {
73  typedef std::vector< std::vector<std::string> > vvstring;
74 
76  input.setParameterList(p);
77 
78  const std::vector<std::string> & names = input.getDofNames();
79  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
80  const vvstring & tangent_field_names = input.getTangentNames();
81 
82  indexerNames_ = input.getIndexerNames();
83  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
84  globalDataKey_ = input.getGlobalDataKey();
85 
86  // allocate fields
87  gatherFields_.resize(names.size());
88  for (std::size_t fd = 0; fd < names.size(); ++fd) {
89  gatherFields_[fd] =
90  PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
91  this->addEvaluatedField(gatherFields_[fd]);
92  }
93 
94  // Setup dependent tangent fields if requested
95  if (tangent_field_names.size()>0) {
96  TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
97 
98  has_tangent_fields_ = true;
99  tangentFields_.resize(gatherFields_.size());
100  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
101  tangentFields_[fd].resize(tangent_field_names[fd].size());
102  for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
103  tangentFields_[fd][i] =
104  PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
105  this->addDependentField(tangentFields_[fd][i]);
106  }
107  }
108  }
109 
110  // figure out what the first active name is
111  std::string firstName = "<none>";
112  if(names.size()>0)
113  firstName = names[0];
114 
115  std::string n = "GatherSolution (BlockedEpetra): "+firstName+" ()";
116  this->setName(n);
117 }
118 
119 // **********************************************************************
120 template<typename TRAITS,typename LO,typename GO>
122 postRegistrationSetup(typename TRAITS::SetupData d,
124 {
125  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
126 
127  indexerIds_.resize(gatherFields_.size());
128  subFieldIds_.resize(gatherFields_.size());
129 
130  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
131  // get field ID from DOF manager
132  const std::string& fieldName = indexerNames_[fd];
133 
134  indexerIds_[fd] = getFieldBlock(fieldName,indexers_);
135  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
136 
137  TEUCHOS_ASSERT(indexerIds_[fd]>=0);
138 
139  // setup the field data object
140  this->utils.setFieldData(gatherFields_[fd],fm);
141  }
142 
143  if (has_tangent_fields_) {
144  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd)
145  for (std::size_t i=0; i<tangentFields_[fd].size(); ++i)
146  this->utils.setFieldData(tangentFields_[fd][i],fm);
147  }
148 
149  indexerNames_.clear(); // Don't need this anymore
150 }
151 
152 // **********************************************************************
153 template<typename TRAITS,typename LO,typename GO>
155 preEvaluate(typename TRAITS::PreEvalData d)
156 {
157  typedef BlockedEpetraLinearObjContainer BLOC;
158 
159  using Teuchos::RCP;
160  using Teuchos::rcp_dynamic_cast;
161 
162  RCP<GlobalEvaluationData> ged;
163 
164  // first try refactored ReadOnly container
165  std::string post = useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X";
166  if(d.gedc.containsDataObject(globalDataKey_+post)) {
167  ged = d.gedc.getDataObject(globalDataKey_+post);
168 
169  RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,true);
170 
171  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(ro_ged->getGhostedVector());
172 
173  // post condition
174  TEUCHOS_TEST_FOR_EXCEPTION(x_==Teuchos::null,std::logic_error,
175  "Gather Residual: Can't find x vector in GEDC \"" << globalDataKey_ << "\" (" << post << "). "
176  "A cast failed for " << ged << ". Type is " << Teuchos::typeName(*ged));
177 
178  return;
179  }
180  else {
181  ged = d.gedc.getDataObject(globalDataKey_);
182 
183  // extract linear object container
184  RCP<const BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<const BlockedVector_ReadOnly_GlobalEvaluationData>(ged);
185  RCP<const BlockedEpetraLinearObjContainer> blockedContainer = rcp_dynamic_cast<const BLOC>(ged);
186 
187  if(ro_ged!=Teuchos::null) {
188  RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,true);
189 
190  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(ro_ged->getGhostedVector());
191  }
192  else if(blockedContainer!=Teuchos::null) {
193  if (useTimeDerivativeSolutionVector_)
194  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_dxdt());
195  else
196  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_x());
197  }
198 
199  // post condition
200  TEUCHOS_TEST_FOR_EXCEPTION(x_==Teuchos::null,std::logic_error,
201  "Gather Residual: Can't find x vector in GEDC \"" << globalDataKey_ << "\" (" << post << "). "
202  "A cast failed for " << ged << ". Type is " << Teuchos::typeName(*ged));
203  }
204 }
205 
206 // **********************************************************************
207 template<typename TRAITS,typename LO,typename GO>
209 evaluateFields(typename TRAITS::EvalData workset)
210 {
211  using Teuchos::RCP;
212  using Teuchos::ArrayRCP;
213  using Teuchos::ptrFromRef;
214  using Teuchos::rcp_dynamic_cast;
215 
216  using Thyra::VectorBase;
217  using Thyra::SpmdVectorBase;
219 
220  // for convenience pull out some objects from workset
221  std::string blockId = this->wda(workset).block_id;
222  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
223 
224  // loop over the fields to be gathered
225  Teuchos::ArrayRCP<const double> local_x;
226  for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
227 
228  PHX::MDField<ScalarT,Cell,NODE> & field = gatherFields_[fieldIndex];
229 
230  int indexerId = indexerIds_[fieldIndex];
231  int subFieldNum = subFieldIds_[fieldIndex];
232 
233  // grab local data for inputing
234  Teuchos::ArrayRCP<const double> local_x;
235  rcp_dynamic_cast<SpmdVectorBase<double> >(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(local_x));
236 
237  auto subRowIndexer = indexers_[indexerId];
238  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
239 
240  // gather operation for each cell in workset
241  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
242  LO cellLocalId = localCellIds[worksetCellIndex];
243 
244  const std::vector<int> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
245 
246  // loop over basis functions and fill the fields
247  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
248  int offset = elmtOffset[basis];
249  int lid = LIDs[offset];
250 
251  // TEUCHOS_ASSERT(indexerId==GIDs[offset].first);
252  // TEUCHOS_ASSERT(lid<local_x.size() && lid>=0);
253 
254  field(worksetCellIndex,basis) = local_x[lid];
255  }
256  }
257  }
258 }
259 
260 // **********************************************************************
261 // Specialization: Tangent
262 // **********************************************************************
263 
264 template<typename TRAITS,typename LO,typename GO>
267  const Teuchos::ParameterList& p)
268  : indexers_(indexers)
269  , has_tangent_fields_(false)
270 {
271  typedef std::vector< std::vector<std::string> > vvstring;
272 
274  input.setParameterList(p);
275 
276  const std::vector<std::string> & names = input.getDofNames();
277  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
278  const vvstring & tangent_field_names = input.getTangentNames();
279 
280  indexerNames_ = input.getIndexerNames();
281  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
282  globalDataKey_ = input.getGlobalDataKey();
283 
284  // allocate fields
285  gatherFields_.resize(names.size());
286  for (std::size_t fd = 0; fd < names.size(); ++fd) {
287  gatherFields_[fd] =
288  PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
289  this->addEvaluatedField(gatherFields_[fd]);
290  }
291 
292  // Setup dependent tangent fields if requested
293  if (tangent_field_names.size()>0) {
294  TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
295 
296  has_tangent_fields_ = true;
297  tangentFields_.resize(gatherFields_.size());
298  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
299  tangentFields_[fd].resize(tangent_field_names[fd].size());
300  for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
301  tangentFields_[fd][i] =
302  PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
303  this->addDependentField(tangentFields_[fd][i]);
304  }
305  }
306  }
307 
308  // figure out what the first active name is
309  std::string firstName = "<none>";
310  if(names.size()>0)
311  firstName = names[0];
312 
313  std::string n = "GatherSolution Tangent (BlockedEpetra): "+firstName+" ()";
314  this->setName(n);
315 }
316 
317 // **********************************************************************
318 template<typename TRAITS,typename LO,typename GO>
320 postRegistrationSetup(typename TRAITS::SetupData d,
322 {
323  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
324 
325  indexerIds_.resize(gatherFields_.size());
326  subFieldIds_.resize(gatherFields_.size());
327 
328  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
329  // get field ID from DOF manager
330  const std::string& fieldName = indexerNames_[fd];
331 
332  indexerIds_[fd] = getFieldBlock(fieldName,indexers_);
333  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
334 
335  TEUCHOS_ASSERT(indexerIds_[fd]>=0);
336 
337  // setup the field data object
338  this->utils.setFieldData(gatherFields_[fd],fm);
339  }
340 
341  if (has_tangent_fields_) {
342  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd)
343  for (std::size_t i=0; i<tangentFields_[fd].size(); ++i)
344  this->utils.setFieldData(tangentFields_[fd][i],fm);
345  }
346 
347  indexerNames_.clear(); // Don't need this anymore
348 }
349 
350 // **********************************************************************
351 template<typename TRAITS,typename LO,typename GO>
353 preEvaluate(typename TRAITS::PreEvalData d)
354 {
355  typedef BlockedEpetraLinearObjContainer BLOC;
356 
357  using Teuchos::RCP;
358  using Teuchos::rcp_dynamic_cast;
359 
360  RCP<GlobalEvaluationData> ged;
361 
362  // first try refactored ReadOnly container
363  std::string post = useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X";
364  if(d.gedc.containsDataObject(globalDataKey_+post)) {
365  ged = d.gedc.getDataObject(globalDataKey_+post);
366 
367  RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,true);
368 
369  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(ro_ged->getGhostedVector());
370 
371  return;
372  }
373  else {
374  ged = d.gedc.getDataObject(globalDataKey_);
375 
376  // extract linear object container
377  RCP<const BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<const BlockedVector_ReadOnly_GlobalEvaluationData>(ged);
378  RCP<const BlockedEpetraLinearObjContainer> blockedContainer = rcp_dynamic_cast<const BLOC>(ged);
379 
380  if(ro_ged!=Teuchos::null) {
381  RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,true);
382 
383  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(ro_ged->getGhostedVector());
384  }
385  else if(blockedContainer!=Teuchos::null) {
386  if (useTimeDerivativeSolutionVector_)
387  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_dxdt());
388  else
389  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_x());
390  }
391  }
392 
393  // post condition
394  TEUCHOS_ASSERT(x_!=Teuchos::null); // someone has to find the x_ vector
395 }
396 
397 // **********************************************************************
398 template<typename TRAITS,typename LO,typename GO>
400 evaluateFields(typename TRAITS::EvalData workset)
401 {
402  using Teuchos::RCP;
403  using Teuchos::ArrayRCP;
404  using Teuchos::ptrFromRef;
405  using Teuchos::rcp_dynamic_cast;
406 
407  using Thyra::VectorBase;
408  using Thyra::SpmdVectorBase;
410 
411  std::vector<std::pair<int,GO> > GIDs;
412  std::vector<int> LIDs;
413 
414  // for convenience pull out some objects from workset
415  std::string blockId = this->wda(workset).block_id;
416  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
417 
418  // loop over the fields to be gathered
419  Teuchos::ArrayRCP<const double> local_x;
420  for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
421 
422  PHX::MDField<ScalarT,Cell,NODE> & field = gatherFields_[fieldIndex];
423 
424  int indexerId = indexerIds_[fieldIndex];
425  int subFieldNum = subFieldIds_[fieldIndex];
426 
427  // grab local data for inputing
428  Teuchos::ArrayRCP<const double> local_x;
429  rcp_dynamic_cast<SpmdVectorBase<double> >(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(local_x));
430 
431  auto subRowIndexer = indexers_[indexerId];
432  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
433 
434  // gather operation for each cell in workset
435  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
436  LO cellLocalId = localCellIds[worksetCellIndex];
437 
438  const std::vector<int> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
439 
440  // loop over basis functions and fill the fields
441  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
442  int offset = elmtOffset[basis];
443  int lid = LIDs[offset];
444 
445  // TEUCHOS_ASSERT(indexerId==GIDs[offset].first);
446  // TEUCHOS_ASSERT(lid<local_x.size() && lid>=0);
447 
448  if (!has_tangent_fields_)
449  field(worksetCellIndex,basis) = local_x[lid];
450  else {
451  field(worksetCellIndex,basis).val() = local_x[lid];
452  for (std::size_t i=0; i<tangentFields_[fieldIndex].size(); ++i)
453  field(worksetCellIndex,basis).fastAccessDx(i) =
454  tangentFields_[fieldIndex][i](worksetCellIndex,basis).val();
455  }
456  }
457  }
458  }
459 }
460 
461 // **********************************************************************
462 // Specialization: Jacobian
463 // **********************************************************************
464 
465 template<typename TRAITS,typename LO,typename GO>
468  const Teuchos::ParameterList& p)
469  : indexers_(indexers)
470 {
471  // typedef std::vector< std::vector<std::string> > vvstring;
472 
474  input.setParameterList(p);
475 
476  const std::vector<std::string> & names = input.getDofNames();
477  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
478  //const vvstring & tangent_field_names = input.getTangentNames();
479 
480  indexerNames_ = input.getIndexerNames();
481  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
482  globalDataKey_ = input.getGlobalDataKey();
483 
484  gatherSeedIndex_ = input.getGatherSeedIndex();
485  sensitivitiesName_ = input.getSensitivitiesName();
486  disableSensitivities_ = !input.firstSensitivitiesAvailable();
487 
488  // allocate fields
489  gatherFields_.resize(names.size());
490  for (std::size_t fd = 0; fd < names.size(); ++fd) {
491  PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->functional);
492  gatherFields_[fd] = f;
493  this->addEvaluatedField(gatherFields_[fd]);
494  }
495 
496  // figure out what the first active name is
497  std::string firstName = "<none>";
498  if(names.size()>0)
499  firstName = names[0];
500 
501  // print out convenience
502  if(disableSensitivities_) {
503  std::string n = "GatherSolution (BlockedEpetra, No Sensitivities): "+firstName+" ("+PHX::typeAsString<EvalT>()+")";
504  this->setName(n);
505  }
506  else {
507  std::string n = "GatherSolution (BlockedEpetra): "+firstName+" ("+PHX::typeAsString<EvalT>()+") ";
508  this->setName(n);
509  }
510 }
511 
512 // **********************************************************************
513 template<typename TRAITS,typename LO,typename GO>
515 postRegistrationSetup(typename TRAITS::SetupData d,
517 {
518  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
519 
520  indexerIds_.resize(gatherFields_.size());
521  subFieldIds_.resize(gatherFields_.size());
522 
523  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
524  // get field ID from DOF manager
525  const std::string& fieldName = indexerNames_[fd];
526 
527  indexerIds_[fd] = getFieldBlock(fieldName,indexers_);
528  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
529 
530  TEUCHOS_ASSERT(indexerIds_[fd]>=0);
531 
532  // setup the field data object
533  this->utils.setFieldData(gatherFields_[fd],fm);
534  }
535 
536  indexerNames_.clear(); // Don't need this anymore
537 }
538 
539 template<typename TRAITS,typename LO,typename GO>
541 preEvaluate(typename TRAITS::PreEvalData d)
542 {
543  typedef BlockedEpetraLinearObjContainer BLOC;
544 
545  using Teuchos::RCP;
546  using Teuchos::rcp_dynamic_cast;
547 
548  // manage sensitivities
550  if(!disableSensitivities_) {
551  if(d.first_sensitivities_name==sensitivitiesName_)
552  applySensitivities_ = true;
553  else
554  applySensitivities_ = false;
555  }
556  else
557  applySensitivities_ = false;
558 
560 
561  RCP<GlobalEvaluationData> ged;
562 
563  // first try refactored ReadOnly container
564  std::string post = useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X";
565  if(d.gedc.containsDataObject(globalDataKey_+post)) {
566  ged = d.gedc.getDataObject(globalDataKey_+post);
567 
568  RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,true);
569 
570  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(ro_ged->getGhostedVector());
571 
572  // post condition
573  TEUCHOS_TEST_FOR_EXCEPTION(x_==Teuchos::null,std::logic_error,
574  "Gather Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ << post << "\""
575  "A cast failed for " << ged << ". Type is " << Teuchos::typeName(*ged) << std::endl);
576  }
577  else {
578  ged = d.gedc.getDataObject(globalDataKey_);
579 
580  // extract linear object container
581  RCP<const BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<const BlockedVector_ReadOnly_GlobalEvaluationData>(ged);
582  RCP<const BlockedEpetraLinearObjContainer> blockedContainer = rcp_dynamic_cast<const BLOC>(ged);
583 
584  if(ro_ged!=Teuchos::null) {
585  RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<BlockedVector_ReadOnly_GlobalEvaluationData>(ged,true);
586 
587  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(ro_ged->getGhostedVector());
588  }
589  else if(blockedContainer!=Teuchos::null) {
590  if (useTimeDerivativeSolutionVector_)
591  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_dxdt());
592  else
593  x_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_x());
594  }
595 
596  // post condition
597  TEUCHOS_TEST_FOR_EXCEPTION(x_==Teuchos::null,std::logic_error,
598  "Gather Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ << "\" (" << post << "). "
599  "A cast failed for " << ged << ". Type is " << Teuchos::typeName(*ged));
600  }
601 }
602 
603 // **********************************************************************
604 template<typename TRAITS,typename LO,typename GO>
606 evaluateFields(typename TRAITS::EvalData workset)
607 {
608  using Teuchos::RCP;
609  using Teuchos::ArrayRCP;
610  using Teuchos::ptrFromRef;
611  using Teuchos::rcp_dynamic_cast;
612 
613  using Thyra::VectorBase;
614  using Thyra::SpmdVectorBase;
616 
617  // for convenience pull out some objects from workset
618  std::string blockId = this->wda(workset).block_id;
619  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
620 
621  double seed_value = 0.0;
622  if(applySensitivities_) {
623  if (useTimeDerivativeSolutionVector_ && gatherSeedIndex_<0) {
624  seed_value = workset.alpha;
625  }
626  else if (gatherSeedIndex_<0) {
627  seed_value = workset.beta;
628  }
629  else if(!useTimeDerivativeSolutionVector_) {
630  seed_value = workset.gather_seeds[gatherSeedIndex_];
631  }
632  else {
633  TEUCHOS_ASSERT(false);
634  }
635  }
636 
637  // turn off sensitivies: this may be faster if we don't expand the term
638  // but I suspect not because anywhere it is used the full complement of
639  // sensitivies will be needed anyway.
640  if(!applySensitivities_)
641  seed_value = 0.0;
642 
643  std::vector<int> blockOffsets;
644  computeBlockOffsets(blockId,indexers_,blockOffsets);
645 
646  // NOTE: A reordering of these loops will likely improve performance
647  // The "getGIDFieldOffsets may be expensive. However the
648  // "getElementGIDs" can be cheaper. However the lookup for LIDs
649  // may be more expensive!
650 
651  int numDerivs = blockOffsets[blockOffsets.size()-1];
652 
653  // loop over the fields to be gathered
654  for(std::size_t fieldIndex=0;
655  fieldIndex<gatherFields_.size();fieldIndex++) {
656 
657  PHX::MDField<ScalarT,Cell,NODE> & field = gatherFields_[fieldIndex];
658 
659  int indexerId = indexerIds_[fieldIndex];
660  int subFieldNum = subFieldIds_[fieldIndex];
661 
662  // grab local data for inputing
663  Teuchos::ArrayRCP<const double> local_x;
664  rcp_dynamic_cast<SpmdVectorBase<double> >(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(local_x));
665 
666  auto subRowIndexer = indexers_[indexerId];
667  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
668 
669  int startBlkOffset = blockOffsets[indexerId];
670 
671  // gather operation for each cell in workset
672  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
673  LO cellLocalId = localCellIds[worksetCellIndex];
674 
675  const std::vector<int> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
676 
677  // loop over basis functions and fill the fields
678  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
679  int offset = elmtOffset[basis];
680  int lid = LIDs[offset];
681 
682  // set the value and seed the FAD object
683  field(worksetCellIndex,basis) = ScalarT(numDerivs, local_x[lid]);
684  field(worksetCellIndex,basis).fastAccessDx(startBlkOffset+offset) = seed_value;
685  }
686  }
687  }
688 }
689 
690 // **********************************************************************
691 
692 #endif
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)
PHX::MDField< const ScalarT > input
PHX::MDField< ScalarT > vector
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
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.