Panzer Version of the Day
Loading...
Searching...
No Matches
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//
75template<typename TRAITS, typename LO, typename GO>
78 const Teuchos::RCP<const panzer::GlobalIndexer>& indexer,
79 const Teuchos::ParameterList& p)
80 :
81 globalIndexer_(indexer),
82 hasTangentFields_(false)
83{
85 using PHX::MDField;
86 using PHX::print;
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 + " (Residual)");
136 this->setName(n);
137} // end of Initializing Constructor (Residual Specialization)
138
140//
141// postRegistrationSetup() (Residual Specialization)
142//
144template<typename TRAITS, typename LO, typename GO>
145void
148 typename TRAITS::SetupData /* d */,
150{
151 using std::logic_error;
152 using std::size_t;
153 using std::string;
154 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size())
155 int numFields(gatherFields_.size());
156 fieldIds_.resize(numFields);
157 for (int fd(0); fd < numFields; ++fd)
158 {
159 // Get the field ID from the DOF manager.
160 const string& fieldName(indexerNames_[fd]);
161 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
162
163 // This is the error return code; raise the alarm.
164 TEUCHOS_TEST_FOR_EXCEPTION(fieldIds_[fd] == -1, logic_error,
165 "GatherSolution_Epetra<Residual>: Could not find field \"" +
166 fieldName + "\" in the global indexer. ")
167 } // end loop over gatherFields_
168 indexerNames_.clear();
169} // end of postRegistrationSetup() (Residual Specialization)
170
172//
173// preEvaluate() (Residual Specialization)
174//
176template<typename TRAITS, typename LO, typename GO>
177void
180 typename TRAITS::PreEvalData d)
181{
182 using std::string;
183 using Teuchos::RCP;
184 using Teuchos::rcp_dynamic_cast;
188 using LOC = panzer::LinearObjContainer;
190
191 // First try the refactored ReadOnly container.
192 RCP<GED> ged;
193 string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
194 if (d.gedc->containsDataObject(globalDataKey_ + post))
195 {
196 ged = d.gedc->getDataObject(globalDataKey_ + post);
197 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
198 return;
199 } // end of the refactored ReadOnly way
200
201 // Now try the old path.
202 ged = d.gedc->getDataObject(globalDataKey_);
203 {
204 // Try to extract the linear object container.
205 auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
206 auto locPair = rcp_dynamic_cast<LPGED>(ged);
207 if (not locPair.is_null())
208 {
209 RCP<LOC> loc = locPair->getGhostedLOC();
210 epetraContainer = rcp_dynamic_cast<ELOC>(loc);
211 } // end if (not locPair.is_null())
212 if (not epetraContainer.is_null())
213 {
214 if (useTimeDerivativeSolutionVector_)
215 x_ = epetraContainer->get_dxdt();
216 else // if (not useTimeDerivativeSolutionVector_)
217 x_ = epetraContainer->get_x();
218 return;
219 } // end if (not epetraContainer.is_null())
220 } // end of the old path
221
222 // As a last resort, try to extract an EpetraVector_ReadOnly object. This
223 // will throw an exception if it doesn't work.
224 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
225} // end of preEvaluate() (Residual Specialization)
226
228//
229// evaluateFields() (Residual Specialization)
230//
232template<typename TRAITS, typename LO, typename GO>
233void
236 typename TRAITS::EvalData workset)
237{
238 using PHX::MDField;
239 using std::size_t;
240 using std::string;
241 using std::vector;
242 using Teuchos::ArrayRCP;
243 using Teuchos::ptrFromRef;
244 using Teuchos::RCP;
245 using Teuchos::rcp_dynamic_cast;
246 using Thyra::SpmdVectorBase;
247
248 // For convenience, pull out some objects from the workset.
249 string blockId(this->wda(workset).block_id);
250 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
251 int numCells(localCellIds.size()), numFields(gatherFields_.size());
252
253 // NOTE: A reordering of these loops will likely improve performance. The
254 // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
255 // can be cheaper. However the lookup for LIDs may be more expensive!
256
257 auto LIDs = globalIndexer_->getLIDs();
258 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
259 Kokkos::deep_copy(LIDs_h, LIDs);
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 auto field_h = Kokkos::create_mirror_view(field.get_static_view());
265 int fieldNum(fieldIds_[fieldInd]);
266 const vector<int>& elmtOffset =
267 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
268 int numBases(elmtOffset.size());
269 // Gather operation for each cell in the workset.
270 for (int cell(0); cell < numCells; ++cell)
271 {
272 size_t cellLocalId(localCellIds[cell]);
273 // Loop over the basis functions and fill the fields.
274 for (int basis(0); basis < numBases; ++basis)
275 {
276 int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
277 if (x_.is_null())
278 field_h(cell, basis) = (*xEvRoGed_)[lid];
279 else
280 field_h(cell, basis) = (*x_)[lid];
281 } // end loop over the basis functions
282 } // end loop over localCellIds
283 Kokkos::deep_copy(field.get_static_view(), field_h);
284 } // end loop over the fields to be gathered
285
286} // end of evaluateFields() (Residual Specialization)
287
289//
290// Initializing Constructor (Tangent Specialization)
291//
293template<typename TRAITS, typename LO, typename GO>
296 const Teuchos::RCP<const panzer::GlobalIndexer>& indexer,
297 const Teuchos::ParameterList& p)
298 :
299 globalIndexer_(indexer),
300 hasTangentFields_(false)
301{
302 using panzer::PureBasis;
303 using PHX::MDField;
304 using PHX::print;
305 using std::size_t;
306 using std::string;
307 using std::vector;
308 using Teuchos::RCP;
309 using vvstring = std::vector<std::vector<std::string>>;
311 input.setParameterList(p);
312 const vector<string>& names = input.getDofNames();
313 RCP<const PureBasis> basis = input.getBasis();
314 const vvstring& tangentFieldNames = input.getTangentNames();
315 indexerNames_ = input.getIndexerNames();
316 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
317 globalDataKey_ = input.getGlobalDataKey();
318
319 // Allocate the fields.
320 int numFields(names.size());
321 gatherFields_.resize(numFields);
322 for (int fd(0); fd < numFields; ++fd)
323 {
324 gatherFields_[fd] =
325 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
326 this->addEvaluatedField(gatherFields_[fd]);
327 } // end loop over names
328
329 // Set up the dependent tangent fields, if requested.
330 if (tangentFieldNames.size() > 0)
331 {
332 TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size())
333 hasTangentFields_ = true;
334 tangentFields_.resize(numFields);
335 for (int fd(0); fd < numFields; ++fd)
336 {
337 int numTangentFields(tangentFieldNames[fd].size());
338 tangentFields_[fd].resize(numTangentFields);
339 for (int i(0); i < numTangentFields; ++i)
340 {
341 tangentFields_[fd][i] =
342 MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
343 basis->functional);
344 this->addDependentField(tangentFields_[fd][i]);
345 } // end loop over tangeng_field_names[fd]
346 } // end loop over gatherFields_
347 } // end if (tangentFieldNames.size() > 0)
348
349 // Figure out what the first active name is.
350 string firstName("<none>");
351 if (numFields > 0)
352 firstName = names[0];
353 string n("GatherSolution (Epetra): " + firstName + " (" +
354 print<EvalT>() + ")");
355 this->setName(n);
356} // end of Initializing Constructor (Tangent Specialization)
357
359//
360// postRegistrationSetup() (Tangent Specialization)
361//
363template<typename TRAITS, typename LO, typename GO>
364void
367 typename TRAITS::SetupData /* d */,
369{
370 using std::logic_error;
371 using std::size_t;
372 using std::string;
373 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size())
374 int numFields(gatherFields_.size());
375 fieldIds_.resize(numFields);
376 for (int fd(0); fd < numFields; ++fd)
377 {
378 // Get the field ID from the DOF manager.
379 const string& fieldName(indexerNames_[fd]);
380 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
381
382 // This is the error return code; raise the alarm.
383 TEUCHOS_TEST_FOR_EXCEPTION(fieldIds_[fd] == -1, logic_error,
384 "GatherSolution_Epetra<Tangent>: Could not find field \"" + fieldName
385 + "\" in the global indexer. ")
386 } // end loop over gatherFields_
387 indexerNames_.clear();
388} // end of postRegistrationSetup() (Tangent Specialization)
389
391//
392// preEvaluate() (Tangent Specialization)
393//
395template<typename TRAITS, typename LO, typename GO>
396void
399 typename TRAITS::PreEvalData d)
400{
401 using std::string;
402 using Teuchos::RCP;
403 using Teuchos::rcp_dynamic_cast;
407 using LOC = panzer::LinearObjContainer;
409
410 // First try the refactored ReadOnly container.
411 RCP<GED> ged;
412 string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
413 if (d.gedc->containsDataObject(globalDataKey_ + post))
414 {
415 ged = d.gedc->getDataObject(globalDataKey_ + post);
416 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
417 return;
418 } // end of the refactored ReadOnly way
419
420 // Now try the old path.
421 ged = d.gedc->getDataObject(globalDataKey_);
422 {
423 // Try to extract the linear object container.
424 auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
425 auto locPair = rcp_dynamic_cast<LPGED>(ged);
426 if (not locPair.is_null())
427 {
428 RCP<LOC> loc = locPair->getGhostedLOC();
429 epetraContainer = rcp_dynamic_cast<ELOC>(loc);
430 } // end if (not locPair.is_null())
431 if (not epetraContainer.is_null())
432 {
433 if (useTimeDerivativeSolutionVector_)
434 x_ = epetraContainer->get_dxdt();
435 else // if (not useTimeDerivativeSolutionVector_)
436 x_ = epetraContainer->get_x();
437 return;
438 } // end if (not epetraContainer.is_null())
439 } // end of the old path
440
441 // As a last resort, try to extract an EpetraVector_ReadOnly object. This
442 // will throw an exception if it doesn't work.
443 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
444} // end of preEvaluate() (Tangent Specialization)
445
447//
448// evaluateFields() (Tangent Specialization)
449//
451template<typename TRAITS, typename LO, typename GO>
452void
455 typename TRAITS::EvalData workset)
456{
457 using PHX::MDField;
458 using std::size_t;
459 using std::string;
460 using std::vector;
461 using Teuchos::ArrayRCP;
462 using Teuchos::ptrFromRef;
463 using Teuchos::RCP;
464 using Teuchos::rcp_dynamic_cast;
465 using Thyra::SpmdVectorBase;
466
467 // For convenience, pull out some objects from the workset.
468 string blockId(this->wda(workset).block_id);
469 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
470 int numCells(localCellIds.size()), numFields(gatherFields_.size());
471
472 // NOTE: A reordering of these loops will likely improve performance. The
473 // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
474 // can be cheaper. However the lookup for LIDs may be more expensive!
475 auto LIDs = globalIndexer_->getLIDs();
476 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
477 Kokkos::deep_copy(LIDs_h, LIDs);
478 // Loop over the fields to be gathered.
479 for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
480 {
481 MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
482 auto field_h = Kokkos::create_mirror_view(field.get_static_view());
483 int fieldNum(fieldIds_[fieldInd]);
484 const vector<int>& elmtOffset =
485 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
486 int numBases(elmtOffset.size());
487 // Gather operation for each cell in the workset.
488 for (int cell(0); cell < numCells; ++cell)
489 {
490 size_t cellLocalId(localCellIds[cell]);
491 // Loop over the basis functions and fill the fields.
492 for (int basis(0); basis < numBases; ++basis)
493 {
494 int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
495 if (x_.is_null())
496 field_h(cell, basis) = (*xEvRoGed_)[lid];
497 else
498 field_h(cell, basis) = (*x_)[lid];
499 } // end loop over the basis functions
500 } // end loop over localCellIds
501
502 // Deal with the tangent fields.
503 if (hasTangentFields_) {
504 // Loop over the tangent fields and fill them in.
505 int numTangentFields(tangentFields_[fieldInd].size());
506 for (int i(0); i < numTangentFields; ++i){
507 auto tf = Kokkos::create_mirror_view(tangentFields_[fieldInd][i].get_static_view());
508 Kokkos::deep_copy(tf, tangentFields_[fieldInd][i].get_static_view());
509 // Loop over the cells in the workset.
510 for (int cell(0); cell < numCells; ++cell) {
511 // Loop over the fields to be gathered.
512 int fieldNum(fieldIds_[fieldInd]);
513 const vector<int>& elmtOffset =
514 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
515 int numBases(elmtOffset.size());
516
517 // Loop over the basis functions.
518 for (int basis(0); basis < numBases; ++basis){
519 field_h(cell, basis).fastAccessDx(i) =
520 tf(cell, basis).val();
521 } // end loop over the basis functions
522 } // end loop over the cells in the workset
523 } // end loop over numTangentFields
524 } // end if (hasTangentFields_)
525 Kokkos::deep_copy(field.get_static_view(), field_h);
526 } // end loop over the fields to be gathered
527} // end of evaluateFields() (Tangent Specialization)
528
530//
531// Initializing Constructor (Jacobian Specialization)
532//
534template<typename TRAITS, typename LO, typename GO>
537 const Teuchos::RCP<const panzer::GlobalIndexer>& indexer,
538 const Teuchos::ParameterList& p)
539 :
540 globalIndexer_(indexer)
541{
542 using panzer::PureBasis;
543 using PHX::MDField;
544 using PHX::print;
545 using std::size_t;
546 using std::string;
547 using std::vector;
548 using Teuchos::RCP;
550 input.setParameterList(p);
551 const vector<string>& names = input.getDofNames();
552 RCP<const PureBasis> basis = input.getBasis();
553 indexerNames_ = input.getIndexerNames();
554 useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
555 globalDataKey_ = input.getGlobalDataKey();
556 gatherSeedIndex_ = input.getGatherSeedIndex();
557 sensitivitiesName_ = input.getSensitivitiesName();
558 disableSensitivities_ = not input.firstSensitivitiesAvailable();
559
560 // Allocate the fields.
561 int numFields(names.size());
562 gatherFields_.resize(numFields);
563 for (int fd(0); fd < numFields; ++fd)
564 {
565 MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
566 gatherFields_[fd] = f;
567 this->addEvaluatedField(gatherFields_[fd]);
568 } // end loop over names
569
570 // Figure out what the first active name is.
571 string firstName("<none>"), n("GatherSolution (Epetra");
572 if (numFields > 0)
573 firstName = names[0];
574 if (disableSensitivities_)
575 n += ", No Sensitivities";
576 n += "): " + firstName + " (Jacobian)";
577 this->setName(n);
578} // end of Initializing Constructor (Jacobian Specialization)
579
581//
582// postRegistrationSetup() (Jacobian Specialization)
583//
585template<typename TRAITS, typename LO, typename GO>
586void
589 typename TRAITS::SetupData /* d */,
591{
592 using std::logic_error;
593 using std::size_t;
594 using std::string;
595 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size())
596 int numFields(gatherFields_.size());
597 fieldIds_.resize(numFields);
598 for (int fd(0); fd < numFields; ++fd)
599 {
600 // Get the field ID from the DOF manager.
601 const string& fieldName(indexerNames_[fd]);
602 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
603
604 // This is the error return code; raise the alarm.
605 TEUCHOS_TEST_FOR_EXCEPTION(fieldIds_[fd] == -1, logic_error,
606 "GatherSolution_Epetra<Jacobian>: Could not find field \"" +
607 fieldName + "\" in the global indexer. ")
608 } // end loop over gatherFields_
609 indexerNames_.clear();
610} // end of postRegistrationSetup() (Jacobian Specialization)
611
613//
614// preEvaluate() (Jacobian Specialization)
615//
617template<typename TRAITS, typename LO, typename GO>
618void
621 typename TRAITS::PreEvalData d)
622{
623 using std::string;
624 using Teuchos::RCP;
625 using Teuchos::rcp_dynamic_cast;
629 using LOC = panzer::LinearObjContainer;
631
632 // Manage sensitivities.
633 applySensitivities_ = false;
634 if ((not disableSensitivities_ ) and
635 (d.first_sensitivities_name == sensitivitiesName_))
636 applySensitivities_ = true;
637
638 // First try the refactored ReadOnly container.
639 RCP<GED> ged;
640 string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
641 if (d.gedc->containsDataObject(globalDataKey_ + post))
642 {
643 ged = d.gedc->getDataObject(globalDataKey_ + post);
644 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
645 return;
646 } // end of the refactored ReadOnly way
647
648 // Now try the old path.
649 ged = d.gedc->getDataObject(globalDataKey_);
650 {
651 // Try to extract the linear object container.
652 auto epetraContainer = rcp_dynamic_cast<ELOC>(ged);
653 auto locPair = rcp_dynamic_cast<LPGED>(ged);
654 if (not locPair.is_null())
655 {
656 RCP<LOC> loc = locPair->getGhostedLOC();
657 epetraContainer = rcp_dynamic_cast<ELOC>(loc);
658 } // end if (not locPair.is_null())
659 if (not epetraContainer.is_null())
660 {
661 if (useTimeDerivativeSolutionVector_)
662 x_ = epetraContainer->get_dxdt();
663 else // if (not useTimeDerivativeSolutionVector_)
664 x_ = epetraContainer->get_x();
665 return;
666 } // end if (not epetraContainer.is_null())
667 } // end of the old path
668
669 // As a last resort, try to extract an EpetraVector_ReadOnly object. This
670 // will throw an exception if it doesn't work.
671 xEvRoGed_ = rcp_dynamic_cast<EVROGED>(ged, true);
672} // end of preEvaluate() (Jacobian Specialization)
673
675//
676// evaluateFields() (Jacobian Specialization)
677//
679template<typename TRAITS, typename LO, typename GO>
680void
683 typename TRAITS::EvalData workset)
684{
685 using PHX::MDField;
686 using std::size_t;
687 using std::string;
688 using std::vector;
689 using Teuchos::ArrayRCP;
690 using Teuchos::ptrFromRef;
691 using Teuchos::RCP;
692 using Teuchos::rcp_dynamic_cast;
693 using Thyra::SpmdVectorBase;
694
695 // For convenience, pull out some objects from the workset.
696 string blockId(this->wda(workset).block_id);
697 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
698 int numFields(gatherFields_.size()), numCells(localCellIds.size());
699
700 // Set a sensitivity seed value.
701 double seedValue(0);
702 if (applySensitivities_)
703 {
704 if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
705 seedValue = workset.alpha;
706 else if (gatherSeedIndex_ < 0)
707 seedValue = workset.beta;
708 else if (not useTimeDerivativeSolutionVector_)
709 seedValue = workset.gather_seeds[gatherSeedIndex_];
710 else
711 TEUCHOS_ASSERT(false)
712 } // end if (applySensitivities_)
713
714 // Turn off sensitivies. This may be faster if we don't expand the term, but
715 // I suspect not, because anywhere it is used the full complement of
716 // sensitivies will be needed anyway.
717 if (not applySensitivities_)
718 seedValue = 0.0;
719
720 // Interface worksets handle DOFs from two element blocks. The derivative
721 // offset for the other element block must be shifted by the derivative side
722 // of my element block.
723 int dos(0);
724 if (this->wda.getDetailsIndex() == 1)
725 {
726 // Get the DOF count for my element block.
727 dos = globalIndexer_->getElementBlockGIDCount(workset.details(0).block_id);
728 } // end if (this->wda.getDetailsIndex() == 1)
729
730 // NOTE: A reordering of these loops will likely improve performance. The
731 // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
732 // can be cheaper. However the lookup for LIDs may be more expensive!
733 auto LIDs = globalIndexer_->getLIDs();
734 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
735 Kokkos::deep_copy(LIDs_h, LIDs);
736 // Loop over the fields to be gathered.
737 for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
738 {
739 MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
740 auto field_h = Kokkos::create_mirror_view(field.get_static_view());
741 int fieldNum(fieldIds_[fieldInd]);
742 const vector<int>& elmtOffset =
743 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
744 int numBases(elmtOffset.size());
745 // Gather operation for each cell in the workset.
746 for (int cell(0); cell < numCells; ++cell)
747 {
748 size_t cellLocalId(localCellIds[cell]);
749 // Loop over the basis functions and fill the fields.
750 for (int basis(0); basis < numBases; ++basis)
751 {
752 int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
753 if (x_.is_null())
754 field_h(cell, basis) = (*xEvRoGed_)[lid];
755 else
756 field_h(cell, basis) = (*x_)[lid];
757 } // end loop over the basis functions
758 } // end loop over localCellIds
759
760 // Deal with the sensitivities.
761 if (applySensitivities_) {
762 int fieldNum(fieldIds_[fieldInd]);
763 const vector<int>& elmtOffset =
764 globalIndexer_->getGIDFieldOffsets(blockId, fieldNum);
765 int numBases(elmtOffset.size());
766
767 // Gather operation for each cell in the workset.
768 for (int cell(0); cell < numCells; ++cell)
769 {
770 // Loop over the basis functions.
771 for (int basis(0); basis < numBases; ++basis)
772 {
773 // Seed the FAD object.
774 int offset(elmtOffset[basis]);
775
776 field_h(cell, basis).fastAccessDx(dos + offset) = seedValue;
777 } // end loop over the basis functions
778 } // end loop over localCellIds
779 } // end if (applySensitivities_)
780 Kokkos::deep_copy(field.get_static_view(), field_h);
781 } // end loop over the fields to be gathered
782
783} // end of evaluateFields() (Jacobian Specialization)
784
785#endif // __Panzer_GatherSolution_Epetra_impl_hpp__
int numFields
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
This class provides a boundary exchange communication mechanism for vectors.
Gathers solution values from the Newton solution vector into the nodal fields of the field manager.
std::string getGlobalDataKey() const
Name of the global evaluation data container to use for the source vector (all types)
void setParameterList(const Teuchos::ParameterList &pl)
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
const std::vector< std::vector< std::string > > & getTangentNames() const
Get the name of the tangent fields (tangent only)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
const std::vector< std::string > & getIndexerNames() const
std::string getSensitivitiesName() const
The name of the sensitivities. Enables sensitivities at "preEvaluate" time (Jacobian and Hessian)
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
Description and data layouts associated with a particular basis.