graphicsCurvedShapes.py

You can view and download this file on Github: graphicsCurvedShapes.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  The curved shapes of GraphicsData, for looking at them: 6-node (quadratic) triangles
  5#           ('triangles6' of a TriangleList) with quadratic edges ('edges3'), quadratic lines
  6#           (Lines with shape 'quadratic') and spheres (type 'Spheres'), each beside the flat
  7#           shape of exudyn.graphics it compares to (#2709).
  8#           Row y=0: cylinder, torus, sphere and vase of 6-node triangles - few, large elements (45 degrees),
  9#           the renderer splits them when it draws; the cylinder and the vase carry red edges3 on their rims.
 10#           Row y=-3: graphics.Cylinder, Torus, Sphere, SolidOfRevolution, LinkedCylinders (with a bore) and a hollow
 11#           Sphere - curved as well, of 6-node
 12#           triangles (the sphere with edges is a TriangleList of flat triangles); row y=-6: the same shapes as
 13#           flat triangles (graphics.Triangles6ToTriangles), with their facets.
 14#           Row y=3: quadratic lines (a circle of 4 segments, a helix), many spheres, a whole graphics.Sphere.
 15#           At x=12: a rotating rigid body with a curved cylinder, its rims and a quadratic ring.
 16#           What to look at:
 17#           - the curved rows are smooth, the flat row shows its facets; the red rims lie on the surfaces;
 18#           - key V, openGL.advanced.curvedTriangleTilingAngle: 90 shows the coarse elements, 5 a fine split,
 19#             and the change shows at once; curvedTriangleMaxTiling limits it;
 20#           - key T (face edges): the edges of the 6-node triangles are curved, not those of the split;
 21#           - CTRL+R (raytracer): the same shapes, the spheres exact; with transparent faces as well;
 22#           - SPACE: the body rotates, its surface, rims and ring move together.
 23#
 24# Author:   Johannes Gerstmayr
 25# Date:     2026-10-01
 26#
 27# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
 28#
 29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 30
 31import exudyn as exu
 32from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 33import exudyn.graphics as graphics
 34import numpy as np
 35
 36SC = exu.SystemContainer()
 37mbs = SC.AddSystem()
 38
 39#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 40#a surface F(u,v) of 6-node triangles: nu x nv quadratic patches, each split into two 6-node triangles;
 41#normalFunction gives the normals at the nodes, None uses those of the geometry; rims=True adds the
 42#boundaries v=v0 and v=v1 as quadratic edges
 43def Surface6(F, uRange, vRange, nu, nv, color, normalFunction=None, rims=False, edgeColor=graphics.color.red):
 44    us = np.linspace(uRange[0], uRange[1], 2*nu+1)
 45    vs = np.linspace(vRange[0], vRange[1], 2*nv+1)
 46    Node = lambda i, j: j*(2*nu+1) + i #node (i,j) of the grid of (2nu+1) x (2nv+1) nodes
 47    points = [F(u, v) for v in vs for u in us]
 48    triangles6 = []
 49    for j in range(nv):
 50        for i in range(nu):
 51            (a, b, c, d) = (Node(2*i, 2*j), Node(2*i+2, 2*j), Node(2*i+2, 2*j+2), Node(2*i, 2*j+2))
 52            (ab, bc, cd, da, center) = (Node(2*i+1, 2*j), Node(2*i+2, 2*j+1), Node(2*i+1, 2*j+2),
 53                                        Node(2*i, 2*j+1), Node(2*i+1, 2*j+1))
 54            triangles6 += [[a, b, c, ab, bc, center], [a, c, d, center, cd, da]] #corners, then mid nodes 01, 12, 20
 55    g = {'type':'TriangleList', 'points':np.array(points).flatten(),
 56         'colors':np.tile(color, len(points)), 'triangles6':np.array(triangles6).flatten()}
 57    if normalFunction is not None:
 58        g['normals'] = np.array([normalFunction(u, v) for v in vs for u in us]).flatten()
 59    if rims:
 60        g['edges3'] = np.array([[Node(2*i, j), Node(2*i+2, j), Node(2*i+1, j)]
 61                                for j in [0, 2*nv] for i in range(nu)]) #end points, then the mid node
 62        g['edgeColor'] = edgeColor
 63    return g
 64
 65r = 0.5
 66gCurved = []
 67#cylinder: 8 elements around, 45 degrees each
 68gCurved += [Surface6(lambda u, v: [r*np.cos(u), r*np.sin(u), v], [0, 2*np.pi], [0, 1.5], 8, 1,
 69                     graphics.color.steelblue, normalFunction=lambda u, v: [np.cos(u), np.sin(u), 0], rims=True)]
 70#torus: 8 x 4 elements
 71(R0, r0) = (0.6, 0.25)
 72gCurved += [graphics.Move(Surface6(lambda u, v: [(R0+r0*np.cos(v))*np.cos(u), (R0+r0*np.cos(v))*np.sin(u), r0*np.sin(v)],
 73                                   [0, 2*np.pi], [0, 2*np.pi], 8, 4, graphics.color.orange,
 74                                   normalFunction=lambda u, v: [np.cos(v)*np.cos(u), np.cos(v)*np.sin(u), np.sin(v)]),
 75                          [3, 0, 0.5])]
 76#sphere: 8 elements around, 4 from pole to pole
 77gCurved += [graphics.Move(Surface6(lambda u, v: [r*np.cos(v)*np.cos(u), r*np.cos(v)*np.sin(u), r*np.sin(v)],
 78                                   [0, 2*np.pi], [-np.pi/2, np.pi/2], 8, 4, graphics.color.lawngreen,
 79                                   normalFunction=lambda u, v: [np.cos(v)*np.cos(u), np.cos(v)*np.sin(u), np.sin(v)]),
 80                          [6, 0, 0.5])]
 81#vase, a solid of revolution with the radius rV(z), without normals: at each point the mean of the normals of the geometry
 82rV = lambda z: 0.3 + 0.2*np.sin(2.5*z)
 83gCurved += [graphics.Move(Surface6(lambda u, v: [rV(v)*np.cos(u), rV(v)*np.sin(u), v], [0, 2*np.pi], [0, 1.5], 8, 3,
 84                                   graphics.color.dodgerblue, rims=True), [9, 0, 0])]
 85
 86#the shapes of exudyn.graphics, curved, and flat for comparison
 87contour = [[z, rV(z)] for z in np.linspace(0, 1.5, 9)]
 88gFlat = [graphics.Cylinder(pAxis=[0, -3, 0], vAxis=[0, 0, 1.5], radius=r, color=graphics.color.steelblue, nTiles=8, addEdges=True),
 89         graphics.Torus(point=[3, -3, 0.5], axis=[0, 0, 1], radiusMajor=R0, radiusMinor=r0, color=graphics.color.orange,
 90                        nTilesMajor=12, nTilesMinor=8),
 91         graphics.Sphere(point=[6, -3, 0.5], radius=r, color=graphics.color.lawngreen, nTiles=8, addEdges=True), #with edges: a TriangleList
 92         graphics.SolidOfRevolution(pAxis=[9, -3, 0], vAxis=[0, 0, 1], contour=contour, color=graphics.color.dodgerblue,
 93                                    nTiles=12, addEdges=True),
 94         graphics.LinkedCylinders(point0=[12, -3, 0], point1=[13, -3, 0], axisCylinder=[0, 0, 0.5], radius0=0.4,
 95                                  radius1=0.25, radiusInner0=0.2, nTiles=16, color=graphics.color.steelblue, addEdges=True),
 96         graphics.Sphere(point=[15, -3, 0.5], radius=r, innerRadius=0.7*r, majorAngleMin=-0.3*np.pi, majorAngleMax=0.2*np.pi,
 97                         color=graphics.color.orange, nTiles=12, addEdges=True)]
 98gFlat += [graphics.Move(graphics.Triangles6ToTriangles(g), [0, -3, 0]) for g in gFlat]
 99
100#lines, spheres
101s = np.sqrt(0.5)
102circle = [[r, 0, 0], [r*s, r*s, 0], [0, r, 0], [-r*s, r*s, 0], [-r, 0, 0], [-r*s, -r*s, 0], [0, -r, 0], [r*s, -r*s, 0], [r, 0, 0]]
103helix = [[0.4*np.cos(t), 0.4*np.sin(t), 0.1*t] for t in np.linspace(0, 4*np.pi, 17)] #2 turns, 8 quadratic segments
104gLines = [graphics.Move(graphics.Lines(circle, color=graphics.color.red, shape='quadratic'), [0, 3, 0.5]),
105          graphics.Move(graphics.Lines(helix, color=graphics.color.blue, shape='quadratic'), [3, 3, 0]),
106          {'type':'Lines', 'points':[[2.4, 3, 0], [3.6, 3, 0]], 'colors':[graphics.color.black]*2}, #given as rows
107          graphics.Spheres(points=[[6+0.3*i, 3+0.3*j, 0.2*(i+j)] for i in range(-2, 3) for j in range(-2, 3)],
108                           radii=[0.05+0.02*(i % 4) for i in range(25)],
109                           colors=[graphics.colorList[i % 10] for i in range(25)], nTiles=16),
110          graphics.Sphere(point=[9, 3, 0.5], radius=r, color=graphics.color.grey, nTiles=16)] #a whole sphere: the type Spheres
111
112oGround = mbs.CreateGround(graphicsDataList=gCurved + gFlat + gLines +
113                           [graphics.CheckerBoard(point=[4.5, 0, -0.5], size=16, nTiles=8)])
114
115#a rotating body with a curved cylinder, its rims and a quadratic ring
116gBody = [graphics.Move(Surface6(lambda u, v: [0.3*np.cos(u), 0.3*np.sin(u), v], [0, 2*np.pi], [-0.6, 0.6], 3, 1,
117                                graphics.color.steelblue, normalFunction=lambda u, v: [np.cos(u), np.sin(u), 0], rims=True),
118                       [0, 0, 0], RotationMatrixY(0.5*np.pi)),
119         graphics.Lines([[0, 0.6*np.cos(t), 0.6*np.sin(t)] for t in np.linspace(0, 2*np.pi, 7)], color=graphics.color.red,
120                        shape='quadratic')]
121mbs.CreateRigidBody(inertia=InertiaCylinder(1000, 1.2, 0.3, axis=0), referencePosition=[12, 0, 0.5],
122                    initialAngularVelocity=[0, 0, 1], gravity=[0, 0, 0], graphicsDataList=gBody)
123mbs.Assemble()
124
125#what is in the graphics data, without a window: the curved shapes as they are, and refined once
126data = SC.renderer.GetGraphicsData()
127exu.Print('native shapes: triangles6 =', len(data['triangles6']['items']), ', lines3 =', len(data['lines3']['items']),
128          ', spheres =', len(data['spheres']['items']))
129dataFlat = SC.renderer.GetGraphicsData(flatShapes=True)
130exu.Print('refined once:  triangles =', len(dataFlat['triangles']['items']), ', lines =', len(dataFlat['lines']['items']))
131
132SC.visualizationSettings.openGL.multiSampling = 4
133SC.visualizationSettings.openGL.light0.shadow = 0.3
134SC.visualizationSettings.openGL.faceEdgesColor = [0, 0, 0, 1]
135SC.visualizationSettings.nodes.show = False
136SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
137
138simulationSettings = exu.SimulationSettings()
139simulationSettings.timeIntegration.numberOfSteps = 1000
140simulationSettings.timeIntegration.endTime = 1
141
142SC.renderer.Start()
143if SC.renderer.IsActive(): #with a window: rotate for a while, in real time
144    simulationSettings.timeIntegration.endTime = 200
145    simulationSettings.timeIntegration.numberOfSteps = 200000
146    simulationSettings.timeIntegration.realtime.active = True
147    SC.renderer.DoIdleTasks()
148mbs.SolveDynamic(simulationSettings)
149SC.renderer.DoIdleTasks()
150SC.renderer.Stop() #safely close rendering window!