Tutorial for `PETScMatrix` and `PETScVector` object

1

Hello everyone,
I am new to FEniCS. I have to calculate modal participation factor of a dynamic system. For this I am following Modal analysis of an elastic structure code as a reference till assemblage of stiffness and mass matrix. The system that I am working on, has very large degrees of freedom. So, calculating with numpy library in python is giving Large memory: bad allocation error. So I am bound to do calculations in PETScMatrix and PETScVector object. Can you please suggest any tutorial or link for understanding basic operations in PETScMatrix and PETScVector object .

2

Consider this excellent example in the numerical tours tutorial by @bleyerj.

3

Hi,
I just added a small part on computing modal participation factors and effective masses at the end of the modal analysis demo
You don’t necessarily need to work with PETSc objects but just UFL forms, although the latter is also possible.

1 Like
4

Thank you Sir,
The demo is working well for a small degrees of freedom system, but for large degrees of freedom system, I am facing this error.

image (46)
image (46)717×355 109 KB

5

I don’t see how changing the number of degrees of freedom would cause such an error. Are you sure that you did not redefine some variable before running again the code.
In general, consider posting the full code, even better a Minimal Working Example, which reproduces your error.

2 Likes
6

Sir, this is the code that I am using.

from dolfin import *
import numpy as np
import json
import matplotlib.pyplot as plt
import scipy as sp
%matplotlib inline
mesh = Mesh()
with XDMFFile("mesh/tetra.xdmf") as infile:
    infile.read(mesh)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
E, nu = (71e9), (0.33) #N/m**2
rho=2820 #Kg/m**3


mu = E/2./(1+nu)
lmbda = 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)
u_ = TrialFunction(V)
du = TestFunction(V)
support = CompiledSubDomain("x[1] == -0.020")
bc = DirichletBC(V, Constant((0.,0.,0.)), support)
k_form = inner(sigma(du),eps(u_))*dx

m_form = rho*dot(du,u_)*dx
K = PETScMatrix()
assemble(k_form, tensor=K)


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

bc.zero(M)
bc.apply(K)

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

N_eig = 40   # 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

# Extraction
eig_v = []
eig_l = []
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)

    # 3D eigenfrequency
    freq_3D = np.sqrt(r)/2/pi
    print("Mode {}:".format(i + 1))
    print("Solid FE: {0:8.5f} [Hz]".format(freq_3D))
    
    
    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector "+str(i))
    eigenmode.vector()[:] = rx
    file_results.write(eigenmode, 0)
    eig_l.append(freq_3D)
    eig_v.append(rx)

u = Function(V, name="Unit displacement")
u.interpolate(Constant((0, 0, 1)))
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: {:.2e}".format(qi))
    print("  Modal mass: {:.4f}".format(mi))
    print("  Effective mass: {:.2e}".format(meff_i))
    print("  Relative contribution: {:.2f} %".format(100 * meff_i / total_mass))

    combined_mass += meff_i
    
print(
    "\nTotal relative mass of the first {} modes: {:.2f} %".format(
        N_eig, 100 * combined_mass / total_mass
    )
)
7

Hi, xi should be a dofin Function not a PETScVector when using action. You should append eigenmode to eig_v, not rx.

1 Like