Problem when trying to solve a Plane Tension problem (2D)

dolfinx
1

Good evening, excuse me, I am an engineering student who wants to learn more about finite elements. To do this, I set out to create a small Python code in Colab. The program will solve a Plane Stress problem (2D). I have a formulation, but I am having trouble finding the solution, I can’t use LinearProblem(a, l, bcs=[bc]). Here’s part of the code:

# Define the function space
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim,)))

# Create the zero vertical displacement constraint at the interface
zero = dolfinx.fem.Function(V)
with zero.vector.localForm() as loc:
    loc.set(0.0)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define stress and strain
def eps(v):
    return sym(grad(v))

def sigma(v, lmbda, mu):
    return lmbda * tr(eps(v)) * Identity(2) + 2 * mu * eps(v)

# Material properties
E_presa = 24000000  # kN/m^2
nu_presa = 0.2
rho_presa = 22  # kN/m^3

E_suelo = 24000000  # kN/m^2
nu_suelo = 0.2
rho_suelo = 26  # kN/m^3

# Define constants for constitutive relations
lmbda_presa = E_presa * nu_presa / ((1 + nu_presa) * (1 - 2 * nu_presa))
mu_presa = E_presa / (2 * (1 + nu_presa))

lmbda_suelo = E_suelo * nu_suelo / ((1 + nu_suelo) * (1 - 2 * nu_suelo))
mu_suelo = E_suelo / (2 * (1 + nu_suelo))

# Boundary conditions for sliding at the dam-soil interface
def interface_pres_suelo(x):
    return np.isclose(x[1], 0) & (x[0] >= 0) & (x[0] <= 40)

# Zero vertical displacement boundary condition
dofs_y = dolfinx.fem.locate_dofs_geometrical((V.sub(1), V), interface_pres_suelo) # Locate dofs in the y direction
bc_vertical = dolfinx.fem.dirichletbc(zero, dofs_y, V.sub(1))

# Clamped boundary conditions on the lateral edges and bottom of the soil
def empotrado(x):
    return np.isclose(x[0], -ancho_suelo) | np.isclose(x[0], 40 + ancho_suelo) | np.isclose(x[1], -profundidad_suelo)

# Clamped boundary conditions on the lateral edges and bottom of the soil
dofs_suelo = dolfinx.fem.locate_dofs_geometrical(V, empotrado) # Locate the dofs for clamping
bc_suelo = dolfinx.fem.dirichletbc(zero, dofs_suelo)

# Boundary conditions for hydrostatic pressure on the left face of the dam
x = SpatialCoordinate(mesh)
presion_hidrost = 1000 * 9.81 * (50 - x[1])  # Hydrostatic pressure distribution

# Define the measure ds for boundary terms
ds = Measure('ds', domain=mesh)

# External forces (self-weight of the dam and water)
f_presa = dolfinx.fem.Constant(mesh, (0, -rho_presa * 9.81))
f_suelo = dolfinx.fem.Constant(mesh, (0, -rho_suelo * 9.81))

n = FacetNormal(mesh)

# Variational formulation with boundary conditions
a = inner(sigma(u, lmbda_presa, mu_presa), eps(v)) * dx + inner(sigma(u, lmbda_suelo, mu_suelo), eps(v)) * dx
l = inner(f_presa, v) * dx + inner(f_suelo, v) * dx + dot(presion_hidrost * n, v) * ds

# Solve the problem
problem = dolfinx.fem.LinearProblem(a, l, bcs=[bc])
u_sol = problem.solve()

Seeking an alternative solution, I used another method, but it does not run directly; instead, it restarts the environment. Code:

# Define the function space
V = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim,)))

# Create the zero vertical displacement constraint at the interface
zero = fem.Function(V)
with zero.vector.localForm() as loc:
    loc.set(0.0)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define stress and strain
def eps(v):
    return sym(grad(v))

def sigma(v, lmbda, mu):
    return lmbda * tr(eps(v)) * Identity(2) + 2 * mu * eps(v)

# Material properties
E_presa = 24000000  # kN/m^2
nu_presa = 0.2
rho_presa = 22  # kN/m^3

E_suelo = 24000000  # kN/m^2
nu_suelo = 0.2
rho_suelo = 26  # kN/m^3

# Define constants for constitutive relations
lmbda_presa = E_presa * nu_presa / ((1 + nu_presa) * (1 - 2 * nu_presa))
mu_presa = E_presa / (2 * (1 + nu_presa))

lmbda_suelo = E_suelo * nu_suelo / ((1 + nu_suelo) * (1 - 2 * nu_suelo))
mu_suelo = E_suelo / (2 * (1 + nu_suelo))

# Boundary conditions for sliding at the dam-soil interface
def interface_pres_suelo(x):
    return np.isclose(x[1], 0) & (x[0] >= 0) & (x[0] <= 40)

# Locate DOFs at the interface to set vertical displacement (y direction)
dofs_y = fem.locate_dofs_geometrical(V, interface_pres_suelo)

# Create Dirichlet condition for vertical displacement (y = 0) at the dam-soil interface
bc_vertical = fem.dirichletbc(PETSc.ScalarType((0.0, 0.0)), dofs_y, V)

# Clamped boundary conditions on the lateral edges and bottom of the soil
ancho_suelo = 80  # Extension of the soil on both sides of the dam
profundidad_suelo = 150  # Depth of the soil

def empotrado(x):
    return np.isclose(x[0], -ancho_suelo) | np.isclose(x[0], 40 + ancho_suelo) | np.isclose(x[1], -profundidad_suelo)

dofs_suelo = fem.locate_dofs_geometrical(V, empotrado)  # Locate the dofs for clamping
bc_suelo = fem.dirichletbc(zero, dofs_suelo)

# Boundary conditions for hydrostatic pressure on the left face of the dam
x = SpatialCoordinate(mesh)
presion_hidrost = 1000 * 9.81 * (50 - x[1])  # Hydrostatic pressure distribution

# Define the measure ds for boundary terms
ds = Measure('ds', domain=mesh)

# Define the normal vector for boundary integration
n = FacetNormal(mesh)

# External forces (self-weight of the dam and water)
f_presa = fem.Constant(mesh, PETSc.ScalarType((0, -rho_presa * 9.81)))
f_suelo = fem.Constant(mesh, PETSc.ScalarType((0, -rho_suelo * 9.81)))

# Variational formulation with boundary conditions
a = inner(sigma(u, lmbda_presa, mu_presa), eps(v)) * dx + inner(sigma(u, lmbda_suelo, mu_suelo), eps(v)) * dx
l = inner(f_presa, v) * dx + inner(f_suelo, v) * dx + dot(presion_hidrost * n, v) * ds

# Create matrices and vectors for the linear system
A = fem.assemble_matrix(fem.form(a), bcs=[bc_vertical, bc_suelo])
b = fem.assemble_vector(fem.form(l))

# Define the initial solution before applying lifting
u_sol = fem.Function(V)

# Convert `b` to a numpy ndarray for apply_lifting
b_array = b.array

# Apply lifting to the matrix and vector of the system
fem.apply_lifting(b_array, [fem.form(a)], bcs=[[bc_vertical, bc_suelo]], x0=[u_sol.vector.array], scale=1.0,
                  constants=[], coeffs=[])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(b, [bc_vertical, bc_suelo])

# Solve the linear system using PETSc
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.solve(b, u_sol.vector)

I know it’s a conceptual error. Could you help me with documentation to better understand the problem I’m having and the best way to approach it? Thanks in advance, greetings from Bolivia.
link Colab: Google Colab

2

Please note that it is preferred that a reproducible example is posted directly at the forum, rather than through external links (as links might disappear in the future if you change the location/ownership of the code.

Secondly, the code you have provided does not explain what doesn’t work. Are you getting an error message when running either of the codes above?

3

Thank you, I will keep this in mind for future questions. The problem was an incorrect handling of syntax; by adding a couple of lines, it was resolved. However, now I have a conceptual problem. I’m not sure if the way I define my boundary conditions is correct. The problem involves a dam resting on the ground, subjected to hydrostatic pressure. I assumed the following:

  • Zero vertical displacement at the dam-soil interface: At the interface between the dam and the soil, a zero vertical displacement is imposed, allowing horizontal sliding.
  • Clamping on the lateral edges and bottom of the soil: Zero displacements (clamping) are imposed on the lateral edges of the soil and at the bottom.
  • Hydrostatic pressure on the left face of the dam: A hydrostatic pressure is applied that varies linearly with the water height on the left face of the dam, affecting only up to the water height.
    Code:
try:
    import dolfinx
except ImportError:
    !wget "https://github.com/fem-on-colab/fem-on-colab.github.io/raw/df47889/releases/fenicsx-install-real.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
    import dolfinx

try:
    import gmsh
except ImportError:
    !wget "https://github.com/fem-on-colab/fem-on-colab.github.io/raw/df47889/releases/gmsh-install.sh" -O "/tmp/gmsh-install.sh" && bash "/tmp/gmsh-install.sh"
    import gmsh

!pip install meshio

MESH GENERATION WITH GMSH

import gmsh
import numpy as np
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx import mesh

# Initialize Gmsh
gmsh.initialize()
gmsh.model.add("Dam and Soil")

# Dam and soil parameters
height_triangle = 42  # meters
base_triangle = 34  # meters
base_rectangle = 40.0  # meters
height_rectangle = 8.0  # meters
soil_depth = 3 * (height_triangle + height_rectangle)  # Three times the total height of the dam
soil_width = 2 * base_rectangle  # Twice the rectangular base

# Definition of dam points
A = gmsh.model.geo.addPoint(0, 0, 0)
B = gmsh.model.geo.addPoint(base_rectangle, 0, 0)
C = gmsh.model.geo.addPoint(base_rectangle, height_rectangle, 0)
D = gmsh.model.geo.addPoint(0, height_rectangle, 0)
E = gmsh.model.geo.addPoint(base_triangle, height_rectangle, 0)
F = gmsh.model.geo.addPoint(0, height_triangle + height_rectangle, 0)

# Definition of dam lines
l1 = gmsh.model.geo.addLine(A, B)  # Base line of the rectangular part
l2 = gmsh.model.geo.addLine(B, C)  # Right vertical line of the rectangular part
l3 = gmsh.model.geo.addLine(C, E)  # Top right line of the rectangular part
l4 = gmsh.model.geo.addLine(E, F)  # Diagonal line of the triangular part
l5 = gmsh.model.geo.addLine(F, D)  # Left vertical line of the triangular part
l6 = gmsh.model.geo.addLine(D, A)  # Left vertical line of the rectangular part

# Definition of dam surfaces
loop_dam = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4, l5, l6])
surface_dam = gmsh.model.geo.addPlaneSurface([loop_dam])

# Definition of soil
P1 = gmsh.model.geo.addPoint(-soil_width, -soil_depth, 0)
P2 = gmsh.model.geo.addPoint(base_rectangle + soil_width, -soil_depth, 0)
P3 = gmsh.model.geo.addPoint(base_rectangle + soil_width, 0, 0)
P4 = gmsh.model.geo.addPoint(-soil_width, 0, 0)

l8 = gmsh.model.geo.addLine(P1, P2)
l9 = gmsh.model.geo.addLine(P2, P3)
l10 = gmsh.model.geo.addLine(P3, P4)
l11 = gmsh.model.geo.addLine(P4, P1)

loop_soil = gmsh.model.geo.addCurveLoop([l8, l9, l10, l11])
surface_soil = gmsh.model.geo.addPlaneSurface([loop_soil])

# Synchronize geometry
gmsh.model.geo.synchronize()

# Define physical groups for applying loads and boundary conditions
gmsh.model.addPhysicalGroup(2, [surface_dam], 1)  # Dam
gmsh.model.addPhysicalGroup(2, [surface_soil], 2)  # Soil

# Additional physical groups
gmsh.model.addPhysicalGroup(1, [l5], 3)  # Left face of the dam (for water pressure)
gmsh.model.addPhysicalGroup(1, [l1], 4)  # Base of the dam (for dam-soil interaction)
gmsh.model.addPhysicalGroup(1, [l8], 5)  # Soil surface at the water interface (water weight)

# Specify global characteristic length
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 10)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 10)

# Generate the 2D mesh
gmsh.model.mesh.generate(2)

# Export the mesh to .msh format for use in DOLFINx
gmsh.write("dam_soil.msh")

# Read the mesh with DOLFINx
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("dam_soil.msh", MPI.COMM_WORLD, gdim=2)

print("Geometry defined and mesh generated successfully.")

Problem Formulation and Solution

import dolfinx
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.mesh import create_rectangle, CellType
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx.fem.petsc import LinearProblem
from ufl import TrialFunction, TestFunction, dx, inner, sym, grad, Identity, tr, dot, FacetNormal, SpatialCoordinate, Measure

print(f"Mesh geometry dimension: {mesh.geometry.dim}")
print(f"Mesh coordinates: {mesh.geometry.x}")

# Define the function space
V = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim,)))

# Check mesh information
print(f"Number of cells in the mesh: {mesh.topology.index_map(mesh.topology.dim).size_local}")
print(f"Number of vertices in the mesh: {mesh.topology.index_map(0).size_local}")
print(f"Number of DOFs per node (Expected for 2D Elasticity): {V.dofmap.index_map_bs}")
print(f"Total number of DOFs: {V.dofmap.index_map.size_local * V.dofmap.index_map_bs}")

# Create zero vertical displacement constraint at the interface
zero = fem.Function(V)
with zero.vector.localForm() as loc:
    loc.set(0.0)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define stress and strain
def eps(v):
    return sym(grad(v))

def sigma(v, lmbda, mu):
    return lmbda * tr(eps(v)) * Identity(2) + 2 * mu * eps(v)

# Material properties
E_dam = 24000000  # kN/m^2
nu_dam = 0.2
rho_dam = 2.2426  # kN/m^3

E_soil = 24000000  # kN/m^2
nu_soil = 0.2
rho_soil = 2.6504  # kN/m^3

# Define constants for constitutive relations
lmbda_dam = E_dam * nu_dam / ((1 + nu_dam) * (1 - 2 * nu_dam))
mu_dam = E_dam / (2 * (1 + nu_dam))

lmbda_soil = E_soil * nu_soil / ((1 + nu_soil) * (1 - 2 * nu_soil))
mu_soil = E_soil / (2 * (1 + nu_soil))

# Boundary conditions for sliding at the dam-soil interface
def interface_dam_soil(x):
    return np.isclose(x[1], 0) & (x[0] >= 0) & (x[0] <= 40)

# Vertical displacement boundary condition
dofs_y = fem.locate_dofs_geometrical((V.sub(1), V), interface_dam_soil)  # Locate dofs in the y direction
bc_vertical = fem.dirichletbc(zero, dofs_y, V.sub(1))

# Clamped boundary conditions on the lateral edges and bottom of the soil
soil_width = 80  # Soil extension on both sides of the dam
soil_depth = 150  # Soil depth

def clamped(x):
    return np.isclose(x[0], -soil_width) | np.isclose(x[0], 40 + soil_width) | np.isclose(x[1], -soil_depth)

dofs_soil = fem.locate_dofs_geometrical(V, clamped)  # Locate the dofs for clamping
bc_soil = fem.dirichletbc(zero, dofs_soil)

from ufl import conditional, le
# Boundary conditions for hydrostatic pressure on the left face of the dam
water_height = 40  # Example: water height
rho_water = 1
g = 9.81  # Acceleration due to gravity m/s2

if water_height > 0:
    x = SpatialCoordinate(mesh)
    # Define hydrostatic pressure applying the condition that it only affects up to the water height
    hydrostatic_pressure = rho_water * g * (water_height - x[1])
else:
    hydrostatic_pressure = 0  # No hydrostatic pressure if water height is zero

# Define the measure ds for boundary terms
ds = Measure('ds', domain=mesh)

# Define the normal vector for boundary integration
n = FacetNormal(mesh)

# External forces (self-weight of the dam and water)
f_dam = fem.Constant(mesh, (0, -rho_dam * 9.81))
f_soil = fem.Constant(mesh, (0, -rho_soil * 9.81))

# Variational formulation with boundary conditions
a = inner(sigma(u, lmbda_dam, mu_dam), eps(v)) * dx + inner(sigma(u, lmbda_soil, mu_soil), eps(v)) * dx
l = inner(f_dam, v) * dx + inner(f_soil, v) * dx + dot(hydrostatic_pressure * n, v) * ds

# PETSc solver options
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps"
}

# Solve the problem using LinearProblem
problem = LinearProblem(a, l, bcs=[bc_vertical, bc_soil], petsc_options=petsc_options)
u_sol = problem.solve()

print(f"Solution vector size u_sol: {u_sol.vector.array.size}")
print(u_sol.vector.array)

The question arises because I was testing the code by changing the water height, and for an empty dam, the results don’t seem reasonable to me. I would appreciate it if you could provide some bibliography, documentation, or an example. Thank you in advance.