1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
|
from .libngpy._geom2d import SplineGeometry, Solid2d, CSG2d, Rectangle, Circle, EdgeInfo, PointInfo
from .meshing import meshsize
import math as math
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 = {}
duplicate_cnt = 0
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 "pml_" + spline_name in normals \
and normals["pml_" + spline_name] != spline.GetNormal(0):
duplicate_cnt += 1
spline_name = spline_name + "_duplicate_" + str(duplicate_cnt)
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
bc = lambda s : EdgeInfo(bc=s)
maxh = lambda h : EdgeInfo(maxh=h)
def cp(p_or_px, py_or_none = None):
if py_or_none is None:
return EdgeInfo(control_point=p)
else:
return EdgeInfo(control_point=(p_or_px,py_or_none))
def Ellipse(center, a, b, bc="ellipse", mat="ellipse"):
"""Creates ellipse centered at point center with principle axis a and b.
Parameters
---------
center : Vec2
center of ellipse
a : Vec2
first principle axis, needs to be perpendicular to b
b : Vec2
second principle axis, needs to be perpendicular to a
bc : string
boundary name
mat : string
material name
"""
if abs(a[0]*b[0] + a[1]*b[1]) > 1e-12:
raise Exception("In Ellipse: principle axis a and b are not perpendicular")
ellipse = Circle( center=(0,0), radius=1.0, mat=mat, bc=bc )
alpha = math.pi/2-math.atan2(a[0],a[1])
r_a = math.sqrt(a[0]**2+a[1]**2)
r_b = math.sqrt(b[0]**2+b[1]**2)
ellipse.Scale( (r_a,r_b) )
ellipse.Rotate( alpha/math.pi*180, center=(0,0) )
ellipse.Move( center )
return ellipse
|
