Panzer  Version of the Day
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_GATHER_SOLUTION_EPETRA_IMPL_HPP
44 #define PANZER_GATHER_SOLUTION_EPETRA_IMPL_HPP
45 
46 #include "Teuchos_Assert.hpp"
47 #include "Phalanx_DataLayout.hpp"
48 
50 #include "Panzer_PureBasis.hpp"
56 
57 #include "Teuchos_FancyOStream.hpp"
58 
59 #include "Epetra_Vector.h"
60 #include "Epetra_Map.h"
61 
62 // **********************************************************************
63 // Specialization: Residual
64 // **********************************************************************
65 
66 template<typename TRAITS,typename LO,typename GO>
69  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & indexer,
70  const Teuchos::ParameterList& p)
71  : globalIndexer_(indexer)
72  , has_tangent_fields_(false)
73 {
74  typedef std::vector< std::vector<std::string> > vvstring;
75 
77  input.setParameterList(p);
78 
79  const std::vector<std::string> & names = input.getDofNames();
80  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
81  const vvstring & tangent_field_names = input.getTangentNames();
82 
83  indexerNames_ = input.getIndexerNames();
84  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
85  globalDataKey_ = input.getGlobalDataKey();
86 
87  // allocate fields
88  gatherFields_.resize(names.size());
89  for (std::size_t fd = 0; fd < names.size(); ++fd) {
90  gatherFields_[fd] =
91  PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
92  this->addEvaluatedField(gatherFields_[fd]);
93  }
94 
95  // Setup dependent tangent fields if requested
96  if (tangent_field_names.size()>0) {
97  TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
98 
99  has_tangent_fields_ = true;
100  tangentFields_.resize(gatherFields_.size());
101  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
102  tangentFields_[fd].resize(tangent_field_names[fd].size());
103  for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
104  tangentFields_[fd][i] =
105  PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
106  this->addDependentField(tangentFields_[fd][i]);
107  }
108  }
109  }
110 
111  // figure out what the first active name is
112  std::string firstName = "<none>";
113  if(names.size()>0)
114  firstName = names[0];
115 
116  std::string n = "GatherSolution (Epetra): "+firstName+" (Residual)";
117  this->setName(n);
118 }
119 
120 // **********************************************************************
121 template<typename TRAITS,typename LO,typename GO>
123 postRegistrationSetup(typename TRAITS::SetupData d,
125 {
126  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
127 
128  fieldIds_.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  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
134 
135  // this is the error return code, raise the alarm
136  if(fieldIds_[fd]==-1) {
137  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
138  "GatherSolution_Epetra<Residual>: Could not find field \"" + fieldName + "\" in the global indexer. ");
139  // wouldn't it be nice to print more information???
140  }
141 
142  // setup the field data object
143  this->utils.setFieldData(gatherFields_[fd],fm);
144  }
145 
146  if (has_tangent_fields_) {
147  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd)
148  for (std::size_t i=0; i<tangentFields_[fd].size(); ++i)
149  this->utils.setFieldData(tangentFields_[fd][i],fm);
150  }
151 
152  indexerNames_.clear(); // Don't need this anymore
153 }
154 
155 // **********************************************************************
156 template<typename TRAITS,typename LO,typename GO>
158 preEvaluate(typename TRAITS::PreEvalData d)
159 {
160  using Teuchos::RCP;
161  using Teuchos::rcp_dynamic_cast;
162 
163  RCP<GlobalEvaluationData> ged;
164 
165  // first try refactored ReadOnly container
166  std::string post = useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X";
167  if(d.gedc.containsDataObject(globalDataKey_+post)) {
168  ged = d.gedc.getDataObject(globalDataKey_+post);
169  RCP<EpetraVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<EpetraVector_ReadOnly_GlobalEvaluationData>(ged,true);
170 
171  x_ = ro_ged->getGhostedVector_Epetra();
172 
173  return;
174  }
175 
176  // now try old path
177  ged = d.gedc.getDataObject(globalDataKey_);
178 
179  // try to extract linear object container
180  {
181  RCP<EpetraLinearObjContainer> epetraContainer = rcp_dynamic_cast<EpetraLinearObjContainer>(ged);
182  RCP<LOCPair_GlobalEvaluationData> loc_pair = rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(ged);
183 
184  if(loc_pair!=Teuchos::null) {
185  Teuchos::RCP<LinearObjContainer> loc = loc_pair->getGhostedLOC();
186  // extract linear object container
187  epetraContainer = rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
188  }
189 
190  if(epetraContainer!=Teuchos::null) {
191  if (useTimeDerivativeSolutionVector_)
192  x_ = epetraContainer->get_dxdt();
193  else
194  x_ = epetraContainer->get_x();
195 
196  return; // epetraContainer was found
197  }
198  }
199 
200  // try to extract an EpetraVector_ReadOnly object (this is the last resort!, it throws if not found)
201  {
202  RCP<EpetraVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<EpetraVector_ReadOnly_GlobalEvaluationData>(ged,true);
203 
204  x_ = ro_ged->getGhostedVector_Epetra();
205  }
206 }
207 
208 // **********************************************************************
209 template<typename TRAITS,typename LO,typename GO>
211 evaluateFields(typename TRAITS::EvalData workset)
212 {
213  std::vector<int> LIDs;
214 
215  // for convenience pull out some objects from workset
216  std::string blockId = this->wda(workset).block_id;
217  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
218 
219  // NOTE: A reordering of these loops will likely improve performance
220  // The "getGIDFieldOffsets may be expensive. However the
221  // "getElementGIDs" can be cheaper. However the lookup for LIDs
222  // may be more expensive!
223 
224  Epetra_Vector & x = *x_;
225 
226  // gather operation for each cell in workset
227  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
228  std::size_t cellLocalId = localCellIds[worksetCellIndex];
229 
230  LIDs = globalIndexer_->getElementLIDs(cellLocalId);
231 
232  // loop over the fields to be gathered
233  for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
234  int fieldNum = fieldIds_[fieldIndex];
235  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
236 
237  // loop over basis functions and fill the fields
238  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
239  int offset = elmtOffset[basis];
240  int lid = LIDs[offset];
241  (gatherFields_[fieldIndex])(worksetCellIndex,basis) = x[lid];
242  }
243  }
244  }
245 }
246 
247 // **********************************************************************
248 // Specialization: Tangent
249 // **********************************************************************
250 
251 template<typename TRAITS,typename LO,typename GO>
254  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & indexer,
255  const Teuchos::ParameterList& p)
256  : globalIndexer_(indexer)
257  , has_tangent_fields_(false)
258 {
259  typedef std::vector< std::vector<std::string> > vvstring;
260 
262  input.setParameterList(p);
263 
264  const std::vector<std::string> & names = input.getDofNames();
265  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
266  const vvstring & tangent_field_names = input.getTangentNames();
267 
268  indexerNames_ = input.getIndexerNames();
269  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
270  globalDataKey_ = input.getGlobalDataKey();
271 
272  // allocate fields
273  gatherFields_.resize(names.size());
274  for (std::size_t fd = 0; fd < names.size(); ++fd) {
275  gatherFields_[fd] =
276  PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
277  this->addEvaluatedField(gatherFields_[fd]);
278  }
279 
280  // Setup dependent tangent fields if requested
281  if (tangent_field_names.size()>0) {
282  TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
283 
284  has_tangent_fields_ = true;
285  tangentFields_.resize(gatherFields_.size());
286  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
287  tangentFields_[fd].resize(tangent_field_names[fd].size());
288  for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
289  tangentFields_[fd][i] =
290  PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
291  this->addDependentField(tangentFields_[fd][i]);
292  }
293  }
294  }
295 
296  // figure out what the first active name is
297  std::string firstName = "<none>";
298  if(names.size()>0)
299  firstName = names[0];
300 
301  std::string n = "GatherSolution (Epetra): "+firstName+" (Tangent)";
302  this->setName(n);
303 }
304 
305 // **********************************************************************
306 template<typename TRAITS,typename LO,typename GO>
308 postRegistrationSetup(typename TRAITS::SetupData d,
310 {
311  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
312 
313  fieldIds_.resize(gatherFields_.size());
314 
315  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
316  // get field ID from DOF manager
317  const std::string& fieldName = indexerNames_[fd];
318  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
319 
320  // this is the error return code, raise the alarm
321  if(fieldIds_[fd]==-1) {
322  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
323  "GatherSolution_Epetra<Tangent>: Could not find field \"" + fieldName + "\" in the global indexer. ");
324  // wouldn't it be nice to print more information???
325  }
326 
327  // setup the field data object
328  this->utils.setFieldData(gatherFields_[fd],fm);
329  }
330 
331  if (has_tangent_fields_) {
332  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd)
333  for (std::size_t i=0; i<tangentFields_[fd].size(); ++i)
334  this->utils.setFieldData(tangentFields_[fd][i],fm);
335  }
336 
337  indexerNames_.clear();
338 }
339 
340 // **********************************************************************
341 template<typename TRAITS,typename LO,typename GO>
343 preEvaluate(typename TRAITS::PreEvalData d)
344 {
345  using Teuchos::RCP;
346  using Teuchos::rcp_dynamic_cast;
347 
348  RCP<GlobalEvaluationData> ged;
349 
350  // first try refactored ReadOnly container
351  std::string post = useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X";
352  if(d.gedc.containsDataObject(globalDataKey_+post)) {
353  ged = d.gedc.getDataObject(globalDataKey_+post);
354  RCP<EpetraVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<EpetraVector_ReadOnly_GlobalEvaluationData>(ged,true);
355 
356  x_ = ro_ged->getGhostedVector_Epetra();
357 
358  return;
359  }
360 
361  ged = d.gedc.getDataObject(globalDataKey_);
362 
363  // try to extract linear object container
364  {
365  RCP<EpetraLinearObjContainer> epetraContainer = rcp_dynamic_cast<EpetraLinearObjContainer>(ged);
366  RCP<LOCPair_GlobalEvaluationData> loc_pair = rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(ged);
367 
368  if(loc_pair!=Teuchos::null) {
369  Teuchos::RCP<LinearObjContainer> loc = loc_pair->getGhostedLOC();
370  // extract linear object container
371  epetraContainer = rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
372  }
373 
374  if(epetraContainer!=Teuchos::null) {
375  if (useTimeDerivativeSolutionVector_)
376  x_ = epetraContainer->get_dxdt();
377  else
378  x_ = epetraContainer->get_x();
379 
380  return; // epetraContainer was found
381  }
382  }
383 
384  // try to extract an EpetraVector_ReadOnly object (this is the last resort!, it throws if not found)
385  {
386  RCP<EpetraVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<EpetraVector_ReadOnly_GlobalEvaluationData>(ged,true);
387 
388  x_ = ro_ged->getGhostedVector_Epetra();
389  }
390 }
391 
392 // **********************************************************************
393 template<typename TRAITS,typename LO,typename GO>
395 evaluateFields(typename TRAITS::EvalData workset)
396 {
397  std::vector<int> LIDs;
398 
399  // for convenience pull out some objects from workset
400  std::string blockId = this->wda(workset).block_id;
401  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
402 
403  // NOTE: A reordering of these loops will likely improve performance
404  // The "getGIDFieldOffsets may be expensive. However the
405  // "getElementGIDs" can be cheaper. However the lookup for LIDs
406  // may be more expensive!
407 
408  Epetra_Vector & x = *x_;
409 
410  // gather operation for each cell in workset
411  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
412  std::size_t cellLocalId = localCellIds[worksetCellIndex];
413 
414  LIDs = globalIndexer_->getElementLIDs(cellLocalId);
415 
416  // loop over the fields to be gathered
417  for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
418  int fieldNum = fieldIds_[fieldIndex];
419  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
420 
421  // loop over basis functions and fill the fields
422  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
423  int offset = elmtOffset[basis];
424  int lid = LIDs[offset];
425  if (!has_tangent_fields_)
426  (gatherFields_[fieldIndex])(worksetCellIndex,basis) = x[lid];
427  else {
428  (gatherFields_[fieldIndex])(worksetCellIndex,basis).val() = x[lid];
429  for (std::size_t i=0; i<tangentFields_[fieldIndex].size(); ++i)
430  (gatherFields_[fieldIndex])(worksetCellIndex,basis).fastAccessDx(i) =
431  tangentFields_[fieldIndex][i](worksetCellIndex,basis).val();
432  }
433  }
434  }
435  }
436 }
437 
438 // **********************************************************************
439 // Specialization: Jacobian
440 // **********************************************************************
441 
442 template<typename TRAITS,typename LO,typename GO>
445  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & indexer,
446  const Teuchos::ParameterList& p)
447  : globalIndexer_(indexer)
448 {
449  // typedef std::vector< std::vector<std::string> > vvstring;
450 
452  input.setParameterList(p);
453 
454  const std::vector<std::string> & names = input.getDofNames();
455  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
456  //const vvstring & tangent_field_names = input.getTangentNames();
457 
458  indexerNames_ = input.getIndexerNames();
459  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
460  globalDataKey_ = input.getGlobalDataKey();
461 
462  gatherSeedIndex_ = input.getGatherSeedIndex();
463  sensitivitiesName_ = input.getSensitivitiesName();
464  disableSensitivities_ = !input.firstSensitivitiesAvailable();
465 
466  // allocate fields
467  gatherFields_.resize(names.size());
468  for (std::size_t fd = 0; fd < names.size(); ++fd) {
469  PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->functional);
470  gatherFields_[fd] = f;
471  this->addEvaluatedField(gatherFields_[fd]);
472  }
473 
474  // figure out what the first active name is
475  std::string firstName = "<none>";
476  if(names.size()>0)
477  firstName = names[0];
478 
479  // print out convenience
480  if(disableSensitivities_) {
481  std::string n = "GatherSolution (Epetra, No Sensitivities): "+firstName+" (Jacobian)";
482  this->setName(n);
483  }
484  else {
485  std::string n = "GatherSolution (Epetra): "+firstName+" ("+PHX::typeAsString<EvalT>()+") ";
486  this->setName(n);
487  }
488 }
489 
490 // **********************************************************************
491 template<typename TRAITS,typename LO,typename GO>
493 postRegistrationSetup(typename TRAITS::SetupData d,
495 {
496  // globalIndexer_ = d.globalIndexer_;
497  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
498 
499  fieldIds_.resize(gatherFields_.size());
500 
501  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
502  // get field ID from DOF manager
503  const std::string& fieldName = indexerNames_[fd];
504  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
505 
506  // this is the error return code, raise the alarm
507  if(fieldIds_[fd]==-1) {
508  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
509  "GatherSolution_Epetra<Jacobian>: Could not find field \"" + fieldName + "\" in the global indexer. ");
510  // wouldn't it be nice to print more information???
511  }
512 
513  // setup the field data object
514  this->utils.setFieldData(gatherFields_[fd],fm);
515  }
516 
517  indexerNames_.clear(); // Don't need this anymore
518 }
519 
520 // **********************************************************************
521 template<typename TRAITS,typename LO,typename GO>
523 preEvaluate(typename TRAITS::PreEvalData d)
524 {
525  using Teuchos::RCP;
526  using Teuchos::rcp_dynamic_cast;
527 
528  // manage sensitivities
530  if(!disableSensitivities_) {
531  if(d.first_sensitivities_name==sensitivitiesName_)
532  applySensitivities_ = true;
533  else
534  applySensitivities_ = false;
535  }
536  else
537  applySensitivities_ = false;
538 
540 
541  RCP<GlobalEvaluationData> ged;
542 
543  // first try refactored ReadOnly container
544  std::string post = useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X";
545  if(d.gedc.containsDataObject(globalDataKey_+post)) {
546  ged = d.gedc.getDataObject(globalDataKey_+post);
547 
548  RCP<EpetraVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<EpetraVector_ReadOnly_GlobalEvaluationData>(ged,true);
549 
550  x_ = ro_ged->getGhostedVector_Epetra();
551 
552  return;
553  }
554 
555  ged = d.gedc.getDataObject(globalDataKey_);
556 
557  // try to extract linear object container
558  {
559  RCP<EpetraLinearObjContainer> epetraContainer = rcp_dynamic_cast<EpetraLinearObjContainer>(ged);
560  RCP<LOCPair_GlobalEvaluationData> loc_pair = rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(ged);
561 
562  if(loc_pair!=Teuchos::null) {
563  Teuchos::RCP<LinearObjContainer> loc = loc_pair->getGhostedLOC();
564  // extract linear object container
565  epetraContainer = rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
566  }
567 
568  if(epetraContainer!=Teuchos::null) {
569  if (useTimeDerivativeSolutionVector_)
570  x_ = epetraContainer->get_dxdt();
571  else
572  x_ = epetraContainer->get_x();
573 
574  return; // epetraContainer was found
575  }
576  }
577 
578  // try to extract an EpetraVector_ReadOnly object (this is the last resort!, it throws if not found)
579  {
580  RCP<EpetraVector_ReadOnly_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<EpetraVector_ReadOnly_GlobalEvaluationData>(ged,true);
581 
582  x_ = ro_ged->getGhostedVector_Epetra();
583  }
584 }
585 
586 // **********************************************************************
587 template<typename TRAITS,typename LO,typename GO>
589 evaluateFields(typename TRAITS::EvalData workset)
590 {
591  // for convenience pull out some objects from workset
592  std::string blockId = this->wda(workset).block_id;
593  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
594 
595  double seed_value = 0.0;
596  if(applySensitivities_) {
597  if (useTimeDerivativeSolutionVector_ && gatherSeedIndex_<0) {
598  seed_value = workset.alpha;
599  }
600  else if (gatherSeedIndex_<0) {
601  seed_value = workset.beta;
602  }
603  else if(!useTimeDerivativeSolutionVector_) {
604  seed_value = workset.gather_seeds[gatherSeedIndex_];
605  }
606  else {
607  TEUCHOS_ASSERT(false);
608  }
609  }
610 
611  Epetra_Vector & x = *x_;
612 
613  // turn off sensitivies: this may be faster if we don't expand the term
614  // but I suspect not because anywhere it is used the full complement of
615  // sensitivies will be needed anyway.
616  if(!applySensitivities_)
617  seed_value = 0.0;
618 
619  // NOTE: A reordering of these loops will likely improve performance
620  // The "getGIDFieldOffsets may be expensive. However the
621  // "getElementGIDs" can be cheaper. However the lookup for LIDs
622  // may be more expensive!
623 
624  // Interface worksets handle DOFs from two element blocks. The derivative
625  // offset for the other element block must be shifted by the derivative side
626  // of my element block.
627  int dos = 0;
628  if (this->wda.getDetailsIndex() == 1) {
629  // Get the DOF count for my element block.
630  dos = globalIndexer_->getElementBlockGIDCount(workset.details(0).block_id);
631  }
632 
633  // loop over the fields to be gathered
634  for(std::size_t fieldIndex=0;
635  fieldIndex<gatherFields_.size();fieldIndex++) {
636  PHX::MDField<ScalarT,Cell,NODE> & field = gatherFields_[fieldIndex];
637  int fieldNum = fieldIds_[fieldIndex];
638  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
639 
640  // gather operation for each cell in workset
641  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
642  std::size_t cellLocalId = localCellIds[worksetCellIndex];
643 
644  const std::vector<int> & LIDs = globalIndexer_->getElementLIDs(cellLocalId);
645 
646  if(!applySensitivities_) {
647  // loop over basis functions and fill the fields
648  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
649  int offset = elmtOffset[basis];
650  int lid = LIDs[offset];
651 
652  // set the value and seed the FAD object
653  field(worksetCellIndex,basis) = x[lid];
654  }
655  }
656  else {
657  // loop over basis functions and fill the fields
658  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
659  int offset = elmtOffset[basis];
660  int lid = LIDs[offset];
661 
662  // set the value and seed the FAD object
663  field(worksetCellIndex,basis).val() = x[lid];
664  field(worksetCellIndex,basis).fastAccessDx(dos + offset) = seed_value;
665  }
666  }
667  }
668  }
669 }
670 
671 // **********************************************************************
672 
673 #endif
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
PHX::MDField< const ScalarT > input
Teuchos::RCP< Epetra_Vector > getGhostedVector_Epetra() const
Get the ghosted vector (Epetra version)
PHX::MDField< const ScalarT, Cell, IP > field
Teuchos::RCP< LinearObjContainer > getGhostedLOC() const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
const Teuchos::RCP< Epetra_Vector > get_dxdt() const