ffc team mailing list archive

[<kirby@xxxxxxxxxxxx>]

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