Right way to evaluate a function at the vertices

Could you attach the code?

Kent

On 9 June 2012 20:15, Maurice van Leeuwen <
<email address hidden>> wrote:

> Question #199961 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/199961
>
> Description changed to:
> 1.Introduction
> For a linear elastic material I want to determine the displacements at the
> vertices.
> Using some custom code geometry (polygons) are generated, then meshed
> using Tetgen.
> This mesh is then converted to a Dolfin XML Mesh (which is undocumented?),
> this is done
> with some custom code (for convenience) but is otherwise identical to
> dolfin-convert.
>
> The displacements are determined using the usual u(mesh.coordinates()[n])
> method that
> is presented in the Fenics book. For some meshes evaluating the
> displacement results
> in an error for one or more vertices. These meshes are all based on the
> same simple geometry,
> errors occur when the mesh is finer (more vertices).
>
>
> As usual the displacements can be found by solving the linear variational
> problem:
>
> mesh = Mesh("mesh.xml") #dolfin xml mesh
> V = VectorFunctionSpace(mesh, "Lagrange", 2)
>
> ## As usual boundary conditions, forces are applied,
> ## equations for the linear elastic problem are specified,
> ## nothing special happening. Omitted for readability.
>
> u = Function(V)
> problem = LinearVariationalProblem(a, L, u, bcs=bc)
> solver = LinearVariationalSolver(problem)
> solver.parameters["symmetric"] = True
> solver.solve()
>
> After solving, the function u contains the displacements for
> functionspace(V).
> Now there are a couple of ways to determine the displacements at the
> vertices.
>
> 2. Function evaluation
> As said, the function evalution method as presented in the Fenics book
> (pg.192)
> can be used. This is done for all vertex coordinates in the mesh:
>
> meshCoordinates = mesh.coordinates()
> for n in range(meshCoordinates.shape[0]):
> displacement = u(meshCoordinates[n])
>
> 3. Errors
> For fine meshes (based on the same simple geometry) with more vertices
> (more than
> 400, but might be different for other geometry) the following error occurs
> when
> evaluating u:
>
> *** Error: Unable to evaluate function at point.
> *** Reason: The point is not inside the domain. Consider setting
> "allow_extrapolation" to allow extrapolation.
> *** Where: This error was encountered inside Function.cpp.
>
> After setting parameters["allow_extrapolation"] = True:
>
> *** Error: Unable to perform call to DOLFIN function %s.
> *** Reason: The function compute has not been implemented (in compute
> line 87118088).
> *** Where: This error was encountered inside log.h.
>
> This problem occurs always for the same vertices, though inspection of the
> mesh doesn't
> indicate any reasons why that vertex is problematic.
>
> 4. Possible PVD / DOF method
> Displacements at the vertices can also be obtained by exporting to pvd:
> pvd = File("displacements.pvd")
> pvd << u
>
> Similar for dolfin xml:
> xml = File("displacements.xml")
> xml << u
>
> Exporting to pvd goes without errors. So I expect the file to contain the
> problematic vertex displacement, although I haven't checked.
>
> Exporting to xml gives a file with displacements at the DOFs, also without
> issuing an error. I believe this is similar to u.vector()[n], but I'm
> unsure
> how to calculate vertex values from DOF values.
>
> 5.
> Can anyone reflect on why u(mesh.coordinates()[n]) isn't working properly,
> and how to properly evaluate the vertex displacements? Also it would be
> equally useful if I know how to convert from the DOF to vertex values.
>
> Thank you!
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I've place the necessary files on pastebin. You can also check the imagebin links to see the difference between the two meshes. For the conflicting mesh-fine.xml I've added red dots to indicate the problematic vertices. Please notice that export to xml and pvd seems to be without problems.

linearelastic.py
http://pastebin.com/RGCnCnB8
uncomment line 3 to enable extrapolation (gives another error)
switch line 6 and 7 to evaluate the coarse mesh (no errors)

mesh-fine.xml
http://pastebin.com/6vQQrrmV
532 vertices
1308 cells
http://imagebin.org/215834
http://imagebin.org/215835

mesh-coarse.xml
http://pastebin.com/72nrbcxm
468 vertices
1165 cells
http://imagebin.org/215837
http://imagebin.org/215838

Additionally I can confirm that this method:
  pvd = File("displacements.pvd")
  pvd << u
does indeed contain correct values for u, for the vertices that throw an error with the other method.

> 3. Errors
> For fine meshes (based on the same simple geometry) with more vertices (more than
> 400, but might be different for other geometry) the following error occurs when
> evaluating u:
>
> *** Error: Unable to evaluate function at point.
> *** Reason: The point is not inside the domain. Consider setting "allow_extrapolation" to allow extrapolation.
> *** Where: This error was encountered inside Function.cpp.
>
> After setting parameters["allow_extrapolation"] = True:
>
> *** Error: Unable to perform call to DOLFIN function %s.
> *** Reason: The function compute has not been implemented (in compute line 87118088).
> *** Where: This error was encountered inside log.h.
>
> This problem occurs always for the same vertices, though inspection of the mesh doesn't
> indicate any reasons why that vertex is problematic.

That error message looks strange. What version of dolfin do you use? The
%s at the end of the first line was removed some 200 commits ago.

Johan

It seems to work fine with dolfin 1.0.

Kent

On 10 June 2012 14:11, Maurice van Leeuwen <
<email address hidden>> wrote:

> Question #199961 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/199961
>
> Status: Answered => Open
>
> Maurice van Leeuwen is still having a problem:
> I've place the necessary files on pastebin. You can also check the
> imagebin links to see the difference between the two meshes. For the
> conflicting mesh-fine.xml I've added red dots to indicate the
> problematic vertices. Please notice that export to xml and pvd seems to
> be without problems.
>
> linearelastic.py
> http://pastebin.com/RGCnCnB8
> uncomment line 3 to enable extrapolation (gives another error)
> switch line 6 and 7 to evaluate the coarse mesh (no errors)
>
> mesh-fine.xml
> http://pastebin.com/6vQQrrmV
> 532 vertices
> 1308 cells
> http://imagebin.org/215834
> http://imagebin.org/215835
>
> mesh-coarse.xml
> http://pastebin.com/72nrbcxm
> 468 vertices
> 1165 cells
> http://imagebin.org/215837
> http://imagebin.org/215838
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I should note I'm using the 1.0 release for windows. I believe the uBLAS solver is used. Maybe the problem is in the used platform or used solver?

I suspected the problem may lie in the uBLAS solver - so I switched to Ubuntu Linux to have the PETSc solver.
Installed with:
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt-get update
sudo apt-get install fenics
sudo apt-get dist-upgrade

So it should be the latest version? Still having the extrapolation error (or the %s error when allowing extrapolation).

It would be interesting to know what configuration you are on Kent! I'll see if I can install a more recent source
like Johan Hake mentioned..

Yes, that should give you the latest stable version. Which version of Ubuntu do you have? It works for me with Ubuntu 12.04.

I'm also on Ubuntu 12.04 (32-bit). I've also tried it with the latest development build (sudo add-apt-repository ppa:fenics-packages/fenics-dev). I'm still having an error running the script. I wonder what differs from our configuration if the error doesn't occur to you?

With extrapolation enabled, running the script for the fine mesh:
*** Error: Unable to perform call to DOLFIN function.
*** Reason: The function compute has not been implemented (in /build/buildd/dolfin-1.0.0/dolfin/intersection/IntersectionOperatorImplementation.h line 294).
*** Where: This error was encountered inside log.h.

So it seems related to this problem: https://answers.launchpad.net/dolfin/+question/199448

Anyway, how does the pvd export make it work? Can I compose the vertex values from u.vector() someway?
Exporting to pvd and then importing again would make it work, but that isn't playing nice..

I got the same error on a 32 bits Ubuntu 12.04 machine.

Thanks Anders Logg, that solved my question.

Still unclear to my why the interpolation failed. But I've come up with three methods to find the vertex values:

(this one didn't actually work for me, probably implemented it wrong. Just posting it for reference:)
  import numpy
  displacements = numpy.zeros(mesh.num_vertices(), dtype='d')
  u.compute_vertex_values(displacements, mesh)

Another way:
  displacements = dolfin.cpp.FunctionPlotData(u, mesh)
  displacements.vertex_values()

Or interpolating u from a higher order Lagrange functionspace to a linear Lagrange space, as seen in 1.1.6, and 1.1.11 in the book:
  V = VectorFunctionSpace(mesh, "Lagrange", 2)
  V2 = VectorFunctionSpace(mesh, "Lagrange", 1)
  u = Function(V)
  #<insert solve stuff here>, then interpolate to linear lagrange space:
  uV2 = interpolate(u, V2)
  uV2v = uV2.vector()
  if(len(uV2v) == 3*mesh.num_vertices()):
 uV2a = numpy.reshape(uV2v, (mesh.num_vertices(), 3), order='F')

uV2a now is a numpy array of shape [num_vertices(), 3]