ffc team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxxxxxxx>]

Quoting 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.

This would be good.


> 
> 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/
> 
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
>