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