Calculation of the steady-state heat field in the part: Julia, Gmsh, Gridap, VTK and Makie¶

Let us apply a chain of tools that allows us to design in Engee, parameterise the model, apply a computational mesh, construct the problem to be solved by finite element method and finally calculate the heat transfer in a complex part.

Introduction¶

For problems requiring the exact solution of partial equations (e.g. heat transfer or deformable body mechanics), specialised tools (like Gridap or FEniCS) give more accurate and flexible results. They can handle arbitrary meshes and complex boundary conditions, which is critical for scientific calculations.

Pkg.add(["Gridap", "GridapGmsh", "CairoMakie", "WGLMakie", "WriteVTK", "ReadVTK", "GeometryBasics"])

For this example, we use a library whose installation automatically adds the Gmsh package to the system.

using GridapGmsh

Let's take a simple model example from the repository of this library and modify it slightly:

the boundary 1 region will combine the inner surface of the wider hole (right)
area boundary 2 - the entire inner surface of the hole with a protrusion (left)

The creation of such a parameterised model can be done in a text editor, e.g. by visualising the results in https://gmsh.info/ or in one of the many CAD systems for solid modelling, including open source CAD systems: OpenCASCADE, OpenSCAD and FreeCAD.

Let's load the example in .geo format (basic format for geometry in Gmsh) and create a mesh:

GridapGmsh.gmsh.initialize()
GridapGmsh.gmsh.clear()
GridapGmsh.gmsh.parser.parse("$(@__DIR__)/demo.geo");
GridapGmsh.gmsh.model.geo.synchronize();
GridapGmsh.gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.01)  # минимальный размер
GridapGmsh.gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)  # максимальный размер
GridapGmsh.gmsh.model.mesh.generate(3);
GridapGmsh.gmsh.write("demo.msh");

Let's output the model of our pattern, which is a set of tetrahedrons (Makie doesn't allow working with tetrahedrons, so we manually converted them to triangles). The colouring of the graph corresponds to the X coordinate of the points.

using Gridap, GridapGmsh, GeometryBasics

import CairoMakie: mesh, NoShading
# import WGLMakie: mesh, NoShading # На замену CairoMakie - для интерактивных графиков

model = GmshDiscreteModel("demo.msh");

grid = get_grid(get_triangulation(model))
points = reinterpret(reshape, Float64, grid.node_coordinates)
cells = hcat(grid.cell_node_ids...)

tetrahedra = [ (c[1], c[2], c[3], c[4]) for c in eachcol(cells) ]
all_faces = [ [TriangleFace{Int}(ta, tb, tc),TriangleFace{Int}(ta, tb, td),TriangleFace{Int}(ta, tc, td),TriangleFace{Int}(tb, tc, td)] for (ta, tb, tc, td) in tetrahedra ]
all_faces = reshape(reinterpret(Int, vcat(all_faces...)), 3, :)
colors = points[1,:]

mesh(points', all_faces', color=colors, shading=NoShading)

Info    : Reading 'demo.msh'...
Info    : 188 entities
Info    : 3302 nodes
Info    : 12065 elements
Info    : Done reading 'demo.msh'

Warning : Gmsh has aleady been initialized

Numerical problem formulation¶

We use the code from basic example for GridapGmsh, where we solve the Laplace's equation with Dirichlet boundary conditions in the weak (variational) formulation by the finite element method (FEM).

Laplace's equation and boundary conditions: We search for the function $u(x)$, satisfying the equation $\Delta u = 0$ in the region $\Omega$ with boundary conditions $u = 0$ at $\Gamma_1$, $u = 1$ at $\Gamma_2$, where $\Omega$ is the region defined by the mesh demo.msh, $\Gamma_1$ and $\Gamma_2$ are the boundaries with tags "boundary1" and "boundary2" respectively.

Slack form (variational problem formulation) is a key method for modern FEM, it relaxes the smoothness requirements of the solution, also with it the boundary conditions can be given directly without approximating them. And, of course, it is the main method of discretisation of the problem, which allows us to specify it in basis functions and solve it as a linear system $Au = b$. So, our problem in the variational formulation is transformed into $a(u, v) = \int_{\Omega} \nabla v \cdot \nabla u \, d\Omega = 0 \quad \forall v \in V_0$, where $\nabla$ is the gradient, $\Omega$ is the region of integration, $V_0$ is the trial space turning to zero at the boundaries.

Note how close the description given in the code ∫( ∇(v)⋅∇(u) )dΩ is to the Laplace equation in the variational formulation $a(u, v) = \int_{\Omega} \nabla v \cdot \nabla u \, d\Omega = 0 \quad \forall v \in V_0$.

using Gridap

model = GmshDiscreteModel("demo.msh");         # Загружаем геометрию сетки из файла *.msh
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)  # Линейные лагранжевы элементы 1-го порядка
V = TestFESpace(model,reffe,dirichlet_tags=["boundary1","boundary2"]) # Пробное пространство
U = TrialFESpace(V,[0,1])                      # Определение пространства решений
                                               # - на boundary1 значение решения u = 0
                                               # - на boundary2 значение u = 1
Ω = Triangulation(model)                       # Триангуляция и мера интегрирования
dΩ = Measure(Ω,2*order)
a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ                      # Слабая форма уравнения Лапласа
l(v) = 0                                       # Линейная форма (правая часть)
op = AffineFEOperator(a,l,U,V)                 # Собираем линейный оператор задачи
uh = solve(op)                                 # Решение системы

writevtk(Ω,"demo",cellfields=["uh"=>uh]);

Info    : Reading 'demo.msh'...
Info    : 188 entities
Info    : 3302 nodes
Info    : 12065 elements
Info    : Done reading 'demo.msh'

We will go a little further and, in addition to on-the-fly visualisation, we will show how to visualise the calculation saved in VTK. But, first, we will show the solution we have just obtained:

using GeometryBasics
nodes = GridapGmsh.get_node_coordinates(Ω)
cell_node_ids = GridapGmsh.get_cell_node_ids(Ω)
uh_values = [uh(node) for node in nodes]
points = [Point3f(x...) for x in nodes]
faces = [TriangleFace(cell_node_ids[i][1], cell_node_ids[i][2], cell_node_ids[i][3]) for i in 1:length(cell_node_ids)]
mesh(points, faces, color=uh_values, colormap=:inferno, shading=NoShading)

Visualisation of the solution from VTK file¶

For viewing VTK files, the most convenient solution is the open source programme ParaView. But it is also valuable for us to have interactive reports that are fully discoverable in Engee. So we will open this file in Julia and output the solution directly in this script.

If we run this part asynchronously with the solution (which can take a long time), we will need to connect these libraries:

using WriteVTK, ReadVTK

They will help us read the VTK file and display it on 3D graphics. There are additional libraries for some of the functions, but we will make do with low-level commands for now:

create a matrix points for vertex coordinates
get the vector uh_values, where the values of the field uh of our solution will lie.
split tetrahedrons cells into a set of triangles all_faces

vtk = VTKFile("demo.vtu")
points = get_points(vtk)
uh_values = get_data(get_point_data(vtk)["uh"])

cells = get_cells(vtk)

# Разобьем массив cells.connectivity на тетраедры, а затем на треугольники
using GeometryBasics
all_faces = TriangleFace{Int}[]
offset = 1
while offset <= length(cells.connectivity)
    ta, tb, tc, td = cells.connectivity[offset:offset+3]
    push!(all_faces, TriangleFace(ta, tb, tc))  # грань 1
    push!(all_faces, TriangleFace(ta, tb, td))  # грань 2
    push!(all_faces, TriangleFace(ta, tc, td))  # грань 3
    push!(all_faces, TriangleFace(tb, tc, td))  # грань 4
    offset += 4
end

# CairoMakie не принимает последний треугольник в модели, временно используем [1:end-1]
mesh(points, all_faces[1:end-1], color=uh_values, colormap=:coolwarm, shading=NoShading)

Conclusion¶

We have implemented the solution of the Laplace equation on a tetrahedral mesh, saved the results in VTK format and visualised them - this is the basis for modelling physical fields (thermal, electrical, mechanical). This approach ensures reproducibility, accuracy and scalability to real engineering problems where detailed calculations, rather than simplified system approximations, are required.