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