Skip to content

Commit

Permalink
bug fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@cgcgcg
cgcgcg committed Nov 27, 2023
1 parent 9d60c74 commit aa50efc
Show file tree
Hide file tree
Showing 2 changed files with 73 additions and 77 deletions.
148 changes: 72 additions & 76 deletions fem/PyNucleus_fem/mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,13 +108,20 @@ def pacman(h=0.1, **kwargs):
return mesh


def uniformSquare(N, M=None, ax=0, ay=0, bx=1, by=1, crossed=False, preserveLinesHorizontal=[], preserveLinesVertical=[]):
if M is None:
M = max(int(np.around((by-ay)/(bx-ax)))*N, 2)
assert N >= 2
assert M >= 2
xVals = np.linspace(ax, bx, N)
yVals = np.linspace(ay, by, M)
def uniformSquare(N=2, M=None, ax=0, ay=0, bx=1, by=1, crossed=False, preserveLinesHorizontal=[], preserveLinesVertical=[], xVals=None, yVals=None):
if xVals is None:
assert N is not None
assert N >= 2
xVals = np.linspace(ax, bx, N)
else:
N = xVals.shape[0]
if yVals is None:
if M is None:
M = max(int(np.around((by-ay)/(bx-ax)))*N, 2)
assert M >= 2
yVals = np.linspace(ay, by, M)
else:
M = yVals.shape[0]
x, y = np.meshgrid(xVals, yVals)
for yVal in preserveLinesHorizontal:
assert (yVals-yVal).min() < 1e-10
Expand Down Expand Up @@ -310,96 +317,85 @@ def squareWithInteractions(ax, ay, bx, by,
if not uniform:
from . meshConstruction import (circularSegment,
line,
polygon,
transformationRestriction)
if h is None:
h = horizon
bottomLeft = np.array([ax, ay])
bottomRight = np.array([bx, ay])
topRight = np.array([bx, by])
topLeft = np.array([ax, by])

horizontalOffset = np.array([horizon, 0.])
verticalOffset = np.array([0., horizon])

center = np.array([(ax+bx)/2, (ay+by)/2])

numPointsPerUnitLength = int(np.ceil(1/h))

assert len(preserveLinesVertical) == 0 or len(preserveLinesHorizontal) == 0
if len(preserveLinesVertical)+len(preserveLinesHorizontal) > 0:
preserve = preserveLinesVertical+preserveLinesHorizontal

c1 = circularSegment(bottomLeft, horizon, np.pi, 3/2*np.pi, numPointsPerUnitLength)

x1 = preserve[0]
c2 = line((ax, ay), (x1, ay))
for k in range(len(preserve)-1):
x1 = preserve[k]
x2 = preserve[k+1]
c2 = c2+line((x1, ay), (x2, ay))
x2 = preserve[-1]
c2 = c2+line((x2, ay), (bx, ay))
c1 = c1 + (c2+(0., -horizon))
else:
c1 = circularSegment(bottomLeft, horizon, np.pi, 3/2*np.pi, numPointsPerUnitLength)
c2 = line((ax, ay), (bx, ay))
c1 = c1 + (c2+(0., -horizon))
c3 = line((ax, ay), (ax, ay-horizon))
c4 = line((ax, ay), (ax-horizon, ay))
c = c1+c2+c3+c4

frame = (c + (c*(center, np.pi/2)) + (c*(center, np.pi)) + (c*(center, -np.pi/2)))

if len(preserveLinesVertical) > 0:
d = line((0, ay-horizon), (0, ay))
x1 = preserve[0]
d = d + line((0, ay), (0, x1))
for k in range(len(preserve)-1):
x1 = preserve[k]
x2 = preserve[k+1]
d += line((0, x1), (0, x2))
x2 = preserve[-1]
d = d + line((0, x2), (0, by))
d = d + line((0, by), (0, by+horizon))
for x in preserveLinesVertical:
assert ax <= x <= bx
frame += (d+(x, 0.))
if len(preserveLinesHorizontal) > 0:
d = line((ax-horizon, 0), (ax, 0))+line((ax, 0), (bx, 0))+line((bx, 0), (bx+horizon, 0))

d = line((ax-horizon, 0), (ax, 0))
x1 = preserve[0]
d = d + line((ax, 0), (x1, 0))
for k in range(len(preserve)-1):
x1 = preserve[k]
x2 = preserve[k+1]
d += line((x1, 0), (x2, 0))
x2 = preserve[-1]
d = d + line((x2, 0), (bx, 0))
d = d + line((bx, 0), (bx+horizon, 0))

for y in preserveLinesHorizontal:
assert ay <= y <= by
frame += (d+(0, y))

lineHorizontal = polygon([(0., 0.)] + [(p-ax, 0.) for p in preserveLinesVertical] + [(bx-ax, 0.)], doClose=False)
lineVertical = polygon([(0., 0.)] + [(0., p-ay) for p in preserveLinesHorizontal] + [(0., by-ay)], doClose=False)

d1 = (circularSegment(bottomLeft, horizon, np.pi, 1.5*np.pi, numPointsPerUnitLength) +
line(bottomLeft, bottomLeft-horizontalOffset) +
line(bottomLeft, bottomLeft-verticalOffset) +
(lineHorizontal+bottomLeft)+
(lineHorizontal+(bottomLeft-verticalOffset)))

d2 = (circularSegment(bottomRight, horizon, -0.5*np.pi, 0., numPointsPerUnitLength) +
line(bottomRight, bottomRight+horizontalOffset) +
line(bottomRight, bottomRight-verticalOffset) +
(lineVertical+(bottomRight+horizontalOffset)) +
(lineVertical+bottomRight))

d3 = (circularSegment(topRight, horizon, 0, 0.5*np.pi, numPointsPerUnitLength) +
line(topRight, topRight+horizontalOffset) +
line(topRight, topRight+verticalOffset) +
(lineHorizontal+topLeft) +
(lineHorizontal+(topLeft+verticalOffset)))

d4 = (circularSegment(topLeft, horizon, 0.5*np.pi, np.pi, numPointsPerUnitLength) +
line(topLeft, topLeft-horizontalOffset) +
line(topLeft, topLeft+verticalOffset) +
(lineVertical+bottomLeft) +
(lineVertical+(bottomLeft-horizontalOffset)))

frame = d1 + d2 + d3 + d4

frame.holes.append(center)

if innerRadius > 0:
frame += transformationRestriction(circularSegment(center, innerRadius, 0, 2*np.pi, numPointsPerUnitLength),
center-(innerRadius, innerRadius),
center+(innerRadius, innerRadius))
mesh = frame.mesh(max_volume=h**2, min_angle=30, **kwargs)
else:
frame.holes.append(center)
mesh = frame.mesh(max_volume=h**2, min_angle=30, **kwargs)

eps = 1e-10
N1 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 0]-ax) < eps,
np.logical_and(mesh.vertices_as_array[:, 1] >= ay-eps,
mesh.vertices_as_array[:, 1] <= by+eps)).sum()
N2 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 0]-bx) < eps,
np.logical_and(mesh.vertices_as_array[:, 1] >= ay-eps,
mesh.vertices_as_array[:, 1] <= by+eps)).sum()
M1 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 1]-ay) < eps,
idx1 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 0]-ax) < eps,
np.logical_and(mesh.vertices_as_array[:, 1] >= ay-eps,
mesh.vertices_as_array[:, 1] <= by+eps))
idx2 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 0]-bx) < eps,
np.logical_and(mesh.vertices_as_array[:, 1] >= ay-eps,
mesh.vertices_as_array[:, 1] <= by+eps))
yVals1 = np.sort(mesh.vertices_as_array[idx1, 1])
yVals2 = np.sort(mesh.vertices_as_array[idx2, 1])
assert np.allclose(yVals1, yVals2), (yVals1, yVals2)

idx3 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 1]-ay) < eps,
np.logical_and(mesh.vertices_as_array[:, 0] >= ax-eps,
mesh.vertices_as_array[:, 0] <= bx+eps)).sum()
M2 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 1]-by) < eps,
mesh.vertices_as_array[:, 0] <= bx+eps))
idx4 = np.logical_and(np.absolute(mesh.vertices_as_array[:, 1]-by) < eps,
np.logical_and(mesh.vertices_as_array[:, 0] >= ax-eps,
mesh.vertices_as_array[:, 0] <= bx+eps)).sum()
assert N1 == N2, (N1, N2)
assert M1 == M2, (M1, M2)
mesh2 = uniformSquare(N=N1, M=M1, ax=ax, ay=ay, bx=bx, by=by)
mesh.vertices_as_array[:, 0] <= bx+eps))
xVals3 = np.sort(mesh.vertices_as_array[idx3, 0])
xVals4 = np.sort(mesh.vertices_as_array[idx4, 0])
assert np.allclose(xVals3, xVals4), (xVals3, xVals4)
mesh2 = uniformSquare(ax=ax, ay=ay, bx=bx, by=by, xVals=xVals3, yVals=yVals1)
mesh2.plot()
mesh = snapMeshes(mesh, mesh2)

location = uninitialized((mesh.num_vertices), dtype=INDEX)
Expand Down
2 changes: 1 addition & 1 deletion fem/PyNucleus_fem/meshConstruction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ def __init__(self, points, facets, holes=[]):
self.meshTransformations = []

def __add__(self, other):
if isinstance(other, tuple):
if isinstance(other, (tuple, np.ndarray)):
newPoints = [(other[0]+p[0], other[1]+p[1]) for p in self.points]
newHoles = [(other[0]+p[0], other[1]+p[1]) for p in self.holes]
newSegment = segment(newPoints, self.facets, newHoles)
Expand Down

0 comments on commit aa50efc

Please sign in to comment.