Panzer  Version of the Day
Panzer_STK_SetupUtilities.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 #include "Panzer_Workset_Builder.hpp"
45 #include "Teuchos_Assert.hpp"
46 
47 #include <stk_mesh/base/Selector.hpp>
48 #include <stk_mesh/base/GetEntities.hpp>
49 
50 namespace panzer_stk {
51 Teuchos::RCP<std::vector<panzer::Workset> >
53  const panzer::PhysicsBlock & pb)
54 {
55  using namespace workset_utils;
56 
57  std::vector<std::string> element_blocks;
58 
59  std::vector<std::size_t> local_cell_ids;
60  Kokkos::DynRankView<double,PHX::Device> cell_vertex_coordinates;
61 
62  getIdsAndVertices(mesh, pb.elementBlockID(), local_cell_ids, cell_vertex_coordinates);
63 
64  // only build workset if there are elements to worry about
65  // this may be processor dependent, so an element block
66  // may not have elements and thus no contribution
67  // on this processor
68  return panzer::buildWorksets(pb, local_cell_ids, cell_vertex_coordinates);
69 }
70 
71 Teuchos::RCP<std::vector<panzer::Workset> >
73  const std::string & eBlock,
74  const panzer::WorksetNeeds & needs)
75 {
76  using namespace workset_utils;
77 
78  std::vector<std::string> element_blocks;
79 
80  std::vector<std::size_t> local_cell_ids;
81  Kokkos::DynRankView<double,PHX::Device> cell_vertex_coordinates;
82 
83  getIdsAndVertices(mesh, eBlock, local_cell_ids, cell_vertex_coordinates);
84 
85  // only build workset if there are elements to worry about
86  // this may be processor dependent, so an element block
87  // may not have elements and thus no contribution
88  // on this processor
89  return panzer::buildWorksets(needs, eBlock, local_cell_ids, cell_vertex_coordinates);
90 }
91 
92 Teuchos::RCP<std::vector<panzer::Workset> >
94  const panzer::PhysicsBlock & pb,
95  const std::string & sideset,
96  bool useCascade)
97 {
98  using namespace workset_utils;
99  using Teuchos::RCP;
100 
101  std::vector<stk::mesh::Entity> sideEntities;
102 
103  try {
104  // grab local entities on this side
105  // ...catch any failure...primarily wrong side set and element block info
106  if(!useCascade)
107  mesh.getMySides(sideset,pb.elementBlockID(),sideEntities);
108  else
109  mesh.getAllSides(sideset,pb.elementBlockID(),sideEntities);
110  }
112  std::stringstream ss;
113  std::vector<std::string> sideSets;
114  mesh.getSidesetNames(sideSets);
115 
116  // build an error message
117  ss << e.what() << "\nChoose one of:\n";
118  for(std::size_t i=0;i<sideSets.size();i++)
119  ss << "\"" << sideSets[i] << "\"\n";
120 
121  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
122  }
124  std::stringstream ss;
125  std::vector<std::string> elementBlocks;
126  mesh.getElementBlockNames(elementBlocks);
127 
128  // build an error message
129  ss << e.what() << "\nChoose one of:\n";
130  for(std::size_t i=0;i<elementBlocks.size();i++)
131  ss << "\"" << elementBlocks[i] << "\"\n";
132 
133  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
134  }
135  catch(std::logic_error & e) {
136  std::stringstream ss;
137  ss << e.what() << "\nUnrecognized logic error.\n";
138 
139  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
140  }
141 
142  std::vector<stk::mesh::Entity> elements;
143  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> > local_cell_ids;
144  if(!useCascade) {
145  unsigned subcell_dim = pb.cellData().baseCellDimension()-1;
146  std::vector<std::size_t> local_side_ids;
147  getSideElements(mesh, pb.elementBlockID(),
148  sideEntities,local_side_ids,elements);
149 
150  // build local cell_ids, mapped by local side id
151  for(std::size_t elm=0;elm<elements.size();++elm) {
152  stk::mesh::Entity element = elements[elm];
153 
154  local_cell_ids[std::make_pair(subcell_dim,local_side_ids[elm])].push_back(mesh.elementLocalId(element));
155  }
156  }
157  else {
158  std::vector<std::size_t> local_subcell_ids, subcell_dim;
160  sideEntities,subcell_dim,local_subcell_ids,elements);
161 
162  // build local cell_ids, mapped by local side id
163  for(std::size_t elm=0;elm<elements.size();++elm) {
164  stk::mesh::Entity element = elements[elm];
165 
166  local_cell_ids[std::make_pair(subcell_dim[elm],local_subcell_ids[elm])].push_back(mesh.elementLocalId(element));
167  }
168  }
169 
170  // only build workset if there are elements to worry about
171  // this may be processor dependent, so a defined boundary
172  // condition may have not elements and thus no contribution
173  // on this processor
174  if(elements.size()!=0) {
175  Teuchos::RCP<const shards::CellTopology> topo = mesh.getCellTopology(pb.elementBlockID());
176 
177  // worksets to be returned
178  Teuchos::RCP<std::vector<panzer::Workset> > worksets = Teuchos::rcp(new std::vector<panzer::Workset>);
179 
180  // loop over each side
181  for(std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> >::const_iterator itr=local_cell_ids.begin();
182  itr!=local_cell_ids.end();++itr) {
183 
184  if(itr->second.size()==0)
185  continue;
186 
187  Kokkos::DynRankView<double,PHX::Device> vertices;
188  mesh.getElementVertices(itr->second,pb.elementBlockID(),vertices);
189 
190  Teuchos::RCP<std::vector<panzer::Workset> > current
191  = panzer::buildWorksets(pb, itr->second, vertices);
192 
193  // correct worksets so the sides are correct
194  for(std::size_t w=0;w<current->size();w++) {
195  (*current)[w].subcell_dim = itr->first.first;
196  (*current)[w].subcell_index = itr->first.second;
197  }
198 
199  // append new worksets
200  worksets->insert(worksets->end(),current->begin(),current->end());
201  }
202 
203  return worksets;
204  }
205 
206  // return Teuchos::null;
207  return Teuchos::rcp(new std::vector<panzer::Workset>());
208 }
209 
210 Teuchos::RCP<std::vector<panzer::Workset> >
212  const panzer::WorksetNeeds & needs,
213  const std::string & sideset,
214  const std::string & eBlock,
215  bool useCascade)
216 {
217  using namespace workset_utils;
218  using Teuchos::RCP;
219 
220  std::vector<stk::mesh::Entity> sideEntities;
221 
222  try {
223  // grab local entities on this side
224  // ...catch any failure...primarily wrong side set and element block info
225  if(!useCascade)
226  mesh.getMySides(sideset,eBlock,sideEntities);
227  else
228  mesh.getAllSides(sideset,eBlock,sideEntities);
229  }
231  std::stringstream ss;
232  std::vector<std::string> sideSets;
233  mesh.getSidesetNames(sideSets);
234 
235  // build an error message
236  ss << e.what() << "\nChoose one of:\n";
237  for(std::size_t i=0;i<sideSets.size();i++)
238  ss << "\"" << sideSets[i] << "\"\n";
239 
240  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
241  }
243  std::stringstream ss;
244  std::vector<std::string> elementBlocks;
245  mesh.getElementBlockNames(elementBlocks);
246 
247  // build an error message
248  ss << e.what() << "\nChoose one of:\n";
249  for(std::size_t i=0;i<elementBlocks.size();i++)
250  ss << "\"" << elementBlocks[i] << "\"\n";
251 
252  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
253  }
254  catch(std::logic_error & e) {
255  std::stringstream ss;
256  ss << e.what() << "\nUnrecognized logic error.\n";
257 
258  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
259  }
260 
261  std::vector<stk::mesh::Entity> elements;
262  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> > local_cell_ids;
263  if(!useCascade) {
264  unsigned subcell_dim = needs.cellData.baseCellDimension()-1;
265  std::vector<std::size_t> local_side_ids;
266  getSideElements(mesh, eBlock,
267  sideEntities,local_side_ids,elements);
268 
269  // build local cell_ids, mapped by local side id
270  for(std::size_t elm=0;elm<elements.size();++elm) {
271  stk::mesh::Entity element = elements[elm];
272 
273  local_cell_ids[std::make_pair(subcell_dim,local_side_ids[elm])].push_back(mesh.elementLocalId(element));
274  }
275  }
276  else {
277  std::vector<std::size_t> local_subcell_ids, subcell_dim;
278  getSideElementCascade(mesh, eBlock,
279  sideEntities,subcell_dim,local_subcell_ids,elements);
280 
281  // build local cell_ids, mapped by local side id
282  for(std::size_t elm=0;elm<elements.size();++elm) {
283  stk::mesh::Entity element = elements[elm];
284 
285  local_cell_ids[std::make_pair(subcell_dim[elm],local_subcell_ids[elm])].push_back(mesh.elementLocalId(element));
286  }
287  }
288 
289  // only build workset if there are elements to worry about
290  // this may be processor dependent, so a defined boundary
291  // condition may have not elements and thus no contribution
292  // on this processor
293  if(elements.size()!=0) {
294  Teuchos::RCP<const shards::CellTopology> topo = mesh.getCellTopology(eBlock);
295 
296  // worksets to be returned
297  Teuchos::RCP<std::vector<panzer::Workset> > worksets = Teuchos::rcp(new std::vector<panzer::Workset>);
298 
299  // loop over each side
300  for(std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> >::const_iterator itr=local_cell_ids.begin();
301  itr!=local_cell_ids.end();++itr) {
302 
303  if(itr->second.size()==0)
304  continue;
305 
306  Kokkos::DynRankView<double,PHX::Device> vertices;
307  mesh.getElementVertices(itr->second,eBlock,vertices);
308 
309  Teuchos::RCP<std::vector<panzer::Workset> > current
310  = panzer::buildWorksets(needs, eBlock, itr->second, vertices);
311 
312  // correct worksets so the sides are correct
313  for(std::size_t w=0;w<current->size();w++) {
314  (*current)[w].subcell_dim = itr->first.first;
315  (*current)[w].subcell_index = itr->first.second;
316  }
317 
318  // append new worksets
319  worksets->insert(worksets->end(),current->begin(),current->end());
320  }
321 
322  return worksets;
323  }
324 
325  // return Teuchos::null;
326  return Teuchos::rcp(new std::vector<panzer::Workset>());
327 }
328 
329 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
331  const panzer::PhysicsBlock & pb_a,
332  const panzer::PhysicsBlock & pb_b,
333  const std::string & sideset)
334 {
335  using namespace workset_utils;
336  using Teuchos::RCP;
337 
338  std::vector<stk::mesh::Entity> sideEntities; // we will reduce a_ and b_ to this vector
339 
340  try {
341  // grab local entities on this side
342  // ...catch any failure...primarily wrong side set and element block info
343 
344  // we can't use getMySides because it only returns locally owned sides
345  // this gurantees all the sides are extracted (element ownership is considered
346  // we we call getSideElements below)
347 
348  stk::mesh::Part * sidePart = mesh.getSideset(sideset);
349  TEUCHOS_TEST_FOR_EXCEPTION(sidePart==0,std::logic_error,
350  "Unknown side set \"" << sideset << "\"");
351 
352  stk::mesh::Selector side = *sidePart;
353  // stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & side;
354 
355  // grab elements
356  stk::mesh::get_selected_entities(side,mesh.getBulkData()->buckets(mesh.getSideRank()),sideEntities);
357  }
359  std::stringstream ss;
360  std::vector<std::string> elementBlocks;
361  mesh.getElementBlockNames(elementBlocks);
362 
363  // build an error message
364  ss << e.what() << "\nChoose one of:\n";
365  for(std::size_t i=0;i<elementBlocks.size();i++)
366  ss << "\"" << elementBlocks[i] << "\"\n";
367 
368  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
369  }
370  catch(std::logic_error & e) {
371  std::stringstream ss;
372  ss << e.what() << "\nUnrecognized logic error.\n";
373 
374  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
375  }
376 
377  std::vector<stk::mesh::Entity> elements_a, elements_b;
378  std::vector<std::size_t> local_cell_ids_a, local_cell_ids_b;
379  std::vector<std::size_t> local_side_ids_a, local_side_ids_b;
380 
381  // this enforces that "a" elements must be owned.
382  getSideElements(mesh, pb_a.elementBlockID(), pb_b.elementBlockID(), sideEntities,
383  local_side_ids_a,elements_a,
384  local_side_ids_b,elements_b);
385 
386  TEUCHOS_TEST_FOR_EXCEPTION(elements_a.size()!=elements_b.size(),std::logic_error,
387  "For a DG type boundary, the number of elements on the \"left\" and \"right\" is not the same.");
388 
389  // only build workset if there are elements to worry about
390  // this may be processor dependent, so a defined boundary
391  // condition may have not elements and thus no contribution
392  // on this processor
393  if(elements_a.size()==0)
394  return Teuchos::rcp(new std::map<unsigned,panzer::Workset>);
395 
396  // loop over elements of this block (note the assures that element_a and element_b
397  // are the same size, the ordering is the same because the order of sideEntities is
398  // the same
399  for(std::size_t elm=0;elm<elements_a.size();++elm) {
400  stk::mesh::Entity element_a = elements_a[elm];
401  stk::mesh::Entity element_b = elements_b[elm];
402 
403  local_cell_ids_a.push_back(mesh.elementLocalId(element_a));
404  local_cell_ids_b.push_back(mesh.elementLocalId(element_b));
405  }
406 
407  Kokkos::DynRankView<double,PHX::Device> vertex_coordinates_a, vertex_coordinates_b;
408  mesh.getElementVertices(local_cell_ids_a,pb_a.elementBlockID(),vertex_coordinates_a);
409  mesh.getElementVertices(local_cell_ids_b,pb_b.elementBlockID(),vertex_coordinates_b);
410 
411  // worksets to be returned
412  return buildBCWorkset(pb_a, local_cell_ids_a, local_side_ids_a, vertex_coordinates_a,
413  pb_b, local_cell_ids_b, local_side_ids_b, vertex_coordinates_b);
414 }
415 
416 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
418  const panzer::PhysicsBlock & pb,
419  const std::string & sidesetID)
420 {
421  using namespace workset_utils;
422  using Teuchos::RCP;
423 
424  std::vector<stk::mesh::Entity> sideEntities;
425 
426  try {
427  // grab local entities on this side
428  // ...catch any failure...primarily wrong side set and element block info
429  mesh.getMySides(sidesetID,pb.elementBlockID(),sideEntities);
430  }
432  std::stringstream ss;
433  std::vector<std::string> sideSets;
434  mesh.getSidesetNames(sideSets);
435 
436  // build an error message
437  ss << e.what() << "\nChoose one of:\n";
438  for(std::size_t i=0;i<sideSets.size();i++)
439  ss << "\"" << sideSets[i] << "\"\n";
440 
441  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
442  }
444  std::stringstream ss;
445  std::vector<std::string> elementBlocks;
446  mesh.getElementBlockNames(elementBlocks);
447 
448  // build an error message
449  ss << e.what() << "\nChoose one of:\n";
450  for(std::size_t i=0;i<elementBlocks.size();i++)
451  ss << "\"" << elementBlocks[i] << "\"\n";
452 
453  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
454  }
455  catch(std::logic_error & e) {
456  std::stringstream ss;
457  ss << e.what() << "\nUnrecognized logic error.\n";
458 
459  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
460  }
461 
462  std::vector<stk::mesh::Entity> elements;
463  std::vector<std::size_t> local_cell_ids;
464  std::vector<std::size_t> local_side_ids;
465  getSideElements(mesh, pb.elementBlockID(),
466  sideEntities,local_side_ids,elements);
467 
468  // loop over elements of this block
469  for(std::size_t elm=0;elm<elements.size();++elm) {
470  stk::mesh::Entity element = elements[elm];
471 
472  local_cell_ids.push_back(mesh.elementLocalId(element));
473  }
474 
475  // only build workset if there are elements to worry about
476  // this may be processor dependent, so a defined boundary
477  // condition may have not elements and thus no contribution
478  // on this processor
479  if(elements.size()!=0) {
480  Teuchos::RCP<const shards::CellTopology> topo
481  = mesh.getCellTopology(pb.elementBlockID());
482 
483  Kokkos::DynRankView<double,PHX::Device> vertices;
484  mesh.getElementVertices(local_cell_ids,pb.elementBlockID(),vertices);
485 
486  return panzer::buildBCWorkset(pb, local_cell_ids, local_side_ids, vertices);
487  }
488 
489  return Teuchos::null;
490 }
491 
492 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
494  const panzer::WorksetNeeds & needs,
495  const std::string & eblockID,
496  const std::string & sidesetID)
497 {
498  using namespace workset_utils;
499  using Teuchos::RCP;
500 
501  std::vector<stk::mesh::Entity> sideEntities;
502 
503  try {
504  // grab local entities on this side
505  // ...catch any failure...primarily wrong side set and element block info
506  mesh.getMySides(sidesetID,eblockID,sideEntities);
507  }
509  std::stringstream ss;
510  std::vector<std::string> sideSets;
511  mesh.getSidesetNames(sideSets);
512 
513  // build an error message
514  ss << e.what() << "\nChoose one of:\n";
515  for(std::size_t i=0;i<sideSets.size();i++)
516  ss << "\"" << sideSets[i] << "\"\n";
517 
518  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
519  }
521  std::stringstream ss;
522  std::vector<std::string> elementBlocks;
523  mesh.getElementBlockNames(elementBlocks);
524 
525  // build an error message
526  ss << e.what() << "\nChoose one of:\n";
527  for(std::size_t i=0;i<elementBlocks.size();i++)
528  ss << "\"" << elementBlocks[i] << "\"\n";
529 
530  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
531  }
532  catch(std::logic_error & e) {
533  std::stringstream ss;
534  ss << e.what() << "\nUnrecognized logic error.\n";
535 
536  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
537  }
538 
539  std::vector<stk::mesh::Entity> elements;
540  std::vector<std::size_t> local_cell_ids;
541  std::vector<std::size_t> local_side_ids;
542  getSideElements(mesh, eblockID,
543  sideEntities,local_side_ids,elements);
544 
545  // loop over elements of this block
546  for(std::size_t elm=0;elm<elements.size();++elm) {
547  stk::mesh::Entity element = elements[elm];
548 
549  local_cell_ids.push_back(mesh.elementLocalId(element));
550  }
551 
552  // only build workset if there are elements to worry about
553  // this may be processor dependent, so a defined boundary
554  // condition may have not elements and thus no contribution
555  // on this processor
556  if(elements.size()!=0) {
557  Teuchos::RCP<const shards::CellTopology> topo
558  = mesh.getCellTopology(eblockID);
559 
560  Kokkos::DynRankView<double,PHX::Device> vertices;
561  mesh.getElementVertices(local_cell_ids,eblockID,vertices);
562 
563  return panzer::buildBCWorkset(needs, eblockID, local_cell_ids, local_side_ids, vertices);
564  }
565 
566  return Teuchos::null;
567 }
568 
569 namespace workset_utils {
570 
572  const std::string & blockId,
573  const std::vector<stk::mesh::Entity> & entities,
574  std::vector<std::size_t> & localEntityIds,
575  std::vector<stk::mesh::Entity> & elements)
576 {
577  // for verifying that an element is in specified block
578  stk::mesh::Part * blockPart = mesh.getElementBlockPart(blockId);
579  stk::mesh::Part * ownedPart = mesh.getOwnedPart();
580  stk::mesh::BulkData& bulkData = *mesh.getBulkData();
581 
582  // loop over each entitiy extracting elements and local entity ID that
583  // are containted in specified block.
584  std::vector<stk::mesh::Entity>::const_iterator entityItr;
585  for(entityItr=entities.begin();entityItr!=entities.end();++entityItr) {
586  stk::mesh::Entity entity = *entityItr;
587 
588  const size_t num_rels = bulkData.num_elements(entity);
589  stk::mesh::Entity const* relations = bulkData.begin_elements(entity);
590  stk::mesh::ConnectivityOrdinal const* ordinals = bulkData.begin_element_ordinals(entity);
591  for(std::size_t e=0; e<num_rels; ++e) {
592  stk::mesh::Entity element = relations[e];
593  std::size_t entityId = ordinals[e];
594 
595  // is this element in requested block
596  stk::mesh::Bucket const& bucket = bulkData.bucket(element);
597  bool inBlock = bucket.member(*blockPart);
598  bool onProc = bucket.member(*ownedPart);
599  if(inBlock && onProc) {
600  // add element and Side ID to output vectors
601  elements.push_back(element);
602  localEntityIds.push_back(entityId);
603  }
604  }
605  }
606 }
607 
609  const std::string & blockId,
610  const std::vector<stk::mesh::Entity> & entities,
611  std::vector<std::size_t> & localEntityIds,
612  std::vector<stk::mesh::Entity> & elements)
613 {
614  // for verifying that an element is in specified block
615  stk::mesh::Part * blockPart = mesh.getElementBlockPart(blockId);
616  stk::mesh::Part * universalPart = &mesh.getMetaData()->universal_part();
617  stk::mesh::BulkData& bulkData = *mesh.getBulkData();
618 
619  // loop over each entitiy extracting elements and local entity ID that
620  // are containted in specified block.
621  std::vector<stk::mesh::Entity>::const_iterator entityItr;
622  for(entityItr=entities.begin();entityItr!=entities.end();++entityItr) {
623  stk::mesh::Entity entity = *entityItr;
624 
625  const size_t num_rels = bulkData.num_elements(entity);
626  stk::mesh::Entity const* element_rels = bulkData.begin_elements(entity);
627  stk::mesh::ConnectivityOrdinal const* ordinals = bulkData.begin_element_ordinals(entity);
628  for(std::size_t e=0; e<num_rels; ++e) {
629  stk::mesh::Entity element = element_rels[e];
630  std::size_t entityId = ordinals[e];
631 
632  // is this element in requested block
633  stk::mesh::Bucket const& bucket = bulkData.bucket(element);
634  bool inBlock = bucket.member(*blockPart);
635  bool onProc = bucket.member(*universalPart);
636  if(inBlock && onProc) {
637  // add element and Side ID to output vectors
638  elements.push_back(element);
639  localEntityIds.push_back(entityId);
640  }
641  }
642  }
643 }
644 
646  const std::string & blockId,
647  const std::vector<stk::mesh::Entity> & sides,
648  std::vector<std::size_t> & localSubcellDim,
649  std::vector<std::size_t> & localSubcellIds,
650  std::vector<stk::mesh::Entity> & elements)
651 {
652  // This is the alogrithm, for computing the side element
653  // cascade. The requirements are that for a particular set of sides
654  // we compute all elements and subcells where they touch the side. Note
655  // that elements can be and will be repeated within this list.
656 
657  std::vector<std::vector<stk::mesh::Entity> > subcells;
658  getSubcellEntities(mesh,sides,subcells);
659  subcells.push_back(sides);
660 
661  // subcells now contains a unique list of faces, edges and nodes that
662  // intersect with the sides
663 
664  for(std::size_t d=0;d<subcells.size();d++) {
665  std::vector<std::size_t> subcellIds;
666  std::vector<stk::mesh::Entity> subcellElements;
667 
668  // find elements connected to the subcells and their local subcell information
669  getSubcellElements(mesh,blockId,subcells[d],subcellIds,subcellElements);
670 
671  // sanity check
672  TEUCHOS_ASSERT(subcellIds.size()==subcellElements.size());
673 
674  // concatenate with found elements
675  localSubcellDim.insert(localSubcellDim.end(),subcellElements.size(),d);
676  localSubcellIds.insert(localSubcellIds.end(),subcellIds.begin(),subcellIds.end());
677  elements.insert(elements.end(),subcellElements.begin(),subcellElements.end());
678  }
679 }
680 
682  const std::string & blockId,
683  const std::vector<stk::mesh::Entity> & sides,
684  std::vector<std::size_t> & localSideIds,
685  std::vector<stk::mesh::Entity> & elements)
686 {
687  getSubcellElements(mesh,blockId,sides,localSideIds,elements);
688 }
689 
691  const std::string & blockId_a,
692  const std::string & blockId_b,
693  const std::vector<stk::mesh::Entity> & sides,
694  std::vector<std::size_t> & localSideIds_a,
695  std::vector<stk::mesh::Entity> & elements_a,
696  std::vector<std::size_t> & localSideIds_b,
697  std::vector<stk::mesh::Entity> & elements_b)
698 {
699  // for verifying that an element is in specified block
700  stk::mesh::Part * blockPart_a = mesh.getElementBlockPart(blockId_a);
701  stk::mesh::Part * blockPart_b = mesh.getElementBlockPart(blockId_b);
702  stk::mesh::Part * ownedPart = mesh.getOwnedPart();
703  stk::mesh::Part * universalPart = &mesh.getMetaData()->universal_part();
704  stk::mesh::BulkData& bulkData = *mesh.getBulkData();
705 
706  // loop over each entitiy extracting elements and local entity ID that
707  // are containted in specified block.
708  std::vector<stk::mesh::Entity>::const_iterator sidesItr;
709  for(sidesItr=sides.begin();sidesItr!=sides.end();++sidesItr) {
710  stk::mesh::Entity side = *sidesItr;
711 
712  // these are used below the loop to insert into the appropriate vectors
713  stk::mesh::Entity element_a = stk::mesh::Entity(), element_b = stk::mesh::Entity();
714  std::size_t entityId_a=0, entityId_b=0;
715 
716  const size_t num_rels = bulkData.num_elements(side);
717  stk::mesh::Entity const* element_rels = bulkData.begin_elements(side);
718  stk::mesh::ConnectivityOrdinal const* ordinals = bulkData.begin_element_ordinals(side);
719  for(std::size_t e=0; e<num_rels; ++e) {
720  stk::mesh::Entity element = element_rels[e];
721  std::size_t entityId = ordinals[e];
722 
723  // is this element in requested block
724  stk::mesh::Bucket const& bucket = bulkData.bucket(element);
725  bool inBlock_a = bucket.member(*blockPart_a);
726  bool inBlock_b = bucket.member(*blockPart_b);
727  bool onProc = bucket.member(*ownedPart);
728  bool unProc = bucket.member(*universalPart);
729 
730  if(inBlock_a && onProc) {
731  TEUCHOS_ASSERT(element_a==stk::mesh::Entity()); // sanity check
732  element_a = element;
733  entityId_a = entityId;
734  }
735  if(inBlock_b && unProc) {
736  TEUCHOS_ASSERT(element_b==stk::mesh::Entity()); // sanity check
737  element_b = element;
738  entityId_b = entityId;
739  }
740  }
741 
742  if(element_a!=stk::mesh::Entity() && element_b!=stk::mesh::Entity()) { // add element and Side ID to output vectors
743  elements_a.push_back(element_a);
744  localSideIds_a.push_back(entityId_a);
745 
746  // add element and Side ID to output vectors
747  elements_b.push_back(element_b);
748  localSideIds_b.push_back(entityId_b);
749  }
750  }
751 }
752 
754  const std::string & blockId,
755  const std::vector<stk::mesh::Entity> & nodes,
756  std::vector<std::size_t> & localNodeIds,
757  std::vector<stk::mesh::Entity> & elements)
758 {
759  getSubcellElements(mesh,blockId,nodes,localNodeIds,elements);
760 }
761 
763  const std::vector<stk::mesh::Entity> & entities,
764  std::vector<std::vector<stk::mesh::Entity> > & subcells)
765 {
766  // exit if there is no work to do
767  if(entities.size()==0) {
768  subcells.clear();
769  return;
770  }
771 
772  stk::mesh::BulkData& bulkData = *mesh.getBulkData();
773  stk::mesh::EntityRank master_rank = bulkData.entity_rank(entities[0]);
774 
775  std::vector<std::set<stk::mesh::Entity> > subcells_set(master_rank);
776 
777  // loop over each entitiy extracting elements and local entity ID that
778  // are containted in specified block.
779  std::vector<stk::mesh::Entity>::const_iterator entityItr;
780  for(entityItr=entities.begin();entityItr!=entities.end();++entityItr) {
781  stk::mesh::Entity entity = *entityItr;
782 
783  // sanity check, enforcing that there is only one rank
784  TEUCHOS_ASSERT(bulkData.entity_rank(entity)==master_rank);
785 
786  for(int i=0; i<master_rank; i++) {
787  stk::mesh::EntityRank const to_rank = static_cast<stk::mesh::EntityRank>(i);
788  if (bulkData.connectivity_map().valid(master_rank, to_rank)) {
789  const size_t num_rels = bulkData.num_connectivity(entity, to_rank);
790  stk::mesh::Entity const* relations = bulkData.begin(entity, to_rank);
791 
792  // for each relation insert the appropriate entity (into the set
793  // which gurantees uniqueness
794  for(std::size_t e=0; e<num_rels; ++e) {
795  stk::mesh::Entity subcell = relations[e];
796 
797  subcells_set[i].insert(subcell);
798  }
799  }
800  }
801  }
802 
803  // copy unique entities into vector
804  subcells.clear();
805  subcells.resize(subcells_set.size());
806  for(std::size_t i=0;i<subcells_set.size();i++)
807  subcells[i].assign(subcells_set[i].begin(),subcells_set[i].end());
808 }
809 
810 }
811 
812 }
Teuchos::RCP< std::map< unsigned, Workset > > buildBCWorkset(const PhysicsBlock &volume_pb, const std::vector< std::size_t > &local_cell_ids, const std::vector< std::size_t > &local_side_ids, const ArrayT &vertex_coordinates)
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block count
void getElementVertices(const std::vector< std::size_t > &localIds, ArrayT &vertices) const
void getMySides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void getIdsAndVertices(const panzer_stk::STK_Interface &mesh, std::string blockId, std::vector< std::size_t > &localIds, ArrayT &vertices)
Object that contains information on the physics and discretization of a block of elements with the SA...
std::size_t elementLocalId(stk::mesh::Entity elmt) const
void getSubcellEntities(const panzer_stk::STK_Interface &mesh, const std::vector< stk::mesh::Entity > &entities, std::vector< std::vector< stk::mesh::Entity > > &subcells)
std::string elementBlockID() const
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
Teuchos::RCP< std::vector< panzer::Workset > > buildWorksets(const panzer_stk::STK_Interface &mesh, const panzer::PhysicsBlock &pb)
void getSubcellElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &entities, std::vector< std::size_t > &localEntityIds, std::vector< stk::mesh::Entity > &elements)
stk::mesh::Part * getOwnedPart() const
Get a pointer to the locally owned part.
stk::mesh::EntityRank getSideRank() const
PHX::MDField< ScalarT > vector
void getElementBlockNames(std::vector< std::string > &names) const
stk::mesh::Part * getSideset(const std::string &name) const
int baseCellDimension() const
Dimension of the base cell. NOT the dimension of the local side, even if the side() method returns tr...
void getSideElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &sides, std::vector< std::size_t > &localSideIds, std::vector< stk::mesh::Entity > &elements)
void getUniversalSubcellElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &entities, std::vector< std::size_t > &localEntityIds, std::vector< stk::mesh::Entity > &elements)
Teuchos::RCP< std::vector< Workset > > buildWorksets(const PhysicsBlock &pb, const std::vector< std::size_t > &local_cell_ids, const ArrayT &vertex_coordinates)
const panzer::CellData & cellData() const
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const
void getAllSides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void getSidesetNames(std::vector< std::string > &name) const
Teuchos::RCP< std::map< unsigned, panzer::Workset > > buildBCWorksets(const panzer_stk::STK_Interface &mesh, const panzer::PhysicsBlock &pb_a, const panzer::PhysicsBlock &pb_b, const std::string &sideset)
void getSideElementCascade(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &sides, std::vector< std::size_t > &localSubcellDim, std::vector< std::size_t > &localSubcellIds, std::vector< stk::mesh::Entity > &elements)
void getNodeElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &nodes, std::vector< std::size_t > &localNodeIds, std::vector< stk::mesh::Entity > &elements)