ferrite-fem / ferritegmsh.jl Goto Github PK
View Code? Open in Web Editor NEWLicense: MIT License
License: MIT License
Hi,
I have problems with generate desired mesh when using spline/bspline.
When using spline/bspline one gives very messy mesh grid:
while the desired mesh from python using dolfinx and gmsh is like
Here attached the code for production of the messy mesh in Julia and desired mesh in python
using Ferrite
using FerriteGmsh
using FerriteViz
gmsh.initialize()
gdim = 2
p1 = gmsh.model.occ.addPoint(0.0, 0.0, 0, 0)
p2 = gmsh.model.occ.addPoint(1.0, 0.0, 0, 0)
p3 = gmsh.model.occ.addPoint(1.0, 0.5, 0, 0)
p4 = gmsh.model.occ.addPoint(1.0, 1.0, 0, 0)
s1 = gmsh.model.occ.addBSpline([p1, p2, p3, p4])
p2 = gmsh.model.occ.addPoint(0.0, 1.0, 0, 0)
p3 = gmsh.model.occ.addPoint(0.5, 1.0, 0, 0)
s2 = gmsh.model.occ.addSpline([p4, p3, p2, p1])
ll = gmsh.model.occ.addCurveLoop([s1, s2])
pl = gmsh.model.occ.addPlaneSurface([ll])
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], 1)
meshSize = 0.1
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",meshSize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",meshSize)
gmsh.model.mesh.generate(gdim)
dim = Int64(gmsh.model.getDimension())
facedim = dim - 1
nodes = tonodes()
elements, gmsh_elementidx = toelements(dim)
cellsets = tocellsets(dim, gmsh_elementidx)
domaincellset = cellsets["1"]
elements = elements[collect(domaincellset)]
boundarydict = toboundary(facedim)
facesets = tofacesets(boundarydict, elements)
gmsh.finalize()
grid = Grid(elements, nodes, facesets=facesets, cellsets=cellsets)
import WGLMakie
FerriteViz.wireframe(grid,markersize=10,strokewidth=2)
import gmsh
from dolfinx.io import gmshio
import pyvista
from mpi4py import MPI
from dolfinx import fem
from dolfinx import plot
gmsh.initialize()
p1 = gmsh.model.occ.addPoint(0.0, 0.0, 0, 0)
p2 = gmsh.model.occ.addPoint(1.0, 0.0, 0, 0)
p3 = gmsh.model.occ.addPoint(1.0, 0.5, 0, 0)
p4 = gmsh.model.occ.addPoint(1.0, 1.0, 0, 0)
s1 = gmsh.model.occ.addBSpline([p1, p2, p3, p4])
p2 = gmsh.model.occ.addPoint(0.0, 1.0, 0, 0)
p3 = gmsh.model.occ.addPoint(0.5, 1.0, 0, 0)
s2 = gmsh.model.occ.addSpline([p4, p3, p2, p1])
ll = gmsh.model.occ.addCurveLoop([s1, s2])
pl = gmsh.model.occ.addPlaneSurface([ll])
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], 1)
meshSize = 0.1
gdim = 2
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",meshSize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",meshSize)
gmsh.model.mesh.generate(gdim)
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
V = fem.functionspace(domain, ("Lagrange", 1))
cells, types, x = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.show()
gmsh.finalize()
Until no one else comes up with a solid approach gmsh.jl should live inside this repo, since I'll never be able to register gmsh.jl
Hello,
I'm new in FEM and in Ferrite. So far I have been able to use this library and Gmsh with different geometry and it worked fine, but this one gives me trouble:
Here is the error message:
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
Here is a minimal geo file that reproduce the problem:
// 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};
And here is a minimal Julia code to reproduce the error:
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
I have tried with different mesh type and size. I also tried feeding the geo file to Ferrite instead of a msh file, but the problem persist. Quad gives a higher detJ, and reducing the mesh size also gives higher detJ, but detJ stays negative. I can't really decrease the mesh size further.
Is there anything else I can try?
Currently, node and edgesets are not supported. Should be not too hard to support them.
Currently parsing of embedded elements (e.g. Triangles in 3D) does not fully work. As a hotfix I have
try
togrid("/home/dogiermann/untitled.msh")
catch
end
nodes = tonodes();
elements, gmsh_elementidx = toelements(2);
grid = Grid(elements, nodes);
but I think it would be nice to have a proper solution.
For me this is important to allow parsing of mixed-dimensional problems later on.
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
For the upcoming Ferrite v0.4 release, an update of compat will be needed: Ferrite-FEM/Ferrite.jl#668
The problem here seems to be the call to gmsh.model.mesh.set_periodic()
, it doesn't happen without that.
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))
I played around a bit with this package, mainly just for reading .geo and .msh files. It is very nice compared to my own little mesh reader I used in the past!
Anyway, I was thinking I could use "Physical Groups" to mark regions in the mesh and have that translate to cellsets in Ferrite. That seems to work if there is just one surface in the group, the remaining ones seems to be forgotten. I attach the .msh file here: test.log (I had to change .msh
extension .log
for GitHub to accept it). Checking e.g.
julia> using Ferrite, FerriteGmsh
julia> g = saved_file_to_grid("test.msh");
julia> vtk_grid("test", g) do vtk
vtk_cellset(vtk, g)
end
only one of the circles are included in the "inclusions"
cellset.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.