closest_cell function throwing exception

On Fri, Dec 10, 2010 at 12:20:04AM -0000, Jack wrote:
> New question #137265 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/137265
>
> This statement
>
> index = mesh.closest_cell(p);
>
> terminates the program with the error:
>
> terminate called after throwing an instance of 'std::runtime_error'
> what(): *** Error: The function 'compute' has not been implemented (in /Users/JJ//FEniCS/src/dolfin/dolfin/mesh/IntersectionOperatorImplementation.h line 265).
> [AttackinhgChess-3:22966] *** Process received signal ***
> [AttackinhgChess-3:22966] Signal: Abort trap (6)
> [AttackinhgChess-3:22966] Signal code: (0)
> [AttackinhgChess-3:22966] [ 0] 2 libSystem.B.dylib 0x00007fff82b0367a _sigtramp + 26
> [AttackinhgChess-3:22966] [ 1] 3 ??? 0x000000014cee0210 0x0 + 5585633808
> [AttackinhgChess-3:22966] [ 2] 4 libstdc++.6.dylib 0x0000000101b3c365 _ZN9__gnu_cxx27__verbose_terminate_handlerEv + 293
> [AttackinhgChess-3:22966] [ 3] 5 libstdc++.6.dylib 0x0000000101ad084a _ZN10__cxxabiv111__terminateEPFvvE + 10
> [AttackinhgChess-3:22966] [ 4] 6 libstdc++.6.dylib 0x0000000101ad08a3 _ZSt9terminatev + 19
> [AttackinhgChess-3:22966] [ 5] 7 libstdc++.6.dylib 0x0000000101ad0926 __cxa_throw + 102
> [AttackinhgChess-3:22966] [ 6] 8 libdolfin.0.dylib 0x0000000100070d34 _ZNK6dolfin6Logger5errorESs + 132
> [AttackinhgChess-3:22966] [ 7] 9 libdolfin.0.dylib 0x0000000100070b7b _ZN6dolfin5errorESsz + 283
> [AttackinhgChess-3:22966] [ 8] 10 libdolfin.0.dylib 0x00000001001a8d83 _ZNK6dolfin36IntersectionOperatorImplementation_dINS_15TetrahedronCellEN4CGAL16Simple_cartesianIdEEE22closest_point_and_cellERKNS_5PointE + 83
> [AttackinhgChess-3:22966] [ 9] 11 libdolfin.0.dylib 0x000000010018e7bb _ZNK6dolfin36IntersectionOperatorImplementation_dINS_15TetrahedronCellEN4CGAL16Simple_cartesianIdEEE12closest_cellERKNS_5PointE + 27
> [AttackinhgChess-3:22966] [10] 12 colour_mesh 0x000000010000f6a7 main + 1351
> [AttackinhgChess-3:22966] [11] 13 colour_mesh 0x000000010000f0c4 start + 52
> [AttackinhgChess-3:22966] [12] 14 ??? 0x0000000000000001 0x0 + 1
> [AttackinhgChess-3:22966] *** End of error message ***
>
> Is it not implemented for 3D?

Correct. Unfortunately, the functionality we depend on is not
implemented for tetrahedra in CGAL.

Andre Massing might be able to give you some more details.

--
Anders

I want to interpolate a function of other functions on a three dimensional function space by subclassing dolfin::Expression.

When I do this in a 3d mesh (created by gmsh or with UnitSphere), I get the error "Unable to evaluate function at given point (not inside domain). Set option 'allow_extrapolation' to allow extrapolation"

If I set "allow_extrapolation" to true, I get the same error as above.

Does your answer mean, we have to wait until the functionality is implemented in CGAL?
Do you know, if they even have plans to implement this feature?

At this moment, I use CGAL 3.7.

Andre Massing is working on this, so stay put.

Martin

On 30 June 2011 16:21, Andreas Naumann
<email address hidden> wrote:
> Question #137265 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/137265
>
> Andreas Naumann requested more information:
> I want to interpolate a function of other functions on a three
> dimensional function space by subclassing dolfin::Expression.
>
> When I do this in a 3d mesh (created by gmsh or with UnitSphere), I get
> the error "Unable to evaluate function at given point (not inside
> domain). Set option 'allow_extrapolation' to allow extrapolation"
>
> If I set "allow_extrapolation" to true, I get the same error as above.
>
> Does your answer mean, we have to wait until the functionality is implemented in CGAL?
> Do you know, if they even have plans to implement this feature?
>
> At this moment, I use CGAL 3.7.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Hi!

On 06/30/2011 04:21 PM, Andreas Naumann wrote:
> Question #137265 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/137265
>
> Andreas Naumann requested more information:
> I want to interpolate a function of other functions on a three
> dimensional function space by subclassing dolfin::Expression.
>
> When I do this in a 3d mesh (created by gmsh or with UnitSphere), I get
> the error "Unable to evaluate function at given point (not inside
> domain). Set option 'allow_extrapolation' to allow extrapolation"
>
> If I set "allow_extrapolation" to true, I get the same error as above.
>
> Does your answer mean, we have to wait until the functionality is implemented in CGAL?
> Do you know, if they even have plans to implement this feature?
>
> At this moment, I use CGAL 3.7.

The resolution of that issue depends on a couple of things.
Maybe you could give us some more information, as in
-which DOLFIN version are you using
-minimal example where the error occurs.

There are 2 comments I can give right away:
1) allow_extrapolation depends on different features, like closest point
and cell computations, which is incomplete in 3D due to some missing
functionality in CGAL we we plan to add, see blueprint

https://blueprints.launchpad.net/dolfin/+spec/closest-point-3d

2) In some recent version of DOLFIN Marie has implemented a workaround
for that functionality in order to do the extrapolation in 3D as well
(IIRC) which works fine. Depending on your DOLFIN version, you may
consider an update. Marie knows the versions where it has been added.

HTH,
Andre

Hi @all:

thank you for the fast answer.

dolfin version: 0.9.10

Example:
ufl-file: Poisson.ufl (from the examples, triangle replaced by tetrahedron)
cpp:
#include <dolfin.h>
#include "Poisson3D.h"

using namespace dolfin;
struct Evaler : Expression {
        Function& f;
        Evaler(Function& f):
                f(f)
        {}
        void eval(Array< double >& vals, const Array<double>& x) const {
                f.eval(vals,x);
        }
};

Mesh mesh = UnitSphere(20);
typedef Poisson3D::FunctionSpace Space;

parameters["allow_extrapolation"] = true;

 Constant c(1.0);
        Function f1(space);
        f1=c;
        File testFile("meshtest.pvd");
        testFile << f1;
        Function f2(space);
        Evaler ev(f1);
        f2 = ev; //will throw "Unable to evaluate function at given point (not inside domain). Set option 'allow_extrapolation' to allow extrapolation'"
        return 0;

Regards,
Andreas

On 06/30/2011 04:56 PM, Andreas Naumann wrote:
> Question #137265 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/137265
>
> Andreas Naumann proposed the following answer:
> Hi @all:
>
> thank you for the fast answer.
>
> dolfin version: 0.9.10
>
> Example:
> ufl-file: Poisson.ufl (from the examples, triangle replaced by tetrahedron)
> cpp:
> #include<dolfin.h>
> #include "Poisson3D.h"
>
> using namespace dolfin;
> struct Evaler : Expression {
> Function& f;
> Evaler(Function& f):
> f(f)
> {}
> void eval(Array< double>& vals, const Array<double>& x) const {
> f.eval(vals,x);
> }
> };
>
>
> Mesh mesh = UnitSphere(20);
> typedef Poisson3D::FunctionSpace Space;
>
> parameters["allow_extrapolation"] = true;
>
> Constant c(1.0);
> Function f1(space);
> f1=c;
> File testFile("meshtest.pvd");
> testFile<< f1;
> Function f2(space);
> Evaler ev(f1);
> f2 = ev; //will throw "Unable to evaluate function at given point (not inside domain). Set option 'allow_extrapolation' to allow extrapolation'"
> return 0;

Is your test example complete? I am missing for example the definition
of the variable "space".

Andre

>
> Regards,
> Andreas
>

sorry, my selection started and ended from the wrong lines. the full c++ main is:

int main() {
        //Mesh mesh("zylinder.xml.gz");
        Mesh mesh = UnitSphere(20);
        typedef Poisson3D::FunctionSpace Space;
        Space space(mesh);

        parameters["allow_extrapolation"] = true;
        Constant c(1.0);
        Function f1(space);
        f1=c;
        File testFile("meshtest.pvd");
        testFile << f1;
        Function f2(space);
        Evaler ev(f1);
        f2 = ev; //will throw "Unable to evaluate function at given point (not inside domain). Set option 'allow_extrapolation' to allow extrapolation'"
        return 0;
}

On 06/30/2011 05:31 PM, Andreas Naumann wrote:
> Question #137265 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/137265
>
> Andreas Naumann proposed the following answer:
> sorry, my selection started and ended from the wrong lines. the full c++
> main is:
>
>
> int main() {
> //Mesh mesh("zylinder.xml.gz");
> Mesh mesh = UnitSphere(20);
> typedef Poisson3D::FunctionSpace Space;
> Space space(mesh);
>
> parameters["allow_extrapolation"] = true;
> Constant c(1.0);
> Function f1(space);
> f1=c;
> File testFile("meshtest.pvd");
> testFile<< f1;
> Function f2(space);
> Evaler ev(f1);
> f2 = ev; //will throw "Unable to evaluate function at given point (not inside domain). Set option 'allow_extrapolation' to allow extrapolation'"
> return 0;
> }
>

Great, last request:) If the cylinder mesh is not to big, could you
supply that as well? Then we have a complete test case at hand. Thanks a
lot!
Andre

i uploaded it to
http://www.file-upload.net/download-3549438/zylinder.xml.gz.html

it's size is approx 174kB.

I hate tests, which depend on files, that's why I chose the UnitSphere, where the error is also thrown

Andreas

I tried the current tarball (dolfin 0.9.11) and got the same error :(
Do I have to set some additional flags?

Andreas

Hi!

On 07/02/2011 09:41 AM, Andreas Naumann wrote:
> Question #137265 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/137265
>
> Andreas Naumann proposed the following answer:
> I tried the current tarball (dolfin 0.9.11) and got the same error :(
> Do I have to set some additional flags?
>
> Andreas
>

I am not aware of any additional flags which might be needed here, but
Anders or Marie will know more about that.
Unfortunately I have not time to look into this before Monday/Tuesday (I
am moving this week-end), but just one quick question:
At which points to you want to evaluate the function?
Could you give an example?
If you want to evaluate the function outside the region (as in the point
is not within the mesh domain modulus some precision errors) then I
don't know whether Marie's extension was aiming at that case as well. I
roughly remember that is was meant to evaluate functions at points which
were just slightly outside because of finite precision.

HTH
Andre

I don't even want to evaluate outside, simply interpolate a given expression which depends on a (dolfin) function. The points, dolfin evaluates the expression are slightly outside the domain.

After installing the newer dolfin, I hoped to get at least another behavior/another error message. That's why, I thought the extension has to be enabled in some manner.

I'll ask Marie and wish you a good weekend.

Andreas

Sorry for the late response.

Could you try overloading the

  void eval(Array<double>& values, const Array<double>& x,
            const ufc::cell& ufc_cell) const

method and call f.eval(values, x, cell) instead? The mixed-poisson demo provides an example for this method.

@Marie:

Thank you for the hint.
In my small test, this method works.
I'll test it in the whole program later.

I think, the problem is solved now.