I am solving equation lineare piezoelectric and simulate it on paraweiv but since I have a problem on my result and I don’t know where I made mistake I would like you to help me. Here is the code and also the result on paraweiv
equation I have to solve 
my code
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
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_31 = -6.03
e_33 = 15.49
alpha_M= 1.267e5
alpha_K= 6.259e-10
rho_Piezo = Constant(7700.)
c = as_matrix([[c_11, c_12, c_13, 0., 0., 0.],
[0., 0., 0., 0., 0., 0.],
[0., 0., c_33, 0., 0., 0.],
[0., 0., 0., c_44, 0., 0],
[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., (c_11-c_12)/2]])
e_piezo = as_matrix([[0. , 0. , 0., 0. , e_15, 0.],
[0., 0., 0., e_15, 0., 0.],
[e_31, e_31, e_33, 0., 0., 0.]])
# transpose form of piezo tensor
et_piezo = as_matrix([[0., 0., e_31],
[0., 0., e_31],
[0., 0., e_33],
[0., e_15, 0.],
[e_15, 0., 0.],
[0., 0., 0.]])
eps_s = as_matrix([[eps_11, 0., e_31],
[0., eps_11, e_31],
[0., 0., eps_11]])
resolution = 5
cylinder = Cylinder(Point(0, 0, 0), Point(0, 0, 1), 0.5, 0.5)
mesh= generate_mesh (cylinder, resolution)
dof_coordinates= mesh.coordinates()
dof_coordinates[:, 0] *= 0.1
dof_coordinates[:, 1] *= 0.1
dof_coordinates[:, 2] *= 0.01
mesh.bounding_box_tree().build(mesh)
V = VectorFunctionSpace(mesh, "CG",3)
VV = VectorElement("CG", mesh.ufl_cell(), 3)
WW = FiniteElement("CG", mesh.ufl_cell(), 1)
M = FunctionSpace(mesh, MixedElement([VV,VV, WW]))
#initialize boundary condition
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
top = Top()
bot = Bottom()
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_parts.set_all(0)
top.mark(boundary_parts, 1)
bot.mark(boundary_parts, 2)
ds = Measure("ds", subdomain_data=boundary_parts)
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())
# build the trial and test function spaces over a domain
#Defining the "Trial" functions
u1,z1,phi1=TrialFunctions(M)
#Defining the "Test" functions
u,z,phi=TestFunctions(M)
#The e-form can thus be reformulated using Voigt notation
# rewrite the tensor into a vector using Voigt notation
def as_voigt_vector(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],ten[1,2],ten[0,2],ten[0,1]])
# rewrite a vector into a tensor using Voigt notation
def from_voigt_vector(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 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) + np.transpose(grad(u)))
tmps = Function(M)
dt = 0.01
# the bilineare equation
'''
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 + (eBu 1 − ε S ∇φ 1 , ∇w) L 2
+(u 1 , yi) L 2 − Θ · dt(z 1 , yi) L 2
'''
L2 = norm(tmps, 'L2')
a = rho_Piezo*inner(z1,u)*dx
+ dot(0.5,dt*inner((c*(as_voigt_vector(epsilon(u1))) + et_piezo*grad(phi1)),as_voigt_vector(epsilon(u)))*L2)
+ dot(0.5,dot(dt,dot(alpha_M,rho_Piezo*inner(z1,u)*L2)))
+ dot(0.5,dot(dt,alpha_K*inner(c*as_voigt_vector(epsilon(z1)), as_voigt_vector(epsilon(u)))*L2 ))
+ inner(e_piezo*as_voigt_vector(epsilon(u1))-eps_s*grad(phi1), grad(phi))*L2
+ inner(u1,z)*L2
-dot(0.5, dt*inner(z1,z)*L2)
'''
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,u)*dx
+ dot(0.5,dt*inner((c*(as_voigt_vector(epsilon(u0))) + et_piezo*grad(phi0)),as_voigt_vector(epsilon(u)))*L2)
+ dot(0.5,dot(dt,dot(alpha_M,rho_Piezo*inner(z0,u)*L2)))
+ dot(0.5,dot(dt,alpha_K*inner(c*as_voigt_vector(epsilon(z0)), as_voigt_vector(epsilon(u)))*L2 ))
- inner(u0,z)*L2
- dot(0.5, dt*inner(z0,z)*L2)
phi_top = 0.1
phi_bot = 0.0
bc_top = DirichletBC(M.sub(2),phi_top,boundary_parts,1)
bc_bot=DirichletBC(M.sub(2),phi_bot,boundary_parts,2)
bcp = [bc_top, bc_bot]
A=assemble(a)
L=assemble(-l)
for bc in bcp:
bc.apply(A, L)
solve(A,tmps.vector(),L)
u_result,z_result,phi_result=tmps.split()
vtkfile = File('piezo_phi1.pvd')
vtkfile << phi_result my result on Paraweiv is not good I would like to have help on this point
