This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Current example of how to solve a PDE on a manifold

0 votes

I'm looking for an example of how to use Fenics to solve a PDE on a manifold. The only examples that I've been able to find are on

http://www.geosci-model-dev.net/6/2099/2013/

and from the question on this site

https://fenicsproject.org/qa/2897/inhomogeneous-neumann-boundary-condition-poisson-manifold?show=2897#q2897

Unfortunately, none of these codes will run with the modern version of Fenics. I'm new to Fenics so I'm not sure how to interpret the errors to be able to fix them.

If anybody could fix one of the examples to be able to run with the current version of Fenics, or (preferably) produce a minimum working example of how to solve say the poisson equation on a manifold I would greatly appreciate it.

asked May 1, 2017 by rviertel FEniCS Novice (700 points)

1 Answer

0 votes
 
Best answer

This should work, it is Arnold's example. The only things that have changed since then are:
- You have to explicitly set the Expression's degree
- You must create Mixed Spaces with a MixedElement (at least that's how it works for me)

from dolfin import *
from numpy import log

# forcing function, high degree for better approximation
f = Expression("exp(-10.*(pow(x[0], 2) + pow(x[1]-1/sqrt(2.), 2) + pow(x[2]-1/sqrt(2.), 2)))", 
degree = 4)

mesh = Mesh("hemisphere_mesh.xml")
global_normal = Expression(("x[0]", "x[1]", "x[2]"), degree = 1)
mesh.init_cell_orientations(global_normal)

# finite element mesh and spaces
V0 = FunctionSpace(mesh, "RT", 3)
V1 = FunctionSpace(mesh, "DG", 2)
R = FunctionSpace(mesh, "Real", 0)
Mel = MixedElement([V0.ufl_element(), V1.ufl_element(), R.ufl_element()])
V = FunctionSpace(mesh, Mel)
# essential boundary condition
def bdry(x, on_boundary):
    return on_boundary

bc = DirichletBC(V.sub(0), Constant((0., 0., 0.0)), bdry)

# trial and test functions
(u, p, c0) = TrialFunctions(V)
(v, q, c1) = TestFunctions(V)

# bilinear form, right hand side
a = (dot(u, v) - p*div(v) - div(u)*q + c0*q + p*c1) * dx
L = -f*q * dx

# assemble, solve
up = Function(V)
solve(a == L, up, bc)
(u, p, c0) = up.split()

I didn't know about setting the normal, so thanks for the question. Let me know in case of anything. Best regards!

answered May 3, 2017 by nabarnaf FEniCS User (2,940 points)
selected May 3, 2017 by rviertel
...