Archive for November, 2009

Chain model – simple simulation

Sunday, November 15th, 2009

Wrote a simple script that simulates chain models in Microstation. A series of chain segments of a specific length is defined. Some knots are fix others can move.

 

‘Simulation of a chain model
‘Sebastian Gmelin
’15.11.2009

Option Explicit

Const gravity As Double = 1 / 5
Const springforce As Double = 1 / 8
Const attrForce As Double = 3

Public Type sphere
    number As Long
    radius As Double
    center As Point3d
    newCenter As Point3d
    sElements(3) As Element
    numNeighbours As Long
    Neighbours(10) As Long
    connections(10) As Element
    moveable As Boolean
End Type

Dim spheres() As sphere

Function createSphere(pos As Point3d, rad As Double, number As Long, moveable As Boolean) As sphere
   
    Dim elem As Element
    Dim mat As Matrix3d
   
    createSphere.center = pos
    createSphere.newCenter = pos
    createSphere.radius = rad
    createSphere.moveable = moveable
    createSphere.number = number
    createSphere.numNeighbours = 0
   
    Set elem = CreateLineElement2(Nothing, pos, pos)
    elem.Color = 3
    If moveable Then elem.Color = 1
    elem.LineWeight = 8
    elem.LineStyle = ActiveDesignFile.LineStyles(1)
    ActiveModelReference.AddElement elem
    elem.Redraw
    Set createSphere.sElements(0) = elem
   
    mat = Matrix3dFromAxisAndRotationAngle(2, 0)
    Set elem = CreateEllipseElement2(Nothing, pos, rad, rad, mat, msdFillModeNotFilled)
    elem.Color = 6
    elem.LineWeight = 0
    elem.LineStyle = ActiveDesignFile.LineStyles(2)
    ActiveModelReference.AddElement elem
    elem.Redraw
    Set createSphere.sElements(1) = elem
   
    mat = Matrix3dFromAxisAndRotationAngle(1, Pi / 2)
    Set elem = CreateEllipseElement2(Nothing, pos, rad, rad, mat, msdFillModeNotFilled)
    elem.Color = 6
    elem.LineWeight = 0
    elem.LineStyle = ActiveDesignFile.LineStyles(2)
    ActiveModelReference.AddElement elem
    elem.Redraw
    Set createSphere.sElements(2) = elem
   
    mat = Matrix3dFromAxisAndRotationAngle(0, Pi / 2)
    Set elem = CreateEllipseElement2(Nothing, pos, rad, rad, mat, msdFillModeNotFilled)
    elem.Color = 6
    elem.LineWeight = 0
    elem.LineStyle = ActiveDesignFile.LineStyles(2)
    ActiveModelReference.AddElement elem
    elem.Redraw
    Set createSphere.sElements(3) = elem
   
End Function

Sub moveSphere(s As sphere, vec As Point3d)

    Dim i As Long
   
    s.center = Point3dAdd(s.center, vec)
   
    For i = 0 To UBound(s.sElements)
        Call s.sElements(i).Move(vec)
        s.sElements(i).Rewrite
        s.sElements(i).Redraw
        ‘ActiveDesignFile.Views(1).Redraw
    Next i
   
End Sub

Sub moveConnections()

    Dim s As Long
    Dim c As Long
   
    For s = 0 To UBound(spheres)
        For c = 0 To spheres(s).numNeighbours – 1
            Call spheres(s).connections(c).AsLineElement.ModifyVertex(0, spheres(s).center)
            Call spheres(s).connections(c).AsLineElement.ModifyVertex(1, spheres(spheres(s).Neighbours(c)).center)
            spheres(s).connections(c).Rewrite
            spheres(s).connections(c).Redraw
        Next c
    Next s

End Sub

Sub makeSphereFix(s As Long)

    spheres(s).moveable = False
    spheres(s).sElements(0).Color = 3
    spheres(s).sElements(0).Rewrite
    spheres(s).sElements(0).Redraw

End Sub

Sub makeSphereMoveAble(s As Long)

    spheres(s).moveable = True
    spheres(s).sElements(0).Color = 1
    spheres(s).sElements(0).Rewrite
    spheres(s).sElements(0).Redraw

End Sub

Sub drawConnections(s As Long)

    Dim elem As Element
    Dim points(1) As Point3d
    Dim i As Long
   
    For i = 0 To spheres(s).numNeighbours – 1
        points(0) = spheres(s).center
        points(1) = spheres(spheres(s).Neighbours(i)).center
        Set elem = CreateLineElement1(Nothing, points)
        elem.Color = 2
        elem.LineWeight = 2
        elem.LineStyle = ActiveDesignFile.LineStyles(1)
        ActiveModelReference.AddElement elem
        elem.Redraw
        Set spheres(s).connections(i) = elem
    Next i

End Sub

Sub calculateSphereMove(s As Long)

    Dim i As Long
    Dim k As Long
    Dim direction As Point3d
    Dim dist As Double
    Dim length As Double
    Dim grav As Boolean
    Dim vec As Point3d
    Dim tension As Double
   
    If spheres(s).moveable Then
    tension = 0
        For i = 0 To spheres(s).numNeighbours – 1
            direction = Point3dSubtract(spheres(spheres(s).Neighbours(i)).center, spheres(s).newCenter)
            dist = Abs(Point3dMagnitude(direction))
            ‘Debug.Print (“distance ” + CStr(dist))
            length = spheres(s).radius + spheres(spheres(s).Neighbours(i)).radius
            ‘Debug.Print (“length ” + CStr(length))
            tension = tension + dist – length
        Next i
        tension = tension / (spheres(s).numNeighbours)
        grav = True
        If tension > 0 Then
            spheres(s).newCenter = Point3dAdd(spheres(s).newCenter, Point3dFromXYZ(0, 0, -0.1 * gravity))
        Else
            spheres(s).newCenter = Point3dAdd(spheres(s).newCenter, Point3dFromXYZ(0, 0, -1 * gravity))
        End If
        vec = Point3dFromXYZ(0, 0, 0)
        For i = 0 To spheres(s).numNeighbours – 1
            direction = Point3dSubtract(spheres(spheres(s).Neighbours(i)).center, spheres(s).newCenter)
            dist = Abs(Point3dMagnitude(direction))
            ‘Debug.Print (“distance ” + CStr(dist))
            length = spheres(s).radius + spheres(spheres(s).Neighbours(i)).radius
            ‘Debug.Print (“length ” + CStr(length))
            If (dist – length) > 0 Then
                grav = False
                direction = Point3dScale(direction, attrForce * springforce * (dist – length) / (dist))
            Else
                direction = Point3dScale(direction, springforce * (dist – length) / (dist))
            End If
            vec = Point3dAdd(vec, direction)
        Next i
        spheres(s).newCenter = Point3dAdd(spheres(s).newCenter, vec)
    End If

End Sub

Sub moveSpheres()

    Dim i As Long
    Dim s As Long
    Dim vec As Point3d
   
    For s = 0 To UBound(spheres)
        vec = Point3dSubtract(spheres(s).newCenter, spheres(s).center)
        spheres(s).center = spheres(s).newCenter
   
        For i = 0 To UBound(spheres(s).sElements)
            Call spheres(s).sElements(i).Move(vec)
            spheres(s).sElements(i).Rewrite
            spheres(s).sElements(i).Redraw
            ‘ActiveDesignFile.Views(1).Redraw
        Next i
    Next s

End Sub

 

Sub Main()

    Dim i As Long
    Dim k As Long
    Dim s As Long
   
    Randomize

    ReDim spheres(99)
    For i = 0 To 9
        For k = 0 To 9
            spheres(i * 10 + k) = createSphere(Point3dFromXYZ(i, k, 0), 0.65, i * 10 + k, True)
        Next k
    Next i
   
    For i = 0 To 9
        For k = 0 To 9
            ‘spheres(i * 10 + k) = createSphere(Point3dFromXYZ(i, k, 0), rnd * 0.2 + 0.5, i * 10 + k, True)
            If i < 9 Then
                spheres(i * 10 + k).Neighbours(spheres(i * 10 + k).numNeighbours) = i * 10 + k + 10
                spheres(i * 10 + k).numNeighbours = spheres(i * 10 + k).numNeighbours + 1
            End If
            If k < 9 Then
                spheres(i * 10 + k).Neighbours(spheres(i * 10 + k).numNeighbours) = i * 10 + k + 1
                spheres(i * 10 + k).numNeighbours = spheres(i * 10 + k).numNeighbours + 1
            End If
            If i > 0 Then
                spheres(i * 10 + k).Neighbours(spheres(i * 10 + k).numNeighbours) = i * 10 + k – 10
                spheres(i * 10 + k).numNeighbours = spheres(i * 10 + k).numNeighbours + 1
            End If
            If k > 0 Then
                spheres(i * 10 + k).Neighbours(spheres(i * 10 + k).numNeighbours) = i * 10 + k – 1
                spheres(i * 10 + k).numNeighbours = spheres(i * 10 + k).numNeighbours + 1
            End If
            Call drawConnections(i * 10 + k)
        Next k
    Next i
   
    ‘Call makeSphereFix(0)
    Call makeSphereFix(9)
    ‘Call makeSphereFix(55)
    Call makeSphereFix(90)
    Call makeSphereFix(99)
   
    For i = 0 To 75
‘        Call moveSphere(spheres(Round(Rnd * 99)), Point3dFromXYZ(0, 0, Rnd * 0.5 – 0.25))
‘        Call moveConnections
        For s = 0 To UBound(spheres)
            Call calculateSphereMove(s)
        Next s
        Call moveSpheres
        Call moveConnections
        ActiveDesignFile.Views(1).Redraw
    Next i
   
End Sub

Share

Hyperbolic Paraboloids – math model

Wednesday, November 11th, 2009

Hyperbolic paraboloids (HP) are discribed by the formula
\frac{x^2}{a^2} - \frac{y^2}{b^2} - 2z = 0
The below script creates HPs in Microstation within a certain range of  x and y coordinates. The dashed lines indicate the regular grid on coordinate system. These are of parabolic shape. The thick diagonal lines are all straight. They are the reason why this shape is interesting for construction, as it can be discribed as a series of straight lines. Intersecting the surface with a horizontal plane creates hyperbolas (blue curves).

Hyperbolic paraboloid

Hyperbolic paraboloid

‘Hyperbolic Paraboloids
‘Sebastian Gmelin
’11.11.2009

Const a As Double = 2
Const b As Double = 2

Const range As Long = 8

Dim hp() As Point3d

Dim x As Double
Dim y As Double
Dim z As Double
   
Dim point As Point3d

Dim ele As element
Sub drawGrid()

Dim points() As Point3d

ReDim points(range * 2)

    For x = range * -1 To range
        For y = range * -1 To range
            points(y + range) = hp(x + range, y + range)
        Next y
        Set ele = CreateLineElement1(Nothing, points)
        ele.Color = 4
        ele.LineWeight = 0
        ActiveModelReference.AddElement ele
        ele.Redraw
    Next x
   
    For y = range * -1 To range
        For x = range * -1 To range
            points(x + range) = hp(x + range, y + range)
        Next x
        Set ele = CreateLineElement1(Nothing, points)
        ele.Color = 4
        ele.LineWeight = 0
        ActiveModelReference.AddElement ele
        ele.Redraw
    Next y

End Sub
Sub drawDiaGrid()

Dim points() As Point3d
Dim i As Long

    For i = 1 To range * 2
        x = range * -1
        ReDim points(i)
        For y = range – i To range
            points(y – range + i) = hp(x + y + i, y + range)
        Next y
        Set ele = CreateLineElement1(Nothing, points)
        ele.Color = 7
        ele.LineWeight = 1
        ActiveModelReference.AddElement ele
        ele.Redraw
       
        For y = i To 0 Step -1
            points(i – y) = hp(i – y, y)
        Next y
        Set ele = CreateLineElement1(Nothing, points)
        ele.Color = 7
        ele.LineWeight = 1
        ActiveModelReference.AddElement ele
        ele.Redraw
       
        If i < range * 2 Then
            y = range * -1
            ReDim points(range * 2 – i)
            For x = range * -1 + i To range
                points(x + range – i) = hp(x + range, y + range * 2 + x – i)
            Next x
            Set ele = CreateLineElement1(Nothing, points)
            ele.Color = 7
            ele.LineWeight = 1
            ActiveModelReference.AddElement ele
            ele.Redraw
           
            y = range
            ReDim points(range * 2 – i)
            For x = range * -1 + i To range
                points(x + range – i) = hp(x + range, y – x + i)
            Next x
            Set ele = CreateLineElement1(Nothing, points)
            ele.Color = 7
            ele.LineWeight = 1
            ActiveModelReference.AddElement ele
            ele.Redraw
        End If
    Next i
   
   
End Sub

Sub horizontalCut(cz As Double)

Dim points() As Point3d
Dim pointsNeg() As Point3d

ReDim points(range * 2)
ReDim pointsNeg(range * 2)
   
    If cz >= 0 Then
        For y = range * -1 To range
            x = Math.Sqr(a * a * (((y * y) / (b * b)) + 2 * cz))
            points(y + range) = Point3dFromXYZ(x, y, cz)
            pointsNeg(y + range) = Point3dFromXYZ(-1 * x, y, cz)
        Next y
    Else
        For x = range * -1 To range
            y = Math.Sqr(b * b * (((x * x) / (a * a)) – 2 * cz))
            points(x + range) = Point3dFromXYZ(x, y, cz)
            pointsNeg(x + range) = Point3dFromXYZ(x, -1 * y, cz)
        Next x
    End If

    Set ele = CreateLineElement1(Nothing, points)
    ele.Color = 5
    ele.LineWeight = 2
    ActiveModelReference.AddElement ele
    ele.Redraw
   
    Set ele = CreateLineElement1(Nothing, pointsNeg)
    ele.Color = 5
    ele.LineWeight = 2
    ActiveModelReference.AddElement ele
    ele.Redraw

End Sub

 

Sub Main()
   
    ReDim hp(range * 2, range * 2)
   
    For x = range * -1 To range
        For y = range * -1 To range
            z = (((x * x) / (a * a)) – ((y * y) / (b * b))) / 2
            point = Point3dFromXYZ(x, y, z)
            Set ele = CreateLineElement2(Nothing, point, point)
            ele.Color = 3
            ele.LineWeight = 10
            ActiveModelReference.AddElement ele
            ele.Redraw
            hp(x + range, y + range) = point
        Next y
    Next x
   
    drawGrid
    drawDiaGrid
    horizontalCut (range / 2)
    horizontalCut (range / -2)

End Sub

Share

Update to research plan

Monday, November 2nd, 2009

Added diagrams to the research plan page.
http://www.gmelin.li/PhD/research-plan/

Share

Bezier Curves, algorithm of de Casteljau

Sunday, November 1st, 2009

I started reading on Curve models and geometries, different types of Spline algorithms. This is a first attempt to generate Bezier curvers using the algorithm of de Casteljau. The script uses a selected polyline as the control polygon of the curve. The Bezier is created in several steps to illustrate the process. It runs under Microstation MVBA.

Bezier curve of degree 3

Bezier curve of degree 3

Bezier curve of degree 4

Bezier curve of degree 4

Bezier curve of degree 5

Bezier curve of degree 5

Bezier curve of degree 8

Bezier curve of degree 8

‘Linestring Bezier Approximation after Casteljau’s algorithm
‘Sebastian Gmelin
’24.10.2009

‘draw and select linestring and run script

Const cDiv As Long = 5 ‘subdivision at which guidlines are drawn and a copy of the drawing set is created
Const div As Long = 50 ‘subdivision between cDiv for bezier spline points (total subdivision = cDiv * div)

Dim cPoints() As Point3d
Dim bPoints() As Point3d
Dim cSize As Point3d
Dim cLine As Element
Private Sub ScanDesignFile()
‘scan active designfile for selected linestrings
    Dim oElEnum As ElementEnumerator
    Dim oElem As Element

    ‘get collection of selected elements
    Set oElEnum = ActiveModelReference.GetSelectedElements
    ActiveModelReference.UnselectAllElements
   
    ‘go through selection set
    While oElEnum.MoveNext
        Set oElem = oElEnum.Current
        ‘check if active element is linestring
        If (oElem.IsVertexList) Then
            Debug.Print (“Linestring”)
            cPoints = oElem.AsVertexList.GetVertices ‘save collection of vertices
            Set cLine = oElem
        End If
    Wend
End Sub

Private Sub getSize()
‘get linestring size for element copies offset

Dim i As Long
Dim min As Point3d
Dim max As Point3d
Dim lineString As Element

    min = Point3dFromXYZ(1E+17, 1E+17, 1E+17)
    max = Point3dFromXYZ(-1E+17, -1E+17, -1E+17)
   
    For i = 0 To UBound(cPoints)
        If cPoints(i).X > max.X Then max.X = cPoints(i).X
        If cPoints(i).Y > max.Y Then max.Y = cPoints(i).Y
        If cPoints(i).Z > max.Z Then max.Z = cPoints(i).Z
        If cPoints(i).X < min.X Then min.X = cPoints(i).X
        If cPoints(i).Y < min.Y Then min.Y = cPoints(i).Y
        If cPoints(i).Z < min.Z Then min.Z = cPoints(i).Z
    Next
   
    cSize = Point3dFromXYZ(0, max.Y – min.Y, 0)
    cSize.Y = cSize.Y * -1.2
   
End Sub

Sub Casteljau()

Dim i As Long
Dim k As Double
Dim p As Long
Dim points() As Point3d ‘last parent point collection
Dim dPoints() As Point3d ‘division points
Dim mPoints() As Point3d ‘points to draw
Dim lineString As Element
Dim color As Long
Dim count As Long
Dim moveCount As Long

ReDim bPoints(0)
bPoints(0) = cPoints(0)
count = 0
moveCount = 0

For k = 1 / (div * cDiv) To 1 Step 1 / (div * cDiv) ‘k is scale factor of divisions
    points = cPoints
    dPoints = cPoints
    color = 48
    ‘calculate divisions
    While UBound(dPoints) > 1
    ‘points is the collection of vertices of the parent linestring
    ‘dPoints is the child linestring at scalefactor k
        For i = 0 To UBound(points) – 1
            dPoints(i) = Point3dAdd(Point3dScale(Point3dSubtract(points(i + 1), points(i)), k), points(i))
        Next i
        ReDim Preserve dPoints(UBound(points) – 1)
        ‘draw guidlines
        If count = div Then
            mPoints = dPoints
            ‘calculate offset
            For p = 0 To UBound(mPoints)
                mPoints(p) = Point3dAdd(mPoints(p), Point3dScale(cSize, moveCount))
            Next
            ‘draw line
            Set lineString = CreateLineElement1(Nothing, mPoints)
            lineString.color = color
            lineString.LineStyle = ActiveDesignFile.LineStyles(3)
            ActiveModelReference.AddElement lineString
            lineString.Redraw
            color = color + 48
        End If
        points = dPoints
    Wend
    count = count + 1
    ‘draw progress of bezier curve
    If count > div Then
        count = 0
        ‘copy control polygon
        Set lineString = ActiveModelReference.CopyElement(cLine)
        Call lineString.Move(Point3dScale(cSize, moveCount))
        ActiveModelReference.AddElement lineString
        lineString.Redraw
        moveCount = moveCount + 1
        ‘calculate offset of beziers pline
        mPoints = bPoints
        For p = 0 To UBound(mPoints)
            mPoints(p) = Point3dAdd(mPoints(p), Point3dScale(cSize, moveCount – 1))
        Next
        ‘draw bezier progress
        Set lineString = CreateLineElement1(Nothing, mPoints)
        lineString.color = 3
        lineString.LineWeight = 2
        ActiveModelReference.AddElement lineString
        lineString.Redraw
    End If
    ‘add calculated division point to collection of bezier points
    ReDim Preserve bPoints(UBound(bPoints) + 1)
    bPoints(UBound(bPoints)) = Point3dAdd(Point3dScale(Point3dSubtract(points(1), points(0)), k), points(0))
Next k

‘draw final bezier curve
mPoints = bPoints
For p = 0 To UBound(mPoints)
    mPoints(p) = Point3dAdd(mPoints(p), Point3dScale(cSize, moveCount – 1))
Next
Set lineString = CreateLineElement1(Nothing, mPoints)
lineString.color = 3
lineString.LineWeight = 2
ActiveModelReference.AddElement lineString
lineString.Redraw

End Sub
Sub Main()

    If ActiveModelReference.AnyElementsSelected Then
        ScanDesignFile
        getSize
        Casteljau
    End If

End Sub

more Bezier curves

more Bezier curves

more Bezier curves

more Bezier curves

more Bezier curves

more Bezier curves

more Bezier curves

more Bezier curves

more Bezier curves

more Bezier curves

Share