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