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