Subject: [UFC-dev] added higher mesh variable

From: Shawn Walker <[email-masked]>
Date: Wed, 22 Apr 2009 17:28:30 -0400 (EDT)

Here is the changeset that adds a `higher_order_coordinates' variable for 
storing higher order mesh data.  This is a very minor change so please 
push this.

A changeset for DOLFIN is coming immediately after this.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Sat, 25 Apr 2009 00:14:31 +0200

> A changeset for DOLFIN is coming immediately after this.

I'm not sure what to do about this. It's problematic to add
experimental work to UFC since it must be stable. In particular, any
small change to ufc.h means that all forms must be recompiled
everywhere for everyone.

So before we make a change to UFC, we need to know exactly what we
need. Which also means I can't import your DOLFIN patch since it
depends on the UFC patch.

I see you've added

    double** higher_order_coordinates;

to ufc::cell. This is analogous to what is now implemented in
MeshGeometry and the mesh XML format so I think it's good.

The question is what other information we need. As it works now (for
the standard ufc::cell), UFC code generated by a form compiler knows
what to expect from for a ufc::cell argument. If higher order mappings
should work the same way, then the generated code and thus the form
compilers need to know which mapping should be used and also the
length of higher_order_coordinates. Is this what you were thinking?

Before we do much more about it, more people need to weigh in on it as
it affects DOLFIN, UFC, SyFi and FFC.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Fri, 24 Apr 2009 23:04:02 -0400 (EDT)

> small change to ufc.h means that all forms must be recompiled

I see...

> length of higher_order_coordinates. Is this what you were thinking?

Yes, that is what I was thinking.  There needs to be some variable that 
stores the number of DoF's used to represent the geometry.  And I also 
think there should be a string with the "name" of the mapping.  That way 
there can be an error check, because presumably FFC will `assume' a type 
of mapping (based on user input).

For the stuff I did in MeshGeometry, I plan on going back and making a 
single struct variable that has all the pertinent higher order mesh 
variables together.  Depending on what is needed in ufc::cell, the same 
thing may be in order.

> it affects DOLFIN, UFC, SyFi and FFC.

ok.  I say you need:

1.  higher order vertex data that is already in the PROPER
     order/numbering. 
2.  the number of higher order vertices.
3.  a string identifying the type of map i.e. "P2_Lagrange".
4.  AND a boolean variable stating whether the cell is affine or not (i.e. 
straight or curved).  This is useful because then the tabulate_tensor 
routines can have a switch between the standard code, and the probably 
more complicated curved element code.  Remember, in most cases, only the 
boundary triangles need to be curved; the interior elements can be 
straight (IN MOST CASES).

And I would probably put this into a single struct.  I don't think a whole 
separate class is necessary.  But I am open to suggestions.

I hope people comment on this.  This effects many things besides the 
matrix assembly.  It affects interpolation, computing the facet normal, 
etc...

- Shawn

---

From: Kent Andre <[email-masked]>
Date: Mon, 27 Apr 2009 11:26:13 +0200

> it affects DOLFIN, UFC, SyFi and FFC.

But is there any other way around this. It would be nice with higher
order meshes and UFC should not stop this. 

An alternative to changing the cell class would be to make a subclass
of cell. Would this work ? 

Kent

---

From: Anders Logg <[email-masked]>
Date: Mon, 27 Apr 2009 14:52:25 +0200

> of cell. Would this work ? 

How about just using the current ufc::cell data structure as it is but
let coordinates hold all the coordinates?

This could also be the final solution. Then everything that's needed
is an extra argument to tabulate_tensor that tells the generated code
whether the cell is affinely mapped or not. The flag could simply be
an integer: 1 means affine, 2 means quadratic etc.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Mon, 27 Apr 2009 09:21:58 -0400 (EDT)

> an integer: 1 means affine, 2 means quadratic etc.

But you still need to modify the ufc::cell code, I think.  There is also 
an implicit assumption that the higher order coordinates should contain 
the standard mesh vertex coordinates.  Of course, this is true for 
most practical cases.  But for more fancy mappings, maybe this is not the 
case.

Also, in the ufc::cell code, you currently read in the cell coordinates 
using info in MeshTopology.  However, the higher order coordinate info 
resides in MeshGeometry (which is where it belongs).  So you would still 
need to modify ufc.h.   Remember, there is higher order cell data that is 
contained in MeshGeometry.

Is it really that hard to change ufc.h?  Other things have to be 
recompiled, but isn't that automatic?

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Mon, 27 Apr 2009 15:27:43 +0200

On Mon, Apr 27, 2009 at 09:21:58AM -0400, Shawn Walker wrote:
> the standard mesh vertex coordinates.  Of course, this is true for most 
> practical cases.  But for more fancy mappings, maybe this is not the  
> case.

It seems to me that a reasonable assumption would be to limit the
cases to P1, P2, P3, etc, that is, mappings that can be written down
using standard Lagrange bases so then the vertices will always be
included. They would also be first in the list meaning that the code
would actually work (but might not give accurate results) even if it
were generated for affine mappings.

> contained in MeshGeometry.

Where is MeshTopology used for this? I looked in UFCCell.h which is
where the coordinates are copied to ufc::cell and there MeshGeometry
is used.

> recompiled, but isn't that automatic?

Yes, it's easy to change, but a main point with UFC is that we
shouldn't change it.

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 27 Apr 2009 14:32:20 +0100

> Yes, it's easy to change, but a main point with UFC is that we

UFC will need to be extended as time goes on, but it is hard to know 
from the outset how it should be done. What about using some IFDEF's or 
non-pure virtual functions in the development version to allow 
experimentation? These can then either be removed or added to UFC at 
release time.

Garth

> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> UFC-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufc-dev

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 27 Apr 2009 15:43:28 +0200

> experimentation? These can then either be removed or added to UFC at

Or subclasses with non-pure virtual functions:

class experimental_cell_integral: public ufc::cell_integral
{
  void foo() const { throw ...("Experimental feature not implemented."); }
};

or

namespace eufc {
  class cell_integral: public ufc::cell_integral
  {
    void foo() const { throw ...("Experimental feature not implemented."); }
  };
}

We can define these in "experimental_ufc.h" or "eufc.h" to keep the
official header file constant.

Then the DOLFIN code that uses experimental features must be clearly marked:

  ufc::cell_integral *itg = form.create_cell_integral(0)
  eufc::cell_integral *eitg = dynamic_cast<eufc::cell_integral>(itg);

and can then use "if(eitg)" to select between experimental and
non-experimental code.

Martin

---

From: Shawn Walker <[email-masked]>
Date: Mon, 27 Apr 2009 09:46:11 -0400 (EDT)

>> were generated for affine mappings.

I am not totally against limiting it to Lagrange only...

>> where the coordinates are copied to ufc::cell and there MeshGeometry

What I meant was that MeshTopology is used to get the local cell entities, 
and from that you get the local vertex coordinates (e.g. for a standard 
triangle with 3 vertices).  But this is done in the DOLFIN code, so it 
isn't a problem.

> UFC will need to be extended as time goes on, but it is hard to know from the 
> outset how it should be done. What about using some IFDEF's or non-pure 
> virtual functions in the development version to allow experimentation? These 
> can then either be removed or added to UFC at release time.

This might not be bad.  Eventually, UFC will need to handle these higher 
order meshes, and so it should contain some kind of data struct to hold 
that info.

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Mon, 27 Apr 2009 09:49:45 -0400 (EDT)

> and can then use "if(eitg)" to select between experimental and

This seems kind of complicated, but that is just me.  Originally, I was 
only asking for a small modification of ufc.h.  The real code would be 
contained in DOLFIN or the code generated by FFC.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Mon, 27 Apr 2009 17:39:27 +0200

> and can then use "if(eitg)" to select between experimental and

I like having a separate header file for experimental stuff. Maybe it
could be placed in DOLFIN? This will make it easier to experiment (no
need for synchronization) and keep UFC constant.

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 27 Apr 2009 16:42:59 +0100

> need for synchronization) and keep UFC constant.

Having it in DOLFIN would be nice for those of us without write 
permission to the ufc repository.

Garth

> DOLFIN-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/dolfin-dev

---

From: Anders Logg <[email-masked]>
Date: Mon, 27 Apr 2009 17:49:07 +0200

> permission to the ufc repository.

I'd be happy to give you write access if you want it (and Martin does
not object).

I thought the problem was I didn't have time to make the changes, but
maybe the problem is no one else has access.

Anyhow, any changes need to be discussed at some length. UFC is also
fairly well documented (in contrast to other projects we know of) so
any changes to the interface need to be documented in the manual.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Mon, 27 Apr 2009 12:16:14 -0400 (EDT)

> any changes to the interface need to be documented in the manual.

How about I just add the ``double** higher_order_coordinates;'' as a 
public variable to UFCCell.h?   And I do the same with whatever else is 
needed.  Then, once the interface is figured out, we can add this to ufc.h 
properly.  Sorry, I didn't do this to begin with.

- Shawn

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 27 Apr 2009 18:37:40 +0200

> I'd be happy to give you write access if you want it (and Martin does

No objection here.

> any changes to the interface need to be documented in the manual.

Agreed. (When moved from the experimental interface to the official, anyway).

Martin

---

From: Shawn Walker <[email-masked]>
Date: Mon, 27 Apr 2009 23:39:08 -0400 (EDT)

> and can then use "if(eitg)" to select between experimental and

In a similar vein, could we just have another file named eufc.h, and put 
an IFNDEF somewhere that would use eufc.h instead of ufc.h?  That way I 
could modify eufc.h all I want, and people don't have to use it.  But I'm 
not sure how to do this.  Ideally, this would be an option for scons like 
enableExpUFC=true.  Or is it only necessary to include eufc.h in files 
that FFC generates?  We just need something for testing.

Obviously, I cannot just modify UFCCell.h.  I tried that, but FFC cannot 
access variables declared in the sub-class (woops...  :(  ).

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Tue, 28 Apr 2009 07:33:09 +0200

> access variables declared in the sub-class (woops...  :(  ).

One option could be to create a file named ufc.h and put it in

  dolfin/fem/ufc.h

and change all #include <ufc.h> to #include "ufc.h".

Then the DOLFIN version of ufc.h will include either the installed
ufc.h or another file named ufce.h which is placed in

  dolfin/fem/ufce.h

That file contains data structures named the same way as in the
official ufc.h but with modifications (so the rest of the code won't
need to be changed much).

Then in the DOLFIN ufc.h control file we place an #ifdef for whether
to include the official ufc.h or the experimental one:

#ifdef UFC_EXPERIMENTAL
#include "ufce.h"
#else
#include <ufc.h>
#endif

You can set the flag by adding

  customCxxFlags="-DUFC_EXPERIMENTAL"

to scons.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Tue, 28 Apr 2009 09:06:28 -0400 (EDT)

> and change all #include <ufc.h> to #include "ufc.h".

Are there several places to change this?  All the demos?

---

From: Shawn Walker <[email-masked]>
Date: Tue, 28 Apr 2009 22:27:15 -0400 (EDT)

On Tue, 28 Apr 2009, Anders Logg wrote:

>> On Mon, 27 Apr 2009, Martin Sandve Aln?s wrote:
>>
>>>>>>>> On l?., 2009-04-25 at 00:14 +0200, Anders Logg wrote:

There is already a file called UFC.h.  This seems a little confusing to 
use ufc.h for this.  Can we rename this?

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Tue, 28 Apr 2009 22:50:11 -0400 (EDT)

>  customCxxFlags="-DUFC_EXPERIMENTAL"

Wouldn't it make more sense to have this exp_ufc.h be a part of ufc?  That 
way when you install ufc, there can be an `enableExpUFC' option.  In this 
case, scons will copy exp_ufc.h to the build directory and rename it to 
ufc.h.  The exp_ufc.h file will basically be identical to ufc.h, so it 
makes sense to keep it together.  This will also put in a convenient 
`buffer' for experimenting with ufc, and allow for easy moving over of 
additions to ufc.h.  Could someone please put this in?  Please?  :)

- Shawn

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Wed, 29 Apr 2009 08:34:43 +0200

On Wed, Apr 29, 2009 at 4:50 AM, Shawn Walker <[email-masked]> wrote:
>>>>>>> most
> ufc.h.  The exp_ufc.h file will basically be identical to ufc.h, so it makes
> sense to keep it together.  This will also put in a convenient `buffer' for
> experimenting with ufc, and allow for easy moving over of additions to
> ufc.h.  Could someone please put this in?  Please?  :)

I don't see the point with all this build system tweaking and file moving,
my solution was only a few lines C++ and allows you to use both "official"
ufc forms and "extended" versions in the same dolfin build, with
modifications only where things change. Which means a lot less
compilation time, and you can add stuff to the extended version
and implement support in dolfin and add stuff to ffc later.

But it's your time, suit yourself.

Martin

---

From: Anders Logg <[email-masked]>
Date: Wed, 29 Apr 2009 13:36:02 +0200

> > ufc.h.  Could someone please put this in?  Please?  :)

It's better to have it in DOLFIN. It makes it easier to experiment (no
need to change UFC). When features have stabilized, we move it to UFC.

I agree having both ufc.h and UFC.h is confusing. How about this
instead. In DOLFIN, we never include ufc.h, only UFC.h. And at the top
of that file, we put

#ifdef UFC_EXPERIMENTAL
#include "ufce.h"
#else
#include <ufc.h>
#endif

That should work.

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 29 Apr 2009 12:38:36 +0100

Anders Logg wrote:
> instead. In DOLFIN, we never include ufc.h, only UFC.h. 

We do include ufc.h in a number of places. This needs to be fixed so 
that a ufc(e).h file can be placed in DOLFIN.

Garth

And at the top

---

From: Kent Andre <[email-masked]>
Date: Wed, 29 Apr 2009 13:46:40 +0200

> that a ufc(e).h file can be placed in DOLFIN.

What about generated code ? 

Kent

---

From: Anders Logg <[email-masked]>
Date: Wed, 29 Apr 2009 13:49:12 +0200

> > > instead. In DOLFIN, we never include ufc.h, only UFC.h. 

What I meant was, we need to change this so we only include UFC.h in
places where we now include ufc.h. The inclusion of ufc.h will be
handled by UFC.h and it will pick the right one.

> What about generated code ? 

The experimental form compiler or the person hand-writing the
experimental code would need to handle that.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Wed, 29 Apr 2009 10:45:50 -0400 (EDT)

Here is a changeset that now has an option for using an experimental ufc.h 
that is internal to dolfin.  I made a new sub-directory: ufc_exp, and then 
put ufc.h in that.  That way the file name can be the same.  And we don't 
need to change other symbols/names that have `ufc' in them.  And there is 
another added file: dolfin_ufc.h, that does the selection.

Please check this and push it.

I also think it would be good to add an enable option to scons that would 
select the experimental ufc.h file, instead of manually modifying the 
scons script.  But I don't know how to do this.

- Shawn

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 29 Apr 2009 20:32:43 +0100

> experimental code would need to handle that.

There is come generated code in DOLFIN (in ale and mf) which complicates 
things.

Garth

---

From: Shawn Walker <[email-masked]>
Date: Wed, 29 Apr 2009 15:38:06 -0400 (EDT)

> There is come generated code in DOLFIN (in ale and mf) which complicates

Yeah, I saw that, but I wasn't sure what to do.  Should I change all that 
over too?

- Shawn

---

From: Ola Skavhaug <[email-masked]>
Date: Thu, 30 Apr 2009 09:50:39 +0200

Why not go for Martin's suggestion?

Inheritance is a well known concept in C++, and it avoids using the
evil macros that notoriously break compilation when people forget
their ifndefs. Extending ufc with subclasss in a namespace, say eufc,
dynamic casting in experimental assemblers when needed, and finally,
moving functionality up the inheritance chain when settled seems to me
a better solution than what is happening now?

Ola

On Wed, Apr 29, 2009 at 1:36 PM, Anders Logg <[email-masked]> wrote:
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkn4O6IACgkQTuwUCDsYZdG3bQCglt8WQtlFN3U/s73vkEr+rS0E
> lYIAn0tSAnDOgHgyEfnE48ZMFJecjfB8
> =p76f
> -----END PGP SIGNATURE-----
>

-- 
Ola Skavhaug

---

From: Shawn Walker <[email-masked]>
Date: Thu, 30 Apr 2009 08:54:23 -0400 (EDT)

> a better solution than what is happening now?

Why not do it?

Reasons:

1. I'm not that adept at C++ to know how to do what you are saying.  Nor 
do I forsee myself learning this in the near future. (i.e. like all this 
dynamic casting stuff).

2. This `eufc' is only for experimentation and I would assume it will not 
be a permanent piece of the code.

3. The original reason for doing this, was to basically just add some 
freaking variables!  Why do we have to make it so complicated? 
Furthermore, these additions (if made to the real ufc.h) would not have 
changed anything.  The current modules now would just not use those new 
data structs.  It would have been completely backwards compatible.  But I 
understand you want to keep ufc and its doc's up to date.

If someone wants to set this up, be my guest (see #1).  I would certainly 
appreciate it.  I don't know how to do it.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Thu, 30 Apr 2009 18:25:38 +0200

On Thu, Apr 30, 2009 at 08:54:23AM -0400, Shawn Walker wrote:
> freaking variables!

I know, and that's exactly the point. Before we had UFC, the interface
between FFC and DOLFIN changed all the time. Now we have two form
compilers and other projects using UFC or considering using it. So we
need to keep it stable and make sure we know what we want before
making changes. Any change to ufc.h will require changes to DOLFIN, FFC,
SyFi/SFC, and all code that has ever been generated for UFC 1.1.

> Why do we have to make it so complicated?  

I talked to Martin and Ola today and they convinced me we should use
inheritance, not macros (#ifdef).

The reason is that with inheritance, we can freely extend the classes
in ufc.h and the code (DOLFIN) will still work with the official UFC
interface. With #ifdef, one would not be able to assemble forms
compiled for the official UFC interface with DOLFIN when the
UFC_EXPERIMENTAL flag is set.

Here's a summary of what needs to be done:

1. Create a file named ufce.h and place it in dolfin/fem/ as part of
DOLFIN.

2. This file defines subclasses of the UFC classes that need to be
extended. It can look something like this:

#include <ufc.h>

namespace ufce
{
  
  class cell : public ufc::cell
  {
  public:

    /// Array of coordinates for higher order vertices of the cell
    double** higher_order_coordinates;

  };

}

3. UFCCell needs to inherit from ufce::cell, not ufc::cell. That means
it is still a ufc::cell since ufce::cell is a ufc::cell.

4. In places where we assume that we have a ufce::cell and not a
ufc::cell, we make a dynamic_cast and issue an error if the cast fails
(which indicates that the form that comes in is not a ufce form but a
ufc form). A dynamic_cast of a pointer can be done like this:

  Bar* bar = dynamic_cast<Bar*>(foo);

Here, foo is a pointer to a Foo and Bar is a subclass of Foo.
If the cast fails, a zero pointer will be returned indicating that foo
is not a Bar, only a Foo.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Thu, 30 Apr 2009 14:02:31 -0400 (EDT)

> SyFi/SFC, and all code that has ever been generated for UFC 1.1.

ok, ok.

> UFC_EXPERIMENTAL flag is set.

ok.

> 1. Create a file named ufce.h and place it in dolfin/fem/ as part of

This should be named ufcexp.h.  If someone is visually scanning the code, 
this can be missed.

> extended. It can look something like this:

namespace ufcexp

> it is still a ufc::cell since ufce::cell is a ufc::cell.

This seems reasonable...  again ufcexp::cell.  the other way can be easily 
overlooked.

> is not a Bar, only a Foo.

ok, so for example, say we are inside a tabulate_tensor routine.

when we read in `c.coordinates', everything is fine.  No problem.

Now, we want to read in the `higher_order_coordinates'.  So, we cannot do 
`c.higher_order_coordinates' because the routine only sees it as a 
ufc::cell (because we want to be backwards compatible).

So, we have to do this:

ufcexp::cell MYc = dynamic_cast<ufcexp::cell>(c);

Then we can read:

`MYc.higher_order_coordinates'.

Is this right?  If so, that doesn't sound too bad.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Thu, 30 Apr 2009 20:32:31 +0200

> This should be named ufcexp.h.  If someone is visually scanning the code, 

Good point.

> This seems reasonable...  again ufcexp::cell.  the other way can be 
> easily overlooked.

Sure.

> Is this right?  If so, that doesn't sound too bad.

Almost:

  const ufcexp::cell* cexp = dynamic_cast<const ufcexp::cell*>(c);
  dolfin_assert(cexp);

  cexp->higher_order_coordinates[...] 

It needs to be a pointer and it also needs to be const.

You might also get around with

  const ufcexp::cell& cexp = *dynamic_cast<const ufcexp::cell*>(c);
  cexp.higher_order_coordinates[...]

This should throw an exception (bad_cast) if it does not work.

But perhaps the best would be

  const ufcexp::cell* cexp = dynamic_cast<const ufcexp::cell*>(c);
  if (!cexp)
    error("Informative error message.");

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Thu, 30 Apr 2009 15:28:52 -0400 (EDT)

Ok, here is a changeset.  Only two files are modified: UFCCell.h and 
added a ufcexp.h.  It compiles and runs on my machine, so it should be 
fine.  Some of the code for the higher_order_coordinates is commented out 
in the UFCCell.h.  That is because I still have to re-submit the mesh 
geometry class modification.

Anyway, please check this and push it.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Thu, 30 Apr 2009 22:10:55 +0200

> Anyway, please check this and push it.

Looks good. Pushed.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Thu, 30 Apr 2009 16:53:47 -0400 (EDT)

Here is another changeset for reading in the higher order mesh vertices. 
It works on my end.  Please check and push.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Fri, 1 May 2009 09:53:06 +0200

Pushed.

-- 
Anders