It should also be pointed out that you might (depending on circumstances) be making life overly complicated for yourself.

The mesh coordinates are just a Function themselves, so if you have a CG1 Function then u.dat.data_ro is exactly the Function evaluated at the mesh vertices, in the same order as mesh.coordinate.dat.data_ro . Furthermore, you can achieve the same thing for the nodes of any CG or DG space by reinterpolating the coordinates. See the instructions at:

https://www.firedrakeproject.org/interpolation.html#interpolation-from-external-data

(admittedly those instructions are for interpolating into Firedrake rather than out of, but the same principle applies)

Regards,

David

From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk>
Reply-To: firedrake <firedrake@imperial.ac.uk>
Date: Friday, 26 October 2018 at 13:07
To: firedrake <firedrake@imperial.ac.uk>
Subject: Re: [firedrake] .at() and PointNotInDomainError

Ok thanks Lawrence.

---

De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Lawrence Mitchell <wence@gmx.li>
Envoyé : vendredi 26 octobre 2018 11:12
À : firedrake
Objet : Re: [firedrake] .at() and PointNotInDomainError

Dear Floriane,

> On 26 Oct 2018, at 11:06, Floriane Gidel [RPG] <mmfg@leeds.ac.uk> wrote:
>
> Dear all,
>
> When trying to evaluate a function u on some mesh points with
>
> mesh = Mesh('coarse/Omega_coarse.msh')
> V = FunctionSpace(mesh, "CG",1)
> V_Vec = VectorFunctionSpace(mesh, "CG", 1)
> u = Function(V)
> points = Function(V_vec).interpolate(SpatialCoordinate(mesh))
> dat = u.at(points.dat.data_ro)
>
> I get the following error:
>
> File "lin_dyn_unicell.py", line 117, in <module>
>     dat = u.at(points.dat.data_ro)
>   File "/Users/gidel/Documents/Firedrake/firedrake/src/firedrake/firedrake/function.py", line 585, in at
>     raise PointNotInDomainError(self.function_space().mesh(), points[i].reshape(-1))
> firedrake.function.PointNotInDomainError: domain <Mesh #1> does not contain point [200.          66.66666667]
>
>
> How is that possible while both function spaces are defined on the same mesh?

You have run into the delightful world of computational geometry. Unfortunately points might not be determined to be in any cell if they fall close to cell boundaries.

You can reduce the tightness of the tolerance used to check if a point is in a cell by saying:

data = u.at(points, tolerance=1e-6)

Thanks,

Lawrence
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake