OpenCascade geometry Periodic Boundary Conditions

Netgen
,
1

Hello,

I would like to impose periodic boundary conditions on the bottom and top of a space-time cylinder mesh.

To create a 2pi rotation in space-time of a subdomain in the domain, we use here the OpenCascade library of NGSolve, but it seems not to include he attribute PeriodicSurfaces() to impose the boundary condition:

[color=#000000] File “occ_periodic.py”, line 86, in [/color]
geoOCC.PeriodicSurfaces(bottom,top)
AttributeError: ‘netgen.libngpy._NgOCC.OCCGeometry’ object has no attribute ‘PeriodicSurfaces’

Is there a way to do this? We were looking at other solutions as well. For example the function GlobalInterfaceSpace(), but we are not sure of to use the mapping argument for our purposes (and if it is possible).

Thank you very much in advance for your help!

Alessio

2

Hi,
for opencascade the syntax is a little different, here a cube with 2 faces identified:


from netgen.occ import *
from ngsolve import *
from netgen.meshing import IdentificationType

a = Box((0,0,0),(1,1,1))

a.faces.Min(X).Identify(a.faces.Max(X),
                        "periodic",
                        IdentificationType.PERIODIC)

geo = OCCGeometry(a)
mesh = Mesh(geo.GenerateMesh(maxh=0.2))
Draw(mesh)

If your transformation between the faces is not a translation you have to create a netgen.occ.gp_Trsf and give the trafo as a fourth argument to identify (also if you want to identify multiple faces at once).

Look at the output if the identification was successful there should be something like

Map face 1 → 2

Also check identified points in the gui with view->Viewing Options->Mesh->Show identified Points

Best
Christopher

3

Hello Christopher,

thank you very much for your quick answer!

I have tried the command you suggested and it works in our case. We have a full 2*pi rotation and so a translation in Z-axis is fine.

Many thanks,
Alessio

4

Hi,

I am trying to impose periodic boundary conditions on the top and the bottom of a rotating object. Now, using a 360 degree rotation the periodic identification works pretty good. But when using only a 90 degree rotation, I have to create this trafo netgen.occ.gp_Trsf and use it as the fourth argument in the command as mentioned in the last post.

Here a short example, what I did so far:

from ngsolve import *
from netgen.occ import *
import math
import time
from ngsolve.internal import visoptions
from ngsolve.internal import viewoptions
from netgen.meshing import IdentificationType

origin = (0,0,0)
T = 1
desiredRad = 0.1
rotationAngle = pi/2

face1 = WorkPlane().Circle(0.25,0.25,desiredRad).Face()
spine = Segment(origin, (0,0,T))
cyl = Cylinder(origin, Z, r=0.5*desiredRad, h=T).faces[0]
heli = Edge(Segment((0,0), (rotationAngle, T)), cyl)
pipe1 = Pipe(spine, face1, auxspine=heli).mat(“pipe1”)
pipe1.faces.name=“pipe_outer”
pipe1.faces.Min(Z).name =“bottom”

rot = Rotation(Axis((0,0,0), Z), 180)
pipe1.faces.Min(Z).Identify(pipe1.faces.Max(Z), “periodic”, IdentificationType.PERIODIC, rot)

geoOCC = OCCGeometry(pipe1)
mesh = Mesh(geoOCC.GenerateMesh(maxh=0.1))
mesh.Curve(3)

Unfortunately, this does not work, after trying to adapt the rot object in several ways.

Thanks in advance for your help!

Best, Mario

5

Hi Mario,

you have to give the transformation from the first to the second face, which consists of rotation (by 90 degrees) and translation, then it works:

rot = Rotation(Axis((0,0,0), Z), 90)
trans = Translation(TZ)
pipe1.faces.Min(Z).Identify(pipe1.faces.Max(Z), “periodic”, IdentificationType.PERIODIC, rottrans)

best, Joachim

6

Dear Professor Schöberl,

thank you very much for taking the time to check this for us and find the solution, or better said our mistake.

Now I am happy, it doesn’t matter if the solution is or is not unique :wink:

All the best,
Alessio

7

Hi, I am trying to use this function with no success. here is my code so far. It builds a helical coil and bounding box and then I try to make the top and bottom faces of the box periodic. I would like both entire faces to be periodic. what I mean is that part of each has an interior face of the intersection of the coil with the box face.

I’ve determined, with the netgen webgui, which face numbers these are so 4 and 2 are bigger parts of the box faces I need to correspond. 7 and 8 are the respective smaller faces representing the ends of the coil.

from ngsolve import *
from netgen.occ import *
import math
from netgen.meshing import IdentificationType
from ngsolve.webgui import Draw
from netgen.webgui import Draw as DrawGeo

r = 0.01
h = 0.06
extra = h / 2

# cyl = (base pt, axis dir, radius, height)
cyl = Cylinder((0, 0, 0), Z, r, h).faces[0]
extracyl = Cylinder((0, 0, -extra), Z, r, h + 2 * extra).faces[0]
seg = Segment((0, 0), (2 * math.pi, h))
seg2 = Segment((0, 0), (4 * math.pi, 2 * h))
tseg = seg.Edge()
heli = Edge(seg, cyl)
heli2 = Edge(seg2, extracyl)
spiral = Wire([heli])
spiral2 = Wire([heli2])

piperadius = 0.003
circ = Face(Wire([Circle(heli.start, heli.start_tangent, piperadius)]))
coil = Pipe(spiral, circ)
circ2 = Face(Wire([Circle(heli2.start, heli2.start_tangent, piperadius)]))
coil2 = Pipe(spiral2, circ2)

box = Box((r * 3, r * 3, 0), (-r * 3, -r * 3, h))

coilfinal = coil2 - (coil2 - box)
air = box - coilfinal
geo = Glue([air, coilfinal])

# trans = Translation(h*Z)
geo.faces[8].Identify(geo.faces[7], "periodic", IdentificationType.PERIODIC)

geo.faces[4].Identify(geo.faces[2], "periodic", IdentificationType.PERIODIC)

# %%
# geo = OCCGeometry(boxandcoilmain)
geo = OCCGeometry(geo)
mesh = Mesh(geo.GenerateMesh(meshsize.coarse, maxh=0.3)).Curve(3)
Draw(mesh, clipping={"y": 1, "z": 0, "dist": 1.14})

If a python notebook version of this is preferred, I can post that too…

I have tried adding a transformation arg into the call. the one I use is seen commented out above the first Identify call. I don’t see any identified points when I turn that option on in the mesh settings.

also, i have shown that the identification/periodicity works if the pipe is straight but only on outer faces 2 and 4. the inner faces 7 and 8 do not work even in that base case even when those are the only two I try to match. here is the extra code for that which depends on the code above up to including the definition of coilfinal.

# the following is a way to confirm that the box cuts the coil at
# exactly one turn so that the end cap faces end up aligned. I extrude from
# one to the other linearly in z dir
p = coilfinal.faces[1].Extrude(h * Z)
pview = OCCGeometry(p)
straightcoil = p - (p - box)
air2 = box - straightcoil
straightexample = Glue([air2, straightcoil])
straightexample.faces[7].Identify(straightexample.faces[8], "periodic",
                                  IdentificationType.PERIODIC)
straightexample.faces[2].Identify(straightexample.faces[4], "periodic",
                                  IdentificationType.PERIODIC)

geoOCC = OCCGeometry(straightexample)
mesh = Mesh(geoOCC.GenerateMesh(meshsize.coarse, maxh=0.3)).Curve(3)
Draw(mesh, clipping={"y": 1, "z": 0, "dist": 1.14})

Thanks for your help.
Micron Jenkins

8

Hi Micron,

Regarding the internal surface, this should automatically be periodic, since Glue() creates a body that has two regions. Overlapping faces and edges merge into one, if you don’t want this you would have to mesh them separately and then merge them. Therefore, you will probably not find your identification points.

With respect to identification of faces etc.

To see if a face got periodic, simply set a maxh arg only for the master face and check after using .Identify() if the mesh is changing on both sides equally. If not, it did not work.

master_face.maxh=h

master_face.Identify(other=minor_surface)

You can also set a GridFunction() on a periodic FESpace to see if it worked.

In OCC the order of faces in the array body.faces can change when functions like fillet etc. are used. Therefore I would always use the name (or a positional argument like .Min(X)) of the face to determine the array index just before applying Identify() to avoid getting the wrong faces. Furthermore, the periodic property of a face can be lost again if shaping operations as mentioned above are used afterwards. Therefore, always use it only at the end.

All the best,

Leonhard