Free(): invalid next size (normal)

General
1

Hello everyone, i need help with an error i encountered in a project where I am applying finite element method and topology optimization. Firstly, I exported a file in STEP format to MSH format using the GMSH application. Then, I converted this file to XDMF format. I apply the finite element method to the model first, followed by topology optimization. However, I am encountering an error. I am sharing my code and the error message with you, and I would appreciate any help you can provide.
Here my .msh file GitHub - ctaskum/-ekiller

import meshio
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *

msh = meshio.read(filename="/home/zeynep/FEM/square2.msh")

for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data

for key in msh.cell_data_dict["gmsh:geometrical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:geometrical"][key]

triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", triangle_cells)],
                            cell_data={"name_to_read": [triangle_data]})

meshio.write("mf.xdmf", triangle_mesh)


# XDMF 
mesh = Mesh()
xdmf_file = XDMFFile("mf.xdmf")
xdmf_file.read(mesh)

plot(mesh)
plt.show()


# Algorithmic parameters
niternp = 20 
niter = 80 
pmax = 4   
exponent_update_frequency = 4 
tol_mass = 1e-4 
thetamin = 0.001

# Problem parameters
thetamoy = 0.4 
E = Constant(1)  
nu = Constant(0.3)  
lamda = E*nu/(1+nu)/(1-2*nu) 
mu = E/(2*(1+nu)) 
f = Constant((0, -1, 0)) 


def left(x, on_boundary): 
    return near(x[0], -400) and on_boundary 

def load(x, on_boundary): 
    return near(x[0], 650) and near(x[1], 150, 50) 

# Facets
facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(load).mark(facets, 1) 
ds = Measure("ds", subdomain_data=facets) 


# Function space for density field
V0 = FunctionSpace(mesh, "DG", 0) 
# Function space for displacement
V2 = VectorFunctionSpace(mesh, "CG", 3) 
# Fixed boundary condtions
bc = DirichletBC(V2, Constant((0, 0, 0)), left) 

p = Constant(1) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = Constant(1) # Lagrange multiplier for volume constraint


thetaold = Function(V0, name="Density")
thetaold.interpolate(Constant(thetamoy))
coeff = thetaold**p
theta = Function(V0)

volume = assemble(Constant(1.)*dx(domain=mesh))
avg_density_0 = assemble(thetaold*dx)/volume # initial average density
avg_density = 0.


def eps(v):
    return sym(grad(v))
def sigma(v):
    return coeff*(lamda*div(v)*Identity(3)+2*mu*eps(v))
def energy_density(u, v):
    return inner(sigma(u), eps(v))


# Inhomogeneous elastic variational problem
u_ = TestFunction(V2)
du = TrialFunction(V2)
a = inner(sigma(u_), eps(du))*dx
L = dot(f, u_)*ds(1)


def local_project(v, V):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    u = Function(V)
    solver.solve_local_rhs(u)
    return u

def update_theta():
    theta.assign(local_project((p*coeff*energy_density(u, u)/lagrange)**(1/(p+1)), V0))
    thetav = theta.vector().get_local()
    theta.vector().set_local(np.maximum(np.minimum(1, thetav), thetamin))
    theta.vector().apply("insert")
    avg_density = assemble(theta*dx)/volume
    return avg_density



def update_lagrange_multiplier(avg_density):
    avg_density1 = avg_density
    # Initial bracketing of Lagrange multiplier
    if (avg_density1 < avg_density_0):
        lagmin = float(lagrange)
        while (avg_density < avg_density_0):
            lagrange.assign(Constant(lagrange/2))
            avg_density = update_theta()
        lagmax = float(lagrange)
    elif (avg_density1 > avg_density_0):
        lagmax = float(lagrange)
        while (avg_density > avg_density_0):
            lagrange.assign(Constant(lagrange*2))
            avg_density = update_theta()
        lagmin = float(lagrange)
    else:
        lagmin = float(lagrange)
        lagmax = float(lagrange)

    # Dichotomy on Lagrange multiplier
    inddico=0
    while ((abs(1.-avg_density/avg_density_0)) > tol_mass):
        lagrange.assign(Constant((lagmax+lagmin)/2))
        avg_density = update_theta()
        inddico += 1;
        if (avg_density < avg_density_0):
            lagmin = float(lagrange)
        else:
            lagmax = float(lagrange)
    print("   Dichotomy iterations:", inddico)


def update_exponent(exponent_counter):
    exponent_counter += 1
    if (i < niternp):
        p.assign(Constant(1))
    elif (i >= niternp):
        if i == niternp:
            print("\n Starting penalized iterations\n")
        if ((abs(compliance-old_compliance) < 0.01*compliance_history[0]) and
            (exponent_counter > exponent_update_frequency) ):
            # average gray level
            gray_level = assemble((theta-thetamin)*(1.-theta)*dx)*4/volume
            p.assign(Constant(min(float(p)*(1+0.3**(1.+gray_level/2)), pmax)))
            exponent_counter = 0
            print("   Updated SIMP exponent p = ", float(p))
    return exponent_counter

u = Function(V2, name="Displacement")
old_compliance = 1e30
ffile = XDMFFile("topology_optimization.xdmf")
ffile.parameters["flush_output"]=True
ffile.parameters["functions_share_mesh"]=True
compliance_history = []
for i in range(niter):
    solve(a == L, u, bc, solver_parameters={"linear_solver": "cg", "preconditioner": "hypre_amg"})
    ffile.write(thetaold, i)
    ffile.write(u, i)

    compliance = assemble(action(L, u))
    compliance_history.append(compliance)
    print("Iteration {}: compliance =".format(i), compliance)

    avg_density = update_theta()

    update_lagrange_multiplier(avg_density)

    exponent_counter = update_exponent(exponent_counter)

    # Update theta field and compliance
    thetaold.assign(theta)
    old_compliance = compliance


plot(theta, cmap="bone_r")
plt.title("Final density")
plt.show();

plt.figure()
plt.plot(np.arange(1, niter+1), compliance_history)
ax = plt.gca()
ymax = ax.get_ylim()[1]
plt.plot([niternp, niternp], [0, ymax], "--k")
plt.annotate(r"$\leftarrow$ Penalized iterations $\rightarrow$", xy=[niternp+1, ymax*0.02], fontsize=14)
plt.xlabel("Number of iterations")
plt.ylabel("Compliance")
plt.title("Convergence history", fontsize=16)
plt.show();

zeynep@LAPTOP-3KN21JEE:~$  /usr/bin/env /bin/python3 /home/zeynep/.vscode-server/extensions/ms-python.debugpy-2024.14.0-linux-x64/bundled/libs/debugpy/adapter/../../debugpy/launcher 45387 -- /home/zeynep/FEM/fem1.py 

Solving linear variational problem.
free(): invalid next size (normal)
[LAPTOP-3KN21JEE:37362] *** Process received signal ***
[LAPTOP-3KN21JEE:37362] Signal: Aborted (6)
[LAPTOP-3KN21JEE:37362] Signal code:  (-6)
2

First of all, shouldn’t your variational problem be:

# Inhomogeneous elastic variational problem
du = TestFunction(V2)
u_ = TrialFunction(V2)
a = inner(sigma(u_), eps(du)) * dx
L = dot(f, du) * ds(1)

Secondly:

def left(x, on_boundary):
    return near(x[0], -400) and on_boundary

does not match your grid, as the boundary isn’t as -400 but -920, as indicated by the warning:

Solving linear variational problem.
  *** Warning: Found no facets matching domain for boundary condition.

Please also note that you should further debug your error by removing everything that is not needed to reproduce the error. For instance, you never get past the first solve, so you only need the code:

import meshio
from dolfin import *

msh = meshio.read(filename="Catia2.msh")

for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data

for key in msh.cell_data_dict["gmsh:geometrical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:geometrical"][key]

triangle_mesh = meshio.Mesh(
    points=msh.points,
    cells=[("triangle", triangle_cells)],
    cell_data={"name_to_read": [triangle_data]},
)

meshio.write("test.xdmf", triangle_mesh)


# XDMF
mesh = Mesh()
xdmf_file = XDMFFile("test.xdmf")
xdmf_file.read(mesh)
xdmf_file.close()
# Algorithmic parameters
niter = 80

# Problem parameters
V0 = FunctionSpace(mesh, "DG", 0)
thetamoy = 0.4
E = Constant(1)
nu = Constant(0.3)
lamda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / (2 * (1 + nu))
f = Constant((0, -1, 0))


def left(x, on_boundary):
    return near(x[0], -920) and on_boundary


def load(x, on_boundary):
    return near(x[0], 650) and near(x[1], 150, 50)


# Facets
facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(load).mark(facets, 1)
ds = Measure("ds", subdomain_data=facets)


V2 = VectorFunctionSpace(mesh, "CG", 3)
# Fixed boundary condtions
bc = DirichletBC(V2, Constant((0, 0, 0)), left)

p = Constant(1)  # SIMP penalty exponent
thetaold = Function(V0, name="Density")
thetaold.interpolate(Constant(thetamoy))
coeff = thetaold**p


def eps(v):
    return sym(grad(v))


def sigma(v):
    return coeff * (lamda * div(v) * Identity(3) + 2 * mu * eps(v))


# Inhomogeneous elastic variational problem
du = TestFunction(V2)
u_ = TrialFunction(V2)
a = inner(sigma(u_), eps(du)) * dx
L = dot(f, du) * ds(1)


u = Function(V2, name="Displacement")
solve(
    a == L,
    u,
    bc,
    solver_parameters={"linear_solver": "cg", "preconditioner": "hypre_amg"},
)

image
image2243×1867 505 KB

Secondly, removing your solver parameters:

solve(
    a == L,
    u,
    bc,
    # solver_parameters={"linear_solver": "cg", "preconditioner": "hypre_amg"},
)

lets the code run,meaning that there seems to be an issue with your preconditioner.
Changing it to

solve(
    a == L,
    u,
    bc,
    solver_parameters={"linear_solver": "cg", "preconditioner": "jacobi"},
)

Lets it run, meaning that there is an issue with hypre_amg.

As you are using legacy-FEniCS, i wouldn’t expect anyone putting alot of effort into figuring out what is going wrong within PETSc/DOLFIN interface.

3

Hi! Thank you for your help. I fixed my boundaries and ran the debug, but then I got this error:

Exception has occurred: ValueError
cannot create std::vector larger than max_size()
  File "/home/zeynep/FEM/fem1.py", line 35, in <module>
    xdmf_file.read(mesh)
ValueError: cannot create std::vector larger than max_size()

at this line:
xdmf_file.read(mesh)

4

This must be related to the mesh you are reading in. Are you sure you are reading in the correct mesh?

5

Hello again, I solved this error by providing the full path of the XDMF file. Later, I decided to migrate my code from Dolfin to DolfinX. My mesh file is the same, I just changed the name.

from mpi4py import MPI
import pyvista as pv
import meshio
import numpy as np
import matplotlib.pyplot as plt
import dolfinx
from petsc4py import PETSc
from dolfinx.io import XDMFFile
import dolfinx.fem as fem
import ufl
from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem import FunctionSpace, Function, dirichletbc, assemble,assemble_matrix, assemble_vector,locate_dofs_geometrical, assemble_scalar,form
from ufl import TestFunction, TrialFunction, inner, grad, div, sym, Identity, dot, Measure,Constant
import ufl.finiteelement
import basix
from dolfinx.mesh import locate_entities_boundary, MeshTags
from dolfinx import mesh
from petsc4py.PETSc import ScalarType

msh = meshio.read(filename="/home/zeynep/FEM/square2.msh")

for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data

for key in msh.cell_data_dict["gmsh:geometrical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:geometrical"][key]

triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", triangle_cells)],
                            cell_data={"name_to_read": [triangle_data]})

with XDMFFile(MPI.COMM_WORLD, "/home/zeynep/FEM/square2.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")


# Algorithmic parameters
niternp = 20 
niter = 80 
pmax = 4  
exponent_update_frequency = 4 
tol_mass = 1e-4 
thetamin = 0.001 

thetamoy = 0.4
E = 1
nu =0.3
lamda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / (2 * (1 + nu))
f = np.array([0 ,-1])


def left(x): 
    return np.isclose(x[0], -920)

def load(x): 
    return np.logical_and(np.isclose(x[0], 80), np.isclose(x[1], 470))


#FACETS
facet_dim = mesh.topology.dim - 1
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, facet_dim, left)
left_tags = np.full(len(left_facets), 1, dtype=np.int32)

load_facets = locate_entities_boundary(mesh, facet_dim, load)
load_tags = np.full(len(load_facets), 2, dtype=np.int32)

facet_indices = np.concatenate([left_facets, load_facets])
facet_values = np.concatenate([left_tags, load_tags])

facets =  dolfinx.mesh.meshtags(mesh, facet_dim, facet_indices, facet_values)

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



# Problem parameters
element1 = basix.ufl.element("DG", mesh.ufl_cell().cellname(), degree=0)
V0 = fem.functionspace(mesh, element1)

element2= basix.ufl.element("CG", mesh.ufl_cell().cellname(), degree=2)
V2 = fem.functionspace(mesh, element2)


#DIRICHLETBC
value = 0
boundary_dofs = fem.locate_dofs_geometrical(V2, left)
bc = fem.dirichletbc(ScalarType(value), boundary_dofs, V2)


p = 1 # SIMP penalty exponent

exponent_counter = 0
lagrange = 1 

thetaold = fem.Function(V0, name="Density")
thetaold.x.array[:]=thetamoy
coeff = thetaold**p
theta = Function(V0)

dx = Measure("dx", domain=mesh)
volume = assemble_scalar(form(1.0 * dx))
avg_density_0 = assemble_scalar(form(thetaold * dx)) / volume
avg_density = 0.0


def eps(v):
    return sym(grad(v))

def sigma(v):
    return coeff * (lamda * div(v) * Identity(2) + 2 * mu * eps(v))


# Inhomogeneous elastic variational problem
du = TestFunction(V2)
u_ = TrialFunction(V2)
a = inner(sigma(u_), eps(du)) * dx
L = dot(f, u_) * ds(1)

However, I am getting this error, I would be very grateful if you could assist me with it.

Exception has occurred: ValueError
Symmetric part of tensor with rank != 2 is undefined.
  File "/home/zeynep/FEM/fem1.py", line 139, in eps
    return sym(grad(v))
           ^^^^^^^^^^^^
  File "/home/zeynep/FEM/fem1.py", line 142, in sigma
    return coeff * (lamda * div(v) * Identity(2) + 2 * mu * eps(v))
                                                            ^^^^^^
  File "/home/zeynep/FEM/fem1.py", line 148, in <module>
    a = inner(sigma(u_), eps(du)) * dx
              ^^^^^^^^^
ValueError: Symmetric part of tensor with rank != 2 is undefined.
6

Neither of the spaces here are vector spaces.
See for instance: Introduction to the Unified Form Language — FEniCS Workshop

7

Thank you, i will correct this part, but before there was an another error when transitioning my code to Fenicsx. This part seems to be incorrect because when i run it, a new XDMF file is not created. How should this part be? I would be very grateful if you could help.

msh = meshio.read(filename="/home/zeynep/FEM/square2.msh")

for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data

for key in msh.cell_data_dict["gmsh:geometrical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:geometrical"][key]

triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", triangle_cells)],
                            cell_data={"name_to_read": [triangle_data]})

with XDMFFile(MPI.COMM_WORLD, "/home/zeynep/FEM/square2.xdmf", "w") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
8

You never store a mesh file here, so I am not sure what you would expect.