Panzer  Version of the Day
Panzer_FieldAggPattern.cpp
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 
44 
46 
47 using Teuchos::RCP;
48 using Teuchos::rcp;
49 
50 namespace panzer {
51 
53 {
54 }
55 
56 FieldAggPattern::FieldAggPattern(std::vector<std::pair<int,Teuchos::RCP<const FieldPattern> > > & patterns,
57  const Teuchos::RCP<const FieldPattern> & geomAggPattern)
58 {
59  buildPattern(patterns,geomAggPattern);
60 }
61 
62 Teuchos::RCP<const FieldPattern> FieldAggPattern::getGeometricAggFieldPattern() const
63 {
64  return geomAggPattern_;
65 }
66 
67 void FieldAggPattern::buildPattern(const std::vector<std::pair<int,Teuchos::RCP<const FieldPattern> > > & patterns,
68  const Teuchos::RCP<const FieldPattern> & geomAggPattern)
69 {
70  TEUCHOS_ASSERT(patterns.size()>0);
71 
72  // build geometric information
73  if(geomAggPattern==Teuchos::null) {
74  std::vector<FPPtr> patternVec;
75 
76  // convert map into vector
77  std::vector<std::pair<int,FPPtr> >::const_iterator itr;
78  for(itr=patterns.begin();itr!=patterns.end();++itr)
79  patternVec.push_back(itr->second);
80 
81  // build geometric aggregate field pattern
82  geomAggPattern_ = rcp(new GeometricAggFieldPattern(patternVec));
83  }
84  else
85  geomAggPattern_ = geomAggPattern;
86 
87  patterns_ = patterns;
88 
89  buildFieldIdToPatternIdx(); // builds look up from fieldId to index into patterns_ vector
90  buildFieldIdsVector(); // meat of matter: build field pattern information
91  buildFieldPatternData(); // this hanldes the "getSubcell*" information
92 }
93 
95 {
96  FPPtr geomPattern = getGeometricAggFieldPattern();
97  TEUCHOS_TEST_FOR_EXCEPTION(geomPattern==Teuchos::null,std::logic_error,
98  "Geometric field pattern not yet set, call buildPatterns first");
99 
100  return geomPattern->getDimension();
101 }
102 
103 shards::CellTopology FieldAggPattern::getCellTopology() const
104 {
105  FPPtr geomPattern = getGeometricAggFieldPattern();
106  TEUCHOS_TEST_FOR_EXCEPTION(geomPattern==Teuchos::null,std::logic_error,
107  "Geometric field pattern not yet set, call buildPatterns first");
108 
109  return geomPattern->getCellTopology();
110 }
111 
113 {
114  return patternData_[dimension].size();
115 }
116 
117 const std::vector<int> & FieldAggPattern::getSubcellIndices(int dimension, int subcell) const
118 {
119  return patternData_[dimension][subcell];
120 }
121 
122 void FieldAggPattern::getSubcellClosureIndices(int, int, std::vector<int> &) const
123 {
124  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
125  "FieldAggPattern::getSubcellClosureIndices should not be called");
126 }
127 
128 void FieldAggPattern::print(std::ostream & os) const
129 {
131 
132  os << "FieldPattern: FieldAggPattern" << std::endl;
133  os << "FieldPattern: |numFields| = " << numFields_.size() << std::endl;
134  os << "FieldPattern: numFields = [ ";
135  int total=0;
136  for(std::size_t i=0;i<numFields_.size();i++) {
137  os << numFields_[i] << " ";
138  total += numFields_[i];
139  }
140  os << "]" << std::endl;
141  os << "FieldPattern: |fieldIds| = " << fieldIds_.size() << " (" << total << ")" << std::endl;
142  os << "FieldPattern: fieldIds = [ ";
143  for(std::size_t i=0;i<fieldIds_.size();i++)
144  os << fieldIds_[i] << " ";
145  os << "]" << std::endl;
146  os << "FieldPattern: local offsets\n";
147 
148  std::map<int,int>::const_iterator itr;
149  for(itr=fieldIdToPatternIdx_.begin();itr!=fieldIdToPatternIdx_.end();++itr) {
150  int fieldId = itr->first;
151  const std::vector<int> & offsets = localOffsets(fieldId);
152  os << "FieldPattern: field " << itr->first << " = [ ";
153  for(std::size_t i=0;i<offsets.size();i++)
154  os << offsets[i] << " ";
155  os << "]" << std::endl;
156  }
157 }
158 
159 Teuchos::RCP<const FieldPattern> FieldAggPattern::getFieldPattern(int fieldId) const
160 {
161  std::map<int,int>::const_iterator idxIter = fieldIdToPatternIdx_.find(fieldId);
162  TEUCHOS_TEST_FOR_EXCEPTION(idxIter==fieldIdToPatternIdx_.end(),std::logic_error,
163  "FieldID = " << fieldId << " not defined in this pattern");
164 
165  return patterns_[idxIter->second].second;
166 }
167 
169 {
170  // convert map into vector
171  int index = 0;
172  std::vector<std::pair<int,FPPtr> >::const_iterator itr;
173  for(itr=patterns_.begin();itr!=patterns_.end();++itr,++index)
174  fieldIdToPatternIdx_[itr->first] = index;
175 }
176 
178 {
179  FPPtr geomAggPattern = getGeometricAggFieldPattern();
180 
181  // numFields_: stores number of fields per geometric ID
182  // fieldIds_: stores field IDs for each entry field indicated by numFields_
183  numFields_.resize(geomAggPattern->numberIds());
184  fieldIds_.clear();
185 
186  // build numFields_ and fieldIds_ vectors
187  for(int dim=0;dim<getDimension()+1;dim++) {
188  int numSubcell = geomAggPattern->getSubcellCount(dim);
189  for(int sc=0;sc<numSubcell;sc++) {
190  // merge the field patterns for multiple fields
191  // on a specific dimension and subcell
192  mergeFieldPatterns(dim,sc);
193  }
194  }
195 }
196 
197 void FieldAggPattern::mergeFieldPatterns(int dim,int subcell)
198 {
199  // make sure that the geometric IDs count is equal to the maxDOFs
200  const std::vector<int> & geomIndices = getGeometricAggFieldPattern()->getSubcellIndices(dim,subcell);
201 
202  // a single geometric entity is assigned
203  // indices size should be either 1 or 0.
204  TEUCHOS_ASSERT(geomIndices.size()<=1);
205 
206  if (geomIndices.size() > 0) {
207  const int geomIndex = geomIndices[0];
208 
209  numFields_[geomIndex] = 0; // no fields initially
210 
211  std::vector<std::pair<int,FPPtr> >::const_iterator itr;
212  for(itr=patterns_.begin();itr!=patterns_.end();++itr) {
213  std::size_t fieldSize = itr->second->getSubcellIndices(dim,subcell).size();
214 
215  // add field ID if their are enough in current pattern
216  for (std::size_t i=0;i<fieldSize;++i)
217  fieldIds_.push_back(itr->first);
218  numFields_[geomIndex]+=fieldSize;
219  }
220  TEUCHOS_ASSERT(numFields_[geomIndex]>=0);
221  }
222 }
223 
225 {
226  int dimension = getDimension();
227  FPPtr geomIdsPattern = getGeometricAggFieldPattern();
228 
229  // build patternData_ vector for implementation of the FieldPattern
230  // functions (indicies will index into fieldIds_)
231  patternData_.resize(dimension+1);
232  int nextIndex = 0;
233  for(int d=0;d<dimension+1;d++) {
234  int numSubcell = geomIdsPattern->getSubcellCount(d);
235  patternData_[d].resize(numSubcell);
236 
237  // loop through sub cells
238  for(int sc=0;sc<numSubcell;sc++) {
239  // a single geometric entity is assigned to field pattern
240  // geomIds size should be either 0 or 1.
241  const std::vector<int> & geomIds = geomIdsPattern->getSubcellIndices(d,sc);
242  TEUCHOS_ASSERT(geomIds.size()<=1);
243  if (geomIds.size() > 0) {
244  const int geomId = geomIds[0];
245  const int numFields = numFields_[geomId];
246  for(int k=0;k<numFields;k++,nextIndex++)
247  patternData_[d][sc].push_back(nextIndex);
248  }
249  }
250  }
251 }
252 
253 const std::vector<int> & FieldAggPattern::localOffsets(int fieldId) const
254 {
255  // lazy evaluation
256  std::map<int,std::vector<int> >::const_iterator itr = fieldOffsets_.find(fieldId);
257  if(itr!=fieldOffsets_.end())
258  return itr->second;
259 
260  std::vector<int> & offsets = fieldOffsets_[fieldId];
261  localOffsets_build(fieldId,offsets);
262  return offsets;
263 }
264 
265 bool FieldAggPattern::LessThan::operator()(const Teuchos::Tuple<int,3> & a,const Teuchos::Tuple<int,3> & b) const
266 {
267  if(a[0] < b[0]) return true;
268  if(a[0] > b[0]) return false;
269 
270  // a[0]==b[0]
271  if(a[1] < b[1]) return true;
272  if(a[1] > b[1]) return false;
273 
274  // a[1]==b[1] && a[0]==b[0]
275  if(a[2] < b[2]) return true;
276  if(a[2] > b[2]) return false;
277 
278  // a[2]==b[2] && a[1]==b[1] && a[0]==b[0]
279  return false; // these are equal to, but not less than!
280 }
281 
282 //const std::vector<int> &
283 const std::pair<std::vector<int>,std::vector<int> > &
284 FieldAggPattern::localOffsets_closure(int fieldId,int subcellDim,int subcellId) const
285 {
286  // lazy evaluation
287  typedef std::map<Teuchos::Tuple<int,3>, std::pair<std::vector<int>,std::vector<int> >,LessThan> OffsetMap;
288 
289  Teuchos::Tuple<int,3> subcellTuple = Teuchos::tuple<int>(fieldId,subcellDim,subcellId);
290 
291  OffsetMap::const_iterator itr
292  = fieldSubcellOffsets_closure_.find(subcellTuple);
293  if(itr!=fieldSubcellOffsets_closure_.end()) {
294  return itr->second;
295  }
296 
297  TEUCHOS_TEST_FOR_EXCEPTION(subcellDim>=getDimension(),std::logic_error,
298  "FieldAggPattern::localOffsets_closure precondition subcellDim<getDimension() failed");
299  TEUCHOS_TEST_FOR_EXCEPTION(subcellId<0,std::logic_error,
300  "FieldAggPattern::localOffsets_closure precondition subcellId>=0 failed");
301  TEUCHOS_TEST_FOR_EXCEPTION(subcellId>=getSubcellCount(subcellDim),std::logic_error,
302  "FieldAggPattern::localOffsets_closure precondition subcellId<getSubcellCount(subcellDim) failed");
303 
304  // build vector for sub cell closure indices
306 
307  // grab field offsets
308  const std::vector<int> & fieldOffsets = localOffsets(fieldId);
309 
310  // get offsets into field offsets for the closure indices (from field pattern)
311  std::vector<int> closureOffsets;
312  FPPtr fieldPattern = getFieldPattern(fieldId);
313  fieldPattern->getSubcellClosureIndices(subcellDim,subcellId,closureOffsets);
314 
315  // build closure indices into the correct location in lazy evaluation map.
316  std::pair<std::vector<int>,std::vector<int> > & indicesPair
317  = fieldSubcellOffsets_closure_[subcellTuple];
318 
319  std::vector<int> & closureIndices = indicesPair.first;
320  for(std::size_t i=0;i<closureOffsets.size();i++)
321  closureIndices.push_back(fieldOffsets[closureOffsets[i]]);
322 
323  std::vector<int> & basisIndices = indicesPair.second;
324  basisIndices.assign(closureOffsets.begin(),closureOffsets.end());
325 
326  TEUCHOS_ASSERT(fieldSubcellOffsets_closure_[subcellTuple].first.size()==closureIndices.size());
327  return fieldSubcellOffsets_closure_[subcellTuple];
328 }
329 
330 void FieldAggPattern::localOffsets_build(int fieldId,std::vector<int> & offsets) const
331 {
332  // This function makes the assumption that if there are multiple indices
333  // shared by a subcell then the GeometricAggFieldPattern reflects this.
334  // This is a fine assumption on edges and faces because the symmetries require
335  // extra information about ordering. However, on nodes and "volumes" the
336  // assumption appears to be stupid. For consistency we will make it.
337  //
338  // This code needs to be tested with a basis that has multple IDs at
339  // a node or "volume" sub cell.
340 
341  FPPtr fieldPattern = getFieldPattern(fieldId);
342 
343  offsets.clear();
344  offsets.resize(fieldPattern->numberIds(),-111111); // fill with some negative number
345  // for testing
346 
347  // this will offsets for all IDs associated with the field
348  // but using a geometric ordering
349  std::vector<int> fieldIdsGeomOrder;
350  for(std::size_t i=0;i<fieldIds_.size();++i) {
351  if(fieldIds_[i]==fieldId)
352  fieldIdsGeomOrder.push_back(i);
353  }
354  TEUCHOS_ASSERT((int) fieldIdsGeomOrder.size()==fieldPattern->numberIds());
355 
356  // built geometric ordering for this pattern: will index into fieldIdsGeomOrder
357  GeometricAggFieldPattern geomPattern(fieldPattern);
358  int cnt = 0;
359  for(int dim=0;dim<geomPattern.getDimension()+1;dim++) {
360  for(int sc=0;sc<geomPattern.getSubcellCount(dim);sc++) {
361  const std::vector<int> & fIndices = fieldPattern->getSubcellIndices(dim,sc);
362 
363  for(std::size_t i=0;i<fIndices.size();i++)
364  offsets[fIndices[i]] = fieldIdsGeomOrder[cnt++];
365  }
366  }
367 
368  // check for failure/bug
369  for(std::size_t i=0;i<offsets.size();i++) {
370  TEUCHOS_ASSERT(offsets[i]>=0);
371  }
372 }
373 
374 }
void mergeFieldPatterns(int dim, int subcell)
std::map< int, std::vector< int > > fieldOffsets_
std::map< Teuchos::Tuple< int, 3 >, std::pair< std::vector< int >, std::vector< int > >, LessThan > fieldSubcellOffsets_closure_
std::map< int, int > fieldIdToPatternIdx_
std::vector< std::pair< int, Teuchos::RCP< const FieldPattern > > > patterns_
virtual const std::vector< int > & getSubcellIndices(int dimension, int subcell) const
std::vector< std::vector< std::vector< int > > > patternData_
virtual Teuchos::RCP< const FieldPattern > getFieldPattern(int fieldId) const
bool operator()(const Teuchos::Tuple< int, 3 > &a, const Teuchos::Tuple< int, 3 > &b) const
Teuchos::RCP< const FieldPattern > geomAggPattern_
PHX::MDField< ScalarT > vector
int numFields
virtual void buildPattern(const std::vector< std::pair< int, Teuchos::RCP< const FieldPattern > > > &patterns, const Teuchos::RCP< const FieldPattern > &geomAggPattern=Teuchos::null)
void localOffsets_build(int fieldId, std::vector< int > &offsets) const
virtual int getDimension() const
virtual Teuchos::RCP< const FieldPattern > getGeometricAggFieldPattern() const
virtual void print(std::ostream &os) const
virtual void print(std::ostream &os) const
Print this pattern.
const std::vector< int > & localOffsets(int fieldId) const
virtual void getSubcellClosureIndices(int, int, std::vector< int > &) const
virtual shards::CellTopology getCellTopology() const
virtual int getSubcellCount(int dimension) const
Teuchos::RCP< const FieldPattern > FPPtr
Kokkos::View< const int *, PHX::Device > offsets
const std::pair< std::vector< int >, std::vector< int > > & localOffsets_closure(int fieldId, int subcellDim, int subcellId) const