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 <[email protected]> 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 <[email protected]> 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
signature.asc
Description: Digital signature
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
