hello everyone i have a problem with my code if i apply the project function it does not accept mixed functions i would like to know if anyone can help me
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.)
Rayleigh damping coefficients
eta_m = Constant(0.1)
#eta_k = Constant(0.1)
#Generalized-alpha method parameters
alpha_m = Constant(0.5)
alpha_f = Constant(0.5)
gamma = Constant(0.5+alpha_f-alpha_m)
beta = Constant((gamma+0.5)*1/4)
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_11 = 8.15e-9
eps_33 = 6.58e-9
eps_s = as_matrix([[eps_11, 0., e_31],
[0., eps_11, e_31],
[0., 0., eps_11]])
Setup mesh
resolution = 30
cylinder = Cylinder(Point(0, 0, 0), Point(0, 0, 1), 5, 5)
Omega = generate_mesh (cylinder, resolution)
x = Omega.coordinates()
x[:, 0] *= 0.1
x[:, 1] *= 0.1
x[:, 2] *= 0.01
Omega.bounding_box_tree().build(Omega)
V = VectorFunctionSpace(Omega, “CG”,3)
VV = VectorElement(“CG”, Omega.ufl_cell(), 3)
W = FunctionSpace(Omega,“CG”,1)
WW = FiniteElement(“CG”, Omega.ufl_cell(), 1)
M = FunctionSpace(Omega, 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”, Omega, Omega.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=3)
tmps.split()
tmps()
u0=project(zero3d,M.sub(0))
z0=project(zero3d,M.sub(1))
phi0=project(zero1d,M.sub(2))
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)
RuntimeError Traceback (most recent call last)
in
48 zero3d=Expression((‘0.0’,‘0.0’,‘0.0’), degree=3)
49
—> 50 u0=project(zero3d,M.sub(0))
51 z0=project(zero3d,M.sub(1))
52 phi0=project(zero1d,M.sub(2))
~/anaconda3/envs/fenicsfix2/lib/python3.8/site-packages/dolfin/fem/projection.py in project(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)
135 # Solve linear system for projection
136 if function is None:
→ 137 function = Function(V)
138 cpp.la.solve(A, function.vector(), b, solver_type, preconditioner_type)
139
~/anaconda3/envs/fenicsfix2/lib/python3.8/site-packages/dolfin/function/function.py in init(self, *args, **kwargs)
214 if len(args) == 1:
215 # If passing only the FunctionSpace
→ 216 self._cpp_object = cpp.function.Function(V._cpp_object)
217 elif len(args) == 2:
218 if isinstance(args[1], cpp.la.GenericVector):
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to create function.
*** Reason: Cannot be created from subspace. Consider collapsing the function space.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset: 77c758717da1cb8017bf95109363bb1a91a4f8e5
*** -------------------------------------------------------------------------
#The d-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
‘’’
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,dtinner((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,betainner(cas_voigt_vector(epsilon(z1)), as_voigt_vector(epsilon(u)))*L2 ))
- inner(e_piezoas_voigt_vector(epsilon(u1))-eps_sgrad(phi1), grad(phi))*L2
- inner(u1,z)L2
-dot(0.5, dtinner(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,dtinner((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,betainner(cas_voigt_vector(epsilon(z0)), as_voigt_vector(epsilon(u)))*L2 ))
- inner(u0,z)*L2
- dot(0.5, dt*inner(z0,z)*L2)
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
NameError Traceback (most recent call last)
in
12
13
—> 14 a = rho_Piezoinner(z1,u)dx
15 + dot(0.5,dtinner((c(as_voigt_vector(epsilon(u1))) + et_piezo*grad(phi1)),as_voigt_vector(epsilon(u)))L2)
16 + dot(0.5,dot(dt,dot(alpha_m,rho_Piezoinner(z1,u)*L2)))
NameError: name ‘z1’ is not defined
Get x values of the sine wave
dt = 0.01
t = np.arange(0, 10.0, dt);
Amplitude of the sine wave is sine of a variable like time
potential = np.sinc(3np.pi(t-(20.0/2)))
Plot a sine wave using time and amplitude obtained for the sine wave
plt.plot(t, potential)
Give a title for the sine wave plot
plt.title(‘Sine wave’)
Give x axis label for the sine wave plot
plt.xlabel(‘Time’)
Give y axis label for the sine wave plot
plt.ylabel(‘Potential’)
plt.grid(True, which=‘both’)
plt.axhline(y=0, color=‘k’)
Display the sine wave
plt.show()
def potential_result(phi_top, phi_bot):
bc_top = DirichletBC(M.sub(2),phi_top,boundary_parts,1)
bc_bot=DirichletBC(M.sub(2),phi_bot,boundary_parts
bcp = []
bcp.append(bc_top)
bcp.append(bc_bot)
# Assembly the system
A=assemble(a)
L=assemble(-l)
bcp[0].apply(A)
bcp[1].apply(L)
#solve(A, vector, L)
solve(A,tmps.vector(),L)
u_result,z_result,phi_result=tmps.split()
assigner1=FunctionAssigner([V,V,W],tmps.function_space())
assigner1.assign([u0,z0,phi0],tmps)
print(Omega.hmax())
print(Omega.hmin())
plot(Omega)
plot(u_result)
vtkfile = File('inversepiezo_test01.pvd')
vtkfile << u_result