Modal participation factor for mixed elements

1

Hi all,
I am using mixed-elements for doing modal analysis. I want to extend it to compute the modal participation factor. The geo file is attached here: geometry - Google Drive . The MWE is-

from dolfin import *
import numpy as np
import json

mesh = Mesh()
with XDMFFile("mesh/tetra.xdmf") as infile:
    infile.read(mesh)

# Function space, trial and test functions
elem_P = VectorElement('Lagrange', mesh.ufl_cell(), 2)
elem_R = VectorElement('Real', mesh.ufl_cell(), 0)
elem = MixedElement([elem_P, elem_R, elem_R])
V = FunctionSpace(mesh, elem)

u_, U_, C_ = TrialFunctions(V)
du, dU, dC = TestFunctions(V)

# Material properties, constitutive law
E, nu = 200e9, 0.3  
rho = Constant(7850)
mu = Constant(E/2./(1+nu))
lmbda = Constant(E*nu/(1+nu)/(1-2*nu))

def eps(v):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

# Boundary conditions
support1 = CompiledSubDomain("near(x[0], 1.0, tol) && on_boundary", tol=1e-14)
support2 = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14)
support3 = CompiledSubDomain("near(x[1], 1.0, tol) && on_boundary", tol=1e-14)
support4 = CompiledSubDomain("near(x[1], 0.0, tol) && on_boundary", tol=1e-14)

bc_l = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support1)
bc_r = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support2)
bc_t = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support3)
bc_b = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support4)
bc=[bc_l,bc_r,bc_t,bc_b]

# Surface markers
force_face = CompiledSubDomain("(x[0] - 0.41478)*(x[0] - 0.41478) + (x[1] - 0.46787)*(x[1] - 0.46787)  < 0.16*0.16")
line_mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
line_mf.set_all(0)

force_face.mark(line_mf, 2)

# Parameters of the point mass
hole_radius = 0.15
plate_thk = 0.1
mass = 1000   # kg
area = 2*np.pi*hole_radius*plate_thk

M = Constant(mass)
A = Constant(area)

# Mass and stiffness forms
ds = Measure('ds', mesh, subdomain_data=line_mf)
k_form = inner(sigma(du),eps(u_))*dx - (inner(C_, du-dU) + inner(u_-U_, dC))*ds(2) 
m_form = rho*inner(u_,du)*dx + M/A*inner(U_,du)*ds(2)

K = PETScMatrix()
assemble(k_form, tensor=K)

M = PETScMatrix()
assemble(m_form, tensor=M)

for bcs in bc:
    bcs.zero(M)
    bcs.apply(K)

# Solution
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.

n_eig = 5   # number of eigenvalues
print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(n_eig)

# Set up file for exporting results
file_results = XDMFFile("modal_analysis.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
eigenmodes = []
eig_vals = []
eig_v = []
# Extraction
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    
    # 3D eigenfrequency
    freq_3D = np.sign(r)*np.sqrt(np.abs(r))/2/pi
    print("Solid FE ({0:d}): {1:8.5f} [Hz]".format(i+1, freq_3D))
    
    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector "+str(i))
    eigenmode.vector()[:] = rx
    file_results.write(eigenmode.split()[0], 0)
    eigenmodes.append(rx)
    eig_v.append(eigenmode)
    eig_vals.append(freq_3D)

V_new = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
u = Function(V_new, name="Unit displacement")
u.interpolate(Constant((0, 1, 0)))


combined_mass = 0
for i, xi in enumerate(eig_v):
    qi = assemble(action(action(m_form, u), xi))
    mi = assemble(action(action(m_form, xi), xi))
    meff_i = (qi / mi) ** 2
    total_mass = assemble(rho * dx(domain=mesh))

    print("-" * 50)
    print("Mode {}:".format(i + 1))
    print("  Modal participation factor: {:}".format(qi))
    print("  Modal mass: {:}".format(mi))
    print("  Effective mass: {:}".format(meff_i))

The demo for calculation of modal participation factor can be found here