From 63ff4e957b27b9987703e8e6c034e00e1a608ebc Mon Sep 17 00:00:00 2001 From: Yeongdo Park Date: Wed, 28 Dec 2022 19:12:00 +0900 Subject: [PATCH] paraview programmable sources --- annulus-R19.py | 86 ++++++++++++++++++++++++++++++++++++++++++++++++++ annulus.py | 84 ++++++++++++++++++++++++++++++++++++++++++++++++ anulus.py | 69 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 239 insertions(+) create mode 100644 annulus-R19.py create mode 100644 annulus.py create mode 100644 anulus.py diff --git a/annulus-R19.py b/annulus-R19.py new file mode 100644 index 0000000..b6eb518 --- /dev/null +++ b/annulus-R19.py @@ -0,0 +1,86 @@ +import numpy as np + +# Number of circumferential faces (should be even number) +n_faces = 40 + +# Outer Radius +ro = 19. # mm + +# Width +w = 3.5 # mm + +# Circumference / Width +aspect = 2 * np.pi * ro / w +print(aspect, np.round(aspect/2)*2) +n_faces = np.round(aspect/2).astype(np.int)*2 + +golden = (1+np.sqrt(5))/2 +# Ratio - Circumferential faces upper and lower side +ratio = 1/(1+golden) # 1/2 # 2/(1+np.sqrt(5))# 0 ~ 1 +ratio = np.max([0., ratio]) +ratio = np.min([1., ratio]) + +offset = 2 * ratio - 1. + +pdo = self.GetPolyDataOutput() + +dtheta = 2 * np.pi / n_faces + +multiplier = np.arange(n_faces, dtype=np.double) +multiplier[1::2] = multiplier[1::2] + offset + +theta = multiplier * dtheta + +x = ro * np.sin(theta) +y = ro * np.cos(theta) + +torsion = dtheta + 0.0*dtheta + +sin_dtheta = np.sin(torsion) +cos_dtheta = np.cos(torsion) + +# rotation_mat = np.array([cos_dtheta, - sin_dtheta, sin_dtheta, cos_dtheta]).reshape((2,2)) + +x_ = cos_dtheta * x - sin_dtheta * y +y_ = sin_dtheta * x + cos_dtheta * y + +x1 = np.concatenate((x_[1:], x_[:1])) +y1 = np.concatenate((y_[1:], y_[:1])) + +nodes = vtk.vtkPoints() +for i, (xi, yi) in enumerate(zip(x, y)): + nodes.InsertNextPoint(xi, yi, 0.0) + +for i, (xi, yi) in enumerate(zip(x1, y1)): + nodes.InsertNextPoint(xi, yi, w) + +elements = vtk.vtkCellArray() +for i in range(n_faces): + + # Create the object that stores the list of nodeIDs + elementIdList0 = vtk.vtkIdList() + + # Insert the nodeIDs in the appropriate order. + i_next = (i + 1) % n_faces + elementIdList0.InsertNextId(i) + elementIdList0.InsertNextId(i+n_faces) + elementIdList0.InsertNextId(i_next+n_faces) + elementIdList0.InsertNextId(i_next) + + elements.InsertNextCell(elementIdList0) + +elementIdListTop = vtk.vtkIdList() +elementIdListBottom = vtk.vtkIdList() +for i in range(n_faces): + + # Insert the nodeIDs in the appropriate order. + elementIdListTop.InsertNextId(i) + elementIdListBottom.InsertNextId(i+n_faces) + +elements.InsertNextCell(elementIdListTop) +elements.InsertNextCell(elementIdListBottom) + +# Specify that the variable 'nodes' as the points in the polydata +pdo.SetPoints(nodes) +# Similarly, specify that the variable 'elements' as the polys in polydata +pdo.SetPolys(elements) \ No newline at end of file diff --git a/annulus.py b/annulus.py new file mode 100644 index 0000000..1b570d4 --- /dev/null +++ b/annulus.py @@ -0,0 +1,84 @@ +import numpy as np + +# Number of circumferential faces (should be even number) +n_faces = 36 + +# Outer Radius +ro = 22. # mm + +# Width +w = 2.8 # mm + +# Circumference / Width +aspect = 2 * np.pi * ro / w +print(aspect) + +# Ratio - Circumferential faces upper and lower side +ratio = 2/(1+np.sqrt(5))# 0 ~ 1 +ratio = np.max([0., ratio]) +ratio = np.min([1., ratio]) + +offset = 2 * ratio - 1. + +pdo = self.GetPolyDataOutput() + +dtheta = 2 * np.pi / n_faces + +multiplier = np.arange(n_faces, dtype=np.double) +multiplier[1::2] = multiplier[1::2] + offset + +theta = multiplier * dtheta + +x = ro * np.sin(theta) +y = ro * np.cos(theta) + +torsion = dtheta + 0.0*dtheta + +sin_dtheta = np.sin(torsion) +cos_dtheta = np.cos(torsion) + +# rotation_mat = np.array([cos_dtheta, - sin_dtheta, sin_dtheta, cos_dtheta]).reshape((2,2)) + +x_ = cos_dtheta * x - sin_dtheta * y +y_ = sin_dtheta * x + cos_dtheta * y + +x1 = np.concatenate((x_[1:], x_[:1])) +y1 = np.concatenate((y_[1:], y_[:1])) + +nodes = vtk.vtkPoints() +for i, (xi, yi) in enumerate(zip(x, y)): + nodes.InsertNextPoint(xi, yi, 0.0) + +for i, (xi, yi) in enumerate(zip(x1, y1)): + nodes.InsertNextPoint(xi, yi, w) + +elements = vtk.vtkCellArray() +for i in range(n_faces): + + # Create the object that stores the list of nodeIDs + elementIdList0 = vtk.vtkIdList() + + # Insert the nodeIDs in the appropriate order. + i_next = (i + 1) % n_faces + elementIdList0.InsertNextId(i) + elementIdList0.InsertNextId(i+n_faces) + elementIdList0.InsertNextId(i_next+n_faces) + elementIdList0.InsertNextId(i_next) + + elements.InsertNextCell(elementIdList0) + +elementIdListTop = vtk.vtkIdList() +elementIdListBottom = vtk.vtkIdList() +for i in range(n_faces): + + # Insert the nodeIDs in the appropriate order. + elementIdListTop.InsertNextId(i) + elementIdListBottom.InsertNextId(i+n_faces) + +elements.InsertNextCell(elementIdListTop) +elements.InsertNextCell(elementIdListBottom) + +# Specify that the variable 'nodes' as the points in the polydata +pdo.SetPoints(nodes) +# Similarly, specify that the variable 'elements' as the polys in polydata +pdo.SetPolys(elements) \ No newline at end of file diff --git a/anulus.py b/anulus.py new file mode 100644 index 0000000..3d41a3e --- /dev/null +++ b/anulus.py @@ -0,0 +1,69 @@ +import numpy as np + +# Number of circumferential faces (should be even number) +n_faces = 40 + +# Outer Radius +ro = 27. # mm + +# Width +w = 4.2 # mm + +# Ratio +ratio = 2/(1+np.sqrt(5))# 0 ~ 1 +ratio = np.max([0., ratio]) +ratio = np.min([1., ratio]) + +offset = 2 * ratio - 1. + +pdo = self.GetPolyDataOutput() + +dtheta = 2 * np.pi / n_faces + +multiplier = np.arange(n_faces, dtype=np.double) +multiplier[1::2] = multiplier[1::2] + offset + +theta = multiplier * dtheta + +x = ro * np.sin(theta) +y = ro * np.cos(theta) + +torsion = dtheta + 0.2*dtheta + +sin_dtheta = np.sin(torsion) +cos_dtheta = np.cos(torsion) + +# rotation_mat = np.array([cos_dtheta, - sin_dtheta, sin_dtheta, cos_dtheta]).reshape((2,2)) + +x_ = cos_dtheta * x - sin_dtheta * y +y_ = sin_dtheta * x + cos_dtheta * y + +x1 = np.concatenate((x_[1:], x_[:1])) +y1 = np.concatenate((y_[1:], y_[:1])) + +nodes = vtk.vtkPoints() +for i, (xi, yi) in enumerate(zip(x, y)): + nodes.InsertNextPoint(xi, yi, 0.0) + +for i, (xi, yi) in enumerate(zip(x1, y1)): + nodes.InsertNextPoint(xi, yi, w) + +elements = vtk.vtkCellArray() +for i in range(n_faces): + + # Create the object that stores the list of nodeIDs + elementIdList0 = vtk.vtkIdList() + + # Insert the nodeIDs in the appropriate order. + i_next = (i + 1) % n_faces + elementIdList0.InsertNextId(i) + elementIdList0.InsertNextId(i+n_faces) + elementIdList0.InsertNextId(i_next+n_faces) + elementIdList0.InsertNextId(i_next) + + elements.InsertNextCell(elementIdList0) + +# Specify that the variable 'nodes' as the points in the polydata +pdo.SetPoints(nodes) +# Similarly, specify that the variable 'elements' as the polys in polydata +pdo.SetPolys(elements) \ No newline at end of file