fcFEM - FEA from start to finish

About the development of the FEM module/workbench.

Moderator: bernd

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

Re: fcFEM - FEA from start to finish

Postby bernd » Mon Feb 11, 2019 11:23 pm

FreeCAD has some methods isEqual, isPartner, isSame but they all only work on exact the same geomety, but not on equal geometry. Means two boxes could have equal geometry but they are not the same.

It is like two women could wear an equal dress (normally they never would, but they could :mrgreen: ), but never ever the same dress.

There is no method in part wb to check for equal geometry. I did some basic one for FEM in mesh tools. See https://github.com/FreeCAD/FreeCAD/blob ... 1422-L1454 It might not work on a compound because of CenterOfMass. Just try.
User avatar
bernd
Posts: 7321
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland

Re: fcFEM - FEA from start to finish

Postby bernd » Mon Feb 11, 2019 11:26 pm

User avatar
HarryvL
Posts: 941
Joined: Sat Jan 06, 2018 7:38 pm

Re: fcFEM - FEA from start to finish

Postby HarryvL » Tue Feb 12, 2019 5:13 am

bernd wrote:
Mon Feb 11, 2019 11:23 pm
FreeCAD has some methods isEqual, isPartner, isSame but they all only work on exact the same geomety, but not on equal geometry. Means two boxes could have equal geometry but they are not the same.

It is like two women could wear an equal dress (normally they never would, but they could :mrgreen: ), but never ever the same dress.

There is no method in part wb to check for equal geometry. I did some basic one for FEM in mesh tools. See https://github.com/FreeCAD/FreeCAD/blob ... 1422-L1454 It might not work on a compound because of CenterOfMass. Just try.
:D

I tried all three (IsPartner, IsSame and isEqual) and they only identify identical faces in my case. The search continues...
User avatar
bernd
Posts: 7321
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland

Re: fcFEM - FEA from start to finish

Postby bernd » Tue Feb 12, 2019 5:10 pm

Have you tried?
bernd wrote:
Mon Feb 11, 2019 11:23 pm
There is no method in part wb to check for equal geometry. I did some basic one for FEM in mesh tools. See https://github.com/FreeCAD/FreeCAD/blob ... 1422-L1454 It might not work on a compound because of CenterOfMass. Just try.
User avatar
HarryvL
Posts: 941
Joined: Sat Jan 06, 2018 7:38 pm

Re: fcFEM - FEA from start to finish

Postby HarryvL » Tue Feb 12, 2019 6:09 pm

bernd wrote:
Tue Feb 12, 2019 5:10 pm
Have you tried?
bernd wrote:
Mon Feb 11, 2019 11:23 pm
There is no method in part wb to check for equal geometry. I did some basic one for FEM in mesh tools. See https://github.com/FreeCAD/FreeCAD/blob ... 1422-L1454 It might not work on a compound because of CenterOfMass. Just try.
Not yet. I am first trying to figure out how FreeCAD interprets what GMSH returns in terms of element faces and element nodes when using Coherance=FALSE and feeding it with a compound of Boolean fragments.

In fact every contact has two sets of nodes and two sets of element faces, with one set of faces linked to one set of nodes. Both Boolean Fragment Faces at the contact refer to both sets. Quite confusing.
User avatar
HarryvL
Posts: 941
Joined: Sat Jan 06, 2018 7:38 pm

Re: fcFEM - FEA from start to finish

Postby HarryvL » Tue Feb 12, 2019 6:13 pm

Interface.FCStd
(19.57 KiB) Downloaded 3 times

Code: Select all

    elfaces=[]
    numelfaces=0
    for obj in App.ActiveDocument.Compound.Links:
        print("obj: {}\n".format(obj.Name))
        bfcount=0
        for bfface in obj.Shape.Faces:
            bfcount+=1
            print(">> object_face {}: {}\n".format(bfcount,bfface))
            for elface in mesh.getFacesByFace(bfface):
                numelfaces+=1
                elfaces.append(elface)
                print(">>>>object: {} element_face: {} with coordinates: \n>>>>>> node 1: {} {} \n>>>>>> node 2: {} {} \n>>>>>> node 3: {} {}\n".format(
                    obj.Name, elface,
                    mesh.getElementNodes(elface)[0], mesh.Nodes[mesh.getElementNodes(elface)[0]],
                    mesh.getElementNodes(elface)[1], mesh.Nodes[mesh.getElementNodes(elface)[1]],
                    mesh.getElementNodes(elface)[2], mesh.Nodes[mesh.getElementNodes(elface)[2]])
                      )
    print ("number of element faces: {}".format(numelfaces))
    print ("element faces array: {}".format(sorted(elfaces)))
User avatar
bernd
Posts: 7321
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland

Re: fcFEM - FEA from start to finish

Postby bernd » Tue Feb 12, 2019 9:04 pm

have you tried a normal Compound of Solids, not BooleanFragments?
User avatar
HarryvL
Posts: 941
Joined: Sat Jan 06, 2018 7:38 pm

Re: fcFEM - FEA from start to finish

Postby HarryvL » Tue Feb 12, 2019 9:25 pm

Yes, that is actually what the attached example is. However, in general I would need Boolean Fragments for more complex bodies than the solid primitives in part WB?
User avatar
HarryvL
Posts: 941
Joined: Sat Jan 06, 2018 7:38 pm

Re: fcFEM - FEA from start to finish

Postby HarryvL » Tue Feb 12, 2019 9:28 pm

CORRECTED TO RUN STAND ALONE:
HarryvL wrote:
Tue Feb 12, 2019 6:13 pm

Code: Select all

elfaces=[]
numelfaces=0
mesh = App.ActiveDocument.getObject("FEMMeshGmsh").FemMesh
for obj in App.ActiveDocument.Compound.Links:
    print("obj: {}\n".format(obj.Name))
    bfcount=0
    for bfface in obj.Shape.Faces:
        bfcount+=1
        print(">> object_face {}: {}\n".format(bfcount,bfface))
        for elface in mesh.getFacesByFace(bfface):
            numelfaces+=1
            elfaces.append(elface)
            print(">>>>object: {} element_face: {} with coordinates: \n>>>>>> node 1: {} {} \n>>>>>> node 2: {} {} \n>>>>>> node 3: {} {}\n".format(
                obj.Name, elface,
                mesh.getElementNodes(elface)[0], mesh.Nodes[mesh.getElementNodes(elface)[0]],
                mesh.getElementNodes(elface)[1], mesh.Nodes[mesh.getElementNodes(elface)[1]],
                mesh.getElementNodes(elface)[2], mesh.Nodes[mesh.getElementNodes(elface)[2]])
                  )
print ("number of element faces: {}".format(numelfaces))
print ("element faces array: {}".format(sorted(elfaces)))
User avatar
HarryvL
Posts: 941
Joined: Sat Jan 06, 2018 7:38 pm

Re: fcFEM - FEA from start to finish

Postby HarryvL » Wed Feb 13, 2019 9:24 pm

HarryvL wrote:
Sun Feb 10, 2019 7:27 am
HarryvL wrote:
Sat Feb 09, 2019 7:44 am
I am planning to profile the code and report out on where time is spent. Clearly, solving the matrix equations is the most computationally intensive, but here is where NumPy/SciPy comes the rescue. Including smart use of preconditioning with classical techniques, like Cuthill McKee band optimization (also available as a routine in SciPy btw). Assembly of the global stiffness matrix is not normally the most intensive, but involves a lot of loops. Here I Need to focus on maximizing use of masking, matrix multiplication and list comprehension instead of dum loops. Suggestions on cutting out for loops are welcome.
without any optimisation applied and naive application of Numpy I get the following results:



extract information from FreeCAD objects....................... 0.006 seconds
prepare finite element input................................... 0.528 seconds
calculate the global stiffness matrix and global load vector... 1.407 seconds
apply displacement boundary conditions......................... 0.003 seconds
solve the global siffness matrix equation...................... 18.344 seconds
calculate stresses from displacements.......................... 0.000 seconds
paste results in the FEM result object......................... 0.013 seconds


This is with 785 2nd order tetrahedral elements and 1399 nodes.
Simply upgrading to NumPy 1.16.3 reduces the time taken for NumPy.Solve() from 18.3 to 1.1 seconds