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
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”?
Hey,
I tried to read a 3D mesh from gmsh with
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.
See for instance: Accessing and marking imported boundaries - #8 by dokken
On how to use newer versions of meshio
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
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
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:
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.
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
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?
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
Thank you. I have create a new post with all informations. Mesh, Submesh , Torus, Boundary