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_GatherSolution_BlockedEpetra_impl_hpp__
44 #define __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
45 
47 //
48 // Include Files
49 //
51 
52 // Panzer
56 #include "Panzer_PureBasis.hpp"
60 
61 // Phalanx
62 #include "Phalanx_DataLayout.hpp"
63 
64 // Teuchos
65 #include "Teuchos_Assert.hpp"
66 #include "Teuchos_FancyOStream.hpp"
67 
68 // Thyra
69 #include "Thyra_ProductVectorBase.hpp"
70 #include "Thyra_SpmdVectorBase.hpp"
71 
73 //
74 // Initializing Constructor (Residual Specialization)
75 //
77 template<typename TRAITS, typename LO, typename GO>
81  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO, int>>>&
82  indexers,
83  const Teuchos::ParameterList& p)
84  :
85  indexers_(indexers),
86  hasTangentFields_(false)
87 {
88  using panzer::PureBasis;
89  using PHX::MDField;
90  using PHX::typeAsString;
91  using std::size_t;
92  using std::string;
93  using std::vector;
94  using Teuchos::RCP;
95  using vvstring = std::vector<std::vector<std::string>>;
97  input.setParameterList(p);
98  const vector<string>& names = input.getDofNames();
99  RCP<const PureBasis> basis = input.getBasis();
100  const vvstring& tangentFieldNames = input.getTangentNames();
101  indexerNames_ = input.getIndexerNames();
102  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
103  globalDataKey_ = input.getGlobalDataKey();
104 
105  // Allocate the fields.
106  int numFields(names.size());
107  gatherFields_.resize(numFields);
108  for (int fd(0); fd < numFields; ++fd)
109  {
110  gatherFields_[fd] =
111  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
112  this->addEvaluatedField(gatherFields_[fd]);
113  } // end loop over names
114 
115  // Setup the dependent tangent fields, if requested.
116  if (tangentFieldNames.size() > 0)
117  {
118  TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size());
119  hasTangentFields_ = true;
120  tangentFields_.resize(numFields);
121  for (int fd(0); fd < numFields; ++fd)
122  {
123  int numTangentFields(tangentFieldNames[fd].size());
124  tangentFields_[fd].resize(numTangentFields);
125  for (int i(0); i < numTangentFields; ++i)
126  {
127  tangentFields_[fd][i] =
128  MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
129  basis->functional);
130  this->addDependentField(tangentFields_[fd][i]);
131  } // end loop over tangentFieldNames[fd]
132  } // end loop over gatherFields_
133  } // end if we have tangent fields
134 
135  // Figure out what the first active name is.
136  string firstName("<none>");
137  if (numFields > 0)
138  firstName = names[0];
139  string n("GatherSolution (BlockedEpetra): " + firstName + " (" +
140  typeAsString<EvalT>() + ")");
141  this->setName(n);
142 } // end of Initializing Constructor (Residual Specialization)
143 
145 //
146 // postRegistrationSetup() (Residual Specialization)
147 //
149 template<typename TRAITS, typename LO, typename GO>
150 void
151 panzer::
154  typename TRAITS::SetupData /* d */,
155  PHX::FieldManager<TRAITS>& /* fm */)
156 {
157  using std::size_t;
158  using std::string;
159  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
160  int numFields(gatherFields_.size());
161  indexerIds_.resize(numFields);
162  subFieldIds_.resize(numFields);
163  for (int fd(0); fd < numFields; ++fd)
164  {
165  // Get the field ID from the DOF manager.
166  const string& fieldName(indexerNames_[fd]);
167  indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
168  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
169  TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
170  } // end loop over gatherFields_
171  indexerNames_.clear();
172 } // end of postRegistrationSetup() (Residual Specialization)
173 
175 //
176 // preEvaluate() (Residual Specialization)
177 //
179 template<typename TRAITS, typename LO, typename GO>
180 void
181 panzer::
184  typename TRAITS::PreEvalData d)
185 {
186  using std::logic_error;
187  using std::string;
188  using Teuchos::RCP;
189  using Teuchos::rcp_dynamic_cast;
190  using Teuchos::typeName;
194  using GED = panzer::GlobalEvaluationData;
195 
196  // First try the refactored ReadOnly container.
197  RCP<GED> ged;
198  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
199  if (d.gedc->containsDataObject(globalDataKey_ + post))
200  {
201  ged = d.gedc->getDataObject(globalDataKey_ + post);
202  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
203  return;
204  } // end of the refactored ReadOnly way
205 
206  // Now try the old path.
207  {
208  ged = d.gedc->getDataObject(globalDataKey_);
209 
210  // Extract the linear object container.
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())
216  {
217  if (useTimeDerivativeSolutionVector_)
218  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
219  else // if (not useTimeDerivativeSolutionVector_)
220  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
221  TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
222  "Residual: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
223  "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
224  typeName(*ged));
225  } // end if we have a roGed or beLoc
226  } // end of the old path
227 } // end of preEvaluate() (Residual Specialization)
228 
230 //
231 // evaluateFields() (Residual Specialization)
232 //
234 template<typename TRAITS, typename LO, typename GO>
235 void
236 panzer::
239  typename TRAITS::EvalData workset)
240 {
241  using PHX::MDField;
242  using std::size_t;
243  using std::string;
244  using std::vector;
245  using Teuchos::ArrayRCP;
246  using Teuchos::ptrFromRef;
247  using Teuchos::RCP;
248  using Teuchos::rcp_dynamic_cast;
250  using Thyra::SpmdVectorBase;
251  using Thyra::VectorBase;
252 
253  // For convenience, pull out some objects from the workset.
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());
257 
258  if (x_.is_null())
259  {
260  // Loop over the fields to be gathered.
261  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
262  {
263  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
264  int indexerId(indexerIds_[fieldInd]),
265  subFieldNum(subFieldIds_[fieldInd]);
266 
267  // Grab the local data for inputing.
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());
273 
274  // Gather operation for each cell in the workset.
275  for (int cell(0); cell < numCells; ++cell)
276  {
277  LO cellLocalId = localCellIds[cell];
278  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
279 
280  // Loop over the basis functions and fill the fields.
281  for (int basis(0); basis < numBases; ++basis)
282  {
283  int offset(elmtOffset[basis]), lid(LIDs[offset]);
284  field(cell, basis) = (*xEvRoGed)[lid];
285  } // end loop over the basis functions
286  } // end loop over localCellIds
287  } // end loop over the fields to be gathered
288  }
289  else // if (not x_.is_null())
290  {
291  // Loop over the fields to be gathered.
292  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
293  {
294  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
295  int indexerId(indexerIds_[fieldInd]),
296  subFieldNum(subFieldIds_[fieldInd]);
297 
298  // Grab the local data for inputing.
299  ArrayRCP<const double> x;
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());
306 
307  // Gather operation for each cell in the workset.
308  for (int cell(0); cell < numCells; ++cell)
309  {
310  LO cellLocalId = localCellIds[cell];
311  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
312 
313  // Loop over the basis functions and fill the fields.
314  for (int basis(0); basis < numBases; ++basis)
315  {
316  int offset(elmtOffset[basis]), lid(LIDs[offset]);
317  field(cell, basis) = x[lid];
318  } // end loop over the basis functions
319  } // end loop over localCellIds
320  } // end loop over the fields to be gathered
321  } // end if (x_.is_null()) or not
322 } // end of evaluateFields() (Residual Specialization)
323 
325 //
326 // Initializing Constructor (Tangent Specialization)
327 //
329 template<typename TRAITS, typename LO, typename GO>
332  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO, int>>>&
333  indexers,
334  const Teuchos::ParameterList& p)
335  :
336  indexers_(indexers),
337  hasTangentFields_(false)
338 {
339  using panzer::PureBasis;
340  using PHX::MDField;
341  using PHX::typeAsString;
342  using std::size_t;
343  using std::string;
344  using std::vector;
345  using Teuchos::RCP;
346  using vvstring = std::vector<std::vector<std::string>>;
347  GatherSolution_Input input;
348  input.setParameterList(p);
349  const vector<string>& names = input.getDofNames();
350  RCP<const PureBasis> basis = input.getBasis();
351  const vvstring& tangentFieldNames = input.getTangentNames();
352  indexerNames_ = input.getIndexerNames();
353  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
354  globalDataKey_ = input.getGlobalDataKey();
355 
356  // Allocate the fields.
357  int numFields(names.size());
358  gatherFields_.resize(numFields);
359  for (int fd(0); fd < numFields; ++fd)
360  {
361  gatherFields_[fd] =
362  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
363  this->addEvaluatedField(gatherFields_[fd]);
364  } // end loop over names
365 
366  // Set up the dependent tangent fields, if requested.
367  if (tangentFieldNames.size() > 0)
368  {
369  TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size());
370  hasTangentFields_ = true;
371  tangentFields_.resize(numFields);
372  for (int fd(0); fd < numFields; ++fd)
373  {
374  int numTangentFields(tangentFieldNames[fd].size());
375  tangentFields_[fd].resize(numTangentFields);
376  for (int i(0); i < numTangentFields; ++i)
377  {
378  tangentFields_[fd][i] =
379  MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
380  basis->functional);
381  this->addDependentField(tangentFields_[fd][i]);
382  } // end loop over tangentFieldNames
383  } // end loop over gatherFields_
384  } // end if we have tangent fields
385 
386  // Figure out what the first active name is.
387  string firstName("<none>");
388  if (numFields > 0)
389  firstName = names[0];
390  string n("GatherSolution Tangent (BlockedEpetra): " + firstName + " (" +
391  typeAsString<EvalT>() + ")");
392  this->setName(n);
393 } // end of Initializing Constructor (Tangent Specialization)
394 
396 //
397 // postRegistrationSetup() (Tangent Specialization)
398 //
400 template<typename TRAITS, typename LO, typename GO>
401 void
404  typename TRAITS::SetupData /* d */,
405  PHX::FieldManager<TRAITS>& /* fm */)
406 {
407  using std::size_t;
408  using std::string;
409  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
410  int numFields(gatherFields_.size());
411  indexerIds_.resize(numFields);
412  subFieldIds_.resize(numFields);
413  for (int fd(0); fd < numFields; ++fd)
414  {
415  // Get the field ID from the DOF manager.
416  const string& fieldName(indexerNames_[fd]);
417  indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
418  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
419  TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
420  } // end loop over gatherFields_
421  indexerNames_.clear();
422 } // end of postRegistrationSetup() (Tangent Specialization)
423 
425 //
426 // preEvaluate() (Tangent Specialization)
427 //
429 template<typename TRAITS, typename LO, typename GO>
430 void
433  typename TRAITS::PreEvalData d)
434 {
435  using std::logic_error;
436  using std::string;
437  using Teuchos::RCP;
438  using Teuchos::rcp_dynamic_cast;
439  using Teuchos::typeName;
443  using GED = panzer::GlobalEvaluationData;
444 
445  // First try the refactored ReadOnly container.
446  RCP<GED> ged;
447  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
448  if (d.gedc->containsDataObject(globalDataKey_ + post))
449  {
450  ged = d.gedc->getDataObject(globalDataKey_ + post);
451  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
452  return;
453  } // end of the refactored ReadOnly way
454 
455  // Now try the old path.
456  {
457  ged = d.gedc->getDataObject(globalDataKey_);
458 
459  // Extract the linear object container.
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())
465  {
466  if (useTimeDerivativeSolutionVector_)
467  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
468  else // if (not useTimeDerivativeSolutionVector_)
469  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
470  TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
471  "Tangent: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
472  "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
473  typeName(*ged));
474  } // end if we have a roGed or beLoc
475  } // end of the old path
476 } // end of preEvaluate() (Tangent Specialization)
477 
479 //
480 // evaluateFields() (Tangent Specialization)
481 //
483 template<typename TRAITS, typename LO, typename GO>
484 void
487  typename TRAITS::EvalData workset)
488 {
489  using PHX::MDField;
490  using std::pair;
491  using std::size_t;
492  using std::string;
493  using std::vector;
494  using Teuchos::ArrayRCP;
495  using Teuchos::ptrFromRef;
496  using Teuchos::RCP;
497  using Teuchos::rcp_dynamic_cast;
499  using Thyra::SpmdVectorBase;
500  using Thyra::VectorBase;
501 
502  // For convenience, pull out some objects from the workset.
503  vector<pair<int, GO>> GIDs;
504  vector<int> LIDs;
505  string blockId(this->wda(workset).block_id);
506  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
507  int numFields(gatherFields_.size()), numCells(localCellIds.size());
508 
509  if (x_.is_null())
510  {
511  // Loop over the fields to be gathered.
512  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
513  {
514  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
515  int indexerId(indexerIds_[fieldInd]),
516  subFieldNum(subFieldIds_[fieldInd]);
517 
518  // Grab the local data for inputing.
519  auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
520  auto subRowIndexer = indexers_[indexerId];
521  const vector<int>& elmtOffset =
522  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
523  int numBases(elmtOffset.size());
524 
525  // Gather operation for each cell in the workset.
526  for (int cell(0); cell < numCells; ++cell)
527  {
528  LO cellLocalId = localCellIds[cell];
529  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
530 
531  // Loop over the basis functions and fill the fields.
532  for (int basis(0); basis < numBases; ++basis)
533  {
534  int offset(elmtOffset[basis]), lid(LIDs[offset]);
535  field(cell, basis) = (*xEvRoGed)[lid];
536  } // end loop over the basis functions
537  } // end loop over localCellIds
538  } // end loop over the fields to be gathered
539  }
540  else // if (not x_.is_null())
541  {
542  // Loop over the fields to be gathered.
543  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
544  {
545  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
546  int indexerId(indexerIds_[fieldInd]),
547  subFieldNum(subFieldIds_[fieldInd]);
548 
549  // Grab the local data for inputing.
550  ArrayRCP<const double> x;
551  rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
552  getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
553  auto subRowIndexer = indexers_[indexerId];
554  const vector<int>& elmtOffset =
555  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
556  int numBases(elmtOffset.size());
557 
558  // Gather operation for each cell in the workset.
559  for (int cell(0); cell < numCells; ++cell)
560  {
561  LO cellLocalId = localCellIds[cell];
562  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
563 
564  // Loop over the basis functions and fill the fields.
565  for (int basis(0); basis < numBases; ++basis)
566  {
567  int offset(elmtOffset[basis]), lid(LIDs[offset]);
568  field(cell, basis) = x[lid];
569  } // end loop over the basis functions
570  } // end loop over localCellIds
571  } // end loop over the fields to be gathered
572  } // end if (x_.is_null()) or not
573 
574  // Deal with the tangent fields.
575  if (hasTangentFields_)
576  {
577  // Loop over the fields to be gathered.
578  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
579  {
580  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
581  int indexerId(indexerIds_[fieldInd]),
582  subFieldNum(subFieldIds_[fieldInd]);
583  auto subRowIndexer = indexers_[indexerId];
584  const vector<int>& elmtOffset =
585  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
586  int numBases(elmtOffset.size());
587 
588  // Gather operation for each cell in the workset.
589  for (int cell(0); cell < numCells; ++cell)
590  {
591  // Loop over the basis functions and fill the fields.
592  for (int basis(0); basis < numBases; ++basis)
593  {
594  int numTangentFields(tangentFields_[fieldInd].size());
595  for (int i(0); i < numTangentFields; ++i)
596  field(cell, basis).fastAccessDx(i) =
597  tangentFields_[fieldInd][i](cell, basis).val();
598  } // end loop over the basis functions
599  } // end loop over localCellIds
600  } // end loop over the fields to be gathered
601  } // end if (hasTangentFields_)
602 } // end of evaluateFields() (Tangent Specialization)
603 
605 //
606 // Initializing Constructor (Jacobian Specialization)
607 //
609 template<typename TRAITS, typename LO, typename GO>
610 panzer::
613  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO, int>>>&
614  indexers,
615  const Teuchos::ParameterList& p)
616  :
617  indexers_(indexers)
618 {
619  using panzer::PureBasis;
620  using PHX::MDField;
621  using PHX::typeAsString;
622  using std::size_t;
623  using std::string;
624  using std::vector;
625  using Teuchos::RCP;
626  GatherSolution_Input input;
627  input.setParameterList(p);
628  const vector<string>& names = input.getDofNames();
629  RCP<const PureBasis> basis = input.getBasis();
630  indexerNames_ = input.getIndexerNames();
631  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
632  globalDataKey_ = input.getGlobalDataKey();
633  gatherSeedIndex_ = input.getGatherSeedIndex();
634  sensitivitiesName_ = input.getSensitivitiesName();
635  disableSensitivities_ = not input.firstSensitivitiesAvailable();
636 
637  // Allocate the fields.
638  int numFields(names.size());
639  gatherFields_.resize(numFields);
640  for (int fd(0); fd < numFields; ++fd)
641  {
642  MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
643  gatherFields_[fd] = f;
644  this->addEvaluatedField(gatherFields_[fd]);
645  } // end loop over names
646 
647  // Figure out what the first active name is.
648  string firstName("<none>"), n("GatherSolution (BlockedEpetra");
649  if (numFields > 0)
650  firstName = names[0];
651  if (disableSensitivities_)
652  n += ", No Sensitivities";
653  n += "): " + firstName + " (" + typeAsString<EvalT>() + ")";
654  this->setName(n);
655 } // end of Initializing Constructor (Jacobian Specialization)
656 
658 //
659 // postRegistrationSetup() (Jacobian Specialization)
660 //
662 template<typename TRAITS, typename LO, typename GO>
663 void
664 panzer::
667  typename TRAITS::SetupData /* d */,
668  PHX::FieldManager<TRAITS>& /* fm */)
669 {
670  using std::size_t;
671  using std::string;
672  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
673  int numFields(gatherFields_.size());
674  indexerIds_.resize(numFields);
675  subFieldIds_.resize(numFields);
676  for (int fd(0); fd < numFields; ++fd)
677  {
678  // Get the field ID from the DOF manager.
679  const string& fieldName(indexerNames_[fd]);
680  indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
681  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
682  TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
683  } // end loop over gatherFields_
684  indexerNames_.clear();
685 } // end of postRegistrationSetup() (Jacobian Specialization)
686 
688 //
689 // preEvaluate() (Jacobian Specialization)
690 //
692 template<typename TRAITS, typename LO, typename GO>
693 void
694 panzer::
697  typename TRAITS::PreEvalData d)
698 {
699  using std::endl;
700  using std::logic_error;
701  using std::string;
702  using Teuchos::RCP;
703  using Teuchos::rcp_dynamic_cast;
704  using Teuchos::typeName;
708  using GED = panzer::GlobalEvaluationData;
709 
710  // Manage sensitivities.
711  applySensitivities_ = false;
712  if ((not disableSensitivities_ ) and
713  (d.first_sensitivities_name == sensitivitiesName_))
714  applySensitivities_ = true;
715 
716  // First try the refactored ReadOnly container.
717  RCP<GED> ged;
718  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
719  if (d.gedc->containsDataObject(globalDataKey_ + post))
720  {
721  ged = d.gedc->getDataObject(globalDataKey_ + post);
722  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
723  return;
724  } // end of the refactored ReadOnly way
725 
726  // Now try the old path.
727  {
728  ged = d.gedc->getDataObject(globalDataKey_);
729 
730  // Extract the linear object container.
731  auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
732  auto beLoc = rcp_dynamic_cast<const BELOC>(ged);
733  if (not roGed.is_null())
734  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
735  else if (not beLoc.is_null())
736  {
737  if (useTimeDerivativeSolutionVector_)
738  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
739  else // if (not useTimeDerivativeSolutionVector_)
740  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
741  TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
742  "Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ <<
743  "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
744  typeName(*ged));
745  } // end if we have a roGed or beLoc
746  } // end of the old path
747 } // end of preEvaluate() (Jacobian Specialization)
748 
750 //
751 // evaluateFields() (Jacobian Specialization)
752 //
754 template<typename TRAITS, typename LO, typename GO>
755 void
756 panzer::
759  typename TRAITS::EvalData workset)
760 {
761  using PHX::MDField;
762  using std::size_t;
763  using std::string;
764  using std::vector;
765  using Teuchos::ArrayRCP;
766  using Teuchos::ptrFromRef;
767  using Teuchos::RCP;
768  using Teuchos::rcp_dynamic_cast;
770  using Thyra::SpmdVectorBase;
771  using Thyra::VectorBase;
772 
773  // For convenience, pull out some objects from the workset.
774  string blockId(this->wda(workset).block_id);
775  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
776  int numFields(gatherFields_.size()), numCells(localCellIds.size());
777  double seedValue(0);
778  if (applySensitivities_)
779  {
780  if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
781  seedValue = workset.alpha;
782  else if (gatherSeedIndex_ < 0)
783  seedValue = workset.beta;
784  else if (not useTimeDerivativeSolutionVector_)
785  seedValue = workset.gather_seeds[gatherSeedIndex_];
786  else
787  TEUCHOS_ASSERT(false);
788  } // end if (applySensitivities_)
789 
790  // Turn off sensitivies: This may be faster if we don't expand the term, but
791  // I suspect not, because anywhere it is used the full complement of
792  // sensitivies will be needed anyway.
793  if (not applySensitivities_)
794  seedValue = 0.0;
795 
796  vector<int> blockOffsets;
797  computeBlockOffsets(blockId, indexers_, blockOffsets);
798  int numDerivs(blockOffsets[blockOffsets.size() - 1]);
799 
800  // NOTE: A reordering of these loops will likely improve performance. The
801  // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
802  // can be cheaper. However the lookup for LIDs may be more expensive!
803  if (x_.is_null())
804  {
805  // Loop over the fields to be gathered.
806  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
807  {
808  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
809  int indexerId(indexerIds_[fieldInd]),
810  subFieldNum(subFieldIds_[fieldInd]);
811 
812  // Grab the local data for inputing.
813  auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
814  auto subRowIndexer = indexers_[indexerId];
815  const vector<int>& elmtOffset =
816  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
817  int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
818 
819  // Gather operation for each cell in the workset.
820  for (int cell(0); cell < numCells; ++cell)
821  {
822  LO cellLocalId = localCellIds[cell];
823  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
824 
825  // Loop over the basis functions and fill the fields.
826  for (int basis(0); basis < numBases; ++basis)
827  {
828  // Set the value and seed the FAD object.
829  int offset(elmtOffset[basis]), lid(LIDs[offset]);
830  field(cell, basis) = ScalarT(numDerivs, (*xEvRoGed)[lid]);
831  field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
832  } // end loop over the basis functions
833  } // end loop over localCellIds
834  } // end loop over the fields to be gathered
835  }
836  else // if (not x_.is_null())
837  {
838  // Loop over the fields to be gathered.
839  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
840  {
841  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
842  int indexerId(indexerIds_[fieldInd]),
843  subFieldNum(subFieldIds_[fieldInd]);
844 
845  // Grab the local data for inputing.
846  ArrayRCP<const double> x;
847  rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
848  getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
849  auto subRowIndexer = indexers_[indexerId];
850  const vector<int>& elmtOffset =
851  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
852  int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
853 
854  // Gather operation for each cell in the workset.
855  for (int cell(0); cell < numCells; ++cell)
856  {
857  LO cellLocalId = localCellIds[cell];
858  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
859 
860  // Loop over the basis functions and fill the fields.
861  for (int basis(0); basis < numBases; ++basis)
862  {
863  // Set the value and seed the FAD object.
864  int offset(elmtOffset[basis]), lid(LIDs[offset]);
865  field(cell, basis) = ScalarT(numDerivs, x[lid]);
866  field(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
867  } // end loop over the basis functions
868  } // end loop over localCellIds
869  } // end loop over the fields to be gathered
870  } // end if (x_.is_null()) or not
871 } // end of evaluateFields() (Jacobian Specialization)
872 
873 #endif // __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
std::string getGlobalDataKey() const
Name of the global evaluation data container to use for the source vector (all types) ...
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
void setParameterList(const Teuchos::ParameterList &pl)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis, std::vector< int > &blockOffsets)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
int numFields
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Post-Registration Setup.
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
Constructor.
std::string getSensitivitiesName() const
The name of the sensitivities. Enables sensitivities at "preEvaluate" time (Jacobian and Hessian) ...
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&#39;ll contribute, or in which we&#39;ll store, the result of computing this integral...
const std::vector< std::string > & getIndexerNames() const
void evaluateFields(typename TRAITS::EvalData d)
Evaluate Fields.
Description and data layouts associated with a particular basis.
const std::vector< std::vector< std::string > > & getTangentNames() const
Get the name of the tangent fields (tangent only)