Panzer  Version of the Day
Panzer_ScatterDirichletResidual_Epetra_Hessian_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_ScatterDirichletResidual_Epetra_Hessian_impl_hpp__
44 #define __Panzer_ScatterDirichletResidual_Epetra_Hessian_impl_hpp__
45 
46 // only do this if required by the user
47 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
48 
50 
51 // the includes for this file come in as a result of the includes in the main
52 // Epetra scatter dirichlet residual file
53 
54 namespace panzer {
55 
56 // **************************************************************
57 // Hessian Specialization
58 // **************************************************************
59 template<typename TRAITS,typename LO,typename GO>
60 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
61 ScatterDirichletResidual_Epetra(const Teuchos::RCP<const UniqueGlobalIndexer<LO,GO> > & indexer,
62  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & cIndexer,
63  const Teuchos::ParameterList& p)
64  : globalIndexer_(indexer)
65  , colGlobalIndexer_(cIndexer)
66  , globalDataKey_("Residual Scatter Container")
67  , preserveDiagonal_(false)
68 {
69  std::string scatterName = p.get<std::string>("Scatter Name");
70  scatterHolder_ =
71  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
72 
73  // get names to be evaluated
74  const std::vector<std::string>& names =
75  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
76 
77  // grab map from evaluated names to field names
78  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
79 
81  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
82 
83  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
84  local_side_id_ = p.get<int>("Local Side ID");
85 
86  // build the vector of fields that this is dependent on
87  scatterFields_.resize(names.size());
88  for (std::size_t eq = 0; eq < names.size(); ++eq) {
89  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
90 
91  // tell the field manager that we depend on this field
92  this->addDependentField(scatterFields_[eq]);
93  }
94 
95  checkApplyBC_ = p.get<bool>("Check Apply BC");
96  if (checkApplyBC_) {
97  applyBC_.resize(names.size());
98  for (std::size_t eq = 0; eq < names.size(); ++eq) {
99  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
100  this->addDependentField(applyBC_[eq]);
101  }
102  }
103 
104  // this is what this evaluator provides
105  this->addEvaluatedField(*scatterHolder_);
106 
107  if (p.isType<std::string>("Global Data Key"))
108  globalDataKey_ = p.get<std::string>("Global Data Key");
109 
110  if (p.isType<bool>("Preserve Diagonal"))
111  preserveDiagonal_ = p.get<bool>("Preserve Diagonal");
112 
113  this->setName(scatterName+" Scatter Residual (Hessian)");
114 }
115 
116 template<typename TRAITS,typename LO,typename GO>
117 void
118 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
119 postRegistrationSetup(typename TRAITS::SetupData d,
121 {
122  fieldIds_.resize(scatterFields_.size());
123 
124  // load required field numbers for fast use
125  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
126  // get field ID from DOF manager
127  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
128  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
129  // fill field data object
130  this->utils.setFieldData(scatterFields_[fd],fm);
131 
132  if (checkApplyBC_)
133  this->utils.setFieldData(applyBC_[fd],fm);
134  }
135 
136  // get the number of nodes (Should be renamed basis)
137  num_nodes = scatterFields_[0].dimension(1);
138  num_eq = scatterFields_.size();
139 }
140 
141 template<typename TRAITS,typename LO,typename GO>
142 void
143 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
144 preEvaluate(typename TRAITS::PreEvalData d)
145 {
146  // extract linear object container
147  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc.getDataObject(globalDataKey_));
148 
149  if(epetraContainer_==Teuchos::null) {
150  // extract linear object container
151  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),true)->getGhostedLOC();
152  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc,true);
153 
154  dirichletCounter_ = Teuchos::null;
155  }
156  else {
157  // extract dirichlet counter from container
158  Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc.getDataObject("Dirichlet Counter");
159  Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(dataContainer,true);
160 
161  dirichletCounter_ = epetraContainer->get_f();
162  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
163  }
164 }
165 
166 template<typename TRAITS,typename LO,typename GO>
167 void
168 ScatterDirichletResidual_Epetra<panzer::Traits::Hessian,TRAITS,LO,GO>::
169 evaluateFields(typename TRAITS::EvalData workset)
170 {
171  // TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
172  // "ScatterDirichletResidual_Epetra<Hessian> is not yet implemented"); // just in case
173 
174  std::vector<int> cLIDs, rLIDs;
175  std::vector<double> jacRow;
176 
177  bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
178 
179  // for convenience pull out some objects from workset
180  std::string blockId = this->wda(workset).block_id;
181  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
182 
183  Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
184  TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
185  Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
186  // NOTE: A reordering of these loops will likely improve performance
187  // The "getGIDFieldOffsets may be expensive. However the
188  // "getElementGIDs" can be cheaper. However the lookup for LIDs
189  // may be more expensive!
190 
191  // scatter operation for each cell in workset
192  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
193  std::size_t cellLocalId = localCellIds[worksetCellIndex];
194 
195  rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
196  if(useColumnIndexer)
197  cLIDs = colGlobalIndexer_->getElementLIDs(cellLocalId);
198  else
199  cLIDs = rLIDs;
200 
201  // loop over each field to be scattered
202  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
203  int fieldNum = fieldIds_[fieldIndex];
204 
205  // this call "should" get the right ordering according to the Intrepid2 basis
206  const std::pair<std::vector<int>,std::vector<int> > & indicePair
207  = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
208  const std::vector<int> & elmtOffset = indicePair.first;
209  const std::vector<int> & basisIdMap = indicePair.second;
210 
211  // loop over basis functions
212  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
213  int offset = elmtOffset[basis];
214  int row = rLIDs[offset];
215  if(row<0) // not on this processor
216  continue;
217 
218  int basisId = basisIdMap[basis];
219 
220  if (checkApplyBC_)
221  if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
222  continue;
223 
224  // zero out matrix row
225  {
226  int numEntries = 0;
227  int * rowIndices = 0;
228  double * rowValues = 0;
229 
230  Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
231 
232  for(int i=0;i<numEntries;i++) {
233  if(preserveDiagonal_) {
234  if(row!=rowIndices[i])
235  rowValues[i] = 0.0;
236  }
237  else
238  rowValues[i] = 0.0;
239  }
240  }
241 
242  // int gid = GIDs[offset];
243  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
244 
245  if(dirichletCounter_!=Teuchos::null) {
246  // std::cout << "Writing " << row << " " << dirichletCounter_->MyLength() << std::endl;
247  (*dirichletCounter_)[row] = 1.0; // mark row as dirichlet
248  }
249 
250  // loop over the sensitivity indices: all DOFs on a cell
251  jacRow.resize(scatterField.size());
252  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
253  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex).fastAccessDx(0);
254 
255  if(!preserveDiagonal_) {
256  int err = Jac->ReplaceMyValues(row, cLIDs.size(),
257  panzer::ptrFromStlVector(jacRow),
258  panzer::ptrFromStlVector(cLIDs));
259  TEUCHOS_ASSERT(err==0);
260  }
261  }
262  }
263  }
264 }
265 
266 }
267 
268 // **************************************************************
269 #endif
270 
271 #endif
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T * ptrFromStlVector(std::vector< T > &v)
bool isType(const std::string &name) const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
#define TEUCHOS_ASSERT(assertion_test)
const Teuchos::RCP< Epetra_Vector > get_f() const