From 9edcbe2dd58983a4c1dea4c9b682f7c0a2d0d0a4 Mon Sep 17 00:00:00 2001 From: joha2 Date: Fri, 26 May 2017 21:41:13 +0100 Subject: [PATCH] FEM: Fenics meshes, add support for xdmf format, some more improvements --- src/Mod/Fem/App/CMakeLists.txt | 4 + src/Mod/Fem/CMakeLists.txt | 6 +- src/Mod/Fem/Init.py | 4 +- src/Mod/Fem/importFenicsMesh.py | 339 ++------------------------------ src/Mod/Fem/importToolsFem.py | 72 +++++++ src/Mod/Fem/readFenicsXDMF.py | 40 ++++ src/Mod/Fem/readFenicsXML.py | 233 ++++++++++++++++++++++ src/Mod/Fem/writeFenicsXDMF.py | 200 +++++++++++++++++++ src/Mod/Fem/writeFenicsXML.py | 102 ++++++++++ 9 files changed, 674 insertions(+), 326 deletions(-) create mode 100644 src/Mod/Fem/readFenicsXDMF.py create mode 100644 src/Mod/Fem/readFenicsXML.py create mode 100644 src/Mod/Fem/writeFenicsXDMF.py create mode 100644 src/Mod/Fem/writeFenicsXML.py diff --git a/src/Mod/Fem/App/CMakeLists.txt b/src/Mod/Fem/App/CMakeLists.txt index 92fe5126f8..32b30a0478 100644 --- a/src/Mod/Fem/App/CMakeLists.txt +++ b/src/Mod/Fem/App/CMakeLists.txt @@ -72,6 +72,10 @@ SET(FemScripts_SRCS convert2TetGen.py importCcxDatResults.py importCcxFrdResults.py + readFenicsXML.py + writeFenicsXML.py + readFenicsXDMF.py + writeFenicsXDMF.py importFenicsMesh.py importInpMesh.py importToolsFem.py diff --git a/src/Mod/Fem/CMakeLists.txt b/src/Mod/Fem/CMakeLists.txt index 320e019aad..3a0c6c5953 100755 --- a/src/Mod/Fem/CMakeLists.txt +++ b/src/Mod/Fem/CMakeLists.txt @@ -20,7 +20,11 @@ INSTALL( convert2TetGen.py importCcxDatResults.py importCcxFrdResults.py - importFenicsMesh.py + readFenicsXML.py + writeFenicsXML.py + readFenicsXDMF.py + writeFenicsXDMF.py + importFenicsMesh.py importInpMesh.py importToolsFem.py importVTKResults.py diff --git a/src/Mod/Fem/Init.py b/src/Mod/Fem/Init.py index c9d9b1efd0..c7f1c66d04 100644 --- a/src/Mod/Fem/Init.py +++ b/src/Mod/Fem/Init.py @@ -37,8 +37,8 @@ if("BUILD_FEM_VTK" in FreeCAD.__cmake__): 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("Fenics mesh file (*.xml *.xdmf)", "importFenicsMesh") +FreeCAD.addExportType("Fenics mesh file (*.xml *.xdmf)", "importFenicsMesh") FreeCAD.addImportType("Mesh from Calculix/Abaqus input file (*.inp)", "importInpMesh") FreeCAD.addImportType("Z88 mesh file (*.txt)", "importZ88Mesh") FreeCAD.addExportType("Z88 mesh file (*.txt)", "importZ88Mesh") diff --git a/src/Mod/Fem/importFenicsMesh.py b/src/Mod/Fem/importFenicsMesh.py index b8feb85e66..48878e3d0b 100644 --- a/src/Mod/Fem/importFenicsMesh.py +++ b/src/Mod/Fem/importFenicsMesh.py @@ -25,11 +25,6 @@ __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 reader and writer for FEM workbench @@ -37,8 +32,11 @@ __url__ = "http://www.freecadweb.org" import FreeCAD import importToolsFem import os -import itertools -from lxml import etree # parsing xml files and exporting + +import readFenicsXML +import writeFenicsXML +import readFenicsXDMF +import writeFenicsXDMF # Template copied from importZ88Mesh.py. Thanks Bernd! @@ -67,7 +65,7 @@ def insert(filename, docname): import_fenics_mesh(filename) -def export(objectslist, filename): +def export(objectslist, fileString): "called when freecad exports a file" if len(objectslist) != 1: FreeCAD.Console.PrintError("This exporter can only export one object.\n") @@ -77,329 +75,24 @@ def export(objectslist, filename): FreeCAD.Console.PrintError("No FEM mesh object selected.\n") return - write_fenics_mesh(obj, filename) + if fileString != "": + fileName, fileExtension = os.path.splitext(fileString) + if fileExtension.lower() == '.xml': + writeFenicsXML.write_fenics_mesh_xml(obj, fileString) + elif fileExtension.lower() == '.xdmf': + writeFenicsXDMF.write_fenics_mesh_xdmf(obj, fileString) - -########## 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) - """ - dim = None - - if fem_mesh_obj.FemMesh.Nodes != (): - dim = 0 - if fem_mesh_obj.FemMesh.Edges != (): - dim = 1 - if fem_mesh_obj.FemMesh.Faces != (): - dim = 2 - if fem_mesh_obj.FemMesh.Volumes != (): - dim = 3 - - return dim - - -def get_FemMeshObjectElementTypes(fem_mesh_obj, remove_zero_element_entries=True): - """ - Spit out all elements in the mesh with their appropriate dimension. - """ - FreeCAD_element_names = [ - "Node", "Edge", "Hexa", "Polygon", "Polyhedron", - "Prism", "Pyramid", "Quadrangle", "Tetra", "Triangle"] - FreeCAD_element_dims = [0, 1, 3, 2, 3, 3, 3, 2, 3, 2] - - elements_list_with_zero = [(eval("fem_mesh_obj.FemMesh." + s + "Count"), s, d) for (s, d) in zip(FreeCAD_element_names, FreeCAD_element_dims)] - # ugly but necessary - if remove_zero_element_entries: - elements_list = [(num, s, d) for (num, s, d) in elements_list_with_zero if num > 0] - else: - elements_list = elements_list_with_zero - - return elements_list - - -def get_MaxDimElementFromList(elem_list): - """ - Gets element with the maximal dimension in the mesh to determine cells. - """ - elem_list.sort(key=lambda (num, s, d): d) - return elem_list[-1] - - -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.) - """ - - FreeCAD_to_Fenics_dict = { - "Triangle": "triangle", - "Tetra": "tetrahedron", - "Hexa": "hexahedron", - "Edge": "interval", - "Node": "point", - "Quadrangle": "quadrilateral", - - "Polygon": "unknown", "Polyhedron": "unknown", - "Prism": "unknown", "Pyramid": "unknown", - } - - print("Converting " + fem_mesh_obj.Label + " to fenics XML File") - print("Dimension of mesh: %d" % (get_FemMeshObjectDimension(fem_mesh_obj),)) - - elements_in_mesh = get_FemMeshObjectElementTypes(fem_mesh_obj) - print("Elements appearing in mesh: %s" % (str(elements_in_mesh),)) - celltype_in_mesh = get_MaxDimElementFromList(elements_in_mesh) - (num_cells, cellname_fc, dim_cell) = celltype_in_mesh - cellname_fenics = FreeCAD_to_Fenics_dict[cellname_fc] - print("Celltype in mesh -> %s and its Fenics name: %s" % (str(celltype_in_mesh), cellname_fenics)) - - root = etree.Element("dolfin", dolfin="http://fenicsproject.org") - meshchild = etree.SubElement(root, "mesh", celltype=cellname_fenics, dim=str(dim_cell)) - vertices = etree.SubElement(meshchild, "vertices", size=str(fem_mesh_obj.FemMesh.NodeCount)) - - for (nodeind, fc_vec) in fem_mesh_obj.FemMesh.Nodes.iteritems(): # python2 - etree.SubElement( - vertices, "vertex", index=str(nodeind - 1), - # FC starts from 1, fenics starts from 0 to size-1 - x=str(fc_vec[0]), y=str(fc_vec[1]), z=str(fc_vec[2])) - - cells = etree.SubElement(meshchild, "cells", size=str(num_cells)) - if dim_cell == 3: - fc_cells = fem_mesh_obj.FemMesh.Volumes - elif dim_cell == 2: - fc_cells = fem_mesh_obj.FemMesh.Faces - elif dim_cell == 1: - fc_cells = fem_mesh_obj.FemMesh.Edges - else: - fc_cells = () - - for (fen_ind, fc_volume_ind) in enumerate(fc_cells): - # FC starts after all other entities, fenics start from 0 to size-1 - nodeindices = fem_mesh_obj.FemMesh.getElementNodes(fc_volume_ind) - - cell_args = {} - for (vi, ni) in enumerate(nodeindices): - cell_args["v" + str(vi)] = str(ni - 1) - # generate as many v entries in dict as nodes are listed in cell (works only for first order elements) - - etree.SubElement(cells, cellname_fenics, index=str(fen_ind), **cell_args) - - etree.SubElement(meshchild, "data") - - fp = pyopen(outputfile, "w") - fp.write(etree.tostring(root, pretty_print=True)) - fp.close() - -############ Import Section ############ + # write_fenics_mesh(obj, filename) def import_fenics_mesh(filename, analysis=None): '''insert a FreeCAD FEM Mesh object in the ActiveDocument ''' - mesh_data = read_fenics_mesh(filename) + mesh_data = readFenicsXML.read_fenics_mesh_xml(filename) + # xdmf not operational + mesh_name = os.path.basename(os.path.splitext(filename)[0]) femmesh = importToolsFem.make_femmesh(mesh_data) if femmesh: mesh_object = FreeCAD.ActiveDocument.addObject('Fem::FemMeshObject', mesh_name) mesh_object.FemMesh = femmesh - - -def read_fenics_mesh(xmlfilename): - ''' - Returns element dictionary to be evaluated by make_femmesh later - ''' - - 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'] - } diff --git a/src/Mod/Fem/importToolsFem.py b/src/Mod/Fem/importToolsFem.py index a85d0312be..d0b8cf7c5a 100644 --- a/src/Mod/Fem/importToolsFem.py +++ b/src/Mod/Fem/importToolsFem.py @@ -33,6 +33,78 @@ from math import pow, sqrt import numpy as np +def get_FemMeshObjectOrder(fem_mesh_obj): + """ + Gets element order. Element order counting based on number of nodes on + edges. Edge with 2 nodes -> linear elements, Edge with 3 nodes -> + quadratic elements, and so on. No edges in mesh -> not determined. + (Is this possible? Seems to be a very degenerate case.) + If there are edges with different number of nodes appearing, return + list of orders. + """ + presumable_order = None + + edges = fem_mesh_obj.FemMesh.Edges + + if edges != (): + edges_length_set = list({len(fem_mesh_obj.FemMesh.getElementNodes(e)) for e in edges}) + # only need set to eliminate double entries + + if len(edges_length_set) == 1: + presumable_order = edges_length_set[0] - 1 + else: + presumable_order = [el - 1 for el in edges_length_set] + else: + print("Found no edges in mesh: Element order determination does not work without them.") + print(presumable_order) + + return presumable_order + + +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) + """ + dim = None + + if fem_mesh_obj.FemMesh.Nodes != (): + dim = 0 + if fem_mesh_obj.FemMesh.Edges != (): + dim = 1 + if fem_mesh_obj.FemMesh.Faces != (): + dim = 2 + if fem_mesh_obj.FemMesh.Volumes != (): + dim = 3 + + return dim + + +def get_FemMeshObjectElementTypes(fem_mesh_obj, remove_zero_element_entries=True): + """ + Spit out all elements in the mesh with their appropriate dimension. + """ + FreeCAD_element_names_dims = { + "Node": 0, "Edge": 1, "Hexa": 3, "Polygon": 2, "Polyhedron": 3, + "Prism": 3, "Pyramid": 3, "Quadrangle": 2, "Tetra": 3, "Triangle": 2} + + elements_list_with_zero = [(eval("fem_mesh_obj.FemMesh." + s + "Count"), s, d) for (s, d) in FreeCAD_element_names_dims.iteritems()] + # ugly but necessary + if remove_zero_element_entries: + elements_list = [(num, s, d) for (num, s, d) in elements_list_with_zero if num > 0] + else: + elements_list = elements_list_with_zero + + return elements_list + + +def get_MaxDimElementFromList(elem_list): + """ + Gets element with the maximal dimension in the mesh to determine cells. + """ + elem_list.sort(key=lambda (num, s, d): d) + return elem_list[-1] + + def make_femmesh(mesh_data): ''' makes an FreeCAD FEM Mesh object from FEM Mesh data ''' diff --git a/src/Mod/Fem/readFenicsXDMF.py b/src/Mod/Fem/readFenicsXDMF.py new file mode 100644 index 0000000000..d4c35577d1 --- /dev/null +++ b/src/Mod/Fem/readFenicsXDMF.py @@ -0,0 +1,40 @@ +# *************************************************************************** +# * * +# * Copyright (c) 2017 - Johannes Hartung * +# * * +# * This program is free software; you can redistribute it and/or modify * +# * it under the terms of the GNU Lesser General Public License (LGPL) * +# * as published by the Free Software Foundation; either version 2 of * +# * the License, or (at your option) any later version. * +# * for detail see the LICENCE text file. * +# * * +# * This program is distributed in the hope that it will be useful, * +# * but WITHOUT ANY WARRANTY; without even the implied warranty of * +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +# * GNU Library General Public License for more details. * +# * * +# * You should have received a copy of the GNU Library General Public * +# * License along with this program; if not, write to the Free Software * +# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * +# * USA * +# * * +# *************************************************************************** +from __future__ import print_function + +__title__ = "FreeCAD Fenics XDMF mesh reader" +__author__ = "Johannes Hartung" +__url__ = "http://www.freecadweb.org" + +## @package importFenicsXDMF +# \ingroup FEM +# \brief FreeCAD Fenics Mesh XDMF reader for FEM workbench + +def read_fenics_mesh_xdmf(xdmffilename): + + print("Not operational, yet") + + return {'Nodes': {}, + 'Hexa8Elem': {}, 'Penta6Elem': {}, 'Tetra4Elem': {}, 'Tetra10Elem': {}, + 'Penta15Elem': {}, 'Hexa20Elem': {}, 'Tria3Elem': {}, 'Tria6Elem': {}, + 'Quad4Elem': {}, 'Quad8Elem': {}, 'Seg2Elem': {} + } diff --git a/src/Mod/Fem/readFenicsXML.py b/src/Mod/Fem/readFenicsXML.py new file mode 100644 index 0000000000..afd1ee5e6c --- /dev/null +++ b/src/Mod/Fem/readFenicsXML.py @@ -0,0 +1,233 @@ +# *************************************************************************** +# * * +# * Copyright (c) 2017 - Johannes Hartung * +# * * +# * This program is free software; you can redistribute it and/or modify * +# * it under the terms of the GNU Lesser General Public License (LGPL) * +# * as published by the Free Software Foundation; either version 2 of * +# * the License, or (at your option) any later version. * +# * for detail see the LICENCE text file. * +# * * +# * This program is distributed in the hope that it will be useful, * +# * but WITHOUT ANY WARRANTY; without even the implied warranty of * +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +# * GNU Library General Public License for more details. * +# * * +# * You should have received a copy of the GNU Library General Public * +# * License along with this program; if not, write to the Free Software * +# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * +# * USA * +# * * +# *************************************************************************** +from __future__ import print_function + +__title__ = "FreeCAD Fenics XML mesh reader" +__author__ = "Johannes Hartung" +__url__ = "http://www.freecadweb.org" + +## @package importFenicsXML +# \ingroup FEM +# \brief FreeCAD Fenics Mesh XML reader for FEM workbench + + +import FreeCAD +import importToolsFem +from lxml import etree # parsing xml files and exporting +import itertools + + +def read_fenics_mesh_xml(xmlfilename): + ''' + Returns element dictionary to be evaluated by make_femmesh later + ''' + + 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'] + } diff --git a/src/Mod/Fem/writeFenicsXDMF.py b/src/Mod/Fem/writeFenicsXDMF.py new file mode 100644 index 0000000000..e8bf37fab8 --- /dev/null +++ b/src/Mod/Fem/writeFenicsXDMF.py @@ -0,0 +1,200 @@ +# *************************************************************************** +# * * +# * Copyright (c) 2017 - Johannes Hartung * +# * * +# * This program is free software; you can redistribute it and/or modify * +# * it under the terms of the GNU Lesser General Public License (LGPL) * +# * as published by the Free Software Foundation; either version 2 of * +# * the License, or (at your option) any later version. * +# * for detail see the LICENCE text file. * +# * * +# * This program is distributed in the hope that it will be useful, * +# * but WITHOUT ANY WARRANTY; without even the implied warranty of * +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +# * GNU Library General Public License for more details. * +# * * +# * You should have received a copy of the GNU Library General Public * +# * License along with this program; if not, write to the Free Software * +# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * +# * USA * +# * * +# *************************************************************************** +from __future__ import print_function + +__title__ = "FreeCAD Fenics XDMF mesh writer" +__author__ = "Johannes Hartung" +__url__ = "http://www.freecadweb.org" + +## @package exportFenicsXDMF +# \ingroup FEM +# \brief FreeCAD Fenics Mesh XDMF writer for FEM workbench + +from importToolsFem import get_FemMeshObjectDimension, get_FemMeshObjectElementTypes, get_MaxDimElementFromList, get_FemMeshObjectOrder +from lxml import etree # parsing xml files and exporting +import numpy as np + +ENCODING_ASCII = 'ASCII' +ENCODING_HDF5 = 'HDF5' + +# TODO: export mesh functions (to be defined, cell functions, vertex functions, facet functions) +# TODO: integrate cell function +# TODO: check pyopen for other files + +# we need numpy functions to later access and process large data sets in a fast manner +# also the hd5 support better works together with numpy + + +def numpy_array_to_str(npa): + res = "" + dt = str(npa.dtype) + if 'int' in dt: + res = "\n".join([" ".join([("%d" % s) for s in a]) for a in npa.tolist()]) + elif 'float' in dt: + res = "\n".join([" ".join([("%3.6f" % s) for s in a]) for a in npa.tolist()]) + return res + + +def points_to_numpy(pts): + return np.array([[p.x, p.y, p.z] for p in pts]) + + +def tuples_to_numpy(tpls): + return np.array([list(t) for t in tpls]) + + +def write_fenics_mesh_points_xdmf(fem_mesh_obj, geometrynode, encoding=ENCODING_ASCII): + """ + Writes either into hdf5 file or into open mesh file + """ + + numnodes = fem_mesh_obj.FemMesh.NodeCount + + # dim = get_MaxDimElementFromList(get_FemMeshObjectElementTypes(fem_mesh_obj))[2] + # if dim == 2: + # geometrynode.set("GeometryType", "XY") + # elif dim == 3: + # geometrynode.set("GeometryType", "XYZ") + + geometrynode.set("GeometryType", "XYZ") + + # TODO: investigate: real two dimensional geometry. At the moment it is saved as + # flat 3d geometry. + + recalc_nodes_ind_dict = {} + + if encoding == ENCODING_ASCII: + dataitem = etree.SubElement(geometrynode, "DataItem", Dimensions="%d %d" % (numnodes, 3), Format="XML") + nodes = [] + for (ind, (key, node)) in enumerate(fem_mesh_obj.FemMesh.Nodes.iteritems()): + nodes.append(node) + recalc_nodes_ind_dict[key] = ind + + dataitem.text = numpy_array_to_str(points_to_numpy(nodes)) + elif encoding == ENCODING_HDF5: + pass + + return recalc_nodes_ind_dict + + +def write_fenics_mesh_volumes_xdmf(fem_mesh_obj, topologynode, rd, encoding=ENCODING_ASCII): + (num_cells, name_cell, dim_cell) = get_MaxDimElementFromList(get_FemMeshObjectElementTypes(fem_mesh_obj)) + element_order = get_FemMeshObjectOrder(fem_mesh_obj) + + FreeCAD_to_Fenics_XDMF_dict = { + ("Node", 1): ("polyvertex", 1), + ("Edge", 1): ("polyline", 2), + ("Edge", 2): ("edge_3", 3), + ("Triangle", 1): ("triangle", 3), + ("Triangle", 2): ("tri_6", 6), + ("Tetra", 1): ("tetrahedron", 4), + ("Tetra", 2): ("tet_10", 10) + } + + (topology_type, nodes_per_element) = FreeCAD_to_Fenics_XDMF_dict[(name_cell, element_order)] + + topologynode.set("TopologyType", topology_type) + topologynode.set("NumberOfElements", str(num_cells)) + topologynode.set("NodesPerElement", str(nodes_per_element)) + + if dim_cell == 3: + fc_cells = fem_mesh_obj.FemMesh.Volumes + elif dim_cell == 2: + fc_cells = fem_mesh_obj.FemMesh.Faces + elif dim_cell == 1: + fc_cells = fem_mesh_obj.FemMesh.Edges + elif dim_cell == 0: + fc_cells = fem_mesh_obj.FemMesh.Nodes + else: + fc_cells = [] + print("Dimension of mesh incompatible with export XDMF function: %d" % (dim_cell,)) + + nodeindices = [(rd[ind] for ind in fem_mesh_obj.FemMesh.getElementNodes(fc_volume_ind)) for (fen_ind, fc_volume_ind) in enumerate(fc_cells)] + # FC starts after all other entities, fenics start from 0 to size-1 + # write nodeindices into dict to access them later + + if encoding == ENCODING_ASCII: + dataitem = etree.SubElement(topologynode, "DataItem", NumberType="UInt", Dimensions="%d %d" % (num_cells, nodes_per_element), Format="XML") + dataitem.text = numpy_array_to_str(tuples_to_numpy(nodeindices)) + elif encoding == ENCODING_HDF5: + pass + + +def write_fenics_mesh_cellfunctions(fem_mesh_obj, mycellvalues, attributenode, encoding=ENCODING_ASCII): + attributenode.set("AttributeType", "Scalar") + attributenode.set("Center", "Cell") + attributenode.set("Name", "f") + + (num_cells, name_cell, dim_cell) = get_MaxDimElementFromList(get_FemMeshObjectElementTypes(fem_mesh_obj)) + + if encoding == ENCODING_ASCII: + dataitem = etree.SubElement(attributenode, "DataItem", Dimensions="%d %d" % (num_cells, 1), Format="XML") + dataitem.text = numpy_array_to_str(np.random.random((num_cells, 1))) + elif encoding == ENCODING_HDF5: + pass + + +def write_fenics_mesh_xdmf(fem_mesh_obj, outputfile, encoding=ENCODING_ASCII): + """ + For the export of xdmf. + """ + + FreeCAD_to_Fenics_dict = { + "Triangle": "triangle", + "Tetra": "tetrahedron", + "Hexa": "hexahedron", + "Edge": "interval", + "Node": "point", + "Quadrangle": "quadrilateral", + + "Polygon": "unknown", "Polyhedron": "unknown", + "Prism": "unknown", "Pyramid": "unknown", + } + + print("Converting " + fem_mesh_obj.Label + " to fenics XDMF File") + print("Dimension of mesh: %d" % (get_FemMeshObjectDimension(fem_mesh_obj),)) + + elements_in_mesh = get_FemMeshObjectElementTypes(fem_mesh_obj) + print("Elements appearing in mesh: %s" % (str(elements_in_mesh),)) + celltype_in_mesh = get_MaxDimElementFromList(elements_in_mesh) + (num_cells, cellname_fc, dim_cell) = celltype_in_mesh + cellname_fenics = FreeCAD_to_Fenics_dict[cellname_fc] + print("Celltype in mesh -> %s and its Fenics dolfin name: %s" % (str(celltype_in_mesh), cellname_fenics)) + + root = etree.Element("Xdmf", version="3.0") + domain = etree.SubElement(root, "Domain") + grid = etree.SubElement(domain, "Grid", Name="mesh", GridType="Uniform") + topology = etree.SubElement(grid, "Topology") + geometry = etree.SubElement(grid, "Geometry") + + # attribute = etree.SubElement(grid, "Attribute") # for cell functions + + recalc_dict = write_fenics_mesh_points_xdmf(fem_mesh_obj, geometry, encoding=encoding) + write_fenics_mesh_volumes_xdmf(fem_mesh_obj, topology, recalc_dict, encoding=encoding) + + # TODO: improve cell functions support + # write_fenics_mesh_cellfunctions(fem_mesh_obj, {}, attribute, encoding=encoding) + + fp = open(outputfile, "w") + fp.write('''\n\n''') + fp.write(etree.tostring(root, pretty_print=True)) + fp.close() diff --git a/src/Mod/Fem/writeFenicsXML.py b/src/Mod/Fem/writeFenicsXML.py new file mode 100644 index 0000000000..97f5f630f2 --- /dev/null +++ b/src/Mod/Fem/writeFenicsXML.py @@ -0,0 +1,102 @@ +# *************************************************************************** +# * * +# * Copyright (c) 2017 - Johannes Hartung * +# * * +# * This program is free software; you can redistribute it and/or modify * +# * it under the terms of the GNU Lesser General Public License (LGPL) * +# * as published by the Free Software Foundation; either version 2 of * +# * the License, or (at your option) any later version. * +# * for detail see the LICENCE text file. * +# * * +# * This program is distributed in the hope that it will be useful, * +# * but WITHOUT ANY WARRANTY; without even the implied warranty of * +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +# * GNU Library General Public License for more details. * +# * * +# * You should have received a copy of the GNU Library General Public * +# * License along with this program; if not, write to the Free Software * +# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * +# * USA * +# * * +# *************************************************************************** +from __future__ import print_function + +__title__ = "FreeCAD Fenics XML mesh writer" +__author__ = "Johannes Hartung" +__url__ = "http://www.freecadweb.org" + +## @package exportFenicsXML +# \ingroup FEM +# \brief FreeCAD Fenics Mesh XML writer for FEM workbench + + +from importToolsFem import get_FemMeshObjectDimension, get_FemMeshObjectElementTypes, get_MaxDimElementFromList +from lxml import etree # parsing xml files and exporting + + +def write_fenics_mesh_xml(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 (what to do? deny export or reduce element order?) + + FreeCAD_to_Fenics_dict = { + "Triangle": "triangle", + "Tetra": "tetrahedron", + "Hexa": "hexahedron", + "Edge": "interval", + "Node": "point", + "Quadrangle": "quadrilateral", + + "Polygon": "unknown", "Polyhedron": "unknown", + "Prism": "unknown", "Pyramid": "unknown", + } + + print("Converting " + fem_mesh_obj.Label + " to fenics XML File") + print("Dimension of mesh: %d" % (get_FemMeshObjectDimension(fem_mesh_obj),)) + + elements_in_mesh = get_FemMeshObjectElementTypes(fem_mesh_obj) + print("Elements appearing in mesh: %s" % (str(elements_in_mesh),)) + celltype_in_mesh = get_MaxDimElementFromList(elements_in_mesh) + (num_cells, cellname_fc, dim_cell) = celltype_in_mesh + cellname_fenics = FreeCAD_to_Fenics_dict[cellname_fc] + print("Celltype in mesh -> %s and its Fenics name: %s" % (str(celltype_in_mesh), cellname_fenics)) + + root = etree.Element("dolfin", dolfin="http://fenicsproject.org") + meshchild = etree.SubElement(root, "mesh", celltype=cellname_fenics, dim=str(dim_cell)) + vertices = etree.SubElement(meshchild, "vertices", size=str(fem_mesh_obj.FemMesh.NodeCount)) + + for (nodeind, fc_vec) in fem_mesh_obj.FemMesh.Nodes.iteritems(): # python2 + etree.SubElement( + vertices, "vertex", index=str(nodeind - 1), + # FC starts from 1, fenics starts from 0 to size-1 + x=str(fc_vec[0]), y=str(fc_vec[1]), z=str(fc_vec[2])) + + cells = etree.SubElement(meshchild, "cells", size=str(num_cells)) + if dim_cell == 3: + fc_cells = fem_mesh_obj.FemMesh.Volumes + elif dim_cell == 2: + fc_cells = fem_mesh_obj.FemMesh.Faces + elif dim_cell == 1: + fc_cells = fem_mesh_obj.FemMesh.Edges + else: + fc_cells = () + + for (fen_ind, fc_volume_ind) in enumerate(fc_cells): + # FC starts after all other entities, fenics start from 0 to size-1 + nodeindices = fem_mesh_obj.FemMesh.getElementNodes(fc_volume_ind) + + cell_args = {} + for (vi, ni) in enumerate(nodeindices): + cell_args["v" + str(vi)] = str(ni - 1) + # generate as many v entries in dict as nodes are listed in cell (works only for first order elements) + + etree.SubElement(cells, cellname_fenics, index=str(fen_ind), **cell_args) + + etree.SubElement(meshchild, "data") + + fp = open(outputfile, "w") # TODO: what about pyopen? + fp.write(etree.tostring(root, pretty_print=True)) + fp.close()