Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Periodicity for a parallelepiped

More
2 years 4 months ago - 2 years 3 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 3 months ago by s.yang.
More
2 years 3 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 3 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.113 seconds