Defining a parameter in element coordinate system

1

Hi

I am solving a Hyperelasticity problem (compressible neo-Hookean model). I want to add a new term to the strain energy function: C[0,0]*cos(theta)

Where (\theta) is defined in XY plane in global coordinate system:

3

Now I want to evaluate the parameter \theta in elementwise fashion. To be more specific I want to consider the orientation of the XY plane with elements orientation. This figure shows what I am looking at:

2
2897×704 79.6 KB

In general, I want to map the XY plane (in global coordinate system) on the outer facet of each element (e.g. I need to define the \theta in local or element coordinate system) and \theta is constant in each element. That is why I think I need to solve it elementwise (I have never done it before!). I just need to figure out how I should consider the orientation of XY plane with orientation of elements. Does anybody know how I can do that in FEniCS?
This is my minimal work for a simple unit cube (The cube example is for simplicity but I am looking for a general solution working for any type of geometry like a sphere etc):

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(2, 2, 2)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

left = Left()
right = Right()

boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': 2})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 2})

bcl = DirichletBC(V, Constant((0,0,0)), boundaries, 1)

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

force  = Constant((1.0,  0.0, 0.0))

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

theta = 0

# The last term is the new term added to strain energy function that includes theta parameter
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 + C[0,0]*cos(theta)

Pi = psi*dx - dot(force, u)*ds(2)

F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bcl, J=J,)
2

I don’t see a general way to do this that would be insensitive to node numbering (which should not affect the physics). I suspect there is some confusion in your modeling that should be cleared up first. If you’re trying to add some spatially-varying orthotropic behavior to the material, you might look into specifying a unit vector field, \mathbf{N}, and including a term in your energy density involving I_4 = \mathbf{N}\cdot\mathbf{C}\mathbf{N} (i.e., I_4 = dot(N,C*N) in UFL). For \mathbf{N} constant on each element, you could represent it in a DG0 VectorFunctionSpace.

3

Thanks for your response.That is right. I am trying to define a spatially-varying behavior to my material. The unit vectors in my problem define the element plane. Please see this for clarification:

tt
tt561×710 103 KB

The red vectors define the plane of that element. In other words when the element rotates, the vectors also rotate with elements. I am trying to define such a unit vector field! Sorry if I was not able to explain it better.
Regarding what you suggested, how can I specify that N is constant on each element? Could you please provide a simple example?

4

I think the difficulty here is that, for a tetrahedral element in a general unstructured mesh, it is not possible to uniquely define “the plane of that element”. Thus, given a single scalar angle \theta for each element, it is not possible to define \mathbf{N} without additional assumptions or information. (For instance, in the above diagram, what if one picked a different edge as the x-axis of the local coordinate system? How does one know which orthogonal direction to take as the y-axis?)

The most convenient way to set the material direction field \mathbf{N} really depends on what is upstream of the analysis, providing the original problem data on material anisotropy (e.g., SALS data providing a collagen fiber orientation in a soft tissue sample, design drawings specifying the orientation of fiber reinforcement in a composite structure, etc.).

1 Like
5

Thanks for your clarification. It cleared up some of my confusions. Based on what you explained, I think it is better to use Hexahedron element instead of Tetrahedron. At this point I am looking at a structured mesh. So I came up with an idea. I made a structured mesh using Hexahedron element. Then I extracted the node numbers and their coordinates on each element. Then I figured out the nodes ordering on each element. Then I calculated the unit vectors defining the plane of each element. Please look at this figure showing what I did:
S
In the above figure, the red vectors (I call it a vector) correspond to the X direction of each element plane and the blue vectors ( b vector) correspond to the Y direction of each element plane. Then I saved a and b vectors for each element in a txt file. The txt file looks like this: (3 numbers in front of each vector define the direction of each vector in 3D)

element id = 1
	a vector: -0.203,0.023,-0.979
	b vector: -0.992,0.129,0.000

element id = 2
	a vector: -0.329,0.019,-0.944
	b vector: -0.982,0.128,0.139

element id = 3
	a vector: -0.448,0.015,-0.894
	b vector: -0.950,0.126,0.285
.
.
.

If the element plane for all elements was XY plane, then the parameter N could be defined as : N = [cos \theta ,sin \theta,0].
What I am thinking about, is to define each element as a subdomain and then iterate the variational form over all subdomains (e.g. elements) but for each subdomain I use a transformation based on unit vectors of that element to define appropriate N for that element.
As I was able to find the unit vectors of each element defining the element plane, could you please tell me if this idea may work or not? (With assumption of using structured mesh and Hexahedron elements)
Thanks again!

6

Hi,
sorry but is your geometry a portion of a hollow sphere ? If so, it seems to me that your local unit vectors are just the spherical coordinate basis vectors e_\phi and e_\theta which you can easily define analytically.
Otherwise, if you have the outward unit normal N, you could just define you red vector as (e_z \times N)/\|e_z \times N\|.

7

Hi
Yes the geometry is a part of a hollow sphere. The problem of defining the unit vectors for each element (e.g. The unit vectors that define the element plane) has already been solved as I stated above. The issue that still remains is how to define the parameter N that varies spatially. Let me explain it a little more:
As David indicated, I have this term in my strain energy function : N.C.N where C is the right cauchy deformation tensor and N is defined as : N = [cos \theta , sin\theta , 0]. In other words, N defines the orientation of the fiber in the element plane. For example if \theta=0 they are all oriented like the red vectors and if \theta=90 they are all oriented like the blue vectors. Now I have all of the unit vectors (red and blue) for each element. I am just trying to figure out how I can define the term N = [cos \theta , sin\theta , 0] in UFL so it recognizes that N has been defined in the element plane when the unit vectors of that element are already known.
For example when an element is parallel to the XY plane, the red unit vector is: [1,0,0] and blue unit vector is: [0,1,0] but for the next element the red vector are different (For instance the red is: [-0.2,0.3,-0.9] and the blue is: [-0.9,0.1,0.000]).
I hope I could explain my problem properly.
I appreciate any hint or example helping me resolve this issue.

8

Well if you know the unit vectors a and b in each element, then you just have:

theta = Constant(pi/4)
N = cos(theta)*a+sin(theta)*b
1 Like
9

Thank Jeremy. It makes sense.
When it comes to implementation, my solution (As explained above) is to define each element as a subdomain (e.g. each cube is considered and marked as a subdomain) and then with having the a and b vector for that subdomain I define:
N = cos(theta)*a+sin(theta)*b for that subdomain and in my varational form I have to iterate over all subdomains (e.g. dx(1), dx(2) , …)
Do you think this solution may work or there might be a better way to do it?
Thanks!

10

Hi, yes I think that a and b can be represented using a UserExpression or a DG0 Function from reading your txt file.

11

I tried to use DG0 space in the above code by doing :

V = VectorFunctionSpace(mesh, 'DG',0)

However I got this error:

`Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.

I guess I should stick to the CG space and read the vectors from txt file as you suggested.
Anyways, thank you guys for your nice suggestions.`