 # Mesh, Submesh , Torus, Boundary

Hello,

I want to solve the following problem.

\Omega is unit cube, T is torus with T \subset \Omega and T is supposed to be in the middle of \Omega.

find (u,p) \in V \times S, s.t.
\begin{equation} \begin{cases} curl(\alpha(\left| curl u \right| )curl u)+\nabla p = g \text{ in } \Omega \\ div u=0 \text{ in } \Omega \\ u \times n = 0 \text{ auf } \partial\Omega. \end{cases} \end{equation}

with the weak formulation

\begin{equation}\int_{\Omega} \alpha(\left|curlu(x) \right|)curlu(x) \cdot curlv(x) dx + \int_{\Omega} v (x)\cdot \nabla p(x) dx = \int_{\Omega} g \cdot v dx \end{equation}
\begin{equation} \int_{\Omega} u(x)\cdot \nabla q(x) dx = 0 \end{equation}

where V is the finite elemente space with Nedelec elements and S is the space with piecewise linear elements.
u, p should exist in complete \Omega
g= \chi_T *(\frac{z}{\sqrt{x^2+z^2}},0,\frac{x}{\sqrt{x^2+z^2}}) ^T
\begin{align}\alpha (s) = \begin{cases} 0.001 \,on \,T \\ 1- \frac{1}{2(s^2+1)}\end{cases} on \Omega/T \end{align}.

My problem is to create the domains and subdomains respectively the meshes so that I can integrate over the subdomains :
\int_{\Omega} g \cdot v dx = \int_{T} (\frac{z}{\sqrt{x^2+z^2}},0,\frac{x}{\sqrt{x^2+z^2}}) ^T \cdot v dx
and

\int_{\Omega} \alpha(\left|curlu(x) \right|)curlu(x) \cdot curlv(x) dx = \int_{\Omega/T} (1- \frac{1}{2(|curlu(x) |^2+1)})curlu(x) \cdot curlv(x) dx + \int_{T} 0.001curlu(x) \cdot curlv(x) dx

My biggest Problem is how to define T and \Omega / T and then the boundary conditions.

The interior of the torus including the boundary satisfies
\left( \sqrt{x^2 + y^2} - R \right)^2 + z^2 \leq r^2
so for everything outside the torus
\left( \sqrt{x^2 + y^2} - R \right)^2 + z^2 \ge r^2
or? It should not only be outside the torus but also inside the cube. So x,y,z \in [-1,1]^3.

With submesh and thanks to dokken : Define a Torus in the UnitCubemesh I tried

from fenics import*
import numpy as np

# -------------------------------------------------
#Mesh Omega
N = 50
mesh = BoxMesh(Point(-1, -1, -1), Point(1, 1, 1), N, N, N)
print('dim mesh:', mesh.geometry().dim())  # Dimension mesh

# -------------------------------------------------
#Torus
r = Constant(0.1)
R = Constant(0.2)
class Torus(SubDomain):
def inside(self, x, on_boundary):
return (np.sqrt(x * x + x * x) - float(R)) ** 2 + x * x <= float(r) ** 2

torus = Torus()

# -------------------------------------------------
mft = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
mft.set_all(0)
torus.mark(mft, 1)
File("mft.pvd") << mft
dxT = Measure('dx', domain=mesh, subdomain_data=mft)

submesh_torus = SubMesh(mesh, mft, 1)
print('dim submesh T:', submesh_torus.geometry().dim())  # Dimension mesh

# -------------------------------------------------
# Omega/T
class Omega_Torus(SubDomain):
def inside(self, x, on_boundary):
return (np.sqrt(x * x + x * x) - float(R)) ** 2 + x * x >= float(r) ** 2 and -1 <= x <= 1 and -1 <= x <= 1 and -1 <= x <= 1

omega_torus = Omega_Torus()
# Torus T
mf = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
mf.set_all(0)
omega_torus.mark(mf, 1)
File("mf.pvd") << mf

dxOT = Measure('dx', domain=mesh, subdomain_data=mf)
submesh_OT = SubMesh(mesh, mf, 1)
print('dim submesh OT:', submesh_OT.geometry().dim())  # Dimension mesh


Is it possible to do it this way? When I look at the meshes in paraview, they look matching.

But how I have to create the boundary conditions?

Thank you and
best regards,
noya

Hello everybody,
I’m trying to set up the boundary conditions, but don’t know exactly how.

I currently have

from fenics import*
import numpy as np

# -------------------------------------------------
#Mesh Omega
N = 50
mesh = BoxMesh(Point(-1, -1, -1), Point(1, 1, 1), N, N, N)

# -------------------------------------------------
#Torus
r = Constant(0.1)
R = Constant(0.2)
class Torus(SubDomain):
def inside(self, x, on_boundary):
return (np.sqrt(x * x + x * x) - float(R)) ** 2 + x * x <= float(r) ** 2

torus = Torus()

# -------------------------------------------------
# Meshfunction Torus T
mft = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
mft.set_all(0)
torus.mark(mft, 1)
File("mft.pvd") << mft
dxT = Measure('dx', domain=mesh, subdomain_data=mft)

submesh_torus = SubMesh(mesh, mft, 1)

# -------------------------------------------------
# Omega/T
class Omega_Torus(SubDomain):
def inside(self, x, on_boundary):
return (np.sqrt(x * x + x * x) - float(R)) ** 2 + x * x > float(r) ** 2 and -1 <= x <= 1 and -1 <= x <= 1 and -1 <= x <= 1

omega_torus = Omega_Torus()
# Meshfunction Omega/T
mf = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
mf.set_all(0)
omega_torus.mark(mf, 1)
File("mf.pvd") << mf

dxOT = Measure('dx', domain=mesh, subdomain_data=mf)
submesh_OT = SubMesh(mesh, mf, 1)

# -------------------------------------------------

V = FiniteElement("N1curl", mesh.ufl_cell(), 1)
S = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Z = MixedElement([V, S])
W = FunctionSpace(mesh, Z)
W1 = FunctionSpace(submesh_OT, Z)
T = FunctionSpace(submesh_torus, S)

#Boundary \Omega
def O_boundary(x, on_boundary):
return on_boundary
w0 = Function(W)
bc0 = DirichletBC(W, w0, O_boundary)
#Boundary \Omega / T
def OT_boundary(x, on_boundary):
return on_boundary
w1 = Function(W1)
bc1 = DirichletBC(W, w1, OT_boundary)
#Boundary T
def T_boundary(x, on_boundary):
return on_boundary
w2 = Function(T)
bc2 = DirichletBC(W, w2, T_boundary)

bc = [bc0, bc1, bc2]


Is this correct?

Hi @Noya,

You may want to try creating your mesh using an external meshing tool, e.g. Gmsh. The advantage of this approach is that you can mark the subdomains and boundaries using Gmsh, save the markers to files, and load them into FEniCS MeshFunctions, which can be used to define integration measures and boundary conditions. Using your current code (i.e. with BoxMesh), you obtain a mesh that doesn’t respect the internal boundary of the torus, i.e. you will have elements that pass through the torus boundary.

Below is a code using the Gmsh Python API to generate the mesh, and meshio to convert to XDMF files. Running this code produces three files: mesh.xdmf contains only the mesh, domains.xdmf contains markers which identify the domains \Omega and T, and boundaries.xdmf contains markers that identify the boundaries \partial \Omega and \partial T.

import gmsh
import meshio
import numpy as np

gmsh.initialize()

R = 0.5  # Major radius
r = 0.1  # Minor radius
L = 2    # Box dimension

hmin = r/2   # Min mesh size
hmax = L/10  # Max mesh size

# Create geometry
gmsh.model.occ.rotate([(3,2)],0,0,0,0,0,1,np.pi)
gmsh.model.occ.fragment([(3,3)], [(3,1),(3,2)])
gmsh.model.occ.synchronize()

# Create mesh fields to control element size
gmsh.model.mesh.field.setString(d, "F", "Sqrt((Sqrt(x^2+y^2)-{R:f})^2+z^2)".format(R=R))
gmsh.model.mesh.field.setNumber(t, "InField", d)
gmsh.model.mesh.field.setNumber(t, "DistMin", r)
gmsh.model.mesh.field.setNumber(t, "DistMax", 2*r)
gmsh.model.mesh.field.setNumber(t, "SizeMin", hmin)
gmsh.model.mesh.field.setNumber(t, "SizeMax", hmax)
gmsh.model.mesh.field.setAsBackgroundMesh(t)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)

# Generate physical groups to mark domains and boundaries
gmsh.model.addPhysicalGroup(2, [7,8], 3)         # Internal boundary

# Generate mesh and save
gmsh.model.mesh.generate()
gmsh.write("torus_mesh.msh")
gmsh.finalize()

# Convert mesh to XDMF for use in FEniCS
tet_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]

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


You can then load the markers contained in domains.xdmf and boundaries.xdmf, e.g.:

from fenics import*
import numpy as np

# -------------------------------------------------
#Mesh Omega
N = 50
mesh = Mesh()
with XDMFFile("mesh.xdmf") as mshfile:

# -------------------------------------------------
# Meshfunction over domains
mft = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
with XDMFFile("domains.xdmf") as domfile:
File("mft.pvd") << mft
dxT = Measure('dx', domain=mesh, subdomain_data=mft)

submesh_torus = SubMesh(mesh, mft, 1)
submesh_OT = SubMesh(mesh, mft, 2)

# -------------------------------------------------

V = FiniteElement("N1curl", mesh.ufl_cell(), 1)
S = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Z = MixedElement([V, S])
W = FunctionSpace(mesh, Z)
W1 = FunctionSpace(submesh_OT, Z)
T = FunctionSpace(submesh_torus, S)

# -------------------------------------------------
# Meshfunction over facets
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("boundaries.xdmf") as bndfile:
mff = MeshFunction('size_t', mesh, mvc)
for i in range(mff.size()):
if mff[i] not in (3,4):
mff[i] = 0
File("mff.pvd") << mff
dS = Measure('dS', domain=mesh, subdomain_data=mff)

#Boundary \Omega
w0 = Function(W)
bc0 = DirichletBC(W, w0, mff, 4)
#Boundary T
w2 = Function(T)
bc2 = DirichletBC(W, w0, mff, 3)

bc = [bc0, bc2]


With a suitable set of filters in Paraview, I obtain the following for mff.pvd: and the following for mft.pvd: 2 Likes

WOW. Thank you.
Last week I tried gmsh, but my supervisor tells me to do it with submesh. But your work with gmsh is so much nicer. Wow.
I have some questions to understand everything.

# Create geometry
gmsh.model.occ.rotate([(3,2)],0,0,0,0,0,1,np.pi)
gmsh.model.occ.fragment([(3,3)], [(3,1),(3,2)])
gmsh.model.occ.synchronize()


why are two tori added here resp. why two half tori? and not one complete torus?
And what is the “fragment”?

N = 50
mesh = Mesh()
with XDMFFile("mesh.xdmf") as mshfile:


Does I need N? Or where I need it?

Actually I’m a little bit confused about mff and mft.
Can you please explain the difference? And what both do exactly?

dxT = Measure('dx', domain=mesh, subdomain_data=mft)


So I can calculate the integral \int_T f(x) dx for all x in T?

dS = Measure('dS', domain=mesh, subdomain_data=mff)


What does dS exactly do?

And how can I calculate the integral over \Omega/T?

I also have a question about paraview.
Which filters do you use there?

Thank you so much.

Fragment means that it takes multiple geometries, merge them, but keep the boundary between them. For more information see: Gmsh 4.8.4

You should no need N any longer.

mft is the marked cells (tetrahedrons) of the torus and surrounding box, while mff is the marked facets of the boundary of the torus (with value 3) and outer boundary (with value 4).

if you want to compute the integral over the torus, you can use:

dx = Measure("dx", domain=mesh, subdomain_data=mft)
print(assemble(Constant(1)*dx(1)))


which should give you the volume of the torus (note that dx(1) indicates that we are integrating all cells tagged in mft with value 1).
Similarly, for the box surrounding the torus (marked with 2), you would get
print(assemble(Constant(1)*dx(2)))

1 Like

Okay. Thank you.

so mft represents the entire area including the torus T?
And mff with value 3 represents the boundary of the torus and with value 4 the boundary of the cube?

I want to create i.e the integral \int_T f*v*dx

where f = Expression(('-x/sqrt(x*x+x*x)', '0.0', 'x/sqrt(x*x+x*x)'), degree=2) and vh = TestFunction(W)

Is this

L = inner(f, vh) * dx(1)


? And if I take dx(2) I integrate over \Omega/T and for dx I integrate over \Omega?

Yes. I would suggest you sit down and try some simple integrals (Such as the one I suggest), where you can verify the result, before writing your whole code.

Thank you all so much. Now I will try to play with the different integrals and learn more about it. I found that Gmsh was unable to mesh the geometry when using one complete torus. The error message mentioned that Gmsh was unable to create a periodic boundary mesh. I’m not exactly certain why this error is encountered, but splitting the torus into two half tori solved the issue. My guess is that it makes the topology of the toroidal surface more amenable to meshing.

Hello all,
I want the torus to be in the x-z plane. And now I have defined the structure as follows in a geo-file. And my

SetFactory("OpenCASCADE");

L = 2;
R = 0.2;
r = 0.1;
hmin = r/2;
hmax = L/10;

Box(1) = {-L/2,-L/2,-L/2,L,L,L};  //tag=1
Torus(2) = {0, 0, 0, R, r, 2*Pi};  //tag=2
Rotate {{1, 0, 0}, {0, 0, 0}, Pi/2}{
Volume{2};
}
v() = BooleanFragments{ Volume{1, 2}; Delete; }{};
Physical Volume(22) = {v};
Physical Volume(23) = {v};
Physical Surface(3) = {7};    //innerer Rand
Physical Surface(4) = {8,9,10,11,12,13};    //aeusserer Rand
BooleanUnion{ Volume{2}; Delete; }{ Volume{3}; Delete; }
Field = MathEval;
Field.F = "Sqrt((Sqrt(x^2+y^2)-{R})^2+z^2)";

Field = Threshold;
Field.DistMax = 2*r;
Field.DistMin = r;
Field.SizeMin = hmin;
Field.SizeMax = hmax;
Field.InField = Field;
Background Field = 2;
Mesh 3;


from fenics import*
import meshio

tet_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
tri_data = msh.cell_data_dict["gmsh:physical"]["triangle"]

meshio.write("meshtorus.xdmf",
meshio.Mesh(points=msh.points,
cells={"tetra": msh.cells_dict["tetra"]}
)
)
meshio.write("randtorus.xdmf",
meshio.Mesh(points=msh.points,
cells={"triangle": msh.cells_dict["triangle"]},
cell_data={"bnd_marker": [msh.cell_data_dict["gmsh:physical"]["triangle"]]}
)
)
meshio.write("gebiettorus.xdmf",
meshio.Mesh(
points=msh.points,
cells={"tetra": msh.cells_dict["tetra"]},
cell_data={"dom_marker": [msh.get_cell_data("gmsh:physical", "tetra")]}
)
)

mesh = Mesh()
with XDMFFile("meshtorus.xdmf") as mshfile:

mft = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
with XDMFFile("gebiettorus.xdmf") as domfile:
File("mft.pvd") << mft
dx = Measure('dx', domain=mesh, subdomain_data=mft)

submesh_torus = SubMesh(mesh, mft, 22)
submesh_OT = SubMesh(mesh, mft, 23)

V = FiniteElement("N1curl", mesh.ufl_cell(), 1)
S = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Z = MixedElement([V, S])
W = FunctionSpace(mesh, Z)
T = FunctionSpace(submesh_torus, Z)
OT = FunctionSpace(submesh_OT, Z)
print('Dim = DoF: ', W.dim())  # Dimension W

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("randtorus.xdmf") as bndfile:
mff = MeshFunction('size_t', mesh, mvc)
for i in range(mff.size()):
if mff[i] not in (3, 4):
mff[i] = 0
File("mff.pvd") << mff
dS = Measure('dS', domain=mesh, subdomain_data=mff)
ds = Measure("ds", domain=mesh, subdomain_data=mff)
print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(22)))
print(assemble(Constant(1)*dx(23)))
print(assemble(Constant(1)*ds(4)))
print(assemble(Constant(1)*dS(3)))

w0 = Function(W)
bc0 = DirichletBC(W, w0, mff, 4)

w1 = Function(T)
bc1 = DirichletBC(W, w1, mff, 3)

bc = [bc0, bc1]

(vh, qh) = TestFunctions(W)
(yh, ph) = TrialFunctions(W)

z0 = Function(W)
y0, p0 = split(z0)

s = inner(curl(y0), curl(y0))
nu_cube = 1 - (1 / (2 * (s + 1)))
nu_torus = 1.0000064

f = Expression(('-x/sqrt(x*x+x*x)', '0', 'x/sqrt(x*x+x*x)'), degree=2)

# weak formulation
a1 = nu_cube * inner(curl(yh), curl(vh)) * dx(23)
a2 = nu_torus * inner(curl(yh), curl(vh)) * dx(22)
a = a1 + a2
b = inner(vh, grad(ph)) * dx
c = inner(yh, grad(qh)) * dx
# A = a + b + c
A = a + b + c
L = inner(f, vh) * dx(22)

maxiter = 1000
tol = 1e-15

m = Function(W)
for n in range(1, maxiter):
solve(A == L, m, bc)
yh, ph = m.split()
error1 = errornorm(yh, z0.sub(0), 'Hcurl')
error2 = errornorm(ph, z0.sub(1), 'H10')
error = error1 + error2
print(error1)
print(error2)
print(n, error)
if error < tol:
break
z0.assign(m)

yh.rename("yh", "")
ph.rename("ph", "")

with XDMFFile(MPI.comm_world, 'lsgtorus.xdmf') as file:
file.parameters.update(
{
"functions_share_mesh": True,
"rewrite_function_mesh": False
})
file.write(yh, 0)
file.write(ph, 0)


But it seems that the mft is not correct or not as above (s.Mesh, Submesh , Torus, Boundary - #3 by conpierce8)
And the error between p^1-p^0 is 0 and the total error is already 0 in the second step. Furthermore I can’t present p in the visualization of the solution with paraview. Something is still going wrong, but I don’t know exactly what. Could someone give me a hint?
The complete output is

Dim = DoF: 8303
8.000000000000014
0.01920222306126237
7.9807977769387515
23.999999999999957
0.6612439999362151
Solving linear variational problem.
Building point search tree to accelerate distance queries.
Computed bounding box tree with 65 nodes for 33 points.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
0.3431120777953642
0.0
1 0.3431120777953642
Solving linear variational problem.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
0.0
0.0
2 0.0
Process finished with exit code 0

Best regards, noya

Hi @Noya,

Are you sure your file is meshing correctly? When I open your file with the Gmsh GUI on my machine, I get a “Segmentation fault (core dumped)” message in the command line, and the GUI immediately closes.

When I comment out line 19 (the BooleanUnion), the GUI opens and I see the error:

• line 28: syntax error (Field)

and the only elements appear to be triangles on the outer surface of the box. The error on line 28 can be removed by using the syntax found at Field in Gmsh t10.

When you fix that error, you may receive another error (invalid token on expression), because the syntax {R} in your expression for Field is not correct. You should use Sprintf (see Field in Gmsh t10) to insert the value of R into your expression.

When you fix these errors, Gmsh should be able to fully mesh your geometry. As a final note, the expression in Field should be modified for a torus in the xz-plane:

\sqrt{ \left(\sqrt{x^2 + z^2}^2 - R\right)^2 + y^2}
2 Likes

Hey @conpierce8 ,

thank you. I thought it was meshed correctly since I didn’t get an error message and I could see both the 3D and 2D elements.
I have now changed line 28 to :
Field.InField = 1;

Now I have exactly the error you describe.

But with Field I have problems.
I have read the documentation and tutorial 10 (Gmsh 4.8.4) .
I have now

Field.F = Sprintf("Sqrt((Sqrt(x*x+z*z)-%g)^2+y*y)",R);


So I get no error message. Is this correct?

And I get the output

Dim = DoF: 60386
7.999999999999986
0.038275152052199134
7.961724847947793
23.999999999999904
0.7830624867869774
Solving linear variational problem.
Building point search tree to accelerate distance queries.
Computed bounding box tree with 3437 nodes for 1719 points.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
5.16191108202644
1.219863453306725e-05
1 5.1619232806609725
Solving linear variational problem.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
0.0
0.0
2 0.0
Process finished with exit code 0

But it seems that mft is still not correct and my problem with the visualization of the solution is the same as above.

But y is looking very good in the visualization Thanks and best regards,
noya

Hi @Noya,

A few things:

• What result do you expect for mft? When I visualize mft.pvd in Paraview, it appears as I would expect (equal to 22 inside the torus, and equal to 23 outside).
• What are you attempting to accomplish with the for loop? Since this does not appear to be a time-dependent problem, the loop should be unnecessary.
• The reason your error is 0 in the second step is because you have solved exactly the same problem in both the first and second steps. Check out the line solve(A == L, m, bc): all of the arguments to the solve function have no dependence on n. Thus the solution m is the same whether n==1 or n==2, and your error will be 0.
• If (as I suspect is the case) you are trying to monitor the convergence of your nonlinear problem, the for loop is unnecessary. Rather, you should replace solve(A == L, m, bc) with solve(F == 0, m, bc), as described in the Nonlinear Poisson tutorial. The syntax F == 0 signals to the solve method that you are solving a nonlinear problem (as opposed to A == L, which tells the solve method you have a linear problem).
• What precisely do you mean when you say “I can’t present p in the visualization of the solution with paraview”? When I open lsgtorus.xdmf in Paraview, I find two data arrays named yh and ph. If I use an “Extract Cells by Region” filter, I can see that ph is 0 outside the torus, but takes on nonzero values inside the torus.
1 Like

I can guess the torus in the cube, but it’s all one color. Is it possible to change this?

EDIT : rescale helps me here too Thank you

I solve an “own” algorithm instead of newton or similar.
(S1) Set n=1 and choose y_h^{(0)} \in V_h
(S2) solve the linear System for y_h^{(n)} \in V_h
\begin{equation} \int_{\Omega} \nu(x,\left|rot y^{n-1}_h\right|)rot y^{n}_h \cdot rotv_h dx + \int_{\Omega} v_h \cdot \nabla p^{n}_h dx = (f,v_h)_{L^2(\Omega)} \quad\forall v_h\in V_h \end{equation}
\begin{equation} \int_{\Omega} y^{n}_h\cdot \nabla q_hdx = 0 \quad\forall q \in S_h \end{equation}
(S3) If \parallel y_h^{n} -y_h^{n-1} \parallel_{H(rot)} +\parallel p_h^{n} -p_h^{n-1} \parallel_{H^1_0(\Omega)} < \epsilon STOP
(S4) Set n=n+1 and go to (S2)

Hm.I have solved the same problem with another geometry and there I started with y0,p0 and there the solution has been updated at each step and the total error at step 8 is less than tol. And actually nothing should change in the algorithm just because I change the geometry, or?

Yes, exactly, I actually have a nonlinear problem, but by using the initial value y0 in nu, my problem becomes a linear problem with the help of the algorithm.

This is exactly what I would expect, as well as for y.

but it seems that p is 0 everywhere

What do I wrong?
EDIT : I have rescaled ph and now I have

Best regards
Thank you
Noya

Ah, I understand now.

It seems that the only issue remaining is that problem converges in an unreasonably small number of iterations? I don’t see why anything should change in the algorithm. Unfortunately I can’t really think of any reason why this geometry seems to converge faster.

After thinking about this a bit more, have you checked the curl of y_h outside the torus (i.e. in \Omega - T)? On my machine, it appears that \textrm{rot}\,{y_h} is identically 0 outside the torus. Therefore, in your expression

\nu\left(x, |\textrm{rot}\,y_h|\right) = 1 - \frac{1}{2\left(|\textrm{rot}\,y_h|^2 + 1\right)}; \quad \forall x\in(\Omega-T)

the nonlinearity disappears and you have \nu = 1/2. This explains why your problem converges in one step.