diff --git a/Makefile.am b/Makefile.am index 13745fd..709dcb3 100644 --- a/Makefile.am +++ b/Makefile.am @@ -23,6 +23,7 @@ libTempestRemap_la_include_HEADERS = \ src/TriangularQuadrature.h \ src/CommandLine.h \ src/Defines.h \ + src/FunctionTimer.h \ src/GaussLobattoQuadrature.h \ src/kdtree.h \ src/order32.h \ @@ -55,6 +56,7 @@ libTempestRemap_la_include_HEADERS = \ libTempestRemap_la_SOURCES = \ src/Announce.cpp \ src/PolynomialInterp.cpp \ + src/FunctionTimer.cpp \ src/GridElements.cpp \ src/MeshUtilities.cpp \ src/MeshUtilitiesFuzzy.cpp \ @@ -68,8 +70,9 @@ libTempestRemap_la_SOURCES = \ src/GenerateUTMMesh.cpp \ src/GenerateGLLMetaData.cpp \ src/GenerateLambertConfConicMesh.cpp \ - src/GenerateOverlapMesh.cpp \ - src/GenerateOverlapMesh_v1.cpp \ + src/GenerateOverlapMeshEdge.cpp \ + src/GenerateOverlapMeshKdx.cpp \ + src/GenerateOverlapMeshLint.cpp \ src/GaussQuadrature.cpp \ src/GaussLobattoQuadrature.cpp \ src/LegendrePolynomial.cpp \ @@ -110,12 +113,11 @@ GenerateVolumetricMesh_SOURCES = src/GenerateVolumetricMesh.cpp # Overlap mesh generators GenerateOverlapMesh_SOURCES = src/GenerateOverlapMeshExe.cpp -GenerateOverlapMesh_v1_SOURCES = src/GenerateOverlapMeshExe_v1.cpp +GenerateOverlapMeshParallel_SOURCES = src/GenerateOverlapMeshParallelExe.cpp # Compute and apply the offline mapping weights ApplyOfflineMap_SOURCES = src/ApplyOfflineMapExe.cpp GenerateOfflineMap_SOURCES = src/GenerateOfflineMapExe.cpp -# GenerateOfflineMap_v1_SOURCES = src/GenerateOfflineMap_v1.cpp GenerateGLLMetaData_SOURCES = src/GenerateGLLMetaDataExe.cpp GenerateTransposeMap_SOURCES = src/GenerateTransposeMap.cpp @@ -134,7 +136,7 @@ bin_PROGRAMS = GenerateTestData \ GenerateCSMesh GenerateTransectMesh GenerateStereographicMesh GenerateRLLMesh \ GenerateUTMMesh GenerateICOMesh GenerateRectilinearMeshFromFile \ GenerateVolumetricMesh GenerateLambertConfConicMesh \ - GenerateOverlapMesh GenerateOverlapMesh_v1 \ + GenerateOverlapMesh GenerateOverlapMeshParallel \ ApplyOfflineMap GenerateOfflineMap \ CalculateDiffNorms GenerateGLLMetaData \ GenerateTransposeMap CoarsenRectilinearData \ diff --git a/src/Defines.h b/src/Defines.h index 2e7fdce..cb2b6e7 100644 --- a/src/Defines.h +++ b/src/Defines.h @@ -68,7 +68,8 @@ static const Real ReferenceTolerance = 1.0e-12; // is the slowest). // //#define OVERLAPMESH_RETAIN_REPEATED_NODES -#define OVERLAPMESH_USE_UNSORTED_MAP +#define OVERLAPMESH_USE_SORTED_MAP +//#define OVERLAPMESH_USE_UNSORTED_MAP //#define OVERLAPMESH_USE_NODE_MULTIMAP /////////////////////////////////////////////////////////////////////////////// diff --git a/src/FiniteElementTools.cpp b/src/FiniteElementTools.cpp index ce189af..44b6cbf 100644 --- a/src/FiniteElementTools.cpp +++ b/src/FiniteElementTools.cpp @@ -294,7 +294,7 @@ double GenerateMetaData( // Verify face areas are available if (!fNoBubble) { - if (mesh.vecFaceArea.GetRows() != nElements) { + if (mesh.vecFaceArea.size() != nElements) { _EXCEPTIONT("Face area information unavailable or incorrect"); } } diff --git a/src/FunctionTimer.cpp b/src/FunctionTimer.cpp new file mode 100644 index 0000000..92b04ce --- /dev/null +++ b/src/FunctionTimer.cpp @@ -0,0 +1,199 @@ +/////////////////////////////////////////////////////////////////////////////// +/// +/// \file FunctionTimer.cpp +/// \author Paul Ullrich +/// \version November 2, 2021 +/// +/// +/// Copyright 2021 Paul Ullrich +/// +/// This file is distributed as part of the Tempest source code package. +/// Permission is granted to use, copy, modify and distribute this +/// source code and its documentation under the terms of the GNU General +/// Public License. This software is provided "as is" without express +/// or implied warranty. +/// + +#include "FunctionTimer.h" +#include "Exception.h" + +#include +#include + +/////////////////////////////////////////////////////////////////////////////// + +FunctionTimer::GroupDataMap FunctionTimer::m_mapGroupData; + +/////////////////////////////////////////////////////////////////////////////// + +FunctionTimer::FunctionTimer(const char *szGroup) { + + // Start the timer + m_fStopped = false; + + // Assign group name + if (szGroup == NULL) { + m_strGroup = ""; + } else { + m_strGroup = szGroup; + } + + // Assign start time + m_tpStartTime = std::chrono::high_resolution_clock::now(); +} + +/////////////////////////////////////////////////////////////////////////////// + +void FunctionTimer::Reset() { + m_fStopped = false; + m_tpStartTime = std::chrono::high_resolution_clock::now(); +} + +/////////////////////////////////////////////////////////////////////////////// + +std::chrono::microseconds FunctionTimer::Time(bool fDone) { + + if (!m_fStopped) { + m_tpStopTime = std::chrono::high_resolution_clock::now(); + } + + std::chrono::microseconds msTime = + std::chrono::duration_cast( + m_tpStopTime - m_tpStartTime); + + // If no name associated with this timer, ignore fDone. + if (m_strGroup == "") { + return msTime; + } + + // Add the time to the group record + if ((fDone) && (!m_fStopped)) { + m_fStopped = true; + + GroupDataMap::iterator iter; + + iter = m_mapGroupData.find(m_strGroup); + + // Add to existing group record + if (iter != m_mapGroupData.end()) { + iter->second.iTotalTime += msTime; + iter->second.nEntries++; + + // Create new group record + } else { + GroupDataPair gdp; + gdp.first = m_strGroup; + gdp.second.iTotalTime = msTime; + gdp.second.nEntries = 1; + + m_mapGroupData.insert(gdp); + } + } + + return msTime; +} + +/////////////////////////////////////////////////////////////////////////////// + +std::chrono::microseconds FunctionTimer::StopTime() { + return Time(true); +} + +/////////////////////////////////////////////////////////////////////////////// + +const FunctionTimer::TimerGroupData & FunctionTimer::GetGroupTimeRecord( + const char *szName +) { + GroupDataMap::iterator iter; + + iter = m_mapGroupData.find(szName); + + // Retrieve existing group record + if (iter != m_mapGroupData.end()) { + return iter->second; + + // Group record does not exist + } else { + _EXCEPTION1("Group time record %s does not exist.", szName); + } +} + +/////////////////////////////////////////////////////////////////////////////// + +std::chrono::microseconds FunctionTimer::GetTotalGroupTime(const char *szName) { + + GroupDataMap::iterator iter; + + iter = m_mapGroupData.find(szName); + + // Retrieve existing group record + if (iter != m_mapGroupData.end()) { + const TimerGroupData & tgd = iter->second; + + return (tgd.iTotalTime); + + // Group record does not exist + } else { + return std::chrono::microseconds(0); + } +} + +/////////////////////////////////////////////////////////////////////////////// + +std::chrono::microseconds FunctionTimer::GetAverageGroupTime(const char *szName) { + + GroupDataMap::iterator iter; + + iter = m_mapGroupData.find(szName); + + // Retrieve existing group record + if (iter != m_mapGroupData.end()) { + const TimerGroupData & tgd = iter->second; + + return (tgd.iTotalTime / tgd.nEntries); + + // Group record does not exist + } else { + return std::chrono::microseconds(0); + } +} + +/////////////////////////////////////////////////////////////////////////////// + +unsigned long FunctionTimer::GetNumberOfEntries(const char *szName) { + + GroupDataMap::iterator iter; + + iter = m_mapGroupData.find(szName); + + // Retrieve existing group record + if (iter != m_mapGroupData.end()) { + const TimerGroupData & tgd = iter->second; + + return (tgd.nEntries); + + // Group record does not exist + } else { + return 0; + } +} + +/////////////////////////////////////////////////////////////////////////////// + +void FunctionTimer::ResetGroupTimeRecord(const char *szName) { + GroupDataMap::iterator iter; + + iter = m_mapGroupData.find(szName); + + // Retrieve existing group record + if (iter != m_mapGroupData.end()) { + m_mapGroupData.erase(iter); + + // Group record does not exist + } else { + _EXCEPTION1("Group time record %s does not exist.", szName); + } +} + +/////////////////////////////////////////////////////////////////////////////// + diff --git a/src/FunctionTimer.h b/src/FunctionTimer.h new file mode 100644 index 0000000..b4ed0a2 --- /dev/null +++ b/src/FunctionTimer.h @@ -0,0 +1,143 @@ +/////////////////////////////////////////////////////////////////////////////// +/// +/// \file FunctionTimer.h +/// \author Paul Ullrich +/// \version November 2, 2021 +/// +/// +/// Copyright 2021 Paul Ullrich +/// +/// This file is distributed as part of the Tempest source code package. +/// Permission is granted to use, copy, modify and distribute this +/// source code and its documentation under the terms of the GNU General +/// Public License. This software is provided "as is" without express +/// or implied warranty. +/// + +#ifndef _FUNCTIONTIMER_H_ +#define _FUNCTIONTIMER_H_ + +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// FunctionTimer is a class used for timing operations or groups of +/// operations. Timing is provided via the Unix gettimeofday class +/// and is calculated in microseconds. +/// +class FunctionTimer { + +public: + /// + /// A structure for storing group data. + /// + struct TimerGroupData { + std::chrono::microseconds iTotalTime; + unsigned long nEntries; + }; + + /// + /// The map structure for storing group data. + /// + typedef std::map GroupDataMap; + + typedef std::pair GroupDataPair; + +public: + /// + /// Constructor. + /// + /// + /// Group name used for organizing grouped operations. + /// + FunctionTimer(const char *szGroup = NULL); + + /// + /// Destructor. + /// + virtual ~FunctionTimer() { + StopTime(); + } + + /// + /// Reset the timer. + /// + void Reset(); + + /// + /// Return the time elapsed since this timer began. + /// + /// + /// If true stores the elapsed time in the group structure. + /// + std::chrono::microseconds Time(bool fDone = false); + + /// + /// Return the time elapsed since this timer began and store in + /// group data. + /// + std::chrono::microseconds StopTime(); + +public: + /// + /// Retrieve a group data record. + /// + static const TimerGroupData & GetGroupTimeRecord(const char *szName); + + /// + /// Retrieve the total time from a group data record. + /// + static std::chrono::microseconds GetTotalGroupTime(const char *szName); + + /// + /// Retrieve the average time from a group data record. + /// + static std::chrono::microseconds GetAverageGroupTime(const char *szName); + + /// + /// Retrieve the number of entries from a group data record. + /// + static unsigned long GetNumberOfEntries(const char *szName); + + /// + /// Reset the group data record. + /// + static void ResetGroupTimeRecord(const char *szName); + +private: + /// + /// Group data. + /// + static GroupDataMap m_mapGroupData; + +private: + /// + /// Timer is stopped. + /// + bool m_fStopped; + + /// + /// Time at which this timer was constructed. + /// + /// Time at which this timer was stopped. + /// + /// Group name associated with this timer. + /// + std::string m_strGroup; +}; + +/////////////////////////////////////////////////////////////////////////////// + +#endif + diff --git a/src/GenerateOfflineMap.cpp b/src/GenerateOfflineMap.cpp index faaf3be..65b1694 100644 --- a/src/GenerateOfflineMap.cpp +++ b/src/GenerateOfflineMap.cpp @@ -344,10 +344,10 @@ try { _ASSERT(meshOverlap.vecSourceFaceIx.size() == meshOverlap.faces.size()); _ASSERT(meshOverlap.vecTargetFaceIx.size() == meshOverlap.faces.size()); - _ASSERT(meshOverlap.vecFaceArea.GetRows() == meshOverlap.faces.size()); + _ASSERT(meshOverlap.vecFaceArea.size() == meshOverlap.faces.size()); - _ASSERT(meshSource.vecFaceArea.GetRows() == meshSource.faces.size()); - _ASSERT(meshTarget.vecFaceArea.GetRows() == meshTarget.faces.size()); + _ASSERT(meshSource.vecFaceArea.size() == meshSource.faces.size()); + _ASSERT(meshTarget.vecFaceArea.size() == meshTarget.faces.size()); for (int i = 0; i < meshOverlap.faces.size(); i++) { dSourceArea[ meshOverlap.vecSourceFaceIx[i] ] += meshOverlap.vecFaceArea[i]; diff --git a/src/GenerateOverlapMesh_v1.cpp b/src/GenerateOverlapMeshEdge.cpp similarity index 78% rename from src/GenerateOverlapMesh_v1.cpp rename to src/GenerateOverlapMeshEdge.cpp index afcfe87..993d9dc 100644 --- a/src/GenerateOverlapMesh_v1.cpp +++ b/src/GenerateOverlapMeshEdge.cpp @@ -1,6 +1,6 @@ /////////////////////////////////////////////////////////////////////////////// /// -/// \file GenerateOverlapMesh.cpp +/// \file GenerateOverlapMeshEdge.cpp /// \author Paul Ullrich /// \version March 7, 2014 /// @@ -22,17 +22,36 @@ #include "OverlapMesh.h" #include "netcdfcpp.h" +#include "NetCDFUtilities.h" #include /////////////////////////////////////////////////////////////////////////////// extern "C" -int GenerateOverlapMesh_v1(std::string strMeshA, std::string strMeshB, Mesh& meshOverlap, std::string strOverlapMesh, std::string strMethod, bool fNoValidate) { +int GenerateOverlapMeshEdge( + std::string strMeshA, + std::string strMeshB, + Mesh& meshOverlap, + std::string strOverlapMesh, + std::string strOutputFormat, + std::string strMethod, + bool fNoValidate +) { NcError error(NcError::silent_nonfatal); try { + // Check command line parameters (data type arguments) + STLStringHelper::ToLower(strOutputFormat); + + NcFile::FileFormat eOutputFormat = + GetNcFileFormatFromString(strOutputFormat); + if (eOutputFormat == NcFile::BadFormat) { + _EXCEPTION1("Invalid \"out_format\" value (%s), " + "expected [Classic|Offset64Bits|Netcdf4|Netcdf4Classic]", + strOutputFormat.c_str()); + } // Method string OverlapMeshMethod method; @@ -97,12 +116,12 @@ try { AnnounceEndBlock(NULL); AnnounceStartBlock("Construct overlap mesh"); - GenerateOverlapMesh_v1(meshA, meshB, meshOverlap, method); + GenerateOverlapMeshEdge(meshA, meshB, meshOverlap, method); AnnounceEndBlock(NULL); // Write the overlap mesh AnnounceStartBlock("Writing overlap mesh"); - meshOverlap.Write(strOverlapMesh.c_str()); + meshOverlap.Write(strOverlapMesh.c_str(), eOutputFormat); AnnounceEndBlock(NULL); } catch(Exception & e) { diff --git a/src/GenerateOverlapMeshExe.cpp b/src/GenerateOverlapMeshExe.cpp index aeab154..abad701 100644 --- a/src/GenerateOverlapMeshExe.cpp +++ b/src/GenerateOverlapMeshExe.cpp @@ -1,6 +1,6 @@ /////////////////////////////////////////////////////////////////////////////// /// -/// \file GenerateOverlapMesh.cpp +/// \file GenerateOverlapMeshExe.cpp /// \author Paul Ullrich /// \version March 7, 2014 /// @@ -17,6 +17,7 @@ #include "Announce.h" #include "CommandLine.h" #include "GridElements.h" +#include "STLStringHelper.h" #include "TempestRemapAPI.h" @@ -39,6 +40,9 @@ int main(int argc, char** argv) { // Overlap grid generation method std::string strMethod; + // Overlap grid generation algorithm + std::string strAlgorithm; + // No validation of the meshes bool fNoValidate; @@ -62,6 +66,7 @@ int main(int argc, char** argv) { CommandLineString(strOverlapMesh, "out", "overlap.g"); CommandLineString(strOutputFormat, "out_format", "netcdf4"); CommandLineStringD(strMethod, "method", "fuzzy", "(fuzzy|exact|mixed)"); + CommandLineStringD(strAlgorithm, "alg", "kdx", "(edge|kdx|lint)"); CommandLineBool(fNoValidate, "novalidate"); CommandLineBool(fHasConcaveFacesA, "concavea"); CommandLineBool(fHasConcaveFacesB, "concaveb"); @@ -73,18 +78,65 @@ int main(int argc, char** argv) { AnnounceBanner(); + // Change algorithm to lowercase + STLStringHelper::ToLower(strAlgorithm); + // Call the actual mesh generator - Mesh meshOverlap; - int err = - GenerateOverlapMesh( - strMeshA, strMeshB, - meshOverlap, strOverlapMesh, strOutputFormat, - strMethod, fNoValidate, - fHasConcaveFacesA, fHasConcaveFacesB, - fAllowNoOverlap, - fVerbose); - - if (err) exit(err); + Mesh meshOverlap; + + // Edge algorithm + if (strAlgorithm == "edge") { + if (fHasConcaveFacesA) { + Announce("WARNING: --alg \"edge\" does not support concave faces. Argument --concavea ignored."); + } + if (fHasConcaveFacesB) { + Announce("WARNING: --alg \"edge\" does not support concave faces. Argument --concaveb ignored."); + } + if (fAllowNoOverlap) { + Announce("WARNING: Argument --allow_no_overlap has no effect with --alg \"edge\""); + } + + int err = + GenerateOverlapMeshEdge( + strMeshA, strMeshB, + meshOverlap, strOverlapMesh, strOutputFormat, + strMethod, fNoValidate); + + if (err) exit(err); + + // KD-tree search + } else if (strAlgorithm == "kdx" ) { + int err = + GenerateOverlapMeshKdx( + strMeshA, strMeshB, + meshOverlap, strOverlapMesh, strOutputFormat, + strMethod, fNoValidate, + fHasConcaveFacesA, fHasConcaveFacesB, + fAllowNoOverlap, + fVerbose); + + if (err) exit(err); + + // Line search with interval tree + } else if (strAlgorithm == "lint" ) { + if (fAllowNoOverlap) { + Announce("WARNING: Argument --allow_no_overlap has no effect with --alg \"lint\""); + } + + int err = + GenerateOverlapMeshLint( + strMeshA, strMeshB, + meshOverlap, strOverlapMesh, strOutputFormat, + strMethod, fNoValidate, + fHasConcaveFacesA, fHasConcaveFacesB, + false, "/tmp", + fVerbose); + + if (err) exit(err); + + } else { + _EXCEPTIONT("Invalid value of --alg, expected \"edge\", \"kdx\" or \"lint\""); + } AnnounceBanner(); diff --git a/src/GenerateOverlapMeshExe_v1.cpp b/src/GenerateOverlapMeshExe_v1.cpp deleted file mode 100644 index 5e1fbf5..0000000 --- a/src/GenerateOverlapMeshExe_v1.cpp +++ /dev/null @@ -1,66 +0,0 @@ -/////////////////////////////////////////////////////////////////////////////// -/// -/// \file GenerateOverlapMesh.cpp -/// \author Paul Ullrich -/// \version March 7, 2014 -/// -/// -/// Copyright 2000-2014 Paul Ullrich -/// -/// This file is distributed as part of the Tempest source code package. -/// Permission is granted to use, copy, modify and distribute this -/// source code and its documentation under the terms of the GNU General -/// Public License. This software is provided "as is" without express -/// or implied warranty. -/// - -#include "Announce.h" -#include "CommandLine.h" -#include "Exception.h" -#include "GridElements.h" - -#include "TempestRemapAPI.h" - -/////////////////////////////////////////////////////////////////////////////// - -int main(int argc, char** argv) { - - // Input mesh A - std::string strMeshA; - - // Input mesh B - std::string strMeshB; - - // Output mesh file - std::string strOverlapMesh; - - // Overlap grid generation method - std::string strMethod; - - // No validation of the meshes - bool fNoValidate; - - // Parse the command line - BeginCommandLine() - CommandLineString(strMeshA, "a", ""); - CommandLineString(strMeshB, "b", ""); - CommandLineString(strOverlapMesh, "out", "overlap.g"); - CommandLineStringD(strMethod, "method", "fuzzy", "(fuzzy|exact|mixed)"); - CommandLineBool(fNoValidate, "novalidate"); - - ParseCommandLine(argc, argv); - EndCommandLine(argv) - - AnnounceBanner(); - - // Call the actual mesh generator - Mesh meshOverlap; - int err = GenerateOverlapMesh_v1(strMeshA, strMeshB, meshOverlap, strOverlapMesh, strMethod, fNoValidate); - - AnnounceBanner(); - - if (err) exit(err); - else return 0; -} - -/////////////////////////////////////////////////////////////////////////////// diff --git a/src/GenerateOverlapMesh.cpp b/src/GenerateOverlapMeshKdx.cpp similarity index 97% rename from src/GenerateOverlapMesh.cpp rename to src/GenerateOverlapMeshKdx.cpp index c01b14c..7347c58 100644 --- a/src/GenerateOverlapMesh.cpp +++ b/src/GenerateOverlapMeshKdx.cpp @@ -1,6 +1,6 @@ /////////////////////////////////////////////////////////////////////////////// /// -/// \file GenerateOverlapMesh.cpp +/// \file GenerateOverlapMeshKdx.cpp /// \author Paul Ullrich /// \version March 7, 2014 /// @@ -29,7 +29,7 @@ /////////////////////////////////////////////////////////////////////////////// extern "C" -int GenerateOverlapWithMeshes ( +int GenerateOverlapWithMeshesKdx ( Mesh & meshA, Mesh & meshB, Mesh & meshOverlap, @@ -81,7 +81,7 @@ int GenerateOverlapWithMeshes ( meshOverlap.type = Mesh::MeshType_Overlap; AnnounceStartBlock ( "Construct overlap mesh" ); - GenerateOverlapMesh_v2 ( + GenerateOverlapMeshKdx ( meshA, meshB, meshOverlap, method, @@ -145,7 +145,7 @@ int GenerateOverlapWithMeshes ( /////////////////////////////////////////////////////////////////////////////// extern "C" -int GenerateOverlapMesh( +int GenerateOverlapMeshKdx( std::string strMeshA, std::string strMeshB, Mesh & meshOverlap, @@ -226,7 +226,7 @@ int GenerateOverlapMesh( AnnounceEndBlock ( NULL ); int err = - GenerateOverlapWithMeshes ( + GenerateOverlapWithMeshesKdx ( meshA, meshB, meshOverlap, strOverlapMesh, diff --git a/src/GenerateOverlapMeshLint.cpp b/src/GenerateOverlapMeshLint.cpp new file mode 100644 index 0000000..d6f3eaa --- /dev/null +++ b/src/GenerateOverlapMeshLint.cpp @@ -0,0 +1,451 @@ +/////////////////////////////////////////////////////////////////////////////// +/// +/// \file GenerateOverlapMeshLint.cpp +/// \author Paul Ullrich +/// \version August 19, 2021 +/// +/// +/// Copyright 2021 Paul Ullrich +/// +/// This file is distributed as part of the Tempest source code package. +/// Permission is granted to use, copy, modify and distribute this +/// source code and its documentation under the terms of the GNU General +/// Public License. This software is provided "as is" without express +/// or implied warranty. +/// + +#include "Announce.h" +#include "CommandLine.h" +#include "STLStringHelper.h" +#include "Exception.h" +#include "GridElements.h" +#include "OverlapMesh.h" + +#include "netcdfcpp.h" +#include "NetCDFUtilities.h" + +#include + +#if defined(TEMPEST_MPIOMP) +#include +#endif + +/////////////////////////////////////////////////////////////////////////////// + +extern "C" +int GenerateOverlapWithMeshesLint ( + Mesh & meshA, + Mesh & meshB, + Mesh & meshOverlap, + std::string strOverlapMesh, + std::string strOutputFormat, + std::string strMethod, + const bool fNoValidate, + const bool fHasConcaveFacesA, + const bool fHasConcaveFacesB, + const bool fParallel, + std::string strTempDir, + const bool fVerbose +) { + + NcError error ( NcError::silent_nonfatal ); + + try + { + // Check command line parameters (data type arguments) + STLStringHelper::ToLower(strOutputFormat); + + NcFile::FileFormat eOutputFormat = + GetNcFileFormatFromString(strOutputFormat); + if (eOutputFormat == NcFile::BadFormat) { + _EXCEPTION1("Invalid \"out_format\" value (%s), " + "expected [Classic|Offset64Bits|Netcdf4|Netcdf4Classic]", + strOutputFormat.c_str()); + } + + // Method string + OverlapMeshMethod method; + STLStringHelper::ToLower ( strMethod ); + + if (strMethod == "fuzzy") { + method = OverlapMeshMethod_Fuzzy; + + } else if (strMethod == "exact") { + method = OverlapMeshMethod_Exact; + + } else if (strMethod == "mixed") { + method = OverlapMeshMethod_Mixed; + + } else { + _EXCEPTIONT("Invalid \"method\" value, expected \"fuzzy\", \"exact\" or \"mixed\""); + } + + // Convexify Mesh A + if (fHasConcaveFacesA) { + Mesh meshTemp = meshA; + ConvexifyMesh(meshTemp, meshA, fVerbose); + } + + // Validate Mesh A + if (!fNoValidate) { + AnnounceStartBlock("Validate mesh A"); + meshA.Validate(); + AnnounceEndBlock(NULL); + } + + // Convexify Mesh B + if (fHasConcaveFacesB) { + Mesh meshTemp = meshB; + ConvexifyMesh(meshTemp, meshB, fVerbose); + } + + // Validate Mesh B + if (!fNoValidate) { + AnnounceStartBlock("Validate mesh B"); + meshB.Validate(); + AnnounceEndBlock(NULL); + } + + // Generate meshes + meshOverlap.type = Mesh::MeshType_Overlap; + + AnnounceStartBlock ( "Construct overlap mesh" ); + GenerateOverlapMeshLint ( + meshA, meshB, + meshOverlap, + method, + fParallel, + strTempDir, + fVerbose ); + AnnounceEndBlock ( NULL ); + + /* + GenerateOverlapMeshFromFace( + meshA, + meshB, + 0, + meshOverlap, + method); + + // Construct the reverse node array on both meshes + AnnounceStartBlock("Constructing reverse node array on input mesh"); + meshA.ConstructReverseNodeArray(); + AnnounceEndBlock(NULL); + + AnnounceStartBlock("Constructing reverse node array on output mesh"); + meshB.ConstructReverseNodeArray(); + AnnounceEndBlock(NULL); + + // Equalize nearly coincident nodes on these Meshes + AnnounceStartBlock("Equalize coicident Nodes"); + EqualizeCoincidentNodes(meshA, meshB); + AnnounceEndBlock(NULL); + + // Construct the overlap mesh + Mesh meshOverlap; + + AnnounceStartBlock("Construct overlap mesh"); + GenerateOverlapMesh_v1(meshA, meshB, meshOverlap, method); + AnnounceEndBlock(NULL); + */ + + // Write the overlap mesh + if ( strOverlapMesh.size() ) + { + AnnounceStartBlock("Writing overlap mesh"); + meshOverlap.Write(strOverlapMesh.c_str(), eOutputFormat); + AnnounceEndBlock(NULL); + } + + return 0; + + } + catch ( Exception& e ) + { + Announce ( e.ToString().c_str() ); + return ( 0 ); + + } + catch ( ... ) + { + return ( 0 ); + } +} + +/////////////////////////////////////////////////////////////////////////////// + +extern "C" +int GenerateOverlapMeshLint( + std::string strMeshA, + std::string strMeshB, + Mesh & meshOverlap, + std::string strOverlapMesh, + std::string strOutputFormat, + std::string strMethod, + const bool fNoValidate, + const bool fHasConcaveFacesA, + const bool fHasConcaveFacesB, + const bool fParallel, + std::string strTempDir, + const bool fVerbose +) { + + NcError error ( NcError::silent_nonfatal ); + + try + { + // Set output format + STLStringHelper::ToLower(strOutputFormat); + + NcFile::FileFormat eOutputFormat = + GetNcFileFormatFromString(strOutputFormat); + if (eOutputFormat == NcFile::BadFormat) { + _EXCEPTION1("Invalid \"out_format\" value (%s), " + "expected [Classic|Offset64Bits|Netcdf4|Netcdf4Classic]", + strOutputFormat.c_str()); + } + + // Prepare for parallel computation + int nMPIRank = 0; + int nMPISize = 1; + +#if defined(TEMPEST_MPIOMP) + if (fParallel) { + MPI_Comm_rank(MPI_COMM_WORLD, &nMPIRank); + MPI_Comm_size(MPI_COMM_WORLD, &nMPISize); + } +#endif + + // Input meshes + Mesh meshA; + Mesh meshB; + + // Run in serial + if (nMPISize == 1) { + AnnounceStartBlock("Loading mesh A"); + meshA.Read(strMeshA); + meshA.RemoveZeroEdges(); + AnnounceEndBlock(NULL); + + AnnounceStartBlock("Loading mesh B"); + meshB.Read(strMeshB); + meshB.RemoveZeroEdges(); + AnnounceEndBlock(NULL); + + int err = + GenerateOverlapWithMeshesLint ( + meshA, meshB, + meshOverlap, + strOverlapMesh, + strOutputFormat, + strMethod, + fNoValidate, + fHasConcaveFacesA, + fHasConcaveFacesB, + fParallel, strTempDir, + fVerbose); + + return err; + + // Run in parallel + } else { + +#if !defined(TEMPEST_MPIOMP) + _EXCEPTIONT("MPI not enabled: Cannot run in parallel"); +#else + // Number of pieces (used in parallel operation) + int nPiecesA = sqrt(nMPISize); + int nPiecesB = nMPISize / nPiecesA; + int nPiecesTotal = nPiecesA * nPiecesB; + + // Split mesh in serial (on MPI Rank 0) + if (nMPIRank == 0) { + AnnounceStartBlock("Loading mesh A"); + meshA.Read(strMeshA); + meshA.RemoveZeroEdges(); + AnnounceEndBlock(NULL); + + Announce("Splitting mesh A into %i pieces", nPiecesA); + meshA.SplitAndWrite(strTempDir + std::string("/meshA"), nPiecesA); + + AnnounceStartBlock("Loading mesh B"); + meshB.Read(strMeshB); + meshB.RemoveZeroEdges(); + AnnounceEndBlock(NULL); + + AnnounceStartBlock("Splitting mesh B into %i pieces", nPiecesB); + meshB.SplitAndWrite(strTempDir + std::string("/meshB"), nPiecesB); + AnnounceEndBlock(NULL); + } + + // Synchronize tasks + if (fParallel) { + MPI_Barrier(MPI_COMM_WORLD); + } + + // Computer overlap meshes in parallel + AnnounceStartBlock("Computing overlap meshes (in parallel)"); + + int iPieceAB = 0; + for (int iPieceA = 0; iPieceA < nPiecesA; iPieceA++) { + for (int iPieceB = 0; iPieceB < nPiecesB; iPieceB++) { + + if (iPieceAB % nMPISize == nMPIRank) { + + char szIndexA[10]; + snprintf(szIndexA, 10, "%06i", iPieceA); + char szIndexB[10]; + snprintf(szIndexB, 10, "%06i", iPieceB); + + std::string strMeshAFile = strTempDir + std::string("/meshA") + szIndexA + std::string(".g"); + std::string strMeshBFile = strTempDir + std::string("/meshB") + szIndexB + std::string(".g"); + std::string strMeshOutFile = strTempDir + std::string("/meshOvA") + szIndexA + std::string("B") + szIndexB + std::string(".g"); + std::string strMeshLogFile = strTempDir + std::string("/logA") + szIndexA + std::string("B") + szIndexB + std::string(".txt"); + + FILE * fpLog = fopen(strMeshLogFile.c_str(), "w"); + if (fpLog == NULL) { + _EXCEPTION1("Unable to open log file \"%s\" for writing", strMeshLogFile.c_str()); + } + AnnounceSetOutputBuffer(fpLog); + AnnounceOutputOnAllRanks(); + + Announce("Processor %i computing overlap on mesh pieces A%i and B%i\n", nMPIRank, iPieceA, iPieceB); + + Mesh meshA(strMeshAFile); + Mesh meshB(strMeshBFile); + Mesh meshOverlapSub; + + int err = + GenerateOverlapWithMeshesLint ( + meshA, meshB, + meshOverlapSub, + strMeshOutFile, + strOutputFormat, + strMethod, + fNoValidate, + fHasConcaveFacesA, + fHasConcaveFacesB, + fParallel, strTempDir, + fVerbose); + + AnnounceOnlyOutputOnRankZero(); + AnnounceSetOutputBuffer(stdout); + + if (err) { + MPI_Abort(MPI_COMM_WORLD, err); + return err; + } + } + + iPieceAB++; + } + } + + AnnounceEndBlock("Done"); + + // Synchronize tasks + if (fParallel) { + MPI_Barrier(MPI_COMM_WORLD); + } + + + + // Merge overlap meshes in serial (on MPI Rank 0) + if (nMPIRank == 0) { + AnnounceStartBlock("Merging meshes"); + + Mesh meshCombined; + meshCombined.Read(strTempDir + std::string("/meshOvA000000B000000.g")); + meshCombined.BeginAppend(); + + int iPieceAB = 0; + for (int iPieceA = 0; iPieceA < nPiecesA; iPieceA++) { + for (int iPieceB = 0; iPieceB < nPiecesB; iPieceB++) { + if ((iPieceA == 0) && (iPieceB == 0)) { + continue; + } + + char szIndexA[10]; + snprintf(szIndexA, 10, "%06i", iPieceA); + char szIndexB[10]; + snprintf(szIndexB, 10, "%06i", iPieceB); + std::string strMeshOvFile = strTempDir + std::string("/meshOvA") + szIndexA + std::string("B") + szIndexB + std::string(".g"); + + Mesh meshOther; + meshOther.Read(strMeshOvFile); + meshCombined.Append(meshOther); + } + } + + meshCombined.EndAppend(); + + AnnounceEndBlock("Done"); + + AnnounceStartBlock("Writing final overlap mesh"); + meshCombined.Write(strOverlapMesh); + AnnounceEndBlock("Done"); + } + + // Clean-up files + if (nMPIRank == 0) { + AnnounceStartBlock("Cleanup temporary files"); + + for (int iPieceA = 0; iPieceA < nPiecesA; iPieceA++) { + char szIndexA[10]; + snprintf(szIndexA, 10, "%06i", iPieceA); + std::string strMeshAFile = strTempDir + std::string("/meshA") + szIndexA + std::string(".g"); + Announce("\"%s\"", strMeshAFile.c_str()); + remove(strMeshAFile.c_str()); + } + for (int iPieceB = 0; iPieceB < nPiecesB; iPieceB++) { + char szIndexB[10]; + snprintf(szIndexB, 10, "%06i", iPieceB); + std::string strMeshBFile = strTempDir + std::string("/meshB") + szIndexB + std::string(".g"); + Announce("\"%s\"", strMeshBFile.c_str()); + remove(strMeshBFile.c_str()); + } + + for (int iPieceA = 0; iPieceA < nPiecesA; iPieceA++) { + for (int iPieceB = 0; iPieceB < nPiecesB; iPieceB++) { + + char szIndexA[10]; + snprintf(szIndexA, 10, "%06i", iPieceA); + char szIndexB[10]; + snprintf(szIndexB, 10, "%06i", iPieceB); + std::string strMeshOvFile = strTempDir + std::string("/meshOvA") + szIndexA + std::string("B") + szIndexB + std::string(".g"); + std::string strMeshLogFile = strTempDir + std::string("/logA") + szIndexA + std::string("B") + szIndexB + std::string(".txt"); + + Announce("\"%s\"", strMeshOvFile.c_str()); + remove(strMeshOvFile.c_str()); + //Announce("\"%s\"", strMeshLogFile.c_str()); + //remove(strMeshLogFile.c_str()); + } + } + + AnnounceEndBlock("Done"); + } +#endif + } + + return 0; + } + catch ( Exception& e ) + { + Announce ( e.ToString().c_str() ); +#if defined(TEMPEST_MPIOMP) + int flag; + MPI_Initialized(&flag); + if (flag) { + MPI_Abort(MPI_COMM_WORLD, -1); + } +#endif + return ( 0 ); + + } + catch ( ... ) + { + return ( 0 ); + } +} + +/////////////////////////////////////////////////////////////////////////////// diff --git a/src/GenerateOverlapMeshParallelExe.cpp b/src/GenerateOverlapMeshParallelExe.cpp new file mode 100644 index 0000000..f140bee --- /dev/null +++ b/src/GenerateOverlapMeshParallelExe.cpp @@ -0,0 +1,137 @@ +/////////////////////////////////////////////////////////////////////////////// +/// +/// \file GenerateOverlapMeshParallelExe.cpp +/// \author Paul Ullrich +/// \version August 23, 2021 +/// +/// +/// Copyright 2021 Paul Ullrich +/// +/// This file is distributed as part of the Tempest source code package. +/// Permission is granted to use, copy, modify and distribute this +/// source code and its documentation under the terms of the GNU General +/// Public License. This software is provided "as is" without express +/// or implied warranty. +/// + +#include "Announce.h" +#include "CommandLine.h" +#include "GridElements.h" +#include "STLStringHelper.h" + +#include "TempestRemapAPI.h" + +#if defined(TEMPEST_MPIOMP) +#include +#endif + +/////////////////////////////////////////////////////////////////////////////// + +int main(int argc, char** argv) { + +#if defined(TEMPEST_MPIOMP) + MPI_Init(&argc, &argv); + AnnounceOnlyOutputOnRankZero(); +#endif + + // Input mesh A + std::string strMeshA; + + // Input mesh B + std::string strMeshB; + + // Output mesh file + std::string strOverlapMesh; + + // Output format + std::string strOutputFormat; + + // Overlap grid generation method + std::string strMethod; + + // Overlap grid generation algorithm + std::string strAlgorithm; + + // No validation of the meshes + bool fNoValidate; + + // Concave elements may be present in mesh A + bool fHasConcaveFacesA; + + // Concave elements may be present in mesh B + bool fHasConcaveFacesB; + + // Allow for the case of no overlap element found (enabling this may + // cause the mesh generator to generate an incomplete overlap mesh) + bool fAllowNoOverlap; + + // Run in parallel + bool fParallel; + + // Temporary file directory + std::string strTempDir; + + // Verbose + bool fVerbose; + + // Parse the command line + BeginCommandLine() + CommandLineString(strMeshA, "a", ""); + CommandLineString(strMeshB, "b", ""); + CommandLineString(strOverlapMesh, "out", "overlap.g"); + CommandLineString(strOutputFormat, "out_format", "netcdf4"); + CommandLineStringD(strMethod, "method", "fuzzy", "(fuzzy|exact|mixed)"); + CommandLineStringD(strAlgorithm, "alg", "lint", "(lint)"); + CommandLineBool(fNoValidate, "novalidate"); + CommandLineBool(fHasConcaveFacesA, "concavea"); + CommandLineBool(fHasConcaveFacesB, "concaveb"); + CommandLineBool(fAllowNoOverlap, "allow_no_overlap"); + CommandLineString(strTempDir, "tmpdir", "/tmp"); + CommandLineBool(fVerbose, "verbose"); + + ParseCommandLine(argc, argv); + EndCommandLine(argv) + + AnnounceBanner(); + +#if !defined(TEMPEST_MPIOMP) + _EXCEPTIONT("TempestRemap was not compiled with MPI: Parallel operation not supported"); +#endif + + // Change algorithm to lowercase + STLStringHelper::ToLower(strAlgorithm); + + // Call the actual mesh generator + Mesh meshOverlap; + + // Edge algorithm + if (strAlgorithm == "lint" ) { + if (fAllowNoOverlap) { + Announce("WARNING: Argument --allow_no_overlap has no effect with --alg \"lint\""); + } + + int err = + GenerateOverlapMeshLint( + strMeshA, strMeshB, + meshOverlap, strOverlapMesh, strOutputFormat, + strMethod, fNoValidate, + fHasConcaveFacesA, fHasConcaveFacesB, + true, strTempDir, + fVerbose); + + if (err) exit(err); + + } else { + _EXCEPTIONT("Invalid value of --alg, expected \"lint\""); + } + + AnnounceBanner(); + +#if defined(TEMPEST_MPIOMP) + MPI_Finalize(); +#endif + + return 0; +} + +/////////////////////////////////////////////////////////////////////////////// diff --git a/src/GenerateRectilinearMeshFromFile.cpp b/src/GenerateRectilinearMeshFromFile.cpp index 7e0f254..a001373 100644 --- a/src/GenerateRectilinearMeshFromFile.cpp +++ b/src/GenerateRectilinearMeshFromFile.cpp @@ -232,6 +232,8 @@ try { // Build the mesh AnnounceStartBlock("Building mesh"); + long lTotalFaces = 0; + for (long j = 0; j < lDim0Size; j++) { for (long i = 0; i < lDim1Size; i++) { @@ -239,6 +241,8 @@ try { continue; } + lTotalFaces++; + int iGivenValues = ((ixV(j ,i ) != InvalidNode)?(1):(0)) + ((ixV(j+1,i ) != InvalidNode)?(1):(0)) @@ -313,6 +317,17 @@ try { AnnounceEndBlock("Done"); + // Store rectilinear coordinates if all points are available + if (lTotalFaces == lDim0Size * lDim1Size) { + mesh.vecGridDimSize.resize(2); + mesh.vecGridDimSize[0] = lDim0Size; + mesh.vecGridDimSize[1] = lDim1Size; + + mesh.vecGridDimName.resize(2); + mesh.vecGridDimName[0] = "dim0"; + mesh.vecGridDimName[1] = "dim1"; + } + AnnounceEndBlock("Done"); } else { diff --git a/src/GridElements.cpp b/src/GridElements.cpp index 0066a1d..8f5c1c3 100644 --- a/src/GridElements.cpp +++ b/src/GridElements.cpp @@ -28,6 +28,7 @@ #include #include #include +#include #include "netcdfcpp.h" #include "triangle.h" @@ -125,6 +126,37 @@ void Mesh::ConstructReverseNodeArray() { /////////////////////////////////////////////////////////////////////////////// +Real Mesh::SumFaceAreas() const { + + // Calculate accumulated area carefully + static const int Jump = 10; + std::vector vecFaceAreaBak; + vecFaceAreaBak.resize(vecFaceArea.size()); + memcpy(&(vecFaceAreaBak[0]), &(vecFaceArea[0]), + vecFaceArea.size() * sizeof(double)); + + for (;;) { + if (vecFaceAreaBak.size() == 1) { + break; + } + for (int i = 0; i <= (vecFaceAreaBak.size()-1) / Jump; i++) { + int ixRef = Jump * i; + vecFaceAreaBak[i] = vecFaceAreaBak[ixRef]; + for (int j = 1; j < Jump; j++) { + if (ixRef + j >= vecFaceAreaBak.size()) { + break; + } + vecFaceAreaBak[i] += vecFaceAreaBak[ixRef + j]; + } + } + vecFaceAreaBak.resize((vecFaceAreaBak.size()-1) / Jump + 1); + } + + return vecFaceAreaBak[0]; +} + +/////////////////////////////////////////////////////////////////////////////// + Real Mesh::CalculateFaceAreas( bool fContainsConcaveFaces ) { @@ -135,7 +167,7 @@ Real Mesh::CalculateFaceAreas( // Calculate areas int nCount = 0; - vecFaceArea.Allocate(faces.size()); + vecFaceArea.resize(faces.size(), 0); // Calculate the area of each Face if (fContainsConcaveFaces) { @@ -165,31 +197,8 @@ Real Mesh::CalculateFaceAreas( } } - // Calculate accumulated area carefully - static const int Jump = 10; - std::vector vecFaceAreaBak; - vecFaceAreaBak.resize(vecFaceArea.GetRows()); - memcpy(&(vecFaceAreaBak[0]), &(vecFaceArea[0]), - vecFaceArea.GetRows() * sizeof(double)); - - for (;;) { - if (vecFaceAreaBak.size() == 1) { - break; - } - for (int i = 0; i <= (vecFaceAreaBak.size()-1) / Jump; i++) { - int ixRef = Jump * i; - vecFaceAreaBak[i] = vecFaceAreaBak[ixRef]; - for (int j = 1; j < Jump; j++) { - if (ixRef + j >= vecFaceAreaBak.size()) { - break; - } - vecFaceAreaBak[i] += vecFaceAreaBak[ixRef + j]; - } - } - vecFaceAreaBak.resize((vecFaceAreaBak.size()-1) / Jump + 1); - } - - return vecFaceAreaBak[0]; + // Return the sum of face areas + return SumFaceAreas(); } /////////////////////////////////////////////////////////////////////////////// @@ -197,12 +206,12 @@ Real Mesh::CalculateFaceAreas( Real Mesh::CalculateFaceAreasFromOverlap( const Mesh & meshOverlap ) { - if (meshOverlap.vecFaceArea.GetRows() == 0) { + if (meshOverlap.vecFaceArea.size() == 0) { _EXCEPTIONT("MeshOverlap Face Areas have not been calculated"); } // Set all Face areas to zero - vecFaceArea.Allocate(faces.size()); + vecFaceArea.resize(faces.size(), 0); // Loop over all Faces in meshOverlap double dTotalArea = 0.0; @@ -210,7 +219,7 @@ Real Mesh::CalculateFaceAreasFromOverlap( for (int i = 0; i < meshOverlap.faces.size(); i++) { int ixFirstFace = meshOverlap.vecSourceFaceIx[i]; - if (ixFirstFace >= vecFaceArea.GetRows()) { + if (ixFirstFace >= vecFaceArea.size()) { _EXCEPTIONT("Overlap Mesh FirstFaceIx contains invalid " "Face index"); } @@ -223,6 +232,293 @@ Real Mesh::CalculateFaceAreasFromOverlap( /////////////////////////////////////////////////////////////////////////////// +void Mesh::GenerateLatLonBoxes( + std::vector< LatLonBox > & vecLatLonBoxes +) const { + vecLatLonBoxes.resize(faces.size()); + + for (int i = 0; i < faces.size(); i++) { + const Face & face = faces[i]; + + // Check if face contains pole points (if so we should get iWindingNumber != 0) + int iWindingNumber = 0; + for (int j = 0; j < face.edges.size(); j++) { + const Edge & edge = face.edges[j]; + + const Node & n1 = nodes[edge[0]]; + const Node & n2 = nodes[edge[1]]; + + // One of the edges crosses z=0; assume this does not contain the pole point + if ((n1.z <= -ReferenceTolerance) && (n2.z >= ReferenceTolerance)) { + iWindingNumber = 0; + break; + } + if ((n1.z >= ReferenceTolerance) && (n2.z <= -ReferenceTolerance)) { + iWindingNumber = 0; + break; + } + + // One of the points on this edge is a pole point (counts as face including pole points) + if (fabs(fabs(n1.z) - 1.0) < ReferenceTolerance) { + iWindingNumber = 2000; + break; + } + + // n1 is on the line of constant x (winding should be accounted for in previous step) + if (fabs(n1.y) < ReferenceTolerance) { + + // edge passes through the pole + // TODO: This only works for great circle arcs. Fix for lines of constant latitude. + if (fabs(n2.y) < ReferenceTolerance) { + if ((n1.x >= ReferenceTolerance) && (n1.x <= ReferenceTolerance)) { + iWindingNumber = 3000; + break; + } + if ((n1.x <= ReferenceTolerance) && (n1.x >= ReferenceTolerance)) { + iWindingNumber = 4000; + break; + } + } + continue; + } + + // Both endpoints of line lay on same side of y = 0 + if ((n1.y <= ReferenceTolerance) && (n2.y <= ReferenceTolerance)) { + continue; + } + if ((n1.y >= ReferenceTolerance) && (n2.y >= ReferenceTolerance)) { + continue; + } + + // Determine intersection with line (x,y,z)=(x,0,1) + double dDenom = n1.y * n2.z - n2.y * n1.z; + if (fabs(dDenom) < ReferenceTolerance) { + continue; + } + + double dXint = (n1.y * n2.x - n1.x * n2.y) / dDenom; + if (dXint > 0.0) { + if ((n1.y < ReferenceTolerance) && (n2.y > ReferenceTolerance)) { + iWindingNumber++; + } else { + iWindingNumber--; + } + } else { + if ((n1.y < ReferenceTolerance) && (n2.y > ReferenceTolerance)) { + iWindingNumber--; + } else { + iWindingNumber++; + } + } + } + + // Face containing a pole point; latlon box determined by minimum/maximum latitude + if (abs(iWindingNumber) > 1) { + + // TODO: Fix for faces that only contain the pole at one edge endpoint + double dLatExtentRad; + + for (int j = 0; j < face.edges.size(); j++) { + const Edge & edge = face.edges[j]; + + const Node & n1 = nodes[edge[0]]; + const Node & n2 = nodes[edge[1]]; + + double dLatRad; + double dLonRad; + + // Set the latitude extent + if (j == 0) { + if (n1.z < 0.0) { + dLatExtentRad = -0.5 * M_PI; + } else { + dLatExtentRad = 0.5 * M_PI; + } + } + + // Insert edge endpoint into box + nodes[edge[0]].ToLatLonRad(dLatRad, dLonRad); + if (fabs(dLatRad) < fabs(dLatExtentRad)) { + dLatExtentRad = dLatRad; + } + + // Edges that are great circle arcs + if (edge.type == Edge::Type_GreatCircleArc) { + + // Determine if latitude is maximized between endpoints + double dDotN1N2 = DotProduct(n1,n2); + + double dDenom = (n1.z + n2.z) * (dDotN1N2 - 1.0); + + // Either repeated point or Amax out of range + if (fabs(dDenom) < ReferenceTolerance) { + continue; + } + + // Maximum latitude occurs between endpoints of edge + double dAmax = (n1.z * dDotN1N2 - n2.z) / dDenom; + if ((dAmax > 0.0) && (dAmax < 1.0)) { + + Node n3 = n1 * (1.0 - dAmax) + n2 * dAmax; + n3.NormalizeInPlace(); + + dLatRad = n3.z; + if (dLatRad > 1.0) { + dLatRad = 0.5 * M_PI; + } else if (dLatRad < -1.0) { + dLatRad = -0.5 * M_PI; + } else { + dLatRad = asin(dLatRad); + } + if (fabs(dLatRad) < fabs(dLatExtentRad)) { + dLatExtentRad = dLatRad; + } + } + + // Edges that are lines of constant latitude + } else if (edge.type == Edge::Type_ConstantLatitude) { + + } else { + _EXCEPTION1("Unsupported edge type (%i)", edge.type); + } + } + + // Define latlon box + if (dLatExtentRad < 0.0) { + vecLatLonBoxes[i].is_null = false; + vecLatLonBoxes[i].lon[0] = 0.0; + vecLatLonBoxes[i].lon[1] = 2.0 * M_PI; + vecLatLonBoxes[i].lat[0] = -0.5 * M_PI; + vecLatLonBoxes[i].lat[1] = dLatExtentRad; + + } else { + vecLatLonBoxes[i].is_null = false; + vecLatLonBoxes[i].lon[0] = 0.0; + vecLatLonBoxes[i].lon[1] = 2.0 * M_PI; + vecLatLonBoxes[i].lat[0] = dLatExtentRad; + vecLatLonBoxes[i].lat[1] = 0.5 * M_PI; + } +/* + printf("%1.5e %1.5e %1.5e %1.5e\n", + vecLatLonBoxes[i].lat[0], + vecLatLonBoxes[i].lat[1], + vecLatLonBoxes[i].lon[0], + vecLatLonBoxes[i].lon[1]); +*/ + // Normal face + } else { + for (int j = 0; j < face.edges.size(); j++) { + const Edge & edge = face.edges[j]; + + // Edges that are great circle arcs + if (edge.type == Edge::Type_GreatCircleArc) { + + const Node & n1 = nodes[edge[0]]; + const Node & n2 = nodes[edge[1]]; + + double dDotN1N2 = DotProduct(n1,n2); + + double dDenom = (n1.z + n2.z) * (dDotN1N2 - 1.0); + + double dLatRad; + double dLonRad; + + // Insert first end points of the edge + nodes[edge[0]].ToLatLonRad(dLatRad, dLonRad); + vecLatLonBoxes[i].insert(dLatRad, dLonRad); + + // Either repeated point or Amax out of range + if (fabs(dDenom) < ReferenceTolerance) { + continue; + } + + // These points don't need to be inserted + //nodes[edge[1]].ToLatLonRad(dLatRad, dLonRad); + //vecLatLonBoxes[i].insert(dLatRad, dLonRad); + + // Maximum latitude occurs between endpoints of edge + double dAmax = (n1.z * dDotN1N2 - n2.z) / dDenom; + if ((dAmax > 0.0) && (dAmax < 1.0)) { + + Node n3 = n1 * (1.0 - dAmax) + n2 * dAmax; + n3.NormalizeInPlace(); + + dLatRad = n3.z; + if (dLatRad > 1.0) { + dLatRad = 0.5 * M_PI; + } else if (dLatRad < -1.0) { + dLatRad = -0.5 * M_PI; + } else { + dLatRad = asin(dLatRad); + } + vecLatLonBoxes[i].insert(dLatRad, dLonRad); + } + + // Edges that are lines of constant latitude + } else if (edge.type == Edge::Type_ConstantLatitude) { + const Node & n1 = nodes[edge[0]]; + const Node & n2 = nodes[edge[0]]; + + double dLatRad; + double dLonRad; + + n1.ToLatLonRad(dLatRad, dLonRad); + vecLatLonBoxes[i].insert(dLatRad, dLonRad); + + n2.ToLatLonRad(dLatRad, dLonRad); + vecLatLonBoxes[i].insert(dLatRad, dLonRad); + + } else { + _EXCEPTION1("Unsupported edge type (%i)", edge.type); + } + } + } + + //vecLatLonBoxes[i].lon[0] = 0.0; + //vecLatLonBoxes[i].lon[1] = 2.0 * M_PI; + //vecLatLonBoxes[i].lat[0] = -0.5 * M_PI; + //vecLatLonBoxes[i].lat[1] = 0.5 * M_PI; + + // Verify face is of non-zero size + _ASSERT(vecLatLonBoxes[i].lon[0] != vecLatLonBoxes[i].lon[1]); + _ASSERT(vecLatLonBoxes[i].lat[0] != vecLatLonBoxes[i].lat[1]); +/* + // TODO: REMOVE AFTER TESTING + // DEBUG: Check edge points are all in latlon box + for (int j = 0; j < face.edges.size(); j++) { + const Edge & edge = face.edges[j]; + + const Node & n1 = nodes[edge[0]]; + const Node & n2 = nodes[edge[1]]; + + for (int k = 0; k < 11; k++) { + double dA = static_cast(k)/10.0; + + Node n3 = n1 * (1.0 - dA) + n2 * dA; + + n3 = n3.Normalized(); + + double dLatRad; + double dLonRad; + n3.ToLatLonRad(dLatRad, dLonRad); + + if (!vecLatLonBoxes[i].contains_eps(dLatRad, dLonRad, 1.0e-14)) { + printf("\nA: %1.15e lat: %1.15e lon: %1.15e\n", dA, dLatRad, dLonRad); + printf("lat [%1.15e, %1.15e] \nlon [%1.15e, %1.15e]\n", + vecLatLonBoxes[i].lat[0], + vecLatLonBoxes[i].lat[1], + vecLatLonBoxes[i].lon[0], + vecLatLonBoxes[i].lon[1]); + _EXCEPTION(); + } + } + } +*/ + } +} + +/////////////////////////////////////////////////////////////////////////////// + void Mesh::ExchangeFirstAndSecondMesh() { // Verify all vectors are the same size @@ -347,7 +643,7 @@ void Mesh::Write( vecBlockSizes.resize(mapBlockSizes.size()); vecBlockSizeFaces.resize(mapBlockSizes.size()); - AnnounceStartBlock("Nodes per element"); + AnnounceStartBlock("Nodes per element:"); iterBlockSize = mapBlockSizes.begin(); iBlock = 1; for (; iterBlockSize != mapBlockSizes.end(); iterBlockSize++) { @@ -1520,9 +1816,164 @@ void Mesh::Read(const std::string & strFile) { /////////////////////////////////////////////////////////////////////////////// -void Mesh::RemoveZeroEdges() { +void Mesh::SplitAndWrite( + const std::string & strFilePrefix, + size_t sPieces +) { + _ASSERT(sPieces > 0); + + ConstructReverseNodeArray(); + + size_t nSplitMeshSize; + if (faces.size() % sPieces == 0) { + nSplitMeshSize = faces.size() / sPieces; + } else { + nSplitMeshSize = faces.size() / sPieces + 1; + } + _ASSERT(nSplitMeshSize * sPieces >= faces.size()); + + int iRootIndex = 0; + std::vector vecFacesUsed(faces.size(), false); + + for (size_t iCurrentPiece = 0; iCurrentPiece < sPieces; iCurrentPiece++) { + Mesh meshPiece; + std::map mapPieceNodes; - // Remove zero edges from all Faces + // Construct this piece of the mesh + while (meshPiece.faces.size() < nSplitMeshSize) { + if (iRootIndex == faces.size()) { + break; + } + if (vecFacesUsed[iRootIndex]) { + iRootIndex++; + continue; + } + + std::queue queueSearchFaces; + queueSearchFaces.push(iRootIndex); + + while ((!queueSearchFaces.empty()) && (meshPiece.faces.size() < nSplitMeshSize)) { + int iNextFace = queueSearchFaces.front(); + queueSearchFaces.pop(); + if (vecFacesUsed[iNextFace]) { + continue; + } + vecFacesUsed[iNextFace] = true; + + // Add this face to meshPiece, avoiding repeated nodes + const Face & face = faces[iNextFace]; + + Face faceCurrent(face.edges.size()); + for (int iFaceNode = 0; iFaceNode < face.edges.size(); iFaceNode++) { + const Node & nodeCurrent = nodes[face[iFaceNode]]; + auto it = mapPieceNodes.find(nodeCurrent); + if (it == mapPieceNodes.end()) { + int ixMeshNode = meshPiece.nodes.size(); + mapPieceNodes.insert( + std::pair( + nodeCurrent, + ixMeshNode)); + meshPiece.nodes.push_back(nodeCurrent); + faceCurrent.SetNode(iFaceNode, ixMeshNode); + } else { + faceCurrent.SetNode(iFaceNode, it->second); + } + } + meshPiece.faces.push_back(faceCurrent); + + // Add all unused neighbors of this face's vertices to the search queue + for (int iFaceNode = 0; iFaceNode < face.edges.size(); iFaceNode++) { + const std::set & setNodes = revnodearray[face[iFaceNode]]; + for (auto it : setNodes) { + if (!vecFacesUsed[it]) { + queueSearchFaces.push(it); + } + } + } + } + } + + // Write the mesh to disk + char szIndex[10]; + snprintf(szIndex, 10, "%06lu", iCurrentPiece); + std::string strMeshFile = strFilePrefix + szIndex + std::string(".g"); + AnnounceStartBlock("Writing \"%s\"", strMeshFile.c_str()); + meshPiece.Write(strMeshFile, NcFile::Netcdf4); + AnnounceEndBlock(NULL); + } +} + +/////////////////////////////////////////////////////////////////////////////// + +void Mesh::BeginAppend() { + nodemap.clear(); + for (int i = 0; i < nodes.size(); i++) { + nodemap.insert(NodeMapPair(nodes[i], i)); + } +} + +/////////////////////////////////////////////////////////////////////////////// + +void Mesh::Append( + const Mesh & meshOther +) { + vecFaceArea.clear(); + vecMask.clear(); + edgemap.clear(); + revnodearray.clear(); + vecMultiFaceMap.clear(); + + if (nodemap.size() != nodes.size()) { + BeginAppend(); + } + + if ((vecSourceFaceIx.size() != 0) && (meshOther.vecSourceFaceIx.size() != 0)) { + vecSourceFaceIx.reserve(vecSourceFaceIx.size() + meshOther.vecSourceFaceIx.size()); + } + if ((vecTargetFaceIx.size() != 0) && (meshOther.vecTargetFaceIx.size() != 0)) { + vecTargetFaceIx.reserve(vecTargetFaceIx.size() + meshOther.vecTargetFaceIx.size()); + } + + for (int i = 0; i < meshOther.faces.size(); i++) { + const Face & face = meshOther.faces[i]; + Face faceNew(face.size()); + + for (int j = 0; j < face.size(); j++) { + const Node & node = meshOther.nodes[face[j]]; + int ixNode; + auto it = nodemap.find(node); + if (it == nodemap.end()) { + ixNode = nodes.size(); + nodes.push_back(node); + nodemap.insert(NodeMapPair(node, ixNode)); + + } else { + ixNode = it->second; + } + faceNew.SetNode(j, ixNode); + } + + int ixNewFace = face.size(); + faces.push_back(faceNew); + + if ((vecSourceFaceIx.size() != 0) && (meshOther.vecSourceFaceIx.size() != 0)) { + vecSourceFaceIx.push_back(meshOther.vecSourceFaceIx[i]); + } + if ((vecTargetFaceIx.size() != 0) && (meshOther.vecTargetFaceIx.size() != 0)) { + vecTargetFaceIx.push_back(meshOther.vecTargetFaceIx[i]); + } + } +} + +/////////////////////////////////////////////////////////////////////////////// + +void Mesh::EndAppend() { + nodemap.clear(); +} + +/////////////////////////////////////////////////////////////////////////////// + +void Mesh::RemoveZeroEdges() { for (int i = 0; i < faces.size(); i++) { faces[i].RemoveZeroEdges(); } diff --git a/src/GridElements.h b/src/GridElements.h index e10fdfa..2b4fa43 100644 --- a/src/GridElements.h +++ b/src/GridElements.h @@ -21,6 +21,7 @@ #include "Defines.h" #include "CoordTransforms.h" +#include "LatLonBox.h" #include #include @@ -204,6 +205,16 @@ class Node { return (*this)/Magnitude();; } + /// + /// Project node onto the unit sphere + /// + void NormalizeInPlace() { + Real mag = Magnitude(); + x /= mag; + y /= mag; + z /= mag; + } + /// /// Magnitude of this node. /// @@ -211,6 +222,18 @@ class Node { return sqrt(x * x + y * y + z * z); } + /// + /// Lat/lon coordinates of this node in radians. + /// Return value of lat is in the interval [-0.5 * M_PI, 0.5 * M_PI] + /// Return value of lon is in the interval [0.0, 2.0 * M_PI) + /// + void ToLatLonRad( + double & lat, + double & lon + ) const { + XYZtoRLL_Rad(x, y, z, lon, lat); + } + /// /// Output node to stdout. /// @@ -228,7 +251,7 @@ typedef std::vector NodeVector; /// /// A map between Nodes and indices. /// -#if defined(OVERLAPMESH_RETAIN_REPEATED_NODES) +#if defined(OVERLAPMESH_USE_SORTED_MAP) typedef std::map NodeMap; #endif @@ -576,6 +599,13 @@ class Face { return edges[ix][0]; } + /// + /// Face size. + /// + size_t size() const { + return edges.size(); + } + /// /// Set a node. /// @@ -677,13 +707,18 @@ class Mesh { /// /// Vector of Face areas. /// - DataArray1D vecFaceArea; + std::vector vecFaceArea; /// /// Vector storing mask variable for this mesh. /// std::vector vecMask; + /// + /// A NodeMap used for appending Meshes. + /// + NodeMap nodemap; + /// /// EdgeMap for this mesh. /// @@ -742,6 +777,11 @@ class Mesh { /// void ConstructReverseNodeArray(); + /// + /// Sum the Face areas. + /// + Real SumFaceAreas() const; + /// /// Calculate Face areas. /// @@ -756,6 +796,13 @@ class Mesh { const Mesh & meshOverlap ); + /// + /// Generate the lat-lon rectangles for all faces in the mesh. + /// + void GenerateLatLonBoxes( + std::vector< LatLonBox > & vecLatLonBoxes + ) const; + /// /// Sort Faces by the opposite source mesh. /// @@ -789,6 +836,34 @@ class Mesh { /// void Read(const std::string & strFile); + /// + /// Split the mesh into the given number of pieces and write it to + /// multiple files. + /// + void SplitAndWrite( + const std::string & strFilePrefix, + size_t sPieces + ); + + /// + /// Indicate that we are beginning to append other Meshes to this Mesh. + /// It is not mandatory to call this prior to calling Append(), but + /// it does improve performance. + /// + void BeginAppend(); + + /// + /// Append meshOther to this Mesh. + /// + void Append( + const Mesh & meshOther + ); + + /// + /// Indicate that we are done appending other Meshes to this Mesh. + /// + void EndAppend(); + /// /// Remove zero edges from all Faces. /// diff --git a/src/LatLonBox.h b/src/LatLonBox.h new file mode 100644 index 0000000..f206cee --- /dev/null +++ b/src/LatLonBox.h @@ -0,0 +1,344 @@ +/////////////////////////////////////////////////////////////////////////////// +/// +/// \file LatLonBox.h +/// \author Paul Ullrich +/// \version August 19, 2021 +/// +/// +/// Copyright 2020 Paul Ullrich +/// +/// This file is distributed as part of the Tempest source code package. +/// Permission is granted to use, copy, modify and distribute this +/// source code and its documentation under the terms of the GNU General +/// Public License. This software is provided "as is" without express +/// or implied warranty. +/// + +#ifndef _LATLONBOX_H_ +#define _LATLONBOX_H_ + +#include "Exception.h" +#include +#include +#include +#include + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// A structure for storing a bounding box in latitude / longitude space. +/// The box is treated as singularly periodic in the longitudinal direction +/// if lon_periodic is set to true. +/// +template +class LatLonBox { + +public: + /// + /// Flag indicating this is a null box. + /// + bool is_null; + + /// + /// Flag indicating this is a regional. + /// + bool lon_periodic; + + /// + /// Width of the longitude variable. + /// + Type lon_width; + + /// + /// Bounding longitudes (endpoints are included). + /// + Type lon[2]; + + /// + /// Bounding latitudes (endpoints are included). + /// + Type lat[2]; + +public: + /// + /// Constructor. + /// + LatLonBox( + bool a_lon_periodic = true, + Type a_lon_width = static_cast(2.0 * M_PI) + ) : + is_null(true), + lon_periodic(a_lon_periodic), + lon_width(a_lon_width) + { + lon[0] = static_cast(0); + lon[1] = static_cast(0); + lat[0] = static_cast(0); + lat[1] = static_cast(0); + } + + /// + /// Width of this LatLonBox. + /// + Type width() const { + if (is_null) { + return 0; + } + + if (!lon_periodic) { + return (lon[1] - lon[0]); + } + + if (lon[0] == lon[1]) { + return static_cast(0); + } else if (lon[0] <= lon[1]) { + return (lon[1] - lon[0]); + } else { + return (lon[1] - lon[0] + lon_width); + } + } + + /// + /// Insert a point into this LatLonBox. + /// + void insert( + const Type & lat_pt, + const Type & lon_pt + ) { + if (is_null) { + is_null = false; + lat[0] = lat_pt; + lat[1] = lat_pt; + lon[0] = lon_pt; + lon[1] = lon_pt; + return; + } + if (lon_pt < static_cast(0)) { + std::stringstream strError; + strError << "lon_pt out of range (" << lon_pt << " < 0)" << std::endl; + _EXCEPTIONT(strError.str().c_str()); + } + if (lon_pt > lon_width) { + std::stringstream strError; + strError << "lon_pt out of range (" << lon_pt << " > "; + strError << lon_width << ")" << std::endl; + _EXCEPTIONT(strError.str().c_str()); + } + + // Expand latitudes + if (lat_pt > lat[1]) { + lat[1] = lat_pt; + } + if (lat_pt < lat[0]) { + lat[0] = lat_pt; + } + + // Expand longitude, if non-periodic + if (!lon_periodic) { + if (lon_pt > lon[1]) { + lon[1] = lon_pt; + } + if (lon_pt < lon[0]) { + lon[0] = lon_pt; + } + return; + } + + // New longitude lies within existing range + if (lon[0] <= lon[1]) { + if ((lon_pt >= lon[0]) && (lon_pt <= lon[1])) { + return; + } + } else { + if ((lon_pt >= lon[0]) || (lon_pt <= lon[1])) { + return; + } + } + + // New longitude lies outside of existing range + LatLonBox boxA(*this); + boxA.lon[0] = lon_pt; + + LatLonBox boxB(*this); + boxB.lon[1] = lon_pt; + + // The updated box is the box of minimum width + Type dWidthNow = width(); + Type dWidthA = boxA.width(); + Type dWidthB = boxB.width(); + + if ((dWidthA - dWidthNow < -1.0e-14) || (dWidthB - dWidthNow < -1.0e-14)) { + _EXCEPTIONT("Logic error"); + } + + if (dWidthA < dWidthB) { + (*this) = boxA; + } else { + (*this) = boxB; + } + } + + /// + /// Determine if this LatLonBox contains the given point. + /// + bool contains( + const Type & lat_pt, + const Type & lon_pt + ) { + // Check latitudes + if (lat[0] > lat_pt) { + return false; + } + if (lat[1] < lat_pt) { + return false; + } + + // Check longitudes, if non-periodic + if (!lon_periodic) { + if (lon[0] > lon_pt) { + return false; + } + if (lon[1] < lon_pt) { + return false; + } + return true; + } + + // This box crosses lon 360 + if (lon[0] > lon[1]) { + if (lon_pt >= lon[0]) { + return true; + } + if (lon_pt <= lon[1]) { + return true; + } + return false; + } + + // This box does not cross lon 360 + if ((lon_pt >= lon[0]) && (lon_pt <= lon[1])) { + return true; + } + return false; + } + + /// + /// Determine if this LatLonBox contains the given point, up to + /// a prescribed machine epsilon. + /// + bool contains_eps( + const Type & lat_pt, + const Type & lon_pt, + const Type & eps + ) { + // Check latitudes + if (lat[0] > lat_pt + eps) { + return false; + } + if (lat[1] < lat_pt - eps) { + return false; + } + + // Check longitudes, if non-periodic + if (!lon_periodic) { + if (lon[0] > lon_pt + eps) { + return false; + } + if (lon[1] < lon_pt - eps) { + return false; + } + return true; + } + + // This box crosses lon 360 + if (lon[0] > lon[1]) { + if (lon_pt + eps >= lon[0]) { + return true; + } + if (lon_pt - eps <= lon[1]) { + return true; + } + return false; + } + + // This box does not cross lon 360 + if ((lon_pt + eps >= lon[0]) && (lon_pt - eps <= lon[1])) { + return true; + } + return false; + } + + /// + /// Determine if this LatLonBox overlaps with box. + /// + bool overlaps( + const LatLonBox & box + ) const { + + // Check flags + if ((is_null) || (box.is_null)) { + return false; + } + if (lon_periodic != box.lon_periodic) { + _EXCEPTION2("Inconsistent periodicity flag (%s/%s)", + (lon_periodic)?("true"):("false"), + (box.lon_periodic)?("true"):("false")); + } + if (lon_width != box.lon_width) { + _EXCEPTIONT("Inconsistent box lon_width in comparison"); + } + + // Check latitudes + if (lat[0] > box.lat[1]) { + return false; + } + if (lat[1] < box.lat[0]) { + return false; + } + + // Cases when longitude is periodic + if (lon_periodic) { + + // Both boxes cross lon 360 + if ((lon[0] > lon[1]) && (box.lon[0] > box.lon[1])) { + return true; + } + + // This box crosses lon 360 + if (lon[0] > lon[1]) { + if (box.lon[1] >= lon[0]) { + return true; + } + if (box.lon[0] <= lon[1]) { + return true; + } + return false; + } + + // That box crosses lon 360 + if (box.lon[0] > box.lon[1]) { + if (lon[1] >= box.lon[0]) { + return true; + } + if (lon[0] <= box.lon[1]) { + return true; + } + return false; + } + } + + // No boxes cross lon 360 + if (box.lon[1] < lon[0]) { + return false; + } + if (box.lon[0] > lon[1]) { + return false; + } + return true; + } +}; + +/////////////////////////////////////////////////////////////////////////////// + +#endif // _LATLONBOX_H_ + diff --git a/src/Makefile.gmake b/src/Makefile.gmake index 8ea4ec3..1ea9334 100644 --- a/src/Makefile.gmake +++ b/src/Makefile.gmake @@ -15,6 +15,7 @@ include $(TEMPESTREMAPDIR)/mk/framework.make UTIL_FILES= Announce.cpp \ FiniteElementTools.cpp \ FiniteVolumeTools.cpp \ + FunctionTimer.cpp \ GaussLobattoQuadrature.cpp \ GaussQuadrature.cpp \ GridElements.cpp \ @@ -39,8 +40,9 @@ UTIL_FILES= Announce.cpp \ GenerateICOMesh.cpp \ GenerateLambertConfConicMesh.cpp \ GenerateOfflineMap.cpp \ - GenerateOverlapMesh.cpp \ - GenerateOverlapMesh_v1.cpp \ + GenerateOverlapMeshEdge.cpp \ + GenerateOverlapMeshKdx.cpp \ + GenerateOverlapMeshLint.cpp \ GenerateRLLMesh.cpp \ GenerateRectilinearMeshFromFile.cpp \ GenerateUTMMesh.cpp \ @@ -63,7 +65,7 @@ GenerateVolumetricMesh_FILES= GenerateVolumetricMesh.cpp # Overlap mesh generators GenerateOverlapMesh_FILES= GenerateOverlapMeshExe.cpp -GenerateOverlapMesh_v1_FILES= GenerateOverlapMeshExe_v1.cpp +GenerateOverlapMeshParallel_FILES= GenerateOverlapMeshParallelExe.cpp # Compute and apply the offline mapping weights ApplyOfflineMap_FILES= ApplyOfflineMapExe.cpp @@ -99,7 +101,7 @@ EXEC_TARGETS= ApplyOfflineMap \ GenerateLambertConfConicMesh \ GenerateOfflineMap \ GenerateOverlapMesh \ - GenerateOverlapMesh_v1 \ + GenerateOverlapMeshParallel \ GenerateRLLMesh \ GenerateRectilinearMeshFromFile \ GenerateUTMMesh \ @@ -132,7 +134,7 @@ GenerateICOMesh_EXE: $(GenerateICOMesh_FILES:%.cpp=$(BUILDDIR)/%.o) GenerateLambertConfConicMesh_EXE: $(GenerateLambertConfConicMesh_FILES:%.cpp=$(BUILDDIR)/%.o) GenerateVolumetricMesh_EXE: $(GenerateVolumetricMesh_FILES:%.cpp=$(BUILDDIR)/%.o) GenerateOverlapMesh_EXE: $(GenerateOverlapMesh_FILES:%.cpp=$(BUILDDIR)/%.o) -GenerateOverlapMesh_v1_EXE: $(GenerateOverlapMesh_v1_FILES:%.cpp=$(BUILDDIR)/%.o) +GenerateOverlapMeshParallel_EXE: $(GenerateOverlapMeshParallel_FILES:%.cpp=$(BUILDDIR)/%.o) ApplyOfflineMap_EXE: $(ApplyOfflineMap_FILES:%.cpp=$(BUILDDIR)/%.o) GenerateOfflineMap_EXE: $(GenerateOfflineMap_FILES:%.cpp=$(BUILDDIR)/%.o) GenerateGLLMetaData_EXE: $(GenerateGLLMetaData_FILES:%.cpp=$(BUILDDIR)/%.o) diff --git a/src/OfflineMap.h b/src/OfflineMap.h index b0dfb00..eaf2df7 100644 --- a/src/OfflineMap.h +++ b/src/OfflineMap.h @@ -494,15 +494,17 @@ class OfflineMap { /// /// Set the vector of areas associated with the source mesh. /// - void SetSourceAreas(const DataArray1D & dSourceAreas) { - m_dSourceAreas = dSourceAreas; + void SetSourceAreas(const std::vector & dSourceAreas) { + m_dSourceAreas.Allocate(dSourceAreas.size()); + memcpy(&(m_dSourceAreas[0]), &(dSourceAreas[0]), dSourceAreas.size() * sizeof(double)); } /// /// Set the vector of areas associated with the target mesh. /// - void SetTargetAreas(const DataArray1D & dTargetAreas) { - m_dTargetAreas = dTargetAreas; + void SetTargetAreas(const std::vector & dTargetAreas) { + m_dTargetAreas.Allocate(dTargetAreas.size()); + memcpy(&(m_dTargetAreas[0]), &(dTargetAreas[0]), dTargetAreas.size() * sizeof(double)); } public: diff --git a/src/OverlapMesh.cpp b/src/OverlapMesh.cpp index fcfee99..234bfc9 100644 --- a/src/OverlapMesh.cpp +++ b/src/OverlapMesh.cpp @@ -2,10 +2,10 @@ /// /// \file OverlapMesh.cpp /// \author Paul Ullrich -/// \version March 7, 2014 +/// \version August 19, 2021 /// /// -/// Copyright 2000-2014 Paul Ullrich +/// Copyright 2020 Paul Ullrich /// /// This file is distributed as part of the Tempest source code package. /// Permission is granted to use, copy, modify and distribute this @@ -22,10 +22,16 @@ #include "Announce.h" #include "kdtree.h" +#include "interval_tree.hpp" #include #include #include +#include + +#if defined(TEMPEST_MPIOMP) +#include +#endif /////////////////////////////////////////////////////////////////////////////// @@ -1165,7 +1171,7 @@ bool GenerateOverlapFaces( /////////////////////////////////////////////////////////////////////////////// -void GenerateOverlapMesh_v1( +void GenerateOverlapMeshEdge( const Mesh & meshSource, const Mesh & meshTarget, Mesh & meshOverlap, @@ -1395,10 +1401,19 @@ void GenerateOverlapFace( const EdgeVector & evecSource = faceSource.edges; // List outputList = subjectPolygon + nodevecOutput.reserve(evecSource.size() + evecTarget.size()); for (int i = 0; i < evecTarget.size(); i++) { nodevecOutput.push_back(nodesTarget[evecTarget[i][0]]); } + // Allocate space for intersections + std::vector vecIntersections; + vecIntersections.reserve(2); + + // Input node vector + NodeVector nodevecInput; + nodevecInput.reserve(evecSource.size() + evecTarget.size()); + // for (Edge clipEdge in clipPolygon) do for (int i = 0; i < evecSource.size(); i++) { @@ -1408,7 +1423,7 @@ void GenerateOverlapFace( } // List inputList = outputList; - NodeVector nodevecInput = nodevecOutput; + nodevecInput = nodevecOutput; // outputList.clear(); nodevecOutput.clear(); @@ -1446,7 +1461,7 @@ void GenerateOverlapFace( if (iNodeEdgeSideS < 0) { // outputList.add(ComputeIntersection(S,E,clipEdge)); - std::vector vecIntersections; + vecIntersections.clear(); bool fCoincident = utils.CalculateEdgeIntersectionsSemiClip( nodeS, @@ -1486,7 +1501,7 @@ void GenerateOverlapFace( } else if (iNodeEdgeSideS >= 0) { // outputList.add(ComputeIntersection(S,E,clipEdge)); - std::vector vecIntersections; + vecIntersections.clear(); bool fCoincident = utils.CalculateEdgeIntersectionsSemiClip( nodeS, @@ -1597,6 +1612,9 @@ int FindFaceContainingNode( queueTargetFaces.push(ixTargetFaceSeed); setExaminedTargetFaces.insert(ixTargetFaceSeed); + NodeVector nodevecOverlap; + nodevecOverlap.reserve(16); + int ixIterate = 0; while (!queueTargetFaces.empty()) { @@ -1628,7 +1646,7 @@ int FindFaceContainingNode( // Node on boundary of target face; check for overlap if (loc != Face::NodeLocation_Exterior) { - NodeVector nodevecOverlap; + nodevecOverlap.clear(); GenerateOverlapFace( meshSource, @@ -1883,13 +1901,14 @@ void GenerateOverlapMeshFromFace( meshOverlap.vecSourceFaceIx.push_back(ixSourceFace); meshOverlap.vecTargetFaceIx.push_back(ixCurrentTargetFace); + meshOverlap.vecFaceArea.push_back(dArea); } } } /////////////////////////////////////////////////////////////////////////////// -void GenerateOverlapMesh_v2( +void GenerateOverlapMeshKdx( const Mesh & meshSource, const Mesh & meshTarget, Mesh & meshOverlap, @@ -1897,7 +1916,7 @@ void GenerateOverlapMesh_v2( const bool fAllowNoOverlap, const bool fVerbose ) { -#if defined(OVERLAPMESH_RETAIN_REPEATED_NODES) +#if defined(OVERLAPMESH_USE_SORTED_MAP) NodeMap nodemapOverlap; #endif #if defined(OVERLAPMESH_USE_UNSORTED_MAP) @@ -1922,7 +1941,6 @@ void GenerateOverlapMesh_v2( // Generate Overlap mesh for each Face for (int i = 0; i < meshSource.faces.size(); i++) { - //for (int i = 3999555; i < 3999555+1; i++) { if (fVerbose) { std::string strAnnounce = "Source Face " + std::to_string((long long)i); AnnounceStartBlock(strAnnounce.c_str()); @@ -2025,7 +2043,480 @@ void GenerateOverlapMesh_v2( */ // Calculate Face areas //if (fVerbose) { - double dTotalAreaOverlap = meshOverlap.CalculateFaceAreas(false); + double dTotalAreaOverlap = meshOverlap.SumFaceAreas(); + Announce("Overlap Mesh Geometric Area: %1.15e (%1.15e)", dTotalAreaOverlap, 4.0 * M_PI); + //} +} + +/////////////////////////////////////////////////////////////////////////////// + +class RectangleEdge { + +public: + RectangleEdge( + double _value, + int _ix, + bool _begin, + bool _source + ) : + value(_value), + ix(_ix), + begin(_begin), + source(_source) + { } + + bool operator<(const RectangleEdge & redge) const { + // end edges always come before begin edges + if (value == redge.value) { + if (!begin && redge.begin) { + return true; + } + if (begin == redge.begin) { + return (ix < redge.ix); + } + return false; + } + return (value < redge.value); + } + +public: + double value; + int ix; + bool begin; + bool source; +}; + +/////////////////////////////////////////////////////////////////////////////// + +void GenerateOverlapMeshLint( + const Mesh & meshSource, + const Mesh & meshTarget, + Mesh & meshOverlap, + OverlapMeshMethod method, + const bool fParallel, + std::string strTempDir, + const bool fVerbose +) { + +#if defined(OVERLAPMESH_USE_SORTED_MAP) + NodeMap nodemapOverlap; +#endif +#if defined(OVERLAPMESH_USE_UNSORTED_MAP) + NodeMap nodemapOverlap; +#endif +#if defined(OVERLAPMESH_USE_NODE_MULTIMAP) + NodeMap nodemapOverlap(ReferenceTolerance, OVERLAPMESH_BIN_WIDTH); +#endif + + // Generate lat-lon rectangles for each face on the source and target mesh + std::vector< LatLonBox > vecMeshSourceLatLonBox; + std::vector< LatLonBox > vecMeshTargetLatLonBox; + + meshSource.GenerateLatLonBoxes(vecMeshSourceLatLonBox); + meshTarget.GenerateLatLonBoxes(vecMeshTargetLatLonBox); + + // Define interval trees + // TODO: Change to open intervals + lib_interval_tree::interval_tree_t inttreeSourceActive; + lib_interval_tree::interval_tree_t inttreeTargetActive; + + // Sort the longitude boundaries on the source and target mesh in order of increasing longitude + std::vector vecSortedLon; + vecSortedLon.reserve(2 * (vecMeshSourceLatLonBox.size() + vecMeshTargetLatLonBox.size())); + + for (int i = 0; i < vecMeshSourceLatLonBox.size(); i++) { + vecSortedLon.push_back({vecMeshSourceLatLonBox[i].lon[0], i, true, true}); + vecSortedLon.push_back({vecMeshSourceLatLonBox[i].lon[1], i, false, true}); + + _ASSERT(vecMeshSourceLatLonBox[i].lon[0] != vecMeshSourceLatLonBox[i].lon[1]); + if (vecMeshSourceLatLonBox[i].lon[0] > vecMeshSourceLatLonBox[i].lon[1]) { + inttreeSourceActive.insert({vecMeshSourceLatLonBox[i].lat[0], vecMeshSourceLatLonBox[i].lat[1], i}); + } + } + for (int i = 0; i < vecMeshTargetLatLonBox.size(); i++) { + vecSortedLon.push_back({vecMeshTargetLatLonBox[i].lon[0], i, true, false}); + vecSortedLon.push_back({vecMeshTargetLatLonBox[i].lon[1], i, false, false}); + + _ASSERT(vecMeshTargetLatLonBox[i].lon[0] != vecMeshTargetLatLonBox[i].lon[1]); + if (vecMeshTargetLatLonBox[i].lon[0] > vecMeshTargetLatLonBox[i].lon[1]) { + inttreeTargetActive.insert({vecMeshTargetLatLonBox[i].lat[0], vecMeshTargetLatLonBox[i].lat[1], i}); + } + } + + std::sort(vecSortedLon.begin(), vecSortedLon.end()); + _ASSERT(vecSortedLon.size() == 2 * (vecMeshSourceLatLonBox.size() + vecMeshTargetLatLonBox.size())); + + // Accumulated Face data (sort by source face ix) + struct LintFace { + Face face; + int ixSourceFace; + int ixTargetFace; + double dFaceArea; + + // Constructor + LintFace( + const Face & _face, + int _ixSourceFace, + int _ixTargetFace, + double _dFaceArea + ) : + face(_face), + ixSourceFace(_ixSourceFace), + ixTargetFace(_ixTargetFace), + dFaceArea(_dFaceArea) + { } + + // Comparator using source face ix + bool operator<(const LintFace & lf) const { + return (ixSourceFace < lf.ixSourceFace); + } + }; + + std::multiset setAccumulatedLintFaces; + + // Perform line search + std::set< std::pair > setTested; + + size_t iFace = 0; + + for (auto it : vecSortedLon) { + iFace++; + if (!fVerbose && ((iFace % 10000) == 0)) { + double dFrac = 100.0*static_cast(iFace)/static_cast(vecSortedLon.size()); + Announce("Percent Complete %2.1f%%", dFrac); + } + + // Vector containing indices of overlapped faces + std::vector< std::pair > vecOverlapIx; + + // Edge is of type source + if (it.source) { + if (it.begin) { +/* + for (int iTarget = 0; iTarget < vecMeshTargetLatLonBox.size(); iTarget++) { + if (vecMeshTargetLatLonBox[iTarget].overlaps(vecMeshSourceLatLonBox[it.ix])) { + vecOverlapIx.push_back({it.ix, iTarget}); + } + } +*/ + + auto itinterval = + inttreeSourceActive.insert( + {vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1], it.ix}); + _ASSERT(itinterval != inttreeSourceActive.end()); +/* + lib_interval_tree::interval intSource({vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1], it.ix}); + for (auto itTargetInterval : inttreeTargetActive) { + if (itTargetInterval.overlaps(intSource)) { + vecOverlapIx.push_back({it.ix, itTargetInterval.data()}); + } + } +*/ +/* + lib_interval_tree::interval intSource({vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1], it.ix}); + for (auto itTargetInterval : inttreeTargetActive) { + if (itTargetInterval.overlaps(intSource)) { + const LatLonBox & boxTarget = vecMeshTargetLatLonBox[itTargetInterval.data()]; + const LatLonBox & boxSource = vecMeshSourceLatLonBox[it.ix]; + + if (boxSource.overlaps(boxTarget)) { + vecOverlapIx.push_back({it.ix, itTargetInterval.data()}); + } + } + } +*/ +/* + inttreeTargetActive.overlap_find_all( + {vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1], it.ix}, + [&it, &vecOverlapIx](lib_interval_tree::interval_tree_t::iterator ittarget){ + vecOverlapIx.push_back({it.ix, ittarget->data()}); return true;}, + false); +*/ + inttreeTargetActive.overlap_find_all( + {vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1], it.ix}, + [&it, &vecOverlapIx, &vecMeshTargetLatLonBox, &vecMeshSourceLatLonBox]( + lib_interval_tree::interval_tree_t::iterator ittarget + ){ + const LatLonBox & boxSource = vecMeshSourceLatLonBox[it.ix]; + const LatLonBox & boxTarget = vecMeshTargetLatLonBox[ittarget->data()]; + + // Due to periodicity we may have tested this source/target pair already + if ((boxSource.lon[0] <= boxSource.lon[1]) || (boxSource.lon[1] <= boxTarget.lon[0])) { + vecOverlapIx.push_back({it.ix, ittarget->data()}); + } + + return true; + }, + true); +/* + for (auto itxx : vecOverlapIx) { + if ((itxx.first == 18) && (itxx.second == 40)) { + std::cout << "TESTA" << it.value << std::endl; + } + } +*/ +/* + std::sort(vecOverlapIx.begin(), vecOverlapIx.end()); + std::sort(vecOverlapIx2.begin(), vecOverlapIx2.end()); + if (vecOverlapIx != vecOverlapIx2) { + + printf("Contents of interval tree:\n"); + lib_interval_tree::interval_tree_t inttreeTargetCopy; + for (auto itTargetInterval : inttreeTargetActive) { + inttreeTargetCopy.insert(itTargetInterval); + printf("%1.15e %1.15e\n", itTargetInterval.low(), itTargetInterval.high()); + } + + std::vector< std::pair > vecOverlapIx3; + inttreeTargetCopy.overlap_find_all( + {vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1], it.ix}, + [&it, &vecOverlapIx3](lib_interval_tree::interval_tree_t::iterator ittarget){ + vecOverlapIx3.push_back({it.ix, ittarget->data()}); return true;}, + false); + printf("%lu %lu %lu\n", vecOverlapIx.size(), vecOverlapIx2.size(), vecOverlapIx3.size()); + + printf("Search interval:\n"); + printf("%1.15e %1.15e\n", vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1]); + for (auto it : vecOverlapIx) { + printf("A(%i,%i) %1.15e %1.15e\n", it.first, it.second, vecMeshTargetLatLonBox[it.second].lat[0], vecMeshTargetLatLonBox[it.second].lat[1]); + } + for (auto it : vecOverlapIx2) { + printf("B(%i,%i) %1.15e %1.15e\n", it.first, it.second, vecMeshTargetLatLonBox[it.second].lat[0], vecMeshTargetLatLonBox[it.second].lat[1]); + } + for (auto it : vecOverlapIx3) { + printf("C(%i,%i) %1.15e %1.15e\n", it.first, it.second, vecMeshTargetLatLonBox[it.second].lat[0], vecMeshTargetLatLonBox[it.second].lat[1]); + } + + _EXCEPTION(); + } +*/ + } else { + auto itinterval = + inttreeSourceActive.find( + {vecMeshSourceLatLonBox[it.ix].lat[0], vecMeshSourceLatLonBox[it.ix].lat[1], it.ix}); + _ASSERT(itinterval->data() == it.ix); + _ASSERT(itinterval != inttreeSourceActive.end()); + inttreeSourceActive.erase(itinterval); + } + + // Edge is of type target + } else { + if (it.begin) { + auto itinterval = + inttreeTargetActive.insert( + {vecMeshTargetLatLonBox[it.ix].lat[0], vecMeshTargetLatLonBox[it.ix].lat[1], it.ix}); + _ASSERT(itinterval != inttreeTargetActive.end()); +/* + lib_interval_tree::interval intTarget({vecMeshTargetLatLonBox[it.ix].lat[0], vecMeshTargetLatLonBox[it.ix].lat[1], it.ix}); + for (auto itSourceInterval : inttreeSourceActive) { + if (itSourceInterval.overlaps(intTarget)) { + const LatLonBox & boxTarget = vecMeshTargetLatLonBox[it.ix]; + const LatLonBox & boxSource = vecMeshSourceLatLonBox[itSourceInterval.data()]; + + // Due to periodicity we may have tested this source/target pair already + if (boxTarget.lon[0] > boxTarget.lon[1]) { + if (boxTarget.lon[1] > boxSource.lon[0]) { + continue; + } + } + + vecOverlapIx.push_back({itSourceInterval.data(), it.ix}); + } + } +*/ + + inttreeSourceActive.overlap_find_all( + {vecMeshTargetLatLonBox[it.ix].lat[0], vecMeshTargetLatLonBox[it.ix].lat[1], it.ix}, + [&it, &vecOverlapIx, &vecMeshTargetLatLonBox, &vecMeshSourceLatLonBox]( + lib_interval_tree::interval_tree_t::iterator itsource + ){ + const LatLonBox & boxTarget = vecMeshTargetLatLonBox[it.ix]; + const LatLonBox & boxSource = vecMeshSourceLatLonBox[itsource->data()]; + + // Due to periodicity we may have tested this source/target pair already + if ((boxTarget.lon[0] <= boxTarget.lon[1]) || (boxTarget.lon[1] <= boxSource.lon[0])) { + vecOverlapIx.push_back({itsource->data(), it.ix}); + } + + return true; + }, + true); +/* + for (auto itxx : vecOverlapIx) { + if ((itxx.first == 18) && (itxx.second == 40)) { + std::cout << "TESTB " << it.value << std::endl; + } + } +*/ + } else { + auto itinterval = + inttreeTargetActive.find( + {vecMeshTargetLatLonBox[it.ix].lat[0], vecMeshTargetLatLonBox[it.ix].lat[1], it.ix}); + _ASSERT(itinterval->data() == it.ix); + _ASSERT(itinterval != inttreeTargetActive.end()); + inttreeTargetActive.erase(itinterval); + } + } +/* + for (auto itoverlap : vecOverlapIx) { + auto itTested = setTested.find(itoverlap); + if (itTested != setTested.end()) { + printf("\nERROR: Two faces already used for overlaps are being used again:\n"); + printf("A: %i B: %i\n", itTested->first, itTested->second); + + const LatLonBox & boxSource = vecMeshSourceLatLonBox[itTested->first]; + const LatLonBox & boxTarget = vecMeshTargetLatLonBox[itTested->second]; + + printf("A: [%1.5f %1.5f] x [%1.5f %1.5f]\n", + boxSource.lat[0], boxSource.lat[1], boxSource.lon[0], boxSource.lon[1]); + printf("B: [%1.5f %1.5f] x [%1.5f %1.5f]\n", + boxTarget.lat[0], boxTarget.lat[1], boxTarget.lon[0], boxTarget.lon[1]); + + printf("Value (%1.5f) Index (%i) Source (%i) Begin (%i)\n", + it.value, static_cast(it.ix), static_cast(it.source), static_cast(it.begin)); + _EXCEPTION(); + } + setTested.insert(itoverlap); + } +*/ + // Find the overlap polygons + for (auto itoverlap : vecOverlapIx) { + + _ASSERT(vecMeshSourceLatLonBox[itoverlap.first].overlaps(vecMeshTargetLatLonBox[itoverlap.second])); + + NodeVector nodevecOutput; + + GenerateOverlapFace( + meshSource, + meshTarget, + itoverlap.first, + itoverlap.second, + nodevecOutput + ); + + if (nodevecOutput.size() == 0) { + continue; + + } else if (nodevecOutput.size() < 3) { + _EXCEPTIONT("Overlap polygon consists of fewer than 3 nodes"); + } + + // Calculate face area + Face faceTemp(nodevecOutput.size()); + for (int i = 0; i < nodevecOutput.size(); i++) { + faceTemp.SetNode(i, i); + } + Real dArea = CalculateFaceArea(faceTemp, nodevecOutput); + + if (dArea < 1.0e-13) { + continue; + } + + // Insert Face into meshOverlap, using nodemapOverlap to remove redundant nodes + Face faceNew(nodevecOutput.size()); + for (int i = 0; i < nodevecOutput.size(); i++) { +#if defined(OVERLAPMESH_RETAIN_REPEATED_NODES) + faceNew.SetNode(i, meshOverlap.nodes.size()); + meshOverlap.nodes.push_back(nodevecOutput[i]); +#else + NodeMapConstIterator iter + = nodemapOverlap.find(nodevecOutput[i]); + + if (iter != nodemapOverlap.end()) { + faceNew.SetNode(i, iter->second); + } else { + int iNextNodeMapOverlapIx = nodemapOverlap.size(); + faceNew.SetNode(i, iNextNodeMapOverlapIx); + nodemapOverlap.insert( + NodeMapPair(nodevecOutput[i], iNextNodeMapOverlapIx)); + } +#endif + + } + + setAccumulatedLintFaces.insert( + LintFace(faceNew, itoverlap.first, itoverlap.second, dArea)); + + //meshOverlap.faces.push_back(faceNew); + //meshOverlap.vecSourceFaceIx.push_back(itoverlap.first); + //meshOverlap.vecTargetFaceIx.push_back(itoverlap.second); + //meshOverlap.vecFaceArea.push_back(dArea); + } + } + + // Sort mesh faces by source face ix + meshOverlap.faces.reserve(setAccumulatedLintFaces.size()); + meshOverlap.vecSourceFaceIx.reserve(setAccumulatedLintFaces.size()); + meshOverlap.vecTargetFaceIx.reserve(setAccumulatedLintFaces.size()); + meshOverlap.vecFaceArea.reserve(setAccumulatedLintFaces.size()); + for (auto it : setAccumulatedLintFaces) { + meshOverlap.faces.push_back(it.face); + meshOverlap.vecSourceFaceIx.push_back(it.ixSourceFace); + meshOverlap.vecTargetFaceIx.push_back(it.ixTargetFace); + meshOverlap.vecFaceArea.push_back(it.dFaceArea); + } + + // Replace parent indices if meshSource has a MultiFaceMap + if (meshSource.vecMultiFaceMap.size() != 0) { + for (int f = 0; f < meshOverlap.faces.size(); f++) { + meshOverlap.vecSourceFaceIx[f] = + meshSource.vecMultiFaceMap[meshOverlap.vecSourceFaceIx[f]]; + } + } + + // Replace parent indices if meshTarget has a MultiFaceMap + if (meshTarget.vecMultiFaceMap.size() != 0) { + for (int f = 0; f < meshOverlap.faces.size(); f++) { + meshOverlap.vecTargetFaceIx[f] = + meshSource.vecMultiFaceMap[meshOverlap.vecTargetFaceIx[f]]; + } + } + + // Replace parent indices if meshSource has a MultiFaceMap + if (meshSource.vecMultiFaceMap.size() != 0) { + for (int f = 0; f < meshOverlap.faces.size(); f++) { + meshOverlap.vecSourceFaceIx[f] = + meshSource.vecMultiFaceMap[meshOverlap.vecSourceFaceIx[f]]; + } + } + + // Replace parent indices if meshTarget has a MultiFaceMap + if (meshTarget.vecMultiFaceMap.size() != 0) { + for (int f = 0; f < meshOverlap.faces.size(); f++) { + meshOverlap.vecTargetFaceIx[f] = + meshSource.vecMultiFaceMap[meshOverlap.vecTargetFaceIx[f]]; + } + } + +#if !defined(OVERLAPMESH_RETAIN_REPEATED_NODES) + // Insert all Nodes from nodemapOverlap into meshOverlap.nodes + meshOverlap.nodes.resize(nodemapOverlap.size()); + + NodeMapConstIterator iter = nodemapOverlap.begin(); + for (; iter != nodemapOverlap.end(); iter++) { + meshOverlap.nodes[iter->second] = iter->first; + } +#endif + +/* +#if !defined(OVERLAPMESH_RETAIN_REPEATED_NODES) + meshOverlap.RemoveCoincidentNodes(false); +#endif +*/ +/* + // Check concavity of overlap mesh + AnnounceStartBlock("Testing concavity of overlap mesh"); + for (int i = 0; i < meshOverlap.faces.size(); i++) { + bool fIsConcave = meshOverlap.IsFaceConcave(i); + if (fIsConcave) { + _EXCEPTIONT("Concave element detected in overlap mesh"); + } + } + AnnounceEndBlock("Done"); +*/ + // Calculate Face areas + //if (fVerbose) { + double dTotalAreaOverlap = meshOverlap.SumFaceAreas(); Announce("Overlap Mesh Geometric Area: %1.15e (%1.15e)", dTotalAreaOverlap, 4.0 * M_PI); //} } diff --git a/src/OverlapMesh.h b/src/OverlapMesh.h index 3d22384..859c218 100644 --- a/src/OverlapMesh.h +++ b/src/OverlapMesh.h @@ -2,10 +2,10 @@ /// /// \file OverlapMesh.h /// \author Paul Ullrich -/// \version March 7, 2014 +/// \version August 19, 2021 /// /// -/// Copyright 2000-2014 Paul Ullrich +/// Copyright 2021 Paul Ullrich /// /// This file is distributed as part of the Tempest source code package. /// Permission is granted to use, copy, modify and distribute this @@ -38,12 +38,12 @@ enum OverlapMeshMethod { /// Generate the mesh obtained by overlapping meshes meshSource and /// meshTarget. /// -void GenerateOverlapMesh_v1( +void GenerateOverlapMeshEdge( const Mesh & meshSource, const Mesh & meshTarget, Mesh & meshOverlap, OverlapMeshMethod method, - const bool verbose = true + const bool fVerbose = true ); /////////////////////////////////////////////////////////////////////////////// @@ -69,13 +69,29 @@ void GenerateOverlapFace( /// Generate the mesh obtained by overlapping meshes meshSource and /// meshTarget. /// -void GenerateOverlapMesh_v2( +void GenerateOverlapMeshKdx( const Mesh & meshSource, const Mesh & meshTarget, Mesh & meshOverlap, OverlapMeshMethod method, const bool fAllowNoOverlap, - const bool verbose = true + const bool fVerbose = true +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Generate the mesh obtained by overlapping meshes meshSource and +/// meshTarget. +/// +void GenerateOverlapMeshLint( + const Mesh & meshSource, + const Mesh & meshTarget, + Mesh & meshOverlap, + OverlapMeshMethod method, + const bool fParallel, + std::string strTempDir, + const bool fVerbose = true ); /////////////////////////////////////////////////////////////////////////////// diff --git a/src/TempestRemapAPI.h b/src/TempestRemapAPI.h index 9933d55..ac8f22c 100644 --- a/src/TempestRemapAPI.h +++ b/src/TempestRemapAPI.h @@ -172,17 +172,57 @@ extern "C" { bool fAllowNoOverlap = false, bool fVerbose = true ); - // Old version of the implementation to compute the overlap mesh - // given a source and target mesh file names - int GenerateOverlapMesh_v1 ( + /// + /// Old version of the implementation to compute the overlap mesh + /// given a source and target mesh file names. + /// + int GenerateOverlapMeshEdge ( std::string strMeshA, std::string strMeshB, Mesh& meshOverlap, std::string strOverlapMesh, + std::string strOutputFormat, std::string strMethod, const bool fNoValidate = true ); - // Generate the Gauss-Lobatto-Legendre metadata for the given Mesh + /// + /// Compute the overlap mesh given a source and target mesh file names + /// using the kd-tree algorithm. + /// + int GenerateOverlapMeshKdx ( + std::string strMeshA, + std::string strMeshB, + Mesh & meshOverlap, + std::string strOverlapMesh, + std::string strOutputFormat, + std::string strMethod, + bool fNoValidate, + bool fHasConcaveFacesA = false, + bool fHasConcaveFacesB = false, + bool fAllowNoOverlap = false, + bool fVerbose = true ); + + /// + /// Compute the overlap mesh given a source and target mesh file names + /// using the line sweep algorithm. + /// + int GenerateOverlapMeshLint ( + std::string strMeshA, + std::string strMeshB, + Mesh & meshOverlap, + std::string strOverlapMesh, + std::string strOutputFormat, + std::string strMethod, + bool fNoValidate, + bool fHasConcaveFacesA = false, + bool fHasConcaveFacesB = false, + bool fParallel = false, + std::string strTempDir = "/tmp", + bool fVerbose = true ); + + /// + /// Generate the Gauss-Lobatto-Legendre metadata for the given Mesh. + /// int GenerateGLLMetaData ( std::string strMesh, Mesh & meshOut, @@ -438,7 +478,6 @@ extern "C" { int GenerateConnectivityData( const Mesh & meshIn, std::vector< std::set > & vecConnectivity ); - } #endif // TEMPESTREMAP_API_H diff --git a/src/interval_tree.hpp b/src/interval_tree.hpp new file mode 100644 index 0000000..1e20ae6 --- /dev/null +++ b/src/interval_tree.hpp @@ -0,0 +1,1603 @@ +// This is free and unencumbered software released into the public domain. +// +// Anyone is free to copy, modify, publish, use, compile, sell, or +// distribute this software, either in source code form or as a compiled +// binary, for any purpose, commercial or non-commercial, and by any +// means. +// +// In jurisdictions that recognize copyright laws, the author or authors +// of this software dedicate any and all copyright interest in the +// software to the public domain. We make this dedication for the benefit +// of the public at large and to the detriment of our heirs and +// successors. We intend this dedication to be an overt act of +// relinquishment in perpetuity of all present and future rights to this +// software under copyright law. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +// IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR +// OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// +// For more information, please refer to +// +// Code obtained from: +// Ebbeke, Tim (2021) "interval-tree" https://github.com/5cript/interval-tree, +// accessed August 22, 2021 commit f5bff66f88daf47b03a5a612c58f09e2a399c199. +// + +#pragma once + +#include "interval_tree_fwd.hpp" +#include "interval_types.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace lib_interval_tree +{ +//############################################################################################################ + enum class rb_color + { + fail, + red, + black, + double_black + }; +//############################################################################################################ + using default_interval_value_type = int; + using default_interval_data_type = int; +//############################################################################################################ + template + struct interval + { + public: + using value_type = numerical_type; + using interval_kind = interval_kind_; + using data_type = data_type_; + + /** + * Constructs an interval. low MUST be smaller than high. + */ +#ifndef INTERVAL_TREE_SAFE_INTERVALS + interval(value_type low, value_type high, data_type data) + : low_{low} + , high_{high} + , data_{data} + { + assert(low <= high); + } +#else + interval(value_type low, value_type high, data_type data) + : low_{std::min(low, high)} + , high_{std::max(low, high)} + , data_{data} + { + } +#endif + + /** + * Returns if both intervals equal. + */ + friend bool operator==(interval const& lhs, interval const& other) + { + return lhs.low_ == other.low_ && lhs.high_ == other.high_ && lhs.data_ == other.data_; + } + + /** + * Returns if both intervals are different. + */ + friend bool operator!=(interval const& lhs, interval const& other) + { + return lhs.low_ != other.low_ || lhs.high_ != other.high_ || lhs.data_ != other.data_; + } + + /** + * Returns the lower bound of the interval + */ + value_type low() const + { + return low_; + } + + /** + * Returns the upper bound of the interval + */ + value_type high() const + { + return high_; + } + + /** + * Get the data associated with this interval. + */ + const data_type & data() const { + return data_; + } + + /** + * Returns whether the intervals overlap. + * For when both intervals are closed. + */ + bool overlaps(value_type l, value_type h) const + { + return low_ <= h && l <= high_; + } + + /** + * Returns whether the intervals overlap, excluding border. + * For when at least one interval is open (l&r). + */ + bool overlaps_exclusive(value_type l, value_type h) const + { + return low_ < h && l < high_; + } + + /** + * Returns whether the intervals overlap + */ + bool overlaps(interval const& other) const + { + return overlaps(other.low_, other.high_); + } + + /** + * Returns whether the intervals overlap, excluding border. + */ + bool overlaps_exclusive(interval const& other) const + { + return overlaps_exclusive(other.low_, other.high_); + } + + /** + * Returns whether the given value is in this. + */ + bool within(value_type value) const + { + return interval_kind::within(low_, high_, value); + } + + /** + * Returns whether the given interval is in this. + */ + bool within(interval const& other) const + { + return low_ <= other.low_ && high_ >= other.high_; + } + + /** + * Calculates the distance between the two intervals. + * Overlapping intervals have 0 distance. + */ + value_type operator-(interval const& other) const + { + if (overlaps(other)) + return 0; + if (high_ < other.low_) + return other.low_ - high_; + else + return low_ - other.high_; + } + + /** + * Returns the size of the interval. + */ + value_type size() const + { + return high_ - low_; + } + + /** + * Creates a new interval from this and other, that contains both intervals and whatever + * is between. + */ + interval join(interval const& other) const + { + return {std::min(low_, other.low_), std::max(high_, other.high_)}; + } + + private: + value_type low_; + value_type high_; + data_type data_; + }; +//############################################################################################################ + /** + * Creates a safe interval that puts the lower bound left automatically. + */ + template + interval make_safe_interval(numerical_type lhs, numerical_type rhs) + { + return interval {std::min(lhs, rhs), std::max(lhs, rhs)}; + } +//############################################################################################################ + template > + class node + { + private: + using node_type = node ; + + public: + using interval_type = interval_type_; + using value_type = numerical_type; + + public: + friend lib_interval_tree::interval_tree ; + friend lib_interval_tree::const_interval_tree_iterator >; + friend lib_interval_tree::interval_tree_iterator >; + + public: + node(node* parent, interval_type interval) + : interval_{std::move(interval)} + , max_{interval.high()} + , parent_{parent} + , left_{} + , right_{} + , color_{rb_color::fail} + { + } + + ~node() + { + } + + interval_type interval() const + { + return interval_; + } + + value_type max() const + { + return max_; + } + + bool is_left() const noexcept + { + return this == parent_->left_; + } + + bool is_right() const noexcept + { + return this == parent_->right_; + } + + bool is_root() const noexcept + { + return !parent_; + } + + /** + * Returns the color of the node. + */ + rb_color color() const + { + return color_; + } + + /** + * Returns the parent node up the tree. + */ + node const* parent() const + { + return parent_; + } + + /** + * Returns the left node (readonly). + */ + node const* left() const + { + return left_; + } + + /** + * Returns the right node (readonly). + */ + node const* right() const + { + return right_; + } + + /** + * Returns the height of the node in the tree. Where height = how many parents does it have. + * The root has no parents and is therefor has height 0. + */ + int height() const + { + int counter{0}; + for (auto* p = parent_; p != nullptr; p = p->parent_) + ++counter; + return counter; + } + + /** + * Returns the lower bound of the interval of this node + */ + value_type low() const + { + return interval_.low(); + } + + /** + * Returns the upper bound of the interval of this node + */ + value_type high() const + { + return interval_.high(); + } + + /** + * Returns the data of the interval of this node + */ + data_type data() const + { + return interval_.data(); + } + +private: + void set_interval(interval_type const& ival) + { + interval_ = ival; + } + + void kill() const noexcept + { + auto* parent = parent_; + if (is_left()) + { + delete parent_->left_; + parent->left_ = nullptr; + } + else + { + delete parent_->right_; + parent->right_ = nullptr; + } + } + + private: + interval_type interval_; + value_type max_; + node* parent_; + node* left_; + node* right_; + rb_color color_; + }; +//############################################################################################################ + template + class basic_interval_tree_iterator : public std::forward_iterator_tag + { + public: + friend interval_tree ; + + using tree_type = interval_tree ; + using value_type = node_type; + + using node_ptr_t = typename std::conditional < + std::is_const ::type>::value, + node_type const*, + node_type* + >::type; + + public: + constexpr basic_interval_tree_iterator(basic_interval_tree_iterator const&) = default; + basic_interval_tree_iterator& operator=(basic_interval_tree_iterator const&) = default; + + bool operator!=(basic_interval_tree_iterator const& other) const + { + return node_ != other.node_; + } + + bool operator==(basic_interval_tree_iterator const& other) const + { + return node_ == other.node_; + } + + /** + * Returns the max property of the node. + */ + typename node_type::interval_type::value_type max() const + { + return node_->max(); + } + + /** + * Returns the color of the node. + */ + rb_color color() const + { + return node_->color(); + } + + typename tree_type::interval_type interval() const + { + return node_->interval(); + } + + virtual ~basic_interval_tree_iterator() = default; + + protected: + basic_interval_tree_iterator(node_ptr_t node, owner_type owner) + : node_{node} + , owner_{owner} + { + } + + protected: + node_ptr_t node_; + owner_type owner_; + }; +//############################################################################################################ + template + class const_interval_tree_iterator + : public basic_interval_tree_iterator const*> + { + public: + using tree_type = interval_tree ; + using iterator_base = basic_interval_tree_iterator ; + using value_type = typename iterator_base::value_type; + using iterator_base::node_; + using iterator_base::owner_; + + friend tree_type; + + public: + const_interval_tree_iterator& operator++() + { + if (!node_) + { + node_ = owner_->root_; + + if (!node_) + return *this; + + while(node_->left_) + node_ = node_->left_; + } + + if (node_->right_) + { + node_ = node_->right_; + + while (node_->left_) + node_ = node_->left_; + } + else + { + auto* parent = node_->parent_; + while (parent != nullptr && node_ == parent->right_) + { + node_ = parent; + parent = parent->parent_; + } + node_ = parent; + } + + return *this; + } + + const_interval_tree_iterator operator++(int) + { + const_interval_tree_iterator cpy = *this; + operator++(); + return cpy; + } + + typename value_type::interval_type operator*() const + { + if (node_) + return node_->interval(); + else + throw std::out_of_range("dereferencing interval_tree_iterator out of bounds"); + } + + /** + * Returns an iterator to the parent of this node. + * will equal std::end(tree) if there is no parent node. + */ + const_interval_tree_iterator parent() const + { + if (node_) + return {node_->parent_, owner_}; + else + throw std::out_of_range("interval_tree_iterator out of bounds"); + } + + /** + * Continues down the left side of this node. + * will equal std::end(tree) if there is no left node. + */ + const_interval_tree_iterator left() const + { + if (node_) + return {node_->left_, owner_}; + else + throw std::out_of_range("interval_tree_iterator out of bounds"); + } + + /** + * Continues down the right side of this node. + * will equal std::end(tree) if there is no right node. + */ + const_interval_tree_iterator right() const + { + if (node_) + return {node_->right_, owner_}; + else + throw std::out_of_range("interval_tree_iterator out of bounds"); + } + + value_type const* operator->() const + { + return node_; + } + + private: + const_interval_tree_iterator(node_type const* node, tree_type const* owner) + : basic_interval_tree_iterator {node, owner} + { + } + }; +//############################################################################################################ + template + class interval_tree_iterator + : public basic_interval_tree_iterator *> + { + public: + using tree_type = interval_tree ; + using iterator_base = basic_interval_tree_iterator ; + using value_type = typename iterator_base::value_type; + using iterator_base::node_; + using iterator_base::owner_; + + friend tree_type; + + public: + interval_tree_iterator& operator++() + { + if (!node_) + { + node_ = owner_->root_; + + if (!node_) + return *this; + + while(node_->left_) + node_ = node_->left_; + } + + if (node_->right_) + { + node_ = node_->right_; + + while (node_->left_) + node_ = node_->left_; + } + else + { + auto* parent = node_->parent_; + while (parent != nullptr && node_ == parent->right_) + { + node_ = parent; + parent = parent->parent_; + } + node_ = parent; + } + + return *this; + } + + interval_tree_iterator operator++(int) + { + interval_tree_iterator cpy = *this; + operator++(); + return cpy; + } + + /** + * Returns an iterator to the parent of this node. + * will equal std::end(tree) if there is no parent node. + */ + interval_tree_iterator parent() + { + if (node_) + return {node_->parent_, owner_}; + else + throw std::out_of_range("interval_tree_iterator out of bounds"); + } + + /** + * Continues down the left side of this node. + * will equal std::end(tree) if there is no left node. + */ + interval_tree_iterator left() + { + if (node_) + return {node_->left_, owner_}; + else + throw std::out_of_range("interval_tree_iterator out of bounds"); + } + + /** + * Continues down the right side of this node. + * will equal std::end(tree) if there is no right node. + */ + interval_tree_iterator right() + { + if (node_) + return {node_->right_, owner_}; + else + throw std::out_of_range("interval_tree_iterator out of bounds"); + } + + typename value_type::interval_type operator*() const + { + if (node_) + return node_->interval(); + else + throw std::out_of_range("interval_tree_iterator out of bounds"); + } + + value_type* operator->() + { + return node_; + } + + private: + interval_tree_iterator(node_type* node, tree_type* owner) + : basic_interval_tree_iterator {node, owner} + { + } + }; +//############################################################################################################ + template > + class interval_tree + { + public: + using interval_type = IntervalT; + using value_type = typename interval_type::value_type; + using data_type = typename interval_type::data_type; + using node_type = node ; + using iterator = interval_tree_iterator ; + using const_iterator = const_interval_tree_iterator ; + using size_type = long long; + + public: + friend const_interval_tree_iterator ; + friend interval_tree_iterator ; + + interval_tree() + : root_{nullptr} + , size_{0} + { + } + + ~interval_tree() + { + clear(); + } + + interval_tree(interval_tree const& other) + : root_{nullptr} + , size_{0} + { + operator=(other); + } + + public: + interval_tree& operator=(interval_tree const& other) + { + if (!empty()) + clear(); + + if (other.root_ != nullptr) + root_ = copyTreeImpl(other.root_, nullptr); + + size_ = other.size_; + + return *this; + } + + /** + * Removes all from this tree. + */ + void clear() + { + for (auto i = std::begin(*this); i != std::end(*this);) + i = erase(i); + } + + /** + * Returns the root node from this tree. + */ + iterator root() + { + return {root_, this}; + } + + /** + * Returns the root node from this tree. + */ + const_iterator root() const + { + return {root_, this}; + } + + /** + * Inserts an interval into the tree. + */ + template + iterator insert(IntervalType&& ival) + { + node_type* z = new node_type(nullptr, std::forward (ival)); + node_type* y = nullptr; + node_type* x = root_; + while (x) + { + y = x; + if (z->interval_.low() < x->interval_.low()) + x = x->left_; + else + x = x->right_; + } + z->parent_ = y; + if (!y) + root_ = z; + else if (z->interval_.low() < y->interval_.low()) + y->left_ = z; + else + y->right_ = z; + z->color_ = rb_color::red; + + insert_fixup(z); + recalculate_max(z); + ++size_; + return {z, this}; + } + + /** + * Inserts an interval into the tree if no other interval overlaps it. + * Otherwise merge the interval with the being overlapped. + * + * @param ival The interval + * @param exclusive Exclude borders. + */ + iterator insert_overlap(interval_type const& ival, bool exclusive = false) + { + auto iter = overlap_find(ival, exclusive); + if (iter == end()) + return insert(ival); + else + { + auto merged = iter->interval().join(ival); + erase(iter); + return insert(merged); + } + } + + /** + * Erases the element pointed to be iter. + */ + iterator erase(iterator iter) + { + if (!iter.node_) + throw std::out_of_range("cannot erase end iterator"); + + auto next = iter; + ++next; + + node_type* y; + if (!iter.node_->left_ || !iter.node_->right_) + y = iter.node_; + else + y = successor(iter.node_); + + node_type* x; + if (y->left_) + x = y->left_; + else + x = y->right_; + + if (x) + x->parent_ = y->parent_; + + auto* x_parent = y->parent_; + + if (!y->parent_) + root_ = x; + else if (y->is_left()) + y->parent_->left_ = x; + else + y->parent_->right_ = x; + + if (y != iter.node_) + { + iter.node_->interval_ = y->interval_; + iter.node_->max_ = y->max_; + recalculate_max(iter.node_); + } + + if (x && x->color_ == rb_color::red) + { + if (x_parent) + erase_fixup(x, x_parent, y->is_left()); + else + x->color_ = rb_color::black; + } + + delete y; + + --size_; + return next; + } + + /** + * Returns the size of the object. + */ + size_type size() const + { + return size_; + } + + /** + * Finds the first exact match. + * + * @param ival The interval to find an exact match for within the tree. + * @param compare A comparison function to use. + */ + template + iterator find(interval_type const& ival, CompareFunctionT const& compare) + { + if (root_ == nullptr) + return end(); + return iterator{find_i(root_, ival, compare), this}; + } + + /** + * Finds the first exact match. + * + * @param ival The interval to find an exact match for within the tree. + * @param compare A comparison function to use. + */ + template + const_iterator find(interval_type const& ival, CompareFunctionT const& compare) const + { + if (root_ == nullptr) + return end(); + return const_iterator{find_i(root_, ival, compare), this}; + } + + /** + * Finds the first exact match. + * + * @param ival The interval to find an exact match for within the tree. + */ + iterator find(interval_type const& ival) + { + return find(ival, [](interval_type const& lhs, interval_type const& rhs){return lhs == rhs;}); + } + /** + * Finds the first exact match. + * + * @param ival The interval to find an exact match for within the tree. + */ + const_iterator find(interval_type const& ival) const + { + return find(ival, [](interval_type const& lhs, interval_type const& rhs){return lhs == rhs;}); + } + + /** + * Finds all exact matches and returns the amount of intervals found. + */ + template + void find_all(interval_type const& ival, FunctionT const& on_find, CompareFunctionT const& compare) + { + if (root_ == nullptr) + return; + find_all_i(root_, ival, on_find, compare); + } + template + void find_all(interval_type const& ival, FunctionT const& on_find, CompareFunctionT const& compare) const + { + if (root_ == nullptr) + return; + find_all_i(root_, ival, on_find, compare); + } + + template + void find_all(interval_type const& ival, FunctionT const& on_find) + { + find_all(ival, on_find, [](interval_type const& lhs, interval_type const& rhs){return lhs == rhs;}); + } + template + void find_all(interval_type const& ival, FunctionT const& on_find) const + { + find_all(ival, on_find, [](interval_type const& lhs, interval_type const& rhs){return lhs == rhs;}); + } + + /** + * Finds the next exact match EXCLUDING from. + * + * @param from The iterator to search from EXCLUSIVE! + * @param ival The interval to find an exact match for within the tree. + * @param compare A comparison function to use. + */ + template + iterator find_next_in_subtree(iterator from, interval_type const& ival, CompareFunctionT const& compare) + { + if (root_ == nullptr) + return end(); + return iterator{find_i_ex(from.node_, ival, compare), this}; + } + template + const_iterator find_next_in_subtree(iterator from, interval_type const& ival, CompareFunctionT const& compare) const + { + if (root_ == nullptr) + return end(); + return iterator{find_i_ex(from.node_, ival, compare), this}; + } + + /** + * Finds the next exact match EXCLUDING from. + * + * @param from The iterator to search from, EXCLUSIVE! + * @param ival The interval to find an exact match for within the tree. + * @param compare A comparison function to use. + */ + iterator find_next_in_subtree(iterator from, interval_type const& ival) + { + return find_next_in_subtree(from, ival, [](interval_type const& lhs, interval_type const& rhs){return lhs == rhs;}); + } + const_iterator find_next_in_subtree(iterator from, interval_type const& ival) const + { + return find_next_in_subtree(from, ival, [](interval_type const& lhs, interval_type const& rhs){return lhs == rhs;}); + } + + /** + * Finds the first interval that overlaps with ival. + * + * @param ival The interval to find an overlap for within the tree. + * @param exclusive Exclude edges? + */ + iterator overlap_find(interval_type const& ival, bool exclusive = false) + { + if (root_ == nullptr) + return end(); + if (exclusive) + return iterator{overlap_find_i(root_, ival), this}; + else + return iterator{overlap_find_i(root_, ival), this}; + } + const_iterator overlap_find(interval_type const& ival, bool exclusive = false) const + { + if (root_ == nullptr) + return end(); + if (exclusive) + return const_iterator{overlap_find_i(root_, ival), this}; + else + return const_iterator{overlap_find_i(root_, ival), this}; + } + + /** + * Finds all intervals that overlaps with ival. + * + * @param ival The interval to find an overlap for within the tree. + * @param exclusive Exclude edges? + */ + template + void overlap_find_all(interval_type const& ival, FunctionT const& on_find, bool exclusive = false) + { + if (root_ == nullptr) + return; + if (exclusive) + overlap_find_all_i(root_, ival, on_find); + else + overlap_find_all_i(root_, ival, on_find); + } + template + void overlap_find_all(interval_type const& ival, FunctionT const& on_find, bool exclusive = false) const + { + if (root_ == nullptr) + return; + if (exclusive) + overlap_find_all_i(root_, ival, on_find); + else + overlap_find_all_i(root_, ival, on_find); + } + + /** + * Finds the next interval that overlaps with ival + * + * @param from The iterator to start from, EXCLUSIVE! + * @param ival The interval to find an overlap for within the tree. + * @param exclusive Exclude edges? + */ + iterator overlap_find_next_in_subtree(iterator from, interval_type const& ival, bool exclusive = false) + { + if (root_ == nullptr) + return end(); + return iterator{overlap_find_i_ex(from.node_, ival, exclusive), this}; + } + const_iterator overlap_find_next_in_subtree(const_iterator from, interval_type const& ival, bool exclusive = false) const + { + if (root_ == nullptr) + return end(); + return const_iterator {overlap_find_i_ex(from.node_, ival, exclusive), this}; + } + + /** + * Deoverlaps the tree but returns it as a copy. + */ + interval_tree deoverlap_copy() + { + interval_tree fresh; + for (auto i = begin(), e = end(); i != e; ++i) + fresh.insert_overlap(*i); + + return fresh; + } + + /** + * Merges all overlapping intervals by erasing overlapping intervals and reinserting the merged interval. + */ + interval_tree& deoverlap() + { + *this = deoverlap_copy(); + return *this; + } + + /** + * Only works with deoverlapped trees. + * Creates an interval tree that contains all gaps between the intervals as intervals. + */ + interval_tree punch() const + { + if (empty()) + return {}; + auto min = std::begin(*this)->interval().low(); + auto max = root_->max_; + return punch({min, max}); + } + + /** + * Only works with deoverlapped trees. + * Removes all intervals from the given interval and produces a tree that contains the remaining intervals. + * This is basically the other punch overload with ival = [tree_lowest, tree_highest] + */ + interval_tree punch(interval_type const& ival) const + { + if (empty()) + return {}; + + interval_tree result; + auto i = std::begin(*this); + if (ival.low() < i->interval().low()) + result.insert({ival.low(), i->interval().low()}); + + for (auto e = end(); i != e; ++i) + { + auto next = i; ++next; + if (next != e) + result.insert({i->interval().high(), next->interval().low()}); + else + break; + } + + if (i != end() && i->interval().high() < ival.high()) + result.insert({i->interval().high(), ival.high()}); + + return result; + } + + iterator begin() + { + if (!root_) + return {nullptr, this}; + + auto* iter = root_; + + while (iter->left_) + iter = iter->left_; + + return{iter, this}; + } + iterator end() + { + return {nullptr, this}; + } + + const_iterator cbegin() const + { + if (!root_) + return {nullptr, this}; + + auto* iter = root_; + + while (iter->left_) + iter = iter->left_; + + return const_iterator{iter, this}; + } + const_iterator cend() const + { + return const_iterator{nullptr, this}; + } + const_iterator begin() const + { + return cbegin(); + } + const_iterator end() const + { + return cend(); + } + + /** + * Returns wether or not the tree is empty + */ + bool empty() const noexcept + { + return root_ == nullptr; + } + + private: + node_type* copyTreeImpl(node_type* root, node_type* parent) + { + if (root) + { + auto* cpy = new node_type(parent, root->interval()); + cpy->color_ = root->color_; + cpy->max_ = root->max_; + cpy->left_ = copyTreeImpl(root->left_, cpy); + cpy->right_ = copyTreeImpl(root->right_, cpy); + return cpy; + } + return nullptr; + }; + + template + bool find_all_i(node_type* ptr, interval_type const& ival, FunctionT const& on_find, ComparatorFunctionT const& compare) + { + if (compare(ptr->interval(), ival)) + { + if (!on_find(IteratorT{ptr, this})) + return false; + } + if (ptr->left_ && ival.high() <= ptr->left_->max()) + { + // no right? can only continue left + if (!ptr->right_ || ival.low() > ptr->right_->max()) + return find_all_i(ptr->left_, ival, on_find, compare); + + if (!find_all_i(ptr->left_, ival, on_find, compare)) + return false; + } + if (ptr->right_ && ival.high() <= ptr->right_->max()) + { + if (!ptr->left_ || ival.low() > ptr->left_->max()) + return find_all_i(ptr->right_, ival, on_find, compare); + + if (!find_all_i(ptr->right_, ival, on_find, compare)) + return false; + } + return true; + } + + template + node_type* find_i(node_type* ptr, interval_type const& ival, ComparatorFunctionT const& compare) + { + if (compare(ptr->interval(), ival)) + return ptr; + else + return find_i_ex(ptr, ival, compare); + } + + // excludes ptr + template + node_type* find_i_ex(node_type* ptr, interval_type const& ival, ComparatorFunctionT const& compare) + { + if (ptr->left_ && ival.high() <= ptr->left_->max()) + { + // no right? can only continue left + if (!ptr->right_ || ival.low() > ptr->right_->max()) + return find_i(ptr->left_, ival, compare); + + auto* res = find_i(ptr->left_, ival, compare); + if (res != nullptr) + return res; + } + if (ptr->right_ && ival.high() <= ptr->right_->max()) + { + if (!ptr->left_ || ival.low() > ptr->left_->max()) + return find_i(ptr->right_, ival, compare); + + auto* res = find_i(ptr->right_, ival, compare); + if (res != nullptr) + return res; + } + return nullptr; + } + + template + node_type* overlap_find_i(node_type* ptr, interval_type const& ival) + { +#if __cplusplus >= 201703L + if constexpr (Exclusive) +#else + if (Exclusive) +#endif + { + if (ptr->interval().overlaps_exclusive(ival)) + return ptr; + } + else + { + if (ptr->interval().overlaps(ival)) + return ptr; + } + + return overlap_find_i_ex(ptr, ival); + } + + template + bool overlap_find_all_i(node_type* ptr, interval_type const& ival, FunctionT const& on_find) + { +#if __cplusplus >= 201703L + if constexpr (Exclusive) +#else + if (Exclusive) +#endif + { + if (ptr->interval().overlaps_exclusive(ival)) + { + if (!on_find(IteratorT{ptr, this})) + { + return false; + } + } + } + else + { + if (ptr->interval().overlaps(ival)) + { + if (!on_find(IteratorT{ptr, this})) + { + return false; + } + } + } + if (ptr->left_ && ptr->left_->max() >= ival.low()) + { + // no right? can only continue left + // or interval low is bigger than max of right branch. + if (!ptr->right_ || ival.low() > ptr->right_->max()) + return overlap_find_all_i(ptr->left_, ival, on_find); + + if (!overlap_find_all_i(ptr->left_, ival, on_find)) + return false; + } + if (ptr->right_ && ptr->right_->max() >= ival.low()) + { + if (!ptr->left_ || ival.low() > ptr->right_->max()) + return overlap_find_all_i(ptr->right_, ival, on_find); + + if (!overlap_find_all_i(ptr->right_, ival, on_find)) + return false; + } + return true; + } + + // excludes ptr + template + node_type* overlap_find_i_ex(node_type* ptr, interval_type const& ival) + { + if (ptr->left_ && ptr->left_->max() >= ival.low()) + { + // no right? can only continue left + // or upper bounds higher than what is contained right? continue left. + if (!ptr->right_ || ival.low() > ptr->right_->max()) + return overlap_find_i(ptr->left_, ival); + + auto* res = overlap_find_i(ptr->left_, ival); + if (res != nullptr) + return res; + } + if (ptr->right_ && ptr->right_->max() >= ival.low()) + { + if (!ptr->left_ || ival.low() > ptr->left_->max()) + return overlap_find_i(ptr->right_, ival); + + auto* res = overlap_find_i(ptr->right_, ival); + if (res != nullptr) + return res; + } + return nullptr; + } + + node_type* successor(node_type* node) + { + if (node->right_) + return minimum(node->right_); + auto* y = node->parent_; + while (y && node == y->right_) + { + node = y; + y = y->parent_; + } + return y; + } + + bool is_descendant(iterator par, iterator desc) + { + auto p = desc->parent_; + for (; p && p != par.node_; p = p->parent_) {} + return p != nullptr; + } + + /** + * Set v inplace of u. Does not delete u. + * Creates orphaned nodes. A transplant call must be succeeded by delete calls. + */ + void transplant(node_type* u, node_type* v) + { + if (u->is_root()) + root_ = v; + else if (u->is_left()) + u->parent_->left_ = v; + else + u->parent_->right_ = v; + if (v) + v->parent_ = u->parent_; + } + + /** + * Get leftest of x. + */ + node_type* minimum(node_type* x) const + { + while (x->left_) + x = x->left_; + return x; + } + + void left_rotate(node_type* x) + { + auto* y = x->right_; + x->right_ = y->left_; + if (y->left_) + y->left_->parent_ = x; + + y->parent_ = x->parent_; + if (!x->parent_) + root_ = y; + else if (x->is_left()) + x->parent_->left_ = y; + else + x->parent_->right_ = y; + + y->left_ = x; + x->parent_ = y; + + // max fixup + if (x->left_ && x->right_) + x->max_ = std::max(x->interval_.high(), std::max(x->left_->max_, x->right_->max_)); + else if (x->left_) + x->max_ = std::max(x->interval_.high(), x->left_->max_); + else if (x->right_) + x->max_ = std::max(x->interval_.high(), x->right_->max_); + else + x->max_ = x->interval_.high(); + + if (y->right_) + y->max_ = std::max(y->interval_.high(), std::max(y->right_->max_, x->max_)); + else + y->max_ = std::max(y->interval_.high(), x->max_); + } + + void right_rotate(node_type* y) + { + auto* x = y->left_; + y->left_ = x->right_; + + if (x->right_) + x->right_->parent_ = y; + + x->parent_= y->parent_; + if (!y->parent_) + root_ = x; + else if (y->is_left()) + y->parent_->left_ = x; + else + y->parent_->right_ = x; + + x->right_ = y; + y->parent_ = x; + + // max fixup + if (y->left_ && y->right_) + y->max_ = std::max(y->interval_.high(), std::max(y->left_->max_, y->right_->max_)); + else if (y->left_) + y->max_ = std::max(y->interval_.high(), y->left_->max_); + else if (y->right_) + y->max_ = std::max(y->interval_.high(), y->right_->max_); + else + y->max_ = y->interval_.high(); + + if (x->left_) + x->max_ = std::max(x->interval_.high(), std::max(x->left_->max_, y->max_)); + else + x->max_ = std::max(x->interval_.high(), y->max_); + } + + void recalculate_max(node_type* reacalculation_root) + { + auto* p = reacalculation_root; + while (p && p->max_ <= reacalculation_root->max_) + { + if (p->left_ && p->left_->max_ > p->max_) + p->max_ = p->left_->max_; + if (p->right_ && p->right_->max_ > p->max_) + p->max_ = p->right_->max_; + p = p->parent_; + } + } + + void insert_fixup(node_type* z) + { + while (z->parent_ && z->parent_->color_ == rb_color::red) + { + if (!z->parent_->parent_) + break; + if (z->parent_ == z->parent_->parent_->left_) + { + node_type* y = z->parent_->parent_->right_; + if (y && y->color_ == rb_color::red) + { + z->parent_->color_ = rb_color::black; + y->color_ = rb_color::black; + z->parent_->parent_->color_ = rb_color::red; + z = z->parent_->parent_; + } + else + { + if (z == z->parent_->right_) + { + z = z->parent_; + left_rotate(z); + } + z->parent_->color_ = rb_color::black; + z->parent_->parent_->color_ = rb_color::red; + right_rotate(z->parent_->parent_); + } + } + else + { + node_type* y = z->parent_->parent_->left_; + if (y && y->color_ == rb_color::red) + { + z->parent_->color_ = rb_color::black; + y->color_ = rb_color::black; + z->parent_->parent_->color_ = rb_color::red; + z = z->parent_->parent_; + } + else + { + if (z->is_left()) + { + z = z->parent_; + right_rotate(z); + } + z->parent_->color_ = rb_color::black; + z->parent_->parent_->color_ = rb_color::red; + left_rotate(z->parent_->parent_); + } + } + } + root_->color_ = rb_color::black; + } + + void erase_fixup(node_type* x, node_type* x_parent, bool y_is_left) + { + while (x != root_ && x->color_ == rb_color::black) + { + node_type* w; + if (y_is_left) + { + w = x_parent->right_; + if (w->color_ == rb_color::red) + { + w->color_ = rb_color::black; + x_parent->color_ = rb_color::red; + left_rotate(x_parent); + w = x_parent->right_; + } + + if (w->left_->color_ == rb_color::black && w->right_->color_ == rb_color::black) + { + w->color_ = rb_color::red; + x = x_parent; + x_parent = x->parent_; + y_is_left = (x == x_parent->left_); + } + else + { + if (w->right_->color_ == rb_color::black) + { + w->left_->color_ = rb_color::black; + w->color_ = rb_color::red; + right_rotate(w); + w = x_parent->right_; + } + + w->color_ = x_parent->color_; + x_parent->color_ = rb_color::black; + if (w->right_) + w->right_->color_ = rb_color::black; + + left_rotate(x_parent); + x = root_; + x_parent = nullptr; + } + } + else + { + w = x_parent->left_; + if (w->color_ == rb_color::red) + { + w->color_ = rb_color::black; + x_parent->color_ = rb_color::red; + right_rotate(x_parent); + w = x_parent->left_; + } + + if (w->right_->color_ == rb_color::black && w->left_->color_ == rb_color::black) + { + w->color_ = rb_color::red; + x = x_parent; + x_parent = x->parent_; + y_is_left = (x == x_parent->left_); + } + else + { + if (w->left_->color_ == rb_color::black) + { + w->right_->color_ = rb_color::black; + w->color_ = rb_color::red; + left_rotate(w); + w = x_parent->left_; + } + + w->color_ = x_parent->color_; + x_parent->color_ = rb_color::black; + if (w->left_) + w->left_->color_ = rb_color::black; + + right_rotate(x_parent); + x = root_; + x_parent = nullptr; + } + } + } + + x->color_ = rb_color::black; + } + + private: + node_type* root_; + size_type size_; + }; +//############################################################################################################ + template + using interval_tree_t = interval_tree >; +//############################################################################################################ +} diff --git a/src/interval_tree_fwd.hpp b/src/interval_tree_fwd.hpp new file mode 100644 index 0000000..7547b2e --- /dev/null +++ b/src/interval_tree_fwd.hpp @@ -0,0 +1,52 @@ +// This is free and unencumbered software released into the public domain. +// +// Anyone is free to copy, modify, publish, use, compile, sell, or +// distribute this software, either in source code form or as a compiled +// binary, for any purpose, commercial or non-commercial, and by any +// means. +// +// In jurisdictions that recognize copyright laws, the author or authors +// of this software dedicate any and all copyright interest in the +// software to the public domain. We make this dedication for the benefit +// of the public at large and to the detriment of our heirs and +// successors. We intend this dedication to be an overt act of +// relinquishment in perpetuity of all present and future rights to this +// software under copyright law. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +// IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR +// OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// +// For more information, please refer to +// +// Code obtained from: +// Ebbeke, Tim (2021) "interval-tree" https://github.com/5cript/interval-tree, +// accessed August 22, 2021 commit f5bff66f88daf47b03a5a612c58f09e2a399c199. +// + +#pragma once + +namespace lib_interval_tree +{ + template + struct interval; + + template + class interval_tree; + + template + class node; + + template + class basic_interval_tree_iterator; + + template + class const_interval_tree_iterator; + + template + class interval_tree_iterator; +} diff --git a/src/interval_types.hpp b/src/interval_types.hpp new file mode 100644 index 0000000..9610c46 --- /dev/null +++ b/src/interval_types.hpp @@ -0,0 +1,71 @@ +// This is free and unencumbered software released into the public domain. +// +// Anyone is free to copy, modify, publish, use, compile, sell, or +// distribute this software, either in source code form or as a compiled +// binary, for any purpose, commercial or non-commercial, and by any +// means. +// +// In jurisdictions that recognize copyright laws, the author or authors +// of this software dedicate any and all copyright interest in the +// software to the public domain. We make this dedication for the benefit +// of the public at large and to the detriment of our heirs and +// successors. We intend this dedication to be an overt act of +// relinquishment in perpetuity of all present and future rights to this +// software under copyright law. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +// IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR +// OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// +// For more information, please refer to +// +// Code obtained from: +// Ebbeke, Tim (2021) "interval-tree" https://github.com/5cript/interval-tree, +// accessed August 22, 2021 commit f5bff66f88daf47b03a5a612c58f09e2a399c199. +// + +#pragma once + +namespace lib_interval_tree +{ + // (] + struct left_open + { + template + static inline bool within(numerical_type b, numerical_type e, numerical_type p) + { + return (b < p) && (p <= e); + } + }; + // [) + struct right_open + { + template + static inline bool within(numerical_type b, numerical_type e, numerical_type p) + { + return (b <= p) && (p < e); + } + }; + // [] + struct closed + { + template + static inline bool within(numerical_type b, numerical_type e, numerical_type p) + { + return (b <= p) && (p <= e); + } + }; + // () + struct open + { + template + static inline bool within(numerical_type b, numerical_type e, numerical_type p) + { + return (b < p) && (p < e); + } + }; +}