Problems with facet_markers assignment

dolfinx
,
1

Good afternoon, I have the following error.
I am generating a mesh, the problem arises when I want to define a facet_markers to apply a load on that line, at first I tried to Save it along with the rest with the command ‘meshio.write’, however, I got the following error:

ValueError                                Traceback (most recent call last)
<ipython-input-193-d7866242edc3> in <cell line: 102>()
    100 
    101 # Crear el archivo XDMF incluyendo tanto los elementos como las facetas
--> 102 meshio.write("presa_suelo.xdmf", meshio.Mesh(points=points_2d, cells=cells, cell_data=cell_data))

/usr/local/lib/python3.10/dist-packages/meshio/_mesh.py in __init__(self, points, cells, point_data, cell_data, field_data, point_sets, cell_sets, gmsh_periodic, info)
    172         for key, data in self.cell_data.items():
    173             if len(data) != len(cells):
--> 174                 raise ValueError(
    175                     f"Incompatible cell data '{key}'. "
    176                     f"{len(cells)} cell blocks, but '{key}' has {len(data)} blocks."

ValueError: Incompatible cell data 'triangle'. 2 cell blocks, but 'triangle' has 1 blocks.

I tried everything, the problem is only removed when I do not generate physical groups of dimension 1 (Line). Since said physical group is necessary, I don’t know how to proceed, I would appreciate your help, I leave the code in this part:

import gmsh
import meshio
import numpy as np

# Inicializar Gmsh
gmsh.initialize()
gmsh.model.add("Presa y Suelo")

# Parámetros de la presa y el suelo
altura_triangulo = 42  # metros
base_triangulo = 34  # metros
base_rectangulo = 40.0  # metros
altura_rectangulo = 8.0  # metros
profundidad_suelo = 2 * (altura_triangulo + altura_rectangulo)  # Dos veces la altura total de la presa
ancho_suelo = 2 * base_rectangulo  # El doble de la base rectangular

# Definición de puntos de la presa
A = gmsh.model.geo.addPoint(0, 0, 0)
B = gmsh.model.geo.addPoint(base_rectangulo, 0, 0)
C = gmsh.model.geo.addPoint(base_rectangulo, altura_rectangulo, 0)
D = gmsh.model.geo.addPoint(0, altura_rectangulo, 0)
E = gmsh.model.geo.addPoint(base_triangulo, altura_rectangulo, 0)
F = gmsh.model.geo.addPoint(0, altura_triangulo + altura_rectangulo, 0)

# Definición de líneas de la presa
l1 = gmsh.model.geo.addLine(A, B)
l2 = gmsh.model.geo.addLine(B, C)
l3 = gmsh.model.geo.addLine(C, E)
l4 = gmsh.model.geo.addLine(E, F)
l5 = gmsh.model.geo.addLine(F, D)
l6 = gmsh.model.geo.addLine(D, A)

# Definición de superficies de la presa
loop_presa = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4, l5, l6])
superficie_presa = gmsh.model.geo.addPlaneSurface([loop_presa])

# Definición del suelo
P1 = gmsh.model.geo.addPoint(-ancho_suelo, -profundidad_suelo, 0)
P2 = gmsh.model.geo.addPoint(base_rectangulo + ancho_suelo, -profundidad_suelo, 0)
P3 = gmsh.model.geo.addPoint(base_rectangulo + ancho_suelo, 0, 0)
P4 = gmsh.model.geo.addPoint(-ancho_suelo, 0, 0)

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

l12 = gmsh.model.geo.addLine(F,A)  # Línea para la presión del agua

# Corregir el sentido de las líneas para la definición correcta de la superficie del suelo
loop_suelo = gmsh.model.geo.addCurveLoop([l7, l8, l9, -l1, l10, l11])
superficie_suelo = gmsh.model.geo.addPlaneSurface([loop_suelo])

# Sincronización de la geometría
gmsh.model.geo.synchronize()

# Definición de grupos físicos para aplicar cargas y condiciones de frontera
gmsh.model.addPhysicalGroup(2, [superficie_presa], 1)  # Presa
gmsh.model.addPhysicalGroup(2, [superficie_suelo], 2)  # Suelo

## Grupos físicos adicionales SÓLO AL COMETAR ÉSTA LÍNEA EL CÓDIGO FUNCIONA, SIN EMBARGO, NECESITO EL GRUPO FÍSICO PARA LA ASIGNACIÓN CORRECTA DE LA CARGA HIDROSTÁTICA
gmsh.model.addPhysicalGroup(1, [l12], 3)  # Cara izquierda de la presa (para presión del agua)

# Forzar elementos triangulares en la malla
gmsh.option.setNumber("Mesh.ElementOrder", 1)
gmsh.option.setNumber("Mesh.Algorithm", 6)  # Tipo de Mallado
gmsh.option.setNumber("Mesh.RecombineAll", 0)

# Ajuste del tamaño característico
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 5)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 10)

# Genera la malla 2D
gmsh.model.mesh.generate(2)

# Exportar la malla a formato .msh
gmsh.write("presa_suelo.msh")

# Finalizar Gmsh
gmsh.finalize()

# Leer la malla de Gmsh utilizando Meshio
msh = meshio.read("presa_suelo.msh")

# Convertir puntos a 2D eliminando la coordenada z
points_2d = msh.points[:, :2]

# Obtener celdas triangulares y líneas
triangle_cells = msh.cells_dict["triangle"]
cells = [("triangle", triangle_cells)]
cell_data = {"triangle": [msh.cell_data_dict["gmsh:physical"]["triangle"]]}

# Verificar y agregar las facetas (líneas) si existen
if "line" in msh.cells_dict:
    line_cells = msh.cells_dict["line"]
    cells.append(("line", line_cells))
    line_data = msh.cell_data_dict["gmsh:physical"]["line"]
    cell_data["line"] = [line_data]

# Crear el archivo XDMF incluyendo tanto los elementos como las facetas
meshio.write("presa_suelo.xdmf", meshio.Mesh(points=points_2d, cells=cells, cell_data=cell_data))

#HYDROSTATIC LOAD

ds_water = Measure('ds', domain=mesh, subdomain_data=facet_markers, subdomain_id=3)
L += presion_hidrostatica * inner(Constant((0, 1)), v) * ds_water

Thanks in advance