On Tue, Jul 01, 2008 at 04:59:54PM +0200, Andy Ray Terrel wrote:
Okay so I am trying to get the Scott Vogelius Stokes algorithm working
again, but I've hit a bit of a wall with the copying 4th ordered
functions. In the algorithm there is a term that gets updated every
iteration (w = w + r * u). So I could do this with a form as Anders
has suggested before, but it seems silly to use the 10K lines to do
what essentially is a vector axpy. Since some people at simula are
working on the function classes, I pose the question: Should we be
able to do this with functions? For example something like:
...
1 Function w;
2 w = Function(mesh, 0.0);
3 BilinearForm a(-1e3);
4 for (...) {
5 LinearForm L(f,w);
6 LinearPDE pde(a,L,mesh,bcs);
7 pde.solve(u);
8 w += 1e3 * u;
...
9 }
...
I think this should be possible. There could be a simple check to see
whether or not the involved Functions are defined on the same mesh
with the same element and in that case perform the operation on the
underlying Vectors (and otherwise give an error message).
We're about to start reworking the DofMap and Function classes. We
recently came to some kind of conlusion on the linear algebra classes
so this is up next.
Okay so assignment doesn't work for Constant functions (simple fix I
put in my local repo). Okay operaters +, +=, * are not defined, so I
expand line 8 to:
w.vector().axpy(1e3,u.vector());
Okay then line 8 won't work because w was previously a constant
function. I tried adding something like:
if (iter == 0) w = u;
but that didn't work because it is 4th order and the ElementLibrary
doesn't define it.
I have tried monkeying around with creating a better assignment
operator and copy constructor but then I run into other problems. For
example who owns the finite_element and DOF, in my current code they
are getting destroyed with the form and that kills later preprocessing
I do.
I have some hacks that basically work, but I wanted to poll the list
to see what is being planned and put my two cents in about it.
What if you just start with a discrete function which is zero?
The following would work in Python:
w = Function(element, mesh, Vector())
Or start by projecting whatever your initial data is to a discrete
function:
w = project(u0, element)
where u0 is some user-defined function (overloading eval) and
returning for example sin(x[0]).
