From d67308f4f816c5f20efa53d93eab0bc7af61e068 Mon Sep 17 00:00:00 2001 From: Esukhareva Date: Tue, 19 Nov 2024 18:10:56 +0000 Subject: [PATCH 1/2] bos #43483 [CEA 42969] Viscous Layers generation issue --- src/StdMeshers/StdMeshers_ViscousLayers2D.cxx | 40 ++++++++++--------- 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/src/StdMeshers/StdMeshers_ViscousLayers2D.cxx b/src/StdMeshers/StdMeshers_ViscousLayers2D.cxx index c25d0d05b..e6797ca15 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers2D.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers2D.cxx @@ -1867,26 +1867,30 @@ bool _ViscousBuilder2D::shrink() Geom2dAdaptor_Curve edgeCurve( pcurve, Min( uf, ul ), Max( uf, ul )); Geom2dAdaptor_Curve seg2Curve( seg2Line ); Geom2dInt_GInter curveInt( edgeCurve, seg2Curve, 1e-7, 1e-7 ); + + // In the older version length2D was set to this value only inside the !convex if block + // But it seems that length2D can be set here anyway, because if not set valid value of length2D here, + // it will be calculated later using length1D, and it can be not valid in cases if length1D is too large or too small. + length2D = L2->_lEdges[iFSeg2]._length2D; + + /* convex VERTEX + * L seg2 + * | o---o--- + * | / | + * |/ | L2 + * x------x--- */ + /* concave VERTEX + * o-----o--- + * \ | + * \ | L2 + * x--x--- + * / + * L / */ isConvex = ( curveInt.IsDone() && !curveInt.IsEmpty() ); - if ( isConvex ) { - /* convex VERTEX */ + if ( isConvex ) + { length1D = Abs( u - curveInt.Point( 1 ).ParamOnFirst() ); - double maxDist2d = 2 * L2->_lEdges[ iLSeg2 ]._length2D; - isConvex = ( length1D < maxDist2d * len1dTo2dRatio ); - /* |L seg2 - * | o---o--- - * | / | - * |/ | L2 - * x------x--- */ - } - if ( !isConvex ) { /* concave VERTEX */ /* o-----o--- - * \ | - * \ | L2 - * x--x--- - * / - * L / */ - length2D = L2->_lEdges[ iFSeg2 ]._length2D; - //if ( L2->_advancable ) continue; + length2D = Max(length2D, length1D / len1dTo2dRatio); } } else // L2 is advancable but in the face adjacent by L From 63567f57a1452106c71ac9eb9452d6f88118b460 Mon Sep 17 00:00:00 2001 From: mbs Date: Wed, 13 Nov 2024 15:50:59 +0000 Subject: [PATCH 2/2] [bos #43494] fixed Arithmetic1D discretization and added new error compensation function --- doc/gui/input/1d_meshing_hypo.rst | 7 +- src/StdMeshers/StdMeshers_Regular_1D.cxx | 120 +++++++++++++++++++++-- 2 files changed, 116 insertions(+), 11 deletions(-) diff --git a/doc/gui/input/1d_meshing_hypo.rst b/doc/gui/input/1d_meshing_hypo.rst index eba4a83b6..21c2c50c5 100644 --- a/doc/gui/input/1d_meshing_hypo.rst +++ b/doc/gui/input/1d_meshing_hypo.rst @@ -1,4 +1,5 @@ .. _a1d_meshing_hypo_page: +.. |larr| unicode:: U+02190 .. LEFTWARDS ARROW ********************* 1D Meshing Hypotheses @@ -57,7 +58,7 @@ Adaptive hypothesis Arithmetic Progression hypothesis ################################# -**Arithmetic Progression** hypothesis allows to split edges into segments with a length that changes in arithmetic progression (Lk = Lk-1 + d) beginning from a given starting length and up to a given end length. +**Arithmetic Progression** hypothesis allows to split edges into segments with a length that changes in arithmetic progression (L\ :sub:`k`\ |larr| L\ :sub:`k-1`\ + d) beginning from a given starting length and up to a given end length. The splitting direction is defined by the orientation of the underlying geometrical edge. **Reverse Edges** list box allows specifying the edges, for which the splitting should be made in the direction opposite to their orientation. This list box is usable only if a geometry object is selected for meshing. In this case it is possible to select edges to be reversed either directly picking them in the 3D viewer or by selecting the edges or groups of edges in the Object Browser. Use **Add** button to add the selected edges to the list. @@ -82,7 +83,7 @@ The splitting direction is defined by the orientation of the underlying geometri Geometric Progression hypothesis ################################ -**Geometric Progression** hypothesis allows splitting edges into segments with a length that changes in geometric progression (Lk = Lk-1 * d) starting from a given **Start Length** and with a given **Common Ratio**. +**Geometric Progression** hypothesis allows splitting edges into segments with a length that changes in geometric progression (L\ :sub:`k`\ |larr| L\ :sub:`k-1`\ * d) starting from a given **Start Length** and with a given **Common Ratio**. The splitting direction is defined by the orientation of the underlying geometrical edge. **Reverse Edges** list box allows specifying the edges, for which the splitting should be made in the direction opposite to their orientation. This list box is usable only if a geometry object is selected for meshing. In this case it is possible to select edges to be reversed either directly picking them in the 3D viewer or by selecting the edges or groups of edges in the Object Browser. Use **Add** button to add the selected edges to the list. @@ -177,7 +178,7 @@ You can set the type of node distribution for this hypothesis in the **Hypothesi **Scale Distribution** - length of segments gradually changes depending on the **Scale Factor**, which is a ratio of the first segment length to the last segment length. -Length of segments changes in geometric progression with the common ratio (A) depending on the **Scale Factor** (S) and **Number of Segments** (N) as follows: A = S**(1/(N-1)). For an edge of length L, length of the first segment is L * (1 - A)/(1 - A**N) +Length of segments changes in geometric progression with the common ratio (A) depending on the **Scale Factor** (S) and **Number of Segments** (N) as follows: A = S\ :sup:`(1/(N-1))`\ . For an edge of length L, length of the first segment is L * (1 - A)/(1 - A\ :sup:`N`\ ) .. image:: ../images/a-nbsegments2.png :align: center diff --git a/src/StdMeshers/StdMeshers_Regular_1D.cxx b/src/StdMeshers/StdMeshers_Regular_1D.cxx index 13b63bcc6..bd0b99e0b 100644 --- a/src/StdMeshers/StdMeshers_Regular_1D.cxx +++ b/src/StdMeshers/StdMeshers_Regular_1D.cxx @@ -492,6 +492,105 @@ static void compensateError(double a1, double an, } } + +//================================================================================ +/*! + * \brief adjust internal node parameters so that the last segment length == an, + * and by distributing the error for the total length of curve segments + * in relation to the target length computed from the current parameters + * \param a1 - the first segment length + * \param an - the last segment length + * \param U1 - the first edge parameter + * \param Un - the last edge parameter + * \param length - the edge length + * \param C3d - the edge curve + * \param theParams - internal node parameters to adjust + */ +//================================================================================ + +static void distributeError(double a1, double an, + double U1, double Un, + double length, + Adaptor3d_Curve& C3d, + list & theParams) +{ + // Compute the error of the total length based in the current curve parameters + double tol = Min( Precision::Confusion(), 0.01 * Min(a1, an) ); + double totalLength = 0.0; + double prevParam = U1; + list segLengths; + list::iterator itU = theParams.begin(); + for ( ; itU != theParams.end(); ++itU ) + { + // Compute the curve length between two adjacent parameters and sum them up + double curLength = GCPnts_AbscissaPoint::Length(C3d, prevParam, *itU, tol); + segLengths.push_back(curLength); + totalLength += curLength; + prevParam = *itU; + } + // Calculate the error between the total length of all segments based on given parameters + // and the target length of the edge itself + double error = totalLength - length; + // Compute the sum of all internal segments (= total computed length minus the length of + // the start and end segments) + double midLength = totalLength - (a1 + an); + + // We only need to distribute the error, if the current parametrization is not correct, + // and if there are multiple internal segments + smIdType nPar = theParams.size(); + if ( a1 + an <= length && nPar > 1 && fabs(error) > tol ) + { + // Update the length of each internal segment (start and end length are given and not changed) + double newTotalLength = 0.0; + double newLength; + double relError = error / midLength; + list newSegLengths; + list::iterator itL = segLengths.begin(); + for ( ; itL != segLengths.end(); ++itL ) + { + // Do not update, but copy the first and the last segment lengths + newLength = *itL; + if (itL != segLengths.begin() && itL != --segLengths.end()) + { + newLength -= newLength * relError; + } + newSegLengths.push_back(newLength); + newTotalLength += newLength; + } + bool reverse = ( U1 > Un ); + + // Update the parameters of the curve based on the new lengths + double curveLength, tol2, U; + double prevU = U1; + itU = theParams.begin(); + itL = newSegLengths.begin(); + for ( ; itU != theParams.end(); ++itU, ++itL ) + { + curveLength = (reverse ? -(*itL) : *itL); + tol2 = Min( Precision::Confusion(), fabs(curveLength) / 100. ); + GCPnts_AbscissaPoint Discret( tol2, C3d, curveLength, prevU ); + if ( !Discret.IsDone() ) + { + return; + } + U = Discret.Parameter(); + + double sign = reverse ? -1 : 1; + if ( sign*U1 < sign*U && sign*U < sign*Un ) + { + *itU = U; + } + else + { + *itU = (sign*U >= sign*Un ? Un : U1); + break; + } + prevU = U; + } + } +} + + //================================================================================ /*! * \brief Class used to clean mesh on edges when 0D hyp modified. @@ -1042,8 +1141,9 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh & theMesh, return error ( SMESH_Comment("Invalid segment lengths (")< numeric_limits::min() ? ( 1+( an-a1 )/q ) : ( 1+theLength/a1 )); + // Compute first the number of segments and then the arithmetic increment based on that number + int n = static_cast(2 * theLength / ( a1 + an ) + 0.5); + double q = (n > 1 ? ( an - a1 ) / (n - 1) : 0.0); double U1 = theReverse ? l : f; double Un = theReverse ? f : l; @@ -1054,19 +1154,23 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh & theMesh, eltSize = -eltSize; q = -q; } - while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) { + for (int i=0; i at the distance // from the point of parameter . GCPnts_AbscissaPoint Discret( tol, theC3d, eltSize, param ); if ( !Discret.IsDone() ) break; param = Discret.Parameter(); - if ( param > f && param < l ) - theParams.push_back( param ); - else - break; + theParams.push_back( param ); eltSize += q; } - compensateError( a1, an, U1, Un, theLength, theC3d, theParams ); + + distributeError( a1, an, U1, Un, theLength, theC3d, theParams ); + + // Do not include the parameter for the start or end of an edge in the list of parameters + // NOTE: it is required to correctly distribute the error + if (fabs(theParams.front() - U1) < tol) theParams.pop_front(); + if (fabs(theParams.back() - Un) < tol) theParams.pop_back(); + if ( theReverse ) theParams.reverse(); // NPAL18025 return true;