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 Periodicity — NGS-Py 6.2.1907 documentation. 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.5sqrt(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

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)