Using GMSH In Simple Problem

mesh
1

Hello!

I am new to Fenics and am trying to run a very simple nonlinear problem over the Unit Sphere. Since, to my knowledge, the built in mesh for a unit sphere no longer exists in fenics, I decided to generate one with gmsh since I need to learn it anyway. I tried to look into how other people did this, but the solutions were either using an outdated version of meshio or their problem was far more complicated than mine and it was difficult to understand what was going on.

So - imagine I have a unit_sphere.msh in my working directory. How would I go about running this same code, but with the unit_sphere.msh I made in gmsh instead of the unit cube?

Thank you!

from __future__ import print_function
from fenics import *
from dolfin import *
import matplotlib.pyplot as plt
import time
import meshio

T = 2.0            # final time
num_steps = 50     # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx = ny = nz = 30
mesh = UnitCubeMesh(nx, ny, nz)
V = FunctionSpace(mesh, 'P', 1)

# Define initial value - change for cube
u_0 = Expression('exp(1 - 1 / (1 - pow(x[0], 2) + pow(x[1], 2) + pow(x[2], 2), 2))', degree=2)
u_n = interpolate(u_0, V)

# Define variational problem
u = Function(V)
u.vector()[:] = 1000
v = TestFunction(V)

D = Constant((1/1000)/.18)
a = Constant(1/.18)

beta = Constant(0.9)
one  = Constant(1)

# F = u*v*dx + dt*D*dot(grad(u), grad(v))*dx - (u_n + dt*a*(one - beta*u))*u*v*dx
F = u*v*dx + dt*D*dot(grad(u), grad(v))*dx - u_n*u*v*dx
# Create VTK file for saving solution
vtkfile = File('CarTCellModel/solution.pvd')

# Time-stepping
t = 0
for n in range(num_steps):
    # Update current time
    t += dt

    # Compute solution
    J = derivative(F, u)
    solve(F == 0, u, J = J)

    # Save to file and plot solution
    vtkfile << (u, t)

    print(assemble(u*dx))

    # Update previous solution
    u_n.assign(u)

mesh2 = UnitSquareMesh(nx, ny)
V2 = FunctionSpace(mesh2, 'P', 1)
u2 = Function(V2)
u  = interpolate(u, V)

coords = mesh2.coordinates()
dof_map = vertex_to_dof_map(V2)

# Pick your values of z
z = 1.0

for i in range(V2.dim()):
    u2.vector()[dof_map[i]] = u(coords[i][0],coords[i][1],z)

plot(u2)
# interactive()
# # Hold plot
# interactive()
# # hold plot
plt.show()
2

Have you considered the approach of: What is wrong with my mesh? - #8 by dokken
which as been covered in several posts, including: 3D Stokes flow, meshed in gmsh - #13 by dokken