from dolfin import *
from mpi4py import MPI
import meshio
export_path = "/home/your_path/"
mesh_comm = MPI.COMM_WORLD
from solver_gym_fenics import create_mesh
mesh = Mesh()
# The mesh generated with gmsh's python API
mesh_file = "cylinder_mesh.msh"
# xdmf file that will contain the mesh
xdmf_file = "mesh.xdmf"
# xdmf file that will contain the 1D elements mesh (facets)
facet_file = "facet_mesh.xdmf"
# read mesh from msh file
msh = meshio.read(mesh_file)
# create the 1D element mesh
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write(facet_file, line_mesh)
# create the xdmf mesh
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write(xdmf_file, triangle_mesh)
# writing xdmf file into Mesh object
with XDMFFile(xdmf_file) as infile:
infile.read(mesh)
# retrieving the facets
mvc1 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(facet_file) as infile:
infile.read(mvc1, "name_to_read")
facet = cpp.mesh.MeshFunctionSizet(mesh, mvc1)
velocity_element = VectorElement("Lagrange", mesh.ufl_cell(), 2)
pressure_element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, velocity_element * pressure_element)
w = Function(W)
wall_marker, inlet_marker, outlet_marker = range(1, 4)
obstacle_marker = 4
solution_file = XDMFFile(mesh_comm, export_path + "mwe.xdmf")
restart = False
if not restart:
i = 0
bc = DirichletBC(W.sub(0), Constant((1, -1)), facet, wall_marker)
solution_file.write(mesh)
else:
i = 1
bc = DirichletBC(W.sub(0), Constant((-1, 1)), facet, wall_marker)
bc.apply(w.vector())
solution_file.write_checkpoint(w.sub(0), "mixed_function", i, XDMFFile.Encoding.HDF5, True)
# solution_file.write(w, i)
solution_file.close()
Here you can find my mwe. I launch it, it creates the mesh file and a solution file. But there is an error in paraview, and I can’t know why. I hope you understood what I wanted to do in my post, so that you can explain me what I am doing wrong with the XDMFFile
class.
This mwe needs a cylinder_mesh
file, which is created using this script:
############################################
#------------ Mesh Creation ---------------#
############################################
import gmsh
import time
import sys
import numpy
from env_data import r, c_x, c_y, gdim
from solver_gym_fenics import printr0
# Calculating the time that mesh generation will take
tic = time.perf_counter()
gmsh.initialize()
c_x, c_y = c_x[0], c_y[0]
gdim = 2
gmsh.model.add("cylinder")
# target mesh size for any point
lc = 0.05
# lc = 0.1
############################################
#---------- Geometry creation -------------#
############################################
## 0 dim : Points definition
# Rectangle delimiting the flow
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(2.2, 0, 0, lc, 2)
gmsh.model.geo.addPoint(2.2, 0.41, 0, lc, 3)
gmsh.model.geo.addPoint(0, 0.41, 0, lc, 4)
# Points to create the circle for the cylinder (4 points + center)
gmsh.model.geo.addPoint(c_x, c_y, 0, lc, 5) # center
gmsh.model.geo.addPoint(c_x + r, c_y, 0, lc, 6)
gmsh.model.geo.addPoint(c_x, c_y + r, 0, lc, 7)
gmsh.model.geo.addPoint(c_x - r, c_y, 0, lc, 8)
gmsh.model.geo.addPoint(c_x, c_y - r, 0, lc, 9)
# Points to create the polygon (pentagon in this case) that will include the cylinder and the wake
# In that pentagon, the mesh will be refined
gmsh.model.geo.addPoint(c_x - 2*r, c_y, 0, lc, 10)
gmsh.model.geo.addPoint(c_x, c_y + 2*r, 0, lc, 11)
gmsh.model.geo.addPoint(2.2, c_y + 2*r, 0, lc, 12)
gmsh.model.geo.addPoint(2.2, c_y - 2*r, 0, lc, 13)
gmsh.model.geo.addPoint(c_x, c_y - 2*r, 0, lc, 14)
## 1 dim: Lines and Curves
# Lines for the rectangle
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 13, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addLine(12 , 3, 5)
# Lines for the pentagone
gmsh.model.geo.addLine(10, 11, 6)
gmsh.model.geo.addLine(11, 12, 7)
gmsh.model.geo.addLine(12, 13, 8)
gmsh.model.geo.addLine(13, 14, 9)
gmsh.model.geo.addLine(14, 10, 10)
# Creation of the circle for the cylinder (4 CircleArcs joining the 4 points previously defined)
gmsh.model.geo.addCircleArc(6, 5, 7, tag=11)
gmsh.model.geo.addCircleArc(7, 5, 8, tag=12)
gmsh.model.geo.addCircleArc(8, 5, 9, tag=13)
gmsh.model.geo.addCircleArc(9, 5, 6, tag=14)
## 2 dim
# Creation of the plane surface: rectangle - pentagone
gmsh.model.geo.addCurveLoop([1, 2, 9, 10, 6, 7, 5, 3, 4], 1)
gmsh.model.geo.addCurveLoop(list(range(6, 11)), 2)
gmsh.model.geo.addPlaneSurface([1], 1)
############################################
#------ Parameters & Physical groups ------#
############################################
# Creation of a curve for the polygon
curved_loop_list = [2]
# Creation of curves loops for the cylinder
gmsh.model.geo.addCurveLoop([11, 12, 13, 14], 3)
curved_loop_list.append(3)
# Creation of the plane surface: polygon - cylinder
gmsh.model.geo.addPlaneSurface(curved_loop_list, 2)
gmsh.model.geo.synchronize()
## Fields
# Distance field to calculate distance from the three cylinders
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "CurvesList", list(range(11, 15)))
gmsh.model.mesh.field.setNumber(1, "Sampling", 100)
# Threshold field to refine the mesh around the three cylinders
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "SizeMin", lc / 9)
gmsh.model.mesh.field.setNumber(2, "SizeMax", lc)
gmsh.model.mesh.field.setNumber(2, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(2, "DistMax", r)
# Constant mesh size inside polygon
gmsh.model.mesh.field.add("Constant", 3)
gmsh.model.mesh.field.setNumber(3, "VIn", lc/6)
gmsh.model.mesh.field.setNumber(3, "VOut", lc)
gmsh.model.mesh.field.setNumbers(3, "SurfacesList", [2])
# Distance field to calculate distance from the outlet
gmsh.model.mesh.field.add("Distance", 4)
gmsh.model.mesh.field.setNumbers(4, "CurvesList", [2, 8, 5])
gmsh.model.mesh.field.setNumber(4, "Sampling", 100)
# Threshold field to refine the mesh around the three cylinders
gmsh.model.mesh.field.add("Threshold", 5)
gmsh.model.mesh.field.setNumber(5, "InField", 4)
gmsh.model.mesh.field.setNumber(5, "SizeMin", lc/2)
gmsh.model.mesh.field.setNumber(5, "SizeMax", lc)
gmsh.model.mesh.field.setNumber(5, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(5, "DistMax", r)
# Final field as in t10.py, take the minimum mesh size for all Threshold fields.
gmsh.model.mesh.field.add("Min", 6)
gmsh.model.mesh.field.setNumbers(6, "FieldsList", [2, 3, 5])
gmsh.model.mesh.field.setAsBackgroundMesh(6)
## Create Physical Groups
gmsh.model.addPhysicalGroup(1, [1, 3], name="walls") # walls: up and down sides of the rectangle
gmsh.model.addPhysicalGroup(1, [4], name="inlet") # inlet: left side of the rectangle
gmsh.model.addPhysicalGroup(1, [2, 8, 5], name="outlet") # outlet: right side of the rectangle
gmsh.model.addPhysicalGroup(1, list(range(11, 15)), name="cylinder") # obstacel
gmsh.model.addPhysicalGroup(2, [1, 2], name="control_volume", tag=5) # Reunion of the two plane surfaces: where fluid will flow
gmsh.model.geo.synchronize()
# printing how long the mesh creation took
toc = time.perf_counter()
print("Mesh generation took {} secondes !".format(toc-tic))
############################################
#----- Mesh Saving & Post-processing ------#
############################################
# Generate and write to file
gmsh.model.mesh.generate(gdim)
gmsh.write("cylinder_mesh.msh")
gmsh.write("cylinder_mesh.vtk")
# Launch the GUI to see the results:
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
# retrieving points of the fluid domain (our physical group in gmsh)
points = gmsh.model.mesh.getNodesForPhysicalGroup(2, 5)[1]
# getting the coordinates of the mesh points
coordinates = []
for i in range(0, len(points), 3):
coordinates.append(list(points[i:i+3]))
two_D_coord = []
for coord in coordinates:
two_D_coord.append(coord[:-1])
# saving the coordinates a file
with open('mesh_data.npy', 'wb') as f:
numpy.save(f, two_D_coord)
So what I would like to do:
- create the mesh using the second script.
- launch the first script with
restart = False
- launch the first script again with
restart = True
.
And I would like to have a solution_file.xdmf
containing 2 steps of the solution, one created when restart = False
, and the other created when restart = True
.