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

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

@Heda,

I'm sorry I'm a few days late getting responding to your messages in particular!
I'm so far behind reading any of the work you've done that I didn't feel it was right to comment before I had done so.

Embarrassingly, I STILL haven't done so! But I wanted to say that I will get back to you as soon as I can, and I cannot even begin to thank you enough for taking the tiny amount of work I've done and inflating it with 99% of your amazing work to make this piece of art.

Thank you.
You'll be the first name I mention when I finish adding the hex mesher.

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

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Also, @Heda,

With your permission, could I please instantiate your work in github?
I don't want to take credit for work I did not do, but I also am having a lot of trouble building on your work this way.

- Aero
heda
Veteran
Posts: 1348
Joined: Sat Dec 12, 2015 5:49 pm

Re: Creating a bespoke shell mesher in FreeCAD

Post by heda »

aerospaceweeb wrote: Thu Aug 05, 2021 3:54 am With your permission, could I please instantiate your work in github?
absolutely - that was the general idea - to get you kick-started :-), well, that and as a first intro for myself on wb's.

meanwhile I have realized that it should have been made from the "freecad.workbench_starterkit" template (making it pip installable),
ah well, maybe next time around,
in any case it would not have hurt that it at least was mentioned as alternative on Workbench_creation
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 »

heda wrote: Thu Aug 05, 2021 6:35 pm ... meanwhile I have realized that it should have been made from the "freecad.workbench_starterkit" template (making it pip installable),
ah well, maybe next time around, ...
That would be cool and it is much cleaner because none FreeCAD code is separated from FreeCAD related code.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Pardon my absence,

I'm finally on vacation and I've got a bit of time to work on this more.

I've decided that I can't try to solve multiple problems at the same time. These problems are now--

1.) GMSH, Salome, and FreeCAD FEM don't strongly store IDs on finite element entities; instead opting to use indices of the finite element data arrays themelves as the "IDs" until export.
This is in contrast to most commercial finite element oriented meshers and pre-processors that I've used (Ansys, Comsol, Hypermesh, Abaqus, FeMap, Patran, etc) but IS exactly like how I know that a ton of CFD preprocessors do it. I know for certain that those CFD softwares did it because adding an ID to every node would have increased the database sizes by at least 33%, and there was never really a strong reason to have IDs before data export.

2.) Because gmsh is a robsut, top down mesh generator, orphan mesh manipulation capabilities are sparse.
This means that a lot of the things that are high on my list of key features I was mising when trying to make a mesh of a wing-box assembly are things I don't know how to do. The most important feature, by an extremely large margin, is equivolencing nodes together.
The gmsh encapsulation in FreeCAD proper is rock solid as far as I can tell, and its specialty is making contiguous meshes of singular bodies, either surface or solid.
The issue I have, for my specific application coming from the world of aerospace engineering, is that contiguous surface meshes, equivolenced together aren't as common as they used to be for most FEMs. Nowadays most assemblies I see are "fully sprung" where you have each fastener represented as its own 1D element.
A lap joint, or just a basic bolted joint, is typically something like two shell meshes, where there's an RBE3 grabbing a few nodes on the shell mesh of one side. The central node of the RBE3is attached to a CBAR/CBUSH element which represents the fastener, and that fastener crosses over from one shell mesh to the other, where it attaches to another RBE3, which is then attached to that shell mesh.
This leads to the next issue--

3.) Lack of support for RBE2, RBE3, and RBAR esque entities
I'm not aware of any support in gmsh for RBE2 or RBE3 elements, which then leads me to think that there's no way to do this in FreeCAD. The good news is that they are remarkably easy elements to make, which have benefits on the runtimes of models when they're used. I've seen many coarse models where there are more RBE2 and RBE3 elements than regular elements, as they're frequently used to smear masses around using CONM2 elements.

4.) There's no way to export multiple mesh bodies
This is possibly not too much of an issue, as I imagine you could store multiple meshes together in a single mesh entitiy, but coupled with the next issue and issue #1, this becomes a massive issue.

5.) There's no way to control the IDs of a mesh entity on export
If you wanted to mesh two things, and then export them individually, you'd have both of their node ID and element IDs starting at 1; meaning they would conflict with eachother, preventing them from being run in a simulation or easily assembled in any other software.

I think I want to try solving the equivalencing problem.
I've gotten some boilerplate code to work with the KDTree library from scipy. It's very, very, VERY promising, and once I'm done trying to get some primitives to work, I'll have something to go off of once I'm back from my vacation and at my desktop computer.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Github page is up.

https://github.com/zchlrnr/FCMesher

Added a shell mesh thickener to help with an unrelated project. It just fit to add it here.

Gui support/integration is still in work, but standardizing the backend is what I'm focusing.

Open to suggestions for features to add, but here's my list in no particular order
1.) Equivolencing features
2.) The ability to export multiple mesh entities in the model at the same time, without the IDs overlapping with eachother.
3.) Transition all tools that make bdfs to use the build in tools in FreeCAD.
4.) Make a macro with a little popup window that can thicken a selected shell mesh.
5.) Add ability to visualize element normals. I don't know if that exists in free cad yet.
6.) Add ability to reverse element normals / make uniform the element normals of a shell mesh
7.) Add support for solid element thickener for tri elements
8.) Implement error handling for solid_mesh_by_thickened_shell_mesh
9.) Add element, node, and property ID offsetting subroutine. Related to #2.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Does anybody on here know how to get the nodes from an element?

I'm almost done with an equivalence routine, but I can't seem to figure out how to get the nodes that comprise an element, or as freeCAD calls them, "faces" if and only if they are tria or quad elements.
Screenshot from 2021-08-21 16-58-17.png
Screenshot from 2021-08-21 16-58-17.png (236.96 KiB) Viewed 2830 times
I'm sure it's in getAllDerivedFrom, but I can't figure out how to parse what that function creates.

EDIT_01: Found it. It's getElementNodes. All is well.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

All,

I've created a shell mesh thickener along normal vectors.

Code: Select all

# App = FreeCAD, Gui = FreeCADGui
from copy import copy as copy
import FreeCAD, Part, Fem
from PySide import QtGui
import os
import sys
import numpy as np
from scipy.spatial import KDTree

Vector = App.Vector

class Form(QtGui.QDialog): # {{{
    """ Set N_layers and total thickness
    """
    # savior --> https://doc.qt.io/qtforpython/tutorials/basictutorial/dialog.html
    # https://zetcode.com/gui/pysidetutorial/layoutmanagement/
    N_layers = 3
    thickness = 1
    
    def __init__(self): # {{{
        super(Form, self).__init__()
        self.setModal(True)
        self.makeUI()
        # }}}
        
    def makeUI(self): # {{{
        label_layers = QtGui.QLabel('N_layers')
        spin_layers = self.spin_layers = QtGui.QSpinBox()
        spin_layers.setValue(self.N_layers)
        spin_layers.setRange(1, 100)
       
        label_thickness = QtGui.QLabel('total thickness')
        thickness_field = self.thickness_field = QtGui.QLineEdit(str(self.thickness))

        btn = self.btn = QtGui.QPushButton('Thicken Shell Mesh')
        btn.clicked.connect(self.make_mesh)
        
        layout = QtGui.QGridLayout()
        layout.addWidget(label_layers, 0, 0)
        layout.addWidget(spin_layers, 0, 1)
        layout.addWidget(label_thickness, 1, 0)
        layout.addWidget(thickness_field, 1, 1)
        layout.addWidget(btn, 2, 1)
        
        self.setLayout(layout)
        self.show()
    # }}}

    def get_values(self): # {{{
        return self.N_layers, self.thickness
    # }}}

    def make_mesh(self): # {{{
        self.N_layers = self.spin_layers.value()
        self.thickness = float(self.thickness_field.text())
        main(self.N_layers, self.thickness)
        self.close()
    # }}}
# }}}
def get_E2N_nodes_and_E2T(mesh_objects_to_merge): # {{{
    """ takes in FemMesh objects, returns a combined E2N and nodes
    - [ ] Make it also return an E2T once other element types are supported
    """

    mesh_types_expected = []
    # check all mesh entities for pre-coded types # {{{
    for obj in mesh_objects_to_merge:
        # if there are any edge elements, raise error and abort
        if obj.EdgeCount != 0:
            s = "Edges not supported as of 2021.08.22"
            raise ValueError(s)
        elif obj.TriangleCount != 0:
            s = "Triangles not supported as of 2021.08.22"
            raise ValueError(s)
        elif obj.QuadrangleCount !=0:
            mesh_types_expected.append(15)
        elif obj.HexaCount != 0:
            s = "Hexas not supported as of 2021.08.22"
            raise ValueError(s)
        elif obj.TetraCount != 0:
            s = "Tetras not supported as of 2021.08.22"
            raise ValueError(s)
        elif obj.VolumeCount != 0:
            s = "Volumes not supported as of 2021.08.22"
            raise ValueError(s)
        elif obj.PyramidCount != 0:
            s = "Pyramids not supported as of 2021.08.22"
            raise ValueError(s)
        elif obj.PrismCount != 0:
            s = "Prisms not supported as of 2021.08.22"
            raise ValueError(s)
    # }}}
    mesh_types_expected = list(set(mesh_types_expected))

    # get [body ID, node ID old, node ID new]
    nid_transform = []
    body_ID = 1
    NID = 1
    for obj in mesh_objects_to_merge:
        for n in obj.Nodes:
            nid_transform.append([body_ID, n, NID])
            NID += 1
        body_ID += 1
    # get [body ID, EID old, EID new]
    eid_transform = []
    body_ID = 1
    EID = 1
    for obj in mesh_objects_to_merge:
        for e in obj.Faces:
            eid_transform.append([body_ID, e, EID])
            EID += 1
        body_ID += 1

    # create E2N for the ascending order
    current_body_ID = 0
    EID = 1
    E2N = {}
    for obj in mesh_objects_to_merge:
        current_body_ID += 1
        # get the part of eid transform that's in this body
        eid_transform_relevant = []
        for e in eid_transform:
            if e[0] == current_body_ID:
                eid_transform_relevant.append(e)
        # get the part of nid transform that's in this body
        nid_transform_relevant = []
        for n in nid_transform:
            if n[0] == current_body_ID:
                nid_transform_relevant.append(n)
        # making node ID dictionary [old --> new]
        n_rel = {}
        for n in nid_transform_relevant:
            n_rel[n[1]] = n[2]
        # for the relevant elements, make E2N
        for e in eid_transform_relevant:
            EID_old = e[1]
            EID_new = e[2]
            nodes_in_elm = list(obj.getElementNodes(EID_old))
            # replace nodes in elm with its new node IDs
            new_nodes_in_elm = []
            for n in nodes_in_elm:
                NID_new = n_rel[n]
                new_nodes_in_elm.append(NID_new)
            E2N[EID_new] = new_nodes_in_elm
    # create nodes
    nodes = {}
    for n in nid_transform:
        body_index = n[0] - 1
        old_NID = n[1]
        new_NID = n[2]
        obj = mesh_objects_to_merge[body_index]
        nodes[new_NID] = list(obj.Nodes[old_NID])

    # check that there was only one type of element, and it was 15 (CQUAD4)

    E2T = {}
    if len(mesh_types_expected) == 1 and mesh_types_expected[0] == 15:
        # Populate the E2T
        for e in E2N:
            E2T[e] = 15
    else:
        raise ValueError("Unknown elements passed into mesh equivalencer.")

    return [E2N, E2T, nodes]
#}}}
def get_E2NormVec(nodes, E2N): #{{{
    """ Compute the normal vector elements
    """
    E2NormVec = {}
    for EID in list(E2N.keys()):
        NodeIDs = E2N[EID]
        these_nodes = []
        for N in NodeIDs:
            # get the X coord of this node ID
            for i in list(nodes.keys()):
                if i == N:
                    x = nodes[i][0]
                    y = nodes[i][1]
                    z = nodes[i][2]
                    these_nodes.append([N, x, y, z])
        # compute every normal vector possible
        N = these_nodes
        # Get coordinates of all four points
        P1 = np.array([N[0][1], N[0][2], N[0][3]])
        P2 = np.array([N[1][1], N[1][2], N[1][3]])
        P3 = np.array([N[2][1], N[2][2], N[2][3]])
        P4 = np.array([N[3][1], N[3][2], N[3][3]])
        # Get vectors encircling the element
        V12 = P2 - P1
        V23 = P3 - P2
        V34 = P4 - P3
        V41 = P1 - P4
        # Get all the cross products
        NV1 = np.cross(V12, V23)
        NV2 = np.cross(V23, V34)
        NV3 = np.cross(V34, V41)
        NV4 = np.cross(V41, V12)
        # Compute average of all these vectors
        NV = NV1 + NV2 + NV3 + NV4
        # Normalize the magnitude
        mag = (NV[0]**2 + NV[1]**2 + NV[2]**2)**0.5
        NV = NV/mag
        E2NormVec[EID]= [NV[0], NV[1], NV[2]]
    return E2NormVec
    #}}}
def get_N2NormVec(E2NormVec, E2N, nodes): # {{{
    N2NormVec = {}
    N2E = get_N2E(E2N)
    
    for NID in nodes:
        elms_with_this_node = N2E[NID]
        X_comp = 0
        Y_comp = 0
        Z_comp = 0
        for EID in elms_with_this_node:
            X_comp += E2NormVec[EID][0]
            Y_comp += E2NormVec[EID][1]
            Z_comp += E2NormVec[EID][2]
        mag = ((X_comp**2) + (Y_comp**2) + (Z_comp**2))**0.5
        X = X_comp/mag
        Y = Y_comp/mag
        Z = Z_comp/mag
        N2NormVec[NID] = [X, Y, Z]
    return N2NormVec
# }}}
def get_N2E(E2N): # {{{
    """ Turns the E2N around, giving a dict of N2E
    """
    N2E = {}
    for EID, Element in E2N.items():
        # go through nodes in every element
        for NID in Element:
            # if the node's not in N2E yet, prepare for it to be
            if NID not in N2E:
                N2E[NID] = []
            # store the EID with that node ID we're on
            N2E[NID].append(EID)
    return N2E
# }}}
def get_new_nodes(N_layers, thickness, nodes, N2NormVec): # {{{
    # create nodes translated to the correct positions as new_nodes
    new_nodes = {}
    node_ID_offset = 0
    for layer in range(N_layers+1):
        # get vector magnitude
        mag = layer * (thickness / N_layers)
        for node in nodes.keys():
            x = nodes[node][0] + mag * N2NormVec[node][0]
            y = nodes[node][1] + mag * N2NormVec[node][1]
            z = nodes[node][2] + mag * N2NormVec[node][2]
            new_nodes[node + node_ID_offset] = [x, y, z]
        node_ID_offset += len(nodes.keys())
    return new_nodes
# }}}
def get_new_E2N(E2N, nodes, N_layers): # {{{
    new_E2N = {}
    element_ID_offset = 0
    for layer in range(N_layers):
        for element in E2N.keys():
            N1 = E2N[element][0] + layer * len(nodes.keys())
            N2 = E2N[element][1] + layer * len(nodes.keys())
            N3 = E2N[element][2] + layer * len(nodes.keys())
            N4 = E2N[element][3] + layer * len(nodes.keys())
            N5 = E2N[element][0] + (layer + 1) * len(nodes.keys())
            N6 = E2N[element][1] + (layer + 1) * len(nodes.keys())
            N7 = E2N[element][2] + (layer + 1) * len(nodes.keys())
            N8 = E2N[element][3] + (layer + 1) * len(nodes.keys())
            new_EID = element + element_ID_offset
            new_E2N[new_EID] = [N1, N2, N3, N4, N5, N6, N7, N8]
        element_ID_offset += len(E2N.keys())
    return new_E2N
# }}}
def main(N_layers, thickness): # {{{
    # gather FemMeshObject instances from selections
    mesh_objects_to_merge = [] 
    # gather edges from selection
    selected_edges = []
    for obj in Gui.Selection.getSelectionEx():
        if obj.TypeName == "Fem::FemMeshObject":
            mesh_objects_to_merge.append(obj.Object.FemMesh)
        elif obj.TypeName == "Fem::FemMeshObjectPython":
            mesh_objects_to_merge.append(obj.Object.FemMesh)
        elif obj.HasSubObjects:
            for sub in obj.SubObjects:
                if isinstance(sub, Part.Edge):
                    selected_edges.append(sub)
        else:
            obj_edges = obj.Object.Shape.Edges
            if len(obj_edges) == 1:
                selected_edges.append(obj_edges[0])
            
    # if there are no mesh objects selected, raise an error
    if len(mesh_objects_to_merge) == 0:
        raise ValueError("No mesh entities selected.")

    # determine the mode of the shell thickening
    # if there's no selected edges, set mode = 0
    if len(selected_edges) == 0:
        mode = 0
    # if there's one item in selected_edges, set mode = 1
    elif len(selected_edges) == 1:
        mode = 1

    # Throw errors on thicken mode
    if mode == 1:
        raise ValueError("As of 2021.08.29, no support for thicken mode 1")

    [E2N, E2T, nodes] = get_E2N_nodes_and_E2T(mesh_objects_to_merge)

    # Construct element ID to normal vector data structure
    E2NormVec = get_E2NormVec(nodes, E2N)

    # Construct node ID to normal vector data structure
    N2NormVec = get_N2NormVec(E2NormVec, E2N, nodes)

    # create nodes translated to the correct positions as new_nodes
    new_nodes =  get_new_nodes(N_layers, thickness, nodes, N2NormVec)

    # create elements calling out the new_nodes
    new_E2N = get_new_E2N(E2N, nodes, N_layers)

    # Now have new_E2N and new_nodes

    # create new FemMesh container called thickened
    thickened = Fem.FemMesh()

    # add all of the nodes to new container
    for node in new_nodes.keys():
        x = new_nodes[node][0]
        y = new_nodes[node][1]
        z = new_nodes[node][2]
        thickened.addNode(x, y, z, node)

    # create the elements of this container
    for element in new_E2N.keys():
        thickened.addVolume([*new_E2N[element]], element)

    # set graphical object to render correctly
    doc = App.ActiveDocument
    obj = doc.addObject("Fem::FemMeshObject", "thickened")
    obj.FemMesh = thickened 
    obj.Placement.Base = FreeCAD.Vector(0, 0, 0)
    obj.ViewObject.DisplayMode = "Faces, Wireframe & Nodes"
    obj.ViewObject.BackfaceCulling = False

    # render object
    doc.recompute()
# }}}

if __name__ == '__main__':
    form = Form()
Here's a clip of it working.


https://youtu.be/utuIerIpaog


The single most important part of this is the fact that it doesn't need geometry to work.

It is very common to have a legacy shell mesh as a .bdf file, and then have to take it into a solids only multiphysics solver, which necessitates the act of "thickening" it, and then subdividing it later.

My next course of action will be to add support for HEX8 elements into the equivalence tool that I've made.

I still have no way to easily apply properties to any of these meshes, as they don't have geometry, but a workaround for that is in work.

The leading proposal on another thread I can't recall was to make geometry for each element, but I haven't had much luck with that. It basically rendered the feature tree unusable as it filled with as many body objects as there were elements in the model, and needing geometry in the first place still doesn't address the fact that the majority of enterprise FEM work needs to use orphan meshes, as analysis frequently needs to leapfrog past design bottlenecks in order to justify said design changes in the first place, but I digress.

I think the next orders of business could also be to add a "reverse direction" checkbox, as well as a "extrude about mid surface" checkbox. Not sure how "in demand" those would be.

Complementary tools to this would include automatic node set generation for these bricks, as the nodes on the three new topological faces of the thickened solid would be very easy to reference using a closed form algebraic solution.

What do any of you all think? How can this be made better?
heda
Veteran
Posts: 1348
Joined: Sat Dec 12, 2015 5:49 pm

Re: Creating a bespoke shell mesher in FreeCAD

Post by heda »

looks like the fortran is beginning to fade away a bit :-)

as a person having done a bit of node-picking in a distant past, also i felt hampered by the fact that one could not directly operate on the mesh, but also have to say that it is really nice to just select on geometry when that is around.

hopefully someone will come around to open up pandoras box to get node-picking up and running, but it probably would be quite a bit of code, since one also should make convenience functions like selecting all nodes on a face, or on crest edge, etc to make the whole thing somewhat bearable to work with.
some of that functionality is already partly there in mesh wb, have you explored what could be done by just making it a mesh rather than a geometry?
sadly I do not think one can select/highlight on a mesh-wb-mesh, but it might be that it is not much work for a core-dev to make that possible and expose it to python.

meanwhile you should be able to use https://github.com/mwganson/MeshRemodel to create geometry out of selected parts of a mesh, not sure if you first have to make it a mesh-wb-mesh - but that you will figure out quickly.

hang in there, at some point you will have enough momentum so that people will join in and finetune it and start to fill in gaps.
aerospaceweeb
Posts: 118
Joined: Fri Apr 09, 2021 3:26 am

Re: Creating a bespoke shell mesher in FreeCAD

Post by aerospaceweeb »

Thank you Heda!

I could not have done it without you.

On the topic of node picking, I actually believe that everything I need is almost already here! The simple fact of the matter is that an excellent container for node lists exists, which maintains data connecting the label of the originating mesh entity to the node ID list, and that's the "Node Set" primitive, which has so far proven extremely easy to write code to and from.

The only thing I have left that I could possibly want in terms of easy primitive containers would be an element set container, akin to the node set container. The ability to manually add elements and nodes to those containers is an excellent replacement for 80-90% of typical femap or patran esque workflows.

With regards to the mesh workbench, I haven't really looked at it too hard. The concept of the mesh workbench is kind of at odds with how a lot of FEM work ought be done, in my opinion. It's the same reason that I didn't immediately go to Blender to try and make a FEM meshing tool there, despite the fact that it has a TON of advantages over pretty much anything else for specific things, like quad meshing shell structures. The mesh workbench does have a TON of code in it though, so perhaps I'll sit down with a beer or two and peruse it this next weekend.

In geometry related news, I'm actually having a TON of sucess using geometry to pick FEM data WITHOUT explicit association! I thought at first that FreeCAD had ways of doing it, but to my understanding it doesn't! It doesn't, despite the fact that FreeCAD has so many excellent geometry querying tools, that a script to "find all nodes on this surface" was surprisingly easy to write. Its performance isn't good yet though, but I did write it before I knew about KDTREE so perhaps I need to revisit it again.

In the past day, I've tried making a small tool to split a GMSH mesh into different mesh entities by element type, as I think that is a particularly annoying thing that gmsh does by default.

I also think that I can make a pretty excellent hex mesher of solids that have six sides, but that's only pseudocode at the moment. I have no desire to do any of the heavy lifting that GMSH has already done. I think that the best I should aim for, before other folks jump in to assist (why would they, FEM is such a boring development category ;) ) is the ability to do attachments between multiple different FemMeshes. The ID schema are a REAL pain in the ass, and I HATE HATE HATE that Salome defaulted to the destrction of any good node and element numbering schemes I can devise, but before I can get a workaround for that ready, I need the ability to just get an assembly out of the darn software and into my textual file of choice.

I hope it doesn't alienate people when they see it's nastran formats before calculix formats.
Post Reply