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
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],
]
)
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],
]
)
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],
]
)
Eps_s = as_tensor([[eps_11, 0.0, 0.0], [0.0, eps_11, 0.0], [0.0, 0.0, eps_33]])
def strain3voigt(ten):
return as_vector(
[ten[0, 0], ten[1, 1], ten[2, 2], 2 * ten[1, 2], 2 * ten[0, 2], 2 * ten[0, 1]]
)
def voigt3stress(vec):
return as_tensor(
[[vec[0], vec[5], vec[4]], [vec[5], vec[1], vec[3]], [vec[4], vec[3], vec[2]]]
)
def B(u):
return 1.0 / 2.0 * (grad(u) + (grad(u).T))
def E_field_val_vector(Ei):
return as_vector([Ei[0], Ei[1], Ei[2]])
def elect_field_E(phi):
return -grad(phi)
return voigt3stress(
dot(C, strain3voigt(B(u))) + E_t_piezo * E_field_val_vector(elect_field_E(phi))
)
def disp_D(Di):
return as_vector([Di[0], Di[1], Di[2]])
def elect_disp_D(u, phi):
elect_disp = E_piezo * strain3voigt(B(u)) - Eps_s * elect_field_E(phi)
return disp_D(elect_disp)
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two slave edges
return bool ((near(x[0], 0) or near(x[1], 0) or near(x[2], 0)) and
(not ((near(x[0], a) and near(x[2], a)) or
(near(x[0], a) and near(x[1], a)) or
(near(x[1], a) and near(x[2], a)))) and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
#### define mapping for a single point in the box, such that 3 mappings are required
if near(x[0], a) and near(x[1], a) and near(x[2],a):
y[0] = x[0] - a
y[1] = x[1] - a
y[2] = x[2] - a
##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
if near(x[0], a) and near(x[2], a):
y[0] = x[0] - a
y[1] = x[1]
y[2] = x[2] - a
elif near(x[1], a) and near(x[2], a):
y[0] = x[0]
y[1] = x[1] - a
y[2] = x[2] - a
elif near(x[0], a) and near(x[1], a):
y[0] = x[0] - a
y[1] = x[1] - a
y[2] = x[2]
#### right maps to left: left/right is defined as the x-direction
elif near(x[0], a):
y[0] = x[0] - a
y[1] = x[1]
y[2] = x[2]
### back maps to front: front/back is defined as the y-direction
elif near(x[1], a):
y[0] = x[0]
y[1] = x[1] - a
y[2] = x[2]
#### top maps to bottom: top/bottom is defined as the z-direction
else:
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - a
mesh = UnitCubeMesh(20, 20, 20)
pb = PeriodicBoundary()
CGV_elem = VectorElement(“Lagrange”, mesh.ufl_cell(), 1)
CGS_elem = FiniteElement(“CG”, mesh.ufl_cell(), 1)
R_element = FiniteElement(“R”, mesh.ufl_cell(), 0)
W_elem = MixedElement([CGV_elem, CGS_elem, R_element, R_element, R_element, R_element])
W = FunctionSpace(mesh, W_elem, constrained_domain=pb)
(u, v, j0, j1, j2, l) = TrialFunctions(W)
(alpha, chi, q0, q1, q2, m) = TestFunctions(W)
a = inner(stress(u, v), strain3voigt(B(alpha)))dx(mesh) + inner(elect_disp_D(u, v), g(chi))dx(mesh) + j0alpha[0]dx(mesh) + j1alpha[1]dx(mesh) + j2alpha[2]dx(mesh) + q0u[0]dx(mesh) + q1u[1]dx(mesh) + q2u[2]dx(mesh) + lchidx(mesh) + mvdx(mesh)
L = 0*dx(mesh)
w = Function(W)
solve(a == L, w)
(Psi00, psi00, j0, j1, j2, l) = w.split()