Periodic cube mesh

Good morning,

I have a problem with periodic meshes :wink:

Attached, there is a script which defines a periodic cube of length L and a structured cube mesh (Joachim posted this somewhere in the forum).
However, I cannot transfer the periodicity to the cube mesh.

This problem is already in 2D. I generally do not know how to tell a manually generated mesh that it should have periodic facets…

I would be grateful for any help.

Best, Philipp

Attachment: PeriodicCubeMesh.py

Hi Philipp,

The only thing you do have to do is to set the vertex identification for all your periodic points. This is done by the function mesh.AddPointIdentification. It takes the point identifications of the two points, the identification number (1 based, which periodic boundary it’s on if you have multiple ones) and the identification type (which is 2 for periodic). If you change your first loop where you add the points to

pids = [] for i in range(N+1): for j in range(N+1): for k in range(N+1): pids.append (mesh.Add(MeshPoint(Pnt(L*i/N, L*j/N, L*k/N)))) if k is 0: slave = pids[-1] if k is N: mesh.AddPointIdentification(pids[-1],slave,identnr=1,type=2)

it works for me. I’ve attached the file with a simple h1 test.
Another tip: If you select in the gui the mesh in the visualization and then go to View-> Viewing options->Mesh and disable “show filled triangles” and enable “show identified points” you can see if your point identification is correct. This helps a lot when debugging :wink:

Best
Christopher

Attachment: PeriodicCubeMesh.py

Thank you very much, Christopher!

The test also works on my computer and the “show identified points” is also very neat :wink:

If I see this correctly, the mesh is now periodic only in one direction.
Somehow, I cannot identify the points in a second direction.
I always end up with identifying one corner point with all points on the opposite face…

So, could you please tell me how to extend the periodicity in one direction to periodicity in, let’s say, all directions? :slight_smile:

Best, Philipp

Hi Philipp,

if you identify the other directions, make sure you increase the “identnr” when you use

mesh.AddPointIdentification(...,identnr=2,type=2)

This function creates an array with identifications for the vertices of two opposite faces.
If you now use the same “identnr” for the identification of two other faces, you overwrite some data.

Best,
Christoph

Hi Christoph,

Thank you for your help! I finally managed to get a periodic box manually now :slight_smile:

Best, Philipp

Hey guys,

I actually have sort of a follow-up question on manually generated meshes.

In the script in this discussion, each face gets a name via e.g. FaceDescriptor(…, bc=1), right?
So this name in form of [int list] is used to specify Dirichlet boundaries, the Set() operator, or billinear form integrator which act only on part of the boundary.

If I use the integer list, I get the following warning:

 warning: Boundaries( [int list] ) is deprecated, pls generate Region 

So, how can I generate regions (I guess named in form of strings?) for the boundary faces in this example?

Thank you very much for your help!

Best, Philipp

Attachment: PeriodicCubeMesh_2.py

The preferred way is to use names instead of numbers.
If you don’t have names, you can create a region from a BitArray, like

Region(mesh,BND, BitArray([True,False,False,True]))

This region is defined on the first and last boundary.

Thank you, Joachim!

How can I use this region in, for example, “definedon=mesh.Boundaries(…)”?
Is it possible to assign names to regions?

Best, Philipp

You can set boundary names on a Netgen mesh with
mesh.SetBCName(index, name)
BC indices are 1 based.

Thank you, Christopher!

In the attached script, I tried to do it this way.
However, I always get the following error:

illegal bc-number

Can you tell me what I am doing wrong?

Thanks for your help.

Best, Philipp

Attachment: PeriodicCubeMesh_3.py

This is a bit messy… A quickfix is to set names for all boundary conditions. I will discuss a nice solution with Joachim when he’s back in Vienna.
Best
Christopher