- Thank you received: 0
Periodicity for a parallelepiped
2 years 10 months ago - 2 years 10 months ago #4124
by s.yang
Periodicity for a parallelepiped was created by s.yang
Hi everyone,
I am trying to solve a problem in a 3D parallelepiped domain (in particular a right parallelogrammic prism or a right rhombic prism) with periodic boundary conditions on two pairs of faces (left, right and front, back) and Dirichlet boundary conditions on the other two faces (top, bottom).
I know we can build periodic meshes for a cube as in ngsolve.org/docu/v6.2.1907/how_to/periodic.html . I tried to use this in a parallelepiped domain, but the solution didn't have the desired periodicity. I checked that everything is fine for a 2D analogue (parallelogram domain) or a cube domain in 3D situation.
I solve a Laplace equation with f=x^2 on the domain of a parallelepiped generated by vertices (0,0,0), (sqrt(3)/2,1/2,0), (sqrt(3)/2,-1/2,0) and (0,0,5). What I did is copied here in the end of message.
Could this be a bug related to the CSGometry()? Or did I make a mistake? I appreciate any suggestion on this. Many thanks and happy holidays!
from netgen.csg import *
geo = CSGeometry()
left = Plane(Pnt(0,0,0),Vec(-0.5,0.5*sqrt(3),0))
right = Plane(Pnt(0.5*sqrt(3),-0.5,0),Vec(0.5,-0.5*sqrt(3),0))
bot = Plane(Pnt(0,0,0),Vec(0,0,-1)).bc('bc_dirichlet')
top = Plane(Pnt(0,0,5),Vec(0,0,1)).bc('bc_dirichlet')
back = Plane(Pnt(0,0,0),Vec(-0.5,-0.5*sqrt(3),0))
front = Plane(Pnt(0.5*sqrt(3),0.5,0),Vec(0.5,0.5*sqrt(3),0))
parallelepiped = left * right * front * back * bot * top
geo.Add(parallelepiped)
geo.PeriodicSurfaces(left,right)
geo.PeriodicSurfaces(front,back)
ngmesh1 = geo.GenerateMesh(maxh=0.1)
from ngsolve import *
mesh = Mesh(ngmesh1)
fes = Periodic(H1(mesh,order=3,dirichlet="bc_dirichlet"))
u,v = fes.TrialFunction(), fes.TestFunction()
a = BilinearForm(fes,symmetric=True)
a += SymbolicBFI(grad(u) * grad(v))
f = LinearForm(fes)
f += SymbolicLFI(x*x*v)
u = GridFunction(fes,"u")
with TaskManager():
a.Assemble()
f.Assemble()
u.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw(u)
I am trying to solve a problem in a 3D parallelepiped domain (in particular a right parallelogrammic prism or a right rhombic prism) with periodic boundary conditions on two pairs of faces (left, right and front, back) and Dirichlet boundary conditions on the other two faces (top, bottom).
I know we can build periodic meshes for a cube as in ngsolve.org/docu/v6.2.1907/how_to/periodic.html . I tried to use this in a parallelepiped domain, but the solution didn't have the desired periodicity. I checked that everything is fine for a 2D analogue (parallelogram domain) or a cube domain in 3D situation.
I solve a Laplace equation with f=x^2 on the domain of a parallelepiped generated by vertices (0,0,0), (sqrt(3)/2,1/2,0), (sqrt(3)/2,-1/2,0) and (0,0,5). What I did is copied here in the end of message.
Could this be a bug related to the CSGometry()? Or did I make a mistake? I appreciate any suggestion on this. Many thanks and happy holidays!
from netgen.csg import *
geo = CSGeometry()
left = Plane(Pnt(0,0,0),Vec(-0.5,0.5*sqrt(3),0))
right = Plane(Pnt(0.5*sqrt(3),-0.5,0),Vec(0.5,-0.5*sqrt(3),0))
bot = Plane(Pnt(0,0,0),Vec(0,0,-1)).bc('bc_dirichlet')
top = Plane(Pnt(0,0,5),Vec(0,0,1)).bc('bc_dirichlet')
back = Plane(Pnt(0,0,0),Vec(-0.5,-0.5*sqrt(3),0))
front = Plane(Pnt(0.5*sqrt(3),0.5,0),Vec(0.5,0.5*sqrt(3),0))
parallelepiped = left * right * front * back * bot * top
geo.Add(parallelepiped)
geo.PeriodicSurfaces(left,right)
geo.PeriodicSurfaces(front,back)
ngmesh1 = geo.GenerateMesh(maxh=0.1)
from ngsolve import *
mesh = Mesh(ngmesh1)
fes = Periodic(H1(mesh,order=3,dirichlet="bc_dirichlet"))
u,v = fes.TrialFunction(), fes.TestFunction()
a = BilinearForm(fes,symmetric=True)
a += SymbolicBFI(grad(u) * grad(v))
f = LinearForm(fes)
f += SymbolicLFI(x*x*v)
u = GridFunction(fes,"u")
with TaskManager():
a.Assemble()
f.Assemble()
u.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw(u)
Last edit: 2 years 10 months ago by s.yang.
- christopher
- Offline
- Administrator
Less
More
- Thank you received: 101
2 years 10 months ago #4161
by christopher
Replied by christopher on topic Periodicity for a parallelepiped
There is a argument to the PeroidicSurfaces function that lets you specify the mapping. As default it assumes it to be normal to the planes.
Code:
geo.PeriodicSurfaces(left,right, Trafo((0.5*sqrt(3),-0.5,0)))
geo.PeriodicSurfaces(front,back, Trafo((-0.5*sqrt(3),-0.5,0)))
2 years 10 months ago #4163
by s.yang
Replied by s.yang on topic Periodicity for a parallelepiped
Hi Christopher,
Thanks a lot for pointing this out! I have already solved the problem by setting up the meshes manually as in ngsolve.org/forum/ngspy-forum/106-questi...esh-generation-in-3d . It is still great to know this easy solution.
Many thanks again.
Best wishes
Thanks a lot for pointing this out! I have already solved the problem by setting up the meshes manually as in ngsolve.org/forum/ngspy-forum/106-questi...esh-generation-in-3d . It is still great to know this easy solution.
Many thanks again.
Best wishes
Time to create page: 0.100 seconds