Creating a bespoke shell mesher in FreeCAD

About the development of the FEM module/workbench.

Moderator: bernd

aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Hey all,

I've made this macro that can create quad meshes as a ruled surface between curves.
https://youtu.be/GkX37Bhy_p4
Capture.PNG
Capture.PNG (142.11 KiB) Viewed 3457 times
Can anyone help me improve this? Add it to the FEM gui somehow?
Any kind of interactive mesh generation is my goal here.

Anyone have any equivalencing scripts they'd care to contribute?

Did I do anything wrong in this macro? What fixes it?

Code: Select all

# There is a way to create an arc between points.
import FreeCAD
import Part
from FreeCAD import Base
import Fem
from math import pi as pi

def main():
    # get selections
    N_selections = 0
    sel = []
    for o in Gui.Selection.getSelectionEx():
        for s in o.SubElementNames:
            for s in o.SubObjects:
                N_selections = N_selections + 1
                sel.append(s)
    # Check that there are the correct number of selections
    if N_selections != 2:
        print("Select Two Edges\n")
        return
    Curve_01 = sel[0]
    Curve_02 = sel[1]
    # Need to make Curve_01 and Curve_02
    # It should be a priority to get these from a gui
    # Populate, Create, and show Curve_1
    #X1 = Base.Vector(0, 0, 0)
    #Y1 = Base.Vector(1, 1, 0)
    #Z1 = Base.Vector(2, 0.5, 0)
    #arc = Part.Arc(X1, Y1, Z1)
    #Curve_01 = arc.toShape()
    #Part.show(Curve_01)
    Arc_Length_01 = Curve_01.Length

    # Populate, Create, and show Curve_2
    #X2 = Base.Vector(0, 0, 6)
    #Y2 = Base.Vector(1, 1, 6)
    #Z2 = Base.Vector(2, 0, 6)
    #arc = Part.Arc(X2, Y2, Z2)
    #Curve_02 = arc.toShape()
    #Part.show(Curve_02)       
    Arc_Length_02 = Curve_02.Length

    # Can get document name with FreeCAD.ActiveDocument.Name
    # Knowledge From Here
    # https://wiki.freecadweb.org/Topological_data_scripting

	# I can get a list of "N" points by calling 'valueAt' on "Arc_Length"
    N_elms_X = 20
    N_elms_Y = 20

    # get the nodes on curve_01
    Nodes_On_Curve_1 = get_nodes_from_curve(Curve_01, N_elms_X)

	# get the nodes on curve_02
    Nodes_On_Curve_2 = get_nodes_from_curve(Curve_02, N_elms_X)

    # Make nodes by tracing the streamlines of the surface
    nodes = make_nodes_of_ruled_mesh(Nodes_On_Curve_1, Nodes_On_Curve_2, N_elms_Y)
    
    # Now to make the elements
    E2N = make_elements_of_ruled_mesh(N_elms_X, N_elms_Y)
    
    # Now have E2N and Nodes array
    # Add all of the nodes to a container called "a"
    a = Fem.FemMesh()
    
    # add all of the nodes to container "a"
    for i in range(len(nodes)):
        a.addNode(nodes[i][0], nodes[i][1], nodes[i][2], nodes[i][3])
    # add all of the elements to container "a"
    for i in range(len(E2N)):
        a.addFace([E2N[i][1], E2N[i][2], E2N[i][3], E2N[i][4]],E2N[i][0])
    # Making it render correctly
    obj = App.ActiveDocument.addObject("Fem::FemMeshObject", "a")
    obj.FemMesh = a
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

def make_elements_of_ruled_mesh(N_elms_X, N_elms_Y):
    Nx = N_elms_X + 1
    Ny = N_elms_Y + 1
    EID = 0
    E2N = []
    for j in range(Ny):
        for i in range(Nx):
            NID = 1 + Nx*j + i
            # skip rolling over elements
            if i == (Nx - 1):
                continue
            if j == (Ny - 1):
                continue
            EID = EID + 1
            N1 = NID
            N2 = 1 + Nx*(j+1) + i
            N3 = N2 + 1
            N4 = N1 + 1
            # E2N = [EID, N1, N2, N3, N4]
            E2N.append([EID, N1, N2, N3, N4])
    return E2N

def get_nodes_from_curve(Curve_Handle, N_elms):
    Arc_Length = Curve_Handle.Length
    nodes_on_curve = []
    for i in range(N_elms+1):
        dist_along = (i/N_elms)
        # Is it this? This worked originally.
        #dist_along = (i/N_elms)*Arc_Length 
        This_Node = Curve_Handle.valueAt(dist_along)
        x = This_Node.x
        y = This_Node.y
        z = This_Node.z
        nodes_on_curve.append([x, y, z])
    return nodes_on_curve

def make_nodes_of_ruled_mesh(Nodes_1, Nodes_2, N_elms_Y):
    # Make nodes by tracing the streamlines of the surface
    nodes = []
    # first set of nodes will be from Nodes_On_Curve_1
    NID = 0
    for i in range(len(Nodes_1)):
        # Get coordinates of the node on Curve 1
        x1 = Nodes_1[i][0]
        y1 = Nodes_1[i][1]
        z1 = Nodes_1[i][2]
        # Get coordinates of the node on Curve 2
        x2 = Nodes_2[i][0]
        y2 = Nodes_2[i][1]
        z2 = Nodes_2[i][2]
        # current node will be a weighted average of coords
        for j in range(N_elms_Y+1):
            NID = NID + 1
            x = (x2-x1)*(j/N_elms_Y) + x1
            y = (y2-y1)*(j/N_elms_Y) + y1
            z = (z2-z1)*(j/N_elms_Y) + z1
            nodes.append([x,y,z,NID])
    return nodes

if __name__ == '__main__':
    main()
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

It just kind of seems to... idk explode sometimes?

Its behavior appears to be very inconsistent and I don't know why.
For example, when get_nodes_from_curve looks like this:

Code: Select all

def get_nodes_from_curve(Curve_Handle, N_elms):
    Arc_Length = Curve_Handle.Length
    nodes_on_curve = []
    for i in range(N_elms+1):
        #dist_along = (i/N_elms)
        # Is it this? This worked originally.
        dist_along = (i/N_elms)*Arc_Length 
        This_Node = Curve_Handle.valueAt(dist_along)
        x = This_Node.x
        y = This_Node.y
        z = This_Node.z
        print([x, y, z])
        nodes_on_curve.append([x, y, z])
    return nodes_on_curve
This linear example works.
Linear_Example_Perfectly.PNG
Linear_Example_Perfectly.PNG (93.75 KiB) Viewed 3442 times
However, when run in that state, it fails the original test between two bsplines
Failure.PNG
Failure.PNG (103.52 KiB) Viewed 3442 times
I feel like I'm doing something extremely stupid.







Code: Select all

# There is a way to create an arc between points.
import FreeCAD
import Part
from FreeCAD import Base
import Fem
from math import pi as pi

def main():
    # get selections
    N_selections = 0
    sel = []
    for o in Gui.Selection.getSelectionEx():
        for s in o.SubElementNames:
            for s in o.SubObjects:
                N_selections = N_selections + 1
                sel.append(s)
    # Check that there are the correct number of selections
    if N_selections != 2:
        print("Select Two Edges\n")
        return
    Curve_01 = sel[0]
    Curve_02 = sel[1]

    Arc_Length_01 = Curve_01.Length

    Arc_Length_02 = Curve_02.Length

    # Knowledge From Here
    # https://wiki.freecadweb.org/Topological_data_scripting

	# I can get a list of "N" points by calling 'valueAt' on "Arc_Length"
    N_elms_X = 20
    N_elms_Y = 20

    # get the nodes on curve_01
    Nodes_On_Curve_1 = get_nodes_from_curve(Curve_01, N_elms_X)
    
	# get the nodes on curve_02
    Nodes_On_Curve_2 = get_nodes_from_curve(Curve_02, N_elms_X)

    # Make nodes by tracing the streamlines of the surface
    nodes = make_nodes_of_ruled_mesh(Nodes_On_Curve_1, Nodes_On_Curve_2, N_elms_Y)

    # Now to make the elements
    E2N = make_elements_of_ruled_mesh(N_elms_X, N_elms_Y)

    # Now have E2N and Nodes array
    # Add all of the nodes to a container called "a"
    a = Fem.FemMesh()

    # add all of the nodes to container "a"
    for i in range(len(nodes)):
        a.addNode(nodes[i][0], nodes[i][1], nodes[i][2], nodes[i][3])
    # add all of the elements to container "a"
    for i in range(len(E2N)):
        a.addFace([E2N[i][1], E2N[i][2], E2N[i][3], E2N[i][4]],E2N[i][0])
    # Making it render correctly
    obj = App.ActiveDocument.addObject("Fem::FemMeshObject", "a")
    obj.FemMesh = a
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

def make_elements_of_ruled_mesh(N_elms_X, N_elms_Y):
    Nx = N_elms_X + 1
    Ny = N_elms_Y + 1
    EID = 0
    E2N = []
    for j in range(Ny):
        for i in range(Nx):
            NID = 1 + Nx*j + i
            # skip rolling over elements
            if i == (Nx - 1):
                continue
            if j == (Ny - 1):
                continue
            EID = EID + 1
            N1 = NID
            N2 = 1 + Nx*(j+1) + i
            N3 = N2 + 1
            N4 = N1 + 1
            # E2N = [EID, N1, N2, N3, N4]
            E2N.append([EID, N1, N2, N3, N4])
    return E2N

def get_nodes_from_curve(Curve_Handle, N_elms):
    Arc_Length = Curve_Handle.Length
    nodes_on_curve = []
    for i in range(N_elms+1):
        #dist_along = (i/N_elms)
        # Is it this? This worked originally.
        dist_along = (i/N_elms)*Arc_Length 
        This_Node = Curve_Handle.valueAt(dist_along)
        x = This_Node.x
        y = This_Node.y
        z = This_Node.z
        print([x, y, z])
        nodes_on_curve.append([x, y, z])
    return nodes_on_curve

def make_nodes_of_ruled_mesh(Nodes_1, Nodes_2, N_elms_Y):
    # Make nodes by tracing the streamlines of the surface
    nodes = []
    # first set of nodes will be from Nodes_On_Curve_1
    NID = 0
    for i in range(len(Nodes_1)):
        # Get coordinates of the node on Curve 1
        x1 = Nodes_1[i][0]
        y1 = Nodes_1[i][1]
        z1 = Nodes_1[i][2]
        # Get coordinates of the node on Curve 2
        x2 = Nodes_2[i][0]
        y2 = Nodes_2[i][1]
        z2 = Nodes_2[i][2]
        # current node will be a weighted average of coords
        for j in range(N_elms_Y+1):
            NID = NID + 1
            x = (x2-x1)*(j/N_elms_Y) + x1
            y = (y2-y1)*(j/N_elms_Y) + y1
            z = (z2-z1)*(j/N_elms_Y) + z1
            nodes.append([x,y,z,NID])
    return nodes

if __name__ == '__main__':
    main()
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

OKAY I fixed it by doing something stupid, and I have no clue why this ever happens but it does and this has worked so far.
This_works_somehow.PNG
This_works_somehow.PNG (206.88 KiB) Viewed 3426 times
Basically, I had to write get_nodes_from_curve like this

Code: Select all

def get_nodes_from_curve(Curve_Handle, N_elms):
    Arc_Length = Curve_Handle.Length
    nodes_on_curve = []
    for i in range(N_elms+1):
        #dist_along = (i/N_elms)
        # Is it this? This worked originally.
        # Maybe need to check if FirstParameter=0.0 and LastParameter=1.0?
        first_par = Curve_Handle.FirstParameter
        last_par = Curve_Handle.LastParameter
        if first_par == 0.0 and last_par == 1.0:
            dist_along = (i/N_elms)
        else:
            dist_along = (i/N_elms)*Arc_Length
        This_Node = Curve_Handle.valueAt(dist_along)
        x = This_Node.x
        y = This_Node.y
        z = This_Node.z
        nodes_on_curve.append([x, y, z])
    return nodes_on_curve
With a stupid if statement that checks if both FirstParameter and LastParameter are 0.0 and 1.0 respectively. If they are, it uses the first formualation

Code: Select all

dist_along = (i/N_elms)
and if it's not, it uses the second formulation

Code: Select all

dist_along = (i/N_elms)*Arc_Length
No clue why this works, but it appears to.

Code: Select all

# There is a way to create an arc between points.
import FreeCAD
import Part
from FreeCAD import Base
import Fem
from math import pi as pi

def main():
    # get selections
    N_selections = 0
    sel = []
    for o in Gui.Selection.getSelectionEx():
        for s in o.SubElementNames:
            for s in o.SubObjects:
                N_selections = N_selections + 1
                sel.append(s)
    # Check that there are the correct number of selections
    if N_selections != 2:
        print("Select Two Edges\n")
        return
    Curve_01 = sel[0]
    Curve_02 = sel[1]

    Arc_Length_01 = Curve_01.Length

    Arc_Length_02 = Curve_02.Length

    # Knowledge From Here
    # https://wiki.freecadweb.org/Topological_data_scripting

	# I can get a list of "N" points by calling 'valueAt' on "Arc_Length"
    N_elms_X = 10
    N_elms_Y = 10

    # get the nodes on curve_01
    Nodes_On_Curve_1 = get_nodes_from_curve(Curve_01, N_elms_X)
    #return
	# get the nodes on curve_02
    Nodes_On_Curve_2 = get_nodes_from_curve(Curve_02, N_elms_X)

    # Make nodes by tracing the streamlines of the surface
    nodes = make_nodes_of_ruled_mesh(Nodes_On_Curve_1, Nodes_On_Curve_2, N_elms_Y)

    # Now to make the elements
    E2N = make_elements_of_ruled_mesh(N_elms_X, N_elms_Y)

    # Now have E2N and Nodes array
    # Add all of the nodes to a container called "a"
    a = Fem.FemMesh()

    # add all of the nodes to container "a"
    for i in range(len(nodes)):
        a.addNode(nodes[i][0], nodes[i][1], nodes[i][2], nodes[i][3])
    # add all of the elements to container "a"
    for i in range(len(E2N)):
        a.addFace([E2N[i][1], E2N[i][2], E2N[i][3], E2N[i][4]],E2N[i][0])
    # Making it render correctly
    obj = App.ActiveDocument.addObject("Fem::FemMeshObject", "a")
    obj.FemMesh = a
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

def make_elements_of_ruled_mesh(N_elms_X, N_elms_Y):
    Nx = N_elms_X + 1
    Ny = N_elms_Y + 1
    EID = 0
    E2N = []
    for j in range(Ny):
        for i in range(Nx):
            NID = 1 + Nx*j + i
            # skip rolling over elements
            if i == (Nx - 1):
                continue
            if j == (Ny - 1):
                continue
            EID = EID + 1
            N1 = NID
            N2 = 1 + Nx*(j+1) + i
            N3 = N2 + 1
            N4 = N1 + 1
            # E2N = [EID, N1, N2, N3, N4]
            E2N.append([EID, N1, N2, N3, N4])
    return E2N

def get_nodes_from_curve(Curve_Handle, N_elms):
    Arc_Length = Curve_Handle.Length
    nodes_on_curve = []
    for i in range(N_elms+1):
        #dist_along = (i/N_elms)
        # Is it this? This worked originally.
        # Maybe need to check if FirstParameter=0.0 and LastParameter=1.0?
        first_par = Curve_Handle.FirstParameter
        last_par = Curve_Handle.LastParameter
        if first_par == 0.0 and last_par == 1.0:
            dist_along = (i/N_elms)
        else:
            dist_along = (i/N_elms)*Arc_Length
        This_Node = Curve_Handle.valueAt(dist_along)
        x = This_Node.x
        y = This_Node.y
        z = This_Node.z
        nodes_on_curve.append([x, y, z])
    return nodes_on_curve

def make_nodes_of_ruled_mesh(Nodes_1, Nodes_2, N_elms_Y):
    # Make nodes by tracing the streamlines of the surface
    nodes = []
    # first set of nodes will be from Nodes_On_Curve_1
    NID = 0
    for i in range(len(Nodes_1)):
        # Get coordinates of the node on Curve 1
        x1 = Nodes_1[i][0]
        y1 = Nodes_1[i][1]
        z1 = Nodes_1[i][2]
        # Get coordinates of the node on Curve 2
        x2 = Nodes_2[i][0]
        y2 = Nodes_2[i][1]
        z2 = Nodes_2[i][2]
        # current node will be a weighted average of coords
        for j in range(N_elms_Y+1):
            NID = NID + 1
            x = (x2-x1)*(j/N_elms_Y) + x1
            y = (y2-y1)*(j/N_elms_Y) + y1
            z = (z2-z1)*(j/N_elms_Y) + z1
            nodes.append([x,y,z,NID])
    return nodes

if __name__ == '__main__':
    main()
heda
Veteran
Posts: 1348
Joined: Sat Dec 12, 2015 5:49 pm

Re: Creating a bespoke shell mesher in FreeCAD

Post by heda »

cool.

one can add macros to meny/toolbar with Customize_Toolbars for example, or create your own workbench (search forum).

took a glance at the code and pythonized it a bit,
of course a matter of taste, but often makes it easier to read and takes less time to write :-)

the construction of the points of the curve ended up being a one-liner, but also showed that use of "discretize" can throw curved balls.
the first/last parameters are not limited to range [0,1] - suppose you have made that assumption, it is also not the same as curve.Length,
(just to check in the python console for some different examples).

as for the selection, it looks like when one selects for example a spline made with draft in the gui, sometimes the selection object has subobjects, i.e. the edge, and sometimes not - maybe that is what is behind the "unpredictable" behaviour. took the opportunity to expand the acceptance of what can be selected, so now it seemingly works when one selects edges of a shape as well.

btw, the code as it is now is only intended to work with equal number of nodes in x/y, right?

Code: Select all

# App = FreeCAD, Gui = FreeCADGui
import FreeCAD, Part, Fem


def main():
    # gather edges from selection
    sel_edges = []
    for obj in Gui.Selection.getSelectionEx():
        if obj.HasSubObjects:
            for sub in obj.SubObjects:
                if isinstance(sub, Part.Edge):
                    sel_edges.append(sub)
        else:
            obj_edges = obj.Object.Shape.Edges
            if len(obj_edges) == 1:
                sel_edges.append(obj_edges[0])
    # Check that there are the correct number of selections
    if len(sel_edges) != 2:
        print("Select Two Edges", len(sel_edges))
        return
    Curve_01, Curve_02 = sel_edges

    # Knowledge From Here
    # https://wiki.freecadweb.org/Topological_data_scripting

    # I can get a list of "N" points by calling 'valueAt' on "Arc_Length"
    N_elms_X = N_elms_Y = 10

    # get the nodes on curve_01
    Nodes_On_Curve_1 = get_nodes_from_curve(Curve_01, N_elms_X)

	# get the nodes on curve_02
    Nodes_On_Curve_2 = get_nodes_from_curve(Curve_02, N_elms_X)

    # Make nodes by tracing the streamlines of the surface
    nodes = make_nodes_of_ruled_mesh(Nodes_On_Curve_1, Nodes_On_Curve_2, N_elms_Y)

    # Now to make the elements
    E2N = make_elements_of_ruled_mesh(N_elms_X, N_elms_Y)

    # Now have E2N and Nodes array
    # Add all of the nodes to a container called "a"
    a = Fem.FemMesh()

    # add all of the nodes to container "a"
    for node in nodes:
        a.addNode(*node)
    # add all of the elements to container "a"
    for eid, elnodes in E2N.items():
        a.addFace(elnodes, eid)
    # Making it render correctly
    obj = App.ActiveDocument.addObject("Fem::FemMeshObject", "a")
    obj.FemMesh = a
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

def make_elements_of_ruled_mesh(N_elms_X, N_elms_Y):
    Nx, Ny = N_elms_X + 1, N_elms_Y + 1
    EID = 0
    E2N = dict()
    for j in range(Ny):
        for i in range(Nx):
            NID = 1 + Nx*j + i
            # skip rolling over elements
            if i == (Nx - 1) or j == (Ny - 1):
                continue
            EID += 1
            N1 = NID
            N2 = 1 + Nx*(j+1) + i
            N3 = N2 + 1
            N4 = N1 + 1
            # E2N = [EID, N1, N2, N3, N4]
            E2N[EID] = [N1, N2, N3, N4]
    return E2N

def get_nodes_from_curve(Curve_Handle, N_elms):
    return tuple((n.x, n.y, n.z) for n in Curve_Handle.discretize(N_elms+1))

def make_nodes_of_ruled_mesh(Nodes_1, Nodes_2, N_elms_Y):
    """Make nodes by tracing the streamlines of the surface"""
    nodes = []
    # first set of nodes will be from Nodes_On_Curve_1
    NID = 0
    for node1i, node2i in zip(Nodes_1, Nodes_2):
        # Get coordinates of the node for Curve 1 & 2
        x1, y1, z1 = node1i
        x2, y2, z2 = node2i
        # current node will be a weighted average of coords
        for j in range(N_elms_Y+1):
            NID += 1
            x = (x2 - x1) * (j / N_elms_Y) + x1
            y = (y2 - y1) * (j / N_elms_Y) + y1
            z = (z2 - z1) * (j / N_elms_Y) + z1
            nodes.append([x, y, z, NID])
    return nodes

if __name__ == '__main__':
    main()


not directly related, but when doing this I got a weird behaviour using discretize that very much looks like a bug to me.
maybe someone proficient in c-side of things can take a look at that, maybe it deserves a unit-test as well.

the following example gives the difference between using discretize directly on the edge, or on the curve of the edge.

Code: Select all

>>> import Part
>>> box = Part.makeBox(4, 4, 4)
>>> for edge in box.Edges:
...   print('-'*25)
...   print(edge.discretize(5))
...   print(edge.Curve.discretize(5))
... 
-------------------------
[Vector (0.0, 0.0, 0.0), Vector (0.0, 0.0, 1.0), Vector (0.0, 0.0, 2.0), Vector (0.0, 0.0, 3.0), Vector (0.0, 0.0, 4.0)]
[Vector (0.0, 0.0, -2e+100), Vector (0.0, 0.0, -1e+100), Vector (0.0, 0.0, 0.0), Vector (0.0, 0.0, 1e+100), Vector (0.0, 0.0, 2e+100)]
-------------------------
[Vector (0.0, 0.0, 4.0), Vector (0.0, 1.0, 4.0), Vector (0.0, 2.0, 4.0), Vector (0.0, 3.0, 4.0), Vector (0.0, 4.0, 4.0)]
[Vector (0.0, -2e+100, 4.0), Vector (0.0, -1e+100, 4.0), Vector (0.0, 0.0, 4.0), Vector (0.0, 1e+100, 4.0), Vector (0.0, 2e+100, 4.0)]
-------------------------
[Vector (0.0, 4.0, 0.0), Vector (0.0, 4.0, 1.0), Vector (0.0, 4.0, 2.0), Vector (0.0, 4.0, 3.0), Vector (0.0, 4.0, 4.0)]
[Vector (0.0, 4.0, -2e+100), Vector (0.0, 4.0, -1e+100), Vector (0.0, 4.0, 0.0), Vector (0.0, 4.0, 1e+100), Vector (0.0, 4.0, 2e+100)]
-------------------------
[Vector (0.0, 0.0, 0.0), Vector (0.0, 1.0, 0.0), Vector (0.0, 2.0, 0.0), Vector (0.0, 3.0, 0.0), Vector (0.0, 4.0, 0.0)]
[Vector (0.0, -2e+100, 0.0), Vector (0.0, -1e+100, 0.0), Vector (0.0, 0.0, 0.0), Vector (0.0, 1e+100, 0.0), Vector (0.0, 2e+100, 0.0)]
-------------------------
[Vector (4.0, 0.0, 0.0), Vector (4.0, 0.0, 1.0), Vector (4.0, 0.0, 2.0), Vector (4.0, 0.0, 3.0), Vector (4.0, 0.0, 4.0)]
[Vector (4.0, 0.0, -2e+100), Vector (4.0, 0.0, -1e+100), Vector (4.0, 0.0, 0.0), Vector (4.0, 0.0, 1e+100), Vector (4.0, 0.0, 2e+100)]
-------------------------
[Vector (4.0, 0.0, 4.0), Vector (4.0, 1.0, 4.0), Vector (4.0, 2.0, 4.0), Vector (4.0, 3.0, 4.0), Vector (4.0, 4.0, 4.0)]
[Vector (4.0, -2e+100, 4.0), Vector (4.0, -1e+100, 4.0), Vector (4.0, 0.0, 4.0), Vector (4.0, 1e+100, 4.0), Vector (4.0, 2e+100, 4.0)]
-------------------------
[Vector (4.0, 4.0, 0.0), Vector (4.0, 4.0, 1.0), Vector (4.0, 4.0, 2.0), Vector (4.0, 4.0, 3.0), Vector (4.0, 4.0, 4.0)]
[Vector (4.0, 4.0, -2e+100), Vector (4.0, 4.0, -1e+100), Vector (4.0, 4.0, 0.0), Vector (4.0, 4.0, 1e+100), Vector (4.0, 4.0, 2e+100)]
-------------------------
[Vector (4.0, 0.0, 0.0), Vector (4.0, 1.0, 0.0), Vector (4.0, 2.0, 0.0), Vector (4.0, 3.0, 0.0), Vector (4.0, 4.0, 0.0)]
[Vector (4.0, -2e+100, 0.0), Vector (4.0, -1e+100, 0.0), Vector (4.0, 0.0, 0.0), Vector (4.0, 1e+100, 0.0), Vector (4.0, 2e+100, 0.0)]
-------------------------
[Vector (0.0, 0.0, 0.0), Vector (1.0, 0.0, 0.0), Vector (2.0, 0.0, 0.0), Vector (3.0, 0.0, 0.0), Vector (4.0, 0.0, 0.0)]
[Vector (-2e+100, 0.0, 0.0), Vector (-1e+100, 0.0, 0.0), Vector (0.0, 0.0, 0.0), Vector (1e+100, 0.0, 0.0), Vector (2e+100, 0.0, 0.0)]
-------------------------
[Vector (0.0, 0.0, 4.0), Vector (1.0, 0.0, 4.0), Vector (2.0, 0.0, 4.0), Vector (3.0, 0.0, 4.0), Vector (4.0, 0.0, 4.0)]
[Vector (-2e+100, 0.0, 4.0), Vector (-1e+100, 0.0, 4.0), Vector (0.0, 0.0, 4.0), Vector (1e+100, 0.0, 4.0), Vector (2e+100, 0.0, 4.0)]
-------------------------
[Vector (0.0, 4.0, 0.0), Vector (1.0, 4.0, 0.0), Vector (2.0, 4.0, 0.0), Vector (3.0, 4.0, 0.0), Vector (4.0, 4.0, 0.0)]
[Vector (-2e+100, 4.0, 0.0), Vector (-1e+100, 4.0, 0.0), Vector (0.0, 4.0, 0.0), Vector (1e+100, 4.0, 0.0), Vector (2e+100, 4.0, 0.0)]
-------------------------
[Vector (0.0, 4.0, 4.0), Vector (1.0, 4.0, 4.0), Vector (2.0, 4.0, 4.0), Vector (3.0, 4.0, 4.0), Vector (4.0, 4.0, 4.0)]
[Vector (-2e+100, 4.0, 4.0), Vector (-1e+100, 4.0, 4.0), Vector (0.0, 4.0, 4.0), Vector (1e+100, 4.0, 4.0), Vector (2e+100, 4.0, 4.0)]
>>> 

aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Thank you so much heda!!

I love the way this looks. I apologize for my lack of style.

As may have been obvious, I come from a fortran/matlab/octave background. My main hobby nowadays is helping maintain Mystran, and this is hopefully the first step towards a mesher that can be useful for shell optimized solvers like Mystran and Nastran95. I look forward to getting to know python better.

It was actually intended to work with different numbers of nodes in x and y directions. I believe that in the original form I posted it in, it did work that way but I introduced a bug at some point that broke it. I didn't notice it because I was too busy investigating the way that discretizing the edges worked. I'll be fixing that first thing tonight.

Thank you so much for elaborating on that strange splitting behavior, and telling me that it may be a bug. Should I go about trying to report it somewhere?
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Oh dear my sincerest apologies but a lot of this fancy python stuff is flying straight over my head.

The worst part is I can tell you avoided doing anything particularly difficult to follow and it's still a bit foreign to me.

I'll absolutely keep this as a template but if you see my old ugly style stick around, know it's not because I dislike your work. I'm just not good enough at python (YET!) to understand much of this.

Thank you once more.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Here's proof that the algorithms should inherently permit any number of elements on either side.
Proof_The_Algorithm_Should_Be_Right.PNG
Proof_The_Algorithm_Should_Be_Right.PNG (54.77 KiB) Viewed 3268 times
And here's the spreadsheet I used to check it.
Algorithm_Exploration.xlsx
(18.07 KiB) Downloaded 59 times
Pardon it not being in OBS format I had it almost done before realizing.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Rejoice!
Asymetric_Meshes_Work_Once_More.PNG
Asymetric_Meshes_Work_Once_More.PNG (195.92 KiB) Viewed 3247 times
I found the bug.
Classic indexing problem. Basic math. I'm an idiot.

Code: Select all

# App = FreeCAD, Gui = FreeCADGui
import FreeCAD, Part, Fem


def main():
    # gather edges from selection
    sel_edges = []
    for obj in Gui.Selection.getSelectionEx():
        if obj.HasSubObjects:
            for sub in obj.SubObjects:
                if isinstance(sub, Part.Edge):
                    sel_edges.append(sub)
        else:
            obj_edges = obj.Object.Shape.Edges
            if len(obj_edges) == 1:
                sel_edges.append(obj_edges[0])
    # Check that there are the correct number of selections
    if len(sel_edges) != 2:
        print("Select Two Edges", len(sel_edges))
        return
    Curve_01, Curve_02 = sel_edges

    # Knowledge From Here
    # https://wiki.freecadweb.org/Topological_data_scripting

    # I can get a list of "N" points by calling 'valueAt' on "Arc_Length"
    N_elms_X = 2
    N_elms_Y = 1

    # get the nodes on curve_01
    Nodes_On_Curve_1 = get_nodes_from_curve(Curve_01, N_elms_X)

	# get the nodes on curve_02
    Nodes_On_Curve_2 = get_nodes_from_curve(Curve_02, N_elms_X)
    
    # Make nodes by tracing the streamlines of the surface
    nodes = make_nodes_of_ruled_mesh(Nodes_On_Curve_1, Nodes_On_Curve_2, N_elms_Y)
    
    # Now to make the elements
    E2N = make_elements_of_ruled_mesh(N_elms_X, N_elms_Y)
    
    # Now have E2N and Nodes array
    # Add all of the nodes to a container called "a"
    a = Fem.FemMesh()

    # add all of the nodes to container "a"
    for node in nodes:
        a.addNode(*node)
    # add all of the elements to container "a"
    for E in E2N:
        print(E)
        a.addFace([E[1],E[2],E[3],E[4]],E[0])
    #return
    
    # Making it render correctly
    obj = App.ActiveDocument.addObject("Fem::FemMeshObject", "a")
    obj.FemMesh = a
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

def make_elements_of_ruled_mesh(N_elms_X, N_elms_Y):
    Nx = N_elms_X + 1
    Ny = N_elms_Y + 1
    EID = 0
    E2N = []
    for j in range(Ny):
        for i in range(Nx):
            NID = 1 + Nx*j + i
            # skip rolling over elements
            if i == (Nx - 1) or j == (Ny - 1):
                continue
            EID += 1
            N1 = NID
            N2 = 1 + Nx*(j+1) + i
            N3 = N2 + 1
            N4 = N1 + 1
            E2N.append([EID, N1, N2, N3, N4])
    return E2N

def get_nodes_from_curve(Curve_Handle, N_elms):
    return tuple((n.x, n.y, n.z) for n in Curve_Handle.discretize(N_elms+1))

def make_nodes_of_ruled_mesh(Nodes_1, Nodes_2, N_elms_Y):
    # Make nodes by tracing the streamlines of the surface
    nodes = []
    # first set of nodes will be from Nodes_On_Curve_1
    NID = 0
    for i in range(N_elms_Y+1):
        for j in range(len(Nodes_1)):
            NID = NID + 1
            # Get coordinates of the node on Curve 1
            x1 = Nodes_1[j][0]
            y1 = Nodes_1[j][1]
            z1 = Nodes_1[j][2]
            # Get coordinates of the node on Curve 2
            x2 = Nodes_2[j][0]
            y2 = Nodes_2[j][1]
            z2 = Nodes_2[j][2]
            # get location of current node in ruled surf
            x = (x2-x1)*(i/N_elms_Y) + x1
            y = (y2-y1)*(i/N_elms_Y) + y1
            z = (z2-z1)*(i/N_elms_Y) + z1
            nodes.append([x,y,z,NID])
    return nodes

if __name__ == '__main__':
    main()
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Started creating metrics to calculate quantities needed to judge mesh quality. First metric, which will be used to determine if curves should be assumed in opposite directions is the warping constant. To quote the MSC Nastran Reference Manual for 2021


compute the warping of each element according to MSC Nastran definition
According to MSC Nastran 2021 Reference Manual
Warp Test: This test evaluates how far out of plane the four corner
grid points are by measuring the distance of each point from a "mean"
plane passing through the locations of the four points.
The corner points are alternately H units above and H units below this
mean plane. If the lengths of the diagonals of the element are denoted
by D1 and D2, the warping coefficient is obtained from the equation
WC = H/2*(D1+D2)
If this value exceeds the tolerance, an informational message is produced.

The goal of this will be to make it able to tell if a curve needs to be reversed, as I've been trying on and off for almost a year now to overcome the inability to "reverse" curves and am yet to succeed. If all goes well it should be kinda hard to make a criss-crossed mesh with this.

Time will tell. Cue dramatic music.

Code: Select all

# App = FreeCAD, Gui = FreeCADGui
import FreeCAD, Part, Fem
"""
Development Goals:
==================
2021.07.25 || - [X] Add back in support for asymmetric mesh seeds
202X.XX.XX || - [ ] Add element warping function. Store in E2Warp
202X.XX.XX || - [ ] Use E2Warp to contemplate switching order of a curve
202X.XX.XX || - [ ] 
"""

def main():
    ## gather edges from selection
    sel_edges = []
    for obj in Gui.Selection.getSelectionEx():
        if obj.HasSubObjects:
            for sub in obj.SubObjects:
                if isinstance(sub, Part.Edge):
                    sel_edges.append(sub)
        else:
            obj_edges = obj.Object.Shape.Edges
            if len(obj_edges) == 1:
                sel_edges.append(obj_edges[0])
    # Check that there are the correct number of selections
    if len(sel_edges) != 2:
        print("Select Two Edges", len(sel_edges))
        return
    Curve_01, Curve_02 = sel_edges

    N_elms_X = 5
    N_elms_Y = 3

    # get the nodes on curve_01
    Nodes_On_Curve_1 = get_nodes_from_curve(Curve_01, N_elms_X)

	# get the nodes on curve_02
    Nodes_On_Curve_2 = get_nodes_from_curve(Curve_02, N_elms_X)
    
    # Make nodes by tracing the streamlines of the surface
    nodes = make_nodes_of_ruled_mesh(Nodes_On_Curve_1, Nodes_On_Curve_2, N_elms_Y)
    
    # Now to make the elements
    E2N = make_elements_of_ruled_mesh(N_elms_X, N_elms_Y)
    
    
    E2Warp = get_E2Warp(E2N, nodes)
    #return
    
    
    # Now have E2N and Nodes array
    # Add all of the nodes to a container called "a"
    a = Fem.FemMesh()

    # add all of the nodes to container "a"
    for node in nodes:
        a.addNode(*node)
    # add all of the elements to container "a"
    for E in E2N:
        print(E)
        a.addFace([E[1],E[2],E[3],E[4]],E[0])
    
    # Making it render correctly
    obj = App.ActiveDocument.addObject("Fem::FemMeshObject", "a")
    obj.FemMesh = a
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

def get_E2Warp(E2N, nodes):
    """
    compute the warping of each element according to MSC Nastran definition
    According to MSC Nastran 2021 Reference Manual
    Warp Test: This test evaluates how far out of plane the four corner
    grid points are by measuring the distance of each point from a "mean"
    plane passing through the locations of the four points.
    The corner points are alternately H units above and H units below this
    mean plane. If the lengths of the diagonals of the element are denoted
    by D1 and D2, the warping coefficient is obtained from the equation
    WC = H/2*(D1+D2)
    If this value exceeds the tolerance, an informational message is produced.
    """
    # [EID, WarpingConstant]
    E2Warp = []
    # for each element, get the nodes
    for E in E2N:
        EID = E[0]
        NID1 = E[1]
        NID2 = E[2]
        NID3 = E[3]
        NID4 = E[4]
        # scrape out data and get coords of all nodes
        for n in nodes: 
            NID = n[3]
            if NID == NID1:
                n1 = n       # [x1, y1, z1, NID1]
            elif NID == NID2:
                n2 = n       # [x2, y2, z2, NID2]
            elif NID == NID3:
                n3 = n       # [x3, y3, z3, NID3]
            elif NID == NID4:
                n4 = n       # [x4, y4, z4, NID4]
        # need midpoint between n1 and n3
        m13x = (n1[0] + n3[0])/2
        m13y = (n1[1] + n3[1])/2
        m13z = (n1[2] + n3[2])/2
        # need midpoint between n2 and n4
        m24x = (n2[0] + n4[0])/2
        m24y = (n2[1] + n4[1])/2
        m24z = (n2[2] + n4[2])/2
        # "H" is half the distance between these points
        dmx = m24x - m13x
        dmy = m24y - m13y
        dmz = m24z - m13z
        H = 0.5 * (((dmx ** 2) + (dmy ** 2) + (dmz ** 2))**0.5)
        # computing D1
        dD1x = n3[0] - n1[0]
        dD1y = n3[1] - n1[1]
        dD1z = n3[2] - n1[2]
        D1 = ((dD1x ** 2)+(dD1y ** 2)+(dD1z ** 2)) ** 0.5
        # computing D2
        dD2x = n4[0] - n2[0]
        dD2y = n4[1] - n2[1]
        dD2z = n4[2] - n2[2]
        D2 = ((dD2x ** 2)+(dD2y ** 2)+(dD2z ** 2)) ** 0.5
        WC = H / 2 * (D1 + D2)
        # Append EID, WC to E2Warp
        E2Warp.append([EID, WC])
    return E2Warp
    
def make_elements_of_ruled_mesh(N_elms_X, N_elms_Y):
    Nx = N_elms_X + 1
    Ny = N_elms_Y + 1
    EID = 0
    E2N = []
    for j in range(Ny):
        for i in range(Nx):
            NID = 1 + Nx*j + i
            # skip rolling over elements
            if i == (Nx - 1) or j == (Ny - 1):
                continue
            EID += 1
            N1 = NID
            N2 = 1 + Nx*(j+1) + i
            N3 = N2 + 1
            N4 = N1 + 1
            E2N.append([EID, N1, N2, N3, N4])
    return E2N

def get_nodes_from_curve(Curve_Handle, N_elms):
    return tuple((n.x, n.y, n.z) for n in Curve_Handle.discretize(N_elms+1))

def make_nodes_of_ruled_mesh(Nodes_1, Nodes_2, N_elms_Y):
    # Make nodes by tracing the streamlines of the surface
    nodes = []
    # first set of nodes will be from Nodes_On_Curve_1
    NID = 0
    for i in range(N_elms_Y+1):
        for j in range(len(Nodes_1)):
            NID = NID + 1
            # Get coordinates of the node on Curve 1
            x1 = Nodes_1[j][0]
            y1 = Nodes_1[j][1]
            z1 = Nodes_1[j][2]
            # Get coordinates of the node on Curve 2
            x2 = Nodes_2[j][0]
            y2 = Nodes_2[j][1]
            z2 = Nodes_2[j][2]
            # get location of current node in ruled surf
            x = (x2-x1)*(i/N_elms_Y) + x1
            y = (y2-y1)*(i/N_elms_Y) + y1
            z = (z2-z1)*(i/N_elms_Y) + z1
            nodes.append([x,y,z,NID])
    return nodes

if __name__ == '__main__':
    main()
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Code now supports excessive warping prevention.

Uses a quirk of the mathematics to do this without re-calculating the Element to Node array (E2N).
It's a bafflingly good idea.

Code: Select all

# App = FreeCAD, Gui = FreeCADGui
import FreeCAD, Part, Fem
"""
Development Goals:
==================
2021.07.25 || - [X] Add back in support for asymmetric mesh seeds
2021.07.25 || - [X] Add element warping function. Store in E2Warp
202X.XX.XX || - [X] Use E2Warp to contemplate switching order of a curve
202X.XX.XX || - [ ] 
"""

def main():
    ## gather edges from selection
    sel_edges = []
    for obj in Gui.Selection.getSelectionEx():
        if obj.HasSubObjects:
            for sub in obj.SubObjects:
                if isinstance(sub, Part.Edge):
                    sel_edges.append(sub)
        else:
            obj_edges = obj.Object.Shape.Edges
            if len(obj_edges) == 1:
                sel_edges.append(obj_edges[0])
    # Check that there are the correct number of selections
    if len(sel_edges) != 2:
        print("Select Two Edges", len(sel_edges))
        return
    Curve_01, Curve_02 = sel_edges

    N_elms_X = 5
    N_elms_Y = 3

    # get the nodes on curve_01
    N_Curve_1 = get_nodes_from_curve(Curve_01, N_elms_X)

	# get the nodes on curve_02
    N_Curve_2 = get_nodes_from_curve(Curve_02, N_elms_X)
    
    # Make nodes by tracing the streamlines of the surface
    nodes = make_nodes_of_ruled_mesh(N_Curve_1, N_Curve_2, N_elms_Y)
    
    # Now to make the elements
    E2N = make_elements_of_ruled_mesh(N_elms_X, N_elms_Y)
    
    # computing E2Warp
    E2Warp = get_E2Warp(E2N, nodes)
    

    # assemble package to pass into correct_curve_order_for_warp
    #        1.) N_Curve_1
    #        2.) N_Curve_2
    #        3.) N_elms_Y
    #        4.) E2N
    #        5.) E2Warp
    #        6.) nodes
    warp_correction_pack = [N_Curve_1, N_Curve_2, N_elms_Y, E2N, E2Warp, nodes]
    nodes = correct_curve_order_for_warp(warp_correction_pack)

    # Now have E2N and Nodes array
    # Add all of the nodes to a container called "a"
    a = Fem.FemMesh()

    # add all of the nodes to container "a"
    for node in nodes:
        a.addNode(*node)
    # add all of the elements to container "a"
    for E in E2N:
        print(E)
        a.addFace([E[1],E[2],E[3],E[4]],E[0])
    
    # Making it render correctly
    obj = App.ActiveDocument.addObject("Fem::FemMeshObject", "a")
    obj.FemMesh = a
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

def correct_curve_order_for_warp(pack):
    # unpacking input
    N_Curve_1 = pack[0]
    N_Curve_2 = pack[1]
    N_elms_Y = pack[2]
    E2N = pack[3]
    E2Warp = pack[4]
    nodes = pack[5]

    # if the curve needs reversing, it will change nodes
    # Reverse one of the curves, in order to see if improves mesh
    N_Curve_1_rev = N_Curve_1[::-1] # ???? python sucks holy crap
    # how do people think this is better than octave or julia

    # Make nodes of reversed curve 1 mesh
    nodes_rev = make_nodes_of_ruled_mesh(N_Curve_1_rev, N_Curve_2, N_elms_Y)

    # Make elements of reversed curve 1 mesh
    # awesomely, the same E2N holds true! How cool is that?!
    # computing E2Warp_rev
    E2Warp_rev = get_E2Warp(E2N, nodes_rev)

    # need to compare rev_warp to regular warp
    # computing sum of all warp in the elements

    total_rev_warp = 0
    for E2W_rev in E2Warp_rev:
        total_rev_warp += E2W_rev[1]

    total_warp = 0
    for E2W in E2Warp:
        total_warp += E2W[1]

    # if there's more total warp in rev than regular, return nodes_rev
    if total_warp > total_rev_warp:
        return nodes_rev
    else:
        return nodes

def get_E2Warp(E2N, nodes):
    """
    compute the warping of each element according to MSC Nastran definition
    According to MSC Nastran 2021 Reference Manual
    Warp Test: This test evaluates how far out of plane the four corner
    grid points are by measuring the distance of each point from a "mean"
    plane passing through the locations of the four points.
    The corner points are alternately H units above and H units below this
    mean plane. If the lengths of the diagonals of the element are denoted
    by D1 and D2, the warping coefficient is obtained from the equation
    WC = H/2*(D1+D2)
    If this value exceeds the tolerance, an informational message is produced.
    """
    # [EID, WarpingConstant]
    E2Warp = []
    # for each element, get the nodes
    for E in E2N:
        EID = E[0]
        NID1 = E[1]
        NID2 = E[2]
        NID3 = E[3]
        NID4 = E[4]
        # scrape out data and get coords of all nodes
        for n in nodes: 
            NID = n[3]
            if NID == NID1:
                n1 = n       # [x1, y1, z1, NID1]
            elif NID == NID2:
                n2 = n       # [x2, y2, z2, NID2]
            elif NID == NID3:
                n3 = n       # [x3, y3, z3, NID3]
            elif NID == NID4:
                n4 = n       # [x4, y4, z4, NID4]
        # need midpoint between n1 and n3
        m13x = (n1[0] + n3[0]) / 2
        m13y = (n1[1] + n3[1]) / 2
        m13z = (n1[2] + n3[2]) / 2
        # need midpoint between n2 and n4
        m24x = (n2[0] + n4[0]) / 2
        m24y = (n2[1] + n4[1]) / 2
        m24z = (n2[2] + n4[2]) / 2
        # "H" is half the distance between these points
        dmx = m24x - m13x
        dmy = m24y - m13y
        dmz = m24z - m13z
        H = 0.5 * (((dmx ** 2) + (dmy ** 2) + (dmz ** 2))**0.5)
        # computing D1
        dD1x = n3[0] - n1[0]
        dD1y = n3[1] - n1[1]
        dD1z = n3[2] - n1[2]
        D1 = ((dD1x ** 2)+(dD1y ** 2)+(dD1z ** 2)) ** 0.5
        # computing D2
        dD2x = n4[0] - n2[0]
        dD2y = n4[1] - n2[1]
        dD2z = n4[2] - n2[2]
        D2 = ((dD2x ** 2)+(dD2y ** 2)+(dD2z ** 2)) ** 0.5
        WC = H / 2 * (D1 + D2)
        # Append EID, WC to E2Warp
        E2Warp.append([EID, WC])
    return E2Warp
    
def make_elements_of_ruled_mesh(N_elms_X, N_elms_Y):
    Nx = N_elms_X + 1
    Ny = N_elms_Y + 1
    EID = 0
    E2N = []
    for j in range(Ny):
        for i in range(Nx):
            NID = 1 + Nx*j + i
            # skip rolling over elements
            if i == (Nx - 1) or j == (Ny - 1):
                continue
            EID += 1
            N1 = NID
            N2 = 1 + Nx*(j+1) + i
            N3 = N2 + 1
            N4 = N1 + 1
            E2N.append([EID, N1, N2, N3, N4])
    return E2N

def get_nodes_from_curve(Curve_Handle, N_elms):
    return tuple((n.x, n.y, n.z) for n in Curve_Handle.discretize(N_elms+1))

def make_nodes_of_ruled_mesh(Nodes_1, Nodes_2, N_elms_Y):
    # Make nodes by tracing the streamlines of the surface
    nodes = []
    # first set of nodes will be from N_Curve_1
    NID = 0
    for i in range(N_elms_Y+1):
        for j in range(len(Nodes_1)):
            NID = NID + 1
            # Get coordinates of the node on Curve 1
            x1 = Nodes_1[j][0]
            y1 = Nodes_1[j][1]
            z1 = Nodes_1[j][2]
            # Get coordinates of the node on Curve 2
            x2 = Nodes_2[j][0]
            y2 = Nodes_2[j][1]
            z2 = Nodes_2[j][2]
            # get location of current node in ruled surf
            x = (x2-x1)*(i/N_elms_Y) + x1
            y = (y2-y1)*(i/N_elms_Y) + y1
            z = (z2-z1)*(i/N_elms_Y) + z1
            nodes.append([x,y,z,NID])
    return nodes

if __name__ == '__main__':
    main()
Post Reply