mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-27 14:20:34 +05:00
179 lines
7.3 KiB
Python
179 lines
7.3 KiB
Python
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
|