From a5feb0b81ce17f8e333b36486cbac20d86000365 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Mon, 21 Oct 2024 11:33:29 +1000 Subject: [PATCH] Add some testing cases for 3D to MWE --- scripts/Embedded/MWEs/jordi_embedded.jl | 31 +++++++++++++++---------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/scripts/Embedded/MWEs/jordi_embedded.jl b/scripts/Embedded/MWEs/jordi_embedded.jl index b4205af..d929a70 100644 --- a/scripts/Embedded/MWEs/jordi_embedded.jl +++ b/scripts/Embedded/MWEs/jordi_embedded.jl @@ -9,11 +9,11 @@ import Gridap.Geometry: get_node_coordinates, collect1d include("../differentiable_trians.jl") order = 1 -n = 40 +n = 15 N = 8 -_model = CartesianDiscreteModel((0,1,0,1),(n,n)) -#_model = CartesianDiscreteModel((0,1,0,1,0,1),(n,n,n)) +# _model = CartesianDiscreteModel((0,1,0,1),(n,n)) +_model = CartesianDiscreteModel((0,1,0,1,0,1),(n,n,n)) cd = Gridap.Geometry.get_cartesian_descriptor(_model) base_model = UnstructuredDiscreteModel(_model) ref_model = refine(base_model, refinement_method = "barycentric") @@ -31,9 +31,16 @@ V_φ = TestFESpace(model,reffe_scalar) # φh = interpolate(x->(1+0.25)abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25,V_φ) # Straight lines with scaling # φh = interpolate(x->abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25/(1+0.25),V_φ) # Straight lines without scaling # φh = interpolate(x->tan(-pi/4)*(x[1]-0.5)+(x[2]-0.5),V_φ) # Angled interface +# φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.5223,V_φ) # Circle -φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.5223,V_φ) # Circle -# φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.25,V_φ) # Sphere +# φh = interpolate(x->max(abs(x[1]-0.5),abs(x[2]-0.5),abs(x[3]-0.5))-0.25,V_φ) # Square prism +# φh = interpolate(x->((x[1]-0.5)^N+(x[2]-0.5)^N+(x[3]-0.5)^N)^(1/N)-0.25,V_φ) # Curved corner example +# φh = interpolate(x->abs(x[1]-0.5)+abs(x[2]-0.5)+abs(x[3]-0.5)-0.25-0/n/10,V_φ) # Diamond prism +# φh = interpolate(x->(1+0.25)abs(x[1]-0.5)+0abs(x[2]-0.5)+0abs(x[3]-0.5)-0.25,V_φ) # Straight lines with scaling +# φh = interpolate(x->abs(x[1]-0.5)+0abs(x[2]-0.5)+0abs(x[3]-0.5)-0.25/(1+0.25),V_φ) # Straight lines without scaling +# φh = interpolate(x->tan(-pi/3)*(x[1]-0.5)+(x[2]-0.5),V_φ) # Angled interface +# φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.53,V_φ) # Sphere +φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])*cos(2π*x[3])-0.11,V_φ) # "Regular" LSF fh = interpolate(x->cos(x[1]*x[2]),V_φ) @@ -54,13 +61,13 @@ cutgeo = cut(model,geo) diff_Ωin = DifferentiableTriangulation(Ωin) diff_Ωout = DifferentiableTriangulation(Ωout) -dΩin = Measure(diff_Ωin,3*order) +dΩin = Measure(diff_Ωin,2*order) j_in(φ) = ∫(fh)dΩin dj_in = gradient(j_in,φh) dj_vec_in = assemble_vector(dj_in,V_φ) norm(dj_vec_in) -dΩout = Measure(diff_Ωout,3*order) +dΩout = Measure(diff_Ωout,2*order) j_out(φ) = ∫(fh)dΩout dj_out = gradient(j_out,φh) dj_vec_out = -assemble_vector(dj_out,V_φ) @@ -72,8 +79,8 @@ dj_expected(q) = ∫(-fh*q/(norm ∘ (∇(φh))))dΓ dj_exp_vec = assemble_vector(dj_expected,V_φ) norm(dj_exp_vec) -norm(dj_vec_in-dj_exp_vec) -norm(dj_vec_out-dj_exp_vec) +println("dj1 error = $(norm(dj_vec_in-dj_exp_vec))") +println("dj1 error = $(norm(dj_vec_out-dj_exp_vec))") ############################################################################################ @@ -82,7 +89,7 @@ include("../SubFacetSkeletons.jl") Γ = EmbeddedBoundary(cutgeo) Λ = Skeleton(Γ) -dΓ = Measure(Γ,3*order) +dΓ = Measure(Γ,2*order) dΛ = Measure(Λ,2*order) Γpts = get_cell_points(Γ) @@ -106,7 +113,7 @@ J2(φ) = ∫(fh)dΓ_AD dJ2_AD = gradient(J2,φh) dj2_AD = assemble_vector(dJ2_AD,V_φ) # TODO: Why is this +ve but AD on volume is -ve??? -norm(dj2 - dj2_AD) +println("dj2 uncorrected error = $(norm(dj2 - dj2_AD))") ############################################################################################ @@ -121,7 +128,7 @@ dΣ = Measure(Σ,2*order) dJ3(w) = dJ2(w) + ∫(-n_S_Σ ⋅ (fh*m_k_Σ * w / ∇ˢφ_Σ))dΣ dj3 = assemble_vector(dJ3,V_φ) -norm(dj3 - dj2_AD) +println("dj2 error = $(norm(dj3 - dj2_AD))") ############################################################################################