Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_GatherSolution_Epetra_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_Epetra_impl_hpp__
44 #define __Panzer_GatherSolution_Epetra_impl_hpp__
45 
47 //
48 // Include Files
49 //
51 
52 // Epetra
53 #include "Epetra_Vector.h"
54 
55 // Panzer
60 #include "Panzer_PureBasis.hpp"
63 
64 // Teuchos
65 #include "Teuchos_Assert.hpp"
66 
67 // Thyra
68 #include "Thyra_SpmdVectorBase.hpp"
69 
71 //
72 // Initializing Constructor (Residual Specialization)
73 //
75 template<typename TRAITS, typename LO, typename GO>
79  const Teuchos::ParameterList& p)
80  :
81  globalIndexer_(indexer),
82  hasTangentFields_(false)
83 {
84  using panzer::PureBasis;
85  using PHX::MDField;
86  using PHX::typeAsString;
87  using std::size_t;
88  using std::vector;
89  using std::string;
90  using Teuchos::RCP;
91  using vvstring = std::vector<std::vector<std::string>>;
93  input.setParameterList(p);
94  const vector<string>& names = input.getDofNames();
95  RCP<const PureBasis> basis = input.getBasis();
96  const vvstring& tangentFieldNames = input.getTangentNames();
97  indexerNames_ = input.getIndexerNames();
98  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
99  globalDataKey_ = input.getGlobalDataKey();
100 
101  // Allocate the fields.
102  int numFields(names.size());
103  gatherFields_.resize(numFields);
104  for (int fd(0); fd < numFields; ++fd)
105  {
106  gatherFields_[fd] =
107  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
108  this->addEvaluatedField(gatherFields_[fd]);
109  } // end loop over names
110 
111  // Set up the dependent tangent fields, if requested.
112  if (tangentFieldNames.size() > 0)
113  {
114  TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size())
115  hasTangentFields_ = true;
116  tangentFields_.resize(numFields);
117  for (int fd(0); fd < numFields; ++fd)
118  {
119  int numTangentFields(tangentFieldNames[fd].size());
120  tangentFields_[fd].resize(numTangentFields);
121  for (int i(0); i < numTangentFields; ++i)
122  {
123  tangentFields_[fd][i] =
124  MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
125  basis->functional);
126  this->addDependentField(tangentFields_[fd][i]);
127  } // end loop over tangentFieldNames[fd]
128  } // end loop over gatherFields
129  } // end if (tangentFieldNames.size() > 0)
130 
131  // Figure out what the first active name is.
132  string firstName("<none>");
133  if (numFields > 0)
134  firstName = names[0];
135  string n("GatherSolution (Epetra): " + firstName + " (" +
136  typeAsString<EvalT>() + ")");
137  this->setName(n);
138 } // end of Initializing Constructor (Residual Specialization)
139 
141 //
142 // postRegistrationSetup() (Residual Specialization)
143 //
145 template<typename TRAITS, typename LO, typename GO>
146 void
149  typename TRAITS::SetupData /* d */,
150  PHX::FieldManager<TRAITS>& /* fm */)
151 {
152  using std::logic_error;
153  using std::size_t;
154  using std::string;
155  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size())
156  int numFields(gatherFields_.size());
157  fieldIds_.resize(numFields);
158  for (int fd(0); fd < numFields; ++fd)
159  {
160  // Get the field ID from the DOF manager.
161  const string& fieldName(indexerNames_[fd]);
162  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
163 
164  // This is the error return code; raise the alarm.
165  TEUCHOS_TEST_FOR_EXCEPTION(fieldIds_[fd] == -1, logic_error,
166  "GatherSolution_Epetra<Residual>: Could not find field \"" +
167  fieldName + "\" in the global indexer. ")
168  } // end loop over gatherFields_
169  indexerNames_.clear();
170 } // end of postRegistrationSetup() (Residual Specialization)
171 
173 //
174 // preEvaluate() (Residual Specialization)
175 //
177 template<typename TRAITS, typename LO, typename GO>
178 void
181  typename TRAITS::PreEvalData d)
182 {
183  using std::string;
184  using Teuchos::RCP;
185  using Teuchos::rcp_dynamic_cast;
188  using GED = panzer::GlobalEvaluationData;
189  using LOC = panzer::LinearObjContainer;
191 
192  // First try the refactored ReadOnly container.
193  RCP<GED> ged;
194  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
195  if (d.gedc->containsDataObject(globalDataKey_ + post))
196  {
197  ged = d.gedc->getDataObject(globalDataKey_ + post);
198  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
199  return;
200  } // end of the refactored ReadOnly way
201 
202  // Now try the old path.
203  ged = d.gedc->getDataObject(globalDataKey_);
204  {
205  // Try to extract the linear object container.
206  auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
207  auto locPair = rcp_dynamic_cast<LPGED>(ged);
208  if (not locPair.is_null())
209  {
210  RCP<LOC> loc = locPair->getGhostedLOC();
211  epetraContainer = rcp_dynamic_cast<ELOC>(loc);
212  } // end if (not locPair.is_null())
213  if (not epetraContainer.is_null())
214  {
215  if (useTimeDerivativeSolutionVector_)
216  x_ = epetraContainer->get_dxdt();
217  else // if (not useTimeDerivativeSolutionVector_)
218  x_ = epetraContainer->get_x();
219  return;
220  } // end if (not epetraContainer.is_null())
221  } // end of the old path
222 
223  // As a last resort, try to extract an EpetraVector_ReadOnly object. This
224  // will throw an exception if it doesn't work.
225  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
226 } // end of preEvaluate() (Residual Specialization)
227 
229 //
230 // evaluateFields() (Residual Specialization)
231 //
233 template<typename TRAITS, typename LO, typename GO>
234 void
237  typename TRAITS::EvalData workset)
238 {
239  using PHX::MDField;
240  using std::size_t;
241  using std::string;
242  using std::vector;
243  using Teuchos::ArrayRCP;
244  using Teuchos::ptrFromRef;
245  using Teuchos::RCP;
246  using Teuchos::rcp_dynamic_cast;
247  using Thyra::SpmdVectorBase;
248 
249  // For convenience, pull out some objects from the workset.
250  string blockId(this->wda(workset).block_id);
251  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
252  int numCells(localCellIds.size()), numFields(gatherFields_.size());
253 
254  // NOTE: A reordering of these loops will likely improve performance. The
255  // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
256  // can be cheaper. However the lookup for LIDs may be more expensive!
257  if (x_.is_null())
258  {
259  // Gather operation for each cell in the workset.
260  for (int cell(0); cell < numCells; ++cell)
261  {
262  size_t cellLocalId(localCellIds[cell]);
263  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
264 
265  // Loop over the fields to be gathered.
266  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
267  {
268  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
269  int fieldNum(fieldIds_[fieldInd]);
270  const vector<int>& elmtOffset =
271  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
272  int numBases(elmtOffset.size());
273 
274  // Loop over the basis functions and fill the fields.
275  for (int basis(0); basis < numBases; ++basis)
276  {
277  int offset(elmtOffset[basis]), lid(LIDs[offset]);
278  field(cell, basis) = (*xEvRoGed_)[lid];
279  } // end loop over the basis functions
280  } // end loop over the fields to be gathered
281  } // end loop over localCellIds
282  }
283  else // if (not x_.is_null())
284  {
285  // Gather operation for each cell in the workset.
286  for (int cell(0); cell < numCells; ++cell)
287  {
288  size_t cellLocalId(localCellIds[cell]);
289  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
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 fieldNum(fieldIds_[fieldInd]);
296  const vector<int>& elmtOffset =
297  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
298  int numBases(elmtOffset.size());
299 
300  // Loop over the basis functions and fill the fields.
301  for (int basis(0); basis < numBases; ++basis)
302  {
303  int offset(elmtOffset[basis]), lid(LIDs[offset]);
304  field(cell, basis) = (*x_)[lid];
305  } // end loop over the basis functions
306  } // end loop over the fields to be gathered
307  } // end loop over localCellIds
308  } // end if (x_.is_null()) or not
309 } // end of evaluateFields() (Residual Specialization)
310 
312 //
313 // Initializing Constructor (Tangent Specialization)
314 //
316 template<typename TRAITS, typename LO, typename GO>
320  const Teuchos::ParameterList& p)
321  :
322  globalIndexer_(indexer),
323  hasTangentFields_(false)
324 {
325  using panzer::PureBasis;
326  using PHX::MDField;
327  using PHX::typeAsString;
328  using std::size_t;
329  using std::string;
330  using std::vector;
331  using Teuchos::RCP;
332  using vvstring = std::vector<std::vector<std::string>>;
333  GatherSolution_Input input;
334  input.setParameterList(p);
335  const vector<string>& names = input.getDofNames();
336  RCP<const PureBasis> basis = input.getBasis();
337  const vvstring& tangentFieldNames = input.getTangentNames();
338  indexerNames_ = input.getIndexerNames();
339  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
340  globalDataKey_ = input.getGlobalDataKey();
341 
342  // Allocate the fields.
343  int numFields(names.size());
344  gatherFields_.resize(numFields);
345  for (int fd(0); fd < numFields; ++fd)
346  {
347  gatherFields_[fd] =
348  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
349  this->addEvaluatedField(gatherFields_[fd]);
350  } // end loop over names
351 
352  // Set up the dependent tangent fields, if requested.
353  if (tangentFieldNames.size() > 0)
354  {
355  TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size())
356  hasTangentFields_ = true;
357  tangentFields_.resize(numFields);
358  for (int fd(0); fd < numFields; ++fd)
359  {
360  int numTangentFields(tangentFieldNames[fd].size());
361  tangentFields_[fd].resize(numTangentFields);
362  for (int i(0); i < numTangentFields; ++i)
363  {
364  tangentFields_[fd][i] =
365  MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
366  basis->functional);
367  this->addDependentField(tangentFields_[fd][i]);
368  } // end loop over tangeng_field_names[fd]
369  } // end loop over gatherFields_
370  } // end if (tangentFieldNames.size() > 0)
371 
372  // Figure out what the first active name is.
373  string firstName("<none>");
374  if (numFields > 0)
375  firstName = names[0];
376  string n("GatherSolution (Epetra): " + firstName + " (" +
377  typeAsString<EvalT>() + ")");
378  this->setName(n);
379 } // end of Initializing Constructor (Tangent Specialization)
380 
382 //
383 // postRegistrationSetup() (Tangent Specialization)
384 //
386 template<typename TRAITS, typename LO, typename GO>
387 void
390  typename TRAITS::SetupData /* d */,
391  PHX::FieldManager<TRAITS>& /* fm */)
392 {
393  using std::logic_error;
394  using std::size_t;
395  using std::string;
396  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size())
397  int numFields(gatherFields_.size());
398  fieldIds_.resize(numFields);
399  for (int fd(0); fd < numFields; ++fd)
400  {
401  // Get the field ID from the DOF manager.
402  const string& fieldName(indexerNames_[fd]);
403  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
404 
405  // This is the error return code; raise the alarm.
406  TEUCHOS_TEST_FOR_EXCEPTION(fieldIds_[fd] == -1, logic_error,
407  "GatherSolution_Epetra<Tangent>: Could not find field \"" + fieldName
408  + "\" in the global indexer. ")
409  } // end loop over gatherFields_
410  indexerNames_.clear();
411 } // end of postRegistrationSetup() (Tangent Specialization)
412 
414 //
415 // preEvaluate() (Tangent Specialization)
416 //
418 template<typename TRAITS, typename LO, typename GO>
419 void
422  typename TRAITS::PreEvalData d)
423 {
424  using std::string;
425  using Teuchos::RCP;
426  using Teuchos::rcp_dynamic_cast;
429  using GED = panzer::GlobalEvaluationData;
430  using LOC = panzer::LinearObjContainer;
432 
433  // First try the refactored ReadOnly container.
434  RCP<GED> ged;
435  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
436  if (d.gedc->containsDataObject(globalDataKey_ + post))
437  {
438  ged = d.gedc->getDataObject(globalDataKey_ + post);
439  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
440  return;
441  } // end of the refactored ReadOnly way
442 
443  // Now try the old path.
444  ged = d.gedc->getDataObject(globalDataKey_);
445  {
446  // Try to extract the linear object container.
447  auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
448  auto locPair = rcp_dynamic_cast<LPGED>(ged);
449  if (not locPair.is_null())
450  {
451  RCP<LOC> loc = locPair->getGhostedLOC();
452  epetraContainer = rcp_dynamic_cast<ELOC>(loc);
453  } // end if (not locPair.is_null())
454  if (not epetraContainer.is_null())
455  {
456  if (useTimeDerivativeSolutionVector_)
457  x_ = epetraContainer->get_dxdt();
458  else // if (not useTimeDerivativeSolutionVector_)
459  x_ = epetraContainer->get_x();
460  return;
461  } // end if (not epetraContainer.is_null())
462  } // end of the old path
463 
464  // As a last resort, try to extract an EpetraVector_ReadOnly object. This
465  // will throw an exception if it doesn't work.
466  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
467 } // end of preEvaluate() (Tangent Specialization)
468 
470 //
471 // evaluateFields() (Tangent Specialization)
472 //
474 template<typename TRAITS, typename LO, typename GO>
475 void
478  typename TRAITS::EvalData workset)
479 {
480  using PHX::MDField;
481  using std::size_t;
482  using std::string;
483  using std::vector;
484  using Teuchos::ArrayRCP;
485  using Teuchos::ptrFromRef;
486  using Teuchos::RCP;
487  using Teuchos::rcp_dynamic_cast;
488  using Thyra::SpmdVectorBase;
489 
490  // For convenience, pull out some objects from the workset.
491  string blockId(this->wda(workset).block_id);
492  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
493  int numCells(localCellIds.size()), numFields(gatherFields_.size());
494 
495  // NOTE: A reordering of these loops will likely improve performance. The
496  // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
497  // can be cheaper. However the lookup for LIDs may be more expensive!
498  if (x_.is_null())
499  {
500  // Gather operation for each cell in the workset.
501  for (int cell(0); cell < numCells; ++cell)
502  {
503  size_t cellLocalId(localCellIds[cell]);
504  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
505 
506  // Loop over the fields to be gathered.
507  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
508  {
509  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
510  int fieldNum(fieldIds_[fieldInd]);
511  const vector<int>& elmtOffset =
512  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
513  int numBases(elmtOffset.size());
514 
515  // Loop over the basis functions and fill the fields.
516  for (int basis(0); basis < numBases; ++basis)
517  {
518  int offset(elmtOffset[basis]), lid(LIDs[offset]);
519  field(cell, basis) = (*xEvRoGed_)[lid];
520  } // end loop over the basis functions
521  } // end loop over the fields to be gathered
522  } // end loop over localCellIds
523  }
524  else // if (not x_.is_null())
525  {
526  // Gather operation for each cell in the workset.
527  for (int cell(0); cell < numCells; ++cell)
528  {
529  size_t cellLocalId(localCellIds[cell]);
530  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
531 
532  // Loop over the fields to be gathered.
533  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
534  {
535  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
536  int fieldNum(fieldIds_[fieldInd]);
537  const vector<int>& elmtOffset =
538  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
539  int numBases(elmtOffset.size());
540 
541  // Loop over the basis functions and fill the fields.
542  for (int basis(0); basis < numBases; ++basis)
543  {
544  int offset(elmtOffset[basis]), lid(LIDs[offset]);
545  field(cell, basis) = (*x_)[lid];
546  } // end loop over the basis functions
547  } // end loop over the fields to be gathered
548  } // end loop over localCellIds
549  } // end if (x_.is_null()) or not
550 
551  // Deal with the tangent fields.
552  if (hasTangentFields_)
553  {
554  // Loop over the cells in the workset.
555  for (int cell(0); cell < numCells; ++cell)
556  {
557  // Loop over the fields to be gathered.
558  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
559  {
560  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
561  int fieldNum(fieldIds_[fieldInd]);
562  const vector<int>& elmtOffset =
563  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
564  int numBases(elmtOffset.size());
565 
566  // Loop over the basis functions.
567  for (int basis(0); basis < numBases; ++basis)
568  {
569  // Loop over the tangent fields and fill them in.
570  int numTangentFields(tangentFields_[fieldInd].size());
571  for (int i(0); i < numTangentFields; ++i)
572  field(cell, basis).fastAccessDx(i) =
573  tangentFields_[fieldInd][i](cell, basis).val();
574  } // end loop over the basis functions
575  } // end loop over the fields to be gathered
576  } // end loop over the cells in the workset
577  } // end if (hasTangentFields_)
578 } // end of evaluateFields() (Tangent Specialization)
579 
581 //
582 // Initializing Constructor (Jacobian Specialization)
583 //
585 template<typename TRAITS, typename LO, typename GO>
589  const Teuchos::ParameterList& p)
590  :
591  globalIndexer_(indexer)
592 {
593  using panzer::PureBasis;
594  using PHX::MDField;
595  using PHX::typeAsString;
596  using std::size_t;
597  using std::string;
598  using std::vector;
599  using Teuchos::RCP;
600  GatherSolution_Input input;
601  input.setParameterList(p);
602  const vector<string>& names = input.getDofNames();
603  RCP<const PureBasis> basis = input.getBasis();
604  indexerNames_ = input.getIndexerNames();
605  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
606  globalDataKey_ = input.getGlobalDataKey();
607  gatherSeedIndex_ = input.getGatherSeedIndex();
608  sensitivitiesName_ = input.getSensitivitiesName();
609  disableSensitivities_ = not input.firstSensitivitiesAvailable();
610 
611  // Allocate the fields.
612  int numFields(names.size());
613  gatherFields_.resize(numFields);
614  for (int fd(0); fd < numFields; ++fd)
615  {
616  MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
617  gatherFields_[fd] = f;
618  this->addEvaluatedField(gatherFields_[fd]);
619  } // end loop over names
620 
621  // Figure out what the first active name is.
622  string firstName("<none>"), n("GatherSolution (Epetra");
623  if (numFields > 0)
624  firstName = names[0];
625  if (disableSensitivities_)
626  n += ", No Sensitivities";
627  n += "): " + firstName + " (" + typeAsString<EvalT>() + ")";
628  this->setName(n);
629 } // end of Initializing Constructor (Jacobian Specialization)
630 
632 //
633 // postRegistrationSetup() (Jacobian Specialization)
634 //
636 template<typename TRAITS, typename LO, typename GO>
637 void
640  typename TRAITS::SetupData /* d */,
641  PHX::FieldManager<TRAITS>& /* fm */)
642 {
643  using std::logic_error;
644  using std::size_t;
645  using std::string;
646  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size())
647  int numFields(gatherFields_.size());
648  fieldIds_.resize(numFields);
649  for (int fd(0); fd < numFields; ++fd)
650  {
651  // Get the field ID from the DOF manager.
652  const string& fieldName(indexerNames_[fd]);
653  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
654 
655  // This is the error return code; raise the alarm.
656  TEUCHOS_TEST_FOR_EXCEPTION(fieldIds_[fd] == -1, logic_error,
657  "GatherSolution_Epetra<Jacobian>: Could not find field \"" +
658  fieldName + "\" in the global indexer. ")
659  } // end loop over gatherFields_
660  indexerNames_.clear();
661 } // end of postRegistrationSetup() (Jacobian Specialization)
662 
664 //
665 // preEvaluate() (Jacobian Specialization)
666 //
668 template<typename TRAITS, typename LO, typename GO>
669 void
672  typename TRAITS::PreEvalData d)
673 {
674  using std::string;
675  using Teuchos::RCP;
676  using Teuchos::rcp_dynamic_cast;
679  using GED = panzer::GlobalEvaluationData;
680  using LOC = panzer::LinearObjContainer;
682 
683  // Manage sensitivities.
684  applySensitivities_ = false;
685  if ((not disableSensitivities_ ) and
686  (d.first_sensitivities_name == sensitivitiesName_))
687  applySensitivities_ = true;
688 
689  // First try the refactored ReadOnly container.
690  RCP<GED> ged;
691  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
692  if (d.gedc->containsDataObject(globalDataKey_ + post))
693  {
694  ged = d.gedc->getDataObject(globalDataKey_ + post);
695  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
696  return;
697  } // end of the refactored ReadOnly way
698 
699  // Now try the old path.
700  ged = d.gedc->getDataObject(globalDataKey_);
701  {
702  // Try to extract the linear object container.
703  auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
704  auto locPair = rcp_dynamic_cast<LPGED>(ged);
705  if (not locPair.is_null())
706  {
707  RCP<LOC> loc = locPair->getGhostedLOC();
708  epetraContainer = rcp_dynamic_cast<ELOC>(loc);
709  } // end if (not locPair.is_null())
710  if (not epetraContainer.is_null())
711  {
712  if (useTimeDerivativeSolutionVector_)
713  x_ = epetraContainer->get_dxdt();
714  else // if (not useTimeDerivativeSolutionVector_)
715  x_ = epetraContainer->get_x();
716  return;
717  } // end if (not epetraContainer.is_null())
718  } // end of the old path
719 
720  // As a last resort, try to extract an EpetraVector_ReadOnly object. This
721  // will throw an exception if it doesn't work.
722  xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
723 } // end of preEvaluate() (Jacobian Specialization)
724 
726 //
727 // evaluateFields() (Jacobian Specialization)
728 //
730 template<typename TRAITS, typename LO, typename GO>
731 void
734  typename TRAITS::EvalData workset)
735 {
736  using PHX::MDField;
737  using std::size_t;
738  using std::string;
739  using std::vector;
740  using Teuchos::ArrayRCP;
741  using Teuchos::ptrFromRef;
742  using Teuchos::RCP;
743  using Teuchos::rcp_dynamic_cast;
744  using Thyra::SpmdVectorBase;
745 
746  // For convenience, pull out some objects from the workset.
747  string blockId(this->wda(workset).block_id);
748  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
749  int numFields(gatherFields_.size()), numCells(localCellIds.size());
750 
751  // Set a sensitivity seed value.
752  double seedValue(0);
753  if (applySensitivities_)
754  {
755  if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
756  seedValue = workset.alpha;
757  else if (gatherSeedIndex_ < 0)
758  seedValue = workset.beta;
759  else if (not useTimeDerivativeSolutionVector_)
760  seedValue = workset.gather_seeds[gatherSeedIndex_];
761  else
762  TEUCHOS_ASSERT(false)
763  } // end if (applySensitivities_)
764 
765  // Turn off sensitivies. This may be faster if we don't expand the term, but
766  // I suspect not, because anywhere it is used the full complement of
767  // sensitivies will be needed anyway.
768  if (not applySensitivities_)
769  seedValue = 0.0;
770 
771  // Interface worksets handle DOFs from two element blocks. The derivative
772  // offset for the other element block must be shifted by the derivative side
773  // of my element block.
774  int dos(0);
775  if (this->wda.getDetailsIndex() == 1)
776  {
777  // Get the DOF count for my element block.
778  dos = globalIndexer_->getElementBlockGIDCount(workset.details(0).block_id);
779  } // end if (this->wda.getDetailsIndex() == 1)
780 
781  // NOTE: A reordering of these loops will likely improve performance. The
782  // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
783  // can be cheaper. However the lookup for LIDs may be more expensive!
784  if (x_.is_null())
785  {
786  // Loop over the fields to be gathered.
787  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
788  {
789  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
790  int fieldNum(fieldIds_[fieldInd]);
791  const vector<int>& elmtOffset =
792  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
793  int numBases(elmtOffset.size());
794 
795  // Gather operation for each cell in the workset.
796  for (int cell(0); cell < numCells; ++cell)
797  {
798  size_t cellLocalId(localCellIds[cell]);
799  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
800 
801  // Loop over the basis functions and fill the fields.
802  for (int basis(0); basis < numBases; ++basis)
803  {
804  int offset(elmtOffset[basis]), lid(LIDs[offset]);
805  field(cell, basis) = (*xEvRoGed_)[lid];
806  } // end loop over the basis functions
807  } // end loop over localCellIds
808  } // end loop over the fields to be gathered
809  }
810  else // if (not x_.is_null())
811  {
812  // Loop over the fields to be gathered.
813  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
814  {
815  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
816  int fieldNum(fieldIds_[fieldInd]);
817  const vector<int>& elmtOffset =
818  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
819  int numBases(elmtOffset.size());
820 
821  // Gather operation for each cell in the workset.
822  for (int cell(0); cell < numCells; ++cell)
823  {
824  size_t cellLocalId(localCellIds[cell]);
825  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
826 
827  // Loop over the basis functions and fill the fields.
828  for (int basis(0); basis < numBases; ++basis)
829  {
830  int offset(elmtOffset[basis]), lid(LIDs[offset]);
831  field(cell, basis) = (*x_)[lid];
832  } // end loop over the basis functions
833  } // end loop over localCellIds
834  } // end loop over the fields to be gathered
835  } // end if (x_.is_null()) or not
836 
837  // Deal with the sensitivities.
838  if (applySensitivities_)
839  {
840  // Loop over the fields to be gathered.
841  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
842  {
843  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
844  int fieldNum(fieldIds_[fieldInd]);
845  const vector<int>& elmtOffset =
846  globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
847  int numBases(elmtOffset.size());
848 
849  // Gather operation for each cell in the workset.
850  for (int cell(0); cell < numCells; ++cell)
851  {
852  // Loop over the basis functions.
853  for (int basis(0); basis < numBases; ++basis)
854  {
855  // Seed the FAD object.
856  int offset(elmtOffset[basis]);
857  field(cell, basis).fastAccessDx(dos + offset) = seedValue;
858  } // end loop over the basis functions
859  } // end loop over localCellIds
860  } // end loop over the fields to be gathered
861  } // end if (applySensitivities_)
862 } // end of evaluateFields() (Jacobian Specialization)
863 
864 #endif // __Panzer_GatherSolution_Epetra_impl_hpp__
const std::vector< std::vector< std::string > > & getTangentNames() const
Get the name of the tangent fields (tangent only)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
void setParameterList(const Teuchos::ParameterList &pl)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::string getGlobalDataKey() const
Name of the global evaluation data container to use for the source vector (all types) ...
int numFields
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
This class provides a boundary exchange communication mechanism for vectors.
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
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
Description and data layouts associated with a particular basis.
#define TEUCHOS_ASSERT(assertion_test)
std::string getSensitivitiesName() const
The name of the sensitivities. Enables sensitivities at &quot;preEvaluate&quot; time (Jacobian and Hessian) ...