Panzer  Version of the Day
Panzer_ScatterDirichletResidual_BlockedEpetra_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_BLOCEDEPETRA_IMPL_HPP
44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_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 
58 #include "Panzer_PureBasis.hpp"
60 
61 #include "Phalanx_DataLayout_MDALayout.hpp"
62 
63 #include "Thyra_SpmdVectorBase.hpp"
64 #include "Thyra_ProductVectorBase.hpp"
65 #include "Thyra_DefaultProductVector.hpp"
66 #include "Thyra_BlockedLinearOpBase.hpp"
67 #include "Thyra_get_Epetra_Operator.hpp"
68 
69 #include "Teuchos_FancyOStream.hpp"
70 
71 #include <unordered_map>
72 
73 // **********************************************************************
74 // Specialization: Residual
75 // **********************************************************************
76 
77 
78 template<typename TRAITS,typename LO,typename GO>
81  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & cIndexers,
82  const Teuchos::ParameterList& p,
83  bool useDiscreteAdjoint)
84  : rowIndexers_(rIndexers)
85  , colIndexers_(cIndexers)
86  , globalDataKey_("Residual Scatter Container")
87 {
88  std::string scatterName = p.get<std::string>("Scatter Name");
89  scatterHolder_ =
90  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
91 
92  // get names to be evaluated
93  const std::vector<std::string>& names =
94  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
95 
96  // grab map from evaluated names to field names
97  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
98 
99  // determine if we are scattering an initial condition
100  scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
101 
102  Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
104  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
105  if (!scatterIC_) {
106  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
107  local_side_id_ = p.get<int>("Local Side ID");
108  }
109 
110  // build the vector of fields that this is dependent on
111  scatterFields_.resize(names.size());
112  for (std::size_t eq = 0; eq < names.size(); ++eq) {
113  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
114 
115  // tell the field manager that we depend on this field
116  this->addDependentField(scatterFields_[eq]);
117  }
118 
119  checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
120  if (checkApplyBC_) {
121  applyBC_.resize(names.size());
122  for (std::size_t eq = 0; eq < names.size(); ++eq) {
123  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
124  this->addDependentField(applyBC_[eq]);
125  }
126  }
127 
128  // this is what this evaluator provides
129  this->addEvaluatedField(*scatterHolder_);
130 
131  if (p.isType<std::string>("Global Data Key"))
132  globalDataKey_ = p.get<std::string>("Global Data Key");
133 
134  this->setName(scatterName+" Scatter Residual");
135 }
136 
137 // **********************************************************************
138 template<typename TRAITS,typename LO,typename GO>
140 postRegistrationSetup(typename TRAITS::SetupData d,
142 {
143  indexerIds_.resize(scatterFields_.size());
144  subFieldIds_.resize(scatterFields_.size());
145 
146  // load required field numbers for fast use
147  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
148  // get field ID from DOF manager
149  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
150 
151  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
152  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
153 
154  // fill field data object
155  this->utils.setFieldData(scatterFields_[fd],fm);
156 
157  if (checkApplyBC_)
158  this->utils.setFieldData(applyBC_[fd],fm);
159  }
160 
161  // get the number of nodes (Should be renamed basis)
162  num_nodes = scatterFields_[0].dimension(1);
163 }
164 
165 // **********************************************************************
166 template<typename TRAITS,typename LO,typename GO>
168 preEvaluate(typename TRAITS::PreEvalData d)
169 {
170  typedef BlockedEpetraLinearObjContainer BLOC;
171  typedef BlockedEpetraLinearObjContainer ELOC;
172 
173  using Teuchos::rcp_dynamic_cast;
175 
176  // extract dirichlet counter from container
177  Teuchos::RCP<BLOC> blockContainer
178  = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc.getDataObject("Dirichlet Counter"),true);
179 
180  dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
181  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
182 
183  // extract linear object container
184  Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject(globalDataKey_));
185  Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc.getDataObject(globalDataKey_));
186 
187  // if its blocked do this
188  if(blockedContainer!=Teuchos::null)
189  r_ = (!scatterIC_) ?
190  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
191  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
192  else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
193  r_ = (!scatterIC_) ?
194  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
195  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
196 
197  TEUCHOS_ASSERT(r_!=Teuchos::null);
198 }
199 
200 // **********************************************************************
201 template<typename TRAITS,typename LO,typename GO>
203 evaluateFields(typename TRAITS::EvalData workset)
204 {
205  using Teuchos::RCP;
206  using Teuchos::ArrayRCP;
207  using Teuchos::ptrFromRef;
208  using Teuchos::rcp_dynamic_cast;
209 
210  using Thyra::VectorBase;
211  using Thyra::SpmdVectorBase;
213 
214  // for convenience pull out some objects from workset
215  std::string blockId = this->wda(workset).block_id;
216  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
217 
218  // NOTE: A reordering of these loops will likely improve performance
219  // The "getGIDFieldOffsets may be expensive. However the
220  // "getElementGIDs" can be cheaper. However the lookup for LIDs
221  // may be more expensive!
222 
223  // loop over each field to be scattered
225  Teuchos::ArrayRCP<double> local_dc;
226  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
227  int rowIndexer = indexerIds_[fieldIndex];
228  int subFieldNum = subFieldIds_[fieldIndex];
229 
230  rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
231  ->getNonconstLocalData(ptrFromRef(local_dc));
232 
233  // grab local data for inputing
234  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
235  ->getNonconstLocalData(ptrFromRef(local_r));
236 
237  auto subRowIndexer = rowIndexers_[rowIndexer];
238 
239  // scatter operation for each cell in workset
240  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
241  std::size_t cellLocalId = localCellIds[worksetCellIndex];
242 
243  const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
244 
245  if (!scatterIC_) {
246  // this call "should" get the right ordering according to the Intrepid2 basis
247  const std::pair<std::vector<int>,std::vector<int> > & indicePair
248  = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
249  const std::vector<int> & elmtOffset = indicePair.first;
250  const std::vector<int> & basisIdMap = indicePair.second;
251 
252  // loop over basis functions
253  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
254  int offset = elmtOffset[basis];
255  int lid = LIDs[offset];
256  if(lid<0) // not on this processor!
257  continue;
258 
259  int basisId = basisIdMap[basis];
260 
261  if (checkApplyBC_)
262  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
263  continue;
264 
265  local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
266 
267  // record that you set a dirichlet condition
268  local_dc[lid] = 1.0;
269  }
270  } else {
271  // this call "should" get the right ordering according to the Intrepid2 basis
272  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
273 
274  // loop over basis functions
275  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
276  int offset = elmtOffset[basis];
277  int lid = LIDs[offset];
278  if(lid<0) // not on this processor!
279  continue;
280 
281  local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
282 
283  // record that you set a dirichlet condition
284  local_dc[lid] = 1.0;
285  }
286  }
287  }
288  }
289 }
290 
291 // **********************************************************************
292 // Specialization: Tangent
293 // **********************************************************************
294 
295 
296 template<typename TRAITS,typename LO,typename GO>
299  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & cIndexers,
300  const Teuchos::ParameterList& p,
301  bool useDiscreteAdjoint)
302  : rowIndexers_(rIndexers)
303  , colIndexers_(cIndexers)
304  , globalDataKey_("Residual Scatter Container")
305 {
306  std::string scatterName = p.get<std::string>("Scatter Name");
307  scatterHolder_ =
308  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
309 
310  // get names to be evaluated
311  const std::vector<std::string>& names =
312  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
313 
314  // grab map from evaluated names to field names
315  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
316 
317  // determine if we are scattering an initial condition
318  scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
319 
320  Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
322  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
323  if (!scatterIC_) {
324  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
325  local_side_id_ = p.get<int>("Local Side ID");
326  }
327 
328  // build the vector of fields that this is dependent on
329  scatterFields_.resize(names.size());
330  for (std::size_t eq = 0; eq < names.size(); ++eq) {
331  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
332 
333  // tell the field manager that we depend on this field
334  this->addDependentField(scatterFields_[eq]);
335  }
336 
337  checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
338  if (checkApplyBC_) {
339  applyBC_.resize(names.size());
340  for (std::size_t eq = 0; eq < names.size(); ++eq) {
341  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
342  this->addDependentField(applyBC_[eq]);
343  }
344  }
345 
346  // this is what this evaluator provides
347  this->addEvaluatedField(*scatterHolder_);
348 
349  if (p.isType<std::string>("Global Data Key"))
350  globalDataKey_ = p.get<std::string>("Global Data Key");
351 
352  this->setName(scatterName+" Scatter Tangent");
353 }
354 
355 // **********************************************************************
356 template<typename TRAITS,typename LO,typename GO>
358 postRegistrationSetup(typename TRAITS::SetupData d,
360 {
361  indexerIds_.resize(scatterFields_.size());
362  subFieldIds_.resize(scatterFields_.size());
363 
364  // load required field numbers for fast use
365  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
366  // get field ID from DOF manager
367  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
368 
369  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
370  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
371 
372  // fill field data object
373  this->utils.setFieldData(scatterFields_[fd],fm);
374 
375  if (checkApplyBC_)
376  this->utils.setFieldData(applyBC_[fd],fm);
377  }
378 
379  // get the number of nodes (Should be renamed basis)
380  num_nodes = scatterFields_[0].dimension(1);
381 }
382 
383 // **********************************************************************
384 template<typename TRAITS,typename LO,typename GO>
386 preEvaluate(typename TRAITS::PreEvalData d)
387 {
388  typedef BlockedEpetraLinearObjContainer BLOC;
389  typedef BlockedEpetraLinearObjContainer ELOC;
390 
391  using Teuchos::rcp_dynamic_cast;
393 
394  // extract dirichlet counter from container
395  Teuchos::RCP<BLOC> blockContainer
396  = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc.getDataObject("Dirichlet Counter"),true);
397 
398  dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
399  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
400 
401  // extract linear object container
402  Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject(globalDataKey_));
403  Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc.getDataObject(globalDataKey_));
404 
405  // if its blocked do this
406  if(blockedContainer!=Teuchos::null)
407  r_ = (!scatterIC_) ?
408  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
409  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
410  else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
411  r_ = (!scatterIC_) ?
412  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
413  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
414 
415  TEUCHOS_ASSERT(r_!=Teuchos::null);
416 }
417 
418 // **********************************************************************
419 template<typename TRAITS,typename LO,typename GO>
421 evaluateFields(typename TRAITS::EvalData workset)
422 {
423  TEUCHOS_ASSERT(false);
424 
425  using Teuchos::RCP;
426  using Teuchos::ArrayRCP;
427  using Teuchos::ptrFromRef;
428  using Teuchos::rcp_dynamic_cast;
429 
430  using Thyra::VectorBase;
431  using Thyra::SpmdVectorBase;
433 
434  // for convenience pull out some objects from workset
435  std::string blockId = this->wda(workset).block_id;
436  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
437 
438  // NOTE: A reordering of these loops will likely improve performance
439  // The "getGIDFieldOffsets may be expensive. However the
440  // "getElementGIDs" can be cheaper. However the lookup for LIDs
441  // may be more expensive!
442 
443  // loop over each field to be scattered
445  Teuchos::ArrayRCP<double> local_dc;
446  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
447  int rowIndexer = indexerIds_[fieldIndex];
448  int subFieldNum = subFieldIds_[fieldIndex];
449 
450  rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
451  ->getNonconstLocalData(ptrFromRef(local_dc));
452 
453  // grab local data for inputing
454  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
455  ->getNonconstLocalData(ptrFromRef(local_r));
456 
457  auto subRowIndexer = rowIndexers_[rowIndexer];
458 
459  // scatter operation for each cell in workset
460  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
461  std::size_t cellLocalId = localCellIds[worksetCellIndex];
462 
463  const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
464 
465  if (!scatterIC_) {
466  // this call "should" get the right ordering according to the Intrepid2 basis
467  const std::pair<std::vector<int>,std::vector<int> > & indicePair
468  = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
469  const std::vector<int> & elmtOffset = indicePair.first;
470  const std::vector<int> & basisIdMap = indicePair.second;
471 
472  // loop over basis functions
473  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
474  int offset = elmtOffset[basis];
475  int lid = LIDs[offset];
476  if(lid<0) // not on this processor!
477  continue;
478 
479  int basisId = basisIdMap[basis];
480 
481  if (checkApplyBC_)
482  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
483  continue;
484 
485  local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
486 
487  // record that you set a dirichlet condition
488  local_dc[lid] = 1.0;
489  }
490  } else {
491  // this call "should" get the right ordering according to the Intrepid2 basis
492  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
493 
494  // loop over basis functions
495  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
496  int offset = elmtOffset[basis];
497  int lid = LIDs[offset];
498  if(lid<0) // not on this processor!
499  continue;
500 
501  local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
502 
503  // record that you set a dirichlet condition
504  local_dc[lid] = 1.0;
505  }
506  }
507  }
508  }
509 }
510 
511 // **********************************************************************
512 // Specialization: Jacobian
513 // **********************************************************************
514 
515 template<typename TRAITS,typename LO,typename GO>
518  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & cIndexers,
519  const Teuchos::ParameterList& p,
520  bool useDiscreteAdjoint)
521  : rowIndexers_(rIndexers)
522  , colIndexers_(cIndexers)
523  , globalDataKey_("Residual Scatter Container")
524 {
525  std::string scatterName = p.get<std::string>("Scatter Name");
526  scatterHolder_ =
527  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
528 
529  // get names to be evaluated
530  const std::vector<std::string>& names =
531  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
532 
533  // grab map from evaluated names to field names
534  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
535 
537  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
538 
539  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
540  local_side_id_ = p.get<int>("Local Side ID");
541 
542  // build the vector of fields that this is dependent on
543  scatterFields_.resize(names.size());
544  for (std::size_t eq = 0; eq < names.size(); ++eq) {
545  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
546 
547  // tell the field manager that we depend on this field
548  this->addDependentField(scatterFields_[eq]);
549  }
550 
551  checkApplyBC_ = p.get<bool>("Check Apply BC");
552  if (checkApplyBC_) {
553  applyBC_.resize(names.size());
554  for (std::size_t eq = 0; eq < names.size(); ++eq) {
555  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
556  this->addDependentField(applyBC_[eq]);
557  }
558  }
559 
560  // this is what this evaluator provides
561  this->addEvaluatedField(*scatterHolder_);
562 
563  if (p.isType<std::string>("Global Data Key"))
564  globalDataKey_ = p.get<std::string>("Global Data Key");
565 
566  if(colIndexers_.size()==0)
567  colIndexers_ = rowIndexers_;
568 
569  this->setName(scatterName+" Scatter Residual (Jacobian)");
570 }
571 
572 // **********************************************************************
573 template<typename TRAITS,typename LO,typename GO>
575 postRegistrationSetup(typename TRAITS::SetupData d,
577 {
578  indexerIds_.resize(scatterFields_.size());
579  subFieldIds_.resize(scatterFields_.size());
580 
581  // load required field numbers for fast use
582  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
583  // get field ID from DOF manager
584  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
585 
586  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
587  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
588 
589  // fill field data object
590  this->utils.setFieldData(scatterFields_[fd],fm);
591 
592  if (checkApplyBC_)
593  this->utils.setFieldData(applyBC_[fd],fm);
594  }
595 
596  // get the number of nodes (Should be renamed basis)
597  num_nodes = scatterFields_[0].dimension(1);
598  num_eq = scatterFields_.size();
599 }
600 
601 // **********************************************************************
602 template<typename TRAITS,typename LO,typename GO>
604 preEvaluate(typename TRAITS::PreEvalData d)
605 {
606  typedef BlockedEpetraLinearObjContainer BLOC;
607 
608  using Teuchos::rcp_dynamic_cast;
609 
610  // extract dirichlet counter from container
611  Teuchos::RCP<const BLOC> blockContainer
612  = rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject("Dirichlet Counter"),true);
613 
614  dirichletCounter_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
615  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
616 
617  // extract linear object container
618  blockContainer = rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject(globalDataKey_),true);
619  TEUCHOS_ASSERT(!Teuchos::is_null(blockContainer));
620 
621  r_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f());
622  Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockContainer->get_A());
623 }
624 
625 // **********************************************************************
626 template<typename TRAITS,typename LO,typename GO>
628 evaluateFields(typename TRAITS::EvalData workset)
629 {
630  using Teuchos::RCP;
631  using Teuchos::ArrayRCP;
632  using Teuchos::ptrFromRef;
633  using Teuchos::rcp_dynamic_cast;
634 
635  using Thyra::SpmdVectorBase;
636 
637  // for convenience pull out some objects from workset
638  std::string blockId = this->wda(workset).block_id;
639  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
640 
641  int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
642 
643  std::vector<int> blockOffsets;
644  computeBlockOffsets(blockId,colIndexers_,blockOffsets);
645 
646  std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
647 
648  // NOTE: A reordering of these loops will likely improve performance
649  // The "getGIDFieldOffsets may be expensive. However the
650  // "getElementGIDs" can be cheaper. However the lookup for LIDs
651  // may be more expensive!
652 
653  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
654  int rowIndexer = indexerIds_[fieldIndex];
655  int subFieldNum = subFieldIds_[fieldIndex];
656 
657  // loop over each field to be scattered
658  Teuchos::ArrayRCP<double> local_dc;
659  rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
660  ->getNonconstLocalData(ptrFromRef(local_dc));
661 
662  // grab local data for inputing
664  if(r_!=Teuchos::null)
665  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
666  ->getNonconstLocalData(ptrFromRef(local_r));
667 
668  auto subRowIndexer = rowIndexers_[rowIndexer];
669  auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
670  const std::vector<int> & subElmtOffset = subIndicePair.first;
671  const std::vector<int> & subBasisIdMap = subIndicePair.second;
672 
673  // scatter operation for each cell in workset
674  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
675  std::size_t cellLocalId = localCellIds[worksetCellIndex];
676 
677  const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
678 
679  // loop over basis functions
680  for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
681  int offset = subElmtOffset[basis];
682  int lid = rLIDs[offset];
683  if(lid<0) // not on this processor
684  continue;
685 
686  int basisId = subBasisIdMap[basis];
687 
688  if (checkApplyBC_)
689  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
690  continue;
691 
692  // zero out matrix row
693  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
694  int start = blockOffsets[colIndexer];
695  int end = blockOffsets[colIndexer+1];
696 
697  if(end-start<=0)
698  continue;
699 
700  // check hash table for jacobian sub block
701  std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
702  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
703 // if you didn't find one before, add it to the hash table
704  if(subJac==Teuchos::null) {
705  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
706 
707  // block operator is null, don't do anything (it is excluded)
708  if(Teuchos::is_null(tOp))
709  continue;
710 
711  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
712  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
713  jacEpetraBlocks[blockIndex] = subJac;
714  }
715 
716  int numEntries = 0;
717  int * rowIndices = 0;
718  double * rowValues = 0;
719 
720  subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
721 
722  for(int i=0;i<numEntries;i++)
723  rowValues[i] = 0.0;
724  }
725 
726  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
727 
728  if(r_!=Teuchos::null)
729  local_r[lid] = scatterField.val();
730  local_dc[lid] = 1.0; // mark row as dirichlet
731 
732  // loop over the sensitivity indices: all DOFs on a cell
733  std::vector<double> jacRow(scatterField.size(),0.0);
734 
735  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
736  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
737 
738  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
739  int start = blockOffsets[colIndexer];
740  int end = blockOffsets[colIndexer+1];
741 
742  if(end-start<=0)
743  continue;
744 
745  auto subColIndexer = colIndexers_[colIndexer];
746  const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
747 
748  TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
749 
750  // check hash table for jacobian sub block
751  std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
752  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
753 
754  // if you didn't find one before, add it to the hash table
755  if(subJac==Teuchos::null) {
756  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
757 
758  // block operator is null, don't do anything (it is excluded)
759  if(Teuchos::is_null(tOp))
760  continue;
761 
762  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
763  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
764  jacEpetraBlocks[blockIndex] = subJac;
765  }
766 
767  // Sum Jacobian
768  int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
769  if(err!=0) {
770  std::stringstream ss;
771  ss << "Failed inserting row: " << " (" << lid << "): ";
772  for(int i=0;i<end-start;i++)
773  ss << cLIDs[i] << " ";
774  ss << std::endl;
775  ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
776 
777  ss << "scatter field = ";
778  scatterFields_[fieldIndex].print(ss);
779  ss << std::endl;
780 
781  TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
782  }
783 
784  }
785  }
786  }
787  }
788 }
789 
790 // **********************************************************************
791 
792 #endif
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis, std::vector< int > &blockOffsets)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
PHX::MDField< ScalarT > vector
Pushes residual values into the residual vector for a Newton-based solve.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)