Hi,
for reference, look at tests/test_msm_symbolic.py.
Currently, the symbolic derivation of right-hand sides corresponding to manufactured solutions works only for scalar variables - briefly, we can test just Poisson or diffusion problems, where the symbolic description of what each term does is: Laplace term: 'c * div( grad( u ) )' (c is scalar) Diffusion term: 'div( K * grad( u ) )' (K is matrix)
I would like to discuss here how to approach the problems involving vector fields, e.g. the linear elasticity. The linear elastic term is defined (LaTeX notation, summation convention - same indices sum together) as
d/dx_j D_{ijkl} e_{kl}(u), or div( D e ), where e is the small strain tensor.
I am thinking of two ways of describing this in SfePy:
- matrix notation (this is used for scalar terms) expr = """ e = 1/2 * (grad( vec( u ) ) + grad( vec( u ) ).T) D = map( D_sym ) s = D * e out = div( s ) """
... ok, but some more complex terms may be difficult to express in this was. Anyway I will try this first.
- indicial notation """ e[i,j] = 1/2 * (derj + deri) map = ? D[i,j,k,l] = map( D_sym ) s[i,j] = D[i,j,k,l] * e[k,l] out[i] = sum( derj) """
... better in a sense that one sees immediately how the operations look. But as I would like to avoid writing for i in range( 3 ): for j in range( 3 ): e[i,j] = ... i.e. writing the cycles explicitly, some parsing would have to be included.
-> my question (Ondrej) - would it be possible to add to sympy some capabilities of Einstein summation convention? maybe a module that would expand names like e_i_j into proper cycles.
r.
On Wed, Jun 25, 2008 at 2:55 PM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi,
for reference, look at tests/test_msm_symbolic.py.
Currently, the symbolic derivation of right-hand sides corresponding to manufactured solutions works only for scalar variables - briefly, we can test just Poisson or diffusion problems, where the symbolic description of what each term does is: Laplace term: 'c * div( grad( u ) )' (c is scalar) Diffusion term: 'div( K * grad( u ) )' (K is matrix)
I would like to discuss here how to approach the problems involving vector fields, e.g. the linear elasticity. The linear elastic term is defined (LaTeX notation, summation convention - same indices sum together) as
d/dx_j D_{ijkl} e_{kl}(u), or div( D e ), where e is the small strain tensor.
I am thinking of two ways of describing this in SfePy:
- matrix notation (this is used for scalar terms) expr = """ e = 1/2 * (grad( vec( u ) ) + grad( vec( u ) ).T) D = map( D_sym ) s = D * e out = div( s ) """
... ok, but some more complex terms may be difficult to express in this was. Anyway I will try this first.
- indicial notation """ e[i,j] = 1/2 * (derj + deri) map = ? D[i,j,k,l] = map( D_sym ) s[i,j] = D[i,j,k,l] * e[k,l] out[i] = sum( derj) """
... better in a sense that one sees immediately how the operations look. But as I would like to avoid writing for i in range( 3 ): for j in range( 3 ): e[i,j] = ... i.e. writing the cycles explicitly, some parsing would have to be included.
-> my question (Ondrej) - would it be possible to add to sympy some capabilities of Einstein summation convention? maybe a module that would expand names like e_i_j into proper cycles.
Yes, it's on my TODO list. I will implement the same as GiNaC has here:
http://www.ginac.de/tutorial/Indexed-objects.html#Indexed%20objects
That way you will just write something like
var("e D i j k l") Indexed(D, i, j, k, l) * Indexed(e, k, l)
and it will do the right thing. Then you will tell it to expand itself and it will probably call your routine to return something to be substituted for
D_0000 D_0100 ... e_00 e_01 e_11 ...
Plus I will implement Matrix support so that for vectors and 2 rank tensors this works out of the box.
Ondrej
On Thu, Jun 26, 2008 at 4:51 PM, Ondrej Certik <ond...@certik.cz> wrote:
On Wed, Jun 25, 2008 at 2:55 PM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Hi,
for reference, look at tests/test_msm_symbolic.py.
Currently, the symbolic derivation of right-hand sides corresponding to manufactured solutions works only for scalar variables - briefly, we can test just Poisson or diffusion problems, where the symbolic description of what each term does is: Laplace term: 'c * div( grad( u ) )' (c is scalar) Diffusion term: 'div( K * grad( u ) )' (K is matrix)
I would like to discuss here how to approach the problems involving vector fields, e.g. the linear elasticity. The linear elastic term is defined (LaTeX notation, summation convention - same indices sum together) as
d/dx_j D_{ijkl} e_{kl}(u), or div( D e ), where e is the small strain tensor.
I am thinking of two ways of describing this in SfePy:
- matrix notation (this is used for scalar terms) expr = """ e = 1/2 * (grad( vec( u ) ) + grad( vec( u ) ).T) D = map( D_sym ) s = D * e out = div( s ) """
... ok, but some more complex terms may be difficult to express in this was. Anyway I will try this first.
- indicial notation """ e[i,j] = 1/2 * (derj + deri) map = ? D[i,j,k,l] = map( D_sym ) s[i,j] = D[i,j,k,l] * e[k,l] out[i] = sum( derj) """
... better in a sense that one sees immediately how the operations look. But as I would like to avoid writing for i in range( 3 ): for j in range( 3 ): e[i,j] = ... i.e. writing the cycles explicitly, some parsing would have to be included.
-> my question (Ondrej) - would it be possible to add to sympy some capabilities of Einstein summation convention? maybe a module that would expand names like e_i_j into proper cycles.
Yes, it's on my TODO list. I will implement the same as GiNaC has here:
http://www.ginac.de/tutorial/Indexed-objects.html#Indexed%20objects
That way you will just write something like
var("e D i j k l") Indexed(D, i, j, k, l) * Indexed(e, k, l)
and it will do the right thing. Then you will tell it to expand itself and it will probably call your routine to return something to be substituted for
D_0000 D_0100 ... e_00 e_01 e_11 ...
Plus I will implement Matrix support so that for vectors and 2 rank tensors this works out of the box.
This is tracked at:
http://code.google.com/p/sympy/issues/detail?id=16
as you can see, this is a really old issue. :)
Ondrej
Ondrej Certik wrote:
Yes, it's on my TODO list. I will implement the same as GiNaC has here:
http://www.ginac.de/tutorial/Indexed-objects.html#Indexed%20objects
That way you will just write something like
var("e D i j k l") Indexed(D, i, j, k, l) * Indexed(e, k, l)
and it will do the right thing. Then you will tell it to expand itself and it will probably call your routine to return something to be substituted for
D_0000 D_0100 ... e_00 e_01 e_11 ...
Plus I will implement Matrix support so that for vectors and 2 rank tensors this works out of the box.
This is tracked at:
http://code.google.com/p/sympy/issues/detail?id=16
as you can see, this is a really old issue. :)
wow, history breathes on me! I look forward having this in SymPy.
r.
participants (2)
-
Ondrej Certik -
Robert Cimrman