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