Panzer  Version of the Day
Panzer_STK_PeriodicBC_Matcher.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 
45 #include <stk_mesh/base/GetEntities.hpp>
46 #include <stk_mesh/base/GetBuckets.hpp>
47 #include <stk_mesh/base/FieldBase.hpp>
48 
49 #include "Panzer_NodeType.hpp"
50 #include "Tpetra_Map.hpp"
51 #include "Tpetra_Import.hpp"
52 #include "Tpetra_Vector.hpp"
53 
54 #include "Teuchos_FancyOStream.hpp"
55 
56 namespace panzer_stk {
57 
58 namespace periodic_helpers {
59 
60 Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >
61 getGlobalPairing(const std::vector<std::size_t> & locallyRequiredIds,
62  const std::vector<std::pair<std::size_t,std::size_t> > & locallyMatchedIds,
63  const STK_Interface & mesh,bool failure)
64 {
65  using LO = panzer::LocalOrdinal;
66  using GO = panzer::GlobalOrdinal;
68  using Map = Tpetra::Map<LO,GO,NODE>;
69  using Importer = Tpetra::Import<LO,GO,NODE>;
70 
71  auto comm = mesh.getComm();
72 
73  // this is needed to prevent hanging: it is unfortunately expensive
74  // need a better way!
75  int myVal = failure ? 1 : 0;
76  int sumVal = 0;
77  Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_SUM,1,&myVal,&sumVal);
78  TEUCHOS_ASSERT(sumVal==0);
79 
80  std::vector<GO> requiredInts(locallyRequiredIds.size());
81  for(std::size_t i=0;i<requiredInts.size();i++)
82  requiredInts[i] = locallyRequiredIds[i];
83 
84  std::vector<GO> providedInts(locallyMatchedIds.size());
85  for(std::size_t i=0;i<locallyMatchedIds.size();i++)
86  providedInts[i] = locallyMatchedIds[i].first;
87 
88  // maps and communciation all set up
89  auto computeInternally = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
90  Teuchos::RCP<Map> requiredMap = Teuchos::rcp(new Map(computeInternally,Teuchos::ArrayView<const GO>(requiredInts),0,comm));
91  Teuchos::RCP<Map> providedMap = Teuchos::rcp(new Map(computeInternally,Teuchos::ArrayView<const GO>(providedInts),0,comm));
92  Importer importer(providedMap,requiredMap);
93 
94  // this is what to distribute
95  Tpetra::Vector<GO,LO,GO,NODE> providedVector(providedMap);
96  auto pvHost = providedVector.getLocalViewHost(Tpetra::Access::OverwriteAll);
97  for(std::size_t i=0;i<locallyMatchedIds.size();i++)
98  pvHost(i,0) = locallyMatchedIds[i].second;
99 
100  Tpetra::Vector<GO,LO,GO,NODE> requiredVector(requiredMap);
101  requiredVector.doImport(providedVector,importer,Tpetra::INSERT);
102 
103 
104  Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > > result
105  = Teuchos::rcp(new std::vector<std::pair<std::size_t,std::size_t> >(requiredInts.size()));
106 
107  auto rvHost = requiredVector.getLocalViewHost(Tpetra::Access::ReadOnly);
108  for(std::size_t i=0;i<result->size();i++) {
109  (*result)[i].first = requiredInts[i];
110  (*result)[i].second = rvHost(i,0);
111  }
112  return result;
113 }
114 
115 
116 
120 Teuchos::RCP<std::vector<std::size_t> >
122  const std::string & sideName, const std::string type_)
123 {
124  Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
125  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
126 
127  // grab nodes owned by requested side
129  std::stringstream ss;
130  ss << "Can't find part=\"" << sideName << "\"" << std::endl;
131  stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
132  stk::mesh::Selector mySides = *side & (metaData->locally_owned_part() | metaData->globally_shared_part());
133 
134  stk::mesh::EntityRank rank;
135  panzer::GlobalOrdinal offset = 0; // offset to avoid giving nodes, edges, faces the same sideId
136  if(type_ == "coord"){
137  rank = mesh.getNodeRank();
138  } else if(type_ == "edge"){
139  rank = mesh.getEdgeRank();
140  offset = mesh.getMaxEntityId(mesh.getNodeRank());
141  } else if(type_ == "face"){
142  rank = mesh.getFaceRank();
143  offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
144  } else {
145  ss << "Can't do BCs of type " << type_ << std::endl;
146  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
147  }
148 
149  std::vector<stk::mesh::Bucket*> const& nodeBuckets =
150  bulkData->get_buckets(rank, mySides);
151 
152  // build id vector
154  std::size_t nodeCount = 0;
155  for(std::size_t b=0;b<nodeBuckets.size();b++)
156  nodeCount += nodeBuckets[b]->size();
157 
158  Teuchos::RCP<std::vector<std::size_t> > sideIds
159  = Teuchos::rcp(new std::vector<std::size_t>(nodeCount));
160 
161  // loop over node buckets
162  for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
163  stk::mesh::Bucket & bucket = *nodeBuckets[b];
164 
165  for(std::size_t n=0;n<bucket.size();n++,index++)
166  (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
167  }
168 
169  return sideIds;
170 }
171 
172 std::pair<Teuchos::RCP<std::vector<std::size_t> >,
173  Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > >
175  const std::string & sideName, const std::string type_)
176 {
177  unsigned physicalDim = mesh.getDimension();
178 
179  Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
180  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
181 
182  // grab nodes owned by requested side
184  std::stringstream ss;
185  ss << "Can't find part=\"" << sideName << "\"" << std::endl;
186  stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
187  stk::mesh::Selector mySides = (*side) & metaData->locally_owned_part();
188 
189  stk::mesh::EntityRank rank;
191  unsigned int offset = 0;
192  if(type_ == "coord"){
193  rank = mesh.getNodeRank();
194  field = & mesh.getCoordinatesField();
195  } else if(type_ == "edge"){
196  rank = mesh.getEdgeRank();
197  field = & mesh.getEdgesField();
198  offset = mesh.getMaxEntityId(mesh.getNodeRank());
199  } else if(type_ == "face"){
200  rank = mesh.getFaceRank();
201  field = & mesh.getFacesField();
202  offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
203  } else {
204  ss << "Can't do BCs of type " << type_ << std::endl;
205  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
206  }
207 
208  std::vector<stk::mesh::Bucket*> const& nodeBuckets =
209  bulkData->get_buckets(rank, mySides);
210 
211  // build id vector
213  std::size_t nodeCount = 0;
214  for(std::size_t b=0;b<nodeBuckets.size();b++)
215  nodeCount += nodeBuckets[b]->size();
216 
217  Teuchos::RCP<std::vector<std::size_t> > sideIds
218  = Teuchos::rcp(new std::vector<std::size_t>(nodeCount));
219  Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > sideCoords
220  = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(nodeCount));
221 
222  // loop over node buckets
223  for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
224  stk::mesh::Bucket & bucket = *nodeBuckets[b];
225  double const* array = stk::mesh::field_data(*field, bucket);
226 
227  for(std::size_t n=0;n<bucket.size();n++,index++) {
228  (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
229  Teuchos::Tuple<double,3> & coord = (*sideCoords)[index];
230 
231  // copy coordinates into multi vector
232  for(std::size_t d=0;d<physicalDim;d++)
233  coord[d] = array[physicalDim*n + d];
234 
235  // need to ensure that higher dimensions are properly zeroed
236  // required for 1D periodic boundary conditions
237  for(std::size_t d=physicalDim;d<3;d++)
238  coord[d] = 0;
239  }
240  }
241 
242  return std::make_pair(sideIds,sideCoords);
243 }
244 
245 std::pair<Teuchos::RCP<std::vector<std::size_t> >,
246  Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > >
248  const std::string & sideName, const std::string type_)
249 {
250  using Teuchos::RCP;
251  using Teuchos::rcp;
252  using LO = panzer::LocalOrdinal;
253  using GO = panzer::GlobalOrdinal;
255  using Map = Tpetra::Map<LO,GO,NODE>;
256  using Importer = Tpetra::Import<LO,GO,NODE>;
257 
258  // Epetra_MpiComm Comm(mesh.getBulkData()->parallel());
259  auto comm = mesh.getComm();
260 
261  unsigned physicalDim = mesh.getDimension();
262 
263  // grab local IDs and coordinates on this side
264  // and build local epetra vector
266 
267  std::pair<Teuchos::RCP<std::vector<std::size_t> >,
268  Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > > sidePair =
269  getLocalSideIdsAndCoords(mesh,sideName,type_);
270 
271  std::vector<std::size_t> & local_side_ids = *sidePair.first;
272  std::vector<Teuchos::Tuple<double,3> > & local_side_coords = *sidePair.second;
273  int nodeCount = local_side_ids.size();
274 
275  // build local Tpetra objects
276  auto computeInternally = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
277  RCP<Map> idMap_ = rcp(new Map(computeInternally,nodeCount,0,comm));
278  RCP<Tpetra::Vector<GO,LO,GO,NODE>> localIdVec_ = rcp(new Tpetra::Vector<GO,LO,GO,NODE>(idMap_));
279  RCP<Tpetra::MultiVector<double,LO,GO,NODE>> localCoordVec_ = rcp(new Tpetra::MultiVector<double,LO,GO,NODE>(idMap_,physicalDim));
280 
281  // copy local Ids and coords into Tpetra vectors
282  auto lidHost = localIdVec_->getLocalViewHost(Tpetra::Access::OverwriteAll);
283  auto lcoordHost = localCoordVec_->getLocalViewHost(Tpetra::Access::OverwriteAll);
284  for(std::size_t n=0;n<local_side_ids.size();n++) {
285  std::size_t nodeId = local_side_ids[n];
286  Teuchos::Tuple<double,3> & coords = local_side_coords[n];
287 
288  lidHost(n,0) = static_cast<GO>(nodeId);
289  for(unsigned d=0;d<physicalDim;d++)
290  lcoordHost(n,d) = coords[d];
291  }
292 
293  // fully distribute epetra vector across all processors
294  // (these are "distributed" or "dist" objects)
296 
297  int dist_nodeCount = idMap_->getGlobalNumElements();
298 
299  // build global Tpetra objects
300  RCP<Map> distMap_ = rcp(new Map(dist_nodeCount,0,comm,Tpetra::LocallyReplicated));
301  RCP<Tpetra::Vector<GO,LO,GO,NODE>> distIdVec_ = rcp(new Tpetra::Vector<GO,LO,GO,NODE>(distMap_));
302  RCP<Tpetra::MultiVector<double,LO,GO,NODE>> distCoordVec_ = rcp(new Tpetra::MultiVector<double,LO,GO,NODE>(distMap_,physicalDim));
303 
304  // export to the localVec object from the "vector" object
305  Importer importer_(idMap_,distMap_);
306  distIdVec_->doImport(*localIdVec_,importer_,Tpetra::INSERT);
307  distCoordVec_->doImport(*localCoordVec_,importer_,Tpetra::INSERT);
308 
309  // convert back to generic stl vector objects
311 
312  Teuchos::RCP<std::vector<std::size_t> > dist_side_ids
313  = Teuchos::rcp(new std::vector<std::size_t>(dist_nodeCount));
314  Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > dist_side_coords
315  = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(dist_nodeCount));
316 
317  // copy local Ids from Tpetra vector
318  const auto didHost = distIdVec_->getLocalViewHost(Tpetra::Access::ReadOnly);
319  const auto dcoordHost = distCoordVec_->getLocalViewHost(Tpetra::Access::ReadOnly);
320  for(std::size_t n=0;n<dist_side_ids->size();++n) {
321  (*dist_side_ids)[n] = didHost(n,0);
322 
323  Teuchos::Tuple<double,3> & coords = (*dist_side_coords)[n];
324  for(unsigned d=0;d<physicalDim;++d) {
325  coords[d] = dcoordHost(n,d);
326  }
327  // ensure that higher dimensions are zero
328  for(unsigned d=physicalDim;d<3;++d)
329  coords[d] = 0;
330  }
331 
332  return std::make_pair(dist_side_ids,dist_side_coords);
333 }
334 
335 } // end periodic_helpers
336 
337 } // end panzer_stk
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
Teuchos::RCP< std::vector< std::size_t > > getLocalSideIds(const STK_Interface &mesh, const std::string &sideName, const std::string type_)
stk::mesh::EntityRank getFaceRank() const
stk::mesh::EntityRank getNodeRank() const
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
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_)
stk::mesh::EntityRank getEdgeRank() const
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
const VectorFieldType & getFacesField() const
stk::mesh::Field< double, stk::mesh::Cartesian > VectorFieldType
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
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...
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)
const VectorFieldType & getCoordinatesField() const
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 > > > > getSideIdsAndCoords(const STK_Interface &mesh, const std::string &sideName, const std::string type_)
const VectorFieldType & getEdgesField() const