Subject: [Ffc] Status

From: Anders Logg <[email-masked]>
Date: Wed, 6 Jan 2010 11:13:57 +0100

We've come pretty far on the rewrite of FFC and only a few functions
remain. As far as I can see, it remains to implement code generation
for the following functions:

Kristian:

  code["evaluate_basis_all"] = not_implemented
  code["evaluate_basis_derivatives"] = not_implemented
  code["evaluate_basis_derivatives_all"] = not_implemented

Marie:

  code["tabulate_entity_dofs"] = not_implemented # Marie doesn't know what this function should do
  code["tabulate_coordinates"] = not_implemented # Marie doesn't believe in this function

Marie or Kristian:

  code["interpolate_vertex_values"] = not_implemented

In addition, we need to get the code generation for the integral
classes working. It's sort of working for the tensor representation
but not at all (?) for the quadrature representation. I'll be working
on getting the tensor representation completely in place.

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 6 Jan 2010 11:24:04 +0100 (CET)

>  code["evaluate_basis_all"] = not_implemented

This should be easy to implement.
Would it be an idea to let evaluate_basis() call evaluate_basis_all() and then pick the correct value?
This will remove a lot of redundant code at the cost of some overhead introduced by the function call.
However, I don't see why one would call evaluate_basis() if evaluate_basis_all() works, I imagine that the assembler will call evaluate_basis_all() in the future for efficiency?

> on getting the tensor representation completely in place.

I don't think the quadrature code for integrals work at all, but I haven't checked. If you set up the infrastructure for combining code from different integrals, I'll take care of the quadrature part.

Kristian

> --
> Anders
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAktEYmQACgkQTuwUCDsYZdHCvgCcDbc6GpfqnJGi5iUJqxj+1uRq
> WTEAnivVMK4Kw1tW7jxdqH+zO8oAJH/j
> =HEMl
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp
>
>

---

From: Anders Logg <[email-masked]>
Date: Wed, 6 Jan 2010 11:57:25 +0100

> However, I don't see why one would call evaluate_basis() if evaluate_basis_all() works, I imagine that the assembler will call evaluate_basis_all() in the future for efficiency?

I don't think the assembler calls evaluate_basis_all. It's only used
in special cases for people that write their own assemblers (Nitsche
on overlapping meshes, crack propagation, partition of unity etc).

But having evaluate_basis call evaluate_basis_all sounds like it would
result in very much overhead. It would evaluate all basis functions
even if just one is needed. Wouldn't it be more natural to do the opposite?

> I don't think the quadrature code for integrals work at all, but I haven't checked. If you set up the infrastructure for combining code from different integrals, I'll take care of the quadrature pa…

ok. I will look at it. But you can also experiment with just changing
one line in codegeneration.py. Change from

  code = generate_tensor_integrals(ir, options)

to

  code = generate_quadrature_integrals(ir, options)

Note that in the new approach, all code needs to be generated as a
string, not a list of declarations but maybe that was only used by the
tensor representation before anyway.

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 6 Jan 2010 14:31:44 +0100 (CET)

> on overlapping meshes, crack propagation, partition of unity etc).

I think I tracked a call from the assembler to evaluate_basis, but that was in a very old DOLFIN I guess. Now evaluate_basis() is just called from one of the Function::eval() functions.

> even if just one is needed. Wouldn't it be more natural to do the opposite?

But if users will use evaluate_basis_all() (if available) rather than evaluate_basis() it is better to call evaluate_basis from evaluate_basis_all than the other way around.
The fastest will be to just include all the code in both functions.
Getting the value of all basis functions is not that demanding once the values of the expansion polynomial basis are computed (which is the same for all basis functions). Then it is just a matrix vector multiply of the coefficients (tabulated) and the basis values (computed at run-time).

> tensor representation before anyway.

OK

Kristian

> iEYEARECAAYFAktEbJQACgkQTuwUCDsYZdE6JwCePzSA6iEP1Pw9ODwYyN+wc+e7
> kukAoIem9SAzancyFMvKlJn//PWVWTdi
> =vZ7I

---

From: Marie Rognes <[email-masked]>
Date: Wed, 06 Jan 2010 16:25:30 +0100

>   code["tabulate_entity_dofs"] = not_implemented # Marie doesn't know what this function should do

Evidently, my comment was unclear ;) 

Should it look like this for CG_1?

  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
  virtual void tabulate_entity_dofs(unsigned int* dofs,
                                    unsigned int d, unsigned int i) const
  {
       dofs[0] = i  
  }

>   code["tabulate_coordinates"] = not_implemented # Marie doesn't believe in this function

What should this function do for the integral moment dofs, where there 
are more than one point per dof?

>   code["interpolate_vertex_values"] = not_implemented

Can fix.

--
Marie

> ------------------------------------------------------------------------
>

---

From: Marie Rognes <[email-masked]>
Date: Wed, 06 Jan 2010 16:29:09 +0100

Kristian Oelgaard wrote:
>>> Would it be an idea to let evaluate_basis() call 
>>> evaluate_basis_all() and then pick the correct value?
>>> This will remove a lot of redundant code at the cost of some 
>>> overhead introduced by the function call.
>>> However, I don't see why one would call evaluate_basis() if 
>>> evaluate_basis_all() works, I imagine that the assembler will call 
>>> evaluate_basis_all() in the future for efficiency?
>>
> I think I tracked a call from the assembler to evaluate_basis, but 
> that was in a very old DOLFIN I guess. Now evaluate_basis() is just 
> called from one of the Function::eval() functions.
>
> But if users will use evaluate_basis_all() (if available) rather than 
> evaluate_basis() it is better to call evaluate_basis from 
> evaluate_basis_all than the other way around.
> Getting the value of all basis functions is not that demanding once 
> the values of the expansion polynomial basis are computed (which is 
> the same for all basis functions). Then it is just a matrix vector 
> multiply of the coefficients (tabulated) and the basis values 
> (computed at run-time).
>

The same applies for evaluate_dof and evaluate_dofs. I thought that I'd 
implement these two differently,
so that evaluate_dof is as fast as possible for a single dof and 
evaluate_dofs is as fast as possible
for all the dofs.

--
Marie

>>> I don't think the quadrature code for integrals work at all, but I 
>>> haven't checked. If you set up the infrastructure for combining code 
>>> from different integrals, I'll take care of the quadrature part.
>>

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 06 Jan 2010 15:32:49 +0000

Marie Rognes wrote:
>       dofs[0] = i   }

An error?

It is useful to know the dof location for some problems. A simple
example is for visualising higher-order Lagrange functions, which is
something that some libraries do support.

Garth

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Jan 2010 14:43:58 +0100

>                                    unsigned int d, unsigned int i) const

Yes, this function should return an array of the dofs (dof indices) for
the dofs associated with entity number i of dimension d.

But perhaps we should check the d and the i as well to make sure that
d == 0 and i <= d + 1?

Here's another example of what should happen for P3 in 2D.

  switch (d)
  {
  case 0:
    if (i > 2)
      throw std::runtime_error("Sensible error message.");
    dofs[0] = i;
  case 1:
    if (i > 2)
      throw std::runtime_error("Sensible error message.");
    dofs[0] = 3 + 2*i;
    dofs[1] = 3 + 2*i + 1;
  case 2:
    if (i > 0)
      throw std::runtime_error("Sensible error message.");
    dofs[0] = 10;
  default:
    throw std::runtime_error("Sensible error message.");
  }

The corresponding function num_entity_dofs should return the size of
the dofs array for each dimension so it should return 0, 2 and 1 in
the above case for d = 0, 1, 2.

> What should this function do for the integral moment dofs, where
> there are more than one point per dof?

I guess we need to throw an exception and think about possibly extending the
UFC interface with this extra information.

--
Anders

---

From: Marie Rognes <[email-masked]>
Date: Thu, 07 Jan 2010 15:04:35 +0100

>     throw std::runtime_error("Sensible error message.");

Ok!

> the above case for d = 0, 1, 2.

Fixed already.

> UFC interface with this extra information.

Will throw.

--
Marie

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Jan 2010 15:16:59 +0100

> >the above case for d = 0, 1, 2.

Nice! :-)

--
Anders

---

From: Marie Rognes <[email-masked]>
Date: Tue, 12 Jan 2010 19:00:54 +0100

> on getting the tensor representation completely in place.

How are things going?

I'm have a bug or two to fix in evaluate_dof for integral moments and 
will then fix evaluate_dofs.
Need to sort out FIAT normal directions before correct code is generated 
for RT/BDM though.
Otherwise good.

Can we move the prefetching of the various format["..."] into the global 
scope of codegeneration.py
or is that evil?

--
Marie

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 12 Jan 2010 19:48:41 +0100 (CET)

2010/1/12 Marie Rognes <[email-masked]>:
> I'm have a bug or two to fix in evaluate_dof for integral moments and will
> then fix evaluate_dofs.
> Need to sort out FIAT normal directions before correct code is generated for
> RT/BDM though.
> Otherwise good.
>

We need to do something different with the way we handle the formatting. At the moment we just import format from cpp in codegeneration.py, but that is clearly not an option if one wants to generate different code.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Tue, 12 Jan 2010 19:48:44 +0100

> >on getting the tensor representation completely in place.

I haven't had time to program much the last couple of days so
tabulate_tensor is in the same shape as before = somewhat working but
needs testing + bug fixes + formatting fixes.

I'll get back to it tomorrow.

> I'm have a bug or two to fix in evaluate_dof for integral moments
> and will then fix evaluate_dofs.
> Need to sort out FIAT normal directions before correct code is
> generated for RT/BDM though.
> Otherwise good.

Splendid.

> Can we move the prefetching of the various format["..."] into the
> global scope of codegeneration.py
> or is that evil?

Isn't it quite useful to have it in each function so we can see
clearly which formatting rules are used by each function?

--
Anders

---

From: Marie Rognes <[email-masked]>
Date: Tue, 12 Jan 2010 20:05:27 +0100

Kristian Oelgaard wrote:
>>>  code["tabulate_entity_dofs"] = not_implemented # Marie doesn't know 
>>> what
>>> this function should do
>>>  code["tabulate_coordinates"] = not_implemented # Marie doesn't 
>>> believe in
>>> this function
>>>
>>> Marie or Kristian:
>>>
> We need to do something different with the way we handle the 
> formatting. At the moment we just import format from cpp in 
> codegeneration.py, but that is clearly not an option if one wants to 
> generate different code.
>

Do we have any examples of different code that we want to generate?

--
Marie

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 12 Jan 2010 20:32:13 +0100 (CET)

2010/1/12 Marie Rognes <[email-masked]>:
>>>>  code["tabulate_coordinates"] = not_implemented # Marie doesn't believe
>>>> in
>>>> this function
>>>>
>>>> Marie or Kristian:
>>>>
>>> for
>> We need to do something different with the way we handle the formatting.
>> At the moment we just import format from cpp in codegeneration.py, but that
>> is clearly not an option if one wants to generate different code.
>>
>

Currently not.
Another issue is cpp.py importing FFC_OPTIONS which will lead to some difficulties I think w.r.t. format["epsilon"].

Kristian

---

From: Anders Logg <[email-masked]>
Date: Tue, 12 Jan 2010 20:49:15 +0100

> >Do we have any examples of different code that we want to generate?

Initially, I had in mind that one would like to generate different
code as an option (maybe Fortran) but it hasn't happened yet.

So for now it's just a convenient way to make it easy to edit how the
code is formatted in one single place without needing to search for +,
* etc throughout the code.

> Another issue is cpp.py importing FFC_OPTIONS which will lead to some difficulties I think w.r.t. format["epsilon"].

Yes, it's an issue. It's something we need to sort out. I'll get to it
in time.

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 12 Jan 2010 21:00:11 +0100 (CET)

> * etc throughout the code.

Yes, I totally agree. And when/if we decide to generate Fortran I will be very surprised if simply feeding the new format will work out of the box anyway.

Kristian

> iEYEARECAAYFAktM0jsACgkQTuwUCDsYZdGJzwCfUsPOQ5MgiB1+GE9y6xUsEzYy
> JtMAn14KSOisUYkuMq/2wlfzyz/4YpFu
> =IwLF

---

From: Anders Logg <[email-masked]>
Date: Tue, 12 Jan 2010 21:30:31 +0100

On Tue, Jan 12, 2010 at 08:09:25PM +0100, Marie Rognes wrote:
> Ok.
> It is just that some formats, e.g. assign, component and return, are
> used by virtually
> all of the non-trivial codegeneration functions.

ok, maybe it would not be so bad to extract 3-4 commonly used
formatting rules at the top.

--
Anders