Panzer  Version of the Day
Panzer_ScatterResidual_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_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP
44 #define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_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 
57 #include "Panzer_PureBasis.hpp"
60 #include "Panzer_HashUtils.hpp"
61 
62 #include "Thyra_SpmdVectorBase.hpp"
63 #include "Thyra_ProductVectorBase.hpp"
64 #include "Thyra_DefaultProductVector.hpp"
65 #include "Thyra_BlockedLinearOpBase.hpp"
66 #include "Thyra_DefaultBlockedLinearOp.hpp"
67 #include "Thyra_get_Epetra_Operator.hpp"
68 
69 #include "Phalanx_DataLayout_MDALayout.hpp"
70 
71 #include "Teuchos_FancyOStream.hpp"
72 
73 // **********************************************************************
74 // Specialization: Residual
75 // **********************************************************************
76 
77 template<typename TRAITS,typename LO,typename GO>
80  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & cIndexers,
81  const Teuchos::ParameterList& p,
82  bool useDiscreteAdjoint)
83  : rowIndexers_(rIndexers)
84  , colIndexers_(cIndexers)
85  , globalDataKey_("Residual Scatter Container")
86 {
87  std::string scatterName = p.get<std::string>("Scatter Name");
88  scatterHolder_ =
89  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
90 
91  // get names to be evaluated
92  const std::vector<std::string>& names =
93  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
94 
95  // grab map from evaluated names to field names
96  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
97 
99  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
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  // this is what this evaluator provides
111  this->addEvaluatedField(*scatterHolder_);
112 
113  if (p.isType<std::string>("Global Data Key"))
114  globalDataKey_ = p.get<std::string>("Global Data Key");
115  if (p.isType<bool>("Use Discrete Adjoint"))
116  useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
117 
118  this->setName(scatterName+" Scatter Residual");
119 }
120 
121 // **********************************************************************
122 template<typename TRAITS,typename LO,typename GO>
124 postRegistrationSetup(typename TRAITS::SetupData d,
126 {
127  indexerIds_.resize(scatterFields_.size());
128  subFieldIds_.resize(scatterFields_.size());
129 
130  // load required field numbers for fast use
131  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
132  // get field ID from DOF manager
133  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
134 
135  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
136  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
137 
138  // fill field data object
139  this->utils.setFieldData(scatterFields_[fd],fm);
140  }
141 }
142 
143 // **********************************************************************
144 template<typename TRAITS,typename LO,typename GO>
146 preEvaluate(typename TRAITS::PreEvalData d)
147 {
148  typedef BlockedEpetraLinearObjContainer BLOC;
149  typedef BlockedEpetraLinearObjContainer ELOC;
150 
151  // extract linear object container
152  Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject(globalDataKey_));
153  Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc.getDataObject(globalDataKey_));
154 
155  // if its blocked do this
156  if(blockedContainer!=Teuchos::null)
157  r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
158  else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
159  r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
160 
161  TEUCHOS_ASSERT(r_!=Teuchos::null);
162 }
163 
164 // **********************************************************************
165 template<typename TRAITS,typename LO,typename GO>
167 evaluateFields(typename TRAITS::EvalData workset)
168 {
169  using Teuchos::RCP;
170  using Teuchos::ptrFromRef;
171  using Teuchos::rcp_dynamic_cast;
172 
173  using Thyra::VectorBase;
174  using Thyra::SpmdVectorBase;
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  // NOTE: A reordering of these loops will likely improve performance
182  // The "getGIDFieldOffsets may be expensive. However the
183  // "getElementGIDs" can be cheaper. However the lookup for LIDs
184  // may be more expensive!
185 
186  // loop over each field to be scattered
188  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
189  int indexerId = indexerIds_[fieldIndex];
190  int subFieldNum = subFieldIds_[fieldIndex];
191 
192  // grab local data for inputing
193  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
194 
195  auto subRowIndexer = rowIndexers_[indexerId];
196  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
197 
198  // scatter operation for each cell in workset
199  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
200  std::size_t cellLocalId = localCellIds[worksetCellIndex];
201 
202  const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
203 
204  // loop over basis functions
205  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
206  int offset = elmtOffset[basis];
207  int lid = LIDs[offset];
208  local_r[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis);
209  }
210  }
211  }
212 }
213 
214 // **********************************************************************
215 // Specialization: Tangent
216 // **********************************************************************
217 
218 template<typename TRAITS,typename LO,typename GO>
221  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & cIndexers,
222  const Teuchos::ParameterList& p,
223  bool useDiscreteAdjoint)
224  : rowIndexers_(rIndexers)
225  , colIndexers_(cIndexers)
226  , globalDataKey_("Residual Scatter Container")
227 {
228  std::string scatterName = p.get<std::string>("Scatter Name");
229  scatterHolder_ =
230  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
231 
232  // get names to be evaluated
233  const std::vector<std::string>& names =
234  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
235 
236  // grab map from evaluated names to field names
237  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
238 
240  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
241 
242  // build the vector of fields that this is dependent on
243  scatterFields_.resize(names.size());
244  for (std::size_t eq = 0; eq < names.size(); ++eq) {
245  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
246 
247  // tell the field manager that we depend on this field
248  this->addDependentField(scatterFields_[eq]);
249  }
250 
251  // this is what this evaluator provides
252  this->addEvaluatedField(*scatterHolder_);
253 
254  if (p.isType<std::string>("Global Data Key"))
255  globalDataKey_ = p.get<std::string>("Global Data Key");
256 
257  this->setName(scatterName+" Scatter Tangent");
258 }
259 
260 // **********************************************************************
261 template<typename TRAITS,typename LO,typename GO>
263 postRegistrationSetup(typename TRAITS::SetupData d,
265 {
266  indexerIds_.resize(scatterFields_.size());
267  subFieldIds_.resize(scatterFields_.size());
268 
269  // load required field numbers for fast use
270  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
271  // get field ID from DOF manager
272  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
273 
274  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
275  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
276 
277  // fill field data object
278  this->utils.setFieldData(scatterFields_[fd],fm);
279  }
280 }
281 
282 // **********************************************************************
283 template<typename TRAITS,typename LO,typename GO>
285 preEvaluate(typename TRAITS::PreEvalData d)
286 {
287  typedef BlockedEpetraLinearObjContainer BLOC;
288  typedef BlockedEpetraLinearObjContainer ELOC;
289 
290  // extract linear object container
291  Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject(globalDataKey_));
292  Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc.getDataObject(globalDataKey_));
293 
294  // if its blocked do this
295  if(blockedContainer!=Teuchos::null)
296  r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
297  else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
298  r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
299 
300  TEUCHOS_ASSERT(r_!=Teuchos::null);
301 }
302 
303 // **********************************************************************
304 template<typename TRAITS,typename LO,typename GO>
306 evaluateFields(typename TRAITS::EvalData workset)
307 {
308  TEUCHOS_ASSERT(false);
309 
310  using Teuchos::RCP;
311  using Teuchos::ptrFromRef;
312  using Teuchos::rcp_dynamic_cast;
313 
314  using Thyra::VectorBase;
315  using Thyra::SpmdVectorBase;
317 
318  // for convenience pull out some objects from workset
319  std::string blockId = this->wda(workset).block_id;
320  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
321 
322  // NOTE: A reordering of these loops will likely improve performance
323  // The "getGIDFieldOffsets may be expensive. However the
324  // "getElementGIDs" can be cheaper. However the lookup for LIDs
325  // may be more expensive!
326 
327  // loop over each field to be scattered
329  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
330  int indexerId = indexerIds_[fieldIndex];
331  int subFieldNum = subFieldIds_[fieldIndex];
332 
333  // grab local data for inputing
334  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
335 
336  auto subRowIndexer = rowIndexers_[indexerId];
337  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
338 
339  // scatter operation for each cell in workset
340  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
341  std::size_t cellLocalId = localCellIds[worksetCellIndex];
342 
343  const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
344 
345  // loop over basis functions
346  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
347  int offset = elmtOffset[basis];
348  int lid = LIDs[offset];
349  local_r[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
350  }
351  }
352  }
353 }
354 
355 // **********************************************************************
356 // Specialization: Jacobian
357 // **********************************************************************
358 
359 template<typename TRAITS,typename LO,typename GO>
362  const std::vector<Teuchos::RCP<const UniqueGlobalIndexer<LO,int> > > & cIndexers,
363  const Teuchos::ParameterList& p,
364  bool useDiscreteAdjoint)
365  : rowIndexers_(rIndexers)
366  , colIndexers_(cIndexers)
367  , globalDataKey_("Residual Scatter Container")
368  , useDiscreteAdjoint_(useDiscreteAdjoint)
369 {
370  std::string scatterName = p.get<std::string>("Scatter Name");
371  scatterHolder_ =
372  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
373 
374  // get names to be evaluated
375  const std::vector<std::string>& names =
376  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
377 
378  // grab map from evaluated names to field names
379  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
380 
382  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
383 
384  // build the vector of fields that this is dependent on
385  scatterFields_.resize(names.size());
386  for (std::size_t eq = 0; eq < names.size(); ++eq) {
387  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
388 
389  // tell the field manager that we depend on this field
390  this->addDependentField(scatterFields_[eq]);
391  }
392 
393  // this is what this evaluator provides
394  this->addEvaluatedField(*scatterHolder_);
395 
396  if (p.isType<std::string>("Global Data Key"))
397  globalDataKey_ = p.get<std::string>("Global Data Key");
398  if (p.isType<bool>("Use Discrete Adjoint"))
399  useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
400 
401  // discrete adjoint does not work with non-square matrices
402  if(useDiscreteAdjoint)
403  { TEUCHOS_ASSERT(colIndexers_.size()==0); }
404 
405  if(colIndexers_.size()==0)
406  colIndexers_ = rowIndexers_;
407 
408  this->setName(scatterName+" Scatter Residual BlockedEpetra (Jacobian)");
409 }
410 
411 // **********************************************************************
412 template<typename TRAITS,typename LO,typename GO>
414 postRegistrationSetup(typename TRAITS::SetupData d,
416 {
417  indexerIds_.resize(scatterFields_.size());
418  subFieldIds_.resize(scatterFields_.size());
419 
420  // load required field numbers for fast use
421  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
422  // get field ID from DOF manager
423  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
424 
425  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
426  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
427 
428  // fill field data object
429  this->utils.setFieldData(scatterFields_[fd],fm);
430  }
431 }
432 
433 // **********************************************************************
434 template<typename TRAITS,typename LO,typename GO>
436 preEvaluate(typename TRAITS::PreEvalData d)
437 {
438  using Teuchos::RCP;
439  using Teuchos::rcp_dynamic_cast;
440 
441  typedef BlockedEpetraLinearObjContainer BLOC;
442  typedef BlockedEpetraLinearObjContainer ELOC;
443 
444  // extract linear object container
445  RCP<const BLOC> blockedContainer = rcp_dynamic_cast<const BLOC>(d.gedc.getDataObject(globalDataKey_));
446  RCP<const ELOC> epetraContainer = rcp_dynamic_cast<const ELOC>(d.gedc.getDataObject(globalDataKey_));
447 
448  // if its blocked do this
449  if(blockedContainer!=Teuchos::null) {
450  r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
451  Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockedContainer->get_A());
452  }
453  else if(epetraContainer!=Teuchos::null) {
454  // if its straight up epetra do this
455  if(epetraContainer->get_f_th()!=Teuchos::null)
456  r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
457 
458  // convert it into a blocked operator
459  RCP<Thyra::LinearOpBase<double> > J = blockedContainer->get_A_th();
460  Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(Thyra::nonconstBlock1x1(J));
461  }
462 
463  TEUCHOS_ASSERT(Jac_!=Teuchos::null);
464 }
465 
466 // **********************************************************************
467 template<typename TRAITS,typename LO,typename GO>
469 evaluateFields(typename TRAITS::EvalData workset)
470 {
471  using Teuchos::RCP;
472  using Teuchos::ArrayRCP;
473  using Teuchos::ptrFromRef;
474  using Teuchos::rcp_dynamic_cast;
475 
476  using Thyra::VectorBase;
477  using Thyra::SpmdVectorBase;
480 
481  std::vector<double> jacRow;
482 
483  // for convenience pull out some objects from workset
484  std::string blockId = this->wda(workset).block_id;
485  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
486 
487  int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
488 
489  std::vector<int> blockOffsets;
490  computeBlockOffsets(blockId,colIndexers_,blockOffsets);
491 
492  std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
493 
494  // loop over each field to be scattered
496  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
497  int rowIndexer = indexerIds_[fieldIndex];
498  int subFieldNum = subFieldIds_[fieldIndex];
499 
500  // grab local data for inputing
501  if(r_!=Teuchos::null)
502  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))->getNonconstLocalData(ptrFromRef(local_r));
503 
504  auto subRowIndexer = rowIndexers_[rowIndexer];
505  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
506 
507  // scatter operation for each cell in workset
508  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
509  std::size_t cellLocalId = localCellIds[worksetCellIndex];
510 
511  const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
512 
513  // loop over the basis functions (currently they are nodes)
514  for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
515  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
516  int rowOffset = elmtOffset[rowBasisNum];
517  int r_lid = rLIDs[rowOffset];
518 
519  // Sum residual
520  if(local_r!=Teuchos::null)
521  local_r[r_lid] += (scatterField.val());
522 
523  // loop over the sensitivity indices: all DOFs on a cell
524  jacRow.resize(scatterField.size());
525 
526  // For Neumann conditions with no dependence on degrees of freedom, there should be no Jacobian contribution
527  if(scatterField.size() == 0)
528  continue;
529 
530  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
531  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
532 
533  // scatter the row to each block
534  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
535  int start = blockOffsets[colIndexer];
536  int end = blockOffsets[colIndexer+1];
537 
538  if(end-start<=0)
539  continue;
540 
541  auto subColIndexer = colIndexers_[colIndexer];
542  const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
543 
544  TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
545 
546  // check hash table for jacobian sub block
547  std::pair<int,int> blockIndex = useDiscreteAdjoint_ ? std::make_pair(colIndexer,rowIndexer)
548  : std::make_pair(rowIndexer,colIndexer);
549  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
550 
551  // if you didn't find one before, add it to the hash table
552  if(subJac==Teuchos::null) {
553  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
554 
555  // block operator is null, don't do anything (it is excluded)
556  if(Teuchos::is_null(tOp))
557  continue;
558 
559  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
560  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
561  jacEpetraBlocks[blockIndex] = subJac;
562  }
563 
564  // Sum Jacobian
565  if(!useDiscreteAdjoint_) {
566  int err = subJac->SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs[0]);
567  if(err!=0) {
568 
569  std::stringstream ss;
570  ss << "Failed inserting row: " << "LID = " << r_lid << "): ";
571  for(int i=start;i<end;i++)
572  ss << cLIDs[i] << " ";
573  ss << std::endl;
574  ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
575 
576  ss << "scatter field = ";
577  scatterFields_[fieldIndex].print(ss);
578  ss << std::endl;
579 
580  TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
581  }
582  }
583  else {
584  for(std::size_t c=0;c<cLIDs.size();c++) {
585  int err = subJac->SumIntoMyValues(cLIDs[c], 1, &jacRow[start+c],&r_lid);
587  }
588  }
589  }
590  } // end rowBasisNum
591  } // end fieldIndex
592  }
593 }
594 
595 // **********************************************************************
596 
597 #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)
PHX::MDField< ScalarT > vector
void resize(const size_type n, const T &val=T())
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isType(const std::string &name) const
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Pushes residual values into the residual vector for a Newton-based solve.
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)