values of Function at specific points

On 08/28/2012 05:50 PM, Andrew E Slaughter wrote:
> New question #207055 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/207055
>
> I have a simple question, that is due to my inexperience with both Python and FEniCS. The following code causes the error shown. I am expecting the result from eval should be a scalar value for the Function evaluated at the specified point, but the eval function requires a numpy.array(float) according to the documentation, which I can not seem to initialize correctly.
>
> In the end, I want an array of the values of phi for each vertex in the cell, from which I can perform a test to assign the cell to a subdomain. I would appreciate any help or suggestions getting this working.
>
> --- ERROR -----------------------------
> TypeError: contiguous numpy array of 'double' expected. Make sure that the numpy array is contiguous and uses dtype='d'.
>
> --- CODE ------------------------------
> from dolfin import *
> import numpy
>
> # Create mesh and phi Function
> mesh = Rectangle(-1.0, -1.0, 1.0, 1.0, 1, 1)
> phi_space = FunctionSpace(mesh, "DG", 1)
> phi_init = Expression("x[0]")
> phi_old = Function(phi_space)
> phi_old.interpolate(phi_init)
>
> # Locate the interface, label cells as subdomain 1
> for c in cells(mesh):
> for v in vertices(c):
> values = numpy.array([0])
> phi_old.eval(values,v.point())
>

If you instead use:

  values = numpy.array([0], dtype='d') # or dtype=numpy.float_
                                       # or use [0.] instead of [0]

and do not access the coordinates from the vertex, as these are not
correctly wrapped to Python (This is a bug). Instead use:

  phi_old.eval(values, mesh.coordinates()[v.index()])

Also keep in mind that you are evaluating a DG Function which
potentially has different values at each vertex dependent on what cell
you are in. Then you can use eval_cell instead:

  phi_old.eval_cell(values, mesh.coordinates()[v.index()], cell)

Johan

Thanks Johan Hake, that solved my question.

I am also implementing this behavior in C++, which does not have a eval_cell member function, but has an overloaded cell. The overloaded eval functions require a ufc::cell and/or a dolphin::cell. What is the difference? Based on the comments above regarding DG elements, I am assuming I need to include the cell somewhere to get the correct value.

Thanks for the help,
Andrew

On 08/29/2012 04:26 PM, Andrew E Slaughter wrote:
> Question #207055 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/207055
>
> Status: Solved => Open
>
> Andrew E Slaughter is still having a problem:
> I am also implementing this behavior in C++, which does not have a
> eval_cell member function, but has an overloaded cell.

In Python the version with cell argument is renamed to eval_cell.

> The overloaded
> eval functions require a ufc::cell and/or a dolphin::cell. What is the
> difference? Based on the comments above regarding DG elements, I am
> assuming I need to include the cell somewhere to get the correct value.

Yes,

Johan

> Thanks for the help,
> Andrew
>

What is the difference between the ufc::cell and dolphin::cell? The three argument version of eval in C++ requires a ufc::cell, but I am looping over a dolphin::cell. How can I use the three argument version in C++?

On Wed, Aug 29, 2012 at 04:01:11PM -0000, Andrew E Slaughter wrote:
> Question #207055 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/207055
>
> Status: Answered => Open
>
> Andrew E Slaughter is still having a problem:
> What is the difference between the ufc::cell and dolphin::cell? The
> three argument version of eval in C++ requires a ufc::cell, but I am
> looping over a dolphin::cell. How can I use the three argument version
> in C++?

You can instantiate a ufc::cell from a dolfin::cell by using the class
UFCCell. If you require it multiple times, as in a loop. Create one
UFCCell object outside the loop and then call the update function
inside the loop.

--
Anders

I came here with a similar question, and the answers so far do not appear to work with the latest version of FEniCS. I am using the latest docker image. I get the error mentioned in the original question:

TypeError: contiguous numpy array of 'double' expected. Make sure that the numpy array is contiguous, with 1 dimension,
and uses dtype=float_.

when running the following script:

import fenics
import numpy

mesh = fenics.RectangleMesh(fenics.Point((-1.0, -1.0)), fenics.Point((1.0, 1.0)), 1, 1)
phi_space = fenics.FunctionSpace(mesh, "DG", 1)
phi_init = fenics.Expression("x[0]")
phi_old = fenics.Function(phi_space)
phi_old.interpolate(phi_init)

for c in cells(mesh):
    for v in vertices(c):
        values = numpy.array([0.], dtype=numpy.float_)
        phi_old.eval(values, v.point())

This is related to my comment on the question here: https://fenicsproject.org/qa/596/setting-condition-for-mesh-refinement, i.e. this provides some history explaining why I'm trying to use eval or eval_cell.