dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #11188
Re: Function design
On Wed, Dec 17, 2008 at 12:46:53PM +0100, Martin Sandve Alnæs wrote:
> 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.
I think we have more checks now than we used to have, but any
additional checks are welcome. Just add them whenever you encounter
a problem that we should have caught.
> >> 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.
ok, whenever anyone has a good suggestion, let's change it.
--
Anders
Attachment:
signature.asc
Description: Digital signature
References