Panzer Version of the Day
Loading...
Searching...
No Matches
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
56namespace panzer_stk {
57
58namespace periodic_helpers {
59
60Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >
61getGlobalPairing(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;
67 using NODE = panzer::TpetraNodeType;
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 {
97 auto pvHost = providedVector.getLocalViewHost(Tpetra::Access::OverwriteAll);
98 for(std::size_t i=0;i<locallyMatchedIds.size();i++)
99 pvHost(i,0) = locallyMatchedIds[i].second;
100 }
101
102 Tpetra::Vector<GO,LO,GO,NODE> requiredVector(requiredMap);
103 requiredVector.doImport(providedVector,importer,Tpetra::INSERT);
104
105
106 Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > > result
107 = Teuchos::rcp(new std::vector<std::pair<std::size_t,std::size_t> >(requiredInts.size()));
108
109 auto rvHost = requiredVector.getLocalViewHost(Tpetra::Access::ReadOnly);
110 for(std::size_t i=0;i<result->size();i++) {
111 (*result)[i].first = requiredInts[i];
112 (*result)[i].second = rvHost(i,0);
113 }
114 return result;
115}
116
117
118
122Teuchos::RCP<std::vector<std::size_t> >
124 const std::string & sideName, const std::string type_)
125{
126 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
127 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
128
129 // grab nodes owned by requested side
131 std::stringstream ss;
132 ss << "Can't find part=\"" << sideName << "\"" << std::endl;
133 stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
134 stk::mesh::Selector mySides = *side & (metaData->locally_owned_part() | metaData->globally_shared_part());
135
136 stk::mesh::EntityRank rank;
137 panzer::GlobalOrdinal offset = 0; // offset to avoid giving nodes, edges, faces the same sideId
138 if(type_ == "coord"){
139 rank = mesh.getNodeRank();
140 } else if(type_ == "edge"){
141 rank = mesh.getEdgeRank();
142 offset = mesh.getMaxEntityId(mesh.getNodeRank());
143 } else if(type_ == "face"){
144 rank = mesh.getFaceRank();
145 offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
146 } else {
147 ss << "Can't do BCs of type " << type_ << std::endl;
148 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
149 }
150
151 std::vector<stk::mesh::Bucket*> const& nodeBuckets =
152 bulkData->get_buckets(rank, mySides);
153
154 // build id vector
156 std::size_t nodeCount = 0;
157 for(std::size_t b=0;b<nodeBuckets.size();b++)
158 nodeCount += nodeBuckets[b]->size();
159
160 Teuchos::RCP<std::vector<std::size_t> > sideIds
161 = Teuchos::rcp(new std::vector<std::size_t>(nodeCount));
162
163 // loop over node buckets
164 for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
165 stk::mesh::Bucket & bucket = *nodeBuckets[b];
166
167 for(std::size_t n=0;n<bucket.size();n++,index++)
168 (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
169 }
170
171 return sideIds;
172}
173
174std::pair<Teuchos::RCP<std::vector<std::size_t> >,
175 Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > >
177 const std::string & sideName, const std::string type_)
178{
179 unsigned physicalDim = mesh.getDimension();
180
181 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
182 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
183
184 // grab nodes owned by requested side
186 std::stringstream ss;
187 ss << "Can't find part=\"" << sideName << "\"" << std::endl;
188 stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
189 stk::mesh::Selector mySides = (*side) & metaData->locally_owned_part();
190
191 stk::mesh::EntityRank rank;
193 stk::mesh::EntityId offset = 0;
194 if(type_ == "coord"){
195 rank = mesh.getNodeRank();
196 field = & mesh.getCoordinatesField();
197 } else if(type_ == "edge"){
198 rank = mesh.getEdgeRank();
199 field = & mesh.getEdgesField();
200 offset = mesh.getMaxEntityId(mesh.getNodeRank());
201 } else if(type_ == "face"){
202 rank = mesh.getFaceRank();
203 field = & mesh.getFacesField();
204 offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
205 } else {
206 ss << "Can't do BCs of type " << type_ << std::endl;
207 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
208 }
209
210 std::vector<stk::mesh::Bucket*> const& nodeBuckets =
211 bulkData->get_buckets(rank, mySides);
212
213 // build id vector
215 std::size_t nodeCount = 0;
216 for(std::size_t b=0;b<nodeBuckets.size();b++)
217 nodeCount += nodeBuckets[b]->size();
218
219 Teuchos::RCP<std::vector<std::size_t> > sideIds
220 = Teuchos::rcp(new std::vector<std::size_t>(nodeCount));
221 Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > sideCoords
222 = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(nodeCount));
223
224 // loop over node buckets
225 for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
226 stk::mesh::Bucket & bucket = *nodeBuckets[b];
227 double const* array = stk::mesh::field_data(*field, bucket);
228
229 for(std::size_t n=0;n<bucket.size();n++,index++) {
230 (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
231 Teuchos::Tuple<double,3> & coord = (*sideCoords)[index];
232
233 // copy coordinates into multi vector
234 for(std::size_t d=0;d<physicalDim;d++)
235 coord[d] = array[physicalDim*n + d];
236
237 // need to ensure that higher dimensions are properly zeroed
238 // required for 1D periodic boundary conditions
239 for(std::size_t d=physicalDim;d<3;d++)
240 coord[d] = 0;
241 }
242 }
243
244 return std::make_pair(sideIds,sideCoords);
245}
246
247std::pair<Teuchos::RCP<std::vector<std::size_t> >,
248 Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > >
250 const std::string & sideName, const std::string type_)
251{
252 using Teuchos::RCP;
253 using Teuchos::rcp;
254 using LO = panzer::LocalOrdinal;
255 using GO = panzer::GlobalOrdinal;
256 using NODE = panzer::TpetraNodeType;
257 using Map = Tpetra::Map<LO,GO,NODE>;
258 using Importer = Tpetra::Import<LO,GO,NODE>;
259
260 // Epetra_MpiComm Comm(mesh.getBulkData()->parallel());
261 auto comm = mesh.getComm();
262
263 unsigned physicalDim = mesh.getDimension();
264
265 // grab local IDs and coordinates on this side
266 // and build local epetra vector
268
269 std::pair<Teuchos::RCP<std::vector<std::size_t> >,
270 Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > > sidePair =
271 getLocalSideIdsAndCoords(mesh,sideName,type_);
272
273 std::vector<std::size_t> & local_side_ids = *sidePair.first;
274 std::vector<Teuchos::Tuple<double,3> > & local_side_coords = *sidePair.second;
275 std::size_t nodeCount = local_side_ids.size();
276
277 // build local Tpetra objects
278 auto computeInternally = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
279 RCP<Map> idMap_ = rcp(new Map(computeInternally,nodeCount,0,comm));
280 RCP<Tpetra::Vector<GO,LO,GO,NODE>> localIdVec_ = rcp(new Tpetra::Vector<GO,LO,GO,NODE>(idMap_));
281 RCP<Tpetra::MultiVector<double,LO,GO,NODE>> localCoordVec_ = rcp(new Tpetra::MultiVector<double,LO,GO,NODE>(idMap_,physicalDim));
282
283 // copy local Ids and coords into Tpetra vectors
284 {
285 auto lidHost = localIdVec_->getLocalViewHost(Tpetra::Access::OverwriteAll);
286 auto lcoordHost = localCoordVec_->getLocalViewHost(Tpetra::Access::OverwriteAll);
287 for(std::size_t n=0;n<local_side_ids.size();n++) {
288 std::size_t nodeId = local_side_ids[n];
289 Teuchos::Tuple<double,3> & coords = local_side_coords[n];
290
291 lidHost(n,0) = static_cast<GO>(nodeId);
292 for(unsigned d=0;d<physicalDim;d++)
293 lcoordHost(n,d) = coords[d];
294 }
295 }
296
297 // fully distribute epetra vector across all processors
298 // (these are "distributed" or "dist" objects)
300
301 std::size_t dist_nodeCount = idMap_->getGlobalNumElements();
302
303 // build global Tpetra objects
304 RCP<Map> distMap_ = rcp(new Map(dist_nodeCount,0,comm,Tpetra::LocallyReplicated));
305 RCP<Tpetra::Vector<GO,LO,GO,NODE>> distIdVec_ = rcp(new Tpetra::Vector<GO,LO,GO,NODE>(distMap_));
306 RCP<Tpetra::MultiVector<double,LO,GO,NODE>> distCoordVec_ = rcp(new Tpetra::MultiVector<double,LO,GO,NODE>(distMap_,physicalDim));
307
308 // export to the localVec object from the "vector" object
309 Importer importer_(idMap_,distMap_);
310 distIdVec_->doImport(*localIdVec_,importer_,Tpetra::INSERT);
311 distCoordVec_->doImport(*localCoordVec_,importer_,Tpetra::INSERT);
312
313 // convert back to generic stl vector objects
315
316 Teuchos::RCP<std::vector<std::size_t> > dist_side_ids
317 = Teuchos::rcp(new std::vector<std::size_t>(dist_nodeCount));
318 Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > dist_side_coords
319 = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(dist_nodeCount));
320
321 // copy local Ids from Tpetra vector
322 const auto didHost = distIdVec_->getLocalViewHost(Tpetra::Access::ReadOnly);
323 const auto dcoordHost = distCoordVec_->getLocalViewHost(Tpetra::Access::ReadOnly);
324 for(std::size_t n=0;n<dist_side_ids->size();++n) {
325 (*dist_side_ids)[n] = didHost(n,0);
326
327 Teuchos::Tuple<double,3> & coords = (*dist_side_coords)[n];
328 for(unsigned d=0;d<physicalDim;++d) {
329 coords[d] = dcoordHost(n,d);
330 }
331 // ensure that higher dimensions are zero
332 for(unsigned d=physicalDim;d<3;++d)
333 coords[d] = 0;
334 }
335
336 return std::make_pair(dist_side_ids,dist_side_coords);
337}
338
339} // end periodic_helpers
340
341} // end panzer_stk
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we're performing.
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
const VectorFieldType & getEdgesField() const
stk::mesh::Field< double, stk::mesh::Cartesian > VectorFieldType
stk::mesh::EntityRank getNodeRank() const
const VectorFieldType & getFacesField() const
stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const
get max entity ID of type entityRank
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
stk::mesh::EntityRank getFaceRank() const
const VectorFieldType & getCoordinatesField() const
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
unsigned getDimension() const
get the dimension
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
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)
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_)
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_)
Teuchos::RCP< std::vector< std::size_t > > getLocalSideIds(const STK_Interface &mesh, const std::string &sideName, const std::string type_)
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType