Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_ScatterDirichletResidual_Tpetra_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_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP
44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP
45 
46 #include "Teuchos_RCP.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 #include "Phalanx_DataLayout.hpp"
50 
51 #include "Panzer_GlobalIndexer.hpp"
52 #include "Panzer_PureBasis.hpp"
57 
58 #include "Phalanx_DataLayout_MDALayout.hpp"
59 
60 #include "Teuchos_FancyOStream.hpp"
61 
62 // **********************************************************************
63 // Specialization: Residual
64 // **********************************************************************
65 
66 
67 template<typename TRAITS,typename LO,typename GO,typename NodeT>
70  const Teuchos::ParameterList& p)
71  : globalIndexer_(indexer)
72  , globalDataKey_("Residual Scatter Container")
73 {
74  std::string scatterName = p.get<std::string>("Scatter Name");
75  scatterHolder_ =
77 
78  // get names to be evaluated
79  const std::vector<std::string>& names =
80  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
81 
82  // grab map from evaluated names to field names
83  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
84 
85  // determine if we are scattering an initial condition
86  scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
87 
88  Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
89  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
90  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional ;
91  if (!scatterIC_) {
92  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
93  local_side_id_ = p.get<int>("Local Side ID");
94  }
95 
96  // build the vector of fields that this is dependent on
97  scatterFields_.resize(names.size());
98  for (std::size_t eq = 0; eq < names.size(); ++eq) {
99  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
100 
101  // tell the field manager that we depend on this field
102  this->addDependentField(scatterFields_[eq]);
103  }
104 
105  checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
106  if (checkApplyBC_) {
107  applyBC_.resize(names.size());
108  for (std::size_t eq = 0; eq < names.size(); ++eq) {
109  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
110  this->addDependentField(applyBC_[eq]);
111  }
112  }
113 
114  // this is what this evaluator provides
115  this->addEvaluatedField(*scatterHolder_);
116 
117  if (p.isType<std::string>("Global Data Key"))
118  globalDataKey_ = p.get<std::string>("Global Data Key");
119 
120  this->setName(scatterName+" Scatter Residual");
121 }
122 
123 // **********************************************************************
124 template<typename TRAITS,typename LO,typename GO,typename NodeT>
126 postRegistrationSetup(typename TRAITS::SetupData /* d */,
127  PHX::FieldManager<TRAITS>& /* fm */)
128 {
129  fieldIds_.resize(scatterFields_.size());
130 
131  // load required field numbers for fast use
132  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
133  // get field ID from DOF manager
134  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
135  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
136  }
137 
138  // get the number of nodes (Should be renamed basis)
139  num_nodes = scatterFields_[0].extent(1);
140 }
141 
142 // **********************************************************************
143 template<typename TRAITS,typename LO,typename GO,typename NodeT>
145 preEvaluate(typename TRAITS::PreEvalData d)
146 {
147  // extract linear object container
148  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
149 
150  if(tpetraContainer_==Teuchos::null) {
151  // extract linear object container
152  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
153  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
154 
155  dirichletCounter_ = Teuchos::null;
156  }
157  else {
158  // extract dirichlet counter from container
159  Teuchos::RCP<LOC> tpetraContainer
160  = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
161 
162  dirichletCounter_ = tpetraContainer->get_f();
163  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
164  }
165 }
166 
167 // **********************************************************************
168 template<typename TRAITS,typename LO,typename GO,typename NodeT>
170 evaluateFields(typename TRAITS::EvalData workset)
171 {
172  std::vector<GO> GIDs;
173  std::vector<LO> LIDs;
174 
175  // for convenience pull out some objects from workset
176  std::string blockId = this->wda(workset).block_id;
177  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
178 
179  Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
180  tpetraContainer_->get_f() :
181  tpetraContainer_->get_x();
182 
183  Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
184  Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
185 
186  // NOTE: A reordering of these loops will likely improve performance
187  // The "getGIDFieldOffsets may be expensive. However the
188  // "getElementGIDs" can be cheaper. However the lookup for LIDs
189  // may be more expensive!
190 
191 
192 
193  // loop over each field to be scattered
194  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
195  int fieldNum = fieldIds_[fieldIndex];
196  auto scatterFields_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
197  Kokkos::deep_copy(scatterFields_h, scatterFields_[fieldIndex].get_static_view());
198 
199  // scatter operation for each cell in workset
200  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
201  std::size_t cellLocalId = localCellIds[worksetCellIndex];
202 
203  globalIndexer_->getElementGIDs(cellLocalId,GIDs);
204 
205  // caculate the local IDs for this element
206  LIDs.resize(GIDs.size());
207  for(std::size_t i=0;i<GIDs.size();i++)
208  LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
209 
210  if (!scatterIC_) {
211  // this call "should" get the right ordering according to the Intrepid2 basis
212  const std::pair<std::vector<int>,std::vector<int> > & indicePair
213  = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
214  const std::vector<int> & elmtOffset = indicePair.first;
215  const std::vector<int> & basisIdMap = indicePair.second;
216 
217  // loop over basis functions
218  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
219  int offset = elmtOffset[basis];
220  LO lid = LIDs[offset];
221  if(lid<0) // not on this processor!
222  continue;
223 
224  int basisId = basisIdMap[basis];
225 
226  if (checkApplyBC_)
227  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
228  continue;
229 
230  r_array[lid] = scatterFields_h(worksetCellIndex,basisId);
231 
232  // record that you set a dirichlet condition
233  dc_array[lid] = 1.0;
234  }
235  } else {
236  // this call "should" get the right ordering according to the Intrepid2 basis
237  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
238 
239  // loop over basis functions
240  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
241  int offset = elmtOffset[basis];
242  LO lid = LIDs[offset];
243  if(lid<0) // not on this processor!
244  continue;
245 
246  r_array[lid] = scatterFields_h(worksetCellIndex,basis);
247 
248  // record that you set a dirichlet condition
249  dc_array[lid] = 1.0;
250  }
251  }
252  }
253  }
254 }
255 
256 // **********************************************************************
257 // Specialization: Tangent
258 // **********************************************************************
259 
260 
261 template<typename TRAITS,typename LO,typename GO,typename NodeT>
264  const Teuchos::ParameterList& p)
265  : globalIndexer_(indexer)
266  , globalDataKey_("Residual Scatter Container")
267 {
268  std::string scatterName = p.get<std::string>("Scatter Name");
269  scatterHolder_ =
271 
272  // get names to be evaluated
273  const std::vector<std::string>& names =
274  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
275 
276  // grab map from evaluated names to field names
277  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
278 
279  // determine if we are scattering an initial condition
280  scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
281 
282  Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
283  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
284  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional ;
285  if (!scatterIC_) {
286  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
287  local_side_id_ = p.get<int>("Local Side ID");
288  }
289 
290  // build the vector of fields that this is dependent on
291  scatterFields_.resize(names.size());
292  for (std::size_t eq = 0; eq < names.size(); ++eq) {
293  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
294 
295  // tell the field manager that we depend on this field
296  this->addDependentField(scatterFields_[eq]);
297  }
298 
299  checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
300  if (checkApplyBC_) {
301  applyBC_.resize(names.size());
302  for (std::size_t eq = 0; eq < names.size(); ++eq) {
303  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
304  this->addDependentField(applyBC_[eq]);
305  }
306  }
307 
308  // this is what this evaluator provides
309  this->addEvaluatedField(*scatterHolder_);
310 
311  if (p.isType<std::string>("Global Data Key"))
312  globalDataKey_ = p.get<std::string>("Global Data Key");
313 
314  this->setName(scatterName+" Scatter Tangent");
315 }
316 
317 // **********************************************************************
318 template<typename TRAITS,typename LO,typename GO,typename NodeT>
320 postRegistrationSetup(typename TRAITS::SetupData /* d */,
321  PHX::FieldManager<TRAITS>& /* fm */)
322 {
323  fieldIds_.resize(scatterFields_.size());
324 
325  // load required field numbers for fast use
326  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
327  // get field ID from DOF manager
328  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
329  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
330  }
331 
332  // get the number of nodes (Should be renamed basis)
333  num_nodes = scatterFields_[0].extent(1);
334 }
335 
336 // **********************************************************************
337 template<typename TRAITS,typename LO,typename GO,typename NodeT>
339 preEvaluate(typename TRAITS::PreEvalData d)
340 {
341  // extract linear object container
342  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
343 
344  if(tpetraContainer_==Teuchos::null) {
345  // extract linear object container
346  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
347  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
348 
349  dirichletCounter_ = Teuchos::null;
350  }
351  else {
352  // extract dirichlet counter from container
353  Teuchos::RCP<LOC> tpetraContainer
354  = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
355 
356  dirichletCounter_ = tpetraContainer->get_f();
357  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
358  }
359 
360  using Teuchos::RCP;
361  using Teuchos::rcp_dynamic_cast;
362 
363  // this is the list of parameters and their names that this scatter has to account for
364  std::vector<std::string> activeParameters =
365  rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
366 
367  // ETP 02/03/16: This code needs to be updated to properly handle scatterIC_
368  TEUCHOS_ASSERT(!scatterIC_);
369  dfdp_vectors_.clear();
370  for(std::size_t i=0;i<activeParameters.size();i++) {
372  rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
373  Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
374  dfdp_vectors_.push_back(vec_array);
375  }
376 }
377 
378 // **********************************************************************
379 template<typename TRAITS,typename LO,typename GO,typename NodeT>
381 evaluateFields(typename TRAITS::EvalData workset)
382 {
383  std::vector<GO> GIDs;
384  std::vector<LO> LIDs;
385 
386  // for convenience pull out some objects from workset
387  std::string blockId = this->wda(workset).block_id;
388  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
389 
390  Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
391  tpetraContainer_->get_f() :
392  tpetraContainer_->get_x();
393 
394  Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
395  Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
396 
397  // NOTE: A reordering of these loops will likely improve performance
398  // The "getGIDFieldOffsets may be expensive. However the
399  // "getElementGIDs" can be cheaper. However the lookup for LIDs
400  // may be more expensive!
401 
402 
403  // scatter operation for each cell in workset
404  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
405  std::size_t cellLocalId = localCellIds[worksetCellIndex];
406 
407  globalIndexer_->getElementGIDs(cellLocalId,GIDs);
408 
409  // caculate the local IDs for this element
410  LIDs.resize(GIDs.size());
411  for(std::size_t i=0;i<GIDs.size();i++)
412  LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
413 
414  // loop over each field to be scattered
415  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
416  int fieldNum = fieldIds_[fieldIndex];
417 
418  if (!scatterIC_) {
419  // this call "should" get the right ordering according to the Intrepid2 basis
420  const std::pair<std::vector<int>,std::vector<int> > & indicePair
421  = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
422  const std::vector<int> & elmtOffset = indicePair.first;
423  const std::vector<int> & basisIdMap = indicePair.second;
424 
425  // loop over basis functions
426  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
427  int offset = elmtOffset[basis];
428  LO lid = LIDs[offset];
429  if(lid<0) // not on this processor!
430  continue;
431 
432  int basisId = basisIdMap[basis];
433 
434  if (checkApplyBC_)
435  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
436  continue;
437 
438  ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
439  //r_array[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
440 
441  // then scatter the sensitivity vectors
442  if(value.size()==0)
443  for(std::size_t d=0;d<dfdp_vectors_.size();d++)
444  dfdp_vectors_[d][lid] = 0.0;
445  else
446  for(int d=0;d<value.size();d++) {
447  dfdp_vectors_[d][lid] = value.fastAccessDx(d);
448  }
449 
450  // record that you set a dirichlet condition
451  dc_array[lid] = 1.0;
452  }
453  } else {
454  // this call "should" get the right ordering according to the Intrepid2 basis
455  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
456 
457  // loop over basis functions
458  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
459  int offset = elmtOffset[basis];
460  LO lid = LIDs[offset];
461  if(lid<0) // not on this processor!
462  continue;
463 
464  ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
465  //r_array[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
466 
467  // then scatter the sensitivity vectors
468  if(value.size()==0)
469  for(std::size_t d=0;d<dfdp_vectors_.size();d++)
470  dfdp_vectors_[d][lid] = 0.0;
471  else
472  for(int d=0;d<value.size();d++) {
473  dfdp_vectors_[d][lid] = value.fastAccessDx(d);
474  }
475 
476  // record that you set a dirichlet condition
477  dc_array[lid] = 1.0;
478  }
479  }
480  }
481  }
482 }
483 
484 // **********************************************************************
485 // Specialization: Jacobian
486 // **********************************************************************
487 
488 template<typename TRAITS,typename LO,typename GO,typename NodeT>
491  const Teuchos::ParameterList& p)
492  : globalIndexer_(indexer)
493  , globalDataKey_("Residual Scatter Container")
494 {
495  std::string scatterName = p.get<std::string>("Scatter Name");
496  scatterHolder_ =
498 
499  // get names to be evaluated
500  const std::vector<std::string>& names =
501  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
502 
503  // grab map from evaluated names to field names
504  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
505 
507  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
508 
509  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
510  local_side_id_ = p.get<int>("Local Side ID");
511 
512  // build the vector of fields that this is dependent on
513  scatterFields_.resize(names.size());
514  for (std::size_t eq = 0; eq < names.size(); ++eq) {
515  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
516 
517  // tell the field manager that we depend on this field
518  this->addDependentField(scatterFields_[eq]);
519  }
520 
521  checkApplyBC_ = p.get<bool>("Check Apply BC");
522  if (checkApplyBC_) {
523  applyBC_.resize(names.size());
524  for (std::size_t eq = 0; eq < names.size(); ++eq) {
525  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
526  this->addDependentField(applyBC_[eq]);
527  }
528  }
529 
530  // this is what this evaluator provides
531  this->addEvaluatedField(*scatterHolder_);
532 
533  if (p.isType<std::string>("Global Data Key"))
534  globalDataKey_ = p.get<std::string>("Global Data Key");
535 
536  this->setName(scatterName+" Scatter Residual (Jacobian)");
537 }
538 
539 // **********************************************************************
540 template<typename TRAITS,typename LO,typename GO,typename NodeT>
542 postRegistrationSetup(typename TRAITS::SetupData /* d */,
543  PHX::FieldManager<TRAITS>& /* fm */)
544 {
545  fieldIds_.resize(scatterFields_.size());
546  // load required field numbers for fast use
547  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
548  // get field ID from DOF manager
549  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
550  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
551  }
552 
553  // get the number of nodes (Should be renamed basis)
554  num_nodes = scatterFields_[0].extent(1);
555  num_eq = scatterFields_.size();
556 }
557 
558 // **********************************************************************
559 template<typename TRAITS,typename LO,typename GO,typename NodeT>
561 preEvaluate(typename TRAITS::PreEvalData d)
562 {
563  // extract linear object container
564  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
565 
566  if(tpetraContainer_==Teuchos::null) {
567  // extract linear object container
568  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
569  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
570 
571  dirichletCounter_ = Teuchos::null;
572  }
573  else {
574  // extract dirichlet counter from container
575  Teuchos::RCP<LOC> tpetraContainer
576  = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
577 
578  dirichletCounter_ = tpetraContainer->get_f();
579  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
580  }
581 }
582 
583 // **********************************************************************
584 template<typename TRAITS,typename LO,typename GO,typename NodeT>
586 evaluateFields(typename TRAITS::EvalData workset)
587 {
588  std::vector<GO> GIDs;
589 
590  // for convenience pull out some objects from workset
591  std::string blockId = this->wda(workset).block_id;
592  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
593 
594  Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
595  Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
596 
597  Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
598  Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
599 
600  // NOTE: A reordering of these loops will likely improve performance
601  // The "getGIDFieldOffsets may be expensive. However the
602  // "getElementGIDs" can be cheaper. However the lookup for LIDs
603  // may be more expensive!
604 
605  // scatter operation for each cell in workset
606  auto LIDs = globalIndexer_->getLIDs();
607  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
608  Kokkos::deep_copy(LIDs_h, LIDs);
609  // loop over each field to be scattered
610  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
611  int fieldNum = fieldIds_[fieldIndex];
612  auto scatterFields_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
613  Kokkos::deep_copy(scatterFields_h, scatterFields_[fieldIndex].get_static_view());
614  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
615  std::size_t cellLocalId = localCellIds[worksetCellIndex];
616 
617  globalIndexer_->getElementGIDs(cellLocalId,GIDs);
618 
619  // this call "should" get the right ordering according to the Intrepid2 basis
620  const std::pair<std::vector<int>,std::vector<int> > & indicePair
621  = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
622  const std::vector<int> & elmtOffset = indicePair.first;
623  const std::vector<int> & basisIdMap = indicePair.second;
624 
625  // loop over basis functions
626  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
627  int offset = elmtOffset[basis];
628  int lid = LIDs_h(cellLocalId, offset);
629  if(lid<0) // not on this processor
630  continue;
631 
632  int basisId = basisIdMap[basis];
633 
634  if (checkApplyBC_)
635  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
636  continue;
637 
638  // zero out matrix row
639  {
640  std::size_t sz = Jac->getNumEntriesInLocalRow(lid);
641  std::size_t numEntries = 0;
642  typename LOC::CrsMatrixType::nonconst_local_inds_host_view_type rowIndices("indices", sz);
643  typename LOC::CrsMatrixType::nonconst_values_host_view_type rowValues("values", sz);
644 
645  // Jac->getLocalRowView(lid,numEntries,rowValues,rowIndices);
646  Jac->getLocalRowCopy(lid,rowIndices,rowValues,numEntries);
647 
648  for(std::size_t i=0;i<numEntries;i++)
649  rowValues(i) = 0.0;
650 
651  Jac->replaceLocalValues(lid,rowIndices,rowValues);
652  }
653 
654  GO gid = GIDs[offset];
655  const ScalarT scatterField = scatterFields_h(worksetCellIndex,basisId);
656 
657  r_array[lid] = scatterField.val();
658  dc_array[lid] = 1.0; // mark row as dirichlet
659 
660  // loop over the sensitivity indices: all DOFs on a cell
661  std::vector<double> jacRow(scatterField.size(),0.0);
662 
663  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
664  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
665  TEUCHOS_ASSERT(jacRow.size()==GIDs.size());
666 
667  Jac->replaceGlobalValues(gid, GIDs, jacRow);
668  }
669  }
670  }
671 }
672 
673 // **********************************************************************
674 
675 #endif
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void resize(const size_type n, const T &val=T())
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pushes residual values into the residual vector for a Newton-based solve.
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)