Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_STK_PeriodicBC_Matcher.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Panzer: A partial differential equation assembly
4 // engine for strongly coupled complex multiphysics systems
5 //
6 // Copyright 2011 NTESS and the Panzer contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
12 
13 #include <stk_mesh/base/GetEntities.hpp>
14 #include <stk_mesh/base/FieldBase.hpp>
15 
16 #include "Panzer_NodeType.hpp"
17 #include "Tpetra_Map.hpp"
18 #include "Tpetra_Import.hpp"
19 #include "Tpetra_Vector.hpp"
20 
21 #include "Teuchos_FancyOStream.hpp"
22 
23 namespace panzer_stk {
24 
25 namespace periodic_helpers {
26 
28 getGlobalPairing(const std::vector<std::size_t> & locallyRequiredIds,
29  const std::vector<std::pair<std::size_t,std::size_t> > & locallyMatchedIds,
30  const STK_Interface & mesh,bool failure)
31 {
32  using LO = panzer::LocalOrdinal;
33  using GO = panzer::GlobalOrdinal;
35  using Map = Tpetra::Map<LO,GO,NODE>;
36  using Importer = Tpetra::Import<LO,GO,NODE>;
37 
38  auto comm = mesh.getComm();
39 
40  // this is needed to prevent hanging: it is unfortunately expensive
41  // need a better way!
42  int myVal = failure ? 1 : 0;
43  int sumVal = 0;
44  Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_SUM,1,&myVal,&sumVal);
45  TEUCHOS_ASSERT(sumVal==0);
46 
47  std::vector<GO> requiredInts(locallyRequiredIds.size());
48  for(std::size_t i=0;i<requiredInts.size();i++)
49  requiredInts[i] = locallyRequiredIds[i];
50 
51  std::vector<GO> providedInts(locallyMatchedIds.size());
52  for(std::size_t i=0;i<locallyMatchedIds.size();i++)
53  providedInts[i] = locallyMatchedIds[i].first;
54 
55  // maps and communciation all set up
57  Teuchos::RCP<Map> requiredMap = Teuchos::rcp(new Map(computeInternally,Teuchos::ArrayView<const GO>(requiredInts),0,comm));
58  Teuchos::RCP<Map> providedMap = Teuchos::rcp(new Map(computeInternally,Teuchos::ArrayView<const GO>(providedInts),0,comm));
59  Importer importer(providedMap,requiredMap);
60 
61  // this is what to distribute
62  Tpetra::Vector<GO,LO,GO,NODE> providedVector(providedMap);
63  {
64  auto pvHost = providedVector.getLocalViewHost(Tpetra::Access::OverwriteAll);
65  for(std::size_t i=0;i<locallyMatchedIds.size();i++)
66  pvHost(i,0) = locallyMatchedIds[i].second;
67  }
68 
69  Tpetra::Vector<GO,LO,GO,NODE> requiredVector(requiredMap);
70  requiredVector.doImport(providedVector,importer,Tpetra::INSERT);
71 
72 
74  = Teuchos::rcp(new std::vector<std::pair<std::size_t,std::size_t> >(requiredInts.size()));
75 
76  auto rvHost = requiredVector.getLocalViewHost(Tpetra::Access::ReadOnly);
77  for(std::size_t i=0;i<result->size();i++) {
78  (*result)[i].first = requiredInts[i];
79  (*result)[i].second = rvHost(i,0);
80  }
81  return result;
82 }
83 
84 
85 
91  const std::string & sideName, const std::string type_)
92 {
95 
96  // grab nodes owned by requested side
98  std::stringstream ss;
99  ss << "Can't find part=\"" << sideName << "\"" << std::endl;
100  stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
101  stk::mesh::Selector mySides = *side & (metaData->locally_owned_part() | metaData->globally_shared_part());
102 
103  stk::mesh::EntityRank rank;
104  panzer::GlobalOrdinal offset = 0; // offset to avoid giving nodes, edges, faces the same sideId
105  if(type_ == "coord"){
106  rank = mesh.getNodeRank();
107  } else if(type_ == "edge"){
108  rank = mesh.getEdgeRank();
109  offset = mesh.getMaxEntityId(mesh.getNodeRank());
110  } else if(type_ == "face"){
111  rank = mesh.getFaceRank();
112  offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
113  } else {
114  ss << "Can't do BCs of type " << type_ << std::endl;
115  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
116  }
117 
118  std::vector<stk::mesh::Bucket*> const& nodeBuckets =
119  bulkData->get_buckets(rank, mySides);
120 
121  // build id vector
123  std::size_t nodeCount = 0;
124  for(std::size_t b=0;b<nodeBuckets.size();b++)
125  nodeCount += nodeBuckets[b]->size();
126 
128  = Teuchos::rcp(new std::vector<std::size_t>(nodeCount));
129 
130  // loop over node buckets
131  for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
132  stk::mesh::Bucket & bucket = *nodeBuckets[b];
133 
134  for(std::size_t n=0;n<bucket.size();n++,index++)
135  (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
136  }
137 
138  return sideIds;
139 }
140 
141 std::pair<Teuchos::RCP<std::vector<std::size_t> >,
144  const std::string & sideName, const std::string type_)
145 {
146  unsigned physicalDim = mesh.getDimension();
147 
150 
151  // grab nodes owned by requested side
153  std::stringstream ss;
154  ss << "Can't find part=\"" << sideName << "\"" << std::endl;
155  stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
156  stk::mesh::Selector mySides = (*side) & metaData->locally_owned_part();
157 
158  stk::mesh::EntityRank rank;
160  stk::mesh::EntityId offset = 0;
161  if(type_ == "coord"){
162  rank = mesh.getNodeRank();
163  field = & mesh.getCoordinatesField();
164  } else if(type_ == "edge"){
165  rank = mesh.getEdgeRank();
166  field = & mesh.getEdgesField();
167  offset = mesh.getMaxEntityId(mesh.getNodeRank());
168  } else if(type_ == "face"){
169  rank = mesh.getFaceRank();
170  field = & mesh.getFacesField();
171  offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
172  } else {
173  ss << "Can't do BCs of type " << type_ << std::endl;
174  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
175  }
176 
177  std::vector<stk::mesh::Bucket*> const& nodeBuckets =
178  bulkData->get_buckets(rank, mySides);
179 
180  // build id vector
182  std::size_t nodeCount = 0;
183  for(std::size_t b=0;b<nodeBuckets.size();b++)
184  nodeCount += nodeBuckets[b]->size();
185 
187  = Teuchos::rcp(new std::vector<std::size_t>(nodeCount));
189  = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(nodeCount));
190 
191  // loop over node buckets
192  for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
193  stk::mesh::Bucket & bucket = *nodeBuckets[b];
194  double const* array = stk::mesh::field_data(*field, bucket);
195 
196  for(std::size_t n=0;n<bucket.size();n++,index++) {
197  (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
198  Teuchos::Tuple<double,3> & coord = (*sideCoords)[index];
199 
200  // copy coordinates into multi vector
201  for(std::size_t d=0;d<physicalDim;d++)
202  coord[d] = array[physicalDim*n + d];
203 
204  // need to ensure that higher dimensions are properly zeroed
205  // required for 1D periodic boundary conditions
206  for(std::size_t d=physicalDim;d<3;d++)
207  coord[d] = 0;
208  }
209  }
210 
211  return std::make_pair(sideIds,sideCoords);
212 }
213 
214 std::pair<Teuchos::RCP<std::vector<std::size_t> >,
217  const std::string & sideName, const std::string type_)
218 {
219  using Teuchos::RCP;
220  using Teuchos::rcp;
221  using LO = panzer::LocalOrdinal;
222  using GO = panzer::GlobalOrdinal;
224  using Map = Tpetra::Map<LO,GO,NODE>;
225  using Importer = Tpetra::Import<LO,GO,NODE>;
226 
227  // Epetra_MpiComm Comm(mesh.getBulkData()->parallel());
228  auto comm = mesh.getComm();
229 
230  unsigned physicalDim = mesh.getDimension();
231 
232  // grab local IDs and coordinates on this side
233  // and build local epetra vector
235 
236  std::pair<Teuchos::RCP<std::vector<std::size_t> >,
238  getLocalSideIdsAndCoords(mesh,sideName,type_);
239 
240  std::vector<std::size_t> & local_side_ids = *sidePair.first;
241  std::vector<Teuchos::Tuple<double,3> > & local_side_coords = *sidePair.second;
242  std::size_t nodeCount = local_side_ids.size();
243 
244  // build local Tpetra objects
246  RCP<Map> idMap_ = rcp(new Map(computeInternally,nodeCount,0,comm));
247  RCP<Tpetra::Vector<GO,LO,GO,NODE>> localIdVec_ = rcp(new Tpetra::Vector<GO,LO,GO,NODE>(idMap_));
248  RCP<Tpetra::MultiVector<double,LO,GO,NODE>> localCoordVec_ = rcp(new Tpetra::MultiVector<double,LO,GO,NODE>(idMap_,physicalDim));
249 
250  // copy local Ids and coords into Tpetra vectors
251  {
252  auto lidHost = localIdVec_->getLocalViewHost(Tpetra::Access::OverwriteAll);
253  auto lcoordHost = localCoordVec_->getLocalViewHost(Tpetra::Access::OverwriteAll);
254  for(std::size_t n=0;n<local_side_ids.size();n++) {
255  std::size_t nodeId = local_side_ids[n];
256  Teuchos::Tuple<double,3> & coords = local_side_coords[n];
257 
258  lidHost(n,0) = static_cast<GO>(nodeId);
259  for(unsigned d=0;d<physicalDim;d++)
260  lcoordHost(n,d) = coords[d];
261  }
262  }
263 
264  // fully distribute epetra vector across all processors
265  // (these are "distributed" or "dist" objects)
267 
268  std::size_t dist_nodeCount = idMap_->getGlobalNumElements();
269 
270  // build global Tpetra objects
271  RCP<Map> distMap_ = rcp(new Map(dist_nodeCount,0,comm,Tpetra::LocallyReplicated));
272  RCP<Tpetra::Vector<GO,LO,GO,NODE>> distIdVec_ = rcp(new Tpetra::Vector<GO,LO,GO,NODE>(distMap_));
273  RCP<Tpetra::MultiVector<double,LO,GO,NODE>> distCoordVec_ = rcp(new Tpetra::MultiVector<double,LO,GO,NODE>(distMap_,physicalDim));
274 
275  // export to the localVec object from the "vector" object
276  Importer importer_(idMap_,distMap_);
277  distIdVec_->doImport(*localIdVec_,importer_,Tpetra::INSERT);
278  distCoordVec_->doImport(*localCoordVec_,importer_,Tpetra::INSERT);
279 
280  // convert back to generic stl vector objects
282 
284  = Teuchos::rcp(new std::vector<std::size_t>(dist_nodeCount));
286  = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(dist_nodeCount));
287 
288  // copy local Ids from Tpetra vector
289  const auto didHost = distIdVec_->getLocalViewHost(Tpetra::Access::ReadOnly);
290  const auto dcoordHost = distCoordVec_->getLocalViewHost(Tpetra::Access::ReadOnly);
291  for(std::size_t n=0;n<dist_side_ids->size();++n) {
292  (*dist_side_ids)[n] = didHost(n,0);
293 
294  Teuchos::Tuple<double,3> & coords = (*dist_side_coords)[n];
295  for(unsigned d=0;d<physicalDim;++d) {
296  coords[d] = dcoordHost(n,d);
297  }
298  // ensure that higher dimensions are zero
299  for(unsigned d=physicalDim;d<3;++d)
300  coords[d] = 0;
301  }
302 
303  return std::make_pair(dist_side_ids,dist_side_coords);
304 }
305 
306 } // end periodic_helpers
307 
308 } // end panzer_stk
stk::mesh::Field< double > VectorFieldType
Teuchos::RCP< std::vector< std::size_t > > getLocalSideIds(const STK_Interface &mesh, const std::string &sideName, const std::string type_)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const
get max entity ID of type entityRank
std::pair< Teuchos::RCP< std::vector< std::size_t > >, Teuchos::RCP< std::vector< Teuchos::Tuple< double, 3 > > > > getLocalSideIdsAndCoords(const STK_Interface &mesh, const std::string &sideName, const std::string type_)
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we&#39;re performing.
unsigned getDimension() const
get the dimension
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
stk::mesh::EntityRank getFaceRank() const
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
const VectorFieldType & getEdgesField() const
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we&#39;ll contribute, or in which we&#39;ll store, the result of computing this integral...
stk::mesh::EntityRank getEdgeRank() const
Teuchos::RCP< std::vector< std::pair< std::size_t, std::size_t > > > getGlobalPairing(const std::vector< std::size_t > &locallyRequiredIds, const std::vector< std::pair< std::size_t, std::size_t > > &locallyMatchedIds, const STK_Interface &mesh, bool failure)
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
const VectorFieldType & getFacesField() const
#define TEUCHOS_ASSERT(assertion_test)
stk::mesh::EntityRank getNodeRank() const
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType
std::pair< Teuchos::RCP< std::vector< std::size_t > >, Teuchos::RCP< std::vector< Teuchos::Tuple< double, 3 > > > > getSideIdsAndCoords(const STK_Interface &mesh, const std::string &sideName, const std::string type_)
const VectorFieldType & getCoordinatesField() const