Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_PeriodicBC_Search.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#ifdef PANZER_HAVE_STKSEARCH
57namespace panzer_stk {
58
59namespace periodic_helpers {
60
61void fillLocalSearchVector(const STK_Interface & mesh, SphereIdVector & searchVector, const double & error,
62 const std::string & sideName, const std::string & type_, const bool & getGhostedIDs)
63{
64
65 // empty optional arguments for real call below
66 // this is partially for backwards compatability but also recognizes that for the first
67 // pairing, these vectors are not needed
68 // they are also not used for side B in every case
69
70 std::vector<std::string> matchedSides;
71 std::vector<SearchId> potentialIDsToRemap;
72
73 fillLocalSearchVector(mesh,searchVector,error,sideName,type_,getGhostedIDs,matchedSides,potentialIDsToRemap);
74
75 return;
76
77}
78
79void fillLocalSearchVector(const STK_Interface & mesh, SphereIdVector & searchVector, const double & error,
80 const std::string & sideName, const std::string & type_, const bool & getGhostedIDs,
81 const std::vector<std::string> & matchedSides, std::vector<SearchId> & potentialIDsToRemap)
82{
83 unsigned physicalDim = mesh.getDimension();
84
85 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
86 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
87
88 const unsigned parallelRank = bulkData->parallel_rank();
89
90 // grab entities owned by requested side
92 std::stringstream ss;
93 ss << "Can't find a sideset named \"" << sideName << "\" in the mesh" << std::endl;
94 stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
95
96 // if ghosted IDs are requested, add in the shared portion
97 stk::mesh::Selector mySides = getGhostedIDs ?
98 (*side) & (metaData->locally_owned_part() | metaData->globally_shared_part()) :
99 (*side) & metaData->locally_owned_part();
100
101 stk::mesh::EntityRank rank;
103 // the sorting further downstream only uses ids which are offset
104 // based on the entity type
105 stk::mesh::EntityId offset = 0;
106 if(type_ == "coord"){
107 rank = mesh.getNodeRank();
108 field = & mesh.getCoordinatesField();
109 // no offset
110 } else if(type_ == "edge"){
111 rank = mesh.getEdgeRank();
112 field = & mesh.getEdgesField();
113 offset = mesh.getMaxEntityId(mesh.getNodeRank());
114 } else if(type_ == "face"){
115 rank = mesh.getFaceRank();
116 field = & mesh.getFacesField();
117 offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
118 } else {
119 ss << "Can't do BCs of type " << type_ << std::endl;
120 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
121 }
122
123 // remove previously matched nodes
124
125 stk::mesh::Selector intersection; // empty set to start
126
127 if (matchedSides.size()>0) {
128 for (size_t j=0; j<matchedSides.size(); ++j) {
129 auto previouslyMatched = matchedSides[j];
130 // add in the overlap between the requested side and the previously matched side
131 intersection = intersection | (mySides & *(metaData->get_part(previouslyMatched)));
132 }
133 // remove all overlaps
134 mySides = mySides - intersection;
135 }
136
137 // get buckets
138 std::vector<stk::mesh::Bucket*> const& entityBuckets =
139 bulkData->get_buckets(rank, mySides);
140
141 // loop over entity buckets
142
143 for(size_t b=0;b<entityBuckets.size();b++) {
144 stk::mesh::Bucket & bucket = *entityBuckets[b];
145 double const* array = stk::mesh::field_data(*field, bucket);
146
147 // loop over entities
148 for(size_t n=0;n<bucket.size();n++) {
149
150 double coord[3]; // coordinates
151 // copy coordinates into multi vector
152 for(size_t d=0;d<physicalDim;d++)
153 coord[d] = array[physicalDim*n + d];
154
155 // need to ensure that higher dimensions are properly zeroed
156 // required for 1D periodic boundary conditions
157 for(size_t d=physicalDim;d<3;d++)
158 coord[d] = 0;
159
160 // add to the coordinate and id to the search vector
161 // a tolerance can be specified
162 // TODO allow for relative tolerances...
163 stk::search::Point<double> center(coord[0],coord[1],coord[2]);
164 stk::mesh::EntityKey trueKey = bulkData->entity_key(bucket[n]);
165 stk::mesh::EntityKey shiftedKey(trueKey.rank(), trueKey.id()+offset);
166 SearchId search_id(shiftedKey, parallelRank);
167
168 searchVector.emplace_back( Sphere(center, error), search_id);
169 }
170 }
171
172 // for multiperiodic case, populate the potentialIDsToRemap vector with the IDs that have
173 // already been matched and fall on this side
174
175 if (matchedSides.size()>0) {
176 TEUCHOS_ASSERT(potentialIDsToRemap.size()==0);
177 // reset mySides
178 mySides = getGhostedIDs ?
179 (*side) & (metaData->locally_owned_part() | metaData->globally_shared_part()) :
180 (*side) & metaData->locally_owned_part();
181
182 std::vector<stk::mesh::Bucket*> const & intersectionEntityBuckets =
183 bulkData->get_buckets(rank, mySides & intersection);
184
185 // loop over entity buckets
186
187 for(size_t b=0;b<intersectionEntityBuckets.size();b++) {
188 stk::mesh::Bucket & bucket = *intersectionEntityBuckets[b];
189 // loop over entities
190 for(size_t n=0;n<bucket.size();n++) {
191 stk::mesh::EntityKey trueKey = bulkData->entity_key(bucket[n]);
192 stk::mesh::EntityKey shiftedKey(trueKey.rank(), trueKey.id()+offset);
193 potentialIDsToRemap.emplace_back(shiftedKey, parallelRank);
194 }
195 }
196 }
197
198 return;
199}
200
201const std::vector<double> computeGlobalCentroid(const STK_Interface & mesh, const std::string & sideName)
202{
203 // TODO too much replicated code here
204 unsigned physicalDim = mesh.getDimension();
205
206 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
207 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
208
209 // grab requested side
211 std::stringstream ss;
212 ss << "Can't find part=\"" << sideName << "\"" << std::endl;
213 stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
214 stk::mesh::Selector mySides = (*side) & metaData->locally_owned_part();
215
216 // get node buckets
217 std::vector<stk::mesh::Bucket*> const& entityBuckets =
218 bulkData->get_buckets(mesh.getNodeRank(), mySides);
219
220 // loop over node buckets
221 double localCentroid[3] = {0.,0.,0.};
222 int localNodeCount = 0;
223 for(size_t b=0;b<entityBuckets.size();b++) {
224 stk::mesh::Bucket & bucket = *entityBuckets[b];
225 double const* array = stk::mesh::field_data(mesh.getCoordinatesField(), bucket);
226
227 // loop over nodes
228 for(size_t n=0;n<bucket.size();n++) {
229
230 ++localNodeCount;
231 // sum (note that unused dimensions are skipped)
232 for(size_t d=0;d<physicalDim;d++)
233 localCentroid[d] += array[physicalDim*n + d];
234
235 }
236 }
237 int globalNodeCount = 0;
238
239 auto comm = mesh.getComm();
240
241 double globalCentroid[3] = { };
242
243 Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,1,&localNodeCount,&globalNodeCount);
244 Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,3,&localCentroid[0],&globalCentroid[0]);
245
246 std::vector<double> result = {0.,0.,0.};
247 for (size_t d=0;d<physicalDim;d++)
248 result[d] = globalCentroid[d]/globalNodeCount;
249
250 return result;
251}
252
253void updateMapping(Teuchos::RCP<std::vector<std::pair<size_t,size_t> > > & currentMatches,
254 const std::vector<std::pair<size_t,size_t> > & previousMatches,
255 const std::vector<SearchId> & IDsToRemap, const STK_Interface & mesh)
256{
257
258 using LO = panzer::LocalOrdinal;
259 using GO = panzer::GlobalOrdinal;
261 using Map = Tpetra::Map<LO,GO,NODE>;
262 using Importer = Tpetra::Import<LO,GO,NODE>;
263
264 auto comm = mesh.getComm();
265
266 // store maps
267 // this is necessary because of the uniqueness requirements
268 // and convenient to update the map
269
270 std::map<size_t,size_t> myPreviousAtoB,myCurrentAtoB;
271 std::map<size_t,std::vector<size_t> > myPreviousBtoA;
272 for (size_t i=0;i<previousMatches.size();++i) {
273 myPreviousAtoB[previousMatches[i].first] = previousMatches[i].second;
274 // may not be one-to-one
275 myPreviousBtoA[previousMatches[i].second].push_back(previousMatches[i].first);
276 }
277 for (size_t i=0;i<currentMatches->size();++i)
278 myCurrentAtoB[(*currentMatches)[i].first] = (*currentMatches)[i].second;
279
280 // find which IDs we need to query to get THEIR A to B map
281 // this means we need to the the B id of our previous match for the IDsToRemap
282 std::vector<GO> requestedAIDs;
283
284 for (auto & id : IDsToRemap)
285 requestedAIDs.push_back(myPreviousAtoB[id.id().id()]);
286
287 // quick and dirty way to get rid of repeated entries on each process
288 std::set<GO> uniqueAIDs(requestedAIDs.begin(),requestedAIDs.end());
289 requestedAIDs = std::vector<GO>(uniqueAIDs.begin(),uniqueAIDs.end());
290
291 // find which A ids we have in our new mapping
292 std::vector<GO> newAIDs(currentMatches->size());
293
294 for (size_t i=0;i<currentMatches->size();++i) {
295 newAIDs[i] = (*currentMatches)[i].first;
296 }
297
298 // create teuchos maps
299 auto computeInternally = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
300 Teuchos::RCP<const Map> testMap = Teuchos::rcp(new Map(computeInternally,&newAIDs[0],newAIDs.size(),0,comm));
301 Teuchos::RCP<const Map> newAIDsMap;
302 // source must be unique across communicator
303 if (!testMap->isOneToOne()){
304 newAIDsMap = Tpetra::createOneToOne<LO,GO,NODE>(testMap);
305 } else {
306 newAIDsMap = testMap;
307 }
308
309 Teuchos::RCP<Map> requestedAIDsMap = Teuchos::rcp(new Map(computeInternally,&requestedAIDs[0],requestedAIDs.size(),0,comm));
310
311 Importer importer(newAIDsMap,requestedAIDsMap);
312
313 // send out my A to B map
314 Tpetra::Vector<GO,LO,GO,NODE> newBIDs(newAIDsMap);
315 auto newBIDsHost = newBIDs.getLocalViewHost(Tpetra::Access::OverwriteAll);
316 auto myGIDs = newAIDsMap->getMyGlobalIndices();
317 for (size_t i=0;i<myGIDs.size();++i)
318 newBIDsHost(i,0) = myCurrentAtoB[myGIDs[i]];
319
320 Tpetra::Vector<GO,LO,GO,NODE> requestedBIDs(requestedAIDsMap);
321 requestedBIDs.doImport(newBIDs,importer,Tpetra::INSERT);
322 auto requestedBIDsHost = requestedBIDs.getLocalViewHost(Tpetra::Access::ReadOnly);
323
324 // now overwrite previous map where necessary...
325 // what about error checking? what is the default is something is requested but not there?
326 // TODO this assumes that teuchos maps and communication does not
327 // alter the ordering in anyway so that AIDs and IDsToRemap correspond appropriately
328 size_t ind = 0;
329 for (const auto & id : requestedAIDs) {
330 // get the corresponding ids to update in the previous map
331 for (const auto & idToUpdate : myPreviousBtoA[id])
332 // update with the final B id
333 myPreviousAtoB[idToUpdate] = requestedBIDsHost(ind,0);
334 ++ind;
335 }
336
337 // and add to new map...
338 // needs to respect the previous ordering or type_vec in getPeriodicNodePairing will be wrong
339 for (const auto & AB : previousMatches) {
340 // so we get the A ids in order
341 auto id = AB.first;
342 // and use the updated previous A to B map
343 (*currentMatches).emplace_back(std::pair<size_t,size_t>(id,myPreviousAtoB[id]));
344 }
345
346 return;
347}
348
349void appendMapping(Teuchos::RCP<std::vector<std::pair<size_t,size_t> > > & currentMatches,
350 const std::vector<std::pair<size_t,size_t> > & previousMatches)
351{
352 // add previous mapping to new map
353 for (const auto & AB : previousMatches)
354 (*currentMatches).push_back(AB);
355
356 return;
357}
358
359} // end periodic_helpers
360
361} // end panzer_stk
362#endif
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.
stk::mesh::Field< double, stk::mesh::Cartesian > VectorFieldType
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType