Skip to content

Commit

Permalink
remove test plots from module
Browse files Browse the repository at this point in the history
  • Loading branch information
pbrubeck committed Oct 24, 2023
1 parent 58a10c6 commit fab288d
Showing 1 changed file with 0 additions and 63 deletions.
63 changes: 0 additions & 63 deletions FIAT/recursive_points.py
Original file line number Diff line number Diff line change
Expand Up @@ -101,66 +101,3 @@ def make_points(self, ref_el, dim, entity_id, order):
return self.recursive_points(ref_el.get_vertices(), order, 1)
else:
raise ValueError("illegal dimension")


if __name__ == "__main__":
from FIAT import reference_element, quadrature
from matplotlib import pyplot as plt

interval = reference_element.UFCInterval()
cg = lambda n: numpy.linspace(0.0, 1.0, n + 1)
dg = lambda n: numpy.linspace(1.0/(n+2.0), (n+1.0)/(n+2.0), n + 1)
gll = lambda n: quadrature.GaussLegendreQuadratureLineRule(interval, n + 1).get_points()
gl = lambda n: quadrature.GaussLobattoLegendreQuadratureLineRule(interval, n + 1).get_points() if n else None

order = 11
cg_points = RecursivePointSet(gll)
dg_points = RecursivePointSet(gl)

ref_el = reference_element.ufc_simplex(2)
h = numpy.sqrt(3)
s = 2*h/3
ref_el.vertices = [(0, s), (-1.0, s-h), (1.0, s-h)]
x = numpy.array(ref_el.vertices + ref_el.vertices[:1])
plt.plot(x[:, 0], x[:, 1], "k")

for d in range(1, 4):
assert cg_points.make_points(reference_element.ufc_simplex(d), d, 0, d) == []

topology = ref_el.get_topology()
for dim in topology:
for entity in topology[dim]:
pts = cg_points.make_points(ref_el, dim, entity, order)
if len(pts):
x = numpy.array(pts)
for r in range(1, 3):
th = r * (2*numpy.pi)/3
ct = numpy.cos(th)
st = numpy.sin(th)
Q = numpy.array([[ct, -st], [st, ct]])
x = numpy.concatenate([x, numpy.dot(x, Q)])
plt.scatter(x[:, 0], x[:, 1])

x0 = 2.0
h = -h
s = 2*h/3
ref_el = reference_element.ufc_simplex(2)
ref_el.vertices = [(x0, s), (x0-1.0, s-h), (x0+1.0, s-h)]

x = numpy.array(ref_el.vertices + ref_el.vertices[:1])
d = len(ref_el.vertices)
x0 = sum(x[:d])/d
plt.plot(x[:, 0], x[:, 1], "k")

pts = dg_points.recursive_points(ref_el.vertices, order)
x = numpy.array(pts)
for r in range(1, 3):
th = r * (2*numpy.pi)/3
ct = numpy.cos(th)
st = numpy.sin(th)
Q = numpy.array([[ct, -st], [st, ct]])
x = numpy.concatenate([x, numpy.dot(x-x0, Q)+x0])
plt.scatter(x[:, 0], x[:, 1])

plt.gca().set_aspect('equal', 'box')
plt.show()

0 comments on commit fab288d

Please sign in to comment.