From 16afe46b435c1d003cfcd0be7e6384be1bdc917e Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Thu, 11 Jul 2024 15:00:36 +1000
Subject: [PATCH] testing

---
 scripts/Embedded/MWEs/zach_embedded.jl        |  53 +++++++++
 .../differentiable_trians_zachs_ver.jl        | 110 ++++++++++++++++++
 2 files changed, 163 insertions(+)
 create mode 100644 scripts/Embedded/MWEs/zach_embedded.jl
 create mode 100644 scripts/Embedded/differentiable_trians_zachs_ver.jl

diff --git a/scripts/Embedded/MWEs/zach_embedded.jl b/scripts/Embedded/MWEs/zach_embedded.jl
new file mode 100644
index 00000000..e496cc74
--- /dev/null
+++ b/scripts/Embedded/MWEs/zach_embedded.jl
@@ -0,0 +1,53 @@
+
+using Gridap
+
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
+import Gridap.Geometry: get_node_coordinates, collect1d
+
+include("../differentiable_trians.jl")
+
+order = 1
+model = CartesianDiscreteModel((0,1,0,1),(30,30))
+Ω = Triangulation(model)
+dΩ = Measure(Ω,2order)
+
+reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+V = TestFESpace(model,reffe_scalar)
+U = TrialFESpace(V)
+V_φ = TestFESpace(model,reffe_scalar)
+
+φ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,2order)
+
+j(φ) = ∫(1)dΩin
+dj = gradient(j,φh)
+dj_vec = assemble_vector(dj,V_φ)
+
+geo = DiscreteGeometry(φh,model)
+cutgeo = cut(model,geo)
+Γ = EmbeddedBoundary(cutgeo)
+dΓ = Measure(Γ,2order)
+dj_expected(q) = ∫(-q)dΓ
+dj_exp_vec = assemble_vector(dj_expected,V_φ)
+
+writevtk(Ω,"results/test",cellfields=["φh"=>φh,"dj"=>FEFunction(V_φ,dj_vec),"dj_expected"=>FEFunction(V_φ,dj_exp_vec)])
+# writevtk(Ωin.recipe(φh),"results/test_phys",cellfields=["φh"=>φh,"dj"=>FEFunction(V_φ,dj_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?
+
diff --git a/scripts/Embedded/differentiable_trians_zachs_ver.jl b/scripts/Embedded/differentiable_trians_zachs_ver.jl
new file mode 100644
index 00000000..a51bb8d4
--- /dev/null
+++ b/scripts/Embedded/differentiable_trians_zachs_ver.jl
@@ -0,0 +1,110 @@
+
+using Gridap.CellData, Gridap.Geometry, Gridap.Helpers
+
+function CellData.get_contribution(a::DomainContribution,trian::Geometry.AppendedTriangulation)
+  if haskey(a.dict,trian)
+    return a.dict[trian]
+  elseif haskey(a.dict,trian.a)
+    return a.dict[trian.a]
+  elseif haskey(a.dict,trian.b)
+    return a.dict[trian.b]
+  else
+    @unreachable """\n
+    There is not contribution associated with the given mesh in this DomainContribution object.
+    """
+  end
+end
+
+
+# DifferentiableTriangulation
+
+mutable struct DifferentiableTriangulation{Dc,Dp} <: Triangulation{Dc,Dp}
+  recipe   :: Function
+  state    :: Triangulation{Dc,Dp}
+  children :: IdDict{UInt64,Measure}
+end
+
+function DifferentiableTriangulation(f::Function,φh)
+  DifferentiableTriangulation(f,f(φh),IdDict{UInt64,Measure}())
+end
+
+(t::DifferentiableTriangulation)(φh) = update_trian!(t,φh)
+
+function update_trian!(trian::DifferentiableTriangulation,φh)
+  trian.state = trian.recipe(φh)
+  for child in values(trian.children)
+    update_measure!(child,trian.state)
+  end
+  return trian
+end
+
+function Gridap.FESpaces._change_argument(
+  op,f,trian::DifferentiableTriangulation,uh::SingleFieldFEFunction
+)
+  U = get_fe_space(uh)
+  function g(cell_u)
+    cf = CellField(U,cell_u)
+    update_trian!(trian,cf)
+    cell_grad = f(cf)
+    get_contribution(cell_grad,trian)
+  end
+  g
+end
+
+function GridapEmbedded.Interfaces.get_background_model(t::DifferentiableTriangulation)
+  get_background_model(t.state)
+end
+
+Gridap.Geometry.get_glue(ttrian::DifferentiableTriangulation,val::Val{d}) where d = get_glue(ttrian.state,val)
+
+# Mutable measure
+
+mutable struct MutableMeasure <: Measure
+  state :: Measure
+  trian :: Triangulation
+  params
+end
+
+function Measure(trian::DifferentiableTriangulation,args...;kwargs...)
+  state = Measure(trian.state,args...;kwargs...)
+  meas  = MutableMeasure(state, trian,(args,kwargs))
+  push!(trian.children, objectid(meas) => meas)
+  return meas
+end
+
+function update_measure!(meas::MutableMeasure,trian::Triangulation)
+  args, kwargs = meas.params
+  meas.state = Measure(trian,args...;kwargs...)
+  return meas
+end
+
+function CellData.integrate(f,b::MutableMeasure)
+  c = integrate(f,b.state.quad)
+  cont = DomainContribution()
+  add_contribution!(cont,b.trian,c)
+  cont
+end
+
+## Fix for _array_cache used by Connor
+using Gridap.Arrays: IndexItemPair,LazyArray,testitem,return_cache,return_value
+function Gridap.Arrays._array_cache!(dict::Dict,a::LazyArray)
+  @boundscheck begin
+    if ! all(map(isconcretetype, map(eltype, a.args)))
+      for n in 1:length(a.args)
+        @notimplementedif ! all(map(isconcretetype, map(eltype, a.args[n])))
+      end
+    end
+    if ! (eltype(a.maps) <: Function)
+      @notimplementedif ! isconcretetype(eltype(a.maps))
+    end
+  end
+  gi = testitem(a.maps)
+  fi = map(testitem,a.args)
+  cg = array_cache(dict,a.maps)
+  cf = map(fi->array_cache(dict,fi),a.args)
+  cgi = return_cache(gi, fi...)
+  index = -1
+  #item = evaluate!(cgi,gi,testargs(gi,fi...)...)
+  item = return_value(gi,fi...)
+  (cg, cgi, cf), IndexItemPair(index, item)
+end
\ No newline at end of file