Got rid of Connor workaround

scripts/Embedded/MWEs/jordi_embedded.jl

@@ -8,32 +8,45 @@ import Gridap.Geometry: get_node_coordinates, collect1d
 
 include("../differentiable_trians.jl")
 
+order = 1
 model = CartesianDiscreteModel((0,1,0,1),(100,100))
 Ω = Triangulation(model)
-dΩ = Measure(Ω,2)
+dΩ = Measure(Ω,2*order)
 
-reffe_scalar = ReferenceFE(lagrangian,Float64,1)
+reffe_scalar = ReferenceFE(lagrangian,Float64,order)
 V = TestFESpace(model,reffe_scalar)
 U = TrialFESpace(V)
 V_φ = TestFESpace(model,reffe_scalar)
 
-φh = interpolate(x->-cos(4π*x[1])*cos(4*pi*x[2])/4-0.2/4,V_φ)
+φh = interpolate(x->max(abs(x[1]-0.5),abs(x[2]-0.5))-0.1,V_φ)
+# φh = interpolate(x->max(abs(x[1]-0.5),abs(x[2]-0.5))-0.11,V_φ)
+# φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.3,V_φ)
 
 Ωin = DifferentiableTriangulation(φh) do φh
   geo = DiscreteGeometry(φh,model)
   cutgeo = cut(model,geo)
   return Triangulation(cutgeo,PHYSICAL_IN)
 end
 
-dΩin = Measure(Ωin,2)
+dΩin = Measure(Ωin,2*order)
+j(φ) = ∫(1)dΩin
+dj = gradient(j,φh)
+dj_vec = assemble_vector(dj,V_φ)
 
-j(φ) = ∫(φ)dΩin
-gradient(j,φh)
+geo = DiscreteGeometry(φh,model)
+cutgeo = cut(model,geo)
+Γ = EmbeddedBoundary(cutgeo)
+dΓ = Measure(Γ,2*order)
+dj_expected(q) = ∫(-q)dΓ
+dj_exp_vec = assemble_vector(dj_expected,V_φ)
+
+norm(dj_vec-dj_exp_vec)
+
+writevtk(Ω,"results/test",cellfields=["φh"=>φh,"dj"=>FEFunction(V_φ,dj_vec),"dj_expected"=>FEFunction(V_φ,dj_exp_vec)])
 
 # To-Do: 
 #   1. Add updateability condition
 #   2. Add caches so we don't have to recompute everything every time
 #   3. Figure out a way to share the cache between different triangulations created from the 
 #      same cut geometry, i.e PHYS_IN/PHYS_OUT/Boundary
 #   4. Anything else?
-

scripts/Embedded/differentiable_trians.jl

@@ -1,5 +1,5 @@
 
-using Gridap.CellData, Gridap.Geometry, Gridap.Helpers
+using Gridap.CellData, Gridap.Geometry, Gridap.Helpers, Gridap.Arrays
 
 function CellData.get_contribution(a::DomainContribution,trian::Geometry.AppendedTriangulation)
   if haskey(a.dict,trian)
@@ -15,6 +15,11 @@ function CellData.get_contribution(a::DomainContribution,trian::Geometry.Appende
   end
 end
 
+function FESpaces._compute_cell_ids(uh,ttrian::AppendedTriangulation)
+  ids_a = FESpaces._compute_cell_ids(uh,ttrian.a)
+  ids_b = FESpaces._compute_cell_ids(uh,ttrian.b)
+  lazy_append(ids_a,ids_b)
+end
 
 # DifferentiableTriangulation
 
@@ -38,7 +43,7 @@ function update_trian!(trian::DifferentiableTriangulation,φh)
   return trian
 end
 
-function Gridap.FESpaces._change_argument(
+function FESpaces._change_argument(
   op,f,trian::DifferentiableTriangulation,uh::SingleFieldFEFunction
 )
   U = get_fe_space(uh)
@@ -51,16 +56,27 @@ function Gridap.FESpaces._change_argument(
   g
 end
 
+function FESpaces._compute_cell_ids(uh,ttrian::DifferentiableTriangulation)
+  FESpaces._compute_cell_ids(uh,ttrian.state)
+end
+
+function Geometry.get_background_model(t::DifferentiableTriangulation)
+  get_background_model(t.state)
+end
+
+Geometry.get_glue(ttrian::DifferentiableTriangulation,val::Val{d}) where d = get_glue(ttrian.state,val)
+
 # Mutable measure
 
 mutable struct MutableMeasure <: Measure
   state :: Measure
-  params 
+  trian :: Triangulation
+  params
 end
 
 function Measure(trian::DifferentiableTriangulation,args...;kwargs...)
   state = Measure(trian.state,args...;kwargs...)
-  meas  = MutableMeasure(state,(args,kwargs))
+  meas  = MutableMeasure(state, trian,(args,kwargs))
   push!(trian.children, objectid(meas) => meas)
   return meas
 end
@@ -71,4 +87,9 @@ function update_measure!(meas::MutableMeasure,trian::Triangulation)
   return meas
 end
 
-CellData.integrate(f,b::MutableMeasure) = integrate(f,b.state)
+function CellData.integrate(f,b::MutableMeasure)
+  c = integrate(f,b.state.quad)
+  cont = DomainContribution()
+  add_contribution!(cont,b.trian,c)
+  cont
+end