Periodicity for a parallelepiped

More
2 years 10 months ago - 2 years 10 months ago #4124 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)
Last edit: 2 years 10 months ago by s.yang.
More
2 years 10 months ago #4161 by christopher
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)))
More
2 years 10 months ago #4163 by s.yang
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
Time to create page: 0.100 seconds