Possible reuse of LU decomposition when solving linear elasticity problem in same mesh with diferent bcs

General
1

Hi,

I would like to optimize a linear elasticity code to run faster by reusing the LU decomposition used by PETSc.

Below I have a version of my code that runs using the function LinearProblem(), and another version where I assemble the Lu=b system.

  1. What I noticed is that to assemble A using assemble_matrix() you need the boundary conditions, why is this? Is L modified? Once I can move the assembly of L outside the loop, how could I reuse the LU decomposition from PETSc?
  2. In the function solve_fea() there is also a step where I fix the displacements along some directions. From what I remember this means the solution is also Lu=b, but with the rows and columns corresponding to the degrees of freedom are removed. I think this also means I have to remove those rows and columns from the decomposition, and then map them back the full problem. Is this correct? How could I do this in PETSc?

A self contained running code is shown below.

Thank you,

Alex


import numpy as np
import time

from mpi4py import MPI
from petsc4py import PETSc

import ufl
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, set_bc, LinearProblem


r_loading = 2.0e-3





def in_sphere(x, point, r_range):
    
    return ((x[0]-point[0])**2 + (x[1]-point[1])**2 + (x[2]-point[2])**2) < r_range**2


def mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal):
   
    dim = domain.topology.dim
    fdim = dim - 1
    domain.topology.create_connectivity(fdim, dim)

    def left_facets(xq):
        return in_sphere(xq, facet_node, r_loading)

    def right_facets(xq):
        return in_sphere(xq, grasp_node, r_loading)

    facets_left  = mesh.locate_entities(domain, fdim, left_facets)
    facets_right = mesh.locate_entities(domain, fdim, right_facets)

    facet_indices = np.concatenate([facets_left, facets_right])
    facet_markers = np.concatenate([np.full_like(facets_left,  1, dtype=np.int32),
                                    np.full_like(facets_right, 2, dtype=np.int32)])

    # Sort by facet index
    sort_perm = np.argsort(facet_indices)
    facet_indices = facet_indices[sort_perm]
    facet_markers = facet_markers[sort_perm]

    mt = mesh.meshtags(domain, fdim, facet_indices, facet_markers)
    return mt


def elasticity_parameters(E=10.0e4, nu=0.20):
    """Lamé parameters from E, nu."""
    lambda_ = E * nu  / ((1.0 + nu) * (1.0 - 2.0 * nu))
    mu      = E * 0.5 /  (1.0 + nu)
    
    return lambda_, mu


def epsilon(u):
    return ufl.sym(ufl.grad(u))


def sigma(u, lambda_, mu):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2.0 * mu * epsilon(u)


def build_primal_problem(domain, facet_tags, lambda_, mu, node_normal, lock_num,
                         p_traction=1.0/(np.pi*r_loading**2), body_force=(0.0, 0.0, -3941.0*9.8)):
    
    dim = domain.topology.dim
    fdim = dim - 1
    ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)

    V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    f_vec = fem.Constant(domain, default_scalar_type(body_force))

    # Traction direction: push along +x (change sign if needed)
    T = fem.Constant(domain, -p_traction * default_scalar_type((node_normal[0], node_normal[1], node_normal[2])))

    a = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx
    L = ufl.dot(f_vec, v) * ufl.dx + ufl.dot(T, v) * ds(2)

    # Dirichlet BC on left jaw (tag 1)
    left_facets = facet_tags.find(1)
    domain.topology.create_connectivity(fdim, dim)
    left_dofs = fem.locate_dofs_topological(V, fdim, left_facets)

    u_D = np.array((0.0, 0.0, 0.0), dtype=default_scalar_type)
    bc_left = fem.dirichletbc(u_D, left_dofs, V)


    # # Dirichlet BC to stop deflection normal to x (lock_num == 1) or y (lock_num == 2) axis   
    # if lock_num == 1:
    #     right_facets = facet_tags.find(2)
    #     right_dofs1   = fem.locate_dofs_topological(V.sub(1), fdim, right_facets)
    #     bc_1 = fem.dirichletbc(0.0, right_dofs1, V.sub(1))
    #     right_dofs2   = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
    #     bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
    #     return V, a, L, [bc_left, bc_1, bc_2]
    # if lock_num == 2:
    #     right_facets = facet_tags.find(2)
    #     right_dofs0   = fem.locate_dofs_topological(V.sub(0), fdim, right_facets)
    #     bc_1 = fem.dirichletbc(0.0, right_dofs0, V.sub(0))
    #     right_dofs2   = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
    #     bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
    #     return V, a, L, [bc_left, bc_1, bc_2]
    
    return V, a, L, [bc_left]




def array_dist_squared(x, p):
    
    return (x[:,0]-p[0])**2 + (x[:,1]-p[1])**2 + (x[:,2]-p[2])**2



def fea_solve(domain, facet_node, grasp_node, grasp_node_normal, lock_num):
 
    lambda_, mu = elasticity_parameters()
    facet_tags  = mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal)

    Vp, a_p, L_p, bcs_p = build_primal_problem(domain, facet_tags, lambda_, mu, grasp_node_normal, lock_num)
    
    problem = LinearProblem(a_p, L_p, bcs=bcs_p,
                            petsc_options_prefix="linear_problem",
                            petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
    uh = problem.solve()

       
    
    # Get displacement at node where load is applied
    dof_coords = Vp.tabulate_dof_coordinates()
    distances  = np.linalg.norm(dof_coords - grasp_node, axis=1)
    dof_index  = np.argmin(distances)
    u_node     = np.sqrt(uh.x.array[3*dof_index]**2 + uh.x.array[3*dof_index+1]**2 + uh.x.array[3*dof_index+2]**2)
    
    return domain, uh, u_node



     
# Part FEA
##################################################################
length = 20.0e-3
width  = 10.0e-3

grasp_facet_nodes = np.array([[0.0,    width/2, width/2],
                              [0.0,    width/2, width/2 + 1.0e-3],
                              [0.0,    width/2, width/2 + 2.0e-3]])
  
grasp_nodes2      = np.array([[length, width/2, width/2],
                              [length, width/2, width/2 + 1.0e-3],
                              [length, width/2, width/2 + 2.0e-3]])

grasp_node_normal = np.array([1.0,    0.0,     0.0    ])

cell_size = 0.5e-3

n_l = int(length/cell_size)
n_w = int(width /cell_size)

lock_num = 1


domain = mesh.create_box(MPI.COMM_WORLD, ((0.0, 0.0, 0.0), (length,width,width)), (n_l,n_w,n_w), mesh.CellType.tetrahedron)


for n in range(grasp_facet_nodes.shape[0]):    
    
    t_start = time.perf_counter()
    domain_f, uh, u_node = fea_solve(domain, grasp_facet_nodes[n,:], grasp_nodes2[n,:], grasp_node_normal, lock_num)
    t_end   = time.perf_counter()
    print(f"FEA elapsed time: {t_end  - t_start:.6f} seconds")


    with io.XDMFFile(domain_f.comm, f"max_p_stress_part_loading_{n+1}.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain_f)
        uh.name   = "u"
        xdmf.write_function(uh)
        
    print(u_node)


print('\n')





# Common assembly steps
t_start = time.perf_counter()

lambda_, mu = elasticity_parameters()

dim = domain.topology.dim
fdim = dim - 1

Vp = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
u  = ufl.TrialFunction(Vp)
v  = ufl.TestFunction(Vp)

body_force=(0.0, 0.0, -3941.0*9.8)
f_vec = fem.Constant(domain, default_scalar_type(body_force))

# Traction direction: push along +x (change sign if needed)
p_traction=1.0/(np.pi*r_loading**2)
T = fem.Constant(domain, -p_traction * default_scalar_type((grasp_node_normal[0], grasp_node_normal[1], grasp_node_normal[2])))

a_p = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx

domain.topology.create_connectivity(fdim, dim)

u_D   = np.array((0.0, 0.0, 0.0), dtype=default_scalar_type)


t_end   = time.perf_counter()
print(f"Common assembly time: {t_end  - t_start:.6f} seconds")





for n in range(grasp_facet_nodes.shape[0]):    
    
    t_start = time.perf_counter()
    
    
    facet_tags  = mark_boundaries(domain, grasp_facet_nodes[n,:],grasp_nodes2[n,:],grasp_node_normal)
    
    ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)

    L_p = ufl.dot(f_vec, v) * ufl.dx + ufl.dot(T, v) * ds(2)

    # Dirichlet BC on left jaw (tag 1)
    left_facets = facet_tags.find(1)
    left_dofs = fem.locate_dofs_topological(Vp, fdim, left_facets)

    bcs   = fem.dirichletbc(u_D, left_dofs, Vp)
    bcs_p = [bcs]
    
    
    
    """
    Solve linear system A u = b with PETSc LU.
    """
    a_form = fem.form(a_p)
    L_form = fem.form(L_p)

    A = assemble_matrix(a_form, bcs=bcs_p)
    A.assemble()

    b = assemble_vector(L_form)
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, bcs_p)

    uh = fem.Function(Vp)
    solver = PETSc.KSP().create(domain.comm)
    solver.setType(PETSc.KSP.Type.PREONLY)
    solver.getPC().setType(PETSc.PC.Type.LU)
    solver.getPC().setFactorSolverType('mumps')
    
    solver.setOperators(A)
    solver.solve(b, uh.x.petsc_vec)
    uh.x.scatter_forward()
    
    
    
    # Get displacement at node where load is applied
    dof_coords = Vp.tabulate_dof_coordinates()
    distances  = np.linalg.norm(dof_coords - grasp_nodes2[n,:], axis=1)
    dof_index  = np.argmin(distances)
    u_node     = np.sqrt(uh.x.array[3*dof_index]**2 + uh.x.array[3*dof_index+1]**2 + uh.x.array[3*dof_index+2]**2)
    
    
    
    
    
    
   
    t_end   = time.perf_counter()
    print(f"FEA elapsed time: {t_end  - t_start:.6f} seconds")


    with io.XDMFFile(domain_f.comm, f"max_p_stress2_part_loading_{n+1}.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain_f)
        uh.name   = "u"
        xdmf.write_function(uh)


    print(u_node)
2

Here’s a shorter answer without looking too much at your code.

The reason the LU-factorization is not reused in the LinearProblem-solves is that the system matrix is reassembled for each problem.solve()-call, since constants might have been changed etc. Since the matrix is reassembled each time, the LU-factorization is discarded.

You can get around this by subclassing the dolfinx LinearProblem and copy-pasting some code like below. The only change I’ve made is to move the matrix assembly to a different method and comment it out in the solve() method. You then call assemble_matrix() before the first solve and any time the system matrix should be changed.

The reason assemble_matrix() needs the boundary conditions is that the rows corresponding to the boundary conditions are changed to be 1 on the diagonal and 0 elsewhere and then the same is done to the columns corresponding to boundary conditions through lifting. This process is explained very well in this tutorial by @dokken. So, the actual values have no effect on the assemble_matrix()-step, only which dofs are included. Therefore you can do this step as well as defining your solver outside your loop, as long as you do the lifting at each step.

from petsc4py import PETSc

import dolfinx
import dolfinx.fem.petsc

from dolfinx.fem.function import Function as _Function
from typing import Sequence

from dolfinx.fem.bcs import bcs_by_block as _bcs_by_block
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting
from dolfinx.fem.forms import extract_function_spaces as _extract_function_spaces

class MyLinearProblem(dolfinx.fem.petsc.LinearProblem):
    """
    A customized version of :class:`dolfinx.fem.petsc.LinearProblem` that
    exposes the `assemble_matrix` method for assembling the system matrix
    separately from the `solve` method. Saves time when solving multiple
    linear systems with the same system matrix but different right-hand sides.
    """
            
    def assemble_matrix(self) -> None:
        """Assemble the system matrix and preconditioner matrix (if any)."""
        # Assemble lhs
        self.A.zeroEntries()
        assemble_matrix(self.A, self.a, bcs=self.bcs)  # type: ignore[arg-type, misc]
        self.A.assemble()

        # Assemble preconditioner
        if self.P_mat is not None:
            self.P_mat.zeroEntries()
            assemble_matrix(self.P_mat, self.preconditioner, bcs=self.bcs)  # type: ignore[arg-type, misc]
            self.P_mat.assemble()


    def solve(self) -> _Function | Sequence[_Function]:
        """Solve the problem.

        This method updates the solution ``u`` function(s) stored in the
        problem instance.

        Note:
            The user is responsible for asserting convergence of the KSP
            solver e.g. ``problem.solver.getConvergedReason() > 0``.
            Alternatively, pass ``"ksp_error_if_not_converged" : True`` in
            ``petsc_options`` to raise a ``PETScError`` on failure.

        Returns:
            The solution function(s).
        """

        # # Assemble lhs
        # self.A.zeroEntries()
        # assemble_matrix(self.A, self.a, bcs=self.bcs)  # type: ignore[arg-type, misc]
        # self.A.assemble()

        # # Assemble preconditioner
        # if self.P_mat is not None:
        #     self.P_mat.zeroEntries()
        #     assemble_matrix(self.P_mat, self.preconditioner, bcs=self.bcs)  # type: ignore[arg-type, misc]
        #     self.P_mat.assemble()

        # Assemble rhs
        dolfinx.la.petsc._zero_vector(self.b)
        assemble_vector(self.b, self.L)  # type: ignore[arg-type]

        # Apply boundary conditions to the rhs
        if self.bcs is not None:
            if isinstance(self.u, Sequence):  # block or nest
                bcs1 = _bcs_by_block(_extract_function_spaces(self.a, 1), self.bcs)  # type: ignore[arg-type]
                apply_lifting(self.b, self.a, bcs=bcs1)  # type: ignore[arg-type]
                dolfinx.la.petsc._ghost_update(
                    self.b,
                    PETSc.InsertMode.ADD,  # type: ignore[attr-defined]
                    PETSc.ScatterMode.REVERSE,  # type: ignore[attr-defined]
                )
                bcs0 = _bcs_by_block(_extract_function_spaces(self.L), self.bcs)  # type: ignore[arg-type]
                dolfinx.fem.petsc.set_bc(self.b, bcs0)
            else:  # single form
                apply_lifting(self.b, [self.a], bcs=[self.bcs])  # type: ignore[arg-type]
                dolfinx.la.petsc._ghost_update(
                    self.b,
                    PETSc.InsertMode.ADD,  # type: ignore[attr-defined]
                    PETSc.ScatterMode.REVERSE,  # type: ignore[attr-defined]
                )
                for bc in self.bcs:
                    bc.set(self.b.array_w)
        else:
            dolfinx.la.petsc._ghost_update(self.b, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE)  # type: ignore[attr-defined]

        # Solve linear system and update ghost values in the solution
        self.solver.solve(self.b, self.x)
        dolfinx.la.petsc._ghost_update(self.x, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD)  # type: ignore[attr-defined]
        dolfinx.fem.petsc.assign(self.x, self.u)
        return self.u
3

Hi @ottarhellan thank you! That is what I was looking for. This answers both questions as well, even in the bc case for question 1 I am still fixing some dofs so L will still be modified. What I don’t see in your code is how does PETSc know when to reuse L? Does L need to also be passed to this new class?

4

You use the class in the exact same way as the dolfinx LinearProblem class, except you need to call assemble_matrix() before the first solve and any time something has changed that will make the system matrix be different. So you pass it your bilinear form a_p, linear form L_p and boundary conditions bcs_p, just like you do in your fea_solve function, as well as the petsc options like using LU factorization direct solvers.

I guess PETSc knows to reuse the factorization of L since it is never changed unless you call assemble_matrix() again.

5

Hi @ottarhellan , I tried implementing the subclass and its running faster but not giving the right nodal displacements (they should match the ones from the first loop). Do you see anything wrong in the code below?

import numpy as np
import time

from mpi4py import MPI
from petsc4py import PETSc

import ufl
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem

import dolfinx.fem.petsc

from dolfinx.fem.function import Function as _Function
from typing import Sequence

from dolfinx.fem.bcs import bcs_by_block as _bcs_by_block
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting
from dolfinx.fem.forms import extract_function_spaces as _extract_function_spaces



r_loading = 2.0e-3



def in_sphere(x, point, r_range):
    
    return ((x[0]-point[0])**2 + (x[1]-point[1])**2 + (x[2]-point[2])**2) < r_range**2


def mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal):
   
    dim = domain.topology.dim
    fdim = dim - 1
    domain.topology.create_connectivity(fdim, dim)

    def left_facets(xq):
        return in_sphere(xq, facet_node, r_loading)

    def right_facets(xq):
        return in_sphere(xq, grasp_node, r_loading)

    facets_left  = mesh.locate_entities(domain, fdim, left_facets)
    facets_right = mesh.locate_entities(domain, fdim, right_facets)

    facet_indices = np.concatenate([facets_left, facets_right])
    facet_markers = np.concatenate([np.full_like(facets_left,  1, dtype=np.int32),
                                    np.full_like(facets_right, 2, dtype=np.int32)])

    # Sort by facet index
    sort_perm = np.argsort(facet_indices)
    facet_indices = facet_indices[sort_perm]
    facet_markers = facet_markers[sort_perm]

    mt = mesh.meshtags(domain, fdim, facet_indices, facet_markers)
    return mt


def elasticity_parameters(E=10.0e4, nu=0.20):
    """Lamé parameters from E, nu."""
    lambda_ = E * nu  / ((1.0 + nu) * (1.0 - 2.0 * nu))
    mu      = E * 0.5 /  (1.0 + nu)
    
    return lambda_, mu


def epsilon(u):
    return ufl.sym(ufl.grad(u))


def sigma(u, lambda_, mu):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2.0 * mu * epsilon(u)


def build_primal_problem(domain, facet_tags, lambda_, mu, node_normal, lock_num,
                         p_traction=1.0/(np.pi*r_loading**2), body_force=(0.0, 0.0, -3941.0*9.8)):
    
    dim = domain.topology.dim
    fdim = dim - 1
    ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)

    V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    f_vec = fem.Constant(domain, default_scalar_type(body_force))

    # Traction direction: push along +x (change sign if needed)
    T = fem.Constant(domain, -p_traction * default_scalar_type((node_normal[0], node_normal[1], node_normal[2])))

    a = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx
    L = ufl.dot(f_vec, v) * ufl.dx + ufl.dot(T, v) * ds(2)

    # Dirichlet BC on left jaw (tag 1)
    left_facets = facet_tags.find(1)
    domain.topology.create_connectivity(fdim, dim)
    left_dofs = fem.locate_dofs_topological(V, fdim, left_facets)

    u_D = np.array((0.0, 0.0, 0.0), dtype=default_scalar_type)
    bc_left = fem.dirichletbc(u_D, left_dofs, V)


    # # Dirichlet BC to stop deflection normal to x (lock_num == 1) or y (lock_num == 2) axis   
    # if lock_num == 1:
    #     right_facets = facet_tags.find(2)
    #     right_dofs1   = fem.locate_dofs_topological(V.sub(1), fdim, right_facets)
    #     bc_1 = fem.dirichletbc(0.0, right_dofs1, V.sub(1))
    #     right_dofs2   = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
    #     bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
    #     return V, a, L, [bc_left, bc_1, bc_2]
    # if lock_num == 2:
    #     right_facets = facet_tags.find(2)
    #     right_dofs0   = fem.locate_dofs_topological(V.sub(0), fdim, right_facets)
    #     bc_1 = fem.dirichletbc(0.0, right_dofs0, V.sub(0))
    #     right_dofs2   = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
    #     bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
    #     return V, a, L, [bc_left, bc_1, bc_2]
    
    return V, a, L, [bc_left]



def array_dist_squared(x, p):
    
    return (x[:,0]-p[0])**2 + (x[:,1]-p[1])**2 + (x[:,2]-p[2])**2



def fea_solve(domain, facet_node, grasp_node, grasp_node_normal, lock_num):
 
    lambda_, mu = elasticity_parameters()
    facet_tags  = mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal)

    Vp, a_p, L_p, bcs_p = build_primal_problem(domain, facet_tags, lambda_, mu, grasp_node_normal, lock_num)
    
    problem = LinearProblem(a_p, L_p, bcs=bcs_p,
                            petsc_options_prefix="linear_problem",
                            petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
    uh = problem.solve()

       
    
    # Get displacement at node where load is applied
    dof_coords = Vp.tabulate_dof_coordinates()
    distances  = np.linalg.norm(dof_coords - grasp_node, axis=1)
    dof_index  = np.argmin(distances)
    u_node     = np.sqrt(uh.x.array[3*dof_index]**2 + uh.x.array[3*dof_index+1]**2 + uh.x.array[3*dof_index+2]**2)
    
    return domain, uh, u_node



class MyLinearProblem(dolfinx.fem.petsc.LinearProblem):
    """
    A customized version of :class:`dolfinx.fem.petsc.LinearProblem` that
    exposes the `assemble_matrix` method for assembling the system matrix
    separately from the `solve` method. Saves time when solving multiple
    linear systems with the same system matrix but different right-hand sides.
    """
            
    def assemble_matrix(self) -> None:
        """Assemble the system matrix and preconditioner matrix (if any)."""
        # Assemble lhs
        self.A.zeroEntries()
        assemble_matrix(self.A, self.a, bcs=self.bcs)  # type: ignore[arg-type, misc]
        self.A.assemble()

        # Assemble preconditioner
        if self.P_mat is not None:
            self.P_mat.zeroEntries()
            assemble_matrix(self.P_mat, self.preconditioner, bcs=self.bcs)  # type: ignore[arg-type, misc]
            self.P_mat.assemble()


    def solve(self) -> _Function | Sequence[_Function]:
        """Solve the problem.

        This method updates the solution ``u`` function(s) stored in the
        problem instance.

        Note:
            The user is responsible for asserting convergence of the KSP
            solver e.g. ``problem.solver.getConvergedReason() > 0``.
            Alternatively, pass ``"ksp_error_if_not_converged" : True`` in
            ``petsc_options`` to raise a ``PETScError`` on failure.

        Returns:
            The solution function(s).
        """

        # # Assemble lhs
        # self.A.zeroEntries()
        # assemble_matrix(self.A, self.a, bcs=self.bcs)  # type: ignore[arg-type, misc]
        # self.A.assemble()

        # # Assemble preconditioner
        # if self.P_mat is not None:
        #     self.P_mat.zeroEntries()
        #     assemble_matrix(self.P_mat, self.preconditioner, bcs=self.bcs)  # type: ignore[arg-type, misc]
        #     self.P_mat.assemble()

        # Assemble rhs
        dolfinx.la.petsc._zero_vector(self.b)
        assemble_vector(self.b, self.L)  # type: ignore[arg-type]

        # Apply boundary conditions to the rhs
        if self.bcs is not None:
            if isinstance(self.u, Sequence):  # block or nest
                bcs1 = _bcs_by_block(_extract_function_spaces(self.a, 1), self.bcs)  # type: ignore[arg-type]
                apply_lifting(self.b, self.a, bcs=bcs1)  # type: ignore[arg-type]
                dolfinx.la.petsc._ghost_update(
                    self.b,
                    PETSc.InsertMode.ADD,  # type: ignore[attr-defined]
                    PETSc.ScatterMode.REVERSE,  # type: ignore[attr-defined]
                )
                bcs0 = _bcs_by_block(_extract_function_spaces(self.L), self.bcs)  # type: ignore[arg-type]
                dolfinx.fem.petsc.set_bc(self.b, bcs0)
            else:  # single form
                apply_lifting(self.b, [self.a], bcs=[self.bcs])  # type: ignore[arg-type]
                dolfinx.la.petsc._ghost_update(
                    self.b,
                    PETSc.InsertMode.ADD,  # type: ignore[attr-defined]
                    PETSc.ScatterMode.REVERSE,  # type: ignore[attr-defined]
                )
                for bc in self.bcs:
                    bc.set(self.b.array_w)
        else:
            dolfinx.la.petsc._ghost_update(self.b, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE)  # type: ignore[attr-defined]

        # Solve linear system and update ghost values in the solution
        self.solver.solve(self.b, self.x)
        dolfinx.la.petsc._ghost_update(self.x, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD)  # type: ignore[attr-defined]
        dolfinx.fem.petsc.assign(self.x, self.u)
        return self.u




     
# Part FEA
##################################################################
length = 20.0e-3
width  = 10.0e-3

grasp_facet_nodes = np.array([[0.0,    width/2, width/2],
                              [0.0,    width/2, width/2 + 1.0e-3],
                              [0.0,    width/2, width/2 + 2.0e-3]])
  
grasp_nodes2      = np.array([[length, width/2, width/2],
                              [length, width/2, width/2 + 1.0e-3],
                              [length, width/2, width/2 + 2.0e-3]])

grasp_node_normal = np.array([1.0,    0.0,     0.0    ])

cell_size = 0.5e-3

n_l = int(length/cell_size)
n_w = int(width /cell_size)

lock_num = 1


domain = mesh.create_box(MPI.COMM_WORLD, ((0.0, 0.0, 0.0), (length,width,width)), (n_l,n_w,n_w), mesh.CellType.tetrahedron)


for n in range(grasp_facet_nodes.shape[0]):    
    
    t_start = time.perf_counter()
    domain_f, uh, u_node = fea_solve(domain, grasp_facet_nodes[n,:], grasp_nodes2[n,:], grasp_node_normal, lock_num)
    t_end   = time.perf_counter()
    print(f"FEA elapsed time: {t_end  - t_start:.6f} seconds")


    with io.XDMFFile(domain_f.comm, f"max_p_stress_part_loading_{n+1}.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain_f)
        uh.name   = "u"
        xdmf.write_function(uh)
        
    print(u_node)


print('\n')



# Common assembly steps
t_start = time.perf_counter()

lambda_, mu = elasticity_parameters()

dim = domain.topology.dim
fdim = dim - 1

Vp = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
u  = ufl.TrialFunction(Vp)
v  = ufl.TestFunction(Vp)

body_force=(0.0, 0.0, -3941.0*9.8)
f_vec = fem.Constant(domain, default_scalar_type(body_force))

# Traction direction: push along +x (change sign if needed)
p_traction=1.0/(np.pi*r_loading**2)
T = fem.Constant(domain, -p_traction * default_scalar_type((grasp_node_normal[0], grasp_node_normal[1], grasp_node_normal[2])))

a_p = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx

domain.topology.create_connectivity(fdim, dim)

u_D   = np.array((0.0, 0.0, 0.0), dtype=default_scalar_type)


t_end   = time.perf_counter()
print(f"Common assembly time: {t_end  - t_start:.6f} seconds")



for n in range(grasp_facet_nodes.shape[0]):    
    
    t_start = time.perf_counter()
    
    
    facet_tags  = mark_boundaries(domain, grasp_facet_nodes[n,:],grasp_nodes2[n,:],grasp_node_normal)
    
    ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)

    L_p = ufl.dot(f_vec, v) * ufl.dx + ufl.dot(T, v) * ds(2)

    # Dirichlet BC on left jaw (tag 1)
    left_facets = facet_tags.find(1)
    left_dofs = fem.locate_dofs_topological(Vp, fdim, left_facets)

    bcs   = fem.dirichletbc(u_D, left_dofs, Vp)
    bcs_p = [bcs]
    
    if n == 0:
        problem2 = MyLinearProblem(a_p, L_p,
                                   petsc_options_prefix="linear_problem",
                                   petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
        problem2.assemble_matrix()

        
    uh = problem2.solve()
    
    
    
    # Get displacement at node where load is applied
    dof_coords = Vp.tabulate_dof_coordinates()
    distances  = np.linalg.norm(dof_coords - grasp_nodes2[n,:], axis=1)
    dof_index  = np.argmin(distances)
    u_node     = np.sqrt(uh.x.array[3*dof_index]**2 + uh.x.array[3*dof_index+1]**2 + uh.x.array[3*dof_index+2]**2)
    
    
    
    t_end   = time.perf_counter()
    print(f"FEA elapsed time: {t_end  - t_start:.6f} seconds")


    with io.XDMFFile(domain_f.comm, f"max_p_stress2_part_loading_{n+1}.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain_f)
        uh.name   = "u"
        xdmf.write_function(uh)


    print(u_node)
6

Hi Alex,

You need to give the boundary conditions to the MyLinearProblem class, otherwise it won’t know to apply the boundary conditions. For you it would be problem2 = MyLinearProblem(a_p, L_p, bcs=bcs_p, .... Each time you call problem2.solve(), it will use the data in bcs to assemble the boundary condition contributions to the right hand side vector using the current data in u_D.

However, I see now that the degrees of freedom that are fixed using the boundary condition is changing each time you are solving the problem, and in this case reusing the LU-factorization won’t work, since the system matrix won’t be the same between iterations. I’m sorry I did not notice this earlier. I think in this case, creating a new standard dolfinx.fem.petsc.LinearProblem for each solve is pretty much as fast as you can expect, I don’t believe the overhead is very large.

7

Hi @ottarhellan , isnt that solved by lifting like you posted initially? From the point of view of LU decomposition this should be possible, since the decomposition of the unconstrained matrix could be used and then the matrices can be modified by not considering the rows and columns coresponding to the zero displacement dofs.

8

I’ve been talking to gemini a bit and it seems like the LU-decomposition reuse you’re thinking about is impossible, without me really understanding all the details. In any case, the matrix that is passed to the LU-decomposition routine is the assembled petsc matrix, which has rows and columns corresponding to bcs zeroed out outside the diagonal already. In addition, gemini tells me it is impossible to get the explicit L and U factors from mumps.

If you really want to reuse LU factorizations here, I think it will take a lot more effort than making some small modifications to existing dolfinx code.

Best of luck with your work.

9

Hi @ottarhellan, I looked into this further and the method to solve faster with multiple zero displacement BCs and the same mesh may be static condensation:

This is similar to Guyan reduction:

From what I read if the stiffness matrix (K in Ku=f) is pivoted to have the boundary displacements on the lower right of a block matrix like this:

\begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix} u_i \\ u_b \end{pmatrix} = \begin{pmatrix} f_i \\ f_b \end{pmatrix}

Where the i subscript is for internal nodes and b is for boundary nodes (or parts of the force vector that can change?). Then the boundary displacements can be found using:

(D-CA^{-1}B)u_b=f_b-CA^{-1} f_i

I think f_i is zero (no applied forces on internal nodes). The term on the left multiplying u_b is the Schur complement:

This should be faster since a smaller system is solved, and A^{-1} does not change since its related to internal nodes that are not subject to different displacement boundary conditions.

This may take some time to implement. There are some points to work out like:

  1. What is an efficient way to get the inverse of matrix A in fenicsx?
  2. What nodes are affected when using zero displacement BCs, traction BCs, and what extra nodes will need to be included to calculate the stress of the body? This will be needed as well. Is that all the nodes of surface elements, not just surface nodes?
  3. From the L form in fenicsx, I think the body forces stay constant in A, how are they applied here?
10

I don’t really have the experience to help much with this. Regarding point 1 however, I can tell you that you should generally handle A^{-1} implicitly. If you set up a solver for A, you can get the action of A^{-1}, i.e. computing x = A^{-1} b by solving Ax = b using your solver, and you can then use that to solve the schur complement system efficiently. You don’t ever actually need A^{-1} and you can probably not hold it in memory anyway, you just need its action on vectors. As I understand it, the idea is to solve

\begin{pmatrix} A & C^T \\ C & 0 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} f \\ u_D \end{pmatrix} ,

where Ax = f is the unconstrained pde (in this case with the permanently fixed left side bcs applied, so A is invertible) and Cx = u_D applies the added Dirichlet boundary conditions and has much fewer rows than A. This will involve solving sub equations involving A and a smaller dense matrix you build using the action of A^{-1} and C.

I’m not convinced this is worth the effort though and would suggest considering whether just making a new LinearProblem for each new solve is fast enough for your needs.

11

Thanks for the information. I agree I’ll use the LinearSolve() for my current work. It’s taking 4 to 20 seconds for a converged mesh at the moment.