CAN ANY ONE SUGGEST A SIMPLE WAY TO CREATE A SPHERICAL MESH OF CERTAIN RADIUS?
‘’’
python3
import random
from dolfin import *
Create mesh and define function space with periodic boundary conditions
sphere = Sphere(Point(0, 0, 0), 200)
mesh = generate_mesh(sphere, 50)
P1 = FiniteElement(mesh,“CG”, 1)
ME = FunctionSpace(mesh, P1*P1,)
Class representing the intial conditions
class InitialConditions(UserExpression):
def init(self, **kwargs):
random.seed(2 + MPI.rank(MPI.comm_world))
super().init(*kwargs)
def eval(self, values, x):
values[0] = 0.50 + 0.02(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
def init(self, a, L):
NonlinearProblem.init(self)
self.L = L
self.a = a
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
Model parameters
lmbda = 2 # surface parameter
dt = 0.5 # time step
theta = 0.5 # time stepping family, e.g. theta=1 → backward Euler, theta=0.5 → Crank-Nicolson
Form compiler options
parameters[“form_compiler”][“optimize”] = True
parameters[“form_compiler”][“cpp_optimize”] = True
Create mesh and define function space
Define trial and test functions
du = TrialFunction(ME)
q, v = TestFunctions(ME)
Define functions
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
Split mixed functions
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
Compute the chemical potential df/dc
c = variable(c)
f = 5*(c-0.3)*2(0.7-c)**2
dfdc = diff(f, c)
mu_(n+theta)
mu_mid = (1.0-theta)mu0 + thetamu
Weak statement of the equations
L0 = cqdx - c0qdx + dtdot(5grad(mu_mid), grad(q))dx
L1 = muvdx - dfdcvdx - dot(lmbdagrad(c), grad(v))*dx
L = L0 + L1
Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters[“linear_solver”] = “lu”
solver.parameters[“convergence_criterion”] = “incremental”
solver.parameters[“relative_tolerance”] = 1e-6
Output file
file = File(“sol1/output.pvd”, “compressed”)
Step in time
t = 0.0
T = 200*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
‘’’