kindred/create
1
1
Fork 0
Code Issues 37 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity

FEM: Fenics meshes, add support for xdmf format, some more improvements

Browse Source
This commit is contained in:
joha2
2017-05-26 21:41:13 +01:00
committed by Yorik van Havre
parent 62a823854f
commit 9edcbe2dd5
9 changed files with 674 additions and 326 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/Mod/Fem/App/CMakeLists.txt
View File

@@ -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

6
src/Mod/Fem/CMakeLists.txt
View File

@@ -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

4
src/Mod/Fem/Init.py
View File

@@ -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")

339
src/Mod/Fem/importFenicsMesh.py
View File

@@ -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']
}

72
src/Mod/Fem/importToolsFem.py
View File

@@ -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
'''

40
src/Mod/Fem/readFenicsXDMF.py Normal file
View File

@@ -0,0 +1,40 @@
# ***************************************************************************
# * *
# * Copyright (c) 2017 - Johannes Hartung <j.hartung@gmx.net> *
# * *
# * 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': {}
}

233
src/Mod/Fem/readFenicsXML.py Normal file
View File

@@ -0,0 +1,233 @@
# ***************************************************************************
# * *
# * Copyright (c) 2017 - Johannes Hartung <j.hartung@gmx.net> *
# * *
# * 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']
}

200
src/Mod/Fem/writeFenicsXDMF.py Normal file
View File

@@ -0,0 +1,200 @@
# ***************************************************************************
# * *
# * Copyright (c) 2017 - Johannes Hartung <j.hartung@gmx.net> *
# * *
# * 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('''<?xml version="1.0"?>\n<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>\n''')
fp.write(etree.tostring(root, pretty_print=True))
fp.close()

102
src/Mod/Fem/writeFenicsXML.py Normal file
View File

@@ -0,0 +1,102 @@
# ***************************************************************************
# * *
# * Copyright (c) 2017 - Johannes Hartung <j.hartung@gmx.net> *
# * *
# * 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()