# Use the result of a PDE for further calculations

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:

1. 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?

2. 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).

3. 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)), # origin
(1, lambda x: np.isclose(x, Rbase) & np.greater_equal(x, np.pi/2)), # base
(2, lambda x: np.isclose(x, Rcap) & np.less(x, np.pi/2)), # cap
(3, lambda x: np.isclose(x, minq)),
(4, lambda x: np.isclose(x, 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)),fem.Constant(domain, ScalarType(minit)),
fem.Constant(domain, ScalarType(minit))])
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}")


Pamela

Hi again. I’ve been able to solve point 3 by doing fem.assemble_scalar(fem.form(F)) after finding m.

Any suggestions on the other two points? I really need to be able to transform from spherical to Cartesian coordinates.

Hi pamgur,

Just my thoughts on point 1. You probably can evaluate and save a list of the spherical coordinates and the values, then calculate its corresponding value in Cartesian coordinates using numpy or matlab yourself. The transform rules are available on the internet and many textbooks. In dolfinx the evaluation is like

‘function.eval(x0, cells)’

With matplotlib, you can easily visualize a list of coordinates and the function values.

Cheers
Carrot

1 Like

Hi Carrot,
Thanks for your answer. Could you please clarify a bit how to perform the evaluation? I’ve tried the following after solving for m:

x0 = [0.5, 1]
tree = geometry.BoundingBoxTree(domain, gdim)
cells = geometry.compute_colliding_cells(domain, tree, x0)
m.eval(x0, cells)


But I get the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In , in <cell line: 3>()
1 x0 = [0.5, 1]
2 tree = geometry.BoundingBoxTree(domain, gdim)
----> 3 cells = geometry.compute_colliding_cells(domain, tree, x0)
5 m.eval(x0, cells)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/geometry.py:63, in compute_colliding_cells(mesh, candidates, x)
50 def compute_colliding_cells(mesh: Mesh, candidates: AdjacencyList_int32, x: numpy.ndarray):
51     """From a mesh, find which cells collide with a set of points.
52
53     Args:
(...)
61         collide with the ith point
62     """
---> 63     return _cpp.geometry.compute_colliding_cells(mesh, candidates, x)

TypeError: compute_colliding_cells(): incompatible function arguments. The following argument types are supported:
1. (mesh: dolfinx.cpp.mesh.Mesh, candidate_cells: dolfinx.cpp.graph.AdjacencyList_int32, points: numpy.ndarray[numpy.float64]) -> Union[dolfinx.cpp.graph.AdjacencyList_int32, numpy.ndarray[numpy.int32]]

Invoked with: <dolfinx.mesh.Mesh object at 0x7fa3dfd57c90>, [0 0 0 ]->[1 3.14159 0 ]
{[0 0 0 ]->[0.967178 1.50933 0 ]
{[0 0 0 ]->[0.782479 0.846396 0 ]
{[0 0 0 ]->[0.718235 0.438816 0 ]
{[0 0 0 ]->[0.494683 0.438692 0 ]
{[0 0 0 ]->[0.318892 0.438692 0 ]
{[0 0 0 ]->[0.318892 0.208535 0 ]
{[0 0 0 ]->[0.182887 0.208535 0 ]
{[0 0 0 ]->[0.181229 0.0981748 0 ]
{[0 0 0 ]->[0.0937538 0.0981748 0 ]
{[0 0 0 ]->[0.0937538 0.06674 0 ]
leaf containing entity (106),
[0 0 0 ]->[0.0666976 0.0981748 0 ]
leaf containing entity (121)}
,
[...many many more lines of stuff like this...]
}
}
}
}
}
}
}
}
}
}
, [0.5, 1]


Hi pamgur,

Sorry for the late reply.

I am not sure about your case but the general usage of eval function can be found from the tutorial here.

Please have a try and let me know if it works for you.

Best,
Carrot