30318764346

Committed 28 Jul 2026 12:55AM UTC coverage: 76.972% (+1.1%) from 75.909%

Build # 30318764346

Build Type

Pull #166

github

Committed by

web-flow

Commit Message

Merge 3346edfd3 into 955f92dca

Pull Request Pull Request #166: Add L2/H1 boundary optimization

Coverage Stats

1631 of 1964 new or added lines in 40 files covered. (83.04%)

48 existing lines in 4 files now uncovered.

9563 of 12424 relevant lines covered (76.97%)

634449.32 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.7

/Source/BoundaryOptimization/CurveOptimization.f90

! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! 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 OR COPYRIGHT HOLDERS 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.
!
! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
!
! --- End License
!
!////////////////////////////////////////////////////////////////////////
!
!      CurveOptimization.f90
!      Created: November 20, 2025 at 11:10 AM
!      By: David Kopriva
!
!  Compute an optimized approximation to a curve.
!
!  There are three main entry points:
!
!
!  (1) SUBROUTINE OptimizeCurve(curve, polyOrder, cuts, options, newName, newID, optimized)
!
!      Use this if the location of the segment boundaries given in cuts is known.
!      Returns an SMCurve subclass "optimized".
!
!  (2) SUBROUTINE OptimizeCurveWithSubdivisions(curve, polyOrder, options, newName, newID, optimized)
!
!      Use this to automatically compute the segments so that the error is is within the desired
!      tolerance, toler. Returns an SMCurve subclass "optimized". Uses subdivision algorithm
!      *** THIS ROUTINE IS INEFFICIENT COMPARED TO THE NEXT. It IS INCLUDED BUT COMMENTED OUT FOR SOME FUTURE NEED ***
!
!  (3) SUBROUTINE OptimizeCurveByMarching(curve, polyOrder, options, newName, newID, optimized)
!
!      Use this to automatically compute the segments so that the error is is within the desired
!      tolerance, toler. Returns an SMCurve subclass "optimized". Uses the marching algorithm
!
!
!////////////////////////////////////////////////////////////////////////
!
   Module CurveOptimization
   USE SMCurveClass
   USE ConstrainedMultiH1Optimization
   USE MultiSegmentModalCurveClass
   IMPLICIT NONE

   INTEGER, PARAMETER, PRIVATE :: MAX_DEPTH = 10
   INTEGER, PARAMETER, PRIVATE :: LEFT_SIDE = 0, RIGHT_SIDE = 1
   INTEGER, PARAMETER          :: USER_NORM = 1, MAX_NORM   = 2, MAX_NORM_DERIV = 3
!
!  ========
   CONTAINS
!  ========
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE OptimizeCurve(curve, polyOrder, cuts, breakIndices, &
                            options, newName, newID, optimized)
!
!  ----------------------------------------------------
!  Given an SMCurve, return and optimal multi-segment
!  curve that approximates it
!  ----------------------------------------------------
!
      IMPLICIT NONE
!
!     ---------
!     Arguments
!     ---------
!
      CLASS(SMCurve), POINTER :: curve              ! The curve to be approximated
      CLASS(SMCurve), POINTER :: optimized          ! The new approximation to be returned
      INTEGER                 :: polyOrder          ! polynomial order of the approximation
      REAL(KIND=RP)           :: cuts(0:)           ! Segment boundaries of the new approximation
      CHARACTER(LEN=*)        :: newName            ! Name of the new curve
      INTEGER                 :: newID              ! ID of the new curve
      TYPE(OptimizerOptions)  :: options            ! parameters for the approximation
      INTEGER, ALLOCATABLE    :: breakIndices(:)    ! index in optimizedCuts where break boundaries occur
!
!     ---------------
!     Local variables
!     ---------------
!
      TYPE(MultiH1Optimizer)                   :: multiOptimizer
      TYPE(MultiSegmentModalCurve)   , POINTER :: MSMCurve

      INTEGER                                  :: quadratureOrder
      INTEGER                                  :: nSegments
      REAL(KIND=RP), ALLOCATABLE               :: optimalCoefficients(:,:,:)
!
!     ---------------
!     Optimizer setup
!     ---------------
!
      nSegments       = SIZE(cuts) - 1
      quadratureOrder = 2*polyOrder

      CALL multiOptimizer % construct(                                 &
                                      N            = polyOrder,        &
                                      M            = quadratureOrder,  &
                                      nSegments    = nSegments,        &
                                      cuts         = cuts,             &
                                      breakIndices = breakIndices,     &
                                      options      = options)
!
!     ------------------
!     Find optimal curve
!     ------------------
!
      ALLOCATE(optimalCoefficients(0:polyOrder,2,nSegments), source = 0.0_RP)
      CALL multiOptimizer % Optimize(curve, optimalCoefficients)
!
!     --------------------------------------------------------------
!     Create a multiSegmentModalCurve from the computed coefficients
!     --------------------------------------------------------------
!
      ALLOCATE(MSMCurve)
      CALL MSMCurve % constructMultiSegmentModalCurve(curve, cuts, optimalCoefficients, newName, newID)

      optimized => MSMCurve

   END SUBROUTINE OptimizeCurve
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE OptimizeCurveByMarching(curve, polyOrder, breaks, breakIndices, &
                                      options, newName, newID, optimized)
!
!  ----------------------------------------------------------------------------------
!  March along the curve and use bisection and secant to find the length of each
!  segment so that the H^1 error is less than the tolerance allowed for that segment.
!  Clean up at the end by applying the global optimization
!  ----------------------------------------------------------------------------------
!
      IMPLICIT NONE
!
!     ---------
!     Arguments
!     ---------
!
      CLASS(SMCurve), POINTER :: curve              ! The curve to be approximated
      CLASS(SMCurve), POINTER :: optimized          ! The new approximation
      INTEGER                 :: polyOrder          ! polynomial order of the new approximation
      REAL(KIND=RP)           :: breaks(0:)         ! Location of forced cut points, starting with t = 0.0
                                                    ! If none, use breaks = [0.0_RP,1.0_RP]
      CHARACTER(LEN=*)        :: newName            ! Name of the new curve
      INTEGER                 :: newID              ! ID of the new curve
      TYPE(OptimizerOptions)  :: options            ! parameters for the approximation
      INTEGER, ALLOCATABLE    :: breakIndices(:)    ! index in optimizedCuts where break boundaries occur
!
!     ---------------
!     Local variables
!     ---------------
!
      REAL(KIND=RP), ALLOCATABLE             :: optimizedCuts(:)
      CLASS(SMCurve)               , POINTER :: optimizedCurve
      CLASS(MultiSegmentModalCurve), POINTER :: msmCurve
      TYPE(GaussQuadratureType)              :: gQuad
      TYPE(OptimizerOptions)                 :: savedOptions   ! parameters for the marching algorithm
      REAL(KIND=RP), ALLOCATABLE             :: errors(:,:)
      REAL(KIND=RP)                          :: e, errorRatio
      INTEGER                                :: n


      CALL ConstructGaussQuadrature(gQuad, 4*polyOrder) ! For error computation. The 4 is arbitrary
!
!     -----------------------------------------------------------------------
!     Save the options because for marching only end constraints will be used
!     -----------------------------------------------------------------------
!
      savedOptions            = options
      options % endConstraint = FIXED_CONSTRAINT

      DO n = 1,3 ! This is a try/try again to account for constraints not being included in marching.
!
!        ---------------------------
!        Find the optimized segments
!        ---------------------------
!
         CALL FindOptimizedCuts(curve, polyOrder, breaks, options, gQuad, &
                                optimizedCuts, breakIndices)
!
!        ------------------------------------------------------------
!        Construct the global polynomial using the optimized segments
!        ------------------------------------------------------------
!
         CALL OptimizeCurve(curve              = curve,              &
                            polyOrder          = polyOrder,          &
                            cuts               = optimizedCuts,      &
                            breakIndices       = breakIndices,       &
                            options            = savedOptions,       &
                            newName            = newName,            &
                            newID              = newID,              &
                            optimized          = optimizedCurve)
!
!        --------------------------------------------------------------------------
!        See if the error is within tolerance. Since the marching can only use the
!        FIXED_CONSTRAINT option, and since errors increase as the constraint
!        order increases, then the error won't be within the tolerance. If not,
!        create a new tolerance that is small enough so that when done we should
!        be OK
!        --------------------------------------------------------------------------
!
         CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
         CALL SegmentErrors(exact          = curve,    &
                            segmentedCurve = msmCurve, &
                            gQuad          = gQuad,    &
                            errors         = errors)
         e = MAXVAL(errors(:,USER_NORM))
         IF ( e .le. options % toler )     EXIT

         errorRatio = options % toler/e
         options % toler = options % toler*errorRatio
!         IF(errorRatio < 1.0_RP) WRITE(0,*) "Error required adjusting", e, errorRatio, TRIM(curve % curveName()) !DEBUG
!
         IF ( n < 3 )     THEN
            DEALLOCATE(errors)
            DEALLOCATE(optimizedCuts)
            CALL releaseMultiSegmentModalCurve(msmCurve)
         END IF

      END DO

      optimized => optimizedCurve
!
   END SUBROUTINE OptimizeCurveByMarching
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE FindOptimizedCuts(curve, polyOrder, breaks, options, gQuad, &
                                optimizedCuts, breakIndices)
!
!  ----------------------------------------------------------------------------------
!  March along the curve to find the length of each
!  segment so that the H^1 error is less than the tolerance allowed for that segment.
!  ----------------------------------------------------------------------------------
!
      IMPLICIT NONE
!
!     ---------
!     Arguments
!     ---------
!
      CLASS(SMCurve), POINTER   , INTENT(IN)  :: curve              ! The curve to be approximated
      INTEGER                   , INTENT(IN)  :: polyOrder          ! polynomial order of the new approximation
      REAL(KIND=RP)             , INTENT(IN)  :: breaks(0:)         ! t Location of forced break points, starting with t = 0.0
!                                                                     If none, use [0.0_RP,1.0_RP]
      TYPE(OptimizerOptions)   ,  INTENT(IN)  :: options            ! parameters for the approximation
      TYPE(GaussQuadratureType),  INTENT(IN)  :: gQuad
      REAL(KIND=RP), ALLOCATABLE, INTENT(OUT) :: optimizedCuts(:)   ! t location of segment boundaries starting with t = 0
      INTEGER      , ALLOCATABLE, INTENT(OUT) :: breakIndices(:)    ! index in optimizedCuts where break boundaries occur
!
!     ---------------
!     Local variables
!     ---------------
!
      REAL(KIND=RP), ALLOCATABLE :: rangeCuts(:)
      INTEGER                    :: nBreaks
      INTEGER                    :: j
      REAL(KIND=RP)              :: limit = 1.0d-14
      REAL(KIND=RP), ALLOCATABLE :: tmpA(:)

      nBreaks = SIZE(breaks) - 1
!
!     -------------------------------------------------------------
!     There will be at least one, and possibly more forced segments
!     -------------------------------------------------------------
!
      CALL FindOptimizedCutsInRange(curve, polyOrder, [breaks(0), breaks(1)-Limit], &
                                    options, gQuad, tmpA)
      breakIndices = [SIZE(tmpA)]
      tmpA(UBOUND(tmpA)) = breaks(1)
      IF(nBreaks == 1) THEN
         ALLOCATE(optimizedCuts(0:SIZE(tmpA)-1))
         optimizedCuts(0:) = tmpA
         RETURN
      END IF
!
!     ----------------------------
!     Otherwise, continue marching
!     ----------------------------
!
      DO j = 2, nBreaks
         CALL FindOptimizedCutsInRange(curve, polyOrder, [breaks(j-1)+limit, breaks(j)-limit], &
                                       options, gQuad, rangeCuts)
         breakIndices  = [breakIndices,breakIndices(j-1) + SIZE(rangeCuts) - 1]
         rangeCuts(UBOUND(rangeCuts)) = breaks(j)
         tmpA = [tmpA, rangeCuts(1:)]
         DEALLOCATE(rangeCuts)
      END DO
      ALLOCATE(optimizedCuts(0:SIZE(tmpA)-1))
      optimizedCuts(0:) = tmpA

   END SUBROUTINE FindOptimizedCuts
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE FindOptimizedCutsInRange(curve, polyOrder, tRange, options, gQuad, optimizedCuts)
!
!  ----------------------------------------------------------------------------------
!  March along the curve to find the length of each
!  segment so that the H^1 error is less than the tolerance allowed for that segment.
!  ----------------------------------------------------------------------------------
!
      IMPLICIT NONE
!
!     ---------
!     Arguments
!     ---------
!
      CLASS(SMCurve), POINTER    :: curve              ! The curve to be approximated
      INTEGER                    :: polyOrder          ! polynomial order of the new approximation
      REAL(KIND=RP)              :: tRange(2)          ! Location of forced cut points, starting with t = 0.0
      TYPE(OptimizerOptions)     :: options            ! parameters for the approximation
      TYPE(GaussQuadratureType)  :: gQuad
      REAL(KIND=RP), ALLOCATABLE :: optimizedCuts(:)
!
!     ---------------
!     Local variables
!     ---------------
!
      LOGICAL                                :: steppingNotDone
      REAL(KIND=RP)                          :: cuts(0:1)
      REAL(KIND=RP)                          :: iterToler = 1.0d-6
      REAL(KIND=RP)                          :: fR, h, f, h2
      REAL(KIND=RP)                          :: tL, tMid, tR
      REAL(KIND=RP)                          :: e, hStar, cFactor, targetError
      INTEGER                                :: k, nCuts

      ALLOCATE(optimizedCuts(0:1))
      optimizedCuts = tRange  ! will be added to as marching progresses
!
!     ---------------------------------------
!     Loop until hit the end at t = tRange(2)
!     ---------------------------------------
!
      tL              = tRange(1)
      tR              = tRange(2)
      cuts            = [tL, tR] ! = [t_{k-1},t_k]
      h               = tR - tL
!
!     ------------------------------------------------------------------
!     If the approximation already satisfies the approximation tolerance
!     then the marchingfunction will be negative. Positive means
!     the tolerance is not met and a search must be done
!     fR = marchingfunction(curve, polyOrder, cuts, options, gQuad)
!     ------------------------------------------------------------------
!
      fR = marchingfunction(curve, polyOrder, cuts, options, gQuad)
      IF(fr .le. 0.0) RETURN
!
!     --------------------------
!     Find the intervals in turn
!     --------------------------
!
      steppingNotDone = .TRUE.
      k               = 0 ! Segment number

      DO WHILE(steppingNotDone)
         k    = k + 1
         tMid = iterate(tL, tR, curve, polyOrder, options, gQuad, iterToler, RIGHT_SIDE)
!
!        -------------------------------------------------
!        converged... add this point to the segments array
!        and prepare for the next
!        -------------------------------------------------
!
         IF (ABS(tMid - tRange(2)) <= iterToler) EXIT
         CALL addCut(cuts = optimizedCuts, valueToAdd = tMid, atIndex = k)
         tL   = tMid
         tR   = tRange(2)
         cuts = [tL,tR]
!
!        ----------------------------------------------------
!        We're done if the remaining segment has small enough
!        errors
!        ----------------------------------------------------
!
         f = marchingFunction(curve, polyOrder, cuts, options, gQuad)
         IF(f .le. 0.0_RP)    EXIT

      END DO !Stepping
!
!     --------------------------------------------------------------------
!     The last segment may be very small since the marching knows nothing
!     about what lies ahead. Use the last error and the fact that the
!     the H^1 error should vary like h^{p} to estimate the segment size
!     needed to produce the desired error. Choose the smaller of this
!     segment size and half of the last two segments to balance the error
!     in the last two segments. Don't bother with changes that are less
!     than 25%, which is a completely arbitrary number.
!     --------------------------------------------------------------------
!
       nCuts       = SIZE(optimizedCuts) - 1
       h2          = optimizedCuts(nCuts) - optimizedCuts(nCuts-2)
       h           = optimizedCuts(nCuts) - optimizedCuts(nCuts-1)
       cuts        = [optimizedCuts(nCuts-1), optimizedCuts(nCuts)]
       targetError = options % safetyFactor*options % toler
       e           = marchingFunction(curve, polyOrder, cuts, options, gQuad) + targetError
       cFactor     = (targetError/(e + 1.0d-14))**(1.0_RP/polyOrder)

       IF ( cFactor > 1.25_RP )     THEN
          hStar       = h*cFactor
          optimizedCuts(nCuts-1) = optimizedCuts(nCuts) - MIN(hStar, 0.5_RP*h2)
       END IF

   END SUBROUTINE FindOptimizedCutsInRange
!
!////////////////////////////////////////////////////////////////////////
!
!
!////////////////////////////////////////////////////////////////////////
!
   REAL(KIND=RP) FUNCTION iterate(tL, tR, curve, polyOrder, options, gQuad, iterToler, side)
   !
   !  -----------------------------------------------------------------------------------
   !  Iterate to find the point t_{k+1} such that the error on [t_k,t_{k+1}] is less than
   !  iterToler.
   !  -----------------------------------------------------------------------------------
   !
         IMPLICIT NONE
   !
   !     ---------
   !     Arguments
   !     ---------
   !
         REAL(KIND=RP)             :: tL, tR
         CLASS(SMCurve), POINTER   :: curve              ! The curve to be approximated
         INTEGER                   :: polyOrder          ! polynomial order of the new approximation
         TYPE(OptimizerOptions)    :: options            ! parameters for the approximation
         TYPE(GaussQuadratureType) :: gQuad
         REAL(KIND=RP)             :: iterToler
         INTEGER                   :: side
   !
   !     ---------------
   !     Local variables
   !     ---------------
   !
         REAL(KIND=RP)   :: cuts(0:1)
         REAL(KIND=RP)   :: fL, fR
         REAL(KIND=RP)   :: tN, tNm1, fN, fNm1, tNp1, fNp1
         INTEGER         :: maxIterB = 6 ! Should get to 0.016
         INTEGER         :: maxIterS = 8 ! Should get to iterToler
         INTEGER         :: maxIter
         INTEGER         :: n
         LOGICAL         :: done

         cuts = [tL, tR] ! = [t_{k-1},t_k]
         iterate = bisect(maxIterB, side, done)
         IF ( done ) RETURN
   !
   !     -----------------------------------
   !     Pegasus method to finish
   !     (fast modification of Regula Falsi)
   !
   !     See paper
   !     Dowell, M., Jarratt, P. (1972)
   !     The “Pegasus” method for computing
   !       the root of an equation.
   !     https://doi.org/10.1007/BF01932959
   !     -----------------------------------
   !
         tN   = tR; fN   = fR
         tNm1 = tL; fNm1 = fL

       ! Important note. This assumes that
         ! fN and fNm1 have opposite sign, otherwise it won't bracket.
         DO n = 1, maxIterS
          ! Regula Falsi position
            ! (Intersection of secant line with the t axis)
          ! tNp1 = tN - (tN - tNm1) * fN / (fN - fNm1)
            tNp1 = (fN * tNm1 - tN * fNm1) / (fN - fNm1)
            cuts(1) = tNp1
            fNp1    = marchingFunction(curve, polyOrder, cuts, options, gQuad)

            ! Check for convergence
            IF ( ( ABS(fNp1) <= iterToler ) .OR. ( ABS(tN - tNm1) <= iterToler ) )    THEN
               iterate = tNp1
               RETURN
            END IF

            ! Update endpoints
            IF ( fNm1 * fNp1 < 0.0_RP )    THEN
               ! tNm1 stays the same
               ! Pegasus scaling of left function value
               fNm1 = fNm1 * fN / (fN + fNp1)

               ! Replace right endpoint
               tN = tNp1
               fN = fNp1
            ELSE
               ! tN stays the same
               ! Pegasus scaling of right function value
               fN = fNm1 * fN / (fNm1 + fNp1)

               ! Replace left endpoint
               tNm1 = tNp1
               fNm1 = fNp1
            END IF
         END DO
   !
   !     ---------------------------------------------------------------
   !     If Pegasus doesn't converge, we get to here. Just use bisection
   !     to finish the job.
   !     ---------------------------------------------------------------
   !
         maxIter = NINT(LOG((tR - tL)/iterToler)/LOG(2.0_RP))+1
         iterate = bisect(maxIter, side, done)
   !
   !     ========
         CONTAINS
   !     ========
   !
   !////////////////////////////////////////////////////////////////////////
   !
      REAL(KIND=RP) FUNCTION bisect(maxIter, side, done)
   !
   !  -----------------------------------------------------------------------------------
   !  Iterate to find the point t_{k+1} such that the error on [t_k,t_{k+1}] is less than
   !  iterToler.
   !  -----------------------------------------------------------------------------------
   !
         IMPLICIT NONE
   !
   !     ---------
   !     Arguments
   !     ---------
   !
         INTEGER                   :: maxIter
         LOGICAL                   :: done
         INTEGER                   :: side
   !
   !     ---------------
   !     Local variables
   !     ---------------
   !
         REAL(KIND=RP)   :: tMid, h
         REAL(KIND=RP)   :: fMid
         INTEGER         :: n

   !
   !    -----------------------------------------------------
   !    Iterate on error function until tolerance is satisfied
   !    -----------------------------------------------------
   !
         fL = - options % safetyFactor*options % toler ! the best case. See marchingFunction
   !
   !     ---------
   !     Bisection
   !     ---------
   !
         done = .FALSE.
         DO n = 1, maxIter
            h       = tR - tL
            tMid    = tL + 0.5_RP*h
            cuts(side) = tMid
            fMid = marchingFunction(curve, polyOrder, cuts, options, gQuad)

            IF(0.5_RP*ABS(h) .le. iterToler) THEN ! If marching R -> L, h < 0
               bisect = tMid
               done   = .TRUE.
               RETURN
            END IF

            IF ( fMid*fL < 0.0_RP )     THEN
               tR = tMid
               fR = fMid
            ELSE
               tL = tMid
               fL = fMid
            END IF
         END DO

         bisect = tMid

      END FUNCTION bisect

      END FUNCTION iterate
!
!////////////////////////////////////////////////////////////////////////
!
   FUNCTION marchingFunction(curve, polyOrder, cuts, options, gQuad) RESULT(f)
      IMPLICIT NONE
!
!     ----------
!     Arguments
!     ----------
!
      CLASS(SMCurve), POINTER    :: curve              ! The curve to be approximated
      INTEGER                    :: polyOrder          ! polynomial order of the new approximation
      REAL(KIND=RP)              :: cuts(0:1)
      TYPE(OptimizerOptions)     :: options            ! parameters for the marching algorithm
      TYPE(GaussQuadratureType)  :: gQuad
!
!     ------------
!     Return value
!     ------------
!
      REAL(KIND=RP) :: f
!
!     ---------------
!     Local variables
!     ---------------
!
      REAL(KIND=RP), ALLOCATABLE             :: errors(:,:)
      CALL curveSegmentErrors(curve, polyOrder, cuts, options, gQuad, errors)
      f = errors(1,USER_NORM) - options % safetyFactor*options % toler

   END FUNCTION marchingFunction
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE curveSegmentErrors(curve, polyOrder, cuts, options, gQuad, errors)
      IMPLICIT NONE
!
!     ----------
!     Arguments
!     ----------
!
      CLASS(SMCurve), POINTER    :: curve              ! The curve to be approximated
      INTEGER                    :: polyOrder          ! polynomial order of the new approximation
      REAL(KIND=RP)              :: cuts(0:1)
      TYPE(GaussQuadratureType)  :: gQuad
      REAL(KIND=RP), ALLOCATABLE :: errors(:,:)
!
!     ---------------
!     Local variables
!     ---------------
!
      CLASS(SMCurve)               , POINTER :: optimizedCurve
      CLASS(MultiSegmentModalCurve), POINTER :: msmCurve
      TYPE(OptimizerOptions)                 :: options  ! parameters for the marching algorithm
      INTEGER, ALLOCATABLE                   :: breakIndices(:)

      breakIndices = [0]
      CALL OptimizeCurve(curve        = curve,                   &
                         polyOrder    = polyOrder,               &
                         cuts         = cuts,                    &
                         breakIndices = breakIndices,            &
                         options      = options,                 &
                         newName      = "tmp",                   &
                         newID        = 0,                       &
                         optimized    = optimizedCurve)

      CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
      CALL SegmentErrors(exact          = curve,    &
                         segmentedCurve = msmCurve, &
                         gQuad          = gQuad,    &
                         errors         = errors)

      CALL releaseBaseCurve(optimizedCurve)

   END SUBROUTINE curveSegmentErrors
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE curveSegmentLengths(curve, polyOrder, cuts, options, gQuad, lengths)
!
!  -----------------------------------------------------
!  Collect the arc lengths of the optimal curve segments
!  -----------------------------------------------------
!
      IMPLICIT NONE
!
!     ----------
!     Arguments
!     ----------
!
      CLASS(SMCurve), POINTER    :: curve              ! The curve to be approximated
      INTEGER                    :: polyOrder          ! polynomial order of the new approximation
      REAL(KIND=RP)              :: cuts(0:)
      TYPE(GaussQuadratureType)  :: gQuad
      REAL(KIND=RP), ALLOCATABLE :: lengths(:)
!
!     ---------------
!     Local variables
!     ---------------
!
      CLASS(SMCurve)               , POINTER :: optimizedCurve
      CLASS(MultiSegmentModalCurve), POINTER :: msmCurve
      TYPE(OptimizerOptions)                 :: options  ! parameters for the marching algorithm
      INTEGER, ALLOCATABLE                   :: breakIndices(:)
      INTEGER                                :: nSegments
      INTEGER                                :: j
!
!     ------------------------------
!     find the optimal approximation
!     ------------------------------
!
      breakIndices = [0]
      CALL OptimizeCurve(curve        = curve,                   &
                         polyOrder    = polyOrder,               &
                         cuts         = cuts,                    &
                         breakIndices = breakIndices,            &
                         options      = options,                 &
                         newName      = "tmp",                   &
                         newID        = 0,                       &
                         optimized    = optimizedCurve)

      CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
!
!     ---------------------------------------
!     Compute the arc lengths of each segment
!     ---------------------------------------
!
      nSegments = msmCurve % nSegments
      ALLOCATE(lengths(nSegments))
      DO j = 1, nSegments
         lengths(j) = segmentArcLength(msmCurve, j ,gQuad % nodes, gQuad % weights)
      END DO
      CALL releaseBaseCurve(optimizedCurve)

   END SUBROUTINE curveSegmentLengths
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE SegmentErrors(exact, segmentedCurve, gQuad, errors)
!
!  -------------------------------------------------------------------
!  Compute the H^1 errors in each of the segments of a segmented curve
!  -------------------------------------------------------------------
!
      IMPLICIT NONE
!
!     ---------
!     Arguments
!     ---------
!
      CLASS(SMCurve)               , POINTER :: exact
      CLASS(MultiSegmentModalCurve), POINTER :: segmentedCurve
      TYPE(GaussQuadratureType)              :: gQuad
      REAL(KIND=RP), ALLOCATABLE             :: errors(:,:)
!
!     ----------------
!     Local variables
!     ----------------
!
      INTEGER       :: nSegments
      INTEGER       :: k

      nSegments = segmentedCurve % nSegments
      ALLOCATE(errors(nSegments,3), source = 0.0_RP)

      DO k = 1, nSegments
         errors(k,:) =  segmentError(k, exact, segmentedCurve, gQuad)
      END DO

   END SUBROUTINE SegmentErrors
!
!////////////////////////////////////////////////////////////////////////
!
   FUNCTION segmentError(k, exact, segmentedCurve, gQuad) RESULT(er)
!
!  ---------------------------------------------------------------
!  Compute the H^1 errors on the k'th segment of a segmented curve
!  ---------------------------------------------------------------
!
      IMPLICIT NONE
!
!     ---------
!     Arguments
!     ---------
!
      INTEGER                                :: k
      CLASS(SMCurve)               , POINTER :: exact
      CLASS(MultiSegmentModalCurve), POINTER :: segmentedCurve
      TYPE(GaussQuadratureType)              :: gQuad
      REAL(KIND=RP)                          :: er(3)
!
!     ----------------
!     Local variables
!     ----------------
!
      INTEGER       :: qOrder
      INTEGER       :: j
      REAL(KIND=RP) :: t, h, dsdt
      REAL(KIND=RP) :: e, e0(3), e1(3), eMax, eMaxDeriv

      qOrder = gQuad % N

      e         = 0.0_RP
      eMax      = 0.0_RP
      eMaxDeriv = 0.0_RP
      h         = segmentedCurve % cuts(k) - segmentedCurve % cuts(k-1)
      dsdt      = 2.0_RP/h
      DO j = 0, qOrder
         t    = segmentedCurve % cuts(k-1) + h*0.5_RP*(gQuad % nodes(j) + 1.0_RP)
         e0   = (exact % positionAt(t)   - segmentedCurve % valueInSegment(k, t, which = LA_EVALUATE_FUNCTION))**2
         e1   = (exact % derivativeAt(t) - dsdt*segmentedCurve % valueInSegment(k, t, which = LA_EVALUATE_DERIVATIVE))**2

         e         = e + (e0(1) + e0(2) + e1(1) + e1(2))*gQuad % weights(j)
         eMax      = MAX(eMax,e0(1),e0(2))
         eMaxDeriv = MAX(eMaxDeriv,e1(1),e1(2))
      END DO
      er(USER_NORM)       = SQRT(0.5_RP*h*e)
      er(MAX_NORM)        = SQRT(eMax)
      er(MAX_NORM_DERIV)  = SQRT(eMaxDeriv)

   END FUNCTION segmentError
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE addCut(cuts, valueToAdd, atIndex)
!
!     ----------------------------------------------------------
!     Cuts is dimensioned as (0:N), but there seems to be no way
!     with an allocatable to note that in the declaration here.
!     So work with shifted indices, shifted by 1
!     ----------------------------------------------------------
!
      IMPLICIT NONE
!
!     ---------
!     Arguments
!     ---------
!
      REAL(KIND=RP), ALLOCATABLE :: cuts(:)
      REAL(KIND=RP) :: valueToAdd
      INTEGER       :: atIndex
!
!     ---------------
!     Local variables
!     ---------------
!
      REAL(KIND=RP), ALLOCATABLE :: newCuts(:)
      INTEGER                    :: N

      N = SIZE(cuts)
      ALLOCATE(newCuts(0:N))
      newCuts(0:atIndex-1)  = cuts(0:atIndex-1)
      newCuts(atIndex+1:)   = cuts(atIndex:)
      newCuts(atIndex)      = valueToAdd

      CALL MOVE_ALLOC(FROM = newCuts, TO = cuts)

   END SUBROUTINE addCut
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE AppendRealArray(base, arrayToAdd)
      IMPLICIT NONE
      REAL(KIND=RP), ALLOCATABLE :: base(:)
      REAL(KIND=RP)              :: arrayToAdd(:)
      base = [base,arrayToAdd]

   END SUBROUTINE AppendRealArray
!
!////////////////////////////////////////////////////////////////////////
!
!   SUBROUTINE OptimizeCurveWithSubdivisions(curve, polyOrder, options, newName, newID, optimized)
!      IMPLICIT NONE
!!
!!     ---------
!!     Arguments
!!     ---------
!!
!      CLASS(SMCurve), POINTER :: curve              ! The curve to be approximated
!      CLASS(SMCurve), POINTER :: optimized          ! The new approximation
!      INTEGER                 :: polyOrder          ! polynomial order of the new approximation
!      CHARACTER(LEN=*)        :: newName            ! Name of the new curve
!      INTEGER                 :: newID              ! ID of the new curve
!      TYPE(OptimizerOptions)  :: options            ! parameters for the approximation
!!
!!     ---------------
!!     Local variables
!!     ---------------
!!
!      INTEGER                                :: lev, j, i
!      INTEGER                                :: offset
!      INTEGER                                :: nSegments
!      REAL(KIND=RP), ALLOCATABLE             :: cuts(:)
!      CLASS(SMCurve)               , POINTER :: optimizedCurve
!      CLASS(MultiSegmentModalCurve), POINTER :: msmCurve
!      TYPE(GaussQuadratureType)              :: gQuad
!      REAL(KIND=RP), ALLOCATABLE             :: errors(:)
!      REAL(KIND=RP)                          :: tMid
!      LOGICAL                                :: subdivisionIsFinished
!
!      CALL ConstructGaussQuadrature(gQuad, 4*polyOrder)
!!
!!     ------------------
!!     Start at top level
!!     ------------------
!!
!      nSegments = 1
!      ALLOCATE(cuts(0:nSegments))
!      cuts = [0.0_RP, 1.0_RP]
!
!      CALL OptimizeCurve(curve     = curve,                       &
!                         polyOrder = polyOrder,                   &
!                         cuts      = cuts,                        &
!                         options   = options,                     &
!                         newName   = newName,                     &
!                         newID     = newID,                       &
!                         optimized = optimizedCurve)
!
!      CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
!      CALL SegmentErrors(exact          = curve,    &
!                         segmentedCurve = msmCurve, &
!                         gQuad          = gQuad,    &
!                         errors         = errors)
!      DO lev = 1, MAX_DEPTH
!         subdivisionIsFinished = .TRUE.
!!
!!        -------------------------------------
!!        Compute new subdivisions as necessary
!!        -------------------------------------
!!
!         offset = 0
!         DO j = 1, UBOUND(errors,1)
!            i = j + offset
!            IF ( errors(j) > options % toler )     THEN
!               tMid = 0.5_RP*(cuts(i) + cuts(i-1))
!               CALL addCut(cuts,valueToAdd = tMid, atIndex = i)
!               offset = offset + 1
!               subdivisionIsFinished = .FALSE.
!            END IF
!         END DO
!         IF(subdivisionIsFinished) EXIT
!!
!!        ----------------------------------
!!        Compute curves on new subdivisions
!!        ----------------------------------
!!
!         CALL releaseBaseCurve(optimizedCurve)
!         DEALLOCATE(errors)
!         CALL OptimizeCurve(curve              = curve,              &
!                            polyOrder          = polyOrder,          &
!                            cuts               = cuts,               &
!                            options            = options,            &
!                            newName            = newName,            &
!                            newID              = newID,              &
!                            optimized          = optimizedCurve)
!
!         CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
!         CALL SegmentErrors(exact          = curve,    &
!                            segmentedCurve = msmCurve, &
!                            gQuad          = gQuad,    &
!                            errors         = errors)
!
!      END DO
!      optimized => optimizedCurve
!
!   END SUBROUTINE OptimizeCurveWithSubdivisions

   END Module CurveOptimization

trixi-framework / HOHQMesh / 30318764346

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous