Merge branch 'main' into dev_io_interface

docs/src/howto/howto.md

@@ -88,7 +88,7 @@ Open from `Bcube.jl/` a REPL and type:
 
 ```julia
 pkg> activate --temp
-pkg> add PkgBenchmark BenchmarkTools
+pkg> add BenchmarkTools PkgBenchmark StaticArrays WriteVTK UnPack
 pkg> dev .
 using PkgBenchmark
 import Bcube
@@ -101,7 +101,7 @@ Then checkout the `main` branch. Start a fresh REPL and type (almost the same):
 
 ```julia
 pkg> activate --temp
-pkg> add PkgBenchmark BenchmarkTools
+pkg> add BenchmarkTools PkgBenchmark StaticArrays WriteVTK UnPack
 pkg> dev .
 using PkgBenchmark
 import Bcube
@@ -151,7 +151,7 @@ Open from `Bcube.jl/` a REPL and type:
 
 ```julia
 pkg> activate --temp
-pkg> add PkgBenchmark
+pkg> add BenchmarkTools PkgBenchmark StaticArrays WriteVTK UnPack
 pkg> dev .
 using PkgBenchmark
 import Bcube

src/Bcube.jl

@@ -51,18 +51,6 @@ include("./mesh/connectivity.jl")
 include("./mesh/mesh.jl")
 export ncells, nnodes, boundary_names, nboundaries, boundary_tag, get_nodes
 
-include("./mesh/gmsh_utils.jl")
-export read_msh,
-    read_msh_with_cell_names,
-    gen_line_mesh,
-    gen_rectangle_mesh,
-    gen_hexa_mesh,
-    gen_disk_mesh,
-    gen_star_disk_mesh,
-    gen_cylinder_mesh,
-    read_partitions,
-    gen_rectangle_mesh_with_tri_and_quad
-
 include("./mesh/mesh_generator.jl")
 export basic_mesh,
     one_cell_mesh,
@@ -78,6 +66,18 @@ export basic_mesh,
     translate,
     translate!
 
+include("./mesh/gmsh_utils.jl")
+export read_msh,
+    read_msh_with_cell_names,
+    gen_line_mesh,
+    gen_rectangle_mesh,
+    gen_hexa_mesh,
+    gen_disk_mesh,
+    gen_star_disk_mesh,
+    gen_cylinder_mesh,
+    read_partitions,
+    gen_rectangle_mesh_with_tri_and_quad
+
 include("./mesh/domain.jl")
 export AbstractDomain,
     CellDomain, InteriorFaceDomain, BoundaryFaceDomain, get_mesh, get_face_normals

src/feoperator/projection_newapi.jl

@@ -120,6 +120,8 @@ end
 _may_scalar(a::SVector{1}) = a[1]
 _may_scalar(a) = a
 
+cell_mean(u::MeshData{CellData}, ::Measure) = u
+
 function build_mass_matrix(u::AbstractFEFunction, dΩ::AbstractMeasure)
     U = get_fespace(u)
     V = TestFESpace(U)

src/mesh/gmsh_utils.jl

@@ -134,8 +134,7 @@ function read_msh_with_cell_names(path::String, spaceDim = 0; verbose = false)
     _spaceDim = spaceDim > 0 ? spaceDim : _compute_space_dim(verbose)
     mesh = _read_msh(_spaceDim, verbose)
 
-    #----- Read cell names
-    # Read elements
+    # Read volumic physical groups (build a dict tag -> name)
     el_tags = gmsh.model.getPhysicalGroups(_spaceDim)
     _el_names = [gmsh.model.getPhysicalName(_dim, _tag) for (_dim, _tag) in el_tags]
     el = [
@@ -145,34 +144,26 @@ function read_msh_with_cell_names(path::String, spaceDim = 0; verbose = false)
     el_names = Dict(convert(Int, _tag) => _name for (_tag, _name) in el)
     el_names_inv = Dict(_name => convert(Int, _tag) for (_tag, _name) in el)
 
-    #@show _el_names
+    # Read cell indices associated to each volumic physical group
     el_cells = Dict{Int, Array{Int}}()
-    #el_cells_inv = Dict{Int,Int}()
     for (_dim, _tag) in el_tags
         v = Int[]
-        for i in 1:length(gmsh.model.getEntitiesForPhysicalGroup(_dim, _tag))
-            tmpTypes, tmpTags, tmpNodeTags = gmsh.model.mesh.getElements(
-                _dim,
-                gmsh.model.getEntitiesForPhysicalGroup(_dim, _tag)[i],
-            )
-            #@show size(tmpTags[1])
-            #@show tmpTags
-            if length(tmpTags) > 0
-                v = vcat(v, Int[i for i in tmpTags[1]])
+
+        for iEntity in gmsh.model.getEntitiesForPhysicalGroup(_dim, _tag)
+            tmpTypes, tmpTags, tmpNodeTags = gmsh.model.mesh.getElements(_dim, iEntity)
+
+            # Notes : a PhysicalGroup "entity" can contain different types of elements.
+            # So `tmpTags` is an array of the cell indices of each type in the Physical group.
+            for _tmpTags in tmpTags
+                v = vcat(v, Int.(_tmpTags))
             end
         end
         el_cells[_tag] = v
-
-        #for i=1:length(v)
-        #    el_cells_inv[glo2loc_cell_indices[v[i]]] = _tag
-        #end
     end
 
     absolute_cell_indices = absolute_indices(mesh, :cell)
     _, glo2loc_cell_indices = densify(absolute_cell_indices; permute_back = true)
 
-    #el_cells = Dict( convert(Int,_tag) => Int[i for i in gmsh.model.mesh.getElementsByType(spaceDim,_tag)] for (_tag,_name) in el_names)
-
     # free gmsh
     gmsh.finalize()
 

src/mesh/mesh_generator.jl

@@ -814,3 +814,62 @@ function _duplicate_mesh(mesh::AbstractMesh)
         bc_faces = boundary_faces(mesh),
     )
 end
+
+"""
+    _compute_space_dim(topodim, lx, ly, lz, tol, verbose::Bool)
+
+Deduce the number of space dimensions from the mesh boundaries : if one (or more) dimension of the bounding
+box is way lower than the other dimensions, the number of space dimension is decreased.
+
+Currently, having for instance (x,z) is not supported. Only (x), or (x,y), or (x,y,z).
+"""
+function _compute_space_dim(topodim, lx, ly, lz, tol, verbose::Bool)
+
+    # Maximum dimension and default value for number of space dimensions
+    lmax = maximum([lx, ly, lz])
+
+    # Now checking several case (complex `if` cascade for comprehensive warning messages)
+
+    # If the topology is 3D, useless to continue, the space dim must be 3
+    topodim == 3 && (return 3)
+
+    if topodim == 2
+        if (lz / lmax < tol)
+            msg = "Warning : the mesh is flat on the z-axis. It is now considered 2D."
+            msg *= " (use `spacedim` argument of `read_msh` if you want to keep 3D coordinates.)"
+            msg *= " Disable this warning with `verbose = false`"
+            verbose && println(msg)
+
+            return 2
+        else
+            # otherwise, it is a surface in a 3D space, we keep 3 coordinates
+            return 3
+        end
+    end
+
+    if topodim == 1
+        if (ly / lmax < tol && lz / lmax < tol)
+            msg = "Warning : the mesh is flat on the y and z axis. It is now considered 1D."
+            msg *= " (use `spacedim` argument of `read_msh` if you want to keep 2D or 3D coordinates.)"
+            msg *= " Disable this warning with `verbose = false`"
+            verbose && println(msg)
+
+            return 1
+        elseif (ly / lmax < tol && lz / lmax > tol)
+            error(
+                "You have a flat y-axis but a non flat z-axis, this is not supported. Consider rotating your mesh.",
+            )
+
+        elseif (ly / lmax > tol && lz / lmax < tol)
+            msg = "Warning : the mesh is flat on the z-axis. It is now considered as a 1D mesh in a 2D space."
+            msg *= " (use `spacedim` argument of `read_msh` if you want to keep 3D coordinates.)"
+            msg *= " Disable this warning with `verbose = false`"
+            verbose && println(msg)
+
+            return 2
+        else
+            # otherwise, it is a line in a 3D space, we keep 3 coordinates
+            return 3
+        end
+    end
+end

src/utils.jl

@@ -198,60 +198,3 @@ end
 end
 @inline tuplemap(f, ::Tuple{}) = ()
 @inline tuplemap(f, ::Tuple{}, ::Tuple{}) = ()
-
-"""
-Deduce the number of space dimensions from the mesh : if one (or more) dimension of the bounding
-box is way lower than the other dimensions, the number of space dimension is decreased.
-
-Currently, having for instance (x,z) is not supported. Only (x), or (x,y), or (x,y,z).
-"""
-function _compute_space_dim(topodim, lx, ly, lz, tol, verbose::Bool)
-
-    # Maximum dimension and default value for number of space dimensions
-    lmax = maximum([lx, ly, lz])
-
-    # Now checking several case (complex `if` cascade for comprehensive warning messages)
-
-    # If the topology is 3D, useless to continue, the space dim must be 3
-    topodim == 3 && (return 3)
-
-    if topodim == 2
-        if (lz / lmax < tol)
-            msg = "Warning : the mesh is flat on the z-axis. It is now considered 2D."
-            msg *= " (use `spacedim` argument of `read_msh` if you want to keep 3D coordinates.)"
-            msg *= " Disable this warning with `verbose = false`"
-            verbose && println(msg)
-
-            return 2
-        else
-            # otherwise, it is a surface in a 3D space, we keep 3 coordinates
-            return 3
-        end
-    end
-
-    if topodim == 1
-        if (ly / lmax < tol && lz / lmax < tol)
-            msg = "Warning : the mesh is flat on the y and z axis. It is now considered 1D."
-            msg *= " (use `spacedim` argument of `read_msh` if you want to keep 2D or 3D coordinates.)"
-            msg *= " Disable this warning with `verbose = false`"
-            verbose && println(msg)
-
-            return 1
-        elseif (ly / lmax < tol && lz / lmax > tol)
-            error(
-                "You have a flat y-axis but a non flat z-axis, this is not supported. Consider rotating your mesh.",
-            )
-
-        elseif (ly / lmax > tol && lz / lmax < tol)
-            msg = "Warning : the mesh is flat on the z-axis. It is now considered as a 1D mesh in a 2D space."
-            msg *= " (use `spacedim` argument of `read_msh` if you want to keep 3D coordinates.)"
-            msg *= " Disable this warning with `verbose = false`"
-            verbose && println(msg)
-
-            return 2
-        else
-            # otherwise, it is a line in a 3D space, we keep 3 coordinates
-            return 3
-        end
-    end
-end

src/writers/vtk.jl

@@ -321,22 +321,45 @@ function _vtk_coords_from_lagrange(shape::Union{Square, Cube}, degree)
 end
 
 """
+    write_vtk_lagrange(
+        basename::String,
+        vars::Dict{String, F},
+        mesh::AbstractMesh,
+        it::Integer = -1,
+        time::Real = 0.0;
+        mesh_degree::Integer = 1,
+        discontinuous::Bool = true,
+        functionSpaceType::AbstractFunctionSpaceType = Lagrange(),
+        collection_append::Bool = false,
+        vtk_kwargs...,
+    ) where {F <: AbstractLazy}
+
     write_vtk_lagrange(
         basename::String,
         vars::Dict{String, F},
         mesh::AbstractMesh,
         U_export::AbstractFESpace,
         it::Integer = -1,
         time::Real = 0.0;
-        collection_append = false,
+        collection_append::Bool = false,
         vtk_kwargs...,
     ) where {F <: AbstractLazy}
 
-Write the provided FEFunction/MeshCellData/CellFunction on the mesh with
-the precision of the Lagrange FESpace provided.
+Write the provided FEFunction/MeshCellData/CellFunction on a mesh of order `mesh_degree` with the nodes
+defined by the `functionSpaceType` (only `Lagrange{:Uniform}` supported for now). The boolean `discontinuous`
+indicate if the node values should be discontinuous or not.
 
 `vars` is a dictionnary of variable name => Union{FEFunction,MeshCellData,CellFunction} to write.
 
+Cell-centered data should be wrapped as `MeshCellData` in the `vars` dict. To convert an `AbstractLazy` to
+a `MeshCellData`, simply use `Bcube.cell_mean` or `MeshCellData ∘ Bcube.var_on_centers`.
+
+# Remarks:
+* If `mesh` is of degree `d`, the solution will be written on a mesh of degree `mesh_degree`, even if this
+number is different from `d`.
+* The degree of the input FEFunction (P1, P2, P3, ...) is not used to define the nodes where the solution is
+written, only `mesh` and `mesh_degree` matter. The FEFunction is simply evaluated on the aforementionned nodes.
+
 # Example
 ```julia
 mesh = rectangle_mesh(6, 7; xmin = -1, xmax = 1.0, ymin = -1, ymax = 1.0)
@@ -346,13 +369,13 @@ projection_l2!(u, f_u, mesh)
 
 vars = Dict("f_u" => f_u, "u" => u, "grad_u" => ∇(u))
 
-for degree_export in 1:5
-    U_export = TrialFESpace(FunctionSpace(:Lagrange, degree_export), mesh)
+for mesh_degree in 1:5
     Bcube.write_vtk_lagrange(
         joinpath(@__DIR__, "output"),
         vars,
-        mesh,
-        U_export,
+        mesh;
+        mesh_degree,
+        discontinuous = false
     )
 end
 ```
@@ -362,6 +385,35 @@ end
 - `collection_append` is not named `append` to enable passing correct `kwargs` to `vtk_grid`
 - remove (once fully validated) : `write_vtk_discontinuous`
 """
+function write_vtk_lagrange(
+    basename::String,
+    vars::Dict{String, F},
+    mesh::AbstractMesh,
+    it::Integer = -1,
+    time::Real = 0.0;
+    mesh_degree::Integer = 1,
+    discontinuous::Bool = true,
+    functionSpaceType::AbstractFunctionSpaceType = Lagrange(),
+    collection_append::Bool = false,
+    vtk_kwargs...,
+) where {F <: AbstractLazy}
+    U_export = TrialFESpace(
+        FunctionSpace(functionSpaceType, mesh_degree),
+        mesh;
+        isContinuous = !discontinuous,
+    )
+    write_vtk_lagrange(
+        basename,
+        vars,
+        mesh,
+        U_export,
+        it,
+        time;
+        collection_append,
+        vtk_kwargs...,
+    )
+end
+
 function write_vtk_lagrange(
     basename::String,
     vars::Dict{String, F},
@@ -379,21 +431,6 @@ function write_vtk_lagrange(
     dhl_export = Bcube._get_dhl(U_export)
     nd = get_ndofs(dhl_export)
 
-    # If the FunctionSpace is a Lagrange P0, we shortcut the rest of this function
-    # to call a more adequate one.
-    if degree_export == 0
-        vals = map(values(vars)) do var
-            _, dim = get_return_type_and_codim(var, mesh)
-            is_scalar = dim[1] == 1
-            value = var_on_centers(var, mesh)
-            value = is_scalar ? value : transpose(value)
-            value
-        end
-        d = Dict(key => (value, VTKCellData()) for (key, value) in zip(keys(vars), vals))
-        write_vtk(basename, it, time, mesh, d; append = collection_append)
-        return nothing
-    end
-
     # extract all `MeshCellData` in `vars` as this type of variable
     # will not be interpolated to nodes and will be written with
     # the `VTKCellData` attribute.

test/checksums.sha1

@@ -1,4 +1,3 @@
-54104f8915cf522913bbac20a8a670e7e7a36e50  write_vtk_lagrange_deg0.vtu
 bc4229dd75d992f8788a76769fe26d7fffef227c  write_vtk_lagrange_deg1.vtu
 6333082e04dc024aed37545bb50e513b2dfd836d  write_vtk_lagrange_deg2.vtu
 afc0e781b9170934c0a95525291e7ec26e01d7ed  write_vtk_lagrange_deg4.vtu