Panzer Version of the Day
Loading...
Searching...
No Matches
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
47using Teuchos::RCP;
48using Teuchos::rcp;
49
50namespace panzer {
51
53{
54}
55
56FieldAggPattern::FieldAggPattern(std::vector<std::tuple<int,panzer::FieldType,Teuchos::RCP<const FieldPattern> > > & patterns,
57 const Teuchos::RCP<const FieldPattern> & geomAggPattern)
58{
59 buildPattern(patterns,geomAggPattern);
60}
61
62Teuchos::RCP<const FieldPattern> FieldAggPattern::getGeometricAggFieldPattern() const
63{
64 return geomAggPattern_;
65}
66
67void 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
103shards::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
112int FieldAggPattern::getSubcellCount(int dimension) const
113{
114 return patternData_[dimension].size();
115}
116
117const std::vector<int> & FieldAggPattern::getSubcellIndices(int dimension, int subcell) const
118{
119 return patternData_[dimension][subcell];
120}
121
122void 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
128void 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
159Teuchos::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 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
277const 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];
286 return offsets;
287}
288
289const PHX::View<const int*> FieldAggPattern::localOffsetsKokkos(int fieldId) const
290{
291 // lazy evaluation
292 std::map<int,PHX::View<int*> >::const_iterator itr = fieldOffsetsKokkos_.find(fieldId);
293 if(itr!=fieldOffsetsKokkos_.end())
294 return itr->second;
295
296 const auto hostOffsetsStdVector = this->localOffsets(fieldId);
297 PHX::View<int*> 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 fieldOffsetsKokkos_[fieldId] = offsets;
303 return offsets;
304}
305
306bool FieldAggPattern::LessThan::operator()(const Teuchos::Tuple<int,3> & a,const Teuchos::Tuple<int,3> & b) const
307{
308 if(a[0] < b[0]) return true;
309 if(a[0] > b[0]) return false;
310
311 // a[0]==b[0]
312 if(a[1] < b[1]) return true;
313 if(a[1] > b[1]) return false;
314
315 // a[1]==b[1] && a[0]==b[0]
316 if(a[2] < b[2]) return true;
317 if(a[2] > b[2]) return false;
318
319 // a[2]==b[2] && a[1]==b[1] && a[0]==b[0]
320 return false; // these are equal to, but not less than!
321}
322
323//const std::vector<int> &
324const std::pair<std::vector<int>,std::vector<int> > &
325FieldAggPattern::localOffsets_closure(int fieldId,int subcellDim,int subcellId) const
326{
327 // lazy evaluation
328 typedef std::map<Teuchos::Tuple<int,3>, std::pair<std::vector<int>,std::vector<int> >,LessThan> OffsetMap;
329
330 Teuchos::Tuple<int,3> subcellTuple = Teuchos::tuple<int>(fieldId,subcellDim,subcellId);
331
332 OffsetMap::const_iterator itr
333 = fieldSubcellOffsets_closure_.find(subcellTuple);
334 if(itr!=fieldSubcellOffsets_closure_.end()) {
335 return itr->second;
336 }
337
338 TEUCHOS_TEST_FOR_EXCEPTION(subcellDim >= getDimension(),std::logic_error,
339 "FieldAggPattern::localOffsets_closure precondition subcellDim<getDimension() failed");
340 TEUCHOS_TEST_FOR_EXCEPTION(subcellId < 0,std::logic_error,
341 "FieldAggPattern::localOffsets_closure precondition subcellId>=0 failed");
342 TEUCHOS_TEST_FOR_EXCEPTION(subcellId>=getSubcellCount(subcellDim),std::logic_error,
343 "FieldAggPattern::localOffsets_closure precondition subcellId<getSubcellCount(subcellDim) failed");
344
345 // build vector for sub cell closure indices
347 const std::vector<int> & fieldOffsets = localOffsets(fieldId);
348
349 // get offsets into field offsets for the closure indices (from field pattern)
350 std::vector<int> closureOffsets;
351 FPPtr fieldPattern = getFieldPattern(fieldId);
352 fieldPattern->getSubcellClosureIndices(subcellDim,subcellId,closureOffsets);
353
354 // build closure indices into the correct location in lazy evaluation map.
355 std::pair<std::vector<int>,std::vector<int> > & indicesPair
356 = fieldSubcellOffsets_closure_[subcellTuple];
357
358 std::vector<int> & closureIndices = indicesPair.first;
359 for(std::size_t i=0;i<closureOffsets.size();i++)
360 closureIndices.push_back(fieldOffsets[closureOffsets[i]]);
361
362 std::vector<int> & basisIndices = indicesPair.second;
363 basisIndices.assign(closureOffsets.begin(),closureOffsets.end());
364
365 TEUCHOS_ASSERT(fieldSubcellOffsets_closure_[subcellTuple].first.size()==closureIndices.size());
366 return fieldSubcellOffsets_closure_[subcellTuple];
367}
368
369void FieldAggPattern::localOffsets_build(int fieldId,std::vector<int> & offsets) const
370{
371 // This function makes the assumption that if there are multiple indices
372 // shared by a subcell then the GeometricAggFieldPattern reflects this.
373 // This is a fine assumption on edges and faces because the symmetries require
374 // extra information about ordering. However, on nodes and "volumes" the
375 // assumption appears to be stupid. For consistency we will make it.
376 //
377 // This code needs to be tested with a basis that has multple IDs at
378 // a node or "volume" sub cell.
379
380 FPPtr fieldPattern = getFieldPattern(fieldId);
381 FieldType fieldType = getFieldType(fieldId);
382
383 offsets.clear();
384 offsets.resize(fieldPattern->numberIds(),-111111); // fill with negative number for testing
385
386 // this will offsets for all IDs associated with the field
387 // but using a geometric ordering
388 std::vector<int> fieldIdsGeomOrder;
389 for(std::size_t i=0;i<fieldIds_.size();++i) {
390 if(fieldIds_[i]==fieldId)
391 fieldIdsGeomOrder.push_back(i);
392 }
393
394 // Check that number of DOFs line up
395 TEUCHOS_ASSERT((int) fieldIdsGeomOrder.size()==fieldPattern->numberIds());
396
397 // Build geometric ordering for this pattern: will index into fieldIdsGeomOrder
398 GeometricAggFieldPattern geomPattern(fieldType,fieldPattern);
399 int cnt = 0;
400 for(int dim=0;dim<geomPattern.getDimension()+1;dim++) {
401 for(int sc=0;sc<geomPattern.getSubcellCount(dim);sc++) {
402 const std::vector<int> & fIndices = fieldPattern->getSubcellIndices(dim,sc);
403
404 for(std::size_t i=0;i<fIndices.size();i++)
405 offsets[fIndices[i]] = fieldIdsGeomOrder[cnt++];
406 }
407 }
408
409 // check for failure/bug
410 for(std::size_t i=0;i<offsets.size();i++) {
411 TEUCHOS_ASSERT(offsets[i]>=0);
412 }
413}
414
415}
int numFields
PHX::View< const int * > offsets
std::map< int, int > fieldIdToPatternIdx_
virtual void buildPattern(const std::vector< std::tuple< int, panzer::FieldType, Teuchos::RCP< const FieldPattern > > > &patterns, const Teuchos::RCP< const FieldPattern > &geomAggPattern=Teuchos::null)
const std::pair< std::vector< int >, std::vector< int > > & localOffsets_closure(int fieldId, int subcellDim, int subcellId) const
virtual void getSubcellClosureIndices(int, int, std::vector< int > &) const
virtual Teuchos::RCP< const FieldPattern > getFieldPattern(int fieldId) const
void mergeFieldPatterns(const FieldType &fieldType)
virtual int getDimension() const
Teuchos::RCP< const FieldPattern > geomAggPattern_
std::map< Teuchos::Tuple< int, 3 >, std::pair< std::vector< int >, std::vector< int > >, LessThan > fieldSubcellOffsets_closure_
const std::vector< int > & localOffsets(int fieldId) const
std::vector< std::tuple< int, panzer::FieldType, Teuchos::RCP< const FieldPattern > > > patterns_
virtual const std::vector< int > & getSubcellIndices(int dimension, int subcell) const
void localOffsets_build(int fieldId, std::vector< int > &offsets) const
virtual void print(std::ostream &os) const
Print this pattern.
Teuchos::RCP< const FieldPattern > FPPtr
std::map< int, PHX::View< int * > > fieldOffsetsKokkos_
Stores the Field offsets for the fieldId key. Note that the key is the fieldId, not the index into th...
std::vector< std::vector< std::vector< int > > > patternData_
virtual Teuchos::RCP< const FieldPattern > getGeometricAggFieldPattern() const
virtual int getSubcellCount(int dimension) const
virtual FieldType getFieldType(int fieldId) 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...
virtual shards::CellTopology getCellTopology() const
const PHX::View< const int * > localOffsetsKokkos(int fieldId) const
virtual void print(std::ostream &os) const
FieldType
The type of discretization to use for a field pattern.
@ DG
Continuous Galerkin Formulation.
bool operator()(const Teuchos::Tuple< int, 3 > &a, const Teuchos::Tuple< int, 3 > &b) const