Assemble matrix in dolfinx

1

Hi, I am solving a modal analysis problem. I am not able assemble K and M matrices with current dolfin-x docker update. Has the syntax changed because it was working with previous docker version. The MWE is attached.

import sys
import dolfinx
import numpy as np
import slepc4py
from scipy.optimize import root
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

from dolfinx.mesh import CellType, locate_entities_boundary, create_box
from dolfinx.fem import (Constant, dirichletbc, Function, VectorFunctionSpace,
                         locate_dofs_topological)
from ufl import (TrialFunction, TestFunction, Identity, inner, dot, nabla_div,
                 nabla_grad, dx)

slepc4py.init(sys.argv)


L, B, H = 5., 0.3, 0.6

mesh = create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([L, H, B])],
    [25, 3, 2],
    cell_type=CellType.tetrahedron)

E, nu = (2E11), (0.3)
rho = (7850)
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)


def sigma(u):
    return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) \
           + 2*mu*epsilon(u)


V = VectorFunctionSpace(mesh, ("CG", 2))

u = TrialFunction(V)
v = TestFunction(V)


def clamped_boundary(x):
    return np.isclose(x[0], 0)


fdim = mesh.topology.dim - 1
boundary_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)

u_D = Function(V)
with u_D.vector.localForm() as loc:
    loc.set(0)

bc = dirichletbc(u_D, locate_dofs_topological(V, fdim, boundary_facets))

T = Constant(mesh, (PETSc.ScalarType(0), PETSc.ScalarType(0),
                    PETSc.ScalarType(0)))
k_form = inner(sigma(u), epsilon(v))*dx
m_form = rho*dot(u, v)*dx


K = dolfinx.fem.assemble_matrix(k_form, bc)
K.assemble()

M = dolfinx.fem.assemble_matrix(m_form, bcs=[bc])
M.assemble()

shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('sinvert')  # spectral transform
shift.setShift(0.0)  # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=50)
eigensolver.setProblemType(2)
eigensolver.setST(shift)
eigensolver.setType('krylovschur')  # Solver type
eigensolver.setWhichEigenpairs(7)  # Target Real
eigensolver.setConvergenceTest(2)
eigensolver.setOperators(K, M)

eps_type = eigensolver.getType()
print("Solution method: %s" % eps_type)
problem_type = eigensolver.getProblemType()
print("Problem type: %s" % problem_type)
which_eigenpairs = eigensolver.getWhichEigenpairs()
print("Eigenpairs type: %s" % which_eigenpairs)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.createVecs()


# Exact solution computation
def falpha(x):
    return np.cos(x) * np.cosh(x) + 1


def alpha(n):
    return root(falpha, (2*n + 1)*np.pi/2)['x'][0]


ef = 0
print("Number of converged eigenpairs %d" % evs)
if evs > 0:
    for i in range(evs):
        e_val = eigensolver.getEigenpair(i, vr, vi)

        if (~np.isclose(e_val.real, 1.0)):
            # Calculation of eigenfrequency from real part of eigenvalue
            freq_3D = np.sqrt(e_val.real)/2/np.pi

            # Beam eigenfrequency
            if ef % 2 == 0:
                # exact solution should correspond to weak axis bending
                I_bend = H*B**3/12.
            else:
                # exact solution should correspond to strong axis bending
                I_bend = B*H**3/12.

            freq_beam = alpha(ef/2)**2*np.sqrt(E*I_bend/(rho*B*H*L**4))/2/np.pi

            print(
                "Solid FE: {0:8.5f} [Hz] "
                "Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))

            ef += 1
2

See: Issue with dolfinx.fem.create_matrix - #2 by dokken

1 Like