netgen/python/geom2d.py

from netgen.libngpy._geom2d import *
from netgen.libngpy._meshing import *

tmp_generate_mesh = SplineGeometry.GenerateMesh

def geom2d_meshing_func (geom, **args):
    if "mp" in args:
        return tmp_generate_mesh (geom, args["mp"])
    else:
        return tmp_generate_mesh (geom, MeshingParameters (**args))


SplineGeometry.GenerateMesh = geom2d_meshing_func


unit_square = SplineGeometry()
pnts = [ (0,0), (1,0), (1,1), (0,1) ]
lines = [ (0,1,1,"bottom"), (1,2,2,"right"), (2,3,3,"top"), (3,0,4,"left") ]
pnums = [unit_square.AppendPoint(*p) for p in pnts]
for l1,l2,bc,bcname in lines:
    unit_square.Append( ["line", pnums[l1], pnums[l2]], bc=bcname)


def MakeRectangle (geo, p1, p2, bc=None, bcs=None, **args):
    p1x, p1y = p1
    p2x, p2y = p2
    p1x,p2x = min(p1x,p2x), max(p1x, p2x)
    p1y,p2y = min(p1y,p2y), max(p1y, p2y)
    if not bcs: bcs=4*[bc]
    pts = [geo.AppendPoint(*p) for p in [(p1x,p1y), (p2x, p1y), (p2x, p2y), (p1x, p2y)]]
    for p1,p2,bc in [(0,1,bcs[0]), (1, 2, bcs[1]), (2, 3, bcs[2]), (3, 0, bcs[3])]:
        geo.Append( ["line", pts[p1], pts[p2]], bc=bc, **args)

def MakeCircle (geo, c, r, **args):
    cx,cy = c
    pts = [geo.AppendPoint(*p) for p in [(cx,cy-r), (cx+r,cy-r), (cx+r,cy), (cx+r,cy+r), \
                                         (cx,cy+r), (cx-r,cy+r), (cx-r,cy), (cx-r,cy-r)]]
    for p1,p2,p3 in [(0,1,2), (2,3,4), (4, 5, 6), (6, 7, 0)]:
        geo.Append( ["spline3", pts[p1], pts[p2], pts[p3]], **args)

    

def CreatePML(geo, pml_size, tol=1e-12):
    """Create a pml layer around the geometry. This function works only on convex geometries and
the highest existing domain number must be named by using the function geo.SetMaterial(domnr, name).
Points in the geometry are assumed to be the same if (pt1 - pt2).Norm() < tol.
Returned is a dict with information to create the pml layer:
  normals: A dict from the names of the linear pml domains to the normal vectors pointing inside the pml."""

    def Start(spline):
        if spline.rightdom == 0:
            return spline.StartPoint()
        return spline.EndPoint()
    def End(spline):
        if spline.rightdom == 0:
            return spline.EndPoint()
        return spline.StartPoint()

    splines = []
    for i in range(geo.GetNSplines()):
        splines.append(geo.GetSpline(i))
    border = []
    is_closed = False
    current_endpoint = None
    while not is_closed:
        for spline in splines:
            if spline.leftdom == 0 or spline.rightdom == 0:
                if current_endpoint is not None:
                    if (Start(spline)-current_endpoint).Norm() < tol:
                        border.append(spline)
                        current_endpoint = End(spline)
                        if (current_endpoint - startpoint).Norm() < tol:
                            is_closed = True
                        break
                else:
                    startpoint = Start(spline)
                    current_endpoint = End(spline)
                    border.append(spline)
                    break
        else:
            raise Exception("Couldn't find closed spline around domain")
    endpointindex_map = []
    for spline in border:
        pnt = End(spline)
        for i in range(geo.GetNPoints()):
            if (pnt - geo.GetPoint(i)).Norm() < tol:
                endpointindex_map.append(i)
                break
        else:
            raise Exception("Couldn't find endpoint of spline in geometry")
    start_ndoms = ndoms = geo.GetNDomains() + 1
    new_spline_domains = []
    normals = {}
    for i, spline in enumerate(border):
        if i == 0:
            global_start = Start(spline) + pml_size * spline.GetNormal(0)
            global_start_pnt = current_start = geo.AppendPoint(global_start[0], global_start[1])
        next_spline = border[(i+1)%len(border)]
        new_end =  End(spline) + pml_size * spline.GetNormal(1)
        spline_name = geo.GetBCName(spline.bc)
        if (new_end - global_start).Norm() < tol:
            new_spline_domains.append(ndoms)
            geo.Append(["line", current_start, global_start_pnt], bc="outer_" + spline_name, leftdomain = ndoms)
            geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms)
            geo.SetMaterial(ndoms, "pml_" + spline_name)
            normals["pml_" + spline_name] = spline.GetNormal(0)
            ndoms += 1
            break
        end = geo.AppendPoint(new_end[0], new_end[1])
        new_spline_domains.append(ndoms)
        geo.Append(["line", current_start, end], bc="outer_" + spline_name, leftdomain = ndoms)
        geo.Append(["line", end, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1)
        geo.SetMaterial(ndoms, "pml_" + spline_name)
        normals["pml_" + spline_name] = spline.GetNormal(0)
        ndoms += 1
        new_start = Start(next_spline) + pml_size * next_spline.GetNormal(0)
        if (new_start - global_start).Norm() < tol:
            geo.Append(["line", end, global_start_pnt], bc="outer", leftdomain = ndoms)
            geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms)
            geo.SetMaterial(ndoms, "pml_corner")
            ndoms += 1
            break
        if (new_end - new_start).Norm() < tol:
            current_start = end
        else:
            current_start = geo.AppendPoint(new_start[0], new_start[1])
            geo.Append(["line", end, current_start], bc="outer", leftdomain = ndoms)
            geo.Append(["line", current_start, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1)
            geo.SetMaterial(ndoms, "pml_corner")
            ndoms += 1
    for spline, domnr in zip(border, new_spline_domains):
        if spline.leftdom == 0:
            spline.leftdom = domnr
        else:
            spline.rightdom = domnr
    return {"normals" : normals}

SplineGeometry.AddCircle = lambda geo, c, r, **args : MakeCircle(geo, c, r, **args)
SplineGeometry.AddRectangle = lambda geo, p1, p2, **args : MakeRectangle(geo, p1, p2, **args)
SplineGeometry.AddSegment = lambda *args, **kwargs : SplineGeometry.Append(*args, **kwargs)
SplineGeometry.AddPoint = lambda *args, **kwargs : SplineGeometry.AppendPoint(*args, **kwargs)
SplineGeometry.CreatePML = CreatePML

__all__ = ['SplineGeometry', 'unit_square']