x2.jl

# -----------------------------------------------------------------------------
#
#  Gmsh Julia extended tutorial 2
#
#  Mesh import, discrete entities, hybrid models, terrain meshing
#
# -----------------------------------------------------------------------------

import gmsh

# The API can be used to import a mesh without reading it from a file, by
# creating nodes and elements on the fly and storing them in model
# entities. These model entities can be existing CAD entities, or can be
# discrete entities, entirely defined by the mesh.
#
# Discrete entities can be reparametrized (see `t13.jl') so that they can be
# remeshed later on; and they can also be combined with built-in CAD entities to
# produce hybrid models.
#
# We combine all these features in this tutorial to perform terrain meshing,
# where the terrain is described by a discrete surface (that we then
# reparametrize) combined with a CAD representation of the underground.

gmsh.initialize()

gmsh.model.add("x2")

# We will create the terrain surface mesh from N x N input data points:
N = 100


# Helper function to return a node tag given two indices i and j:
function tag(i, j)
    return (N + 1) * i + j + 1
end

# The x, y, z coordinates of all the nodes:
coords = []

# The tags of the corresponding nodes:
nodes = []

# The connectivities of the triangle elements (3 node tags per triangle) on the
# terrain surface:
tris = []

# The connectivities of the line elements on the 4 boundaries (2 node tags
# for each line element):
lin = [[], [], [], []]

# The connectivities of the point elements on the 4 corners (1 node tag for each
# point element):
pnt = [tag(0, 0), tag(N, 0), tag(N, N), tag(0, N)]

for i in 0:N
    for j in 0:N
        push!(nodes, tag(i, j))
        push!(coords, float(i) / N, float(j) / N,
              0.05 * sin(10 * float(i + j) / N))
        if i > 0 && j > 0
            push!(tris, tag(i - 1, j - 1), tag(i, j - 1), tag(i - 1, j))
            push!(tris, tag(i, j - 1), tag(i, j), tag(i - 1, j))
        end
        if (i == 0 || i == N) && j > 0
            push!(lin[i == 0 ? 4 : 2], tag(i, j - 1), tag(i, j))
        end
        if (j == 0 || j == N) && i > 0
            push!(lin[j == 0 ? 1 : 3], tag(i - 1, j), tag(i, j))
        end
    end
end

# Create 4 discrete points for the 4 corners of the terrain surface:
for i in 0:3
    gmsh.model.addDiscreteEntity(0, i + 1)
end
gmsh.model.setCoordinates(1, 0, 0, coords[3 * tag(0, 0)])
gmsh.model.setCoordinates(2, 1, 0, coords[3 * tag(N, 0)])
gmsh.model.setCoordinates(3, 1, 1, coords[3 * tag(N, N)])
gmsh.model.setCoordinates(4, 0, 1, coords[3 * tag(0, N)])

# Create 4 discrete bounding curves, with their boundary points:
for i in 0:3
    gmsh.model.addDiscreteEntity(1, i + 1, [i + 1, i < 3 ? i + 2 : 1])
end

# Create one discrete surface, with its bounding curves:
gmsh.model.addDiscreteEntity(2, 1, [1, 2, -3, -4])

# Add all the nodes on the surface (for simplicity... see below):
gmsh.model.mesh.addNodes(2, 1, nodes, coords)

# Add point elements on the 4 points, line elements on the 4 curves, and
# triangle elements on the surface:
for i in 0:3
    # Type 15 for point elements:
    gmsh.model.mesh.addElementsByType(i + 1, 15, [], [pnt[i + 1]])
    # Type 1 for 2-node line elements:
    gmsh.model.mesh.addElementsByType(i + 1, 1, [], lin[i + 1])
end
# Type 2 for 3-node triangle elements:
gmsh.model.mesh.addElementsByType(1, 2, [], tris)

# Reclassify the nodes on the curves and the points (since we put them all on
# the surface before with `addNodes' for simplicity)
gmsh.model.mesh.reclassifyNodes()

# Create a geometry for the discrete curves and surfaces, so that we can remesh
# them later on:
gmsh.model.mesh.createGeometry()

# Note that for more complicated meshes, e.g. for on input unstructured STL
# mesh, we could use `classifySurfaces()' to automatically create the discrete
# entities and the topology; but we would then have to extract the boundaries
# afterwards.

# Create other build-in CAD entities to form one volume below the terrain
# surface. Beware that only built-in CAD entities can be hybrid, i.e. have
# discrete entities on their boundary: OpenCASCADE does not support this
# feature.
p1 = gmsh.model.geo.addPoint(0, 0, -0.5)
p2 = gmsh.model.geo.addPoint(1, 0, -0.5)
p3 = gmsh.model.geo.addPoint(1, 1, -0.5)
p4 = gmsh.model.geo.addPoint(0, 1, -0.5)
c1 = gmsh.model.geo.addLine(p1, p2)
c2 = gmsh.model.geo.addLine(p2, p3)
c3 = gmsh.model.geo.addLine(p3, p4)
c4 = gmsh.model.geo.addLine(p4, p1)
c10 = gmsh.model.geo.addLine(p1, 1)
c11 = gmsh.model.geo.addLine(p2, 2)
c12 = gmsh.model.geo.addLine(p3, 3)
c13 = gmsh.model.geo.addLine(p4, 4)
ll1 = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
s1 = gmsh.model.geo.addPlaneSurface([ll1])
ll3 = gmsh.model.geo.addCurveLoop([c1, c11, -1, -c10])
s3 = gmsh.model.geo.addPlaneSurface([ll3])
ll4 = gmsh.model.geo.addCurveLoop([c2, c12, -2, -c11])
s4 = gmsh.model.geo.addPlaneSurface([ll4])
ll5 = gmsh.model.geo.addCurveLoop([c3, c13, 3, -c12])
s5 = gmsh.model.geo.addPlaneSurface([ll5])
ll6 = gmsh.model.geo.addCurveLoop([c4, c10, 4, -c13])
s6 = gmsh.model.geo.addPlaneSurface([ll6])
sl1 = gmsh.model.geo.addSurfaceLoop([s1, s3, s4, s5, s6, 1])
v1 = gmsh.model.geo.addVolume([sl1])
gmsh.model.geo.synchronize()

# Set this to true to build a fully hex mesh:
#transfinite = true
transfinite = false
transfiniteAuto = false

if transfinite
    NN = 30
    for c in gmsh.model.getEntities(1)
        gmsh.model.mesh.setTransfiniteCurve(c[2], NN)
    end
    for s in gmsh.model.getEntities(2)
        gmsh.model.mesh.setTransfiniteSurface(s[2])
        gmsh.model.mesh.setRecombine(s[1], s[2])
        gmsh.model.mesh.setSmoothing(s[1], s[2], 100)
    end
    gmsh.model.mesh.setTransfiniteVolume(v1)
elseif transfiniteAuto
    gmsh.option.setNumber("Mesh.MeshSizeMin", 0.5)
    gmsh.option.setNumber("Mesh.MeshSizeMax", 0.5)
    # setTransfiniteAutomatic() uses the sizing constraints to set the number
    # of points
    gmsh.model.mesh.setTransfiniteAutomatic()
else
    gmsh.option.setNumber("Mesh.MeshSizeMin", 0.05)
    gmsh.option.setNumber("Mesh.MeshSizeMax", 0.05)
end

gmsh.model.mesh.generate(3)
gmsh.write("x2.msh")

# Launch the GUI to see the results:
if !("-nopopup" in ARGS)
    gmsh.fltk.run()
end

gmsh.finalize()