This page was generated from unit-2.12-periodicity/periodicity.ipynb.

2.12 Periodic Spaces

To define spaces with periodic constraints, we have to create meshes where the nodes on one side are identified with nodes on the opposite side.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *

shape = Rectangle(1,1).Face()

shape.edges.Max(X).name = "right"
shape.edges.Min(X).name = "left"
shape.edges.Max(Y).name = "top"
shape.edges.Min(Y).name = "bot"

shape.edges.Max(Y).Identify(shape.edges.Min(Y), "bt")
shape.edges.Max(X).Identify(shape.edges.Min(X), "lr")

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.1))

plist = []
for pair in mesh.ngmesh.GetIdentifications():
    plist += list(mesh.vertices[pair[0]-1].point) + [0]
    plist += list(mesh.vertices[pair[1]-1].point) + [0]
Draw(mesh, objects=[{"type" : "lines", "position" : plist, "name": "identification", "color": "purple"}]);

fes = Periodic(H1(mesh,order=3))
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx+u*v*dx).Assemble()
f = LinearForm(exp(-100*( (x-0.8)**2+(y-0.8)**2))*v*dx).Assemble()

gfu = GridFunction(fes,"u")
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

Draw (gfu);

A piece of cake

f = WorkPlane(Axes((0,0,0), Y,X)).MoveTo(0.3,0).Rectangle(3,1).Face()
ax = Axis ((0,0,0), Z)
cake = f.Revolve(ax, 30)
cake.faces.Min(Y).name="f1"
cake.faces.Max(Y-0.5*X).name="f2"
cake.faces.Min(Z).name="bot"

cake.faces["f1"][0].Identify(cake.faces["f2"][0], "id",
                            trafo=Rotation(ax, 30))
Draw (cake);

NGSolve does not support elements having dofs on the primary and secondary side. To avoid them we refine the mesh once:

mesh = Mesh(OCCGeometry(cake).GenerateMesh(maxh=0.5)) # .Curve(3)
mesh.ngmesh.Refine()

plist = []
for pair in mesh.ngmesh.GetIdentifications():
    plist += list(mesh.vertices[pair[0]-1].point)
    plist += list(mesh.vertices[pair[1]-1].point)
Draw(mesh, objects=[{"type" : "lines", "position" : plist, "name": "identification", "color": "purple"}]);

fes = Periodic(H1(mesh, order=2, dirichlet="bot"))

u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
f = LinearForm(100*exp(-9*( (x-2.5)**2+y**2+(z-0.5)**2))*v*dx).Assemble()

gfu = GridFunction(fes,"u")
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

Draw (gfu);

Application: computing band diagrams

Computing band diagrams