Generalize optimality criteria.

ElmerCSC · Dec 12, 2024 · d3d9e32 · d3d9e32
1 parent 007ac00
commit d3d9e32
Showing 1 changed file with 7 additions and 4 deletions.
diff --git a/fem/src/modules/TopoOpt.F90 b/fem/src/modules/TopoOpt.F90
@@ -503,7 +503,7 @@ SUBROUTINE UpdateDensities(x,dc,dv,g)
     REAL(KIND=dp) :: g
 
     REAL(KIND=dp), ALLOCATABLE :: xnew(:)
-    REAL(KIND=dp) :: Vi, l1, l2, lmid, move, err, tol, V0
+    REAL(KIND=dp) :: Vi, l1, l2, lmid, move, err, tol, V0, damp
     INTEGER :: k
     LOGICAL :: Visited = .FALSE.
 
@@ -518,6 +518,9 @@ SUBROUTINE UpdateDensities(x,dc,dv,g)
 
     tol = ListGetCReal(Params,'Bisection search tolerance',Found )
     IF(.NOT. Found) tol = 1.0e-6
+
+    damp = ListGetCReal(Params,'Bisection search damping exponent',Found )
+    IF(.NOT. Found) damp = 0.5_dp
 
     ! Desired total volume
     V0 = volFrac * SUM(dv)
@@ -531,7 +534,7 @@ SUBROUTINE UpdateDensities(x,dc,dv,g)
 
       ! Note: xnew in [0,1]
       ! Suggested new density
-      xnew = x*SQRT(-dc/(dv*lmid))
+      xnew = x*(MAX(1.0e-10,-dc/(dv*lmid)))**damp
 
       ! Regulators and limiters
       xnew = MAX(0.0_dp,MAX(x-move,MIN(1.0_dp,MIN(x+move,xnew))))
@@ -547,11 +550,11 @@ SUBROUTINE UpdateDensities(x,dc,dv,g)
       END IF
 
       err = (l2-l1)/(l1+l2)
-      IF( err < tol ) EXIT      
-
       IF( InfoActive(20)) THEN
         PRINT *,'Bisection:',k,Vi,l1,l2,err
       END IF
+
+      IF( err < tol ) EXIT      
     END DO
 
     x = xnew