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?
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, 2) + pow(x, 2) + pow(x, 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],coords[i],z) plot(u2) # interactive() # # Hold plot # interactive() # # hold plot plt.show()