Merge pull request #57 from gridap/gridap_v0.17

Gridap v0.17

.github/workflows/ci.yml

@@ -8,7 +8,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
+          - '1.6'
         os:
           - ubuntu-latest
         arch:
@@ -42,7 +42,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
+          - '1.6'
         os:
           - ubuntu-latest
         arch:
@@ -77,7 +77,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.5'
+          version: '1.6'
       - run: |
           julia --project=docs -e '
             using Pkg
@@ -91,4 +91,4 @@ jobs:
       - run: julia --project=docs docs/make.jl
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

NEWS.md

@@ -4,6 +4,12 @@ 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.8.0] - 2021-11-03
+
+### Changed
+
+ - Switch to Gridap 0.17. This required some changes in the user API. In particular, one needs to specify now if you want to build a triangulation of the "physical" domain or of the "active" domain. E.g, `Triangulation(cutgeo,PHYSICAL)` or `Triangulation(cutgeo,ACTIVE)`. Taking advantage that one can define FE spaces on "body-fitted" triangulation in Gridap v0.17, the user builds the FE space in embedded simulation by using the "active" triangulation.
+
 ## [0.7.0] - 2021-06-29
 
 ### Changed

Project.toml

@@ -1,7 +1,7 @@
 name = "GridapEmbedded"
 uuid = "8838a6a3-0006-4405-b874-385995508d5d"
 authors = ["Francesc Verdugo <[email protected]>"]
-version = "0.7.0"
+version = "0.8.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -18,7 +18,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 AbstractTrees = "0.3.3"
 Combinatorics = "1"
 FillArrays = "0.10, 0.11"
-Gridap = "0.16.3"
+Gridap = "0.17"
 LightGraphs = "1.3.3"
 MiniQhull = "0.1.0, 0.2, 0.3"
 julia = "1.3"

examples/BimaterialLinElastCutFEM/BimaterialLinElastCutFEM.jl

@@ -18,6 +18,7 @@ function main(;n,outputfile=nothing)
   pmin = Point(0.,0.,0.)
   pmax = pmin + L
   bgmodel = CartesianDiscreteModel(pmin,pmax,partition)
+  Ω_bg = Triangulation(bgmodel)
 
   # Identify Dirichlet boundaries
   labeling = get_face_labeling(bgmodel)
@@ -54,13 +55,13 @@ function main(;n,outputfile=nothing)
   # Cut the background model
   cutgeo = cut(bgmodel,union(geo3,geo4))
 
-  # Setup models
-  model1 = DiscreteModel(cutgeo,"steel")
-  model2 = DiscreteModel(cutgeo,"concrete")
+  # Setup interpolation mesh
+  Ω1_act = Triangulation(cutgeo,ACTIVE,"steel")
+  Ω2_act = Triangulation(cutgeo,ACTIVE,"concrete")
 
   # Setup integration meshes
-  Ω1 = Triangulation(cutgeo,"steel")
-  Ω2 = Triangulation(cutgeo,"concrete")
+  Ω1 = Triangulation(cutgeo,PHYSICAL,"steel")
+  Ω2 = Triangulation(cutgeo,PHYSICAL,"concrete")
   Γ = EmbeddedBoundary(cutgeo,"steel","concrete")
 
   # Setup normal vectors
@@ -75,11 +76,11 @@ function main(;n,outputfile=nothing)
 
   # Setup FESpace
 
-  V1 = TestFESpace(model1,
+  V1 = TestFESpace(Ω1_act,
                    ReferenceFE(lagrangian,VectorValue{3,Float64},order),
                    conformity=:H1,
                    dirichlet_tags=["support0","support1"])
-  V2 = TestFESpace(model2,
+  V2 = TestFESpace(Ω2_act,
                    ReferenceFE(lagrangian,VectorValue{3,Float64},order),
                    conformity=:H1,
                    dirichlet_tags=["support0","support1"])
@@ -92,14 +93,14 @@ function main(;n,outputfile=nothing)
 
   # Setup stabilization parameters
 
-  meas_K1 = get_cell_measure(Ω1)
-  meas_K2 = get_cell_measure(Ω2)
-  meas_KΓ = get_cell_measure(Γ)
+  meas_K1 = get_cell_measure(Ω1, Ω_bg)
+  meas_K2 = get_cell_measure(Ω2, Ω_bg)
+  meas_KΓ = get_cell_measure(Γ, Ω_bg)
 
   γ_hat = 2
-  κ1 = (E2*meas_K1) ./ (E2*meas_K1 .+ E1*meas_K2)
-  κ2 = (E1*meas_K2) ./ (E2*meas_K1 .+ E1*meas_K2)
-  β = (γ_hat*meas_KΓ) ./ ( meas_K1/E1 .+ meas_K2/E2 )
+  κ1 = CellField( (E2*meas_K1) ./ (E2*meas_K1 .+ E1*meas_K2), Ω_bg)
+  κ2 = CellField( (E1*meas_K2) ./ (E2*meas_K1 .+ E1*meas_K2), Ω_bg)
+  β  = CellField( (γ_hat*meas_KΓ) ./ ( meas_K1/E1 .+ meas_K2/E2 ), Ω_bg)
 
   # Jump and mean operators for this formulation
 

examples/PoissonCSGCutFEM/PoissonCSGCutFEM.jl

@@ -30,7 +30,7 @@ function main(;n,outputfile=nothing)
   cutgeo = cut(bgmodel,geo8)
 
   # Setup integration meshes
-  Ω = Triangulation(cutgeo,"csg")
+  Ω = Triangulation(cutgeo,PHYSICAL,"csg")
   Γd = EmbeddedBoundary(cutgeo,"csg","source")
   Γg = GhostSkeleton(cutgeo,"csg")
 
@@ -51,8 +51,8 @@ function main(;n,outputfile=nothing)
   dΓg = Measure(Γg,degree)
 
   # Setup FESpace
-  model = DiscreteModel(cutgeo,"csg")
-  V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
+  Ω_act = Triangulation(cutgeo,ACTIVE,"csg")
+  V = TestFESpace(Ω_act,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
   U = TrialFESpace(V)
 
   # Weak form

examples/StokesTubeWithObstacleCutFEM/StokesTubeWithObstacleCutFEM.jl

@@ -50,11 +50,11 @@ function main(;n,outputfile=nothing)
   # Cut the background model
   cutgeo = cut(bgmodel,union(geo1,geo5))
 
-  # Generate the "active" model
-  model = DiscreteModel(cutgeo,"fluid")
+  # Generate the "active" mesh
+  Ω_act = Triangulation(cutgeo,ACTIVE,"fluid")
 
   # Setup integration meshes
-  Ω = Triangulation(cutgeo,"fluid")
+  Ω = Triangulation(cutgeo,PHYSICAL,"fluid")
   Γi = EmbeddedBoundary(cutgeo,"fluid","inlet")
   Γw = EmbeddedBoundary(cutgeo,"fluid","solid")
   Γg = GhostSkeleton(cutgeo,"fluid")
@@ -82,8 +82,8 @@ function main(;n,outputfile=nothing)
   reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P)
   reffe_p = ReferenceFE(lagrangian,Float64,order,space=:P)
 
-  V = TestFESpace(model,reffe_u,conformity=:H1)
-  Q = TestFESpace(model,reffe_p,conformity=:H1)
+  V = TestFESpace(Ω_act,reffe_u,conformity=:H1)
+  Q = TestFESpace(Ω_act,reffe_p,conformity=:H1)
 
   U = TrialFESpace(V)
   P = TrialFESpace(Q)

src/AgFEM/AgFEMSpaces.jl

@@ -1,60 +1,44 @@
 
-function AgFEMSpace(f::SingleFieldFESpace,cell_to_cellin::AbstractVector,g::SingleFieldFESpace=f)
-  AgFEMSpace(f,cell_to_cellin,get_fe_basis(g),get_fe_dof_basis(g))
+function AgFEMSpace(f::SingleFieldFESpace,bgcell_to_bgcellin::AbstractVector,g::SingleFieldFESpace=f)
+  @assert get_triangulation(f) === get_triangulation(g)
+  AgFEMSpace(f,bgcell_to_bgcellin,get_fe_basis(g),get_fe_dof_basis(g))
 end
 
 # Note: cell is in fact bgcell in this function since f will usually be an ExtendedFESpace
 function AgFEMSpace(
   f::SingleFieldFESpace,
-  cell_to_cellin::AbstractVector,
-  cell_shapefuns_g::CellField,
-  cell_dof_basis_g::CellDof)
-
-  # Prepare maps between different cell ids
-  cell_to_isactive = lazy_map(i->(i>0),cell_to_cellin)
-  acell_to_cell = findall( cell_to_isactive  )
-  acell_to_cellin = cell_to_cellin[acell_to_cell]
-  cell_to_acell = zeros(Int32,length(cell_to_cellin))
-  cell_to_acell[cell_to_isactive] .= 1:length(acell_to_cell)
+  bgcell_to_bgcellin::AbstractVector,
+  shfns_g::CellField,
+  dofs_g::CellDof)
 
   # Triangulation made of active cells
-  trian = get_triangulation(f)
-  trian_a = Triangulation(trian,acell_to_cell)
-
-  # Celldata Defined on trian
+  trian_a = get_triangulation(f)
+
+  # Build root cell map (i.e. aggregates) in terms of active cell ids
+  D = num_cell_dims(trian_a)
+  glue = get_glue(trian_a,Val(D))
+  acell_to_bgcell = glue.tface_to_mface
+  bgcell_to_acell = glue.mface_to_tface
+  acell_to_bgcellin = lazy_map(Reindex(bgcell_to_bgcellin),acell_to_bgcell)
+  acell_to_acellin = collect(lazy_map(Reindex(bgcell_to_acell),acell_to_bgcellin))
+
+  # Build shape funs of g by replacing local funs in cut cells by the ones at the root
+  # This needs to be done with shape functions in the physical domain
+  # otherwise shape funs in cut and root cells are the same
+  acell_phys_shapefuns_g = get_array(change_domain(shfns_g,PhysicalDomain()))
+  acell_phys_root_shapefuns_g = lazy_map(Reindex(acell_phys_shapefuns_g),acell_to_acellin)
+  root_shfns_g = GenericCellField(acell_phys_root_shapefuns_g,trian_a,PhysicalDomain())
+
+  # Compute data needed to compute the constraints
   dofs_f = get_fe_dof_basis(f)
   shfns_f = get_fe_basis(f)
-  dofs_g = cell_dof_basis_g
-  shfns_g = cell_shapefuns_g
-
-  # Celldata Defined on trian_a
-  dofs_f_a = change_domain(dofs_f,trian_a,DomainStyle(dofs_f))
-  cell_phys_shapefuns_g = get_array(change_domain(shfns_g,PhysicalDomain()))
-  acell_phys_shapefuns_g = lazy_map(Reindex(cell_phys_shapefuns_g),acell_to_cellin)
-  shfns_g_a = GenericCellField(acell_phys_shapefuns_g,trian_a,PhysicalDomain())
-
-  acell_to_coeffs = dofs_f_a(shfns_g_a)
-  cell_to_proj = dofs_g(shfns_f)
-  acell_to_proj = lazy_map(Reindex(cell_to_proj),acell_to_cellin)
-  acell_to_dof_ids = lazy_map(Reindex(get_cell_dof_ids(f)),acell_to_cell)
-
-  #acell_to_dofs = reindex(get_cell_dofs(f),acell_to_cell)
-  #n_fdofs =
-  #acell_to_fbasis = reindex(get_cell_basis(f),acell_to_cellin)
-  #acell_to_gbasis = reindex(cell_shapefuns_g,acell_to_cellin)
-  #acell_to_dof_fbasis = reindex(get_fe_dof_basis(f),acell_to_cell)
-  #acell_to_dof_gbasis = reindex(cell_dof_basis_g,acell_to_cellin)
-  #@notimplementedif is_in_ref_space(acell_to_dof_fbasis)
-  #@notimplementedif is_in_ref_space(acell_to_gbasis)
-  #@assert RefStyle(acell_to_fbasis) == RefStyle(acell_to_dof_gbasis)
-  #acell_to_coeffs = evaluate(acell_to_dof_fbasis,acell_to_gbasis)
-  #acell_to_proj = evaluate(acell_to_dof_gbasis,acell_to_fbasis)
+  acell_to_coeffs = dofs_f(root_shfns_g)
+  acell_to_proj = dofs_g(shfns_f)
+  acell_to_dof_ids = get_cell_dof_ids(f)
 
   aggdof_to_fdof, aggdof_to_dofs, aggdof_to_coeffs = _setup_agfem_constraints(
     num_free_dofs(f),
-    acell_to_cellin,
-    acell_to_cell,
-    cell_to_acell,
+    acell_to_acellin,
     acell_to_dof_ids,
     acell_to_coeffs,
     acell_to_proj)
@@ -64,22 +48,20 @@ end
 
 function _setup_agfem_constraints(
   n_fdofs,
-  acell_to_cellin,
-  acell_to_cell,
-  cell_to_acell,
+  acell_to_acellin,
   acell_to_dof_ids,
   acell_to_coeffs,
   acell_to_proj)
 
-  n_acells = length(acell_to_cell)
+  n_acells = length(acell_to_acellin)
   fdof_to_isagg = fill(true,n_fdofs)
   fdof_to_acell = zeros(Int32,n_fdofs)
   fdof_to_ldof = zeros(Int8,n_fdofs)
   cache = array_cache(acell_to_dof_ids)
   for acell in 1:n_acells
-    iscut = acell_to_cell[acell] != acell_to_cellin[acell]
+    acellin = acell_to_acellin[acell]
+    iscut = acell != acellin
     dofs = getindex!(cache,acell_to_dof_ids,acell)
-    cellin = acell_to_cellin[acell]
     for (ldof,dof) in enumerate(dofs)
       if dof > 0
         fdof = dof
@@ -98,9 +80,8 @@ function _setup_agfem_constraints(
   for aggdof in 1:n_aggdofs
     fdof = aggdof_to_fdof[aggdof]
     acell = fdof_to_acell[fdof]
-    cellin = acell_to_cellin[acell]
-    acell = cell_to_acell[cellin]
-    dofs = getindex!(cache,acell_to_dof_ids,acell)
+    acellin = acell_to_acellin[acell]
+    dofs = getindex!(cache,acell_to_dof_ids,acellin)
     aggdof_to_dofs_ptrs[aggdof+1] = length(dofs)
   end
 
@@ -111,9 +92,8 @@ function _setup_agfem_constraints(
   for aggdof in 1:n_aggdofs
     fdof = aggdof_to_fdof[aggdof]
     acell = fdof_to_acell[fdof]
-    cellin = acell_to_cellin[acell]
-    acell = cell_to_acell[cellin]
-    dofs = getindex!(cache,acell_to_dof_ids,acell)
+    acellin = acell_to_acellin[acell]
+    dofs = getindex!(cache,acell_to_dof_ids,acellin)
     p = aggdof_to_dofs_ptrs[aggdof]-1
     for (i,dof) in enumerate(dofs)
       aggdof_to_dofs_data[p+i] = dof

src/AgFEM/CellAggregation.jl

@@ -104,11 +104,11 @@ end
 function _aggregate_by_threshold(threshold,cut,geo,loc,facet_to_inoutcut)
   @assert loc in (IN,OUT)
 
-  cutinorout = loc == IN ? (CUTIN,IN) : (CUTOUT,OUT)
-  trian = Triangulation(cut,geo,cutinorout)
+  cutinorout = loc == IN ? (CUT_IN,IN) : (CUT_OUT,OUT)
+  trian = Triangulation(cut,cutinorout,geo)
   model = cut.bgmodel
-  cell_to_cut_meas = get_cell_measure(trian)
   bgtrian = get_triangulation(model)
+  cell_to_cut_meas = get_cell_measure(trian,bgtrian)
   cell_to_meas = get_cell_measure(bgtrian)
   cell_to_unit_cut_meas = lazy_map(/,cell_to_cut_meas,cell_to_meas)
 

src/Exports.jl

@@ -18,6 +18,14 @@ end
 @publish Interfaces GhostSkeleton
 @publish Interfaces IN
 @publish Interfaces OUT
+@publish Interfaces CUT_IN
+@publish Interfaces CUT_OUT
+@publish Interfaces ACTIVE
+@publish Interfaces ACTIVE_IN
+@publish Interfaces ACTIVE_OUT
+@publish Interfaces PHYSICAL
+@publish Interfaces PHYSICAL_IN
+@publish Interfaces PHYSICAL_OUT
 
 @publish LevelSetCutters LevelSetCutter
 @publish LevelSetCutters AnalyticalGeometry