# Problem with mesh in .xdmf format

**URL:** <https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678>\
**Category:** mesh\
**Created:** [February 14, 2022, 5:06pm UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678 "2022-02-14T17:06:42Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![tfkonopka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tfkonopka/32/2788_2.png) [@tfkonopka](https://fenicsproject.discourse.group/u/tfkonopka)\
**Post date:** [February 14, 2022, 5:06pm UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678/1 "2022-02-14T17:06:42Z")

</div>

Hi everyone,

I have built and tested my code using UnitSquareMesh( ) available in Fenics library and it is working properly. However, I need to work with more complex domains. So, I have chosen gmsh software for this purpose. I have saved the mesh in format .msh and convert it to .xdmf file or converting it directly from .geo to .xdmf .

When I use the mesh input with \*.xdmf files, two error message shows up according to the tests that I have made:

```auto
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0

```

Or this other error message :

```auto
HDF5-DIAG: Error detected in HDF5 (1.10.4) MPI-process 0:
  #000: ../../../src/H5F.c line 444 in H5Fcreate(): unable to create file
    major: File accessibilty
    minor: Unable to open file
  #001: ../../../src/H5Fint.c line 1364 in H5F__create(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #002: ../../../src/H5Fint.c line 1615 in H5F_open(): unable to lock the file
    major: File accessibilty
    minor: Unable to open file
  #003: ../../../src/H5FD.c line 1640 in H5FD_lock(): driver lock request failed
    major: Virtual File Layer
    minor: Can't update object
  #004: ../../../src/H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
    major: File accessibilty
    minor: Bad file ID accessed
HDF5-DIAG: Error detected in HDF5 (1.10.4) MPI-process 0:
  #000: ../../../src/H5F.c line 444 in H5Fcreate(): unable to create file
    major: File accessibilty
    minor: Unable to open file
  #001: ../../../src/H5Fint.c line 1364 in H5F__create(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #002: ../../../src/H5Fint.c line 1615 in H5F_open(): unable to lock the file
    major: File accessibilty
    minor: Unable to open file
  #003: ../../../src/H5FD.c line 1640 in H5FD_lock(): driver lock request failed
    major: Virtual File Layer
    minor: Can't update object
  #004: ../../../src/H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
    major: File accessibilty
    minor: Bad file ID accessed
HDF5-DIAG: Error detected in HDF5 (1.10.4) MPI-process 0:
  #000: ../../../src/H5F.c line 444 in H5Fcreate(): unable to create file
    major: File accessibilty
    minor: Unable to open file
  #001: ../../../src/H5Fint.c line 1364 in H5F__create(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #002: ../../../src/H5Fint.c line 1615 in H5F_open(): unable to lock the file
    major: File accessibilty
    minor: Unable to open file
  #003: ../../../src/H5FD.c line 1640 in H5FD_lock(): driver lock request failed
    major: Virtual File Layer
    minor: Can't update object
  #004: ../../../src/H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
    major: File accessibilty
    minor: Bad file ID accessed

```

I am proving the code with UnitSquareMesh( ) code (working) just for checking if necessary, the code with xdmf input (not working), and the code used to convert from .geo file to .xdmf.

Any suggestions would be very helpful.

xdmf Mesh:

```auto
def lmbdainv(s, mu_rel, no, nw):
    return 1.0 / ((1.0 / mu_rel) * s **nw + (1.0 - s)** no)

# Fractional flow function
def F(s, mu_rel, no, nw):
    return s **nw / (s** nw + mu_rel * (1.0 - s) ** no)

def tensor_jump(v, n):
    return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")

class PiecewiseConstant(UserExpression):
    def __init__ (self, values, markers, **kwargs):
        self._values = values
        self._markers = markers
        super(). __init__ (**kwargs)

    def eval_cell(self, values, x, cell):
        values[0] = self._values[self._markers[cell.index]]

    def value_shape(self):
        return tuple()

### Relative viscosity of water w.r.t. crude oil
Kinv = Constant(1 / (1e-6))
# Time step
dt = Constant(0.01)

mu = 1
mu_b = 1
phi = Constant(0.4)
pin = 6000
pout = 5000
sbar = Constant(1)

t = 0.0
T = 1000 * float(dt)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("domains.xdmf") as infile:
    infile.read(mvc)
Markers = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc1 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("boundaries.xdmf") as infile:
    infile.read(mvc1)
boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc1)

order = 1
BDM = FiniteElement("BDM", mesh.ufl_cell(), order)
DG = FiniteElement("DG", mesh.ufl_cell(), order - 1)
# DG1 = FiniteElement("DG", mesh.ufl_cell(), order)
Element = MixedElement([BDM, DG, DG])
mixed_space = FunctionSpace(mesh, Element)
n = FacetNormal(mesh)

# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)

v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)

s_mid = 0.5 * (s0 + s)

# ============= DEFINITION OF SPATIALLY-VARYING PARAMETERS =====================
marker_inner = 1
marker_outer = 0

no_outer = 2
nw_outer = 2
no_inner = 1
nw_inner = 1

no = {marker_inner: no_inner, marker_outer: no_outer}
nw = {marker_inner: nw_inner, marker_outer: nw_outer}

VVV = FunctionSpace(mesh, "DG", 0)

noo = PiecewiseConstant(no, Markers)
noo_proj = project(noo, VVV)
nww = PiecewiseConstant(nw, Markers)
nww_proj = project(nww, VVV)

# =========== END DEFINITION OF SPATIALLY-VARYING PARAMETERS ===================

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=Markers)

bc2 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 2)

bc4 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 4)

bcs = [bc2, bc4] # velocity BC

alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))

stab = (
    mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)

L1 = (
    mu_b * inner(grad(u), grad(v)) * dx(1)
    + inner(v, lmbdainv(s_mid, mu, no_outer, nw_outer) * Kinv * u) * dx(0)
    - div(v) * p * dx
    + pin * dot(v, n) * ds(1)
    + pout * dot(v, n) * ds(3)
    + stab
)

L2 = q * div(u) * dx(0) + q * div(u) * dx(1)

un = 0.5 * (inner(u0, n) + abs(inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - abs(inner(u0, n)))

stabilisation = (
    dt("+")
    * inner(
        jump(r),
        un("+") * F(s_mid, mu, noo_proj, nww_proj)("+")
        - un("-") * F(s_mid, mu, noo_proj, nww_proj)("-"),
    )
    * dS
)

L3 = (
    phi * r * (s - s0) * dx(0)
    + r * (s - s0) * dx(1)
    - dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(0)
    - dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(1)
    + dt * r * F(s_mid, mu, no_outer, nw_outer) * un * ds
    + stabilisation
    + dt * r * un_h * sbar * ds(1)
)

# Total L
L = L1 + L2 + L3

J = derivative(L, U, dU)

problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)

solver.parameters["newton_solver"]["relative_tolerance"] = 1e-6
solver.parameters["newton_solver"]["absolute_tolerance"] = 1e-6
solver.parameters["newton_solver"]["maximum_iterations"] = 10

u_file = XDMFFile("dir2" + "/velocity.xdmf")
p_file = XDMFFile("dir2" + "/pressure.xdmf")
s_file = XDMFFile("dir2" + "/saturation.xdmf")

while t < T:
    t += float(dt)
    U0.assign(U)
    solver.solve()
    u, p, s = U.split()
    p_file.write(p, t)
    s_file.write(s, t)
    u_file.write(u, t)

```

UnitSquareMesh( ) code:

```auto
# Total mobility
def lmbdainv(s, mu_rel, no, nw):
    return 1.0 / ((1.0 / mu_rel) * s **nw + (1.0 - s)** no)

# Fractional flow function
def F(s, mu_rel, no, nw):
    return s **nw / (s** nw + mu_rel * (1.0 - s) ** no)

def tensor_jump(v, n):
    return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")

class PiecewiseConstant(UserExpression):
    def __init__ (self, values, markers, **kwargs):
        self._values = values
        self._markers = markers
        super(). __init__ (**kwargs)

    def eval_cell(self, values, x, cell):
        values[0] = self._values[self._markers[cell.index]]

    def value_shape(self):
        return tuple()

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.3, 0.7)) and between(x[0], (0.1, 0.5))

# Relative viscosity of water w.r.t. crude oil

Kinv = Constant(1 / (1e-6))
# Time step
dt = Constant(0.01)

mu = 1
mu_b = 1
phi = Constant(0.4)
pin = 6000
pout = 5000
sbar = Constant(1)

t = 0.0
T = 100 * float(dt)

mesh = UnitSquareMesh(90, 90)
n = FacetNormal(mesh)

order = 1
BDM = FiniteElement("BDM", mesh.ufl_cell(), order)
DG = FiniteElement("DG", mesh.ufl_cell(), order - 1)
Element = MixedElement([BDM, DG, DG])
mixed_space = FunctionSpace(mesh, Element)

# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)

v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)

s_mid = 0.5 * (s0 + s)

# ============= DEFINITION OF SPATIALLY-VARYING PARAMETERS =====================
marker_inner = 1
marker_outer = 0

no_outer = 2
nw_outer = 2
no_inner = 1
nw_inner = 1

obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(marker_outer)
obstacle.mark(domains, marker_inner)

no = {marker_inner: no_inner, marker_outer: no_outer}
nw = {marker_inner: nw_inner, marker_outer: nw_outer}

VVV = FunctionSpace(mesh, "DG", 0)

noo = PiecewiseConstant(no, domains)
noo_proj = project(noo, VVV)
nww = PiecewiseConstant(nw, domains)
nww_proj = project(nww, VVV)

# =========== END DEFINITION OF SPATIALLY-VARYING PARAMETERS ===================
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))

# Define boundary markers
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=domains)

bc2 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 2)
# # bc3 = DirichletBC(VQ.sub(0), Constant((0.0, 0.0)), boundaries, 3)
bc4 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 4)

bcs = [bc2, bc4] # velocity BC

alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))

stab = (
    mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)

L1 = (
    mu_b * inner(grad(u), grad(v)) * dx(1)
    + inner(v, lmbdainv(s_mid, mu, no_outer, nw_outer) * Kinv * u) * dx(0)
    - div(v) * p * dx
    + pin * dot(v, n) * ds(1)
    + pout * dot(v, n) * ds(3)
    + stab
)

L2 = q * div(u) * dx(0) + q * div(u) * dx(1)

un = 0.5 * (inner(u0, n) + abs(inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - abs(inner(u0, n)))

stabilisation = (
    dt("+")
    * inner(
        jump(r),
        un("+") * F(s_mid, mu, noo_proj, nww_proj)("+")
        - un("-") * F(s_mid, mu, noo_proj, nww_proj)("-"),
    )
    * dS
)

L3 = (
    phi * r * (s - s0) * dx(0)
    + r * (s - s0) * dx(1)
    - dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(0)
    - dt * inner(grad(r), F(s_mid, mu, no_outer, nw_outer) * u) * dx(1)
    + dt * r * F(s_mid, mu, no_outer, nw_outer) * un * ds
    + stabilisation
    + dt * r * un_h * sbar * ds(1)
)

# Total L
L = L1 + L2 + L3

J = derivative(L, U, dU)

problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)

u_file = XDMFFile("dir2" + "/velocity.xdmf")
p_file = XDMFFile("dir2" + "/pressure.xdmf")
s_file = XDMFFile("dir2" + "/saturation.xdmf")

while t < T:
    t += float(dt)
    U0.assign(U)
    solver.solve()
    u, p, s = U.split()
    p_file.write(p, t)
    s_file.write(s, t)
    u_file.write(u, t)

```

Convert from .geo to .xdmf

```auto
import os
from dolfin import *
import meshio

file_in = "holesuCT"
file_out = "holesuCT"

os.system("gmsh " + file_in + ".geo -2 -o " + file_out + ".msh")
os.system("meshio-convert " + file_out + ".msh " + file_out + ".xdmf")

msh = meshio.read("holesuCT.msh")
tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
line_data = msh.cell_data_dict["gmsh:physical"]["line"]

meshio.write(
    "mesh.xdmf",
    meshio.Mesh(
        points=msh.points[:, :2], cells={"triangle": msh.cells_dict["triangle"]}
    ),
)

meshio.write(
    "boundaries.xdmf",
    meshio.Mesh(
        points=msh.points,
        cells={"line": msh.cells_dict["line"]},
        cell_data={"bnd_marker": [msh.cell_data_dict["gmsh:physical"]["line"]]},
    ),
)

meshio.write(
    "domains.xdmf",
    meshio.Mesh(
        points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
        cell_data={"dom_marker": [msh.get_cell_data("gmsh:physical", "triangle")]},
    ),
)

```

---

<div class="post-metadata">

**Author:** ![Kei\_Yamamoto](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kei_yamamoto/32/2904_2.png) [@Kei\_Yamamoto](https://fenicsproject.discourse.group/u/Kei_Yamamoto)\
**Post date:** [February 14, 2022, 6:08pm UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678/2 "2022-02-14T18:08:48Z")

</div>

Hi,

I’m not sure what you want to do with `mesh.xdmf` and `domains.xdmf` as you typically need one mesh for the entire domain and the other mesh for the boundary part (`boundaries.xdmf` in your case)

Have you looked this [tutorial](http://jsdokken.com/converted_files/tutorial_pygmsh.html#second) how to conver `.msh` into `.xdmf ` ??

---

<div class="post-metadata">

**Author:** ![tfkonopka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tfkonopka/32/2788_2.png) [@tfkonopka](https://fenicsproject.discourse.group/u/tfkonopka)\
**Post date:** [February 14, 2022, 9:36pm UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678/3 "2022-02-14T21:36:52Z")

</div>

> [@Kei\_Yamamoto](#):
>
> the entire domain and the other mesh for the boundary part ( `boundaries.xdmf` in your case)

Dear @Kei_Yamamoto ,

I need the " domains.xdmf " because in this problem I have different properties along the domain. So, I need to mark the different domains. For example, like the figure bellow. Here I need to calculate the function `def F(s, mu_rel, no, nw)` using no = nw = 2 in blue domain and no = nw = 1 in red one.  
Moreover, I will construct much more complex domains than this one. So, this is the reason that I am using gmsh.

![16_vuggs](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/3/3813a34cde5f37a1ade196fe6b16fc77e3536093.png)

Sorry about not mentioning this figure previously.

best regards  
Thiago Konopka

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 15, 2022, 9:00am UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678/4 "2022-02-15T09:00:44Z")

</div>

The domains can be saved in the saved file as the mesh, as both of them relate to cells.

Secondly, you have not supplied the code for generating the mesh, and thus it is very hard to give you any further guidance.

---

<div class="post-metadata">

**Author:** ![tfkonopka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tfkonopka/32/2788_2.png) [@tfkonopka](https://fenicsproject.discourse.group/u/tfkonopka)\
**Post date:** [February 15, 2022, 2:21pm UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678/5 "2022-02-15T14:21:35Z")

</div>

Dear @dokken,

I did not know that I could use just one file for mesh and domain. I base this choice on other examples in this forum.

The .geo file bellow represents a very simple case like the figure bellow:

Note: this .geo code is generated from an image processing code for complex geometries with multiple subdomains like the red square.

Best regards  
Thiago Konopka

![vugg_forum](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/7e6d3479fb240db067fbff0b4a74051d62e5b248.png)

```auto
// mesh size description

cl_1 = 0.1;
cl_2 = 0.05;

// boundary points
Point(1) = {0, 0, 0, cl_1};
Point(2) = {1.000, 0, 0, cl_1};
Point(3) = {1.000, 1.000, 0, cl_1};
Point(4) = {0, 1.000, 0, cl_1};

// lines that connect boundary
Line(1) = {1, 2};
Line(2) = {3, 2};
Line(3) = {4, 3};
Line(4) = {1, 4};

  // Mesh Parameters
  Mesh.CharacteristicLengthExtendFromBoundary = 0;
  Mesh.CharacteristicLengthMax = 0.01;
 
  // Define Segment coordinates

// Hole 2

X2 = {0.300, 0.304, 0.308, 0.312, 0.316, 0.320, 0.324, 0.328, 0.332, 0.336, 0.340, 0.344, 0.348, 0.351, 0.355, 0.359, 0.363, 0.367, 0.371, 0.375, 0.379, 0.383, 0.387, 0.391, 0.395, 0.399, 0.403, 0.407, 0.411, 0.415, 0.419, 0.423, 0.427, 0.431, 0.435, 0.439, 0.443, 0.447, 0.450, 0.454, 0.458, 0.462, 0.466, 0.470, 0.474, 0.478, 0.482, 0.486, 0.490, 0.494, 0.498, 0.502, 0.506, 0.510, 0.514, 0.518, 0.522, 0.526, 0.530, 0.534, 0.538, 0.542, 0.546, 0.550, 0.553, 0.557, 0.561, 0.565, 0.569, 0.573, 0.577, 0.581, 0.585, 0.589, 0.593, 0.597, 0.601, 0.605, 0.609, 0.613, 0.617, 0.621, 0.625, 0.629, 0.633, 0.637, 0.641, 0.645, 0.649, 0.652, 0.656, 0.660, 0.664, 0.668, 0.672, 0.676, 0.680, 0.684, 0.688, 0.692, 0.696, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.696, 0.692, 0.688, 0.684, 0.680, 0.676, 0.672, 0.668, 0.664, 0.660, 0.656, 0.652, 0.649, 0.645, 0.641, 0.637, 0.633, 0.629, 0.625, 0.621, 0.617, 0.613, 0.609, 0.605, 0.601, 0.597, 0.593, 0.589, 0.585, 0.581, 0.577, 0.573, 0.569, 0.565, 0.561, 0.557, 0.553, 0.550, 0.546, 0.542, 0.538, 0.534, 0.530, 0.526, 0.522, 0.518, 0.514, 0.510, 0.506, 0.502, 0.498, 0.494, 0.490, 0.486, 0.482, 0.478, 0.474, 0.470, 0.466, 0.462, 0.458, 0.454, 0.450, 0.447, 0.443, 0.439, 0.435, 0.431, 0.427, 0.423, 0.419, 0.415, 0.411, 0.407, 0.403, 0.399, 0.395, 0.391, 0.387, 0.383, 0.379, 0.375, 0.371, 0.367, 0.363, 0.359, 0.355, 0.351, 0.348, 0.344, 0.340, 0.336, 0.332, 0.328, 0.324, 0.320, 0.316, 0.312, 0.308, 0.304, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300}; 
Y2 = {0.3, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.304, 0.308, 0.312, 0.316, 0.320, 0.324, 0.328, 0.332, 0.336, 0.340, 0.344, 0.348, 0.351, 0.355, 0.359, 0.363, 0.367, 0.371, 0.375, 0.379, 0.383, 0.387, 0.391, 0.395, 0.399, 0.403, 0.407, 0.411, 0.415, 0.419, 0.423, 0.427, 0.431, 0.435, 0.439, 0.443, 0.447, 0.450, 0.454, 0.458, 0.462, 0.466, 0.470, 0.474, 0.478, 0.482, 0.486, 0.490, 0.494, 0.498, 0.502, 0.506, 0.510, 0.514, 0.518, 0.522, 0.526, 0.530, 0.534, 0.538, 0.542, 0.546, 0.550, 0.553, 0.557, 0.561, 0.565, 0.569, 0.573, 0.577, 0.581, 0.585, 0.589, 0.593, 0.597, 0.601, 0.605, 0.609, 0.613, 0.617, 0.621, 0.625, 0.629, 0.633, 0.637, 0.641, 0.645, 0.649, 0.652, 0.656, 0.660, 0.664, 0.668, 0.672, 0.676, 0.680, 0.684, 0.688, 0.692, 0.696, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.700, 0.696, 0.692, 0.688, 0.684, 0.680, 0.676, 0.672, 0.668, 0.664, 0.660, 0.656, 0.652, 0.649, 0.645, 0.641, 0.637, 0.633, 0.629, 0.625, 0.621, 0.617, 0.613, 0.609, 0.605, 0.601, 0.597, 0.593, 0.589, 0.585, 0.581, 0.577, 0.573, 0.569, 0.565, 0.561, 0.557, 0.553, 0.550, 0.546, 0.542, 0.538, 0.534, 0.530, 0.526, 0.522, 0.518, 0.514, 0.510, 0.506, 0.502, 0.498, 0.494, 0.490, 0.486, 0.482, 0.478, 0.474, 0.470, 0.466, 0.462, 0.458, 0.454, 0.450, 0.447, 0.443, 0.439, 0.435, 0.431, 0.427, 0.423, 0.419, 0.415, 0.411, 0.407, 0.403, 0.399, 0.395, 0.391, 0.387, 0.383, 0.379, 0.375, 0.371, 0.367, 0.363, 0.359, 0.355, 0.351, 0.348, 0.344, 0.340, 0.336, 0.332, 0.328, 0.324, 0.320, 0.316, 0.312, 0.308, 0.304, 0.300};

// Define spline surface

LN = 90;

nR = #X2[];
p0 = newp;
p = p0;
For i In {0:nR-1}
Point(newp) = {X2[i], Y2[i], 0, cl_2};
EndFor
p2 = newp-1;
BSpline(90) = {p:p2,p};
Line Loop(91) = {90};
Plane Surface(92) = {91};

// Define all surfaces
Line Loop(5) = {1, -2, -3, -4};

Physical Line(1) = {4};
Physical Line(2) = {3};
Physical Line(3) = {2};
Physical Line(4) = {1};

Plane Surface(93) ={5,91};
Physical Surface(0) = {93};
Physical Surface(1) ={92};

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 16, 2022, 9:25am UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678/6 "2022-02-16T09:25:25Z")

</div>

I reverted the post to contain the original question.  
Solution was to change the initial guess:

Adding a non-zero initial guess to your non-linear problem seems to help:

```auto
U.vector()[:] = 0.00001

```

prior to

> [@tfkonopka](#):
>
> ````auto
> > problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
> solver = NonlinearVariationalSolver(problem)
> > ```
> 
> ````

makes the solution converge

---

<div class="post-metadata">

**Author:** ![tfkonopka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tfkonopka/32/2788_2.png) [@tfkonopka](https://fenicsproject.discourse.group/u/tfkonopka)\
**Post date:** [February 16, 2022, 12:38pm UTC](https://fenicsproject.discourse.group/t/problem-with-mesh-in-xdmf-format/7678/7 "2022-02-16T12:38:11Z")

</div>

It is working @dokken,

Thank you