From 196c76c55d3caf872a97c559693bac967d73a278 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 13 Oct 2023 15:06:08 +0200 Subject: [PATCH 1/3] improving `makie` support --- docs/src/geo_and_meshes.md | 18 +++++--- ext/IntiMakieExt.jl | 87 ++------------------------------------ src/domain.jl | 2 + src/mesh.jl | 67 ++++++++++++++++++++++++++++- 4 files changed, 84 insertions(+), 90 deletions(-) diff --git a/docs/src/geo_and_meshes.md b/docs/src/geo_and_meshes.md index 933e6957..bcc9090d 100644 --- a/docs/src/geo_and_meshes.md +++ b/docs/src/geo_and_meshes.md @@ -11,8 +11,10 @@ using Gmsh gmsh.initialize() gmsh.option.setNumber("General.Verbosity", 2) gmsh.model.add("Sphere") +# set mesh size gmsh.model.occ.addSphere(0,0,0,1) gmsh.model.occ.synchronize() +gmsh.option.setNumber("Mesh.MeshSizeMax", 0.2) gmsh.model.mesh.generate(3) ents = gmsh.model.getEntities() Ω = Inti.gmsh_import_domain(;dim=3) @@ -24,16 +26,17 @@ gmsh.finalize() - You can plot a `msh` if you have a `Makie` backend (e.g. `GLMakie`) ```@example gmsh-sphere -using CairoMakie # or GLMakie +using CairoMakie # or GLMakie for interactivity if supported Γ = Inti.external_boundary(Ω) -poly(view(msh,Γ);strokewidth=1,color=:lightgray, transparency=true) +color = [cos(10*x[1]) for x in Inti.nodes(msh)] +poly(view(msh,Γ);strokewidth=0.2,color, transparency=true) ``` You can also plot the volume mesh (see the Documentation of `Makie.poly` for more details on possible arguments): ```@example gmsh-sphere -poly(view(msh,Ω);strokewidth=1,color=:lightgray, transparency=true) +poly(view(msh,Ω);color,transparency=true,strokecolor=:lightgray, alpha = 0.1) ``` Two-dimensional meshes are very similar: @@ -44,12 +47,17 @@ Two-dimensional meshes are very similar: using CairoMakie gmsh.initialize() gmsh.option.setNumber("General.Verbosity", 2) + gmsh.option.setNumber("Mesh.MeshSizeMax", 0.1) gmsh.model.add("Disk") - gmsh.model.occ.addDisk(0,0,0,1,1) + gmsh.model.occ.addDisk(0,0,0,3,1) gmsh.model.occ.synchronize() gmsh.model.mesh.generate(2) Ω = Inti.gmsh_import_domain(;dim=2) msh = Inti.gmsh_import_mesh(Ω;dim=2) gmsh.finalize() - poly(view(msh,Ω);strokewidth=2) + color = [cos(20*x[1]) for x in Inti.nodes(msh)] + fig,ax,p = poly(view(msh,Ω);strokewidth=1,color) + colsize!(fig.layout, 1, Aspect(1, 3)) + resize_to_layout!(fig) + fig ``` diff --git a/ext/IntiMakieExt.jl b/ext/IntiMakieExt.jl index d2b6e231..bc628458 100644 --- a/ext/IntiMakieExt.jl +++ b/ext/IntiMakieExt.jl @@ -27,96 +27,15 @@ function tomakie_dim1(msh::Inti.AbstractMesh{N,T}) where {N,T} return coords end -function tomakie_dim2(msh::Inti.AbstractMesh{N,T}) where {N,T} - coords = Makie.Point{N,T}[] - connec = Int[] - for E in Inti.element_types(msh) - iter = Inti.elements(msh, E) - D = Inti.domain(E) - if D isa Inti.ReferenceTriangle - vtxs = Inti.vertices(D) - for el in iter - for vtx in vtxs - push!(coords, el(vtx)) - push!(connec, length(coords)) - end - end - elseif D isa ReferenceSquare - # split square in two triangles for visualization - vtxs_down = Inti.vertices(Inti.ReferenceTriangle()) - vtxs_up = map(v -> -v .+ 1, vtx_down) - for el in iter - for vtxs in (vtxs_down, vtxs_up) # the two triangles - for vtx in vtxs - push!(coords, el(vtx)) - push!(connec, length(coords)) - end - end - end - end - end - return coords, connec -end - -function tomakie_dim3(msh::Inti.AbstractMesh{N,T}) where {N,T} - coords = Makie.Point{N,T}[] - connec = Int[] - for E in Inti.element_types(msh) - iter = Inti.elements(msh, E) - D = Inti.domain(E) - @assert D isa Inti.ReferenceTetrahedron - vtxs = Inti.vertices(D) - for el in iter - for nf in 1:4 # four faces - for (i, vtx) in enumerate(vtxs) - i == nf && continue # i-th face exclude the i-th vertex - push!(coords, el(vtx)) - push!(connec, length(coords)) - end - end - end - end - return coords, connec -end - -function Makie.convert_arguments(P::Type{<:Makie.Lines}, msh::Inti.AbstractMesh) +function Makie.convert_arguments(::Type{<:Makie.Lines}, msh::Inti.AbstractMesh) @assert Inti.geometric_dimension(msh) == 1 "Lines only supported for meshes of geometric dimension 1" coords = tomakie_dim1(msh) return (coords,) end function Makie.convert_arguments(P::Type{<:Makie.Poly}, msh::Inti.AbstractMesh) - gdim = Inti.geometric_dimension(msh) - if gdim == 2 - coords, connec = tomakie_dim2(msh) - elseif gdim == 3 - coords, connec = tomakie_dim3(msh) - else - error("Poly only supported for meshes of geometric dimension 2 or 3") - end - return Makie.convert_arguments(P, coords, connec) -end - -function Makie.convert_arguments( - P::Type{<:Makie.Arrows}, - msh::Inti.AbstractMesh{N,T}, -) where {N,T} - gdim = Inti.geometric_dimension(msh) - adim = Inti.ambient_dimension(msh) - codim = adim - gdim - @assert codim == 1 "Arrows only supported for meshes of codimension 1" - coords = Makie.Point{N,T}[] - normals = Makie.Point{N,T}[] - for E in Inti.element_types(msh) - iter = Inti.elements(msh, E) - dom = Inti.domain(E) - xc = Inti.center(dom) - for el in iter - push!(coords, el(xc)) - push!(normals, normal(el, xc)) - end - end - return Makie.convert_arguments(P, coords, normals) + connec = Inti.triangle_connectivity(msh) + return Makie.convert_arguments(P,Inti.nodes(msh), connec) end end # module diff --git a/src/domain.jl b/src/domain.jl index 562b3201..670035e2 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -93,3 +93,5 @@ Iterating over a domain means iterating over its entities. Base.iterate(Ω::Domain, state = 1) = iterate(entities(Ω), state) Base.isempty(Ω::Domain) = isempty(entities(Ω)) + +Base.in(ent::AbstractEntity, Ω::Domain) = in(ent, entities(Ω)) diff --git a/src/mesh.jl b/src/mesh.jl index 35043e4a..ce30dab2 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -84,8 +84,17 @@ nodes(msh::LagrangeMesh) = msh.nodes elements(msh::LagrangeMesh) = msh.etype2mat ent2tags(msh::LagrangeMesh) = msh.ent2tags +entities(msh::LagrangeMesh) = keys(msh.ent2tags) + +""" + domain(msh::LagrangeMesh) + +Set of all entities covered by the mesh. """ - dom2elt(m::LagrangeMesh,Ω,E) +domain(msh::LagrangeMesh) = Domain(entities(msh)) + +""" + dom2elt(m::LagrangeMesh,Ω,E)::Vector{Int} Compute the element indices `idxs` of the elements of type `E` composing `Ω`, so that `elements(m)[idxs]` gives all the elements of type `E` meshing `Ω`. @@ -181,6 +190,9 @@ ambient_dimension(::SubMesh{N}) where {N} = N geometric_dimension(msh::SubMesh) = geometric_dimension(msh.domain) +nodes(msh::SubMesh) = nodes(msh.parent) +domain(msh::SubMesh) = msh.domain + element_types(msh::SubMesh) = keys(msh.etype2etags) # ElementIterator for submesh @@ -206,3 +218,56 @@ function Base.iterate(iter::ElementIterator{<:LagrangeElement,<:SubMesh}, state state > length(iter) && (return nothing) return iter[state], state + 1 end + + +triangle_connectivity(msh::SubMesh) = _triangle_connectivity(msh.parent, msh.domain) +triangle_connectivity(msh::LagrangeMesh) = _triangle_connectivity(msh, domain(msh)) +function _triangle_connectivity(msh::Inti.LagrangeMesh{N,T},Ω::Inti.Domain) where {N,T} + connec = Int[] + for E in Inti.element_types(msh) + el_idxs = Inti.dom2elt(msh, Ω, E)::Vector{Int} + isempty(el_idxs) && continue + tags = msh.etype2mat[E]::Matrix{Int} + if E <: Inti.LagrangeTriangle + # extract the first three tags + for n in el_idxs + push!(connec,tags[1,n]) + push!(connec,tags[2,n]) + push!(connec,tags[3,n]) + end + elseif E <: Inti.LagrangeSquare + for n in el_idxs + # lower triangle + push!(connec,tags[1,n]) + push!(connec,tags[2,n]) + push!(connec,tags[3,n]) + # upper triangle + push!(connec,tags[3,n]) + push!(connec,tags[4,n]) + push!(connec,tags[1,n]) + end + elseif E <: Inti.LagrangeTetrahedron + for n in el_idxs + # four faces + push!(connec,tags[1,n]) + push!(connec,tags[2,n]) + push!(connec,tags[3,n]) + # + push!(connec,tags[1,n]) + push!(connec,tags[2,n]) + push!(connec,tags[4,n]) + # + push!(connec,tags[1,n]) + push!(connec,tags[3,n]) + push!(connec,tags[4,n]) + # + push!(connec,tags[2,n]) + push!(connec,tags[3,n]) + push!(connec,tags[4,n]) + end + else + error("element type $E not supported") + end + end + return connec +end From e7631d4d49c07e387bbebe83ef26efaf358a3f63 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Wed, 18 Oct 2023 14:40:23 +0200 Subject: [PATCH 2/3] add docstrings to `triangle_connectivity` --- src/mesh.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/mesh.jl b/src/mesh.jl index ce30dab2..b0fe070c 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -219,7 +219,18 @@ function Base.iterate(iter::ElementIterator{<:LagrangeElement,<:SubMesh}, state return iter[state], state + 1 end +""" + triangle_connectivity(msh::AbstractMesh) + +Return a vector of integers `connec` providing a triangulation of the elements +in the `mesh`. The `connec` vector should be read in groups of three to get the +triangles. +For elements of `LagrangeTriangle` type, the first three nodes are returned. For +elements of `LagrangeSquare` type, the triangulation is done by splitting the +element in two triangles. For elements of `LagrangeTetrahedron` type, the four +faces are returned. This method will error if the element type is not supported. +""" triangle_connectivity(msh::SubMesh) = _triangle_connectivity(msh.parent, msh.domain) triangle_connectivity(msh::LagrangeMesh) = _triangle_connectivity(msh, domain(msh)) function _triangle_connectivity(msh::Inti.LagrangeMesh{N,T},Ω::Inti.Domain) where {N,T} From ac048dd1ff78228a47ea950efe570523611bdcba Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Wed, 18 Oct 2023 14:41:28 +0200 Subject: [PATCH 3/3] autoformat --- ext/IntiMakieExt.jl | 2 +- src/mesh.jl | 48 ++++++++++++++++++++++----------------------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/ext/IntiMakieExt.jl b/ext/IntiMakieExt.jl index 15d45ab3..b4b48b77 100644 --- a/ext/IntiMakieExt.jl +++ b/ext/IntiMakieExt.jl @@ -35,7 +35,7 @@ end function Makie.convert_arguments(P::Type{<:Makie.Poly}, msh::Inti.AbstractMesh) connec = Inti.triangle_connectivity(msh) - return Makie.convert_arguments(P,Inti.nodes(msh), connec) + return Makie.convert_arguments(P, Inti.nodes(msh), connec) end end # module diff --git a/src/mesh.jl b/src/mesh.jl index b0fe070c..7907d38c 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -91,7 +91,7 @@ entities(msh::LagrangeMesh) = keys(msh.ent2tags) Set of all entities covered by the mesh. """ -domain(msh::LagrangeMesh) = Domain(entities(msh)) +domain(msh::LagrangeMesh) = Domain(entities(msh)) """ dom2elt(m::LagrangeMesh,Ω,E)::Vector{Int} @@ -231,9 +231,9 @@ elements of `LagrangeSquare` type, the triangulation is done by splitting the element in two triangles. For elements of `LagrangeTetrahedron` type, the four faces are returned. This method will error if the element type is not supported. """ -triangle_connectivity(msh::SubMesh) = _triangle_connectivity(msh.parent, msh.domain) +triangle_connectivity(msh::SubMesh) = _triangle_connectivity(msh.parent, msh.domain) triangle_connectivity(msh::LagrangeMesh) = _triangle_connectivity(msh, domain(msh)) -function _triangle_connectivity(msh::Inti.LagrangeMesh{N,T},Ω::Inti.Domain) where {N,T} +function _triangle_connectivity(msh::Inti.LagrangeMesh{N,T}, Ω::Inti.Domain) where {N,T} connec = Int[] for E in Inti.element_types(msh) el_idxs = Inti.dom2elt(msh, Ω, E)::Vector{Int} @@ -242,39 +242,39 @@ function _triangle_connectivity(msh::Inti.LagrangeMesh{N,T},Ω::Inti.Domain) whe if E <: Inti.LagrangeTriangle # extract the first three tags for n in el_idxs - push!(connec,tags[1,n]) - push!(connec,tags[2,n]) - push!(connec,tags[3,n]) + push!(connec, tags[1, n]) + push!(connec, tags[2, n]) + push!(connec, tags[3, n]) end elseif E <: Inti.LagrangeSquare for n in el_idxs # lower triangle - push!(connec,tags[1,n]) - push!(connec,tags[2,n]) - push!(connec,tags[3,n]) + push!(connec, tags[1, n]) + push!(connec, tags[2, n]) + push!(connec, tags[3, n]) # upper triangle - push!(connec,tags[3,n]) - push!(connec,tags[4,n]) - push!(connec,tags[1,n]) + push!(connec, tags[3, n]) + push!(connec, tags[4, n]) + push!(connec, tags[1, n]) end elseif E <: Inti.LagrangeTetrahedron for n in el_idxs # four faces - push!(connec,tags[1,n]) - push!(connec,tags[2,n]) - push!(connec,tags[3,n]) + push!(connec, tags[1, n]) + push!(connec, tags[2, n]) + push!(connec, tags[3, n]) # - push!(connec,tags[1,n]) - push!(connec,tags[2,n]) - push!(connec,tags[4,n]) + push!(connec, tags[1, n]) + push!(connec, tags[2, n]) + push!(connec, tags[4, n]) # - push!(connec,tags[1,n]) - push!(connec,tags[3,n]) - push!(connec,tags[4,n]) + push!(connec, tags[1, n]) + push!(connec, tags[3, n]) + push!(connec, tags[4, n]) # - push!(connec,tags[2,n]) - push!(connec,tags[3,n]) - push!(connec,tags[4,n]) + push!(connec, tags[2, n]) + push!(connec, tags[3, n]) + push!(connec, tags[4, n]) end else error("element type $E not supported")