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
|
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']
|
