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

"Unable to evaluate function at point" error

0 votes

Hi
I want to evaluate a variable at a particular point of a domain. My domain is a rectangle. In a part of my code where I have defined my mesh (inside a class) I have moved the mesh coordinates upward (in the y-direction):

def mesh(self):
    .
    .
    m = RectangleMesh(Point(0,0), Point(1, 1),self.Nx, self.Ny)
    x = m.coordinates()
    .
    .
    for x in m.coordinates(): #moving mesh in the y-direction
      x[1] += 0.000002
    return m

I want to evaluate the variable "c" (in my problem) at the upper edge (left corner).

print c(0.0,202.0 * (10**(-6)))

I get this error:

*** Reason:  The point is not inside the domain...

While the y=202.0 * (10**(-6)) is obviously inside the domain!! (It is the upper edge of the domain)!
In addition when I do:

print c(0.0,201.9999 * (10**(-6)))

I do not get the error anymore and it gives me the value of "c" at this point but I need to get this value at the coordinate:
(0.0,202.0 * (10**(-6)))
Here is the complete code:

from dolfin import *
from numpy import cos, pi
import numpy as np

D = 1. * (10**(-9))

R = 8.31

C_0 = 1.

epsilon = 0.025

Temp = 293.0

mu = D / (R*Temp)

F_N = 96485.34

Charge_num=Constant(1.)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.03)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.000002)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], (202. * (10**(-6))))

left = Left()
top = Top()
right = Right()
bottom = Bottom()

class MyMesh:
    def __init__(self, lx,ly, Nx,Ny):
        self.lx = lx
        self.ly = ly
        self.Nx = Nx
        self.Ny = Ny

    def mesh(self):
        m = RectangleMesh(Point(0,0), Point(1, 1),self.Nx, self.Ny)
        x = m.coordinates()

        x[:,:] = (x[:,:] - 0.5) * 2.
        x[:,:] = 0.5 * (cos(pi * (x[:,:] - 1.) / 2.) + 1.)

        x[:,0] = x[:,0]*self.lx
        x[:,1] = x[:,1]*self.ly
        for x in m.coordinates():
          x[1] += 0.000002
        return m

mesh0=MyMesh(0.03,200.*(10 ** (-6)), 500,30)
mesh1 = mesh0.mesh()

boundaries = FacetFunction("size_t", mesh1)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)

CG1_elem = FiniteElement("CG", mesh1.ufl_cell(), 1)

V_c = FunctionSpace(mesh1, CG1_elem)

V_phi = FunctionSpace(mesh1, CG1_elem)

W_elem = MixedElement([CG1_elem, CG1_elem])

W = FunctionSpace(mesh1, W_elem)

z = Function(W)
dz=TrialFunction(W)

c, phi = split(z)

(v_1, v_2) = TestFunctions(W)

dt = 0.01
t = 0
T = 0.1

# Previous solution
W_const = Constant(1.0)
C_previous = interpolate(W_const, V_c)

bc_top = DirichletBC(W.sub(1), Constant(2.0), boundaries, 1)

bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)

bcs = [bc_top, bc_bottom]

# Variational form For nonlinear solver
F = dot(grad(phi), grad(v_2)) * dx \
    + c * v_1 * dx \
    - (F_N / epsilon) * c * v_2 * dx \
    + dt * D * dot(grad(c), grad(v_1)) * dx \
    + dt * Charge_num * mu * F_N * c * dot(grad(phi), grad(v_1)) * dx \
    - C_previous * v_1 * dx + (F_N / epsilon) * C_0 * v_2 * dx\

while t <= T:

    J = derivative(F, z, dz)
    problem = NonlinearVariationalProblem(F, z, bcs,J)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()

    (c, phi) = z.split(True)
    C_previous.assign(c)

    t += dt
    plot(c,key="c", interactive=False)

print c(0.0,202.0 * (10**(-6)))

Any idea how to fix this issue?
Thanks!

asked Jun 20, 2017 by Leo FEniCS Novice (840 points)
edited Jun 20, 2017 by Leo

I don't know if this will help you but you can use DOLFIN_EPS:

print c(0.0,202.0 * (10**(-6)) - DOLFIN_EPS)

Thank you so much. It worked like a charm!
The thing is, my problem is bigger than that! I thought it would be better to open a new topic to discuss it.

...