# Define a Torus in the UnitCubemesh

**URL:** <https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350>\
**Category:** mesh\
**Created:** [August 4, 2021, 2:00pm UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350 "2021-08-04T14:00:29Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 4, 2021, 2:00pm UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/1 "2021-08-04T14:00:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [August 4, 2021, 2:05pm UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/2 "2021-08-04T14:05:13Z")

</div>

You could probably generate the mesh with [gmsh](https://gmsh.info/doc/texinfo/gmsh.html) with a [torus inclusion](https://gmsh.info/doc/texinfo/gmsh.html#index-gmsh_002fmodel_002focc_002faddTorus) using [OpenCASCADE](https://gmsh.info/doc/texinfo/gmsh.html#t16).

---

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 5, 2021, 11:11am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/3 "2021-08-05T11:11:33Z")

</div>

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.

```auto
// Gmsh project created on Wed Aug 4 22:50:25 2021
SetFactory("OpenCASCADE");
//+
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”?

---

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 6, 2021, 7:49am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/4 "2021-08-06T07:49:42Z")

</div>

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

```auto
from fenics import*
import meshio

msh = meshio.read("mesh.msh")
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.

---

<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:** [August 6, 2021, 8:13am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/5 "2021-08-06T08:13:10Z")

</div>

See for instance: [Accessing and marking imported boundaries - #8 by dokken](https://fenicsproject.discourse.group/t/accessing-and-marking-imported-boundaries/5753/8)  
On how to use newer versions of meshio

---

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 6, 2021, 10:42am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/6 "2021-08-06T10:42:19Z")

</div>

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

```auto
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

---

<div class="post-metadata">

**Author:** ![conpierce8](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/conpierce8/32/2464_2.png) [@conpierce8](https://fenicsproject.discourse.group/u/conpierce8)\
**Post date:** [August 6, 2021, 1:58pm UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/7 "2021-08-06T13:58:58Z")

</div>

Hi @Noya,

Note that per [Wikipedia](https://en.wikipedia.org/wiki/Torus), 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

```python
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

```

---

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 7, 2021, 3:50am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/8 "2021-08-07T03:50:39Z")

</div>

Thank you.  
But I have know the following problem

```auto
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

---

<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:** [August 7, 2021, 7:19am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/9 "2021-08-07T07:19:35Z")

</div>

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:

```python
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

```

---

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 7, 2021, 4:31pm UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/10 "2021-08-07T16:31:25Z")

</div>

Thank you.  
The

```auto
torus.mark(mf, 1)

```

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

```auto
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](https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html)?

```auto
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.

```auto
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

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

```

with the Functionspace W is defined as

```auto
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

---

<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:** [August 8, 2021, 11:19am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/11 "2021-08-08T11:19:27Z")

</div>

> [@Noya](#):
>
> This creates a cube of the type [-1,1]^3 [−1,1]3[-1,1]^3 with mesh with the refinement 50 in each direction??

This is correct.

> [@Noya](#):
>
> Now there is createt an \Omega=[-1,1]^3 Ω=[−1,1]3\Omega=[-1,1]^3 and the torus T TT with T \in \Omega T∈ΩT \in \Omega . Can I create the set \Omega/T Ω/T\Omega/T ? Respectively I want to create a function which is constant in T and which is a nonlinear function on \Omega/T Ω/T\Omega/T .

You can for instance use `MeshView`, as described in detail in [[1911.01166] Abstractions and automated algorithms for mixed domain finite element methods](https://arxiv.org/abs/1911.01166)  
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)](https://github.com/cdaversin/mixed-dimensional-examples)

---

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 9, 2021, 8:31am UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/12 "2021-08-09T08:31:14Z")

</div>

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](https://github.com/cdaversin/mixed-dimensional-examples/blob/master/Poisson/poisson_lm_3D.py)

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.

```auto

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?

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

```

```auto

    ## 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?

---

<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:** [August 9, 2021, 4:55pm UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/13 "2021-08-09T16:55:13Z")

</div>

Please make a minimal example.  
As you have strayed quite far from the initial topic of the post, I would suggest one of the following

- Create a new post, where you formulate the mathematical problem you are trying to solve, and a minimal code illustrating how far you have gotten in your implementaton.
- This minimal code should be possible for anyone to execute, and should only contain code that is relevant to the example you are trying to solve for (no post processing for instance).
- For a list of good minimal examples that you can try to replicate, see: [Read before posting: How do I get my question answered? - #4 by nate](https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21/4)

---

<div class="post-metadata">

**Author:** ![Noya](https://avatars.discourse-cdn.com/v4/letter/n/ccd318/32.png) [@Noya](https://fenicsproject.discourse.group/u/Noya)\
**Post date:** [August 10, 2021, 2:09pm UTC](https://fenicsproject.discourse.group/t/define-a-torus-in-the-unitcubemesh/6350/14 "2021-08-10T14:09:17Z")

</div>

Thank you. I have create a new post with all informations. [Mesh, Submesh , Torus, Boundary](https://fenicsproject.discourse.group/t/mesh-submesh-torus-boundary/6378)