Merge pull request #112 from zjwegert/fix-diff_trian-multifield

Add missing update_trian! method for MultiField

Author: JordiManyer

NEWS.md

@@ -4,6 +4,16 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.9.8] - 2025-7-10
+
+### Added
+
+- Added missing update_trian! method for MultiField. Since PR[#112](https://github.com/gridap/GridapEmbedded.jl/pull/112).
+
+### Fixed
+
+- Fixed failed precompilation due to Gridap v0.19.2. Since PR[#112](https://github.com/gridap/GridapEmbedded.jl/pull/112).
+
 ## [0.9.7] - 2025-6-11
 
 ### Added

Project.toml

@@ -1,7 +1,7 @@
 name = "GridapEmbedded"
 uuid = "8838a6a3-0006-4405-b874-385995508d5d"
 authors = ["Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Eric Neiva <eric.neiva@college-de-france.fr>", "Pere Antoni Martorell <pere.antoni.martorell@upc.edu>", "Santiago Badia <santiago.badia@monash.edu>"]
-version = "0.9.7"
+version = "0.9.8"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

src/Interfaces/CutFaceBoundaryTriangulations.jl

@@ -395,38 +395,38 @@ end
 
 ############################################################################################
 # This will go to Gridap
-
-function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField})
-  plus = k(a.plus)
-  minus = k(a.minus)
-  SkeletonPair(plus,minus)
-end
-
-function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField},b::SkeletonPair{<:CellField})
-  plus = k(a.plus,b.plus)
-  minus = k(a.minus,b.minus)
-  SkeletonPair(plus,minus)
-end
-
-import Gridap.TensorValues: inner, outer
-import LinearAlgebra: dot
-import Base: abs, *, +, -, /
-
-for op in (:/,)
-  @eval begin
-    ($op)(a::CellField,b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
-    ($op)(a::SkeletonPair{<:CellField},b::CellField) = Operation($op)(a,b)
-  end
-end
-
-for op in (:outer,:*,:dot,:/)
-  @eval begin
-    ($op)(a::SkeletonPair{<:CellField},b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
-  end
-end
-
-function CellData.change_domain(a::SkeletonPair, ::ReferenceDomain, ::PhysicalDomain)
-  plus = change_domain(a.plus,ReferenceDomain(),PhysicalDomain())
-  minus = change_domain(a.minus,ReferenceDomain(),PhysicalDomain())
-  return SkeletonPair(plus,minus)
-end
+#
+# function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField})
+#   plus = k(a.plus)
+#   minus = k(a.minus)
+#   SkeletonPair(plus,minus)
+# end
+
+# function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField},b::SkeletonPair{<:CellField})
+#   plus = k(a.plus,b.plus)
+#   minus = k(a.minus,b.minus)
+#   SkeletonPair(plus,minus)
+# end
+
+# import Gridap.TensorValues: inner, outer
+# import LinearAlgebra: dot
+# import Base: abs, *, +, -, /
+
+# for op in (:/,)
+#   @eval begin
+#     ($op)(a::CellField,b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
+#     ($op)(a::SkeletonPair{<:CellField},b::CellField) = Operation($op)(a,b)
+#   end
+# end
+
+# for op in (:outer,:*,:dot,:/)
+#   @eval begin
+#     ($op)(a::SkeletonPair{<:CellField},b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
+#   end
+# end
+
+# function CellData.change_domain(a::SkeletonPair, ::ReferenceDomain, ::PhysicalDomain)
+#   plus = change_domain(a.plus,ReferenceDomain(),PhysicalDomain())
+#   minus = change_domain(a.minus,ReferenceDomain(),PhysicalDomain())
+#   return SkeletonPair(plus,minus)
+# end

src/LevelSetCutters/DifferentiableTriangulations.jl

@@ -10,11 +10,11 @@ methods to compute derivatives w.r.t. deformations of the embedded mesh.
 To do so, it propagates dual numbers into the geometric maps mapping cut subcells/subfacets
 to the background mesh.
 
-## Constructor: 
+## Constructor:
 
     DifferentiableTriangulation(trian::Triangulation,fe_space::FESpace)
 
-where `trian` must be an embedded triangulation and `fe_space` is the `FESpace` where 
+where `trian` must be an embedded triangulation and `fe_space` is the `FESpace` where
 the level-set function lives.
 
 """
@@ -63,6 +63,18 @@ function update_trian!(trian::DifferentiableTriangulation,::FESpace,::Nothing)
   return trian
 end
 
+const MultiFieldSpaceTypes = Union{<:MultiFieldFESpace,<:DistributedMultiFieldFESpace}
+
+function update_trian!(trian::DifferentiableTriangulation,space::MultiFieldSpaceTypes,φh)
+  map((Ui,φi)->update_trian!(trian,Ui,φi),space,φh)
+  return trian
+end
+
+function update_trian!(trian::DifferentiableTriangulation,::MultiFieldSpaceTypes,::Nothing)
+  trian.cell_values = nothing
+  return trian
+end
+
 # Autodiff
 
 function FESpaces._change_argument(
@@ -424,8 +436,8 @@ function extract_dualized_cell_values(
 end
 
 # TriangulationView
-# This is mostly used in distributed, where we remove ghost cells by taking a view 
-# of the local triangulations. 
+# This is mostly used in distributed, where we remove ghost cells by taking a view
+# of the local triangulations.
 
 const DifferentiableTriangulationView{Dc,Dp} = Geometry.TriangulationView{Dc,Dp,<:DifferentiableTriangulation}
 
@@ -465,7 +477,7 @@ end
 # We only need to propagate the dual numbers to the CUT cells, which is what the
 # following implementation does:
 
-const DifferentiableAppendedTriangulation{Dc,Dp,A} = 
+const DifferentiableAppendedTriangulation{Dc,Dp,A} =
   AppendedTriangulation{Dc,Dp,<:Union{<:DifferentiableTriangulation,<:DifferentiableTriangulationView{Dc,Dp}}}
 
 function DifferentiableTriangulation(

src/LevelSetCutters/LevelSetCutters.jl

@@ -31,6 +31,9 @@ using Gridap.CellData
 using Gridap.Polynomials
 using Gridap.Visualization
 using Gridap.FESpaces
+using Gridap.MultiField
+
+using GridapDistributed: DistributedMultiFieldFESpace
 
 export LevelSetCutter
 export AnalyticalGeometry

test/DistributedTests/GeometricalDifferentiationTests.jl

@@ -1,7 +1,7 @@
 module DistributedGeometricalDifferentitationTests
 ############################################################################################
-# These tests are meant to verify the correctness differentiation of functionals w.r.t the 
-# level set defining the cut domain. 
+# These tests are meant to verify the correctness differentiation of functionals w.r.t the
+# level set defining the cut domain.
 # They are based on the following work:
 # "Level-set topology optimisation with unfitted finite elements and automatic shape differentiation"
 #   by Z. J. Wegert, J. Manyer, C. Mallon, S. Badia, V. J. Challis (2025)
@@ -117,6 +117,19 @@ function main_generic(
 
   @test norm(dJ_bulk_1_AD_in_u_vec - dJ_bulk_1_exact_in_u_vec) < 1e-10
 
+  # Multifield case
+  U = TestFESpace(model,reffe)
+  V_φu = MultiFieldFESpace([V_φ,U])
+  φuh = interpolate([φh,x->x[1]],V_φu)
+  J_bulk_mult((φ,u)) = ∫(fh)dΩ
+  dJ_bulk_AD = gradient(J_bulk_mult,φuh)
+  dJ_bulk_AD_vec = assemble_vector(dJ_bulk_AD,V_φu)
+
+  dJ_bulk_exact_mult((dφ,du)) = ∫(-fh*dφ/(abs(n_Γ ⋅ ∇(φh))))dΓ
+  dJ_bulk_exact_vec = assemble_vector(dJ_bulk_exact_mult,V_φu)
+
+  @test norm(dJ_bulk_AD_vec - dJ_bulk_exact_vec) < 1e-10
+
   # A.2) Volume integral
 
   g(fh) = ∇(fh)⋅∇(fh)

test/LevelSetCuttersTests/GeometricalDifferentiationTests.jl

@@ -1,14 +1,14 @@
 module GeometricalDifferentiationTests
 ############################################################################################
-# These tests are meant to verify the correctness differentiation of functionals w.r.t the 
-# level set defining the cut domain. 
+# These tests are meant to verify the correctness differentiation of functionals w.r.t the
+# level set defining the cut domain.
 # They are based on the following work:
 # "Level-set topology optimisation with unfitted finite elements and automatic shape differentiation"
 #   by Z. J. Wegert, J. Manyer, C. Mallon, S. Badia, V. J. Challis (2025)
 ############################################################################################
 using Test, FiniteDiff
 
-using Gridap, Gridap.Geometry, Gridap.Adaptivity, Gridap.Arrays
+using Gridap, Gridap.Geometry, Gridap.Adaptivity, Gridap.Arrays, Gridap.MultiField
 using GridapEmbedded, GridapEmbedded.LevelSetCutters, GridapEmbedded.Interfaces
 
 using GridapEmbedded.Interfaces: get_conormal_vector
@@ -17,7 +17,7 @@ using GridapEmbedded.Interfaces: get_ghost_normal_vector
 
 using GridapEmbedded.LevelSetCutters: DifferentiableTriangulation
 
-# We general a simplicial model where the simplices are created in a symmetric way using 
+# We general a simplicial model where the simplices are created in a symmetric way using
 # varycentric refinement of QUADs and HEXs.
 function generate_model(D,n)
   domain = (D==2) ? (0,1,0,1) : (0,1,0,1,0,1)
@@ -124,6 +124,23 @@ function main(
 
   @test abs_error < 1e-10
 
+  # Multifield case
+  U = TestFESpace(model,reffe)
+  V_φu = MultiFieldFESpace([V_φ,U])
+  φuh = interpolate([φh,x->x[1]],V_φu)
+  J_bulk_mult((φ,u)) = ∫(fh)dΩ
+  dJ_bulk_AD = gradient(J_bulk_mult,φuh)
+  dJ_bulk_AD_vec = assemble_vector(dJ_bulk_AD,V_φu)
+
+  dJ_bulk_exact_mult((dφ,du)) = ∫(-fh*dφ/(abs(n_Γ ⋅ ∇(φh))))dΓ
+  dJ_bulk_exact_vec = assemble_vector(dJ_bulk_exact_mult,V_φu)
+
+  abs_error = norm(dJ_bulk_AD_vec - dJ_bulk_exact_vec,Inf)
+  if verbose
+    println("  - norm(dJ_AD - dJ_exact,Inf) = ",abs_error,"    | Multifield test")
+  end
+  @test abs_error < 1e-10
+
   # A.1.1) Volume integral with another field
 
   J_bulk_1(u,φ) = ∫(u+fh)dΩ