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