ferrite-fem / ferritegmsh.jl Goto Github PK

using Gmsh, FerriteGmsh, Ferrite

# Initialize gmsh
Gmsh.initialize()

# Add the points
h = 0.05
o = gmsh.model.geo.add_point(0.0, 0.0, 0.0, h)
p1 = gmsh.model.geo.add_point(0.5, 0.0, 0.0, h)
p2 = gmsh.model.geo.add_point(1.0, 0.0, 0.0, h)
p3 = gmsh.model.geo.add_point(0.0, 1.0, 0.0, h)
p4 = gmsh.model.geo.add_point(0.0, 0.5, 0.0, h)

# Add the lines
l1 = gmsh.model.geo.add_line(p1, p2)
l2 = gmsh.model.geo.add_circle_arc(p2, o, p3)
l3 = gmsh.model.geo.add_line(p3, p4)
l4 = gmsh.model.geo.add_circle_arc(p4, o, p1)

# Create the closed curve loop and the surface
loop = gmsh.model.geo.add_curve_loop([l1, l2, l3, l4])
surf = gmsh.model.geo.add_plane_surface([loop])

# Synchronize the model
gmsh.model.geo.synchronize()

# Create the physical domains
gmsh.model.add_physical_group(1, [l1], -1, "right")
gmsh.model.add_physical_group(1, [l2], -1, "outer")
gmsh.model.add_physical_group(1, [l3], -1, "left")
gmsh.model.add_physical_group(1, [l4], -1, "inner")
gmsh.model.add_physical_group(2, [surf])

# # Add the periodicity constraint using 4x4 affine transformation matrix,
# # see https://en.wikipedia.org/wiki/Transformation_matrix#Affine_transformations
transformation_matrix = zeros(4, 4)
transformation_matrix[1, 2] = 1  # -sin(-pi/2)
transformation_matrix[2, 1] = -1 #  cos(-pi/2)
transformation_matrix[3, 3] = 1
transformation_matrix[4, 4] = 1
transformation_matrix = vec(transformation_matrix')
gmsh.model.mesh.set_periodic(1, [l1], [l3], transformation_matrix)

# Generate a 2D mesh
gmsh.model.mesh.generate(2)

# Convert the mesh to Ferrite Grid
grid_bad = FerriteGmsh.togrid()
vtk_save(vtk_grid("grid_bad", grid_bad))
grid_good = mktempdir() do dir
    path = joinpath(dir, "mesh.msh")
    gmsh.write(path)
    togrid(path)
end
vtk_save(vtk_grid("grid_good", grid_good))

Info    : Reading 'folder\file.msh'...
Info    : 18 entities
Info    : 507 nodes
Info    : 1021 elements
Info    : Done reading 'folder\file.msh'
ERROR: LoadError: ArgumentError: det(J) is not positive: det(J) = -1.1052236851259822
Stacktrace:
 [1] throw_detJ_not_pos(detJ::Float64)
   @ Ferrite C:\some_folder\.julia\packages\Ferrite\qYVTc\src\FEValues\common_values.jl:17
 [2] reinit!(cv::CellScalarValues{2, Float64, RefTetrahedron}, x::Vector{Vec{2, Float64}})
   @ Ferrite C:\some_folder\.julia\packages\Ferrite\qYVTc\src\FEValues\cell_values.jl:163
 [3] reinit!(cv::CellScalarValues{2, Float64, RefTetrahedron}, ci::CellIterator{2, Triangle, Float64, DofHandler{2, Triangle, Float64}})
   @ Ferrite C:\some_folder\.julia\packages\Ferrite\qYVTc\src\iterators.jl:130
 [4] top-level scope
   @ d:\some\folder\temp.jl:17
in expression starting at d:\some\folder\temp.jl:16

// Gmsh project created on Wed Mar 02 09:48:56 2022
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0, 10, 0, 1.0};
//+
Point(3) = {10, 10, 0, 1.0};
//+
Point(4) = {10, 20, 0, 1.0};
//+
Point(5) = {20, 20, 0, 1.0};
//+
Point(6) = {20, 5, 0, 1.0};
//+
Point(7) = {30, 5, 0, 1.0};
//+
Point(8) = {30, 0, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 5};
//+
Line(5) = {5, 6};
//+
Line(6) = {7, 7};
//+
Line(7) = {6, 7};
//+
Line(8) = {7, 8};
//+
Line(9) = {8, 1};
//+
Curve Loop(1) = {2, 3, 4, 5, 7, 8, 9, 1};
//+
Plane Surface(1) = {1};

using Ferrite
using FerriteGmsh

folder = "some folder here"
grid = saved_file_to_grid(folder*"file.msh")

dim = 2
ip = Lagrange{dim, RefTetrahedron, 1}() # or RefCube depending on mesh type
qr = QuadratureRule{dim, RefTetrahedron}(2) # or RefCube depending on mesh type
cellvalues = CellScalarValues(qr, ip);

dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh);

@inbounds for cell in CellIterator(dh)
    reinit!(cellvalues, cell)
end

julia> using Ferrite, FerriteGmsh

julia> g = saved_file_to_grid("test.msh");

julia> vtk_grid("test", g) do vtk
           vtk_cellset(vtk, g)
       end