SkdTreeUtils.cpp
1 /*---------------------------------------------------------------------------*\
2  *
3  * mimmo
4  *
5  * Copyright (C) 2015-2021 OPTIMAD engineering Srl
6  *
7  * -------------------------------------------------------------------------
8  * License
9  * This file is part of mimmo.
10  *
11  * mimmo is free software: you can redistribute it and/or modify it
12  * under the terms of the GNU Lesser General Public License v3 (LGPL)
13  * as published by the Free Software Foundation.
14  *
15  * mimmo is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public License
21  * along with mimmo. If not, see <http://www.gnu.org/licenses/>.
22  *
23 \*---------------------------------------------------------------------------*/
24 
25 # include "SkdTreeUtils.hpp"
26 # include "MimmoCGUtils.hpp"
27 # include <bitpit_surfunstructured.hpp>
28 # include <surface_skd_tree.hpp>
29 # include <CG.hpp>
30 # include <queue>
31 
32 namespace mimmo{
33 
34 namespace skdTreeUtils{
35 
49 double distance(const std::array<double,3> *point, const bitpit::PatchSkdTree *tree, long &id, double r)
50 {
51 
52  double h = std::numeric_limits<double>::max();
53  if(!tree ){
54  throw std::runtime_error("Invalid use of skdTreeUtils::distance method: a void tree is detected.");
55  }
56  if(!dynamic_cast<const bitpit::SurfUnstructured*>(&(tree->getPatch()))){
57  throw std::runtime_error("Invalid use of skdTreeUtils::distance method: a not surface patch tree is detected.");
58  }
59 
60  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestCell(*point, r, &id, &h);
61  return h;
62 
63 }
64 
84 double signedDistance(const std::array<double,3> *point, const bitpit::PatchSkdTree *tree, long &id, std::array<double,3> &normal, double r)
85 {
86 
87  double h = std::numeric_limits<double>::max();
88 
89  if(!tree ){
90  throw std::runtime_error("Invalid use of skdTreeUtils::signedDistance method: a void tree is detected.");
91  }
92  const bitpit::SurfUnstructured *spatch = dynamic_cast<const bitpit::SurfUnstructured*>(&(tree->getPatch()));
93  if(!spatch){
94  throw std::runtime_error("Invalid use of skdTreeUtils::signedDistance method: a not surface patch tree is detected.");
95  }
96 
97  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestCell(*point, r, &id, &h);
98 
99  double s = computePseudoNormal(*point, spatch, id, normal);
100  return s*h;
101 
102 }
103 
119 void distance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, double *distances, double r)
120 {
121  std::vector<double> rs(nP, r);
122  distance(nP, points, tree, ids, distances, rs.data());
123 }
124 
140 void distance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, double *distances, double *r)
141 {
142 
143  // Initialize distances
144  for (int ip = 0; ip < nP; ip++){
145  distances[ip] = std::numeric_limits<double>::max();
146  }
147 
148  if(!tree ){
149  throw std::runtime_error("Invalid use of skdTreeUtils::distance method: a void tree is detected.");
150  }
151  if(!dynamic_cast<const bitpit::SurfUnstructured*>(&(tree->getPatch()))){
152  throw std::runtime_error("Invalid use of skdTreeUtils::distance method: a not surface patch tree is detected.");
153  }
154 
155  for (int ip = 0; ip < nP; ip++){
156  const std::array<double,3> *point = &points[ip];
157  long *id = &ids[ip];
158  double *distance = &distances[ip];
159  double rpoint = r[ip];
160  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestCell(*point, rpoint, id, distance);
161  }
162 
163 }
164 
186 void signedDistance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, double *distances, std::array<double,3> *normals, double r)
187 {
188  std::vector<double> rs(nP, r);
189  signedDistance(nP, points, tree, ids, distances, normals, rs.data());
190 }
191 
213 void signedDistance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, double *distances, std::array<double,3> *normals, double *r)
214 {
215 
216  // Initialize distances
217  for (int ip = 0; ip < nP; ip++){
218  distances[ip] = std::numeric_limits<double>::max();
219  }
220 
221  if(!tree ){
222  throw std::runtime_error("Invalid use of skdTreeUtils::signedDistance method: a void tree is detected.");
223  }
224  const bitpit::SurfUnstructured *spatch = dynamic_cast<const bitpit::SurfUnstructured*>(&(tree->getPatch()));
225  if(!spatch){
226  throw std::runtime_error("Invalid use of skdTreeUtils::signedDistance method: a not surface patch tree is detected.");
227  }
228 
229  for (int ip = 0; ip < nP; ip++){
230  const std::array<double,3> *point = &points[ip];
231  long *id = &ids[ip];
232  double *distance = &distances[ip];
233  double rpoint = r[ip];
234  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestCell(*point, rpoint, id, distance);
235  double s = computePseudoNormal(*point, spatch, *id, normals[ip]);
236  *distance *= s;
237  }
238 
239 }
240 
251 std::vector<long> selectByPatch(bitpit::PatchSkdTree *selection, bitpit::PatchSkdTree *target, double tol){
252 
253  int nleafs = selection->getLeafCount();
254  int nnodess = selection->getNodeCount();
255  std::vector<const bitpit::SkdNode*> leafSelection(nleafs);
256  int count = 0;
257  for (int i=0; i<nnodess; i++){
258  if (selection->getNode(i).isLeaf()){
259  if (bitpit::CGElem::intersectBoxBox(selection->getNode(i).getBoxMin(),
260  selection->getNode(i).getBoxMax(),
261  target->getNode(0).getBoxMin(),
262  target->getNode(0).getBoxMax(),
263  3,
264  tol ) ){
265  leafSelection[count] = &(selection->getNode(i));
266  count++;
267 
268  }
269  }
270  }
271  leafSelection.resize(count);
272 
273 
274  std::vector<long> extracted;
275  extractTarget(target, leafSelection, extracted, tol);
276 
277  return extracted;
278 
279 }
280 
296 void extractTarget(bitpit::PatchSkdTree *target, const std::vector<const bitpit::SkdNode*> & leafSelection, std::vector<long> &extracted, double tol){
297 
298  if(leafSelection.empty()) return;
299  std::size_t rootId = 0;
300 
301  std::vector<long> candidates;
302  std::vector<std::pair< std::size_t, std::vector<const bitpit::SkdNode*> > > nodeStack;
303 
304  std::vector<const bitpit::SkdNode*> tocheck;
305  for (int i=0; i<(int)leafSelection.size(); i++){
306  if (bitpit::CGElem::intersectBoxBox(leafSelection[i]->getBoxMin(),
307  leafSelection[i]->getBoxMax(),
308  target->getNode(rootId).getBoxMin(),
309  target->getNode(rootId).getBoxMax(),
310  3,
311  tol ) )
312  {
313  tocheck.push_back(leafSelection[i]);
314  }
315  }
316  nodeStack.push_back(std::make_pair(rootId, tocheck) );
317 
318  while(!nodeStack.empty()){
319 
320  std::pair<long, std::vector<const bitpit::SkdNode*> > touple = nodeStack.back();
321  const bitpit::SkdNode & node = target->getNode(touple.first);
322  nodeStack.pop_back();
323 
324 
325  bool isLeaf = true;
326  for (int i = bitpit::SkdNode::CHILD_BEGIN; i != bitpit::SkdNode::CHILD_END; ++i) {
327  bitpit::SkdNode::ChildLocation childLocation = static_cast<bitpit::SkdNode::ChildLocation>(i);
328  std::size_t childId = node.getChildId(childLocation);
329  if (childId != bitpit::SkdNode::NULL_ID) {
330  isLeaf = false;
331  tocheck.clear();
332  for (int i=0; i<(int)touple.second.size(); i++){
333  if (bitpit::CGElem::intersectBoxBox(touple.second[i]->getBoxMin(),
334  touple.second[i]->getBoxMax(),
335  target->getNode(childId).getBoxMin(),
336  target->getNode(childId).getBoxMax(),
337  3,
338  tol ) )
339  {
340  tocheck.push_back(touple.second[i]);
341  }
342  }
343  if(!tocheck.empty()) nodeStack.push_back(std::make_pair(childId, tocheck) );
344  }
345  }
346 
347  if (isLeaf) {
348  candidates.push_back(touple.first);
349  }
350  }
351 
352  std::unordered_set<long> cellExtracted;
353  for(const auto & nodeId: candidates){
354  const bitpit::SkdNode & node = target->getNode(nodeId);
355  std::vector<long> cellids = node.getCells();
356  cellExtracted.insert(cellids.begin(), cellids.end());
357  }
358 
359  extracted.insert(extracted.end(), cellExtracted.begin(), cellExtracted.end());
360 }
361 
375 darray3E projectPoint(const std::array<double,3> *point, const bitpit::PatchSkdTree *skdtree, double r )
376 {
377  darray3E projected_point;
378  long id;
379  projectPoint(1, point, skdtree, &projected_point, &id, r);
380  return projected_point;
381 }
382 
398 void projectPoint(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, std::array<double,3> *projected_points, long *ids, double r )
399 {
400  std::vector<double> rs(nP, r);
401  projectPoint(nP, points, tree, projected_points, ids, rs.data());
402 }
403 
419 void projectPoint(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, std::array<double,3> *projected_points, long *ids, double *r )
420 {
421  if(nP == 0) return;
422  // Initialize ids and ranks
423  for (int ip = 0; ip < nP; ip++){
424  ids[ip] = bitpit::Cell::NULL_ID;
425  r[ip] = std::max(r[ip], tree->getPatch().getTol());
426  }
427  std::vector<darray3E> normals(nP);
428 
429  std::vector<double> dist(nP, std::numeric_limits<double>::max());
430  double max_dist = std::numeric_limits<double>::max();
431  while (max_dist == std::numeric_limits<double>::max()){
432  //use method sphere by default
433  signedDistance(nP, points, tree, ids, dist.data(), normals.data(), r);
434  for (int ip = 0; ip < nP; ip++){
435  if (std::abs(dist[ip]) >= max_dist){
436  r[ip] *= 1.5;
437  }
438  }
439  max_dist = std::abs(*std::max_element(dist.begin(), dist.end()));
440  max_dist = std::max(max_dist, std::abs(*std::min_element(dist.begin(), dist.end())));
441  }
442 
443  for (int ip = 0; ip < nP; ip++){
444  projected_points[ip] = points[ip] - dist[ip] * normals[ip];
445  }
446 }
447 
456 long locatePointOnPatch(const std::array<double, 3> &point, const bitpit::PatchSkdTree *tree)
457 {
458  if(!dynamic_cast<const bitpit::SurfUnstructured*>(&(tree->getPatch()))){
459  throw std::runtime_error("Invalid use of skdTreeUtils::locatePointOnPatch method: a non surface patch or void patch was detected.");
460  }
461 
462  // Initialize the cell id and distance
463  long id = bitpit::Cell::NULL_ID;
464  double distance = std::numeric_limits<double>::max();
465 
466  // Find the closest cell to the input point
467  long cellId = bitpit::Cell::NULL_ID;
468  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestCell(point, &cellId, &distance);
469 
470  // Check if the point belongs to the closest cell
471  bool checkBelong = false;
472 
473  bitpit::ElementType eltype = tree->getPatch().getCell(cellId).getType();
474  auto vList = tree->getPatch().getCell(cellId).getVertexIds();
475  dvecarr3E vertCoords(vList.size());
476  int counterList = 0;
477  for(const auto & idV : vList){
478  vertCoords[counterList] = tree->getPatch().getVertexCoords(idV);
479  ++counterList;
480  }
481 
482  switch(eltype){
483  case bitpit::ElementType::LINE:
484  checkBelong = mimmoCGUtils::isPointInsideSegment(point, vertCoords[0], vertCoords[1], tree->getPatch().getTol());
485  break;
486  case bitpit::ElementType::TRIANGLE:
487  checkBelong = mimmoCGUtils::isPointInsideTriangle(point, vertCoords[0], vertCoords[1], vertCoords[2], tree->getPatch().getTol());
488  break;
489  case bitpit::ElementType::PIXEL:
490  case bitpit::ElementType::QUAD:
491  case bitpit::ElementType::POLYGON:
492  checkBelong = mimmoCGUtils::isPointInsidePolygon(point, vertCoords, tree->getPatch().getTol());
493  break;
494  default:
495  throw std::runtime_error("Not supported cell type");
496  break;
497  }
498  if (checkBelong) {
499  id = cellId;
500  }
501 
502  return id;
503 }
504 
516 double
517 computePseudoNormal(const std::array<double, 3> &point, const bitpit::SurfUnstructured *surface_mesh, long id, std::array<double, 3> &pseudo_normal)
518 {
519 
520  pseudo_normal.fill(0.);
521  double h, s(std::numeric_limits<double>::max());
522 
523  //Pseudo normal only for 2D element patches
524  if (id != bitpit::Cell::NULL_ID){
525 
526  const bitpit::Cell & cell = surface_mesh->getCell(id);
527  bitpit::ConstProxyVector<long> vertIds = cell.getVertexIds();
528  dvecarr3E VS(vertIds.size());
529  int count = 0;
530  for (const auto & iV: vertIds){
531  VS[count] = surface_mesh->getVertexCoords(iV);
532  ++count;
533  }
534 
535  darray3E xP = {{0.0,0.0,0.0}};
536  darray3E normal= {{0.0,0.0,0.0}};
537 
538  if ( vertIds.size() == 3 ){ //TRIANGLE
539  darray3E lambda;
540  h = bitpit::CGElem::distancePointTriangle(point, VS[0], VS[1], VS[2],lambda);
541  int count = 0;
542  for(const auto &val: lambda){
543  normal += val * surface_mesh->evalVertexNormal(id,count) ;
544  xP += val * VS[count];
545  ++count;
546  }
547  }else if ( vertIds.size() == 2 ){ //LINE/SEGMENT
548  darray2E lambda;
549  h = bitpit::CGElem::distancePointSegment(point, VS[0], VS[1], lambda);
550  int count = 0;
551  for(const auto &val: lambda){
552  normal += val * surface_mesh->evalVertexNormal(id,count) ;
553  xP += val * VS[count];
554  ++count;
555  }
556  }else{ //GENERAL POLYGON
557  std::vector<double> lambda;
558  h = bitpit::CGElem::distancePointPolygon(point, VS,lambda);
559  int count = 0;
560  for(const auto &val: lambda){
561  normal += val * surface_mesh->evalVertexNormal(id,count) ;
562  xP += val * VS[count];
563  ++count;
564  }
565  }
566 
567  s = sign( dotProduct(normal, point - xP) );
568  if(s == 0.0) s =1.0;
569  h = s * h;
570  //pseudo-normal (direction P and xP closest point on triangle)
571  pseudo_normal = s * (point - xP);
572  double normX = norm2(pseudo_normal);
573  if(normX < surface_mesh->getTol()){
574  pseudo_normal = normal/norm2(normal);
575  }else{
576  pseudo_normal /= normX;
577  }
578  }//end if not id null
579 
580  return s;
581 }
582 
591 bool
592 checkPointBelongsToCell(const std::array<double, 3> &point, const bitpit::SurfUnstructured *surface_mesh, long cellId)
593 {
594  bool checkBelong = false;
595  bitpit::ElementType eltype = surface_mesh->getCell(cellId).getType();
596  auto vList = surface_mesh->getCell(cellId).getVertexIds();
597  dvecarr3E vertCoords(vList.size());
598  int counterList = 0;
599  for(const auto & idV : vList){
600  vertCoords[counterList] = surface_mesh->getVertexCoords(idV);
601  ++counterList;
602  }
603 
604  switch(eltype){
605  case bitpit::ElementType::LINE:
606  checkBelong = mimmoCGUtils::isPointInsideSegment(point, vertCoords[0], vertCoords[1], surface_mesh->getTol());
607  break;
608  case bitpit::ElementType::TRIANGLE:
609  checkBelong = mimmoCGUtils::isPointInsideTriangle(point, vertCoords[0], vertCoords[1], vertCoords[2], surface_mesh->getTol());
610  break;
611  case bitpit::ElementType::PIXEL:
612  case bitpit::ElementType::QUAD:
613  case bitpit::ElementType::POLYGON:
614  checkBelong = mimmoCGUtils::isPointInsidePolygon(point, vertCoords, surface_mesh->getTol());
615  break;
616  default:
617  throw std::runtime_error("Not supported cell type");
618  break;
619  }
620  return checkBelong;
621 }
622 
623 #if MIMMO_ENABLE_MPI
624 
637 void globalDistance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, int *ranks, double *distances, double r, bool shared)
638 {
639  std::vector<double> rs(nP, r);
640  globalDistance(nP, points, tree, ids, ranks, distances, rs.data(), shared);
641 }
642 
656 void globalDistance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, int *ranks, double *distances, double* r, bool shared)
657 {
658 
659  if(!tree ){
660  throw std::runtime_error("Invalid use of skdTreeUtils::distance method: a void tree is detected.");
661  }
662  if(!dynamic_cast<const bitpit::SurfUnstructured*>(&(tree->getPatch()))){
663  throw std::runtime_error("Invalid use of skdTreeUtils::distance method: a not surface patch tree is detected.");
664  }
665 
666  if (shared){
667  findSharedPointClosestGlobalCell(nP, points, tree, ids, ranks, distances, r);
668  } else {
669  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestGlobalCell(nP, points, r, ids, ranks, distances);
670  }
671 
672 }
673 
690 void signedGlobalDistance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, int *ranks, std::array<double,3> *normals, double *distances, double r, bool shared)
691 {
692  std::vector<double> rs(nP, r);
693  signedGlobalDistance(nP, points, tree, ids, ranks, normals, distances, rs.data(), shared);
694 }
695 
712 void signedGlobalDistance(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, int *ranks, std::array<double,3> *normals, double *distances, double *r, bool shared)
713 {
714 
715  if(!tree){
716  throw std::runtime_error("Invalid use of skdTreeUtils::signedGlobalDistance method: a void tree is detected.");
717  }
718  const bitpit::SurfUnstructured & spatch = static_cast<const bitpit::SurfUnstructured&>(tree->getPatch());
719 
720  if (shared){
721  // All processes share the points
722  findSharedPointClosestGlobalCell(nP, points, tree, ids, ranks, distances, r);
723  } else {
724  // Distributed points
725  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestGlobalCell(nP, points, r, ids, ranks, distances);
726  }
727 
728 
729  int myrank = spatch.getRank();
730 
731 
732  // Map with rank to ask info in case of distributed points
733  std::map<int,std::vector<long>> ask_to_rank;
734  std::map<int,std::vector<std::array<double,3>>> point_to_rank;
735  std::map<int,std::vector<int>> point_index_received_from_rank;
736 
737  //Global signs container
738  std::vector<double> signs;
739  if(shared) signs.resize(nP, std::numeric_limits<double>::max());
740 
741  // Loop on points
742  for (int ip = 0; ip < nP; ip++){
743 
744  const std::array<double,3> & point = points[ip];
745  const long & cellId = ids[ip];
746  const int & cellRank = ranks[ip];
747  double & distance = distances[ip];
748  std::array<double,3> & pseudo_normal = normals[ip];
749 
750  // If the current rank is the the owner of the cell compute normal
751  if (cellRank == myrank){
752 
753  double s = computePseudoNormal(point, &spatch, cellId, pseudo_normal);
754  if(cellId != bitpit::Cell::NULL_ID){
755  if (shared) {
756  signs[ip] = s;
757  }else{
758  distance *= s;
759  }
760  }
761 
762  } else {
763 
764  // If not owner fix the normal to maximum value
765  pseudo_normal = std::array<double,3>({std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()});
766 
767  if (!shared){
768  // Push to rank map to ask normal
769  ask_to_rank[cellRank].push_back(cellId);
770 
771  // Push to rank map the related points to ask normal
772  point_to_rank[cellRank].push_back(point);
773 
774  // Push to rank map the index of the points to fill output normals
775  point_index_received_from_rank[cellRank].push_back(ip);
776  }
777 
778  } // end if owner rank
779 
780  } // end loop on points
781 
782  int nProcessors = spatch.getProcessorCount();
783 
784  if (nProcessors == 1 && shared){
785  for (int ip = 0; ip < nP; ip++){
786  const long & cellId = ids[ip];
787  double & distance = distances[ip];
788  double & s = signs[ip];
789  if(cellId != bitpit::Cell::NULL_ID){
790  distance *= s;
791  }
792  }
793  }
794 
795  if (nProcessors > 1 && spatch.isPartitioned()){
796  if (shared){
797 
798  // If all points are shared between processes all reduce on normal by minimum operation
799  MPI_Allreduce(MPI_IN_PLACE, normals, 3*nP, MPI_DOUBLE, MPI_MIN, tree->getPatch().getCommunicator());
800  MPI_Allreduce(MPI_IN_PLACE, signs.data(), nP, MPI_DOUBLE, MPI_MIN, tree->getPatch().getCommunicator());
801  for (int ip = 0; ip < nP; ip++){
802  const long & cellId = ids[ip];
803  double & distance = distances[ip];
804  double & s = signs[ip];
805  if(cellId != bitpit::Cell::NULL_ID){
806  distance *= s;
807  }
808  }
809  } else {
810 
811  // Send/receive number of normals to send to other processes
812  std::vector<int> n_send_to_rank(nProcessors, 0);
813  for (int irank = 0; irank < nProcessors; irank++){
814  int n_ask_to_current_rank = ask_to_rank[irank].size();
815  MPI_Gather(&n_ask_to_current_rank, 1, MPI_INT, n_send_to_rank.data(), 1, MPI_INT, irank, tree->getPatch().getCommunicator());
816  }
817 
818  // Map with rank to send info for cellId and point coordinates
819  std::map<int,std::vector<long>> send_to_rank;
820  std::map<int,std::vector<std::array<double,3>>> point_from_rank;
821  for (int irank = 0; irank < nProcessors; irank++){
822  send_to_rank[irank].resize(n_send_to_rank[irank], bitpit::Cell::NULL_ID);
823  point_from_rank[irank].resize(n_send_to_rank[irank], std::array<double,3>({{0.,0.,0.}}));
824  }
825 
826  // Recover id of the cell and points on which compute the pseudonormal
827  for (int irank = 0; irank < nProcessors; irank++){
828  if (myrank == irank){
829  for (int jrank = 0; jrank < nProcessors; jrank++){
830  if (jrank != irank){
831  int n_send_to_current_rank = n_send_to_rank[jrank];
832  MPI_Recv(send_to_rank[jrank].data(), n_send_to_current_rank, MPI_LONG, jrank, 100, tree->getPatch().getCommunicator(), MPI_STATUS_IGNORE);
833  MPI_Recv(point_from_rank[jrank].data(), n_send_to_current_rank*3, MPI_DOUBLE, jrank, 101, tree->getPatch().getCommunicator(), MPI_STATUS_IGNORE);
834  }
835  }
836  } else {
837  int n_ask_to_current_rank = ask_to_rank[irank].size();
838  MPI_Send(ask_to_rank[irank].data(), n_ask_to_current_rank, MPI_LONG, irank, 100, tree->getPatch().getCommunicator());
839  MPI_Send(point_to_rank[irank].data(), n_ask_to_current_rank*3, MPI_DOUBLE, irank, 101, tree->getPatch().getCommunicator());
840  }
841  }
842 
843  // Compute pseudo-normal for each received point from each rank
844  std::map<int, std::vector<std::array<double,3>>> normal_to_rank;
845  std::map<int, std::vector<double>> sign_to_rank;
846 
847  for (int irank = 0; irank < nProcessors; irank++){
848  normal_to_rank[irank].resize(n_send_to_rank[irank], std::array<double,3>({std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}));
849  sign_to_rank[irank].resize(n_send_to_rank[irank], std::numeric_limits<double>::max());
850  for (int ip = 0; ip < n_send_to_rank[irank]; ip++){
851 
852  const std::array<double,3> & point = point_from_rank[irank][ip];
853  const long & cellId = send_to_rank[irank][ip];
854  std::array<double,3> & pseudo_normal = normal_to_rank[irank][ip];
855  double & s = sign_to_rank[irank][ip];
856 
857  s = computePseudoNormal(point, &spatch, cellId, pseudo_normal);
858 
859  } // end loop on points received
860  } // end loop on ranks
861 
862  // Communicate back the computed pseudonormals
863  std::map<int, std::vector<std::array<double,3>>> normal_from_rank;
864  std::map<int, std::vector<double>> sign_from_rank;
865 
866  for (int irank = 0; irank < nProcessors; irank++){
867  if (myrank == irank){
868  for (int jrank = 0; jrank < nProcessors; jrank++){
869  if (jrank != irank){
870  int n_ask_to_current_rank = ask_to_rank[jrank].size();
871  normal_from_rank[jrank].resize(n_ask_to_current_rank);
872  MPI_Recv(normal_from_rank[jrank].data(), n_ask_to_current_rank*3, MPI_DOUBLE, jrank, 102, tree->getPatch().getCommunicator(), MPI_STATUS_IGNORE);
873 
874  sign_from_rank[jrank].resize(n_ask_to_current_rank);
875  MPI_Recv(sign_from_rank[jrank].data(), n_ask_to_current_rank, MPI_DOUBLE, jrank, 103, tree->getPatch().getCommunicator(), MPI_STATUS_IGNORE);
876  // Fill normals output with the correct received pseudonormals (TODO move out of communications scope?)
877 
878  for (int ip = 0; ip < n_ask_to_current_rank; ip++){
879  int point_index = point_index_received_from_rank[jrank][ip];
880  normals[point_index] = normal_from_rank[jrank][ip];
881  distances[point_index] *= sign_from_rank[jrank][ip];
882  }
883 
884  }
885  }
886  } else {
887  int n_send_to_current_rank = n_send_to_rank[irank];
888  MPI_Send(normal_to_rank[irank].data(), n_send_to_current_rank*3, MPI_DOUBLE, irank, 102, tree->getPatch().getCommunicator());
889  MPI_Send(sign_to_rank[irank].data(), n_send_to_current_rank, MPI_DOUBLE, irank, 103, tree->getPatch().getCommunicator());
890  }
891  }
892 
893  } // end if not shared points
894  }
895 }
896 
910 void projectPointGlobal(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, std::array<double,3> *projected_points, long *ids, int *ranks, double r, bool shared )
911 {
912  std::vector<double> rs(nP, r);
913  projectPointGlobal(nP, points, tree, projected_points, ids, ranks, rs.data(), shared);
914 }
915 
930 void projectPointGlobal(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, std::array<double,3> *projected_points, long *ids, int *ranks, double *r, bool shared )
931 {
932  // Initialize ids and ranks
933  for (int ip = 0; ip < nP; ip++){
934  ids[ip] = bitpit::Cell::NULL_ID;
935  ranks[ip] = -1;
936  r[ip] = std::max(r[ip], tree->getPatch().getTol());
937  }
938 
939  bool checkEmptyList = (nP == 0);
940  MPI_Allreduce(MPI_IN_PLACE, &checkEmptyList, 1, MPI_C_BOOL, MPI_LAND, tree->getPatch().getCommunicator());
941 
942  std::vector<darray3E> workpoints(points, points+nP);
943  std::vector<long> workids(ids, ids+nP);
944  std::vector<int> workranks(ranks, ranks+nP);
945  std::vector<double> workradius(r, r+nP);
946  std::unordered_map<std::size_t,std::size_t> mapPosIndex;
947  for(std::size_t i=0; i<std::size_t(nP); ++i){
948  mapPosIndex[i] = i;
949  }
950  std::vector<darray3E> normals(nP), workNormals(nP);
951  std::vector<double> dist(nP, std::numeric_limits<double>::max()), workDist(nP);
952 
953  int kmax(1000);
954  int kiter(0);
955 
956 
957  while (!checkEmptyList && kiter < kmax){
958  //use method sphere by default
959  signedGlobalDistance(workpoints.size(), workpoints.data(), tree, workids.data(), workranks.data(), workNormals.data(), workDist.data(), workradius.data(), shared);
960 
961  //get all points with distances not calculated.
962  std::vector<std::array<double,3>> failedPoints;
963  std::vector<double> failedRadii;
964  std::unordered_map<std::size_t, std::size_t> failedPosIndex;
965  failedPoints.reserve(workpoints.size());
966  failedRadii.reserve(workpoints.size());
967  int count(0);
968  for(long & val : workids){
969  if(val == bitpit::Cell::NULL_ID){
970  failedPoints.push_back(workpoints[count]);
971  failedPosIndex[failedPoints.size()-1] = mapPosIndex[count];
972  failedRadii.push_back(workradius[count]);
973  }else{
974  dist[mapPosIndex[count]] = workDist[count];
975  normals[mapPosIndex[count]] = workNormals[count];
976  ids[mapPosIndex[count]] = workids[count];
977  ranks[mapPosIndex[count]] = workranks[count];
978  r[mapPosIndex[count]] = workradius[count];
979  }
980  ++count;
981  }
982  workDist.resize(failedPoints.size());
983  workNormals.resize(failedPoints.size());
984  workids.resize(failedPoints.size());
985  workranks.resize(failedPoints.size());
986 
987  std::swap(failedPoints, workpoints);
988  std::swap(failedPosIndex, mapPosIndex);
989  std::swap(failedRadii, workradius);
990 
991 
992  checkEmptyList = workpoints.empty();
993  MPI_Allreduce(MPI_IN_PLACE, &checkEmptyList, 1, MPI_C_BOOL, MPI_LAND, tree->getPatch().getCommunicator());
994 
995  //augment workradius of 1.5 factor
996  if(kiter >= kmax-2){
997  std::vector<double> temp(workradius.size(), std::numeric_limits<double>::max());
998  std::swap(temp, workradius);
999  }else{
1000  workradius *= 1.5;
1001  }
1002  // increment iteration
1003  ++kiter;
1004  }//end while.
1005 
1006  for (int ip = 0; ip < nP; ip++){
1007  projected_points[ip] = points[ip] - dist[ip] * normals[ip];
1008  }
1009 }
1010 
1021 void locatePointOnGlobalPatch(int nP, const std::array<double,3> *points, const bitpit::PatchSkdTree *tree, long *ids, int *ranks, bool shared)
1022 {
1023  if(!dynamic_cast<const bitpit::SurfUnstructured*>(&(tree->getPatch()))){
1024  throw std::runtime_error("Invalid use of skdTreeUtils::locatePointOnPatch method: a non surface patch or void patch was detected.");
1025  }
1026  const bitpit::SurfUnstructured & spatch = static_cast<const bitpit::SurfUnstructured&>(tree->getPatch());
1027 
1028  // Initialize the cell id and distance
1029  for (int ip = 0; ip < nP; ip++){
1030  ids[ip] = bitpit::Cell::NULL_ID;
1031  }
1032  std::vector<double> distances(nP, std::numeric_limits<double>::max());
1033 
1034  // Find the closest cell to the input point
1035  std::vector<long> cellIds(nP, bitpit::Cell::NULL_ID);
1036  std::vector<int> cellRanks(nP, -1);
1037  if (shared){
1038  // All processes share the points
1039  findSharedPointClosestGlobalCell(nP, points, tree, cellIds.data(), ranks, distances.data());
1040  } else {
1041  // Distributed points
1042  static_cast<const bitpit::SurfaceSkdTree*>(tree)->findPointClosestGlobalCell(nP, points, cellIds.data(), ranks, distances.data());
1043  }
1044 
1045  int myrank = spatch.getRank();
1046 
1047  // Map with rank to ask info in case of distributed points
1048  std::map<int,std::vector<long>> ask_to_rank;
1049  std::map<int,std::vector<std::array<double,3>>> point_to_rank;
1050  std::map<int,std::vector<int>> point_index_received_from_rank;
1051 
1052  // Loop on points
1053  for (int ip = 0; ip < nP; ip++){
1054 
1055  const std::array<double,3> & point = points[ip];
1056  const long & cellId = cellIds[ip];
1057  const int & cellRank = ranks[ip];
1058 
1059  // If the current rank is the the owner of the cell check if the point belongs to the closest cell
1060  if (cellRank == myrank){
1061 
1062  bool checkBelong = checkPointBelongsToCell(point, &spatch, cellId);
1063 
1064  if (checkBelong) {
1065  ids[ip] = cellId;
1066  }
1067 
1068  } else {
1069 
1070  // If not owner maintain the Cell::NULL_ID value, equal to minimum long int for bitpit::Cell
1071 
1072  if (!shared){
1073 
1074  // Push to rank map to ask normal
1075  ask_to_rank[cellRank].push_back(cellId);
1076 
1077  // Push to rank map the related points to ask normal
1078  point_to_rank[cellRank].push_back(point);
1079 
1080  // Push to rank map the index of the points to fill output normals
1081  point_index_received_from_rank[cellRank].push_back(ip);
1082 
1083  }
1084 
1085  } // end if owner rank
1086 
1087  } // End loop on points
1088 
1089  int nProcessors = spatch.getProcessorCount();
1090 
1091  if (nProcessors > 1 && spatch.isPartitioned()){
1092  if (shared){
1093 
1094  // If all points are shared between processes all reduce on ids by maximum operation
1095  MPI_Allreduce(MPI_IN_PLACE, ids, nP, MPI_LONG, MPI_MAX, tree->getPatch().getCommunicator());
1096 
1097  } else {
1098 
1099  // Send/receive number of check belong to send to other processes
1100  std::vector<int> n_send_to_rank(nProcessors, 0);
1101  for (int irank = 0; irank < nProcessors; irank++){
1102  int n_ask_to_current_rank = ask_to_rank[irank].size();
1103  MPI_Gather(&n_ask_to_current_rank, 1, MPI_INT, n_send_to_rank.data(), 1, MPI_INT, irank, tree->getPatch().getCommunicator());
1104  }
1105 
1106  // Map with rank to send info for cellId and point coordinates
1107  std::map<int,std::vector<long>> send_to_rank;
1108  std::map<int,std::vector<std::array<double,3>>> point_from_rank;
1109  for (int irank = 0; irank < nProcessors; irank++){
1110  send_to_rank[irank].resize(n_send_to_rank[irank], bitpit::Cell::NULL_ID);
1111  point_from_rank[irank].resize(n_send_to_rank[irank], std::array<double,3>({{0.,0.,0.}}));
1112  }
1113 
1114  // Recover id of the cell and points on which check the belonging
1115  for (int irank = 0; irank < nProcessors; irank++){
1116  if (myrank == irank){
1117  for (int jrank = 0; jrank < nProcessors; jrank++){
1118  if (jrank != irank){
1119  int n_send_to_current_rank = n_send_to_rank[jrank];
1120  MPI_Recv(send_to_rank[jrank].data(), n_send_to_current_rank, MPI_LONG, jrank, 100, tree->getPatch().getCommunicator(), MPI_STATUS_IGNORE);
1121  MPI_Recv(point_from_rank[jrank].data(), n_send_to_current_rank*3, MPI_DOUBLE, jrank, 101, tree->getPatch().getCommunicator(), MPI_STATUS_IGNORE);
1122  }
1123  }
1124  } else {
1125  int n_ask_to_current_rank = ask_to_rank[irank].size();
1126  MPI_Send(ask_to_rank[irank].data(), n_ask_to_current_rank, MPI_LONG, irank, 100, tree->getPatch().getCommunicator());
1127  MPI_Send(point_to_rank[irank].data(), n_ask_to_current_rank*3, MPI_DOUBLE, irank, 101, tree->getPatch().getCommunicator());
1128  }
1129  }
1130 
1131  // Check if the point belong to cell for each received point from each rank
1132  std::map<int, std::vector<int>> checkBelong_to_rank;
1133  for (int irank = 0; irank < nProcessors; irank++){
1134  checkBelong_to_rank[irank].resize(n_send_to_rank[irank], false);
1135  for (int ip = 0; ip < n_send_to_rank[irank]; ip++){
1136 
1137  const std::array<double,3> & point = point_from_rank[irank][ip];
1138  const long & cellId = send_to_rank[irank][ip];
1139  int & checkBelong = checkBelong_to_rank[irank][ip];
1140 
1141  checkBelong = int(checkPointBelongsToCell(point, &spatch, cellId));
1142 
1143  } // end loop on points received
1144  } // end loop on ranks
1145 
1146  // Communicate back the computed cehck belong
1147  std::map<int, std::vector<int>> checkBelong_from_rank;
1148 
1149  for (int irank = 0; irank < nProcessors; irank++){
1150  if (myrank == irank){
1151  for (int jrank = 0; jrank < nProcessors; jrank++){
1152  if (jrank != irank){
1153  int n_ask_to_current_rank = ask_to_rank[jrank].size();
1154  checkBelong_from_rank[jrank].resize(n_ask_to_current_rank);
1155  MPI_Recv(checkBelong_from_rank[jrank].data(), n_ask_to_current_rank, MPI_INT, jrank, 102, tree->getPatch().getCommunicator(), MPI_STATUS_IGNORE);
1156 
1157  // Fill cell ids output if the received check belong is true (TODO move out of communications scope?)
1158  for (int ip = 0; ip < n_ask_to_current_rank; ip++){
1159  int point_index = point_index_received_from_rank[jrank][ip];
1160  long cellId = ask_to_rank[jrank][ip];
1161  // Check received int value (zero or one) casting to boolean
1162  if (bool(checkBelong_from_rank[jrank][ip])) {
1163  ids[point_index] = cellId;
1164  }
1165  }
1166 
1167  }
1168  }
1169  } else {
1170  int n_send_to_current_rank = n_send_to_rank[irank];
1171  MPI_Send(checkBelong_to_rank[irank].data(), n_send_to_current_rank, MPI_INT, irank, 102, tree->getPatch().getCommunicator());
1172  }
1173  }
1174 
1175  } // end if not shared points
1176  }
1177 }
1178 
1190 std::vector<long> selectByGlobalPatch(bitpit::PatchSkdTree *selection, bitpit::PatchSkdTree *target, double tol){
1191 
1192  if (!selection->getPatch().isPartitioned()){
1193  throw std::runtime_error("Error: selection PatchSkdTree communicators not set in selectByGlobalPatch");
1194  }
1195  if (!target->getPatch().isPartitioned()){
1196  throw std::runtime_error("Error: target PatchSkdTree communicators not set in selectByGlobalPatch");
1197  }
1198 
1199  //TODO USE THE REAL IS DISTRIBUTED PATCH IN BITPIT WHEN READY !!!!
1200 // if (!selection->getPatch().isPartitioned() && !target->getPatch().isPartitioned()){
1201 // return(selectByPatch(selection, target, tol));
1202 // }
1203 
1204  // Leaf selection nodes of current rank initialized for each rank of target bounding box
1205  int nprocs = selection->getPatch().getProcessorCount();
1206  int myRank = selection->getPatch().getRank();
1207  int nleafs = selection->getLeafCount();
1208  int nnodess = selection->getNodeCount();
1209  std::unordered_map<int, std::vector<const bitpit::SkdNode*>> rankLeafSelection;
1210  for (int irank = 0; irank < nprocs; irank++){
1211  std::vector<const bitpit::SkdNode*> leafSelection(nleafs);
1212  int count = 0;
1213  for (int i=0; i<nnodess; i++){
1214  if (selection->getNode(i).isLeaf()){
1215  if (bitpit::CGElem::intersectBoxBox(selection->getNode(i).getBoxMin(),
1216  selection->getNode(i).getBoxMax(),
1217  target->getPartitionBox(irank).getBoxMin(),
1218  target->getPartitionBox(irank).getBoxMax(),
1219  3,
1220  tol ) ){
1221  leafSelection[count] = &(selection->getNode(i));
1222  count++;
1223  }
1224  }
1225  }
1226  leafSelection.resize(count);
1227  rankLeafSelection[irank].swap(leafSelection);
1228  }
1229 
1230  // Collect all the leaf selection SkdBox bounding boxes from each process in the map structure
1231  // Collect and store from other processes the leaf selection bounding boxes for the current target rank
1232  std::vector<bitpit::SkdBox> leafSelectionBoxes;
1233 
1234  // Recover the number of selection elements of each process for each target process
1235  // This is the number of elements sent and received to/from each process from/to each process
1236  std::vector<int> nLeafSend(nprocs, 0);
1237  std::vector<int> nLeafRecv(nprocs, 0);
1238  for (int irank = 0; irank < nprocs; irank++){
1239  nLeafSend[irank] = rankLeafSelection[irank].size();
1240  }
1241 
1242  for (int irank = 0; irank < nprocs; irank++){
1243  // Scatter send to all process of irank process data
1244  int * recvbuff = &nLeafRecv[irank];
1245  MPI_Scatter(nLeafSend.data(), 1, MPI_INT, recvbuff, 1, MPI_INT, irank, selection->getPatch().getCommunicator());
1246  }
1247 
1248  // Recover global leaf received to resize local received boxes container
1249  int nGlobalReceived = 0;
1250  for (int val : nLeafRecv){
1251  nGlobalReceived += val;
1252  }
1253 
1254  if (nGlobalReceived > 0){
1255  leafSelectionBoxes.reserve(nGlobalReceived);
1256  }
1257 
1258  // Communicate minimum and maximum point of each leaf bounding box to the right process
1259  // Gather bounding boxes for each process by communicating minimum and maximum point
1260  // coordinates as vector of 6 coordinates for each box
1261  std::vector<double> recvBoxes;
1262  for (int irank = 0; irank < nprocs; irank++){
1263  int nSendData = 0;
1264  std::vector<double> sendBoxes;
1265  int nRecvData = 0;
1266  std::vector<int> recvDispls(nprocs, 0);
1267  std::vector<int> nRecvDataPerProc(nprocs, 0);
1268  if (irank == myRank){
1269  // If root is current rank set receive buffer
1270  nRecvData = 6 * nGlobalReceived;
1271  recvBoxes.resize(nRecvData, std::numeric_limits<double>::max());
1272  // Set receiving displacements to gather bounding boxes
1273  nRecvDataPerProc[0] = nLeafRecv[0] * 6;
1274  for (int jrank = 1; jrank < nprocs; jrank++){
1275  recvDispls[jrank] = recvDispls[jrank - 1] + nLeafRecv[jrank - 1] * 6;
1276  nRecvDataPerProc[jrank] = nLeafRecv[jrank] * 6;
1277  }
1278  }
1279  // Set send buffer for all processes even for root
1280  nSendData = 6 * nLeafSend[irank];
1281  sendBoxes.reserve(nSendData);
1282  for (const bitpit::SkdNode* selectNode : rankLeafSelection[irank]){
1283  for (const double val : selectNode->getBoxMin()){
1284  sendBoxes.push_back(val);
1285  }
1286  for (const double val : selectNode->getBoxMax()){
1287  sendBoxes.push_back(val);
1288  }
1289  }
1290 
1291  // Gather bounding boxes points to the root process
1292  MPI_Gatherv(sendBoxes.data(), nSendData, MPI_DOUBLE, recvBoxes.data(), nRecvDataPerProc.data(), recvDispls.data(), MPI_DOUBLE, irank, selection->getPatch().getCommunicator());
1293 
1294  } // End loop on root processes
1295 
1296  // Fill leaf selection boxes container
1297  std::vector<double>::iterator itBox = recvBoxes.begin();
1298  for (int ibox = 0; ibox < nGlobalReceived; ibox++){
1299  std::array<double, 3> minPoint;
1300  minPoint[0] = *(itBox++);
1301  minPoint[1] = *(itBox++);
1302  minPoint[2] = *(itBox++);
1303  std::array<double, 3> maxPoint;
1304  maxPoint[0] = *(itBox++);
1305  maxPoint[1] = *(itBox++);
1306  maxPoint[2] = *(itBox++);
1307  leafSelectionBoxes.emplace_back(bitpit::SkdBox(minPoint, maxPoint));
1308  }
1309 
1310  std::vector<long> extracted;
1311  extractTarget(target, leafSelectionBoxes, extracted, tol);
1312 
1313  return extracted;
1314 
1315 }
1316 
1331 void extractTarget(bitpit::PatchSkdTree *target, const std::vector<bitpit::SkdBox> &leafSelectionBoxes, std::vector<long> &extracted, double tol)
1332 {
1333 
1334  if(leafSelectionBoxes.empty()) return;
1335  std::size_t rootId = 0;
1336 
1337  // Candidates saved in set to avoid duplicate
1338  std::unordered_set<std::size_t> candidates;
1339 
1340  // Loop on leaf selection boxes and go across the target tree with a node stack
1341  for (const bitpit::SkdBox & selectionBox : leafSelectionBoxes){
1342 
1343  // Initialize stack with target root node
1344  std::queue<std::size_t> nodeStack;
1345  nodeStack.push(rootId);
1346 
1347  // Pop and push target tree nodes checked with current selection boxe
1348  // Fill candidates structure with leaf target nodes intersected by selection box
1349  // Stop when stack is empty
1350  while(!nodeStack.empty()){
1351 
1352  std::size_t nodeId = nodeStack.front();
1353  nodeStack.pop();
1354  const bitpit::SkdNode & node = target->getNode(nodeId);
1355 
1356  bool isLeaf = true;
1357  for (int i = bitpit::SkdNode::CHILD_BEGIN; i != bitpit::SkdNode::CHILD_END; ++i) {
1358  bitpit::SkdNode::ChildLocation childLocation = static_cast<bitpit::SkdNode::ChildLocation>(i);
1359  std::size_t childId = node.getChildId(childLocation);
1360  if (childId != bitpit::SkdNode::NULL_ID) {
1361  isLeaf = false;
1362  if (bitpit::CGElem::intersectBoxBox(selectionBox.getBoxMin(), selectionBox.getBoxMax(),
1363  target->getNode(childId).getBoxMin(), target->getNode(childId).getBoxMax(),3, tol)){
1364  nodeStack.push(childId);
1365  }
1366  }
1367  }
1368  if (isLeaf) {
1369  candidates.insert(nodeId);
1370  }
1371  }
1372  } // End loop on selection boxes
1373 
1374  // Extract cell from candidates nodes
1375  // Do not check for intersection -> exctraction performed only by using bounding boxes of leaf nodes of the skd-trees
1376  std::unordered_set<long> cellExtracted;
1377  for(long nodeId : candidates){
1378  const bitpit::SkdNode & node = target->getNode(nodeId);
1379  std::vector<long> cellids = node.getCells();
1380  cellExtracted.insert(cellids.begin(), cellids.end());
1381  }
1382 
1383  extracted.insert(extracted.end(), cellExtracted.begin(), cellExtracted.end());
1384 }
1385 
1386 
1407 long findSharedPointClosestGlobalCell(int nPoints, const std::array<double, 3> *points, const bitpit::PatchSkdTree *tree,
1408  long *ids, int *ranks, double *distances, double r)
1409 {
1410  std::vector<double> rs(nPoints, r);
1411  return findSharedPointClosestGlobalCell(nPoints, points, tree, ids, ranks, distances, rs.data());
1412 }
1413 
1434 long findSharedPointClosestGlobalCell(int nPoints, const std::array<double, 3> *points, const bitpit::PatchSkdTree *tree,
1435  long *ids, int *ranks, double *distances, double *r)
1436 {
1437  //TODO evaluate if to use a provided r or always a numeric_limits::max here as maxDistance. Clean it up useless stuffs then.
1438  BITPIT_UNUSED(r); // just to suppress warnings.
1439 
1440  long nDistanceEvaluations = 0;
1441  std::vector<double> maxDistances(nPoints, std::numeric_limits<double>::max());
1442 
1443  // Early return is the patch is not partitioned
1444  const bitpit::PatchKernel &patch = tree->getPatch();
1445  if (!patch.isPartitioned()) {
1446  for (int i = 0; i < nPoints; ++i) {
1447  // Evaluate distance
1448  nDistanceEvaluations += findPointClosestCell(points[i], tree, maxDistances[i], ids + i, distances + i);
1449 
1450  // The patch is not partitioned, all cells are local
1451  ranks[i] = patch.getRank();
1452 
1453  }
1454 
1455  return nDistanceEvaluations;
1456 
1457  }
1458 
1459  // Get MPI communicator
1460  // Patch is partitioned call the parallel method
1461  if (!tree->getPatch().isPartitioned()){
1462  throw std::runtime_error ("There is no communicator set for the patch.");
1463  }
1464 
1465  int nprocs = tree->getPatch().getProcessorCount();
1466 
1467  // Initialize distance information
1468  std::vector<bitpit::SkdGlobalCellDistance> globalCellDistances(nPoints);
1469 
1470  // Call local find point closest cell for each point
1471  for (int i = 0; i < nPoints; ++i) {
1472  // Get point information
1473  const std::array<double, 3> &point = points[i];
1474 
1475  // Use a maximum distance for each point given by an estimation
1476  // based on partition bounding boxes. The distance will be lesser
1477  // than or equal to the point maximum distance.
1478  double pointMaxDistance = maxDistances[i];
1479  for (int rank = 0; rank < nprocs; ++rank) {
1480  pointMaxDistance = std::min(tree->getPartitionBox(rank).evalPointMaxDistance(point), pointMaxDistance);
1481  }
1482 
1483  // Get cell distance information
1484  bitpit::SkdGlobalCellDistance &globalCellDistance = globalCellDistances[i];
1485  int &cellRank = globalCellDistance.getRank();
1486  long &cellId = globalCellDistance.getId();
1487  double &cellDistance = globalCellDistance.getDistance();
1488 
1489  // Evaluate local distance from the point
1490  bool interiorOnly = true;
1491  nDistanceEvaluations += findPointClosestCell(point, tree, pointMaxDistance, interiorOnly, &cellId, &cellDistance);
1492 
1493  // Set cell rank
1494  if (cellId != bitpit::Cell::NULL_ID) {
1495  cellRank = patch.getCellRank(cellId);
1496  }
1497  }
1498 
1499  // Exchange distance information
1500  // Communicate the closest cells distances, ranks and ids by reduce with custom minimum distance operation
1501  // Reduce only the right portion of data to the right process
1502  MPI_Datatype globalCellDistanceDatatype = bitpit::SkdGlobalCellDistance::getMPIDatatype();
1503  MPI_Op globalCellDistanceMinOp = bitpit::SkdGlobalCellDistance::getMPIMinOperation();
1504  bitpit::SkdGlobalCellDistance *globalCellDistance = globalCellDistances.data();
1505  for (int irank = 0; irank < nprocs; irank++){
1506  MPI_Allreduce(MPI_IN_PLACE, globalCellDistance, nPoints, globalCellDistanceDatatype, globalCellDistanceMinOp, tree->getPatch().getCommunicator());
1507  }
1508 
1509  // Update output arguments
1510  for (int i = 0; i < nPoints; ++i) {
1511  int globalIndex = i;
1512  bitpit::SkdGlobalCellDistance &globalCellDistance = globalCellDistances[globalIndex];
1513  globalCellDistance.exportData(ranks + i, ids + i, distances + i);
1514  }
1515 
1516  return nDistanceEvaluations;
1517 }
1518 #endif
1519 
1538 long findPointClosestCell(const std::array<double, 3> &point, const bitpit::PatchSkdTree *tree, bool interiorOnly, long *id, double *distance)
1539 {
1540  return findPointClosestCell(point, tree, std::numeric_limits<double>::max(), interiorOnly, id, distance);
1541 
1542 }
1543 
1564 long findPointClosestCell(const std::array<double, 3> &point, const bitpit::PatchSkdTree *tree,
1565  double maxDistance, bool interiorOnly, long *id, double *distance)
1566 {
1567 
1568  // Initialize the cell id
1569  *id = bitpit::Cell::NULL_ID;
1570 
1571  // Initialize the distance with an estimate
1572  //
1573  // The real distance will be lesser than or equal to the estimate.
1574  std::size_t rootId = 0;
1575  const bitpit::SkdNode &root = tree->getNode(rootId);
1576  *distance = std::min(root.evalPointMaxDistance(point), maxDistance);
1577 
1578  // Get a list of candidates nodes
1579  //
1580  // Initialize list of candidates
1581  std::vector<std::size_t> candidateIds;
1582  std::vector<double> candidateMinDistances;
1583 
1584  std::vector<std::size_t> nodeStack;
1585  nodeStack.push_back(rootId);
1586  while (!nodeStack.empty()) {
1587  std::size_t nodeId = nodeStack.back();
1588  const bitpit::SkdNode &node = tree->getNode(nodeId);
1589  nodeStack.pop_back();
1590 
1591  // Do not consider nodes with a minimum distance greater than
1592  // the distance estimate
1593  double nodeMinDistance = node.evalPointMinDistance(point);
1594  if (nodeMinDistance > *distance) {
1595  continue;
1596  }
1597 
1598  // Update the distance estimate
1599  //
1600  // The real distance will be lesser than or equal to the
1601  // estimate.
1602  double nodeMaxDistance = node.evalPointMaxDistance(point);
1603  *distance = std::min(nodeMaxDistance, *distance);
1604 
1605  // If the node is a leaf add it to the candidates, otherwise
1606  // add its children to the stack.
1607  bool isLeaf = true;
1608  for (int i = bitpit::SkdNode::CHILD_BEGIN; i != bitpit::SkdNode::CHILD_END; ++i) {
1609  bitpit::SkdNode::ChildLocation childLocation = static_cast<bitpit::SkdNode::ChildLocation>(i);
1610  std::size_t childId = node.getChildId(childLocation);
1611  if (childId != bitpit::SkdNode::NULL_ID) {
1612  isLeaf = false;
1613  nodeStack.push_back(childId);
1614  }
1615  }
1616 
1617  if (isLeaf) {
1618  candidateIds.push_back(nodeId);
1619  candidateMinDistances.push_back(nodeMinDistance);
1620  }
1621  }
1622 
1623  // Process the candidates and find the closest cell
1624  long nDistanceEvaluations = 0;
1625  for (std::size_t k = 0; k < candidateIds.size(); ++k) {
1626  // Do not consider nodes with a minimum distance greater than
1627  // the distance estimate
1628  if (candidateMinDistances[k] > *distance) {
1629  continue;
1630  }
1631 
1632  // Evaluate the distance
1633  std::size_t nodeId = candidateIds[k];
1634  const bitpit::SkdNode &node = tree->getNode(nodeId);
1635 
1636  node.updatePointClosestCell(point, interiorOnly, id, distance);
1637  ++nDistanceEvaluations;
1638  }
1639 
1640  // If no closest cell was found set the distance to the maximum
1641  // representable distance.
1642  if (*id == bitpit::Cell::NULL_ID) {
1643  *distance = std::numeric_limits<double>::max();
1644  }
1645 
1646  return nDistanceEvaluations;
1647 }
1648 
1649 #if MIMMO_ENABLE_MPI
1650 
1664 long findPointClosestGlobalCell(int nPoints, const std::array<double, 3> *points,
1665  const bitpit::PatchSkdTree *tree, long *ids, int *ranks, double *distances)
1666 {
1667 
1668  long nDistanceEvaluations = 0;
1669 
1670  std::vector<double> maxDistances(nPoints, std::numeric_limits<double>::max());
1671 
1672  // Early return is the patch is not partitioned
1673  const bitpit::PatchKernel &patch = tree->getPatch();
1674  if (!patch.isPartitioned()) {
1675  for (int i = 0; i < nPoints; ++i) {
1676  // Evaluate distance
1677  nDistanceEvaluations += findPointClosestCell(points[i], tree, maxDistances[i], ids + i, distances + i);
1678 
1679  // The patch is not partitioned, all cells are local
1680  ranks[i] = patch.getRank();
1681 
1682  }
1683 
1684  return nDistanceEvaluations;
1685 
1686  }
1687 
1688  // Get MPI communicator
1689  // Patch is partitioned call the parallel method
1690  if (!tree->getPatch().isPartitioned()){
1691  throw std::runtime_error ("There is no communicator set for the patch.");
1692  }
1693 
1694  MPI_Comm communicator = tree->getPatch().getCommunicator();
1695  int nprocs = tree->getPatch().getProcessorCount();
1696  int myRank = tree->getPatch().getRank();
1697 
1698  // Gather the number of points associated to each process
1699  std::vector<int> pointsCount(nprocs);
1700  MPI_Allgather(&nPoints, 1, MPI_INT, pointsCount.data(), 1, MPI_INT, communicator);
1701 
1702  // Evaluate information for data communications
1703  std::vector<int> globalPointsDispls(nprocs, 0);
1704  std::vector<int> globalPointsOffsets(nprocs, 0);
1705  std::vector<int> globalPointsDataCount(nprocs, 0);
1706 
1707  globalPointsDataCount[0] = 3 * pointsCount[0];
1708  for (int i = 1; i < nprocs; ++i) {
1709  globalPointsDispls[i] = globalPointsDispls[i - 1] + 3 * pointsCount[i - 1];
1710  globalPointsOffsets[i] = globalPointsOffsets[i - 1] + pointsCount[i - 1];
1711  globalPointsDataCount[i] = 3 * pointsCount[i];
1712  }
1713 
1714  int nGlobalPoints = globalPointsDispls.back() + pointsCount.back();
1715 
1716  // Gather point coordinates
1717  std::vector<std::array<double,3>> globalPoints(nGlobalPoints);
1718  int pointsDataCount = 3 * nPoints;
1719  MPI_Allgatherv(points, pointsDataCount, MPI_DOUBLE, globalPoints.data(),
1720  globalPointsDataCount.data(), globalPointsDispls.data(), MPI_DOUBLE, communicator);
1721 
1722  // Gather vector with all maximum distances
1723  std::vector<double> globalMaxDistances(nGlobalPoints);
1724  MPI_Allgatherv(maxDistances.data(), nPoints, MPI_DOUBLE, globalMaxDistances.data(),
1725  pointsCount.data(), globalPointsOffsets.data(), MPI_DOUBLE, communicator);
1726 
1727  // Initialize distance information
1728  std::vector<bitpit::SkdGlobalCellDistance> globalCellDistances(nGlobalPoints);
1729 
1730  // Call local find point closest cell for each global point collected
1731  for (int i = 0; i < nGlobalPoints; ++i) {
1732  // Get point information
1733  const std::array<double, 3> &point = globalPoints[i];
1734 
1735  // Use a maximum distance for each point given by an estimation
1736  // based on partition bounding boxes. The distance will be lesser
1737  // than or equal to the point maximum distance.
1738  double pointMaxDistance = globalMaxDistances[i];
1739  for (int rank = 0; rank < nprocs; ++rank) {
1740  pointMaxDistance = std::min(tree->getPartitionBox(rank).evalPointMaxDistance(point), pointMaxDistance);
1741  }
1742 
1743  // Get cell distance information
1744  bitpit::SkdGlobalCellDistance &globalCellDistance = globalCellDistances[i];
1745  int &cellRank = globalCellDistance.getRank();
1746  long &cellId = globalCellDistance.getId();
1747  double &cellDistance = globalCellDistance.getDistance();
1748 
1749  // Evaluate local distance from the point
1750  bool interiorOnly = true;
1751  nDistanceEvaluations += findPointClosestCell(point, tree, pointMaxDistance, interiorOnly, &cellId, &cellDistance);
1752 
1753  // Set cell rank
1754  if (cellId != bitpit::Cell::NULL_ID) {
1755  cellRank = patch.getCellRank(cellId);
1756  }
1757  }
1758 
1759  // Exchange distance information
1760  MPI_Datatype globalCellDistanceDatatype = bitpit::SkdGlobalCellDistance::getMPIDatatype();
1761  MPI_Op globalCellDistanceMinOp = bitpit::SkdGlobalCellDistance::getMPIMinOperation();
1762  for (int rank = 0; rank < nprocs; ++rank) {
1763  bitpit::SkdGlobalCellDistance *globalCellDistance = globalCellDistances.data() + globalPointsOffsets[rank];
1764  if (myRank == rank) {
1765  MPI_Reduce(MPI_IN_PLACE, globalCellDistance, pointsCount[rank], globalCellDistanceDatatype, globalCellDistanceMinOp, rank, communicator);
1766  } else {
1767  MPI_Reduce(globalCellDistance, globalCellDistance, pointsCount[rank], globalCellDistanceDatatype, globalCellDistanceMinOp, rank, communicator);
1768  }
1769  }
1770 
1771  // Update output arguments
1772  for (int i = 0; i < nPoints; ++i) {
1773  int globalIndex = i + globalPointsOffsets[myRank];
1774  bitpit::SkdGlobalCellDistance &globalCellDistance = globalCellDistances[globalIndex];
1775  globalCellDistance.exportData(ranks + i, ids + i, distances + i);
1776  }
1777 
1778  return nDistanceEvaluations;
1779 }
1780 #endif
1781 }
1782 
1783 }; // end namespace mimmo
mimmo::skdTreeUtils::findPointClosestCell
long findPointClosestCell(const std::array< double, 3 > &point, const bitpit::PatchSkdTree *tree, bool interiorOnly, long *id, double *distance)
Definition: SkdTreeUtils.cpp:1538
mimmo::skdTreeUtils::signedDistance
double signedDistance(const std::array< double, 3 > *point, const bitpit::PatchSkdTree *tree, long &id, std::array< double, 3 > &normal, double r)
Definition: SkdTreeUtils.cpp:84
mimmo::skdTreeUtils::computePseudoNormal
double computePseudoNormal(const std::array< double, 3 > &point, const bitpit::SurfUnstructured *surface_mesh, long id, std::array< double, 3 > &pseudo_normal)
Definition: SkdTreeUtils.cpp:517
dvecarr3E
std::vector< darray3E > dvecarr3E
Definition: mimmoTypeDef.hpp:58
mimmo::skdTreeUtils::checkPointBelongsToCell
bool checkPointBelongsToCell(const std::array< double, 3 > &point, const bitpit::SurfUnstructured *surface_mesh, long cellId)
Definition: SkdTreeUtils.cpp:592
mimmo::mimmoCGUtils::isPointInsideSegment
bool isPointInsideSegment(const darray3E &p, const darray3E &V0, const darray3E &V1, double tol)
Definition: MimmoCGUtils.cpp:67
mimmo::skdTreeUtils::extractTarget
void extractTarget(bitpit::PatchSkdTree *target, const std::vector< const bitpit::SkdNode * > &leafSelection, std::vector< long > &extracted, double tol)
Definition: SkdTreeUtils.cpp:296
mimmo::skdTreeUtils::locatePointOnPatch
long locatePointOnPatch(const std::array< double, 3 > &point, const bitpit::PatchSkdTree *tree)
Definition: SkdTreeUtils.cpp:456
mimmo::skdTreeUtils::selectByPatch
std::vector< long > selectByPatch(bitpit::PatchSkdTree *selection, bitpit::PatchSkdTree *target, double tol)
Definition: SkdTreeUtils.cpp:251
darray2E
std::array< double, 2 > darray2E
Definition: mimmoTypeDef.hpp:49
mimmo::mimmoCGUtils::isPointInsidePolygon
bool isPointInsidePolygon(const darray3E &p, const dvecarr3E &vertCoords, double tol)
Definition: MimmoCGUtils.cpp:109
darray3E
std::array< double, 3 > darray3E
Definition: mimmoTypeDef.hpp:50
mimmo::skdTreeUtils::projectPoint
darray3E projectPoint(const std::array< double, 3 > *point, const bitpit::PatchSkdTree *skdtree, double r)
Definition: SkdTreeUtils.cpp:375
mimmo::mimmoCGUtils::isPointInsideTriangle
bool isPointInsideTriangle(const darray3E &p, const darray3E &V0, const darray3E &V1, const darray3E &V2, double tol)
Definition: MimmoCGUtils.cpp:85
mimmo::skdTreeUtils::distance
double distance(const std::array< double, 3 > *point, const bitpit::PatchSkdTree *tree, long &id, double r)
Definition: SkdTreeUtils.cpp:49
Generated on Tue Mar 16 2021 16:26:40 for mimmo by   doxygen 1.8.17