1from .libngpy._geom2d import SplineGeometry, Solid2d, CSG2d, Rectangle, Circle, EdgeInfo, PointInfo
2from .meshing import meshsize
3import math as math
4
5unit_square = SplineGeometry()
6_pnts = [ (0,0), (1,0), (1,1), (0,1) ]
7_lines = [ (0,1,1,"bottom"), (1,2,2,"right"), (2,3,3,"top"), (3,0,4,"left") ]
8_pnums = [unit_square.AppendPoint(*p) for p in _pnts]
9for l1,l2,bc,bcname in _lines:
10    unit_square.Append( ["line", _pnums[l1], _pnums[l2]], bc=bcname)
11
12
13def MakeRectangle (geo, p1, p2, bc=None, bcs=None, **args):
14    p1x, p1y = p1
15    p2x, p2y = p2
16    p1x,p2x = min(p1x,p2x), max(p1x, p2x)
17    p1y,p2y = min(p1y,p2y), max(p1y, p2y)
18    if not bcs: bcs=4*[bc]
19    pts = [geo.AppendPoint(*p) for p in [(p1x,p1y), (p2x, p1y), (p2x, p2y), (p1x, p2y)]]
20    for p1,p2,bc in [(0,1,bcs[0]), (1, 2, bcs[1]), (2, 3, bcs[2]), (3, 0, bcs[3])]:
21        geo.Append( ["line", pts[p1], pts[p2]], bc=bc, **args)
22
23def MakeCircle (geo, c, r, **args):
24    cx,cy = c
25    pts = [geo.AppendPoint(*p) for p in [(cx,cy-r), (cx+r,cy-r), (cx+r,cy), (cx+r,cy+r), \
26                                         (cx,cy+r), (cx-r,cy+r), (cx-r,cy), (cx-r,cy-r)]]
27    for p1,p2,p3 in [(0,1,2), (2,3,4), (4, 5, 6), (6, 7, 0)]:
28        geo.Append( ["spline3", pts[p1], pts[p2], pts[p3]], **args)
29
30
31
32def CreatePML(geo, pml_size, tol=1e-12):
33    """Create a pml layer around the geometry. This function works only on convex geometries and
34the highest existing domain number must be named by using the function geo.SetMaterial(domnr, name).
35Points in the geometry are assumed to be the same if (pt1 - pt2).Norm() < tol.
36Returned is a dict with information to create the pml layer:
37  normals: A dict from the names of the linear pml domains to the normal vectors pointing inside the pml."""
38
39    def Start(spline):
40        if spline.rightdom == 0:
41            return spline.StartPoint()
42        return spline.EndPoint()
43    def End(spline):
44        if spline.rightdom == 0:
45            return spline.EndPoint()
46        return spline.StartPoint()
47
48    splines = []
49    for i in range(geo.GetNSplines()):
50        splines.append(geo.GetSpline(i))
51    border = []
52    is_closed = False
53    current_endpoint = None
54    while not is_closed:
55        for spline in splines:
56            if spline.leftdom == 0 or spline.rightdom == 0:
57                if current_endpoint is not None:
58                    if (Start(spline)-current_endpoint).Norm() < tol:
59                        border.append(spline)
60                        current_endpoint = End(spline)
61                        if (current_endpoint - startpoint).Norm() < tol:
62                            is_closed = True
63                        break
64                else:
65                    startpoint = Start(spline)
66                    current_endpoint = End(spline)
67                    border.append(spline)
68                    break
69        else:
70            raise Exception("Couldn't find closed spline around domain")
71    endpointindex_map = []
72    for spline in border:
73        pnt = End(spline)
74        for i in range(geo.GetNPoints()):
75            if (pnt - geo.GetPoint(i)).Norm() < tol:
76                endpointindex_map.append(i)
77                break
78        else:
79            raise Exception("Couldn't find endpoint of spline in geometry")
80    start_ndoms = ndoms = geo.GetNDomains() + 1
81    new_spline_domains = []
82    normals = {}
83    duplicate_cnt = 0
84
85    for i, spline in enumerate(border):
86        if i == 0:
87            global_start = Start(spline) + pml_size * spline.GetNormal(0)
88            global_start_pnt = current_start = geo.AppendPoint(global_start[0], global_start[1])
89        next_spline = border[(i+1)%len(border)]
90        new_end =  End(spline) + pml_size * spline.GetNormal(1)
91        spline_name = geo.GetBCName(spline.bc)
92
93        if "pml_" + spline_name in normals \
94        and normals["pml_" + spline_name] != spline.GetNormal(0):
95            duplicate_cnt += 1
96            spline_name = spline_name + "_duplicate_" + str(duplicate_cnt)
97
98        if (new_end - global_start).Norm() < tol:
99            new_spline_domains.append(ndoms)
100            geo.Append(["line", current_start, global_start_pnt], bc="outer_" + spline_name, leftdomain = ndoms)
101            geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms)
102            geo.SetMaterial(ndoms, "pml_" + spline_name)
103            normals["pml_" + spline_name] = spline.GetNormal(0)
104            ndoms += 1
105            break
106        end = geo.AppendPoint(new_end[0], new_end[1])
107        new_spline_domains.append(ndoms)
108        geo.Append(["line", current_start, end], bc="outer_" + spline_name, leftdomain = ndoms)
109        geo.Append(["line", end, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1)
110        geo.SetMaterial(ndoms, "pml_" + spline_name)
111        normals["pml_" + spline_name] = spline.GetNormal(0)
112        ndoms += 1
113        new_start = Start(next_spline) + pml_size * next_spline.GetNormal(0)
114        if (new_start - global_start).Norm() < tol:
115            geo.Append(["line", end, global_start_pnt], bc="outer", leftdomain = ndoms)
116            geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms)
117            geo.SetMaterial(ndoms, "pml_corner")
118            ndoms += 1
119            break
120        if (new_end - new_start).Norm() < tol:
121            current_start = end
122        else:
123            current_start = geo.AppendPoint(new_start[0], new_start[1])
124            geo.Append(["line", end, current_start], bc="outer", leftdomain = ndoms)
125            geo.Append(["line", current_start, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1)
126            geo.SetMaterial(ndoms, "pml_corner")
127            ndoms += 1
128    for spline, domnr in zip(border, new_spline_domains):
129        if spline.leftdom == 0:
130            spline.leftdom = domnr
131        else:
132            spline.rightdom = domnr
133    return {"normals" : normals}
134
135SplineGeometry.AddCircle = lambda geo, c, r, **args : MakeCircle(geo, c, r, **args)
136SplineGeometry.AddRectangle = lambda geo, p1, p2, **args : MakeRectangle(geo, p1, p2, **args)
137SplineGeometry.AddSegment = lambda *args, **kwargs : SplineGeometry.Append(*args, **kwargs)
138SplineGeometry.AddPoint = lambda *args, **kwargs : SplineGeometry.AppendPoint(*args, **kwargs)
139SplineGeometry.CreatePML = CreatePML
140
141bc = lambda s : EdgeInfo(bc=s)
142maxh = lambda h : EdgeInfo(maxh=h)
143def cp(p_or_px, py_or_none = None):
144    if py_or_none is None:
145        return EdgeInfo(control_point=p)
146    else:
147        return EdgeInfo(control_point=(p_or_px,py_or_none))
148
149
150def Ellipse(center, a, b, bc="ellipse", mat="ellipse"):
151    """Creates ellipse centered at point center with principle axis a and b.
152
153    Parameters
154    ---------
155    center : Vec2
156      center of ellipse
157    a : Vec2
158      first principle axis, needs to be perpendicular to b
159    b : Vec2
160      second principle axis, needs to be perpendicular to a
161    bc : string
162      boundary name
163    mat : string
164      material name
165    """
166    if abs(a[0]*b[0] + a[1]*b[1]) > 1e-12:
167        raise Exception("In Ellipse: principle axis a and b are not perpendicular")
168
169    ellipse = Circle( center=(0,0), radius=1.0, mat=mat, bc=bc )
170
171    alpha = math.pi/2-math.atan2(a[0],a[1])
172    r_a = math.sqrt(a[0]**2+a[1]**2)
173    r_b = math.sqrt(b[0]**2+b[1]**2)
174    ellipse.Scale( (r_a,r_b) )
175    ellipse.Rotate( alpha/math.pi*180, center=(0,0) )
176    ellipse.Move( center )
177
178    return ellipse
179