coords = [0.0 1.0 2.0 0.1 0.9 1.9; 0.0 0.1 0.0 0.9 1.0 0.9]
elts = [[1, 2, 5, 4], [2, 3, 6], [2, 6, 5]]
m = Mesh(coords, elts, 2)
mplot(m) |> mconf()
2 The Mesh interface
2.1 Making meshes
2.1.1 Mesh by coordinates and connectivity.
2.1.2 Build mesh incrementally
l = 80
h1 = 2
h2 = 12
nn = 21
ne = nn - 1
ne2 = ne ÷ 2
m = Mesh(1, 2)
is1 = addnodes!(m, [-l / 2, 0], [l / 2, 0], nn)
is2 = addnodes!(m, polynomialinterpolation([-l / 2, h1], [0, h2], [l / 2, h1]), nn)
addelements!(m, is1, is1 .+ 1, n=ne, group=:chord)
addelements!(m, is2, is2 .+ 1, n=ne, group=:chord)
addelements!(m, is1, is2, group=:vert)
addelements!(m, is2, is1 .+ 1, n=ne2, group=:diag)
addelements!(m, is1 .+ ne2, is2 .+ (ne2 + 1), n=ne2, group=:diag)
mplot(m, edgecolor=:groups, nodesvisible=true) |> mconf()
2.1.3 Read mesh from Gmsh file
Specify meshfile generated by Gmsh:
m = Mesh(meshpath("complex-g1.msh"))
mplot(m, edgesvisible=true) |> mconf()
If an Gmsh input file (extension .geo) is given, the meshfile is generated automatically:
m = Mesh(meshpath("complex-g1.geo"))
mplot(m, edgesvisible=true) |> mconf()
As an alternative, one can pass the Gmsh input as string:
m = Mesh(
"""
Point(1) = {0, 0, 0, 0};
Point(2) = {1, 0, 0, 0};
Point(3) = {0.75, 1, 0, 0};
Point(4) = {0.25, 1, 0, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
MeshSize{:} = 0.25;
Mesh 2;
"""
)
mplot(m, edgesvisible=true) |> mconf()
2.2 Working with meshes
A mesh consists of mesh entities:
| Entity Type | Parametric dimension |
|---|---|
| Node | 0 |
| Edge | 1 |
| Face | 2 |
| Solid | 3 |
| Element | Maximum dimension in mesh |
The Element type simply refers to the parts of the highest dimension present in the mesh. It is introduced in order to be compatible with common finite element terminology.
For the nodes of mesh or mesh entity m, the functions
nnodes(m): Number of nodesnode(m, n): Nodennodes(m): All nodesnodeindices(m): Indices of nodes
exist. There are similar functions for edges, faces, solids and elements. In addition, the functions length(edge) and area(face) exist.
Coordinates of the mesh or mesh entities are accessed using
coordinates(m): All node coordinates in a matrixcoordinates(m, n): Coordinates of nodencoordinates(m, indices): Coordinates of nodes withindices
Here is an example
coords = [0.0 1.0 2.0 0.1 0.9 1.9; 0.0 0.1 0.0 0.9 1.0 0.9]
elts = [[1, 2, 5, 4], [2, 3, 6], [2, 6, 5]]
m = Mesh(coords, elts, 2)
e = element(m, 2)
println("Total node count: ", nnodes(m))
println(" Total area: ", sum(area.(faces(m))))
println(" Node 3: ", node(m, 3))
println(" Edge 8: ", edge(m, 8))
println(" Element: ", e)
println(" Edge of element: ", edge(e, 2))
println(" Node indices: ", nodeindices(e))
println(" Coordinates: ", coordinates(e))
println(" Coordinates: ", coordinates(e, 3))Total node count: 6
Total area: 1.7
Node 3: MMJMesh.Meshes.Node(3)[2.0, 0.0]
Edge 8: MMJMesh.Meshes.Edge(8)[6, 5]
Element: MMJMesh.Meshes.Face(2)[2, 3, 6]
Edge of element: MMJMesh.Meshes.Edge(6)[3, 6]
Node indices: [2, 3, 6]
Coordinates: [1.0 2.0 1.9; 0.1 0.0 0.9]
Coordinates: [1.9, 0.9]A loop can be used to process specific entity types
for e = faces(m)
println(e)
endMMJMesh.Meshes.Face(1)[1, 2, 5, 4]
MMJMesh.Meshes.Face(2)[2, 3, 6]
MMJMesh.Meshes.Face(3)[2, 6, 5]or all entities associated with the mesh
for e = entities(m)
println(e)
endMMJMesh.Meshes.Node(1)[0.0, 0.0]
MMJMesh.Meshes.Node(2)[1.0, 0.1]
MMJMesh.Meshes.Node(3)[2.0, 0.0]
MMJMesh.Meshes.Node(4)[0.1, 0.9]
MMJMesh.Meshes.Node(5)[0.9, 1.0]
MMJMesh.Meshes.Node(6)[1.9, 0.9]
MMJMesh.Meshes.Edge(1)[1, 2]
MMJMesh.Meshes.Edge(2)[2, 5]
MMJMesh.Meshes.Edge(3)[4, 5]
MMJMesh.Meshes.Edge(4)[1, 4]
MMJMesh.Meshes.Edge(5)[2, 3]
MMJMesh.Meshes.Edge(6)[3, 6]
MMJMesh.Meshes.Edge(7)[2, 6]
MMJMesh.Meshes.Edge(8)[6, 5]
MMJMesh.Meshes.Face(1)[1, 2, 5, 4]
MMJMesh.Meshes.Face(2)[2, 3, 6]
MMJMesh.Meshes.Face(3)[2, 6, 5]or a mesh entity
for e = entities(face(m, 3))
println(e)
endMMJMesh.Meshes.Node(2)[1.0, 0.1]
MMJMesh.Meshes.Node(6)[1.9, 0.9]
MMJMesh.Meshes.Node(5)[0.9, 1.0]
MMJMesh.Meshes.Edge(7)[2, 6]
MMJMesh.Meshes.Edge(8)[6, 5]
MMJMesh.Meshes.Edge(2)[2, 5]
MMJMesh.Meshes.Face(1)[1, 2, 5, 4]
MMJMesh.Meshes.Face(2)[2, 3, 6]Nodes can be selected based on a predicate.
m = Mesh(meshpath("complex-g1.msh"))
n1 = nodeindices(m, on(Point(0.75, 1)))
n2 = nodeindices(m, on(Segment([0.75, 0], [2, 0])))
n3 = nodeindices(m, on(VLine(0.5)))
p = mplot(m) |> mconf()
scatter!(coordinates(m, n1), color=:green)
scatter!(coordinates(m, n2), color=:hotpink)
scatter!(coordinates(m, n3), color=:orange)
p
In addition, a single can be identified as well.
nodeindex(m, on(Point(0.75, 1)))132.3 Collections of Entities: Groups
Mesh entities can be collected in groups some of which are predefined.
We take a look at a mesh generated by Gmsh. The input file complex-g1.geo defines the physical curve c1 andthe three physical surfaces s1, s2, s3. When plotted, these are visualized by colors.
m = Mesh(meshpath("complex-g1.msh"))
mplot(m) |> mconf()
In the MMJMesh mesh, groups are represented by symbols. In addition to the group :c1 of edges, MMJMesh defines the group :c10 of nodes on the corresponding curve by appending a 0 to the name.
groupnames(m)5-element Vector{Symbol}:
:c1
:c10
:s1
:s2
:s3Furthermore, some groups are predefined by MMJMesh.
groupnames(m, predefined=true)13-element Vector{Symbol}:
:boundaryedges
:boundaryfaces
:boundarynodes
:c1
:c10
:edges
:elements
:faces
:nodes
:s1
:s2
:s3
:solidsA group collects indices of entities in the group
show(group(m, :boundarynodes))MMJMesh.Meshes.NodeGroup[1:19, 31:34, 46:106, 116:157]A group behaves like a vector of integer values.
g = group(m, :boundarynodes)
g[17:23]7-element Vector{Integer}:
17
18
19
31
32
33
34Groups can be used to filter access to coordinates.
m = Mesh(meshpath("complex-g1.msh"))
p = mplot(m) |> mconf()
scatter!(coordinates(m, group(m, :boundarynodes)), color=:hotpink)
p
Additional groups are defined by specifying a dimension, a name and indices.
definegroup!(m, 2, :f1, 100:300)
definegroup!(m, 2, :f2, 200:600)
definegroup!(m, 1, :e1, 600:800)
definegroup!(m, 1, :e2, 700:900)
mplot(m) |> mconf()
2.4 Information associated with Entities: Data
Example how to attach data to mesh entities:
m = Mesh(meshpath("complex-g1.msh"))
n1 = node(m, 1)
n2 = node(m, 2)
e22 = element(m, 22)
setdata!(n1, :d1, 42)
setdata!(e22, :d2, 61)
println("data( n1, :d1) = ", data(n1, :d1))
println("data( n2, :d1) = ", data(n2, :d1))
println("data(e22, :d2) = ", data(e22, :d2))data( n1, :d1) = 42
data( n2, :d1) = nothing
data(e22, :d2) = 61Use hasdata to test, if data is defined:
println("hasdata( n1, :d1) = ", hasdata(n1, :d1))
println("hasdata( n2, :d1) = ", hasdata(n2, :d1))hasdata( n1, :d1) = true
hasdata( n2, :d1) = falseExample how to attach data to groups:
m = Mesh(meshpath("complex-g1.msh"))
g1 = group(m, :s1)
g2 = group(m, :s2)
g3 = group(m, :s3)
println("s1: ", g1)
println("s2: ", g2)
println("s3: ", g3)
setdata!(g1, :foo, "f1")
setdata!(g2, :foo, "f2")
setdata!(g3, :foo, "f3")
e2 = element(m, 5)
e490 = element(m, 490)
e1132 = element(m, 1132)
println("data(e2, :foo) = ", data(e2, :foo))
println("data(e490, :foo) = ", data(e490, :foo))
println("data(e1132, :foo) = ", data(e1132, :foo))s1: MMJMesh.Meshes.FaceGroup[1:488]
s2: MMJMesh.Meshes.FaceGroup[489:1130]
s3: MMJMesh.Meshes.FaceGroup[1131:1374]
data(e2, :foo) = f1
data(e490, :foo) = f2
data(e1132, :foo) = f3Example how to attach data to elements:
m = Mesh(0 .. 10, 20)
setdata!(group(m, :elements), :foo, 99)
println(data(element(m, 2), :foo))99Set data for nodes and access it on an element level
using Statistics # For mean - function
m = Mesh(:quadtri)
setdata!(m, :a, 100 * (1:nnodes(m)))
println(" Value for node 3: ", data(node(m, 3), :a))
println(" Nodes of element 1: ", nodeindices(element(m, 1)))
println("Values of element 1: ", data(element(m, 1), :a))
println(" Mean for element 1: ", mean(data(element(m, 1), :a))) Value for node 3: 300
Nodes of element 1: [1, 2, 5, 4]
Values of element 1: [100, 200, 500, 400]
Mean for element 1: 300.0