FEM: fenics xml mesh format, add mesh import

src/Mod/Fem/Init.py

@@ -37,6 +37,7 @@
 
 FreeCAD.addExportType("FEM formats (*.unv *.med *.dat *.inp)", "Fem")
 FreeCAD.addImportType("CalculiX result (*.frd)", "importCcxFrdResults")
+FreeCAD.addImportType("Fenics mesh file (*.xml)", "importFenicsMesh")
 FreeCAD.addExportType("Fenics mesh file (*.xml)", "importFenicsMesh")
 FreeCAD.addImportType("Mesh from Calculix/Abaqus input file (*.inp)", "importInpMesh")
 FreeCAD.addImportType("Z88 mesh file (*.txt)", "importZ88Mesh")

src/Mod/Fem/importFenicsMesh.py

@@ -21,16 +21,23 @@
 # ***************************************************************************
 from __future__ import print_function
 
-__title__ = "FreeCAD Fenics mesh writer"
+__title__ = "FreeCAD Fenics mesh reader and writer"
 __author__ = "Johannes Hartung"
 __url__ = "http://www.freecadweb.org"
 
+
+# TODO: check for second order elements
+# TODO: export mesh functions (to be defined, cell functions, vertex functions, facet functions)
+
+
 ## @package importFenicsMesh
 #  \ingroup FEM
-#  \brief FreeCAD Fenics Mesh writer for FEM workbench
+#  \brief FreeCAD Fenics Mesh reader and writer for FEM workbench
 
 import FreeCAD
+import FemMeshTools
 import os
+import itertools
 from lxml import etree  # parsing xml files and exporting
 
 
@@ -57,8 +64,7 @@ def insert(filename, docname):
     except NameError:
         doc = FreeCAD.newDocument(docname)
     FreeCAD.ActiveDocument = doc
-    FreeCAD.Console.PrintError("Not yet implemented")
-    # import_fenics_mesh(filename)
+    import_fenics_mesh(filename)
 
 
 def export(objectslist, filename):
@@ -76,6 +82,8 @@ def export(objectslist, filename):
 
 ########## module specific methods ##########
 # Helper
+
+########## Export Section ###################
 def get_FemMeshObjectDimension(fem_mesh_obj):
     """ Count all entities in an abstract sense, to distinguish which dimension the mesh is
         (i.e. linemesh, facemesh, volumemesh)
@@ -126,8 +134,6 @@ def write_fenics_mesh(fem_mesh_obj, outputfile):
         For the export, we only have to use the highest dimensional entities and their
         vertices to be exported. (For second order elements, we have to delete the mid element nodes.)
     """
-    # TODO: check for second order elements
-    # TODO: export mesh functions (to be defined, cell functions, vertex functions, facet functions)
 
     FreeCAD_to_Fenics_dict = {
         "Triangle": "triangle",
@@ -187,3 +193,210 @@ def write_fenics_mesh(fem_mesh_obj, outputfile):
     fp = pyopen(outputfile, "w")
     fp.write(etree.tostring(root, pretty_print=True))
     fp.close()
+
+############ Import Section ############
+
+
+def import_fenics_mesh(filename, analysis=None):
+    '''insert a FreeCAD FEM Mesh object in the ActiveDocument
+    '''
+    mesh_data = read_fenics_mesh(filename)
+    mesh_name = os.path.basename(os.path.splitext(filename)[0])
+    femmesh = FemMeshTools.make_femmesh(mesh_data)
+    if femmesh:
+        mesh_object = FreeCAD.ActiveDocument.addObject('Fem::FemMeshObject', mesh_name)
+        mesh_object.FemMesh = femmesh
+
+
+def read_fenics_mesh(xmlfilename):
+
+    Fenics_to_FreeCAD_dict = {
+        "triangle": "tria3",
+        "tetrahedron": "tetra4",
+        "hexahedron": "hexa8",
+        "interval": "seg2",
+        "quadrilateral": "quad4",
+    }
+
+    def read_mesh_block(mesh_block):
+        '''
+            Reading mesh block from XML file.
+            The mesh block only contains cells and vertices.
+        '''
+        dim = int(mesh_block.get("dim"))
+        cell_type = mesh_block.get("celltype")
+
+        vertex_size = 0
+
+        print("Mesh dimension: %d" % (dim,))
+        print("Mesh cell type: %s" % (cell_type,))
+
+        cells_parts_dim = {'point': {0: 1},
+                           'interval': {0: 2, 1: 1},
+                           'triangle': {0: 3, 1: 3, 2: 1},
+                           'tetrahedron': {0: 4, 1: 6, 2: 4, 3: 1},
+                           'quadrilateral': {0: 4, 1: 4, 2: 1},
+                           'hexahedron': {0: 8, 1: 12, 2: 6, 3: 1}}
+
+        find_vertices = mesh_block.find("vertices")
+        find_cells = mesh_block.find("cells")
+
+        nodes_dict = {}
+        cell_dict = {}
+
+        if find_vertices is None:
+            print("No vertices found!")
+        else:
+            vertex_size = int(find_vertices.attrib.get("size"))
+            print("Reading %d vertices" % (vertex_size,))
+
+            for vertex in find_vertices:
+                ind = int(vertex.get("index"))
+
+                if vertex.tag.lower() == 'vertex':
+                    [node_x, node_y, node_z] = [float(vertex.get(coord, 0.)) for coord in ["x", "y", "z"]]
+
+                    nodes_dict[ind+1] = FreeCAD.Vector(node_x, node_y, node_z)
+                    # increase node index by one, since fenics starts at 0, FreeCAD at 1
+                    # print("%d %f %f %f" % (ind, node_x, node_y, node_z))
+                else:
+                    print("found strange vertex tag: %s" % (vertex.tag,))
+
+        if find_cells is None:
+            print("No cells found!")
+        else:
+            print("Reading %d cells" % (int(find_cells.attrib.get("size")),))
+            for cell in find_cells:
+                ind = int(cell.get("index"))
+
+                if cell.tag.lower() != cell_type.lower():
+                    print("Strange mismatch between cell type %s and cell tag %s" % (cell_type, cell.tag.lower()))
+                num_vertices = cells_parts_dim[cell_type][0]
+
+                vtupel = tuple([int(cell.get("v"+str(vnum)))+1 for vnum in range(num_vertices)])
+                # generate "v0", "v1", ... from dimension lookup table
+                # increase numbers by one to match FC numbering convention
+
+                cell_dict[ind+1] = vtupel
+
+                # valtupel = tuple([ind] + list(vtupel))
+                # print(("%d " + ("%d "*len(vtupel))) % valtupel)
+
+        return (nodes_dict, cell_dict, cell_type, dim)
+
+    def generate_lower_dimensional_structures(nodes, cell_dict, cell_type, dim):
+
+        def correct_volume_det(element_dict):
+            '''
+                Checks whether the cell elements
+                all have the same volume (<0?)
+                sign (is necessary to avoid negative
+                Jacobian errors).
+                Works only with tet4 and tri3 elements at the moment
+            '''
+            if dim == 3:
+                for (ind, tet) in element_dict['tetra4'].iteritems():
+                    v0 = nodes[tet[0]]
+                    v1 = nodes[tet[1]]
+                    v2 = nodes[tet[2]]
+                    v3 = nodes[tet[3]]
+                    a = v1 - v0
+                    b = v2 - v0
+                    c = v3 - v0
+                    if a.dot(b.cross(c)) > 0:
+                        element_dict['tetra4'][ind] = (tet[1], tet[0], tet[2], tet[3])
+            if dim == 2:
+                nz = FreeCAD.Vector(0., 0., 1.)
+                for (ind, tria) in element_dict['tria3'].iteritems():
+                    v0 = nodes[tria[0]]
+                    v1 = nodes[tria[1]]
+                    v2 = nodes[tria[2]]
+                    a = v1 - v0
+                    b = v2 - v0
+                    if nz.dot(a.cross(b)) < 0:
+                        element_dict['tria3'][ind] = (tria[1], tria[0], tria[2])
+
+        element_dict = {}
+        element_counter = {}
+
+        # TODO: remove upper level lookup
+        for (key, val) in Fenics_to_FreeCAD_dict.iteritems():
+            element_dict[val] = {}
+            element_counter[key] = 0  # count every distinct element and sub element type
+
+        def addtupletodict(di, tpl, counter):
+            sortedtpl = tuple(sorted(tpl))
+            if di.get(sortedtpl) is None:
+                di[sortedtpl] = counter
+                counter += 1
+            return counter
+
+        def invertdict(dic):
+            invdic = {}
+            for (key, it) in dic.iteritems():
+                invdic[it] = key
+            return invdic
+
+        num_vert_dict = {'interval': 2,
+                         'triangle': 3,
+                         'tetrahedron': 4,
+                         'hexahedron': 8,
+                         'quadrilateral': 4}
+        lower_dims_dict = {'interval': [],
+                           'triangle': ['interval'],
+                           'tetrahedron': ['triangle', 'interval'],
+                           'hexahedron': ['quadrilateral', 'interval'],
+                           'quadrilateral': ['interval']}
+
+        for (cell_index, cell) in cell_dict.iteritems():
+            cell_lower_dims = lower_dims_dict[cell_type]
+            element_counter[cell_type] += 1
+            element_dict[Fenics_to_FreeCAD_dict[cell_type]][cell] = element_counter[cell_type]
+            for ld in cell_lower_dims:
+                for vertextuple in itertools.combinations(cell, num_vert_dict[ld]):
+                    element_counter[ld] = addtupletodict(
+                                                        element_dict[Fenics_to_FreeCAD_dict[ld]],
+                                                        vertextuple,
+                                                        element_counter[ld])
+
+        length_counter = len(nodes)
+        for (key, val_dict) in element_dict.iteritems():
+            # to ensure distinct indices for FreeCAD
+            for (vkey, it) in val_dict.iteritems():
+                val_dict[vkey] = it + length_counter
+            length_counter += len(val_dict)
+            # inverse of the dict (dict[key] = val -> dict[val] = key)
+            element_dict[key] = invertdict(val_dict)
+
+        correct_volume_det(element_dict)
+
+        return element_dict
+
+    nodes = {}
+    element_dict = {}
+    # TODO: remove two times initialization
+    for val in Fenics_to_FreeCAD_dict.itervalues():
+        element_dict[val] = {}
+
+    tree = etree.parse(xmlfilename)
+    root = tree.getroot()
+
+    if root.tag.lower() != "dolfin":
+        print("Strange root tag, should be dolfin!")
+
+    find_mesh = root.find("mesh")
+    if find_mesh is not None:  # these are consistency checks of the XML structure
+        print("Mesh found")
+        (nodes, cells_dict, cell_type, dim) = read_mesh_block(find_mesh)
+        element_dict = generate_lower_dimensional_structures(nodes, cells_dict, cell_type, dim)
+    else:
+        print("No mesh found")
+
+    if root.find("data") is not None:
+        print("Internal mesh data found")
+
+    return {'Nodes': nodes,
+            'Hexa8Elem': {}, 'Penta6Elem': {}, 'Tetra4Elem': element_dict['tetra4'], 'Tetra10Elem': {},
+            'Penta15Elem': {}, 'Hexa20Elem': {}, 'Tria3Elem': element_dict['tria3'], 'Tria6Elem': {},
+            'Quad4Elem': element_dict['quad4'], 'Quad8Elem': {}, 'Seg2Elem': element_dict['seg2']
+            }