Linear elasticity with an anisotropic material

1

I want to solve a problem in linear elasticity with an anisotropic material:
example ∇⋅S=0.Here, S denotes the stress, with Sij=Cijkl:εkl. ε is the linearized Green-Lagrange strain ε(u)=1/2(∇u+∇uT) u is the displacement and C is the tensor of elasticity.When i try to apply the formula in the case of my equation it shows me this error i would like to know how to resolve this. your help would be welcome.Because I’m stuck on this point.

# --------------------
# Parameters
# --------------------
c_11 = 1.19e11
c_12 = 0.84e11
c_13 = 0.83e11
c_33 = 1.17e11
c_44 = 0.21e11

eps_11 = 8.15e-9
eps_33 = 6.58e-9

e_15 = 12.09
e_13 = -6.03
e_33 = 15.49

alpha= 1.267e5
beta= 6.259e-10
theta  = 0.5 
rho_Piezo = Constant(7700.)
# different tensor of  material 

# Elasticity tensor c^E
C = as_tensor([[c_11,c_12,c_13,0.,0.,0.],[c_12,c_11,c_13,0.,0.,0.],[c_13,c_13,c_33,0.,0.,0.],[0.,0.,0.,c_44,0.,0],[0.,0.,0.,0.,c_44,0.],[0.,0.,0.,0.,0.,(c_11-c_12)/2]])


#  coupling tensor e
E_piezo = as_tensor([[0. , 0. , 0., 0. , e_15, 0.],[0., 0., 0., e_15, 0., 0.],[e_13, e_13, e_33, 0., 0., 0.]])


# transpose form  tensor
E_t_piezo =as_tensor([[0., 0., e_13],[0., 0., e_13],[0., 0., e_33],[0., e_15, 0.],[e_15, 0., 0.],[0., 0., 0.]])

# Permittivitatstensor  epsilon^S
Eps_s = as_tensor([[eps_11, 0., 0.],[0., eps_11, 0.],[0., 0., eps_33]])
# Setup mesh
resolution = 30

geometry = Cylinder(Point(0, 0, 0), Point(0, 0, 0.001), 0.005, 0.005)

# Making Mesh (30 corresponds to the mesh density)
mesh = generate_mesh(geometry, resolution)
# --------------------
# Functions and classes
# --------------------

V = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)

W = FiniteElement("CG", mesh.ufl_cell(), 1)

M = FunctionSpace(mesh, MixedElement([V,V, W]))


# Get the needed boundaries: 
# 1. To apply a predefined stress on it We keep it simple: the top face sees the same stress vector. 
# 2. For zero Displacement we choose the bottom face of the disc
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.001)

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.)
# first define the function for the strain tensor in fenics notation
def epsilon(u):
    # note that it would be shorter to use sym(grad(u)) but the way used now
    # represents the formula
    return 1./2.*(grad(u) + (grad(u).T))
def Sigma(u):
    return dot(C, epsilon(u))
# define the potential
dt = 0.01
t = np.arange(0.00199578, 20.0, dt)

# Amplitude of the sine wave is sine of a variable like time
potential = np.sin(3*np.pi*(t-(20.0/2)))
phi = Expression("t", degree=0, t=potential[0])

# --------------------
# Boundary conditions
# --------------------
top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()

# Create the boundary condition for displacement
boundary = []
boundary.append(DirichletBC(M.sub(2), 0., bottom_boundary))
boundary.append(DirichletBC(M.sub(2), phi, top_boundary))

for i in range(len(t)):
    phi.t = t[i]
    print(assemble(phi*dx(domain=mesh)))
# As the stress is defined over the outer surface we need to create a 
# integral over that surface
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
top_boundary.mark(boundaries, 1)
bottom_boundary.mark(boundaries, 2)

dx = Measure("dx", domain=mesh)
ds = Measure("ds", subdomain_data=boundaries)
mf =  MeshFunction("size_t", mesh, mesh.topology().dim()-1)
mf.set_all(0)
top_boundary.mark(mf, 1)
bottom_boundary.mark(mf, 2)
# the initial values
zero1d=Expression('0.0', degree=1)
zero3d=Expression(('0.0','0.0','0.0'), degree=1)

u0=interpolate(zero3d,M.sub(0).collapse())
z0=interpolate(zero3d,M.sub(1).collapse())
phi0=interpolate(zero1d,M.sub(2).collapse())
# --------------------
# Function spaces
# --------------------

#Defining the "Trial" functions
u,z,phi=TrialFunctions(M)

#Defining the "Test" functions
v,y,w=TestFunctions(M)

a = rho_Piezo*inner(z, v)*dx
+ theta*dt*inner((Sigma(u)+dot(E_t_piezo, gad(phi) )),epsilon(v) )* dx
+ theta*dt*alpha*rho_Piezo*inner(z1, v) * dx
+ theta*dt*beta*inner(dot(C, epsilon(z)), epsilonB(v)) * dx
+ inner((dot(E_piezo, epsilon(u))-dot(Eps_s,grad(phi))), grad(w)) *dx
+ inner(u, y)*dx
- theta*dt*(inner(z, y))* dx

l = -rho_Piezo*inner(z0, v)*dx
+ theta*dt*inner((sigma(u0)+dot(E_t_piezo, gad(phi0))),epsilon(v)) * dx
+ theta*dt*alpha*rho_Piezo*inner(z0, v)*dx
+ theta*dt*Beta*inner(dot(C, epsilon(z0)),epsilon(v)) * dx
- inner(u0, z)*dx
- theta*dt*(inner(z0, y)) * dx

Dimension mismatch in dot product.

---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-64-0058a608ed80> in <module>
     12 
     13 a = rho_Piezo*inner(z, v)*dx
---> 14 + theta*dt*inner((Sigma(u)+dot(E_t_piezo, gad(phi) )),epsilon(v) )* dx
     15 + theta*dt*alpha*rho_Piezo*inner(z1, v) * dx
     16 + theta*dt*beta*inner(dot(C, epsilon(z)), epsilonB(v)) * dx

<ipython-input-59-b7ff8cb334a9> in Sigma(u)
      1 def Sigma(u):
----> 2     return dot(C, epsilon(u))

~/anaconda3/envs/fenicsfix2/lib/python3.8/site-packages/ufl/operators.py in dot(a, b)
    185     if a.ufl_shape == () and b.ufl_shape == ():
    186         return a * b
--> 187     return Dot(a, b)
    188 
    189 

~/anaconda3/envs/fenicsfix2/lib/python3.8/site-packages/ufl/tensoralgebra.py in __new__(cls, a, b)
    206                   "got arguments with ranks %d and %d." % (ar, br))
    207         if not (scalar or ash[-1] == bsh[0]):
--> 208             error("Dimension mismatch in dot product.")
    209 
    210         # Simplification

~/anaconda3/envs/fenicsfix2/lib/python3.8/site-packages/ufl/log.py in error(self, *message)
    170         "Write error message and raise an exception."
    171         self._log.error(*message)
--> 172         raise self._exception_type(self._format_raw(*message))
    173 
    174     def begin(self, *message):

UFLException: Dimension mismatch in dot product.
2 Likes
2

Hi, you are using a matrix representation for the 4th order tensor while still using a tensor representation for the strain. You should convert the latter to a vector notation to be able to compute the dot product.

3

thank you for your answer.
I have done this but I receive errors here is the form in voigt notation. That’s why I wonder where is the problem.Because I do not really know what to do anymore

from fenics import *
from mshr import *
import numpy as np
from math import *
from matplotlib import pyplot as plt
from dolfin import *
import matplotlib.cm as cm
import matplotlib.tri as tri
from mpl_toolkits.axes_grid1 import make_axes_locatable
import ufl
c_11 = 1.19e11
c_12 = 0.84e11
c_13 = 0.83e11
c_33 = 1.17e11
c_44 = 0.21e11

eps_11 = 8.15e-9
eps_33 = 6.58e-9

e_15 = 12.09
e_13 = -6.03
e_33 = 15.49

alpha= 1.267e5
beta= 6.259e-10
theta  = 0.5 
rho_Piezo = Constant(7700.)
# different tensor of  material 

# Elasticity tensor c^E
C = as_tensor([[c_11,c_12,c_13,0.,0.,0.],[c_12,c_11,c_13,0.,0.,0.],[c_13,c_13,c_33,0.,0.,0.],[0.,0.,0.,c_44,0.,0],[0.,0.,0.,0.,c_44,0.],[0.,0.,0.,0.,0.,(c_11-c_12)/2]])


#  coupling tensor e
E_piezo = as_tensor([[0. , 0. , 0., 0. , e_15, 0.],[0., 0., 0., e_15, 0., 0.],[e_13, e_13, e_33, 0., 0., 0.]])


# transpose form  tensor
E_t_piezo =as_tensor([[0., 0., e_13],[0., 0., e_13],[0., 0., e_33],[0., e_15, 0.],[e_15, 0., 0.],[0., 0., 0.]])

# Permittivitatstensor  epsilon^S
Eps_s = as_tensor([[eps_11, 0., 0.],[0., eps_11, 0.],[0., 0., eps_33]])
# Setup mesh
resolution = 30

geometry = Cylinder(Point(0, 0, 0), Point(0, 0, 0.001), 0.005, 0.005)

# Making Mesh (30 corresponds to the mesh density)
mesh = generate_mesh(geometry, resolution)
# --------------------
# Functions and classes
# --------------------

V = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)

W = FiniteElement("CG", mesh.ufl_cell(), 1)

M = FunctionSpace(mesh, MixedElement([V,V, W]))


# Get the needed boundaries: 
# 1. To apply a predefined stress on it We keep it simple: the top face sees the same stress vector. 
# 2. For zero Displacement we choose the bottom face of the disc
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.001)

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.)
# first define the function for the strain tensor in fenics notation
#definition of the different components of my Stress tensor equation (sigma=C^E:Bu-e_transpose.E)

# rewrite the tensor into a vector using Voigt notation
def strain3voigt(ten):
    # FEniCS does not know anything about Voigt notation, 
    # so one need to access the components directly as eps[0, 0] etc.   
    return as_vector([ten[0,0],ten[1,1],ten[2,2],2*ten[1,2],2*ten[0,2],2*ten[0,1]])

# rewrite a vector into a tensor using Voigt notation
def voigt3stress(vec):
    # FEniCS does not know anything about Voigt notation, 
    # so one need to rewrite the vector by hand into a tensor.   
    return as_tensor([[vec[0], vec[5], vec[4]], 
                      [vec[5], vec[1], vec[3]],
                      [vec[4], vec[3], vec[2]]])

# first define the function for the strain tensor in fenics notation
def B(u):
    # note that it would be shorter to use sym(grad(u)) but the way used now
    # represents the formula
    return 1./2.*(grad(u) + (grad(u).T))
 
#define the electric field vector                 
def E_field_val_vector(Ei): 
    return as_vector([Ei[0], Ei[1], Ei[2]])
                  
#define the electric field in relation to potential
def elect_field_E(phi):
    return -grad(phi)
                   
# sigma for coupled linear elasicity and electrostatics                  
def sigma(u,phi):
                  
    return voigt3stress(dot(C,strain3voigt(B(u)))+ E_t_piezo*E_field_val_vector(elect_field_E(phi)))
                  
#the electric displacements vector
def disp_D(Di): 
                  
    return as_vector(Di[0], Di[1], Di[2])

#definition of the different components of my electric displacements equation (D=e^s:Bu+eps^S.E) 
# elect_disp_D for coupled linear elasicity and electrostatics   
def elect_disp_D(u,phi):
                  
    elect_disp=(E_piezo*strain3voigt(B(u)) - Eps_s*_elect_field_E(phi))
                  
    return disp_D(elect_disp)
# define the potential
dt = 0.01
t = np.arange(0.00199578, 20.0, dt)

# Amplitude of the sine wave is sine of a variable like time
potential = np.sin(3*np.pi*(t-(20.0/2)))
phi = Expression("t", degree=0, t=potential[0])

# --------------------
# Boundary conditions
# --------------------
top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()

# Create the boundary condition for displacement
boundary = []
boundary.append(DirichletBC(M.sub(2), 0., bottom_boundary))
boundary.append(DirichletBC(M.sub(2), phi, top_boundary))

for i in range(len(t)):
    phi.t = t[i]
    print(assemble(phi*dx(domain=mesh)))
# As the stress is defined over the outer surface we need to create a 
# integral over that surface
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
top_boundary.mark(boundaries, 1)
bottom_boundary.mark(boundaries, 2)

dx = Measure("dx", domain=mesh)
ds = Measure("ds", subdomain_data=boundaries)
mf =  MeshFunction("size_t", mesh, mesh.topology().dim()-1)
mf.set_all(0)
top_boundary.mark(mf, 1)
bottom_boundary.mark(mf, 2)
# the initial values
zero1d=Expression('0.0', degree=1)
zero3d=Expression(('0.0','0.0','0.0'), degree=1)

u0=interpolate(zero3d,M.sub(0).collapse())
z0=interpolate(zero3d,M.sub(1).collapse())
phi0=interpolate(zero1d,M.sub(2).collapse())
# --------------------
# Function spaces
# --------------------

#Defining the "Trial" functions
u1,z1,phi1=TrialFunctions(M)

#Defining the "Test" functions
v,y,w=TestFunctions(M)
# --------------------
# Weak form
# --------------------
tmps = Function(M)


'''
#a := ρ(z 1 , v)L 2 + Θ · dt((c E Bu1 + eT ∇φ1 ),Bv) L 2 + Θ · dt · α · ρ(z 1 , v)L 2
#+Θ · dt · β(cE Bz 1 , Bv) L 2 + (eBu1 − ε S ∇φ 1 , ∇w) L 2
#+(u 1 , y) L 2 − Θ · dt(z 1 , y) L 2
'''

a = rho_Piezo*inner(z1, v)*dx
+ theta*dt*inner(sigma(u1, phi1),strain3voigt(B(v)) * dx
+ theta*dt*alpha*rho_Piezo*inner(z1, v) * dx
+ theta*dt*beta*inner(dot(C, strain3voigt(B(z1))), strain3voigt(B(v))) * dx
+ inner(elect_disp_D(u1, phi1), grad(w)) *dx
+ inner(u1, y)*dx
- theta*dt*(inner(z1, y))* dx

'''
l := −ρ(z 0 , v) L 2 + Θ · dt((cE Bu 0 + eT ∇φ 0 ), Bv) L 2 + Θ · dt · α · ρ(z 0 , v) L 2
+Θ · dt · β(c E Bz 0 , Bv) L 2 − (u 0 , y) L 2 − Θ · dt(z 0 , y) L 2 ,
'''

                 
l = -rho_Piezo*inner(z0, v)*dx
+ theta*dt*inner(sigma(u0, phi0), strain3voigt(B(v))) * dx
+ theta*dt*alpha*rho_Piezo*inner(z0, v)*dx
+ theta*dt*Beta*inner(dot(C, strain3voigt(B(z0))), strain3voigt(B(v))) * dx
- inner(u0, z)*dx
- theta*dt*(inner(z0, y))) * dx

 File "<ipython-input-38-393ce687fbac>", line 24
    '''
    ^
SyntaxError: invalid syntax

#assemble only A and L 
A=assemble(a)
L=assemble(-l)
#apply A and L
for bc in boundary:
    bc.apply(A,L)
#solve the weak equation

solve(A,tmps.vector(),L)

u_result,z_result,phi_result=tmps.split()
FunctionAssigner([M.sub(0).collapse(),M.sub(1).collapse(),M.sub(2).collapse()],tmps.function_space())
assign(tmps.sub(0), u0)
assign(tmps.sub(1), z0)
assign(tmps.sub(2),phi0)
vtkfile = File('phi_rult.pvd')
vtkfile << phi_result
4

In a you compute the inner product between sigma and epsilon with tensor representation for sigma and vector Voigt representation for epsilon, that’s why.

5

Please it’s actually this ((c E Bu1 + eT ∇φ1 ),Bv) L 2 scalar product that I want to do that’s why I took the inner. Can you please propose me an alternative to do this.

6

@stefan: Would ufl.indices be useful in this situation? It is a function that generates index objects which can be used to define any indicial product, e.g. C_{ijkl} u_{k,l} v_{i,j} could be written as:

from fenics import *
from ufl import indices
import numpy as np

# Create four free indices for defining the scalar product
i,j,k,l = indices(4)

# Fourth-order stiffness tensor
C = Constant(np.ones((3,3,3,3)))

# Variational form
mesh = UnitCubeMesh(4,4,4)
V = VectorFunctionSpace(mesh, "Lagrange", 2)
u = TrialFunction(V)
v = TestFunction(V)

a = (C[i,j,k,l] * sym(grad(u))[k,l] * sym(grad(v))[i,j])*dx

Note that instead of np.ones((3,3,3,3)) you will need to populate a numpy array with the components of your anisotropic stiffness tensor.

7

Just replace the strain3voigt(B(v)) by B(v). There were many left-over typos in the code, the following one runs for me. But please check that it’s indeed what you want:

from fenics import *
from mshr import *
import numpy as np
from math import *
from matplotlib import pyplot as plt
from dolfin import *
import matplotlib.cm as cm
import matplotlib.tri as tri
from mpl_toolkits.axes_grid1 import make_axes_locatable
import ufl

c_11 = 1.19e11
c_12 = 0.84e11
c_13 = 0.83e11
c_33 = 1.17e11
c_44 = 0.21e11

eps_11 = 8.15e-9
eps_33 = 6.58e-9

e_15 = 12.09
e_13 = -6.03
e_33 = 15.49

alpha = 1.267e5
beta = 6.259e-10
theta = 0.5
rho_Piezo = Constant(7700.0)
# different tensor of  material

# Elasticity tensor c^E
C = as_tensor(
    [
        [c_11, c_12, c_13, 0.0, 0.0, 0.0],
        [c_12, c_11, c_13, 0.0, 0.0, 0.0],
        [c_13, c_13, c_33, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, c_44, 0.0, 0],
        [0.0, 0.0, 0.0, 0.0, c_44, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0, (c_11 - c_12) / 2],
    ]
)


#  coupling tensor e
E_piezo = as_tensor(
    [
        [0.0, 0.0, 0.0, 0.0, e_15, 0.0],
        [0.0, 0.0, 0.0, e_15, 0.0, 0.0],
        [e_13, e_13, e_33, 0.0, 0.0, 0.0],
    ]
)


# transpose form  tensor
E_t_piezo = as_tensor(
    [
        [0.0, 0.0, e_13],
        [0.0, 0.0, e_13],
        [0.0, 0.0, e_33],
        [0.0, e_15, 0.0],
        [e_15, 0.0, 0.0],
        [0.0, 0.0, 0.0],
    ]
)

# Permittivitatstensor  epsilon^S
Eps_s = as_tensor([[eps_11, 0.0, 0.0], [0.0, eps_11, 0.0], [0.0, 0.0, eps_33]])
# Setup mesh
resolution = 10

geometry = Cylinder(Point(0, 0, 0), Point(0, 0, 0.001), 0.005, 0.005)

# Making Mesh (30 corresponds to the mesh density)
mesh = generate_mesh(geometry, resolution)
# --------------------
# Functions and classes
# --------------------

V = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)

W = FiniteElement("CG", mesh.ufl_cell(), 1)

M = FunctionSpace(mesh, MixedElement([V, V, W]))


# Get the needed boundaries:
# 1. To apply a predefined stress on it We keep it simple: the top face sees the same stress vector.
# 2. For zero Displacement we choose the bottom face of the disc
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.001)


class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.0)


# first define the function for the strain tensor in fenics notation
# definition of the different components of my Stress tensor equation (sigma=C^E:Bu-e_transpose.E)

# rewrite the tensor into a vector using Voigt notation
def strain3voigt(ten):
    # FEniCS does not know anything about Voigt notation,
    # so one need to access the components directly as eps[0, 0] etc.
    return as_vector(
        [ten[0, 0], ten[1, 1], ten[2, 2], 2 * ten[1, 2], 2 * ten[0, 2], 2 * ten[0, 1]]
    )


# rewrite a vector into a tensor using Voigt notation
def voigt3stress(vec):
    # FEniCS does not know anything about Voigt notation,
    # so one need to rewrite the vector by hand into a tensor.
    return as_tensor(
        [[vec[0], vec[5], vec[4]], [vec[5], vec[1], vec[3]], [vec[4], vec[3], vec[2]]]
    )


# first define the function for the strain tensor in fenics notation
def B(u):
    # note that it would be shorter to use sym(grad(u)) but the way used now
    # represents the formula
    return 1.0 / 2.0 * (grad(u) + (grad(u).T))


# define the electric field vector
def E_field_val_vector(Ei):
    return as_vector([Ei[0], Ei[1], Ei[2]])


# define the electric field in relation to potential
def elect_field_E(phi):
    return -grad(phi)


# sigma for coupled linear elasicity and electrostatics
def sigma(u, phi):

    return voigt3stress(
        dot(C, strain3voigt(B(u))) + E_t_piezo * E_field_val_vector(elect_field_E(phi))
    )


# the electric displacements vector
def disp_D(Di):

    return as_vector([Di[0], Di[1], Di[2]])


# definition of the different components of my electric displacements equation (D=e^s:Bu+eps^S.E)
# elect_disp_D for coupled linear elasicity and electrostatics
def elect_disp_D(u, phi):

    elect_disp = E_piezo * strain3voigt(B(u)) - Eps_s * elect_field_E(phi)

    return disp_D(elect_disp)


# define the potential
dt = 0.01
t = np.arange(0.00199578, 20.0, dt)

# Amplitude of the sine wave is sine of a variable like time
potential = np.sin(3 * np.pi * (t - (20.0 / 2)))
phi = Expression("t", degree=0, t=potential[0])

# --------------------
# Boundary conditions
# --------------------
top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()

# Create the boundary condition for displacement
boundary = []
boundary.append(DirichletBC(M.sub(2), 0.0, bottom_boundary))
boundary.append(DirichletBC(M.sub(2), phi, top_boundary))

for i in range(len(t)):
    phi.t = t[i]
    print(assemble(phi * dx(domain=mesh)))
# As the stress is defined over the outer surface we need to create a
# integral over that surface
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top_boundary.mark(boundaries, 1)
bottom_boundary.mark(boundaries, 2)

dx = Measure("dx", domain=mesh)
ds = Measure("ds", subdomain_data=boundaries)
mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
mf.set_all(0)
top_boundary.mark(mf, 1)
bottom_boundary.mark(mf, 2)
# the initial values
zero1d = Expression("0.0", degree=1)
zero3d = Expression(("0.0", "0.0", "0.0"), degree=1)

u0 = interpolate(zero3d, M.sub(0).collapse())
z0 = interpolate(zero3d, M.sub(1).collapse())
phi0 = interpolate(zero1d, M.sub(2).collapse())
# --------------------
# Function spaces
# --------------------

# Defining the "Trial" functions
u1, z1, phi1 = TrialFunctions(M)

# Defining the "Test" functions
v, y, w = TestFunctions(M)
# --------------------
# Weak form
# --------------------
tmps = Function(M)


a = (
    rho_Piezo * inner(z1, v) * dx
    + theta * dt * inner(sigma(u1, phi1), B(v)) * dx
    + theta * dt * alpha * rho_Piezo * inner(z1, v) * dx
    + theta * dt * beta * inner(dot(C, strain3voigt(B(z1))), strain3voigt(B(v))) * dx
    + inner(elect_disp_D(u1, phi1), grad(w)) * dx
    + inner(u1, y) * dx
    - theta * dt * (inner(z1, y)) * dx
)


l = (
    -rho_Piezo * inner(z0, v) * dx
    + theta * dt * inner(sigma(u0, phi0), B(v)) * dx
    + theta * dt * alpha * rho_Piezo * inner(z0, v) * dx
    + theta * dt * beta * inner(dot(C, strain3voigt(B(z0))), strain3voigt(B(v))) * dx
    - inner(u0, y) * dx
    - theta * dt * (inner(z0, y)) * dx
)
# assemble only A and L
A = assemble(a)
L = assemble(-l)
# apply A and L
for bc in boundary:
    bc.apply(A, L)
# solve the weak equation

solve(A, tmps.vector(), L)

u_result, z_result, phi_result = tmps.split()
FunctionAssigner(
    [M.sub(0).collapse(), M.sub(1).collapse(), M.sub(2).collapse()],
    tmps.function_space(),
)
assign(tmps.sub(0), u0)
assign(tmps.sub(1), z0)
assign(tmps.sub(2), phi0)
vtkfile = File("phi_rult.pvd")
vtkfile << phi_result
4 Likes
8

I sincerely do not understand in fact I apply at the beginning a sinusiodal voltage as defined in the program on the upper side of my material and normally I had to have that by visualizing on fenics but it appeared all blue that is to say it receives no voltage. do you have an explanation of this phenomenon.because here is what I would like to have

havegetsame
havegetsame859×452 7.13 KB

but here is what I obtain by simulating on paraview.
inversesolve
inversesolve764×384 5.06 KB

9

Hi, have you figured out why? I encountered the same issue when trying to write a code based on this demo.
Thank you

10

Hi all, I was having a similar problem. and the solution provided by @bleyerj solved it.

In the beginning, it produced unexpected results. Then I refined the mesh and got the expected results. I guess with Anisotropic problems the mesh has to be appropriately refined.

Thanks, @bleyerj for sharing your knowledge with us, and thanks @stefan for explaining the problem appropriately.