Add some testing cases for 3D to MWE

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))")
 
 ############################################################################################