How to convert the output mesh to a .geo file?

avishmj · July 20, 2021, 1:09pm

In the below MWE, I am considering an annulus and applying traction on the inner boundary. Now I want to introduce a mesh on the inner hole for further computations. According to me, this can be done in two ways:

Create the mesh in-script by collecting the nodes on the inner boundary of the deformed geometry

or

Create a .geo file of the deformed mesh and introduce the mesh in gmsh.

Here is the deformed geometry I am working with followed by an MWE:

MWE:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from ufl import *
from ufl import cofac
from dolfin import *
import numpy as np
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C = C2 - C1

mesh = generate_mesh(C, 100)

class inner(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class outer(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary
class annulus(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4


boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
inner().mark(boundary_markers, 1)
outer().mark(boundary_markers, 2)
annulus().mark(surface_markers, 3)



plot(mesh)
plt.show()

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(Constant(1)*ds(1)),assemble(Constant(1)*ds(2)))
print(assemble(Constant(1)*dx(3)))


V = VectorFunctionSpace(mesh, "Lagrange", 1)
n=FacetNormal(mesh)

du = TrialFunction(V)
v  = TestFunction(V)             
u  = Function(V)                 
T = Constant(5.0)
bcs=[]

# Kinematics
d = mesh.geometry().dim()
I = Identity(d)             
F = I + grad(u)             
C = F.T*F                   
M = cofac(F)
Finvt = (inv(F)).T
Ic = tr(C)
J  = det(F)

E, nu = 10.0, 0.499   
mu, lmbda = 27.9, Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) + ((nu)*(mu)/(1-2*nu))*(J-1)**2 - mu*(ln(J))

Coords = Expression(('x[0]','x[1]'),degree=0)
coords = project(Coords, V)
phi = coords+u

Pi = psi*dx - (0.5)*dot(M*T*(-n), phi)*ds(1)
ddPI = derivative(Pi, u, v)
solve(ddPI == 0, u, bcs,
     form_compiler_parameters=ffc_options)
     
     
mesh1 = Mesh(mesh) 
X = mesh1.coordinates()
X += np.vstack(map(u, X))

mesh_moved = Mesh(mesh1)
outer_nodes = BoundaryMesh(mesh_moved, "exterior", True)
translated_center = np.mean(outer_nodes.coordinates(), axis=0)

n = np.shape(mesh_moved.coordinates())[0]
translation = np.zeros([n, 2])
translation[:, 0] = translated_center[0]; translation[:, 1] = translated_center[1]
mesh_moved.coordinates()[:] = mesh_moved.coordinates() - translation
trans = np.zeros([n,1])
trans[:,0] = sqrt(translated_center[0]**2 + translated_center[1]**2)



file = File("MWE.pvd");
file << u;

I did search some literature on how to do this but couldn’t find a suitable one. Also, I haven’t tried this before. So any suggestion would be much appreciated.

dokken · July 20, 2021, 4:06pm

Please note that a .geo file is a file containing GMSH language specific commands, and not usually vertices and cells (as their msh format does).

I would suggest the following:

Write your current mesh to file (for instance xdmf).
Use meshio to convert it to the msh format
Use the Merge command in gmsh to add the mesh to geo.

Minimal example:

from dolfin import UnitSquareMesh, XDMFFile

mesh = UnitSquareMesh(2, 2)
xdmf = XDMFFile("mesh.xdmf")
xdmf.write(mesh)
xdmf.close()

import meshio
meshio_mesh = meshio.read("mesh.xdmf")
meshio.write("meshio_mesh.msh", meshio_mesh)

Geo file (for instance called mesh.geo):

Merge "meshio_mesh.msh";

avishmj · July 21, 2021, 1:04pm

Hi @dokken , Thanks for the suggestion.

Do I also need to pass the boundary info to meshio.write? Currently, it seems to be only giving me the mesh. So when I am trying to fill in the hole, gmsh is not letting me create a plane surface on the hole because the boundaries are not registered. Or maybe I am doing something wrong, which is most probable.

Here is the MWE again, with your suggestion plugged in:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from ufl import *
from ufl import cofac
from dolfin import *
import numpy as np
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C = C2 - C1

mesh = generate_mesh(C, 100)

class inner(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class outer(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary
class annulus(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4


boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
inner().mark(boundary_markers, 1)
outer().mark(boundary_markers, 2)
annulus().mark(surface_markers, 3)



#plot(mesh)
#plt.show()

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(Constant(1)*ds(1)),assemble(Constant(1)*ds(2)))
print(assemble(Constant(1)*dx(3)))


V = VectorFunctionSpace(mesh, "Lagrange", 1)
n=FacetNormal(mesh)

du = TrialFunction(V)
v  = TestFunction(V)             
u  = Function(V)                 
T = Constant(5.0)
bcs=[]

# Kinematics
d = mesh.geometry().dim()
I = Identity(d)             
F = I + grad(u)             
C = F.T*F                   
M = cofac(F)
Finvt = (inv(F)).T
Ic = tr(C)
J  = det(F)

E, nu = 10.0, 0.499   
mu, lmbda = 27.9, Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) + ((nu)*(mu)/(1-2*nu))*(J-1)**2 - mu*(ln(J))

Coords = Expression(('x[0]','x[1]'),degree=0)
coords = project(Coords, V)
phi = coords+u

Pi = psi*dx - (0.5)*dot(M*T*(-n), phi)*ds(1)
ddPI = derivative(Pi, u, v)
solve(ddPI == 0, u, bcs,
     form_compiler_parameters=ffc_options)
     
     
mesh1 = Mesh(mesh) 
X = mesh1.coordinates()
X += np.vstack(map(u, X))

mesh_moved = Mesh(mesh1)
outer_nodes = BoundaryMesh(mesh_moved, "exterior", True)
translated_center = np.mean(outer_nodes.coordinates(), axis=0)

n = np.shape(mesh_moved.coordinates())[0]
translation = np.zeros([n, 2])
translation[:, 0] = translated_center[0]; translation[:, 1] = translated_center[1]
mesh_moved.coordinates()[:] = mesh_moved.coordinates() - translation
trans = np.zeros([n,1])
trans[:,0] = sqrt(translated_center[0]**2 + translated_center[1]**2)

xdmf = XDMFFile("meshmvd.xdmf")
xdmf.write(mesh_moved)
xdmf.close()

import meshio
meshio_mesh = meshio.read("meshmvd.xdmf")
meshio.write("annulus.msh", meshio_mesh)

file = File("MWE.pvd");
file << u;

Then I followed the following steps.

Downloaded the annulus.msh file
Created an empty .geo file, named it annulusesp.geo
File → merge → annulus.msh

This gives me the following geometry in gmsh:

Now since no boundaries are passed on to the annulus.msh file (I assume), it is not letting me create a plane surface on the hole by selecting the inner boundary of the mesh.

Now to fill in the hole, do I need to create a spline by selecting the nodes on the edge (my original geometry is messy than the one in the MWE) or is there a way I can pass the boundaries to meshio.write("annulus.msh", meshio_mesh)?

Topic		Replies	Views
Exact Procedure to work with a .geo file on FEniCS (with meshio)? mesh	11	3597	July 18, 2019
Need some help in converting gmsh file into the form which can be used in fenics code mesh	5	398	April 3, 2023
Error in using meshio.write() and create_mesh() for a mesh imported from gmsh	7	1843	September 22, 2021
Remeshing (using gmsh?) mesh	5	3001	July 25, 2019
Working with different subdomains: from meshio to dolfin mesh	10	2329	March 11, 2021

How to convert the output mesh to a .geo file?

Related topics