45 #include <stk_mesh/base/GetEntities.hpp> 46 #include <stk_mesh/base/GetBuckets.hpp> 47 #include <stk_mesh/base/FieldBase.hpp> 50 #include "Tpetra_Map.hpp" 51 #include "Tpetra_Import.hpp" 52 #include "Tpetra_Vector.hpp" 54 #include "Teuchos_FancyOStream.hpp" 58 namespace periodic_helpers {
60 Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >
62 const std::vector<std::pair<std::size_t,std::size_t> > & locallyMatchedIds,
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>;
75 int myVal = failure ? 1 : 0;
77 Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_SUM,1,&myVal,&sumVal);
78 TEUCHOS_ASSERT(sumVal==0);
80 std::vector<GO> requiredInts(locallyRequiredIds.size());
81 for(std::size_t i=0;i<requiredInts.size();i++)
82 requiredInts[i] = locallyRequiredIds[i];
84 std::vector<GO> providedInts(locallyMatchedIds.size());
85 for(std::size_t i=0;i<locallyMatchedIds.size();i++)
86 providedInts[i] = locallyMatchedIds[i].first;
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);
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;
100 Tpetra::Vector<GO,LO,GO,NODE> requiredVector(requiredMap);
101 requiredVector.doImport(providedVector,importer,Tpetra::INSERT);
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()));
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);
120 Teuchos::RCP<std::vector<std::size_t> >
122 const std::string & sideName,
const std::string type_)
124 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.
getMetaData();
125 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
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());
134 stk::mesh::EntityRank rank;
135 panzer::GlobalOrdinal offset = 0;
136 if(type_ ==
"coord"){
138 }
else if(type_ ==
"edge"){
141 }
else if(type_ ==
"face"){
145 ss <<
"Can't do BCs of type " << type_ << std::endl;
146 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error, ss.str())
149 std::vector<stk::mesh::Bucket*>
const& nodeBuckets =
150 bulkData->get_buckets(rank, mySides);
154 std::size_t nodeCount = 0;
155 for(std::size_t b=0;b<nodeBuckets.size();b++)
156 nodeCount += nodeBuckets[b]->size();
158 Teuchos::RCP<std::vector<std::size_t> > sideIds
159 = Teuchos::rcp(
new std::vector<std::size_t>(nodeCount));
162 for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
163 stk::mesh::Bucket & bucket = *nodeBuckets[b];
165 for(std::size_t n=0;n<bucket.size();n++,index++)
166 (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
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_)
179 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.
getMetaData();
180 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
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();
189 stk::mesh::EntityRank rank;
191 unsigned int offset = 0;
192 if(type_ ==
"coord"){
195 }
else if(type_ ==
"edge"){
199 }
else if(type_ ==
"face"){
204 ss <<
"Can't do BCs of type " << type_ << std::endl;
205 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error, ss.str())
208 std::vector<stk::mesh::Bucket*>
const& nodeBuckets =
209 bulkData->get_buckets(rank, mySides);
213 std::size_t nodeCount = 0;
214 for(std::size_t b=0;b<nodeBuckets.size();b++)
215 nodeCount += nodeBuckets[b]->size();
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));
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);
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];
232 for(std::size_t d=0;d<physicalDim;d++)
233 coord[d] = array[physicalDim*n + d];
237 for(std::size_t d=physicalDim;d<3;d++)
242 return std::make_pair(sideIds,sideCoords);
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_)
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>;
267 std::pair<Teuchos::RCP<std::vector<std::size_t> >,
268 Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > > sidePair =
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();
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));
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];
288 lidHost(n,0) =
static_cast<GO
>(nodeId);
289 for(
unsigned d=0;d<physicalDim;d++)
290 lcoordHost(n,d) = coords[d];
297 int dist_nodeCount = idMap_->getGlobalNumElements();
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));
305 Importer importer_(idMap_,distMap_);
306 distIdVec_->doImport(*localIdVec_,importer_,Tpetra::INSERT);
307 distCoordVec_->doImport(*localCoordVec_,importer_,Tpetra::INSERT);
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));
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);
323 Teuchos::Tuple<double,3> & coords = (*dist_side_coords)[n];
324 for(
unsigned d=0;d<physicalDim;++d) {
325 coords[d] = dcoordHost(n,d);
328 for(
unsigned d=physicalDim;d<3;++d)
332 return std::make_pair(dist_side_ids,dist_side_coords);
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'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'll contribute, or in which we'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