How to apply point forces by modifying the assembled matrix for static linear elasticity

dolfinx
1

Hi Fenicsx community,

I’m trying to modify the assembled matrix from the linear part in order to apply a point force. I saw from here : Dirac delta distribution (dolfinx) that the Dirac delta distribution is not yet implemented in fenicsx, unless I misunderstood.

But like said in https://fenicsproject.org/qa/14258/how-to-apply-pointsource-to-linearvariationalproblem/, there is another option, add the force on assembled vectors.

So far, I’ve tried to do the second solution, it works, but only if I solve the problem with numpy, by transferring the dense matrix to python, but it is not an efficient method. Here is what I did :

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import dx, ds, Measure, TrialFunction, TestFunction, sym, grad, div, Identity, inner, dot
from dolfinx import mesh, fem, plot, io
import numpy as np

domain = mesh.create_rectangle(MPI.COMM_SELF, ((-2, 0), (2, 1)), (50,30), mesh.CellType.quadrilateral)
V2 = fem.VectorFunctionSpace(domain, ("CG", 2))

E = fem.Constant(domain, ScalarType(1))
nu = fem.Constant(domain, ScalarType(0.3))
lamda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / (2 * (1 + nu))
T = fem.Constant(domain, ScalarType((0, 0)))

boundary_dofs = fem.locate_dofs_topological(V2, f_dim, lambda x: np.isclose(x[0], -2))
bc = fem.dirichletbc(np.zeros([2], dtype=ScalarType), boundary_dofs, V2)

dofs = fem.locate_dofs_geometrical(V2, lambda x: np.isclose(x[0], 2) & np.isclose(x[1], 0.5))

u = TrialFunction(V2)
v = TestFunction(V2)

def eps(u):
    return sym(grad(u))

def sigma(u):
    return (lamda * div(u) * Identity(2) + 2 * mu * eps(u))

a = inner(sigma(u), eps(v)) * dx
L = dot(T, v) * ds

problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.b.setValue(dofs, -1)  # Assign the force on the good dof

uh = problem.solve()
print((uh.x.array == 0).all())  # return True

print(problem.b.getArray()[dofs])  # After solving the problem it is = 0 instead of -1 :(

Thanks for your help

1 Like
2

solve calls assemble vector which would overwrite any entries in vector:

You would need to use the PETSc KSP interface directly, or make a modified version of the solve command, that applies this entry to the vector after assembling and applying lifting.

1 Like
3

Thank you very much @dokken it works perfectly.

Here the code that I use :

class MyLinearProblem(fem.petsc.LinearProblem):
    def __init__(self, *args, **kwargs):
        if 'point_force_dofs' in kwargs and 'point_force_values' in kwargs:
            point_force_dofs = kwargs['point_force_dofs']
            point_force_values = kwargs['point_force_values']
            if not isinstance(type(point_force_dofs), np.ndarray):
                point_force_dofs = np.array([point_force_dofs], dtype='int64')
            if not isinstance(type(point_force_values), np.ndarray):
                point_force_values = np.array([point_force_values])
            point_force_dofs = point_force_dofs.flatten()
            point_force_values = point_force_values.flatten()

            if len(point_force_dofs) != len(point_force_values):
                raise f"The length of point_force_dofs {len(point_force_values)} is not equal to the length of " \
                      f"point_force_values"

            del kwargs['point_force_dofs'], kwargs['point_force_values']
        else:
            point_force_dofs = None
            point_force_values = None
        self._point_force_dofs = point_force_dofs
        self._point_force_values = point_force_values

        super().__init__(*args, **kwargs)

    def solve(self):
        # Assemble lhs
        self._A.zeroEntries()
        from dolfinx.fem.petsc import _assemble_matrix_mat
        _assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
        self._A.assemble()

        # Assemble rhs
        with self._b.localForm() as b_loc:
            b_loc.set(0)
        from dolfinx.fem.petsc import assemble_vector
        assemble_vector(self._b, self._L)

        # Apply boundary conditions to the rhs
        from dolfinx.fem.petsc import apply_lifting
        apply_lifting(self._b, [self._a], bcs=[self.bcs])
        self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        from dolfinx.fem.petsc import set_bc
        set_bc(self._b, self.bcs)
        
        # Add point loads
        if self._point_force_values is not None and self._point_force_dofs is not None:
            self._b.setValue(self._point_force_dofs, self._point_force_values)

        # Solve linear system and update ghost values in the solution
        self._solver.solve(self._b, self._x)
        self.u.x.scatter_forward()

        return self.u

Have a nice evening !

1 Like