Anders Logg wrote:
> 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]).
>
I do stuff like line 8 all the time, but usually much more complex.
You can initialise w as a discrete function, then to perform line 8
access the vector associated with each Function, w.vector(), u.vector().
You can then get the vector values from u in an array, loop of the array
elements, and then put the result back into w.
You can do this efficiently without using set() and get() if you use
uBlasVector.
These types of operations will definitely be a consideration in the
re-design on Function. I use far more projections for trivial operations
than I would like to.
Garth
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> DOLFIN-dev mailing list
> [email protected]
> http://www.fenics.org/mailman/listinfo/dolfin-dev
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
