It’s not clear what you’ve done for the plots in your previous post, but a correct modification of the example code in poisson-polar.py to solve on the physical annular domain would look something like the following:
from fenics import *
import math
M = 16; N = 32
mesh = RectangleMesh(Point(1,0),
Point(2,math.pi/2),
M,N)
xi = SpatialCoordinate(mesh)
r = xi[0]
theta = xi[1]
x = as_vector([r*cos(theta),
r*sin(theta)])
# Deform the mesh to the physical domain
# using a projection of the analytical
# mapping's displacement onto a piecewise
# linear vector-valued function space.
Vvec = VectorFunctionSpace(mesh,"CG",1)
u_hat = project(x-xi,Vvec)
ALE.move(mesh,u_hat)
# Solve with no change of variables.
kappa = Constant(2)
f = Constant(2)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
a = kappa*dot(grad(u),grad(v))*dx
L = f*v*dx
bc = DirichletBC(V,Constant(0),"on_boundary")
u = Function(V)
solve(a==L,u,bc)
u.rename("u","u")
File("u.pvd") << u
The resulting solution is nearly identical to the one on the reference domain (up to only slight differences from the piecewise-linear approximation of the mapping needed to mesh the curved domain geometry in physical space).
For your formulation in reference config.(or polar coordinates), you used only the determinant of deformation gradient for replacing
dx(which for this 2D problem is the area integral and not volume integral).
Note that I also transformed spatial derivatives, using the chain rule.
However, according to Nanson’s formula, we need to have other terms as well. I did try to account for those terms, but was not able to solve and was getting error in formulating them.
Nanson’s formula would only be necessary in formulations with boundary integrals (i.e., using the ds measure). It looks like you’re still solving with strongly-enforced Dirichlet boundary conditions, so there is no need to use Nanson’s formula.