Subject: [FFC-dev] mixed element issue?

From: "Robert Kirby" <[email-masked]>
Date: Wed, 16 Apr 2008 13:57:53 -0500

Hi all, I'm trying to do a survey of all the different ways of enforcing
incompressibility in the Navier-Stokes equations, and FFC/FIAT of course has
the most advanced element capabilities for doing this...

Unfortunately, I have hit a snag (either in my ignorance or through a bug)
in my grand scheme, starting with the Oseen equations (fix the convective
velocity for fixed points, don't do the full Newton linearization).

It seems that if you put a vector element into a mixed element and then take
it back out, you can no longer index it.

Putting the following in a .form (everything that doesn't involve pressure
spaces) works beatifully
#start .form file here
vel_element = VectorElement("Lagrange","triangle",2)
u = TrialFunction( vel_element )
v = TestFunction( vel_element )
u0 = Function( vel_element )

a = D(v[i],j)*D(u[i],j)*dx + v[i]*u0[j]*D(u[i],j)*dx
#end .form file

On the other hand, if I do:
# start new .form file here
vel_element = VectorElement("Lagrange","triangle",2)
pre_element = FiniteElement("Lagrange","triangle",1)
element = vel_element + pre_element
(u,p) = TrialFunctions( element )
(v,q) = TestFunctions( element )
u0 = Function( vel_element )

a = D(v[i],j)*D(u[i],j)*dx + v[i]*u0[j]*D(u[i],j)*dx
# end new .form file

I get an error:

Preprocessing form file: Oseen.form --> Oseen.py

*** object cannot be interpreted as an index
*** To get more information about this error, rerun ffc with the option -d1.

I've attached the output using -d1.

Any corrections either of my stupidity or of FFC would be appreciated :)

Rob

---

From: "Johan Hoffman" <[email-masked]>
Date: Wed, 16 Apr 2008 21:24:02 +0200 (MEST)

> Any corrections either of my stupidity or of FFC would be appreciated :)

Hi Rob!

I don't think you have to question your sanity, we have experienced
similar problems when dealing with the compressible Euler equations: with
scalar density and energy together with a momentum vector. We have been
forced to use built in operators instead, which is non optimal since the
definitions of the action of some operators is a bit unclear.

I do not know what the underlying problem is, but we have brought this up
before and it seems to be non-trivial? Or maybe (hopefully) I am wrong?

This is also something I want to stress for the coming UFL; that to me in
terms of operators, the most important is to have a fully reliable system
for indexing. Then to add a nice collection of well defined operators on
top is great, but secondary (to me).

I wasn't sure I was supposed to edit the UFL page directly with this
wishes? What do you say Anders?

/Johan

>
> Rob
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: "Murtazo Nazarov" <[email-masked]>
Date: Wed, 16 Apr 2008 21:38:31 +0200 (CEST)

> wishes? What do you say Anders?

Having index notation would be the best in general. Because, you know that
you have full control of the system you are solving. A similar issue I
posted a month ago. It seems there is no index notation defined in mixed
elements. As Johan pointed it out we use mix elements in compressible
solver. Now, instead of short and clear definition using index notation I
split out the equations to the several one: like for momentum I have three
separate equations which makes difficulties if I need to change some part
of model.

In the future, since the aim of Unicorn is to have a general solver for
continuous mechanics, having index notation would be preferable. I hope it
will be included in coming releases.

/murtazo

---

From: Marie Rognes <[email-masked]>
Date: Wed, 16 Apr 2008 21:41:38 +0200

Robert Kirby wrote:
>
> It seems that if you put a vector element into a mixed element and 
> then take it back out, you can no longer index it.
>

What happens is that TestFunction (and TrialFunction) returns an object 
(BasisFunction) that has a well-defined operator for taking components 
with arbitrary indices.

TestFunctions (and TrialFunctions), however, returns "all the 
components" structured as a tuple of lists. Taking "an arbitrary index" 
of a python list doesn't work that well.

Will think some more about whether this is 
trivial/easy/difficult/impossible to change.

PS: If you want something like this to work *right* now try iterating 
over i and j manually ;) Simplify will contract it back for you.

-- 
Marie E. Rognes
Ph.D Fellow, 
Centre of Mathematics for Applications, 
University of Oslo
http://folk.uio.no/meg

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Apr 2008 22:14:42 +0200

> over i and j manually ;) Simplify will contract it back for you.

Yes. In some more detail, when you do

FFC first creates a vector-valued element (mixed element) with three
components, the first two for u and the last for p. Then when you do

  (u, p) = TrialFunctions(element)

FFC picks out the three components of the vector-valued function
w = (u, p), indexes w and returns ([w[0], w[1]], w[2]).

This breaks the index notation since Python lists obviously can't
be index by an FFC Index.

I'd recommend doing something like

  a = sum([D(v[i],j)*D(u[i],j)*dx + v[i]*u0[j]*D(u[i],j)*dx for i in range(2) for j in range(2)])

Does that work for you?

I'm pretty sure we won't invest much time in improving the FFC form
language before switching to UFL. We (mostly Martin) have done some
work on it. You can follow the progress here:

  http://www.fenics.org/dev/ufl

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Apr 2008 22:18:50 +0200

> top is great, but secondary (to me).

Indexing is one of many priorities for UFL.

> wishes? What do you say Anders?

Yes, just edit the page with suggestions. That's what it's for.

-- 
Anders

---

From: "Robert Kirby" <[email-masked]>
Date: Wed, 16 Apr 2008 15:20:04 -0500

Thanks for all the feedback.  I will try the sum + list comprehension.

Rob

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 16 Apr 2008 23:49:39 +0200

> Yes, just edit the page with suggestions. That's what it's for.

All suggestions are welcome!

Perhaps the best way to add suggestions is to write an example form
with wanted syntax and explain with verbose comments what you mean by
it. Then I can use the code as a unit test as well.

Note that the page may be slightly out of sync with what's
implemented, I'll get around to cleaning this up in a few weeks (other
pressing matters right now), and can merge your additions in then.

--
Martin

---

From: "Robert Kirby" <[email-masked]>
Date: Tue, 22 Apr 2008 21:49:18 -0500

A question as I get back up the learning curve:
We now have H(div) capabilities in FFC and also DG capabilities thanks to
the hard work of many people.

Has anyone been so perverse as to try to put these together?

I ask because Bangerth & Kanschat have a neat technique for Navier-Stokes
wherein
they do the velocity with Raviart-Thomas elements.  These are in H(div) but
not H^1, so the
put DG-type penalty terms in the tangential direction.

Any thoughts on whether this will work?  If it is known not to work, I would
like to request the feature of DG + Piola-mapped elements...  If it is
unknown, I will be happy to be a guinea pug.

Thanks,
Rob

On Wed, Apr 16, 2008 at 4:49 PM, Martin Sandve Alnæs <[email-masked]>
wrote:

> >  > > incompressibility in the Navier-Stokes equations, and FFC/FIAT of
> course
> >  > > has
> then
> >  > > take
> >  > > *** To get more information about this error, rerun ffc with the
> option
> >  > > -d1.
> >  > >

---

From: "Marie E. Rognes" <[email-masked]>
Date: Wed, 23 Apr 2008 11:24:36 +0200

Robert Kirby wrote:
> We now have H(div) capabilities in FFC and also DG capabilities thanks 
> to the hard work of many people.
>

Yes, I think so. Take a look at the thread

    "DG operators for BDM" at the ffc-dev archives from September 30th 
2007.

Was it something like this you were thinking of?

--
Marie