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

FEM: code formating, max line length < 100, fem in out modules

Browse Source
This commit is contained in:
Bernd Hahnebach
2019-05-27 07:59:02 +02:00
parent dbc0c489cd
commit 53a2e1e59d
11 changed files with 545 additions and 183 deletions
Download Patch File Download Diff File Expand all files Collapse all files

77
src/Mod/Fem/feminout/convert2TetGen.py
View File

@@ -46,12 +46,19 @@ def exportMeshToTetGenPoly(meshToExport, filePath, beVerbose=1):
f = open(filePath, 'w')
f.write("# This file was generated from FreeCAD geometry\n")
f.write("# Part 1 - node list\n")
'''
f.write("%(TotalNumOfPoints)i %(NumOfDimensions)i %(NumOfProperties)i %(BoundaryMarkerExists)i\n" % {
'TotalNumOfPoints': len(allVertices),
'NumOfDimensions': 3,
'NumOfProperties': 0,
'BoundaryMarkerExists': 0
})
'''
f.write(
"TotalNumOfPoints: {}, NumOfDimensions; {}, "
"NumOfProperties: {}, BoundaryMarkerExists: {}\n"
.format(len(allVertices), 3, 0, 0)
)
for PointIndex in range(len(allVertices)):
f.write("%(PointIndex)5i %(x) e %(y) e %(z) e\n" % {
'PointIndex': PointIndex,
@@ -93,10 +100,13 @@ def exportMeshToTetGenPoly(meshToExport, filePath, beVerbose=1):
EdgeKeys = EdgeFacets.keys()
# disconnectedEdges = len(EdgeKeys)
if beVerbose == 1:
FreeCAD.Console.PrintMessage('\nBoundaryMarker:' + repr(BoundaryMarker) + ' ' + repr(len(EdgeFacets)))
FreeCAD.Console.PrintMessage(
'\nBoundaryMarker:' + repr(BoundaryMarker) + ' ' + repr(len(EdgeFacets))
)
searchForPair = 1
# Main loop: first search for all complementary facets, then fill one branch and repeat while edges are available
# Main loop: first search for all complementary facets
# then fill one branch and repeat while edges are available
while len(EdgeFacets) > 0:
removeEdge = 0
for EdgeIndex in EdgeKeys:
@@ -142,7 +152,9 @@ def exportMeshToTetGenPoly(meshToExport, filePath, beVerbose=1):
searchForPair = 0
# End of main loop
if beVerbose == 1:
FreeCAD.Console.PrintMessage('\nNew BoundaryMarker:' + repr(BoundaryMarker) + ' ' + repr(len(EdgeFacets)))
FreeCAD.Console.PrintMessage(
'\nNew BoundaryMarker:' + repr(BoundaryMarker) + ' ' + repr(len(EdgeFacets))
)
## Part 2 - write all facets to *.poly file
f.write("# Part 2 - facet list\n")
@@ -194,7 +206,8 @@ def createMesh():
# Init objects
if beVerbose == 1:
FreeCAD.Console.PrintMessage("\nInit Objects...")
# App.closeDocument(App.ActiveDocument.Label) #closeDocument after restart of macro. Needs any ActiveDocument.
# closeDocument after restart of macro. Needs any ActiveDocument.
# App.closeDocument(App.ActiveDocument.Label)
AppPyDoc = App.newDocument(PyDocumentName)
NSideBox = AppPyDoc.addObject("Part::Box", NSideBoxName)
PSideBox = AppPyDoc.addObject("Part::Box", PSideBoxName)
@@ -204,17 +217,33 @@ def createMesh():
AdsorbtionBox = AppPyDoc.addObject("Part::Box", AdsorbtionBoxName)
pnMesh = AppPyDoc.addObject("Mesh::Feature", pnMeshName)
BoxList = [NSideBox, DepletionBox, PSideBox, OxideBox, AdsorbtionBox, SurfDepletionBox]
BoxList = [
NSideBox,
DepletionBox,
PSideBox,
OxideBox,
AdsorbtionBox,
SurfDepletionBox
]
NSideBoxMesh = Mesh.Mesh()
PSideBoxMesh = Mesh.Mesh()
DepletionBoxMesh = Mesh.Mesh()
SurfDepletionBoxMesh = Mesh.Mesh()
OxideBoxMesh = Mesh.Mesh()
AdsorbtionBoxMesh = Mesh.Mesh()
BoxMeshList = [NSideBoxMesh, DepletionBoxMesh, PSideBoxMesh, OxideBoxMesh, AdsorbtionBoxMesh, SurfDepletionBoxMesh]
BoxMeshList = [
NSideBoxMesh,
DepletionBoxMesh,
PSideBoxMesh,
OxideBoxMesh,
AdsorbtionBoxMesh,
SurfDepletionBoxMesh
]
if beVerbose == 1:
if len(BoxList) != len(BoxMeshList):
FreeCAD.Console.PrintMessage("\n ERROR! Input len() of BoxList and BoxMeshList is not the same! ")
FreeCAD.Console.PrintMessage(
"\n ERROR! Input len() of BoxList and BoxMeshList is not the same! "
)
## Set sizes in nanometers
if beVerbose == 1:
@@ -253,12 +282,30 @@ def createMesh():
# Object placement
Rot = App.Rotation(0, 0, 0, 1)
NSideBox.Placement = App.Placement(App.Vector(0, 0, -BulkHeight), Rot)
PSideBox.Placement = App.Placement(App.Vector(DepletionSize * 2 + BulkLength, 0, -BulkHeight), Rot)
DepletionBox.Placement = App.Placement(App.Vector(BulkLength, 0, -BulkHeight), Rot)
SurfDepletionBox.Placement = App.Placement(App.Vector(0, 0, 0), Rot)
OxideBox.Placement = App.Placement(App.Vector(0, 0, DepletionSize), Rot)
AdsorbtionBox.Placement = App.Placement(App.Vector(0, 0, DepletionSize + OxideThickness), Rot)
NSideBox.Placement = App.Placement(
App.Vector(0, 0, -BulkHeight),
Rot
)
PSideBox.Placement = App.Placement(
App.Vector(DepletionSize * 2 + BulkLength, 0, -BulkHeight),
Rot
)
DepletionBox.Placement = App.Placement(
App.Vector(BulkLength, 0, -BulkHeight),
Rot
)
SurfDepletionBox.Placement = App.Placement(
App.Vector(0, 0, 0),
Rot
)
OxideBox.Placement = App.Placement(
App.Vector(0, 0, DepletionSize),
Rot
)
AdsorbtionBox.Placement = App.Placement(
App.Vector(0, 0, DepletionSize + OxideThickness),
Rot
)
## Unite
if beVerbose == 1:
@@ -271,7 +318,9 @@ def createMesh():
# for index in range(len(BoxList)):
for index in range(len(BoxList) - 1): # Manual hack
BoxMeshList[index].addFacets(BoxList[index].Shape.tessellate(tessellationTollerance))
BoxMeshList[index].addFacets(
BoxList[index].Shape.tessellate(tessellationTollerance)
)
nmesh.addMesh(BoxMeshList[index])
nmesh.removeDuplicatedPoints()

18
src/Mod/Fem/feminout/importCcxDatResults.py
View File

@@ -44,13 +44,18 @@ elif open.__module__ == 'io':
pyopen = open
def open(filename):
def open(
filename
):
"called when freecad opens a file"
docname = os.path.splitext(os.path.basename(filename))[0]
insert(filename, docname)
def insert(filename, docname):
def insert(
filename,
docname
):
"called when freecad wants to import a file"
try:
doc = FreeCAD.getDocument(docname)
@@ -61,14 +66,19 @@ def insert(filename, docname):
# ********* module specific methods *********
def import_dat(filename, Analysis=None):
def import_dat(
filename,
Analysis=None
):
r = readResult(filename)
# print("Results {}".format(r))
return r
# read a calculix result file and extract the data
def readResult(dat_input):
def readResult(
dat_input
):
FreeCAD.Console.PrintMessage('Read ccx results from dat file: {}\n'.format(dat_input))
dat_file = pyopen(dat_input, "r")
eigenvalue_output_section_found = False

167
src/Mod/Fem/feminout/importCcxFrdResults.py
View File

@@ -49,7 +49,10 @@ def open(filename):
insert(filename, docname)
def insert(filename, docname):
def insert(
filename,
docname
):
"called when freecad wants to import a file"
try:
doc = FreeCAD.getDocument(docname)
@@ -60,7 +63,11 @@ def insert(filename, docname):
# ********* module specific methods *********
def importFrd(filename, analysis=None, result_name_prefix=None):
def importFrd(
filename,
analysis=None,
result_name_prefix=None
):
from . import importToolsFem
import ObjectsFem
if result_name_prefix is None:
@@ -72,11 +79,16 @@ def importFrd(filename, analysis=None, result_name_prefix=None):
analysis_object = analysis
mesh = importToolsFem.make_femmesh(m)
result_mesh_object = ObjectsFem.makeMeshResult(FreeCAD.ActiveDocument, 'Result_mesh')
result_mesh_object = ObjectsFem.makeMeshResult(
FreeCAD.ActiveDocument,
'Result_mesh'
)
result_mesh_object.FemMesh = mesh
number_of_increments = len(m['Results'])
FreeCAD.Console.PrintLog('Increments: ' + str(number_of_increments) + '\n')
FreeCAD.Console.PrintLog(
'Increments: ' + str(number_of_increments) + '\n'
)
if len(m['Results']) > 0:
for result_set in m['Results']:
if 'number' in result_set:
@@ -86,11 +98,20 @@ def importFrd(filename, analysis=None, result_name_prefix=None):
step_time = result_set['time']
step_time = round(step_time, 2)
if eigenmode_number > 0:
results_name = result_name_prefix + 'mode_' + str(eigenmode_number) + '_results'
results_name = (
'{}mode_{}_results'
.format(result_name_prefix, eigenmode_number)
)
elif number_of_increments > 1:
results_name = result_name_prefix + 'time_' + str(step_time) + '_results'
results_name = (
'{}time_{}_results'
.format(result_name_prefix, step_time)
)
else:
results_name = result_name_prefix + 'results'
results_name = (
'{}results'
.format(result_name_prefix)
)
res_obj = ObjectsFem.makeResultMechanical(FreeCAD.ActiveDocument, results_name)
res_obj.Mesh = result_mesh_object
@@ -101,16 +122,23 @@ def importFrd(filename, analysis=None, result_name_prefix=None):
import femresult.resulttools as restools
if not res_obj.MassFlowRate:
# only compact result if not Flow 1D results
# compact result object, workaround for bug 2873, https://www.freecadweb.org/tracker/view.php?id=2873
# compact result object, workaround for bug 2873
# https://www.freecadweb.org/tracker/view.php?id=2873
res_obj = restools.compact_result(res_obj)
res_obj = restools.add_disp_apps(res_obj) # fill DisplacementLengths
res_obj = restools.add_von_mises(res_obj) # fill StressValues
res_obj = restools.add_principal_stress(res_obj) # fill PrincipalMax, PrincipalMed, PrincipalMin, MaxShear
res_obj = restools.fill_femresult_stats(res_obj) # fill Stats
# fill DisplacementLengths
res_obj = restools.add_disp_apps(res_obj)
# fill StressValues
res_obj = restools.add_von_mises(res_obj)
# fill PrincipalMax, PrincipalMed, PrincipalMin, MaxShear
res_obj = restools.add_principal_stress(res_obj)
# fill Stats
res_obj = restools.fill_femresult_stats(res_obj)
else:
error_message = (
"We have nodes but no results in frd file, which means we only have a mesh in frd file. "
"Usually this happens for analysis type 'NOANALYSIS' or if CalculiX returned no results because "
"We have nodes but no results in frd file, "
"which means we only have a mesh in frd file. "
"Usually this happens for analysis type 'NOANALYSIS' "
"or if CalculiX returned no results because "
"of nonpositive jacobian determinant in at least one element.\n"
)
FreeCAD.Console.PrintMessage(error_message)
@@ -124,13 +152,21 @@ def importFrd(filename, analysis=None, result_name_prefix=None):
FreeCAD.ActiveDocument.recompute()
else:
FreeCAD.Console.PrintError('Problem on frd file import. No nodes found in frd file.\n')
FreeCAD.Console.PrintError(
'Problem on frd file import. No nodes found in frd file.\n'
)
return res_obj
# read a calculix result file and extract the nodes, displacement vectors and stress values.
def read_frd_result(frd_input):
FreeCAD.Console.PrintMessage('Read ccx results from frd file: {}\n'.format(frd_input))
# read a calculix result file and extract the nodes
# displacement vectors and stress values.
def read_frd_result(
frd_input
):
FreeCAD.Console.PrintMessage(
'Read ccx results from frd file: {}\n'
.format(frd_input)
)
inout_nodes = []
inout_nodes_file = frd_input.rsplit('.', 1)[0] + '_inout_nodes.txt'
if os.path.exists(inout_nodes_file):
@@ -247,7 +283,8 @@ def read_frd_result(frd_input):
elif elemType == 4 and input_continues is False:
# first line
# C3D20 Calculix --> hexa20 FreeCAD
# N6, N7, N8, N5, N2, N3, N4, N1, N14, N15, N16, N13, N10, N11, N12, N9, N18, N19, N20, N17
# N6, N7, N8, N5, N2, N3, N4, N1, N14, N15
# N16, N13, N10, N11, N12, N9, N18, N19, N20, N17
nd1 = int(line[3:13])
nd2 = int(line[13:23])
nd3 = int(line[23:33])
@@ -272,18 +309,31 @@ def read_frd_result(frd_input):
nd19 = int(line[83:93])
nd20 = int(line[93:103])
input_continues = False
# CalculiX uses a different node order in input file *.inp and result file *.frd for hexa20 (C3D20)
# according to Guido (the developer of ccx), see note in in first line of cgx manuel part element types
# ccx (and thus the *.inp) follows the ABAQUS convention (documented in the ccx-documentation)
# cgx (and thus the *.frd) follows the FAM2 convention (documented in the cgx-documentation)
# FAM32 is from the company FEGS limited, maybe this company does not exist any more)
# elements_hexa20[elem] = (nd6, nd7, nd8, nd5, nd2, nd3, nd4, nd1, nd14, nd15,
# nd16, nd13, nd10, nd11, nd12, nd9, nd18, nd19, nd20, nd17)
# elements_hexa20[elem] = (nd6, nd7, nd8, nd5, nd2, nd3, nd4, nd1, nd14, nd15,
# nd16, nd13, nd18, nd19, nd20, nd17, nd10, nd11, nd12, nd9)
# hexa20 import works with the following frd file node assignment
elements_hexa20[elem] = (nd8, nd5, nd6, nd7, nd4, nd1, nd2, nd3, nd20, nd17,
nd18, nd19, nd12, nd9, nd10, nd11, nd16, nd13, nd14, nd15)
'''
CalculiX uses a different node order in
input file *.inp and result file *.frd for hexa20 (C3D20)
according to Guido (the developer of ccx):
see note in in first line of cgx manuel part element types
ccx (and thus the *.inp) follows the ABAQUS convention
documented in the ccx-documentation
cgx (and thus the *.frd) follows the FAM2 convention
documented in the cgx-documentation
FAM32 is from the company FEGS limited
maybe this company does not exist any more
elements_hexa20[elem] = (
nd6, nd7, nd8, nd5, nd2, nd3, nd4, nd1, nd14, nd15,
nd16, nd13, nd10, nd11, nd12, nd9, nd18, nd19, nd20, nd17
)
elements_hexa20[elem] = (
nd6, nd7, nd8, nd5, nd2, nd3, nd4, nd1, nd14, nd15,
nd16, nd13, nd18, nd19, nd20, nd17, nd10, nd11, nd12, nd9
)
hexa20 import works with the following frd file node assignment
'''
elements_hexa20[elem] = (
nd8, nd5, nd6, nd7, nd4, nd1, nd2, nd3, nd20, nd17,
nd18, nd19, nd12, nd9, nd10, nd11, nd16, nd13, nd14, nd15
)
# print(elements_hexa20[elem])
elif elemType == 5 and input_continues is False:
# first line
@@ -308,11 +358,19 @@ def read_frd_result(frd_input):
nd14 = int(line[33:43])
nd15 = int(line[43:53])
input_continues = False
# CalculiX uses a different node order in input file *.inp and result file *.frd for penta15 (C3D15), see notes at hexa20
# elements_penta15[elem] = (nd5, nd6, nd4, nd2, nd3, nd1, nd11, nd12, nd10, nd8,
# nd9, nd7, nd14, nd15, nd13) # order of the *.inp file
elements_penta15[elem] = (nd5, nd6, nd4, nd2, nd3, nd1, nd14, nd15, nd13, nd8,
nd9, nd7, nd11, nd12, nd10)
'''
CalculiX uses a different node order in
input file *.inp and result file *.frd for penta15 (C3D15)
see notes at hexa20
elements_penta15[elem] = (
nd5, nd6, nd4, nd2, nd3, nd1, nd11, nd12, nd10, nd8,
nd9, nd7, nd14, nd15, nd13
) # order of the *.inp file
'''
elements_penta15[elem] = (
nd5, nd6, nd4, nd2, nd3, nd1, nd14, nd15, nd13, nd8,
nd9, nd7, nd11, nd12, nd10
)
elif elemType == 6:
# C3D10 Calculix --> tetra10 FreeCAD
# N2, N1, N3, N4, N5, N7, N6, N9, N8, N10
@@ -373,7 +431,9 @@ def read_frd_result(frd_input):
elif elemType == 12:
# B32 CalculiX --> seg3 FreeCAD
# Also D element element number
# CalculiX uses a different node order in input file *.inp and result file *.frd for seg3 (B32), see notes at hexa20
# CalculiX uses a different node order in
# input file *.inp and result file *.frd for seg3 (B32)
# see notes at hexa20
# N1, N2 ,N3
nd1 = int(line[3:13])
nd2 = int(line[13:23])
@@ -381,11 +441,14 @@ def read_frd_result(frd_input):
if inout_nodes:
for i in range(len(inout_nodes)):
if nd1 == int(inout_nodes[i][1]):
elements_seg3[elem] = (int(inout_nodes[i][2]), nd3, nd1) # fluid inlet node numbering
# fluid inlet node numbering
elements_seg3[elem] = (int(inout_nodes[i][2]), nd3, nd1)
elif nd3 == int(inout_nodes[i][1]):
elements_seg3[elem] = (nd1, int(inout_nodes[i][2]), nd3) # fluid outlet node numbering
# fluid outlet node numbering
elements_seg3[elem] = (nd1, int(inout_nodes[i][2]), nd3)
else:
elements_seg3[elem] = (nd1, nd2, nd3) # normal node numbering for D, B32 elements
# normal node numbering for D, B32 elements
elements_seg3[elem] = (nd1, nd2, nd3)
# Check if we found new eigenmode line
if line[5:10] == "PMODE":
@@ -399,7 +462,8 @@ def read_frd_result(frd_input):
mode_time_found = True
if mode_time_found and (line[2:7] == "100CL"):
# we found the new time step line
# !!! be careful here, there is timetemp and timestep! TODO: use more differ names
# !!! be careful here, there is timetemp and timestep!
# TODO: use more differ names
timetemp = float(line[13:25])
if timetemp > timestep:
timestep = timetemp
@@ -480,7 +544,8 @@ def read_frd_result(frd_input):
for i in range(len(inout_nodes)):
if elem == int(inout_nodes[i][1]):
node = int(inout_nodes[i][2])
mode_massflow[node] = (massflow * 1000) # convert units to kg/s from t/s
# convert units to kg/s from t/s
mode_massflow[node] = (massflow * 1000)
# Check if we found a network pressure section
if line[5:11] == "STPRES":
@@ -565,7 +630,9 @@ def read_frd_result(frd_input):
if line[1:5] == "9999":
end_of_frd_data_found = True
if (mode_eigen_changed or mode_time_changed or end_of_frd_data_found) and end_of_section_found and not node_element_section:
if (mode_eigen_changed or mode_time_changed or end_of_frd_data_found) \
and end_of_section_found \
and not node_element_section:
'''
print('\n\n----Append mode_results to results')
@@ -581,11 +648,13 @@ def read_frd_result(frd_input):
# append mode_results to results and reset mode_result
results.append(mode_results)
mode_results = {}
mode_results['number'] = float('NaN') # https://forum.freecadweb.org/viewtopic.php?f=18&t=32649&start=10#p274686
# https://forum.freecadweb.org/viewtopic.php?f=18&t=32649&start=10#p274686
mode_results['number'] = float('NaN')
mode_results['time'] = float('NaN')
end_of_section_found = False
# on changed --> write changed values in mode_result --> will be the first to do on an empty mode_result
# on changed --> write changed values in mode_result
# will be the first to do on an empty mode_result
if mode_eigen_changed:
mode_results['number'] = eigenmode
mode_eigen_changed = False
@@ -596,7 +665,8 @@ def read_frd_result(frd_input):
mode_time_found = False
mode_time_changed = False
# here we are in the indent of loop for every line in frd file, do not add a print here :-)
# here we are in the indent of loop for every line in frd file
# do not add a print here :-)
# close frd file if loop over all lines is finished
frd_file.close()
@@ -613,9 +683,12 @@ def read_frd_result(frd_input):
if not inout_nodes:
if results:
if 'mflow' in results[0] or 'npressure' in results[0]:
FreeCAD.Console.PrintError('We have mflow or npressure, but no inout_nodes file.\n')
FreeCAD.Console.PrintError(
'We have mflow or npressure, but no inout_nodes file.\n'
)
if not nodes:
FreeCAD.Console.PrintError('FEM: No nodes found in Frd file.\n')
return {
'Nodes': nodes,
'Seg2Elem': elements_seg2,

87
src/Mod/Fem/feminout/importToolsFem.py
View File

@@ -31,7 +31,9 @@ __url__ = "http://www.freecadweb.org"
import FreeCAD
def get_FemMeshObjectMeshGroups(fem_mesh_obj):
def get_FemMeshObjectMeshGroups(
fem_mesh_obj
):
"""
Get mesh groups from mesh. This also throws no exception if there
is no Groups property at all (e.g. Netgen meshes).
@@ -45,7 +47,9 @@ def get_FemMeshObjectMeshGroups(fem_mesh_obj):
return gmshgroups
def get_FemMeshObjectOrder(fem_mesh_obj):
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 ->
@@ -72,7 +76,9 @@ def get_FemMeshObjectOrder(fem_mesh_obj):
return presumable_order
def get_FemMeshObjectDimension(fem_mesh_obj):
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)
"""
@@ -90,7 +96,10 @@ def get_FemMeshObjectDimension(fem_mesh_obj):
return dim
def get_FemMeshObjectElementTypes(fem_mesh_obj, remove_zero_element_entries=True):
def get_FemMeshObjectElementTypes(
fem_mesh_obj,
remove_zero_element_entries=True
):
"""
Spit out all elements in the mesh with their appropriate dimension.
"""
@@ -99,7 +108,9 @@ def get_FemMeshObjectElementTypes(fem_mesh_obj, remove_zero_element_entries=True
"Prism": 3, "Pyramid": 3, "Quadrangle": 2, "Tetra": 3, "Triangle": 2}
eval_dict = locals() # to access local variables from eval
elements_list_with_zero = [(eval("fem_mesh_obj.FemMesh." + s + "Count", eval_dict), s, d) for (s, d) in FreeCAD_element_names_dims.items()]
elements_list_with_zero = [(
eval("fem_mesh_obj.FemMesh." + s + "Count", eval_dict), s, d
) for (s, d) in FreeCAD_element_names_dims.items()]
# 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]
@@ -109,7 +120,9 @@ def get_FemMeshObjectElementTypes(fem_mesh_obj, remove_zero_element_entries=True
return elements_list
def get_MaxDimElementFromList(elem_list):
def get_MaxDimElementFromList(
elem_list
):
"""
Gets element with the maximal dimension in the mesh to determine cells.
"""
@@ -117,7 +130,9 @@ def get_MaxDimElementFromList(elem_list):
return elem_list[-1]
def make_femmesh(mesh_data):
def make_femmesh(
mesh_data
):
''' makes an FreeCAD FEM Mesh object from FEM Mesh data
'''
import Fem
@@ -169,8 +184,10 @@ def make_femmesh(mesh_data):
elms_hexa20 = m['Hexa20Elem']
for i in elms_hexa20:
e = elms_hexa20[i]
mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
e[10], e[11], e[12], e[13], e[14], e[15], e[16], e[17], e[18], e[19]], i)
mesh.addVolume([
e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
e[10], e[11], e[12], e[13], e[14], e[15], e[16], e[17], e[18], e[19]
], i)
elms_tria3 = m['Tria3Elem']
for i in elms_tria3:
e = elms_tria3[i]
@@ -195,12 +212,29 @@ def make_femmesh(mesh_data):
for i in elms_seg3:
e = elms_seg3[i]
mesh.addEdge([e[0], e[1], e[2]], i)
FreeCAD.Console.PrintLog("imported mesh: {} nodes, {} HEXA8, {} PENTA6, {} TETRA4, {} TETRA10, {} PENTA15".format(
len(nds), len(elms_hexa8), len(elms_penta6), len(elms_tetra4), len(elms_tetra10), len(elms_penta15)
))
FreeCAD.Console.PrintLog("imported mesh: {} HEXA20, {} TRIA3, {} TRIA6, {} QUAD4, {} QUAD8, {} SEG2, {} SEG3".format(
len(elms_hexa20), len(elms_tria3), len(elms_tria6), len(elms_quad4), len(elms_quad8), len(elms_seg2), len(elms_seg3)
))
FreeCAD.Console.PrintLog(
"imported mesh: {} nodes, {} HEXA8, {} PENTA6, {} TETRA4, {} TETRA10, {} PENTA15"
.format(
len(nds),
len(elms_hexa8),
len(elms_penta6),
len(elms_tetra4),
len(elms_tetra10),
len(elms_penta15)
)
)
FreeCAD.Console.PrintLog(
"imported mesh: {} HEXA20, {} TRIA3, {} TRIA6, {} QUAD4, {} QUAD8, {} SEG2, {} SEG3"
.format(
len(elms_hexa20),
len(elms_tria3),
len(elms_tria6),
len(elms_quad4),
len(elms_quad8),
len(elms_seg2),
len(elms_seg3)
)
)
else:
FreeCAD.Console.PrintError("No Elements found!\n")
else:
@@ -208,14 +242,18 @@ def make_femmesh(mesh_data):
return mesh
def fill_femresult_mechanical(res_obj, result_set):
def fill_femresult_mechanical(
res_obj,
result_set
):
''' fills a FreeCAD FEM mechanical result object with result data
'''
if 'time' in result_set:
step_time = result_set['time']
step_time = round(step_time, 2)
# if disp exists, fill res_obj.NodeNumbers and res_obj.DisplacementVectors as well as stress and strain
# if disp exists, fill res_obj.NodeNumbers and
# res_obj.DisplacementVectors as well as stress and strain
if 'disp' in result_set:
disp = result_set['disp']
res_obj.DisplacementVectors = list(map((lambda x: x), disp.values()))
@@ -223,7 +261,8 @@ def fill_femresult_mechanical(res_obj, result_set):
# fill res_obj.NodeStressXX etc if they exist in result_set
# list values are just added
# Should we check if the key in stress and strain dict is the same as the number in NodeNumbers?
# Should we check if the key in stress and strain dict
# is the same as the number in NodeNumbers?
if 'stress' in result_set:
stress = result_set['stress']
Sxx = []
@@ -232,7 +271,8 @@ def fill_femresult_mechanical(res_obj, result_set):
Sxy = []
Sxz = []
Syz = []
for i, values_S in enumerate(stress.values()): # values_S .. stress_tensor .. (Sxx, Syy, Szz, Sxy, Sxz, Syz)
# values_S .. stress_tensor .. (Sxx, Syy, Szz, Sxy, Sxz, Syz)
for i, values_S in enumerate(stress.values()):
Sxx.append(values_S[0])
Syy.append(values_S[1])
Szz.append(values_S[2])
@@ -255,7 +295,8 @@ def fill_femresult_mechanical(res_obj, result_set):
Exy = []
Exz = []
Eyz = []
for i, values_E in enumerate(strain.values()): # values_E .. straintuple .. (Exx, Eyy, Ezz, Exy, Exz, Eyz)
# values_E .. straintuple .. (Exx, Eyy, Ezz, Exy, Exz, Eyz)
for i, values_E in enumerate(strain.values()):
Exx.append(values_E[0])
Eyy.append(values_E[1])
Ezz.append(values_E[2])
@@ -287,7 +328,8 @@ def fill_femresult_mechanical(res_obj, result_set):
res_obj.Peeq = list(Peeq.values())
# fill res_obj.Temperature if they exist
# TODO, check if it is possible to have Temperature without disp, we would need to set NodeNumbers than
# TODO, check if it is possible to have Temperature without disp
# we would need to set NodeNumbers than
if 'temp' in result_set:
Temperature = result_set['temp']
if len(Temperature) > 0:
@@ -311,7 +353,8 @@ def fill_femresult_mechanical(res_obj, result_set):
if len(MassFlow) > 0:
res_obj.MassFlowRate = list(map((lambda x: x), MassFlow.values()))
res_obj.Time = step_time
res_obj.NodeNumbers = list(MassFlow.keys()) # disp does not exist, res_obj.NodeNumbers needs to be set
# disp does not exist, res_obj.NodeNumbers needs to be set
res_obj.NodeNumbers = list(MassFlow.keys())
# fill res_obj.NetworkPressure, disp does not exist, see MassFlow
if 'npressure' in result_set:

65
src/Mod/Fem/feminout/importVTKResults.py
View File

@@ -44,13 +44,18 @@ elif open.__module__ == 'io':
pyopen = open
def open(filename):
def open(
filename
):
"called when freecad opens a file"
docname = os.path.splitext(os.path.basename(filename))[0]
insert(filename, docname)
def insert(filename, docname):
def insert(
filename,
docname
):
"called when freecad wants to import a file"
try:
doc = FreeCAD.getDocument(docname)
@@ -60,30 +65,47 @@ def insert(filename, docname):
importVtk(filename)
def export(objectslist, filename):
def export(
objectslist,
filename
):
"called when freecad exports an object to vtk"
if len(objectslist) > 1: # the case of no selected obj is caught by FreeCAD already
FreeCAD.Console.PrintError("This exporter can only export one object at once\n")
FreeCAD.Console.PrintError(
"This exporter can only export one object at once\n"
)
return
obj = objectslist[0]
if obj.isDerivedFrom("Fem::FemPostPipeline"):
FreeCAD.Console.PrintError('Export of a VTK post object to vtk is not yet implemented !\n')
FreeCAD.Console.PrintError(
'Export of a VTK post object to vtk is not yet implemented !\n'
)
return
elif obj.isDerivedFrom("Fem::FemMeshObject"):
FreeCAD.Console.PrintError('Use export to FEM mesh formats to export a FEM mesh object to vtk!\n')
FreeCAD.Console.PrintError(
'Use export to FEM mesh formats to export a FEM mesh object to vtk!\n'
)
return
elif obj.isDerivedFrom("Fem::FemResultObject"):
Fem.writeResult(filename, obj)
else:
FreeCAD.Console.PrintError('Selected object is not supported by export to VTK.\n')
FreeCAD.Console.PrintError(
'Selected object is not supported by export to VTK.\n'
)
return
# ********* module specific methods *********
def importVtk(filename, object_name=None, object_type=None):
def importVtk(
filename,
object_name=None,
object_type=None
):
if not object_type:
vtkinout_prefs = FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Fem/InOutVtk")
vtkinout_prefs = FreeCAD.ParamGet(
"User parameter:BaseApp/Preferences/Mod/Fem/InOutVtk"
)
object_type = vtkinout_prefs.GetInt("ImportObject", 0)
if not object_name:
object_name = os.path.splitext(os.path.basename(filename))[0]
@@ -97,10 +119,16 @@ def importVtk(filename, object_name=None, object_type=None):
# FreeCAD result object
importVtkFCResult(filename, object_name)
else:
FreeCAD.Console.PrintError('Error, wrong parameter in VTK import pref: {}\n'.format(object_type))
FreeCAD.Console.PrintError(
'Error, wrong parameter in VTK import pref: {}\n'
.format(object_type)
)
def importVtkVtkResult(filename, resultname):
def importVtkVtkResult(
filename,
resultname
):
vtk_result_obj = FreeCAD.ActiveDocument.addObject("Fem::FemPostPipeline", resultname)
vtk_result_obj.read(filename)
vtk_result_obj.touch()
@@ -108,7 +136,10 @@ def importVtkVtkResult(filename, resultname):
return vtk_result_obj
def importVtkFemMesh(filename, meshname):
def importVtkFemMesh(
filename,
meshname
):
meshobj = FreeCAD.ActiveDocument.addObject("Fem::FemMeshObject", meshname)
meshobj.FemMesh = Fem.read(filename)
meshobj.touch()
@@ -116,7 +147,12 @@ def importVtkFemMesh(filename, meshname):
return meshobj
def importVtkFCResult(filename, resultname, analysis=None, result_name_prefix=None):
def importVtkFCResult(
filename,
resultname,
analysis=None,
result_name_prefix=None
):
# only fields from vtk are imported if they exactly named as the FreeCAD result properties
# See _getFreeCADMechResultProperties() in FemVTKTools.cpp for the supported names
@@ -128,7 +164,8 @@ def importVtkFCResult(filename, resultname, analysis=None, result_name_prefix=No
results_name = result_name_prefix + 'results'
result_obj = ObjectsFem.makeResultMechanical(FreeCAD.ActiveDocument, results_name)
Fem.readResult(filename, result_obj.Name) # readResult always creates a new femmesh named ResultMesh
# readResult always creates a new femmesh named ResultMesh
Fem.readResult(filename, result_obj.Name)
# add missing DisplacementLengths (They should have been added by Fem.readResult)
if not result_obj.DisplacementLengths:

140
src/Mod/Fem/feminout/importZ88Mesh.py
View File

@@ -41,13 +41,18 @@ elif open.__module__ == 'io':
pyopen = open
def open(filename):
def open(
filename
):
"called when freecad opens a file"
docname = os.path.splitext(os.path.basename(filename))[0]
insert(filename, docname)
def insert(filename, docname):
def insert(
filename,
docname
):
"called when freecad wants to import a file"
try:
doc = FreeCAD.getDocument(docname)
@@ -57,7 +62,10 @@ def insert(filename, docname):
import_z88_mesh(filename)
def export(objectslist, filename):
def export(
objectslist,
filename
):
"called when freecad exports a file"
if len(objectslist) != 1:
FreeCAD.Console.PrintError("This exporter can only export one object.\n")
@@ -76,7 +84,10 @@ def export(objectslist, filename):
# ********* module specific methods *********
def write(fem_mesh, filename):
def write(
fem_mesh,
filename
):
'''directly write a FemMesh to a Z88 mesh file format
fem_mesh: a FemMesh'''
@@ -92,7 +103,9 @@ def write(fem_mesh, filename):
f.close()
def read(filename):
def read(
filename
):
'''read a FemMesh from a Z88 mesh file and return the FemMesh
'''
# no document object is created, just the FemMesh is returned
@@ -101,8 +114,12 @@ def read(filename):
return importToolsFem.make_femmesh(mesh_data)
def import_z88_mesh(filename, analysis=None):
'''read a FEM mesh from a Z88 mesh file and insert a FreeCAD FEM Mesh object in the ActiveDocument
def import_z88_mesh(
filename,
analysis=None
):
'''read a FEM mesh from a Z88 mesh file and
insert a FreeCAD FEM Mesh object in the ActiveDocument
'''
femmesh = read(filename)
mesh_name = os.path.basename(os.path.splitext(filename)[0])
@@ -111,7 +128,9 @@ def import_z88_mesh(filename, analysis=None):
mesh_object.FemMesh = femmesh
def read_z88_mesh(z88_mesh_input):
def read_z88_mesh(
z88_mesh_input
):
''' reads a z88 mesh file z88i1.txt (Z88OSV14) or z88structure.txt (Z88AuroraV3)
and extracts the nodes and elements
'''
@@ -140,8 +159,11 @@ def read_z88_mesh(z88_mesh_input):
nodes_count = int(mesh_info[1])
elements_count = int(mesh_info[2])
kflag = int(mesh_info[4])
if kflag: # for non rotational elements ist --> kflag = 0 --> cartesian, kflag = 1 polar coordinates
FreeCAD.Console.PrintError("KFLAG = 1, Rotational coordinates not supported at the moment\n")
# for non rotational elements ist --> kflag = 0 --> cartesian, kflag = 1 polar coordinates
if kflag:
FreeCAD.Console.PrintError(
"KFLAG = 1, Rotational coordinates not supported at the moment\n"
)
return {}
nodes_first_line = 2 # first line is mesh_info
nodes_last_line = nodes_count + 1
@@ -186,33 +208,57 @@ def read_z88_mesh(z88_mesh_input):
# not supported elements
if z88_element_type == 8:
# torus8
FreeCAD.Console.PrintError("Z88 Element No. 8, torus8\n")
FreeCAD.Console.PrintError("Rotational elements are not supported at the moment\n")
FreeCAD.Console.PrintError(
"Z88 Element No. 8, torus8\n"
)
FreeCAD.Console.PrintError(
"Rotational elements are not supported at the moment\n"
)
return {}
elif z88_element_type == 12:
# torus12
FreeCAD.Console.PrintError("Z88 Element No. 12, torus12\n")
FreeCAD.Console.PrintError("Rotational elements are not supported at the moment\n")
FreeCAD.Console.PrintError(
"Z88 Element No. 12, torus12\n"
)
FreeCAD.Console.PrintError(
"Rotational elements are not supported at the moment\n"
)
return {}
elif z88_element_type == 15:
# torus6
FreeCAD.Console.PrintError("Z88 Element No. 15, torus6\n")
FreeCAD.Console.PrintError("Rotational elements are not supported at the moment\n")
FreeCAD.Console.PrintError(
"Z88 Element No. 15, torus6\n"
)
FreeCAD.Console.PrintError(
"Rotational elements are not supported at the moment\n"
)
return {}
elif z88_element_type == 19:
# platte16
FreeCAD.Console.PrintError("Z88 Element No. 19, platte16\n")
FreeCAD.Console.PrintError("Not supported at the moment\n")
FreeCAD.Console.PrintError(
"Z88 Element No. 19, platte16\n"
)
FreeCAD.Console.PrintError(
"Not supported at the moment\n"
)
return {}
elif z88_element_type == 21:
# schale16, mixture made from hexa8 and hexa20 (thickness is linear)
FreeCAD.Console.PrintError("Z88 Element No. 21, schale16\n")
FreeCAD.Console.PrintError("Not supported at the moment\n")
FreeCAD.Console.PrintError(
"Z88 Element No. 21, schale16\n"
)
FreeCAD.Console.PrintError(
"Not supported at the moment\n"
)
return {}
elif z88_element_type == 22:
# schale12, mixtrue made from prism6 and prism15 (thickness is linear)
FreeCAD.Console.PrintError("Z88 Element No. 22, schale12\n")
FreeCAD.Console.PrintError("Not supported at the moment\n")
FreeCAD.Console.PrintError(
"Z88 Element No. 22, schale12\n"
)
FreeCAD.Console.PrintError(
"Not supported at the moment\n"
)
return {}
# supported elements
@@ -263,7 +309,8 @@ def read_z88_mesh(z88_mesh_input):
input_continues = False
elif z88_element_type == 16:
# volume16 Z88 --> tetra10 FreeCAD
# N1, N2, N4, N3, N5, N8, N10, N7, N6, N9, , Z88 to FC is different as FC to Z88
# N1, N2, N4, N3, N5, N8, N10, N7, N6, N9
# Z88 to FC is different as FC to Z88
nd1 = int(linecolumns[0])
nd2 = int(linecolumns[1])
nd3 = int(linecolumns[2])
@@ -274,7 +321,9 @@ def read_z88_mesh(z88_mesh_input):
nd8 = int(linecolumns[7])
nd9 = int(linecolumns[8])
nd10 = int(linecolumns[9])
elements_tetra10[elem_no] = (nd1, nd2, nd4, nd3, nd5, nd8, nd10, nd7, nd6, nd9)
elements_tetra10[elem_no] = (
nd1, nd2, nd4, nd3, nd5, nd8, nd10, nd7, nd6, nd9
)
input_continues = False
elif z88_element_type == 1:
# volume1 Z88 --> hexa8 FreeCAD
@@ -291,9 +340,11 @@ def read_z88_mesh(z88_mesh_input):
input_continues = False
elif z88_element_type == 10:
# volume10 Z88 --> hexa20 FreeCAD
# N2, N3, N4, N1, N6, N7, N8, N5, N10, N11, N12, N9, N14, N15, N16, N13, N18, N19, N20, N17
# N2, N3, N4, N1, N6, N7, N8, N5, N10, N11
# N12, N9, N14, N15, N16, N13, N18, N19, N20, N17
# or turn by 90 degree and they match !
# N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, N11, N12, N13, N14, N15, N16, N17, N18, N19, N20
# N1, N2, N3, N4, N5, N6, N7, N8, N9, N10
# N11, N12, N13, N14, N15, N16, N17, N18, N19, N20
nd1 = int(linecolumns[0])
nd2 = int(linecolumns[1])
nd3 = int(linecolumns[2])
@@ -314,11 +365,14 @@ def read_z88_mesh(z88_mesh_input):
nd18 = int(linecolumns[17])
nd19 = int(linecolumns[18])
nd20 = int(linecolumns[19])
elements_hexa20[elem_no] = (nd1, nd2, nd3, nd4, nd5, nd6, nd7, nd8, nd9, nd10,
nd11, nd12, nd13, nd14, nd15, nd16, nd17, nd18, nd19, nd20)
elements_hexa20[elem_no] = (
nd1, nd2, nd3, nd4, nd5, nd6, nd7, nd8, nd9, nd10,
nd11, nd12, nd13, nd14, nd15, nd16, nd17, nd18, nd19, nd20
)
input_continues = False
# unknown elements, some examples have -1 for some teaching reasons to show some other stuff
# unknown elements
# some examples have -1 for some teaching reasons to show some other stuff
else:
FreeCAD.Console.PrintError("Unknown element\n")
return {}
@@ -348,7 +402,12 @@ def read_z88_mesh(z88_mesh_input):
# write z88 Mesh
def write_z88_mesh_to_file(femnodes_mesh, femelement_table, z88_element_type, f):
def write_z88_mesh_to_file(
femnodes_mesh,
femelement_table,
z88_element_type,
f
):
node_dimension = 3 # 2 for 2D not supported
if (
z88_element_type == 4
@@ -431,24 +490,33 @@ def write_z88_mesh_to_file(femnodes_mesh, femelement_table, z88_element_type, f)
n[0], n[1], n[2], n[3], n[4], n[5], n[6], n[7]))
elif z88_element_type == 10:
# hexa20 FreeCAD --> volume10 Z88
# N2, N3, N4, N1, N6, N7, N8, N5, N10, N11, N12, N9, N14, N15, N16, N13, N18, N19, N20, N17
# N2, N3, N4, N1, N6, N7, N8, N5, N10, N11
# N12, N9, N14, N15, N16, N13, N18, N19, N20, N17
# or turn by 90 degree and they match !
# N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, N11, N12, N13, N14, N15, N16, N17, N18, N19, N20
# N1, N2, N3, N4, N5, N6, N7, N8, N9, N10
# N11, N12, N13, N14, N15, N16, N17, N18, N19, N20
f.write("{0} {1}\n".format(element, z88_element_type, element))
f.write(
"{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11} {12} {13} {14} {15} {16} {17} {18} {19}\n"
"{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} "
"{10} {11} {12} {13} {14} {15} {16} {17} {18} {19}\n"
.format(
n[0], n[1], n[2], n[3], n[4], n[5], n[6], n[7], n[8], n[9], n[10], n[11], n[12], n[13], n[14], n[15], n[16], n[17], n[18], n[19]
n[0], n[1], n[2], n[3], n[4], n[5], n[6], n[7], n[8], n[9],
n[10], n[11], n[12], n[13], n[14], n[15], n[16], n[17], n[18], n[19]
)
)
else:
FreeCAD.Console.PrintError("Writing of Z88 elementtype {0} not supported.\n".format(z88_element_type))
FreeCAD.Console.PrintError(
"Writing of Z88 elementtype {0} not supported.\n".format(z88_element_type)
)
# TODO support schale12 (made from prism15) and schale16 (made from hexa20)
return
# Helper
def get_z88_element_type(femmesh, femelement_table=None):
def get_z88_element_type(
femmesh,
femelement_table=None
):
import femmesh.meshtools as FemMeshTools
if not femmesh:
print("Error: No femmesh!")

31
src/Mod/Fem/feminout/importZ88O2Results.py
View File

@@ -44,13 +44,18 @@ elif open.__module__ == 'io':
pyopen = open
def open(filename):
def open(
filename
):
"called when freecad opens a file"
docname = os.path.splitext(os.path.basename(filename))[0]
insert(filename, docname)
def insert(filename, docname):
def insert(
filename,
docname
):
"called when freecad wants to import a file"
try:
doc = FreeCAD.getDocument(docname)
@@ -61,14 +66,19 @@ def insert(filename, docname):
# ********* module specific methods *********
def import_z88_disp(filename, analysis=None, result_name_prefix=None):
def import_z88_disp(
filename,
analysis=None,
result_name_prefix=None
):
'''insert a FreeCAD FEM mechanical result object in the ActiveDocument
pure usage:
import feminout.importZ88O2Results as importZ88O2Results
disp_file = '/pathtofile/z88o2.txt'
importZ88O2Results.import_z88_disp(disp_file)
the z888i1.txt FEMMesh file needs to be in the same directory as z88o2.txt (ahh, make a new document first ;-))
the z888i1.txt FEMMesh file needs to be in the same directory as z88o2.txt
# ahh, make a new document first ;-)
'''
from . import importZ88Mesh
from . import importToolsFem
@@ -86,7 +96,10 @@ def import_z88_disp(filename, analysis=None, result_name_prefix=None):
mesh_file = filename.replace('o2', 'i1')
mesh_data = importZ88Mesh.read_z88_mesh(mesh_file)
femmesh = importToolsFem.make_femmesh(mesh_data)
result_mesh_object = ObjectsFem.makeMeshResult(FreeCAD.ActiveDocument, 'Result_mesh')
result_mesh_object = ObjectsFem.makeMeshResult(
FreeCAD.ActiveDocument,
'Result_mesh'
)
result_mesh_object.FemMesh = femmesh
else:
FreeCAD.Console.PrintError('Z88 mesh file z88i1.txt not found!')
@@ -108,10 +121,14 @@ def import_z88_disp(filename, analysis=None, result_name_prefix=None):
FreeCAD.ActiveDocument.recompute()
else:
FreeCAD.Console.PrintError('Problem on Z88 result file import. No nodes found in Z88 result file.\n')
FreeCAD.Console.PrintError(
'Problem on Z88 result file import. No nodes found in Z88 result file.\n'
)
def read_z88_disp(z88_disp_input):
def read_z88_disp(
z88_disp_input
):
'''
read a z88 disp file and extract the nodes and elements
z88 Displacement output file is z88o2.txt

19
src/Mod/Fem/feminout/readFenicsXDMF.py
View File

@@ -35,8 +35,17 @@ def read_fenics_mesh_xdmf(xdmffilename):
FreeCAD.Console.PrintMessage("Not operational, yet\n")
return {'Nodes': {},
'Hexa8Elem': {}, 'Penta6Elem': {}, 'Tetra4Elem': {}, 'Tetra10Elem': {},
'Penta15Elem': {}, 'Hexa20Elem': {}, 'Tria3Elem': {}, 'Tria6Elem': {},
'Quad4Elem': {}, 'Quad8Elem': {}, 'Seg2Elem': {}
}
return {
'Nodes': {},
'Hexa8Elem': {},
'Penta6Elem': {},
'Tetra4Elem': {},
'Tetra10Elem': {},
'Penta15Elem': {},
'Hexa20Elem': {},
'Tria3Elem': {},
'Tria6Elem': {},
'Quad4Elem': {},
'Quad8Elem': {},
'Seg2Elem': {}
}

13
src/Mod/Fem/feminout/readFenicsXML.py
View File

@@ -86,7 +86,9 @@ def read_fenics_mesh_xml(xmlfilename):
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"]]
[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
@@ -102,10 +104,15 @@ def read_fenics_mesh_xml(xmlfilename):
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()))
print(
"Strange mismatch between cell type {} and cell tag {}"
.format(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)])
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

97
src/Mod/Fem/feminout/writeFenicsXDMF.py
View File

@@ -63,7 +63,9 @@ FreeCAD_to_Fenics_XDMF_dict = {
# also the hd5 support works better together with numpy
def numpy_array_to_str(npa):
def numpy_array_to_str(
npa
):
res = ""
dt = str(npa.dtype)
if 'int' in dt:
@@ -73,15 +75,25 @@ def numpy_array_to_str(npa):
return res
def points_to_numpy(pts, dim=3):
def points_to_numpy(
pts,
dim=3
):
return np.array([[p.x, p.y, p.z] for p in pts])[:, :dim]
def tuples_to_numpy(tpls, numbers_per_line):
def tuples_to_numpy(
tpls,
numbers_per_line
):
return np.array([list(t) for t in tpls])[:, :numbers_per_line]
def write_fenics_mesh_points_xdmf(fem_mesh_obj, geometrynode, encoding=ENCODING_ASCII):
def write_fenics_mesh_points_xdmf(
fem_mesh_obj,
geometrynode,
encoding=ENCODING_ASCII
):
"""
Writes either into hdf5 file or into open mesh file
"""
@@ -117,19 +129,23 @@ def write_fenics_mesh_points_xdmf(fem_mesh_obj, geometrynode, encoding=ENCODING_
return recalc_nodes_ind_dict
def write_fenics_mesh_codim_xdmf(fem_mesh_obj,
topologynode,
nodes_dict,
codim=0,
encoding=ENCODING_ASCII):
def write_fenics_mesh_codim_xdmf(
fem_mesh_obj,
topologynode,
nodes_dict,
codim=0,
encoding=ENCODING_ASCII
):
mesh_dimension = get_FemMeshObjectDimension(fem_mesh_obj)
element_types = get_FemMeshObjectElementTypes(fem_mesh_obj, remove_zero_element_entries=True)
element_order = get_FemMeshObjectOrder(fem_mesh_obj)
# we get all elements from mesh to decide which one to write by selection of codim
# nodeindices = [
# (nodes_dict[ind] for ind in fem_mesh_obj.FemMesh.getElementNodes(fc_volume_ind)) for (fen_ind, fc_volume_ind) in enumerate(fc_cells)
# ]
'''
nodeindices = [(
nodes_dict[ind] for ind in fem_mesh_obj.FemMesh.getElementNodes(fc_volume_ind)
) for (fen_ind, fc_volume_ind) in enumerate(fc_cells)]
'''
writeout_element_dimension = mesh_dimension - codim
(num_topo, name_topo, dim_topo) = (0, "", 0)
@@ -155,9 +171,9 @@ def write_fenics_mesh_codim_xdmf(fem_mesh_obj,
fc_topo = []
print("Dimension of mesh incompatible with export XDMF function: %d" % (dim_topo,))
nodeindices = [
(nodes_dict[ind] for ind in fem_mesh_obj.FemMesh.getElementNodes(fc_topo_ind)) for (fen_ind, fc_topo_ind) in enumerate(fc_topo)
]
nodeindices = [(
nodes_dict[ind] for ind in fem_mesh_obj.FemMesh.getElementNodes(fc_topo_ind)
) for (fen_ind, fc_topo_ind) in enumerate(fc_topo)]
if encoding == ENCODING_ASCII:
dataitem = ET.SubElement(
@@ -173,7 +189,11 @@ def write_fenics_mesh_codim_xdmf(fem_mesh_obj,
return fc_topo
def write_fenics_mesh_scalar_cellfunctions(name, cell_array, attributenode, encoding=ENCODING_ASCII):
def write_fenics_mesh_scalar_cellfunctions(
name, cell_array,
attributenode,
encoding=ENCODING_ASCII
):
attributenode.set("AttributeType", "Scalar")
attributenode.set("Center", "Cell")
attributenode.set("Name", name)
@@ -230,7 +250,12 @@ Example: mesh with two topologies and one mesh function for the facet one
"""
def write_fenics_mesh_xdmf(fem_mesh_obj, outputfile, group_values_dict={}, encoding=ENCODING_ASCII):
def write_fenics_mesh_xdmf(
fem_mesh_obj,
outputfile,
group_values_dict={},
encoding=ENCODING_ASCII
):
"""
For the export of xdmf.
"""
@@ -255,7 +280,10 @@ def write_fenics_mesh_xdmf(fem_mesh_obj, outputfile, group_values_dict={}, encod
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))
print(
"Celltype in mesh -> {} and its Fenics dolfin name: {}"
.format(celltype_in_mesh, cellname_fenics)
)
root = ET.Element("Xdmf", version="3.0")
domain = ET.SubElement(root, "Domain")
@@ -267,11 +295,20 @@ def write_fenics_mesh_xdmf(fem_mesh_obj, outputfile, group_values_dict={}, encod
# TODO: for every marked group write own grid node with topology (ref if cells)
# geometry ref, attribute
#####################################
# ***********************************
# write base topo and geometry
nodes_dict = write_fenics_mesh_points_xdmf(fem_mesh_obj, base_geometry, encoding=encoding)
write_fenics_mesh_codim_xdmf(fem_mesh_obj, base_topology, nodes_dict, codim=0, encoding=encoding)
#####################################
nodes_dict = write_fenics_mesh_points_xdmf(
fem_mesh_obj,
base_geometry,
encoding=encoding
)
write_fenics_mesh_codim_xdmf(
fem_mesh_obj, base_topology,
nodes_dict,
codim=0,
encoding=encoding
)
# ***********************************
fem_mesh = fem_mesh_obj.FemMesh
gmshgroups = get_FemMeshObjectMeshGroups(fem_mesh_obj)
@@ -310,16 +347,22 @@ def write_fenics_mesh_xdmf(fem_mesh_obj, outputfile, group_values_dict={}, encod
# TODO: is it better to save all groups each at once or collect all codim equal
# groups to put them into one function?
# TODO: nevertheless there has to be a dialog which fixes the default value and the mark value
# TODO: nevertheless there has to be a dialog
# which fixes the default value and the mark value
for e in fem_mesh.getGroupElements(g):
elem_dict[e] = elem_mark_group
val_array = np.array([elem_dict.get(e, elem_mark_default) for e in mesh_function_topology_description])
val_array = np.array([
elem_dict.get(e, elem_mark_default) for e in mesh_function_topology_description
])
topo_array = np.vstack((val_array,)).T
write_fenics_mesh_scalar_cellfunctions(mesh_function_name,
topo_array,
mesh_function_attribute, encoding=ENCODING_ASCII)
write_fenics_mesh_scalar_cellfunctions(
mesh_function_name,
topo_array,
mesh_function_attribute,
encoding=ENCODING_ASCII
)
# TODO: improve cell functions support

14
src/Mod/Fem/feminout/writeFenicsXML.py
View File

@@ -30,14 +30,17 @@ __url__ = "http://www.freecadweb.org"
# \brief FreeCAD Fenics Mesh XML writer for FEM workbench
from .importToolsFem import get_FemMeshObjectDimension, get_FemMeshObjectElementTypes, get_MaxDimElementFromList
from .importToolsFem import get_FemMeshObjectDimension
from .importToolsFem import get_FemMeshObjectElementTypes
from .importToolsFem import get_MaxDimElementFromList
from xml.etree import ElementTree as ET # 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.)
vertices to be exported.
For second order elements, we have to delete the mid element nodes.
"""
# TODO: check for second order elements
@@ -74,7 +77,9 @@ def write_fenics_mesh_xml(fem_mesh_obj, outputfile):
(num_cells, cellname_fc, dim_cell) = celltype_in_mesh
cellname_fenics = FreeCAD_to_Fenics_dict[cellname_fc]
num_verts_cell = XML_Number_of_Nodes_dict[cellname_fenics]
print("Celltype in mesh -> %s and its Fenics name: %s" % (str(celltype_in_mesh), cellname_fenics))
print(
"Celltype in mesh -> %s and its Fenics name: %s" % (str(celltype_in_mesh), cellname_fenics)
)
root = ET.Element("dolfin", dolfin="http://fenicsproject.org")
meshchild = ET.SubElement(root, "mesh", celltype=cellname_fenics, dim=str(dim_cell))
@@ -104,7 +109,8 @@ def write_fenics_mesh_xml(fem_mesh_obj, outputfile):
for (vi, ni) in enumerate(nodeindices):
if vi < num_verts_cell: # XML only supports first order meshs
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)
# generate as many v entries in dict as nodes are listed in cell
# works only for first order elements
ET.SubElement(cells, cellname_fenics, index=str(fen_ind), **cell_args)