diff --git a/scripts/Embedded/Examples/2d_elastic_compliance.jl b/scripts/Embedded/Examples/2d_elastic_compliance.jl index bb5e299f..3499fb96 100644 --- a/scripts/Embedded/Examples/2d_elastic_compliance.jl +++ b/scripts/Embedded/Examples/2d_elastic_compliance.jl @@ -1,7 +1,7 @@ using Pkg; Pkg.activate() using Gridap,GridapTopOpt -include("../embedded_measures.jl") +include("../embedded_measures_AD_DISABLED.jl") function main(path="./results/UnfittedFEM_elastic_compliance_ALM/") ## Parameters diff --git a/scripts/Embedded/Examples/2d_inverse_homenisation.jl b/scripts/Embedded/Examples/2d_inverse_homenisation.jl index f5a0aff6..c6e2d328 100644 --- a/scripts/Embedded/Examples/2d_inverse_homenisation.jl +++ b/scripts/Embedded/Examples/2d_inverse_homenisation.jl @@ -1,7 +1,7 @@ using Pkg; Pkg.activate() using Gridap,GridapTopOpt -include("../embedded_measures.jl") +include("../embedded_measures_AD_DISABLED.jl") function main(path="./results/UnfittedFEM_inverse_homenisation_ALM/") ## Parameters diff --git a/scripts/Embedded/Examples/2d_thermal.jl b/scripts/Embedded/Examples/2d_thermal.jl index 086d14fa..44c2efa1 100644 --- a/scripts/Embedded/Examples/2d_thermal.jl +++ b/scripts/Embedded/Examples/2d_thermal.jl @@ -1,7 +1,7 @@ using Pkg; Pkg.activate() using Gridap,GridapTopOpt -include("../embedded_measures.jl") +include("../embedded_measures_AD_DISABLED.jl") function main() path="./results/UnfittedFEM_thermal_compliance_ALM/" diff --git a/scripts/Embedded/Examples/2d_thermal_optimised.jl b/scripts/Embedded/Examples/2d_thermal_optimised.jl index 8a6536a2..4c5f5f3d 100644 --- a/scripts/Embedded/Examples/2d_thermal_optimised.jl +++ b/scripts/Embedded/Examples/2d_thermal_optimised.jl @@ -1,7 +1,7 @@ using Pkg; Pkg.activate() using Gridap,GridapTopOpt,GridapEmbedded -include("../embedded_measures.jl") +include("../embedded_measures_AD_DISABLED.jl") function main() path="./results/UnfittedFEM_thermal_compliance_ALM_Optimised_test/" diff --git a/scripts/Embedded/Examples/2d_thermal_optimised_MPI.jl b/scripts/Embedded/Examples/2d_thermal_optimised_MPI.jl index 7939b04d..71854bee 100644 --- a/scripts/Embedded/Examples/2d_thermal_optimised_MPI.jl +++ b/scripts/Embedded/Examples/2d_thermal_optimised_MPI.jl @@ -1,7 +1,7 @@ using Gridap,GridapTopOpt,GridapEmbedded using GridapDistributed, GridapPETSc, PartitionedArrays -include("../embedded_measures.jl") +include("../embedded_measures_AD_DISABLED.jl") function main(parts,distribute) ranks = distribute(LinearIndices((prod(parts),))) diff --git a/scripts/Embedded/MWEs/gridap_embedded_grad_testing_MWE.jl b/scripts/Embedded/MWEs/gridap_embedded_grad_testing_MWE.jl new file mode 100644 index 00000000..83d33e42 --- /dev/null +++ b/scripts/Embedded/MWEs/gridap_embedded_grad_testing_MWE.jl @@ -0,0 +1,30 @@ +using Pkg; Pkg.activate() + +using Gridap + +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData +import Gridap.Geometry: get_node_coordinates, collect1d + +include("../embedded_measures.jl") + +model = CartesianDiscreteModel((0,1,0,1),(100,100)); +Ω = Triangulation(model) +dΩ = Measure(Ω,2) + +reffe_scalar = ReferenceFE(lagrangian,Float64,1) +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_φ) +embedded_meas = EmbeddedMeasureCache(φh,V_φ) +update_meas(args...) = update_embedded_measures!(embedded_meas,args...) +get_meas(args...) = get_embedded_measures(embedded_meas,args...) + +_j(φ,dΩ,dΩ1,dΩ2,dΓ) = ∫(φ)dΩ1 + ∫(φ)dΩ2 + ∫(φ)dΓ + +j_iem = IntegrandWithEmbeddedMeasure(_j,(dΩ,),update_meas) +j_iem(φh) +gradient(j_iem,φh) \ No newline at end of file diff --git a/scripts/Embedded/MWEs/gridap_embedded_grad_testing_dist_MWE.jl b/scripts/Embedded/MWEs/gridap_embedded_grad_testing_dist_MWE.jl new file mode 100644 index 00000000..4115145a --- /dev/null +++ b/scripts/Embedded/MWEs/gridap_embedded_grad_testing_dist_MWE.jl @@ -0,0 +1,37 @@ +using Pkg; Pkg.activate() + +using Gridap + +using GridapDistributed, GridapPETSc, PartitionedArrays + +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData +import Gridap.Geometry: get_node_coordinates, collect1d + +include("../embedded_measures.jl") + +parts = (2,2); +ranks = with_debug() do distribute + distribute(LinearIndices((prod(parts),))) +end + +model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(100,100)); +Ω = Triangulation(model) +dΩ = Measure(Ω,2) + +reffe_scalar = ReferenceFE(lagrangian,Float64,1) +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_φ) +embedded_meas = EmbeddedMeasureCache(φh,V_φ) +update_meas(args...) = update_embedded_measures!(embedded_meas,args...) +get_meas(args...) = get_embedded_measures(embedded_meas,args...) + +_j(φ,dΩ,dΩ1,dΩ2,dΓ) = ∫(φ)dΩ1 + ∫(φ)dΩ2 + ∫(φ)dΓ + +j_iem = IntegrandWithEmbeddedMeasure(_j,(dΩ,),update_meas) +j_iem(φh) +gradient(j_iem,φh) \ No newline at end of file diff --git a/scripts/Embedded/MWEs/gridap_embedded_testing_MWE.jl b/scripts/Embedded/MWEs/gridap_embedded_testing_MWE.jl deleted file mode 100644 index 29a0cfc8..00000000 --- a/scripts/Embedded/MWEs/gridap_embedded_testing_MWE.jl +++ /dev/null @@ -1,106 +0,0 @@ -using Pkg; Pkg.activate() - -using Gridap - -using GridapEmbedded -using GridapEmbedded.LevelSetCutters -using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData -import Gridap.Geometry: get_node_coordinates, collect1d - -function cv_to_dof(cv,V) - fv=zeros(eltype(eltype(cv)),num_free_dofs(V)) - gather_free_values!(fv,V,cv) -end - -function field_to_cv(uh::FEFunction) - get_cell_dof_values(uh) -end - -function field_to_cv(cf::CellField) - cv=cf.cell_field.args[1] -end - -function get_geo_params(ϕₕ::FEFunction,Vbg) - Ωbg = get_triangulation(Vbg) - bgmodel = get_background_model(Ωbg) - point_to_coords = collect1d(get_node_coordinates(bgmodel)) - ls_to_point_to_value_unmasked = field_to_cv(ϕₕ) - p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg) - geo1 = DiscreteGeometry(p0,point_to_coords) - geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="") - get_geo_params(geo1,geo2,bgmodel) -end - -function get_geo_params(ϕₕ::CellField,Vbg) - Ωbg = get_triangulation(Vbg) - bgmodel = get_background_model(Ωbg) - point_to_coords = collect1d(get_node_coordinates(bgmodel)) - ls_to_point_to_value_unmasked = field_to_cv(ϕₕ) - p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg) - geo1 = DiscreteGeometry(p0,point_to_coords) - geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="") - get_geo_params(geo1,geo2,bgmodel) -end - -function get_geo_params(ϕ::AbstractVector,Vbg) - Ωbg = get_triangulation(Vbg) - bgmodel = get_background_model(Ωbg) - point_to_coords = collect1d(get_node_coordinates(bgmodel)) - geo1 = DiscreteGeometry(ϕ,point_to_coords,name="") - geo2 = DiscreteGeometry(-ϕ,point_to_coords,name="") - get_geo_params(geo1,geo2,bgmodel) -end - -function get_geo_params(geo1,geo2,bgmodel) - cutgeo1= cut(bgmodel,geo1) - cutgeo2= cut(bgmodel,geo2) - # Setup interpolation meshes - Ω1_act = Triangulation(cutgeo1,ACTIVE) - Ω2_act = Triangulation(cutgeo2,ACTIVE) - # Setup integration meshes - Ω1 = Triangulation(cutgeo1,PHYSICAL) - Ω2 = Triangulation(cutgeo2,PHYSICAL) - Γ = EmbeddedBoundary(cutgeo1) - # Setup Lebesgue measures - order = 1 - degree = 2*order - dΩ1 = Measure(Ω1,degree) - dΩ2 = Measure(Ω2,degree) - dΓ = Measure(Γ,degree) - (;dΩ1,dΩ2,dΓ)#,debug=(;cutgeo1,cutgeo2,order, Ω1,Ω2,Ω1_act,Ω2_act,Γ)) -end - -model = CartesianDiscreteModel((0,1,0,1),(100,100)); -Ω = Triangulation(model) -dΩ = Measure(Ω,2) -Γ_N = BoundaryTriangulation(model,tags=6) -dΓ_N = Measure(Γ_N,2) - -reffe_scalar = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,1) # ReferenceFE(lagrangian,Float64,1) -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_φ) - -# _Γ = get_triangulation(Gridap.CellData.get_cell_quadrature(get_geo_params(φh,V_φ).dΓ)) -# writevtk(_Γ,"./results/Gamma",cellfields=["normal"=>get_normal_vector(_Γ)]) - -function a(u,v,φ) - geo = get_geo_params(φ,V_φ) - dΩ1 = geo.dΩ1; dΩ2 = geo.dΩ2 - ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-3*∇(u)⋅∇(v))dΩ2 -end - -l(v,φ) = ∫(v)dΓ_N - -op = AffineFEOperator((u,v) -> a(u,v,φh),v -> l(v,φh),U,V) - -function j(u,φ) - geo = get_geo_params(φ,V_φ) - dΩ1 = geo.dΩ1; dΩ2 = geo.dΩ2 - ∫(u)dΩ1 + ∫(u)dΩ2 -end - -∇(φ -> a(zero(U),zero(U),φ))(φh) # error calling _gradient_nd! -∇(φ -> j(zero(U),φ))(φh) # domain contribution not found error \ No newline at end of file diff --git a/scripts/Embedded/MWEs/gridap_embedded_testing_MWE_periodic.jl b/scripts/Embedded/MWEs/gridap_embedded_testing_MWE_periodic.jl deleted file mode 100644 index 85cf9b17..00000000 --- a/scripts/Embedded/MWEs/gridap_embedded_testing_MWE_periodic.jl +++ /dev/null @@ -1,34 +0,0 @@ -using Pkg; Pkg.activate() - -using Gridap -using GridapEmbedded,GridapEmbedded.LevelSetCutters -import Gridap.Geometry: get_node_coordinates,collect1d - -model = CartesianDiscreteModel((0,1,0,1),(100,100),isperiodic=(true,true)); -V_φ = TestFESpace(model,ReferenceFE(lagrangian,Float64,1)) - -## Discrete -φh = interpolate(x->-cos(4π*(x[1]))*cos(4*pi*x[2])/4-0.2/4,V_φ) -point_to_coords = collect1d(get_node_coordinates(model)) -geo = DiscreteGeometry(φh(point_to_coords),point_to_coords,name="") -cutgeo = cut(model,geo) - -Ω = Triangulation(cutgeo,PHYSICAL) -Ω_act = Triangulation(cutgeo,ACTIVE) - -# writevtk(Ω,"./results/discrete_geo_periodic",cellfields=["φh"=>φh]) - -V_φ_test = TestFESpace(Ω_act,ReferenceFE(lagrangian,Float64,1)) -φh_test = interpolate(x->sqrt(x[1]^2+x[2]^2)-0.5,V_φ_test) -writevtk(Ω,"./results/discrete_geo_periodic",cellfields=["φh"=>φh_test]) - -## Analytic -analytic_geo=AnalyticalGeometry(x->-cos(4π*(x[1]))*cos(4*pi*x[2])/4-0.2/4) -cutgeo_analytic = cut(model,analytic_geo) -Ω_analytic = Triangulation(cutgeo_analytic,PHYSICAL) -Ω_analytic_act = Triangulation(cutgeo_analytic,ACTIVE) -# writevtk(Ω_analytic,"./results/analytic_geo_periodic",cellfields=["φh"=>φh]) - -V_φ_test = TestFESpace(Ω_analytic_act,ReferenceFE(lagrangian,Float64,1)) -φh_test = interpolate(x->sqrt(x[1]^2+x[2]^2)-0.5,V_φ_test) -writevtk(Ω_analytic,"./results/analytic_geo_periodic",cellfields=["φh"=>φh_test]) \ No newline at end of file diff --git a/scripts/Embedded/embedded_measures.jl b/scripts/Embedded/embedded_measures.jl index 96e514b5..f2d0f904 100644 --- a/scripts/Embedded/embedded_measures.jl +++ b/scripts/Embedded/embedded_measures.jl @@ -7,74 +7,20 @@ using GridapEmbedded.Distributed: DistributedDiscreteGeometry using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Helpers import Gridap.Geometry: get_node_coordinates, collect1d -# function cv_to_dof(cv,V) -# fv=zeros(eltype(eltype(cv)),num_free_dofs(V)) -# gather_free_values!(fv,V,cv) -# end - -# function field_to_cv(uh::FEFunction) -# get_cell_dof_values(uh) -# end - -# function field_to_cv(cf::CellField) -# cf.cell_field.args[1] -# end - -# function get_geo_params(φh::FEFunction,Vbg) -# Ωbg = get_triangulation(Vbg) -# bgmodel = get_background_model(Ωbg) -# point_to_coords = collect1d(get_node_coordinates(bgmodel)) -# ls_to_point_to_value_unmasked = field_to_cv(φh) -# p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg) -# geo1 = DiscreteGeometry(p0,point_to_coords) -# geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="") -# get_geo_params(geo1,geo2,bgmodel) -# end - -# function get_geo_params(φh::CellField,Vbg) -# Ωbg = get_triangulation(Vbg) -# bgmodel = get_background_model(Ωbg) -# point_to_coords = collect1d(get_node_coordinates(bgmodel)) -# ls_to_point_to_value_unmasked = field_to_cv(φh) -# p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg) -# geo1 = DiscreteGeometry(p0,point_to_coords) -# geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="") -# get_geo_params(geo1,geo2,bgmodel) -# end - -# function get_geo_params(φ::AbstractVector,Vbg) -# Ωbg = get_triangulation(Vbg) -# bgmodel = get_background_model(Ωbg) -# point_to_coords = collect1d(get_node_coordinates(bgmodel)) -# geo1 = DiscreteGeometry(φ,point_to_coords,name="") -# geo2 = DiscreteGeometry(-φ,point_to_coords,name="") -# get_geo_params(geo1,geo2,bgmodel) -# end - -# The above is more efficent for serial problems but does not work for periodic problems or distirbuted -# problems. This is subject to change. - -_DiscreteGeometry(φ,model::DistributedDiscreteModel;name::String="") = - DistributedDiscreteGeometry(φ,model;name) - -_DiscreteGeometry(φh::CellField,model::CartesianDiscreteModel;name::String="") = - DiscreteGeometry(φh(collect1d(get_node_coordinates(model))),collect1d(get_node_coordinates(model));name) - -function get_geo_params(ϕh::CellField,Vbg)#::GridapDistributed.DistributedFESpace) +function get_geo_params(ϕh::CellField,Vbg) Ωbg = get_triangulation(Vbg) bgmodel = get_background_model(Ωbg) - ϕhminus = FEFunction(Vbg,-get_free_dof_values(ϕh)) - geo1 = _DiscreteGeometry(ϕh,bgmodel) - geo2 = _DiscreteGeometry(ϕhminus,bgmodel,name="") + geo1 = DiscreteGeometry(ϕh,bgmodel) + geo2 = DiscreteGeometry(-ϕh,bgmodel,name="") get_geo_params(geo1,geo2,bgmodel) end -function get_geo_params(φ::AbstractVector,Vbg)#::GridapDistributed.DistributedFESpace) +function get_geo_params(φ::AbstractVector,Vbg) Ωbg = get_triangulation(Vbg) bgmodel = get_background_model(Ωbg) ϕh = FEFunction(Vbg,φ); ϕhminus = FEFunction(Vbg,-φ) - geo1 = _DiscreteGeometry(ϕh,bgmodel,name="") - geo2 = _DiscreteGeometry(ϕhminus,bgmodel,name="") + geo1 = DiscreteGeometry(ϕh,bgmodel,name="") + geo2 = DiscreteGeometry(ϕhminus,bgmodel,name="") get_geo_params(geo1,geo2,bgmodel) end @@ -108,14 +54,21 @@ mutable struct EmbeddedMeasureCache end end -function update_embedded_measures!(φ,s::EmbeddedMeasureCache) +function update_embedded_measures!(s::EmbeddedMeasureCache,φ) trians, measures = get_geo_params(φ,s.space) s.measures = measures s.trians = trians return s.measures end -function get_embedded_measures(φ,s::EmbeddedMeasureCache) +function update_embedded_measures!(s::EmbeddedMeasureCache,φ,space) + trians, measures = get_geo_params(φ,space) + s.measures = measures + s.trians = trians + return s.measures +end + +function get_embedded_measures(s::EmbeddedMeasureCache,args...) return s.measures end @@ -136,43 +89,35 @@ end (F::IntegrandWithEmbeddedMeasure)(args...) = F.F(args...,F.dΩ...,F.get_embedded_dΩ(args[end])...) function Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh::Vector{<:FEFunction},K::Int) - @warn "Automatic differentation is currently disabled for `IntegrandWithEmbeddedMeasure` types" maxlog=1 @check 0 < K <= length(uh) _f(uk) = if K == length(uh) F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...,F.get_embedded_dΩ(uk)...) else F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...,F.get_embedded_dΩ(uh[end])...) end - # return Gridap.gradient(_f,uh[K]) # AD is currently disabled due to error (under investigation) - return Gridap.gradient(u -> ∑(∫(0)F.dΩ[i] for i = 1:length(F.dΩ)),uh[K]) + return Gridap.gradient(_f,uh[K]) end -# This doesn't currently work, we need a nice way to differentiate local_embedded_measures -# like the above. function Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh::Vector,K::Int) - @warn "Automatic differentation is currently disabled for `IntegrandWithEmbeddedMeasure` types" maxlog=1 @check 0 < K <= length(uh) local_fields = map(local_views,uh) |> to_parray_of_arrays local_measures = map(local_views,F.dΩ) |> to_parray_of_arrays - # if K == length(uh) - # # Not sure how to do this just yet... - # else - # local_embedded_measures = map(local_views,F.get_embedded_dΩ(uh[end])) |> to_parray_of_arrays - # contribs = map(local_measures,local_embedded_measures,local_fields) do dΩ,dΩ_embedded,lf - # _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,dΩ_embedded...) - # return Gridap.Fields.gradient(_f,lf[K]) - # end - # end - - # Placeholder - local_embedded_measures = map(local_views,F.get_embedded_dΩ(uh[end])) |> to_parray_of_arrays - contribs = map(local_measures,local_embedded_measures,local_fields) do dΩ,dΩ_embedded,lf - _f = u -> ∑(∫(0)dΩ[i] for i = 1:length(dΩ)) - return Gridap.Fields.gradient(_f,lf[K]) + if K == length(uh) + # Need to test. My hope is that locally this will work - ghosts will be wrong but + # will be thrown away later. + contribs = map(local_measures,local_views(uh[end]),local_fields) do dΩ,φh,lf + _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,F.get_embedded_dΩ(φh,get_fe_space(φh))...) + return Gridap.Fields.gradient(_f,lf[K]) + end + else + local_embedded_measures = map(local_views,F.get_embedded_dΩ(uh[end])) |> to_parray_of_arrays + contribs = map(local_measures,local_embedded_measures,local_fields) do dΩ,dΩ_embedded,lf + _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,dΩ_embedded...) + return Gridap.Fields.gradient(_f,lf[K]) + end end return DistributedDomainContribution(contribs) end - Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh) = Gridap.gradient(F,[uh],1) \ No newline at end of file diff --git a/scripts/Embedded/embedded_measures_AD_DISABLED.jl b/scripts/Embedded/embedded_measures_AD_DISABLED.jl new file mode 100644 index 00000000..d60af406 --- /dev/null +++ b/scripts/Embedded/embedded_measures_AD_DISABLED.jl @@ -0,0 +1,180 @@ +using GridapTopOpt: to_parray_of_arrays +using GridapDistributed +using GridapDistributed: DistributedDiscreteModel, DistributedDomainContribution +using GridapEmbedded.LevelSetCutters +# using GridapEmbedded.LevelSetCutters: DiscreteGeometry +using GridapEmbedded.Distributed: DistributedDiscreteGeometry +using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Helpers +import Gridap.Geometry: get_node_coordinates, collect1d + +# function cv_to_dof(cv,V) +# fv=zeros(eltype(eltype(cv)),num_free_dofs(V)) +# gather_free_values!(fv,V,cv) +# end + +# function field_to_cv(uh::FEFunction) +# get_cell_dof_values(uh) +# end + +# function field_to_cv(cf::CellField) +# cf.cell_field.args[1] +# end + +# function get_geo_params(φh::FEFunction,Vbg) +# Ωbg = get_triangulation(Vbg) +# bgmodel = get_background_model(Ωbg) +# point_to_coords = collect1d(get_node_coordinates(bgmodel)) +# ls_to_point_to_value_unmasked = field_to_cv(φh) +# p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg) +# geo1 = DiscreteGeometry(p0,point_to_coords) +# geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="") +# get_geo_params(geo1,geo2,bgmodel) +# end + +# function get_geo_params(φh::CellField,Vbg) +# Ωbg = get_triangulation(Vbg) +# bgmodel = get_background_model(Ωbg) +# point_to_coords = collect1d(get_node_coordinates(bgmodel)) +# ls_to_point_to_value_unmasked = field_to_cv(φh) +# p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg) +# geo1 = DiscreteGeometry(p0,point_to_coords) +# geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="") +# get_geo_params(geo1,geo2,bgmodel) +# end + +# function get_geo_params(φ::AbstractVector,Vbg) +# Ωbg = get_triangulation(Vbg) +# bgmodel = get_background_model(Ωbg) +# point_to_coords = collect1d(get_node_coordinates(bgmodel)) +# geo1 = DiscreteGeometry(φ,point_to_coords,name="") +# geo2 = DiscreteGeometry(-φ,point_to_coords,name="") +# get_geo_params(geo1,geo2,bgmodel) +# end + +# The above is more efficent for serial problems but does not work for periodic problems or distirbuted +# problems. This is subject to change. + +# No longer required +# _DiscreteGeometry(φ,model::DistributedDiscreteModel;name::String="") = +# DistributedDiscreteGeometry(φ,model;name) + +# _DiscreteGeometry(φh::CellField,model::CartesianDiscreteModel;name::String="") = +# DiscreteGeometry(φh(collect1d(get_node_coordinates(model))),collect1d(get_node_coordinates(model));name) + +function get_geo_params(ϕh::CellField,Vbg)#::GridapDistributed.DistributedFESpace) + Ωbg = get_triangulation(Vbg) + bgmodel = get_background_model(Ωbg) + geo1 = DiscreteGeometry(ϕh,bgmodel) + geo2 = DiscreteGeometry(-ϕh,bgmodel,name="") + get_geo_params(geo1,geo2,bgmodel) +end + +function get_geo_params(φ::AbstractVector,Vbg)#::GridapDistributed.DistributedFESpace) + Ωbg = get_triangulation(Vbg) + bgmodel = get_background_model(Ωbg) + ϕh = FEFunction(Vbg,φ); ϕhminus = FEFunction(Vbg,-φ) + geo1 = DiscreteGeometry(ϕh,bgmodel,name="") + geo2 = DiscreteGeometry(ϕhminus,bgmodel,name="") + get_geo_params(geo1,geo2,bgmodel) +end + +function get_geo_params(geo1,geo2,bgmodel) + cutgeo1= cut(bgmodel,geo1) + cutgeo2= cut(bgmodel,geo2) + # Setup integration meshes + Ω1 = Triangulation(cutgeo1,PHYSICAL) + Ω2 = Triangulation(cutgeo2,PHYSICAL) + Γ = EmbeddedBoundary(cutgeo1) + # Setup Lebesgue measures + order = 1 + degree = 2*order + dΩ1 = Measure(Ω1,degree) + dΩ2 = Measure(Ω2,degree) + dΓ = Measure(Γ,degree) + (Ω1,Ω2,Γ),(dΩ1,dΩ2,dΓ) +end + +####################### + +# For problems where the derivatives are known, we only want to update measures once +mutable struct EmbeddedMeasureCache + const space + trians + measures + + function EmbeddedMeasureCache(φ,space) + measures = get_geo_params(φ,space) + new(space,measures) + end +end + +function update_embedded_measures!(φ,s::EmbeddedMeasureCache) + trians, measures = get_geo_params(φ,s.space) + s.measures = measures + s.trians = trians + return s.measures +end + +function get_embedded_measures(φ,s::EmbeddedMeasureCache) + return s.measures +end + +function get_embedded_triangulations(s::EmbeddedMeasureCache) + return s.trians +end + +####################### + +import GridapTopOpt: AbstractIntegrandWithMeasure + +struct IntegrandWithEmbeddedMeasure{A,B,C} <: AbstractIntegrandWithMeasure + F :: A + dΩ :: B + get_embedded_dΩ :: C +end + +(F::IntegrandWithEmbeddedMeasure)(args...) = F.F(args...,F.dΩ...,F.get_embedded_dΩ(args[end])...) + +function Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh::Vector{<:FEFunction},K::Int) + @warn "Automatic differentation is currently disabled for `IntegrandWithEmbeddedMeasure` types" maxlog=1 + @check 0 < K <= length(uh) + _f(uk) = if K == length(uh) + F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...,F.get_embedded_dΩ(uk)...) + else + F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...,F.get_embedded_dΩ(uh[end])...) + end + # return Gridap.gradient(_f,uh[K]) # AD is currently disabled due to error (under investigation) + return Gridap.gradient(u -> ∑(∫(0)F.dΩ[i] for i = 1:length(F.dΩ)),uh[K]) +end + +function Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh::Vector,K::Int) + @warn "Automatic differentation is currently disabled for `IntegrandWithEmbeddedMeasure` types" maxlog=1 + @check 0 < K <= length(uh) + local_fields = map(local_views,uh) |> to_parray_of_arrays + local_measures = map(local_views,F.dΩ) |> to_parray_of_arrays + + # if K == length(uh) + # # Need to test. My hope is that locally this will work - ghosts will be wrong but + # # will be thrown away later. + # contribs = map(local_measures,local_views(uh[end]),local_fields) do dΩ,φh,lf + # _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,F.get_embedded_dΩ(φh)...) + # return Gridap.Fields.gradient(_f,lf[K]) + # end + # else + # local_embedded_measures = map(local_views,F.get_embedded_dΩ(uh[end])) |> to_parray_of_arrays + # contribs = map(local_measures,local_embedded_measures,local_fields) do dΩ,dΩ_embedded,lf + # _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,dΩ_embedded...) + # return Gridap.Fields.gradient(_f,lf[K]) + # end + # end + + # Placeholder + local_embedded_measures = map(local_views,F.get_embedded_dΩ(uh[end])) |> to_parray_of_arrays + contribs = map(local_measures,local_embedded_measures,local_fields) do dΩ,dΩ_embedded,lf + _f = u -> ∑(∫(0)dΩ[i] for i = 1:length(dΩ)) + return Gridap.Fields.gradient(_f,lf[K]) + end + + return DistributedDomainContribution(contribs) +end +Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh) = Gridap.gradient(F,[uh],1) \ No newline at end of file