Imposing a given force on a part of the domain boundary

variational formulation
1

Good morning,

I try to adapt the example "Setting multiple Dirichlet, Neumann and Robin condition of the Fenicsx documentation to my problem.

I want to apply a uniform force ((1, 0) on the boundary x[0] = 0 of a rectangle domain (L x W).
To do so I mark the factes as explained in “robin_neumann_dirichlet.html”.
Then I define the linear part of the variationnal formulation using ufl.ds(4)
Unfortunately I obtained a null displacement solution, which means the prescribed force are not taken into account.
If I use ufl.ds, (ie integration over all the boundary) the solution is no more zero.
I tried to integrate over ufl.ds(1), ufl.ds(2), … but I always obtained a naught displacement solution.

I tied to find a solution on the forum but I did succeed.
I apologize if this topics as already been addressed.
Here is my script.

Thank you in advance for any answer

Regards

Xavier

External libraries

------------------

import os
import numpy as np
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import ufl
from petsc4py.PETSc import ScalarType
from dolfinx import io

Problem data

L = 1.
W = 0.2
mu = 1.
lmbda = 1.25
p0 = 10.
pL = 0.
nx = 20
ny = 6

deg_u = 1

Def strain tensor

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

Def stress tensor

def sigma(u):
return lmbda * ufl.nabla_div(u) * ufl.Identity(2) + 2. * mu * epsilon(u)

Generating the mesh

Omega = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array([0., 0.]), np.array([L, W])],
[nx, ny], cell_type=mesh.CellType.triangle)

Identifying the facets to apply boundary conditions

boundaries = [(1, lambda x : np.isclose(x[1], 0.)),
(2, lambda x : np.isclose(x[0], L )),
(3, lambda x : np.isclose(x[1], W )),
(4, lambda x : np.isclose(x[0], 0.))]
facet_indices, facet_markers = ,
fdim = Omega.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(Omega, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(Omega, fdim, facet_indices[sorted_facets],
facet_markers[sorted_facets])
ds = ufl.Measure(“ds”, domain=Omega, subdomain_data=facet_tag)

Function spaces and functions

V = fem.VectorFunctionSpace(Omega, (“CG”, deg_u))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

Problem definition

Define boundary conditions

ud = np.array([0., 0.], dtype = ScalarType)
facets = facet_tag.find(1)
dofs = fem.locate_dofs_topological(V, fdim, facets)
bc1 = fem.dirichletbc(ud, dofs, V)
facets = facet_tag.find(3)
dofs = fem.locate_dofs_topological(V, fdim, facets)
bc3 = fem.dirichletbc(ud, dofs, V)

Define variational problem

a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
f = fem.Constant(Omega, ScalarType((1., 0.)))
L = ufl.dot(f, v) * ufl.ds(4)

Computing the solution

problem = fem.petsc.LinearProblem(a, L, bcs=[bc1, bc3],
petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})
uh = problem.solve()
uh.name = “Displacement”

Computing the norm of the displacement

norm_uh = np.sqrt(Omega.comm.allreduce(fem.assemble_scalar(fem.form(ufl.inner(uh, uh) * ufl.dx)), op=MPI.SUM))
print("norm uh = "+str(norm_uh))

Writing XDMFFile for Paraview

with io.XDMFFile(Omega.comm, “marking_boundary.xdmf”, “w”) as xdmf:
xdmf.write_mesh(Omega)
xdmf.write_function(uh)

2

It seems I find the solution.
I just replace ufl.ds(4) by ds(4) and it works.
I just have to use the measure I defined with the line

ds = ufl.Measure(“ds”, domain=Omega, subdomain_data=facet_tag)

I close the topics

1 Like