Where does spatial coordinate x come from?

General
1

Hi All,
I’m confused by the spatial coordinate x in the following general boundary definition in python:

def boundary(x, on_boundary):
    return on_boundary

on_boundary is a boolean value provided by DOLFIN, but where does x come from? For example, in Section 2.2 of book ‘Solving PDEs in Python: The FEniCS Tutorial I’:

from fenics import *

mesh = UnitSquareMesh(8, 8)

V = FunctionSpace(mesh,'P', 1)

u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):

    return on_boundary

bc = DirichletBC(V,u_D,boundary)

u= TrialFunction(V)

u = TrialFunction(V)

v = TestFunction(V)

f = Constant(-6.0)

a = dot(grad(u), grad(v))*dx

L = f*v*dx

# Compute solution

u = Function(V)

solve(a == L, u, bc)

# Plot solution and mesh

plot(u)

plot(mesh)

x may come from mesh = UnitSquareMesh(8, 8), but somehow implicitly.

Now if I have a complicated formula (use formula here to differ from key word Expression) in the following for a field initialization which depends on spatial coordinate x, how can I access x?

x = SpatialCoordinate(mesh)
...
def Tini(aspect, RaT):
    q0 = aspect**(7.0/3.0)*(RaT/(2.0*math.sqrt(math.pi)))**(2.0/3.0)/(1.0+aspect**4.0)**(2.0/3.0)
    Q0 = 2.0*math.sqrt(aspect/(math.pi*q0))
    Tu = 0.5*math.erf(0.5*(1.0-x[1])*math.sqrt(q0/x[0]))
    Tl = 1.0 - 0.5*math.erf(0.5*x[1]*math.sqrt(q0/(aspect-x[0])))
    Tr = 0.5 + 0.5*Q0*math.sqrt(q0/(math.pi+math.pi*x[1]))*math.exp(-x[0]**2*q0/(4.0*x[1]+4.0))
    Ts = 0.5 - 0.5*Q0*math.sqrt(q0/(2.0*math.pi-x[1]*math.pi))*math.exp(-(aspect-x[0])**2*q0/(8.0-4.0*x[1]))

    T0 = Tu + Tl + Tr + Ts - 1.5

    return T0

This snippet will give me an error Division.__float__ returned non-float (type NotImplementedType) at Tu line just because of x[0] (verified by neglecting x[0]). Also I understand x is just a symbol, so how do I access x in this python function?

This is something fundamental that I don’t know much about, I’m still new to FEniCS. Really appreciate your help and suggestion!

2

Depending on where you useboundary, x wil differ. It will always be an array of (x,y,z) coordinates.
When used in

It will be the coordinates of the degrees of freedom in V. Each tabulated dof coordinate is sent into the boundary function to determine if it should have a DirichletBC applied to it.

Here you are mixing ufl and the math module. Depending on how you want to use Tini you need to use either only ufl or only math/numpy modules.

3

What a simple mistake! Honestly I’m still struggling with ufl and python syntax, I don’t have a clear way to tell their difference when setting up a demo. Thanks Dokken, really appreciated!