Hi all.
I’ve been working with FEniCSx for some days and after a lot of help from this forum I’ve been able to write a code that successfully derives and solves my system of PDEs. Now I would like to use the result for further calculations, but I have several doubts:
-
The problem is such that spherical coordinates and basis are the natural set to work on. However, for easy visualisation and further work I would need to transform the result (a vector field) and the mesh to Cartesian coordinates and basis. How can I achieve this?
-
Can I somehow take the result and the mesh (or their Cartesian versions) and plot them with matplotlib? I’ve been using pyvista, but it doesn’t seem so versatile (truth be told, I also I don’t really understand all the steps of the script).
-
The system of PDEs comes from minimising the energy functional
F[\mathbf{m}]=\int_{\Omega}\mathrm dx\,r^2\sin\theta f_B(\mathbf m)+\int_{\partial\Omega}\mathrm ds\,r^2\sin\theta f_S(\mathbf m) . Solving it yields the vector field \mathbf m^*. I would now like to evaluate F(\mathbf m^*). Any hints on how I can do this?
Here’s the code in case it helps:
import numpy as np
import matplotlib.pyplot as plt
import ufl
from dolfinx import mesh, fem, nls, log, plot
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
dim = 3 # dimension of the OP
Rbase = 1 # radius of the base (sphere)
Rcap = Rbase # radius of the cap
minq = 0
bins = 10
domain = mesh.create_rectangle(MPI.COMM_WORLD, np.array([[0,minq],[Rbase,np.pi]]),[bins,int(bins*(np.pi-minq))])
x = ufl.SpatialCoordinate(domain)
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, dim)
V = fem.FunctionSpace(domain, element) # space of vector functions with dim components
boundaries = [(0, lambda x: np.isclose(x[0], 0)), # origin
(1, lambda x: np.isclose(x[0], Rbase) & np.greater_equal(x[1], np.pi/2)), # base
(2, lambda x: np.isclose(x[0], Rcap) & np.less(x[1], np.pi/2)), # cap
(3, lambda x: np.isclose(x[1], minq)),
(4, lambda x: np.isclose(x[1], np.pi))]
facet_indices, facet_markers = [], []
tdim = domain.topology.dim; fdim = tdim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.measure.Measure("ds", domain=domain, subdomain_data=facet_tag)
a = 0.1; c = 0.01 # thermotropic parameters
d = 0.01 # stiffness
w = 0.1 # anchoring
m0base = [1,0,0] # preferred value of the OP in the base
m0cap = [0,1,0] # preferred value of the OP in the cap
m = fem.Function(V, name='m') # dim-dimensional vector OP
dm = ufl.TrialFunction(V) # variation of the OP
phi = ufl.TestFunction(V) # test function
def Jacob(x):
"""Jacobian of the transformation from Cartesian to spherical coordinates."""
xr, xq = x
return xr**2*ufl.sin(xq)
def curl_sph(A, x):
"""Curl of a vector field depending only on r and theta in spherical coordinates."""
xr, xq = x; jr, jq = 0, 1
Ar, Aq, Af = A
curl_r = (ufl.Dx(Af*ufl.sin(xq),jq))/(xr*ufl.sin(xq))
curl_q = (-ufl.Dx(xr*Af,jr))/xr
curl_f = (ufl.Dx(xr*Aq,jr)-ufl.Dx(Ar,jq))/xr
return ufl.as_vector([curl_r, curl_q, curl_f])
def div_sph(A, x):
"""Divergence of a vector field depending only on r and theta in spherical coordinates."""
xr, xq = x; jr, jq = 0, 1
Ar, Aq, Af = A
div_r = ufl.Dx(xr**2*Ar,jr)/xr**2
div_q = ufl.Dx(Aq*ufl.sin(xq),jq)/(xr*ufl.sin(xq))
div_f = 0
return div_r + div_q + div_f
def fB(u, x):
"""Bulk energy density considering the thermotropic and elastic contributions."""
fT = a/2.*ufl.dot(u,u) + c/4.*ufl.dot(u,u)**2
fD = d*( (div_sph(u,x))**2 + (ufl.dot(u,curl_sph(u,x)))**2 + (ufl.cross(u,curl_sph(u,x)))**2 )
return fT+fD
def fS(u,u0):
"""Surface energy density. Quadratic anchoring to a preferred value of the OP in the boundary."""
u0 = ufl.as_vector(u0)
return w/2*ufl.dot(u-u0,u-u0)
F = fB(m,x)*Jacob(x)*ufl.dx + fS(m,m0base)*Jacob(x)*ds(1) + fS(m,m0cap)*Jacob(x)*ds(2) # total energy
dF = ufl.derivative(F, m, phi) # variation of the energy
JF = ufl.derivative(dF, m, dm) # Jacobian of the energy
problem = fem.petsc.NonlinearProblem(dF, m, J=JF)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-12
solver.max_it = 100
solver.report = True
# set initial guess
minit = [1,0,0]
minit = ufl.as_vector([fem.Constant(domain, ScalarType(minit[0])),fem.Constant(domain, ScalarType(minit[1])),
fem.Constant(domain, ScalarType(minit[2]))])
expr = fem.Expression(minit, V.element.interpolation_points())
m.interpolate(expr)
solver.nonzero_initial_guess = True
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(m)
assert(converged)
print(f"Number of interations: {n:d}")
Thanks in advance!
Pamela