##
Thanks for the suggestion, I tried to follow your supply code and i have implemented the above equations, But I get an error stating that :

*** Error: Unable to assemble system.

*** Reason: expected a bilinear form for a.

*** Where: This error was encountered inside SystemAssembler.cpp.

*** Process: 0

*** DOLFIN version: 2019.2.0.dev0

*** Git changeset: unknown

I have attached my code below, It would be appreciable to point out the mistake:

from dolfin import *

import matplotlib.pyplot as plt

import numpy as np

import timeit

#units: m, kg, s, V, K

################################## Mesh part ##########################

“”"

$PhysicalNames

2D 2 “top”

2D 3 “bottom”

3D 1 “cube”

“”"

#
load the mesh data from Gmsh

mesh_file = ‘cube.h5’

mesh = Mesh()

hdf = HDF5File(mesh.mpi_comm(), mesh_file, ‘r’)

hdf.read(mesh, ‘/mesh’, False)

cells = MeshFunction(‘size_t’, mesh, dim=3)

hdf.read(cells, ‘/subdomains’)

facets = MeshFunction(‘size_t’, mesh, dim=2)

hdf.read(facets, ‘/facets’)

#
JZ: use the Measure function

dx = Measure(‘dx’, domain=mesh, subdomain_data=cells)

ds = Measure(‘ds’, domain=mesh, subdomain_data=facets)

#
Function spaces

N = FacetNormal(mesh)

Ue = VectorElement(“CG”, mesh.ufl_cell(), 2) # displacement vector element

Ve = FiniteElement(“CG”, mesh.ufl_cell(), 1) # voltage finite element

#u,phi

W = FunctionSpace(mesh, MixedElement([Ue, Ve]))

#
Parameters of Piezomaterial (PIC-255)

#
--------------------

c_11 = 1.23e11

c_12 = 7.67e10

c_13 = 7.025e10

c_33 = 9.711e10

c_44 = 2.226e10

eps_11 = 8.234

eps_33 = 7.588

e_15 = 11.91

e_13 = -7.15

e_33 = 13.75

rho_Piezo = Constant(7800) #in kg/m3

#
Order of matrices similar to ANSYS notation: x, y, z, xy, yz,xz

#
IEEE notation: x, y, z, yz, xz, xy

#
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]])

#
piezoelectric coupling tensor e

et_piezo = as_tensor([[0. , 0. , 0., 0. , 0., e_15],

[0., 0., 0., 0., e_15, 0.],

[e_13, e_13, e_33, 0., 0., 0.]])

#
transpose form of piezo tensor

e_piezo = as_tensor([[0., 0., e_13],

[0., 0., e_13],

[0., 0., e_33],

[0., 0, 0.],

[0., e_15, 0.],

[e_15, 0., 0.]])

#
Permittivitatstensor epsilon^S

eps_s = as_tensor([[eps_11, 0., 0.],

[0., eps_11, 0.],

[0., 0., eps_33]])

#
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],ten[0,1],ten[1,2],ten[0,2]])

#
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[3], vec[5]],

[vec[3], vec[1], vec[4]],

[vec[5], vec[4], vec[2]]])

#
first define the function for the strain tensor in fenics notation

def B(u):

return 0.5*(nabla_grad(u) + nabla_grad(u).T)

#define the electric field vector

def E_field_vector(Ei):

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

#define the electric field in relation to potential

def E_field(v):

return -grad(v)

#definition of the different components of my Stress tensor equation (sigma=C^E:Bu-e_piezo.E)

def sigma(u, v):

return voigt3stress(dot(C, strain3voigt(B(u))) + e_piezo * E_field_vector(E_field(v)))

#
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 D_field(u, v):

D = et_piezo * strain3voigt(B(u)) - eps_s * E_field(v)

return disp_D(D)

#Defining the “Trial” functions

w = Function(W)

(U, V) = split(w)

#Defining the test functions

(utest, vtest) = TestFunctions(W)

#
dirichlet boundary condition for fixed support

bcs = [DirichletBC(W.sub(0), Constant((0.,0.,0.)), facets, 2)]

#
Forms & matrices

U_form = inner(sigma(U, V), B(utest)) * dx

V_form = inner(D_field(U,V), grad(vtest)) * dx

K_form = U_form + V_form

M_form = rho_Piezo*inner(U, utest)*dx

L = utest[0]*dx

K, M = PETScMatrix(), PETScMatrix()

b = PETScVector()

assemble_system(K_form,L,bcs,A_tensor=K)

assemble_system(M_form,L,bcs,A_tensor=M)

[bc.zero(M) for bc in bcs]

#
--------------------Eigen-solver--------------------

eigensolver = SLEPcEigenSolver(K, M)

eigensolver.parameters[‘problem_type’] = ‘gen_hermitian’

eigensolver.parameters[‘spectral_transform’] = ‘shift-and-invert’

eigensolver.parameters[‘spectral_shift’] = 0.

N_eig = 6 # number of eigenvalues

print(“Computing {} first eigenvalues…”.format(N_eig))

eigensolver.solve(N_eig)