Local coordinate system changes

Asked by Nick Davies on 2013-02-22

Hi Everyone,
I have a problem where I need to use different orientations of coordinates at each element. It is a 3D problem with the mesh the shape of a cone, and tetrahedral elements. The global coordinate system is Cartesian with (0,0,0) at the centre of the bottom polygon. The problem is the youngs and shear modulus’s which can be obtained for the material are orientated with what are usually called ‘radial’, ‘tangential’ and ‘vertical’ axis. The radial axis points toward the point (0,0,z) and the vertical axis points in the direction of the cells (not necessarily vertical in the global sense) with the tangential axis being normal to them. So these change throughout the mesh so that the radial is always pointing toward the centre.

In two dimensions this would be equitant to if we had say a pentagon with six verities (one at 0,0 and each of the other 5 at the corners in Cartesian coordinates) and 5 triangular cells each with the vertex 0,0 and two corners of the pentagon. Now we want to compress the pentagon but we only have modulus’s of elasticity for along an axis running from 0,0 to the corners (radial) and normal to the radial axis is the tangential axis, also with a modulus.

My question is, how can I implement a map in dolfin which will allow for the global coordinate system to be used to position the nodes and for the force to be applied to a boundary, but then when the deformation is calculated it is calculated using the directions the modulus’s of elasticity for? (ie radial/tangential directions for the local coordinates)
I hope this makes sense.
Also can I upload images to Launchpad? That would make the explanation clearer
Thanks for any help
Nick

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Marie Rognes
Solved:
2013-03-16
Last query:
2013-03-16
Last reply:
2013-03-12
Johan Hake (johan-hake) said : #1

This should be possible. Try defining your local coordinate system as a
transformation tensor which you multiply with your modulus vector.

Johan

On 02/22/2013 01:20 AM, Nick Davies wrote:
> New question #222534 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/222534
>
> Hi Everyone,
> I have a problem where I need to use different orientations of coordinates at each element. It is a 3D problem with the mesh the shape of a cone, and tetrahedral elements. The global coordinate system is Cartesian with (0,0,0) at the centre of the bottom polygon. The problem is the youngs and shear modulus’s which can be obtained for the material are orientated with what are usually called ‘radial’, ‘tangential’ and ‘vertical’ axis. The radial axis points toward the point (0,0,z) and the vertical axis points in the direction of the cells (not necessarily vertical in the global sense) with the tangential axis being normal to them. So these change throughout the mesh so that the radial is always pointing toward the centre.
>
> In two dimensions this would be equitant to if we had say a pentagon with six verities (one at 0,0 and each of the other 5 at the corners in Cartesian coordinates) and 5 triangular cells each with the vertex 0,0 and two corners of the pentagon. Now we want to compress the pentagon but we only have modulus’s of elasticity for along an axis running from 0,0 to the corners (radial) and normal to the radial axis is the tangential axis, also with a modulus.
>
> My question is, how can I implement a map in dolfin which will allow for the global coordinate system to be used to position the nodes and for the force to be applied to a boundary, but then when the deformation is calculated it is calculated using the directions the modulus’s of elasticity for? (ie radial/tangential directions for the local coordinates)
> I hope this makes sense.
> Also can I upload images to Launchpad? That would make the explanation clearer
> Thanks for any help
> Nick
>
>

Nick Davies (ntd14) said : #2

Thanks for that Johan,
Is there any reason anyone can think of that this same approach wouldn't work applied directly to the stiffness matrix? This means there is no need to transform the Poisson ratios separately.

something such as:

D(global) = G(transpose) D(local) G
Where D is the the 6x6 stiffness matrix and G is a 6x6 transformation matrix

Is dolfin able to take the transpose of a matrix consisting of expressions of x[0..2]?

Marie Rognes (meg-simula) said : #3

I don't quite understand your follow-up question: could you please provide a small amount of code to make this more concrete?

Nick Davies (ntd14) said : #4

Hi Marie,
Sorry your right that wasn't very clear. Currently I have implemented the solution the way Johan suggested, however I have not been able to work out how to make dolfin create a matrix whose elements are a function of x[0..2] here is a shorted version of how I am currently implementing the solution, the problem is how slow it makes the code.

Er = 2.1
Et = 2.3
… defining all youngs and shear mods along with the Poisson ratios

### creating expressions for the entries of the stiffness matrix, (although constant currently these will need to vary depending on their global position later, this is why they are set as expressions)
class CMrrExpression(Expression):
 def eval(self, values, x):
  values[0] = (-1+Vtv*Vvt)*Er/sig
CM11 = CMrrExpression()

class CMttExpression(Expression):
 def eval(self, values, x):
  values[0] = (-1+Vrv*Vvr)*Et/sig
CM22 = CMttExpression()

 …......... (for each stiffness matrix entry)

#calculating angles to transform from rtv tp xyz
class trans_angles(Expression):
 def eval(self, values, x):
  Ar = math.atan2(x[0],x[1])
  ## cockwise radial to cartesian
  values[0] = Ar
  values[1] = Ar + math.pi/2
.... for the 9 angles all of these angles will potentially need to very with global position in the future
 def value_shape(self):
 return (9,)

trans = trans_angles()

Axr = trans[0]
Axt = trans[1]
...for all 9 angles

G11 = cos(Axr)*cos(Axr)
G12 = cos(Ayr)*cos(Ayr)
….. for each entry of the 36 element transformation matrix

GT11 = G11
GT12 = G21
….for each element of the transpose of the transformation matrix

CMT11 = CM11*GT11+CM12*GT21+CM13*GT31+CM14*GT41+CM15*GT51+CM16*GT61
CMT12 = CM11*GT12+CM12*GT22+CM13*GT32+CM14*GT42+CM15*GT52+CM16*GT62
….....for each element of the new stiffness matrix

CMF11 = CMT11*G11+CMT12*G21+CMT13*G31+CMT14*G41+CMT15*G51+CMT16*G61
CMF12 = CMT11*G12+CMT12*G22+CMT13*G32+CMT14*G42+CMT15*G52+CMT16*G62
…. for each element to convert from rtv coordinates to xyz

…...a bunch of stuff to define forces and set up the cauchy green tensors to calc strain energy etc.

#calc strain density function, anisotropic

psi = 0.5*CMF11*E[0,0]*E[0,0] + 0.5*CMF12*E[0,0]*E[1,1] + 0.5*CMF13*E[0,0]*E[2,2] + 0.5*CMF14*E[0,0]*2*E[0,1] + 0.5*CMF15*E[0,0]*2*E[0,2] + 0.5*CMF16*E[0,0]*2*E[1,2]\
…....... the rest of the terms to calculate the strain energy density

Pi = psi*dx - dot(T, u)*ds(1) #total potential energy
F = derivative(Pi, u, v)
J = derivative(F, u, du)

I would like to change it all into matrices rather than the manual element by element implementation that I am currently using to increase the speed. Is there any way simple way of doing this? When I have tried to do this then I couldn't get matrices to have elements which were expressions.

Thanks for any help

Best Marie Rognes (meg-simula) said : #5

First, for speed, subclassing the Expression class in Python is quite a slow approach. Take a look at the documentation for Expression for alternatives.

Second, here is a code snippet example that illustrates how to transform local tensors via transformation tensors, all which can depend on spatial coordinates, using matrix notation. Hope this might help.
--

from dolfin import *

# Define components of transformation matrix: these can be Functions,
# Expressions, Constants as you like it
g00 = Expression("x[0]")
g01 = g00
g10 = g00
g11 = g00

# Define transformation matrix via as_matrix function in UFL
G = as_matrix(((g00, g01), (g10, g11)))

# Some modulus tensor defined via as_vector and diag functions in UFL
m0 = Expression("x[1]")
m1 = 0.2
M = diag(as_vector((m0, m1)))

# Define some expression in terms of transformed modulus tensor
N = G*M*G.T

Third, if you find the jit compiler is slow there is a big improvement on
the way in development. It is missing some features still, but will be
announced soon.

Martin
Den 12. mars 2013 09:11 skrev "Marie Rognes" <
<email address hidden>> følgende:

> Question #222534 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/222534
>
> Status: Open => Answered
>
> Marie Rognes proposed the following answer:
> First, for speed, subclassing the Expression class in Python is quite a
> slow approach. Take a look at the documentation for Expression for
> alternatives.
>
> Second, here is a code snippet example that illustrates how to transform
> local tensors via transformation tensors, all which can depend on spatial
> coordinates, using matrix notation. Hope this might help.
> --
>
> from dolfin import *
>
> # Define components of transformation matrix: these can be Functions,
> # Expressions, Constants as you like it
> g00 = Expression("x[0]")
> g01 = g00
> g10 = g00
> g11 = g00
>
> # Define transformation matrix via as_matrix function in UFL
> G = as_matrix(((g00, g01), (g10, g11)))
>
> # Some modulus tensor defined via as_vector and diag functions in UFL
> m0 = Expression("x[1]")
> m1 = 0.2
> M = diag(as_vector((m0, m1)))
>
> # Define some expression in terms of transformed modulus tensor
> N = G*M*G.T
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Nick Davies (ntd14) said : #7

Thanks Marie Rognes, that solved my question.

Nick Davies (ntd14) said : #8

Thanks for your help. Looking forward to the update in the JIT compiler