diff --git a/src/norm.jl b/src/norm.jl index aedd1b8..751910b 100644 --- a/src/norm.jl +++ b/src/norm.jl @@ -427,12 +427,21 @@ function penalty(root::Cell{Data, Dim, T, L}, xc, xc_init, dist_ref, mu, degree for d = 1:Dim dist += (xc[d,i] - xc_init[d,i])^2 end - phi += dist/(dist_ref[i]^2) + phi += 0.5*mu*dist/(dist_ref[i]^2) end # compute the diagonal norm based on x H = diagonal_norm(root, xc, degree) + # # compute the discrete KS function + # rho = 100.0 + # minH = minimum(real.(H)) + # sum_exp = 0.0 + # for i = 1:num_nodes + # sum_exp += exp(rho*(minH - H[i])) + # end + # phi += -minH + log(sum_exp/num_nodes)/rho + # add the penalties tol = 1e-5 for i = 1:num_nodes @@ -443,7 +452,6 @@ function penalty(root::Cell{Data, Dim, T, L}, xc, xc_init, dist_ref, mu, degree end end - phi *= 0.5 return phi end @@ -454,11 +462,31 @@ function penalty_grad!(g::AbstractVector{T}, root::Cell{Data, Dim, T, L}, # need the diagonal norm for the reverse sweep H = diagonal_norm(root, xc, degree) + # rho = 100.0 + # minH = minimum(real.(H)) + # sum_exp = 0.0 + # for i = 1:num_nodes + # sum_exp += exp(rho*(minH - H[i])) + # end + # start the reverse sweep fill!(g, zero(T)) + # return phi # phi *= 0.5 - phi_bar = 0.5 + # phi_bar = 0.5 + + # return phi + phi_bar = 1.0 + + # # phi += -minH + log(sum_exp/num_nodes)/rho + # sum_exp_bar = phi_bar/(sum_exp * rho) + # H_bar = zero(H) + # for i = 1:num_nodes + # # sum_exp += exp(rho*(minH - H[i])) + # H_bar[i] -= sum_exp_bar*exp(rho*(minH - H[i]))*rho + # end + # add the penalties tol = 1e-5 H_bar = zero(H) @@ -478,8 +506,8 @@ function penalty_grad!(g::AbstractVector{T}, root::Cell{Data, Dim, T, L}, # compute the norm part of the penalty for i = 1:num_nodes - # phi += dist/(dist_ref[i]^2) - dist_bar = phi_bar/(dist_ref[i]^2) + # phi += 0.5*mu*dist/(dist_ref[i]^2) + dist_bar = 0.5*mu*phi_bar/(dist_ref[i]^2) for d = 1:Dim # dist += (xc[d,i] - xc_init[d,i])^2 xc_bar[d,i] += dist_bar*2.0*(xc[d,i] - xc_init[d,i]) diff --git a/test/opt-norm.jl b/test/opt-norm.jl index d09ae7f..28125ee 100644 --- a/test/opt-norm.jl +++ b/test/opt-norm.jl @@ -14,7 +14,7 @@ using SparseArrays #Random.seed!(42) Dim = 2 -degree = 5 +degree = 2 # use a unit HyperRectangle root = Cell(SVector(ntuple(i -> 0.0, Dim)), @@ -24,7 +24,7 @@ root = Cell(SVector(ntuple(i -> 0.0, Dim)), # DGD dof locations num_basis = binomial(Dim + degree, Dim) -num_nodes = 10*num_basis +num_nodes = 100*num_basis points = rand(Dim, num_nodes) # refine mesh, build sentencil, and evaluate particular quad and nullspace @@ -46,7 +46,7 @@ CutDGD.build_nn_stencils!(root, points, degree) #end #println(dist_ref) dist_ref = ones(num_nodes) -mu = 100.0 +mu = 0.0001 g = zeros(num_nodes*Dim) @@ -71,7 +71,7 @@ for n = 1:max_iter obj0 = obj dxc = reshape(g, (Dim, num_nodes)) - alpha = 1.0 + alpha = 10.0 for k = 1:max_line global points -= alpha*dxc obj = CutDGD.penalty(root, points, points_init, dist_ref, mu, degree)