# Define a Torus in the UnitCubemesh

Hey,
I have a domain \Omega in \mathbb{R}^3, which is the unit cube. And in this domain there should be a torus T, so that my right side L is only on the torus (L=0 outside the torus). Is there any documentation where I can read something about this? Or does anyone else have ideas and could give me some help?
I have never worked with something like this before and so i am a bit confused.

best regards,
noya

You could probably generate the mesh with gmsh with a torus inclusion using OpenCASCADE.

Thank you.

I think I have created the right domain.

Can you please take a look and check it? I am not so sure about the boundaries.

// Gmsh project created on Wed Aug  4 22:50:25 2021
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {1, 1, 0, 1.0};
//+
Point(4) = {0, 1, 0, 1.0};
//+
Point(5) = {1, 0, 1, 1.0};
//+
Point(6) = {0, 0, 1, 1.0};
//+
Point(7) = {0, 1, 1, 1.0};
//+
Point(8) = {1, 1, 1, 1.0};
//+
Line(1) = {6, 5};
//+
Line(2) = {2, 3};
//+
Line(3) = {5, 2};
//+
Line(4) = {2, 1};
//+
Line(5) = {1, 4};
//+
Line(6) = {4, 3};
//+
Line(7) = {3, 8};
//+
Line(8) = {8, 7};
//+
Line(9) = {7, 4};
//+
Line(10) = {1, 6};
//+
Line(11) = {6, 7};
//+
Line(12) = {8, 5};
//+
Curve Loop(1) = {10, 11, 9, -5};
//+
Physical Curve(1) = {9, 5, 10, 3, 2, 7, 8, 12, 11, 1, 6, 4};
//+
Curve Loop(2) = {9, -5, 10, 11};
//+
Plane Surface(1) = {2};
//+
Curve Loop(3) = {8, -11, 1, -12};
//+
Plane Surface(2) = {3};
//+
Curve Loop(4) = {7, 8, 9, 6};
//+
Plane Surface(3) = {4};
//+
Curve Loop(5) = {5, 6, -2, 4};
//+
Plane Surface(4) = {5};
//+
Curve Loop(6) = {4, 10, 1, 3};
//+
Plane Surface(5) = {6};
//+
Curve Loop(7) = {12, 3, 2, 7};
//+
Plane Surface(6) = {7};
//+
Physical Surface(2) = {2, 3, 4, 5, 1, 6};
//+
Torus(1) = {0.5, 0.5, 0.5, 0.2, 0.1, 2*Pi};
//+
Rotate {{1, 0, 0}, {0, 0, 0}, -Pi/2} {
Volume{1};
}
//+
Translate {0, 0, 1} {
Volume{1};
}
//+
Physical Curve(3) = {13, 14};
//+
Physical Surface(4) = {7};


And if it’s ok I have to convert the mesh so I can use it in my code.

But how can I say : “the right side L is in the torus and the rest lives in the cube”?

Hey,
I tried to read a 3D mesh from gmsh with

from fenics import*
import meshio

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


but I get the error

usr/bin/python3.8 /home/noya/donut.py
Traceback (most recent call last):
File “/home/noya/donut.py”, line 5, in
meshio.write(“mesh.xdmf”, meshio.Mesh(points=msh.points, cells={“tetra”: msh.cells[“tetra”]}))
TypeError: list indices must be integers or slices, not str

What do I wrong? Of sure it can be that my mesh is wrong.

See for instance: Accessing and marking imported boundaries - #8 by dokken
On how to use newer versions of meshio

1 Like

Hello again,

I have now been advised by my supervisor to solve the problem via “submesh”.
But for “submesh” I need a “subdomain”.
This “subdomain” should be my torus.
The representation of a torus is
\left(\begin{array}{c} x \\ y \\ z \end{array}\right)= \left(\begin{array}{c} (R+r\cdot cos(\phi))sin(\sigma) \\ (R+r\cdot cos(\phi))sin(\sigma) \\ r\cdot sin(\phi) \end{array}\right)
or in Cartesian coordinates
(x^2+y^2+z^2+R^2-r^2)^2=4R^2(x^2+y^2)

The torus should lie in the center with for example R=0,2 and r =0,01.

How can I define it?
I tried

from fenics import*

mesh = UnitcubeMesh(5, 5, 5)

r = Constant(0.01)
R = Constant(0.2)

class torus(SubDomain):
def inside(self, x, on_boundary):
return ((R+r*cos(alpha))*cos(beta)<=x[0] <= (R+r*cos(alpha))*cos(beta)and (R+r*cos(alpha))*cos(beta)<=x[1] <= (R+r*cos(alpha))*cos(beta) and r*sin(alpha)<=x[2]<= r*sin(alpha))


but I don’t know how to define the return argument so that I get the torus.
How can i define that alpha, beta \in [0,2\pi]
Can someone help me and explain this?

Thanks,

Noya

Hi @Noya,

Note that per Wikipedia, the interior of the torus (including the boundary) satisfies

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

Therefore, you should be able to define your SubDomain using

class torus(SubDomain):
def inside(self, x, on_boundary):
return (sqrt(x[0]*x[0] + x[1]*x[1]) - R)**2 + x[2]*x[2] <= r**2

1 Like

Thank you.
But I have know the following problem

from fenics import*

mesh = UnitCubeMesh(10, 10, 10)
r = Constant(0.01)
R = Constant(0.2)
class Torus(SubDomain):
def inside(self, x, on_boundary):
return (sqrt(x[0] * x[0] + x[1] * x[1]) - R) ** 2 + x[2] * x[2] <= r ** 2

torus = Torus()

# Initialize mesh function
mf = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
mf.set_all(0)

dx1 = Measure('dx', domain=mesh, subdomain_data=mf)

submesh_torus = SubMesh(mesh, mf, 1)
# -------------------------------------------------
torus_file = File('torus.pvd')
torus_file << submesh_torus


/usr/bin/python3.8 /home/donut.py
*** Warning: Mesh is empty, unable to create entities of dimension 0.
Process finished with exit code 0

Why is the mesh empty? Can somebody explain it to me?

Thanks,
Noya

There are quite a few issues with your code:

1. You do not actually use the torus variable to mark the MeshFunction.

After correcting for this, there are quite a few errors due to the usage of dolfin variables, such as Constant and sqrt. These can be corrected as in the code below.

1. Most crucially, the torus you are creating is way to narrow, i.e. r=0.01 is smaller than your mesh size of 0.1. Also, the UnitCube does not encapsulate the whole torus, and thus you need to use another mesh.

I have addressed all these issues in the code below:

from fenics import*
import numpy as np
N = 50
mesh = BoxMesh(Point(-1, -1, -1), Point(1, 1, 1), N, N, N)
r = Constant(0.1)
R = Constant(0.2)

class Torus(SubDomain):
def inside(self, x, on_boundary):
return (np.sqrt(x[0] * x[0] + x[1] * x[1]) - float(R)) ** 2 + x[2] * x[2] <= float(r) ** 2

torus = Torus()

# Initialize mesh function
mf = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
mf.set_all(0)
torus.mark(mf, 1)
File("mf.pvd") << mf
dx1 = Measure('dx', domain=mesh, subdomain_data=mf)

submesh_torus = SubMesh(mesh, mf, 1)
# -------------------------------------------------
torus_file = File('torus.pvd')
torus_file << submesh_torus

1 Like

Thank you.
The

torus.mark(mf, 1)


always caused me errors before and that’s why I commented it out…

N = 50
mesh = BoxMesh(Point(-1, -1, -1), Point(1, 1, 1), N, N, N)


This creates a cube of the type [-1,1]^3 with mesh with the refinement 50 in each direction??

Now there is createt an \Omega=[-1,1]^3 and the torus T with T \in \Omega. Can I create the set \Omega/T ? Respectively I want to create a function which is constant in T and which is a nonlinear function on \Omega/T.

Can I do it like in the tutorial Subdomains and boundary conditions?

z0 = Function(W)
y0, p0 = split(z0)
tol = 1E-14
s = inner(curl(y0), curl(y0))
a_c = 1 - (1 / (2 * (s + 1)))
a_t = 1
a = Expression('np.sqrt(x[0] * x[0] + x[1] * x[1]) - float(R)) ** 2 + x[2] * x[2] <= float(r) ** 2 ?a_t: a_c ', degree=0,
tol=tol, a_t =a_t, a_c=a_c)


And I have to define the boundary conditions for my torus and \Omega, but I’m not sure how to do it.

class TBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and **???**
t_boundary = TBoundary()


I have set there questionsmarks cause I think there is something missing.
and for \Omega i need Dirichletboundary conditions. Normally I have done it with

def boundary(x, on_boundary):
return on_boundary
w = Function(W)
bc = DirichletBC(W, w, boundary)


with the Functionspace W is defined as

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


It would be nice if someone can answer on my questions and give me the link to documentation or something else.

Best regards,
Noya

This is correct.

You can for instance use MeshView, as described in detail in [1911.01166] Abstractions and automated algorithms for mixed domain finite element methods
with the code available at: GitHub - cdaversin/mixed-dimensional-examples: Code to reproduce numerical examples presented in "Abstractions and automated algorithms for mixed-dimensional finite element methods" (2019)

Thank you.
I am still having trouble understanding how MeshView works.
The following code is from mixed-dimensional-examples/poisson_lm_3D.py at master · cdaversin/mixed-dimensional-examples · GitHub

There a unitcubemesh is created and then for “marker” a meshfunction with the mark facets. And now, as long as f is in facets(mesh), marker[f] is set equal to 0.5-DOLFIN_EPS smaller than the center of f for all x smaller than 0.5+DOLFIN_EPS.


from dolfin import *

n_values = [2, 4, 8, 16, 32]
exact = Expression("x[0]*(1-x[0])", degree=2)
order= 1

def poisson(n):
mesh = UnitCubeMesh(n, n, n)

marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
marker[f] = 0.5 - DOLFIN_EPS < f.midpoint().x() < 0.5 + DOLFIN_EPS


Now the submesh is generated here, that I have a slice? Or which one is chosen?

    ## Create submesh ##
submesh = MeshView.create(marker, 1)


## Setup formulation ##
V = FunctionSpace(mesh, "CG", order)
LM = FunctionSpace(submesh, "CG", order)
W = MixedFunctionSpace(V,LM)

def boundary(x):
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
bc = DirichletBC(V, Constant(0.0), boundary)

(u,l) = TrialFunctions(W)
(v,e) = TestFunctions(W)

dV = Measure("dx", domain=W.sub_space(0).mesh())
dL = Measure("dx", domain=W.sub_space(1).mesh())


And how can I create with this my \Omega/T?