Creating a bespoke shell mesher in FreeCAD

About the development of the FEM module/workbench.

Moderator: bernd

User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Creating a bespoke shell mesher in FreeCAD

Post by bernd »

very cool stuff ...

two remarks in the regard of FreeCAD.

- At the beginning define "doc = App.ActiveDocument" than use doc, this way later it would even be possible to use the code not on the active document.

- In FEM everything is separated in Gui and App. FreeCAD can run headless withcout Gui, but only if the developer has taken this into account.

Every Gui code should be inside a

Code: Select all

if App.GuiUp:
    # Gui related imports and code
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Creating a bespoke shell mesher in FreeCAD

Post by bernd »

- we would need a icon for this cool stuff ...
User avatar
johnwang
Veteran
Posts: 1345
Joined: Sun Jan 27, 2019 12:41 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by johnwang »

Very nice.
hfc series CAE workbenches for FreeCAD (hfcNastran95, hfcMystran, hfcFrame3DD, hfcSU2 and more)
heda
Veteran
Posts: 1348
Joined: Sat Dec 12, 2015 5:49 pm

Re: Creating a bespoke shell mesher in FreeCAD

Post by heda »

"trying to report it" - dunno - have a feeling if it is a real bug and spotted by some core-dev, a report other than this is probably not needed.
first I thought it was some jibberish in there with overflow, but looking at the sequence of number it does not look like total jibberish, there is a pattern in there, at least on the printout.

on the "fortran shines through" - indeed, however no harm in that as such - any code that works is good code in my book, and the unpacking in python does feel alien if one is used to index into everything - but nowadays, once I got used to a bit of python, all this indexing is hurting my eyes :-), and dicts are pretty handy, also as a list-like-containers now that they are ordered.

small quick guide on unpacking, function calls in python - one just have to let go of the "one 2 one" notion ;)

Code: Select all


def fcn1(a, b, c):
    print(sum((a, b, c)))

fcn1(1,2,3) # prints 6


def fcn2(seq):
    print(sum(seq))

fcn2((1,2,3)) # prints 6, input could as well be [1,2,3]


def fcn3(args):
    print(sum(args))

fcn3((1,2,3)) # call is with sequence


def fcn4(*args):
    print(sum(args))

fcn4(1,2,3) # call is with 3 args
fcn4(*(1,2,3)) # still a call with 3 args, but now unpacked within the call



a, b = (1, 2)
c = (1, 2)
a, b = 1, 2

a = iter((1, 2))

next(a)

for i, j in ((1, 2), (3, 4)):
    print(i, j)
    
# prints
#1 2
#3 4

for i, j in zip((1, 2), (3, 4)):
    print(i, j)

#prints
#1 3
#2 4


yeah, it might not be the pretties reverse, using the [::-1] trick, but it is handy, for one it works on tuples as well as lists, and it makes a copy of the list, if one has a list one could use .reverse(), but that one is in-place, so if one wants the original around one have to copy it first.

on the equal nbr elements, it must have been me reading to code wrong to begin with - it looked like that to me (and I never really tested it), so then I did the zip, and with the zip it will limit it to equal length (unless one goes through some hoops). in any case, great that it works on different node counts in x and y.

impressive to get that warp algo working.

could not keep my fingers away from the code once again :-),
mainly because when posting yesterday I considered if I should have mentioned the thing with edge flip, but clearly you were on top of that anyway.

this time less intrusive, did an alternative way to check for edge-flip, suppose that way works as long as one only considers a 2 edged ruled surface, also put in a small dialogue there to select number of nodes as basic example.

Code: Select all


# App = FreeCAD, Gui = FreeCADGui
import FreeCAD, Part, Fem
from PySide import QtGui


Vector = App.Vector


"""
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 || - [ ] 
"""


class Form(QtGui.QDialog):
    """flip the true in __name__ == '__main__' to bypass the dialogue"""
    N_elms_X = 5
    N_elms_Y = 3
    
    def __init__(self):
        super(Form, self).__init__()
        self.setModal(True)
        self.makeUI()
        
    def makeUI(self):
        
        labelx = QtGui.QLabel('X elements')
        spinx = self.spinx = QtGui.QSpinBox()
        spinx.setValue(self.N_elms_X)
        spinx.setRange(2, 100)
        
        labely = QtGui.QLabel('Y elements')
        spiny = self.spiny = QtGui.QSpinBox()
        spiny.setValue(self.N_elms_Y)
        spiny.setRange(2, 100)

        btn = self.btn = QtGui.QPushButton('Mesh')
        btn.clicked.connect(self.make_mesh)
        
        layout = QtGui.QGridLayout()
        layout.addWidget(labelx, 0, 0)
        layout.addWidget(spinx, 1, 0)
        layout.addWidget(labely, 0, 1)
        layout.addWidget(spiny, 1, 1)
        layout.addWidget(btn, 2, 1)
        
        self.setLayout(layout)
        self.show()

    def get_values(self):
        return self.N_elms_X, self.N_elms_Y

    def make_mesh(self):
        self.N_elms_X = self.spinx.value()
        self.N_elms_Y = self.spiny.value()
        main(self.N_elms_X, self.N_elms_Y)
        self.close()





def main(N_elms_X, N_elms_Y):
    ## 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)
    
    # check if curve 2 should be reversed
    # use the fact that an edge has two vertexes
    vert2vec = lambda vert: Vector(*(getattr(vert, xyz) for xyz in 'XYZ'))
    e1p1 = vert2vec(Curve_01.Vertexes[0])
    e2p1, e2p2 = (vert2vec(vert) for vert in Curve_02.Vertexes)
    reverse = e1p1.distanceToPoint(e2p1) > e1p1.distanceToPoint(e2p2)
    if reverse:
        print('curve should be flipped')
        N_Curve_2 = N_Curve_2[::-1]
    
    
    # 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[0])
    
    # Making it render correctly
    doc = App.ActiveDocument
    obj = doc.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
    
    doc.recompute()

def correct_curve_order_for_warp(N_Curve_1, N_Curve_2, N_elms_Y,
                                 E2N, E2Warp, nodes):

    # 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:
        print('reversed edge')
        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 += 1
            # Get coordinates of the node on Curve 1
            x1, y1, z1 = Nodes_1[j]
            # Get coordinates of the node on Curve 2
            x2, y2, z2 = Nodes_2[j]
            # 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__':
    N_elms_X, N_elms_Y = 5, 3

    if True:
        form = Form()
    else:
        main(N_elms_X, N_elms_Y)
hope it is of some use.
keithsloan52
Veteran
Posts: 2756
Joined: Mon Feb 27, 2012 5:31 pm

Re: Creating a bespoke shell mesher in FreeCAD

Post by keithsloan52 »

Curious - Gmsh has options to create Quad Meshes for any shape. I don't think the FreeCAD options for using Gmsh as used by the Mesh Workbench and FEM offer/expose this option but Gmsh DOES have the facility. So my question why reinvent the wheel? Surely it would be better to add the facility to create Quad meshes to the existing FreeCAD facilities for invoking Gmsh.
User avatar
johnwang
Veteran
Posts: 1345
Joined: Sun Jan 27, 2019 12:41 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by johnwang »

keithsloan52 wrote: Mon Jul 26, 2021 8:41 pm Curious - Gmsh has options to create Quad Meshes for any shape. I don't think the FreeCAD options for using Gmsh as used by the Mesh Workbench and FEM offer/expose this option but Gmsh DOES have the facility. So my question why reinvent the wheel? Surely it would be better to add the facility to create Quad meshes to the existing FreeCAD facilities for invoking Gmsh.
Add ability to generate Quad mesh with GMSH is good.
Without using gmsh, adds more options.
hfc series CAE workbenches for FreeCAD (hfcNastran95, hfcMystran, hfcFrame3DD, hfcSU2 and more)
User avatar
johnwang
Veteran
Posts: 1345
Joined: Sun Jan 27, 2019 12:41 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by johnwang »

aerospaceweeb wrote: Sun Jul 25, 2021 1:53 am Add it to the FEM gui somehow?
Any kind of interactive mesh generation is my goal here.
Do you want to make a workbench?

FEM gui needs to touch cpp file and no guarantee to be accepted.
hfc series CAE workbenches for FreeCAD (hfcNastran95, hfcMystran, hfcFrame3DD, hfcSU2 and more)
User avatar
johnwang
Veteran
Posts: 1345
Joined: Sun Jan 27, 2019 12:41 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by johnwang »

Hope we can select an edge , then set it as a hinge. It should tell which two meshes it connected and which constrain to release.

Edit:
Maybe just after selected the edge, select the two meshes and save these info into a Hinge object.
Or select two meshes and find out the connection by code.
hinge2.jpg
hinge2.jpg (52.09 KiB) Viewed 2537 times
hinge.FCStd
(9.68 KiB) Downloaded 48 times
hfc series CAE workbenches for FreeCAD (hfcNastran95, hfcMystran, hfcFrame3DD, hfcSU2 and more)
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

keithsloan52 wrote: Mon Jul 26, 2021 8:41 pm So my question why reinvent the wheel? Surely it would be better to add the facility to create Quad meshes to the existing FreeCAD facilities for invoking Gmsh.
GMSH is no ordinary wheel. Gmsh is a pennyfarthing with afterburning turbojets that only reach peak efficiency above speeds of Mach 3.

I've spent over 150 hours of my life trying to learn how to make the GMSH gui do what I want it to, and I hated very second of it. It's frequently easier to write matlab code that writes a nastran bdf file than it is to make gmsh do what I want it to. In fifty years, when quantum computers let us auto-tetmesh everything and have it run in a reasonable amount of time, tools like gmsh will be the bees knees. For coarse meshes though, there's really, honestly nothing great in the realm of open source pre-processors, especially for extensibility.

FreeCADs python based scripting interface is a match made in heaven when it comes to interactive mesh generation and easily accessed geometric entity data-scraping.

Meshing algorithms are generally super easy, if your solver supports the right kinds of primitives, and if you're not foolish enough to trust stress results read directly from an assembly level finite element model.

IN FACT, I hypothesize that nearly all of the "good" open source pre-processors out there aren't as good as my crappy and juvenile workflow in FreeCAD. This is for the simple reason that the ability to easily manipulate geometry is more important than the ability to mesh obscure and poorly realized geometry.

The bleeding edge pre-processors nowadays like the newer Simcenter-NX, 3DX-Simulia, MSC Apex, and Hyperworks all sell themselves based on things like their geometry manipulation tools. Hell, even back in the day, MSC Patran marketed itself on its geometric kernel, which was WAY ahead of its time in many respects (compared to the competition). Patran is STILL the standard pre-processor against which all others are compared; usually subconsciously. Patran is a bottom up modeler. Make points. Make Curves. Make Surfaces. Make Solids. Mesh. Profit. Go bankrupt at the end of the year because Patran is 10 times more expensive than it has any business being.

So things like the Curves workbench make FreeCAD an abhorrently more powerful tool than basically anything else like gmsh.

You know how much of a pain in the ass it is to simply split a surface and remesh it? It's nuts!! And what's EVEN MORE infuriating about gmsh is the fact that there are people out there who are so good at using it that they could make me eat my words with a single youtube video!! Its features are so powerful I honestly wonder how in the heck I've screwed up the act of learning it for so long.

Buy hey, if me insulting gmsh gets some rockstar developer to make it build better into FreeCAD, I consider this a win-win.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

heda wrote: Mon Jul 26, 2021 6:14 pm "trying to report it" - dunno - have a feeling if it is a real bug and spotted by some core-dev, a report other than this is probably not needed.
first I thought it was some jibberish in there with overflow, but looking at the sequence of number it does not look like total jibberish, there is a pattern in there, at least on the printout.
...
I swear Heda you are a damn rockstar.

Reading code like this and seeing how subtley different it is to the garbage I pump out is more educational than anything I could possibly practice on my own. Thank you so much for your help.

Home from work. Let's get down to business.
Post Reply