Periodicity for a parallelepiped

5 months 3 days ago - 4 months 3 weeks 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 . 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


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():
f.Assemble() = a.mat.Inverse(fes.FreeDofs()) * f.vec

Last edit: 4 months 3 weeks ago by s.yang.

Please Log in or Create an account to join the conversation.

4 months 1 week 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.
geo.PeriodicSurfaces(left,right, Trafo((0.5*sqrt(3),-0.5,0))) geo.PeriodicSurfaces(front,back, Trafo((-0.5*sqrt(3),-0.5,0)))

Please Log in or Create an account to join the conversation.

4 months 1 week 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 . It is still great to know this easy solution.

Many thanks again.

Best wishes

Please Log in or Create an account to join the conversation.

Time to create page: 0.128 seconds