2. Meshing from Neper geometries

Neper is a free / open source software package for polycrystal generation and meshing, see https://neper.info

To run this notebook you first have to install neper.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import subprocess, random

n = 500
subprocess.run(["neper", "-T", "-n", str(n), "-format", "tess"])

========================    N   e   p   e   r    =======================
Info   : A software package for polycrystal generation and meshing.
Info   : Version 4.9.1-6
Info   : Built with: gsl|muparser|opengjk|openmp|nlopt|libscotch (full)
Info   : Running on 8 threads.
Info   : <https://neper.info>
Info   : Copyright (C) 2003-2022, and GNU GPL'd, by Romain Quey.
Info   : No initialization file found (`/Users/joachim/.neperrc').
Info   : ---------------------------------------------------------------
Info   : MODULE  -T loaded with arguments:
Info   : [ini file] (none)
Info   : [com line] -n 500 -format tess
Info   : ---------------------------------------------------------------
Info   : Reading input data...
Info   : Creating domain...
Info   : Creating tessellation...
Info   :   - Setting seeds... 
Info   :   - Running tessellation...
Info   : Generating crystal orientations...
Info   : Writing results...
Info   :     [o] Writing file `n500-id1.tess'...
Info   :     [o] Wrote file `n500-id1.tess'.
Info   : Elapsed time: 0.058 secs.
========================================================================

CompletedProcess(args=['neper', '-T', '-n', '500', '-format', 'tess'], returncode=0)

verts = {}
edges = {}
faces = {} 
solids = []

f = open("n" + str(n) + "-id1.tess")

def randomcol():
    r = random.uniform(0, 1)
    g = random.uniform(0, 1)
    b = random.uniform(0, 1)
    return (r,g,b)


while True:
    line = f.readline()
    if not line: break

    if line.split()[0] == "**vertex":
        num = int(f.readline())
        print ("found", num, "vertices")

        for i in range(num):
            line = f.readline()            
            nr,x,y,z,hhh = line.split()
            verts[int(nr)] = Vertex(Pnt(float(x), float(y), float(z)))

    if line.split()[0] == "**edge":
        num = int(f.readline())
        print ("found", num, "edges")

        for i in range(num):
            line = f.readline()            
            nr,i1,i2,hhh = line.split()
            edge = Edge(verts[int(i1)], verts[int(i2)])
            edges[int(nr)] = edge
            edges[-int(nr)] = edge.Reversed()
            
    if line.split()[0] == "**face":
        num = int(f.readline())
        print ("found", num, "faces")

        for i in range(num):
            l1 = f.readline()            
            l2 = f.readline()            
            l3 = f.readline()            
            l4 = f.readline()
            nr,*_ = l1.split()
            nume,*fedges = l2.split()
            
            face = Face(Wire( [edges[int(enr)]  for enr in fedges] ))
            faces[int(nr)] = face
            faces[-int(nr)] = face.Reversed()

            
    if line.split()[0] == "**polyhedron":
        num = int(f.readline())
        print ("found", num, "polyhedra")

        for i in range(num):
            nr,nrf,*polyfaces = f.readline().split()
            solids.append (Solid(Glue( [faces[int(fnr)] for fnr in polyfaces ] )))
            solids[-1].faces.col = randomcol()                        
            
shape = Glue(solids)
geo = OCCGeometry(shape)

Draw (shape);

found 3000 vertices
found 5996 edges
found 3497 faces
found 500 polyhedra

mp = netgen.meshing.MeshingParameters(segmentsperedge=0.01)  # don't restrict mesh-size for short edges
with TaskManager(): # pajetrace=10**9):
    mesh = geo.GenerateMesh(maxh=0.05, mp=mp)
print ("ne =", mesh.ne)

ne = 61315

Draw (mesh); 

# ngsglobals.msg_level=3
for t in Timers():
    if t["time"] > 0.5 or "MeshSmoothing" in t["name"] or "CreateSorted" in t["name"]:
        print (t)

{'name': 'CreateSortedTable', 'time': 0.012723039726541333, 'counts': 9, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'CreateSortedTable', 'time': 0.02327775230313894, 'counts': 18, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'OptimizeVolume', 'time': 0.83072191421068, 'counts': 1, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'CreateSortedTable', 'time': 0.004418909341025568, 'counts': 130, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'CreateSortedTable', 'time': 0.018465976499903704, 'counts': 130, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'CreateSortedTable', 'time': 0.06590248851542903, 'counts': 1306, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'Meshing3::Delaunay', 'time': 0.8495747045590126, 'counts': 130, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'CreateSortedTable', 'time': 0.04904309777539196, 'counts': 786, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'MeshVolume', 'time': 0.5174805119012881, 'counts': 1, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'Combineimprove 2D start', 'time': 0.5200971695360285, 'counts': 31473, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'Combineimprove 2D', 'time': 0.6557157864020968, 'counts': 31473, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'MeshSmoothing 2D', 'time': 0.6876867731043612, 'counts': 62946, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'Optimization 2D', 'time': 1.6466368335487098, 'counts': 1, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'rest of loop', 'time': 1.6711248072033023, 'counts': 3497, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'Surface Meshing', 'time': 4.542678435949811, 'counts': 1, 'flops': 0.0, 'Gflop/s': 0.0}