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

90.14

/Source/BoundaryOptimization/OptimizerTests.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
!
!////////////////////////////////////////////////////////////////////////
!
!      OptimizerTests.f90
!      Created: November 22, 2025 at 8:58 AM
!      By: David Kopriva
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE polynomialApproximationTest
!
!  -------------------------------------------------------------------
!  Approximate simple basis functions with a three segment polynomial
!  This gives an example of how to use the optimizer and multiSegment
!  modal curve classes. The return value fail can be used for testing.
!  -------------------------------------------------------------------
!
      USE CurveOptimization
      USE SMParametricEquationCurveClass
      USE FTAssertions
      IMPLICIT NONE

      INTEGER                                  :: optType = H1_NORM! L2_NORM or H1_NORM

      CLASS(SMCurve)                 , POINTER :: crv
      TYPE(SMParametricEquationCurve), POINTER :: simpleCurve
      CLASS(SMCurve)                 , POINTER :: optimizedCurve
      CLASS(MultiSegmentModalCurve)  , POINTER :: msmCurve
      TYPE(GaussQuadratureType)                :: gQuad
      TYPE(OptimizerOptions)                   :: options

      INTEGER                    :: nSegments
      INTEGER                    :: polyOrder
      REAL(KIND=RP), ALLOCATABLE :: cuts(:)
      REAL(KIND=RP), ALLOCATABLE :: errors(:,:)
      INTEGER      , ALLOCATABLE :: breakIndices(:)
      CHARACTER(LEN=64)          :: xEqn, yEqn, zEqn
      REAL(KIND=RP)              :: tol = 1.0d-7

      CALL SetDefaultOptions(options)
      options % whichNorm = optType
!
!     --------------------
!     Curve to approximate
!     --------------------
!
      xEqn = "x(t) = 2*t-1"                   ! mapped L_0
      yEqn = "y(t) = 0.5*(3*(2*t-1)^2 - 1)"   ! mapped L_1
      zEqn = "z(t) = 0.0"

      ALLOCATE(simpleCurve)
      CALL simpleCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "simpleCurve", id = 1)
      crv => simpleCurve
!
!     ----------
!     Optimizing
!     ----------
!
      nSegments       = 3
      polyOrder       = 4
      ALLOCATE(cuts(0:nSegments))
      cuts         = [0.0_RP, 0.3_RP, 0.6_RP, 1.0_RP] ! for SMCurves, t\in [0,1]
      breakIndices = [0]

      CALL OptimizeCurve(curve              = crv,                   &
                         polyOrder          = polyOrder,             &
                         cuts               = cuts,                  &
                         breakIndices       = breakIndices,          &
                         options            = options,               &
                         newName            = "OptimizedSimpleCurve",&
                         newID              = simpleCurve % id() + 1,&
                         optimized          = optimizedCurve)
!
!     --------------
!     Compute errors
!     --------------
!
      CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
      CALL ConstructGaussQuadrature(gQuad, 2*polyOrder)
      CALL SegmentErrors(exact          = crv, &
                         segmentedCurve = msmCurve, &
                         gQuad          = gQuad, &
                         errors         = errors)
      CALL FTAssert(MAXVAL(errors(:,USER_NORM)) .le. tol, msg = "PolynomialApproximationTest")
!
!     --------
!     Clean up
!     --------
!
      CALL releaseBaseCurve(optimizedCurve)
      CALL releasePECurve(simpleCurve)

   END SUBROUTINE PolynomialApproximationTest
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE BlobCurveTest()
!
!  -------------------------------------------------------------------
!  Ensure that the marching optimization is consistent with previous
!  results.
!  -------------------------------------------------------------------
!
      USE CurveOptimization
      USE SMParametricEquationCurveClass
      USE FTAssertions
      IMPLICIT NONE
!
!     ---------------
!     Test parameters
!     ---------------
!
      INTEGER                                :: optType    = H1_NORM
      INTEGER                                :: polyOrder  = 6
      INTEGER                                :: smoothness = 2 ! C^2 smoothness
      REAL(KIND=RP)                          :: toler      = 0.1_RP
      REAL(KIND=RP)                          :: testTol    = 2.0d-9
!
!     ---------------
!     Local Variables
!     ---------------
!
      CLASS(SMCurve)                 , POINTER :: crv
      TYPE(SMParametricEquationCurve), POINTER :: blobCurve
      CLASS(SMCurve)                 , POINTER :: optimizedCurve
      CLASS(MultiSegmentModalCurve)  , POINTER :: msmCurve
      TYPE(GaussQuadratureType)                :: gQuad
      TYPE(OptimizerOptions)                   :: options
      CHARACTER(LEN=64)                        :: xEqn, yEqn, zEqn
      REAL(KIND=RP), ALLOCATABLE               :: errors(:,:)
      INTEGER      , ALLOCATABLE               :: breakIndices(:)
!
!     ------------------------------------------------------
!     If algorithms change, these values will have to change
!     ------------------------------------------------------
!
      REAL(KIND=RP), DIMENSION(0:14) :: targetCuts = [0.0000000000000000d0, 8.1265069873694071d-002, 0.15045565047752610d0, &
                                                      0.22077110332899968d0, 0.29160293158998757d0 , 0.36571518148058291d0, &
                                                      0.43544720244301666d0, 0.50460363752397241d0 , 0.58656575126088406d0, &
                                                      0.66015636956282953d0, 0.72937540809332357d0 , 0.79962651199484203d0, &
                                                      0.87133816457876168d0, 0.93566908228937584d0 , 1.0000000000000000d0]
      REAL(KIND=RP), DIMENSION(14) :: targetErrs   = [6.5702408339004884d-002, 9.2756204392096850d-002, 9.4318406193022147d-002, &
                                                      8.9851311172061341d-002, 8.4840693387210162d-002, 9.1113841663074166d-002, &
                                                      8.9981484563422739d-002, 6.6436074374897389d-002, 8.5125646477200778d-002, &
                                                      9.4312464884582489d-002, 9.2294291838550244d-002, 7.9708305374154340d-002, &
                                                      6.9615346724282365d-002, 6.0416555870896144d-002]
!
!     --------------
!     Set up options
!     --------------
!
      CALL SetDefaultOptions(options)
      options % internalConstraint = smoothness
      options % toler              = toler
      options % whichNorm          = optType
!
!     -------------------
!     Create Winters Blob
!     -------------------
!
      xEqn = "x(t) = 4*cos(2*pi*t) - 3/5*cos(8*pi*t)^3"
      yEqn = "y(t) = 4*sin(2*pi*t) - 0.5*sin(11*pi*t)^2"
      zEqn = "z(t) = 0.0"
      ALLOCATE(blobCurve)
      CALL blobCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "blobCurve", id = 1)
      crv => blobCurve
!
!     ------------------------------------
!     Find the optimal curve approximation
!     ------------------------------------
!
      CALL OptimizeCurveByMarching(curve              = crv,             &
                                   polyOrder          = polyOrder,       &
                                   breaks             = [0.0_RP,1.0_RP], &
                                   breakIndices       = breakIndices,    &
                                   options            = options,         &
                                   newName            = "BlobTest",      &
                                   newID              = crv % id() + 1,  &
                                   optimized          = optimizedCurve)
!
!     --------------
!     Compute errors
!     --------------
!
      CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
      CALL ConstructGaussQuadrature(gQuad, 2*polyOrder)
      CALL SegmentErrors(exact          = crv, &
                         segmentedCurve = msmCurve, &
                         gQuad          = gQuad, &
                         errors         = errors)
!
!     ------------
!     Check errors
!     ------------
!
      CALL FTAssert(MAXVAL(ABS(errors(:,USER_NORM) - targetErrs)) .le. testTol, msg = "Errors dont match")
      CALL FTAssert(MAXVAL(ABS(targetCuts - msmCurve % cuts))     .le. testTol, msg = "Segments dont match")
!
!     --------
!     Clean up
!     --------
!
      CALL releaseBaseCurve(optimizedCurve)
      CALL releasePECurve(blobCurve)

   END SUBROUTINE BlobCurveTest
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE BlobBreaksTest()
!
!  -------------------------------------------------------------------
!  Ensure that the marching optimization is consistent with previous
!  results.
!  -------------------------------------------------------------------
!
      USE CurveOptimization
      USE SMParametricEquationCurveClass
      USE FTAssertions
      IMPLICIT NONE
!
!     ---------------
!     Test parameters
!     ---------------
!
      INTEGER                                :: optType    = H1_NORM
      INTEGER                                :: polyOrder  = 6
      INTEGER                                :: smoothness = 2 ! C^2 smoothness
      REAL(KIND=RP)                          :: toler      = 0.1_RP, testTol = 2.0d-9
      REAL(KIND=RP)                          :: testBreaks(0:3) = [0.0_RP, 0.4_RP, 0.7_RP, 1.0_RP]
!
!     ---------------
!     Local Variables
!     ---------------
!
      CLASS(SMCurve)                 , POINTER :: crv
      TYPE(SMParametricEquationCurve), POINTER :: blobCurve
      CLASS(SMCurve)                 , POINTER :: optimizedCurve
      CLASS(MultiSegmentModalCurve)  , POINTER :: msmCurve
      TYPE(GaussQuadratureType)                :: gQuad
      TYPE(OptimizerOptions)                   :: options
      CHARACTER(LEN=64)                        :: xEqn, yEqn, zEqn
      REAL(KIND=RP), ALLOCATABLE               :: optimizedCuts(:)
      INTEGER      , ALLOCATABLE               :: breakIndices(:)
!
!     ------------------------------------------------------
!     If algorithms change, these values will have to change
!     ------------------------------------------------------
!
      REAL(KIND=RP) :: targetCuts(0:13) = [0.0000000000000000d0 , 9.7013000086126913d-002, 0.18597895538838266d0, &
                                           0.27478075822849279d0, 0.33739037911424141d0  , 0.40000000000000002d0, &
                                           0.48775051048802515d0, 0.57664143712361982d0  , 0.63832071856180495d0, &
                                           0.69999999999999996d0, 0.78923224968415717d0  , 0.87992823095351413d0, &
                                           0.93996411547675207d0, 1.0000000000000000d0]
      INTEGER       :: targetIndices(3) = [6, 10, 14]
!
!     --------------
!     Set up options
!     --------------
!
      CALL SetDefaultOptions(options)
      options % internalConstraint = smoothness
      options % toler              = toler
      options % whichNorm          = optType
!
!     -------------------
!     Create Winters Blob
!     -------------------
!
      xEqn = "x(t) = 4*cos(2*pi*t) - 3/5*cos(8*pi*t)^3"
      yEqn = "y(t) = 4*sin(2*pi*t) - 0.5*sin(11*pi*t)^2"
      zEqn = "z(t) = 0.0"
      ALLOCATE(blobCurve)
      CALL blobCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "blobCurve", id = 1)
      crv => blobCurve
!
!     ------------------------------------
!     Find the optimal curve approximation
!     ------------------------------------
!
      CALL ConstructGaussQuadrature(gQuad, 4*polyOrder) ! For error computation. The 4 is arbitrary
      CALL FindOptimizedCuts(crv, polyOrder, testBreaks, options, gQuad, optimizedCuts, breakIndices)
!
!     ------------
!     Check errors
!     ------------
!
      CALL FTAssert(MAXVAL(ABS(optimizedCuts - targetCuts)) .le. testTol, msg = "Segments dont match")
      CALL FTAssert(MAXVAL(ABS(targetIndices - breakIndices)) == 0, msg = "Break Indices dont match")
!
!     --------
!     Clean up
!     --------
!
      CALL releasePECurve(blobCurve)

   END SUBROUTINE BlobBreaksTest
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE CircleTest()
!
!  -------------------------------------------------------------------
!  Ensure that the marching optimization is consistent with previous
!  results.
!  -------------------------------------------------------------------
!
      USE CurveOptimization
      USE SMParametricEquationCurveClass
      USE FTAssertions
      IMPLICIT NONE
!
!     ---------------
!     Test parameters
!     ---------------
!
      INTEGER                                :: optType    = H1_NORM
      INTEGER                                :: polyOrder  = 6
      INTEGER                                :: smoothness = 2 ! C^2 smoothness
      REAL(KIND=RP)                          :: toler      = 0.10_RP
!
!     ---------------
!     Local Variables
!     ---------------
!
      CLASS(SMCurve)                 , POINTER :: crv
      TYPE(SMParametricEquationCurve), POINTER :: circleCurve
      CLASS(SMCurve)                 , POINTER :: optimizedCurve
      CLASS(MultiSegmentModalCurve)  , POINTER :: msmCurve
      TYPE(OptimizerOptions)                   :: options
      TYPE(GaussQuadratureType)                :: gQuad
      CHARACTER(LEN=64)                        :: xEqn, yEqn, zEqn
      REAL(KIND=RP)                            :: arcL
      INTEGER                    , ALLOCATABLE :: breakIndices(:)
      REAL(KIND=RP)              , ALLOCATABLE :: aLengths(:)

!
!     --------------
!     Set up options
!     --------------
!
      CALL SetDefaultOptions(options)
      options % internalConstraint = smoothness
      options % toler              = toler
      options % whichNorm          = optType
!
!     ---------------
!     Create a Circle
!     ---------------
!
      xEqn = "x(t) = 4*cos(2*pi*t)"
      yEqn = "y(t) = 4*sin(2*pi*t)"
      zEqn = "z(t) = 0.0"
      ALLOCATE(circleCurve)
      CALL circleCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "circleCurve", id = 1)
      crv => circleCurve
!
!     ------------------------------------
!     Find the optimal curve approximation
!     ------------------------------------
!
      CALL OptimizeCurveByMarching(curve              = crv,             &
                                   polyOrder          = polyOrder,       &
                                   breaks             = [0.0_RP,1.0_RP], &
                                   breakindices       = breakIndices,    &
                                   options            = options,         &
                                   newName            = "CircleTest",    &
                                   newID              = crv % id() + 1,  &
                                   optimized          = optimizedCurve)
!
!     ------------------------
!     Compute total arc length
!     ------------------------
!
      CALL ConstructGaussQuadrature(gQuad, 4*polyOrder) ! The 4 is arbitrary
      CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)
      arcL = msmCurve % arcLength(gQuad)
      CALL FTAssertEqual(expectedValue = 8.0_RP*PI, &
                         actualValue   = arcL,  &
                         relTol        = 1.0d-4, &
                         msg           = "Arc Length estimation")
!
!     -------------------------------
!     Compute arc lengths of segments
!     -------------------------------
!
      CALL curveSegmentLengths(crv, polyOrder, [0.0_RP, 0.5_RP, 1.0_RP], &
                               options, gQuad, aLengths)
      CALL FTAssertEqual(expectedValue = 4.0_RP*PI, &
                         actualValue   = aLengths(1),  &
                         relTol        = 1.0d-4, &
                         msg           = "Arc Length estimation 1")
      CALL FTAssertEqual(expectedValue = 4.0_RP*PI, &
                         actualValue   = aLengths(2),  &
                         relTol        = 1.0d-4, &
                         msg           = "Arc Length estimation 2")
!
!     --------
!     Clean up
!     --------
!
      CALL releaseBaseCurve(optimizedCurve)
      CALL releasePECurve(circleCurve)

   END SUBROUTINE CircleTest
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE SmoothnessCheck
      USE FTAssertions
      USE ConstrainedMultiH1Optimization
      IMPLICIT NONE

      CALL FTAssert(spotCheckSmoothnessConditionsIsOK(),msg = "Smoothness Continuity array")

   END SUBROUTINE SmoothnessCheck
!
!////////////////////////////////////////////////////////////////////////
!
   SUBROUTINE SegmentedCurveCheck
!
!  -------------------------------------------------
!  Compute optimal segments for a half circle domain
!  -------------------------------------------------
!
      USE SMConstants
      USE SMLineClass
      USE SMChainedCurveClass
      USE SMEllipticArcClass
      USE CurveOptimization
      USE FTExceptionClass
      USE SharedExceptionManagerModule
      USE FTAssertions
      USE, INTRINSIC :: iso_fortran_env, only : stderr => ERROR_UNIT
      IMPLICIT NONE
!
!     --------
!     Geometry
!     --------
!
      CLASS(SMLine)        , POINTER :: line
      CLASS(SMEllipticArc) , POINTER :: circle
      TYPE (SMChainedCurve), POINTER :: chain
      CLASS(SMCurve)       , POINTER :: curvePtr => NULL()
      CLASS(FTObject)      , POINTER :: obj
!
!     ------------
!     Optimization
!     ------------
!
      TYPE(GaussQuadratureType)  :: gQuad
      TYPE(OptimizerOptions)     :: options
      REAL(KIND=RP), ALLOCATABLE :: optimizedCuts(:)
      INTEGER      , ALLOCATABLE :: breakIndices(:)
      REAL(KIND=RP)              :: testBreaks(0:2)   = [0.0_RP, 0.5_RP, 1.0_RP]
      REAL(KIND=RP)              :: expectedCuts(0:7) = [0.0000000000000000d0 , 8.5660589005650267d-002, 0.17132074027521632d0, &
                                                         0.25698133016711511d0, 0.34264190234426545d0  , 0.42829047404568898d0, &
                                                         0.50000000000000000d0, 1.0000000000000000d0]
      INTEGER                    :: expectedIndices(2) = [7, 8]

!
!     ---------------
!     Test parameters
!     ---------------
!
      INTEGER       :: optType    = H1_NORM
      INTEGER       :: polyOrder  = 3
      INTEGER       :: smoothness = 0 ! C^0 smoothness
      REAL(KIND=RP) :: toler      = 0.0010_RP, testTol = 2.0d-9
!
!     -----
!     Other
!     -----
!
      TYPE(FTException), POINTER :: exception
      INTEGER                    :: errorSeverity = FT_ERROR_NONE

      CALL SetDefaultOptions(options)
      options % internalConstraint = smoothness
      options % toler              = toler
      options % whichNorm          = optType
!
!     --------
!     Geometry
!     --------
!
      ALLOCATE(line)
      CALL line % initWithStartEndNameAndID(xStart = [-1.0_RP,0.0_RP,0.0_RP], &
                                            xEnd   = [1.0_RP,0.0_RP,0.0_RP],  &
                                            cName  = "line", id = 1)
      ALLOCATE(circle)
      CALL circle % initWithParametersNameAndID(center     = [0.0_RP, 0.0_RP, 0.0_RP], &
                                                radius     = 1.0_RP,                   &
                                                startAngle = 0.0_RP,                   &
                                                endAngle   = PI,                       &
                                                cName      = "circle",                 &
                                                id         = 2)
      ALLOCATE(chain)
      CALL chain % initChainWithNameAndID(chainName = "chain",id = 3)
!
      curvePtr => circle
      CALL chain % addCurve(curvePtr)
      obj => circle
      CALL release(obj)

      curvePtr => line
      CALL chain % addCurve(curvePtr)
      obj => line
      CALL release(obj)
!
      CALL chain % complete(innerOrOuterCurve = OUTER,chainMustClose = .TRUE.)

      IF ( catch() )     THEN
         WRITE(stderr,*)
         WRITE(stderr,*)  "------------------------------------------------------------------"
         WRITE(stderr,*)
         WRITE(stderr,*)  "The following errors were found when constructing the check:"

         DO
            exception => popLastException()
            IF ( .NOT.ASSOCIATED(exception) )     EXIT
            CALL exception % printDescription(stderr)
            errorSeverity = MAX(errorSeverity, exception % severity())
         END DO
         WRITE(stderr,*)
         WRITE(stderr,*)  "------------------------------------------------------------------"
         WRITE(stderr,*)

         IF ( errorSeverity > FT_ERROR_WARNING )     THEN
            ERROR STOP "The Errors were Fatal. Cannot generate optimized curve."
         END IF
      END IF
!
!     -------------
!     Optimizations
!     -------------
!
      curvePtr => chain
      CALL ConstructGaussQuadrature(gQuad, 4*polyOrder) ! For error computation. The 4 is arbitrary
      CALL FindOptimizedCuts(curvePtr, polyOrder, testBreaks, options, gQuad, optimizedCuts, breakIndices)
!
!     ------------
!     Check errors
!     ------------
!
      CALL FTAssert(MAXVAL(ABS(optimizedCuts - expectedCuts)) .le. testTol, msg = "Segments dont match")
      CALL FTAssert(MAXVAL(ABS(expectedIndices - breakIndices)) == 0      , msg = "Break Indices dont match")

      CALL releaseChainedCurve(chain)

   END SUBROUTINE SegmentedCurveCheck


1	! MIT License
2	!
3	! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
4	!
5	! Permission is hereby granted, free of charge, to any person obtaining a copy
6	! of this software and associated documentation files (the "Software"), to deal
7	! in the Software without restriction, including without limitation the rights
8	! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9	! copies of the Software, and to permit persons to whom the Software is
10	! furnished to do so, subject to the following conditions:
11	!
12	! The above copyright notice and this permission notice shall be included in all
13	! copies or substantial portions of the Software.
14	!
15	! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16	! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17	! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18	! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19	! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20	! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21	! SOFTWARE.
22	!
23	! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
24	! public domain software:
25	! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
26	! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
27	!
28	! --- End License
29	!
30	!////////////////////////////////////////////////////////////////////////
31	!
32	! OptimizerTests.f90
33	! Created: November 22, 2025 at 8:58 AM
34	! By: David Kopriva
35	!
36	!////////////////////////////////////////////////////////////////////////
37	!
38	SUBROUTINE polynomialApproximationTest	1✔
39	!
40	! -------------------------------------------------------------------
41	! Approximate simple basis functions with a three segment polynomial
42	! This gives an example of how to use the optimizer and multiSegment
43	! modal curve classes. The return value fail can be used for testing.
44	! -------------------------------------------------------------------
45	!
46	USE CurveOptimization
47	USE SMParametricEquationCurveClass
48	USE FTAssertions
49	IMPLICIT NONE
50
51	INTEGER :: optType = H1_NORM! L2_NORM or H1_NORM
52
53	CLASS(SMCurve) , POINTER :: crv
54	TYPE(SMParametricEquationCurve), POINTER :: simpleCurve
55	CLASS(SMCurve) , POINTER :: optimizedCurve
56	CLASS(MultiSegmentModalCurve) , POINTER :: msmCurve
57	TYPE(GaussQuadratureType) :: gQuad	1✔
58	TYPE(OptimizerOptions) :: options
59
60	INTEGER :: nSegments
61	INTEGER :: polyOrder
62	REAL(KIND=RP), ALLOCATABLE :: cuts(:)	1✔
63	REAL(KIND=RP), ALLOCATABLE :: errors(:,:)	1✔
64	INTEGER , ALLOCATABLE :: breakIndices(:)	1✔
65	CHARACTER(LEN=64) :: xEqn, yEqn, zEqn
66	REAL(KIND=RP) :: tol = 1.0d-7
67
68	CALL SetDefaultOptions(options)	1✔
69	options % whichNorm = optType	1✔
70	!
71	! --------------------
72	! Curve to approximate
73	! --------------------
74	!
75	xEqn = "x(t) = 2*t-1" ! mapped L_0	1✔
76	yEqn = "y(t) = 0.5(3(2*t-1)^2 - 1)" ! mapped L_1	1✔
77	zEqn = "z(t) = 0.0"	1✔
78
79	ALLOCATE(simpleCurve)	1✔
80	CALL simpleCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "simpleCurve", id = 1)	1✔
81	crv => simpleCurve	1✔
82	!
83	! ----------
84	! Optimizing
85	! ----------
86	!
87	nSegments = 3	1✔
88	polyOrder = 4	1✔
89	ALLOCATE(cuts(0:nSegments))	1✔
90	cuts = [0.0_RP, 0.3_RP, 0.6_RP, 1.0_RP] ! for SMCurves, t\in [0,1]	6✔
91	breakIndices = [0]	2✔
92
93	CALL OptimizeCurve(curve = crv, &
94	polyOrder = polyOrder, &
95	cuts = cuts, &
96	breakIndices = breakIndices, &
97	options = options, &
98	newName = "OptimizedSimpleCurve",&
99	newID = simpleCurve % id() + 1,&
100	optimized = optimizedCurve)	1✔
101	!
102	! --------------
103	! Compute errors
104	! --------------
105	!
106	CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)	1✔
107	CALL ConstructGaussQuadrature(gQuad, 2*polyOrder)	1✔
108	CALL SegmentErrors(exact = crv, &
109	segmentedCurve = msmCurve, &
110	gQuad = gQuad, &
111	errors = errors)	1✔
112	CALL FTAssert(MAXVAL(errors(:,USER_NORM)) .le. tol, msg = "PolynomialApproximationTest")	4✔
113	!
114	! --------
115	! Clean up
116	! --------
117	!
118	CALL releaseBaseCurve(optimizedCurve)	1✔
119	CALL releasePECurve(simpleCurve)	1✔
120
121	END SUBROUTINE PolynomialApproximationTest	1✔
122	!
123	!////////////////////////////////////////////////////////////////////////
124	!
125	SUBROUTINE BlobCurveTest()	1✔
126	!
127	! -------------------------------------------------------------------
128	! Ensure that the marching optimization is consistent with previous
129	! results.
130	! -------------------------------------------------------------------
131	!
132	USE CurveOptimization
133	USE SMParametricEquationCurveClass
134	USE FTAssertions
135	IMPLICIT NONE
136	!
137	! ---------------
138	! Test parameters
139	! ---------------
140	!
141	INTEGER :: optType = H1_NORM
142	INTEGER :: polyOrder = 6
143	INTEGER :: smoothness = 2 ! C^2 smoothness
144	REAL(KIND=RP) :: toler = 0.1_RP
145	REAL(KIND=RP) :: testTol = 2.0d-9
146	!
147	! ---------------
148	! Local Variables
149	! ---------------
150	!
151	CLASS(SMCurve) , POINTER :: crv
152	TYPE(SMParametricEquationCurve), POINTER :: blobCurve
153	CLASS(SMCurve) , POINTER :: optimizedCurve
154	CLASS(MultiSegmentModalCurve) , POINTER :: msmCurve
155	TYPE(GaussQuadratureType) :: gQuad	1✔
156	TYPE(OptimizerOptions) :: options
157	CHARACTER(LEN=64) :: xEqn, yEqn, zEqn
158	REAL(KIND=RP), ALLOCATABLE :: errors(:,:)	1✔
159	INTEGER , ALLOCATABLE :: breakIndices(:)	1✔
160	!
161	! ------------------------------------------------------
162	! If algorithms change, these values will have to change
163	! ------------------------------------------------------
164	!
165	REAL(KIND=RP), DIMENSION(0:14) :: targetCuts = [0.0000000000000000d0, 8.1265069873694071d-002, 0.15045565047752610d0, &
166	0.22077110332899968d0, 0.29160293158998757d0 , 0.36571518148058291d0, &
167	0.43544720244301666d0, 0.50460363752397241d0 , 0.58656575126088406d0, &
168	0.66015636956282953d0, 0.72937540809332357d0 , 0.79962651199484203d0, &
169	0.87133816457876168d0, 0.93566908228937584d0 , 1.0000000000000000d0]
170	REAL(KIND=RP), DIMENSION(14) :: targetErrs = [6.5702408339004884d-002, 9.2756204392096850d-002, 9.4318406193022147d-002, &
171	8.9851311172061341d-002, 8.4840693387210162d-002, 9.1113841663074166d-002, &
172	8.9981484563422739d-002, 6.6436074374897389d-002, 8.5125646477200778d-002, &
173	9.4312464884582489d-002, 9.2294291838550244d-002, 7.9708305374154340d-002, &
174	6.9615346724282365d-002, 6.0416555870896144d-002]
175	!
176	! --------------
177	! Set up options
178	! --------------
179	!
180	CALL SetDefaultOptions(options)	1✔
181	options % internalConstraint = smoothness	1✔
182	options % toler = toler	1✔
183	options % whichNorm = optType	1✔
184	!
185	! -------------------
186	! Create Winters Blob
187	! -------------------
188	!
189	xEqn = "x(t) = 4cos(2pit) - 3/5cos(8pit)^3"	1✔
190	yEqn = "y(t) = 4sin(2pit) - 0.5sin(11pit)^2"	1✔
191	zEqn = "z(t) = 0.0"	1✔
192	ALLOCATE(blobCurve)	1✔
193	CALL blobCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "blobCurve", id = 1)	1✔
194	crv => blobCurve	1✔
195	!
196	! ------------------------------------
197	! Find the optimal curve approximation
198	! ------------------------------------
199	!
200	CALL OptimizeCurveByMarching(curve = crv, &
201	polyOrder = polyOrder, &
202	breaks = [0.0_RP,1.0_RP], &
203	breakIndices = breakIndices, &
204	options = options, &
205	newName = "BlobTest", &
206	newID = crv % id() + 1, &
207	optimized = optimizedCurve)	1✔
208	!
209	! --------------
210	! Compute errors
211	! --------------
212	!
213	CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)	1✔
214	CALL ConstructGaussQuadrature(gQuad, 2*polyOrder)	1✔
215	CALL SegmentErrors(exact = crv, &
216	segmentedCurve = msmCurve, &
217	gQuad = gQuad, &
218	errors = errors)	1✔
219	!
220	! ------------
221	! Check errors
222	! ------------
223	!
224	CALL FTAssert(MAXVAL(ABS(errors(:,USER_NORM) - targetErrs)) .le. testTol, msg = "Errors dont match")	15✔
225	CALL FTAssert(MAXVAL(ABS(targetCuts - msmCurve % cuts)) .le. testTol, msg = "Segments dont match")	16✔
226	!
227	! --------
228	! Clean up
229	! --------
230	!
231	CALL releaseBaseCurve(optimizedCurve)	1✔
232	CALL releasePECurve(blobCurve)	1✔
233
234	END SUBROUTINE BlobCurveTest	1✔
235	!
236	!////////////////////////////////////////////////////////////////////////
237	!
238	SUBROUTINE BlobBreaksTest()	1✔
239	!
240	! -------------------------------------------------------------------
241	! Ensure that the marching optimization is consistent with previous
242	! results.
243	! -------------------------------------------------------------------
244	!
245	USE CurveOptimization
246	USE SMParametricEquationCurveClass
247	USE FTAssertions
248	IMPLICIT NONE
249	!
250	! ---------------
251	! Test parameters
252	! ---------------
253	!
254	INTEGER :: optType = H1_NORM
255	INTEGER :: polyOrder = 6
256	INTEGER :: smoothness = 2 ! C^2 smoothness
257	REAL(KIND=RP) :: toler = 0.1_RP, testTol = 2.0d-9
258	REAL(KIND=RP) :: testBreaks(0:3) = [0.0_RP, 0.4_RP, 0.7_RP, 1.0_RP]
259	!
260	! ---------------
261	! Local Variables
262	! ---------------
263	!
264	CLASS(SMCurve) , POINTER :: crv
265	TYPE(SMParametricEquationCurve), POINTER :: blobCurve
266	CLASS(SMCurve) , POINTER :: optimizedCurve
267	CLASS(MultiSegmentModalCurve) , POINTER :: msmCurve
268	TYPE(GaussQuadratureType) :: gQuad	1✔
269	TYPE(OptimizerOptions) :: options
270	CHARACTER(LEN=64) :: xEqn, yEqn, zEqn
271	REAL(KIND=RP), ALLOCATABLE :: optimizedCuts(:)	1✔
272	INTEGER , ALLOCATABLE :: breakIndices(:)	1✔
273	!
274	! ------------------------------------------------------
275	! If algorithms change, these values will have to change
276	! ------------------------------------------------------
277	!
278	REAL(KIND=RP) :: targetCuts(0:13) = [0.0000000000000000d0 , 9.7013000086126913d-002, 0.18597895538838266d0, &
279	0.27478075822849279d0, 0.33739037911424141d0 , 0.40000000000000002d0, &
280	0.48775051048802515d0, 0.57664143712361982d0 , 0.63832071856180495d0, &
281	0.69999999999999996d0, 0.78923224968415717d0 , 0.87992823095351413d0, &
282	0.93996411547675207d0, 1.0000000000000000d0]
283	INTEGER :: targetIndices(3) = [6, 10, 14]
284	!
285	! --------------
286	! Set up options
287	! --------------
288	!
289	CALL SetDefaultOptions(options)	1✔
290	options % internalConstraint = smoothness	1✔
291	options % toler = toler	1✔
292	options % whichNorm = optType	1✔
293	!
294	! -------------------
295	! Create Winters Blob
296	! -------------------
297	!
298	xEqn = "x(t) = 4cos(2pit) - 3/5cos(8pit)^3"	1✔
299	yEqn = "y(t) = 4sin(2pit) - 0.5sin(11pit)^2"	1✔
300	zEqn = "z(t) = 0.0"	1✔
301	ALLOCATE(blobCurve)	1✔
302	CALL blobCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "blobCurve", id = 1)	1✔
303	crv => blobCurve	1✔
304	!
305	! ------------------------------------
306	! Find the optimal curve approximation
307	! ------------------------------------
308	!
309	CALL ConstructGaussQuadrature(gQuad, 4*polyOrder) ! For error computation. The 4 is arbitrary	1✔
310	CALL FindOptimizedCuts(crv, polyOrder, testBreaks, options, gQuad, optimizedCuts, breakIndices)	1✔
311	!
312	! ------------
313	! Check errors
314	! ------------
315	!
316	CALL FTAssert(MAXVAL(ABS(optimizedCuts - targetCuts)) .le. testTol, msg = "Segments dont match")	15✔
317	CALL FTAssert(MAXVAL(ABS(targetIndices - breakIndices)) == 0, msg = "Break Indices dont match")	4✔
318	!
319	! --------
320	! Clean up
321	! --------
322	!
323	CALL releasePECurve(blobCurve)	1✔
324
325	END SUBROUTINE BlobBreaksTest	1✔
326	!
327	!////////////////////////////////////////////////////////////////////////
328	!
329	SUBROUTINE CircleTest()	1✔
330	!
331	! -------------------------------------------------------------------
332	! Ensure that the marching optimization is consistent with previous
333	! results.
334	! -------------------------------------------------------------------
335	!
336	USE CurveOptimization
337	USE SMParametricEquationCurveClass
338	USE FTAssertions
339	IMPLICIT NONE
340	!
341	! ---------------
342	! Test parameters
343	! ---------------
344	!
345	INTEGER :: optType = H1_NORM
346	INTEGER :: polyOrder = 6
347	INTEGER :: smoothness = 2 ! C^2 smoothness
348	REAL(KIND=RP) :: toler = 0.10_RP
349	!
350	! ---------------
351	! Local Variables
352	! ---------------
353	!
354	CLASS(SMCurve) , POINTER :: crv
355	TYPE(SMParametricEquationCurve), POINTER :: circleCurve
356	CLASS(SMCurve) , POINTER :: optimizedCurve
357	CLASS(MultiSegmentModalCurve) , POINTER :: msmCurve
358	TYPE(OptimizerOptions) :: options
359	TYPE(GaussQuadratureType) :: gQuad	1✔
360	CHARACTER(LEN=64) :: xEqn, yEqn, zEqn
361	REAL(KIND=RP) :: arcL
362	INTEGER , ALLOCATABLE :: breakIndices(:)	1✔
363	REAL(KIND=RP) , ALLOCATABLE :: aLengths(:)	1✔
364
365	!
366	! --------------
367	! Set up options
368	! --------------
369	!
370	CALL SetDefaultOptions(options)	1✔
371	options % internalConstraint = smoothness	1✔
372	options % toler = toler	1✔
373	options % whichNorm = optType	1✔
374	!
375	! ---------------
376	! Create a Circle
377	! ---------------
378	!
379	xEqn = "x(t) = 4cos(2pi*t)"	1✔
380	yEqn = "y(t) = 4sin(2pi*t)"	1✔
381	zEqn = "z(t) = 0.0"	1✔
382	ALLOCATE(circleCurve)	1✔
383	CALL circleCurve % initWithEquationsNameAndID(xEqn, yEqn, zEqn, curveName = "circleCurve", id = 1)	1✔
384	crv => circleCurve	1✔
385	!
386	! ------------------------------------
387	! Find the optimal curve approximation
388	! ------------------------------------
389	!
390	CALL OptimizeCurveByMarching(curve = crv, &
391	polyOrder = polyOrder, &
392	breaks = [0.0_RP,1.0_RP], &
393	breakindices = breakIndices, &
394	options = options, &
395	newName = "CircleTest", &
396	newID = crv % id() + 1, &
397	optimized = optimizedCurve)	1✔
398	!
399	! ------------------------
400	! Compute total arc length
401	! ------------------------
402	!
403	CALL ConstructGaussQuadrature(gQuad, 4*polyOrder) ! The 4 is arbitrary	1✔
404	CALL castToMultiSegmentModalCurve(optimizedCurve, msmCurve)	1✔
405	arcL = msmCurve % arcLength(gQuad)	1✔
406	CALL FTAssertEqual(expectedValue = 8.0_RP*PI, &
407	actualValue = arcL, &
408	relTol = 1.0d-4, &
409	msg = "Arc Length estimation")	1✔
410	!
411	! -------------------------------
412	! Compute arc lengths of segments
413	! -------------------------------
414	!
415	CALL curveSegmentLengths(crv, polyOrder, [0.0_RP, 0.5_RP, 1.0_RP], &
416	options, gQuad, aLengths)	1✔
417	CALL FTAssertEqual(expectedValue = 4.0_RP*PI, &
418	actualValue = aLengths(1), &
419	relTol = 1.0d-4, &
420	msg = "Arc Length estimation 1")	1✔
421	CALL FTAssertEqual(expectedValue = 4.0_RP*PI, &
422	actualValue = aLengths(2), &
423	relTol = 1.0d-4, &
424	msg = "Arc Length estimation 2")	1✔
425	!
426	! --------
427	! Clean up
428	! --------
429	!
430	CALL releaseBaseCurve(optimizedCurve)	1✔
431	CALL releasePECurve(circleCurve)	1✔
432
433	END SUBROUTINE CircleTest	1✔
434	!
435	!////////////////////////////////////////////////////////////////////////
436	!
437	SUBROUTINE SmoothnessCheck	1✔
438	USE FTAssertions
439	USE ConstrainedMultiH1Optimization
440	IMPLICIT NONE
441
442	CALL FTAssert(spotCheckSmoothnessConditionsIsOK(),msg = "Smoothness Continuity array")	1✔
443
444	END SUBROUTINE SmoothnessCheck	1✔
445	!
446	!////////////////////////////////////////////////////////////////////////
447	!
448	SUBROUTINE SegmentedCurveCheck	1✔
449	!
450	! -------------------------------------------------
451	! Compute optimal segments for a half circle domain
452	! -------------------------------------------------
453	!
454	USE SMConstants
455	USE SMLineClass
456	USE SMChainedCurveClass
457	USE SMEllipticArcClass
458	USE CurveOptimization
459	USE FTExceptionClass
460	USE SharedExceptionManagerModule
461	USE FTAssertions
462	USE, INTRINSIC :: iso_fortran_env, only : stderr => ERROR_UNIT
463	IMPLICIT NONE
464	!
465	! --------
466	! Geometry
467	! --------
468	!
469	CLASS(SMLine) , POINTER :: line
470	CLASS(SMEllipticArc) , POINTER :: circle
471	TYPE (SMChainedCurve), POINTER :: chain
472	CLASS(SMCurve) , POINTER :: curvePtr => NULL()
473	CLASS(FTObject) , POINTER :: obj
474	!
475	! ------------
476	! Optimization
477	! ------------
478	!
479	TYPE(GaussQuadratureType) :: gQuad	1✔
480	TYPE(OptimizerOptions) :: options
481	REAL(KIND=RP), ALLOCATABLE :: optimizedCuts(:)	1✔
482	INTEGER , ALLOCATABLE :: breakIndices(:)	1✔
483	REAL(KIND=RP) :: testBreaks(0:2) = [0.0_RP, 0.5_RP, 1.0_RP]
484	REAL(KIND=RP) :: expectedCuts(0:7) = [0.0000000000000000d0 , 8.5660589005650267d-002, 0.17132074027521632d0, &
485	0.25698133016711511d0, 0.34264190234426545d0 , 0.42829047404568898d0, &
486	0.50000000000000000d0, 1.0000000000000000d0]
487	INTEGER :: expectedIndices(2) = [7, 8]
488
489	!
490	! ---------------
491	! Test parameters
492	! ---------------
493	!
494	INTEGER :: optType = H1_NORM
495	INTEGER :: polyOrder = 3
496	INTEGER :: smoothness = 0 ! C^0 smoothness
497	REAL(KIND=RP) :: toler = 0.0010_RP, testTol = 2.0d-9
498	!
499	! -----
500	! Other
501	! -----
502	!
503	TYPE(FTException), POINTER :: exception
504	INTEGER :: errorSeverity = FT_ERROR_NONE
505
506	CALL SetDefaultOptions(options)	1✔
507	options % internalConstraint = smoothness	1✔
508	options % toler = toler	1✔
509	options % whichNorm = optType	1✔
510	!
511	! --------
512	! Geometry
513	! --------
514	!
515	ALLOCATE(line)	1✔
516	CALL line % initWithStartEndNameAndID(xStart = [-1.0_RP,0.0_RP,0.0_RP], &
517	xEnd = [1.0_RP,0.0_RP,0.0_RP], &
518	cName = "line", id = 1)	1✔
519	ALLOCATE(circle)	1✔
520	CALL circle % initWithParametersNameAndID(center = [0.0_RP, 0.0_RP, 0.0_RP], &
521	radius = 1.0_RP, &
522	startAngle = 0.0_RP, &
523	endAngle = PI, &
524	cName = "circle", &
525	id = 2)	1✔
526	ALLOCATE(chain)	1✔
527	CALL chain % initChainWithNameAndID(chainName = "chain",id = 3)	1✔
528	!
529	curvePtr => circle	1✔
530	CALL chain % addCurve(curvePtr)	1✔
531	obj => circle	1✔
532	CALL release(obj)	1✔
533
534	curvePtr => line	1✔
535	CALL chain % addCurve(curvePtr)	1✔
536	obj => line	1✔
537	CALL release(obj)	1✔
538	!
539	CALL chain % complete(innerOrOuterCurve = OUTER,chainMustClose = .TRUE.)	1✔
540
541	IF ( catch() ) THEN	1✔
NEW 542	WRITE(stderr,*)	×
NEW 543	WRITE(stderr,*) "------------------------------------------------------------------"	×
NEW 544	WRITE(stderr,*)	×
NEW 545	WRITE(stderr,*) "The following errors were found when constructing the check:"	×
546
NEW 547	DO	×
NEW 548	exception => popLastException()	×
NEW 549	IF ( .NOT.ASSOCIATED(exception) ) EXIT	×
NEW 550	CALL exception % printDescription(stderr)	×
NEW 551	errorSeverity = MAX(errorSeverity, exception % severity())	×
552	END DO
NEW 553	WRITE(stderr,*)	×
NEW 554	WRITE(stderr,*) "------------------------------------------------------------------"	×
NEW 555	WRITE(stderr,*)	×
556
NEW 557	IF ( errorSeverity > FT_ERROR_WARNING ) THEN	×
NEW 558	ERROR STOP "The Errors were Fatal. Cannot generate optimized curve."	×
559	END IF
560	END IF
561	!
562	! -------------
563	! Optimizations
564	! -------------
565	!
566	curvePtr => chain	1✔
567	CALL ConstructGaussQuadrature(gQuad, 4*polyOrder) ! For error computation. The 4 is arbitrary	1✔
568	CALL FindOptimizedCuts(curvePtr, polyOrder, testBreaks, options, gQuad, optimizedCuts, breakIndices)	1✔
569	!
570	! ------------
571	! Check errors
572	! ------------
573	!
574	CALL FTAssert(MAXVAL(ABS(optimizedCuts - expectedCuts)) .le. testTol, msg = "Segments dont match")	9✔
575	CALL FTAssert(MAXVAL(ABS(expectedIndices - breakIndices)) == 0 , msg = "Break Indices dont match")	3✔
576
577	CALL releaseChainedCurve(chain)	1✔
578
579	END SUBROUTINE SegmentedCurveCheck	1✔
580

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