dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #09274
Re: pydolfin
It works if you replace
Function(element, mesh, h_e)
by
Function(element, mesh, Vector())
You should get an error message when trying to create a Function from
a numpy array (unless Martin has done some tricks for initialization
of vector constants that I have missed).
Martin?
--
Anders
On Thu, Aug 21, 2008 at 08:44:32PM +0200, Evan Lezar wrote:
> Hi there
>
> I have implemented a simple test case and have attached it. In writing it, I
> made a discovery. The vector that is passed as a list of coefficients to the
> function constructor (h_e in the example) needs to be a multidimensional array
> with the first dimension 1. If one swaps the dimensions (1, num_dofs) on line
> 26, then the test crashes when field.eval is called.
>
> However, according to me, the eval function does not work as expected. Since
> the list of coefficients is zero and the function should be evaluated as the
> weighted sum of the basis functions using the coefficients as weights, this
> means that the function value returned for any position within the mesh should
> be zero. However, the first component (let us say x) of H is always -1 and the
> second component is not changed and retains the value that was assigned to it
> before the call to evaluate.
>
> I will see about getting a C++ test up and running to test the same thing.
>
> Please let me know if there are any other questions.
> Evan
>
> On Thu, Aug 21, 2008 at 6:06 PM, Evan Lezar <evanlezar@xxxxxxxxx> wrote:
>
>
>
> On Thu, Aug 21, 2008 at 5:49 PM, Ola Skavhaug <skavhaug@xxxxxxxxx> wrote:
>
> Evan Lezar skrev den 21/08-2008 f lgende:
> > Thanks for the speedy replies from all of you.
> >
> > Firstly platform information
> >
> > I am running ubuntu hardy, gcc version 4.2.3, swig 1.3.33, python 2.5
> with
> > release versions of numpy and scipy. I am building dolfin directly
> from
> > source with mpi, but have not enabled petsc or any of the other
> optional
> > packages.
> >
> > What other information do you require?
> >
> > In terms of my specific problem, see the following code snippet:
> >
> > from dolfin import *
> > import numpy as N
> >
> > mesh = UnitSquare(4,4)
> > element = FiniteElement("Nedelec", "triangle", 0)
> >
> > ### Code to assemble system matrices and solve an eigenvalue problem
> >
> > # h_e is an array of coefficients associated with the degrees of
> freedom in
> > the system obtained from the eigenvectors V
> > h_e = V[:,sorted_index[first_index]]
> >
> > # setup a function to evaluate - the solution of the magnetic field
> > field = Function(element, self.mesh, h_e)
> >
> > # arbitrary point to evaluate the field
> > points = N.array((0.5,0.25))
> > # array to store the field value at a point
> > H = N.array((0.0,0.0))
> >
> > print "pre_eval"
> >
> > #evaluate the function
> > field.eval(H, points)
> >
> > # Ouput the field at the point
> > print "H",
> > print H
> >
> >
> > When I run this code, I get the following output:
> >
> > pre_eval
> > [labby:25160] *** Process received signal ***
> > [labby:25160] Signal: Segmentation fault (11)
> > [labby:25160] Signal code: Address not mapped (1)
> > [labby:25160] Failing at address: 0xde8f58ec
> > [labby:25160] [ 0] [0xb7f81440]
> > [labby:25160] [ 1] /usr/bin/python(PyEval_EvalCodeEx+0x6e7)
> [0x80cb0d7]
> > [labby:25160] [ 2] /usr/bin/python(PyEval_EvalFrameEx+0x565e)
> [0x80c92de]
> > [labby:25160] [ 3] /usr/bin/python(PyEval_EvalCodeEx+0x6e7)
> [0x80cb0d7]
> > [labby:25160] [ 4] /usr/bin/python(PyEval_EvalFrameEx+0x565e)
> [0x80c92de]
> > [labby:25160] [ 5] /usr/bin/python(PyEval_EvalCodeEx+0x6e7)
> [0x80cb0d7]
> > [labby:25160] [ 6] /usr/bin/python(PyEval_EvalCode+0x57) [0x80cb227]
> > [labby:25160] [ 7] /usr/bin/python(PyRun_FileExFlags+0xf8)
> [0x80ea6d8]
> > [labby:25160] [ 8] /usr/bin/python(PyRun_SimpleFileExFlags+0x199)
> > [0x80ea979]
> > [labby:25160] [ 9] /usr/bin/python(Py_Main+0xa35) [0x8059335]
> > [labby:25160] [10] /usr/bin/python(main+0x22) [0x80587f2]
> > [labby:25160] [11] /lib/tls/i686/cmov/libc.so.6
> (__libc_start_main+0xe0)
> > [0xb7de8450]
> > [labby:25160] [12] /usr/bin/python [0x8058761]
> > [labby:25160] *** End of error message ***
> >
> > With a previous build I have gotten different output and was able to
> > actually display the contents of H, but it would shill crash once the
> method
> > where the code was being run was terminated. I logged a bug against
> it at
> > http://www.fenics.org/cgi-bin/bugzilla/show_bug.cgi?id=55
> >
> > As I said in my mail, I see myself using fenics in the long run and
> would
> > like to get involved as much as possible, but some pointers as to
> where to
> > start would be much apprciated.
>
> Have you implemented the C++ version as well, such that you know that
> what you
> do is correct?
>
>
> Clearly, we are not catching what's going wrong here (missing check of
> some
> sort).
>
> My guess is that the error occurs in site-packages/dolfin/function.py,
> when
> the Function is initialized. Can you boil down the example to a shorter
> script
> and post it to this list? I'll have a look at it.
>
> Ola
>
>
> Hi
>
> I have not implemented a C++ version, but I will look into it - I am a
> little more comfortable with python at the moment, but I suppose I need to
> jump into the deep end :)
>
> I will write a simple example script and then attach that.
>
> Evan
>
>
>
> # A script to test the evaluation of a vector function from pydoflin
>
> __author__ = "Evan Lezar (evanlezar@xxxxxxxxx)"
> __date__ = "2008-08-21"
> __copyright__ = "Copyright (C) 2008 Evan Lezar"
> __license__ = "GNU LGPL Version 2.1"
>
> from dolfin import *
> import numpy as N
>
> def test_zero_function(order = 1):
> mesh = UnitSquare(2,2)
> element = FiniteElement("Nedelec", "triangle", 0)
>
> v = TestFunction(element)
> u = TrialFunction(element)
>
> L = dot(v, u)
> # assemble T to get the size of the system
> T = assemble(L*dx, mesh)
>
> num_dofs = T.size(1)
>
>
> # set the coeficients of all the dofs to 0
> h_e = N.zeros((1,num_dofs), dtype=N.float64)
>
> # initialize the function
> field = Function(element, mesh, h_e)
>
> point = N.array((0.0,0.0))
> H = N.array((0.0,0.0))
>
> print "About to evaluate the function"
>
> field.eval(H, point)
>
> print "H",
> print H
>
> assert ( H[0] == 0 and H[1] == 0 )
>
> if __name__ == "__main__":
> order = 0
> test_zero_function(order)
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
Attachment:
signature.asc
Description: Digital signature
Follow ups
References