draw curvature "m" and "C" in meshplot

pymech/meshplot.py

@@ -3,6 +3,7 @@
 import OpenGL
 from OpenGL.GL import *
 from OpenGL.GLU import *
+from math import sqrt, atan2, asin, cos, sin
 
 class MeshFrame(wx.Frame):
     """
@@ -38,6 +39,9 @@ def __init__(self,
         # create a menu bar
         #self.makeMenuBar()
 
+        # mesh display properties
+        self.curve_points = 12
+
         # data to be drawn
         self.vertices = []
         self.edges = []
@@ -155,33 +159,106 @@ def OnDraw(self, *args, **kwargs):
         self.SwapBuffers()
 
     def buildMesh(self, mesh):
-        k = 0
+        """
+        Builds the edges to be drawn based on the mesh representation.
+        """
+
+        # gives the indices of the vertices of an element in a position array
+        def vertex_indices(iface):
+            if iface == 0:
+                return (0, 0)
+            elif iface == 1:
+                return (0, -1)
+            elif iface == 2:
+                return (-1, -1)
+            else:
+                return (-1, 0)
+
+        current_point = 0
+        first_point = 0
         for el in mesh.elem:
-            self.vertices.append((
-                el.pos[0, 0, 0, 0],
-                el.pos[1, 0, 0, 0],
-                0.,
-            ))
-            self.vertices.append((
-                el.pos[0, 0, 0, -1],
-                el.pos[1, 0, 0, -1],
-                0.,
-            ))
-            self.vertices.append((
-                el.pos[0, 0, -1, -1],
-                el.pos[1, 0, -1, -1],
-                0.,
-            ))
-            self.vertices.append((
-                el.pos[0, 0, -1, 0],
-                el.pos[1, 0, -1, 0],
-                0.,
-            ))
-            self.edges.append((k, k + 1))
-            self.edges.append((k + 1, k + 2))
-            self.edges.append((k + 2, k + 3))
-            self.edges.append((k + 3, k))
-            k += 4
+            first_point = current_point
+            for iface in range(4):
+                j0, i0 = vertex_indices(iface)
+                if el.ccurv[iface] == '':
+                    self.vertices.append((
+                        el.pos[0, 0, j0, i0],
+                        el.pos[1, 0, j0, i0],
+                        0.,
+                    ))
+                    if iface < 3:
+                        next_point = current_point + 1
+                    else:
+                        next_point = first_point
+                    self.edges.append((current_point, next_point))
+                    current_point += 1
+                elif el.ccurv[iface] == 'm':
+                    # we should draw a parabola passing through the current vertex, the midpoint, and the next vertex.
+                    x0, y0 = el.pos[0:2, 0, j0, i0]
+                    xm, ym = el.curv[iface][0:2]
+                    j1, i1 = vertex_indices((iface + 1) % 4)
+                    x1, y1 = el.pos[0:2, 0, j1, i1]
+                    # quadratic Lagrange interpolation between points
+                    for ipt in range(self.curve_points):
+                        # tp varies between 0 and 1
+                        tp = ipt / self.curve_points
+                        xp = (
+                            x0 * 2 * (tp - 0.5) * (tp - 1)
+                            - xm * 4 * tp * (tp - 1)
+                            + x1 * 2 * tp * (tp - 0.5)
+                        )
+                        yp = (
+                            y0 * 2 * (tp - 0.5) * (tp - 1)
+                            - ym * 4 * tp * (tp - 1)
+                            + y1 * 2 * tp * (tp - 0.5)
+                        )
+                        self.vertices.append((xp, yp, 0))
+                        if iface == 3 and ipt == self.curve_points - 1:
+                            next_point = first_point
+                        else:
+                            next_point = current_point + 1
+                        self.edges.append((current_point, next_point))
+                        current_point += 1
+                elif el.ccurv[iface] == 'C':
+                    # draw a circle of given radius passing through the next vertex and the current one
+                    # first, find the distance between the midpoint of the segment ((x0, y0), (x1, y1)) and the center (xc, yc) of the circle
+                    radius = el.curv[iface][0]  # this can be positive or negative depending on direction
+                    x0, y0 = el.pos[0:2, 0, j0, i0]
+                    j1, i1 = vertex_indices((iface + 1) % 4)
+                    x1, y1 = el.pos[0:2, 0, j1, i1]
+                    # length of the segment
+                    ls2 = (x1 - x0) ** 2 + (y1 - y0) ** 2
+                    try:
+                        dist = radius * sqrt(1 - ls2 / (4 * radius ** 2))
+                    except ValueError:
+                        raise ValueError("the radius of the curved edge is too small")
+                    # midpoint of the edge
+                    xm = 0.5 * (x0 + x1)
+                    ym = 0.5 * (y0 + y1)
+                    # outward normal direction
+                    ls = sqrt(ls2)
+                    nx = (y1 - y0) / ls
+                    ny = -(x1 - x0) / ls
+                    # position of the centre
+                    xc = xm - nx * dist
+                    yc = ym - ny * dist
+                    # now find the range of arguments spanned by the circle arc
+                    # argument to the centre of the edge
+                    theta0 = atan2(ym - yc, xm - xc)
+                    dtheta = asin(ls / (2 * radius))
+                    theta_min = theta0 - dtheta
+                    # Now, add the points
+                    for itheta in range(self.curve_points):
+                        theta = theta_min + 2 * dtheta * itheta / self.curve_points
+                        xp = xc + abs(radius) * cos(theta)
+                        yp = yc + abs(radius) * sin(theta)
+                        self.vertices.append((xp, yp, 0))
+                        if iface == 3 and itheta == self.curve_points - 1:
+                            next_point = first_point
+                        else:
+                            next_point = current_point + 1
+                        self.edges.append((current_point, next_point))
+                        current_point += 1
 
     def setLimits(self, mesh):
         """