Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
57  const Teuchos::RCP<const FieldPattern> & geomAggPattern)
58 {
59  buildPattern(patterns,geomAggPattern);
60 }
61 
63 {
64  return geomAggPattern_;
65 }
66 
67 void FieldAggPattern::buildPattern(const std::vector<std::tuple<int,panzer::FieldType,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<std::pair<FieldType,FPPtr>> patternVec;
75 
76  // convert map into vector
77  auto itr = patterns.cbegin();
78  for(;itr!=patterns.cend();++itr)
79  patternVec.push_back(std::make_pair(std::get<1>(*itr),std::get<2>(*itr)));
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 handles 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 
112 int FieldAggPattern::getSubcellCount(int dimension) const
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 
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 std::get<2>(patterns_[idxIter->second]);
166 }
167 
169 {
170  std::map<int,int>::const_iterator idxIter = fieldIdToPatternIdx_.find(fieldId);
171  TEUCHOS_TEST_FOR_EXCEPTION(idxIter==fieldIdToPatternIdx_.end(),std::logic_error,
172  "FieldID = " << fieldId << " not defined in this pattern");
173 
174  return std::get<1>(patterns_[idxIter->second]);
175 }
176 
178 {
179  // convert map into vector
180  int index = 0;
181  auto itr = patterns_.cbegin();
182  for(;itr!=patterns_.cend();++itr,++index)
183  fieldIdToPatternIdx_[std::get<0>(*itr)] = index;
184 }
185 
187 {
188  FPPtr geomAggPattern = getGeometricAggFieldPattern();
189 
190  // numFields_: stores number of fields per geometric ID
191  // fieldIds_: stores field IDs for each entry field indicated by numFields_
192  numFields_.resize(geomAggPattern->numberIds(),0);
193  fieldIds_.clear();
194 
195  // build numFields_ and fieldIds_ vectors
196 
197  // Merge the field patterns for multiple fields for each dimension
198  // and subcell. We do all the CG first, then all DG. This allows us
199  // to use one offset for mapping DOFs to subcells later on.
202 }
203 
205 {
206  auto geomAggPattern = getGeometricAggFieldPattern();
207 
208  // For DG, all DOFs are added to the internal cell
209  const int cellDim = this->getDimension();
210  const int numDimensions = getDimension()+1;
211 
212  for(int dim=0;dim<numDimensions;dim++) {
213  int numSubcell = geomAggPattern->getSubcellCount(dim);
214  for(int subcell=0;subcell<numSubcell;subcell++) {
215 
216  // Get the geometric index to add the DOF indices to. CG adds to
217  // the subcell we are iterating over, (dim,subcel), while DG
218  // adds to the internal cell (cellDim,0)
219  const std::vector<int> * geomIndices = nullptr;
220  if (fieldType == FieldType::CG)
221  geomIndices = &(getGeometricAggFieldPattern()->getSubcellIndices(dim,subcell));
222  else
223  geomIndices = &(getGeometricAggFieldPattern()->getSubcellIndices(cellDim,0));
224 
225  if (geomIndices->size() > 0) {
226  const int geomIndex = (*geomIndices)[0];
227 
228  auto itr = patterns_.cbegin();
229  for(;itr!=patterns_.cend();++itr) {
230  // Only merge in if pattern matches the FieldType.
231  if (std::get<1>(*itr) == fieldType) {
232  const std::size_t fieldSize = std::get<2>(*itr)->getSubcellIndices(dim,subcell).size();
233 
234  // add field ID if there are enough in current pattern
235  for (std::size_t i=0;i<fieldSize;++i)
236  fieldIds_.push_back(std::get<0>(*itr));
237  numFields_[geomIndex]+=fieldSize;
238  }
239  }
240  TEUCHOS_ASSERT(numFields_[geomIndex]>=0);
241  }
242 
243  }
244  }
245 
246 }
247 
249 {
250  int dimension = getDimension();
251  FPPtr geomIdsPattern = getGeometricAggFieldPattern();
252 
253  // build patternData_ vector for implementation of the FieldPattern
254  // functions (indicies will index into fieldIds_)
255  patternData_.resize(dimension+1);
256  int nextIndex = 0;
257  for(int d=0;d<dimension+1;d++) {
258  int numSubcell = geomIdsPattern->getSubcellCount(d);
259  patternData_[d].resize(numSubcell);
260 
261  // loop through sub cells
262  for(int sc=0;sc<numSubcell;sc++) {
263  // a single geometric entity is assigned to field pattern
264  // geomIds size should be either 0 or 1.
265  const std::vector<int> & geomIds = geomIdsPattern->getSubcellIndices(d,sc);
266  TEUCHOS_ASSERT(geomIds.size()<=1);
267  if (geomIds.size() > 0) {
268  const int geomId = geomIds[0];
269  const int numFields = numFields_[geomId];
270  for(int k=0;k<numFields;k++,nextIndex++)
271  patternData_[d][sc].push_back(nextIndex);
272  }
273  }
274  }
275 }
276 
277 const std::vector<int> & FieldAggPattern::localOffsets(int fieldId) const
278 {
279  // lazy evaluation
280  std::map<int,std::vector<int> >::const_iterator itr = fieldOffsets_.find(fieldId);
281  if(itr!=fieldOffsets_.end())
282  return itr->second;
283 
284  std::vector<int> & offsets = fieldOffsets_[fieldId];
285  localOffsets_build(fieldId,offsets);
286  return offsets;
287 }
288 
289 const Kokkos::View<const int*,PHX::Device> FieldAggPattern::localOffsetsKokkos(int fieldId) const
290 {
291  // lazy evaluation
292  std::map<int,Kokkos::View<int*,PHX::Device> >::const_iterator itr = fieldOffsetsKokkos_.find(fieldId);
293  if(itr!=fieldOffsetsKokkos_.end())
294  return itr->second;
295 
296  const auto hostOffsetsStdVector = this->localOffsets(fieldId);
297  Kokkos::View<int*,PHX::Device> offsets("panzer::FieldAggPattern::localOffsetsKokkos",hostOffsetsStdVector.size());
298  auto hostOffsets = Kokkos::create_mirror_view(offsets);
299  for (size_t i=0; i < hostOffsetsStdVector.size(); ++i)
300  hostOffsets(i) = hostOffsetsStdVector[i];
301  Kokkos::deep_copy(offsets,hostOffsets);
302  typename PHX::Device().fence();
303 
304  fieldOffsetsKokkos_[fieldId] = offsets;
305  return offsets;
306 }
307 
309 {
310  if(a[0] < b[0]) return true;
311  if(a[0] > b[0]) return false;
312 
313  // a[0]==b[0]
314  if(a[1] < b[1]) return true;
315  if(a[1] > b[1]) return false;
316 
317  // a[1]==b[1] && a[0]==b[0]
318  if(a[2] < b[2]) return true;
319  if(a[2] > b[2]) return false;
320 
321  // a[2]==b[2] && a[1]==b[1] && a[0]==b[0]
322  return false; // these are equal to, but not less than!
323 }
324 
325 //const std::vector<int> &
326 const std::pair<std::vector<int>,std::vector<int> > &
327 FieldAggPattern::localOffsets_closure(int fieldId,int subcellDim,int subcellId) const
328 {
329  // lazy evaluation
330  typedef std::map<Teuchos::Tuple<int,3>, std::pair<std::vector<int>,std::vector<int> >,LessThan> OffsetMap;
331 
332  Teuchos::Tuple<int,3> subcellTuple = Teuchos::tuple<int>(fieldId,subcellDim,subcellId);
333 
334  OffsetMap::const_iterator itr
335  = fieldSubcellOffsets_closure_.find(subcellTuple);
336  if(itr!=fieldSubcellOffsets_closure_.end()) {
337  return itr->second;
338  }
339 
340  TEUCHOS_TEST_FOR_EXCEPTION(subcellDim >= getDimension(),std::logic_error,
341  "FieldAggPattern::localOffsets_closure precondition subcellDim<getDimension() failed");
342  TEUCHOS_TEST_FOR_EXCEPTION(subcellId < 0,std::logic_error,
343  "FieldAggPattern::localOffsets_closure precondition subcellId>=0 failed");
344  TEUCHOS_TEST_FOR_EXCEPTION(subcellId>=getSubcellCount(subcellDim),std::logic_error,
345  "FieldAggPattern::localOffsets_closure precondition subcellId<getSubcellCount(subcellDim) failed");
346 
347  // build vector for sub cell closure indices
349  const std::vector<int> & fieldOffsets = localOffsets(fieldId);
350 
351  // get offsets into field offsets for the closure indices (from field pattern)
352  std::vector<int> closureOffsets;
353  FPPtr fieldPattern = getFieldPattern(fieldId);
354  fieldPattern->getSubcellClosureIndices(subcellDim,subcellId,closureOffsets);
355 
356  // build closure indices into the correct location in lazy evaluation map.
357  std::pair<std::vector<int>,std::vector<int> > & indicesPair
358  = fieldSubcellOffsets_closure_[subcellTuple];
359 
360  std::vector<int> & closureIndices = indicesPair.first;
361  for(std::size_t i=0;i<closureOffsets.size();i++)
362  closureIndices.push_back(fieldOffsets[closureOffsets[i]]);
363 
364  std::vector<int> & basisIndices = indicesPair.second;
365  basisIndices.assign(closureOffsets.begin(),closureOffsets.end());
366 
367  TEUCHOS_ASSERT(fieldSubcellOffsets_closure_[subcellTuple].first.size()==closureIndices.size());
368  return fieldSubcellOffsets_closure_[subcellTuple];
369 }
370 
371 void FieldAggPattern::localOffsets_build(int fieldId,std::vector<int> & offsets) const
372 {
373  // This function makes the assumption that if there are multiple indices
374  // shared by a subcell then the GeometricAggFieldPattern reflects this.
375  // This is a fine assumption on edges and faces because the symmetries require
376  // extra information about ordering. However, on nodes and "volumes" the
377  // assumption appears to be stupid. For consistency we will make it.
378  //
379  // This code needs to be tested with a basis that has multple IDs at
380  // a node or "volume" sub cell.
381 
382  FPPtr fieldPattern = getFieldPattern(fieldId);
383  FieldType fieldType = getFieldType(fieldId);
384 
385  offsets.clear();
386  offsets.resize(fieldPattern->numberIds(),-111111); // fill with negative number for testing
387 
388  // this will offsets for all IDs associated with the field
389  // but using a geometric ordering
390  std::vector<int> fieldIdsGeomOrder;
391  for(std::size_t i=0;i<fieldIds_.size();++i) {
392  if(fieldIds_[i]==fieldId)
393  fieldIdsGeomOrder.push_back(i);
394  }
395 
396  // Check that number of DOFs line up
397  TEUCHOS_ASSERT((int) fieldIdsGeomOrder.size()==fieldPattern->numberIds());
398 
399  // Build geometric ordering for this pattern: will index into fieldIdsGeomOrder
400  GeometricAggFieldPattern geomPattern(fieldType,fieldPattern);
401  int cnt = 0;
402  for(int dim=0;dim<geomPattern.getDimension()+1;dim++) {
403  for(int sc=0;sc<geomPattern.getSubcellCount(dim);sc++) {
404  const std::vector<int> & fIndices = fieldPattern->getSubcellIndices(dim,sc);
405 
406  for(std::size_t i=0;i<fIndices.size();i++)
407  offsets[fIndices[i]] = fieldIdsGeomOrder[cnt++];
408  }
409  }
410 
411  // check for failure/bug
412  for(std::size_t i=0;i<offsets.size();i++) {
413  TEUCHOS_ASSERT(offsets[i]>=0);
414  }
415 }
416 
417 }
const std::pair< std::vector< int >, std::vector< int > > & localOffsets_closure(int fieldId, int subcellDim, int subcellId) const
std::map< int, std::vector< int > > fieldOffsets_
Stores the Field offsets for the fieldId key. Note that the key is the fieldId, not the index into th...
std::map< int, Kokkos::View< int *, PHX::Device > > fieldOffsetsKokkos_
Stores the Field offsets for the fieldId key. Note that the key is the fieldId, not the index into th...
std::map< Teuchos::Tuple< int, 3 >, std::pair< std::vector< int >, std::vector< int > >, LessThan > fieldSubcellOffsets_closure_
virtual Teuchos::RCP< const FieldPattern > getFieldPattern(int fieldId) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual Teuchos::RCP< const FieldPattern > getGeometricAggFieldPattern() const
std::map< int, int > fieldIdToPatternIdx_
virtual int getSubcellCount(int dimension) const
std::vector< std::vector< std::vector< int > > > patternData_
FieldType
The type of discretization to use for a field pattern.
const Kokkos::View< const int *, PHX::Device > localOffsetsKokkos(int fieldId) const
virtual FieldType getFieldType(int fieldId) const
virtual int getDimension() const
Teuchos::RCP< const FieldPattern > geomAggPattern_
int numFields
virtual void getSubcellClosureIndices(int, int, std::vector< int > &) const
virtual void buildPattern(const std::vector< std::tuple< int, panzer::FieldType, Teuchos::RCP< const FieldPattern > > > &patterns, const Teuchos::RCP< const FieldPattern > &geomAggPattern=Teuchos::null)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void localOffsets_build(int fieldId, std::vector< int > &offsets) const
Continuous Galerkin Formulation.
virtual void print(std::ostream &os) const
virtual void print(std::ostream &os) const
Print this pattern.
std::vector< std::tuple< int, panzer::FieldType, Teuchos::RCP< const FieldPattern > > > patterns_
bool operator()(const Teuchos::Tuple< int, 3 > &a, const Teuchos::Tuple< int, 3 > &b) const
virtual const std::vector< int > & getSubcellIndices(int dimension, int subcell) const
#define TEUCHOS_ASSERT(assertion_test)
virtual shards::CellTopology getCellTopology() const
Kokkos::View< const int *, PHX::Device > offsets
void mergeFieldPatterns(const FieldType &fieldType)
const std::vector< int > & localOffsets(int fieldId) const