Issue with Mesh Study in Elasticity Problem: Inconsistent Results with Increased Elements

Hey all, I’m working on a simple elasticity problem and encountering an issue during my mesh study. The code runs fine with a low number of elements (39904), but when I increase the number of elements (158564 elements), I get constant results across the entire domain. Here are the two meshes I’m using


The first mesh works well but the second mesh doesn’t.

Here are the outputs for both meshes:


I don’t know what the problem is. Is that because the direct solver that I am using?

And here is the code:

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
                         locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
                 div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type

# .........................................  Non-Dimensionalization Parameters
Fl = 1e-3
Fs = 1e-5
Fp = 1e-4
# # # ....................
C_star = Fl ** 2 / Fs
h_star = 1 / Fl
p_star = Fl ** 3
F_star = Fs
# -------------------------Material Properties
E = 1 * 130e9 * C_star
nu = 0.28
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# ......................................... Reading mesh file
h = 0.76e-3 / Fl
d1 = 1.13e-3 / Fl
tb = 0.1e-3 / Fl
########################
H = h + tb
mesh, cell_markers, facet_markers = gmshio.read_from_msh("Pyramid.msh", MPI.COMM_WORLD,
                                                         gdim=3)  # gmshio.read_from_msh always needs physical group in .msh file
# ......................................... Define Loading
Force = 200.00000 / Fs
Area2 = 1.2769000000000007e-06 / Fl ** 2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
MIXEL = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
Space = functionspace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
u = TrialFunction(Space)
del_u = TestFunction(Space)
# ....................................... Mark boundary regions
tol = 1 * 1e-16


def Bottom(x):
    return np.isclose(x[2], 0, atol=tol)


# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, 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 = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# # ...................................... BCs
# Clamp BC for left side
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
clamp_dofs = locate_dofs_topological(Space, mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(np.array([0, 0, 0], dtype=default_scalar_type), clamp_dofs, Space)
bc = [bc0]  # BC 4
# ......................................... Define tensor
delta = Identity(3)
i, j, k, l = indices(4)


# Strain Tensor
def StrainT(u):
    return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])


# Stress Tensor
def Elas(la, mu):
    return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
                     (i, j, k, l))  # Elasticity tensor


def StressT(u):
    return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])


# ......................................... Define Variational Formulation
VarForm = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
a = VarForm
L = tr_S2[i] * del_u[i] * ds(2)
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = assemble_vector(linear_form)
apply_lifting(b, [bilinear_form], [bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bc)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

problem = LinearProblem(a, L, bcs=bc,
                        petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()

with VTXWriter(mesh.comm, ("Pyramid/DIS_Pyramid.bp"), [Result]) as vtx:
    vtx.write(0.0)

That’s possible. Please try using superlu or superlu_dist and check if you still get the same result.

Apart from this, it’s impossible for anyone to help you any further, because we don’t have access to the meshes and thus we can’t run your code.

Thanks for the reply, I forgot to provide the mesh. The following code generates the mesh in GMSH.


import gmsh
import sys

##################
gmsh.initialize()
gmsh.model.add("Model")
# =====================================    Goccmetry      ======================
import gmsh
import sys
import math

gmsh.initialize()
Fl = 1e-3
lc = 0.0
h = 0.76e-3 / Fl
d1 = 1.13e-3 / Fl
d2 = d1
lb = 0.1e-3 / Fl
tb = 0.1e-3 / Fl
theta = math.atan(h / (0.5 * (0.00272 / Fl - d1)))  #54.74
Points = [[0, 0, 0], [2 * lb + d1 + 2 * h / math.tan(theta), 0, 0],
          [2 * lb + d1 + 2 * h / math.tan(theta), 2 * lb + d2 + 2 * h / math.tan(theta), 0],
          [0, 2 * lb + d2 + 2 * h / math.tan(theta), 0]]
print(d1 + 2 * h / math.tan(theta))
PN = 1  # point Number
LN = 1  # Line Number
SN = 1  # Surface Number
VN = 1  # Volum Number
Vol_list = []
Sur_list = []
Sur_list2 = []
Sur_list3 = []
Sur_list4 = []
Sur_list5 = []
Sur_list6 = []
Edge_list_1 = []  # top edges
Edge_list_2 = []  # Bottom edges
Edge_list_3 = []  # Side edges
Edge_list_4 = []  # in-plane edges
Edge_list_5 = []  # corner edges
Edge_list_6 = []  # Base edges
i=0
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta), Points[i][1] + lb + h / math.tan(theta), tb + h,
                        lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta) + d1, Points[i][1] + lb + h / math.tan(theta),
                        tb + h, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_1.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta) + d1, Points[i][1] + lb + h / math.tan(theta) + d2,
                        tb + h, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_1.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + h / math.tan(theta), Points[i][1] + lb + h / math.tan(theta) + d2,
                        tb + h, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_1.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
Edge_list_1.append(LN)
LN += 1
#########
gmsh.model.geo.addPoint(Points[i][0] + lb, Points[i][1] + lb, tb, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1, Points[i][1] + lb, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_2.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1,
                        Points[i][1] + lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_2.append(LN)
LN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb, Points[i][1] + lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
Edge_list_2.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
Edge_list_2.append(LN)
LN += 1
####################
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
Edge_list_3.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
Edge_list_3.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
Edge_list_3.append(LN)
LN += 1
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
Edge_list_3.append(LN)
LN += 1
####################
gmsh.model.geo.addCurveLoop([LN - 12, LN - 11, LN - 10, LN - 9], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 12, (LN - 3), -(LN - 8), -(LN - 4)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 11, (LN - 2), -(LN - 7), -(LN - 3)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 10, (LN - 1), -(LN - 6), -(LN - 2)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 9, (LN - 4), -(LN - 5), -(LN - 1)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 8, LN - 7, LN - 6, LN - 5], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list.append(SN - 1)
print(SN)
gmsh.model.geo.addSurfaceLoop([SN - 1, SN - 2, SN - 3, SN - 4, SN - 5, SN - 6], VN)
gmsh.model.geo.addVolume([VN], VN)
Vol_list.append(VN)
VN += 1
gmsh.model.geo.addPoint(Points[i][0], Points[i][1], tb, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1 + lb, Points[i][1], tb, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0] + lb + 2 * h / math.tan(theta) + d1 + lb,
                        Points[i][1] + 2 * lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Points[i][0], Points[i][1] + 2 * lb + 2 * h / math.tan(theta) + d2, tb, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 4, PN - 3, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 2, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addCurveLoop([LN - 16, LN - 7, -(LN - 4), -(LN - 8)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 15, LN - 6, -(LN - 3), -(LN - 7)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list3.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 14, LN - 5, -(LN - 2), -(LN - 6)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list4.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 13, LN - 8, -(LN - 1), -(LN - 5)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list5.append(SN - 1)

gmsh.model.geo.synchronize()
Corner_points = [PN - 4, PN - 3, PN - 2, PN - 1]
Left_down_corner = gmsh.model.getValue(0, Corner_points[0], [])
Right_down_corner = gmsh.model.getValue(0, Corner_points[1], [])
Right_up_corner = gmsh.model.getValue(0, Corner_points[2], [])
Left_up_corner = gmsh.model.getValue(0, Corner_points[3], [])

gmsh.model.geo.addPoint(Left_down_corner[0], Left_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_down_corner[0], Right_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_up_corner[0], Right_up_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Left_up_corner[0], Left_up_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
LN += 1
Edge_list_6.append(LN-1)
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 4, PN - 3, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 2, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addCurveLoop([LN - 12, LN - 7, -(LN - 4), -(LN - 8)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 11, LN - 6, -(LN - 3), -(LN - 7)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list3.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 10, LN - 5, -(LN - 2), -(LN - 6)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list4.append(SN - 1)
gmsh.model.geo.addCurveLoop([LN - 9, LN - 8, -(LN - 1), -(LN - 5)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list5.append(SN - 1)
gmsh.model.geo.synchronize()
Corner_points = [PN - 12, PN - 11, PN - 10, PN - 9]
Left_down_corner = gmsh.model.getValue(0, Corner_points[0], [])
Right_down_corner = gmsh.model.getValue(0, Corner_points[1], [])
Right_up_corner = gmsh.model.getValue(0, Corner_points[2], [])
Left_up_corner = gmsh.model.getValue(0, Corner_points[3], [])
gmsh.model.geo.addPoint(Left_down_corner[0], Left_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_down_corner[0], Right_down_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Right_up_corner[0], Right_up_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addPoint(Left_up_corner[0], Left_up_corner[1], 0, lc, PN)
PN += 1
gmsh.model.geo.addLine(PN - 8, PN - 4, LN)
LN += 1
Edge_list_4.append(LN-1)
gmsh.model.geo.addLine(PN - 7, PN - 3, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 6, PN - 2, LN)
LN += 1
Edge_list_4.append(LN - 1)
gmsh.model.geo.addLine(PN - 5, PN - 1, LN)
LN += 1
Edge_list_4.append(LN - 1)
# ################3
gmsh.model.geo.addLine(PN - 4, PN - 16, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 15, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 14, LN)
LN += 1
Edge_list_6.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 13, LN)
LN += 1
Edge_list_6.append(LN - 1)
################
gmsh.model.geo.addLine(PN - 4, PN - 3, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 3, PN - 2, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 2, PN - 1, LN)
LN += 1
Edge_list_5.append(LN - 1)
gmsh.model.geo.addLine(PN - 1, PN - 4, LN)
LN += 1
Edge_list_5.append(LN - 1)
# ################################
gmsh.model.geo.addCurveLoop([LN - 28, LN - 20, (LN - 12), (LN - 8)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)
Sur_list5.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 16, LN - 11, -(LN - 4), -(LN - 12)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 27, LN - 19, (LN - 11), (LN - 7)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)
Sur_list3.append(SN - 1)


gmsh.model.geo.addCurveLoop([LN - 4, LN - 7, -(LN - 36), -(LN - 8)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list2.append(SN - 1)

gmsh.model.geo.addSurfaceLoop(Sur_list2, VN)
gmsh.model.geo.addVolume([VN], VN)
Vol_list.append(VN)
VN += 1
############################
gmsh.model.geo.addCurveLoop([LN - 3, LN - 6, -(LN - 35), -(LN - 7)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list3.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 15, LN - 10, -(LN - 3), -(LN - 11)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list3.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 10, LN - 6, (LN - 26), (LN - 18)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list3.append(SN - 1)
gmsh.model.geo.addSurfaceLoop(Sur_list3, VN)
gmsh.model.geo.addVolume([VN], VN)
Vol_list.append(VN)
VN += 1
Sur_list4.append(SN - 1)
################################
gmsh.model.geo.addCurveLoop([LN - 14, LN - 9, -(LN - 2), -(LN - 10)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list4.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 2, LN - 5, -(LN - 34), -(LN - 6)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list4.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 9, LN - 5, (LN - 25), (LN - 17)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list4.append(SN - 1)
gmsh.model.geo.addSurfaceLoop(Sur_list4, VN)
gmsh.model.geo.addVolume([VN], VN)
Vol_list.append(VN)
VN += 1
Sur_list5.append(SN - 1)
###############################
gmsh.model.geo.addCurveLoop([LN - 13, LN - 12, -(LN - 1), -(LN - 9)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list5.append(SN - 1)

gmsh.model.geo.addCurveLoop([LN - 1, LN - 8, -(LN - 33), -(LN - 5)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list5.append(SN - 1)

gmsh.model.geo.addSurfaceLoop(Sur_list5, VN)
gmsh.model.geo.addVolume([VN], VN)
Vol_list.append(VN)
VN += 1
######################################
Sur_list6 = [SN-21, SN-9, SN-8, SN-4, SN-1]
gmsh.model.geo.addCurveLoop([LN - 4, LN - 3, (LN - 2), (LN - 1)], SN)
gmsh.model.geo.addPlaneSurface([SN], SN)
SN += 1
Sur_list6.append(SN - 1)

gmsh.model.geo.addSurfaceLoop(Sur_list6, VN)
gmsh.model.geo.addVolume([VN], VN)
Vol_list.append(VN)
VN += 1
# =====================================    Mesh Option      ======================
## UP
N_up = 50
type_up = "Bump"
R_up = 0.5
## Bottom
N_bottom = N_up
type_bottom = "Bump"
R_bottom = 0.5
## Side
N_side = 60
type_side = "Bump"
R_side = 1
####
N_inPlane = 6
###
N_base= 6
# ----------------------- UP
for i in Edge_list_1:
    gmsh.model.geo.mesh.setTransfiniteCurve(i, N_up, type_up, R_up)
for i in Edge_list_2:
    gmsh.model.geo.mesh.setTransfiniteCurve(i, N_bottom, type_bottom, R_bottom)
for i in Edge_list_3:
    gmsh.model.geo.mesh.setTransfiniteCurve(i, N_side, type_side, R_side)

for i in Edge_list_4:
    gmsh.model.geo.mesh.setTransfiniteCurve(i, N_inPlane)#, type_side, R_side)

for i in Edge_list_5:
    gmsh.model.geo.mesh.setTransfiniteCurve(i, N_bottom)#, type_side, R_side)

for i in Edge_list_6:
    gmsh.model.geo.mesh.setTransfiniteCurve(i, N_base)#, type_side, R_side)


for i in Sur_list + Sur_list2 +Sur_list3 + Sur_list4 + Sur_list5 + Sur_list6:
    gmsh.model.geo.mesh.setTransfiniteSurface(i)

for i in Sur_list + Sur_list2 +Sur_list3 + Sur_list4 + Sur_list5 + Sur_list6:
    gmsh.model.geo.mesh.setRecombine(2, i)


for i in Vol_list:
    gmsh.model.geo.mesh.setTransfiniteVolume(i)


gmsh.model.geo.remove_all_duplicates()
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(3, Vol_list, name="Pyramid")
gmsh.option.setNumber('Mesh.MeshSizeMax', 0.005)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3)
# ##############
gmsh.write("Pyramid.msh")
######################
if '-nopopup' not in sys.argv:
    gmsh.fltk.run()
#################################
gmsh.finalize()


@francesco-ballarin , When I use the superlu_dist the code terminates due to the low memory error, even with the first mesh that has lower number of elements. And when I use superlu, I get the following

Thanks

Then I would guess you are running out of RAM. Is there any way you can try the same code on a workstation with more RAM?

Thanks @francesco-ballarin for your reply. I can’t upgrade the RAM on my PC at the moment, but I can plan to get a system with higher RAM in the future. However, the issue is that it doesn’t work for a mesh with a lower number of elements (39904), which I can run with mumps. Does this mean that using superlu_dist or superlu works better for large systems (high DOF) but requires more RAM than using mumps?

Do you think there’s no problem with my code, and that the strange results in my sumulation are just something that happens in large models when using mumps?

I have one more question. I’m currently using a direct solver in this code. Should I try using an iterative solver instead? Could that potentially solve the issue?

Thanks,

If the issue is the direct solver and/or the RAM, then using an iterative solver could help. Unfortunately I can’t help much further, since neither solid mechanics nor iterative solvers are my field of expertise. You may want to wait for further comments by someone more knowledgeable than me in those fields.

Thanks @francesco-ballarin for your time. I will wait to see if I get more comments.

Hey all,
I still need help with this issue and I am not sure if going after setting up an iterative solver is the best way to solve this problem. I would appreciate any comments on that.
Thanks,

Could you print the output of problem.solver.getConvergedReason() after the linear solver has given you a constant solution.

As a follow-up, how much RAM and how many processors do you have, and how many processes are you using with mpirun?

the out put is : 4
This is a linear problem. Shouldn’t that be 1?

I have 128 GB RAM and If i use one processor it is much faster (I don’t know why), but I can use up to 8.

Why is that super slow when I use more processors?

And here is the out put I get from
dolfinx.common.list_timings(mesh.comm,[dolfinx.common.TimingType.wall])

Thanks

I do not have the sufficient memory on my laptop to run the code, so this would have to wait for me until I am on a bigger system.

How did you install dolfinx? using conda, docker, spack or something else?

Btw. your problem has 3 922 347 dofs, which is really quite alot to throw at a direct solver.

Thanks for your help,
I installed it using conda. And, currently I am working on setting an iterative solver for my problem.