Panzer  Version of the Day
Panzer_WorksetContainer.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 
48 #include "Panzer_Dimension.hpp"
49 
50 namespace panzer {
51 
54  : worksetSize_(1)
55 {}
56 
57 WorksetContainer::WorksetContainer(const Teuchos::RCP<const WorksetFactoryBase> & factory,
58  const std::vector<Teuchos::RCP<PhysicsBlock> > & physicsBlocks,
59  std::size_t wkstSz)
60  : wkstFactory_(factory), worksetSize_(wkstSz)
61 {
62  setPhysicsBlockVector(physicsBlocks);
63 }
64 
65 WorksetContainer::WorksetContainer(const Teuchos::RCP<const WorksetFactoryBase> & factory,
66  const std::map<std::string,WorksetNeeds> & needs)
67  : wkstFactory_(factory), worksetSize_(-1)
68 {
69  // thats all!
70  ebToNeeds_ = needs;
71 }
72 
77  : wkstFactory_(wc.wkstFactory_)
78  , worksetSize_(wc.worksetSize_)
79 {
80 }
81 
82 void WorksetContainer::setPhysicsBlockVector(const std::vector<Teuchos::RCP<PhysicsBlock> > & physicsBlocks)
83 {
84  using Teuchos::RCP;
85 
86  for(std::size_t i=0;i<physicsBlocks.size();i++) {
87  ebToPb_[physicsBlocks[i]->elementBlockID()] = physicsBlocks[i];
88  }
89 
90  for(std::size_t i=0;i<physicsBlocks.size();i++) {
91  WorksetNeeds & needs = ebToNeeds_[physicsBlocks[i]->elementBlockID()];
92 
93  needs.cellData = physicsBlocks[i]->cellData();
94 
95  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules = physicsBlocks[i]->getIntegrationRules();
96  for(std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules.begin();
97  ir_itr != int_rules.end(); ++ir_itr)
98  needs.int_rules.push_back(ir_itr->second);
99 
100  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases = physicsBlocks[i]->getBases();
101  const std::vector<StrPureBasisPair>& fieldToBasis = physicsBlocks[i]->getProvidedDOFs();
102  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases.begin();
103  b_itr != bases.end(); ++b_itr) {
104  needs.bases.push_back(b_itr->second);
105 
106  bool found = false;
107  for(std::size_t d=0;d<fieldToBasis.size();d++) {
108  if(fieldToBasis[d].second->name()==b_itr->second->name()) {
109  // add representative basis for this field
110  needs.rep_field_name.push_back(fieldToBasis[d].first);
111  found = true;
112 
113  break;
114  }
115  }
116 
117  // this should always work if physics blocks are correctly constructed
118  TEUCHOS_ASSERT(found);
119  }
120  }
121 }
122 
127 {
128  volWorksets_.clear();
129  sideWorksets_.clear();
130 }
131 
133 const PhysicsBlock & WorksetContainer::lookupPhysicsBlock(const std::string & eBlock) const
134 {
135  std::map<std::string,Teuchos::RCP<PhysicsBlock> >::const_iterator itr = ebToPb_.find(eBlock);
136 
137  TEUCHOS_TEST_FOR_EXCEPTION(itr==ebToPb_.end(),std::logic_error,
138  "WorksetContainer::lookupPhysicsBlock no PhysicsBlock object is associated "
139  "with the element block \""+eBlock+"\".");
140 
141  return *itr->second;
142 }
143 
145 const WorksetNeeds & WorksetContainer::lookupNeeds(const std::string & eBlock) const
146 {
147  std::map<std::string,WorksetNeeds>::const_iterator itr = ebToNeeds_.find(eBlock);
148 
149  TEUCHOS_TEST_FOR_EXCEPTION(itr==ebToNeeds_.end(),std::logic_error,
150  "WorksetContainer::lookupNeeds no WorksetNeeds object is associated "
151  "with the element block \""+eBlock+"\".");
152 
153  return itr->second;
154 }
155 
157 Teuchos::RCP<std::vector<Workset> >
158 WorksetContainer::getVolumeWorksets(const std::string & eBlock)
159 {
160  const WorksetDescriptor wd = blockDescriptor(eBlock);
161 
162  return getWorksets(wd);
163 }
164 
165 Teuchos::RCP<std::vector<Workset> >
167 {
168  Teuchos::RCP<std::vector<Workset> > worksetVector;
169  VolumeMap::iterator itr = volWorksets_.find(wd);
170  if(itr==volWorksets_.end()) {
171  // couldn't find workset, build it!
172  const WorksetNeeds & needs = lookupNeeds(wd.getElementBlock());
173  worksetVector = wkstFactory_->getWorksets(wd,needs);
174 
175  // apply orientations to the just constructed worksets
176  if(worksetVector!=Teuchos::null)
177  applyOrientations(wd.getElementBlock(),*worksetVector);
178 
179  // store vector for reuse in the future
180  volWorksets_[wd] = worksetVector;
181  }
182  else
183  worksetVector = itr->second;
184 
185  return worksetVector;
186 }
187 
189 Teuchos::RCP<std::map<unsigned,Workset> >
191 {
192  Teuchos::RCP<std::map<unsigned,Workset> > worksetMap;
193  SideId side(bc);
194  SideMap::iterator itr = sideWorksets_.find(side);
195  if(itr==sideWorksets_.end()) {
196  // couldn't find workset, build it!
197  if (bc.bcType() == BCT_Interface)
198  worksetMap = wkstFactory_->getSideWorksets(bc, lookupPhysicsBlock(bc.elementBlockID()),
200  else {
201  const std::string & eBlock = side.eblk_id;
202  //const PhysicsBlock & pb = lookupPhysicsBlock(eBlock);
203  const WorksetNeeds & needs = lookupNeeds(eBlock);
204  worksetMap = wkstFactory_->getSideWorksets(bc,needs);
205  }
206 
207  // apply orientations to the worksets for this side
208  if(worksetMap!=Teuchos::null)
209  applyOrientations(side,*worksetMap);
210 
211  // store map for reuse in the future
212  sideWorksets_[side] = worksetMap;
213  }
214  else {
215  worksetMap = itr->second;
216  }
217 
218  return worksetMap;
219 }
220 
221 void WorksetContainer::allocateVolumeWorksets(const std::vector<std::string> & eBlocks)
222 {
223  for(std::size_t i=0;i<eBlocks.size();i++) {
224  // couldn't find workset, build it!
225  const std::string & eBlock = eBlocks[i];
226  const WorksetNeeds & needs = lookupNeeds(eBlock);
227 
228  // store vector for reuse in the future
229  const WorksetDescriptor wd = blockDescriptor(eBlock);
230  volWorksets_[eBlock] = wkstFactory_->getWorksets(wd,needs);
231 
232  // apply orientations to the worksets for this side
233  if(volWorksets_[eBlock]!=Teuchos::null)
234  applyOrientations(eBlock,*volWorksets_[eBlock]);
235  }
236 }
237 
238 void WorksetContainer::allocateSideWorksets(const std::vector<BC> & bcs)
239 {
240  for(std::size_t i=0;i<bcs.size();i++) {
241  // couldn't find workset, build it!
242  const BC & bc = bcs[i];
243  SideId side(bc);
244  const std::string & eBlock = bc.elementBlockID();
245  const WorksetNeeds & needs = lookupNeeds(eBlock);
246 
247  // store map for reuse in the future
248  sideWorksets_[side] = wkstFactory_->getSideWorksets(bc,needs);
249 
250  // apply orientations to the worksets for this side
251  if(sideWorksets_[side]!=Teuchos::null)
252  applyOrientations(side,*sideWorksets_[side]);
253  }
254 }
255 
257 setGlobalIndexer(const Teuchos::RCP<const panzer::UniqueGlobalIndexerBase> & ugi)
258 {
259  // apply the orientations for stored worksets
260  applyOrientations(ugi);
261 }
262 
264 addBasis(const std::string & type,int order,const std::string & rep_field)
265 {
266  using Teuchos::RCP;
267  using Teuchos::rcp;
268 
269  for(auto itr=ebToNeeds_.begin();itr!=ebToNeeds_.end();++itr) {
270  WorksetNeeds & needs = itr->second;
271  RCP<PureBasis> basis = rcp(new PureBasis(type,order,needs.cellData));
272 
273  // add in the new basis
274  needs.bases.push_back(basis);
275  needs.rep_field_name.push_back(rep_field);
276  }
277 
278  // clear all arrays, lazy evaluation means it will be rebuilt
279  clear();
280 }
281 
283 applyOrientations(const Teuchos::RCP<const panzer::UniqueGlobalIndexerBase> & ugi)
284 {
285  // this gurantees orientations won't accidently be applied twice.
286  TEUCHOS_ASSERT(globalIndexer_==Teuchos::null);
287 
288  globalIndexer_ = ugi;
289 
290  // loop over volume worksets, apply orientations to each
291  for(VolumeMap::iterator itr=volWorksets_.begin();
292  itr!=volWorksets_.end();++itr) {
293  std::string eBlock = itr->first.getElementBlock();
294 
295  applyOrientations(eBlock,*itr->second);
296  }
297 
298  // loop over side worksets, apply orientations to each
299  for(SideMap::iterator itr=sideWorksets_.begin();
300  itr!=sideWorksets_.end();itr++) {
301  SideId sideId = itr->first;
302 
303  applyOrientations(sideId,*itr->second);
304  }
305 }
306 
308 applyOrientations(const std::string & eBlock,std::vector<Workset> & worksets) const
309 {
310  using Teuchos::RCP;
311 
312  typedef double Scalar; // orientation container scalar type
313  typedef PHX::MDField<Scalar,Cell,BASIS> Array; // orientation container array type
314 
316  // this is for volume worksets //
318 
319  // short circuit if no global indexer exists
320  if(globalIndexer_==Teuchos::null) {
321  Teuchos::FancyOStream fout(Teuchos::rcpFromRef(std::cout));
322  fout.setOutputToRootOnly(0);
323 
324  fout << "Panzer Warning: No global indexer assigned to a workset container. "
325  << "Orientation of the basis for edge basis functions cannot be applied, "
326  << "if those basis functions are used, there will be problems!" << std::endl;
327  return;
328  }
329 
330  // loop over each basis requiring orientations, then apply them
332 
333  // Note: It may be faster to loop over the basis pairs on the inside (not really sure)
334 
335  const WorksetNeeds & needs = lookupNeeds(eBlock);
336  TEUCHOS_ASSERT(needs.bases.size()==needs.rep_field_name.size());
337 
338  for(std::size_t i=0;i<needs.bases.size();i++) {
339  const std::string & fieldName = needs.rep_field_name[i];
340  const PureBasis & basis = *needs.bases[i];
341 
342  // no need for this if orientations are not required!
343  if(!basis.requiresOrientations()) continue;
344 
345  // build accessors for orientation fields
346  RCP<const OrientationContainerBase<Scalar,Array> > orientationContainer
347  = buildOrientationContainer<Scalar,Array>(globalIndexer_,fieldName);
348 
349  MDFieldArrayFactory fc_factory("",true);
350 
351  // loop over worksets compute and apply orientations
352  for(std::size_t i=0;i<worksets.size();i++) {
353  for(std::size_t j=0;j<worksets[i].numDetails();j++) {
354 
355  // break out of the workset loop
356  if(worksets[i].num_cells<=0) continue;
357 
358  // int array0_sz = worksets[i].num_cells;
359  // int array0_sz = getWorksetSize();
360  int array0_sz = Teuchos::as<int>(needs.cellData.numCells());
361  int array1_sz = basis.functional->dimension(1);
362  Array orientations = fc_factory.buildStaticArray<double,panzer::Cell,panzer::BASIS>("orientations",array0_sz,array1_sz);
363 
364  WorksetDetails & details = worksets[i](j);
365 
366  // compute orientations using the orientation container (and global indexer eventually)
367  orientationContainer->getOrientations(eBlock,details.cell_local_ids,orientations);
368 
369  for(std::size_t basis_index=0;basis_index<details.bases.size();basis_index++) {
370  Teuchos::RCP<const BasisIRLayout> layout = details.bases[basis_index]->basis_layout;
371  TEUCHOS_ASSERT(layout!=Teuchos::null);
372  TEUCHOS_ASSERT(layout->getBasis()!=Teuchos::null);
373  if(layout->getBasis()->name()==basis.name()) {
374  // apply orientations for this basis
375  details.bases[basis_index]->applyOrientations(orientations);
376  }
377  }
378  }
379  }
380  }
381 }
382 
384 applyOrientations(const SideId & sideId,std::map<unsigned,Workset> & worksets) const
385 {
386  using Teuchos::RCP;
387 
388  typedef double Scalar; // orientation container scalar type
389  typedef PHX::MDField<Scalar,Cell,BASIS> Array; // orientation container array type
390 
392  // this is for side worksets //
394 
395  // short circuit if no global indexer exists
396  if(globalIndexer_==Teuchos::null) {
397  Teuchos::FancyOStream fout(Teuchos::rcpFromRef(std::cout));
398  fout.setOutputToRootOnly(0);
399 
400  fout << "Panzer Warning: No global indexer assigned to a workset container. "
401  << "Orientation of the basis for edge basis functions cannot be applied, "
402  << "if those basis functions are used, there will be problems!";
403  return;
404  }
405 
406  // loop over each basis requiring orientations, then apply them
408 
409  // Note: It may be faster to loop over the basis pairs on the inside (not really sure)
410 
411  const WorksetNeeds & needs = lookupNeeds(sideId.eblk_id);
412  TEUCHOS_ASSERT(needs.bases.size()==needs.rep_field_name.size());
413 
414  for(std::size_t i=0;i<needs.bases.size();i++) {
415  const std::string & fieldName = needs.rep_field_name[i];
416  const PureBasis & basis = *needs.bases[i];
417 
418  // no need for this if orientations are not required!
419  if(!basis.requiresOrientations()) continue;
420 
421  // build accessors for orientation fields
422  RCP<const OrientationContainerBase<Scalar,Array> > orientationContainer
423  = buildOrientationContainer<Scalar,Array>(globalIndexer_,fieldName);
424 
425  // loop over worksets compute and apply orientations
426  for(std::map<unsigned,Workset>::iterator itr=worksets.begin();
427  itr!=worksets.end();++itr) {
428 
429  // break out of the workset loop
430  if(itr->second.num_cells<=0) continue;
431 
432  int array0_sz = itr->second.num_cells;
433  int array1_sz = basis.functional->dimension(1);
434  // Intrepid2FieldContainerFactory fc_factory;
435  MDFieldArrayFactory fc_factory("",true);
436  Array orientations = fc_factory.buildStaticArray<double,panzer::Cell,panzer::BASIS>("orientations",array0_sz,array1_sz);
437 
438  for(std::size_t j=0;j<itr->second.numDetails();j++) {
439 
440  WorksetDetails & details = itr->second(j);
441 
442  // compute orientations using the orientation container (and global indexer eventually)
443  orientationContainer->getOrientations(sideId.eblk_id,details.cell_local_ids,orientations);
444 
445  for(std::size_t basis_index=0;basis_index<details.bases.size();basis_index++) {
446  Teuchos::RCP<const BasisIRLayout> layout = details.bases[basis_index]->basis_layout;
447  TEUCHOS_ASSERT(layout!=Teuchos::null);
448  TEUCHOS_ASSERT(layout->getBasis()!=Teuchos::null);
449  if(layout->getBasis()->name()==basis.name()) {
450  // apply orientations for this basis
451  details.bases[basis_index]->applyOrientations(orientations);
452  }
453  }
454  }
455  }
456  }
457 }
458 
460  const std::vector<std::string> & elementBlockNames,
461  std::map<std::string,Teuchos::RCP<std::vector<Workset> > > & volumeWksts)
462 {
463  for(std::size_t i=0;i<elementBlockNames.size();i++)
464  volumeWksts[elementBlockNames[i]] = wc.getVolumeWorksets(elementBlockNames[i]);
465 }
466 
468  const std::vector<BC> & bcs,
469  std::map<BC,Teuchos::RCP<std::map<unsigned,Workset> >,LessBC> & sideWksts)
470 {
471  for(std::size_t i=0;i<bcs.size();i++) {
472  Teuchos::RCP<std::map<unsigned,Workset> > wksts = wc.getSideWorksets(bcs[i]);
473  if(wksts!=Teuchos::null)
474  sideWksts[bcs[i]] = wksts;
475  }
476 }
477 
478 }
std::size_t basis_index
Teuchos::RCP< std::map< unsigned, Workset > > getSideWorksets(const BC &bc)
Access, and construction of side worksets.
std::vector< Teuchos::RCP< const PureBasis > > bases
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
const PhysicsBlock & lookupPhysicsBlock(const std::string &eBlock) const
Look up an input physics block, throws an exception if it can not be found.
Object that contains information on the physics and discretization of a block of elements with the SA...
void setGlobalIndexer(const Teuchos::RCP< const panzer::UniqueGlobalIndexerBase > &ugi)
WorksetContainer()
Default contructor, starts with no workset factory objects.
std::string elementBlockID2() const
Returns the second element block id associated with this sideset.
Definition: Panzer_BC.cpp:205
BCType bcType() const
Returns the boundary condition type (Dirichlet or Neumann or Interface).
Definition: Panzer_BC.cpp:184
Teuchos::RCP< const WorksetFactoryBase > wkstFactory_
std::vector< Teuchos::RCP< const IntegrationRule > > int_rules
std::vector< std::string > rep_field_name
Class that provides access to worksets on each element block and side set.
PHX::MDField< ScalarT > vector
void allocateVolumeWorksets(const std::vector< std::string > &eBlocks)
void applyOrientations(const Teuchos::RCP< const panzer::UniqueGlobalIndexerBase > &ugi)
Teuchos::RCP< std::vector< Workset > > getWorksets(const WorksetDescriptor &wd)
Access to volume worksets.
std::string getElementBlock() const
Get element block.
const WorksetNeeds & lookupNeeds(const std::string &eBlock) const
Look up an input physics block, throws an exception if it can not be found.
Teuchos::RCP< const panzer::UniqueGlobalIndexerBase > globalIndexer_
std::map< std::string, Teuchos::RCP< PhysicsBlock > > ebToPb_
How to construct worksets.
std::string elementBlockID() const
Returns the element block id associated with this sideset.
Definition: Panzer_BC.cpp:198
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
VolumeMap volWorksets_
Maps element blocks to input physics block objects.
void setPhysicsBlockVector(const std::vector< Teuchos::RCP< PhysicsBlock > > &physicsBlocks)
The physics block vector.
void addBasis(const std::string &type, int order, const std::string &rep_field)
void allocateSideWorksets(const std::vector< BC > &bcs)
WorksetDescriptor blockDescriptor(const std::string &eBlock)
void getSideWorksetsFromContainer(WorksetContainer &wc, const std::vector< BC > &bcs, std::map< BC, Teuchos::RCP< std::map< unsigned, Workset > >, LessBC > &sideWksts)
Stores input information for a boundary condition.
Definition: Panzer_BC.hpp:80
Description and data layouts associated with a particular basis.
std::size_t numCells() const
Teuchos::RCP< std::vector< Workset > > getVolumeWorksets(const std::string &eBlock)
Access to volume worksets.
std::map< std::string, WorksetNeeds > ebToNeeds_
Maps element blocks to input physics block objects.
void getVolumeWorksetsFromContainer(WorksetContainer &wc, const std::vector< std::string > &elementBlockNames, std::map< std::string, Teuchos::RCP< std::vector< Workset > > > &volumeWksts)