ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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

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

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

On Tue, Aug 30, 2005 at 05:48:24AM -1100, kirby@xxxxxxxxxxxx wrote:

> 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 <johanjan@xxxxxxxxxxxxxxxx>   
> >Subject: Re: [FFC-dev] Tensors   
> >To: "Garth N. Wells" <g.n.wells@xxxxxxxxxxxxxxx> 
> >Cc: ffc-dev@xxxxxxxxxx 
> > 
> >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 (logg@xxxxxxxxx) 
> ># 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 
> >FFC-dev@xxxxxxxxxx 
> >http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev 
> 
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
> 

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