dolfin team mailing list archive

["Evan Lezar" <evanlezar@xxxxxxxxx>]

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)