Re: [firedrake] Function at() PointNotInDomainError

24 May 2017

      Hi Lawrence,

Setting the tolerance worked. Thanks!

Best

Wei

________________________________
From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk>
Sent: 24 May 2017 13:21:22
To: firedrake
Subject: Re: [firedrake] Function at() PointNotInDomainError

On 24/05/17 12:12, Pan, Wei wrote:
...
Hey guys,
I'm having problems with the /Function.at()/ function. Below is an
example code
from __future__ import division
    from firedrake import *
    import numpy as np
mesh_base = SquareMesh(8, 8, 1.)
    mesh_hierarchy = MeshHierarchy(mesh_base, 6)
fine_msh = mesh_hierarchy[6]       /# 512 x 512 mesh/
    fine_dg_space = FunctionSpace(fine_msh, "DG", 1)
x, y = SpatialCoordinate(fine_dg_space)
    q0 = Function(fine_dg_space).interpolate(x)
x_array = np.array([_x / 512.0 for _x in range(513)])  # coordinates
    # x_array[0] += 1. / 1000000.0 # a slight nudge away from boundary
    also fails
for i in x_array:
        for j in x_array:
            print [i, j], q0.at(np.array([i, j]))
My domain is [0, 1], mesh size 512x512. I keep on getting the
following PointNotInDomainError when I run the above, where 0.43164062
correspond to 221 of 512. What am I doing wrong?
[0.0, 0.427734375] 0.0
    [0.0, 0.4296875] 0.0
    [0.0, 0.431640625]
    Traceback (most recent call last):
      File "/home/wpan1/PycharmProjects/MyFirstProject/test_at.py",
    line 22, in <module>
        print [i, j], q0.at(np.array([i, j]))
      File
    "/home/wpan1/Documents/firedrake/src/firedrake/firedrake/function.py",
    line 574, in at
        raise PointNotInDomainError(self.function_space().mesh(),
    points[i].reshape(-1))
    firedrake.function.PointNotInDomainError: domain <Mesh #7> does
    not contain point [ 0.          0.43164062]
The geometric predicates used to decide if a point is in a cell are
all done using floating point, so what might be happening is that this
point lies on the boundary between two cells, but to the specified
tolerance, it lies outside from both sides.

Does it help to relax the tolerance with which we check?

Say, for example:

q0.at(..., tolerance=1e-6)

Cheers,

Lawrence