From b002058933bc5f97d49df5bbce7afd9077884b92 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Thu, 28 Jan 2021 19:32:01 +0000 Subject: [PATCH] EgtGeomKernel 2.3a4 : - in Zmap parallelizzata creazione da superficie TriMesh. --- EgtGeomKernel.rc | Bin 11710 -> 11710 bytes IntersLineSurfTm.cpp | 4 +- VolZmap.cpp | 105 ++++++++++++++++--- VolZmap.h | 7 ++ VolZmapCreation.cpp | 234 +++++++++++++++++++++++++++---------------- 5 files changed, 244 insertions(+), 106 deletions(-) diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index c672b2cdb5c6de453cf1bdd04dd60e43f031c535..1264ee347539f9c3f13f44917d00a8c4daf8d15b 100644 GIT binary patch delta 94 zcmdlNy)SyhFE&P#&A-_cnHfzcD{|{@_Trkr0u;H;XNwSVW8B;$>;>dw2zN+>g;Df- LFmBFL4&ed-WjGrO delta 94 zcmdlNy)SyhFE&Qw&A-_cnHh~ID{|{@_Trkr0u;H;XNwSVW8B;$>;>dw2zN+>g;Df- LFmBFL4&ed-W11TR diff --git a/IntersLineSurfTm.cpp b/IntersLineSurfTm.cpp index 59e9312..c479398 100644 --- a/IntersLineSurfTm.cpp +++ b/IntersLineSurfTm.cpp @@ -155,7 +155,7 @@ IntersParLinesSurfTm::IntersParLinesSurfTm( const Frame3d& frLines, const ISurfT //---------------------------------------------------------------------------- bool -IntersParLinesSurfTm::GetInters( const Point3d& ptL, double dLen, ILSIVECTOR& vInfo, bool bFinite) +IntersParLinesSurfTm::GetInters( const Point3d& ptL, double dLen, ILSIVECTOR& vInfo, bool bFinite) const { // verifico validità if ( ! m_bOk) @@ -189,4 +189,4 @@ IntersParLinesSurfTm::GetInters( const Point3d& ptL, double dLen, ILSIVECTOR& vI OrderInfoIntersLineSurfTm( vInfo) ; return true ; -} \ No newline at end of file +} diff --git a/VolZmap.cpp b/VolZmap.cpp index aeae9e3..7b7575a 100644 --- a/VolZmap.cpp +++ b/VolZmap.cpp @@ -1561,18 +1561,16 @@ VolZmap::SetToModifyDexelBlocks( int nGrid, int nDex, int nInt) } //---------------------------------------------------------------------------- -bool -VolZmap::IsBox( void) +bool +VolZmap::IsMapPartABox( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, double& dMinZ, double& dMaxZ) { - // se non tridexel, non posso stabilire con il metodo seguente se è un box - if ( m_nMapNum == 1) - return false ; - // analisi dei tridexel - bool bFound = true ; - for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) { - double dMinZ = m_dMaxZ[nMap] ; - double dMaxZ = m_dMinZ[nMap] ; - for ( int n = 0 ; n < int( m_nDim[nMap]) ; ++ n) { + dMinZ = m_dMaxZ[nMap] ; + dMaxZ = m_dMinZ[nMap] ; + for ( int i = nInfI ; i < nSupI ; ++ i) { + for ( int j = nInfJ ; j < nSupJ ; ++ j) { + if ( ! m_bIsBox) + return true ; + int n = j * m_nNx[nMap] + i ; int nSize = int( m_Values[nMap][n].size()) ; if ( nSize > 1) return false ; @@ -1586,12 +1584,89 @@ VolZmap::IsBox( void) abs( m_Values[nMap][n][0].dMax - dMaxZ) > EPS_SMALL) return false ; } - } + } } - if ( dMaxZ < dMinZ) - bFound = false ; } - return bFound ; + return ( dMaxZ >= dMinZ) ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::IsBox( void) +{ + // Se non tridexel, non posso stabilire con il metodo seguente se è un box + if ( m_nMapNum == 1) + return false ; + // Numero massimo di thread per il calcolo parallelo. + int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ; + // Disponibile un solo thread + if ( nThreadMax == 1) { + for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) { + double dMinZ, dMaxZ ; + if ( ! IsMapPartABox( nMap, 0, m_nNx[nMap], 0, m_nNy[nMap], dMinZ, dMaxZ)) + return false ; + } + return true ; + } + + // Caso di più thread + m_bIsBox = true ; + for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) { + vector< future> vRes ; + vRes.resize( nThreadMax) ; + std::vector vMinZ, vMaxZ ; + vMinZ.resize( nThreadMax) ; + vMaxZ.resize( nThreadMax) ; + if ( m_nNx[nMap] > m_nNy[nMap]) { + int nDexNum = m_nNx[nMap] / nThreadMax ; + int nRemainder = m_nNx[nMap] % nThreadMax ; + int nInfI = 0 ; + int nSupI = 0 ; + for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) { + nInfI = nSupI ; + nSupI = nInfI + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ; + vRes[nThread] = async( launch::async, &VolZmap::IsMapPartABox, this, nMap, + nInfI, nSupI, 0, m_nNy[nMap], ref( vMinZ[nThread]), ref( vMaxZ[nThread])) ; + } + } + else { + int nDexNum = m_nNy[nMap] / nThreadMax ; + int nRemainder = m_nNy[nMap] % nThreadMax ; + int nInfJ = 0 ; + int nSupJ = 0 ; + for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) { + nInfJ = nSupJ ; + nSupJ = nInfJ + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ; + vRes[nThread] = async( launch::async, &VolZmap::IsMapPartABox, this, nMap, + 0, m_nNx[nMap], nInfJ, nSupJ, ref(vMinZ[nThread]), ref(vMaxZ[nThread])) ; + } + } + // Ciclo per attendere che tutti gli async abbiano terminato. + int nTerminated = 0 ; + while ( nTerminated < nThreadMax) { + for ( int nL = 0 ; nL < nThreadMax ; ++ nL) { + // Async terminato + if ( vRes[nL].valid() && vRes[nL].wait_for(chrono::microseconds{ 1 }) == future_status::ready) { + ++ nTerminated ; + if ( ! vRes[nL].get()) { + m_bIsBox = false ; + } + } + } + } + // Se uno dei thread trova che la sua porzione non è un box, non lo può essere il solido intero. + if ( ! m_bIsBox) + return false ; + // Controllo che gli estremi Z siano uguali. + for ( int nT = 1 ; nT < nThreadMax ; ++ nT) { + if ( abs( vMinZ[nT] - vMinZ[0]) > EPS_SMALL) + return false ; + if ( abs( vMaxZ[nT] - vMaxZ[0]) > EPS_SMALL) + return false ; + } + } + + return m_bIsBox ; } //---------------------------------------------------------------------------- diff --git a/VolZmap.h b/VolZmap.h index 8e16b51..f2218b1 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -17,6 +17,7 @@ #include "GeoObjRW.h" #include "Tool.h" #include "/EgtDev/Include/EGkVolZmap.h" +#include "/EgtDev/Include/EGkIntersLineSurfTm.h" #include #include #include @@ -374,6 +375,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool ExpandFromZInterval( IntContainer& IntCont) ; bool FirstExpansionFromZ( int nNumThread, IntervalIndexes IntSt, IntContainerVec& IntervalsToProcessStackVec) ; bool ProcessIntervals( IntContainer& IntervalsToProcess) ; + bool IsMapPartABox( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, double& dMinZ, double& dMaxZ) ; bool IsBox( void) ; // Avoid semplici bool AvoidSimpleBox( const Frame3d& frBox, const Vector3d& vtDiag, bool bPrecise = false) const ; @@ -386,6 +388,10 @@ class VolZmap : public IVolZmap, public IGeoObjRW // Funzioni ausiliarie per metodi avoid bool SingleMapDexelConeCollision( int nStI, int nEnI, int nStJ, int nEnJ, const Point3d& ptRefPoint, const Vector3d& vtRefAx, double dMinRad, double dMaxRad, double dHeight, double dMinBoxH, double dMaxBoxH) const ; + // Funzione per crezione solido in parallelo + bool CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig, + const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM) ; + private : enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ; enum Shape { GENERIC = 0, BOX = 1, EXTRUSION = 2} ; @@ -443,6 +449,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW mutable std::mutex SliceMutex ; mutable std::atomic m_bBreak ; + std::atomic m_bIsBox ; Tool m_Tool ; } ; diff --git a/VolZmapCreation.cpp b/VolZmapCreation.cpp index 0b80a21..583df4c 100644 --- a/VolZmapCreation.cpp +++ b/VolZmapCreation.cpp @@ -17,8 +17,8 @@ #include "CurveLine.h" #include "VolZmap.h" #include "GeoConst.h" -#include "/EgtDev/Include/EGkIntersLineSurfTm.h" #include "/EgtDev/Include/EgtNumUtils.h" +#include using namespace std ; @@ -417,6 +417,111 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double return true ; } +//---------------------------------------------------------------------------- +bool +VolZmap::CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig, + const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM) +{ + if ( nMap < 0 || nMap > 2 || + nInfI < 0 || nInfI > m_nNx[nMap] || + nSupI < 0 || nSupI > m_nNx[nMap] || + nInfJ < 0 || nInfJ > m_nNy[nMap] || + nSupJ < 0 || nSupJ > m_nNy[nMap]) + return false ; + + // Determinazione e ridimensionamento dei dexel interni alla trimesh + for ( int i = nInfI ; i < nSupI ; ++ i) { + for ( int j = nInfJ ; j < nSupJ ; ++ j) { + + // Definisco la retta da intersecare con la trimesh + double dX = ( i + 0.5) * m_dStep ; + double dY = ( j + 0.5) * m_dStep ; + Point3d ptP0( dX, dY, 0) ; + + // Determino le intersezioni della retta con la TriMesh + ILSIVECTOR IntersectionResults ; + intPLSTM.GetInters( ptP0, vtLen.v[(nMap+2)%3], IntersectionResults) ; + + for ( int nI = 0 ; nI < int( IntersectionResults.size()) - 3 ; ++ nI) { + int nJ = nI + 1 ; + int nK = nJ + 1 ; + int nT = nK + 1 ; + int nSgnI = IntersectionResults[nI].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nI].dCosDN > -EPS_SMALL ? 0 : - 1 ; + int nSgnJ = IntersectionResults[nJ].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nJ].dCosDN > -EPS_SMALL ? 0 : - 1 ; + int nSgnK = IntersectionResults[nK].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nK].dCosDN > -EPS_SMALL ? 0 : - 1 ; + int nSgnT = IntersectionResults[nT].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nT].dCosDN > -EPS_SMALL ? 0 : - 1 ; + double dUJ = IntersectionResults[nJ].dU ; + double dUK = IntersectionResults[nK].dU ; + if ( nSgnI != 0 && nSgnI == nSgnJ && nSgnK != 0 && nSgnK == nSgnT && nSgnI == - nSgnT && abs( dUJ - dUK) < EPS_SMALL) { + IntersectionResults.erase( IntersectionResults.begin() + nK) ; + IntersectionResults.erase( IntersectionResults.begin() + nJ) ; + } + } + + int nInt = int( IntersectionResults.size()) ; + + int nPos = j * m_nNx[nMap] + i ; + + bool bInside = false ; + Point3d ptIn ; + Vector3d vtInN ; + + for ( int k = 0 ; k < nInt ; ++ k) { + + int nIntType = IntersectionResults[k].nILTT ; + + // Se c'è intersezione + if ( nIntType != ILTT_NO) { + + double dCos = IntersectionResults[k].dCosDN ; + + // entro nella superficie trimesh + if ( dCos < - EPS_SMALL) { + + ptIn = IntersectionResults[k].ptI ; + + int nT = IntersectionResults[k].nT ; + int nF = Surf.GetFacetFromTria( nT) ; + + Surf.GetFacetNormal( nF, vtInN) ; + + bInside = true ; + } + + // esco dalla superficie trimesh + else if ( dCos > EPS_SMALL && bInside) { + + Point3d ptOut = IntersectionResults[k].ptI ; + + int nT = IntersectionResults[k].nT ; + int nF = Surf.GetFacetFromTria( nT) ; + + Vector3d vtOutN ; + Surf.GetFacetNormal( nF, vtOutN) ; + + int nCurrentSize = int( m_Values[nMap][nPos].size()) ; + + // Aggiungo un tratto al dexel + m_Values[nMap][nPos].resize( nCurrentSize + 1) ; + + // Aggiorno dati del tratto di dexel + m_Values[nMap][nPos][nCurrentSize].dMin = ptIn.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3] ; + m_Values[nMap][nPos][nCurrentSize].dMax = ptOut.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3] ; + m_Values[nMap][nPos][nCurrentSize].vtMinN = vtInN ; + m_Values[nMap][nPos][nCurrentSize].vtMaxN = vtOutN ; + m_Values[nMap][nPos][nCurrentSize].nToolMin = 0 ; + m_Values[nMap][nPos][nCurrentSize].nToolMax = 0 ; + m_Values[nMap][nPos][nCurrentSize].nCompo = 0 ; + + bInside = false ; + } + } + } + } + } + return true ; +} + //---------------------------------------------------------------------------- bool VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex) @@ -487,6 +592,7 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex return false ; // ciclo sulle griglie + bool bCompleted = true ; for ( int g = 0 ; g < m_nMapNum ; ++ g) { // Definisco dei sistemi di riferimento ausiliari @@ -501,93 +607,43 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex // Oggetto per calcolo massivo intersezioni IntersParLinesSurfTm intPLSTM( frMapFrame, Surf) ; - // Determinazione e ridimensionamento dei dexel interni alla trimesh - for ( int i = 0 ; i < m_nNx[g] ; ++ i) { - for ( int j = 0 ; j < m_nNy[g] ; ++ j) { - - // Definisco la retta da intersecare con la trimesh - double dX = ( i + 0.5) * m_dStep ; - double dY = ( j + 0.5) * m_dStep ; - Point3d ptP0( dX, dY, 0) ; - - // Determino le intersezioni della retta con la TriMesh - ILSIVECTOR IntersectionResults ; - intPLSTM.GetInters( ptP0, vtLen.v[(g+2)%3], IntersectionResults) ; - - for ( int nI = 0 ; nI < int( IntersectionResults.size()) - 3 ; ++ nI) { - int nJ = nI + 1 ; - int nK = nJ + 1 ; - int nT = nK + 1 ; - int nSgnI = IntersectionResults[nI].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nI].dCosDN > -EPS_SMALL ? 0 : - 1 ; - int nSgnJ = IntersectionResults[nJ].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nJ].dCosDN > -EPS_SMALL ? 0 : - 1 ; - int nSgnK = IntersectionResults[nK].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nK].dCosDN > -EPS_SMALL ? 0 : - 1 ; - int nSgnT = IntersectionResults[nT].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nT].dCosDN > -EPS_SMALL ? 0 : - 1 ; - double dUJ = IntersectionResults[nJ].dU ; - double dUK = IntersectionResults[nK].dU ; - if ( nSgnI != 0 && nSgnI == nSgnJ && nSgnK != 0 && nSgnK == nSgnT && nSgnI == - nSgnT && abs( dUJ - dUK) < EPS_SMALL) { - IntersectionResults.erase( IntersectionResults.begin() + nK) ; - IntersectionResults.erase( IntersectionResults.begin() + nJ) ; - } - } - - int nInt = int( IntersectionResults.size()) ; - - int nPos = j * m_nNx[g] + i ; - - bool bInside = false ; - Point3d ptIn ; - Vector3d vtInN ; - - for ( int k = 0 ; k < nInt ; ++ k) { - - int nIntType = IntersectionResults[k].nILTT ; - - // Se c'è intersezione - if ( nIntType != ILTT_NO) { - - double dCos = IntersectionResults[k].dCosDN ; - - // entro nella superficie trimesh - if ( dCos < - EPS_SMALL) { - - ptIn = IntersectionResults[k].ptI ; - - int nT = IntersectionResults[k].nT ; - int nF = Surf.GetFacetFromTria( nT) ; - - Surf.GetFacetNormal( nF, vtInN) ; - - bInside = true ; - } - - // esco dalla superficie trimesh - else if ( dCos > EPS_SMALL && bInside) { - - Point3d ptOut = IntersectionResults[k].ptI ; - - int nT = IntersectionResults[k].nT ; - int nF = Surf.GetFacetFromTria( nT) ; - - Vector3d vtOutN ; - Surf.GetFacetNormal( nF, vtOutN) ; - - int nCurrentSize = int( m_Values[g][nPos].size()) ; - - // Aggiungo un tratto al dexel - m_Values[g][nPos].resize( nCurrentSize + 1) ; - - // Aggiorno dati del tratto di dexel - m_Values[g][nPos][nCurrentSize].dMin = ptIn.v[(g+2)%3] - ptMapOrig.v[(g+2)%3] ; - m_Values[g][nPos][nCurrentSize].dMax = ptOut.v[(g+2)%3] - ptMapOrig.v[(g+2)%3] ; - m_Values[g][nPos][nCurrentSize].vtMinN = vtInN ; - m_Values[g][nPos][nCurrentSize].vtMaxN = vtOutN ; - m_Values[g][nPos][nCurrentSize].nToolMin = 0 ; - m_Values[g][nPos][nCurrentSize].nToolMax = 0 ; - m_Values[g][nPos][nCurrentSize].nCompo = 0 ; - - bInside = false ; - } - } + // Numero massimo di thread + int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ; + vector< future> vRes ; + vRes.resize( nThreadMax) ; + if ( m_nNx[g] > m_nNy[g]) { + int nDexNum = m_nNx[g] / nThreadMax ; + int nRemainder = m_nNx[g] % nThreadMax ; + int nInfI = 0 ; + int nSupI = 0 ; + for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) { + nInfI = nSupI ; + nSupI = nInfI + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ; + vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, g, + nInfI, nSupI, 0, m_nNy[g], ref( vtLen), ref( ptMapOrig), ref( Surf), ref( intPLSTM)) ; + } + } + else { + int nDexNum = m_nNy[g] / nThreadMax ; + int nRemainder = m_nNy[g] % nThreadMax ; + int nInfJ = 0 ; + int nSupJ = 0 ; + for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) { + nInfJ = nSupJ ; + nSupJ = nInfJ + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ; + vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, g, + 0, m_nNx[g], nInfJ, nSupJ, ref( vtLen), ref( ptMapOrig), ref( Surf),ref( intPLSTM)) ; + } + } + + // Ciclo per attendere che tutti gli async abbiano terminato. + int nTerminated = 0 ; + while ( nTerminated < nThreadMax) { + for ( int nL = 0 ; nL < nThreadMax ; ++ nL) { + // Async terminato + if ( vRes[nL].valid() && vRes[nL].wait_for( chrono::microseconds{ 1}) == future_status::ready) { + ++ nTerminated ; + bCompleted = bCompleted && vRes[nL].get() ; } } } @@ -607,5 +663,5 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex // Aggiornamento dello stato m_nStatus = OK ; - return true ; + return bCompleted ; }