From 5271d3e32a1d0588dc9eb56742925b2f42877029 Mon Sep 17 00:00:00 2001 From: vdweide Date: Sat, 8 Feb 2025 09:45:44 +0100 Subject: [PATCH 1/9] Created some fem toolbox classes --- Common/include/fem/fem_geometry_structure.hpp | 146 +------- Common/include/geometry/CGeometry.hpp | 1 - Common/include/geometry/CPhysicalGeometry.hpp | 2 + .../geometry/meshreader/CMeshReaderBase.hpp | 2 +- .../toolboxes/classes_multiple_integers.hpp | 97 +++++ .../fem/CFaceOfElement.hpp} | 123 +------ .../include/toolboxes/fem/CMatchingFace.hpp | 76 ++++ .../toolboxes/fem/CReorderElements.hpp | 101 ++++++ .../toolboxes/fem/CSortBoundaryFaces.hpp | 48 +++ Common/include/toolboxes/fem/CSortFaces.hpp | 71 ++++ Common/src/CConfig.cpp | 3 +- Common/src/fem/fem_geometry_structure.cpp | 138 +------- Common/src/fem/meson.build | 3 +- .../CPhysicalGeometryFEM.cpp} | 330 +----------------- Common/src/geometry/meson.build | 1 + Common/src/toolboxes/fem/CFaceOfElement.cpp | 243 +++++++++++++ Common/src/toolboxes/fem/CMatchingFace.cpp | 128 +++++++ Common/src/toolboxes/fem/CReorderElements.cpp | 65 ++++ .../src/toolboxes/fem/CSortBoundaryFaces.cpp | 42 +++ Common/src/toolboxes/fem/CSortFaces.cpp | 108 ++++++ Common/src/toolboxes/fem/meson.build | 5 + Common/src/toolboxes/meson.build | 1 + 22 files changed, 1011 insertions(+), 723 deletions(-) create mode 100644 Common/include/toolboxes/classes_multiple_integers.hpp rename Common/include/{fem/geometry_structure_fem_part.hpp => toolboxes/fem/CFaceOfElement.hpp} (54%) create mode 100644 Common/include/toolboxes/fem/CMatchingFace.hpp create mode 100644 Common/include/toolboxes/fem/CReorderElements.hpp create mode 100644 Common/include/toolboxes/fem/CSortBoundaryFaces.hpp create mode 100644 Common/include/toolboxes/fem/CSortFaces.hpp rename Common/src/{fem/geometry_structure_fem_part.cpp => geometry/CPhysicalGeometryFEM.cpp} (92%) create mode 100644 Common/src/toolboxes/fem/CFaceOfElement.cpp create mode 100644 Common/src/toolboxes/fem/CMatchingFace.cpp create mode 100644 Common/src/toolboxes/fem/CReorderElements.cpp create mode 100644 Common/src/toolboxes/fem/CSortBoundaryFaces.cpp create mode 100644 Common/src/toolboxes/fem/CSortFaces.cpp create mode 100644 Common/src/toolboxes/fem/meson.build diff --git a/Common/include/fem/fem_geometry_structure.hpp b/Common/include/fem/fem_geometry_structure.hpp index ab4cc457f40..864db542caa 100644 --- a/Common/include/fem/fem_geometry_structure.hpp +++ b/Common/include/fem/fem_geometry_structure.hpp @@ -32,152 +32,10 @@ #include "fem_standard_element.hpp" #include "../wall_model.hpp" #include "../linear_algebra/blas_structure.hpp" +#include "../toolboxes/fem/CFaceOfElement.hpp" using namespace std; -/*! - * \class CLong3T - * \brief Help class used to store three longs as one entity. - * \version 8.1.0 "Harrier" - */ -struct CLong3T { - long long0 = 0; /*!< \brief First long to store in this class. */ - long long1 = 0; /*!< \brief Second long to store in this class. */ - long long2 = 0; /*!< \brief Third long to store in this class. */ - - CLong3T() = default; - - CLong3T(const long a, const long b, const long c) { - long0 = a; - long1 = b; - long2 = c; - } - - bool operator<(const CLong3T& other) const; -}; - -/*! - * \class CReorderElements - * \brief Class, used to reorder the owned elements after the partitioning. - * \author E. van der Weide - * \version 8.1.0 "Harrier" - */ -class CReorderElements { - private: - unsigned long globalElemID; /*!< \brief Global element ID of the element. */ - unsigned short timeLevel; /*!< \brief Time level of the element. Only relevant - for time accurate local time stepping. */ - bool commSolution; /*!< \brief Whether or not the solution must be - communicated to other ranks. */ - unsigned short elemType; /*!< \brief Short hand for the element type, Which - stored info of the VTK_Type, polynomial - degree of the solution and whether or - not the Jacobian is constant. */ - public: - /*! - * \brief Constructor of the class, set the member variables to the arguments. - */ - CReorderElements(const unsigned long val_GlobalElemID, const unsigned short val_TimeLevel, - const bool val_CommSolution, const unsigned short val_VTK_Type, const unsigned short val_nPolySol, - const bool val_JacConstant); - - /*! - * \brief Default constructor of the class. Disabled. - */ - CReorderElements(void) = delete; - - /*! - * \brief Less than operator of the class. Needed for the sorting. - */ - bool operator<(const CReorderElements& other) const; - - /*! - * \brief Function to make available the variable commSolution. - * \return Whether or not the solution of the element must be communicated. - */ - inline bool GetCommSolution(void) const { return commSolution; } - - /*! - * \brief Function to make available the element type of the element. - * \return The value of elemType, which stores the VTK type, polynomial degree - and whether or not the Jacobian is constant. - */ - inline unsigned short GetElemType(void) const { return elemType; } - - /*! - * \brief Function to make available the global element ID. - * \return The global element ID of the element. - */ - inline unsigned long GetGlobalElemID(void) const { return globalElemID; } - - /*! - * \brief Function to make available the time level. - * \return The time level of the element. - */ - inline unsigned short GetTimeLevel(void) const { return timeLevel; } - - /*! - * \brief Function, which sets the value of commSolution. - * \param[in] val_CommSolution - value to which commSolution must be set. - */ - inline void SetCommSolution(const bool val_CommSolution) { commSolution = val_CommSolution; } -}; - -/*! - * \class CSortFaces - * \brief Functor, used for a different sorting of the faces than the < operator - * of CFaceOfElement. - * \author E. van der Weide - * \version 8.1.0 "Harrier" - */ -class CVolumeElementFEM; // Forward declaration to avoid problems. -class CSortFaces { - private: - unsigned long nVolElemOwned; /*!< \brief Number of locally owned volume elements. */ - unsigned long nVolElemTot; /*!< \brief Total number of local volume elements . */ - - const CVolumeElementFEM* volElem; /*!< \brief The locally stored volume elements. */ - - public: - /*! - * \brief Constructor of the class. Set the values of the member variables. - */ - CSortFaces(unsigned long val_nVolElemOwned, unsigned long val_nVolElemTot, const CVolumeElementFEM* val_volElem) { - nVolElemOwned = val_nVolElemOwned; - nVolElemTot = val_nVolElemTot; - volElem = val_volElem; - } - - /*! - * \brief Default constructor of the class. Disabled. - */ - CSortFaces(void) = delete; - - /*! - * \brief Operator used for the comparison. - * \param[in] f0 - First face in the comparison. - * \param[in] f1 - Second face in the comparison. - */ - bool operator()(const CFaceOfElement& f0, const CFaceOfElement& f1); -}; - -/*! - * \class CSortBoundaryFaces - * \brief Functor, used for a different sorting of the faces than the < operator - * of CSurfaceElementFEM. - * \author E. van der Weide - * \version 8.1.0 "Harrier" - */ -struct CSurfaceElementFEM; // Forward declaration to avoid problems. -struct CSortBoundaryFaces { - /*! - * \brief Operator used for the comparison. - * \param[in] f0 - First boundary face in the comparison. - * \param[in] f1 - Second boundary face in the comparison. - */ - bool operator()(const CSurfaceElementFEM& f0, const CSurfaceElementFEM& f1); -}; - /*! * \class CVolumeElementFEM * \brief Class to store a volume element for the FEM solver. @@ -271,7 +129,7 @@ class CVolumeElementFEM { /*! * \class CPointFEM - * \brief Class to a point for the FEM solver. + * \brief Class to store a point for the FEM solver. * \author E. van der Weide * \version 8.1.0 "Harrier" */ diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index dd79546f56b..3a05f446e09 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -60,7 +60,6 @@ extern "C" { #include "dual_grid/CTurboVertex.hpp" #include "../CConfig.hpp" -#include "../fem/geometry_structure_fem_part.hpp" #include "../toolboxes/graph_toolbox.hpp" #include "../adt/CADTElemClass.hpp" diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index e6a871bb977..7d2f7d31f09 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/Common/include/geometry/CPhysicalGeometry.hpp @@ -30,6 +30,8 @@ #include "CGeometry.hpp" #include "meshreader/CMeshReaderBase.hpp" #include "../containers/C2DContainer.hpp" +#include "../toolboxes/classes_multiple_integers.hpp" +#include "../toolboxes/fem/CFaceOfElement.hpp" /*! * \class CPhysicalGeometry diff --git a/Common/include/geometry/meshreader/CMeshReaderBase.hpp b/Common/include/geometry/meshreader/CMeshReaderBase.hpp index 622b48db4c8..01a2687ba85 100644 --- a/Common/include/geometry/meshreader/CMeshReaderBase.hpp +++ b/Common/include/geometry/meshreader/CMeshReaderBase.hpp @@ -32,7 +32,7 @@ #include #include "../primal_grid/CPrimalGridFEM.hpp" -#include "../../fem/geometry_structure_fem_part.hpp" +#include "../../toolboxes/fem/CFaceOfElement.hpp" #include "../../parallelization/mpi_structure.hpp" #include "../../CConfig.hpp" diff --git a/Common/include/toolboxes/classes_multiple_integers.hpp b/Common/include/toolboxes/classes_multiple_integers.hpp new file mode 100644 index 00000000000..605cf046505 --- /dev/null +++ b/Common/include/toolboxes/classes_multiple_integers.hpp @@ -0,0 +1,97 @@ +/*! + * \file classes_multiple_integers.hpp + * \brief Header file for the classes that consists of multiple integer types. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#pragma once + +#include + +/*! + * \struct CUnsignedLong2T + * \brief Helper struct used to store two integral types as one entity. + */ +struct CUnsignedLong2T { + unsigned long long0; /*!< \brief First integer to store in this class. */ + unsigned long long1; /*!< \brief Second integer to store in this class. */ + + CUnsignedLong2T(unsigned long a = 0, unsigned long b = 0) : long0(a), long1(b) {} + + inline bool operator<(const CUnsignedLong2T& other) const { + if (long0 != other.long0) return (long0 < other.long0); + return (long1 < other.long1); + } + + inline bool operator==(const CUnsignedLong2T& other) const { + return (long0 == other.long0) && (long1 == other.long1); + } +}; + +/*! + * \struct CUnsignedShort2T + * \brief Help struct used to store two integral types as one entity. + */ +struct CUnsignedShort2T { + unsigned short short0; /*!< \brief First integer to store in this class. */ + unsigned short short1; /*!< \brief Second integer to store in this class. */ + + CUnsignedShort2T(unsigned short a = 0, unsigned short b = 0) : short0(a), short1(b) {} + + inline bool operator<(const CUnsignedShort2T& other) const { + if (short0 != other.short0) return (short0 < other.short0); + return (short1 < other.short1); + } + + inline bool operator==(const CUnsignedShort2T& other) const { + return (short0 == other.short0) && (short1 == other.short1); + } +}; + +/*! + * \class CLong3T + * \brief Help class used to store three longs as one entity. + * \version 8.1.0 "Harrier" + */ +struct CLong3T { + long long0 = 0; /*!< \brief First long to store in this class. */ + long long1 = 0; /*!< \brief Second long to store in this class. */ + long long2 = 0; /*!< \brief Third long to store in this class. */ + + CLong3T() = default; + + CLong3T(const long a, const long b, const long c) { + long0 = a; + long1 = b; + long2 = c; + } + + bool operator<(const CLong3T& other) const { + if (long0 != other.long0) return (long0 < other.long0); + if (long1 != other.long1) return (long1 < other.long1); + if (long2 != other.long2) return (long2 < other.long2); + + return false; + } +}; diff --git a/Common/include/fem/geometry_structure_fem_part.hpp b/Common/include/toolboxes/fem/CFaceOfElement.hpp similarity index 54% rename from Common/include/fem/geometry_structure_fem_part.hpp rename to Common/include/toolboxes/fem/CFaceOfElement.hpp index 4cf2ab12d1f..57f57ac7110 100644 --- a/Common/include/fem/geometry_structure_fem_part.hpp +++ b/Common/include/toolboxes/fem/CFaceOfElement.hpp @@ -1,6 +1,7 @@ /*! - * \file geometry_structure_fem_part.hpp - * \brief Helper classes for the Fluid FEM solver. + * \file CFaceOfElement.hpp + * \brief Header file for the class CFaceOfElement. + * The implementations are in the CFaceOfElement.cpp file. * \author E. van der Weide * \version 8.1.0 "Harrier" * @@ -27,50 +28,13 @@ #pragma once -#include "../basic_types/datatype_structure.hpp" - -#include +#include "../../code_config.hpp" +#include "../../option_structure.hpp" +#include #include +#include -/*! - * \struct CUnsignedLong2T - * \brief Helper struct used to store two integral types as one entity. - */ -struct CUnsignedLong2T { - unsigned long long0; /*!< \brief First integer to store in this class. */ - unsigned long long1; /*!< \brief Second integer to store in this class. */ - - CUnsignedLong2T(unsigned long a = 0, unsigned long b = 0) : long0(a), long1(b) {} - - inline bool operator<(const CUnsignedLong2T& other) const { - if (long0 != other.long0) return (long0 < other.long0); - return (long1 < other.long1); - } - - inline bool operator==(const CUnsignedLong2T& other) const { - return (long0 == other.long0) && (long1 == other.long1); - } -}; - -/*! - * \struct CUnsignedShort2T - * \brief Help struct used to store two integral types as one entity. - */ -struct CUnsignedShort2T { - unsigned short short0; /*!< \brief First integer to store in this class. */ - unsigned short short1; /*!< \brief Second integer to store in this class. */ - - CUnsignedShort2T(unsigned short a = 0, unsigned short b = 0) : short0(a), short1(b) {} - - inline bool operator<(const CUnsignedShort2T& other) const { - if (short0 != other.short0) return (short0 < other.short0); - return (short1 < other.short1); - } - - inline bool operator==(const CUnsignedShort2T& other) const { - return (short0 == other.short0) && (short1 == other.short1); - } -}; +using namespace std; /*! * \class CFaceOfElement @@ -137,74 +101,3 @@ class CFaceOfElement { /*--- Copy function, which copies the data of the given object into the current object. ---*/ void Copy(const CFaceOfElement& other); }; - -/*! - * \class CBoundaryFace - * \brief Help class used in the partitioning of the FEM grid. - It stores a boundary element. - */ -class CBoundaryFace { - public: - unsigned short VTK_Type, nPolyGrid, nDOFsGrid; - unsigned long globalBoundElemID, domainElementID; - std::vector Nodes; - - /* Standard constructor and destructor. Nothing to be done. */ - CBoundaryFace() {} - ~CBoundaryFace() {} - - /* Copy constructor and assignment operator. */ - inline CBoundaryFace(const CBoundaryFace& other) { Copy(other); } - - inline CBoundaryFace& operator=(const CBoundaryFace& other) { - Copy(other); - return (*this); - } - - /* Less than operator. Needed for the sorting. */ - inline bool operator<(const CBoundaryFace& other) const { return (globalBoundElemID < other.globalBoundElemID); } - - private: - /*--- Copy function, which copies the data of the given object into the current object. ---*/ - void Copy(const CBoundaryFace& other); -}; - -/*! - * \class CMatchingFace - * \brief Help class used to determine whether or not (periodic) faces match. - */ -class CMatchingFace { - public: - unsigned short nCornerPoints; /*!< \brief Number of corner points of the face. */ - unsigned short nDim; /*!< \brief Number of spatial dimensions. */ - unsigned short nPoly; /*!< \brief Polynomial degree of the face. */ - unsigned short nDOFsElem; /*!< \brief Number of DOFs of the relevant adjacent element. */ - unsigned short elemType; /*!< \brief Type of the adjacent element. */ - unsigned long elemID; /*!< \brief The relevant adjacent element ID. */ - su2double cornerCoor[4][3]; /*!< \brief Coordinates of the corner points of the face. */ - su2double tolForMatching; /*!< \brief Tolerance for this face for matching points. */ - - /* Standard constructor. */ - CMatchingFace(); - - /* Destructor, nothing to be done. */ - ~CMatchingFace() {} - - /* Copy constructor and assignment operator. */ - inline CMatchingFace(const CMatchingFace& other) { Copy(other); } - - inline CMatchingFace& operator=(const CMatchingFace& other) { - Copy(other); - return (*this); - } - - /* Less than operator. Needed for the sorting and searching. */ - bool operator<(const CMatchingFace& other) const; - - /*--- Member function, which sorts the coordinates of the face. ---*/ - void SortFaceCoordinates(void); - - private: - /*--- Copy function, which copies the data of the given object into the current object. ---*/ - void Copy(const CMatchingFace& other); -}; diff --git a/Common/include/toolboxes/fem/CMatchingFace.hpp b/Common/include/toolboxes/fem/CMatchingFace.hpp new file mode 100644 index 00000000000..9f68bbddd5a --- /dev/null +++ b/Common/include/toolboxes/fem/CMatchingFace.hpp @@ -0,0 +1,76 @@ +/*! + * \file CMatchingFace.hpp + * \brief Header file for the class CMatchingFace. + * The implementations are in the CMatchingFace.cpp file. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#pragma once + +#include "../../code_config.hpp" +#include +#include +#include + +using namespace std; + +/*! + * \class CMatchingFace + * \brief Help class used to determine whether or not (periodic) faces match. + */ +class CMatchingFace { + public: + unsigned short nCornerPoints; /*!< \brief Number of corner points of the face. */ + unsigned short nDim; /*!< \brief Number of spatial dimensions. */ + unsigned short nPoly; /*!< \brief Polynomial degree of the face. */ + unsigned short nDOFsElem; /*!< \brief Number of DOFs of the relevant adjacent element. */ + unsigned short elemType; /*!< \brief Type of the adjacent element. */ + unsigned long elemID; /*!< \brief The relevant adjacent element ID. */ + su2double cornerCoor[4][3]; /*!< \brief Coordinates of the corner points of the face. */ + su2double tolForMatching; /*!< \brief Tolerance for this face for matching points. */ + + /* Standard constructor. */ + CMatchingFace(); + + /* Destructor, nothing to be done. */ + ~CMatchingFace() {} + + /* Copy constructor and assignment operator. */ + inline CMatchingFace(const CMatchingFace& other) { Copy(other); } + + inline CMatchingFace& operator=(const CMatchingFace& other) { + Copy(other); + return (*this); + } + + /* Less than operator. Needed for the sorting and searching. */ + bool operator<(const CMatchingFace& other) const; + + /*--- Member function, which sorts the coordinates of the face. ---*/ + void SortFaceCoordinates(void); + + private: + /*--- Copy function, which copies the data of the given object into the current object. ---*/ + void Copy(const CMatchingFace& other); +}; diff --git a/Common/include/toolboxes/fem/CReorderElements.hpp b/Common/include/toolboxes/fem/CReorderElements.hpp new file mode 100644 index 00000000000..a7408aa3d37 --- /dev/null +++ b/Common/include/toolboxes/fem/CReorderElements.hpp @@ -0,0 +1,101 @@ +/*! + * \file CReorderElements.hpp + * \brief Header file for the class CReorderElements. + * The implementations are in the CReorderElements.cpp file. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#pragma once + + +#include + +using namespace std; + +/*! + * \class CReorderElements + * \brief Class, used to reorder the owned elements after the partitioning. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + */ +class CReorderElements { + private: + unsigned long globalElemID; /*!< \brief Global element ID of the element. */ + unsigned short timeLevel; /*!< \brief Time level of the element. Only relevant + for time accurate local time stepping. */ + bool commSolution; /*!< \brief Whether or not the solution must be + communicated to other ranks. */ + unsigned short elemType; /*!< \brief Short hand for the element type, Which + stored info of the VTK_Type, polynomial + degree of the solution and whether or + not the Jacobian is constant. */ + public: + /*! + * \brief Constructor of the class, set the member variables to the arguments. + */ + CReorderElements(const unsigned long val_GlobalElemID, const unsigned short val_TimeLevel, + const bool val_CommSolution, const unsigned short val_VTK_Type, const unsigned short val_nPolySol, + const bool val_JacConstant); + + /*! + * \brief Default constructor of the class. Disabled. + */ + CReorderElements(void) = delete; + + /*! + * \brief Less than operator of the class. Needed for the sorting. + */ + bool operator<(const CReorderElements& other) const; + + /*! + * \brief Function to make available the variable commSolution. + * \return Whether or not the solution of the element must be communicated. + */ + inline bool GetCommSolution(void) const { return commSolution; } + + /*! + * \brief Function to make available the element type of the element. + * \return The value of elemType, which stores the VTK type, polynomial degree + and whether or not the Jacobian is constant. + */ + inline unsigned short GetElemType(void) const { return elemType; } + + /*! + * \brief Function to make available the global element ID. + * \return The global element ID of the element. + */ + inline unsigned long GetGlobalElemID(void) const { return globalElemID; } + + /*! + * \brief Function to make available the time level. + * \return The time level of the element. + */ + inline unsigned short GetTimeLevel(void) const { return timeLevel; } + + /*! + * \brief Function, which sets the value of commSolution. + * \param[in] val_CommSolution - value to which commSolution must be set. + */ + inline void SetCommSolution(const bool val_CommSolution) { commSolution = val_CommSolution; } +}; diff --git a/Common/include/toolboxes/fem/CSortBoundaryFaces.hpp b/Common/include/toolboxes/fem/CSortBoundaryFaces.hpp new file mode 100644 index 00000000000..8271a084a99 --- /dev/null +++ b/Common/include/toolboxes/fem/CSortBoundaryFaces.hpp @@ -0,0 +1,48 @@ +/*! + * \file CSortBoundaryFaces.hpp + * \brief Header file for the class CSortBoundaryFaces. + * The implementations are in the CSortBoundaryFaces.cpp file. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#pragma once + +#include + +/*! + * \class CSortBoundaryFaces + * \brief Functor, used for a different sorting of the faces than the < operator + * of CSurfaceElementFEM. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + */ +struct CSurfaceElementFEM; // Forward declaration to avoid problems. +struct CSortBoundaryFaces { + /*! + * \brief Operator used for the comparison. + * \param[in] f0 - First boundary face in the comparison. + * \param[in] f1 - Second boundary face in the comparison. + */ + bool operator()(const CSurfaceElementFEM& f0, const CSurfaceElementFEM& f1); +}; diff --git a/Common/include/toolboxes/fem/CSortFaces.hpp b/Common/include/toolboxes/fem/CSortFaces.hpp new file mode 100644 index 00000000000..721472b60a0 --- /dev/null +++ b/Common/include/toolboxes/fem/CSortFaces.hpp @@ -0,0 +1,71 @@ +/*! + * \file CSortFaces.hpp + * \brief Header file for the class CSortFaces. + * The implementations are in the CSortFaces.cpp file. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#pragma once + +#include "CFaceOfElement.hpp" + +using namespace std; + +/*! + * \class CSortFaces + * \brief Functor, used for a different sorting of the faces than the < operator + * of CFaceOfElement. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + */ +class CVolumeElementFEM; // Forward declaration to avoid problems. +class CSortFaces { + private: + unsigned long nVolElemOwned; /*!< \brief Number of locally owned volume elements. */ + unsigned long nVolElemTot; /*!< \brief Total number of local volume elements . */ + + const CVolumeElementFEM* volElem; /*!< \brief The locally stored volume elements. */ + + public: + /*! + * \brief Constructor of the class. Set the values of the member variables. + */ + CSortFaces(unsigned long val_nVolElemOwned, unsigned long val_nVolElemTot, const CVolumeElementFEM* val_volElem) { + nVolElemOwned = val_nVolElemOwned; + nVolElemTot = val_nVolElemTot; + volElem = val_volElem; + } + + /*! + * \brief Default constructor of the class. Disabled. + */ + CSortFaces(void) = delete; + + /*! + * \brief Operator used for the comparison. + * \param[in] f0 - First face in the comparison. + * \param[in] f1 - Second face in the comparison. + */ + bool operator()(const CFaceOfElement& f0, const CFaceOfElement& f1); +}; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index d54d1dc3a6c..c61fc1acb2a 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -27,12 +27,13 @@ #define ENABLE_MAPS #include +#include #include "../include/CConfig.hpp" #undef ENABLE_MAPS #include "../include/fem/fem_gauss_jacobi_quadrature.hpp" -#include "../include/fem/fem_geometry_structure.hpp" +#include "../include/toolboxes/classes_multiple_integers.hpp" #include "../include/basic_types/ad_structure.hpp" #include "../include/toolboxes/printing_toolbox.hpp" diff --git a/Common/src/fem/fem_geometry_structure.cpp b/Common/src/fem/fem_geometry_structure.cpp index c6e52dbea81..183d65d28ef 100644 --- a/Common/src/fem/fem_geometry_structure.cpp +++ b/Common/src/fem/fem_geometry_structure.cpp @@ -26,6 +26,10 @@ */ #include "../../include/toolboxes/CLinearPartitioner.hpp" +#include "../../include/toolboxes/classes_multiple_integers.hpp" +#include "../../include/toolboxes/fem/CReorderElements.hpp" +#include "../../include/toolboxes/fem/CSortFaces.hpp" +#include "../../include/toolboxes/fem/CSortBoundaryFaces.hpp" #include "../../include/fem/fem_geometry_structure.hpp" #include "../../include/geometry/primal_grid/CPrimalGridFEM.hpp" #include "../../include/geometry/primal_grid/CPrimalGridBoundFEM.hpp" @@ -38,140 +42,6 @@ extern "C" void dpotrf_(char*, int*, passivedouble*, int*, int*); extern "C" void dpotri_(char*, int*, passivedouble*, int*, int*); #endif -bool CLong3T::operator<(const CLong3T& other) const { - if (long0 != other.long0) return (long0 < other.long0); - if (long1 != other.long1) return (long1 < other.long1); - if (long2 != other.long2) return (long2 < other.long2); - - return false; -} - -CReorderElements::CReorderElements(const unsigned long val_GlobalElemID, const unsigned short val_TimeLevel, - const bool val_CommSolution, const unsigned short val_VTK_Type, - const unsigned short val_nPolySol, const bool val_JacConstant) { - /* Copy the global elment ID, time level and whether or not this element - must be communicated. */ - globalElemID = val_GlobalElemID; - timeLevel = val_TimeLevel; - commSolution = val_CommSolution; - - /* Create the element type used in this class, which stores information of - the VTK type, the polynomial degree of the solution and whether or not the - Jacobian of the transformation is constant. As it is possible that the - polynomial degree of the solution is zero, this convention is different - from the convention used in the SU2 grid file. */ - elemType = val_VTK_Type + 100 * val_nPolySol; - if (!val_JacConstant) elemType += 50; -} - -bool CReorderElements::operator<(const CReorderElements& other) const { - /* Elements with the lowest time level are stored first. */ - if (timeLevel != other.timeLevel) return timeLevel < other.timeLevel; - - /* Next comparison is whether or not the element must communicate its - solution data to other ranks. Elements which do not need to do this - are stored first. */ - if (commSolution != other.commSolution) return other.commSolution; - - /* Elements of the same element type must be stored as contiguously as - possible to allow for the simultaneous treatment of elements in the - matrix multiplications. */ - if (elemType != other.elemType) return elemType < other.elemType; - - /* The final comparison is based on the global element ID. */ - return globalElemID < other.globalElemID; -} - -bool CSortFaces::operator()(const CFaceOfElement& f0, const CFaceOfElement& f1) { - /*--- Comparison in case both faces are boundary faces. ---*/ - if (f0.faceIndicator >= 0 && f1.faceIndicator >= 0) { - /* Both faces are boundary faces. The first comparison is the boundary - marker, which is stored in faceIndicator. */ - if (f0.faceIndicator != f1.faceIndicator) return f0.faceIndicator < f1.faceIndicator; - - /* Both faces belong to the same boundary marker. The second comparison is - based on the on the local volume ID's of the adjacent elements. As the - volumes are sorted according to the time levels for time accurate - local stepping, there is no need to do this check seperately here. */ - unsigned long ind0 = f0.elemID0 < nVolElemTot ? f0.elemID0 : f0.elemID1; - unsigned long ind1 = f1.elemID0 < nVolElemTot ? f1.elemID0 : f1.elemID1; - - return ind0 < ind1; - } - - /*--- Comparison in case both faces are internal faces. ---*/ - if (f0.faceIndicator == -1 && f1.faceIndicator == -1) { - /* Both faces are internal faces. First determine the minimum and maximum - ID of its adjacent elements. */ - unsigned long elemIDMin0 = min(f0.elemID0, f0.elemID1); - unsigned long elemIDMax0 = max(f0.elemID0, f0.elemID1); - - unsigned long elemIDMin1 = min(f1.elemID0, f1.elemID1); - unsigned long elemIDMax1 = max(f1.elemID0, f1.elemID1); - - /* Determine the situation. */ - if (elemIDMax0 < nVolElemTot && elemIDMax1 < nVolElemTot) { - /* Both faces are matching internal faces. Determine whether or not these - faces are local faces, i.e. faces between locally owned elements. */ - const bool face0IsLocal = elemIDMax0 < nVolElemOwned; - const bool face1IsLocal = elemIDMax1 < nVolElemOwned; - - /* Check if both faces have the same status, i.e. either local or - not local. */ - if (face0IsLocal == face1IsLocal) { - /* Both faces are either local or not local. Determine the time level - of the faces, which is the minimum value of the adjacent volume - elements. */ - const unsigned short timeLevel0 = min(volElem[elemIDMin0].timeLevel, volElem[elemIDMax0].timeLevel); - const unsigned short timeLevel1 = min(volElem[elemIDMin1].timeLevel, volElem[elemIDMax1].timeLevel); - - /* Internal faces with the same status are first sorted according to - their time level. Faces with the smallest time level are numbered - first. Note this is only relevant for time accurate local time - stepping. */ - if (timeLevel0 != timeLevel1) return timeLevel0 < timeLevel1; - - /* The faces belong to the same time level. They are sorted according - to their element ID's in order to increase cache performance. */ - if (elemIDMin0 != elemIDMin1) return elemIDMin0 < elemIDMin1; - return elemIDMax0 < elemIDMax1; - } /* One face is a local face and the other is not. Make sure that - the local faces are numbered first. */ - return face0IsLocal; - - } else if (elemIDMax0 >= nVolElemTot && elemIDMax1 >= nVolElemTot) { - /* Both faces are non-matching internal faces. Sort them according to - their relevant element ID. The time level is not taken into account - yet, because non-matching faces are not possible at the moment with - time accurate local time stepping. */ - return elemIDMin0 < elemIDMin1; - } else { - /* One face is a matching internal face and the other face is a - non-matching internal face. Make sure that the non-matching face - is numbered after the matching face. This is accomplished by comparing - the maximum element ID's. */ - return elemIDMax0 < elemIDMax1; - } - } - - /*--- One face is a boundary face and the other face is an internal face. - Make sure that the boundary face is numbered first. This can be - accomplished by using the > operator for faceIndicator. ---*/ - return f0.faceIndicator > f1.faceIndicator; -} - -bool CSortBoundaryFaces::operator()(const CSurfaceElementFEM& f0, const CSurfaceElementFEM& f1) { - /* First sorting criterion is the index of the standard element. The - boundary faces should be sorted per standard element. Note that the - time level is not taken into account here, because it is assumed that - the surface elements to be sorted belong to one time level. */ - if (f0.indStandardElement != f1.indStandardElement) return f0.indStandardElement < f1.indStandardElement; - - /* The standard elements are the same. The second criterion is the - corresponding volume IDs of the surface elements. */ - return f0.volElemID < f1.volElemID; -} - bool CPointFEM::operator<(const CPointFEM& other) const { if (periodIndexToDonor != other.periodIndexToDonor) return periodIndexToDonor < other.periodIndexToDonor; return globalID < other.globalID; diff --git a/Common/src/fem/meson.build b/Common/src/fem/meson.build index 263fccc0e84..414483d5108 100644 --- a/Common/src/fem/meson.build +++ b/Common/src/fem/meson.build @@ -1,5 +1,4 @@ -common_src += files(['geometry_structure_fem_part.cpp', - 'fem_geometry_structure.cpp', +common_src += files(['fem_geometry_structure.cpp', 'fem_integration_rules.cpp', 'fem_work_estimate_metis.cpp', 'fem_standard_element.cpp', diff --git a/Common/src/fem/geometry_structure_fem_part.cpp b/Common/src/geometry/CPhysicalGeometryFEM.cpp similarity index 92% rename from Common/src/fem/geometry_structure_fem_part.cpp rename to Common/src/geometry/CPhysicalGeometryFEM.cpp index f4734eed85d..e6a30c34fb2 100644 --- a/Common/src/fem/geometry_structure_fem_part.cpp +++ b/Common/src/geometry/CPhysicalGeometryFEM.cpp @@ -1,6 +1,6 @@ /*! - * \file geometry_structure_fem_part.cpp - * \brief Main subroutines for distributin the grid for the Fluid FEM solver. + * \file CPhysicalGeometryFEM.cpp + * \brief Implementation of the FEM related functions of CPhysicalGeometry. * \author F. Palacios, T. Economon * \version 8.1.0 "Harrier" * @@ -26,6 +26,9 @@ */ #include "../../include/toolboxes/CLinearPartitioner.hpp" +#include "../../include/toolboxes/classes_multiple_integers.hpp" +#include "../../include/toolboxes/fem/CFaceOfElement.hpp" +#include "../../include/toolboxes/fem/CMatchingFace.hpp" #include "../../include/geometry/CPhysicalGeometry.hpp" #include "../../include/fem/fem_standard_element.hpp" #include "../../include/geometry/primal_grid/CPrimalGridFEM.hpp" @@ -37,329 +40,6 @@ #include #include #include -/*--- Epsilon definition ---*/ - -CFaceOfElement::CFaceOfElement() { - nCornerPoints = 0; - cornerPoints[0] = cornerPoints[1] = cornerPoints[2] = cornerPoints[3] = ULONG_MAX; - elemID0 = elemID1 = ULONG_MAX; - nPolyGrid0 = nPolyGrid1 = 0; - nPolySol0 = nPolySol1 = 0; - nDOFsElem0 = nDOFsElem1 = 0; - elemType0 = elemType1 = 0; - faceID0 = faceID1 = 0; - periodicIndex = periodicIndexDonor = 0; - faceIndicator = 0; - - JacFaceIsConsideredConstant = false; - elem0IsOwner = false; -} - -CFaceOfElement::CFaceOfElement(const unsigned short VTK_Type, const unsigned short nPoly, const unsigned long* Nodes) { - /* Set the default values of the member variables. */ - nCornerPoints = 0; - cornerPoints[0] = cornerPoints[1] = cornerPoints[2] = cornerPoints[3] = ULONG_MAX; - elemID0 = elemID1 = ULONG_MAX; - nPolyGrid0 = nPolyGrid1 = 0; - nPolySol0 = nPolySol1 = 0; - nDOFsElem0 = nDOFsElem1 = 0; - elemType0 = elemType1 = 0; - faceID0 = faceID1 = 0; - periodicIndex = periodicIndexDonor = 0; - faceIndicator = 0; - - JacFaceIsConsideredConstant = false; - elem0IsOwner = false; - - /* Determine the face element type and set the corner points accordingly. */ - switch (VTK_Type) { - case LINE: { - nCornerPoints = 2; - cornerPoints[0] = Nodes[0]; - cornerPoints[1] = Nodes[nPoly]; - break; - } - - case TRIANGLE: { - const unsigned short ind2 = (nPoly + 1) * (nPoly + 2) / 2 - 1; - nCornerPoints = 3; - cornerPoints[0] = Nodes[0]; - cornerPoints[1] = Nodes[nPoly]; - cornerPoints[2] = Nodes[ind2]; - break; - } - - case QUADRILATERAL: { - const unsigned short ind2 = nPoly * (nPoly + 1); - const unsigned short ind3 = (nPoly + 1) * (nPoly + 1) - 1; - nCornerPoints = 4; - cornerPoints[0] = Nodes[0]; - cornerPoints[1] = Nodes[nPoly]; - cornerPoints[2] = Nodes[ind2]; - cornerPoints[3] = Nodes[ind3]; - break; - } - - default: { - ostringstream message; - message << "Unknown VTK surface element type, " << VTK_Type; - SU2_MPI::Error(message.str(), CURRENT_FUNCTION); - } - } -} - -bool CFaceOfElement::operator<(const CFaceOfElement& other) const { - if (nCornerPoints != other.nCornerPoints) return nCornerPoints < other.nCornerPoints; - - for (unsigned short i = 0; i < nCornerPoints; ++i) - if (cornerPoints[i] != other.cornerPoints[i]) return cornerPoints[i] < other.cornerPoints[i]; - - return false; -} - -bool CFaceOfElement::operator==(const CFaceOfElement& other) const { - if (nCornerPoints != other.nCornerPoints) return false; - for (unsigned short i = 0; i < nCornerPoints; ++i) - if (cornerPoints[i] != other.cornerPoints[i]) return false; - - return true; -} - -void CFaceOfElement::Copy(const CFaceOfElement& other) { - nCornerPoints = other.nCornerPoints; - for (unsigned short i = 0; i < nCornerPoints; ++i) cornerPoints[i] = other.cornerPoints[i]; - for (unsigned short i = nCornerPoints; i < 4; ++i) cornerPoints[i] = ULONG_MAX; - - elemID0 = other.elemID0; - elemID1 = other.elemID1; - - nPolyGrid0 = other.nPolyGrid0; - nPolyGrid1 = other.nPolyGrid1; - nPolySol0 = other.nPolySol0; - nPolySol1 = other.nPolySol1; - - nDOFsElem0 = other.nDOFsElem0; - nDOFsElem1 = other.nDOFsElem1; - elemType0 = other.elemType0; - elemType1 = other.elemType1; - faceID0 = other.faceID0; - faceID1 = other.faceID1; - - periodicIndex = other.periodicIndex; - periodicIndexDonor = other.periodicIndexDonor; - faceIndicator = other.faceIndicator; - - JacFaceIsConsideredConstant = other.JacFaceIsConsideredConstant; - elem0IsOwner = other.elem0IsOwner; -} - -void CFaceOfElement::CreateUniqueNumberingWithOrientation() { - /*--- Determine the element type and create the unique numbering accordingly. ---*/ - bool swapElements = false; - switch (nCornerPoints) { - case 2: { - /* Element is a line. Check if the node numbering must be swapped. If so - also the element information must be swapped, because element 0 is to - the left of the face and element 1 to the right. */ - if (cornerPoints[1] < cornerPoints[0]) { - swap(cornerPoints[0], cornerPoints[1]); - swapElements = true; - } - break; - } - - case 3: { - /* Element is a triangle. The vertices are sorted in increasing order. - If the sequence of the new numbering is opposite to the current - numbering, the element information must be exchanged, because - element 0 is to the left of the face and element 1 to the right. */ - unsigned long nn[] = {cornerPoints[0], cornerPoints[1], cornerPoints[2]}; - unsigned short ind = 0; - if (nn[1] < nn[ind]) ind = 1; - if (nn[2] < nn[ind]) ind = 2; - - unsigned short indm1 = ind == 0 ? 2 : ind - 1; // Next lower index. - unsigned short indp1 = ind == 2 ? 0 : ind + 1; // Next upper index. - - if (nn[indp1] < nn[indm1]) { - /* The orientation of the triangle remains the same. - Store the new sorted node numbering. */ - cornerPoints[0] = nn[ind]; - cornerPoints[1] = nn[indp1]; - cornerPoints[2] = nn[indm1]; - } else { - /* The orientation of the triangle changes. Store the new - sorted node numbering and set swapElements to true. */ - cornerPoints[0] = nn[ind]; - cornerPoints[1] = nn[indm1]; - cornerPoints[2] = nn[indp1]; - swapElements = true; - } - - break; - } - - case 4: { - /* Element is a quadrilateral. The vertices are sorted in increasing order - under the condition neighboring vertices remain neighbors. If the - sequence of the new numbering is opposite to the current - numbering, the element information must be exchanged, because - element 0 is to the left of the face and element 1 to the right. */ - unsigned long nn[] = {cornerPoints[0], cornerPoints[1], cornerPoints[2], cornerPoints[3]}; - unsigned short ind = 0; - if (nn[1] < nn[ind]) ind = 1; - if (nn[2] < nn[ind]) ind = 2; - if (nn[3] < nn[ind]) ind = 3; - - unsigned short indm1 = ind == 0 ? 3 : ind - 1; // Next lower index. - unsigned short indp1 = ind == 3 ? 0 : ind + 1; // Next upper index. - unsigned short indp2 = ind >= 2 ? ind - 2 : ind + 2; // Opposite index. - - if (nn[indp1] < nn[indm1]) { - /* The orientation of the quadrilateral remains the same. - Store the new sorted node numbering. */ - cornerPoints[0] = nn[ind]; - cornerPoints[1] = nn[indp1]; - cornerPoints[2] = nn[indp2]; - cornerPoints[3] = nn[indm1]; - } else { - /* The orientation of the quadrilateral changes. Store the new - sorted node numbering and set swapElements to true. */ - cornerPoints[0] = nn[ind]; - cornerPoints[1] = nn[indm1]; - cornerPoints[2] = nn[indp2]; - cornerPoints[3] = nn[indp1]; - swapElements = true; - } - - break; - } - - default: { - ostringstream message; - message << "Unknown surface element type with " << nCornerPoints << " corners." << endl; - SU2_MPI::Error(message.str(), CURRENT_FUNCTION); - } - } - - /* Swap the element information, if needed. */ - if (swapElements) { - swap(elemID0, elemID1); - swap(nPolyGrid0, nPolyGrid1); - swap(nPolySol0, nPolySol1); - swap(nDOFsElem0, nDOFsElem1); - swap(elemType0, elemType1); - swap(faceID0, faceID1); - } -} - -void CBoundaryFace::Copy(const CBoundaryFace& other) { - VTK_Type = other.VTK_Type; - nPolyGrid = other.nPolyGrid; - nDOFsGrid = other.nDOFsGrid; - globalBoundElemID = other.globalBoundElemID; - domainElementID = other.domainElementID; - Nodes = other.Nodes; -} - -CMatchingFace::CMatchingFace() { - nCornerPoints = 0; - nDim = 0; - nPoly = 0; - nDOFsElem = 0; - elemType = 0; - elemID = 0; - - cornerCoor[0][0] = cornerCoor[0][1] = cornerCoor[0][2] = 0.0; - cornerCoor[1][0] = cornerCoor[1][1] = cornerCoor[1][2] = 0.0; - cornerCoor[2][0] = cornerCoor[2][1] = cornerCoor[2][2] = 0.0; - cornerCoor[3][0] = cornerCoor[3][1] = cornerCoor[3][2] = 0.0; - - tolForMatching = 0.0; -} - -bool CMatchingFace::operator<(const CMatchingFace& other) const { - /* First compare the number of corner points. ---*/ - if (nCornerPoints != other.nCornerPoints) return nCornerPoints < other.nCornerPoints; - - /*--- Determine the tolerance for comparing both objects. ---*/ - const su2double tol = min(tolForMatching, other.tolForMatching); - - /*--- Loop over the number of corner points and dimensions and compare the - coordinates. If considered different, return true if the current face - is considered smaller and false otherwise. ---*/ - for (unsigned short k = 0; k < nCornerPoints; ++k) { - for (unsigned short l = 0; l < nDim; ++l) { - if (fabs(cornerCoor[k][l] - other.cornerCoor[k][l]) > tol) return cornerCoor[k][l] < other.cornerCoor[k][l]; - } - } - - /*--- Both objects are considered the same. Return false. ---*/ - return false; -} - -void CMatchingFace::SortFaceCoordinates() { - /*--- Determine the tolerance for a matching point for this face. This is - accomplished by computing the minimum distance between the points of - the face, multiplied by a relative tolerance. ---*/ - for (unsigned short k = 0; k < nCornerPoints; ++k) { - for (unsigned short j = (k + 1); j < nCornerPoints; ++j) { - su2double dist = 0.0; - for (unsigned short l = 0; l < nDim; ++l) { - su2double ds = cornerCoor[k][l] - cornerCoor[j][l]; - dist += ds * ds; - } - dist = sqrt(dist); - - if (k == 0 && j == 1) - tolForMatching = dist; - else - tolForMatching = min(tolForMatching, dist); - } - } - - tolForMatching *= 1.e-6; - - /*--- Sort the points in increasing order based on the coordinates. - An insertion sort algorithm is used, which is quite efficient - for at most four corner points. ---*/ - for (unsigned short k = 1; k < nCornerPoints; ++k) { - for (unsigned short j = k; j > 0; --j) { - /* Check if cornerCoor[j] is considered less than cornerCoor[j-1]. */ - bool lessThan = false; - for (unsigned short l = 0; l < nDim; ++l) { - if (fabs(cornerCoor[j][l] - cornerCoor[j - 1][l]) > tolForMatching) { - lessThan = cornerCoor[j][l] < cornerCoor[j - 1][l]; - break; - } - } - - /* If cornerCoor[j] is less than cornerCoor[j-1] they must be swapped. - Otherwise an exit can be made from the j-loop. */ - if (lessThan) { - for (unsigned short l = 0; l < nDim; ++l) swap(cornerCoor[j][l], cornerCoor[j - 1][l]); - } else - break; - } - } -} - -void CMatchingFace::Copy(const CMatchingFace& other) { - nCornerPoints = other.nCornerPoints; - nDim = other.nDim; - nPoly = other.nPoly; - nDOFsElem = other.nDOFsElem; - elemType = other.elemType; - elemID = other.elemID; - - for (unsigned short k = 0; k < nCornerPoints; ++k) { - for (unsigned l = 0; l < nDim; ++l) { - cornerCoor[k][l] = other.cornerCoor[k][l]; - } - } - - tolForMatching = other.tolForMatching; -} void CPhysicalGeometry::LoadLinearlyPartitionedPointsFEM(CConfig* config, CMeshReaderBase* mesh) { /*--- Get the partitioned coordinates and their global IDs from the mesh object. ---*/ diff --git a/Common/src/geometry/meson.build b/Common/src/geometry/meson.build index b67e2b98658..aeaabdb1a26 100644 --- a/Common/src/geometry/meson.build +++ b/Common/src/geometry/meson.build @@ -1,5 +1,6 @@ common_src += files(['CGeometry.cpp', 'CPhysicalGeometry.cpp', + 'CPhysicalGeometryFEM.cpp', 'CMultiGridGeometry.cpp', 'CDummyGeometry.cpp', 'CMultiGridQueue.cpp']) diff --git a/Common/src/toolboxes/fem/CFaceOfElement.cpp b/Common/src/toolboxes/fem/CFaceOfElement.cpp new file mode 100644 index 00000000000..cb46aad39a6 --- /dev/null +++ b/Common/src/toolboxes/fem/CFaceOfElement.cpp @@ -0,0 +1,243 @@ +/*! + * \file CFaceOfElement.cpp + * \brief Helper class used in distributing the surface elements and + * creating the surface elements for the DG solver. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#include "../../../include/toolboxes/fem/CFaceOfElement.hpp" + +CFaceOfElement::CFaceOfElement() { + nCornerPoints = 0; + cornerPoints[0] = cornerPoints[1] = cornerPoints[2] = cornerPoints[3] = ULONG_MAX; + elemID0 = elemID1 = ULONG_MAX; + nPolyGrid0 = nPolyGrid1 = 0; + nPolySol0 = nPolySol1 = 0; + nDOFsElem0 = nDOFsElem1 = 0; + elemType0 = elemType1 = 0; + faceID0 = faceID1 = 0; + periodicIndex = periodicIndexDonor = 0; + faceIndicator = 0; + + JacFaceIsConsideredConstant = false; + elem0IsOwner = false; +} + +CFaceOfElement::CFaceOfElement(const unsigned short VTK_Type, const unsigned short nPoly, const unsigned long* Nodes) { + /* Set the default values of the member variables. */ + nCornerPoints = 0; + cornerPoints[0] = cornerPoints[1] = cornerPoints[2] = cornerPoints[3] = ULONG_MAX; + elemID0 = elemID1 = ULONG_MAX; + nPolyGrid0 = nPolyGrid1 = 0; + nPolySol0 = nPolySol1 = 0; + nDOFsElem0 = nDOFsElem1 = 0; + elemType0 = elemType1 = 0; + faceID0 = faceID1 = 0; + periodicIndex = periodicIndexDonor = 0; + faceIndicator = 0; + + JacFaceIsConsideredConstant = false; + elem0IsOwner = false; + + /* Determine the face element type and set the corner points accordingly. */ + switch (VTK_Type) { + case LINE: { + nCornerPoints = 2; + cornerPoints[0] = Nodes[0]; + cornerPoints[1] = Nodes[nPoly]; + break; + } + + case TRIANGLE: { + const unsigned short ind2 = (nPoly + 1) * (nPoly + 2) / 2 - 1; + nCornerPoints = 3; + cornerPoints[0] = Nodes[0]; + cornerPoints[1] = Nodes[nPoly]; + cornerPoints[2] = Nodes[ind2]; + break; + } + + case QUADRILATERAL: { + const unsigned short ind2 = nPoly * (nPoly + 1); + const unsigned short ind3 = (nPoly + 1) * (nPoly + 1) - 1; + nCornerPoints = 4; + cornerPoints[0] = Nodes[0]; + cornerPoints[1] = Nodes[nPoly]; + cornerPoints[2] = Nodes[ind2]; + cornerPoints[3] = Nodes[ind3]; + break; + } + + default: { + ostringstream message; + message << "Unknown VTK surface element type, " << VTK_Type; + SU2_MPI::Error(message.str(), CURRENT_FUNCTION); + } + } +} + +bool CFaceOfElement::operator<(const CFaceOfElement& other) const { + if (nCornerPoints != other.nCornerPoints) return nCornerPoints < other.nCornerPoints; + + for (unsigned short i = 0; i < nCornerPoints; ++i) + if (cornerPoints[i] != other.cornerPoints[i]) return cornerPoints[i] < other.cornerPoints[i]; + + return false; +} + +bool CFaceOfElement::operator==(const CFaceOfElement& other) const { + if (nCornerPoints != other.nCornerPoints) return false; + for (unsigned short i = 0; i < nCornerPoints; ++i) + if (cornerPoints[i] != other.cornerPoints[i]) return false; + + return true; +} + +void CFaceOfElement::Copy(const CFaceOfElement& other) { + nCornerPoints = other.nCornerPoints; + for (unsigned short i = 0; i < nCornerPoints; ++i) cornerPoints[i] = other.cornerPoints[i]; + for (unsigned short i = nCornerPoints; i < 4; ++i) cornerPoints[i] = ULONG_MAX; + + elemID0 = other.elemID0; + elemID1 = other.elemID1; + + nPolyGrid0 = other.nPolyGrid0; + nPolyGrid1 = other.nPolyGrid1; + nPolySol0 = other.nPolySol0; + nPolySol1 = other.nPolySol1; + + nDOFsElem0 = other.nDOFsElem0; + nDOFsElem1 = other.nDOFsElem1; + elemType0 = other.elemType0; + elemType1 = other.elemType1; + faceID0 = other.faceID0; + faceID1 = other.faceID1; + + periodicIndex = other.periodicIndex; + periodicIndexDonor = other.periodicIndexDonor; + faceIndicator = other.faceIndicator; + + JacFaceIsConsideredConstant = other.JacFaceIsConsideredConstant; + elem0IsOwner = other.elem0IsOwner; +} + +void CFaceOfElement::CreateUniqueNumberingWithOrientation() { + /*--- Determine the element type and create the unique numbering accordingly. ---*/ + bool swapElements = false; + switch (nCornerPoints) { + case 2: { + /* Element is a line. Check if the node numbering must be swapped. If so + also the element information must be swapped, because element 0 is to + the left of the face and element 1 to the right. */ + if (cornerPoints[1] < cornerPoints[0]) { + swap(cornerPoints[0], cornerPoints[1]); + swapElements = true; + } + break; + } + + case 3: { + /* Element is a triangle. The vertices are sorted in increasing order. + If the sequence of the new numbering is opposite to the current + numbering, the element information must be exchanged, because + element 0 is to the left of the face and element 1 to the right. */ + unsigned long nn[] = {cornerPoints[0], cornerPoints[1], cornerPoints[2]}; + unsigned short ind = 0; + if (nn[1] < nn[ind]) ind = 1; + if (nn[2] < nn[ind]) ind = 2; + + unsigned short indm1 = ind == 0 ? 2 : ind - 1; // Next lower index. + unsigned short indp1 = ind == 2 ? 0 : ind + 1; // Next upper index. + + if (nn[indp1] < nn[indm1]) { + /* The orientation of the triangle remains the same. + Store the new sorted node numbering. */ + cornerPoints[0] = nn[ind]; + cornerPoints[1] = nn[indp1]; + cornerPoints[2] = nn[indm1]; + } else { + /* The orientation of the triangle changes. Store the new + sorted node numbering and set swapElements to true. */ + cornerPoints[0] = nn[ind]; + cornerPoints[1] = nn[indm1]; + cornerPoints[2] = nn[indp1]; + swapElements = true; + } + + break; + } + + case 4: { + /* Element is a quadrilateral. The vertices are sorted in increasing order + under the condition neighboring vertices remain neighbors. If the + sequence of the new numbering is opposite to the current + numbering, the element information must be exchanged, because + element 0 is to the left of the face and element 1 to the right. */ + unsigned long nn[] = {cornerPoints[0], cornerPoints[1], cornerPoints[2], cornerPoints[3]}; + unsigned short ind = 0; + if (nn[1] < nn[ind]) ind = 1; + if (nn[2] < nn[ind]) ind = 2; + if (nn[3] < nn[ind]) ind = 3; + + unsigned short indm1 = ind == 0 ? 3 : ind - 1; // Next lower index. + unsigned short indp1 = ind == 3 ? 0 : ind + 1; // Next upper index. + unsigned short indp2 = ind >= 2 ? ind - 2 : ind + 2; // Opposite index. + + if (nn[indp1] < nn[indm1]) { + /* The orientation of the quadrilateral remains the same. + Store the new sorted node numbering. */ + cornerPoints[0] = nn[ind]; + cornerPoints[1] = nn[indp1]; + cornerPoints[2] = nn[indp2]; + cornerPoints[3] = nn[indm1]; + } else { + /* The orientation of the quadrilateral changes. Store the new + sorted node numbering and set swapElements to true. */ + cornerPoints[0] = nn[ind]; + cornerPoints[1] = nn[indm1]; + cornerPoints[2] = nn[indp2]; + cornerPoints[3] = nn[indp1]; + swapElements = true; + } + + break; + } + + default: { + ostringstream message; + message << "Unknown surface element type with " << nCornerPoints << " corners." << endl; + SU2_MPI::Error(message.str(), CURRENT_FUNCTION); + } + } + + /* Swap the element information, if needed. */ + if (swapElements) { + swap(elemID0, elemID1); + swap(nPolyGrid0, nPolyGrid1); + swap(nPolySol0, nPolySol1); + swap(nDOFsElem0, nDOFsElem1); + swap(elemType0, elemType1); + swap(faceID0, faceID1); + } +} diff --git a/Common/src/toolboxes/fem/CMatchingFace.cpp b/Common/src/toolboxes/fem/CMatchingFace.cpp new file mode 100644 index 00000000000..3f01b67e8a3 --- /dev/null +++ b/Common/src/toolboxes/fem/CMatchingFace.cpp @@ -0,0 +1,128 @@ +/*! + * \file CMatchingFace.cpp + * \brief Helper class used to determine whether or not faces of + * surface elements for the DG solver are matching faces. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#include "../../../include/toolboxes/fem/CMatchingFace.hpp" + +CMatchingFace::CMatchingFace() { + nCornerPoints = 0; + nDim = 0; + nPoly = 0; + nDOFsElem = 0; + elemType = 0; + elemID = 0; + + cornerCoor[0][0] = cornerCoor[0][1] = cornerCoor[0][2] = 0.0; + cornerCoor[1][0] = cornerCoor[1][1] = cornerCoor[1][2] = 0.0; + cornerCoor[2][0] = cornerCoor[2][1] = cornerCoor[2][2] = 0.0; + cornerCoor[3][0] = cornerCoor[3][1] = cornerCoor[3][2] = 0.0; + + tolForMatching = 0.0; +} + +bool CMatchingFace::operator<(const CMatchingFace& other) const { + /* First compare the number of corner points. ---*/ + if (nCornerPoints != other.nCornerPoints) return nCornerPoints < other.nCornerPoints; + + /*--- Determine the tolerance for comparing both objects. ---*/ + const su2double tol = min(tolForMatching, other.tolForMatching); + + /*--- Loop over the number of corner points and dimensions and compare the + coordinates. If considered different, return true if the current face + is considered smaller and false otherwise. ---*/ + for (unsigned short k = 0; k < nCornerPoints; ++k) { + for (unsigned short l = 0; l < nDim; ++l) { + if (fabs(cornerCoor[k][l] - other.cornerCoor[k][l]) > tol) return cornerCoor[k][l] < other.cornerCoor[k][l]; + } + } + + /*--- Both objects are considered the same. Return false. ---*/ + return false; +} + +void CMatchingFace::SortFaceCoordinates() { + /*--- Determine the tolerance for a matching point for this face. This is + accomplished by computing the minimum distance between the points of + the face, multiplied by a relative tolerance. ---*/ + for (unsigned short k = 0; k < nCornerPoints; ++k) { + for (unsigned short j = (k + 1); j < nCornerPoints; ++j) { + su2double dist = 0.0; + for (unsigned short l = 0; l < nDim; ++l) { + su2double ds = cornerCoor[k][l] - cornerCoor[j][l]; + dist += ds * ds; + } + dist = sqrt(dist); + + if (k == 0 && j == 1) + tolForMatching = dist; + else + tolForMatching = min(tolForMatching, dist); + } + } + + tolForMatching *= 1.e-6; + + /*--- Sort the points in increasing order based on the coordinates. + An insertion sort algorithm is used, which is quite efficient + for at most four corner points. ---*/ + for (unsigned short k = 1; k < nCornerPoints; ++k) { + for (unsigned short j = k; j > 0; --j) { + /* Check if cornerCoor[j] is considered less than cornerCoor[j-1]. */ + bool lessThan = false; + for (unsigned short l = 0; l < nDim; ++l) { + if (fabs(cornerCoor[j][l] - cornerCoor[j - 1][l]) > tolForMatching) { + lessThan = cornerCoor[j][l] < cornerCoor[j - 1][l]; + break; + } + } + + /* If cornerCoor[j] is less than cornerCoor[j-1] they must be swapped. + Otherwise an exit can be made from the j-loop. */ + if (lessThan) { + for (unsigned short l = 0; l < nDim; ++l) swap(cornerCoor[j][l], cornerCoor[j - 1][l]); + } else + break; + } + } +} + +void CMatchingFace::Copy(const CMatchingFace& other) { + nCornerPoints = other.nCornerPoints; + nDim = other.nDim; + nPoly = other.nPoly; + nDOFsElem = other.nDOFsElem; + elemType = other.elemType; + elemID = other.elemID; + + for (unsigned short k = 0; k < nCornerPoints; ++k) { + for (unsigned l = 0; l < nDim; ++l) { + cornerCoor[k][l] = other.cornerCoor[k][l]; + } + } + + tolForMatching = other.tolForMatching; +} diff --git a/Common/src/toolboxes/fem/CReorderElements.cpp b/Common/src/toolboxes/fem/CReorderElements.cpp new file mode 100644 index 00000000000..a8152170652 --- /dev/null +++ b/Common/src/toolboxes/fem/CReorderElements.cpp @@ -0,0 +1,65 @@ +/*! + * \file CReorderElements.cpp + * \brief Helper class used to reorder the owned elements + * after the partitioning + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#include "../../../include/toolboxes/fem/CReorderElements.hpp" + +CReorderElements::CReorderElements(const unsigned long val_GlobalElemID, const unsigned short val_TimeLevel, + const bool val_CommSolution, const unsigned short val_VTK_Type, + const unsigned short val_nPolySol, const bool val_JacConstant) { + /* Copy the global elment ID, time level and whether or not this element + must be communicated. */ + globalElemID = val_GlobalElemID; + timeLevel = val_TimeLevel; + commSolution = val_CommSolution; + + /* Create the element type used in this class, which stores information of + the VTK type, the polynomial degree of the solution and whether or not the + Jacobian of the transformation is constant. As it is possible that the + polynomial degree of the solution is zero, this convention is different + from the convention used in the SU2 grid file. */ + elemType = val_VTK_Type + 100 * val_nPolySol; + if (!val_JacConstant) elemType += 50; +} + +bool CReorderElements::operator<(const CReorderElements& other) const { + /* Elements with the lowest time level are stored first. */ + if (timeLevel != other.timeLevel) return timeLevel < other.timeLevel; + + /* Next comparison is whether or not the element must communicate its + solution data to other ranks. Elements which do not need to do this + are stored first. */ + if (commSolution != other.commSolution) return other.commSolution; + + /* Elements of the same element type must be stored as contiguously as + possible to allow for the simultaneous treatment of elements in the + matrix multiplications. */ + if (elemType != other.elemType) return elemType < other.elemType; + + /* The final comparison is based on the global element ID. */ + return globalElemID < other.globalElemID; +} diff --git a/Common/src/toolboxes/fem/CSortBoundaryFaces.cpp b/Common/src/toolboxes/fem/CSortBoundaryFaces.cpp new file mode 100644 index 00000000000..945b110a011 --- /dev/null +++ b/Common/src/toolboxes/fem/CSortBoundaryFaces.cpp @@ -0,0 +1,42 @@ +/*! + * \file CSortFaces.cpp + * \brief Functor, used for a different sorting of the surface elements + * than the < operator of CSurfaceElementFEM. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#include "../../../include/toolboxes/fem/CSortBoundaryFaces.hpp" +#include "../../../include/fem/fem_geometry_structure.hpp" + +bool CSortBoundaryFaces::operator()(const CSurfaceElementFEM& f0, const CSurfaceElementFEM& f1) { + /* First sorting criterion is the index of the standard element. The + boundary faces should be sorted per standard element. Note that the + time level is not taken into account here, because it is assumed that + the surface elements to be sorted belong to one time level. */ + if (f0.indStandardElement != f1.indStandardElement) return f0.indStandardElement < f1.indStandardElement; + + /* The standard elements are the same. The second criterion is the + corresponding volume IDs of the surface elements. */ + return f0.volElemID < f1.volElemID; +} diff --git a/Common/src/toolboxes/fem/CSortFaces.cpp b/Common/src/toolboxes/fem/CSortFaces.cpp new file mode 100644 index 00000000000..68f9c8fbbc6 --- /dev/null +++ b/Common/src/toolboxes/fem/CSortFaces.cpp @@ -0,0 +1,108 @@ +/*! + * \file CSortFaces.cpp + * \brief Functor, used for a different sorting of the surface elements + * than the < operator of CFaceOfElement. + * \author E. van der Weide + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 . + */ + +#include "../../../include/toolboxes/fem/CSortFaces.hpp" +#include "../../../include/fem/fem_geometry_structure.hpp" + +bool CSortFaces::operator()(const CFaceOfElement& f0, const CFaceOfElement& f1) { + /*--- Comparison in case both faces are boundary faces. ---*/ + if (f0.faceIndicator >= 0 && f1.faceIndicator >= 0) { + /* Both faces are boundary faces. The first comparison is the boundary + marker, which is stored in faceIndicator. */ + if (f0.faceIndicator != f1.faceIndicator) return f0.faceIndicator < f1.faceIndicator; + + /* Both faces belong to the same boundary marker. The second comparison is + based on the on the local volume ID's of the adjacent elements. As the + volumes are sorted according to the time levels for time accurate + local stepping, there is no need to do this check seperately here. */ + unsigned long ind0 = f0.elemID0 < nVolElemTot ? f0.elemID0 : f0.elemID1; + unsigned long ind1 = f1.elemID0 < nVolElemTot ? f1.elemID0 : f1.elemID1; + + return ind0 < ind1; + } + + /*--- Comparison in case both faces are internal faces. ---*/ + if (f0.faceIndicator == -1 && f1.faceIndicator == -1) { + /* Both faces are internal faces. First determine the minimum and maximum + ID of its adjacent elements. */ + unsigned long elemIDMin0 = min(f0.elemID0, f0.elemID1); + unsigned long elemIDMax0 = max(f0.elemID0, f0.elemID1); + + unsigned long elemIDMin1 = min(f1.elemID0, f1.elemID1); + unsigned long elemIDMax1 = max(f1.elemID0, f1.elemID1); + + /* Determine the situation. */ + if (elemIDMax0 < nVolElemTot && elemIDMax1 < nVolElemTot) { + /* Both faces are matching internal faces. Determine whether or not these + faces are local faces, i.e. faces between locally owned elements. */ + const bool face0IsLocal = elemIDMax0 < nVolElemOwned; + const bool face1IsLocal = elemIDMax1 < nVolElemOwned; + + /* Check if both faces have the same status, i.e. either local or + not local. */ + if (face0IsLocal == face1IsLocal) { + /* Both faces are either local or not local. Determine the time level + of the faces, which is the minimum value of the adjacent volume + elements. */ + const unsigned short timeLevel0 = min(volElem[elemIDMin0].timeLevel, volElem[elemIDMax0].timeLevel); + const unsigned short timeLevel1 = min(volElem[elemIDMin1].timeLevel, volElem[elemIDMax1].timeLevel); + + /* Internal faces with the same status are first sorted according to + their time level. Faces with the smallest time level are numbered + first. Note this is only relevant for time accurate local time + stepping. */ + if (timeLevel0 != timeLevel1) return timeLevel0 < timeLevel1; + + /* The faces belong to the same time level. They are sorted according + to their element ID's in order to increase cache performance. */ + if (elemIDMin0 != elemIDMin1) return elemIDMin0 < elemIDMin1; + return elemIDMax0 < elemIDMax1; + } /* One face is a local face and the other is not. Make sure that + the local faces are numbered first. */ + return face0IsLocal; + + } else if (elemIDMax0 >= nVolElemTot && elemIDMax1 >= nVolElemTot) { + /* Both faces are non-matching internal faces. Sort them according to + their relevant element ID. The time level is not taken into account + yet, because non-matching faces are not possible at the moment with + time accurate local time stepping. */ + return elemIDMin0 < elemIDMin1; + } else { + /* One face is a matching internal face and the other face is a + non-matching internal face. Make sure that the non-matching face + is numbered after the matching face. This is accomplished by comparing + the maximum element ID's. */ + return elemIDMax0 < elemIDMax1; + } + } + + /*--- One face is a boundary face and the other face is an internal face. + Make sure that the boundary face is numbered first. This can be + accomplished by using the > operator for faceIndicator. ---*/ + return f0.faceIndicator > f1.faceIndicator; +} diff --git a/Common/src/toolboxes/fem/meson.build b/Common/src/toolboxes/fem/meson.build new file mode 100644 index 00000000000..a37f901fffb --- /dev/null +++ b/Common/src/toolboxes/fem/meson.build @@ -0,0 +1,5 @@ +common_src += files(['CFaceOfElement.cpp', + 'CMatchingFace.cpp', + 'CReorderElements.cpp', + 'CSortBoundaryFaces.cpp', + 'CSortFaces.cpp']) diff --git a/Common/src/toolboxes/meson.build b/Common/src/toolboxes/meson.build index d95ffad5992..fa5e3b234c3 100644 --- a/Common/src/toolboxes/meson.build +++ b/Common/src/toolboxes/meson.build @@ -5,3 +5,4 @@ common_src += files(['CLinearPartitioner.cpp', 'CSymmetricMatrix.cpp']) subdir('MMS') +subdir('fem') From 9c16769c20884a7eb155108a1fdd35a117ebf7ec Mon Sep 17 00:00:00 2001 From: vdweide Date: Mon, 10 Feb 2025 17:24:50 +0100 Subject: [PATCH 2/9] Split the function SetColorFEMGrid_Parallel in several smaller functions to improve readability --- Common/include/geometry/CPhysicalGeometry.hpp | 80 +- .../toolboxes/fem/CReorderElements.hpp | 1 - Common/src/geometry/CPhysicalGeometryFEM.cpp | 890 ++++++++++-------- 3 files changed, 550 insertions(+), 421 deletions(-) diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index 7d2f7d31f09..92dd306b34f 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/Common/include/geometry/CPhysicalGeometry.hpp @@ -499,6 +499,35 @@ class CPhysicalGeometry final : public CGeometry { */ void SetColorFEMGrid_Parallel(CConfig* config) override; + /*! + * \brief Determine the donor elements for the boundary elements on viscous + wall boundaries when wall functions are used. + * \param[in] config - Definition of the particular problem. + */ + void DetermineDonorElementsWallFunctions(CConfig* config); + + #ifdef HAVE_MPI + #ifdef HAVE_PARMETIS + /*! + * \brief Function, which converts the input the format for ParMETIS and calls + * ParMETIS to determine the actual colors of the elements. + * \param[in] adjacency - Adjacency information of the elements. + * \param[in] vwgt - Weights of the vertices of the graph, which are the elements. + * \param[in] adjwgt - Weights of the adjacencies of the graph. + */ + void DetermineFEMColorsViaParMETIS(vector > &adjacency, + vector &vwgt, + vector > &adjwgt); + #endif + #endif + + /*! + * \brief Determine whether or not the Jacobians of the elements and faces + are constant and a length scale of the elements. + * \param[in] config - Definition of the particular problem. + */ + void DetermineFEMConstantJacobiansAndLenScale(CConfig* config); + /*! * \brief Compute the weights of the FEM graph for ParMETIS. * \param[in] config - Definition of the particular problem. @@ -509,24 +538,44 @@ class CPhysicalGeometry final : public CGeometry { * \param[out] vwgt - Weights of the vertices of the graph, i.e. the elements. * \param[out] adjwgt - Weights of the edges of the graph. */ - void ComputeFEMGraphWeights(CConfig* config, const vector& localFaces, - const vector >& adjacency, - const map& mapExternalElemIDToTimeLevel, - vector& vwgt, vector >& adjwgt); + void DetermineFEMGraphWeights(CConfig* config, const vector& localFaces, + const vector >& adjacency, + const map& mapExternalElemIDToTimeLevel, + vector& vwgt, vector >& adjwgt); /*! - * \brief Determine the donor elements for the boundary elements on viscous - wall boundaries when wall functions are used. - * \param[in] config - Definition of the particular problem. + * \brief Function, which determines the adjacency information of the graph + * representation of the grid. + * \param[in] localFaces - Vector containing the local matching faces of the FEM grid. + * \param[out] adjacency - Vector of vectors to store the adjacency. */ - void DetermineDonorElementsWallFunctions(CConfig* config); + void DetermineGraphAdjacency(const vector& localFaces, + vector >& adjacency); /*! - * \brief Determine whether or not the Jacobians of the elements and faces - are constant and a length scale of the elements. - * \param[in] config - Definition of the particular problem. + * \brief Function, which determines the matching faces of a FEM grid. + * \param[in] config - Definition of the particular problem. + * \param[out] localFaces - Vector containing the local faces of the FEM grid. + * On output the matching faces are stored as one face. */ - void DetermineFEMConstantJacobiansAndLenScale(CConfig* config); + void DetermineMatchingFacesFEMGrid(const CConfig* config, vector& localFaces); + + /*! + * \brief Function, which determines the non-matching faces of a FEM grid. + * \param[in] config - Definition of the particular problem. + * \param[in,out] localMatchingFaces - Vector containing the local faces of the FEM grid. + * On output the non-matching faces are removed. + */ + void DetermineNonMatchingFacesFEMGrid(const CConfig* config, vector& localMatchingFaces); + + /*! + * \brief Function, which determines the owner of the internal faces, i.e. which element + * is responsible for computing the fluxes through the face. + * \param[in] localFaces - Vector, which contains the element faces of this rank. + * \param[out] mapExternalElemIDToTimeLevel - Map from the external element ID's to their time level and number of DOFs. + */ + void DetermineOwnershipInternalFaces(vector &localFaces, + map &mapExternalElemIDToTimeLevel); /*! * \brief Determine the neighboring information for periodic faces of a FEM grid. @@ -545,6 +594,13 @@ class CPhysicalGeometry final : public CGeometry { void DetermineTimeLevelElements(CConfig* config, const vector& localFaces, map& mapExternalElemIDToTimeLevel); + /*! + * \brief Function, which stores the information of the local matching faces in the + * data structures of the local elements. + * \param[in] localFaces - Vector, which contains the internal matching faces of this rank. + */ + void StoreFaceInfoInLocalElements(const vector &localFaces); + /*! * \brief Compute 3 grid quality metrics: orthogonality angle, dual cell aspect ratio, and dual cell volume ratio. * \param[in] config - Definition of the particular problem. diff --git a/Common/include/toolboxes/fem/CReorderElements.hpp b/Common/include/toolboxes/fem/CReorderElements.hpp index a7408aa3d37..ac2cd6c6c27 100644 --- a/Common/include/toolboxes/fem/CReorderElements.hpp +++ b/Common/include/toolboxes/fem/CReorderElements.hpp @@ -28,7 +28,6 @@ #pragma once - #include using namespace std; diff --git a/Common/src/geometry/CPhysicalGeometryFEM.cpp b/Common/src/geometry/CPhysicalGeometryFEM.cpp index e6a30c34fb2..ee2857a269b 100644 --- a/Common/src/geometry/CPhysicalGeometryFEM.cpp +++ b/Common/src/geometry/CPhysicalGeometryFEM.cpp @@ -167,11 +167,249 @@ void CPhysicalGeometry::LoadLinearlyPartitionedSurfaceElementsFEM(CConfig* confi } void CPhysicalGeometry::SetColorFEMGrid_Parallel(CConfig* config) { + /*--- Initialize the color vector of the elements. ---*/ for (unsigned long i = 0; i < nElem; ++i) elem[i]->SetColor(0); - /*--- Determine the faces of the elements. ---*/ + /*--- Determine the matching faces of the local elements. ---*/ vector localFaces; + DetermineMatchingFacesFEMGrid(config, localFaces); + + /*--- In the procedure above the periodic boundaries are not found. + A different treatment must be used in order to find these. ---*/ + DeterminePeriodicFacesFEMGrid(config, localFaces); + + /*--- So far only the matching faces have been found. Now find the + non-matching faces of the local elements. ---*/ + DetermineNonMatchingFacesFEMGrid(config, localFaces); + + /*--- Determine whether or not the Jacobians of the elements and faces + can be considered constant. ---*/ + DetermineFEMConstantJacobiansAndLenScale(config); + + /*--- Determine the donor elements for the wall function treatment. ---*/ + DetermineDonorElementsWallFunctions(config); + + /*--- Determine the time levels of the elements. This is only relevant + when time accurate local time stepping is employed. ---*/ + map mapExternalElemIDToTimeLevel; + DetermineTimeLevelElements(config, localFaces, mapExternalElemIDToTimeLevel); + + /*--- Determine the ownership of the internal faces, i.e. which adjacent + element is responsible for computing the fluxes through the face. ---*/ + DetermineOwnershipInternalFaces(localFaces, mapExternalElemIDToTimeLevel); + + /*--- All the matching face information is known now, including periodic + faces. Store the information of the neighbors in the data structure + for the local elements. ---*/ + StoreFaceInfoInLocalElements(localFaces); + + /*--- Create the vector of vectors that describe the connectivity + of the graph. ---*/ + vector > adjacency; + DetermineGraphAdjacency(localFaces, adjacency); + + /* Determine the weigts of the graph. */ + vector vwgt; + vector > adjwgt; + DetermineFEMGraphWeights(config, localFaces, adjacency, mapExternalElemIDToTimeLevel, vwgt, adjwgt); + + /*--- The remainder of this function should only be called if we have parallel + support with MPI. ---*/ +#ifdef HAVE_MPI + + /*--- If the ParMETIS library is compiled and linked, the colors + or determined via ParMETIS. Otherwise a linear distribution + is used, which is usually inefficient. ---*/ +#ifdef HAVE_PARMETIS + DetermineFEMColorsViaParMETIS(adjacency, vwgt, adjwgt); +#else + if (size > SINGLE_NODE) { + if (rank == MASTER_NODE) { + cout << endl; + cout << "--------------------- WARNING -------------------------------" << endl; + cout << "SU2 compiled without PARMETIS. A linear distribution is used." << endl; + cout << "This is very inefficient" << endl; + cout << "-------------------------------------------------------------" << endl; + cout << endl; + } + + /* Set the color to the current rank. */ + for (unsigned long i = 0; i < nElem; ++i) elem[i]->SetColor(rank); + } +#endif /* HAVE_PARMETIS */ +#endif /* HAVE_MPI */ +} + +void CPhysicalGeometry::DetermineGraphAdjacency(const vector& localFaces, + vector > &adjacency) { + + /*--- Define the linear partitioning of the elements. ---*/ + CLinearPartitioner elemPartitioner(Global_nElem, 0); + + /*--- Allocate the memory for the first index of adjacency. ---*/ + adjacency.resize(nElem); + + /*--- Loop over the matching faces to create the adjacency coming from internal faces. ---*/ + for(auto FI=localFaces.begin(); FI!=localFaces.end(); ++FI) { + /*--- Determine the local index of elem0, which is always stored locally, + and add elemID1 to the adjacency list. ---*/ + const unsigned long elem0 = FI->elemID0 - elemPartitioner.GetFirstIndexOnRank(rank); + adjacency[elem0].push_back(FI->elemID1); + + /*--- Check if this is not a periodic face and if the second element is + also a local element. If so, add elemID0 to the adjacency list ---*/ + if (FI->periodicIndex == 0) { + if (elemPartitioner.GetRankContainingIndex(FI->elemID1) == static_cast(rank)) { + const unsigned long elem1 = FI->elemID1 - elemPartitioner.GetFirstIndexOnRank(rank); + adjacency[elem1].push_back(FI->elemID0); + } + } + } + + /*--- It is possible that some neighbors appear multiple times due to e.g. + periodic boundary conditions. ParMETIS is not able to deal with this + situation, hence these multiple entries must be removed. ---*/ + for (unsigned long i = 0; i < nElem; ++i) { + sort(adjacency[i].begin(), adjacency[i].end()); + auto lastEntry = unique(adjacency[i].begin(), adjacency[i].end()); + adjacency[i].erase(lastEntry, adjacency[i].end()); + } + + /*--- Due to periodic boundary conditions it is also possible that self entries + are present. ParMETIS is not able to deal with self entries, hence + they must be removed as well. ---*/ + for (unsigned long i = 0; i < nElem; ++i) { + const unsigned long globalElemID = i + elemPartitioner.GetFirstIndexOnRank(rank); + unsigned long nEntriesNew = adjacency[i].size(); + + for (auto &adj : adjacency[i]) { + if (adj == globalElemID) { + adj = ULONG_MAX; + --nEntriesNew; + } + } + + sort(adjacency[i].begin(), adjacency[i].end()); + adjacency[i].resize(nEntriesNew); + } + + /*--- Possibly add the connectivities in the graph from the wall function + treatment. As these connectivities are one-sided, they must be stored + and communicated later. ---*/ + vector additionalExternalEntriesGraph; + + for (unsigned short iMarker = 0; iMarker < nMarker; ++iMarker) { + for (unsigned long l = 0; l < nElem_Bound[iMarker]; ++l) { + /* Get the global and local element ID adjacent to this boundary face. */ + const unsigned long globalElemID = bound[iMarker][l]->GetDomainElement(); + const unsigned long elemID = globalElemID - elemPartitioner.GetFirstIndexOnRank(rank); + + /* Get the number of donor elements for the wall function treatment + and the pointer to the array which stores this info. */ + const unsigned short nDonors = bound[iMarker][l]->GetNDonorsWallFunctions(); + const unsigned long* donors = bound[iMarker][l]->GetDonorsWallFunctions(); + + /* Loop over the number of donors and add the entry in the graph, + if not already present. */ + for (unsigned short i = 0; i < nDonors; ++i) { + if (donors[i] != globalElemID) { + if (!binary_search(adjacency[elemID].begin(), adjacency[elemID].end(), donors[i])) { + /* Donor not present in the graph for elemID. Add it and sort it + afterwards, such that a binary search can be applied later. */ + adjacency[elemID].push_back(donors[i]); + sort(adjacency[elemID].begin(), adjacency[elemID].end()); + + /* Check if the donor element is stored locally. */ + if (elemPartitioner.GetRankContainingIndex(donors[i]) == static_cast(rank)) { + /* Donor is stored locally. Add the entry to the graph + and sort it afterwards. */ + const unsigned long localDonorID = donors[i] - elemPartitioner.GetFirstIndexOnRank(rank); + adjacency[localDonorID].push_back(globalElemID); + sort(adjacency[localDonorID].begin(), adjacency[localDonorID].end()); + } else { + /* Donor is stored externally. Store the graph entry in + additionalExternalEntriesGraph. */ + additionalExternalEntriesGraph.push_back(donors[i]); + additionalExternalEntriesGraph.push_back(globalElemID); + } + } + } + } + } + } + +#ifdef HAVE_MPI + /*--- Create the send buffers with the additional graph data and determine + to which ranks data is sent. ---*/ + vector > sendBufsGraphData(size, vector(0)); + vector sendToRank(size, 0); + + for (unsigned long i = 0; i < additionalExternalEntriesGraph.size(); i += 2) { + /*--- Determine the rank where this external is stored and update + the corresponding communication buffers accordingly. ---*/ + const unsigned long rankElem = elemPartitioner.GetRankContainingIndex(additionalExternalEntriesGraph[i]); + + sendBufsGraphData[rankElem].push_back(additionalExternalEntriesGraph[i]); + sendBufsGraphData[rankElem].push_back(additionalExternalEntriesGraph[i + 1]); + sendToRank[rankElem] = 1; + } + + /*-- Determine to how many ranks this rank will send data and from how + many ranks it will receive data. ---*/ + int nRankSend = 0; + for (int i = 0; i < size; ++i) { + if (sendToRank[i]) ++nRankSend; + } + + int nRankRecv; + vector sizeSend(size, 1); + SU2_MPI::Reduce_scatter(sendToRank.data(), &nRankRecv, sizeSend.data(), MPI_INT, MPI_SUM, SU2_MPI::GetComm()); + + /* Send the data using non-blocking sends. */ + vector sendReqs(nRankSend); + nRankSend = 0; + for (int i = 0; i < size; ++i) { + if (sendToRank[i]) + SU2_MPI::Isend(sendBufsGraphData[i].data(), sendBufsGraphData[i].size(), MPI_UNSIGNED_LONG, i, i, + SU2_MPI::GetComm(), &sendReqs[nRankSend++]); + } + + /* Loop over the number of ranks from which this rank receives data. */ + for (int i = 0; i < nRankRecv; ++i) { + /* Block until a message with unsigned longs arrives from any processor. + Determine the source and the size of the message. */ + SU2_MPI::Status status; + SU2_MPI::Probe(MPI_ANY_SOURCE, rank, SU2_MPI::GetComm(), &status); + int source = status.MPI_SOURCE; + + int sizeMess; + SU2_MPI::Get_count(&status, MPI_UNSIGNED_LONG, &sizeMess); + + /* Allocate the memory for the receive buffer and receive the message + using a blocking receive. */ + vector recvBuf(sizeMess); + SU2_MPI::Recv(recvBuf.data(), sizeMess, MPI_UNSIGNED_LONG, source, rank, SU2_MPI::GetComm(), &status); + + /* Loop over the contents of the receive buffer and update the + graph accordingly. */ + for (int j = 0; j < sizeMess; j += 2) { + const unsigned long elemID = recvBuf[j] - elemPartitioner.GetFirstIndexOnRank(rank); + adjacency[elemID].push_back(recvBuf[j + 1]); + sort(adjacency[elemID].begin(), adjacency[elemID].end()); + } + } + + /* Complete the non-blocking sends and synchronize the ranks, because + wild cards have been used in the above communication. */ + SU2_MPI::Waitall(nRankSend, sendReqs.data(), MPI_STATUSES_IGNORE); + SU2_MPI::Barrier(SU2_MPI::GetComm()); +#endif +} + +void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig *config, vector &localFaces) { + + /*--- Loop over all elements to determine the faces of the elements. ---*/ for (unsigned long k = 0; k < nElem; k++) { /*--- Get the global IDs of the corner points of all the faces of this element. ---*/ unsigned short nFaces; @@ -471,430 +709,134 @@ void CPhysicalGeometry::SetColorFEMGrid_Parallel(CConfig* config) { SU2_MPI::Recv(recvBuf.data(), sizeMess, MPI_UNSIGNED_LONG, status.MPI_SOURCE, rank + 1, SU2_MPI::GetComm(), &status); - sizeMess /= 9; - unsigned long jj = 0; - for (unsigned long j = 0; j < static_cast(sizeMess); ++j, jj += 9) { - CFaceOfElement thisFace; - thisFace.nCornerPoints = recvBuf[jj]; - thisFace.cornerPoints[0] = recvBuf[jj + 1]; - thisFace.cornerPoints[1] = recvBuf[jj + 2]; - thisFace.cornerPoints[2] = recvBuf[jj + 3]; - thisFace.cornerPoints[3] = recvBuf[jj + 4]; - - vector::iterator low; - low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); - low->elemID1 = recvBuf[jj + 5]; - low->nPolySol1 = recvBuf[jj + 6]; - low->nDOFsElem1 = recvBuf[jj + 7]; - low->elemType1 = recvBuf[jj + 8]; - } - } - - /*--- Complete the second round of non-blocking sends. ---*/ - SU2_MPI::Waitall(nMessRecv, commReqs.data(), MPI_STATUSES_IGNORE); - - /*--- Wild cards have been used in the communication, so - synchronize the ranks to avoid problems. ---*/ - SU2_MPI::Barrier(SU2_MPI::GetComm()); - -#endif - - /*--- In the procedure above the periodic boundaries are not found. - A different treatment must be used in order to find these. ---*/ - DeterminePeriodicFacesFEMGrid(config, localFaces); - - /*--- Determine the total number of non-matching faces in the grid. ---*/ - nFacesLoc = localFaces.size(); - nFacesLocOr = nFacesLoc; - for (unsigned long i = 0; i < nFacesLocOr; ++i) { - if (localFaces[i].elemID1 > Global_nElem) { - localFaces[i].nCornerPoints = 4; - localFaces[i].cornerPoints[0] = Global_nPoint; - localFaces[i].cornerPoints[1] = Global_nPoint; - localFaces[i].cornerPoints[2] = Global_nPoint; - localFaces[i].cornerPoints[3] = Global_nPoint; - --nFacesLoc; - } - } - - nFacesLocOr -= nFacesLoc; - if (nFacesLocOr) { - sort(localFaces.begin(), localFaces.end()); - localFaces.resize(nFacesLoc); - } - - unsigned long nNonMatchingFaces = nFacesLocOr; - -#ifdef HAVE_MPI - SU2_MPI::Reduce(&nFacesLocOr, &nNonMatchingFaces, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm()); -#endif - if (rank == MASTER_NODE && nNonMatchingFaces) { - cout << "There are " << nNonMatchingFaces << " non-matching faces in the grid. " - << "These are ignored in the partitioning." << endl; - } - - /*--- Determine whether or not the Jacobians of the elements and faces - can be considered constant. ---*/ - DetermineFEMConstantJacobiansAndLenScale(config); - - /*--- Determine the donor elements for the wall function treatment. ---*/ - DetermineDonorElementsWallFunctions(config); - - /*--- Determine the time levels of the elements. This is only relevant - when time accurate local time stepping is employed. ---*/ - map mapExternalElemIDToTimeLevel; - DetermineTimeLevelElements(config, localFaces, mapExternalElemIDToTimeLevel); - - /*--- Define the linear partitioning of the elements. ---*/ - CLinearPartitioner elemPartitioner(Global_nElem, 0); - - /*--- Determine the ownership of the internal faces, i.e. which adjacent - element is responsible for computing the fluxes through the face. ---*/ - for (unsigned long i = 0; i < nFacesLoc; ++i) { - /* Check for a matching face. */ - if (localFaces[i].elemID1 < Global_nElem) { - /*--- Determine the time level of Elem0, which is always owned. ---*/ - const unsigned long elemID0 = localFaces[i].elemID0 - elemPartitioner.GetFirstIndexOnRank(rank); - const unsigned short timeLevel0 = elem[elemID0]->GetTimeLevel(); - - /*--- Determine the time level of Elem1, which is either owned or - external. Hence a distinction must be made. ---*/ - unsigned short timeLevel1; - if (elemPartitioner.GetRankContainingIndex(localFaces[i].elemID1) == static_cast(rank)) { - const unsigned long elemID1 = localFaces[i].elemID1 - elemPartitioner.GetFirstIndexOnRank(rank); - timeLevel1 = elem[elemID1]->GetTimeLevel(); - } else { - const auto MI = mapExternalElemIDToTimeLevel.find(localFaces[i].elemID1); - if (MI == mapExternalElemIDToTimeLevel.end()) - SU2_MPI::Error("Entry not found in mapExternalElemIDToTimeLevel", CURRENT_FUNCTION); - timeLevel1 = MI->second.short0; - } - - /* Check if both elements have the same time level. */ - if (timeLevel0 == timeLevel1) { - /* Same time level, hence both elements can own the face. First check whether - elemID0 == elemID1 (which happens for periodic problems with only one - element in the periodic direction), because this is a special case. */ - if (localFaces[i].elemID0 == localFaces[i].elemID1) { - /* This face occurs twice, but should be owned only once. Base this - decision on the periodic index. */ - localFaces[i].elem0IsOwner = localFaces[i].periodicIndex < localFaces[i].periodicIndexDonor; - } else { - /* Different elements on both sides of the face. The ownership decision - below makes an attempt to spread the workload evenly. */ - const unsigned long sumElemID = localFaces[i].elemID0 + localFaces[i].elemID1; - if (sumElemID % 2) - localFaces[i].elem0IsOwner = localFaces[i].elemID0 < localFaces[i].elemID1; - else - localFaces[i].elem0IsOwner = localFaces[i].elemID0 > localFaces[i].elemID1; - } - } else { - /* The time level of both elements differ. The element with the smallest - time level must be the owner of the element. */ - localFaces[i].elem0IsOwner = timeLevel0 < timeLevel1; - } - } else { - /* Non-matching face. Give the ownership to element 0. */ - localFaces[i].elem0IsOwner = true; - } - } - - /*--- All the matching face information is known now, including periodic - faces. Store the information of the neighbors in the data structure - for the local elements. ---*/ - for (unsigned long k = 0; k < nElem; ++k) { - unsigned short nFaces; - unsigned short nPointsPerFace[6]; - unsigned long faceConn[6][4]; - - elem[k]->GetCornerPointsAllFaces(nFaces, nPointsPerFace, faceConn); - elem[k]->InitializeNeighbors(nFaces); - - for (unsigned short i = 0; i < nFaces; ++i) { - CFaceOfElement thisFace; - thisFace.nCornerPoints = nPointsPerFace[i]; - for (unsigned short j = 0; j < nPointsPerFace[i]; ++j) thisFace.cornerPoints[j] = faceConn[i][j]; - thisFace.elemID0 = k + elemPartitioner.GetFirstIndexOnRank(rank); - thisFace.nPolySol0 = elem[k]->GetNPolySol(); - thisFace.nDOFsElem0 = elem[k]->GetNDOFsSol(); - - thisFace.CreateUniqueNumbering(); - - vector::iterator low; - low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); - - if (low != localFaces.end()) { - if (!(thisFace < *low)) { - if (low->elemID0 == thisFace.elemID0) { - elem[k]->SetNeighbor_Elements(low->elemID1, i); - elem[k]->SetOwnerFace(low->elem0IsOwner, i); - } else { - elem[k]->SetNeighbor_Elements(low->elemID0, i); - elem[k]->SetOwnerFace(!(low->elem0IsOwner), i); - } - - if (low->periodicIndex > 0) elem[k]->SetPeriodicIndex(low->periodicIndex - 1, i); - } - } - } - } - - /*--- Create the vector of vectors that describe the connectivity - of the graph. First the faces. ---*/ - vector > adjacency(nElem, vector(0)); - for (unsigned long i = 0; i < nFacesLoc; ++i) { - /*--- Determine the local index of elem0, which is always stored locally, - and add elemID1 to the adjacency list. ---*/ - const unsigned long elem0 = localFaces[i].elemID0 - elemPartitioner.GetFirstIndexOnRank(rank); - adjacency[elem0].push_back(localFaces[i].elemID1); - - /*--- Check if this is not a periodic face and if the second element is - also a local element. If so, add elemID0 to the adjacency list ---*/ - if (localFaces[i].periodicIndex == 0) { - if (elemPartitioner.GetRankContainingIndex(localFaces[i].elemID1) == static_cast(rank)) { - const unsigned long elem1 = localFaces[i].elemID1 - elemPartitioner.GetFirstIndexOnRank(rank); - adjacency[elem1].push_back(localFaces[i].elemID0); - } - } - } - - /* It is possible that some neighbors appear multiple times due to e.g. - periodic boundary conditions. ParMETIS is not able to deal with this - situation, hence these multiple entries must be removed. */ - for (unsigned long i = 0; i < nElem; ++i) { - sort(adjacency[i].begin(), adjacency[i].end()); - vector::iterator lastEntry; - lastEntry = unique(adjacency[i].begin(), adjacency[i].end()); - adjacency[i].erase(lastEntry, adjacency[i].end()); - } - - /* Due to periodic boundary conditions it is also possible that self entries - are present. ParMETIS is not able to deal with self entries, hence - they must be removed as well. */ - for (unsigned long i = 0; i < nElem; ++i) { - const unsigned long globalElemID = i + +elemPartitioner.GetFirstIndexOnRank(rank); - unsigned long nEntriesNew = adjacency[i].size(); - - for (unsigned long j = 0; j < adjacency[i].size(); ++j) { - if (adjacency[i][j] == globalElemID) { - adjacency[i][j] = ULONG_MAX; - --nEntriesNew; - } - } - - sort(adjacency[i].begin(), adjacency[i].end()); - adjacency[i].resize(nEntriesNew); - } - - /*--- Possibly add the connectivities in the graph from the wall function - treatment. As these connectivities are one-sided, they must be stored - and communicated later. ---*/ - vector additionalExternalEntriesGraph; - - for (unsigned short iMarker = 0; iMarker < nMarker; ++iMarker) { - for (unsigned long l = 0; l < nElem_Bound[iMarker]; ++l) { - /* Get the global and local element ID adjacent to this boundary face. */ - const unsigned long globalElemID = bound[iMarker][l]->GetDomainElement(); - const unsigned long elemID = globalElemID - elemPartitioner.GetFirstIndexOnRank(rank); - - /* Get the number of donor elements for the wall function treatment - and the pointer to the array which stores this info. */ - const unsigned short nDonors = bound[iMarker][l]->GetNDonorsWallFunctions(); - const unsigned long* donors = bound[iMarker][l]->GetDonorsWallFunctions(); - - /* Loop over the number of donors and add the entry in the graph, - if not already present. */ - for (unsigned short i = 0; i < nDonors; ++i) { - if (donors[i] != globalElemID) { - if (!binary_search(adjacency[elemID].begin(), adjacency[elemID].end(), donors[i])) { - /* Donor not present in the graph for elemID. Add it and sort it - afterwards, such that a binary search can be applied later. */ - adjacency[elemID].push_back(donors[i]); - sort(adjacency[elemID].begin(), adjacency[elemID].end()); - - /* Check if the donor element is stored locally. */ - if (elemPartitioner.GetRankContainingIndex(donors[i]) == static_cast(rank)) { - /* Donor is stored locally. Add the entry to the graph - and sort it afterwards. */ - const unsigned long localDonorID = donors[i] - elemPartitioner.GetFirstIndexOnRank(rank); - adjacency[localDonorID].push_back(globalElemID); - sort(adjacency[localDonorID].begin(), adjacency[localDonorID].end()); - } else { - /* Donor is stored externally. Store the graph entry in - additionalExternalEntriesGraph. */ - additionalExternalEntriesGraph.push_back(donors[i]); - additionalExternalEntriesGraph.push_back(globalElemID); - } - } - } - } - } - } - -#ifdef HAVE_MPI - /*--- Create the send buffers with the additional graph data and determine - to which ranks data is sent. ---*/ - vector > sendBufsGraphData(size, vector(0)); - vector sendToRank(size, 0); - - for (unsigned long i = 0; i < additionalExternalEntriesGraph.size(); i += 2) { - /*--- Determine the rank where this external is stored and update - the corresponding communication buffers accordingly. ---*/ - const unsigned long rankElem = elemPartitioner.GetRankContainingIndex(additionalExternalEntriesGraph[i]); - - sendBufsGraphData[rankElem].push_back(additionalExternalEntriesGraph[i]); - sendBufsGraphData[rankElem].push_back(additionalExternalEntriesGraph[i + 1]); - sendToRank[rankElem] = 1; - } - - /*-- Determine to how many ranks this rank will send data and from how - many ranks it will receive data. ---*/ - int nRankSend = 0; - for (int i = 0; i < size; ++i) { - if (sendToRank[i]) ++nRankSend; - } - - int nRankRecv; - vector sizeSend(size, 1); - SU2_MPI::Reduce_scatter(sendToRank.data(), &nRankRecv, sizeSend.data(), MPI_INT, MPI_SUM, SU2_MPI::GetComm()); - - /* Send the data using non-blocking sends. */ - vector sendReqs(nRankSend); - nRankSend = 0; - for (int i = 0; i < size; ++i) { - if (sendToRank[i]) - SU2_MPI::Isend(sendBufsGraphData[i].data(), sendBufsGraphData[i].size(), MPI_UNSIGNED_LONG, i, i, - SU2_MPI::GetComm(), &sendReqs[nRankSend++]); - } - - /* Loop over the number of ranks from which this rank receives data. */ - for (int i = 0; i < nRankRecv; ++i) { - /* Block until a message with unsigned longs arrives from any processor. - Determine the source and the size of the message. */ - SU2_MPI::Status status; - SU2_MPI::Probe(MPI_ANY_SOURCE, rank, SU2_MPI::GetComm(), &status); - int source = status.MPI_SOURCE; - - int sizeMess; - SU2_MPI::Get_count(&status, MPI_UNSIGNED_LONG, &sizeMess); - - /* Allocate the memory for the receive buffer and receive the message - using a blocking receive. */ - vector recvBuf(sizeMess); - SU2_MPI::Recv(recvBuf.data(), sizeMess, MPI_UNSIGNED_LONG, source, rank, SU2_MPI::GetComm(), &status); - - /* Loop over the contents of the receive buffer and update the - graph accordingly. */ - for (int j = 0; j < sizeMess; j += 2) { - const unsigned long elemID = recvBuf[j] - elemPartitioner.GetFirstIndexOnRank(rank); - adjacency[elemID].push_back(recvBuf[j + 1]); - sort(adjacency[elemID].begin(), adjacency[elemID].end()); + sizeMess /= 9; + unsigned long jj = 0; + for (unsigned long j = 0; j < static_cast(sizeMess); ++j, jj += 9) { + CFaceOfElement thisFace; + thisFace.nCornerPoints = recvBuf[jj]; + thisFace.cornerPoints[0] = recvBuf[jj + 1]; + thisFace.cornerPoints[1] = recvBuf[jj + 2]; + thisFace.cornerPoints[2] = recvBuf[jj + 3]; + thisFace.cornerPoints[3] = recvBuf[jj + 4]; + + vector::iterator low; + low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); + low->elemID1 = recvBuf[jj + 5]; + low->nPolySol1 = recvBuf[jj + 6]; + low->nDOFsElem1 = recvBuf[jj + 7]; + low->elemType1 = recvBuf[jj + 8]; } } - /* Complete the non-blocking sends amd synchronize the ranks, because - wild cards have been used in the above communication. */ - SU2_MPI::Waitall(nRankSend, sendReqs.data(), MPI_STATUSES_IGNORE); + /*--- Complete the second round of non-blocking sends. ---*/ + SU2_MPI::Waitall(nMessRecv, commReqs.data(), MPI_STATUSES_IGNORE); + + /*--- Wild cards have been used in the communication, so + synchronize the ranks to avoid problems. ---*/ SU2_MPI::Barrier(SU2_MPI::GetComm()); #endif +} - /* Allocate the memory to store the weights of the graph. */ - vector vwgt(2 * nElem); - vector > adjwgt(nElem, vector(0)); - - for (unsigned long i = 0; i < nElem; ++i) adjwgt[i].resize(adjacency[i].size()); - - /* Compute the weigts of the graph. */ - ComputeFEMGraphWeights(config, localFaces, adjacency, mapExternalElemIDToTimeLevel, vwgt, adjwgt); - - /*--- The remainder of this function should only be called if we have parallel - support with MPI and have the ParMETIS library compiled and linked. ---*/ - -#ifdef HAVE_MPI -#ifdef HAVE_PARMETIS +void CPhysicalGeometry::DetermineNonMatchingFacesFEMGrid(const CConfig* config, vector& localMatchingFaces) { - /*--- Only call ParMETIS if we have more than one rank to avoid errors ---*/ - if (size > SINGLE_NODE) { - /*--- The scalar variables and the options array for the call to ParMETIS. ---*/ - idx_t wgtflag = 3; // Weights on both the vertices and edges. - idx_t numflag = 0; // C-numbering. - idx_t ncon = 2; // Number of constraints. - real_t ubvec[] = {1.05, 1.05}; // Tolerances for the vertex weights, recommended value is 1.05. - auto nparts = (idx_t)size; // Number of subdomains. Must be number of MPI ranks. - idx_t options[METIS_NOPTIONS]; // Just use the default options. - METIS_SetDefaultOptions(options); - options[1] = 0; + /*--- If there are single faces still present in localMatchingFaces, these are + non-matching faces. Store these faces in singleFaces and remove them + from localMatchingFaces. ---*/ + vector singleFaces; - /*--- Determine the array, which stores the distribution of the graph nodes - over the ranks. ---*/ - vector vtxdist(size + 1); - for (int i = 0; i <= size; ++i) vtxdist[i] = static_cast(elemPartitioner.GetCumulativeSizeBeforeRank(i)); + /*--- Determine the total number of non-matching faces in the grid. ---*/ + unsigned long nFacesLoc = localMatchingFaces.size(); + for (unsigned long i = 0; i < localMatchingFaces.size(); ++i) { + if (localMatchingFaces[i].elemID1 > Global_nElem) { - /* Create the array xadjPar, which contains the number of edges for each - vertex of the graph in ParMETIS format. */ - vector xadjPar(nElem + 1); - xadjPar[0] = 0; - for (unsigned long i = 0; i < nElem; ++i) xadjPar[i + 1] = xadjPar[i] + (idx_t)adjacency[i].size(); + singleFaces.push_back(localMatchingFaces[i]); - /* Create the adjacency information in ParMETIS format. */ - vector adjacencyPar(xadjPar[nElem]); - unsigned long ii = 0; - for (unsigned long i = 0; i < nElem; ++i) { - for (unsigned long j = 0; j < adjacency[i].size(); ++j, ++ii) adjacencyPar[ii] = adjacency[i][j]; + localMatchingFaces[i].nCornerPoints = 4; + localMatchingFaces[i].cornerPoints[0] = Global_nPoint; + localMatchingFaces[i].cornerPoints[1] = Global_nPoint; + localMatchingFaces[i].cornerPoints[2] = Global_nPoint; + localMatchingFaces[i].cornerPoints[3] = Global_nPoint; + --nFacesLoc; } + } - /* Create the vertex weights in ParMETIS format. */ - vector vwgtPar(nElem * ncon); - for (unsigned long i = 0; i < nElem * ncon; ++i) vwgtPar[i] = (idx_t)ceil(vwgt[i]); + if( singleFaces.size() ) { + sort(localMatchingFaces.begin(), localMatchingFaces.end()); + localMatchingFaces.resize(nFacesLoc); + } - /* Create the adjacency weight in ParMETIS format. */ - vector adjwgtPar(xadjPar[nElem]); - ii = 0; - for (unsigned long i = 0; i < nElem; ++i) { - for (unsigned long j = 0; j < adjwgt[i].size(); ++j, ++ii) adjwgtPar[ii] = (idx_t)ceil(adjwgt[i][j]); - } + /*--- Determine the global number of non-matching faces. ---*/ + nFacesLoc = singleFaces.size(); + unsigned long nNonMatchingFaces = nFacesLoc; - /* Make sure that an equal distribution is obtained. */ - vector tpwgts(size * ncon, 1.0 / ((real_t)size)); +#ifdef HAVE_MPI + SU2_MPI::Allreduce(&nFacesLoc, &nNonMatchingFaces, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); +#endif - /*--- Calling ParMETIS ---*/ - vector part(nElem); - if (rank == MASTER_NODE) cout << "Calling ParMETIS..."; + /*--- If there are no non-matching faces, return, such that in the rest of this + function it can be assumed that non-matching faces are present. ---*/ + if(nNonMatchingFaces == 0) return; - idx_t edgecut; - MPI_Comm comm = SU2_MPI::GetComm(); - ParMETIS_V3_PartKway(vtxdist.data(), xadjPar.data(), adjacencyPar.data(), vwgtPar.data(), adjwgtPar.data(), - &wgtflag, &numflag, &ncon, &nparts, tpwgts.data(), ubvec, options, &edgecut, part.data(), - &comm); - if (rank == MASTER_NODE) { - cout << " graph partitioning complete ("; - cout << edgecut << " edge cuts)." << endl; - } + SU2_MPI::Error(string("Handling of non-matching faces not implemented yet"), CURRENT_FUNCTION); +} - /*--- Set the color of the elements to the outcome of ParMETIS. ---*/ - for (unsigned long i = 0; i < nElem; ++i) elem[i]->SetColor(part[i]); - } +void CPhysicalGeometry::DetermineOwnershipInternalFaces(vector& localFaces, + map &mapExternalElemIDToTimeLevel) { -#else /* HAVE_PARMETIS */ + /*--- Define the linear partitioning of the elements. ---*/ + CLinearPartitioner elemPartitioner(Global_nElem, 0); - if (size > SINGLE_NODE) { - if (rank == MASTER_NODE) { - cout << endl; - cout << "--------------------- WARNING -------------------------------" << endl; - cout << "SU2 compiled without PARMETIS. A linear distribution is used." << endl; - cout << "This is very inefficient" << endl; - cout << "-------------------------------------------------------------" << endl; - cout << endl; - } + /*--- Loop over the locally stored faces. ---*/ + for(auto FI =localFaces.begin(); FI!=localFaces.end(); ++FI) { + /* Check for a matching face. */ + if (FI->elemID1 < Global_nElem) { - /* Set the color to the current rank. */ - for (unsigned long i = 0; i < nElem; ++i) elem[i]->SetColor(rank); - } + /*--- Determine the time level of Elem0, which is always owned. ---*/ + const unsigned long elemID0 = FI->elemID0 - elemPartitioner.GetFirstIndexOnRank(rank); + const unsigned short timeLevel0 = elem[elemID0]->GetTimeLevel(); -#endif /* HAVE_PARMETIS */ + /*--- Determine the time level of Elem1, which is either owned or + external. Hence a distinction must be made. ---*/ + unsigned short timeLevel1; + if (elemPartitioner.GetRankContainingIndex(FI->elemID1) == static_cast(rank)) { + const unsigned long elemID1 = FI->elemID1 - elemPartitioner.GetFirstIndexOnRank(rank); + timeLevel1 = elem[elemID1]->GetTimeLevel(); + } else { + const auto MI = mapExternalElemIDToTimeLevel.find(FI->elemID1); + if (MI == mapExternalElemIDToTimeLevel.end()) + SU2_MPI::Error("Entry not found in mapExternalElemIDToTimeLevel", CURRENT_FUNCTION); + timeLevel1 = MI->second.short0; + } -#endif /* HAVE_MPI */ + /* Check if both elements have the same time level. */ + if (timeLevel0 == timeLevel1) { + /* Same time level, hence both elements can own the face. First check whether + elemID0 == elemID1 (which happens for periodic problems with only one + element in the periodic direction), because this is a special case. */ + if (FI->elemID0 == FI->elemID1) { + /* This face occurs twice, but should be owned only once. Base this + decision on the periodic index. */ + FI->elem0IsOwner = FI->periodicIndex < FI->periodicIndexDonor; + } else { + /* Different elements on both sides of the face. The ownership decision + below makes an attempt to spread the workload evenly. */ + const unsigned long sumElemID = FI->elemID0 + FI->elemID1; + if (sumElemID % 2) + FI->elem0IsOwner = FI->elemID0 < FI->elemID1; + else + FI->elem0IsOwner = FI->elemID0 > FI->elemID1; + } + } else { + /* The time level of both elements differ. The element with the smallest + time level must be the owner of the element. */ + FI->elem0IsOwner = timeLevel0 < timeLevel1; + } + } else { + /* Non-matching face. Give the ownership to element 0. */ + FI->elem0IsOwner = true; + } + } } void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector& localFaces) { @@ -1195,6 +1137,81 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector > &adjacency, + vector &vwgt, + vector > &adjwgt) { + + /*--- Only call ParMETIS if we have more than one rank to avoid errors ---*/ + if (size > SINGLE_NODE) { + + /*--- Define the linear partitioning of the elements. ---*/ + CLinearPartitioner elemPartitioner(Global_nElem, 0); + + /*--- The scalar variables and the options array for the call to ParMETIS. ---*/ + idx_t wgtflag = 3; // Weights on both the vertices and edges. + idx_t numflag = 0; // C-numbering. + idx_t ncon = 2; // Number of constraints. + real_t ubvec[] = {1.05, 1.05}; // Tolerances for the vertex weights, recommended value is 1.05. + auto nparts = (idx_t)size; // Number of subdomains. Must be number of MPI ranks. + idx_t options[METIS_NOPTIONS]; // Just use the default options. + METIS_SetDefaultOptions(options); + options[1] = 0; + + /*--- Determine the array, which stores the distribution of the graph nodes + over the ranks. ---*/ + vector vtxdist(size + 1); + for (int i = 0; i <= size; ++i) vtxdist[i] = static_cast(elemPartitioner.GetCumulativeSizeBeforeRank(i)); + + /* Create the array xadjPar, which contains the number of edges for each + vertex of the graph in ParMETIS format. */ + vector xadjPar(nElem + 1); + xadjPar[0] = 0; + for (unsigned long i = 0; i < nElem; ++i) xadjPar[i + 1] = xadjPar[i] + (idx_t)adjacency[i].size(); + + /* Create the adjacency information in ParMETIS format. */ + vector adjacencyPar(xadjPar[nElem]); + unsigned long ii = 0; + for (unsigned long i = 0; i < nElem; ++i) { + for (unsigned long j = 0; j < adjacency[i].size(); ++j, ++ii) adjacencyPar[ii] = adjacency[i][j]; + } + + /* Create the vertex weights in ParMETIS format. */ + vector vwgtPar(nElem * ncon); + for (unsigned long i = 0; i < nElem * ncon; ++i) vwgtPar[i] = (idx_t)ceil(vwgt[i]); + + /* Create the adjacency weight in ParMETIS format. */ + vector adjwgtPar(xadjPar[nElem]); + ii = 0; + for (unsigned long i = 0; i < nElem; ++i) { + for (unsigned long j = 0; j < adjwgt[i].size(); ++j, ++ii) adjwgtPar[ii] = (idx_t)ceil(adjwgt[i][j]); + } + + /* Make sure that an equal distribution is obtained. */ + vector tpwgts(size * ncon, 1.0 / ((real_t)size)); + + /*--- Calling ParMETIS ---*/ + vector part(nElem); + if (rank == MASTER_NODE) cout << "Calling ParMETIS..."; + + idx_t edgecut; + MPI_Comm comm = SU2_MPI::GetComm(); + ParMETIS_V3_PartKway(vtxdist.data(), xadjPar.data(), adjacencyPar.data(), vwgtPar.data(), adjwgtPar.data(), + &wgtflag, &numflag, &ncon, &nparts, tpwgts.data(), ubvec, options, &edgecut, part.data(), + &comm); + if (rank == MASTER_NODE) { + cout << " graph partitioning complete ("; + cout << edgecut << " edge cuts)." << endl; + } + + /*--- Set the color of the elements to the outcome of ParMETIS. ---*/ + for (unsigned long i = 0; i < nElem; ++i) elem[i]->SetColor(part[i]); + } +} +#endif +#endif + void CPhysicalGeometry::DetermineFEMConstantJacobiansAndLenScale(CConfig* config) { /* Definition of the object that is used to carry out the BLAS calls. */ CBlasStructure blasFunctions; @@ -2681,10 +2698,61 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector } } -void CPhysicalGeometry::ComputeFEMGraphWeights(CConfig* config, const vector& localFaces, - const vector >& adjacency, - const map& mapExternalElemIDToTimeLevel, - vector& vwgt, vector >& adjwgt) { +void CPhysicalGeometry::StoreFaceInfoInLocalElements(const vector &localFaces) { + + /*--- Define the linear partitioning of the elements. ---*/ + CLinearPartitioner elemPartitioner(Global_nElem, 0); + + /*--- Loop over the locally stored elements. ---*/ + for (unsigned long k = 0; k < nElem; ++k) { + + /*--- Get the number of faces and the corner points of these + faces of this element. ---*/ + unsigned short nFaces; + unsigned short nPointsPerFace[6]; + unsigned long faceConn[6][4]; + elem[k]->GetCornerPointsAllFaces(nFaces, nPointsPerFace, faceConn); + + /*--- Initialize the data structures for the neighbors. ---*/ + elem[k]->InitializeNeighbors(nFaces); + + /*--- Loop over the faces of the element. ---*/ + for (unsigned short i = 0; i < nFaces; ++i) { + + /*--- Create an object of the class CFaceOfElement for this face. ---*/ + CFaceOfElement thisFace; + thisFace.nCornerPoints = nPointsPerFace[i]; + for (unsigned short j = 0; j < nPointsPerFace[i]; ++j) thisFace.cornerPoints[j] = faceConn[i][j]; + thisFace.elemID0 = k + elemPartitioner.GetFirstIndexOnRank(rank); + thisFace.nPolySol0 = elem[k]->GetNPolySol(); + thisFace.nDOFsElem0 = elem[k]->GetNDOFsSol(); + + thisFace.CreateUniqueNumbering(); + + /*--- Search for this face in localFaces. Store the neighboring element and + the ownership of this face. Store the periodic index for a periodic face. ---*/ + auto low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); + if (low != localFaces.end()) { + if (!(thisFace < *low)) { + if (low->elemID0 == thisFace.elemID0) { + elem[k]->SetNeighbor_Elements(low->elemID1, i); + elem[k]->SetOwnerFace(low->elem0IsOwner, i); + } else { + elem[k]->SetNeighbor_Elements(low->elemID0, i); + elem[k]->SetOwnerFace(!(low->elem0IsOwner), i); + } + + if (low->periodicIndex > 0) elem[k]->SetPeriodicIndex(low->periodicIndex - 1, i); + } + } + } + } +} + +void CPhysicalGeometry::DetermineFEMGraphWeights(CConfig* config, const vector& localFaces, + const vector >& adjacency, + const map& mapExternalElemIDToTimeLevel, + vector& vwgt, vector >& adjwgt) { /*--- Define the linear partitioning of the elements. ---*/ CLinearPartitioner elemPartitioner(Global_nElem, 0); @@ -2697,6 +2765,13 @@ void CPhysicalGeometry::ComputeFEMGraphWeights(CConfig* config, const vector::const_iterator low; - low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); + auto low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); bool thisFaceFound = false; if (low != localFaces.end()) { From 1773588767b8444eb6027fbaffd0e98562e8aa1e Mon Sep 17 00:00:00 2001 From: vdweide Date: Mon, 10 Feb 2025 17:46:31 +0100 Subject: [PATCH 3/9] Replaced the explicit iterators by auto to improve readability --- Common/include/geometry/CPhysicalGeometry.hpp | 39 +++--- Common/src/geometry/CPhysicalGeometryFEM.cpp | 117 +++++++----------- 2 files changed, 64 insertions(+), 92 deletions(-) diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index 92dd306b34f..7257f546d40 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/Common/include/geometry/CPhysicalGeometry.hpp @@ -506,20 +506,19 @@ class CPhysicalGeometry final : public CGeometry { */ void DetermineDonorElementsWallFunctions(CConfig* config); - #ifdef HAVE_MPI - #ifdef HAVE_PARMETIS - /*! - * \brief Function, which converts the input the format for ParMETIS and calls - * ParMETIS to determine the actual colors of the elements. - * \param[in] adjacency - Adjacency information of the elements. - * \param[in] vwgt - Weights of the vertices of the graph, which are the elements. - * \param[in] adjwgt - Weights of the adjacencies of the graph. - */ - void DetermineFEMColorsViaParMETIS(vector > &adjacency, - vector &vwgt, - vector > &adjwgt); - #endif - #endif +#ifdef HAVE_MPI +#ifdef HAVE_PARMETIS + /*! + * \brief Function, which converts the input the format for ParMETIS and calls + * ParMETIS to determine the actual colors of the elements. + * \param[in] adjacency - Adjacency information of the elements. + * \param[in] vwgt - Weights of the vertices of the graph, which are the elements. + * \param[in] adjwgt - Weights of the adjacencies of the graph. + */ + void DetermineFEMColorsViaParMETIS(vector >& adjacency, vector& vwgt, + vector >& adjwgt); +#endif +#endif /*! * \brief Determine whether or not the Jacobians of the elements and faces @@ -549,8 +548,7 @@ class CPhysicalGeometry final : public CGeometry { * \param[in] localFaces - Vector containing the local matching faces of the FEM grid. * \param[out] adjacency - Vector of vectors to store the adjacency. */ - void DetermineGraphAdjacency(const vector& localFaces, - vector >& adjacency); + void DetermineGraphAdjacency(const vector& localFaces, vector >& adjacency); /*! * \brief Function, which determines the matching faces of a FEM grid. @@ -572,10 +570,11 @@ class CPhysicalGeometry final : public CGeometry { * \brief Function, which determines the owner of the internal faces, i.e. which element * is responsible for computing the fluxes through the face. * \param[in] localFaces - Vector, which contains the element faces of this rank. - * \param[out] mapExternalElemIDToTimeLevel - Map from the external element ID's to their time level and number of DOFs. + * \param[out] mapExternalElemIDToTimeLevel - Map from the external element ID's to their time level and number of + * DOFs. */ - void DetermineOwnershipInternalFaces(vector &localFaces, - map &mapExternalElemIDToTimeLevel); + void DetermineOwnershipInternalFaces(vector& localFaces, + map& mapExternalElemIDToTimeLevel); /*! * \brief Determine the neighboring information for periodic faces of a FEM grid. @@ -599,7 +598,7 @@ class CPhysicalGeometry final : public CGeometry { * data structures of the local elements. * \param[in] localFaces - Vector, which contains the internal matching faces of this rank. */ - void StoreFaceInfoInLocalElements(const vector &localFaces); + void StoreFaceInfoInLocalElements(const vector& localFaces); /*! * \brief Compute 3 grid quality metrics: orthogonality angle, dual cell aspect ratio, and dual cell volume ratio. diff --git a/Common/src/geometry/CPhysicalGeometryFEM.cpp b/Common/src/geometry/CPhysicalGeometryFEM.cpp index ee2857a269b..f8c18b66dea 100644 --- a/Common/src/geometry/CPhysicalGeometryFEM.cpp +++ b/Common/src/geometry/CPhysicalGeometryFEM.cpp @@ -167,7 +167,6 @@ void CPhysicalGeometry::LoadLinearlyPartitionedSurfaceElementsFEM(CConfig* confi } void CPhysicalGeometry::SetColorFEMGrid_Parallel(CConfig* config) { - /*--- Initialize the color vector of the elements. ---*/ for (unsigned long i = 0; i < nElem; ++i) elem[i]->SetColor(0); @@ -242,8 +241,7 @@ void CPhysicalGeometry::SetColorFEMGrid_Parallel(CConfig* config) { } void CPhysicalGeometry::DetermineGraphAdjacency(const vector& localFaces, - vector > &adjacency) { - + vector >& adjacency) { /*--- Define the linear partitioning of the elements. ---*/ CLinearPartitioner elemPartitioner(Global_nElem, 0); @@ -251,7 +249,7 @@ void CPhysicalGeometry::DetermineGraphAdjacency(const vector& lo adjacency.resize(nElem); /*--- Loop over the matching faces to create the adjacency coming from internal faces. ---*/ - for(auto FI=localFaces.begin(); FI!=localFaces.end(); ++FI) { + for (auto FI = localFaces.begin(); FI != localFaces.end(); ++FI) { /*--- Determine the local index of elem0, which is always stored locally, and add elemID1 to the adjacency list. ---*/ const unsigned long elem0 = FI->elemID0 - elemPartitioner.GetFirstIndexOnRank(rank); @@ -283,7 +281,7 @@ void CPhysicalGeometry::DetermineGraphAdjacency(const vector& lo const unsigned long globalElemID = i + elemPartitioner.GetFirstIndexOnRank(rank); unsigned long nEntriesNew = adjacency[i].size(); - for (auto &adj : adjacency[i]) { + for (auto& adj : adjacency[i]) { if (adj == globalElemID) { adj = ULONG_MAX; --nEntriesNew; @@ -407,8 +405,7 @@ void CPhysicalGeometry::DetermineGraphAdjacency(const vector& lo #endif } -void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig *config, vector &localFaces) { - +void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig* config, vector& localFaces) { /*--- Loop over all elements to determine the faces of the elements. ---*/ for (unsigned long k = 0; k < nElem; k++) { /*--- Get the global IDs of the corner points of all the faces of this element. ---*/ @@ -454,8 +451,7 @@ void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig *config, vec for (unsigned short j = 0; j < nPointsPerFace[0]; ++j) thisFace.cornerPoints[j] = faceConn[0][j]; thisFace.CreateUniqueNumbering(); - vector::iterator low; - low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); + auto low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); /*--- Store the corresponding element ID for the boundary element and invalidate the face by setting the element ID to an invalid value. ---*/ @@ -542,8 +538,7 @@ void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig *config, vec does not have to be repeated below. ---*/ vector nFacesComm(size, 0); for (unsigned long i = 0; i < nFacesLocComm; ++i) { - vector::iterator low; - low = lower_bound(facePointsProc.begin(), facePointsProc.end(), localFacesComm[i].cornerPoints[0]); + auto low = lower_bound(facePointsProc.begin(), facePointsProc.end(), localFacesComm[i].cornerPoints[0]); unsigned long rankFace = low - facePointsProc.begin(); if (*low > localFacesComm[i].cornerPoints[0]) --rankFace; @@ -676,8 +671,7 @@ void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig *config, vec sendBufFace[ii + 3] = facesRecv[j].cornerPoints[2]; sendBufFace[ii + 4] = facesRecv[j].cornerPoints[3]; - vector::iterator low; - low = lower_bound(localFacesComm.begin(), localFacesComm.end(), facesRecv[j]); + auto low = lower_bound(localFacesComm.begin(), localFacesComm.end(), facesRecv[j]); if (facesRecv[j].elemID0 == low->elemID0) { sendBufFace[ii + 5] = low->elemID1; sendBufFace[ii + 6] = low->nPolySol1; @@ -719,8 +713,7 @@ void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig *config, vec thisFace.cornerPoints[2] = recvBuf[jj + 3]; thisFace.cornerPoints[3] = recvBuf[jj + 4]; - vector::iterator low; - low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); + auto low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); low->elemID1 = recvBuf[jj + 5]; low->nPolySol1 = recvBuf[jj + 6]; low->nDOFsElem1 = recvBuf[jj + 7]; @@ -738,8 +731,8 @@ void CPhysicalGeometry::DetermineMatchingFacesFEMGrid(const CConfig *config, vec #endif } -void CPhysicalGeometry::DetermineNonMatchingFacesFEMGrid(const CConfig* config, vector& localMatchingFaces) { - +void CPhysicalGeometry::DetermineNonMatchingFacesFEMGrid(const CConfig* config, + vector& localMatchingFaces) { /*--- If there are single faces still present in localMatchingFaces, these are non-matching faces. Store these faces in singleFaces and remove them from localMatchingFaces. ---*/ @@ -749,7 +742,6 @@ void CPhysicalGeometry::DetermineNonMatchingFacesFEMGrid(const CConfig* config, unsigned long nFacesLoc = localMatchingFaces.size(); for (unsigned long i = 0; i < localMatchingFaces.size(); ++i) { if (localMatchingFaces[i].elemID1 > Global_nElem) { - singleFaces.push_back(localMatchingFaces[i]); localMatchingFaces[i].nCornerPoints = 4; @@ -761,7 +753,7 @@ void CPhysicalGeometry::DetermineNonMatchingFacesFEMGrid(const CConfig* config, } } - if( singleFaces.size() ) { + if (singleFaces.size()) { sort(localMatchingFaces.begin(), localMatchingFaces.end()); localMatchingFaces.resize(nFacesLoc); } @@ -776,22 +768,20 @@ void CPhysicalGeometry::DetermineNonMatchingFacesFEMGrid(const CConfig* config, /*--- If there are no non-matching faces, return, such that in the rest of this function it can be assumed that non-matching faces are present. ---*/ - if(nNonMatchingFaces == 0) return; + if (nNonMatchingFaces == 0) return; SU2_MPI::Error(string("Handling of non-matching faces not implemented yet"), CURRENT_FUNCTION); } -void CPhysicalGeometry::DetermineOwnershipInternalFaces(vector& localFaces, - map &mapExternalElemIDToTimeLevel) { - +void CPhysicalGeometry::DetermineOwnershipInternalFaces( + vector& localFaces, map& mapExternalElemIDToTimeLevel) { /*--- Define the linear partitioning of the elements. ---*/ CLinearPartitioner elemPartitioner(Global_nElem, 0); /*--- Loop over the locally stored faces. ---*/ - for(auto FI =localFaces.begin(); FI!=localFaces.end(); ++FI) { + for (auto FI = localFaces.begin(); FI != localFaces.end(); ++FI) { /* Check for a matching face. */ if (FI->elemID1 < Global_nElem) { - /*--- Determine the time level of Elem0, which is always owned. ---*/ const unsigned long elemID0 = FI->elemID0 - elemPartitioner.GetFirstIndexOnRank(rank); const unsigned short timeLevel0 = elem[elemID0]->GetTimeLevel(); @@ -836,7 +826,7 @@ void CPhysicalGeometry::DetermineOwnershipInternalFaces(vector& /* Non-matching face. Give the ownership to element 0. */ FI->elem0IsOwner = true; } - } + } } void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector& localFaces) { @@ -909,8 +899,7 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector::iterator low; - low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); + auto low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); /*--- Store the relevant data in facesDonor. ---*/ facesDonor[k].nDim = nDim; @@ -921,8 +910,7 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vectorelemType0; for (unsigned short j = 0; j < nPointsPerFace[0]; ++j) { - map::const_iterator MI; - MI = globalPointIDToLocalInd.find(faceConn[0][j]); + const auto MI = globalPointIDToLocalInd.find(faceConn[0][j]); unsigned long ind = MI->second; for (unsigned l = 0; l < nDim; ++l) facesDonor[k].cornerCoor[j][l] = nodes->GetCoord(ind, l); @@ -965,7 +953,7 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector doubleLocBuf(13 * sizeLocal); unsigned long cS = 0, cL = 0, cD = 0; - for (vector::const_iterator fIt = facesDonor.begin(); fIt != facesDonor.end(); ++fIt) { + for (auto fIt = facesDonor.begin(); fIt != facesDonor.end(); ++fIt) { shortLocBuf[cS++] = fIt->nCornerPoints; shortLocBuf[cS++] = fIt->nDim; shortLocBuf[cS++] = fIt->nPoly; @@ -1079,8 +1067,7 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector::iterator low; - low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); + auto low = lower_bound(localFaces.begin(), localFaces.end(), thisFace); /*--- Indicate that this face is a periodic face. This is accomplished by setting periodicIndex to iMarker + 1. ---*/ @@ -1098,8 +1085,7 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vectorelemType0; for (unsigned short j = 0; j < nPointsPerFace[0]; ++j) { - map::const_iterator MI; - MI = globalPointIDToLocalInd.find(faceConn[0][j]); + const auto MI = globalPointIDToLocalInd.find(faceConn[0][j]); unsigned long ind = MI->second; const su2double* coor = nodes->GetCoord(ind); @@ -1120,8 +1106,7 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector::const_iterator donorLow; - donorLow = lower_bound(facesDonor.begin(), facesDonor.end(), thisMatchingFace); + const auto donorLow = lower_bound(facesDonor.begin(), facesDonor.end(), thisMatchingFace); if (donorLow != facesDonor.end()) { if (!(thisMatchingFace < *donorLow)) { @@ -1139,13 +1124,11 @@ void CPhysicalGeometry::DeterminePeriodicFacesFEMGrid(CConfig* config, vector > &adjacency, - vector &vwgt, - vector > &adjwgt) { - +void CPhysicalGeometry::DetermineFEMColorsViaParMETIS(vector >& adjacency, + vector& vwgt, + vector >& adjwgt) { /*--- Only call ParMETIS if we have more than one rank to avoid errors ---*/ if (size > SINGLE_NODE) { - /*--- Define the linear partitioning of the elements. ---*/ CLinearPartitioner elemPartitioner(Global_nElem, 0); @@ -1265,7 +1248,7 @@ void CPhysicalGeometry::DetermineFEMConstantJacobiansAndLenScale(CConfig* config for (unsigned short j = 0; j < nDOFs; ++j) { unsigned long nodeID = elem[i]->GetNode(j); - map::const_iterator MI = globalPointIDToLocalInd.find(nodeID); + const auto MI = globalPointIDToLocalInd.find(nodeID); unsigned long ind = MI->second; for (unsigned short k = 0; k < nDim; ++k, ++jj) vecRHS[jj] = nodes->GetCoord(ind, k); } @@ -1614,8 +1597,7 @@ void CPhysicalGeometry::DetermineDonorElementsWallFunctions(CConfig* config) { for (unsigned short k = 0; k < nDOFsPerSubElem[i]; ++k, ++kk) { unsigned long nodeID = elem[l]->GetNode(connSubElems[i][kk]); - map::const_iterator MI; - MI = globalPointIDToLocalInd.find(nodeID); + const auto MI = globalPointIDToLocalInd.find(nodeID); elemConn.push_back(MI->second); } @@ -1771,8 +1753,7 @@ void CPhysicalGeometry::DetermineDonorElementsWallFunctions(CConfig* config) { ii = 0; for (unsigned short j = 0; j < nDOFs; ++j) { unsigned long nodeID = bound[iMarker][l]->GetNode(j); - map::const_iterator MI; - MI = globalPointIDToLocalInd.find(nodeID); + const auto MI = globalPointIDToLocalInd.find(nodeID); nodeID = MI->second; for (unsigned short k = 0; k < nDim; ++k, ++ii) coorBoundFace[ii] = nodes->GetCoord(nodeID, k); } @@ -1902,8 +1883,7 @@ void CPhysicalGeometry::DetermineDonorElementsWallFunctions(CConfig* config) { /* Sort donorElementsFace in increasing order and remove the the double entities. */ sort(donorElementsFace.begin(), donorElementsFace.end()); - vector::iterator lastEntry; - lastEntry = unique(donorElementsFace.begin(), donorElementsFace.end()); + auto lastEntry = unique(donorElementsFace.begin(), donorElementsFace.end()); donorElementsFace.erase(lastEntry, donorElementsFace.end()); /* Store the donor elements in the data structure for @@ -2128,8 +2108,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector if (UMI == Global_to_Local_Elem.end()) { /* This element is an external element. Store it in the map mapExternalElemIDToTimeLevel if not already done so. */ - map::iterator MI; - MI = mapExternalElemIDToTimeLevel.find(FI->elemID1); + auto MI = mapExternalElemIDToTimeLevel.find(FI->elemID1); if (MI == mapExternalElemIDToTimeLevel.end()) mapExternalElemIDToTimeLevel[FI->elemID1] = CUnsignedShort2T(0, 0); } @@ -2220,8 +2199,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector /* Loop over the entries of recvBuf and add them to mapExternalElemIDToTimeLevel, if not present already. */ for (int j = 0; j < sizeMess; ++j) { - map::iterator MI; - MI = mapExternalElemIDToTimeLevel.find(recvBuf[j]); + auto MI = mapExternalElemIDToTimeLevel.find(recvBuf[j]); if (MI == mapExternalElemIDToTimeLevel.end()) mapExternalElemIDToTimeLevel[recvBuf[j]] = CUnsignedShort2T(0, 0); } } @@ -2319,14 +2297,13 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector /*--- of the external element data. ---*/ /*--------------------------------------------------------------------------*/ - map::iterator MI; #ifdef HAVE_MPI /*--- Determine the ranks from which I receive element data during the actual exchange. ---*/ recvFromRank.assign(size, 0); - for (MI = mapExternalElemIDToTimeLevel.begin(); MI != mapExternalElemIDToTimeLevel.end(); ++MI) { + for (auto MI = mapExternalElemIDToTimeLevel.begin(); MI != mapExternalElemIDToTimeLevel.end(); ++MI) { /*--- Determine the rank where this external is stored and set the corresponding index of recvFromRank to 1. ---*/ const unsigned long rankElem = elemPartitioner.GetRankContainingIndex(MI->first); @@ -2350,7 +2327,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector will be received from other ranks. ---*/ vector > recvElem(nRankRecv, vector(0)); - for (MI = mapExternalElemIDToTimeLevel.begin(); MI != mapExternalElemIDToTimeLevel.end(); ++MI) { + for (auto MI = mapExternalElemIDToTimeLevel.begin(); MI != mapExternalElemIDToTimeLevel.end(); ++MI) { const unsigned long elemID = MI->first; const unsigned long rankElem = elemPartitioner.GetRankContainingIndex(elemID); @@ -2362,7 +2339,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector exchange and send over the global element ID's. In order to avoid unnecessary communication, multiple entries are filtered out. ---*/ sendReqs.resize(nRankRecv); - map::const_iterator MRI = mapRankToIndRecv.begin(); + auto MRI = mapRankToIndRecv.begin(); for (int i = 0; i < nRankRecv; ++i, ++MRI) { sort(recvElem[i].begin(), recvElem[i].end()); @@ -2440,7 +2417,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector &status); for (unsigned long j = 0; j < recvElem[i].size(); ++j) { - MI = mapExternalElemIDToTimeLevel.find(recvElem[i][j]); + auto MI = mapExternalElemIDToTimeLevel.find(recvElem[i][j]); if (MI == mapExternalElemIDToTimeLevel.end()) SU2_MPI::Error("Entry not found in mapExternalElemIDToTimeLevel", CURRENT_FUNCTION); MI->second.short0 = returnBuf[i][2 * j]; @@ -2517,7 +2494,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector Retrieve its time level from mapExternalElemIDToTimeLevel and determine the minimum time level. */ const unsigned short timeLevelB = elem[elemID]->GetTimeLevel(); - MI = mapExternalElemIDToTimeLevel.find(donors[i]); + auto MI = mapExternalElemIDToTimeLevel.find(donors[i]); if (MI == mapExternalElemIDToTimeLevel.end()) SU2_MPI::Error("Entry not found in mapExternalElemIDToTimeLevel", CURRENT_FUNCTION); const unsigned short timeLevel = min(timeLevelB, MI->second.short0); @@ -2574,7 +2551,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector /* The second element is stored on a different processor. Retrieve its time level from mapExternalElemIDToTimeLevel and determine the minimum time level. */ - MI = mapExternalElemIDToTimeLevel.find(FI->elemID1); + auto MI = mapExternalElemIDToTimeLevel.find(FI->elemID1); if (MI == mapExternalElemIDToTimeLevel.end()) SU2_MPI::Error("Entry not found in mapExternalElemIDToTimeLevel", CURRENT_FUNCTION); const unsigned short timeLevel = min(timeLevel0, MI->second.short0); @@ -2635,7 +2612,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector &status); for (unsigned long j = 0; j < recvElem[i].size(); ++j) { - MI = mapExternalElemIDToTimeLevel.find(recvElem[i][j]); + auto MI = mapExternalElemIDToTimeLevel.find(recvElem[i][j]); if (MI == mapExternalElemIDToTimeLevel.end()) SU2_MPI::Error("Entry not found in mapExternalElemIDToTimeLevel", CURRENT_FUNCTION); MI->second.short0 = min(returnBuf[i][j], MI->second.short0); @@ -2698,14 +2675,12 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector } } -void CPhysicalGeometry::StoreFaceInfoInLocalElements(const vector &localFaces) { - +void CPhysicalGeometry::StoreFaceInfoInLocalElements(const vector& localFaces) { /*--- Define the linear partitioning of the elements. ---*/ CLinearPartitioner elemPartitioner(Global_nElem, 0); /*--- Loop over the locally stored elements. ---*/ for (unsigned long k = 0; k < nElem; ++k) { - /*--- Get the number of faces and the corner points of these faces of this element. ---*/ unsigned short nFaces; @@ -2718,7 +2693,6 @@ void CPhysicalGeometry::StoreFaceInfoInLocalElements(const vector& localFaces, - const vector >& adjacency, - const map& mapExternalElemIDToTimeLevel, - vector& vwgt, vector >& adjwgt) { +void CPhysicalGeometry::DetermineFEMGraphWeights( + CConfig* config, const vector& localFaces, const vector >& adjacency, + const map& mapExternalElemIDToTimeLevel, vector& vwgt, + vector >& adjwgt) { /*--- Define the linear partitioning of the elements. ---*/ CLinearPartitioner elemPartitioner(Global_nElem, 0); @@ -2766,11 +2740,10 @@ void CPhysicalGeometry::DetermineFEMGraphWeights(CConfig* config, const vector Date: Wed, 12 Feb 2025 16:50:33 +0100 Subject: [PATCH 4/9] Removed the option to use LIBXSMM. It is outdated --- .../include/linear_algebra/blas_structure.hpp | 7 +----- Common/src/linear_algebra/blas_structure.cpp | 22 +++++-------------- SU2_CFD/src/SU2_CFD.cpp | 15 ------------- 3 files changed, 6 insertions(+), 38 deletions(-) diff --git a/Common/include/linear_algebra/blas_structure.hpp b/Common/include/linear_algebra/blas_structure.hpp index cf4ec42ef2a..0097847c4fd 100644 --- a/Common/include/linear_algebra/blas_structure.hpp +++ b/Common/include/linear_algebra/blas_structure.hpp @@ -29,11 +29,6 @@ #pragma once -/* LIBXSMM include files, if supported. */ -#ifdef HAVE_LIBXSMM -#include "libxsmm.h" -#endif - class CConfig; /*! @@ -494,7 +489,7 @@ class CBlasStructure { } private: -#if !(defined(HAVE_LIBXSMM) || defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ +#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) /* Blocking parameters for the outer kernel. We multiply mc x kc blocks of the matrix A with kc x nc panels of the matrix B (this approach is referred diff --git a/Common/src/linear_algebra/blas_structure.cpp b/Common/src/linear_algebra/blas_structure.cpp index 6548888bba9..eab96af3c40 100644 --- a/Common/src/linear_algebra/blas_structure.cpp +++ b/Common/src/linear_algebra/blas_structure.cpp @@ -43,7 +43,7 @@ extern "C" void dgemv_(char*, const int*, const int*, const passivedouble*, cons /* Constructor. Initialize the const member variables, if needed. */ CBlasStructure::CBlasStructure() -#if !(defined(HAVE_LIBXSMM) || defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ +#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) : mc(256), kc(128), @@ -62,37 +62,25 @@ void CBlasStructure::gemm(const int M, const int N, const int K, const su2double #endif #if (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) || \ - !(defined(HAVE_LIBXSMM) || defined(HAVE_MKL) || defined(HAVE_BLAS)) + !(defined(HAVE_MKL) || defined(HAVE_BLAS)) /* Native implementation of the matrix product. This optimized implementation assumes that the matrices are in column major order. This can be accomplished by swapping N and M and A and B. This implementation is based on https://github.com/flame/how-to-optimize-gemm. */ gemm_imp(N, M, K, B, A, C); -#else -#ifdef HAVE_LIBXSMM - - /* The gemm function of libxsmm is used to carry out the multiplication. - Note that libxsmm_gemm expects the matrices in column major order. That's - why the in the calling sequence A and B and M and N are reversed. */ - su2double alpha = 1.0; - su2double beta = 0.0; - char trans = 'N'; - - libxsmm_dgemm(&trans, &trans, &N, &M, &K, &alpha, B, &N, A, &K, &beta, C, &N); - #else // MKL and BLAS /* The standard blas routine dgemm is used for the multiplication. Call dgemm without transposing the matrices. In that case dgemm expects - the matrices in column major order, see the comments for libxsmm. */ + the matrices in column major order. That's why the in the calling + sequence A and B and M and N are reversed. */ su2double alpha = 1.0; su2double beta = 0.0; char trans = 'N'; dgemm_(&trans, &trans, &N, &M, &K, &alpha, B, &N, A, &K, &beta, C, &N); -#endif #endif /* Store the profiling information, if needed. */ @@ -131,7 +119,7 @@ void CBlasStructure::gemv(const int M, const int N, const su2double* A, const su #endif } -#if !(defined(HAVE_LIBXSMM) || defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ +#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) /* Macros for accessing submatrices of a matmul using the leading dimension. */ diff --git a/SU2_CFD/src/SU2_CFD.cpp b/SU2_CFD/src/SU2_CFD.cpp index 33bb1a93d15..148071f5694 100644 --- a/SU2_CFD/src/SU2_CFD.cpp +++ b/SU2_CFD/src/SU2_CFD.cpp @@ -27,11 +27,6 @@ #include "../include/SU2_CFD.hpp" -/* LIBXSMM include files, if supported. */ -#ifdef HAVE_LIBXSMM -#include "libxsmm.h" -#endif - /* Include file, needed for the runtime NaN catching. You also have to include feenableexcept(...) below. */ //#include @@ -77,11 +72,6 @@ int main(int argc, char *argv[]) { /*--- Uncomment the following line if runtime NaN catching is desired. ---*/ // feenableexcept(FE_INVALID | FE_OVERFLOW | FE_DIVBYZERO ); - /*--- Initialize libxsmm, if supported. ---*/ -#ifdef HAVE_LIBXSMM - libxsmm_init(); -#endif - /*--- Create a pointer to the main SU2 Driver ---*/ CDriver* driver = nullptr; @@ -154,11 +144,6 @@ int main(int argc, char *argv[]) { delete driver; - /*---Finalize libxsmm, if supported. ---*/ -#ifdef HAVE_LIBXSMM - libxsmm_finalize(); -#endif - /*--- Finalize AD. ---*/ AD::Finalize(); From 2fc8d69dd5f403a291e4f6a2cd53dbd081b44b93 Mon Sep 17 00:00:00 2001 From: vdweide Date: Wed, 12 Feb 2025 17:16:12 +0100 Subject: [PATCH 5/9] Some cosmetic changes --- Common/include/linear_algebra/blas_structure.hpp | 3 +-- Common/src/linear_algebra/blas_structure.cpp | 9 +++------ 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/Common/include/linear_algebra/blas_structure.hpp b/Common/include/linear_algebra/blas_structure.hpp index 0097847c4fd..14d42cf386a 100644 --- a/Common/include/linear_algebra/blas_structure.hpp +++ b/Common/include/linear_algebra/blas_structure.hpp @@ -489,8 +489,7 @@ class CBlasStructure { } private: -#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ - (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) +#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) /* Blocking parameters for the outer kernel. We multiply mc x kc blocks of the matrix A with kc x nc panels of the matrix B (this approach is referred to as `gebp` in the literature). */ diff --git a/Common/src/linear_algebra/blas_structure.cpp b/Common/src/linear_algebra/blas_structure.cpp index eab96af3c40..4f52fe7b8fe 100644 --- a/Common/src/linear_algebra/blas_structure.cpp +++ b/Common/src/linear_algebra/blas_structure.cpp @@ -43,8 +43,7 @@ extern "C" void dgemv_(char*, const int*, const int*, const passivedouble*, cons /* Constructor. Initialize the const member variables, if needed. */ CBlasStructure::CBlasStructure() -#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ - (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) +#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) : mc(256), kc(128), nc(128) @@ -61,8 +60,7 @@ void CBlasStructure::gemm(const int M, const int N, const int K, const su2double if (config) config->GEMM_Tick(&timeGemm); #endif -#if (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) || \ - !(defined(HAVE_MKL) || defined(HAVE_BLAS)) +#if (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) || !(defined(HAVE_MKL) || defined(HAVE_BLAS)) /* Native implementation of the matrix product. This optimized implementation assumes that the matrices are in column major order. This can be accomplished by swapping N and M and A and B. This implementation is based @@ -119,8 +117,7 @@ void CBlasStructure::gemv(const int M, const int N, const su2double* A, const su #endif } -#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || \ - (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) +#if !(defined(HAVE_BLAS) || defined(HAVE_MKL)) || (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) /* Macros for accessing submatrices of a matmul using the leading dimension. */ #define A(i, j) a[(j)*lda + (i)] From 7f56cb14f008efe68bc184047d5d7475651154ed Mon Sep 17 00:00:00 2001 From: vdweide Date: Thu, 13 Feb 2025 08:15:00 +0100 Subject: [PATCH 6/9] Changed the su2doubles to passivedoubles in the implementation of DetermineFEMGraphWeights --- Common/src/geometry/CPhysicalGeometryFEM.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Common/src/geometry/CPhysicalGeometryFEM.cpp b/Common/src/geometry/CPhysicalGeometryFEM.cpp index f8c18b66dea..68de215fa76 100644 --- a/Common/src/geometry/CPhysicalGeometryFEM.cpp +++ b/Common/src/geometry/CPhysicalGeometryFEM.cpp @@ -2725,8 +2725,8 @@ void CPhysicalGeometry::StoreFaceInfoInLocalElements(const vector& localFaces, const vector >& adjacency, - const map& mapExternalElemIDToTimeLevel, vector& vwgt, - vector >& adjwgt) { + const map& mapExternalElemIDToTimeLevel, vector& vwgt, + vector >& adjwgt) { /*--- Define the linear partitioning of the elements. ---*/ CLinearPartitioner elemPartitioner(Global_nElem, 0); From c7dc04da2467ef6080ca2b795ec588848c5d9ba9 Mon Sep 17 00:00:00 2001 From: vdweide Date: Thu, 13 Feb 2025 09:47:55 +0100 Subject: [PATCH 7/9] Changed the work estimates for ParMetis to passivedouble --- Common/include/fem/fem_standard_element.hpp | 8 ++++---- Common/src/fem/fem_work_estimate_metis.cpp | 8 ++++---- Common/src/geometry/CPhysicalGeometryFEM.cpp | 11 ++++++----- 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/Common/include/fem/fem_standard_element.hpp b/Common/include/fem/fem_standard_element.hpp index b074dac7d0b..244a394ee54 100644 --- a/Common/include/fem/fem_standard_element.hpp +++ b/Common/include/fem/fem_standard_element.hpp @@ -990,7 +990,7 @@ class CFEMStandardElement : public CFEMStandardElementBase { type. This information is used to determine a well balanced partition. * \param[in] config - Object, which contains the input parameters. */ - su2double WorkEstimateMetis(CConfig* config); + passivedouble WorkEstimateMetis(CConfig* config); private: /*! @@ -1420,7 +1420,7 @@ class CFEMStandardInternalFace : public CFEMStandardElementBase { type. This information is used to determine a well balanced partition. * \param[in] config - Object, which contains the input parameters. */ - su2double WorkEstimateMetis(CConfig* config); + passivedouble WorkEstimateMetis(CConfig* config); private: /*! @@ -1636,7 +1636,7 @@ class CFEMStandardBoundaryFace : public CFEMStandardElementBase { type. This information is used to determine a well balanced partition. * \param[in] config - Object, which contains the input parameters. */ - su2double WorkEstimateMetis(CConfig* config); + passivedouble WorkEstimateMetis(CConfig* config); /*! * \brief Function, which estimates the additional amount of work for an element @@ -1645,7 +1645,7 @@ class CFEMStandardBoundaryFace : public CFEMStandardElementBase { * \param[in] config - Object, which contains the input parameters. * \param[in] nPointsWF - Number of points to discretize the wall model. */ - su2double WorkEstimateMetisWallFunctions(CConfig* config, const unsigned short nPointsWF); + passivedouble WorkEstimateMetisWallFunctions(CConfig* config, const unsigned short nPointsWF); private: /*! diff --git a/Common/src/fem/fem_work_estimate_metis.cpp b/Common/src/fem/fem_work_estimate_metis.cpp index 02841cc459a..787fe613573 100644 --- a/Common/src/fem/fem_work_estimate_metis.cpp +++ b/Common/src/fem/fem_work_estimate_metis.cpp @@ -28,22 +28,22 @@ #include "../../include/fem/fem_standard_element.hpp" -su2double CFEMStandardElement::WorkEstimateMetis(CConfig* config) { +passivedouble CFEMStandardElement::WorkEstimateMetis(CConfig* config) { /* TEMPORARY IMPLEMENTATION. */ return nIntegration + 0.1 * nDOFs; } -su2double CFEMStandardInternalFace::WorkEstimateMetis(CConfig* config) { +passivedouble CFEMStandardInternalFace::WorkEstimateMetis(CConfig* config) { /* TEMPORARY IMPLEMENTATION. */ return 2.0 * nIntegration + 0.05 * (nDOFsFaceSide0 + nDOFsFaceSide1); } -su2double CFEMStandardBoundaryFace::WorkEstimateMetis(CConfig* config) { +passivedouble CFEMStandardBoundaryFace::WorkEstimateMetis(CConfig* config) { /* TEMPORARY IMPLEMENTATION. */ return nIntegration + 0.05 * nDOFsFace; } -su2double CFEMStandardBoundaryFace::WorkEstimateMetisWallFunctions(CConfig* config, const unsigned short nPointsWF) { +passivedouble CFEMStandardBoundaryFace::WorkEstimateMetisWallFunctions(CConfig* config, const unsigned short nPointsWF) { /* TEMPORARY IMPLEMENTATION. */ return 0.25 * nIntegration * nPointsWF; } diff --git a/Common/src/geometry/CPhysicalGeometryFEM.cpp b/Common/src/geometry/CPhysicalGeometryFEM.cpp index 68de215fa76..5b1b5a5fea1 100644 --- a/Common/src/geometry/CPhysicalGeometryFEM.cpp +++ b/Common/src/geometry/CPhysicalGeometryFEM.cpp @@ -1162,13 +1162,13 @@ void CPhysicalGeometry::DetermineFEMColorsViaParMETIS(vector vwgtPar(nElem * ncon); - for (unsigned long i = 0; i < nElem * ncon; ++i) vwgtPar[i] = (idx_t)ceil(vwgt[i]); + for (unsigned long i = 0; i < nElem * ncon; ++i) vwgtPar[i] = static_cast(ceil(vwgt[i])); /* Create the adjacency weight in ParMETIS format. */ vector adjwgtPar(xadjPar[nElem]); ii = 0; for (unsigned long i = 0; i < nElem; ++i) { - for (unsigned long j = 0; j < adjwgt[i].size(); ++j, ++ii) adjwgtPar[ii] = (idx_t)ceil(adjwgt[i][j]); + for (unsigned long j = 0; j < adjwgt[i].size(); ++j, ++ii) adjwgtPar[ii] = static_cast(ceil(adjwgt[i][j])); } /* Make sure that an equal distribution is obtained. */ @@ -2960,12 +2960,13 @@ void CPhysicalGeometry::DetermineFEMGraphWeights( /*--- Determine the minimum of the workload of the elements, i.e. 1st vertex weight, over the entire domain. ---*/ - su2double minvwgt = vwgt[0]; + passivedouble minvwgt = vwgt[0]; for (unsigned long i = 0; i < nElem; ++i) minvwgt = min(minvwgt, vwgt[2 * i]); #ifdef HAVE_MPI - su2double locminvwgt = minvwgt; - SU2_MPI::Allreduce(&locminvwgt, &minvwgt, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); + su2double locminvwgt = minvwgt, globminvwgt; + SU2_MPI::Allreduce(&locminvwgt, &globminvwgt, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm()); + minvwgt = SU2_TYPE::GetValue(globminvwgt); #endif /*--- Scale the workload of the elements, the 1st vertex weight, with the From 31f6ed81bc89df4112cfa7e6046c6e602b1f4495 Mon Sep 17 00:00:00 2001 From: vdweide Date: Thu, 13 Feb 2025 09:49:01 +0100 Subject: [PATCH 8/9] Forgot the clang-format --- Common/src/fem/fem_work_estimate_metis.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Common/src/fem/fem_work_estimate_metis.cpp b/Common/src/fem/fem_work_estimate_metis.cpp index 787fe613573..d9161e15365 100644 --- a/Common/src/fem/fem_work_estimate_metis.cpp +++ b/Common/src/fem/fem_work_estimate_metis.cpp @@ -43,7 +43,8 @@ passivedouble CFEMStandardBoundaryFace::WorkEstimateMetis(CConfig* config) { return nIntegration + 0.05 * nDOFsFace; } -passivedouble CFEMStandardBoundaryFace::WorkEstimateMetisWallFunctions(CConfig* config, const unsigned short nPointsWF) { +passivedouble CFEMStandardBoundaryFace::WorkEstimateMetisWallFunctions(CConfig* config, + const unsigned short nPointsWF) { /* TEMPORARY IMPLEMENTATION. */ return 0.25 * nIntegration * nPointsWF; } From 23c3d09d562fc33bfdd78cc0fd125ba6fa71dbdf Mon Sep 17 00:00:00 2001 From: vdweide Date: Sun, 16 Feb 2025 09:01:45 +0100 Subject: [PATCH 9/9] Addressing comments --- .../include/toolboxes/fem/CFaceOfElement.hpp | 9 ++-- .../include/toolboxes/fem/CMatchingFace.hpp | 7 +-- .../toolboxes/fem/CReorderElements.hpp | 10 ++-- .../toolboxes/fem/CSortBoundaryFaces.hpp | 48 ------------------- Common/include/toolboxes/fem/CSortFaces.hpp | 2 +- Common/src/fem/fem_geometry_structure.cpp | 14 +++++- Common/src/toolboxes/fem/CFaceOfElement.cpp | 20 ++++---- .../src/toolboxes/fem/CSortBoundaryFaces.cpp | 42 ---------------- Common/src/toolboxes/fem/meson.build | 1 - 9 files changed, 33 insertions(+), 120 deletions(-) delete mode 100644 Common/include/toolboxes/fem/CSortBoundaryFaces.hpp delete mode 100644 Common/src/toolboxes/fem/CSortBoundaryFaces.cpp diff --git a/Common/include/toolboxes/fem/CFaceOfElement.hpp b/Common/include/toolboxes/fem/CFaceOfElement.hpp index 57f57ac7110..e1f17f67c2b 100644 --- a/Common/include/toolboxes/fem/CFaceOfElement.hpp +++ b/Common/include/toolboxes/fem/CFaceOfElement.hpp @@ -34,8 +34,6 @@ #include #include -using namespace std; - /*! * \class CFaceOfElement * \brief Class used in the partitioning of the FEM grid as well as the building of @@ -68,9 +66,8 @@ class CFaceOfElement { bool elem0IsOwner; /*!< \brief Whether or not the neighboring element 0 is the owner of the face. If false, element 1 is the owner. */ - /* Standard constructor and destructor. */ + /* Constructor. Initialize the member variables. */ CFaceOfElement(); - ~CFaceOfElement() {} /* Alternative constructor to set the corner points. */ CFaceOfElement(const unsigned short VTK_Type, const unsigned short nPoly, const unsigned long* Nodes); @@ -91,11 +88,11 @@ class CFaceOfElement { /*--- Member function, which creates a unique numbering for the corner points. A sort in increasing order is OK for this purpose. ---*/ - inline void CreateUniqueNumbering(void) { std::sort(cornerPoints, cornerPoints + nCornerPoints); } + inline void CreateUniqueNumbering() { std::sort(cornerPoints, cornerPoints + nCornerPoints); } /*--- Member function, which creates a unique numbering for the corner points while the orientation is taken into account. ---*/ - void CreateUniqueNumberingWithOrientation(void); + void CreateUniqueNumberingWithOrientation(); private: /*--- Copy function, which copies the data of the given object into the current object. ---*/ diff --git a/Common/include/toolboxes/fem/CMatchingFace.hpp b/Common/include/toolboxes/fem/CMatchingFace.hpp index 9f68bbddd5a..c8a91d188bc 100644 --- a/Common/include/toolboxes/fem/CMatchingFace.hpp +++ b/Common/include/toolboxes/fem/CMatchingFace.hpp @@ -50,12 +50,9 @@ class CMatchingFace { su2double cornerCoor[4][3]; /*!< \brief Coordinates of the corner points of the face. */ su2double tolForMatching; /*!< \brief Tolerance for this face for matching points. */ - /* Standard constructor. */ + /* Constructor. Initialize the member variables to zero. */ CMatchingFace(); - /* Destructor, nothing to be done. */ - ~CMatchingFace() {} - /* Copy constructor and assignment operator. */ inline CMatchingFace(const CMatchingFace& other) { Copy(other); } @@ -68,7 +65,7 @@ class CMatchingFace { bool operator<(const CMatchingFace& other) const; /*--- Member function, which sorts the coordinates of the face. ---*/ - void SortFaceCoordinates(void); + void SortFaceCoordinates(); private: /*--- Copy function, which copies the data of the given object into the current object. ---*/ diff --git a/Common/include/toolboxes/fem/CReorderElements.hpp b/Common/include/toolboxes/fem/CReorderElements.hpp index ac2cd6c6c27..ad66453f443 100644 --- a/Common/include/toolboxes/fem/CReorderElements.hpp +++ b/Common/include/toolboxes/fem/CReorderElements.hpp @@ -60,7 +60,7 @@ class CReorderElements { /*! * \brief Default constructor of the class. Disabled. */ - CReorderElements(void) = delete; + CReorderElements() = delete; /*! * \brief Less than operator of the class. Needed for the sorting. @@ -71,26 +71,26 @@ class CReorderElements { * \brief Function to make available the variable commSolution. * \return Whether or not the solution of the element must be communicated. */ - inline bool GetCommSolution(void) const { return commSolution; } + inline bool GetCommSolution() const { return commSolution; } /*! * \brief Function to make available the element type of the element. * \return The value of elemType, which stores the VTK type, polynomial degree and whether or not the Jacobian is constant. */ - inline unsigned short GetElemType(void) const { return elemType; } + inline unsigned short GetElemType() const { return elemType; } /*! * \brief Function to make available the global element ID. * \return The global element ID of the element. */ - inline unsigned long GetGlobalElemID(void) const { return globalElemID; } + inline unsigned long GetGlobalElemID() const { return globalElemID; } /*! * \brief Function to make available the time level. * \return The time level of the element. */ - inline unsigned short GetTimeLevel(void) const { return timeLevel; } + inline unsigned short GetTimeLevel() const { return timeLevel; } /*! * \brief Function, which sets the value of commSolution. diff --git a/Common/include/toolboxes/fem/CSortBoundaryFaces.hpp b/Common/include/toolboxes/fem/CSortBoundaryFaces.hpp deleted file mode 100644 index 8271a084a99..00000000000 --- a/Common/include/toolboxes/fem/CSortBoundaryFaces.hpp +++ /dev/null @@ -1,48 +0,0 @@ -/*! - * \file CSortBoundaryFaces.hpp - * \brief Header file for the class CSortBoundaryFaces. - * The implementations are in the CSortBoundaryFaces.cpp file. - * \author E. van der Weide - * \version 8.1.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2024, 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 . - */ - -#pragma once - -#include - -/*! - * \class CSortBoundaryFaces - * \brief Functor, used for a different sorting of the faces than the < operator - * of CSurfaceElementFEM. - * \author E. van der Weide - * \version 8.1.0 "Harrier" - */ -struct CSurfaceElementFEM; // Forward declaration to avoid problems. -struct CSortBoundaryFaces { - /*! - * \brief Operator used for the comparison. - * \param[in] f0 - First boundary face in the comparison. - * \param[in] f1 - Second boundary face in the comparison. - */ - bool operator()(const CSurfaceElementFEM& f0, const CSurfaceElementFEM& f1); -}; diff --git a/Common/include/toolboxes/fem/CSortFaces.hpp b/Common/include/toolboxes/fem/CSortFaces.hpp index 721472b60a0..71d98b62661 100644 --- a/Common/include/toolboxes/fem/CSortFaces.hpp +++ b/Common/include/toolboxes/fem/CSortFaces.hpp @@ -60,7 +60,7 @@ class CSortFaces { /*! * \brief Default constructor of the class. Disabled. */ - CSortFaces(void) = delete; + CSortFaces() = delete; /*! * \brief Operator used for the comparison. diff --git a/Common/src/fem/fem_geometry_structure.cpp b/Common/src/fem/fem_geometry_structure.cpp index 183d65d28ef..0c2468e806d 100644 --- a/Common/src/fem/fem_geometry_structure.cpp +++ b/Common/src/fem/fem_geometry_structure.cpp @@ -29,7 +29,6 @@ #include "../../include/toolboxes/classes_multiple_integers.hpp" #include "../../include/toolboxes/fem/CReorderElements.hpp" #include "../../include/toolboxes/fem/CSortFaces.hpp" -#include "../../include/toolboxes/fem/CSortBoundaryFaces.hpp" #include "../../include/fem/fem_geometry_structure.hpp" #include "../../include/geometry/primal_grid/CPrimalGridFEM.hpp" #include "../../include/geometry/primal_grid/CPrimalGridBoundFEM.hpp" @@ -2974,7 +2973,18 @@ void CMeshFEM_DG::CreateFaces(CConfig* config) { each time level. */ for (unsigned short i = 0; i < nTimeLevels; ++i) sort(surfElem.begin() + boundaries[iMarker].nSurfElem[i], - surfElem.begin() + boundaries[iMarker].nSurfElem[i + 1], CSortBoundaryFaces()); + surfElem.begin() + boundaries[iMarker].nSurfElem[i + 1], + [](const CSurfaceElementFEM& f0, const CSurfaceElementFEM& f1) { + /* First sorting criterion is the index of the standard element. The + boundary faces should be sorted per standard element. Note that the + time level is not taken into account here, because it is assumed that + the surface elements to be sorted belong to one time level. */ + if (f0.indStandardElement != f1.indStandardElement) return f0.indStandardElement < f1.indStandardElement; + + /* The standard elements are the same. The second criterion is the + corresponding volume IDs of the surface elements. */ + return f0.volElemID < f1.volElemID; + }); } } } diff --git a/Common/src/toolboxes/fem/CFaceOfElement.cpp b/Common/src/toolboxes/fem/CFaceOfElement.cpp index cb46aad39a6..439e2985ab8 100644 --- a/Common/src/toolboxes/fem/CFaceOfElement.cpp +++ b/Common/src/toolboxes/fem/CFaceOfElement.cpp @@ -90,7 +90,7 @@ CFaceOfElement::CFaceOfElement(const unsigned short VTK_Type, const unsigned sho } default: { - ostringstream message; + std::ostringstream message; message << "Unknown VTK surface element type, " << VTK_Type; SU2_MPI::Error(message.str(), CURRENT_FUNCTION); } @@ -151,7 +151,7 @@ void CFaceOfElement::CreateUniqueNumberingWithOrientation() { also the element information must be swapped, because element 0 is to the left of the face and element 1 to the right. */ if (cornerPoints[1] < cornerPoints[0]) { - swap(cornerPoints[0], cornerPoints[1]); + std::swap(cornerPoints[0], cornerPoints[1]); swapElements = true; } break; @@ -225,19 +225,19 @@ void CFaceOfElement::CreateUniqueNumberingWithOrientation() { } default: { - ostringstream message; - message << "Unknown surface element type with " << nCornerPoints << " corners." << endl; + std::ostringstream message; + message << "Unknown surface element type with " << nCornerPoints << " corners." << std::endl; SU2_MPI::Error(message.str(), CURRENT_FUNCTION); } } /* Swap the element information, if needed. */ if (swapElements) { - swap(elemID0, elemID1); - swap(nPolyGrid0, nPolyGrid1); - swap(nPolySol0, nPolySol1); - swap(nDOFsElem0, nDOFsElem1); - swap(elemType0, elemType1); - swap(faceID0, faceID1); + std::swap(elemID0, elemID1); + std::swap(nPolyGrid0, nPolyGrid1); + std::swap(nPolySol0, nPolySol1); + std::swap(nDOFsElem0, nDOFsElem1); + std::swap(elemType0, elemType1); + std::swap(faceID0, faceID1); } } diff --git a/Common/src/toolboxes/fem/CSortBoundaryFaces.cpp b/Common/src/toolboxes/fem/CSortBoundaryFaces.cpp deleted file mode 100644 index 945b110a011..00000000000 --- a/Common/src/toolboxes/fem/CSortBoundaryFaces.cpp +++ /dev/null @@ -1,42 +0,0 @@ -/*! - * \file CSortFaces.cpp - * \brief Functor, used for a different sorting of the surface elements - * than the < operator of CSurfaceElementFEM. - * \author E. van der Weide - * \version 8.1.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2024, 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 . - */ - -#include "../../../include/toolboxes/fem/CSortBoundaryFaces.hpp" -#include "../../../include/fem/fem_geometry_structure.hpp" - -bool CSortBoundaryFaces::operator()(const CSurfaceElementFEM& f0, const CSurfaceElementFEM& f1) { - /* First sorting criterion is the index of the standard element. The - boundary faces should be sorted per standard element. Note that the - time level is not taken into account here, because it is assumed that - the surface elements to be sorted belong to one time level. */ - if (f0.indStandardElement != f1.indStandardElement) return f0.indStandardElement < f1.indStandardElement; - - /* The standard elements are the same. The second criterion is the - corresponding volume IDs of the surface elements. */ - return f0.volElemID < f1.volElemID; -} diff --git a/Common/src/toolboxes/fem/meson.build b/Common/src/toolboxes/fem/meson.build index a37f901fffb..3193907b669 100644 --- a/Common/src/toolboxes/fem/meson.build +++ b/Common/src/toolboxes/fem/meson.build @@ -1,5 +1,4 @@ common_src += files(['CFaceOfElement.cpp', 'CMatchingFace.cpp', 'CReorderElements.cpp', - 'CSortBoundaryFaces.cpp', 'CSortFaces.cpp'])