How to apply many forces when solving elasticity problem

dolfinx
3

Hello @francesco-ballarin
Thanks for your quick reply.

For sure I did some forum research before directly asking, but unfortunately there were no help for me at all.

First of all I found this thread and tried to do something similar

class MyLinearProblem(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='int32')
            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 ValueError(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) -> _Function:
        """Solve the problem."""

        # Assemble lhs
        self._A.zeroEntries()
        dolfinx.fem.petsc.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)
            b_loc.setValues(self._point_force_dofs, self._point_force_values)

        dolfinx.fem.petsc.assemble_vector(self._b, self._L)

        # Apply boundary conditions to the rhs
        dolfinx.fem.petsc.apply_lifting(self._b, [self._a], bcs=[self.bcs])
        self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        dolfinx.fem.petsc.set_bc(self._b, self.bcs)

        # 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

The result I got I assume is good enough but there were some things which I can’t explain (I attached the picture see below)

Screenshot 2023-12-30 at 14.12.13
Screenshot 2023-12-30 at 14.12.131920×974 80.8 KB

Also I do not understand the part in code when I just assigned some scalar value into b, but how can I apply the force vector I defined? This is beyond my brain power it seems

My guess is you were pointing to this topic but I just was unable to figure out what to do with that information. I see an example of getting some cells or constructing some matrix but what should I do with that?
I’m sorry if I may look dumb to you at this point but I’m just doing my personal project with 0 knowledge related to this stuff (I’m just started to learn),so, any small hint of how to use it would be really appreciated.

In the end, to be more precise, I guess I do not interested in applying a force to a specific point I think it’s fine to me if I will be able to apply the force onto a cell triangle as well

But if I do that:

...
# selecting desired cell
p_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 2, point))
pt = dolfinx.mesh.meshtags(msh, 2, p_facets, 3)
...
# adding desired cell measure
ds = ufl.Measure('ds', domain=msh, subdomain_data=mt)
ds_p = ufl.Measure('ds', domain=msh, subdomain_data=pt) # subdomain is a desired tetra
...
# making equations
a = ufl.inner(sigma(u_tr), epsilon(u_test))*ufl.dx(domain=msh)
l = ufl.dot(b, u_test) * ufl.dx(domain=msh) + ufl.dot(f, u_test) * ds(1) + ufl.dot(fp, u_test) * ds_p(1) # I hope I just added a part then the force is applied onto top triangle of my cell

problem = LinearProblem(a, l, bcs=[bc_bot], petsc_options = petsc_options)
u = problem.solve()

I got the next error:

File ~/miniconda3/envs/dolfinx/lib/python3.11/site-packages/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
...
           # Assuming single domain
           # Check that subdomain data for each integral type is the same
           for data in sd.get(domain).values():
-->136     assert all([d is data[0] for d in data])
    138 mesh = domain.ufl_cargo()
    139 if mesh is None:

AssertionError:

I assume there is an explanation of that is somewhere but outside my head

To sum up, I didn’t succeed to find any definitive example of my problem. I’m positive to share mine if I managed it in the end.