13 #include <stk_mesh/base/GetEntities.hpp>
14 #include <stk_mesh/base/FieldBase.hpp>
17 #include "Tpetra_Map.hpp"
18 #include "Tpetra_Import.hpp"
19 #include "Tpetra_Vector.hpp"
21 #include "Teuchos_FancyOStream.hpp"
23 #ifdef PANZER_HAVE_STKSEARCH
24 namespace panzer_stk {
26 namespace periodic_helpers {
28 void fillLocalSearchVector(
const STK_Interface & mesh, SphereIdVector & searchVector,
const double & error,
29 const std::string & sideName,
const std::string & type_,
const bool & getGhostedIDs)
37 std::vector<std::string> matchedSides;
38 std::vector<SearchId> potentialIDsToRemap;
40 fillLocalSearchVector(mesh,searchVector,error,sideName,type_,getGhostedIDs,matchedSides,potentialIDsToRemap);
46 void fillLocalSearchVector(
const STK_Interface & mesh, SphereIdVector & searchVector,
const double & error,
47 const std::string & sideName,
const std::string & type_,
const bool & getGhostedIDs,
48 const std::vector<std::string> & matchedSides, std::vector<SearchId> & potentialIDsToRemap)
50 unsigned physicalDim = mesh.getDimension();
55 const unsigned parallelRank = bulkData->parallel_rank();
60 ss <<
"Can't find a sideset named \"" << sideName <<
"\" in the mesh" << std::endl;
61 stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
64 stk::mesh::Selector mySides = getGhostedIDs ?
65 (*side) & (metaData->locally_owned_part() | metaData->globally_shared_part()) :
66 (*side) & metaData->locally_owned_part();
68 stk::mesh::EntityRank rank;
72 stk::mesh::EntityId offset = 0;
74 rank = mesh.getNodeRank();
75 field = & mesh.getCoordinatesField();
77 }
else if(type_ ==
"edge"){
78 rank = mesh.getEdgeRank();
79 field = & mesh.getEdgesField();
80 offset = mesh.getMaxEntityId(mesh.getNodeRank());
81 }
else if(type_ ==
"face"){
82 rank = mesh.getFaceRank();
83 field = & mesh.getFacesField();
84 offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
86 ss <<
"Can't do BCs of type " << type_ << std::endl;
92 stk::mesh::Selector intersection;
94 if (matchedSides.size()>0) {
95 for (
size_t j=0; j<matchedSides.size(); ++j) {
96 auto previouslyMatched = matchedSides[j];
98 intersection = intersection | (mySides & *(metaData->get_part(previouslyMatched)));
101 mySides = mySides - intersection;
105 std::vector<stk::mesh::Bucket*>
const& entityBuckets =
106 bulkData->get_buckets(rank, mySides);
110 for(
size_t b=0;b<entityBuckets.size();b++) {
111 stk::mesh::Bucket & bucket = *entityBuckets[b];
112 double const* array = stk::mesh::field_data(*field, bucket);
115 for(
size_t n=0;n<bucket.size();n++) {
119 for(
size_t d=0;d<physicalDim;d++)
120 coord[d] = array[physicalDim*n + d];
124 for(
size_t d=physicalDim;d<3;d++)
130 stk::search::Point<double> center(coord[0],coord[1],coord[2]);
131 stk::mesh::EntityKey trueKey = bulkData->entity_key(bucket[n]);
132 stk::mesh::EntityKey shiftedKey(trueKey.rank(), trueKey.id()+offset);
133 SearchId search_id(shiftedKey, parallelRank);
135 searchVector.emplace_back( Sphere(center, error), search_id);
142 if (matchedSides.size()>0) {
145 mySides = getGhostedIDs ?
146 (*side) & (metaData->locally_owned_part() | metaData->globally_shared_part()) :
147 (*side) & metaData->locally_owned_part();
149 std::vector<stk::mesh::Bucket*>
const & intersectionEntityBuckets =
150 bulkData->get_buckets(rank, mySides & intersection);
154 for(
size_t b=0;b<intersectionEntityBuckets.size();b++) {
155 stk::mesh::Bucket & bucket = *intersectionEntityBuckets[b];
157 for(
size_t n=0;n<bucket.size();n++) {
158 stk::mesh::EntityKey trueKey = bulkData->entity_key(bucket[n]);
159 stk::mesh::EntityKey shiftedKey(trueKey.rank(), trueKey.id()+offset);
160 potentialIDsToRemap.emplace_back(shiftedKey, parallelRank);
168 const std::vector<double> computeGlobalCentroid(
const STK_Interface & mesh,
const std::string & sideName)
171 unsigned physicalDim = mesh.getDimension();
178 std::stringstream ss;
179 ss <<
"Can't find part=\"" << sideName <<
"\"" << std::endl;
180 stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
181 stk::mesh::Selector mySides = (*side) & metaData->locally_owned_part();
184 std::vector<stk::mesh::Bucket*>
const& entityBuckets =
185 bulkData->get_buckets(mesh.getNodeRank(), mySides);
188 double localCentroid[3] = {0.,0.,0.};
189 int localNodeCount = 0;
190 for(
size_t b=0;b<entityBuckets.size();b++) {
191 stk::mesh::Bucket & bucket = *entityBuckets[b];
192 double const* array = stk::mesh::field_data(mesh.getCoordinatesField(), bucket);
195 for(
size_t n=0;n<bucket.size();n++) {
199 for(
size_t d=0;d<physicalDim;d++)
200 localCentroid[d] += array[physicalDim*n + d];
204 int globalNodeCount = 0;
206 auto comm = mesh.getComm();
208 double globalCentroid[3] = { };
213 std::vector<double>
result = {0.,0.,0.};
214 for (
size_t d=0;d<physicalDim;d++)
215 result[d] = globalCentroid[d]/globalNodeCount;
220 void updateMapping(
Teuchos::RCP<std::vector<std::pair<size_t,size_t> > > & currentMatches,
221 const std::vector<std::pair<size_t,size_t> > & previousMatches,
222 const std::vector<SearchId> & IDsToRemap,
const STK_Interface & mesh)
225 using LO = panzer::LocalOrdinal;
226 using GO = panzer::GlobalOrdinal;
228 using Map = Tpetra::Map<LO,GO,NODE>;
229 using Importer = Tpetra::Import<LO,GO,NODE>;
231 auto comm = mesh.getComm();
237 std::map<size_t,size_t> myPreviousAtoB,myCurrentAtoB;
238 std::map<size_t,std::vector<size_t> > myPreviousBtoA;
239 for (
size_t i=0;i<previousMatches.size();++i) {
240 myPreviousAtoB[previousMatches[i].first] = previousMatches[i].second;
242 myPreviousBtoA[previousMatches[i].second].push_back(previousMatches[i].first);
244 for (
size_t i=0;i<currentMatches->size();++i)
245 myCurrentAtoB[(*currentMatches)[i].first] = (*currentMatches)[i].second;
249 std::vector<GO> requestedAIDs;
251 for (
auto &
id : IDsToRemap)
252 requestedAIDs.push_back(myPreviousAtoB[
id.
id().
id()]);
255 std::set<GO> uniqueAIDs(requestedAIDs.begin(),requestedAIDs.end());
256 requestedAIDs = std::vector<GO>(uniqueAIDs.begin(),uniqueAIDs.end());
259 std::vector<GO> newAIDs(currentMatches->size());
261 for (
size_t i=0;i<currentMatches->size();++i) {
262 newAIDs[i] = (*currentMatches)[i].first;
270 if (!testMap->isOneToOne()){
271 newAIDsMap = Tpetra::createOneToOne<LO,GO,NODE>(testMap);
273 newAIDsMap = testMap;
278 Importer importer(newAIDsMap,requestedAIDsMap);
281 Tpetra::Vector<GO,LO,GO,NODE> newBIDs(newAIDsMap);
282 auto newBIDsHost = newBIDs.getLocalViewHost(Tpetra::Access::OverwriteAll);
283 auto myGIDs = newAIDsMap->getMyGlobalIndices();
284 for (
size_t i=0;i<myGIDs.size();++i)
285 newBIDsHost(i,0) = myCurrentAtoB[myGIDs[i]];
287 Tpetra::Vector<GO,LO,GO,NODE> requestedBIDs(requestedAIDsMap);
288 requestedBIDs.doImport(newBIDs,importer,Tpetra::INSERT);
289 auto requestedBIDsHost = requestedBIDs.getLocalViewHost(Tpetra::Access::ReadOnly);
296 for (
const auto &
id : requestedAIDs) {
298 for (
const auto & idToUpdate : myPreviousBtoA[
id])
300 myPreviousAtoB[idToUpdate] = requestedBIDsHost(ind,0);
306 for (
const auto & AB : previousMatches) {
310 (*currentMatches).emplace_back(std::pair<size_t,size_t>(
id,myPreviousAtoB[
id]));
316 void appendMapping(
Teuchos::RCP<std::vector<std::pair<size_t,size_t> > > & currentMatches,
317 const std::vector<std::pair<size_t,size_t> > & previousMatches)
320 for (
const auto & AB : previousMatches)
321 (*currentMatches).push_back(AB);
stk::mesh::Field< double > VectorFieldType
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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...
#define TEUCHOS_ASSERT(assertion_test)
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType