CGeometry.cpp

/*!
 * \file CGeometry.cpp
 * \brief Implementation of the base geometry class.
 * \author F. Palacios, T. Economon
 * \version 8.6.0 "Harrier"
 *
 * SU2 Project Website: https://su2code.github.io
 *
 * The SU2 Project is maintained by the SU2 Foundation
 * (http://su2foundation.org)
 *
 * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
 *
 * SU2 is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * SU2 is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
 */

#include <unordered_set>

#include "../../include/geometry/CGeometry.hpp"
#include "../../include/geometry/elements/CElement.hpp"
#include "../../include/parallelization/omp_structure.hpp"
#include "../../include/toolboxes/geometry_toolbox.hpp"
#include "../../include/toolboxes/ndflattener.hpp"

CGeometry::CGeometry() : size(SU2_MPI::GetSize()), rank(SU2_MPI::GetRank()) {}

CGeometry::~CGeometry() {
  unsigned long iElem, iElem_Bound, iVertex;
  unsigned short iMarker;

  if (elem != nullptr) {
    for (iElem = 0; iElem < nElem; iElem++) delete elem[iElem];
    delete[] elem;
  }

  if (bound != nullptr) {
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      for (iElem_Bound = 0; iElem_Bound < nElem_Bound[iMarker]; iElem_Bound++) {
        delete bound[iMarker][iElem_Bound];
      }
      delete[] bound[iMarker];
    }
    delete[] bound;
  }

  delete nodes;

  delete edges;

  if (vertex != nullptr) {
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
        delete vertex[iMarker][iVertex];
      }
      delete[] vertex[iMarker];
    }
    delete[] vertex;
  }

  delete[] nElem_Bound;
  delete[] nVertex;
  delete[] Marker_All_SendRecv;
  delete[] Tag_to_Marker;

  delete[] beg_node;
  delete[] end_node;
  delete[] nPointLinear;
  delete[] nPointCumulative;

  if (CustomBoundaryHeatFlux != nullptr) {
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      delete[] CustomBoundaryHeatFlux[iMarker];
    }
    delete[] CustomBoundaryHeatFlux;
  }

  if (CustomBoundaryTemperature != nullptr) {
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      delete[] CustomBoundaryTemperature[iMarker];
    }
    delete[] CustomBoundaryTemperature;
  }

  /*--- Delete structures for MPI point-to-point communication. ---*/

  delete[] bufD_P2PRecv;
  delete[] bufD_P2PSend;
#ifdef CODI_REVERSE_TYPE
  delete[] bufPD_P2PRecv;
  delete[] bufPD_P2PSend;
#endif
#ifdef USE_MIXED_PRECISION
  delete[] bufF_P2PRecv;
  delete[] bufF_P2PSend;
#endif
  delete[] bufS_P2PRecv;
  delete[] bufS_P2PSend;

  delete[] req_P2PSend;
  delete[] req_P2PRecv;

  delete[] reqP_P2PSend;
  delete[] reqP_P2PRecv;

  delete[] nPoint_P2PRecv;
  delete[] nPoint_P2PSend;

  delete[] Neighbors_P2PSend;
  delete[] Neighbors_P2PRecv;

  delete[] Local_Point_P2PSend;
  delete[] Local_Point_P2PRecv;

  /*--- Delete structures for MPI periodic communication. ---*/

  delete[] bufD_PeriodicRecv;
  delete[] bufD_PeriodicSend;

  delete[] bufS_PeriodicRecv;
  delete[] bufS_PeriodicSend;

  delete[] req_PeriodicSend;
  delete[] req_PeriodicRecv;

  delete[] nPoint_PeriodicRecv;
  delete[] nPoint_PeriodicSend;

  delete[] Neighbors_PeriodicSend;
  delete[] Neighbors_PeriodicRecv;

  delete[] Local_Point_PeriodicSend;
  delete[] Local_Point_PeriodicRecv;

  delete[] Local_Marker_PeriodicSend;
  delete[] Local_Marker_PeriodicRecv;
}

void CGeometry::PreprocessP2PComms(CGeometry* geometry, CConfig* config) {
  /*--- We start with the send and receive lists already available in
   the form of SEND_RECEIVE boundary markers. We will loop through
   these markers and establish the neighboring ranks and number of
   send/recv points per pair. We will store this information and set
   up persistent data structures so that we can reuse them throughout
   the calculation for any point-to-point communications. The goal
   is to break the non-blocking comms into InitiateComms() and
   CompleteComms() in separate routines so that we can overlap the
   communication and computation to hide the communication latency. ---*/

  /*--- Local variables. ---*/

  unsigned short iMarker;
  unsigned long nVertexS, nVertexR, iVertex, MarkerS, MarkerR;

  int iRank, iSend, iRecv, count;

  /*--- Create some temporary structures for tracking sends/recvs. ---*/

  int* nPoint_Send_All = new int[size + 1];
  nPoint_Send_All[0] = 0;
  int* nPoint_Recv_All = new int[size + 1];
  nPoint_Recv_All[0] = 0;
  int* nPoint_Flag = new int[size];

  for (iRank = 0; iRank < size; iRank++) {
    nPoint_Send_All[iRank] = 0;
    nPoint_Recv_All[iRank] = 0;
    nPoint_Flag[iRank] = -1;
  }
  nPoint_Send_All[size] = 0;
  nPoint_Recv_All[size] = 0;

  /*--- Loop through all of our SEND_RECEIVE markers and track
   our sends with each rank. ---*/

  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
    if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0)) {
      /*--- Get the destination rank and number of points to send. ---*/

      iRank = config->GetMarker_All_SendRecv(iMarker) - 1;
      nVertexS = geometry->nVertex[iMarker];

      /*--- If we have not visited this element yet, increment our
       number of elements that must be sent to a particular proc. ---*/

      if ((nPoint_Flag[iRank] != static_cast<int>(iMarker))) {
        nPoint_Flag[iRank] = static_cast<int>(iMarker);
        nPoint_Send_All[iRank + 1] += nVertexS;
      }
    }
  }

  delete[] nPoint_Flag;

  /*--- Communicate the number of points to be sent/recv'd amongst
   all processors. After this communication, each proc knows how
   many cells it will receive from each other processor. ---*/

  SU2_MPI::Alltoall(&(nPoint_Send_All[1]), 1, MPI_INT, &(nPoint_Recv_All[1]), 1, MPI_INT, SU2_MPI::GetComm());

  /*--- Prepare to send connectivities. First check how many
   messages we will be sending and receiving. Here we also put
   the counters into cumulative storage format to make the
   communications simpler. ---*/

  nP2PSend = 0;
  nP2PRecv = 0;

  for (iRank = 0; iRank < size; iRank++) {
    if ((iRank != rank) && (nPoint_Send_All[iRank + 1] > 0)) nP2PSend++;
    if ((iRank != rank) && (nPoint_Recv_All[iRank + 1] > 0)) nP2PRecv++;

    nPoint_Send_All[iRank + 1] += nPoint_Send_All[iRank];
    nPoint_Recv_All[iRank + 1] += nPoint_Recv_All[iRank];
  }

  /*--- Allocate only as much memory as we need for the P2P neighbors. ---*/

  nPoint_P2PSend = new int[nP2PSend + 1];
  nPoint_P2PSend[0] = 0;
  nPoint_P2PRecv = new int[nP2PRecv + 1];
  nPoint_P2PRecv[0] = 0;

  Neighbors_P2PSend = new int[nP2PSend];
  Neighbors_P2PRecv = new int[nP2PRecv];

  iSend = 0;
  iRecv = 0;
  for (iRank = 0; iRank < size; iRank++) {
    if ((nPoint_Send_All[iRank + 1] > nPoint_Send_All[iRank]) && (iRank != rank)) {
      Neighbors_P2PSend[iSend] = iRank;
      nPoint_P2PSend[iSend + 1] = nPoint_Send_All[iRank + 1];
      iSend++;
    }

    if ((nPoint_Recv_All[iRank + 1] > nPoint_Recv_All[iRank]) && (iRank != rank)) {
      Neighbors_P2PRecv[iRecv] = iRank;
      nPoint_P2PRecv[iRecv + 1] = nPoint_Recv_All[iRank + 1];
      iRecv++;
    }
  }

  /*--- Create a reverse mapping of the message to the rank so that we
   can quickly access the correct data in the buffers when receiving
   messages dynamically. ---*/

  P2PSend2Neighbor.clear();
  for (iSend = 0; iSend < nP2PSend; iSend++) P2PSend2Neighbor[Neighbors_P2PSend[iSend]] = iSend;

  P2PRecv2Neighbor.clear();
  for (iRecv = 0; iRecv < nP2PRecv; iRecv++) P2PRecv2Neighbor[Neighbors_P2PRecv[iRecv]] = iRecv;

  delete[] nPoint_Send_All;
  delete[] nPoint_Recv_All;

  /*--- Allocate the memory that we need for receiving the conn
   values and then cue up the non-blocking receives. Note that
   we do not include our own rank in the communications. We will
   directly copy our own data later. ---*/

  Local_Point_P2PSend = nullptr;
  Local_Point_P2PSend = new unsigned long[nPoint_P2PSend[nP2PSend]];
  for (iSend = 0; iSend < nPoint_P2PSend[nP2PSend]; iSend++) Local_Point_P2PSend[iSend] = 0;

  Local_Point_P2PRecv = nullptr;
  Local_Point_P2PRecv = new unsigned long[nPoint_P2PRecv[nP2PRecv]];
  for (iRecv = 0; iRecv < nPoint_P2PRecv[nP2PRecv]; iRecv++) Local_Point_P2PRecv[iRecv] = 0;

  /*--- We allocate the memory for communicating values in a later step
   once we know the maximum packet size that we need to communicate. This
   memory is deallocated and reallocated automatically in the case that
   the previously allocated memory is not sufficient. ---*/

  bufD_P2PSend = nullptr;
  bufD_P2PRecv = nullptr;

  bufS_P2PSend = nullptr;
  bufS_P2PRecv = nullptr;

  /*--- Allocate memory for the MPI requests if we need to communicate. ---*/

  if (nP2PSend > 0) req_P2PSend = new SU2_MPI::Request[nP2PSend];
  if (nP2PRecv > 0) req_P2PRecv = new SU2_MPI::Request[nP2PRecv];

  if (nP2PSend > 0) reqP_P2PSend = new PassiveRequest[nP2PSend];
  if (nP2PRecv > 0) reqP_P2PRecv = new PassiveRequest[nP2PRecv];

  /*--- Build lists of local index values for send. ---*/

  count = 0;
  for (iSend = 0; iSend < nP2PSend; iSend++) {
    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
      if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0)) {
        MarkerS = iMarker;
        nVertexS = geometry->nVertex[MarkerS];
        iRank = config->GetMarker_All_SendRecv(MarkerS) - 1;

        if (iRank == Neighbors_P2PSend[iSend]) {
          for (iVertex = 0; iVertex < nVertexS; iVertex++) {
            Local_Point_P2PSend[count] = geometry->vertex[MarkerS][iVertex]->GetNode();
            count++;
          }
        }
      }
    }
  }

  /*--- Build lists of local index values for receive. ---*/

  count = 0;
  for (iRecv = 0; iRecv < nP2PRecv; iRecv++) {
    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
      if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0)) {
        MarkerR = iMarker + 1;
        nVertexR = geometry->nVertex[MarkerR];
        iRank = abs(config->GetMarker_All_SendRecv(MarkerR)) - 1;

        if (iRank == Neighbors_P2PRecv[iRecv]) {
          for (iVertex = 0; iVertex < nVertexR; iVertex++) {
            Local_Point_P2PRecv[count] = geometry->vertex[MarkerR][iVertex]->GetNode();
            count++;
          }
        }
      }
    }
  }

  /*--- In the future, some additional data structures could be created
   here to separate the interior and boundary nodes in order to help
   further overlap computation and communication. ---*/
}

void CGeometry::AllocateP2PComms(unsigned short countPerPoint) {
  /*--- This routine is activated whenever we attempt to perform
   a point-to-point MPI communication with our neighbors but the
   memory buffer allocated is not large enough for the packet size.
   Therefore, we deallocate the previously allocated space and
   reallocate a large enough array. Note that after the first set
   communications, this routine will not need to be called again. ---*/

  if (countPerPoint <= maxCountPerPoint) return;

  BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
    /*--- Store the larger packet size to the class data. ---*/

    maxCountPerPoint = countPerPoint;
    const auto send_size = maxCountPerPoint * nPoint_P2PSend[nP2PSend];
    const auto recv_size = maxCountPerPoint * nPoint_P2PRecv[nP2PSend];

    /*-- Deallocate and reallocate our su2double cummunication memory. ---*/
#define SU2_ALLOC_COMM_BUFFERS(SEND, RECV, TYPE) \
  delete[] SEND;                                 \
  SEND = new TYPE[send_size]();                  \
  delete[] RECV;                                 \
  RECV = new TYPE[recv_size]();

    SU2_ALLOC_COMM_BUFFERS(bufD_P2PSend, bufD_P2PRecv, su2double)
    SU2_ALLOC_COMM_BUFFERS(bufS_P2PSend, bufS_P2PRecv, unsigned short)
#ifdef CODI_REVERSE_TYPE
    SU2_ALLOC_COMM_BUFFERS(bufPD_P2PSend, bufPD_P2PRecv, passivedouble)
#endif
#ifdef USE_MIXED_PRECISION
    SU2_ALLOC_COMM_BUFFERS(bufF_P2PSend, bufF_P2PRecv, su2mixedfloat)
#endif
#undef SU2_ALLOC_COMM_BUFFERS
  }
  END_SU2_OMP_SAFE_GLOBAL_ACCESS
}

void CGeometry::PostP2PRecvs(CGeometry* geometry, const CConfig* config, COMM_TYPE commType,
                             unsigned short countPerPoint, bool val_reverse) const {
  /*--- Launch the non-blocking recv's first. Note that we have stored
   the counts and sources, so we can launch these before we even load
   the data and send from the neighbor ranks. ---*/

  SU2_OMP_MASTER
  for (int iRecv = 0; iRecv < nP2PRecv; iRecv++) {
    const auto iMessage = iRecv;

    /*--- In some instances related to the adjoint solver, we need
     to reverse the direction of communications such that the normal
     send nodes become the recv nodes and vice-versa. ---*/

    if (val_reverse) {
      /*--- Compute our location in the buffer using the send data
       structure since we are reversing the comms. ---*/

      auto offset = countPerPoint * nPoint_P2PSend[iRecv];

      /*--- Take advantage of cumulative storage format to get the number
       of elems that we need to recv. Note again that we select the send
       points here as the recv points. ---*/

      auto nPointP2P = nPoint_P2PSend[iRecv + 1] - nPoint_P2PSend[iRecv];

      /*--- Total count can include multiple pieces of data per element. ---*/

      auto count = countPerPoint * nPointP2P;

      /*--- Get the rank from which we receive the message. Note again
       that we use the send rank as the source instead of the recv rank. ---*/

      auto source = Neighbors_P2PSend[iRecv];
      auto tag = source + 1;

      /*--- Post non-blocking recv for this proc. Note that we use the
       send buffer here too. This is important to make sure the arrays
       are the correct size. ---*/

      switch (commType) {
        case COMM_TYPE::DOUBLE:
          SU2_MPI::Irecv(&(bufD_P2PSend[offset]), count, MPI_DOUBLE, source, tag, SU2_MPI::GetComm(),
                         &(req_P2PRecv[iRecv]));
          break;
#ifdef CODI_REVERSE_TYPE
        case COMM_TYPE::PASSIVE_DOUBLE:
          SelectMPIWrapper<passivedouble>::W::Irecv(&(bufPD_P2PSend[offset]), count, MPI_DOUBLE, source, tag,
                                                    SU2_MPI::GetComm(), &GetP2PRecvReq<passivedouble>()[iRecv]);
          break;
#endif
#ifdef USE_MIXED_PRECISION
        case COMM_TYPE::FLOAT:
          SelectMPIWrapper<su2mixedfloat>::W::Irecv(&(bufF_P2PSend[offset]), count, MPI_FLOAT, source, tag,
                                                    SU2_MPI::GetComm(), &GetP2PRecvReq<su2mixedfloat>()[iRecv]);
          break;
#endif
        case COMM_TYPE::UNSIGNED_SHORT:
          SU2_MPI::Irecv(&(bufS_P2PSend[offset]), count, MPI_UNSIGNED_SHORT, source, tag, SU2_MPI::GetComm(),
                         &(req_P2PRecv[iRecv]));
          break;
        default:
          SU2_MPI::Error("Unrecognized data type for point-to-point MPI comms.", CURRENT_FUNCTION);
          break;
      }

    } else {
      /*--- Compute our location in the recv buffer. ---*/

      auto offset = countPerPoint * nPoint_P2PRecv[iRecv];

      /*--- Take advantage of cumulative storage format to get the number
       of elems that we need to recv. ---*/

      auto nPointP2P = nPoint_P2PRecv[iRecv + 1] - nPoint_P2PRecv[iRecv];

      /*--- Total count can include multiple pieces of data per element. ---*/

      auto count = countPerPoint * nPointP2P;

      /*--- Get the rank from which we receive the message. ---*/

      auto source = Neighbors_P2PRecv[iRecv];
      auto tag = source + 1;

      /*--- Post non-blocking recv for this proc. ---*/

      switch (commType) {
        case COMM_TYPE::DOUBLE:
          SU2_MPI::Irecv(&(bufD_P2PRecv[offset]), count, MPI_DOUBLE, source, tag, SU2_MPI::GetComm(),
                         &(req_P2PRecv[iMessage]));
          break;
#ifdef CODI_REVERSE_TYPE
        case COMM_TYPE::PASSIVE_DOUBLE:
          SelectMPIWrapper<passivedouble>::W::Irecv(&(bufPD_P2PRecv[offset]), count, MPI_DOUBLE, source, tag,
                                                    SU2_MPI::GetComm(), &GetP2PRecvReq<passivedouble>()[iMessage]);
          break;
#endif
#ifdef USE_MIXED_PRECISION
        case COMM_TYPE::FLOAT:
          SelectMPIWrapper<su2mixedfloat>::W::Irecv(&(bufF_P2PRecv[offset]), count, MPI_FLOAT, source, tag,
                                                    SU2_MPI::GetComm(), &GetP2PRecvReq<su2mixedfloat>()[iMessage]);
          break;
#endif
        case COMM_TYPE::UNSIGNED_SHORT:
          SU2_MPI::Irecv(&(bufS_P2PRecv[offset]), count, MPI_UNSIGNED_SHORT, source, tag, SU2_MPI::GetComm(),
                         &(req_P2PRecv[iMessage]));
          break;
        default:
          SU2_MPI::Error("Unrecognized data type for point-to-point MPI comms.", CURRENT_FUNCTION);
          break;
      }
    }
  }
  END_SU2_OMP_MASTER
}

void CGeometry::PostP2PSends(CGeometry* geometry, const CConfig* config, COMM_TYPE commType,
                             unsigned short countPerPoint, int val_iSend, bool val_reverse) const {
  /*--- Post the non-blocking send as soon as the buffer is loaded. ---*/

  /*--- In some instances related to the adjoint solver, we need
   to reverse the direction of communications such that the normal
   send nodes become the recv nodes and vice-versa. ---*/

  SU2_OMP_MASTER
  if (val_reverse) {
    /*--- Compute our location in the buffer using the recv data
     structure since we are reversing the comms. ---*/

    auto offset = countPerPoint * nPoint_P2PRecv[val_iSend];

    /*--- Take advantage of cumulative storage format to get the number
     of points that we need to send. Note again that we select the recv
     points here as the send points. ---*/

    auto nPointP2P = nPoint_P2PRecv[val_iSend + 1] - nPoint_P2PRecv[val_iSend];

    /*--- Total count can include multiple pieces of data per element. ---*/

    auto count = countPerPoint * nPointP2P;

    /*--- Get the rank to which we send the message. Note again
     that we use the recv rank as the dest instead of the send rank. ---*/

    auto dest = Neighbors_P2PRecv[val_iSend];
    auto tag = rank + 1;

    /*--- Post non-blocking send for this proc. Note that we use the
     send buffer here too. This is important to make sure the arrays
     are the correct size. ---*/

    switch (commType) {
      case COMM_TYPE::DOUBLE:
        SU2_MPI::Isend(&(bufD_P2PRecv[offset]), count, MPI_DOUBLE, dest, tag, SU2_MPI::GetComm(),
                       &(req_P2PSend[val_iSend]));
        break;
#ifdef CODI_REVERSE_TYPE
      case COMM_TYPE::PASSIVE_DOUBLE:
        SelectMPIWrapper<passivedouble>::W::Isend(&(bufPD_P2PRecv[offset]), count, MPI_DOUBLE, dest, tag,
                                                  SU2_MPI::GetComm(), &GetP2PSendReq<passivedouble>()[val_iSend]);
        break;
#endif
#ifdef USE_MIXED_PRECISION
      case COMM_TYPE::FLOAT:
        SelectMPIWrapper<su2mixedfloat>::W::Isend(&(bufF_P2PRecv[offset]), count, MPI_FLOAT, dest, tag,
                                                  SU2_MPI::GetComm(), &GetP2PSendReq<su2mixedfloat>()[val_iSend]);
        break;
#endif
      case COMM_TYPE::UNSIGNED_SHORT:
        SU2_MPI::Isend(&(bufS_P2PRecv[offset]), count, MPI_UNSIGNED_SHORT, dest, tag, SU2_MPI::GetComm(),
                       &(req_P2PSend[val_iSend]));
        break;
      default:
        SU2_MPI::Error("Unrecognized data type for point-to-point MPI comms.", CURRENT_FUNCTION);
        break;
    }

  } else {
    /*--- Compute our location in the send buffer. ---*/

    auto offset = countPerPoint * nPoint_P2PSend[val_iSend];

    /*--- Take advantage of cumulative storage format to get the number
     of points that we need to send. ---*/

    auto nPointP2P = nPoint_P2PSend[val_iSend + 1] - nPoint_P2PSend[val_iSend];

    /*--- Total count can include multiple pieces of data per element. ---*/

    auto count = countPerPoint * nPointP2P;

    /*--- Get the rank to which we send the message. ---*/

    auto dest = Neighbors_P2PSend[val_iSend];
    auto tag = rank + 1;

    /*--- Post non-blocking send for this proc. ---*/

    switch (commType) {
      case COMM_TYPE::DOUBLE:
        SU2_MPI::Isend(&(bufD_P2PSend[offset]), count, MPI_DOUBLE, dest, tag, SU2_MPI::GetComm(),
                       &(req_P2PSend[val_iSend]));
        break;
#ifdef CODI_REVERSE_TYPE
      case COMM_TYPE::PASSIVE_DOUBLE:
        SelectMPIWrapper<passivedouble>::W::Isend(&(bufPD_P2PSend[offset]), count, MPI_DOUBLE, dest, tag,
                                                  SU2_MPI::GetComm(), &GetP2PSendReq<passivedouble>()[val_iSend]);
        break;
#endif
#ifdef USE_MIXED_PRECISION
      case COMM_TYPE::FLOAT:
        SelectMPIWrapper<su2mixedfloat>::W::Isend(&(bufF_P2PSend[offset]), count, MPI_FLOAT, dest, tag,
                                                  SU2_MPI::GetComm(), &GetP2PSendReq<su2mixedfloat>()[val_iSend]);
        break;
#endif
      case COMM_TYPE::UNSIGNED_SHORT:
        SU2_MPI::Isend(&(bufS_P2PSend[offset]), count, MPI_UNSIGNED_SHORT, dest, tag, SU2_MPI::GetComm(),
                       &(req_P2PSend[val_iSend]));
        break;
      default:
        SU2_MPI::Error("Unrecognized data type for point-to-point MPI comms.", CURRENT_FUNCTION);
        break;
    }
  }
  END_SU2_OMP_MASTER
}

void CGeometry::GetCommCountAndType(const CConfig* config, MPI_QUANTITIES commType, unsigned short& COUNT_PER_POINT,
                                    COMM_TYPE& MPI_TYPE) const {
  switch (commType) {
    case MPI_QUANTITIES::COORDINATES:
      COUNT_PER_POINT = nDim;
      MPI_TYPE = COMM_TYPE::DOUBLE;
      break;
    case MPI_QUANTITIES::GRID_VELOCITY:
      COUNT_PER_POINT = nDim;
      MPI_TYPE = COMM_TYPE::DOUBLE;
      break;
    case MPI_QUANTITIES::COORDINATES_OLD:
      if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)
        COUNT_PER_POINT = nDim * 2;
      else
        COUNT_PER_POINT = nDim;
      MPI_TYPE = COMM_TYPE::DOUBLE;
      break;
    case MPI_QUANTITIES::MAX_LENGTH:
      COUNT_PER_POINT = 1;
      MPI_TYPE = COMM_TYPE::DOUBLE;
      break;
    case MPI_QUANTITIES::NEIGHBORS:
      COUNT_PER_POINT = 1;
      MPI_TYPE = COMM_TYPE::UNSIGNED_SHORT;
      break;
    default:
      SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION);
      break;
  }
}

void CGeometry::InitiateComms(CGeometry* geometry, const CConfig* config, MPI_QUANTITIES commType) const {
  if (nP2PSend == 0) return;

  /*--- Local variables ---*/

  unsigned short iDim;
  unsigned short COUNT_PER_POINT = 0;
  COMM_TYPE MPI_TYPE{};

  unsigned long iPoint, msg_offset, buf_offset;

  int iMessage, iSend, nSend;

  /*--- Set the size of the data packet and type depending on quantity. ---*/

  GetCommCountAndType(config, commType, COUNT_PER_POINT, MPI_TYPE);

  /*--- Check to make sure we have created a large enough buffer
   for these comms during preprocessing. This is only for the su2double
   buffer. It will be reallocated whenever we find a larger count
   per point. After the first cycle of comms, this should be inactive. ---*/

  geometry->AllocateP2PComms(COUNT_PER_POINT);

  /*--- Set some local pointers to make access simpler. ---*/

  su2double* bufDSend = geometry->bufD_P2PSend;
  unsigned short* bufSSend = geometry->bufS_P2PSend;

  su2double* vector = nullptr;

  /*--- Load the specified quantity from the solver into the generic
   communication buffer in the geometry class. ---*/

  /*--- Post all non-blocking recvs first before sends. ---*/

  geometry->PostP2PRecvs(geometry, config, MPI_TYPE, COUNT_PER_POINT, false);

  for (iMessage = 0; iMessage < nP2PSend; iMessage++) {
    /*--- Get the offset in the buffer for the start of this message. ---*/

    msg_offset = nPoint_P2PSend[iMessage];

    /*--- Total count can include multiple pieces of data per element. ---*/

    nSend = (nPoint_P2PSend[iMessage + 1] - nPoint_P2PSend[iMessage]);

    SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
    for (iSend = 0; iSend < nSend; iSend++) {
      /*--- Get the local index for this communicated data. ---*/

      iPoint = geometry->Local_Point_P2PSend[msg_offset + iSend];

      /*--- Compute the offset in the recv buffer for this point. ---*/

      buf_offset = (msg_offset + iSend) * COUNT_PER_POINT;

      switch (commType) {
        case MPI_QUANTITIES::COORDINATES:
          vector = nodes->GetCoord(iPoint);
          for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset + iDim] = vector[iDim];
          break;
        case MPI_QUANTITIES::GRID_VELOCITY:
          vector = nodes->GetGridVel(iPoint);
          for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset + iDim] = vector[iDim];
          break;
        case MPI_QUANTITIES::COORDINATES_OLD:
          vector = nodes->GetCoord_n(iPoint);
          for (iDim = 0; iDim < nDim; iDim++) {
            bufDSend[buf_offset + iDim] = vector[iDim];
          }
          if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
            vector = nodes->GetCoord_n1(iPoint);
            for (iDim = 0; iDim < nDim; iDim++) {
              bufDSend[buf_offset + nDim + iDim] = vector[iDim];
            }
          }
          break;
        case MPI_QUANTITIES::MAX_LENGTH:
          bufDSend[buf_offset] = nodes->GetMaxLength(iPoint);
          break;
        case MPI_QUANTITIES::NEIGHBORS:
          bufSSend[buf_offset] = geometry->nodes->GetnNeighbor(iPoint);
          break;
        default:
          SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION);
          break;
      }
    }
    END_SU2_OMP_FOR

    /*--- Launch the point-to-point MPI send for this message. ---*/

    geometry->PostP2PSends(geometry, config, MPI_TYPE, COUNT_PER_POINT, iMessage, false);
  }
}

void CGeometry::CompleteComms(CGeometry* geometry, const CConfig* config, MPI_QUANTITIES commType) {
  if (nP2PRecv == 0) return;

  /*--- Local variables ---*/

  unsigned short iDim, COUNT_PER_POINT = 0;
  COMM_TYPE MPI_TYPE{};
  unsigned long iPoint, iRecv, nRecv, msg_offset, buf_offset;

  int ind, source, iMessage, jRecv;

  /*--- Global status so all threads can see the result of Waitany. ---*/
  static SU2_MPI::Status status;

  /*--- Set the size of the data packet and type depending on quantity. ---*/

  GetCommCountAndType(config, commType, COUNT_PER_POINT, MPI_TYPE);

  /*--- Set some local pointers to make access simpler. ---*/

  const su2double* bufDRecv = geometry->bufD_P2PRecv;
  const unsigned short* bufSRecv = geometry->bufS_P2PRecv;

  /*--- Store the data that was communicated into the appropriate
   location within the local class data structures. Note that we
   recv and store the data in any order to take advantage of the
   non-blocking comms. ---*/

  for (iMessage = 0; iMessage < nP2PRecv; iMessage++) {
    /*--- For efficiency, recv the messages dynamically based on
     the order they arrive. ---*/

    SU2_OMP_SAFE_GLOBAL_ACCESS(SU2_MPI::Waitany(nP2PRecv, req_P2PRecv, &ind, &status);)

    /*--- Once we have recv'd a message, get the source rank. ---*/

    source = status.MPI_SOURCE;

    /*--- We know the offsets based on the source rank. ---*/

    jRecv = P2PRecv2Neighbor[source];

    /*--- Get the offset in the buffer for the start of this message. ---*/

    msg_offset = nPoint_P2PRecv[jRecv];

    /*--- Get the number of packets to be received in this message. ---*/

    nRecv = nPoint_P2PRecv[jRecv + 1] - nPoint_P2PRecv[jRecv];

    SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
    for (iRecv = 0; iRecv < nRecv; iRecv++) {
      /*--- Get the local index for this communicated data. ---*/

      iPoint = geometry->Local_Point_P2PRecv[msg_offset + iRecv];

      /*--- Compute the total offset in the recv buffer for this point. ---*/

      buf_offset = (msg_offset + iRecv) * COUNT_PER_POINT;

      /*--- Store the data correctly depending on the quantity. ---*/

      switch (commType) {
        case MPI_QUANTITIES::COORDINATES:
          for (iDim = 0; iDim < nDim; iDim++) nodes->SetCoord(iPoint, iDim, bufDRecv[buf_offset + iDim]);
          break;
        case MPI_QUANTITIES::GRID_VELOCITY:
          for (iDim = 0; iDim < nDim; iDim++) nodes->SetGridVel(iPoint, iDim, bufDRecv[buf_offset + iDim]);
          break;
        case MPI_QUANTITIES::COORDINATES_OLD:
          nodes->SetCoord_n(iPoint, &bufDRecv[buf_offset]);
          if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)
            nodes->SetCoord_n1(iPoint, &bufDRecv[buf_offset + nDim]);
          break;
        case MPI_QUANTITIES::MAX_LENGTH:
          nodes->SetMaxLength(iPoint, bufDRecv[buf_offset]);
          break;
        case MPI_QUANTITIES::NEIGHBORS:
          nodes->SetnNeighbor(iPoint, bufSRecv[buf_offset]);
          break;
        default:
          SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION);
          break;
      }
    }
    END_SU2_OMP_FOR
  }

  /*--- Verify that all non-blocking point-to-point sends have finished.
   Note that this should be satisfied, as we have received all of the
   data in the loop above at this point. ---*/

#ifdef HAVE_MPI
  SU2_OMP_SAFE_GLOBAL_ACCESS(SU2_MPI::Waitall(nP2PSend, req_P2PSend, MPI_STATUS_IGNORE);)
#endif
}

void CGeometry::MatchPeriodic(const CConfig* config, unsigned short val_periodic) {
  unsigned short iMarker, iDim, jMarker, pMarker = 0;
  unsigned short iPeriodic, nPeriodic;

  unsigned long iVertex, iPoint, iPointGlobal, index;
  unsigned long jVertex, jVertex_, jPoint, jPointGlobal;
  unsigned long pVertex = 0, pPoint = 0, pPointGlobal = 0;
  unsigned long nLocalVertex_Periodic = 0, MaxLocalVertex_Periodic = 0;
  unsigned long nPointMatch = 0;

  int iProcessor, pProcessor = 0, nProcessor = size;

  bool isBadMatch = false;

  string Marker_Tag;

  su2double Coord_j[3], dist, mindist, maxdist_local, maxdist_global;
  const su2double *center, *angles, *trans;
  su2double translation[3] = {0.0, 0.0, 0.0}, dx, dy, dz;
  su2double rotMatrix[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
  su2double rotCoord[3] = {0.0, 0.0, 0.0};

  bool pointOnAxis = false;
  bool chkSamePoint = false;
  su2double distToAxis = 0.0;

  /*--- Tolerance for distance-based match to report warning. ---*/

  su2double epsilon = 1e-6;

  /*--- Evaluate the number of periodic boundary conditions ---*/

  nPeriodic = config->GetnMarker_Periodic();

  /*--- Send an initial message to the console. ---*/

  if (rank == MASTER_NODE) {
    cout << "Matching the periodic boundary points for marker pair ";
    cout << val_periodic << "." << endl;
  }

  /*--- Compute the total number of vertices that sit on a periodic
   boundary on our local rank. We only include our "owned" nodes. ---*/

  nLocalVertex_Periodic = 0;
  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
    if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
      iPeriodic = config->GetMarker_All_PerBound(iMarker);
      if ((iPeriodic == val_periodic) || (iPeriodic == val_periodic + nPeriodic / 2)) {
        for (iVertex = 0; iVertex < GetnVertex(iMarker); iVertex++) {
          iPoint = vertex[iMarker][iVertex]->GetNode();
          if (nodes->GetDomain(iPoint)) nLocalVertex_Periodic++;
        }
      }
    }
  }

  /*--- Communicate our local periodic point count globally
   and receive the counts of periodic points from all other ranks.---*/

  auto* Buffer_Send_nVertex = new unsigned long[1];
  auto* Buffer_Recv_nVertex = new unsigned long[nProcessor];

  Buffer_Send_nVertex[0] = nLocalVertex_Periodic;

  /*--- Copy our own count in serial or use collective comms with MPI. ---*/

  SU2_MPI::Allreduce(&nLocalVertex_Periodic, &MaxLocalVertex_Periodic, 1, MPI_UNSIGNED_LONG, MPI_MAX,
                     SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_nVertex, 1, MPI_UNSIGNED_LONG, Buffer_Recv_nVertex, 1, MPI_UNSIGNED_LONG,
                     SU2_MPI::GetComm());

  /*--- Prepare buffers to send the information for each
   periodic point to all ranks so that we can match pairs. ---*/

  auto* Buffer_Send_Coord = new su2double[MaxLocalVertex_Periodic * nDim];
  auto* Buffer_Send_Point = new unsigned long[MaxLocalVertex_Periodic];
  auto* Buffer_Send_GlobalIndex = new unsigned long[MaxLocalVertex_Periodic];
  auto* Buffer_Send_Vertex = new unsigned long[MaxLocalVertex_Periodic];
  auto* Buffer_Send_Marker = new unsigned long[MaxLocalVertex_Periodic];

  auto* Buffer_Recv_Coord = new su2double[nProcessor * MaxLocalVertex_Periodic * nDim];
  auto* Buffer_Recv_Point = new unsigned long[nProcessor * MaxLocalVertex_Periodic];
  auto* Buffer_Recv_GlobalIndex = new unsigned long[nProcessor * MaxLocalVertex_Periodic];
  auto* Buffer_Recv_Vertex = new unsigned long[nProcessor * MaxLocalVertex_Periodic];
  auto* Buffer_Recv_Marker = new unsigned long[nProcessor * MaxLocalVertex_Periodic];

  unsigned long nBuffer_Coord = MaxLocalVertex_Periodic * nDim;
  unsigned long nBuffer_Point = MaxLocalVertex_Periodic;
  unsigned long nBuffer_GlobalIndex = MaxLocalVertex_Periodic;
  unsigned long nBuffer_Vertex = MaxLocalVertex_Periodic;
  unsigned long nBuffer_Marker = MaxLocalVertex_Periodic;

  for (iVertex = 0; iVertex < MaxLocalVertex_Periodic; iVertex++) {
    Buffer_Send_Point[iVertex] = 0;
    Buffer_Send_GlobalIndex[iVertex] = 0;
    Buffer_Send_Vertex[iVertex] = 0;
    Buffer_Send_Marker[iVertex] = 0;
    for (iDim = 0; iDim < nDim; iDim++) Buffer_Send_Coord[iVertex * nDim + iDim] = 0.0;
  }

  /*--- Store the local index, global index, local boundary index,
   marker index, and point coordinates in the buffers for sending.
   Note again that this is only for the current pair of periodic
   markers and for only the "owned" points on each rank. ---*/

  nLocalVertex_Periodic = 0;
  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
    if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
      iPeriodic = config->GetMarker_All_PerBound(iMarker);
      if ((iPeriodic == val_periodic) || (iPeriodic == val_periodic + nPeriodic / 2)) {
        for (iVertex = 0; iVertex < GetnVertex(iMarker); iVertex++) {
          iPoint = vertex[iMarker][iVertex]->GetNode();
          iPointGlobal = nodes->GetGlobalIndex(iPoint);
          if (nodes->GetDomain(iPoint)) {
            Buffer_Send_Point[nLocalVertex_Periodic] = iPoint;
            Buffer_Send_GlobalIndex[nLocalVertex_Periodic] = iPointGlobal;
            Buffer_Send_Vertex[nLocalVertex_Periodic] = iVertex;
            Buffer_Send_Marker[nLocalVertex_Periodic] = iMarker;
            for (iDim = 0; iDim < nDim; iDim++)
              Buffer_Send_Coord[nLocalVertex_Periodic * nDim + iDim] = nodes->GetCoord(iPoint, iDim);
            nLocalVertex_Periodic++;
          }
        }
      }
    }
  }

  /*--- Gather the data for all points on each rank with MPI. ---*/

  SU2_MPI::Allgather(Buffer_Send_Coord, nBuffer_Coord, MPI_DOUBLE, Buffer_Recv_Coord, nBuffer_Coord, MPI_DOUBLE,
                     SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_Point, nBuffer_Point, MPI_UNSIGNED_LONG, Buffer_Recv_Point, nBuffer_Point,
                     MPI_UNSIGNED_LONG, SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_GlobalIndex, nBuffer_GlobalIndex, MPI_UNSIGNED_LONG, Buffer_Recv_GlobalIndex,
                     nBuffer_GlobalIndex, MPI_UNSIGNED_LONG, SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_Vertex, nBuffer_Vertex, MPI_UNSIGNED_LONG, Buffer_Recv_Vertex, nBuffer_Vertex,
                     MPI_UNSIGNED_LONG, SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_Marker, nBuffer_Marker, MPI_UNSIGNED_LONG, Buffer_Recv_Marker, nBuffer_Marker,
                     MPI_UNSIGNED_LONG, SU2_MPI::GetComm());

  /*--- Now that all ranks have the data for all periodic points for
   this pair of periodic markers, we match the individual points
   based on the translation / rotation specified for the marker pair. ---*/

  maxdist_local = 0.0;
  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
    if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
      iPeriodic = config->GetMarker_All_PerBound(iMarker);
      if ((iPeriodic == val_periodic) || (iPeriodic == val_periodic + nPeriodic / 2)) {
        /*--- Retrieve the supplied periodic information. ---*/

        Marker_Tag = config->GetMarker_All_TagBound(iMarker);
        center = config->GetPeriodicRotCenter(Marker_Tag);
        angles = config->GetPeriodicRotAngles(Marker_Tag);
        trans = config->GetPeriodicTranslation(Marker_Tag);

        /*--- Store (center+trans) as it is constant and will be added. ---*/

        translation[0] = center[0] + trans[0];
        translation[1] = center[1] + trans[1];
        translation[2] = center[2] + trans[2];

        /*--- Store angles separately for clarity. Compute sines/cosines. ---*/

        const su2double Theta = angles[0];
        const su2double Phi = angles[1];
        const su2double Psi = angles[2];
        const su2double cosTheta = cos(Theta), sinTheta = sin(Theta);
        const su2double cosPhi = cos(Phi), sinPhi = sin(Phi);
        const su2double cosPsi = cos(Psi), sinPsi = sin(Psi);

        /*--- Compute the rotation matrix. Note that the implicit
         ordering is rotation about the x-axis, y-axis, then z-axis. ---*/

        rotMatrix[0][0] = cosPhi * cosPsi;
        rotMatrix[1][0] = cosPhi * sinPsi;
        rotMatrix[2][0] = -sinPhi;

        rotMatrix[0][1] = sinTheta * sinPhi * cosPsi - cosTheta * sinPsi;
        rotMatrix[1][1] = sinTheta * sinPhi * sinPsi + cosTheta * cosPsi;
        rotMatrix[2][1] = sinTheta * cosPhi;

        rotMatrix[0][2] = cosTheta * sinPhi * cosPsi + sinTheta * sinPsi;
        rotMatrix[1][2] = cosTheta * sinPhi * sinPsi - sinTheta * cosPsi;
        rotMatrix[2][2] = cosTheta * cosPhi;

        /*--- Loop over each point on the periodic marker that this rank
         holds locally and find the matching point from the donor marker. ---*/

        for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
          iPoint = vertex[iMarker][iVertex]->GetNode();
          iPointGlobal = nodes->GetGlobalIndex(iPoint);

          if (nodes->GetDomain(iPoint)) {
            /*--- Coordinates of the current boundary point ---*/

            const su2double* Coord_i = nodes->GetCoord(iPoint);

            /*--- Get the position vector from rotation center to point. ---*/

            dx = Coord_i[0] - center[0];
            dy = Coord_i[1] - center[1];
            dz = su2double(0.0);
            if (nDim == 3) dz = Coord_i[2] - center[2];

            /*--- Compute transformed point coordinates. ---*/

            rotCoord[0] = rotMatrix[0][0] * dx + rotMatrix[0][1] * dy + rotMatrix[0][2] * dz + translation[0];
            rotCoord[1] = rotMatrix[1][0] * dx + rotMatrix[1][1] * dy + rotMatrix[1][2] * dz + translation[1];
            rotCoord[2] = rotMatrix[2][0] * dx + rotMatrix[2][1] * dy + rotMatrix[2][2] * dz + translation[2];

            /*--- Check if the point lies on the axis of rotation. ---*/

            pointOnAxis = false;
            distToAxis = 0.0;
            for (iDim = 0; iDim < nDim; iDim++)
              distToAxis += (rotCoord[iDim] - Coord_i[iDim]) * (rotCoord[iDim] - Coord_i[iDim]);
            distToAxis = sqrt(distToAxis);

            if (distToAxis < epsilon) pointOnAxis = true;

            /*--- Initialize the distance to a large value. ---*/

            mindist = 1E6;
            pProcessor = 0;
            pPoint = 0;

            /*--- Loop over all of the periodic data that was gathered from
             all ranks in order to find the matching periodic point. ---*/

            for (iProcessor = 0; iProcessor < nProcessor; iProcessor++)
              for (jVertex = 0; jVertex < Buffer_Recv_nVertex[iProcessor]; jVertex++) {
                index = iProcessor * MaxLocalVertex_Periodic + jVertex;

                jPoint = Buffer_Recv_Point[index];
                jPointGlobal = Buffer_Recv_GlobalIndex[index];
                jVertex_ = Buffer_Recv_Vertex[index];
                jMarker = Buffer_Recv_Marker[index];

                /*--- Avoid matching the point to itself. Use local point
                 index + processor rather than global index, because global
                 indices are not set on coarse multigrid levels. ---*/

                if ((iProcessor != rank || jPoint != iPoint) || (pointOnAxis)) {
                  dist = 0.0;
                  for (iDim = 0; iDim < nDim; iDim++) {
                    Coord_j[iDim] = Buffer_Recv_Coord[index * nDim + iDim];
                    dist += pow(Coord_j[iDim] - rotCoord[iDim], 2.0);
                  }
                  dist = sqrt(dist);

                  chkSamePoint = (((dist < mindist) && (iProcessor != rank)) ||
                                  ((dist < mindist) && (iProcessor == rank) && (jPoint != iPoint)));

                  if (chkSamePoint || ((dist < mindist) && (pointOnAxis))) {
                    mindist = dist;
                    pProcessor = iProcessor;
                    pPoint = jPoint;
                    pPointGlobal = jPointGlobal;
                    pVertex = jVertex_;
                    pMarker = jMarker;
                  }
                }
              }

            /*--- Store the data for the best match found. ---*/

            vertex[iMarker][iVertex]->SetDonorPoint(pPoint, pPointGlobal, pVertex, pMarker, pProcessor);
            maxdist_local = max(maxdist_local, mindist);
            nPointMatch++;

            if (mindist > epsilon) {
              cout.precision(10);
              cout << endl;
              cout << "   Bad match for point " << iPointGlobal << ".\tNearest";
              cout << " donor distance: " << scientific << mindist << ".";
              maxdist_local = min(maxdist_local, 0.0);
              isBadMatch = true;
            }
          }
        }
      }
    }
  }

  /*--- Communicate the final count of matched points and max distance. ---*/

  unsigned long nPointMatch_Local = nPointMatch;
  SU2_MPI::Reduce(&nPointMatch_Local, &nPointMatch, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm());
  SU2_MPI::Reduce(&maxdist_local, &maxdist_global, 1, MPI_DOUBLE, MPI_MAX, MASTER_NODE, SU2_MPI::GetComm());

  if (rank == MASTER_NODE) {
    if (nPointMatch > 0) {
      cout << " Matched " << nPointMatch << " points with a max distance of: ";
      cout << maxdist_global << "." << endl;
    } else {
      cout << " No matching points for periodic marker pair ";
      cout << val_periodic << " in current zone." << endl;
    }

    if (isBadMatch) {
      cout << endl;
      cout << "\n !!! Warning !!!" << endl;
      cout << "Bad matches found. Computation will continue, but be cautious.\n";
    }
  }

  /*--- Free local memory for communications. ---*/

  delete[] Buffer_Send_Coord;
  delete[] Buffer_Send_Point;
  delete[] Buffer_Recv_Coord;
  delete[] Buffer_Recv_Point;
  delete[] Buffer_Send_nVertex;
  delete[] Buffer_Recv_nVertex;
  delete[] Buffer_Send_GlobalIndex;
  delete[] Buffer_Send_Vertex;
  delete[] Buffer_Send_Marker;
  delete[] Buffer_Recv_GlobalIndex;
  delete[] Buffer_Recv_Vertex;
  delete[] Buffer_Recv_Marker;
}

void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) {
  /*--- We start with the send and receive lists already available in
   the form of stored periodic point-donor pairs. We will loop through
   these markers and establish the neighboring ranks and number of
   send/recv points per pair. We will store this information and set
   up persistent data structures so that we can reuse them throughout
   the calculation for any periodic boundary communications. The goal
   is to break the non-blocking comms into InitiatePeriodicComms() and
   CompletePeriodicComms() in separate routines so that we can overlap the
   communication and computation to hide the communication latency. ---*/

  /*--- Local variables. ---*/

  unsigned short iMarker;
  unsigned long iPoint, iVertex, iPeriodic;

  int iRank, iSend, iRecv, ii, jj;

  /*--- Create some temporary structures for tracking sends/recvs. ---*/

  int* nPoint_Send_All = new int[size + 1];
  nPoint_Send_All[0] = 0;
  int* nPoint_Recv_All = new int[size + 1];
  nPoint_Recv_All[0] = 0;

  /*--- Store a set of unique (point, marker) pairs per destination rank. ---*/
  using PointMarkerPair = std::pair<unsigned long, unsigned long>;

  auto pairHash = [](const PointMarkerPair& p) -> std::size_t {
    // Use the golden ratio constant and bit mixing to generate pair hash
    std::size_t h1 = std::hash<unsigned long>{}(p.first);
    std::size_t h2 = std::hash<unsigned long>{}(p.second);

    return h1 ^ (h2 + 0x9e3779b9 + (h1 << 6) + (h1 >> 2));
  };

  using PointMarkerSet = std::unordered_set<PointMarkerPair, decltype(pairHash)>;
  std::vector<PointMarkerSet> Points_Send_All(size, PointMarkerSet(0, pairHash));

  /*--- Points on a rotation axis are their own periodic match. The other copies of their control
   volume are obtained by applying the rotation 1, 2, ... N-1 times (N = 360 deg / angle), instead
   of once in each direction. Therefore, the vertex on the first marker of the pair sends N-1 times
   to itself, and the vertex on the second marker does not send. ---*/

  const unsigned long nMarker_All = config->GetnMarker_All();
  PeriodicAxisPoints = false;

  /*--- Number of copies of the domain around the rotation axis of a periodic marker (N). ---*/

  auto nCopiesAroundAxis = [&](unsigned short jMarker) {
    /*--- Tolerance, in radians, to accept the rotation angle as an integer fraction of 360 degrees. ---*/
    constexpr passivedouble angleTol = 1e-4;

    /*--- The trace of a rotation matrix is 1 + 2 cos(angle). ---*/
    const auto markerTag = config->GetMarker_All_TagBound(jMarker);
    const su2double* angles = config->GetPeriodicRotAngles(markerTag);
    su2double rotMatrix[MAXNDIM][MAXNDIM];
    GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotMatrix);
    su2double cosAngle = 0.5 * (rotMatrix[0][0] + rotMatrix[1][1] + rotMatrix[2][2] - 1.0);
    if (nDim == 2) cosAngle = cos(angles[2]);
    cosAngle = min(cosAngle, 1.0);
    cosAngle = max(cosAngle, -1.0);
    const su2double angle = acos(cosAngle);

    const int nCopy = (angle > angleTol) ? SU2_TYPE::Int(2 * PI_NUMBER / angle + 0.5) : 0;

    if (fabs(nCopy * angle - 2 * PI_NUMBER) > angleTol) {
      SU2_MPI::Error(
          "Periodic points on the rotation axis require an angle that divides 360 degrees (marker " + markerTag + ").",
          CURRENT_FUNCTION);
    }
    return nCopy;
  };

  /*--- Number of sends of a vertex (1 for the points that are not on a rotation axis),
   and whether it is on a rotation axis. ---*/

  auto nSendOfVertex = [&](unsigned short jMarker, unsigned long jVertex, bool& onAxis) -> unsigned long {
    const auto* vert = geometry->vertex[jMarker][jVertex];
    onAxis = (static_cast<int>(vert->GetDonorProcessor()) == rank) &&
             (static_cast<unsigned long>(vert->GetDonorPoint()) == vert->GetNode());
    if (!onAxis) return 1;
    if (config->GetMarker_All_PerBound(jMarker) > config->GetnMarker_Periodic() / 2) return 0;

    PeriodicAxisPoints = true;
    return nCopiesAroundAxis(jMarker) - 1;
  };
  bool onAxis = false;

  /*--- Loop through all of our periodic markers and track
   our sends with each rank. ---*/

  for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
    if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
      iPeriodic = config->GetMarker_All_PerBound(iMarker);
      for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
        /*--- Get the current periodic point index. We only communicate
         the owned nodes on a rank, as the MPI comms will take care of
         the halos after completing the periodic comms. ---*/

        iPoint = geometry->vertex[iMarker][iVertex]->GetNode();

        if (geometry->nodes->GetDomain(iPoint)) {
          /*--- Get the rank that holds the matching periodic point
           on the other marker in the periodic pair. ---*/

          iRank = static_cast<int>(geometry->vertex[iMarker][iVertex]->GetDonorProcessor());

          /*--- Store the (point, marker) pair in the set for the destination rank,
           once for each time the vertex sends (see above). ---*/
          const auto nCopy = nSendOfVertex(iMarker, iVertex, onAxis);
          for (auto iCopy = 0ul; iCopy < nCopy; iCopy++)
            Points_Send_All[iRank].insert(std::make_pair(iPoint, iMarker + iCopy * nMarker_All));
        }
      }
    }
  }

  for (iRank = 0; iRank < size; iRank++) {
    nPoint_Send_All[iRank + 1] = Points_Send_All[iRank].size();
    nPoint_Recv_All[iRank + 1] = 0;
  }

  /*--- Communicate the number of points to be sent/recv'd amongst
   all processors. After this communication, each proc knows how
   many periodic points it will receive from each other processor. ---*/

  SU2_MPI::Alltoall(&(nPoint_Send_All[1]), 1, MPI_INT, &(nPoint_Recv_All[1]), 1, MPI_INT, SU2_MPI::GetComm());

  /*--- Check how many messages we will be sending and receiving.
   Here we also put the counters into cumulative storage format to
   make the communications simpler. Note that we are allowing each
   rank to communicate to themselves in these counters, although
   it will not be done through MPI. ---*/

  nPeriodicSend = 0;
  nPeriodicRecv = 0;

  for (iRank = 0; iRank < size; iRank++) {
    if ((nPoint_Send_All[iRank + 1] > 0)) nPeriodicSend++;
    if ((nPoint_Recv_All[iRank + 1] > 0)) nPeriodicRecv++;

    nPoint_Send_All[iRank + 1] += nPoint_Send_All[iRank];
    nPoint_Recv_All[iRank + 1] += nPoint_Recv_All[iRank];
  }

  /*--- Allocate only as much memory as needed for the periodic neighbors. ---*/

  nPoint_PeriodicSend = new int[nPeriodicSend + 1];
  nPoint_PeriodicSend[0] = 0;
  nPoint_PeriodicRecv = new int[nPeriodicRecv + 1];
  nPoint_PeriodicRecv[0] = 0;

  Neighbors_PeriodicSend = new int[nPeriodicSend];
  Neighbors_PeriodicRecv = new int[nPeriodicRecv];

  iSend = 0;
  iRecv = 0;
  for (iRank = 0; iRank < size; iRank++) {
    if ((nPoint_Send_All[iRank + 1] > nPoint_Send_All[iRank])) {
      Neighbors_PeriodicSend[iSend] = iRank;
      nPoint_PeriodicSend[iSend + 1] = nPoint_Send_All[iRank + 1];
      iSend++;
    }
    if ((nPoint_Recv_All[iRank + 1] > nPoint_Recv_All[iRank])) {
      Neighbors_PeriodicRecv[iRecv] = iRank;
      nPoint_PeriodicRecv[iRecv + 1] = nPoint_Recv_All[iRank + 1];
      iRecv++;
    }
  }

  /*--- Create a reverse mapping of the message to the rank so that we
   can quickly access the correct data in the buffers when receiving
   messages dynamically later during the iterations. ---*/

  PeriodicSend2Neighbor.clear();
  for (iSend = 0; iSend < nPeriodicSend; iSend++) PeriodicSend2Neighbor[Neighbors_PeriodicSend[iSend]] = iSend;

  PeriodicRecv2Neighbor.clear();
  for (iRecv = 0; iRecv < nPeriodicRecv; iRecv++) PeriodicRecv2Neighbor[Neighbors_PeriodicRecv[iRecv]] = iRecv;

  delete[] nPoint_Send_All;
  delete[] nPoint_Recv_All;

  /*--- Allocate the memory to store the local index values for both
   the send and receive periodic points and periodic index. ---*/

  Local_Point_PeriodicSend = nullptr;
  Local_Point_PeriodicSend = new unsigned long[nPoint_PeriodicSend[nPeriodicSend]];
  for (iSend = 0; iSend < nPoint_PeriodicSend[nPeriodicSend]; iSend++) Local_Point_PeriodicSend[iSend] = 0;

  Local_Marker_PeriodicSend = nullptr;
  Local_Marker_PeriodicSend = new unsigned long[nPoint_PeriodicSend[nPeriodicSend]];
  for (iSend = 0; iSend < nPoint_PeriodicSend[nPeriodicSend]; iSend++) Local_Marker_PeriodicSend[iSend] = 0;

  Local_Point_PeriodicRecv = nullptr;
  Local_Point_PeriodicRecv = new unsigned long[nPoint_PeriodicRecv[nPeriodicRecv]];
  for (iRecv = 0; iRecv < nPoint_PeriodicRecv[nPeriodicRecv]; iRecv++) Local_Point_PeriodicRecv[iRecv] = 0;

  Local_Marker_PeriodicRecv = nullptr;
  Local_Marker_PeriodicRecv = new unsigned long[nPoint_PeriodicRecv[nPeriodicRecv]];
  for (iRecv = 0; iRecv < nPoint_PeriodicRecv[nPeriodicRecv]; iRecv++) Local_Marker_PeriodicRecv[iRecv] = 0;

  Local_Copy_PeriodicSend.assign(nPoint_PeriodicSend[nPeriodicSend], 0);
  Local_Copy_PeriodicRecv.assign(nPoint_PeriodicRecv[nPeriodicRecv], 0);

  /*--- We allocate the buffers for communicating values in a later step
   once we know the maximum packet size that we need to communicate. This
   memory is deallocated and reallocated automatically in the case that
   the previously allocated memory is not sufficient. ---*/

  bufD_PeriodicSend = nullptr;
  bufD_PeriodicRecv = nullptr;

  bufS_PeriodicSend = nullptr;
  bufS_PeriodicRecv = nullptr;

  /*--- Allocate memory for the MPI requests if we need to communicate. ---*/

  if (nPeriodicSend > 0) {
    req_PeriodicSend = new SU2_MPI::Request[nPeriodicSend];
  }
  if (nPeriodicRecv > 0) {
    req_PeriodicRecv = new SU2_MPI::Request[nPeriodicRecv];
  }

  /*--- Allocate arrays for sending the periodic point index, the marker index,
   and the copy index (for points on a rotation axis) to the recv rank so that
   it can store the local values. Therefore, the recv rank can quickly loop
   through the buffers to unpack the data. ---*/

  const unsigned short nPackets = 3;
  auto* idSend = new unsigned long[nPoint_PeriodicSend[nPeriodicSend] * nPackets];
  for (iSend = 0; iSend < nPoint_PeriodicSend[nPeriodicSend] * nPackets; iSend++) idSend[iSend] = 0;

  /*--- Re-use set of unique points per destination rank. ---*/
  for (iRank = 0; iRank < size; iRank++) Points_Send_All[iRank].clear();

  /*--- Build the lists of local index and periodic marker index values. ---*/

  ii = 0;
  jj = 0;
  for (iSend = 0; iSend < nPeriodicSend; iSend++) {
    int destRank = Neighbors_PeriodicSend[iSend];
    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
      if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
        iPeriodic = config->GetMarker_All_PerBound(iMarker);
        for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
          /*--- Get the current periodic point index. We only communicate
           the owned nodes on a rank, as the MPI comms will take care of
           the halos after completing the periodic comms. ---*/

          iPoint = geometry->vertex[iMarker][iVertex]->GetNode();

          if (geometry->nodes->GetDomain(iPoint)) {
            /*--- Get the rank that holds the matching periodic point
             on the other marker in the periodic pair. ---*/

            iRank = static_cast<int>(geometry->vertex[iMarker][iVertex]->GetDonorProcessor());

            /*--- If the rank for the current periodic point matches the
             rank of the current send message, then store the local point
             index on the matching periodic point and the periodic marker
             index to be communicated to the recv rank. ---*/

            if (iRank == destRank) {
              /*--- Check if we have already added this (point, marker) pair
               for this destination rank. Use the result if insert(), which is
               a pair whose second element is success. Points on a rotation
               axis send more than once (see above). ---*/
              const auto nCopy = nSendOfVertex(iMarker, iVertex, onAxis);
              for (auto iCopy = 0ul; iCopy < nCopy; iCopy++) {
                const auto pointMarkerPair = std::make_pair(iPoint, iMarker + iCopy * nMarker_All);
                const auto insertResult = Points_Send_All[destRank].insert(pointMarkerPair);
                if (insertResult.second) {
                  Local_Point_PeriodicSend[ii] = iPoint;
                  Local_Marker_PeriodicSend[ii] = static_cast<unsigned long>(iMarker);
                  Local_Copy_PeriodicSend[ii] = onAxis ? iCopy + 1 : 0;
                  jj = ii * nPackets;
                  idSend[jj] = geometry->vertex[iMarker][iVertex]->GetDonorPoint();
                  jj++;
                  idSend[jj] = static_cast<unsigned long>(iPeriodic);
                  jj++;
                  idSend[jj] = Local_Copy_PeriodicSend[ii];
                  ii++;
                }
              }
            }
          }
        }
      }
    }
  }

  /*--- Allocate arrays for receiving the periodic point index and marker
   index to the recv rank so that it can store the local values. ---*/

  auto* idRecv = new unsigned long[nPoint_PeriodicRecv[nPeriodicRecv] * nPackets];
  for (iRecv = 0; iRecv < nPoint_PeriodicRecv[nPeriodicRecv] * nPackets; iRecv++) idRecv[iRecv] = 0;

#ifdef HAVE_MPI

  int iMessage, offset, count, source, dest, tag;

  /*--- Launch the non-blocking recv's first. Note that we have stored
   the counts and sources, so we can launch these before we even load
   the data and send from the periodically matching ranks. ---*/

  iMessage = 0;
  for (iRecv = 0; iRecv < nPeriodicRecv; iRecv++) {
    /*--- Compute our location in the recv buffer. ---*/

    offset = nPackets * nPoint_PeriodicRecv[iRecv];

    /*--- Take advantage of cumulative storage format to get the number
     of elems that we need to recv. ---*/

    count = nPackets * (nPoint_PeriodicRecv[iRecv + 1] - nPoint_PeriodicRecv[iRecv]);

    /*--- Get the rank from which we receive the message. ---*/

    source = Neighbors_PeriodicRecv[iRecv];
    tag = source + 1;

    /*--- Post non-blocking recv for this proc. ---*/

    SU2_MPI::Irecv(&(static_cast<unsigned long*>(idRecv)[offset]), count, MPI_UNSIGNED_LONG, source, tag,
                   SU2_MPI::GetComm(), &(req_PeriodicRecv[iMessage]));

    /*--- Increment message counter. ---*/

    iMessage++;
  }

  /*--- Post the non-blocking sends. ---*/

  iMessage = 0;
  for (iSend = 0; iSend < nPeriodicSend; iSend++) {
    /*--- Compute our location in the send buffer. ---*/

    offset = nPackets * nPoint_PeriodicSend[iSend];

    /*--- Take advantage of cumulative storage format to get the number
     of points that we need to send. ---*/

    count = nPackets * (nPoint_PeriodicSend[iSend + 1] - nPoint_PeriodicSend[iSend]);

    /*--- Get the rank to which we send the message. ---*/

    dest = Neighbors_PeriodicSend[iSend];
    tag = rank + 1;

    /*--- Post non-blocking send for this proc. ---*/

    SU2_MPI::Isend(&(static_cast<unsigned long*>(idSend)[offset]), count, MPI_UNSIGNED_LONG, dest, tag,
                   SU2_MPI::GetComm(), &(req_PeriodicSend[iMessage]));

    /*--- Increment message counter. ---*/

    iMessage++;
  }

  /*--- Wait for the non-blocking comms to complete. ---*/

  SU2_MPI::Waitall(nPeriodicSend, req_PeriodicSend, MPI_STATUS_IGNORE);
  SU2_MPI::Waitall(nPeriodicRecv, req_PeriodicRecv, MPI_STATUS_IGNORE);

#else

  /*--- Copy my own rank's data into the recv buffer directly in serial. ---*/

  int myStart, myFinal;
  for (int val_iSend = 0; val_iSend < nPeriodicSend; val_iSend++) {
    iRank = geometry->PeriodicRecv2Neighbor[rank];
    iRecv = geometry->nPoint_PeriodicRecv[iRank] * nPackets;
    myStart = nPoint_PeriodicSend[val_iSend] * nPackets;
    myFinal = nPoint_PeriodicSend[val_iSend + 1] * nPackets;
    for (iSend = myStart; iSend < myFinal; iSend++) {
      idRecv[iRecv] = idSend[iSend];
      iRecv++;
    }
  }

#endif

  /*--- Store the local periodic point and marker index values in our
   data structures so we can quickly unpack data during the iterations. ---*/

  ii = 0;
  for (iRecv = 0; iRecv < nPoint_PeriodicRecv[nPeriodicRecv]; iRecv++) {
    Local_Point_PeriodicRecv[iRecv] = idRecv[ii];
    ii++;
    Local_Marker_PeriodicRecv[iRecv] = idRecv[ii];
    ii++;
    Local_Copy_PeriodicRecv[iRecv] = idRecv[ii];
    ii++;
  }

  delete[] idSend;
  delete[] idRecv;
  ComputeModifiedSymmetryNormals(config);
}

void CGeometry::SumPeriodicGeometry(const CConfig* config, su2activematrix& values, int vectorIndex) {
  if (!nPeriodicSend && !nPeriodicRecv) return;
  const auto count = values.cols();
  AllocatePeriodicComms(count);
  for (auto iPair = 1u; iPair <= config->GetnMarker_Periodic() / 2; ++iPair) {
    PostPeriodicRecvs(this, config, COMM_TYPE::DOUBLE, count);
    for (auto iMessage = 0; iMessage < nPeriodicSend; ++iMessage) {
      for (auto iSend = nPoint_PeriodicSend[iMessage]; iSend < nPoint_PeriodicSend[iMessage + 1]; ++iSend) {
        const auto iPoint = Local_Point_PeriodicSend[iSend];
        auto* buffer = bufD_PeriodicSend + count * iSend;
        for (auto iVar = 0u; iVar < count; ++iVar) buffer[iVar] = values(iPoint, iVar);
        if (vectorIndex >= 0) {
          const auto* angles =
              config->GetPeriodicRotAngles(config->GetMarker_All_TagBound(Local_Marker_PeriodicSend[iSend]));
          su2double rotation[3][3], vector[3] = {};
          GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation);
          /*--- Axis points contribute the rotated normal of each additional copy. ---*/
          const auto nCopies = max(1ul, Local_Copy_PeriodicSend[iSend]);
          for (auto iCopy = 0ul; iCopy < nCopies; ++iCopy) {
            std::copy(buffer + vectorIndex, buffer + vectorIndex + nDim, vector);
            for (auto iDim = 0u; iDim < nDim; ++iDim) {
              buffer[vectorIndex + iDim] = 0;
              for (auto jDim = 0u; jDim < nDim; ++jDim)
                buffer[vectorIndex + iDim] += rotation[iDim][jDim] * vector[jDim];
            }
          }
        }
      }
#ifdef HAVE_MPI
      PostPeriodicSends(this, config, COMM_TYPE::DOUBLE, count, iMessage);
#else
      const auto message = PeriodicRecv2Neighbor[rank];
      const auto start = count * nPoint_PeriodicSend[iMessage];
      const auto end = count * nPoint_PeriodicSend[iMessage + 1];
      std::copy(bufD_PeriodicSend + start, bufD_PeriodicSend + end,
                bufD_PeriodicRecv + count * nPoint_PeriodicRecv[message]);
#endif
    }
    for (auto iMessage = 0; iMessage < nPeriodicRecv; ++iMessage) {
      auto message = iMessage;
#ifdef HAVE_MPI
      SU2_MPI::Status status;
      int index;
      SU2_MPI::Waitany(nPeriodicRecv, req_PeriodicRecv, &index, &status);
      message = PeriodicRecv2Neighbor[status.MPI_SOURCE];
#endif
      for (auto iRecv = nPoint_PeriodicRecv[message]; iRecv < nPoint_PeriodicRecv[message + 1]; ++iRecv) {
        const auto periodic = Local_Marker_PeriodicRecv[iRecv];
        if (periodic != iPair && periodic != iPair + config->GetnMarker_Periodic() / 2) continue;
        const auto iPoint = Local_Point_PeriodicRecv[iRecv];
        for (auto iVar = 0u; iVar < count; ++iVar) values(iPoint, iVar) += bufD_PeriodicRecv[count * iRecv + iVar];
      }
    }
    SU2_MPI::Waitall(nPeriodicSend, req_PeriodicSend, MPI_STATUS_IGNORE);
  }
}

void CGeometry::UpdatePeriodicVolumes(const CConfig* config) {
  if (!config->GetnMarker_Periodic() || (!nPeriodicSend && !nPeriodicRecv)) return;
  AllocatePeriodicComms(1);
  BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
    su2activematrix volumes(nPoint, 1);
    for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) volumes(iPoint, 0) = nodes->GetVolume(iPoint);
    SumPeriodicGeometry(config, volumes, -1);
    for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint)
      nodes->SetPeriodicVolume(iPoint, volumes(iPoint, 0) - nodes->GetVolume(iPoint));
  }
  END_SU2_OMP_SAFE_GLOBAL_ACCESS
}

void CGeometry::AllocatePeriodicComms(unsigned short countPerPeriodicPoint) {
  /*--- This routine is activated whenever we attempt to perform
   a periodic MPI communication with our neighbors but the
   memory buffer allocated is not large enough for the packet size.
   Therefore, we deallocate the previously allocated arrays and
   reallocate a large enough array. Note that after the first set
   communications, this routine will not need to be called again. ---*/

  if (countPerPeriodicPoint <= maxCountPerPeriodicPoint) return;

  BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
    /*--- Store the larger packet size to the class data. ---*/

    maxCountPerPeriodicPoint = countPerPeriodicPoint;

    /*--- Store the total size of the send/recv arrays for clarity. ---*/

    auto nSend = countPerPeriodicPoint * nPoint_PeriodicSend[nPeriodicSend];
    auto nRecv = countPerPeriodicPoint * nPoint_PeriodicRecv[nPeriodicRecv];

    /*-- Deallocate and reallocate our cummunication memory. ---*/

    delete[] bufD_PeriodicSend;
    bufD_PeriodicSend = new su2double[nSend]();

    delete[] bufD_PeriodicRecv;
    bufD_PeriodicRecv = new su2double[nRecv]();

    delete[] bufS_PeriodicSend;
    bufS_PeriodicSend = new unsigned short[nSend]();

    delete[] bufS_PeriodicRecv;
    bufS_PeriodicRecv = new unsigned short[nRecv]();
  }
  END_SU2_OMP_SAFE_GLOBAL_ACCESS
}

void CGeometry::PostPeriodicRecvs(CGeometry* geometry, const CConfig* config, COMM_TYPE commType,
                                  unsigned short countPerPeriodicPoint) {
  /*--- In parallel, communicate the data with non-blocking send/recv. ---*/

#ifdef HAVE_MPI

  /*--- Launch the non-blocking recv's first. Note that we have stored
   the counts and sources, so we can launch these before we even load
   the data and send from the neighbor ranks. ---*/

  SU2_OMP_MASTER
  for (int iRecv = 0; iRecv < nPeriodicRecv; iRecv++) {
    /*--- Compute our location in the recv buffer. ---*/

    auto offset = countPerPeriodicPoint * nPoint_PeriodicRecv[iRecv];

    /*--- Take advantage of cumulative storage format to get the number
     of elems that we need to recv. ---*/

    auto nPointPeriodic = nPoint_PeriodicRecv[iRecv + 1] - nPoint_PeriodicRecv[iRecv];

    /*--- Total count can include multiple pieces of data per element. ---*/

    auto count = countPerPeriodicPoint * nPointPeriodic;

    /*--- Get the rank from which we receive the message. ---*/

    auto source = Neighbors_PeriodicRecv[iRecv];
    auto tag = source + 1;

    /*--- Post non-blocking recv for this proc. ---*/

    switch (commType) {
      case COMM_TYPE::DOUBLE:
        SU2_MPI::Irecv(&(static_cast<su2double*>(bufD_PeriodicRecv)[offset]), count, MPI_DOUBLE, source, tag,
                       SU2_MPI::GetComm(), &(req_PeriodicRecv[iRecv]));
        break;
      case COMM_TYPE::UNSIGNED_SHORT:
        SU2_MPI::Irecv(&(static_cast<unsigned short*>(bufS_PeriodicRecv)[offset]), count, MPI_UNSIGNED_SHORT, source,
                       tag, SU2_MPI::GetComm(), &(req_PeriodicRecv[iRecv]));
        break;
      default:
        SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION);
        break;
    }
  }
  END_SU2_OMP_MASTER

#endif
}

void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, COMM_TYPE commType,
                                  unsigned short countPerPeriodicPoint, int val_iSend) const {
  /*--- In parallel, communicate the data with non-blocking send/recv. ---*/

#ifdef HAVE_MPI
  SU2_OMP_MASTER {
    /*--- Post the non-blocking send as soon as the buffer is loaded. ---*/

    /*--- Compute our location in the send buffer. ---*/

    auto offset = countPerPeriodicPoint * nPoint_PeriodicSend[val_iSend];

    /*--- Take advantage of cumulative storage format to get the number
     of points that we need to send. ---*/

    auto nPointPeriodic = (nPoint_PeriodicSend[val_iSend + 1] - nPoint_PeriodicSend[val_iSend]);

    /*--- Total count can include multiple pieces of data per element. ---*/

    auto count = countPerPeriodicPoint * nPointPeriodic;

    /*--- Get the rank to which we send the message. ---*/

    auto dest = Neighbors_PeriodicSend[val_iSend];
    auto tag = rank + 1;

    /*--- Post non-blocking send for this proc. ---*/

    switch (commType) {
      case COMM_TYPE::DOUBLE:
        SU2_MPI::Isend(&(static_cast<su2double*>(bufD_PeriodicSend)[offset]), count, MPI_DOUBLE, dest, tag,
                       SU2_MPI::GetComm(), &(req_PeriodicSend[val_iSend]));
        break;
      case COMM_TYPE::UNSIGNED_SHORT:
        SU2_MPI::Isend(&(static_cast<unsigned short*>(bufS_PeriodicSend)[offset]), count, MPI_UNSIGNED_SHORT, dest, tag,
                       SU2_MPI::GetComm(), &(req_PeriodicSend[val_iSend]));
        break;
      default:
        SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION);
        break;
    }
  }
  END_SU2_OMP_MASTER
#else

  /*--- Copy my own rank's data into the recv buffer directly in serial. ---*/

  int myStart, myFinal, iRecv, iRank;
  iRank = geometry->PeriodicRecv2Neighbor[rank];
  iRecv = geometry->nPoint_PeriodicRecv[iRank] * countPerPeriodicPoint;
  myStart = nPoint_PeriodicSend[val_iSend] * countPerPeriodicPoint;
  myFinal = nPoint_PeriodicSend[val_iSend + 1] * countPerPeriodicPoint;

  switch (commType) {
    case COMM_TYPE::DOUBLE:
      parallelCopy(myFinal - myStart, &bufD_PeriodicSend[myStart], &bufD_PeriodicRecv[iRecv]);
      break;
    case COMM_TYPE::UNSIGNED_SHORT:
      parallelCopy(myFinal - myStart, &bufS_PeriodicSend[myStart], &bufS_PeriodicRecv[iRecv]);
      break;
    default:
      SU2_MPI::Error("Unrecognized data type for periodic MPI comms.", CURRENT_FUNCTION);
      break;
  }

#endif
}

void CGeometry::SetEdges() {
  nEdge = 0;
  for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
    for (auto iNode = 0u; iNode < nodes->GetnPoint(iPoint); iNode++) {
      auto jPoint = nodes->GetPoint(iPoint, iNode);
      for (auto jNode = 0u; jNode < nodes->GetnPoint(jPoint); jNode++) {
        if (nodes->GetPoint(jPoint, jNode) == iPoint) {
          auto TestEdge = nodes->GetEdge(jPoint, jNode);
          if (TestEdge == -1) {
            nodes->SetEdge(iPoint, nEdge, iNode);
            nodes->SetEdge(jPoint, nEdge, jNode);
            nEdge++;
          }
          break;
        }
      }
    }
  }

  edges = new CEdge(nEdge, nDim);

  for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
    for (auto jPoint : nodes->GetPoints(iPoint)) {
      if (iPoint < jPoint) {
        auto iEdge = FindEdge(iPoint, jPoint);
        edges->SetNodes(iEdge, iPoint, jPoint);
      }
    }
  }
  edges->SetPaddingNodes();
}

void CGeometry::SetFaces() {
  //  unsigned long iPoint, jPoint, iFace;
  //  unsigned short jNode, iNode;
  //  long TestFace = 0;
  //
  //  nFace = 0;
  //  for (iPoint = 0; iPoint < nPoint; iPoint++)
  //    for (iNode = 0; iNode < nodes->GetnPoint(iPoint); iNode++) {
  //      jPoint = nodes->GetPoint(iPoint, iNode);
  //      for (jNode = 0; jNode < nodes->GetnPoint(jPoint); jNode++)
  //        if (nodes->GetPoint(jPoint, jNode) == iPoint) {
  //          TestFace = nodes->GetFace(jPoint, jNode);
  //          break;
  //        }
  //      if (TestFace == -1) {
  //        nodes->SetFace(iPoint, nFace, iNode);
  //        nodes->SetFace(jPoint, nFace, jNode);
  //        nFace++;
  //      }
  //    }
  //
  //  face = new CFace*[nFace];
  //
  //  for (iPoint = 0; iPoint < nPoint; iPoint++)
  //    for (iNode = 0; iNode < nodes->GetnPoint(iPoint); iNode++) {
  //      jPoint = nodes->GetPoint(iPoint, iNode);
  //      iFace = FindFace(iPoint, jPoint);
  //      if (iPoint < jPoint) face[iFace] = new CFace(iPoint, jPoint, nDim);
  //    }
}

void CGeometry::TestGeometry() const {
  ofstream para_file;

  para_file.open("test_geometry.dat", ios::out);

  auto* Normal = new su2double[nDim];

  for (unsigned long iEdge = 0; iEdge < nEdge; iEdge++) {
    para_file << "Edge index: " << iEdge << endl;
    para_file << "   Point index: " << edges->GetNode(iEdge, 0) << "\t" << edges->GetNode(iEdge, 1) << endl;
    edges->GetNormal(iEdge, Normal);
    para_file << "      Face normal : ";
    for (unsigned short iDim = 0; iDim < nDim; iDim++) para_file << Normal[iDim] << "\t";
    para_file << endl;
  }

  para_file << endl;
  para_file << endl;
  para_file << endl;
  para_file << endl;

  for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) {
    para_file << "Marker index: " << iMarker << endl;
    for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
      para_file << "   Vertex index: " << iVertex << endl;
      para_file << "      Point index: " << vertex[iMarker][iVertex]->GetNode() << endl;
      para_file << "      Point coordinates : ";
      for (unsigned short iDim = 0; iDim < nDim; iDim++) {
        para_file << nodes->GetCoord(vertex[iMarker][iVertex]->GetNode(), iDim) << "\t";
      }
      para_file << endl;
      vertex[iMarker][iVertex]->GetNormal(Normal);
      para_file << "         Face normal : ";
      for (unsigned short iDim = 0; iDim < nDim; iDim++) para_file << Normal[iDim] << "\t";
      para_file << endl;
    }
  }

  delete[] Normal;
}

bool CGeometry::SegmentIntersectsPlane(const su2double* Segment_P0, const su2double* Segment_P1, su2double Variable_P0,
                                       su2double Variable_P1, const su2double* Plane_P0, const su2double* Plane_Normal,
                                       su2double* Intersection, su2double& Variable_Interp) {
  su2double u[3], v[3], Denominator, Numerator, Aux, ModU;
  su2double epsilon =
      1E-6;  // An epsilon is added to eliminate, as much as possible, the posibility of a line that intersects a point
  unsigned short iDim;

  for (iDim = 0; iDim < 3; iDim++) {
    u[iDim] = Segment_P1[iDim] - Segment_P0[iDim];
    v[iDim] = (Plane_P0[iDim] + epsilon) - Segment_P0[iDim];
  }

  ModU = sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);

  Numerator =
      (Plane_Normal[0] + epsilon) * v[0] + (Plane_Normal[1] + epsilon) * v[1] + (Plane_Normal[2] + epsilon) * v[2];
  Denominator =
      (Plane_Normal[0] + epsilon) * u[0] + (Plane_Normal[1] + epsilon) * u[1] + (Plane_Normal[2] + epsilon) * u[2];

  if (fabs(Denominator) <= 0.0) return (false);  // No intersection.

  Aux = Numerator / Denominator;

  if (Aux < 0.0 || Aux > 1.0) return (false);  // No intersection.

  for (iDim = 0; iDim < 3; iDim++) Intersection[iDim] = Segment_P0[iDim] + Aux * u[iDim];

  /*--- Check that the intersection is in the segment ---*/

  for (iDim = 0; iDim < 3; iDim++) {
    u[iDim] = Segment_P0[iDim] - Intersection[iDim];
    v[iDim] = Segment_P1[iDim] - Intersection[iDim];
  }

  Variable_Interp = Variable_P0 + (Variable_P1 - Variable_P0) * sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]) / ModU;

  Denominator =
      (Plane_Normal[0] + epsilon) * u[0] + (Plane_Normal[1] + epsilon) * u[1] + (Plane_Normal[2] + epsilon) * u[2];
  Numerator =
      (Plane_Normal[0] + epsilon) * v[0] + (Plane_Normal[1] + epsilon) * v[1] + (Plane_Normal[2] + epsilon) * v[2];

  Aux = Numerator * Denominator;

  if (Aux > 0.0) return (false);  // Intersection outside the segment.

  return (true);
}

bool CGeometry::RayIntersectsTriangle(const su2double orig[3], const su2double dir[3], const su2double vert0[3],
                                      const su2double vert1[3], const su2double vert2[3], su2double* intersect) {
  const passivedouble epsilon = 0.000001;
  su2double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
  su2double det, inv_det, t, u, v;

  /*--- Find vectors for two edges sharing vert0 ---*/

  GeometryToolbox::Distance(3, vert1, vert0, edge1);
  GeometryToolbox::Distance(3, vert2, vert0, edge2);

  /*--- Begin calculating determinant - also used to calculate U parameter ---*/

  GeometryToolbox::CrossProduct(dir, edge2, pvec);

  /*--- If determinant is near zero, ray lies in plane of triangle ---*/

  det = GeometryToolbox::DotProduct(3, edge1, pvec);

  if (fabs(det) < epsilon) return (false);

  inv_det = 1.0 / det;

  /*--- Calculate distance from vert0 to ray origin ---*/

  GeometryToolbox::Distance(3, orig, vert0, tvec);

  /*--- Calculate U parameter and test bounds ---*/

  u = inv_det * GeometryToolbox::DotProduct(3, tvec, pvec);

  if (u < 0.0 || u > 1.0) return (false);

  /*--- prepare to test V parameter ---*/

  GeometryToolbox::CrossProduct(tvec, edge1, qvec);

  /*--- Calculate V parameter and test bounds ---*/

  v = inv_det * GeometryToolbox::DotProduct(3, dir, qvec);

  if (v < 0.0 || u + v > 1.0) return (false);

  /*--- Calculate t, ray intersects triangle ---*/

  t = inv_det * GeometryToolbox::DotProduct(3, edge2, qvec);

  /*--- Compute the intersection point in cartesian coordinates ---*/

  intersect[0] = orig[0] + (t * dir[0]);
  intersect[1] = orig[1] + (t * dir[1]);
  intersect[2] = orig[2] + (t * dir[2]);

  return (true);
}

bool CGeometry::SegmentIntersectsLine(const su2double point0[2], const su2double point1[2], const su2double vert0[2],
                                      const su2double vert1[2]) {
  su2double det, diff0_A, diff0_B, diff1_A, diff1_B, intersect[2];

  diff0_A = point0[0] - point1[0];
  diff1_A = point0[1] - point1[1];

  diff0_B = vert0[0] - vert1[0];
  diff1_B = vert0[1] - vert1[1];

  det = (diff0_A) * (diff1_B) - (diff1_A) * (diff0_B);

  if (det == 0) return false;

  /*--- Compute point of intersection ---*/

  intersect[0] = ((point0[0] * point1[1] - point0[1] * point1[0]) * diff0_B -
                  (vert0[0] * vert1[1] - vert0[1] * vert1[0]) * diff0_A) /
                 det;

  intersect[1] = ((point0[0] * point1[1] - point0[1] * point1[0]) * diff1_B -
                  (vert0[0] * vert1[1] - vert0[1] * vert1[0]) * diff1_A) /
                 det;

  /*--- Check that the point is between the two surface points ---*/

  su2double dist0, dist1, length;

  dist0 =
      (intersect[0] - point0[0]) * (intersect[0] - point0[0]) + (intersect[1] - point0[1]) * (intersect[1] - point0[1]);

  dist1 =
      (intersect[0] - point1[0]) * (intersect[0] - point1[0]) + (intersect[1] - point1[1]) * (intersect[1] - point1[1]);

  length = diff0_A * diff0_A + diff1_A * diff1_A;

  return (dist0 <= length) && (dist1 <= length);
}

bool CGeometry::SegmentIntersectsTriangle(su2double point0[3], const su2double point1[3], su2double vert0[3],
                                          su2double vert1[3], su2double vert2[3]) {
  su2double dir[3], intersect[3], u[3], v[3], edge1[3], edge2[3], Plane_Normal[3], Denominator, Numerator, Aux;

  GeometryToolbox::Distance(3, point1, point0, dir);

  if (RayIntersectsTriangle(point0, dir, vert0, vert1, vert2, intersect)) {
    /*--- Check that the intersection is in the segment ---*/

    GeometryToolbox::Distance(3, point0, intersect, u);
    GeometryToolbox::Distance(3, point1, intersect, v);

    GeometryToolbox::Distance(3, vert1, vert0, edge1);
    GeometryToolbox::Distance(3, vert2, vert0, edge2);
    GeometryToolbox::CrossProduct(edge1, edge2, Plane_Normal);

    Denominator = GeometryToolbox::DotProduct(3, Plane_Normal, u);
    Numerator = GeometryToolbox::DotProduct(3, Plane_Normal, v);

    Aux = Numerator * Denominator;

    /*--- Intersection outside the segment ---*/

    if (Aux > 0.0) return (false);

  } else {
    /*--- No intersection with the ray ---*/

    return (false);
  }

  /*--- Intersection inside the segment ---*/

  return (true);
}

void CGeometry::ComputeAirfoil_Section(su2double* Plane_P0, su2double* Plane_Normal, su2double MinXCoord,
                                       su2double MaxXCoord, su2double MinYCoord, su2double MaxYCoord,
                                       su2double MinZCoord, su2double MaxZCoord, const su2double* FlowVariable,
                                       vector<su2double>& Xcoord_Airfoil, vector<su2double>& Ycoord_Airfoil,
                                       vector<su2double>& Zcoord_Airfoil, vector<su2double>& Variable_Airfoil,
                                       bool original_surface, CConfig* config) {
  const bool wasActive = AD::BeginPassive();

  unsigned short iMarker, iNode, jNode, iDim, Index = 0;
  bool intersect;
  long Next_Edge = 0;
  unsigned long iPoint, jPoint, iElem, Trailing_Point, Airfoil_Point, iVertex, iEdge, PointIndex, jEdge;
  su2double Segment_P0[3] = {0.0, 0.0, 0.0}, Segment_P1[3] = {0.0, 0.0, 0.0}, Variable_P0 = 0.0, Variable_P1 = 0.0,
            Intersection[3] = {0.0, 0.0, 0.0}, Trailing_Coord, *VarCoord = nullptr, Variable_Interp,
            v1[3] = {0.0, 0.0, 0.0}, v3[3] = {0.0, 0.0, 0.0}, CrossProduct = 1.0;
  bool Found_Edge;
  passivedouble Dist_Value;
  vector<su2double> Xcoord_Index0, Ycoord_Index0, Zcoord_Index0, Variable_Index0, Xcoord_Index1, Ycoord_Index1,
      Zcoord_Index1, Variable_Index1;
  vector<unsigned long> IGlobalID_Index0, JGlobalID_Index0, IGlobalID_Index1, JGlobalID_Index1, IGlobalID_Airfoil,
      JGlobalID_Airfoil;
  vector<unsigned short> Conection_Index0, Conection_Index1;
  vector<unsigned long> Duplicate;
  su2double** Coord_Variation = nullptr;
  vector<su2double> XcoordExtra, YcoordExtra, ZcoordExtra, VariableExtra;
  vector<unsigned long> IGlobalIDExtra, JGlobalIDExtra;
  vector<bool> AddExtra;
  unsigned long EdgeDonor;
  bool FoundEdge;

#ifdef HAVE_MPI
  unsigned long nLocalEdge, MaxLocalEdge, *Buffer_Send_nEdge, *Buffer_Receive_nEdge, nBuffer_Coord, nBuffer_Variable,
      nBuffer_GlobalID;
  int nProcessor, iProcessor;
  su2double *Buffer_Send_Coord, *Buffer_Receive_Coord;
  su2double *Buffer_Send_Variable, *Buffer_Receive_Variable;
  unsigned long *Buffer_Send_GlobalID, *Buffer_Receive_GlobalID;
#endif

  Xcoord_Airfoil.clear();
  Ycoord_Airfoil.clear();
  Zcoord_Airfoil.clear();
  Variable_Airfoil.clear();
  IGlobalID_Airfoil.clear();
  JGlobalID_Airfoil.clear();

  /*--- Set the right plane in 2D (note the change in Y-Z plane) ---*/

  if (nDim == 2) {
    Plane_P0[0] = 0.0;
    Plane_P0[1] = 0.0;
    Plane_P0[2] = 0.0;
    Plane_Normal[0] = 0.0;
    Plane_Normal[1] = 1.0;
    Plane_Normal[2] = 0.0;
  }

  /*--- Grid movement is stored using a vertices information,
   we should go from vertex to points ---*/

  if (!original_surface) {
    Coord_Variation = new su2double*[nPoint];
    for (iPoint = 0; iPoint < nPoint; iPoint++) Coord_Variation[iPoint] = new su2double[nDim];

    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
      if (config->GetMarker_All_GeoEval(iMarker) == YES) {
        for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
          VarCoord = vertex[iMarker][iVertex]->GetVarCoord();
          iPoint = vertex[iMarker][iVertex]->GetNode();
          for (iDim = 0; iDim < nDim; iDim++) Coord_Variation[iPoint][iDim] = VarCoord[iDim];
        }
      }
    }
  }

  for (iMarker = 0; iMarker < nMarker; iMarker++) {
    if (config->GetMarker_All_GeoEval(iMarker) == YES) {
      for (iElem = 0; iElem < nElem_Bound[iMarker]; iElem++) {
        PointIndex = 0;

        /*--- To decide if an element is going to be used or not should be done element based,
         The first step is to compute and average coordinate for the element ---*/

        su2double AveXCoord = 0.0;
        su2double AveYCoord = 0.0;
        su2double AveZCoord = 0.0;

        for (iNode = 0; iNode < bound[iMarker][iElem]->GetnNodes(); iNode++) {
          iPoint = bound[iMarker][iElem]->GetNode(iNode);
          AveXCoord += nodes->GetCoord(iPoint, 0);
          AveYCoord += nodes->GetCoord(iPoint, 1);
          if (nDim == 3) AveZCoord += nodes->GetCoord(iPoint, 2);
        }

        AveXCoord /= su2double(bound[iMarker][iElem]->GetnNodes());
        AveYCoord /= su2double(bound[iMarker][iElem]->GetnNodes());
        AveZCoord /= su2double(bound[iMarker][iElem]->GetnNodes());

        /*--- To only cut one part of the nacelle based on the cross product
         of the normal to the plane and a vector that connect the point
         with the center line ---*/

        CrossProduct = 1.0;

        if (config->GetGeo_Description() == NACELLE) {
          su2double Tilt_Angle = config->GetNacelleLocation(3) * PI_NUMBER / 180;
          su2double Toe_Angle = config->GetNacelleLocation(4) * PI_NUMBER / 180;

          /*--- Translate to the origin ---*/

          su2double XCoord_Trans = AveXCoord - config->GetNacelleLocation(0);
          su2double YCoord_Trans = AveYCoord - config->GetNacelleLocation(1);
          su2double ZCoord_Trans = AveZCoord - config->GetNacelleLocation(2);

          /*--- Apply tilt angle ---*/

          su2double XCoord_Trans_Tilt = XCoord_Trans * cos(Tilt_Angle) + ZCoord_Trans * sin(Tilt_Angle);
          su2double YCoord_Trans_Tilt = YCoord_Trans;
          su2double ZCoord_Trans_Tilt = ZCoord_Trans * cos(Tilt_Angle) - XCoord_Trans * sin(Tilt_Angle);

          /*--- Apply toe angle ---*/

          su2double YCoord_Trans_Tilt_Toe = XCoord_Trans_Tilt * sin(Toe_Angle) + YCoord_Trans_Tilt * cos(Toe_Angle);
          su2double ZCoord_Trans_Tilt_Toe = ZCoord_Trans_Tilt;

          /*--- Undo plane rotation, we have already rotated the nacelle ---*/

          /*--- Undo tilt angle ---*/

          su2double XPlane_Normal_Tilt = Plane_Normal[0] * cos(-Tilt_Angle) + Plane_Normal[2] * sin(-Tilt_Angle);
          su2double YPlane_Normal_Tilt = Plane_Normal[1];
          su2double ZPlane_Normal_Tilt = Plane_Normal[2] * cos(-Tilt_Angle) - Plane_Normal[0] * sin(-Tilt_Angle);

          /*--- Undo toe angle ---*/

          su2double YPlane_Normal_Tilt_Toe =
              XPlane_Normal_Tilt * sin(-Toe_Angle) + YPlane_Normal_Tilt * cos(-Toe_Angle);
          su2double ZPlane_Normal_Tilt_Toe = ZPlane_Normal_Tilt;

          v1[1] = YCoord_Trans_Tilt_Toe - 0.0;
          v1[2] = ZCoord_Trans_Tilt_Toe - 0.0;
          v3[0] = v1[1] * ZPlane_Normal_Tilt_Toe - v1[2] * YPlane_Normal_Tilt_Toe;
          CrossProduct = v3[0] * 1.0;
        }

        for (unsigned short iFace = 0; iFace < bound[iMarker][iElem]->GetnFaces(); iFace++) {
          iNode = bound[iMarker][iElem]->GetFaces(iFace, 0);
          jNode = bound[iMarker][iElem]->GetFaces(iFace, 1);
          iPoint = bound[iMarker][iElem]->GetNode(iNode);
          jPoint = bound[iMarker][iElem]->GetNode(jNode);

          if ((CrossProduct >= 0.0) && ((AveXCoord > MinXCoord) && (AveXCoord < MaxXCoord)) &&
              ((AveYCoord > MinYCoord) && (AveYCoord < MaxYCoord)) &&
              ((AveZCoord > MinZCoord) && (AveZCoord < MaxZCoord))) {
            Segment_P0[0] = 0.0;
            Segment_P0[1] = 0.0;
            Segment_P0[2] = 0.0;
            Variable_P0 = 0.0;
            Segment_P1[0] = 0.0;
            Segment_P1[1] = 0.0;
            Segment_P1[2] = 0.0;
            Variable_P1 = 0.0;

            for (iDim = 0; iDim < nDim; iDim++) {
              if (original_surface) {
                Segment_P0[iDim] = nodes->GetCoord(iPoint, iDim);
                Segment_P1[iDim] = nodes->GetCoord(jPoint, iDim);
              } else {
                Segment_P0[iDim] = nodes->GetCoord(iPoint, iDim) + Coord_Variation[iPoint][iDim];
                Segment_P1[iDim] = nodes->GetCoord(jPoint, iDim) + Coord_Variation[jPoint][iDim];
              }
            }

            if (FlowVariable != nullptr) {
              Variable_P0 = FlowVariable[iPoint];
              Variable_P1 = FlowVariable[jPoint];
            }

            /*--- In 2D add the points directly (note the change between Y and Z coordinate) ---*/

            if (nDim == 2) {
              Xcoord_Index0.push_back(Segment_P0[0]);
              Xcoord_Index1.push_back(Segment_P1[0]);
              Ycoord_Index0.push_back(Segment_P0[2]);
              Ycoord_Index1.push_back(Segment_P1[2]);
              Zcoord_Index0.push_back(Segment_P0[1]);
              Zcoord_Index1.push_back(Segment_P1[1]);
              Variable_Index0.push_back(Variable_P0);
              Variable_Index1.push_back(Variable_P1);
              IGlobalID_Index0.push_back(nodes->GetGlobalIndex(iPoint));
              IGlobalID_Index1.push_back(nodes->GetGlobalIndex(jPoint));
              JGlobalID_Index0.push_back(nodes->GetGlobalIndex(iPoint));
              JGlobalID_Index1.push_back(nodes->GetGlobalIndex(jPoint));
              PointIndex++;
            }

            /*--- In 3D compute the intersection ---*/

            else if (nDim == 3) {
              intersect = SegmentIntersectsPlane(Segment_P0, Segment_P1, Variable_P0, Variable_P1, Plane_P0,
                                                 Plane_Normal, Intersection, Variable_Interp);
              if (intersect) {
                if (PointIndex == 0) {
                  Xcoord_Index0.push_back(Intersection[0]);
                  Ycoord_Index0.push_back(Intersection[1]);
                  Zcoord_Index0.push_back(Intersection[2]);
                  Variable_Index0.push_back(Variable_Interp);
                  IGlobalID_Index0.push_back(nodes->GetGlobalIndex(iPoint));
                  JGlobalID_Index0.push_back(nodes->GetGlobalIndex(jPoint));
                }
                if (PointIndex == 1) {
                  Xcoord_Index1.push_back(Intersection[0]);
                  Ycoord_Index1.push_back(Intersection[1]);
                  Zcoord_Index1.push_back(Intersection[2]);
                  Variable_Index1.push_back(Variable_Interp);
                  IGlobalID_Index1.push_back(nodes->GetGlobalIndex(iPoint));
                  JGlobalID_Index1.push_back(nodes->GetGlobalIndex(jPoint));
                }
                PointIndex++;
              }
            }
          }
        }
      }
    }
  }

  if (!original_surface) {
    for (iPoint = 0; iPoint < nPoint; iPoint++) delete[] Coord_Variation[iPoint];
    delete[] Coord_Variation;
  }

#ifdef HAVE_MPI

  /*--- Copy the coordinates of all the points in the plane to the master node ---*/

  nLocalEdge = 0, MaxLocalEdge = 0;
  nProcessor = size;

  Buffer_Send_nEdge = new unsigned long[1];
  Buffer_Receive_nEdge = new unsigned long[nProcessor];

  nLocalEdge = Xcoord_Index0.size();

  Buffer_Send_nEdge[0] = nLocalEdge;

  SU2_MPI::Allreduce(&nLocalEdge, &MaxLocalEdge, 1, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_nEdge, 1, MPI_UNSIGNED_LONG, Buffer_Receive_nEdge, 1, MPI_UNSIGNED_LONG,
                     SU2_MPI::GetComm());

  Buffer_Send_Coord = new su2double[MaxLocalEdge * 6];
  Buffer_Receive_Coord = new su2double[nProcessor * MaxLocalEdge * 6];

  Buffer_Send_Variable = new su2double[MaxLocalEdge * 2];
  Buffer_Receive_Variable = new su2double[nProcessor * MaxLocalEdge * 2];

  Buffer_Send_GlobalID = new unsigned long[MaxLocalEdge * 4];
  Buffer_Receive_GlobalID = new unsigned long[nProcessor * MaxLocalEdge * 4];

  nBuffer_Coord = MaxLocalEdge * 6;
  nBuffer_Variable = MaxLocalEdge * 2;
  nBuffer_GlobalID = MaxLocalEdge * 4;

  for (iEdge = 0; iEdge < nLocalEdge; iEdge++) {
    Buffer_Send_Coord[iEdge * 6 + 0] = Xcoord_Index0[iEdge];
    Buffer_Send_Coord[iEdge * 6 + 1] = Ycoord_Index0[iEdge];
    Buffer_Send_Coord[iEdge * 6 + 2] = Zcoord_Index0[iEdge];
    Buffer_Send_Coord[iEdge * 6 + 3] = Xcoord_Index1[iEdge];
    Buffer_Send_Coord[iEdge * 6 + 4] = Ycoord_Index1[iEdge];
    Buffer_Send_Coord[iEdge * 6 + 5] = Zcoord_Index1[iEdge];

    Buffer_Send_Variable[iEdge * 2 + 0] = Variable_Index0[iEdge];
    Buffer_Send_Variable[iEdge * 2 + 1] = Variable_Index1[iEdge];

    Buffer_Send_GlobalID[iEdge * 4 + 0] = IGlobalID_Index0[iEdge];
    Buffer_Send_GlobalID[iEdge * 4 + 1] = JGlobalID_Index0[iEdge];
    Buffer_Send_GlobalID[iEdge * 4 + 2] = IGlobalID_Index1[iEdge];
    Buffer_Send_GlobalID[iEdge * 4 + 3] = JGlobalID_Index1[iEdge];
  }

  SU2_MPI::Allgather(Buffer_Send_Coord, nBuffer_Coord, MPI_DOUBLE, Buffer_Receive_Coord, nBuffer_Coord, MPI_DOUBLE,
                     SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_Variable, nBuffer_Variable, MPI_DOUBLE, Buffer_Receive_Variable, nBuffer_Variable,
                     MPI_DOUBLE, SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_GlobalID, nBuffer_GlobalID, MPI_UNSIGNED_LONG, Buffer_Receive_GlobalID,
                     nBuffer_GlobalID, MPI_UNSIGNED_LONG, SU2_MPI::GetComm());

  /*--- Clean the vectors before adding the new vertices only to the master node ---*/

  Xcoord_Index0.clear();
  Xcoord_Index1.clear();
  Ycoord_Index0.clear();
  Ycoord_Index1.clear();
  Zcoord_Index0.clear();
  Zcoord_Index1.clear();
  Variable_Index0.clear();
  Variable_Index1.clear();
  IGlobalID_Index0.clear();
  IGlobalID_Index1.clear();
  JGlobalID_Index0.clear();
  JGlobalID_Index1.clear();

  /*--- Copy the boundary to the master node vectors ---*/

  if (rank == MASTER_NODE) {
    for (iProcessor = 0; iProcessor < nProcessor; iProcessor++) {
      for (iEdge = 0; iEdge < Buffer_Receive_nEdge[iProcessor]; iEdge++) {
        Xcoord_Index0.push_back(Buffer_Receive_Coord[iProcessor * MaxLocalEdge * 6 + iEdge * 6 + 0]);
        Ycoord_Index0.push_back(Buffer_Receive_Coord[iProcessor * MaxLocalEdge * 6 + iEdge * 6 + 1]);
        Zcoord_Index0.push_back(Buffer_Receive_Coord[iProcessor * MaxLocalEdge * 6 + iEdge * 6 + 2]);
        Xcoord_Index1.push_back(Buffer_Receive_Coord[iProcessor * MaxLocalEdge * 6 + iEdge * 6 + 3]);
        Ycoord_Index1.push_back(Buffer_Receive_Coord[iProcessor * MaxLocalEdge * 6 + iEdge * 6 + 4]);
        Zcoord_Index1.push_back(Buffer_Receive_Coord[iProcessor * MaxLocalEdge * 6 + iEdge * 6 + 5]);

        Variable_Index0.push_back(Buffer_Receive_Variable[iProcessor * MaxLocalEdge * 2 + iEdge * 2 + 0]);
        Variable_Index1.push_back(Buffer_Receive_Variable[iProcessor * MaxLocalEdge * 2 + iEdge * 2 + 1]);

        IGlobalID_Index0.push_back(Buffer_Receive_GlobalID[iProcessor * MaxLocalEdge * 4 + iEdge * 4 + 0]);
        JGlobalID_Index0.push_back(Buffer_Receive_GlobalID[iProcessor * MaxLocalEdge * 4 + iEdge * 4 + 1]);
        IGlobalID_Index1.push_back(Buffer_Receive_GlobalID[iProcessor * MaxLocalEdge * 4 + iEdge * 4 + 2]);
        JGlobalID_Index1.push_back(Buffer_Receive_GlobalID[iProcessor * MaxLocalEdge * 4 + iEdge * 4 + 3]);
      }
    }
  }

  delete[] Buffer_Send_Coord;
  delete[] Buffer_Receive_Coord;
  delete[] Buffer_Send_Variable;
  delete[] Buffer_Receive_Variable;
  delete[] Buffer_Send_GlobalID;
  delete[] Buffer_Receive_GlobalID;
  delete[] Buffer_Send_nEdge;
  delete[] Buffer_Receive_nEdge;

#endif

  if ((rank == MASTER_NODE) && (!Xcoord_Index0.empty())) {
    /*--- Remove singular edges ---*/

    bool Remove;

    do {
      Remove = false;
      for (iEdge = 0; iEdge < Xcoord_Index0.size(); iEdge++) {
        if (((IGlobalID_Index0[iEdge] == IGlobalID_Index1[iEdge]) &&
             (JGlobalID_Index0[iEdge] == JGlobalID_Index1[iEdge])) ||
            ((IGlobalID_Index0[iEdge] == JGlobalID_Index1[iEdge]) &&
             (JGlobalID_Index0[iEdge] == IGlobalID_Index1[iEdge]))) {
          Xcoord_Index0.erase(Xcoord_Index0.begin() + iEdge);
          Ycoord_Index0.erase(Ycoord_Index0.begin() + iEdge);
          Zcoord_Index0.erase(Zcoord_Index0.begin() + iEdge);
          Variable_Index0.erase(Variable_Index0.begin() + iEdge);
          IGlobalID_Index0.erase(IGlobalID_Index0.begin() + iEdge);
          JGlobalID_Index0.erase(JGlobalID_Index0.begin() + iEdge);

          Xcoord_Index1.erase(Xcoord_Index1.begin() + iEdge);
          Ycoord_Index1.erase(Ycoord_Index1.begin() + iEdge);
          Zcoord_Index1.erase(Zcoord_Index1.begin() + iEdge);
          Variable_Index1.erase(Variable_Index1.begin() + iEdge);
          IGlobalID_Index1.erase(IGlobalID_Index1.begin() + iEdge);
          JGlobalID_Index1.erase(JGlobalID_Index1.begin() + iEdge);

          Remove = true;
          break;
        }
        if (Remove) break;
      }
    } while (Remove);

    /*--- Remove repeated edges computing distance, this could happend because the MPI ---*/

    do {
      Remove = false;
      for (iEdge = 0; iEdge < Xcoord_Index0.size() - 1; iEdge++) {
        for (jEdge = iEdge + 1; jEdge < Xcoord_Index0.size(); jEdge++) {
          /*--- Edges with the same orientation ---*/

          if ((((IGlobalID_Index0[iEdge] == IGlobalID_Index0[jEdge]) &&
                (JGlobalID_Index0[iEdge] == JGlobalID_Index0[jEdge])) ||
               ((IGlobalID_Index0[iEdge] == JGlobalID_Index0[jEdge]) &&
                (JGlobalID_Index0[iEdge] == IGlobalID_Index0[jEdge]))) &&
              (((IGlobalID_Index1[iEdge] == IGlobalID_Index1[jEdge]) &&
                (JGlobalID_Index1[iEdge] == JGlobalID_Index1[jEdge])) ||
               ((IGlobalID_Index1[iEdge] == JGlobalID_Index1[jEdge]) &&
                (JGlobalID_Index1[iEdge] == IGlobalID_Index1[jEdge])))) {
            Xcoord_Index0.erase(Xcoord_Index0.begin() + jEdge);
            Ycoord_Index0.erase(Ycoord_Index0.begin() + jEdge);
            Zcoord_Index0.erase(Zcoord_Index0.begin() + jEdge);
            Variable_Index0.erase(Variable_Index0.begin() + jEdge);
            IGlobalID_Index0.erase(IGlobalID_Index0.begin() + jEdge);
            JGlobalID_Index0.erase(JGlobalID_Index0.begin() + jEdge);

            Xcoord_Index1.erase(Xcoord_Index1.begin() + jEdge);
            Ycoord_Index1.erase(Ycoord_Index1.begin() + jEdge);
            Zcoord_Index1.erase(Zcoord_Index1.begin() + jEdge);
            Variable_Index1.erase(Variable_Index1.begin() + jEdge);
            IGlobalID_Index1.erase(IGlobalID_Index1.begin() + jEdge);
            JGlobalID_Index1.erase(JGlobalID_Index1.begin() + jEdge);

            Remove = true;
            break;
          }

          /*--- Edges with oposite orientation ---*/

          if ((((IGlobalID_Index0[iEdge] == IGlobalID_Index1[jEdge]) &&
                (JGlobalID_Index0[iEdge] == JGlobalID_Index1[jEdge])) ||
               ((IGlobalID_Index0[iEdge] == JGlobalID_Index1[jEdge]) &&
                (JGlobalID_Index0[iEdge] == IGlobalID_Index1[jEdge]))) &&
              (((IGlobalID_Index1[iEdge] == IGlobalID_Index0[jEdge]) &&
                (JGlobalID_Index1[iEdge] == JGlobalID_Index0[jEdge])) ||
               ((IGlobalID_Index1[iEdge] == JGlobalID_Index0[jEdge]) &&
                (JGlobalID_Index1[iEdge] == IGlobalID_Index0[jEdge])))) {
            Xcoord_Index0.erase(Xcoord_Index0.begin() + jEdge);
            Ycoord_Index0.erase(Ycoord_Index0.begin() + jEdge);
            Zcoord_Index0.erase(Zcoord_Index0.begin() + jEdge);
            Variable_Index0.erase(Variable_Index0.begin() + jEdge);
            IGlobalID_Index0.erase(IGlobalID_Index0.begin() + jEdge);
            JGlobalID_Index0.erase(JGlobalID_Index0.begin() + jEdge);

            Xcoord_Index1.erase(Xcoord_Index1.begin() + jEdge);
            Ycoord_Index1.erase(Ycoord_Index1.begin() + jEdge);
            Zcoord_Index1.erase(Zcoord_Index1.begin() + jEdge);
            Variable_Index1.erase(Variable_Index1.begin() + jEdge);
            IGlobalID_Index1.erase(IGlobalID_Index1.begin() + jEdge);
            JGlobalID_Index1.erase(JGlobalID_Index1.begin() + jEdge);

            Remove = true;
            break;
          }
          if (Remove) break;
        }
        if (Remove) break;
      }

    } while (Remove);

    if (Xcoord_Index0.size() != 1) {
      /*--- Rotate from the Y-Z plane to the X-Z plane to reuse the rest of subroutines  ---*/

      if (config->GetGeo_Description() == FUSELAGE) {
        su2double Angle = -0.5 * PI_NUMBER;
        for (iEdge = 0; iEdge < Xcoord_Index0.size(); iEdge++) {
          su2double XCoord = Xcoord_Index0[iEdge] * cos(Angle) - Ycoord_Index0[iEdge] * sin(Angle);
          su2double YCoord = Ycoord_Index0[iEdge] * cos(Angle) + Xcoord_Index0[iEdge] * sin(Angle);
          su2double ZCoord = Zcoord_Index0[iEdge];
          Xcoord_Index0[iEdge] = XCoord;
          Ycoord_Index0[iEdge] = YCoord;
          Zcoord_Index0[iEdge] = ZCoord;
          XCoord = Xcoord_Index1[iEdge] * cos(Angle) - Ycoord_Index1[iEdge] * sin(Angle);
          YCoord = Ycoord_Index1[iEdge] * cos(Angle) + Xcoord_Index1[iEdge] * sin(Angle);
          ZCoord = Zcoord_Index1[iEdge];
          Xcoord_Index1[iEdge] = XCoord;
          Ycoord_Index1[iEdge] = YCoord;
          Zcoord_Index1[iEdge] = ZCoord;
        }
      }

      /*--- Rotate nacelle secction to a X-Z plane to reuse the rest of subroutines  ---*/

      if (config->GetGeo_Description() == NACELLE) {
        su2double Tilt_Angle = config->GetNacelleLocation(3) * PI_NUMBER / 180;
        su2double Toe_Angle = config->GetNacelleLocation(4) * PI_NUMBER / 180;
        su2double Theta_deg = atan2(Plane_Normal[1], -Plane_Normal[2]) / PI_NUMBER * 180 + 180;
        su2double Roll_Angle = 0.5 * PI_NUMBER - Theta_deg * PI_NUMBER / 180;

        su2double XCoord_Trans, YCoord_Trans, ZCoord_Trans, XCoord_Trans_Tilt, YCoord_Trans_Tilt, ZCoord_Trans_Tilt,
            XCoord_Trans_Tilt_Toe, YCoord_Trans_Tilt_Toe, ZCoord_Trans_Tilt_Toe, XCoord, YCoord, ZCoord;

        for (iEdge = 0; iEdge < Xcoord_Index0.size(); iEdge++) {
          /*--- First point of the edge ---*/

          /*--- Translate to the origin ---*/

          XCoord_Trans = Xcoord_Index0[iEdge] - config->GetNacelleLocation(0);
          YCoord_Trans = Ycoord_Index0[iEdge] - config->GetNacelleLocation(1);
          ZCoord_Trans = Zcoord_Index0[iEdge] - config->GetNacelleLocation(2);

          /*--- Apply tilt angle ---*/

          XCoord_Trans_Tilt = XCoord_Trans * cos(Tilt_Angle) + ZCoord_Trans * sin(Tilt_Angle);
          YCoord_Trans_Tilt = YCoord_Trans;
          ZCoord_Trans_Tilt = ZCoord_Trans * cos(Tilt_Angle) - XCoord_Trans * sin(Tilt_Angle);

          /*--- Apply toe angle ---*/

          XCoord_Trans_Tilt_Toe = XCoord_Trans_Tilt * cos(Toe_Angle) - YCoord_Trans_Tilt * sin(Toe_Angle);
          YCoord_Trans_Tilt_Toe = XCoord_Trans_Tilt * sin(Toe_Angle) + YCoord_Trans_Tilt * cos(Toe_Angle);
          ZCoord_Trans_Tilt_Toe = ZCoord_Trans_Tilt;

          /*--- Rotate to X-Z plane (roll) ---*/

          XCoord = XCoord_Trans_Tilt_Toe;
          YCoord = YCoord_Trans_Tilt_Toe * cos(Roll_Angle) - ZCoord_Trans_Tilt_Toe * sin(Roll_Angle);
          ZCoord = YCoord_Trans_Tilt_Toe * sin(Roll_Angle) + ZCoord_Trans_Tilt_Toe * cos(Roll_Angle);

          /*--- Update coordinates ---*/

          Xcoord_Index0[iEdge] = XCoord;
          Ycoord_Index0[iEdge] = YCoord;
          Zcoord_Index0[iEdge] = ZCoord;

          /*--- Second point of the edge ---*/

          /*--- Translate to the origin ---*/

          XCoord_Trans = Xcoord_Index1[iEdge] - config->GetNacelleLocation(0);
          YCoord_Trans = Ycoord_Index1[iEdge] - config->GetNacelleLocation(1);
          ZCoord_Trans = Zcoord_Index1[iEdge] - config->GetNacelleLocation(2);

          /*--- Apply tilt angle ---*/

          XCoord_Trans_Tilt = XCoord_Trans * cos(Tilt_Angle) + ZCoord_Trans * sin(Tilt_Angle);
          YCoord_Trans_Tilt = YCoord_Trans;
          ZCoord_Trans_Tilt = ZCoord_Trans * cos(Tilt_Angle) - XCoord_Trans * sin(Tilt_Angle);

          /*--- Apply toe angle ---*/

          XCoord_Trans_Tilt_Toe = XCoord_Trans_Tilt * cos(Toe_Angle) - YCoord_Trans_Tilt * sin(Toe_Angle);
          YCoord_Trans_Tilt_Toe = XCoord_Trans_Tilt * sin(Toe_Angle) + YCoord_Trans_Tilt * cos(Toe_Angle);
          ZCoord_Trans_Tilt_Toe = ZCoord_Trans_Tilt;

          /*--- Rotate to X-Z plane (roll) ---*/

          XCoord = XCoord_Trans_Tilt_Toe;
          YCoord = YCoord_Trans_Tilt_Toe * cos(Roll_Angle) - ZCoord_Trans_Tilt_Toe * sin(Roll_Angle);
          ZCoord = YCoord_Trans_Tilt_Toe * sin(Roll_Angle) + ZCoord_Trans_Tilt_Toe * cos(Roll_Angle);

          /*--- Update coordinates ---*/

          Xcoord_Index1[iEdge] = XCoord;
          Ycoord_Index1[iEdge] = YCoord;
          Zcoord_Index1[iEdge] = ZCoord;
        }
      }

      /*--- Identify the extreme of the curve and close it ---*/

      Conection_Index0.reserve(Xcoord_Index0.size() + 1);
      Conection_Index1.reserve(Xcoord_Index0.size() + 1);

      for (iEdge = 0; iEdge < Xcoord_Index0.size(); iEdge++) {
        Conection_Index0[iEdge] = 0;
        Conection_Index1[iEdge] = 0;
      }

      for (iEdge = 0; iEdge < Xcoord_Index0.size() - 1; iEdge++) {
        for (jEdge = iEdge + 1; jEdge < Xcoord_Index0.size(); jEdge++) {
          if (((IGlobalID_Index0[iEdge] == IGlobalID_Index0[jEdge]) &&
               (JGlobalID_Index0[iEdge] == JGlobalID_Index0[jEdge])) ||
              ((IGlobalID_Index0[iEdge] == JGlobalID_Index0[jEdge]) &&
               (JGlobalID_Index0[iEdge] == IGlobalID_Index0[jEdge]))) {
            Conection_Index0[iEdge]++;
            Conection_Index0[jEdge]++;
          }

          if (((IGlobalID_Index0[iEdge] == IGlobalID_Index1[jEdge]) &&
               (JGlobalID_Index0[iEdge] == JGlobalID_Index1[jEdge])) ||
              ((IGlobalID_Index0[iEdge] == JGlobalID_Index1[jEdge]) &&
               (JGlobalID_Index0[iEdge] == IGlobalID_Index1[jEdge]))) {
            Conection_Index0[iEdge]++;
            Conection_Index1[jEdge]++;
          }

          if (((IGlobalID_Index1[iEdge] == IGlobalID_Index0[jEdge]) &&
               (JGlobalID_Index1[iEdge] == JGlobalID_Index0[jEdge])) ||
              ((IGlobalID_Index1[iEdge] == JGlobalID_Index0[jEdge]) &&
               (JGlobalID_Index1[iEdge] == IGlobalID_Index0[jEdge]))) {
            Conection_Index1[iEdge]++;
            Conection_Index0[jEdge]++;
          }

          if (((IGlobalID_Index1[iEdge] == IGlobalID_Index1[jEdge]) &&
               (JGlobalID_Index1[iEdge] == JGlobalID_Index1[jEdge])) ||
              ((IGlobalID_Index1[iEdge] == JGlobalID_Index1[jEdge]) &&
               (JGlobalID_Index1[iEdge] == IGlobalID_Index1[jEdge]))) {
            Conection_Index1[iEdge]++;
            Conection_Index1[jEdge]++;
          }
        }
      }

      /*--- Connect extremes of the curves ---*/

      /*--- First: Identify the extremes of the curve in the extra vector  ---*/

      for (iEdge = 0; iEdge < Xcoord_Index0.size(); iEdge++) {
        if (Conection_Index0[iEdge] == 0) {
          XcoordExtra.push_back(Xcoord_Index0[iEdge]);
          YcoordExtra.push_back(Ycoord_Index0[iEdge]);
          ZcoordExtra.push_back(Zcoord_Index0[iEdge]);
          VariableExtra.push_back(Variable_Index0[iEdge]);
          IGlobalIDExtra.push_back(IGlobalID_Index0[iEdge]);
          JGlobalIDExtra.push_back(JGlobalID_Index0[iEdge]);
          AddExtra.push_back(true);
        }
        if (Conection_Index1[iEdge] == 0) {
          XcoordExtra.push_back(Xcoord_Index1[iEdge]);
          YcoordExtra.push_back(Ycoord_Index1[iEdge]);
          ZcoordExtra.push_back(Zcoord_Index1[iEdge]);
          VariableExtra.push_back(Variable_Index1[iEdge]);
          IGlobalIDExtra.push_back(IGlobalID_Index1[iEdge]);
          JGlobalIDExtra.push_back(JGlobalID_Index1[iEdge]);
          AddExtra.push_back(true);
        }
      }

      /*--- Second, if it is an open curve then find the closest point to an extreme to close it  ---*/

      if (XcoordExtra.size() > 1) {
        for (iEdge = 0; iEdge < XcoordExtra.size() - 1; iEdge++) {
          su2double MinDist = 1E6;
          FoundEdge = false;
          EdgeDonor = 0;
          for (jEdge = iEdge + 1; jEdge < XcoordExtra.size(); jEdge++) {
            Dist_Value =
                sqrt(pow(SU2_TYPE::GetValue(XcoordExtra[iEdge]) - SU2_TYPE::GetValue(XcoordExtra[jEdge]), 2.0));
            if ((Dist_Value < MinDist) && (AddExtra[iEdge]) && (AddExtra[jEdge])) {
              EdgeDonor = jEdge;
              FoundEdge = true;
            }
          }

          if (FoundEdge) {
            /*--- Add first point of the new edge ---*/

            Xcoord_Index0.push_back(XcoordExtra[iEdge]);
            Ycoord_Index0.push_back(YcoordExtra[iEdge]);
            Zcoord_Index0.push_back(ZcoordExtra[iEdge]);
            Variable_Index0.push_back(VariableExtra[iEdge]);
            IGlobalID_Index0.push_back(IGlobalIDExtra[iEdge]);
            JGlobalID_Index0.push_back(JGlobalIDExtra[iEdge]);
            AddExtra[iEdge] = false;

            /*--- Add second (closest)  point of the new edge ---*/

            Xcoord_Index1.push_back(XcoordExtra[EdgeDonor]);
            Ycoord_Index1.push_back(YcoordExtra[EdgeDonor]);
            Zcoord_Index1.push_back(ZcoordExtra[EdgeDonor]);
            Variable_Index1.push_back(VariableExtra[EdgeDonor]);
            IGlobalID_Index1.push_back(IGlobalIDExtra[EdgeDonor]);
            JGlobalID_Index1.push_back(JGlobalIDExtra[EdgeDonor]);
            AddExtra[EdgeDonor] = false;
          }
        }

      }

      else if (XcoordExtra.size() == 1) {
        cout << "There cutting system has failed, there is an incomplete curve (not used)." << endl;
      }

      /*--- Find and add the trailing edge to to the list
       and the contect the first point to the trailing edge ---*/

      Trailing_Point = 0;
      Trailing_Coord = Xcoord_Index0[0];
      for (iEdge = 1; iEdge < Xcoord_Index0.size(); iEdge++) {
        if (Xcoord_Index0[iEdge] > Trailing_Coord) {
          Trailing_Point = iEdge;
          Trailing_Coord = Xcoord_Index0[iEdge];
        }
      }

      Xcoord_Airfoil.push_back(Xcoord_Index0[Trailing_Point]);
      Ycoord_Airfoil.push_back(Ycoord_Index0[Trailing_Point]);
      Zcoord_Airfoil.push_back(Zcoord_Index0[Trailing_Point]);
      Variable_Airfoil.push_back(Variable_Index0[Trailing_Point]);
      IGlobalID_Airfoil.push_back(IGlobalID_Index0[Trailing_Point]);
      JGlobalID_Airfoil.push_back(JGlobalID_Index0[Trailing_Point]);

      Xcoord_Airfoil.push_back(Xcoord_Index1[Trailing_Point]);
      Ycoord_Airfoil.push_back(Ycoord_Index1[Trailing_Point]);
      Zcoord_Airfoil.push_back(Zcoord_Index1[Trailing_Point]);
      Variable_Airfoil.push_back(Variable_Index1[Trailing_Point]);
      IGlobalID_Airfoil.push_back(IGlobalID_Index1[Trailing_Point]);
      JGlobalID_Airfoil.push_back(JGlobalID_Index1[Trailing_Point]);

      Xcoord_Index0.erase(Xcoord_Index0.begin() + Trailing_Point);
      Ycoord_Index0.erase(Ycoord_Index0.begin() + Trailing_Point);
      Zcoord_Index0.erase(Zcoord_Index0.begin() + Trailing_Point);
      Variable_Index0.erase(Variable_Index0.begin() + Trailing_Point);
      IGlobalID_Index0.erase(IGlobalID_Index0.begin() + Trailing_Point);
      JGlobalID_Index0.erase(JGlobalID_Index0.begin() + Trailing_Point);

      Xcoord_Index1.erase(Xcoord_Index1.begin() + Trailing_Point);
      Ycoord_Index1.erase(Ycoord_Index1.begin() + Trailing_Point);
      Zcoord_Index1.erase(Zcoord_Index1.begin() + Trailing_Point);
      Variable_Index1.erase(Variable_Index1.begin() + Trailing_Point);
      IGlobalID_Index1.erase(IGlobalID_Index1.begin() + Trailing_Point);
      JGlobalID_Index1.erase(JGlobalID_Index1.begin() + Trailing_Point);

      /*--- Algorithm for adding the rest of the points ---*/

      do {
        /*--- Last added point in the list ---*/

        Airfoil_Point = Xcoord_Airfoil.size() - 1;

        /*--- Find the closest point  ---*/

        Found_Edge = false;

        for (iEdge = 0; iEdge < Xcoord_Index0.size(); iEdge++) {
          if (((IGlobalID_Index0[iEdge] == IGlobalID_Airfoil[Airfoil_Point]) &&
               (JGlobalID_Index0[iEdge] == JGlobalID_Airfoil[Airfoil_Point])) ||
              ((IGlobalID_Index0[iEdge] == JGlobalID_Airfoil[Airfoil_Point]) &&
               (JGlobalID_Index0[iEdge] == IGlobalID_Airfoil[Airfoil_Point]))) {
            Next_Edge = iEdge;
            Found_Edge = true;
            Index = 0;
            break;
          }

          if (((IGlobalID_Index1[iEdge] == IGlobalID_Airfoil[Airfoil_Point]) &&
               (JGlobalID_Index1[iEdge] == JGlobalID_Airfoil[Airfoil_Point])) ||
              ((IGlobalID_Index1[iEdge] == JGlobalID_Airfoil[Airfoil_Point]) &&
               (JGlobalID_Index1[iEdge] == IGlobalID_Airfoil[Airfoil_Point]))) {
            Next_Edge = iEdge;
            Found_Edge = true;
            Index = 1;
            break;
          }
        }

        /*--- Add and remove the next point to the list and the next point in the edge ---*/

        if (Found_Edge) {
          if (Index == 0) {
            Xcoord_Airfoil.push_back(Xcoord_Index1[Next_Edge]);
            Ycoord_Airfoil.push_back(Ycoord_Index1[Next_Edge]);
            Zcoord_Airfoil.push_back(Zcoord_Index1[Next_Edge]);
            Variable_Airfoil.push_back(Variable_Index1[Next_Edge]);
            IGlobalID_Airfoil.push_back(IGlobalID_Index1[Next_Edge]);
            JGlobalID_Airfoil.push_back(JGlobalID_Index1[Next_Edge]);
          }

          if (Index == 1) {
            Xcoord_Airfoil.push_back(Xcoord_Index0[Next_Edge]);
            Ycoord_Airfoil.push_back(Ycoord_Index0[Next_Edge]);
            Zcoord_Airfoil.push_back(Zcoord_Index0[Next_Edge]);
            Variable_Airfoil.push_back(Variable_Index0[Next_Edge]);
            IGlobalID_Airfoil.push_back(IGlobalID_Index0[Next_Edge]);
            JGlobalID_Airfoil.push_back(JGlobalID_Index0[Next_Edge]);
          }

          Xcoord_Index0.erase(Xcoord_Index0.begin() + Next_Edge);
          Ycoord_Index0.erase(Ycoord_Index0.begin() + Next_Edge);
          Zcoord_Index0.erase(Zcoord_Index0.begin() + Next_Edge);
          Variable_Index0.erase(Variable_Index0.begin() + Next_Edge);
          IGlobalID_Index0.erase(IGlobalID_Index0.begin() + Next_Edge);
          JGlobalID_Index0.erase(JGlobalID_Index0.begin() + Next_Edge);

          Xcoord_Index1.erase(Xcoord_Index1.begin() + Next_Edge);
          Ycoord_Index1.erase(Ycoord_Index1.begin() + Next_Edge);
          Zcoord_Index1.erase(Zcoord_Index1.begin() + Next_Edge);
          Variable_Index1.erase(Variable_Index1.begin() + Next_Edge);
          IGlobalID_Index1.erase(IGlobalID_Index1.begin() + Next_Edge);
          JGlobalID_Index1.erase(JGlobalID_Index1.begin() + Next_Edge);

        } else {
          break;
        }

      } while (!Xcoord_Index0.empty());

      /*--- Clean the vector before using them again for storing the upper or the lower side ---*/

      Xcoord_Index0.clear();
      Ycoord_Index0.clear();
      Zcoord_Index0.clear();
      Variable_Index0.clear();
      IGlobalID_Index0.clear();
      JGlobalID_Index0.clear();
      Xcoord_Index1.clear();
      Ycoord_Index1.clear();
      Zcoord_Index1.clear();
      Variable_Index1.clear();
      IGlobalID_Index1.clear();
      JGlobalID_Index1.clear();
    }
  }

  AD::EndPassive(wasActive);
}

void CGeometry::RegisterCoordinates() const {
  const bool input = true;

  SU2_OMP_FOR_STAT(roundUpDiv(nPoint, omp_get_num_threads()))
  for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
    nodes->RegisterCoordinates(iPoint, input);
  }
  END_SU2_OMP_FOR
}

void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) {
  geometry_container[MESH_0]->InitiateComms(geometry_container[MESH_0], config, MPI_QUANTITIES::COORDINATES);
  geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, MPI_QUANTITIES::COORDINATES);

  geometry_container[MESH_0]->SetControlVolume(config, UPDATE);
  geometry_container[MESH_0]->SetBoundControlVolume(config, UPDATE);
  geometry_container[MESH_0]->SetMaxLength(config);
  geometry_container[MESH_0]->UpdatePeriodicVolumes(config);

  for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) {
    /*--- Update the control volume structures ---*/

    geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh - 1], UPDATE);
    geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], config, UPDATE);
    geometry_container[iMesh]->SetCoord(geometry_container[iMesh - 1]);
    geometry_container[iMesh]->UpdatePeriodicVolumes(config);
  }

  /*--- Compute the global surface areas for all markers. ---*/
  geometry_container[MESH_0]->ComputeSurfaceAreaCfgFile(config);
}

void CGeometry::SetCustomBoundary(CConfig* config) {
  unsigned short iMarker;
  unsigned long iVertex;
  string Marker_Tag;

  /*--- Initialize quantities for customized boundary conditions.
   * Custom values are initialized with the default values specified in the config (avoiding non physical values) ---*/
  CustomBoundaryTemperature = new su2double*[nMarker];
  CustomBoundaryHeatFlux = new su2double*[nMarker];

  for (iMarker = 0; iMarker < nMarker; iMarker++) {
    Marker_Tag = config->GetMarker_All_TagBound(iMarker);
    CustomBoundaryHeatFlux[iMarker] = nullptr;
    CustomBoundaryTemperature[iMarker] = nullptr;
    if (config->GetMarker_All_PyCustom(iMarker)) {
      switch (config->GetMarker_All_KindBC(iMarker)) {
        case HEAT_FLUX:
          CustomBoundaryHeatFlux[iMarker] = new su2double[nVertex[iMarker]];
          for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
            CustomBoundaryHeatFlux[iMarker][iVertex] = config->GetWall_HeatFlux(Marker_Tag);
          }
          break;
        case ISOTHERMAL:
          CustomBoundaryTemperature[iMarker] = new su2double[nVertex[iMarker]];
          for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
            CustomBoundaryTemperature[iMarker][iVertex] = config->GetIsothermal_Temperature(Marker_Tag);
          }
          break;
        case INLET_FLOW:
          // This case is handled in the solver class.
          break;
        default:
          cout << "WARNING: Marker " << Marker_Tag << " is not customizable. Using default behavior." << endl;
          break;
      }
    }
  }
}

void CGeometry::UpdateCustomBoundaryConditions(CGeometry** geometry_container, CConfig* config) {
  for (auto iMGlevel = 1u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) {
    const auto iMGfine = iMGlevel - 1;
    for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
      if (!config->GetMarker_All_PyCustom(iMarker)) continue;
      switch (config->GetMarker_All_KindBC(iMarker)) {
        case HEAT_FLUX:
          geometry_container[iMGlevel]->SetMultiGridWallHeatFlux(geometry_container[iMGfine], iMarker);
          break;
        case ISOTHERMAL:
          geometry_container[iMGlevel]->SetMultiGridWallTemperature(geometry_container[iMGfine], iMarker);
          break;
        // Inlet flow handled in solver class.
        default:
          break;
      }
    }
  }
}

void CGeometry::ComputeSurfaceAreaCfgFile(const CConfig* config){
    SU2_OMP_MASTER{const auto nMarker_Global = config->GetnMarker_CfgFile();
SurfaceAreaCfgFile.resize(nMarker_Global);
vector<su2double> LocalSurfaceArea(nMarker_Global, 0.0);

/*--- Loop over all local markers ---*/
for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) {
  const auto Local_TagBound = config->GetMarker_All_TagBound(iMarker);

  /*--- Loop over all global markers, and find the local-global pair via
        matching unique string tags. ---*/
  for (unsigned short iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
    const auto Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);
    if (Local_TagBound == Global_TagBound) {
      for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
        const auto iPoint = vertex[iMarker][iVertex]->GetNode();

        if (!nodes->GetDomain(iPoint)) continue;

        const auto AreaNormal = vertex[iMarker][iVertex]->GetNormal();
        const auto Area = GeometryToolbox::Norm(nDim, AreaNormal);

        LocalSurfaceArea[iMarker_Global] += Area;
      }  // for iVertex
    }    // if Local == Global
  }      // for iMarker_Global
}  // for iMarker

SU2_MPI::Allreduce(LocalSurfaceArea.data(), SurfaceAreaCfgFile.data(), SurfaceAreaCfgFile.size(), MPI_DOUBLE, MPI_SUM,
                   SU2_MPI::GetComm());
}
END_SU2_OMP_MASTER
}

su2double CGeometry::GetSurfaceArea(const CConfig* config, unsigned short val_marker) const {
  /*---Find the precomputed marker surface area by local-global string-matching. ---*/
  const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker);

  for (unsigned short iMarker_Global = 0; iMarker_Global < config->GetnMarker_CfgFile(); iMarker_Global++) {
    const auto Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);

    if (Marker_Tag == Global_TagBound) return SurfaceAreaCfgFile[iMarker_Global];
  }

  SU2_MPI::Error("Unable to match local-marker with cfg-marker for Surface Area.", CURRENT_FUNCTION);
  return 0.0;
}

void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
  const su2double MIN_AREA = 1e-12;            // minimum area to consider a normal valid
  const su2double PARALLEL_TOLERANCE = 0.001;  // min. angle between symm. planes to consider them non-parallel.

  /* Check how many symmetry planes there are and use the first (lowest ID) as the basis to orthogonalize against.
   * All nodes that are shared by multiple symmetries have to get a corrected normal. */

  /*--- Compute if markers are straight lines or planes. ---*/
  ComputeSurfStraightness(config, false);

  symmetryNormals.clear();
  symmetryNormals.resize(nMarker);
  std::vector<unsigned short> symMarkers, curvedSymMarkers;

  /*--- Check which markers are symmetry or Euler wall and store in a list. ---*/
  for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) {
    if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE) ||
        (config->GetMarker_All_KindBC(iMarker) == EULER_WALL)) {
      symMarkers.push_back(iMarker);
      if (!boundIsStraight[iMarker]) curvedSymMarkers.push_back(iMarker);
    }
  }

  std::vector<std::unordered_map<unsigned long, std::array<su2double, MAXNDIM>>> periodicNormals(nMarker);
  std::vector<bool> periodicPoints(nPoint, false);
  if (nPeriodicRecv)
    for (auto i = 0; i < nPoint_PeriodicRecv[nPeriodicRecv]; ++i) periodicPoints[Local_Point_PeriodicRecv[i]] = true;

  /*--- A slip wall meeting a periodic face needs the normal of the complete
   * wall patch. Keep markers separate to preserve distinct corner constraints. ---*/
  if (config->GetnMarker_Periodic() && (nPeriodicSend || nPeriodicRecv)) {
    for (auto iCfgMarker = 0u; iCfgMarker < config->GetnMarker_CfgFile(); ++iCfgMarker) {
      const auto tag = config->GetMarker_CfgFile_TagBound(iCfgMarker);
      const auto kind = config->GetMarker_CfgFile_KindBC(tag);
      if (kind != SYMMETRY_PLANE && kind != EULER_WALL) continue;
      su2activematrix normals(nPoint, nDim);
      normals = su2double(0);
      for (const auto iMarker : symMarkers) {
        if (config->GetMarker_All_TagBound(iMarker) != tag) continue;
        for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) {
          const auto iPoint = vertex[iMarker][iVertex]->GetNode();
          for (auto iDim = 0u; iDim < nDim; ++iDim) normals(iPoint, iDim) = vertex[iMarker][iVertex]->GetNormal(iDim);
        }
      }
      SumPeriodicGeometry(config, normals, 0);
      for (const auto iMarker : symMarkers) {
        if (config->GetMarker_All_TagBound(iMarker) != tag) continue;
        for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) {
          const auto iPoint = vertex[iMarker][iVertex]->GetNode();
          if (!periodicPoints[iPoint]) continue;
          const auto area = GeometryToolbox::Norm(nDim, normals[iPoint]);
          if (area <= MIN_AREA) continue;
          auto& normal = symmetryNormals[iMarker][iVertex];
          auto& areaNormal = periodicNormals[iMarker][iVertex];
          for (auto iDim = 0u; iDim < nDim; ++iDim) {
            areaNormal[iDim] = normals(iPoint, iDim);
            normal[iDim] = areaNormal[iDim] / area;
          }
        }
      }
    }
  }

  /*--- Merge the normals of curved markers to be independent of how surfaces are divided. ---*/

  std::unordered_map<unsigned long, std::array<su2double, MAXNDIM>> mergedNormals;

  for (const auto iMarker : curvedSymMarkers) {
    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
      const auto iPoint = vertex[iMarker][iVertex]->GetNode();

      /*--- Determine if this point is shared with other curved symmetries. ---*/
      int count = 0;
      for (const auto jMarker : curvedSymMarkers) {
        count += static_cast<int>(nodes->GetVertex(iPoint, jMarker) >= 0);
      }
      if (count < 2) continue;

      std::array<su2double, MAXNDIM> normal = {};
      const auto corrected = periodicNormals[iMarker].find(iVertex);
      if (corrected != periodicNormals[iMarker].end())
        normal = corrected->second;
      else
        vertex[iMarker][iVertex]->GetNormal(normal.data());

      auto result = mergedNormals.emplace(iPoint, normal);
      const auto inserted = result.second;
      auto it = result.first;
      if (!inserted) {
        for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) {
          it->second[iDim] += normal[iDim];
        }
      }
    }
  }

  for (const auto& item : mergedNormals) {
    const auto iPoint = item.first;
    auto normal = item.second;
    const su2double area = GeometryToolbox::Norm(int(MAXNDIM), normal.data());
    for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) normal[iDim] /= area;

    for (const auto iMarker : curvedSymMarkers) {
      const auto iVertex = nodes->GetVertex(iPoint, iMarker);
      if (iVertex >= 0) symmetryNormals[iMarker][iVertex] = normal;
    }
  }

  /*--- Now do Gramm-Schmidt Process for symmetries that are not both curved lines or planes. ---*/

  /*--- Loop over all markers and find nodes on symmetry planes that are shared with other symmetries. ---*/
  /*--- The first symmetry does not need a corrected normal vector, hence start at 1. ---*/
  for (size_t i = 1; i < symMarkers.size(); ++i) {
    const auto iMarker = symMarkers[i];

    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
      const auto iPoint = vertex[iMarker][iVertex]->GetNode();

      /*--- Halo points do not need to be considered. ---*/
      if (!nodes->GetDomain(iPoint)) continue;

      /*--- Get the vertex normal on the current symmetry, which may be a merged normal. ---*/
      auto GetNormal = [&](unsigned long iMarker, unsigned long iVertex, std::array<su2double, MAXNDIM>& normal) {
        const auto it = symmetryNormals[iMarker].find(iVertex);
        if (it != symmetryNormals[iMarker].end()) {
          normal = it->second;
        } else {
          vertex[iMarker][iVertex]->GetNormal(normal.data());
          const su2double jarea = GeometryToolbox::Norm(int(MAXNDIM), normal.data());
          for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) normal[iDim] /= jarea;
        }
      };
      std::array<su2double, MAXNDIM> iNormal = {};
      GetNormal(iMarker, iVertex, iNormal);

      /*--- Loop over previous symmetries and if this point shares them, make this normal orthogonal to them.
       * It's ok if we normalize merged normals against themselves, we get 0 area and this becomes a no-op. ---*/

      std::vector<size_t> parallelMarkers;  // Track markers with nearly-parallel normals

      for (size_t j = 0; j < i; ++j) {
        const auto jMarker = symMarkers[j];
        const auto jVertex = nodes->GetVertex(iPoint, jMarker);
        if (jVertex < 0) continue;

        std::array<su2double, MAXNDIM> jNormal = {};
        GetNormal(jMarker, jVertex, jNormal);

        const su2double proj = GeometryToolbox::DotProduct(int(MAXNDIM), jNormal.data(), iNormal.data());
        const su2double angleDiff = std::abs(1.0 - std::abs(proj));

        // Check if normals are nearly parallel (within ~2.5 degrees)
        // cos(2.5°) ≈ 0.999, so (1 - cos(2.5°)) ≈ 0.001
        if (angleDiff < PARALLEL_TOLERANCE) {
          // These normals are nearly parallel - average them instead of orthogonalizing
          parallelMarkers.push_back(j);
          for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] += jNormal[iDim];
          continue;
        }

        for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] -= proj * jNormal[iDim];
      }

      /*--- If we found parallel markers, average and store the result for all involved markers ---*/
      if (!parallelMarkers.empty()) {
        // Normalize the averaged normal
        const su2double avgArea = GeometryToolbox::Norm(int(MAXNDIM), iNormal.data());
        if (avgArea > MIN_AREA) {
          for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] /= avgArea;

          // Store the averaged normal for the current marker
          symmetryNormals[iMarker][iVertex] = iNormal;

          // Also update all parallel markers with the same averaged normal
          for (const auto j : parallelMarkers) {
            const auto jMarker = symMarkers[j];
            const auto jVertex = nodes->GetVertex(iPoint, jMarker);
            if (jVertex >= 0) {
              symmetryNormals[jMarker][jVertex] = iNormal;
            }
          }
        }
        continue;  // Skip the normal orthogonalization path below
      }

      /*--- Normalize. If the norm is close to zero it means the normal is a linear combination of previous
       * normals, in this case we don't need to store the corrected normal, using the original in the gradient
       * correction will have no effect since previous corrections will remove components in this direction). ---*/
      const su2double area = GeometryToolbox::Norm(int(MAXNDIM), iNormal.data());

      if (area > MIN_AREA) {
        for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] /= area;
        symmetryNormals[iMarker][iVertex] = iNormal;
      }
    }
  }
}

void CGeometry::ComputeSurfStraightness(const CConfig* config, bool print_on_screen) {
  bool RefUnitNormal_defined;
  unsigned short iDim, iMarker, iMarker_Global, nMarker_Global = config->GetnMarker_CfgFile();
  unsigned long iVertex;
  constexpr passivedouble epsilon = 1.0e-6;
  su2double Area;
  string Local_TagBound, Global_TagBound;

  vector<su2double> Normal(nDim), UnitNormal(nDim), RefUnitNormal(nDim);
  /*--- Assume now that this boundary marker is straight. As soon as one
        AreaElement is found that is not aligned with a Reference then it
        is certain that the boundary marker is not straight and one can
        stop searching. Another possibility is that this process doesn't own
        any nodes of that boundary, in that case we also have to assume the
        boundary is straight. Any boundary type other than SYMMETRY_PLANE or
        EULER_WALL gets the value false (or see cases specified in the
        conditional below) which could be wrong. ---*/
  boundIsStraight.resize(nMarker);
  fill(boundIsStraight.begin(), boundIsStraight.end(), true);

  /*--- Loop over all local markers ---*/
  for (iMarker = 0; iMarker < nMarker; iMarker++) {
    Local_TagBound = config->GetMarker_All_TagBound(iMarker);

    /*--- Marker has to be Symmetry or Euler. Additionally marker can't be a
          moving surface and Grid Movement Elasticity is forbidden as well. All
          other GridMovements are rigid. ---*/
    if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE ||
         config->GetMarker_All_KindBC(iMarker) == EULER_WALL) &&
        !config->GetMarker_Moving_Bool(Local_TagBound) && !config->GetMarker_Deform_Mesh_Bool(Local_TagBound)) {
      /*--- Loop over all global markers, and find the local-global pair via
            matching unique string tags. ---*/
      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
        Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);
        if (Local_TagBound == Global_TagBound) {
          RefUnitNormal_defined = false;
          iVertex = 0;

          while (boundIsStraight[iMarker] && iVertex < nVertex[iMarker]) {
            vertex[iMarker][iVertex]->GetNormal(Normal.data());
            UnitNormal = Normal;

            /*--- Compute unit normal. ---*/
            Area = 0.0;
            for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim] * Normal[iDim];
            Area = sqrt(Area);

            /*--- Negate for outward convention. ---*/
            for (iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] /= -Area;

            /*--- Check if unit normal is within tolerance of the Reference unit normal.
                  Reference unit normal = first unit normal found. ---*/
            if (RefUnitNormal_defined) {
              for (iDim = 0; iDim < nDim; iDim++) {
                if (abs(RefUnitNormal[iDim] - UnitNormal[iDim]) > epsilon) {
                  boundIsStraight[iMarker] = false;
                  break;
                }
              }
            } else {
              RefUnitNormal = UnitNormal;  // deep copy of values
              RefUnitNormal_defined = true;
            }

            iVertex++;
          }  // while iVertex
        }    // if Local == Global
      }      // for iMarker_Global
    } else {
      /*--- Enforce default value: false ---*/
      boundIsStraight[iMarker] = false;
    }  // if sym or euler ...
  }    // for iMarker

  /*--- Communicate results and print on screen. ---*/
  if (print_on_screen) {
    /*--- Additional vector which can later be MPI::Allreduce(d) to pring the results
          on screen as nMarker (local) can vary across ranks. Default 'true' as it can
          happen that a local rank does not contain an element of each surface marker.  ---*/
    vector<bool> boundIsStraightGlobal(nMarker_Global, true);
    /*--- Match local with global tag bound and fill a Global Marker vector. ---*/
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      Local_TagBound = config->GetMarker_All_TagBound(iMarker);
      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
        Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);

        if (Local_TagBound == Global_TagBound) boundIsStraightGlobal[iMarker_Global] = boundIsStraight[iMarker];

      }  // for iMarker_Global
    }    // for iMarker

    vector<int> Buff_Send_isStraight(nMarker_Global), Buff_Recv_isStraight(nMarker_Global);

    /*--- Cast to int as std::vector<boolean> can be a special construct. MPI handling using <int>
          is more straight-forward. ---*/
    for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++)
      Buff_Send_isStraight[iMarker_Global] = static_cast<int>(boundIsStraightGlobal[iMarker_Global]);

    /*--- Product of type <int>(bool) is equivalnt to a 'logical and' ---*/
    SU2_MPI::Allreduce(Buff_Send_isStraight.data(), Buff_Recv_isStraight.data(), nMarker_Global, MPI_INT, MPI_PROD,
                       SU2_MPI::GetComm());

    /*--- Print results on screen. ---*/
    if (rank == MASTER_NODE) {
      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
        if (config->GetMarker_CfgFile_KindBC(config->GetMarker_CfgFile_TagBound(iMarker_Global)) == SYMMETRY_PLANE ||
            config->GetMarker_CfgFile_KindBC(config->GetMarker_CfgFile_TagBound(iMarker_Global)) == EULER_WALL) {
          cout << "Boundary marker " << config->GetMarker_CfgFile_TagBound(iMarker_Global) << " is";
          if (!static_cast<bool>(Buff_Recv_isStraight[iMarker_Global])) cout << " NOT";
          if (nDim == 2) cout << " a single straight." << endl;
          if (nDim == 3) cout << " a single plane." << endl;
        }  // if sym or euler
      }    // for iMarker_Global
    }      // if rank==MASTER
  }        // if print_on_scren
}

void CGeometry::ComputeSurf_Curvature(CConfig* config) {
  unsigned short iMarker, iNeigh_Point, iDim, iNode, iNeighbor_Nodes, Neighbor_Node;
  unsigned long Neighbor_Point, iVertex, iPoint, jPoint, iElem_Bound, iEdge, nLocalVertex, MaxLocalVertex,
      *Buffer_Send_nVertex, *Buffer_Receive_nVertex, TotalnPointDomain;
  vector<unsigned long> Point_NeighborList, Elem_NeighborList, Point_Triangle, Point_Edge, Point_Critical;
  su2double U[3] = {0.0}, V[3] = {0.0}, W[3] = {0.0}, Length_U, Length_V, Length_W, CosValue, Angle_Value, *K,
            *Angle_Defect, *Area_Vertex, *Angle_Alpha, *Angle_Beta, **NormalMeanK, MeanK, GaussK, MaxPrinK, cot_alpha,
            cot_beta, delta, X1, X2, X3, Y1, Y2, Y3, radius, *Buffer_Send_Coord, *Buffer_Receive_Coord, *Coord, Dist,
            MinDist, MaxK, MinK, SigmaK;
  bool* Check_Edge;

  /*--- Allocate surface curvature ---*/
  K = new su2double[nPoint];
  for (iPoint = 0; iPoint < nPoint; iPoint++) K[iPoint] = 0.0;

  if (nDim == 2) {
    /*--- Loop over all the markers ---*/
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
        /*--- Loop through all marker vertices again, this time also
         finding the neighbors of each node.---*/
        for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
          iPoint = vertex[iMarker][iVertex]->GetNode();

          if (nodes->GetDomain(iPoint)) {
            /*--- Loop through neighbors. In 2-D, there should be 2 nodes on either
             side of this vertex that lie on the same surface. ---*/
            Point_Edge.clear();

            for (iNeigh_Point = 0; iNeigh_Point < nodes->GetnPoint(iPoint); iNeigh_Point++) {
              Neighbor_Point = nodes->GetPoint(iPoint, iNeigh_Point);

              /*--- Check if this neighbor lies on the surface. If so,
               add to the list of neighbors. ---*/
              if (nodes->GetPhysicalBoundary(Neighbor_Point)) {
                Point_Edge.push_back(Neighbor_Point);
              }
            }

            if (Point_Edge.size() == 2) {
              /*--- Compute the curvature using three points ---*/
              X1 = nodes->GetCoord(iPoint, 0);
              X2 = nodes->GetCoord(Point_Edge[0], 0);
              X3 = nodes->GetCoord(Point_Edge[1], 0);
              Y1 = nodes->GetCoord(iPoint, 1);
              Y2 = nodes->GetCoord(Point_Edge[0], 1);
              Y3 = nodes->GetCoord(Point_Edge[1], 1);

              radius = sqrt(((X2 - X1) * (X2 - X1) + (Y2 - Y1) * (Y2 - Y1)) *
                            ((X2 - X3) * (X2 - X3) + (Y2 - Y3) * (Y2 - Y3)) *
                            ((X3 - X1) * (X3 - X1) + (Y3 - Y1) * (Y3 - Y1))) /
                       (2.0 * fabs(X1 * Y2 + X2 * Y3 + X3 * Y1 - X1 * Y3 - X2 * Y1 - X3 * Y2) + EPS);

              K[iPoint] = 1.0 / radius;
              nodes->SetCurvature(iPoint, K[iPoint]);
            }
          }
        }
      }
    }

  }

  else {
    Angle_Defect = new su2double[nPoint];
    Area_Vertex = new su2double[nPoint];
    for (iPoint = 0; iPoint < nPoint; iPoint++) {
      Angle_Defect[iPoint] = 2 * PI_NUMBER;
      Area_Vertex[iPoint] = 0.0;
    }

    Angle_Alpha = new su2double[nEdge];
    Angle_Beta = new su2double[nEdge];
    Check_Edge = new bool[nEdge];
    for (iEdge = 0; iEdge < nEdge; iEdge++) {
      Angle_Alpha[iEdge] = 0.0;
      Angle_Beta[iEdge] = 0.0;
      Check_Edge[iEdge] = true;
    }

    NormalMeanK = new su2double*[nPoint];
    for (iPoint = 0; iPoint < nPoint; iPoint++) {
      NormalMeanK[iPoint] = new su2double[nDim];
      for (iDim = 0; iDim < nDim; iDim++) {
        NormalMeanK[iPoint][iDim] = 0.0;
      }
    }

    /*--- Loop over all the markers ---*/
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
        /*--- Loop over all the boundary elements ---*/
        for (iElem_Bound = 0; iElem_Bound < nElem_Bound[iMarker]; iElem_Bound++) {
          /*--- Only triangles ---*/
          if (bound[iMarker][iElem_Bound]->GetVTK_Type() == TRIANGLE) {
            /*--- Loop over all the nodes of the boundary element ---*/
            for (iNode = 0; iNode < bound[iMarker][iElem_Bound]->GetnNodes(); iNode++) {
              iPoint = bound[iMarker][iElem_Bound]->GetNode(iNode);

              Point_Triangle.clear();

              for (iNeighbor_Nodes = 0; iNeighbor_Nodes < bound[iMarker][iElem_Bound]->GetnNeighbor_Nodes(iNode);
                   iNeighbor_Nodes++) {
                Neighbor_Node = bound[iMarker][iElem_Bound]->GetNeighbor_Nodes(iNode, iNeighbor_Nodes);
                Neighbor_Point = bound[iMarker][iElem_Bound]->GetNode(Neighbor_Node);
                Point_Triangle.push_back(Neighbor_Point);
              }

              iEdge = FindEdge(Point_Triangle[0], Point_Triangle[1]);

              for (iDim = 0; iDim < nDim; iDim++) {
                U[iDim] = nodes->GetCoord(Point_Triangle[0], iDim) - nodes->GetCoord(iPoint, iDim);
                V[iDim] = nodes->GetCoord(Point_Triangle[1], iDim) - nodes->GetCoord(iPoint, iDim);
              }

              W[0] = 0.5 * (U[1] * V[2] - U[2] * V[1]);
              W[1] = -0.5 * (U[0] * V[2] - U[2] * V[0]);
              W[2] = 0.5 * (U[0] * V[1] - U[1] * V[0]);

              Length_U = 0.0;
              Length_V = 0.0;
              Length_W = 0.0;
              CosValue = 0.0;
              for (iDim = 0; iDim < nDim; iDim++) {
                Length_U += U[iDim] * U[iDim];
                Length_V += V[iDim] * V[iDim];
                Length_W += W[iDim] * W[iDim];
              }
              Length_U = sqrt(Length_U);
              Length_V = sqrt(Length_V);
              Length_W = sqrt(Length_W);
              for (iDim = 0; iDim < nDim; iDim++) {
                U[iDim] /= Length_U;
                V[iDim] /= Length_V;
                CosValue += U[iDim] * V[iDim];
              }
              if (CosValue >= 1.0) CosValue = 1.0;
              if (CosValue <= -1.0) CosValue = -1.0;

              Angle_Value = acos(CosValue);
              Area_Vertex[iPoint] += Length_W;
              Angle_Defect[iPoint] -= Angle_Value;
              if (Angle_Alpha[iEdge] == 0.0)
                Angle_Alpha[iEdge] = Angle_Value;
              else
                Angle_Beta[iEdge] = Angle_Value;
            }
          }
        }
      }
    }

    /*--- Compute mean curvature ---*/
    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
        for (iElem_Bound = 0; iElem_Bound < nElem_Bound[iMarker]; iElem_Bound++) {
          if (bound[iMarker][iElem_Bound]->GetVTK_Type() == TRIANGLE) {
            for (iNode = 0; iNode < bound[iMarker][iElem_Bound]->GetnNodes(); iNode++) {
              iPoint = bound[iMarker][iElem_Bound]->GetNode(iNode);

              for (iNeighbor_Nodes = 0; iNeighbor_Nodes < bound[iMarker][iElem_Bound]->GetnNeighbor_Nodes(iNode);
                   iNeighbor_Nodes++) {
                Neighbor_Node = bound[iMarker][iElem_Bound]->GetNeighbor_Nodes(iNode, iNeighbor_Nodes);
                jPoint = bound[iMarker][iElem_Bound]->GetNode(Neighbor_Node);

                iEdge = FindEdge(iPoint, jPoint);

                if (Check_Edge[iEdge]) {
                  Check_Edge[iEdge] = false;

                  if (tan(Angle_Alpha[iEdge]) != 0.0)
                    cot_alpha = 1.0 / tan(Angle_Alpha[iEdge]);
                  else
                    cot_alpha = 0.0;
                  if (tan(Angle_Beta[iEdge]) != 0.0)
                    cot_beta = 1.0 / tan(Angle_Beta[iEdge]);
                  else
                    cot_beta = 0.0;

                  /*--- iPoint, and jPoint ---*/
                  for (iDim = 0; iDim < nDim; iDim++) {
                    if (Area_Vertex[iPoint] != 0.0)
                      NormalMeanK[iPoint][iDim] += 3.0 * (cot_alpha + cot_beta) *
                                                   (nodes->GetCoord(iPoint, iDim) - nodes->GetCoord(jPoint, iDim)) /
                                                   Area_Vertex[iPoint];
                    if (Area_Vertex[jPoint] != 0.0)
                      NormalMeanK[jPoint][iDim] += 3.0 * (cot_alpha + cot_beta) *
                                                   (nodes->GetCoord(jPoint, iDim) - nodes->GetCoord(iPoint, iDim)) /
                                                   Area_Vertex[jPoint];
                  }
                }
              }
            }
          }
        }
      }
    }

    /*--- Compute Gauss, mean, max and min principal curvature,
     and set the list of critical points ---*/

    for (iMarker = 0; iMarker < nMarker; iMarker++) {
      if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
        for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
          iPoint = vertex[iMarker][iVertex]->GetNode();

          if (nodes->GetDomain(iPoint)) {
            if (Area_Vertex[iPoint] != 0.0)
              GaussK = 3.0 * Angle_Defect[iPoint] / Area_Vertex[iPoint];
            else
              GaussK = 0.0;

            MeanK = 0.0;
            for (iDim = 0; iDim < nDim; iDim++) MeanK += NormalMeanK[iPoint][iDim] * NormalMeanK[iPoint][iDim];
            MeanK = sqrt(MeanK);

            delta = max((MeanK * MeanK - GaussK), 0.0);

            MaxPrinK = MeanK + sqrt(delta);

            /*--- Store the curvature value ---*/
            K[iPoint] = MaxPrinK;
            nodes->SetCurvature(iPoint, K[iPoint]);
          }
        }
      }
    }

    delete[] Angle_Defect;
    delete[] Area_Vertex;
    delete[] Angle_Alpha;
    delete[] Angle_Beta;
    delete[] Check_Edge;

    for (iPoint = 0; iPoint < nPoint; iPoint++) delete[] NormalMeanK[iPoint];
    delete[] NormalMeanK;
  }

  /*--- Sharp edge detection is based in the statistical
   distribution of the curvature ---*/

  MaxK = K[0];
  MinK = K[0];
  MeanK = 0.0;
  TotalnPointDomain = 0;
  for (iMarker = 0; iMarker < nMarker; iMarker++) {
    if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
      for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
        iPoint = vertex[iMarker][iVertex]->GetNode();
        if (nodes->GetDomain(iPoint)) {
          MaxK = max(MaxK, fabs(K[iPoint]));
          MinK = min(MinK, fabs(K[iPoint]));
          MeanK += fabs(K[iPoint]);
          TotalnPointDomain++;
        }
      }
    }
  }

  su2double MyMeanK = MeanK;
  MeanK = 0.0;
  su2double MyMaxK = MaxK;
  MaxK = 0.0;
  unsigned long MynPointDomain = TotalnPointDomain;
  TotalnPointDomain = 0;
  SU2_MPI::Allreduce(&MyMeanK, &MeanK, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
  SU2_MPI::Allreduce(&MyMaxK, &MaxK, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());
  SU2_MPI::Allreduce(&MynPointDomain, &TotalnPointDomain, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());

  /*--- Compute the mean ---*/
  MeanK /= su2double(TotalnPointDomain);

  /*--- Compute the standard deviation ---*/
  SigmaK = 0.0;
  for (iMarker = 0; iMarker < nMarker; iMarker++) {
    if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
      for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
        iPoint = vertex[iMarker][iVertex]->GetNode();
        if (nodes->GetDomain(iPoint)) {
          SigmaK += (fabs(K[iPoint]) - MeanK) * (fabs(K[iPoint]) - MeanK);
        }
      }
    }
  }

  su2double MySigmaK = SigmaK;
  SigmaK = 0.0;
  SU2_MPI::Allreduce(&MySigmaK, &SigmaK, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());

  SigmaK = sqrt(SigmaK / su2double(TotalnPointDomain));

  if (rank == MASTER_NODE)
    cout << "Max K: " << MaxK << ". Mean K: " << MeanK << ". Standard deviation K: " << SigmaK << "." << endl;

  Point_Critical.clear();

  for (iMarker = 0; iMarker < nMarker; iMarker++) {
    if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
      for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
        iPoint = vertex[iMarker][iVertex]->GetNode();
        if (nodes->GetDomain(iPoint)) {
          if (fabs(K[iPoint]) > MeanK + config->GetRefSharpEdges() * SigmaK) {
            Point_Critical.push_back(iPoint);
          }
        }
      }
    }
  }

  /*--- Variables and buffers needed for MPI ---*/

  Buffer_Send_nVertex = new unsigned long[1];
  Buffer_Receive_nVertex = new unsigned long[size];

  /*--- Count the total number of critical edge nodes. ---*/

  nLocalVertex = Point_Critical.size();
  Buffer_Send_nVertex[0] = nLocalVertex;

  /*--- Communicate to all processors the total number of critical edge nodes. ---*/

  MaxLocalVertex = 0;
  SU2_MPI::Allreduce(&nLocalVertex, &MaxLocalVertex, 1, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm());
  SU2_MPI::Allgather(Buffer_Send_nVertex, 1, MPI_UNSIGNED_LONG, Buffer_Receive_nVertex, 1, MPI_UNSIGNED_LONG,
                     SU2_MPI::GetComm());

  /*--- Create and initialize to zero some buffers to hold the coordinates
   of the boundary nodes that are communicated from each partition (all-to-all). ---*/

  const unsigned long nBuffer = MaxLocalVertex * nDim;
  Buffer_Send_Coord = new su2double[nBuffer]();
  Buffer_Receive_Coord = new su2double[size * nBuffer];

  /*--- Retrieve and store the coordinates of the sharp edges boundary nodes on
   the local partition and broadcast them to all partitions. ---*/

  for (iVertex = 0; iVertex < Point_Critical.size(); iVertex++) {
    iPoint = Point_Critical[iVertex];
    for (iDim = 0; iDim < nDim; iDim++) Buffer_Send_Coord[iVertex * nDim + iDim] = nodes->GetCoord(iPoint, iDim);
  }

  SU2_MPI::Allgather(Buffer_Send_Coord, nBuffer, MPI_DOUBLE, Buffer_Receive_Coord, nBuffer, MPI_DOUBLE,
                     SU2_MPI::GetComm());

  /*--- Loop over all interior mesh nodes on the local partition and compute
   the distances to each of the no-slip boundary nodes in the entire mesh.
   Store the minimum distance to the wall for each interior mesh node. ---*/

  for (iPoint = 0; iPoint < GetnPoint(); iPoint++) {
    Coord = nodes->GetCoord(iPoint);

    MinDist = 1E20;
    for (int iProcessor = 0; iProcessor < size; iProcessor++) {
      for (iVertex = 0; iVertex < Buffer_Receive_nVertex[iProcessor]; iVertex++) {
        Dist = 0.0;
        for (iDim = 0; iDim < nDim; iDim++) {
          Dist += (Coord[iDim] - Buffer_Receive_Coord[(iProcessor * MaxLocalVertex + iVertex) * nDim + iDim]) *
                  (Coord[iDim] - Buffer_Receive_Coord[(iProcessor * MaxLocalVertex + iVertex) * nDim + iDim]);
        }
        if (Dist != 0.0)
          Dist = sqrt(Dist);
        else
          Dist = 0.0;
        if (Dist < MinDist) MinDist = Dist;
      }
    }
    nodes->SetSharpEdge_Distance(iPoint, MinDist);
  }

  /*--- Deallocate Max curvature ---*/
  delete[] K;

  /*--- Deallocate the buffers needed for the MPI communication. ---*/
  delete[] Buffer_Send_Coord;
  delete[] Buffer_Receive_Coord;
  delete[] Buffer_Send_nVertex;
  delete[] Buffer_Receive_nVertex;
}

void CGeometry::FilterValuesAtElementCG(const vector<su2double>& filter_radius,
                                        const vector<pair<ENUM_FILTER_KERNEL, su2double>>& kernels,
                                        const unsigned short search_limit, su2double* values) const {
  /*--- Apply a filter to "input_values". The filter is an averaging process over the neighbourhood
  of each element, which is a circle in 2D and a sphere in 3D of radius "filter_radius".
  The filter is characterized by its kernel, i.e. how the weights are computed. Multiple kernels
  can be specified in which case they are applied sequentially (each one being applied to the
  output values of the previous filter. ---*/

  /*--- Check if we need to do any work. ---*/
  if (kernels.empty()) return;

  /*--- FIRST: Gather the adjacency matrix, element centroids, volumes, and values on every
  processor, this is required because the filter reaches far into adjacent partitions. ---*/

  /*--- Adjacency matrix ---*/
  vector<unsigned long> neighbour_start;
  long* neighbour_idx = nullptr;
  GetGlobalElementAdjacencyMatrix(neighbour_start, neighbour_idx);

  /*--- Element centroids and volumes. ---*/
  auto *cg_elem = new su2double[Global_nElemDomain * nDim], *vol_elem = new su2double[Global_nElemDomain];
#ifdef HAVE_MPI
  /*--- Number of subdomain each point is part of. ---*/
  vector<char> halo_detect(Global_nElemDomain);
#endif

  /*--- Inputs of a filter stage, like with CG and volumes, each processor needs to see everything. ---*/
  auto* work_values = new su2double[Global_nElemDomain];

  /*--- When gathering the neighborhood of each element we use a vector of booleans to indicate
  whether an element is already added to the list of neighbors (one vector per thread). ---*/
  vector<vector<bool>> is_neighbor(omp_get_max_threads());

  /*--- Begin OpenMP parallel section, count total number of searches for which
  the recursion limit is reached and the full neighborhood is not considered. ---*/
  unsigned long limited_searches = 0;

  SU2_OMP_PARALLEL_(reduction(+ : limited_searches)) {
    /*--- Initialize ---*/
    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < Global_nElemDomain; ++iElem) {
      for (unsigned short iDim = 0; iDim < nDim; ++iDim) cg_elem[nDim * iElem + iDim] = 0.0;
      vol_elem[iElem] = 0.0;
    }
    END_SU2_OMP_FOR

    /*--- Populate ---*/
    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < nElem; ++iElem) {
      auto iElem_global = elem[iElem]->GetGlobalIndex();
      for (unsigned short iDim = 0; iDim < nDim; ++iDim) cg_elem[nDim * iElem_global + iDim] = elem[iElem]->GetCG(iDim);
      vol_elem[iElem_global] = elem[iElem]->GetVolume();
    }
    END_SU2_OMP_FOR

#ifdef HAVE_MPI
    /*--- Account for the duplication introduced by the halo elements and the
    reduction using MPI_SUM, which is required to maintain differentiabillity. ---*/
    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < Global_nElemDomain; ++iElem) halo_detect[iElem] = 0;
    END_SU2_OMP_FOR

    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < nElem; ++iElem) halo_detect[elem[iElem]->GetGlobalIndex()] = 1;
    END_SU2_OMP_FOR

    /*--- Share with all processors ---*/
    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
      auto* dbl_buffer = new su2double[Global_nElemDomain * nDim];
      SU2_MPI::Allreduce(cg_elem, dbl_buffer, Global_nElemDomain * nDim, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
      swap(dbl_buffer, cg_elem);
      delete[] dbl_buffer;

      dbl_buffer = new su2double[Global_nElemDomain];
      SU2_MPI::Allreduce(vol_elem, dbl_buffer, Global_nElemDomain, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
      swap(dbl_buffer, vol_elem);
      delete[] dbl_buffer;

      vector<char> char_buffer(Global_nElemDomain);
      MPI_Allreduce(halo_detect.data(), char_buffer.data(), Global_nElemDomain, MPI_CHAR, MPI_SUM, SU2_MPI::GetComm());
      halo_detect.swap(char_buffer);
    }
    END_SU2_OMP_SAFE_GLOBAL_ACCESS

    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < Global_nElemDomain; ++iElem) {
      su2double numRepeat = halo_detect[iElem];
      for (unsigned short iDim = 0; iDim < nDim; ++iDim) cg_elem[nDim * iElem + iDim] /= numRepeat;
      vol_elem[iElem] /= numRepeat;
    }
    END_SU2_OMP_FOR
#endif

    /*--- SECOND: Each processor performs the average for its elements. For each
    element we look for neighbours of neighbours of... until the distance to the
    closest newly found one is greater than the filter radius.  ---*/

    is_neighbor[omp_get_thread_num()].resize(Global_nElemDomain, false);

    for (unsigned long iKernel = 0; iKernel < kernels.size(); ++iKernel) {
      auto kernel_type = kernels[iKernel].first;
      su2double kernel_param = kernels[iKernel].second;
      su2double kernel_radius = filter_radius[iKernel];

      /*--- Synchronize work values ---*/
      /*--- Initialize ---*/
      SU2_OMP_FOR_STAT(256)
      for (auto iElem = 0ul; iElem < Global_nElemDomain; ++iElem) work_values[iElem] = 0.0;
      END_SU2_OMP_FOR

      /*--- Populate ---*/
      SU2_OMP_FOR_STAT(256)
      for (auto iElem = 0ul; iElem < nElem; ++iElem) work_values[elem[iElem]->GetGlobalIndex()] = values[iElem];
      END_SU2_OMP_FOR

#ifdef HAVE_MPI
      /*--- Share with all processors ---*/
      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
        auto* buffer = new su2double[Global_nElemDomain];
        SU2_MPI::Allreduce(work_values, buffer, Global_nElemDomain, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
        swap(buffer, work_values);
        delete[] buffer;
      }
      END_SU2_OMP_SAFE_GLOBAL_ACCESS

      /*--- Account for duplication ---*/
      SU2_OMP_FOR_STAT(256)
      for (auto iElem = 0ul; iElem < Global_nElemDomain; ++iElem) {
        su2double numRepeat = halo_detect[iElem];
        work_values[iElem] /= numRepeat;
      }
      END_SU2_OMP_FOR
#endif

      /*--- Filter ---*/
      SU2_OMP_FOR_DYN(128)
      for (auto iElem = 0ul; iElem < nElem; ++iElem) {
        int thread = omp_get_thread_num();

        /*--- Center of the search ---*/
        auto iElem_global = elem[iElem]->GetGlobalIndex();

        /*--- Find the neighbours of iElem ---*/
        vector<long> neighbours;
        limited_searches +=
            !GetRadialNeighbourhood(iElem_global, SU2_TYPE::GetValue(kernel_radius), search_limit, neighbour_start,
                                    neighbour_idx, cg_elem, neighbours, is_neighbor[thread]);
        /*--- Apply the kernel ---*/
        su2double weight = 0.0, numerator = 0.0, denominator = 0.0;

        switch (kernel_type) {
          /*--- distance-based kernels (weighted averages) ---*/
          case ENUM_FILTER_KERNEL::CONSTANT_WEIGHT:
          case ENUM_FILTER_KERNEL::CONICAL_WEIGHT:
          case ENUM_FILTER_KERNEL::GAUSSIAN_WEIGHT:

            for (auto idx : neighbours) {
              su2double distance = 0.0;
              for (unsigned short iDim = 0; iDim < nDim; ++iDim)
                distance += pow(cg_elem[nDim * iElem_global + iDim] - cg_elem[nDim * idx + iDim], 2);
              distance = sqrt(distance);

              switch (kernel_type) {
                case ENUM_FILTER_KERNEL::CONSTANT_WEIGHT:
                  weight = 1.0;
                  break;
                case ENUM_FILTER_KERNEL::CONICAL_WEIGHT:
                  weight = kernel_radius - distance;
                  break;
                case ENUM_FILTER_KERNEL::GAUSSIAN_WEIGHT:
                  weight = exp(-0.5 * pow(distance / kernel_param, 2));
                  break;
                default:
                  break;
              }
              weight *= vol_elem[idx];
              numerator += weight * work_values[idx];
              denominator += weight;
            }
            values[iElem] = numerator / denominator;
            break;

          /*--- morphology kernels (image processing) ---*/
          case ENUM_FILTER_KERNEL::DILATE_MORPH:
          case ENUM_FILTER_KERNEL::ERODE_MORPH:

            for (auto idx : neighbours) {
              switch (kernel_type) {
                case ENUM_FILTER_KERNEL::DILATE_MORPH:
                  numerator += exp(kernel_param * work_values[idx]);
                  break;
                case ENUM_FILTER_KERNEL::ERODE_MORPH:
                  numerator += exp(kernel_param * (1.0 - work_values[idx]));
                  break;
                default:
                  break;
              }
              denominator += 1.0;
            }
            values[iElem] = log(numerator / denominator) / kernel_param;
            if (kernel_type == ENUM_FILTER_KERNEL::ERODE_MORPH) values[iElem] = 1.0 - values[iElem];
            break;

          default:
            SU2_MPI::Error("Unknown type of filter kernel", CURRENT_FUNCTION);
        }
      }
      END_SU2_OMP_FOR
    }
  }
  END_SU2_OMP_PARALLEL

  limited_searches /= kernels.size();

  unsigned long tmp = limited_searches;
  SU2_MPI::Reduce(&tmp, &limited_searches, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm());

  if (rank == MASTER_NODE && limited_searches > 0)
    cout << "Warning: The filter radius was limited for " << limited_searches << " elements ("
         << limited_searches / (0.01 * Global_nElemDomain) << "%).\n";

  delete[] neighbour_idx;
  delete[] cg_elem;
  delete[] vol_elem;
  delete[] work_values;
}

void CGeometry::GetGlobalElementAdjacencyMatrix(vector<unsigned long>& neighbour_start, long*& neighbour_idx) const {
  if (neighbour_idx != nullptr)
    SU2_MPI::Error("neighbour_idx is expected to be NULL, stopping to avoid a potential memory leak", CURRENT_FUNCTION);

  /*--- Determine how much space we need for the adjacency matrix by counting the
  neighbours of each element, i.e. its number of faces---*/
  auto* nFaces_elem = new unsigned short[Global_nElemDomain];

  SU2_OMP_PARALLEL {
    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < Global_nElemDomain; ++iElem) nFaces_elem[iElem] = 0;
    END_SU2_OMP_FOR

    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < nElem; ++iElem) {
      auto iElem_global = elem[iElem]->GetGlobalIndex();
      nFaces_elem[iElem_global] = elem[iElem]->GetnFaces();
    }
    END_SU2_OMP_FOR
  }
  END_SU2_OMP_PARALLEL
#ifdef HAVE_MPI
  /*--- Share with all processors ---*/
  {
    auto* buffer = new unsigned short[Global_nElemDomain];
    MPI_Allreduce(nFaces_elem, buffer, Global_nElemDomain, MPI_UNSIGNED_SHORT, MPI_MAX, SU2_MPI::GetComm());
    /*--- swap pointers and delete old data to keep the same variable name after reduction ---*/
    swap(buffer, nFaces_elem);
    delete[] buffer;
  }
#endif

  /*--- Vector with the addresses of the start of the neighbours of a given element.
  This is generated by a cumulative sum of the neighbour count. ---*/
  neighbour_start.resize(Global_nElemDomain + 1);

  neighbour_start[0] = 0;
  for (auto iElem = 0ul; iElem < Global_nElemDomain; ++iElem) {
    neighbour_start[iElem + 1] = neighbour_start[iElem] + nFaces_elem[iElem];
  }
  delete[] nFaces_elem;

  /*--- Allocate ---*/
  unsigned long matrix_size = neighbour_start[Global_nElemDomain];
  neighbour_idx = new long[matrix_size];

  SU2_OMP_PARALLEL {
    /*--- Initialize ---*/
    SU2_OMP_FOR_STAT(256)
    for (auto iElem = 0ul; iElem < matrix_size; ++iElem) neighbour_idx[iElem] = -1;
    END_SU2_OMP_FOR

    /*--- Populate ---*/
    SU2_OMP_FOR_STAT(128)
    for (auto iElem = 0ul; iElem < nElem; ++iElem) {
      auto iElem_global = elem[iElem]->GetGlobalIndex();
      auto start_pos = neighbour_start[iElem_global];

      for (unsigned short iFace = 0; iFace < elem[iElem]->GetnFaces(); ++iFace) {
        long neighbour = elem[iElem]->GetNeighbor_Elements(iFace);

        if (neighbour >= 0) {
          neighbour_idx[start_pos + iFace] = elem[neighbour]->GetGlobalIndex();
        }
      }
    }
    END_SU2_OMP_FOR
  }
  END_SU2_OMP_PARALLEL
#ifdef HAVE_MPI
  /*--- Share with all processors ---*/
  {
    long* buffer = new long[matrix_size];
    MPI_Allreduce(neighbour_idx, buffer, matrix_size, MPI_LONG, MPI_MAX, SU2_MPI::GetComm());
    swap(buffer, neighbour_idx);
    delete[] buffer;
  }
#endif
}

bool CGeometry::GetRadialNeighbourhood(const unsigned long iElem_global, const passivedouble radius,
                                       size_t search_limit, const vector<unsigned long>& neighbour_start,
                                       const long* neighbour_idx, const su2double* cg_elem, vector<long>& neighbours,
                                       vector<bool>& is_neighbor) const {
  /*--- Validate inputs if we are debugging. ---*/
  assert(neighbour_start.size() == Global_nElemDomain + 1 && neighbour_idx != nullptr && cg_elem != nullptr &&
         is_neighbor.size() == Global_nElemDomain && "invalid inputs");

  /*--- 0 search_limit means "unlimited" (it will probably
   stop once it gathers the entire domain, probably). ---*/
  if (!search_limit) search_limit = numeric_limits<size_t>::max();

  /*--- Center of the search ---*/
  neighbours.clear();
  neighbours.push_back(iElem_global);
  is_neighbor[iElem_global] = true;

  passivedouble X0[3] = {0.0, 0.0, 0.0};
  for (unsigned short iDim = 0; iDim < nDim; ++iDim) X0[iDim] = SU2_TYPE::GetValue(cg_elem[nDim * iElem_global + iDim]);

  /*--- Loop stops when "neighbours" stops changing size, or degree reaches limit. ---*/
  bool finished = false;
  for (size_t degree = 0, start = 0; degree < search_limit && !finished; ++degree) {
    /*--- For each element of the last degree added consider its immediate
     neighbours, that are not already neighbours, as candidates. ---*/
    vector<long> candidates;

    for (auto it = neighbours.begin() + start; it != neighbours.end(); ++it) {
      /*--- scan row of the adjacency matrix of element *it ---*/
      for (auto i = neighbour_start[*it]; i < neighbour_start[(*it) + 1]; ++i) {
        auto idx = neighbour_idx[i];
        if (idx >= 0)
          if (!is_neighbor[idx]) {
            candidates.push_back(idx);
            /*--- mark as neighbour for now to avoid duplicate candidates. ---*/
            is_neighbor[idx] = true;
          }
      }
    }
    /*--- update start position to fetch next degree candidates. ---*/
    start = neighbours.size();

    /*--- Add candidates within "radius" of X0, if none qualifies we are "finished". ---*/
    finished = true;
    for (auto idx : candidates) {
      /*--- passivedouble as we only need to compare "distance". ---*/
      passivedouble distance = 0.0;
      for (unsigned short iDim = 0; iDim < nDim; ++iDim)
        distance += pow(X0[iDim] - SU2_TYPE::GetValue(cg_elem[nDim * idx + iDim]), 2);

      if (distance < pow(radius, 2)) {
        neighbours.push_back(idx);
        finished = false;
      }
      /*--- not a neighbour in the end. ---*/
      else
        is_neighbor[idx] = false;
    }
  }
  /*--- Restore the state of the working vector for next call. ---*/
  for (auto idx : neighbours) is_neighbor[idx] = false;

  return finished;
}

void CGeometry::SetElemVolume() {
  SU2_OMP_PARALLEL {
    /*--- Create a bank of elements to avoid instantiating inside loop. ---*/
    CElement* elements[4] = {nullptr, nullptr, nullptr, nullptr};

    if (nDim == 2) {
      elements[0] = new CTRIA1();
      elements[1] = new CQUAD4();
    } else {
      elements[0] = new CTETRA1();
      elements[1] = new CPYRAM5();
      elements[2] = new CPRISM6();
      elements[3] = new CHEXA8();
    }

    /*--- Compute and store the volume of each "elem". ---*/
    SU2_OMP_FOR_DYN(128)
    for (unsigned long iElem = 0; iElem < nElem; ++iElem) {
      /*--- Get the appropriate type of element. ---*/
      CElement* element = nullptr;
      switch (elem[iElem]->GetVTK_Type()) {
        case TRIANGLE:
          element = elements[0];
          break;
        case QUADRILATERAL:
          element = elements[1];
          break;
        case TETRAHEDRON:
          element = elements[0];
          break;
        case PYRAMID:
          element = elements[1];
          break;
        case PRISM:
          element = elements[2];
          break;
        case HEXAHEDRON:
          element = elements[3];
          break;
        default:
          SU2_MPI::Error("Cannot compute the area/volume of a 1D element.", CURRENT_FUNCTION);
      }
      /*--- Set the nodal coordinates of the element. ---*/
      for (unsigned short iNode = 0; iNode < elem[iElem]->GetnNodes(); ++iNode) {
        unsigned long node_idx = elem[iElem]->GetNode(iNode);
        for (unsigned short iDim = 0; iDim < nDim; ++iDim) {
          su2double coord = nodes->GetCoord(node_idx, iDim);
          element->SetRef_Coord(iNode, iDim, coord);
        }
      }
      /*--- Compute ---*/
      if (nDim == 2)
        elem[iElem]->SetVolume(element->ComputeArea());
      else
        elem[iElem]->SetVolume(element->ComputeVolume());
    }
    END_SU2_OMP_FOR

    delete elements[0];
    delete elements[1];
    if (nDim == 3) {
      delete elements[2];
      delete elements[3];
    }
  }
  END_SU2_OMP_PARALLEL
}

void CGeometry::SetRotationalVelocity(const CConfig* config, bool print) {
  unsigned long iPoint;
  unsigned short iDim;

  su2double GridVel[3] = {0.0, 0.0, 0.0}, Distance[3] = {0.0, 0.0, 0.0}, Center[3] = {0.0, 0.0, 0.0},
            Omega[3] = {0.0, 0.0, 0.0}, xDot[3] = {0.0, 0.0, 0.0};

  /*--- Center of rotation & angular velocity vector from config ---*/

  for (iDim = 0; iDim < 3; iDim++) {
    Center[iDim] = config->GetMotion_Origin(iDim);
    Omega[iDim] = config->GetRotation_Rate(iDim) / config->GetOmega_Ref();
    xDot[iDim] = config->GetTranslation_Rate(iDim) / config->GetVelocity_Ref();
  }

  su2double L_Ref = config->GetLength_Ref();

  /*--- Print some information to the console ---*/

  if (rank == MASTER_NODE && print) {
    cout << " Rotational origin (x, y, z): ( " << Center[0] << ", " << Center[1];
    cout << ", " << Center[2] << " )\n";
    cout << " Angular velocity about x, y, z axes: ( " << Omega[0] << ", ";
    cout << Omega[1] << ", " << Omega[2] << " ) rad/s" << endl;
    cout << " Translational velocity in x, y, z direction: (" << xDot[0] << ", " << xDot[1] << ", " << xDot[2] << ")."
         << endl;
  }

  /*--- Loop over all nodes and set the rotational and translational velocity ---*/

  for (iPoint = 0; iPoint < nPoint; iPoint++) {
    /*--- Get the coordinates of the current node ---*/

    const su2double* Coord = nodes->GetCoord(iPoint);

    /*--- Calculate the non-dim. distance from the rotation center ---*/

    for (iDim = 0; iDim < nDim; iDim++) Distance[iDim] = (Coord[iDim] - Center[iDim]) / L_Ref;

    /*--- Calculate the angular velocity as omega X r and add translational velocity ---*/

    GridVel[0] = Omega[1] * (Distance[2]) - Omega[2] * (Distance[1]) + xDot[0];
    GridVel[1] = Omega[2] * (Distance[0]) - Omega[0] * (Distance[2]) + xDot[1];
    GridVel[2] = Omega[0] * (Distance[1]) - Omega[1] * (Distance[0]) + xDot[2];

    /*--- Store the grid velocity at this node ---*/

    nodes->SetGridVel(iPoint, GridVel);
  }
}

void CGeometry::SetShroudVelocity(const CConfig* config) {
  unsigned long iPoint, iVertex;
  unsigned short iMarker, iMarkerShroud;
  su2double RotVel[3] = {0.0, 0.0, 0.0};

  /*--- Loop over all vertex in the shroud marker and set the rotational velocity to 0.0 ---*/
  for (iMarker = 0; iMarker < nMarker; iMarker++) {
    for (iMarkerShroud = 0; iMarkerShroud < config->GetnMarker_Shroud(); iMarkerShroud++) {
      if (config->GetMarker_Shroud(iMarkerShroud) == config->GetMarker_All_TagBound(iMarker)) {
        for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
          iPoint = vertex[iMarker][iVertex]->GetNode();
          nodes->SetGridVel(iPoint, RotVel);
        }
      }
    }
  }
}

void CGeometry::SetTranslationalVelocity(const CConfig* config, bool print) {
  su2double xDot[3] = {0.0, 0.0, 0.0};

  /*--- Get the translational velocity vector from config ---*/

  for (unsigned short iDim = 0; iDim < nDim; iDim++)
    xDot[iDim] = config->GetTranslation_Rate(iDim) / config->GetVelocity_Ref();

  /*--- Print some information to the console ---*/

  if (rank == MASTER_NODE && print) {
    cout << " Non-dim. translational velocity: (" << xDot[0] << ", " << xDot[1] << ", " << xDot[2] << ")." << endl;
  }

  /*--- Loop over all nodes and set the translational velocity ---*/

  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) nodes->SetGridVel(iPoint, xDot);
}

void CGeometry::SetWallVelocity(const CConfig* config, bool print) {
  const su2double L_Ref = config->GetLength_Ref();
  const su2double Omega_Ref = config->GetOmega_Ref();
  const su2double Vel_Ref = config->GetVelocity_Ref();

  /*--- Store grid velocity for each node on the moving surface(s).
   Sum and store the x, y, & z velocities due to translation and rotation. ---*/

  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
    if (config->GetMarker_All_Moving(iMarker) != YES) continue;

    /*--- Identify iMarker from the list of those under MARKER_MOVING
     Get prescribed wall speed from config for this marker. ---*/

    const auto Marker_Tag = config->GetMarker_All_TagBound(iMarker);
    const auto jMarker = config->GetMarker_Moving(Marker_Tag);

    su2double xDot[MAXNDIM], Center[MAXNDIM], Omega[MAXNDIM];

    for (auto iDim = 0u; iDim < MAXNDIM; iDim++) {
      Center[iDim] = config->GetMarkerMotion_Origin(jMarker, iDim);
      Omega[iDim] = config->GetMarkerRotationRate(jMarker, iDim) / Omega_Ref;
      xDot[iDim] = config->GetMarkerTranslationRate(jMarker, iDim) / Vel_Ref;
    }

    if (rank == MASTER_NODE && print) {
      cout << " Storing grid velocity for marker: ";
      cout << Marker_Tag << ".\n";
      cout << " Translational velocity: (" << xDot[0] * config->GetVelocity_Ref() << ", "
           << xDot[1] * config->GetVelocity_Ref();
      cout << ", " << xDot[2] * config->GetVelocity_Ref();
      if (config->GetSystemMeasurements() == SI)
        cout << ") m/s.\n";
      else
        cout << ") ft/s.\n";
      cout << " Angular velocity: (" << Omega[0] << ", " << Omega[1];
      cout << ", " << Omega[2] << ") rad/s about origin: (" << Center[0];
      cout << ", " << Center[1] << ", " << Center[2] << ")." << endl;
    }

    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
      const auto iPoint = vertex[iMarker][iVertex]->GetNode();

      /*--- Calculate non-dim. position from rotation center ---*/
      su2double r[MAXNDIM] = {0.0};
      for (auto iDim = 0u; iDim < nDim; iDim++) r[iDim] = (nodes->GetCoord(iPoint, iDim) - Center[iDim]) / L_Ref;

      /*--- Cross Product of angular velocity and distance from center to
       get the rotational velocity. Note that we are adding on the velocity
       due to pure translation as well. ---*/

      su2double GridVel[MAXNDIM];
      GeometryToolbox::CrossProduct(Omega, r, GridVel);

      for (auto iDim = 0u; iDim < nDim; iDim++) nodes->SetGridVel(iPoint, iDim, xDot[iDim] + GridVel[iDim]);
    }
  }
}

void CGeometry::SetGridVelocity(const CConfig* config) {
  /*--- Get timestep and whether to use 1st or 2nd order backward finite differences ---*/

  bool FirstOrder = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST);
  bool SecondOrder = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);

  su2double TimeStep = config->GetDelta_UnstTimeND();

  /*--- Compute the velocity of each node in the volume mesh ---*/

  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
    /*--- Coordinates of the current point at n+1, n, & n-1 time levels ---*/

    const su2double* Coord_nM1 = nodes->GetCoord_n1(iPoint);
    const su2double* Coord_n = nodes->GetCoord_n(iPoint);
    const su2double* Coord_nP1 = nodes->GetCoord(iPoint);

    /*--- Compute and store mesh velocity with 1st or 2nd-order approximation ---*/

    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
      su2double GridVel = 0.0;

      if (FirstOrder) GridVel = (Coord_nP1[iDim] - Coord_n[iDim]) / TimeStep;

      if (SecondOrder) GridVel = (1.5 * Coord_nP1[iDim] - 2.0 * Coord_n[iDim] + 0.5 * Coord_nM1[iDim]) / TimeStep;

      nodes->SetGridVel(iPoint, iDim, GridVel);
    }
  }
}

const CGeometry::LDUSparsePattern& CGeometry::GetSparsePattern(ConnectivityType type, unsigned long fillLvl) {
  bool fvm = (type == ConnectivityType::FiniteVolume);
  auto& grp = fillLvl == 0 ? (fvm ? finiteVolumePatternFill0 : finiteElementPatternFill0)
                           : (fvm ? finiteVolumePatternFillN : finiteElementPatternFillN);
  if (grp.empty()) {
    grp.csr = buildCSRPattern(*this, type, static_cast<su2uint>(fillLvl));
    grp.csr.buildDiagPtr();
    grp.l = buildLowerPattern(grp.csr);
    grp.u = buildUpperPattern(grp.csr);
  }
  return grp;
}

const su2vector<su2uint>& CGeometry::GetLToUTransposeSparsePatternMap(ConnectivityType type) {
  bool fvm = (type == ConnectivityType::FiniteVolume);
  auto& l_to_u = fvm ? finiteVolumeLToUTranspMap : finiteElementLToUTranspMap;
  if (l_to_u.empty()) {
    auto& u_to_l = fvm ? finiteVolumeUToLTranspMap : finiteElementUToLTranspMap;
    const auto& pat = GetSparsePattern(type);
    buildLUTransposeMaps(pat.l, pat.u, l_to_u, u_to_l);
  }
  return l_to_u;
}

const su2vector<su2uint>& CGeometry::GetUToLTransposeSparsePatternMap(ConnectivityType type) {
  bool fvm = (type == ConnectivityType::FiniteVolume);
  auto& u_to_l = fvm ? finiteVolumeUToLTranspMap : finiteElementUToLTranspMap;
  if (u_to_l.empty()) {
    auto& l_to_u = fvm ? finiteVolumeLToUTranspMap : finiteElementLToUTranspMap;
    const auto& pat = GetSparsePattern(type);
    buildLUTransposeMaps(pat.l, pat.u, l_to_u, u_to_l);
  }
  return u_to_l;
}

const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(su2double* efficiency, bool maximizeEdgeColorGroupSize) {
  /*--- Check for dry run mode with dummy geometry. ---*/
  if (nEdge == 0) return edgeColoring;

  /*--- Build if required. ---*/
  if (edgeColoring.empty()) {
    /*--- When not using threading use the natural coloring to reduce overhead. ---*/
    if (omp_get_max_threads() == 1) {
      SetNaturalEdgeColoring();
      if (efficiency != nullptr) *efficiency = 1.0;  // by definition
      return edgeColoring;
    }

    /*--- Create a temporary sparse pattern from the edges. ---*/
    su2vector<unsigned long> outerPtr(nEdge + 1);
    su2vector<unsigned long> innerIdx(nEdge * 2);

    for (unsigned long iEdge = 0; iEdge < nEdge; ++iEdge) {
      outerPtr(iEdge) = 2 * iEdge;
      innerIdx(iEdge * 2 + 0) = edges->GetNode(iEdge, 0);
      innerIdx(iEdge * 2 + 1) = edges->GetNode(iEdge, 1);
    }
    outerPtr(nEdge) = 2 * nEdge;

    CCompressedSparsePatternUL pattern(move(outerPtr), move(innerIdx));

    /*--- Color the edges. ---*/
    constexpr bool balanceColors = true;

    /*--- If requested, find an efficient coloring with maximum color group size (up to edgeColorGroupSize). ---*/
    if (maximizeEdgeColorGroupSize) {
      auto upperEdgeColorGroupSize = edgeColorGroupSize + 1; /* upper bound that is deemed too large */
      auto nextEdgeColorGroupSize = edgeColorGroupSize;      /* next value that we are going to try */
      auto lowerEdgeColorGroupSize = 1ul;                    /* lower bound that is known to work */

      bool admissibleColoring = false; /* keep track wether the last tested coloring is admissible */

      while (true) {
        edgeColoring = colorSparsePattern(pattern, nextEdgeColorGroupSize, false, balanceColors);

        /*--- If the coloring fails, reduce the color group size. ---*/
        if (edgeColoring.empty()) {
          upperEdgeColorGroupSize = nextEdgeColorGroupSize;
          admissibleColoring = false;
        }
        /*--- If the coloring succeeds, check the efficiency. ---*/
        else {
          const su2double currentEfficiency =
              coloringEfficiency(edgeColoring, omp_get_max_threads(), nextEdgeColorGroupSize);

          /*--- If the coloring is not efficient, reduce the color group size. ---*/
          if (currentEfficiency < COLORING_EFF_THRESH) {
            upperEdgeColorGroupSize = nextEdgeColorGroupSize;
            admissibleColoring = false;
          }
          /*--- Otherwise, enlarge the color group size. ---*/
          else {
            lowerEdgeColorGroupSize = nextEdgeColorGroupSize;
            admissibleColoring = true;
          }
        }

        const auto increment = (upperEdgeColorGroupSize - lowerEdgeColorGroupSize) / 2;
        nextEdgeColorGroupSize = lowerEdgeColorGroupSize + increment;

        /*--- Terminating condition. ---*/
        if (increment == 0) {
          break;
        }
      }

      edgeColorGroupSize = nextEdgeColorGroupSize;

      /*--- If the last tested coloring was not admissible, recompute the final coloring. ---*/
      if (!admissibleColoring) {
        edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, false, balanceColors);
      }
    }
    /*--- No adaptivity. ---*/
    else {
      edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, false, balanceColors);
    }

    /*--- If the coloring fails use the natural coloring. This is a
     *    "soft" failure as this "bad" coloring should be detected
     *    downstream and a fallback strategy put in place. ---*/
    if (edgeColoring.empty()) SetNaturalEdgeColoring();
  }

  if (efficiency != nullptr) {
    *efficiency = coloringEfficiency(edgeColoring, omp_get_max_threads(), edgeColorGroupSize);
  }
  return edgeColoring;
}

void CGeometry::SetNaturalEdgeColoring() {
  if (nEdge == 0) return;
  edgeColoring = createNaturalColoring(nEdge);
  /*--- In parallel, set the group size to nEdge to protect client code. ---*/
  if (omp_get_max_threads() > 1) edgeColorGroupSize = nEdge;
}

const CCompressedSparsePatternUL& CGeometry::GetElementColoring(su2double* efficiency) {
  /*--- Check for dry run mode with dummy geometry. ---*/
  if (nElem == 0) return elemColoring;

  /*--- Build if required. ---*/
  if (elemColoring.empty()) {
    /*--- When not using threading use the natural coloring. ---*/
    if (omp_get_max_threads() == 1) {
      SetNaturalElementColoring();
      if (efficiency != nullptr) *efficiency = 1.0;  // by definition
      return elemColoring;
    }

    /*--- Create a temporary sparse pattern from the elements. ---*/
    vector<unsigned long> outerPtr(nElem + 1);
    vector<unsigned long> innerIdx;
    innerIdx.reserve(nElem);

    for (unsigned long iElem = 0; iElem < nElem; ++iElem) {
      outerPtr[iElem] = innerIdx.size();

      for (unsigned short iNode = 0; iNode < elem[iElem]->GetnNodes(); ++iNode) {
        innerIdx.push_back(elem[iElem]->GetNode(iNode));
      }
    }
    outerPtr[nElem] = innerIdx.size();

    CCompressedSparsePatternUL pattern(outerPtr, innerIdx);

    /*--- Color the elements. ---*/
    constexpr bool balanceColors = true;
    elemColoring = colorSparsePattern(pattern, elemColorGroupSize, false, balanceColors);

    /*--- Same as for the edge coloring. ---*/
    if (elemColoring.empty()) SetNaturalElementColoring();
  }

  if (efficiency != nullptr) {
    *efficiency = coloringEfficiency(elemColoring, omp_get_max_threads(), elemColorGroupSize);
  }
  return elemColoring;
}

void CGeometry::SetNaturalElementColoring() {
  if (nElem == 0) return;
  elemColoring = createNaturalColoring(nElem);
  /*--- In parallel, set the group size to nElem to protect client code. ---*/
  if (omp_get_max_threads() > 1) elemColorGroupSize = nElem;
}

void CGeometry::ColorMGLevels(unsigned short nMGLevels, const CGeometry* const* geometry) {
  using tColor = uint8_t;
  constexpr auto nColor = numeric_limits<tColor>::max();

  if (nMGLevels) CoarseGridColor_.resize(nPoint, nMGLevels) = 0;

  for (auto iMesh = nMGLevels; iMesh >= 1; --iMesh) {
    /*--- Color the coarse points. ---*/
    vector<tColor> color;
    const auto& adjacency = geometry[iMesh]->nodes->GetPoints();
    if (colorSparsePattern<tColor, nColor>(adjacency, 1, true, false, &color).empty()) continue;

    /*--- Propagate colors to fine mesh. ---*/
    for (auto step = 0u; step < iMesh; ++step) {
      auto coarseMesh = geometry[iMesh - 1 - step];
      if (step)
        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint) {
          CoarseGridColor_(iPoint, step) = CoarseGridColor_(coarseMesh->nodes->GetParent_CV(iPoint), step - 1);
        }
      else
        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint) {
          CoarseGridColor_(iPoint, step) = color[coarseMesh->nodes->GetParent_CV(iPoint)];
        }
    }
  }
}

const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config) const {
  auto& li = lineletInfo;
  /*--- Do not use "linelets.empty()" as the check, some ranks may have no linelets, and this
   * function contains collective communication that all ranks must call the same number of times. ---*/
  if (li.computed) return li;
  li.computed = true;

  li.lineletIdx.resize(nPoint, CLineletInfo::NO_LINELET);

  /*--- Estimate number of linelets. ---*/

  unsigned long nLinelet = 0;
  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
    if (config->GetSolid_Wall(iMarker) || config->GetMarker_All_KindBC(iMarker) == DISPLACEMENT_BOUNDARY) {
      nLinelet += nVertex[iMarker];
    }
  }

  /*--- If the domain contains well defined Linelets ---*/

  unsigned long maxNPoints = 0, sumNPoints = 0;

  if (nLinelet != 0) {
    /*--- Define the basic linelets, starting from each vertex, preventing duplication of points. ---*/

    li.linelets.resize(nLinelet);
    nLinelet = 0;
    for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
      if (config->GetSolid_Wall(iMarker) || config->GetMarker_All_KindBC(iMarker) == DISPLACEMENT_BOUNDARY) {
        for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
          const auto iPoint = vertex[iMarker][iVertex]->GetNode();
          if (li.lineletIdx[iPoint] == CLineletInfo::NO_LINELET && nodes->GetDomain(iPoint)) {
            li.linelets[nLinelet].push_back(iPoint);
            li.lineletIdx[iPoint] = nLinelet;
            ++nLinelet;
          }
        }
      }
    }
    li.linelets.resize(nLinelet);

    /*--- Create the linelet structure. ---*/

    nLinelet = 0;
    for (auto& linelet : li.linelets) {
      while (linelet.size() < CLineletInfo::MAX_LINELET_POINTS) {
        const auto iPoint = linelet.back();

        /*--- Compute the value of the max and min weights to detect if this region is isotropic. ---*/

        su2double max_weight = 0.0, min_weight = std::numeric_limits<su2double>::max();
        for (auto iNode = 0u; iNode < nodes->GetnPoint(iPoint); iNode++) {
          const auto jPoint = nodes->GetPoint(iPoint, iNode);
          const auto iEdge = nodes->GetEdge(iPoint, iNode);
          const auto* normal = edges->GetNormal(iEdge);
          const su2double area = GeometryToolbox::Norm(nDim, normal);
          const su2double volume_iPoint = nodes->GetVolume(iPoint);
          const su2double volume_jPoint = nodes->GetVolume(jPoint);
          const su2double weight = 0.5 * area * (1.0 / volume_iPoint + 1.0 / volume_jPoint);
          max_weight = max(max_weight, weight);
          min_weight = min(min_weight, weight);
        }

        /*--- Isotropic, stop this linelet. ---*/
        if (min_weight / max_weight > CLineletInfo::ALPHA_ISOTROPIC()) break;

        /*--- Otherwise, add the closest valid neighbor. ---*/

        su2double min_dist2 = std::numeric_limits<su2double>::max();
        auto next_Point = iPoint;
        const auto* iCoord = nodes->GetCoord(iPoint);

        for (const auto jPoint : nodes->GetPoints(iPoint)) {
          if (li.lineletIdx[jPoint] == CLineletInfo::NO_LINELET && nodes->GetDomain(jPoint)) {
            const auto* jCoord = nodes->GetCoord(jPoint);
            const su2double d2 = GeometryToolbox::SquaredDistance(nDim, iCoord, jCoord);
            su2double cosTheta = 1;
            if (linelet.size() > 1) {
              const auto* kCoord = nodes->GetCoord(linelet[linelet.size() - 2]);
              su2double dij[3] = {0.0}, dki[3] = {0.0};
              GeometryToolbox::Distance(nDim, iCoord, kCoord, dki);
              GeometryToolbox::Distance(nDim, jCoord, iCoord, dij);
              cosTheta = GeometryToolbox::DotProduct(3, dki, dij) / sqrt(d2 * GeometryToolbox::SquaredNorm(nDim, dki));
            }
            if (d2 < min_dist2 && cosTheta > 0.7071) {
              next_Point = jPoint;
              min_dist2 = d2;
            }
          }
        }

        /*--- Did not find a suitable point. ---*/
        if (next_Point == iPoint) break;

        linelet.push_back(next_Point);
        li.lineletIdx[next_Point] = nLinelet;
      }
      ++nLinelet;

      maxNPoints = max<unsigned long>(maxNPoints, linelet.size());
      sumNPoints += linelet.size();
    }
  }

  /*--- Average linelet size over all ranks. ---*/

  unsigned long globalNPoints, globalNLineLets;
  SU2_MPI::Allreduce(&sumNPoints, &globalNPoints, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
  SU2_MPI::Allreduce(&nLinelet, &globalNLineLets, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());

  if (rank == MASTER_NODE && globalNLineLets > 0) {
    std::cout << "Computed linelet structure, "
              << static_cast<unsigned long>(passivedouble(globalNPoints) / globalNLineLets)
              << " points in each line (average)." << std::endl;
  }

  /*--- Color the linelets for OpenMP parallelization and visualization. ---*/

  /*--- Adjacency between lines, computed from point neighbors. ---*/
  std::vector<std::vector<unsigned long>> adjacency(nLinelet);
  for (auto iLine = 0ul; iLine < nLinelet; ++iLine) {
    std::unordered_set<unsigned long> neighbors;
    for (const auto iPoint : li.linelets[iLine]) {
      neighbors.insert(iPoint);
      for (const auto jPoint : nodes->GetPoints(iPoint)) {
        neighbors.insert(jPoint);
      }
    }
    adjacency[iLine].reserve(neighbors.size());
    for (const auto iPoint : neighbors) {
      adjacency[iLine].push_back(iPoint);
    }
  }

  const auto coloring = colorSparsePattern<uint8_t, std::numeric_limits<uint8_t>::max()>(
      CCompressedSparsePatternUL(adjacency), 1, false, true);
  /*--- unsigned long because it is communicated as MPI_UNSIGNED_LONG below. ---*/
  const unsigned long nColors = coloring.getOuterSize();

  /*--- Sort linelets by color. ---*/
  std::vector<std::vector<unsigned long>> sortedLinelets;
  sortedLinelets.reserve(nLinelet);
  li.colorOffsets.reserve(nColors + 1);
  li.colorOffsets.push_back(0);

  for (auto iColor = 0ul; iColor < nColors; ++iColor) {
    for (const auto iLine : coloring.getInnerIter(iColor)) {
      /*--- Store the new linelet index for its points. ---*/
      for (const auto iPoint : li.linelets[iLine]) {
        li.lineletIdx[iPoint] = sortedLinelets.size();
      }
      sortedLinelets.push_back(std::move(li.linelets[iLine]));
    }
    li.colorOffsets.push_back(sortedLinelets.size());
  }
  li.linelets = std::move(sortedLinelets);

  /*--- For visualization, offset colors to avoid coloring across ranks. ---*/
  std::vector<unsigned long> allNColors(size);
  SU2_MPI::Allgather(&nColors, 1, MPI_UNSIGNED_LONG, allNColors.data(), 1, MPI_UNSIGNED_LONG, SU2_MPI::GetComm());
  unsigned long offset = 0;
  for (int i = 0; i < rank; ++i) offset += allNColors[i];

  /*--- Finally, transfer colors to points, using "0" as "no linelet". ---*/
  li.lineletColor.resize(nPoint, 0);
  for (auto iColor = 0ul; iColor < nColors; ++iColor) {
    for (auto iLine = li.colorOffsets[iColor]; iLine < li.colorOffsets[iColor + 1]; ++iLine) {
      for (const auto iPoint : li.linelets[iLine]) {
        li.lineletColor[iPoint] = 1 + offset + iColor;
      }
    }
  }

  return li;
}

namespace {
su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, const unsigned long iPoint) {
  const su2double max = std::numeric_limits<su2double>::max();
  su2double distance = max;
  for (const auto jPoint : geometry->nodes->GetPoints(iPoint)) {
    const su2double dist =
        geometry->nodes->GetViscousBoundary(jPoint) ? 0.0 : geometry->nodes->GetWall_Distance(jPoint);
    if (dist > EPS) distance = fmin(distance, dist);
  }
  if (distance > 0 && distance < max) return distance;

  /*--- The point only has wall neighbors, which all have 0 wall distance.
   *    Compute an alternative distance based on volume and wall area. ---*/

  const auto nDim = geometry->GetnDim();
  su2double Normal[3] = {};
  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
    if (!config->GetViscous_Wall(iMarker)) continue;

    const auto iVertex = geometry->nodes->GetVertex(iPoint, iMarker);
    if (iVertex < 0) continue;

    for (auto iDim = 0u; iDim < nDim; iDim++) {
      Normal[iDim] += geometry->vertex[iMarker][iVertex]->GetNormal(iDim);
    }
  }
  const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint);
  return 2 * Vol / GeometryToolbox::Norm(3, Normal);
}
}  // namespace

void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container,
                                    const int record_zone) {
  int nZone = config_container[ZONE_0]->GetnZone();
  bool allEmpty = true;
  vector<bool> wallDistanceNeeded(nZone, false);

  for (int iInst = 0; iInst < config_container[ZONE_0]->GetnTimeInstances(); iInst++) {
    for (int iZone = 0; iZone < nZone; iZone++) {
      /*--- Check if a zone needs the wall distance and store a boolean ---*/

      MAIN_SOLVER kindSolver = config_container[iZone]->GetKind_Solver();
      if (kindSolver == MAIN_SOLVER::RANS || kindSolver == MAIN_SOLVER::INC_RANS ||
          kindSolver == MAIN_SOLVER::DISC_ADJ_RANS || kindSolver == MAIN_SOLVER::DISC_ADJ_INC_RANS ||
          kindSolver == MAIN_SOLVER::FEM_LES || kindSolver == MAIN_SOLVER::FEM_RANS) {
        wallDistanceNeeded[iZone] = true;
      }

      /*--- Set the wall distances in all zones to the numerical limit.
       * This is necessary, because before a computed distance is set, it will be checked
       * whether the new distance is smaller than the currently stored one. ---*/
      CGeometry* geometry = geometry_container[iZone][iInst][MESH_0];
      if (wallDistanceNeeded[iZone]) geometry->SetWallDistance(numeric_limits<su2double>::max());
    }

    /*--- Loop over all zones and compute the ADT based on the viscous walls in that zone ---*/
    for (int iZone = 0; iZone < nZone; iZone++) {
      unique_ptr<CADTElemClass> WallADT =
          geometry_container[iZone][iInst][MESH_0]->ComputeViscousWallADT(config_container[iZone]);
      if (WallADT && !WallADT->IsEmpty()) {
        allEmpty = false;
        /*--- Inner loop over all zones to update the wall distances.
         * It might happen that there is a closer viscous wall in zone iZone for points in zone jZone. ---*/
        for (int jZone = 0; jZone < nZone; jZone++) {
          if (wallDistanceNeeded[jZone])
            geometry_container[jZone][iInst][MESH_0]->SetWallDistance(WallADT.get(), config_container[jZone], iZone);
        }
      }
    }

    /*--- If there are no viscous walls in the entire domain, set distances to zero ---*/
    if (allEmpty) {
      for (int iZone = 0; iZone < nZone; iZone++) {
        CGeometry* geometry = geometry_container[iZone][iInst][MESH_0];
        geometry->SetWallDistance(0.0);
      }
      continue;
    }

    /*--- Otherwise, set wall roughnesses, storing them in a common data structure. ---*/
    // [iZone][iMarker] -> roughness, for this rank
    auto roughness_f = make_pair(nZone, [config_container, geometry_container, iInst](unsigned long iZone) {
      const CConfig* config = config_container[iZone];
      const auto nMarker = geometry_container[iZone][iInst][MESH_0]->GetnMarker();

      return make_pair(nMarker, [config](unsigned long iMarker) {
        return config->GetWallRoughnessProperties(config->GetMarker_All_TagBound(iMarker)).second;
      });
    });
    NdFlattener<2> roughness_local(roughness_f);
    // [rank][iZone][iMarker] -> roughness
    NdFlattener<3> roughness_global(Nd_MPI_Environment(), roughness_local);
    // use it to update roughnesses
    for (int jZone = 0; jZone < nZone; jZone++) {
      if (wallDistanceNeeded[jZone] && config_container[jZone]->GetnRoughWall() > 0) {
        geometry_container[jZone][iInst][MESH_0]->nodes->SetWallRoughness(roughness_global);
      }
    }

    for (int iZone = 0; iZone < nZone; iZone++) {
      /*--- When recording for a specific zone, only compute nearest-neighbor distances (which read vertex
       *    normals) for that zone. Reading normals from other zones at this tape position would create
       *    cross-zone AD dependencies before those zones have had their geometry updated. ---*/
      if (record_zone >= 0 && iZone != record_zone) continue;

      /*--- For the FEM solver, we use a different mesh structure ---*/
      MAIN_SOLVER kindSolver = config_container[iZone]->GetKind_Solver();
      if (!wallDistanceNeeded[iZone] || kindSolver == MAIN_SOLVER::FEM_LES || kindSolver == MAIN_SOLVER::FEM_RANS) {
        continue;
      }
      const auto* config = config_container[iZone];
      auto* geometry = geometry_container[iZone][iInst][MESH_0];

      for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); ++iMarker) {
        const auto viscous = config->GetViscous_Wall(iMarker);

        SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
        for (auto iVertex = 0u; iVertex < geometry->nVertex[iMarker]; iVertex++) {
          const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
          su2double dist = 0;
          if (viscous && geometry->nodes->GetDomain(iPoint)) {
            dist = NearestNeighborDistance(geometry, config, iPoint);
          } else {
            dist = geometry->nodes->GetWall_Distance(iPoint);
          }
          geometry->vertex[iMarker][iVertex]->SetNearestNeighborDistance(dist);
        }
        END_SU2_OMP_FOR
      }
    }
  }
}