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