How to apply many forces when solving elasticity problem

dolfinx
1

Hello everyone,

I’m trying to solve linear elasticity problem. I managed to solve a problem when a box is under pressure from top, but how I also want to add another force which is applied to a specific area (point) to the left-top edge of the box

Here I have 3 conditions of finding bot facets, top facets and point facets

V = dolfinx.fem.VectorFunctionSpace(msh, ("CG", 1)) # displacement function space

# --------- Boundary Conditions --------- #
def bot(x):
    return np.isclose(x[2], 0.0, atol = 0.0006) # def of which points will be included in bc at bottom surface

def top(x):
    return np.isclose(x[2], height, atol=0.001)

def point(x):
    return np.logical_and(np.logical_and(np.isclose(x[0], length, atol = 0.0001), np.isclose(x[1], width / 2 + 0.0003, atol = 0.0006)), np.isclose(x[2], height, atol = 0.0006)) # def on which points additional force will be applied

bot_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 0, bot))
bdofs = dolfinx.fem.locate_dofs_topological(V, 0, bot_facets)

top_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 2, top))
mt = dolfinx.mesh.meshtags(msh, 2, top_facets, 1)

p_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 2, point))
pdofs = dolfinx.fem.locate_dofs_topological(V, 0, p_facets)
pt = dolfinx.mesh.meshtags(msh, 0, top_facets, 3)


bc_bot = dolfinx.fem.dirichletbc(dolfinx.fem.function.Constant(msh, (0.0, 0.0, 0.0)), bdofs, V)     # fixed bottom

Here is my problem definition

from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem.function import Function as _Function
from petsc4py import PETSc


# Body force = mass
grav = 9.8# [N/kg] gravitational force
b = dolfinx.fem.function.Constant(msh, (0.0, 0.0, -rho*grav))

f_mag = 3e5
# pressure
f = dolfinx.fem.function.Constant(msh,(0.0, 0.0, -f_mag))
# additional force
fp = dolfinx.fem.function.Constant(msh,(0.0, 0.0, -1e4))

# --------- Vector Function Space Definition --------- #
# a function space is created from mesh for the displacements to be canculated

u_tr = ufl.TrialFunction(V)
u_test = ufl.TestFunction(V)

ds = ufl.Measure('ds', domain=msh, subdomain_data=mt)

# --------- Stress and Strain Equations --------- #
# Strain function
def epsilon(u):
    return ufl.sym(ufl.grad(u))

# Stress function
def sigma(u):
    return lambda_*ufl.div(u)*ufl.Identity(3) + 2*mu*epsilon(u)

# --------- Weak Formulation of the Problem --------- #
a = ufl.inner(sigma(u_tr), epsilon(u_test))*ufl.dx(domain=msh)

# !!! how to add fp = dolfinx.fem.function.Constant(msh,(0.0, 0.0, -1e4)) applied on point here
l = ufl.dot(b, u_test) * ufl.dx(domain=msh) + ufl.dot(f, u_test) * ds(1)
# --------- Solver --------- #
petsc_options={}
problem = LinearProblem(a, l, bcs=[bc_bot], petsc_options = petsc_options)
u = problem.solve()

so I added measures dx of entire box and apply self box mass, ds for top facets and apply pressure but have how idea how to add desired point force

I tried “dP” meause but ended up with Attrubute error “vertex”

Thanks in advance

2

Please make a search with a keyword like “point source”. There were a couple of posts earlier this month about point loads, so you may want to sort by most recent post. Once you have found the post, make sure to link it here.

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.

4

I would suggest having a look at: How to construct a sparse matrix used for calculating function values at multiple points - #4 by dokken
and subsequent posts.

For your problem, you need to insert the point source contributions in your rhs vector.
For your other questions, it is not easy to help you, as you haven’t provided a minimal reproducible example.
My best guess is that you are using a vector element, and that the point source has not taken the block size of the dofmap into account.

5

Hi, @dokken

Thanks for the reply, will look at it.

Regarding providing an example… Will the colab notebook1
notebook2 links
work?

It won’t work directly in it unfortunately but should work on default dolfinx env (conda, docker)

6

The code is quite far from a minimal example.
A minimal example consists of:

  • The minimal amount of code to reproduce your error

This means that all variables, attributes,classes not required to reproduce the error should be removed.

The reason for this is that it is otherwise quite a lot of work that goes into understanding all the things done in your code.
See for instance:
Read before posting: How do I get my question answered? - #4 by nate for examples of MWEs.

7

Hope I cropped enough those two error when I tried to add another ds measure

Extended LinearProblem class and results are concerning

Not “that” minimial as in your examples, but only necessary code to reproduce my problems I believe

8

Could you instead of gmsh use a built in mesh (I.e. dolfinx.mesh.create_box) and instead of adding links to external storage, post your code directly on the forum using
3x` encapsulation, i.e.

```python
# Add code here

```

For each of the codes, please highlight if:

  1. There is an error, if so add the full traceback with:
    ```bash
Add error traceback here
```
  2. If something is unexpected in visualization, please add plots of the visualization.
9

Understood.

This is the code version in which I tried to add the 2nd force applied onto another boundary (despite it’s called “point” it’s more like an area)

# Example when errorr occurred when I added another ds

import dolfinx
import ufl
import numpy as np
import matplotlib.pyplot as plt
import dolfinx.fem.petsc as ptsc
from mpi4py import MPI

from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem.function import Function as _Function
from petsc4py import PETSc    


# Create a cylinder part with radius 20mm, height: 20mm

height = 0.100 #dimensions in m
length = 0.040 #dimensions in m
width = 0.020 #dimensions in m
 

# create the mesh from the points of the tetrahedral cells 
msh = dolfinx.mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]), np.array([length, width, height])], [100, 100, 100])

# --------- Material Parameters --------- # 

# Density (kg/m3)
rho = 2698.9

# Young's modulus (Pa) and Poisson's ratio
E = 68e9 
nu = 0.30

# Lame's constants
lambda_ = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

V = dolfinx.fem.VectorFunctionSpace(msh, ("CG", 1)) # displacement function space

# --------- Boundary Conditions --------- #
def bot(x):
    return np.isclose(x[2], 0.0, atol = 0.0006) # def of which points will be included in bc at bottom surface

def top(x):
    return np.isclose(x[2], height, atol=0.001)

def point(x):
    return np.logical_and(np.logical_and(np.isclose(x[0], length, atol = 0.001), np.isclose(x[1], width / 2, atol = 0.001)), np.isclose(x[2], height, atol = 0.001)) # def of which points will be included in bc at top surface


# Define surface for traction boundary condition
# Definition of Neumann boundary condition domain

top_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 2, top))
mt = dolfinx.mesh.meshtags(msh, 2, top_facets, 1)

p_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 2, point))
pdofs = dolfinx.fem.locate_dofs_topological(V, 2, p_facets)
pt = dolfinx.mesh.meshtags(msh, 2, p_facets, 3)

bot_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 0, bot))
bdofs = dolfinx.fem.locate_dofs_topological(V, 0, bot_facets)
bc_bot = dolfinx.fem.dirichletbc(dolfinx.fem.function.Constant(msh, (0.0, 0.0, 0.0)), bdofs, V)     # fixed bottom


# Body force = mass
grav = 9.8# [N/kg] gravitational force
b = dolfinx.fem.function.Constant(msh, (0.0, 0.0, -rho*grav))

f_mag = 3e5

f = dolfinx.fem.function.Constant(msh,(0.0, 0.0, -f_mag))
fp = dolfinx.fem.function.Constant(msh,(0.0, 0.0, -1e4))

bc_p = dolfinx.fem.dirichletbc(fp, pdofs, V)    

# --------- Vector Function Space Definition --------- #
# a function space is created from mesh for the displacements to be canculated

u_tr = ufl.TrialFunction(V)
u_test = ufl.TestFunction(V)

ds = ufl.Measure('ds', domain=msh, subdomain_data=mt)
ds_p = ufl.Measure('ds', domain=msh, subdomain_data=pt) # subdomain is a desired tetra

# --------- Stress and Strain Equations --------- #
# Strain function
def epsilon(u):
    return ufl.sym(ufl.grad(u))

# Stress function
def sigma(u):
    return lambda_*ufl.div(u)*ufl.Identity(3) + 2*mu*epsilon(u)

# --------- Weak Formulation of the Problem --------- #
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)
# --------- Solver --------- #
petsc_options={}
print(pdofs)
problem = LinearProblem(a, l, bcs=[bc_bot], petsc_options = petsc_options)
u = problem.solve()

And the exception is:

AssertionError                            Traceback (most recent call last)
Cell In[1], line 101
     99 petsc_options={}
    100 print(pdofs)
--> 101 problem = LinearProblem(a, l, bcs=[bc_bot], petsc_options = petsc_options)
    102 u = problem.solve()

File ~/miniconda3/envs/dolfinx/lib/python3.11/site-packages/dolfinx/fem/petsc.py:590, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
    588 self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options)
    589 self._A = create_matrix(self._a)
--> 590 self._L = _create_form(L, form_compiler_options=form_compiler_options, jit_options=jit_options)
    591 self._b = create_vector(self._L)
    593 if u is None:
    594     # Extract function space from TrialFunction (which is at the
    595     # end of the argument list as it is numbered as 1, while the
    596     # Test function is numbered as 0)

File ~/miniconda3/envs/dolfinx/lib/python3.11/site-packages/dolfinx/fem/forms.py:188, in form(form, dtype, form_compiler_options, jit_options)
    185         return list(map(lambda sub_form: _create_form(sub_form), form))
    186     return form
--> 188 return _create_form(form)

File ~/miniconda3/envs/dolfinx/lib/python3.11/site-packages/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
    180 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
...
--> 136     assert all([d is data[0] for d in data])
    138 mesh = domain.ufl_cargo()
    139 if mesh is None:

AssertionError:
10

You cannot have multiple MeshTags in the same form (at the moment).
However, why don’t you split this top domain into three separate integer tags:

  1. facets in top but not point (tag 1)
  2. facets in top and point (tag 2)
  3. facets in point but not top (tag 3)

Ie.

top_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 2, top))

p_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 2, point))
pdofs = dolfinx.fem.locate_dofs_topological(V, 2, p_facets)
pt = dolfinx.mesh.meshtags(msh, 2, p_facets, 3)

bot_facets = np.sort(dolfinx.mesh.locate_entities_boundary(msh, 0, bot))
bdofs = dolfinx.fem.locate_dofs_topological(V, 0, bot_facets)
bc_bot = dolfinx.fem.dirichletbc(dolfinx.fem.function.Constant(
    msh, (0.0, 0.0, 0.0)), bdofs, V)     # fixed bottom


top_not_p = np.setdiff1d(top_facets, p_facets)
p_not_top = np.setdiff1d(p_facets, top_facets)
both = np.intersect1d(top_facets, p_facets)
facets = [top_not_p, p_not_top, both]
values = [np.full_like(top_not_p, 1), np.full_like(
    p_not_top, 2), np.full_like(both, 3)]
all_facets = np.hstack(facets).astype(np.int32)
all_values = np.hstack(values).astype(np.int32)
sort_order = np.argsort(all_facets)
mt = dolfinx.mesh.meshtags(msh, msh.topology.dim-1, all_facets[sort_order],
                           all_values[sort_order])

Which you can then create appropriate integrals for:
ds((1,2,3)), ds((2,3)) etc

11

Many thanks @dokken

That’s exactly what I was looking for