|
| 1 | +```@meta |
| 2 | +CurrentModule = MultiGridBarrier |
| 3 | +``` |
| 4 | + |
| 5 | +# Meshes from Gmsh |
| 6 | + |
| 7 | +MultiGridBarrier's own constructors (`fem2d_P1`, `fem2d`, `fem3d`, …) build |
| 8 | +meshes of simple domains. For real geometry — CAD shapes, holes, local |
| 9 | +refinement, named boundary parts — use [Gmsh](https://gmsh.info), the standard |
| 10 | +open-source mesh generator. The `MultiGridBarrierGmshExt` extension loads |
| 11 | +automatically when both packages are imported and provides |
| 12 | +[`gmsh_import`](@ref), which converts the current Gmsh mesh (or a `.msh` file) |
| 13 | +into a `Geometry` plus named node sets: |
| 14 | + |
| 15 | +- 3/6-node triangles → `fem2d_P1` / `fem2d_P2` (curved P2 edges supported), |
| 16 | +- 4/9-node quadrilaterals → tensor `fem2d` with `k = 1/2` (curved quads |
| 17 | + supported; non-planar quad meshes become embedded surfaces), |
| 18 | +- 8/27-node hexahedra → `fem3d` with `k = 1/2`. |
| 19 | + |
| 20 | +Gmsh **physical groups** come back as `(vertex, element)` node-pair lists — the |
| 21 | +same format as [`find_boundary`](@ref) — so named boundary parts plug directly |
| 22 | +into `amg`'s `dirichlet_nodes`, and named subdomains into the JuMP front end's |
| 23 | +[`On`](@ref). |
| 24 | + |
| 25 | +!!! note "Requires Gmsh" |
| 26 | + Add the [`Gmsh`](https://github.com/JuliaFEM/Gmsh.jl) package |
| 27 | + (`pkg> add Gmsh`) and load it (`using Gmsh: gmsh`). Gmsh itself is |
| 28 | + GPL-licensed; as an opt-in weak dependency it is only ever loaded if you |
| 29 | + load it. |
| 30 | + |
| 31 | +## Example: square with a hole, mixed boundary conditions |
| 32 | + |
| 33 | +Script the geometry through the `gmsh` API (or `gmsh.open` a `.geo`/`.msh` |
| 34 | +file), name the boundary parts, mesh, and import: |
| 35 | + |
| 36 | +```@example gmsh |
| 37 | +using MultiGridBarrier, PyPlot |
| 38 | +using Gmsh: gmsh |
| 39 | +
|
| 40 | +gmsh.initialize() |
| 41 | +gmsh.option.setNumber("General.Terminal", 0) |
| 42 | +sq = gmsh.model.occ.addRectangle(-1.0, -1.0, 0.0, 2.0, 2.0) |
| 43 | +hole = gmsh.model.occ.addDisk(0.3, 0.2, 0.0, 0.4, 0.4) |
| 44 | +gmsh.model.occ.cut([(2, sq)], [(2, hole)]) |
| 45 | +gmsh.model.occ.synchronize() |
| 46 | +# name the two boundary parts: the circle and the outer square |
| 47 | +circle = Int[]; outer = Int[] |
| 48 | +for (d, t) in gmsh.model.getEntities(1) |
| 49 | + xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(d, t) |
| 50 | + if xmin > -0.9 && xmax < 0.9 && ymin > -0.9 && ymax < 0.9 |
| 51 | + push!(circle, t) |
| 52 | + else |
| 53 | + push!(outer, t) |
| 54 | + end |
| 55 | +end |
| 56 | +gmsh.model.addPhysicalGroup(1, outer, -1, "outer") |
| 57 | +gmsh.model.addPhysicalGroup(1, circle, -1, "hole") |
| 58 | +gmsh.option.setNumber("Mesh.MeshSizeMax", 0.12) |
| 59 | +gmsh.model.mesh.generate(2) |
| 60 | +gm = gmsh_import() |
| 61 | +gmsh.finalize() |
| 62 | +length(gm.geometry.w) # number of quadrature nodes |
| 63 | +``` |
| 64 | + |
| 65 | +The named parts drive per-region Dirichlet conditions exactly like hand-built |
| 66 | +node sets: here `u = 0` on the outer square and `u = 1` on the hole, with the |
| 67 | +p-Laplace energy (`p = 1.5`): |
| 68 | + |
| 69 | +```@example gmsh |
| 70 | +mg = amg(gm.geometry; dirichlet_nodes = Dict(:dirichlet => sort(vcat(gm.regions["outer"], gm.regions["hole"])))) |
| 71 | +onhole = Set(gm.regions["hole"]) |
| 72 | +V, N = size(gm.geometry.x, 1), size(gm.geometry.x, 2) |
| 73 | +gvals = zeros(V * N) |
| 74 | +for (v, e) in gm.regions["hole"] |
| 75 | + gvals[v + (e - 1) * V] = 1.0 |
| 76 | +end |
| 77 | +xf = reshape(gm.geometry.x, :, 2) |
| 78 | +sol = mgb_solve(assemble(mg; p = 1.5, |
| 79 | + g_grid = [gvals fill(100.0, V * N)], |
| 80 | + f = x -> (0.0, 0.0, 0.0, 1.0)); verbose = false) |
| 81 | +plot(sol); savefig("gmsh_hole.svg"); nothing # hide |
| 82 | +close() # hide |
| 83 | +``` |
| 84 | + |
| 85 | + |
| 86 | +## Curved elements |
| 87 | + |
| 88 | +Second-order meshes import as isoparametric elements: `setOrder(2)` before |
| 89 | +importing gives curved 6-node triangles (`fem2d_P2`) or curved 9-node quads |
| 90 | +(tensor `fem2d`, `k = 2`) whose boundary nodes lie on the true geometry. For |
| 91 | +all-quad meshes, set `Mesh.RecombineAll = 1` (with |
| 92 | +`Mesh.SubdivisionAlgorithm = 1` to guarantee no leftover triangles); for |
| 93 | +hexahedra use transfinite/swept volumes, or subdivide a tet mesh with |
| 94 | +`Mesh.SubdivisionAlgorithm = 2`. Orders ≥ 3, serendipity elements, and |
| 95 | +tetrahedra are rejected with actionable error messages. |
| 96 | + |
| 97 | +```@example gmsh |
| 98 | +gmsh.initialize() |
| 99 | +gmsh.option.setNumber("General.Terminal", 0) # hide |
| 100 | +gmsh.model.occ.addDisk(0.0, 0.0, 0.0, 1.0, 1.0) |
| 101 | +gmsh.model.occ.synchronize() |
| 102 | +gmsh.model.addPhysicalGroup(1, [t for (d, t) in gmsh.model.getEntities(1)], -1, "circle") |
| 103 | +gmsh.option.setNumber("Mesh.MeshSizeMax", 0.25) |
| 104 | +gmsh.option.setNumber("Mesh.RecombineAll", 1) |
| 105 | +gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) |
| 106 | +gmsh.model.mesh.generate(2) |
| 107 | +gmsh.model.mesh.setOrder(2) # curved Q2 quads |
| 108 | +gmd = gmsh_import() |
| 109 | +gmsh.finalize() |
| 110 | +sum(gmd.geometry.w) - π # curved quadrature: area error of the disk |
| 111 | +``` |
| 112 | + |
| 113 | +## API reference |
| 114 | + |
| 115 | +```@docs |
| 116 | +gmsh_import |
| 117 | +``` |
0 commit comments