dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

On Wed, Dec 17, 2008 at 11:07 AM, Anders Logg <logg@xxxxxxxxx> wrote:
> On Wed, Dec 17, 2008 at 09:18:42AM +0100, Martin Sandve Alnæs wrote:
>> On Wed, Dec 17, 2008 at 8:44 AM, Anders Logg <logg@xxxxxxxxx> wrote:
>> > On Wed, Dec 17, 2008 at 12:15:38AM +0000, Garth N. Wells wrote:
>> >> I'm running into some problems with the new Function design when trying
>> >> to update a solver. I'll try to sketch the issue as simply as possible.
>> >>
>> >>      Function v;
>> >>      Function u, p;
>> >>
>> >>      // First form for U = (u, p)
>> >>      FirstBilinearForm a_1(V1, V1);
>> >>      FirstLinearForm   L_1(V1);
>> >>      LinearPDE pde_first(a_1, L_1);
>> >>
>> >>      // Second form which depends on u
>> >>      SecondBilinearForm a_2(V2, V2);
>> >>      SecondLinearForm   L_2(V2, u);
>> >>      LinearPDE pde_second(a_2, L_2);
>> >>
>> >>      Function U;
>> >>      for(uint i = 0; i < 10; ++i)
>> >>      {
>> >>        pde_first.solve(U);
>> >>
>> >>        // This step breaks down because FunctionSpaces don't match
>> >>        u = U[0];
>> >>        p = U[1];
>> >>
>> >>        pde_second.solve(v);
>> >>      }
>> >>
>> >> The problem is in assigning Functions since we now check the FunctionSpace.
>> >
>> > Which check is it that fails?
>> >
>> >> To get around this, I tried
>> >>
>> >>      FunctionSpace V(mesh);
>> >>      Function U(V);
>> >>      Function u = U[0];
>> >>      Function p = U[1];
>> >>
>> >> which throws an exception because U does not have a vector.
>> >
>> > What if we remove the check in operator[] and then in the constructor
>> > that gets a SubFunction we rely on v.v.vector() which (after Martin's
>> > fix from yesterday) will always return a vector with the dofs
>> > (possibly interpolated if there were no dofs before).
>>
>> Depends...
>>
>> Function u;
>> Vector & v = u.vector(); // error, no function space, ok
>>
>> MyFunction u(V);
>> Vector & v = u.vector(); // calls interpolate, using MyFunction::eval
>>
>> Function u(V);
>> Vector & v = u.vector(); // What happens here?
>
> I looked more closely at this, and as far as I can tell, your two last
> examples will result in the same thing: a zero vector will be created:
>
> if (!_vector)
> {
>  init();
>  interpolate(*_vector, *_function_space);
> }
>
> The call to init() first creates a zero vector. Then the function will
> be interpolated and the coefficients placed in the vector. This
> ultimately results in a call to
>
>  Function::interpolate(double* coefficients,
>                        const FunctionSpace& V,
>                        const ufc::cell& ufc_cell,
>                        int local_facet) const
>
> This function check if there is a vector of dofs (which there now is)
> and in that case picks the dofs from the vector (without calling
> evaluate_dof or eval).

Ok, so my "fix" didn't work at all. I'll remove it then.

I've spent some time debugging my sfc element code only to find that I was
using the wrong function spaces in a way that dolfin should have seen.
Seems like these parts of the new code needs some more work.
I'll try to make a short example.

>> Also, I had another function space error yesterday.
>> Trying to interpolate a user function into a discrete
>> function space fails because the user function doesn't
>> have a function space, but that shouldn't be necessary.
>> (btw, all the "interpolate" variants makes it difficult to
>> discuss and remember exact signatures etc...)
>
> Should we rename some of them?

I think so.

Martin