Subject: Re: [FFC-dev] Tensors

From: <[email-masked]>
Date: Tue, 30 Aug 2005 05:48:24 -1100

This leads to maintaing two different kinds of tensors -- one 
for FFC and one in Numeric.  In my opinion (but perhaps not 
Anders' :)), FFC is still very experimental.  Your email 
leads one to speculate whether it would be a better practice 
to make everything in FFC a Numeric.array of objects.  Then, 
so long as we overload operators appropriately for the scalar 
objects, we're set for vectors and tensors. 
 
Rob 

---- Original message ---- 
>Date: Mon, 29 Aug 2005 15:56:56 +0200 
>From: Johan Jansson <[email-masked]>   
>Subject: Re: [FFC-dev] Tensors   
>To: "Garth N. Wells" <[email-masked]> 
>Cc: [email-masked] 
> 
>On Wed, Aug 24, 2005 at 01:32:19PM +0200, Garth N. Wells 
wrote: 
>> Hi Johan, 
>>  
>> For me, tensor notation would be really useful.  
>>  
> 
>Hi! 
> 
>I've experimented a bit with different approaches to tensor 
notation, 
>and I've found one that works pretty well, some feedback 
would be 
>nice. The approach I've gotten to work is to simply define 
matrices of 
>objects in the FFC algebra, I've done it using multiarrays 
from the 
>Numeric module, but any capable multiarray/matrix class 
could probably 
>be used, or we could write our own if we need special 
capabilities. 
> 
>What I would like to be able to do is write: 
> 
>epsilonu = 0.5 * (u[i].dx(i) * v[j].dx(j)) 
> 
>and then treat epsilonu as a tensor. But as far as I can 
see, epsilonu 
>is not a tensor in the FFC algebra. At least indexing is not 
yet 
>supported for such objects, so I cannot write epsilonu[1, 2] 
or even 
>epsilonu[1]. I'm not sure epsilonu is supposed to fit the 
tensor 
>concept here. 
> 
>What I did instead is to explicitly construct a tensor 
(matrix): 
> 
>from Numeric import * 
>from math import * 
> 
>def epsilon(u): 
>    eu = ffczeros(3, 3) 
> 
>    for ii in range(3): 
>        for jj in range(3): 
>            eu[ii, jj] = 0.5 * (u[ii].dx(jj) + u[jj].dx(ii)) 
> 
>    return eu 
> 
>also defining the function ffczeros(): 
> 
>global zero 
>zero = 0.0 * Constant() 
> 
>def ffczeros(m, n): 
>    A = array((zero)) 
>    A = resize(A, (m, n)) 
>    A *= 0.0 
> 
>    return A 
> 
>array() and resize() are Numeric.array() and 
Numeric.resize(). The 
>array in Numeric is not really made for arbitrary types, it 
seems to 
>only officially support builtin types, but this works well 
>enough. There is an array class in numarray which does 
support general 
>types, but it is not well developed. As I mentioned, we can 
write our 
>own matrix/tensor class if needed. 
> 
>I can now write: 
> 
>epsilonu = epsilon(u) 
> 
>to achieve what I want. 
> 
>It's then very easy to write forms using the same notation 
as in the 
>original PDEs. For example, here is the E() linear 
transformation from 
>linear elasticity: 
> 
>def E(e, lmbda, mu): 
>    Ee = 2.0 * multiply(mu, e) 
>    for ii in range(3): 
>        for jj in range(3): 
>            if ii == jj: 
>                Ee[ii, jj] += lmbda * trace(e) 
> 
>    return Ee 
> 
>I would have liked to write it like so: 
> 
>def E(e, lmbda, mu): 
>    Ee = 2.0 * multiply(mu, e) + lmbda * trace(e) * 
ffcidentity(3) 
>    return Ee 
> 
>but I have not yet succeeded in constructing an identity 
matrix of FFC 
>algebra objects. It would be nice to have a zero and 
identity element 
>in the FFC algebra, they should just be a Constant() with a 
0.0 or 1.0 
>value. I'm not really clear on the Python operators yet 
either, there 
>seems to be some ambiguity between when an operator is used 
from FFC 
>or Numeric.array. 
> 
>I have attached the elasticity form from the DOLFIN module 
written in 
>the new notation (Elasticity-alternative.form), I've 
verified that it 
>produces the same result, so it works. FFC seems to have 
some problems 
>with namespaces though, so you need to compile it like so: 
> 
>ffc Elasticity-alternative.form 
>python Elasticity-alternative.py 
> 
>> Related to the issue of taking slices and indexing, I had 
some 
>> exchanges with Anders a while ago about the possibility of 
being 
>> able to specify the range of the summation when using the 
index 
>> summation convention. This would be a big help when 
examining 
>> coupled problems. At the moment I'm not able (or don't 
know how) to 
>> use the index notation for these problems, so I'm writing 
out the 
>> bilinear form for explicitly. As an example, I'm looking 
at a 
>> coupled elasticity-diffusion equation at the moment. In 2D 
(3 dof 
>> per node), part of th bilinear form looks like 
> 
>>  
>> a = ( lmbda * (u[0].dx(0) + u[1].dx(1))*(v[0].dx(0) + 
v[1].dx(1))  + \ 
>>        mu * ( \ 
>>         (u[0].dx(0)+u[0].dx(0))*(v[0].dx(0)) + \ 
>>         (u[0].dx(1)+u[1].dx(0))*(v[0].dx(1)) + \ 
>>         (u[1].dx(0)+u[0].dx(1))*(v[1].dx(0)) + \ 
>>         (u[1].dx(1)+u[1].dx(1))*(v[1].dx(1)) ) )* dx + \ 
>>      -Beta*( u[2]*v[0].dx(0) + u[2]*v[1].dx(1) )* dx 
>>  
> 
>> For 3D, I need to add a bunch of terms. What I would like 
to able to 
>> do is use something like: 
> 
> 
>> a = ( lmbda * u[i:0,n].dx(i:0,n) * v[j:0,n].dx(j:0,n) + mu 
* ( 
>> u[i:0,n].dx(j:0,n) + u[j:0,n].dx(i:0,n) ) * 
v[i:0,n].dx(j:0,n) )* dx 
>> + -Beta * u[n+1]*v[i:0,n].dx(i:0,n) * dx 
> 
>> where n is the range of the summation. Setting n = 1 would 
then be 
>> for 2D, and n = 2 for 3D. It might then be an idea to make 
the 
>> spatial dimension available within ffc.  Garth 
> 
>I think this could be done much simpler with the separate 
>tensor/matrix class approach. You would do: 
> 
>def epsilon(u, n): 
>    eu = ffczeros(n, n) 
> 
>    for ii in range(n): 
>        for jj in range(n): 
>            eu[ii, jj] = 0.5 * (u[ii].dx(jj) + u[jj].dx(ii)) 
> 
>    return eu 
> 
>and then: 
> 
>a = epsilon(u, n) * ... -Beta * u[n+1] * ... 
> 
>Since slicing is supported by Numeric.array, and would be 
supported by 
>any other tensor/matrix representation we would choose, it 
would not 
>be necessary to add it to the FFC algebra. It would make the 
matrix 
>initialization much more compact though, if we could just 
write: 
> 
>epsilonu = matrix(u[i:0,n].dx(j:0,n) + u[i:0,n].dx(j:0,n)) 
> 
>(or some other slicing convention), but perhaps it's best to 
keep the 
>FFC algebra as simple as possible. 
> 
>Please tell me what you think about this, if it seems to 
cover the 
>necessary functionality etc. I'm going to try to use it for 
some other 
>problems and see if it's usable. 
> 
>  Johan 
>________________ 
># Copyright (c) 2004 Anders Logg ([email-masked]) 
># Licensed under the GNU GPL Version 2 
># 
># The bilinear form for a mass matrix. 
># 
># Compile this form with FFC: ffc Mass.form 
> 
>from Numeric import * 
>from math import * 
> 
>element = FiniteElement("Vector Lagrange", "tetrahedron", 1) 
> 
>v = BasisFunction(element) 
>u = BasisFunction(element) 
> 
>lmbda = Constant() # Lame coefficient 
>mu    = Constant() # Lame coefficient 
> 
>f = Function(element) # Source 
> 
>global zero 
>zero = 0.0 * Constant() 
> 
>def ffczeros(m, n): 
>    A = array((zero)) 
>    A = resize(A, (m, n)) 
>    A *= 0.0 
>     
>    return A 
> 
>def ffcidentity(m): 
>    A = myzeros(m, m) 
> 
>    for ii in range(m): 
>        for jj in range(m): 
>            if ii == jj: 
>                A[ii, jj] = 1.0 
>                     
>    return A 
> 
>def epsilon(u): 
>    eu = ffczeros(3, 3) 
> 
>    for ii in range(3): 
>        for jj in range(3): 
>            eu[ii, jj] = 0.5 * (u[ii].dx(jj) + u[jj].dx(ii)) 
> 
>    return eu 
> 
> 
>def E(e, lmbda, mu): 
>    Ee = 2.0 * multiply(mu, e) 
>    for ii in range(3): 
>        for jj in range(3): 
>            if ii == jj: 
>                Ee[ii, jj] += lmbda * trace(e) 
> 
>    return Ee 
> 
>sigma = E(epsilon(u), lmbda, mu) 
> 
>a = vdot(sigma, epsilon(v)) * dx 
>L = f[i] * v[i] * dx 
>________________ 
>_______________________________________________ 
>FFC-dev mailing list 
>[email-masked] 
>http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev

---

From: Anders Logg <[email-masked]>
Date: Wed, 31 Aug 2005 10:35:04 -0500

Of course, FFC is still very experimental but I don't think it would
be a good idea to base the FFC algebra on Numeric.array.

1. I think the best solution to Garth's problem would be to let the
constructor of an Index accept one optional extra argument to allow
specification of ranges. Then one can do stuff like

    i = Index(0, d)
    j = Index(0, d)
    a = u[i].dx(i) * v[j].dx(j)

This would force corresponding Indices to match and I think it would
allow Garth to do what he likes. Correct me if I'm wrong.

2. Then Johan's suggestion. If you write

then this is a scalar since repeating an Index twice means summation
over that Index. I think the tensor notation has proven to be quite
powerful so I would like to keep it.

On the other hand, you can do everything you need by indexing basis
functions and derivatives by integers (not Index) as you have noted.
This allows a user to define arbitrary operators as you have done:

def epsilon(u):
    return [ [u[i].dx(j) + u[j].dx(i) for i in range(d)] for j in range(d) ]

The only thing that needs to be added to FFC is then the identity as
far as I can see.

I'm planning to add a list of predefined operators and operands to FFC
to simplify and beautify the notation anyway, so this would fit well
into that category.

Operators that will be added include grad(), rot() and div().

Operands that will be added include h (mesh size) and d (space dimension).

Any other suggestions?

All these can be added on top of the existing machinery (I think). As
far as possible, I would like to keep the basic layer simple and add
complexity by combining the simple primitives.

/Anders

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/

---

From: [email-masked]
Date: Wed, 31 Aug 2005 17:57:50 +0200 (MEST)

> Operands that will be added include h (mesh size) and d (space dimension).

I find it useful to be able use some more geometry information (in
particular when we start integrating on surfaces):

vol - volume of cell
area(i) - area for local face i
hface(i) - h for local face i
...
same for local edges?

It may also be useful to be able to identify a local face from a
collection of vertices; that is, while visiting a cell in the assembly one
might want to get information about the local faces.

On the other hand, we may avoid a lot of this since we intend to assemble
separately over faces (and edges), and then vol and h may be reduced to
the corresponding functions for the face when looping over the faces.

/Johan

---

From: Anders Logg <[email-masked]>
Date: Wed, 31 Aug 2005 11:08:42 -0500

> hface(i) - h for local face i

I can understand if you need h (mesh size) and v (volume), but do you
really need area of a specific local face?

> might want to get information about the local faces.

I think all this messy stuff should be hidden. We need to do this when
we implement integrals over boundaries, but everything that needs to
be visible in FFC is the new operator

    *ds

Karin Kraft and Johan Jansson have already implemented a prototype
version of boundary integrals in DOLFIN as far as I understand, so I
only need to implement the remaining portions in FFC. I will probably
also need to request some new features from FIAT so I can tabulate
basis functions on faces/edges of the reference cells.

> the corresponding functions for the face when looping over the faces.

Do you need both vol and h at the same time in the stabilization?

/Anders

---

From: [email-masked]
Date: Wed, 31 Aug 2005 18:15:43 +0200 (MEST)

> Do you need both vol and h at the same time in the stabilization?

No, this is not for the stabilization. It is for applying conditions on
surfaces, such as weak boundary conditions. But if we do the assembly over
faces, we can take "h" to mean the "face h", so maybe this can remain
hidden. But typically we may want to formulate a form including both
volume and surface integrals, and then one would like to include an h_face
in the surface integral. But as I said; this may be avoided by hiding
"h_face" behind the regular "h". On the other hand, we have both dx and
ds.

/Johan

---

From: Anders Logg <[email-masked]>
Date: Wed, 31 Aug 2005 11:25:48 -0500

> "h_face" behind the regular "h". On the other hand, we have both dx and

ok, I see your point. We could probably hide it.

Does the h on the face really have to be specific to that face? Is it
not enough to use the same h (like the diameter)?

/Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 31 Aug 2005 18:22:15 +0200

> allow Garth to do what he likes. Correct me if I'm wrong.

This would be good.

---

From: Anders Logg <[email-masked]>
Date: Wed, 31 Aug 2005 11:31:34 -0500

> > allow Garth to do what he likes. Correct me if I'm wrong.

Fine. I'll get to it as soon as I can.

/Anders

---

From: Johan Jansson <[email-masked]>
Date: Wed, 31 Aug 2005 19:01:02 +0200

On Wed, Aug 31, 2005 at 10:35:04AM -0500, Anders Logg wrote:

...

> powerful so I would like to keep it.

I'm very sorry, but my expression was totally wrong. I copied it from
the wrong place, and I didn't catch it because it looked pretty
similar. What I meant was:

epsilonu = 0.5 * (u[i].dx(j) + u[j].dx(i))

Here the indices are not repeated, so it should not be a scalar. But
it's still not indexable in FFC. Is it meant to be indexable? Are
Products and Sums also supposed to be tensors?

Here's a minimal example:

# Tensor test

element = FiniteElement("Vector Lagrange", "tetrahedron", 1, 3)

gradu = u[i].dx(j)

tmp = gradu[0]

a = tmp * v[0] * dx

When I compile this, I get:

  File "Foo.py", line 36, in ?
    tmp = gradu[0]
  File "/usr/lib/python2.4/site-packages/ffc/compiler/algebra.py", line 383, in __getitem__
    w.basisfunctions[0] = w.basisfunctions[0][component]
  File "/usr/lib/python2.4/site-packages/ffc/compiler/algebra.py", line 209, in __getitem__
    raise RuntimeError, "Cannot pick component of scalar BasisFunction."
RuntimeError: Cannot pick component of scalar BasisFunction.

Perhaps it's just a bug in FFC? I haven't been able to follow the
error clearly. I see now that Product and Sum do have an indexing
operator, so perhaps I was just confused about this. If they are
indeed tensors, then we can just as well define tensor operators for
them.

>     return [ [u[i].dx(j) + u[j].dx(i) for i in range(d)] for j in range(d) ]

Here you construct a nested list structure which is a sort of tensor
representation though (like Numeric.array). Shouldn't this work (i and
j are Index)?:

def epsilon(u):
    return u[i].dx(j) + u[j].dx(i)

> Operands that will be added include h (mesh size) and d (space dimension).

It would also be good to have matrix operators, but that requires
matrix-valued algebra objects. Matrix multiplication, gradient of a
vector etc. But as I mentioned, as soon as indexing works, then the
rest of the operators should be trivial to implement.

> complexity by combining the simple primitives.

Yes, I agree with this.

  Johan

---

From: Johan Jansson <[email-masked]>
Date: Wed, 31 Aug 2005 21:20:03 +0200

On Wed, Aug 31, 2005 at 11:08:42AM -0500, Anders Logg wrote:

...

> basis functions on faces/edges of the reference cells.

Yes, we have a DOLFIN implementation in 2D with a manually generated
FFC-style header file for a reference case. I can send you the
relevant files with some explanatory comments if/when you start
working on this. Or if someone else wants to work on it perhaps, maybe
it's a good idea to try to have several developers working in every
component so that a single person doesn't get swamped?

I think Rob has already implemented tabulation over boundary entities,
check out this message:

http://www.fenics.org/pipermail/fiat-dev/2005-August/000058.html

  Johan

---

From: Anders Logg <[email-masked]>
Date: Wed, 31 Aug 2005 14:43:20 -0500

> epsilonu = 0.5 * (u[i].dx(j) + u[j].dx(i))

ok, now I see. Yes, this is a tensor (or more correctly a particular
component of a tensor). But it's not indexable since the indices are
already in there (i and j).

You can then declare epsilonv in the same way and do the tensor
contraction by

	    epsilonu * epsilonv

> Products and Sums also supposed to be tensors?

Everything in FFC is a scalar, except for vector-valued basis
functions which need to be indexed in the declaration of a form to
produce something that is scalar. More precisely: everything in FFC is
an already indexed tensor of arbitrary rank. So if you do

   u[i].dx(j).dx(k)

then you have a scalar component of a rank 3 tensor, but with an index
(i, j, k) that need not be fixed.

> RuntimeError: Cannot pick component of scalar BasisFunction.

When you do tmp = u[i].dx(j), you first pick one component of the
vector-valued basis function u, so then you have a scalar. When you
then take the derivative in the j direction you obtain a new scalar
function.

It's like when you write the same thing on a piece of paper. If you
have a vector x, then x_i is not a vector but component i on the
vector. Sometimes you use the notation (x_i) to denote the vector with
components x_i for i = 1,2,.... Maybe we could add something like this
to FFC, possibly a keyword vec so that

    vec(u[i])

is the thing you want.

> indeed tensors, then we can just as well define tensor operators for

You can take two vector-valued BasisFunctions u and v and add them to
get a vector-valued Sum:

    u + v

and then you can pick a component:

    (u + v)[i] = u[i] + v[i]

>     return u[i].dx(j) + u[j].dx(i)

No, same explanation as above.

> rest of the operators should be trivial to implement.

They are in a sense matrix-valued already (like u[i].dx(j)) but you
always work with the scalar components.

I hope things are a little more clear now. I need to write this up in
the FFC manual to avoid confusion.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Wed, 31 Aug 2005 14:56:32 -0500

> http://www.fenics.org/pipermail/fiat-dev/2005-August/000058.html

Great, I must have missed that.

Everything is prepared for boundary integrals in FFC. You can already
now write *ds, boundary terms are separated from interior terms, and
the DOLFIN format generates code for assembly of both terms.

The missing part is to call FIAT from FFC to tabulate on boundaries
and do the integration properly. I figure it will take me about a week
(from the starting point, which is not today) to figure out all the
details and get it running.

You can post the files on the list and then I (or anyone else) know
where to look.

/Anders

---

From: [email-masked]
Date: Thu, 1 Sep 2005 08:11:37 +0200 (MEST)

> not enough to use the same h (like the diameter)?

I guess it is possible to use the "volume h" as an approximation for "face
h", for now. Although, there could be a significant difference if we have
anisotropic meshes, for example. I have used this in boundary layer
computations, and then the h in one direction could be 10 times (or 100
times) the h in another direction. So we really also would like to be able
to decompose the "volume h" in different directions (h_x,h_y,h_z).

But the first priority now as I see it would be to get the regular "volume
h" = cell diameter to work.

/Johan

---

From: Johan Jansson <[email-masked]>
Date: Mon, 5 Sep 2005 00:07:51 +0200

...

> complexity by combining the simple primitives.

I think the tensor notation is converging to a natural choice of
representation. I've tested the Numeric approach a bit more, and it's
not actually necessary to have a zero and identity member in the FFC
algebra, the values 0.0 and 1.0 seem to be fine (so it's possible to
use the Numeric functions zeros() and identity() to construct zero and
identity matrices). Some languages have a dynamic type system, but it
seems Python goes even further than that, which is good for us in this
case at least.

I can now write the elasticity form like this (element definitions
etc. not included):

def grad(u):
    return [ [u[i].dx(j) for i in range(d)] for j in range(d) ]

def epsilon(u):
    return 0.5 * (grad(u) + transpose(grad(u)))

def E(e, lmbda, mu):
    Ee = 2.0 * multiply(mu, e) + \
         multiply(lmbda, multiply(trace(e), identity(3)))

    return Ee

All that's left for the notation to be identical (as identical as
ASCII can be) to written form is to move to operators instead of
functions for multiply(), vdot() etc. I'm not convinced that's a good
idea though, I prefer this notation, since it's unambigous. A "*"
operator can mean any kind of product. Since there are at least two
products which make sense for matrices in this context (and probably
more in other contexts), I don't think we have enough symbols. It
always takes me a few minutes to figure out what the dot operator
means in an article, perhaps I have finally found a remedy in this
notation convention. It will lead to an excessive amount of
parantheses though, but functional language programmers seem to be
able to handle that just fine (i.e. Lisp).

It also doesn't seem necessary to explicitly construct a Numeric.array
object using the array() function, the list representation works
fine. It does seem necessary for the Numeric.array operators to work
though, so perhaps it would be good to define grad() as:

def grad(u):
    return array(( [ [u[i].dx(j) for i in range(d)] for j in range(d) ] ))

Then we can use the slicing operators (and other operators) from
Numeric.array. For example, to extract the first 2 columns of grad(u),
we can write:

grad(u)[:, 0:2]

Perhaps this can be generalized even further by simply defining a
vec() function Anders mentioned before, and extending the concept to
matrix() and tensor(), which would simply construct a Numeric.array
from an FFC algebra object, i.e.:

i = Index()
j = Index()

gradu = matrix(u[i].dx(j))

which would be a much cleaner fit with matrix-valued elements. I.e. a
matrix-valued function "stress" could be manipulated in the same way:

i = Index()
j = Index()

sigma = matrix(stress[i, j])

This would mean fewer special cases, which I think is almost always
good. I still haven't grasped the Index concept fully though, so
perhaps there are other better ways to accomplish this.

If we need tensors of rank > 2, then we would have to replace
Numeric.array with something else. Is this relevant? Is there anyone
who would actually need that functionality?

  Johan

---

From: Anders Logg <[email-masked]>
Date: Tue, 6 Sep 2005 09:23:17 -0500

> a = vdot(sigma, epsilon(v)) * dx

This looks very good, but I'm not completely convinced... The problem
is that when you explicitly construct arrays and do scalar products
manually, by actually summing up the terms rather than working through
an Index, FFC has no idea that it should treat a series of terms as a
scalar products.

Consider the following examples with Poisson written in two different
ways:

1. Using the standard index notation:

    a = u.dx(i)*v.dx(i)*dx

Generated code:

    real G0_0_0 = map.det*(map.g00*map.g00 + map.g01*map.g01);
    ...

    block[0] = 4.999999999999998e-01*G0_0_0 + 4.999999999999997e-01*G0_0_1
             + 4.999999999999997e-01*G0_1_0 + 4.999999999999996e-01*G0_1_1;

2. Writing out the scalar product explicitly:

    a = (u.dx(0)*v.dx(0) + u.dx(1)*v.dx(1))*dx

    real G0_0_0 = map.det*map.g00*map.g00;
    ...

    block[0] = 4.999999999999998e-01*G0_0_0 + 4.999999999999997e-01*G0_0_1
             + 4.999999999999997e-01*G0_1_0 + 4.999999999999996e-01*G0_1_1
             + 4.999999999999998e-01*G1_0_0 + 4.999999999999997e-01*G1_0_1
             + 4.999999999999997e-01*G1_1_0 + 4.999999999999996e-01*G1_1_1;

In the first case, FFC generates a tensor representation of the
element matrix as one tensor product. In the second case, FFC
generates two tensor products, which is potentially less optimal.

There are two ways to solve this: either force a user to use the index
notation and extend the index notation with operators grad, rot, div
etc, or build in some functionality into FFC to automatically find
common factors,

> able to handle that just fine (i.e. Lisp).

This is one advantage of the index notation, since then everything is
spelled out in detail.

> perhaps there are other better ways to accomplish this.

I'm still thinking about the proper way to do this. The approach you
suggest would work very well, allowing maximum flexibility for the
user to define new operators as you have made, but there might be a
penalty in the form of less optimal code.

> who would actually need that functionality?

Numeric.array does tensors with arbitrary rank.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 6 Sep 2005 16:53:58 -0500

Rob and I discussed this earlier today and came to the conclusion that
new operators should be defined using Python built-in lists. This seems
to be both the most flexible approach and the most transparent.

With new operators defined in terms of Python lists, it is the
responsibility of FFC to process the given form and find the proper
tensor representation. In the case of Poisson's equation as discussed
in a previous post, this means that FFC will have to look at the two
terms and factor out the common reference tensor. I have started to
work out the details (computing signatures to identify common
tensors).

A clarification might be necessary. The syntax remains the same as
before, in particular the index notation with implicit summation.
The only thing that will change is that new operators (such as grad,
div, rot) will be added and they will be defined in terms of Python
lists. This already works (as Johan has noted), but is non-optimal.

/Anders

---

From: Johan Jansson <[email-masked]>
Date: Wed, 7 Sep 2005 14:50:09 +0200

> lists. This already works (as Johan has noted), but is non-optimal.

Ok, I didn't realize that it would break up the products. Of course we
shouldn't accept sub-optimal code generation. We don't have any real
time constraints at compile time, so let's make sure the operators
lead to the same code as writing them out manually explicitly would
do.

Great that you seem to have found a solution.

Also nice that Numeric.array supports tensors, I must have missed this
detail somehow in the documentation.

  Johan

---

From: Anders Logg <[email-masked]>
Date: Wed, 7 Sep 2005 08:31:01 -0500

On Wed, Sep 07, 2005 at 02:50:09PM +0200, Johan Jansson wrote:

> [...]
> [...]

I don't think it will add anything significant to the compile time. It
might even reduce it, since terms with common reference tensors will
be combined before the reference tensors are computed and actually
computing all the integrals is what takes most of the time (not the
form algebra).

/Anders