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[0] * x[0] + x[1] * x[1]) - float(R)) ** 2 + x[2] * x[2] <= 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[0] * x[0] + x[1] * x[1]) - float(R)) ** 2 + x[2] * x[2] >= float(r) ** 2 and -1 <= x[0] <= 1 and -1 <= x[1] <= 1 and -1 <= x[0] <= 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
Thank you
