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

FEM: code formating, max line length < 150, all FEM python

Browse Source
This commit is contained in:
Bernd Hahnebach
2019-05-20 12:31:36 +02:00
parent 508a65f115
commit 708d1735fd
6 changed files with 469 additions and 148 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/Mod/Fem/ObjectsFem.py
View File

@@ -323,7 +323,8 @@ def makeMeshNetgen(doc, name="FEMMeshNetgen"):
def makeMeshRegion(doc, base_mesh, element_length=0.0, name="FEMMeshRegion"):
'''makeMeshRegion(document, base_mesh, [element_length], [name]): creates a FEM mesh region object to define properties for a region of a FEM mesh'''
'''makeMeshRegion(document, base_mesh, [element_length], [name]):
creates a FEM mesh region object to define properties for a region of a FEM mesh'''
obj = doc.addObject("Fem::FeaturePython", name)
from femobjects import _FemMeshRegion
_FemMeshRegion._FemMeshRegion(obj)

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

@@ -216,7 +216,12 @@ def read_z88_mesh(z88_mesh_input):
return {}
# supported elements
elif z88_element_type == 2 or z88_element_type == 4 or z88_element_type == 5 or z88_element_type == 9 or z88_element_type == 13 or z88_element_type == 25:
elif z88_element_type == 2 \
or z88_element_type == 4 \
or z88_element_type == 5 \
or z88_element_type == 9 \
or z88_element_type == 13 \
or z88_element_type == 25:
# stab4 or stab5 or welle5 or beam13 or beam25 Z88 --> seg2 FreeCAD
# N1, N2
nd1 = int(linecolumns[0])
@@ -368,16 +373,27 @@ def write_z88_mesh_to_file(femnodes_mesh, femelement_table, z88_element_type, f)
written_by = "written by FreeCAD"
# first line, some z88 specific stuff
f.write("{0} {1} {2} {3} {4} {5}\n".format(node_dimension, node_count, element_count, dofs, unknown_flag, written_by))
f.write("{0} {1} {2} {3} {4} {5}\n".format(
node_dimension, node_count, element_count, dofs, unknown_flag, written_by)
)
# nodes
for node in femnodes_mesh:
vec = femnodes_mesh[node]
f.write("{0} {1} {2:.6f} {3:.6f} {4:.6f}\n".format(node, node_dof, vec.x, vec.y, vec.z, node))
f.write(
"{0} {1} {2:.6f} {3:.6f} {4:.6f}\n"
.format(node, node_dof, vec.x, vec.y, vec.z, node)
)
# elements
for element in femelement_table:
# z88_element_type is checked for every element, but mixed elements are not supported up to date
# z88_element_type is checked for every element
# but mixed elements are not supported up to date
n = femelement_table[element]
if z88_element_type == 2 or z88_element_type == 4 or z88_element_type == 5 or z88_element_type == 9 or z88_element_type == 13 or z88_element_type == 25:
if z88_element_type == 2 \
or z88_element_type == 4 \
or z88_element_type == 5 \
or z88_element_type == 9 \
or z88_element_type == 13 \
or z88_element_type == 25:
# seg2 FreeCAD --> stab4 Z88
# N1, N2
f.write("{0} {1}\n".format(element, z88_element_type, element))
@@ -419,8 +435,12 @@ def write_z88_mesh_to_file(femnodes_mesh, femelement_table, z88_element_type, f)
# 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
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".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]))
f.write(
"{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]
)
)
else:
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)

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

@@ -99,7 +99,12 @@ def write_fenics_mesh_points_xdmf(fem_mesh_obj, geometrynode, encoding=ENCODING_
recalc_nodes_ind_dict = {}
if encoding == ENCODING_ASCII:
dataitem = ET.SubElement(geometrynode, "DataItem", Dimensions="%d %d" % (numnodes, effective_dim), Format="XML")
dataitem = ET.SubElement(
geometrynode,
"DataItem",
Dimensions="%d %d" % (numnodes, effective_dim),
Format="XML"
)
nodes = []
for (ind, (key, node)) in enumerate(list(fem_mesh_obj.FemMesh.Nodes.items())):
nodes.append(node)
@@ -122,7 +127,9 @@ def write_fenics_mesh_codim_xdmf(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)
@@ -148,10 +155,17 @@ 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(topologynode, "DataItem", NumberType="UInt", Dimensions="%d %d" % (num_topo, nodes_per_element), Format="XML")
dataitem = ET.SubElement(
topologynode, "DataItem",
NumberType="UInt",
Dimensions="%d %d" % (num_topo, nodes_per_element),
Format="XML"
)
dataitem.text = numpy_array_to_str(tuples_to_numpy(nodeindices, nodes_per_element))
elif encoding == ENCODING_HDF5:
pass
@@ -167,7 +181,11 @@ def write_fenics_mesh_scalar_cellfunctions(name, cell_array, attributenode, enco
(num_cells, num_dims) = np.shape(cell_array)
if encoding == ENCODING_ASCII:
dataitem = ET.SubElement(attributenode, "DataItem", Dimensions="%d %d" % (num_cells, num_dims), Format="XML")
dataitem = ET.SubElement(
attributenode, "DataItem",
Dimensions="%d %d" % (num_cells, num_dims),
Format="XML"
)
dataitem.text = numpy_array_to_str(cell_array)
elif encoding == ENCODING_HDF5:
pass
@@ -269,13 +287,19 @@ def write_fenics_mesh_xdmf(fem_mesh_obj, outputfile, group_values_dict={}, encod
print('group id: %d (label: %s) with element type %s and codim %d'
% (g, mesh_function_name, mesh_function_type, mesh_function_codim))
mesh_function_grid = ET.SubElement(domain, "Grid", Name=mesh_function_name + "_mesh", GridType="Uniform")
mesh_function_grid = ET.SubElement(
domain, "Grid",
Name=mesh_function_name + "_mesh",
GridType="Uniform"
)
mesh_function_topology = ET.SubElement(mesh_function_grid, "Topology")
mesh_function_topology_description = write_fenics_mesh_codim_xdmf(fem_mesh_obj,
mesh_function_topology,
nodes_dict,
codim=mesh_function_codim, encoding=encoding)
mesh_function_topology_description = write_fenics_mesh_codim_xdmf(
fem_mesh_obj,
mesh_function_topology,
nodes_dict,
codim=mesh_function_codim, encoding=encoding
)
mesh_function_geometry = ET.SubElement(mesh_function_grid, "Geometry", Reference="XML")
mesh_function_geometry.text = "/Xdmf/Domain/Grid/Geometry"

90
src/Mod/Fem/femmesh/gmshtools.py
View File

@@ -139,7 +139,10 @@ class GmshTools():
def create_mesh(self):
print("\nWe are going to start Gmsh FEM mesh run!")
print(' Part to mesh: Name --> ' + self.part_obj.Name + ', Label --> ' + self.part_obj.Label + ', ShapeType --> ' + self.part_obj.Shape.ShapeType)
print(
' Part to mesh: Name --> {}, Label --> {}, ShapeType --> {}'
.format(self.part_obj.Name, self.part_obj.Label, self.part_obj.Shape.ShapeType)
)
print(' CharacteristicLengthMax: ' + str(self.clmax))
print(' CharacteristicLengthMin: ' + str(self.clmin))
print(' ElementOrder: ' + self.order)
@@ -226,12 +229,18 @@ class GmshTools():
output = output.decode('utf-8')
gmsh_path = output.split('\n')[0]
elif p1.wait() == 1:
error_message = "Gmsh binary gmsh not found in standard system binary path. Please install Gmsh or set path to binary in FEM preferences tab Gmsh.\n"
error_message = (
"Gmsh binary gmsh not found in standard system binary path. "
"Please install Gmsh or set path to binary in FEM preferences tab Gmsh.\n"
)
FreeCAD.Console.PrintError(error_message)
raise Exception(error_message)
self.gmsh_bin = gmsh_path
else:
error_message = "No standard location implemented for your operating system. Set GMHS binary path in FEM preferences.\n"
error_message = (
"No standard location implemented for your operating system. "
"Set GMHS binary path in FEM preferences.\n"
)
FreeCAD.Console.PrintError(error_message)
raise Exception(error_message)
else:
@@ -288,17 +297,21 @@ class GmshTools():
pass
else:
print(' Mesh regions, we need to get the elements.')
# by the use of MeshRegion object and a BooleanSplitCompound there could be problems with node numbers see
# by the use of MeshRegion object and a BooleanSplitCompound
# there could be problems with node numbers see
# http://forum.freecadweb.org/viewtopic.php?f=18&t=18780&start=40#p149467
# http://forum.freecadweb.org/viewtopic.php?f=18&t=18780&p=149520#p149520
part = self.part_obj
if self.mesh_obj.MeshRegionList:
if part.Shape.ShapeType == "Compound" and hasattr(part, "Proxy"): # other part obj might not have a Proxy, thus an exception would be raised
# other part obj might not have a Proxy, thus an exception would be raised
if part.Shape.ShapeType == "Compound" and hasattr(part, "Proxy"):
if (part.Proxy.Type == "FeatureBooleanFragments" or part.Proxy.Type == "FeatureSlice" or part.Proxy.Type == "FeatureXOR"):
error_message = (
' The mesh to shape is a boolean split tools Compound and the mesh has mesh region list. '
' The mesh to shape is a boolean split tools Compound '
'and the mesh has mesh region list. '
'Gmsh could return unexpected meshes in such circumstances. '
'It is strongly recommended to extract the shape to mesh from the Compound and use this one.'
'It is strongly recommended to extract the shape to mesh '
'from the Compound and use this one.'
)
FreeCAD.Console.PrintError(error_message + "\n")
# TODO: no gui popup because FreeCAD will be in a endless output loop
@@ -312,7 +325,8 @@ class GmshTools():
if mr_obj.References:
for sub in mr_obj.References:
# print(sub[0]) # Part the elements belongs to
# check if the shape of the mesh region is an element of the Part to mesh, if not try to find the element in the shape to mesh
# check if the shape of the mesh region is an element of the Part to mesh,
# if not try to find the element in the shape to mesh
search_ele_in_shape_to_mesh = False
if not self.part_obj.Shape.isSame(sub[0].Shape):
# print(" One element of the meshregion " + mr_obj.Name + " is not an element of the Part to mesh.")
@@ -322,21 +336,32 @@ class GmshTools():
# print(elems) # elems --> element
if search_ele_in_shape_to_mesh:
# we're going to try to find the element in the Shape to mesh and use the found element as elems
ele_shape = FemMeshTools.get_element(sub[0], elems) # the method getElement(element) does not return Solid elements
# the method getElement(element) does not return Solid elements
ele_shape = FemMeshTools.get_element(sub[0], elems)
found_element = FemMeshTools.find_element_in_shape(self.part_obj.Shape, ele_shape)
if found_element:
elems = found_element
else:
FreeCAD.Console.PrintError("One element of the meshregion " + mr_obj.Name + " could not be found in the Part to mesh. It will be ignored.\n")
FreeCAD.Console.PrintError(
"One element of the meshregion {} could not be found "
"in the Part to mesh. It will be ignored.\n"
.format(mr_obj.Name)
)
# print(elems) # element
if elems not in self.ele_length_map:
self.ele_length_map[elems] = Units.Quantity(mr_obj.CharacteristicLength).Value
else:
FreeCAD.Console.PrintError("The element " + elems + " of the meshregion " + mr_obj.Name + " has been added to another mesh region.\n")
FreeCAD.Console.PrintError(
"The element " + elems + " of the meshregion " + mr_obj.Name + " has been added to another mesh region.\n"
)
else:
FreeCAD.Console.PrintError("The meshregion: " + mr_obj.Name + " is not used to create the mesh because the reference list is empty.\n")
FreeCAD.Console.PrintError(
"The meshregion: " + mr_obj.Name + " is not used to create the mesh because the reference list is empty.\n"
)
else:
FreeCAD.Console.PrintError("The meshregion: " + mr_obj.Name + " is not used to create the mesh because the CharacteristicLength is 0.0 mm.\n")
FreeCAD.Console.PrintError(
"The meshregion: " + mr_obj.Name + " is not used to create the mesh because the CharacteristicLength is 0.0 mm.\n"
)
for eleml in self.ele_length_map:
ele_shape = FemMeshTools.get_element(self.part_obj, eleml) # the method getElement(element) does not return Solid elements
ele_vertexes = FemMeshTools.get_vertexes_by_element(self.part_obj.Shape, ele_shape)
@@ -354,8 +379,13 @@ class GmshTools():
else:
print(' Mesh boundary layers, we need to get the elements.')
if self.part_obj.Shape.ShapeType == 'Compound':
# see http://forum.freecadweb.org/viewtopic.php?f=18&t=18780&start=40#p149467 and http://forum.freecadweb.org/viewtopic.php?f=18&t=18780&p=149520#p149520
err = "Gmsh could return unexpected meshes for a boolean split tools Compound. It is strongly recommended to extract the shape to mesh from the Compound and use this one."
# see http://forum.freecadweb.org/viewtopic.php?f=18&t=18780&start=40#p149467 and
# http://forum.freecadweb.org/viewtopic.php?f=18&t=18780&p=149520#p149520
err = (
"Gmsh could return unexpected meshes for a boolean split tools Compound. "
"It is strongly recommended to extract the shape to mesh "
"from the Compound and use this one."
)
FreeCAD.Console.PrintError(err + "\n")
for mr_obj in self.mesh_obj.MeshBoundaryLayerList:
if mr_obj.MinimumThickness and Units.Quantity(mr_obj.MinimumThickness).Value > 0:
@@ -363,7 +393,9 @@ class GmshTools():
belem_list = []
for sub in mr_obj.References:
# print(sub[0]) # Part the elements belongs to
# check if the shape of the mesh boundary_layer is an element of the Part to mesh, if not try to find the element in the shape to mesh
# check if the shape of the mesh boundary_layer is an
# element of the Part to mesh
# if not try to find the element in the shape to mesh
search_ele_in_shape_to_mesh = False
if not self.part_obj.Shape.isSame(sub[0].Shape):
# print(" One element of the mesh boundary layer " + mr_obj.Name + " is not an element of the Part to mesh.")
@@ -373,19 +405,28 @@ class GmshTools():
# print(elems) # elems --> element
if search_ele_in_shape_to_mesh:
# we try to find the element it in the Shape to mesh and use the found element as elems
ele_shape = FemMeshTools.get_element(sub[0], elems) # the method getElement(element) does not return Solid elements
# the method getElement(element) does not return Solid elements
ele_shape = FemMeshTools.get_element(sub[0], elems)
found_element = FemMeshTools.find_element_in_shape(self.part_obj.Shape, ele_shape)
if found_element: # also
elems = found_element
else:
FreeCAD.Console.PrintError("One element of the mesh boundary layer " + mr_obj.Name + " could not be found in the Part to mesh. It will be ignored.\n")
FreeCAD.Console.PrintError(
"One element of the mesh boundary layer {} could not be found "
"in the Part to mesh. It will be ignored.\n"
.format(mr_obj.Name)
)
# print(elems) # element
if elems not in self.bl_boundary_list:
# fetch settings in DocumentObject, fan setting is not implemented
belem_list.append(elems)
self.bl_boundary_list.append(elems)
else:
FreeCAD.Console.PrintError("The element " + elems + " of the mesh boundary layer " + mr_obj.Name + " has been added to another mesh boundary layer.\n")
FreeCAD.Console.PrintError(
"The element {} of the mesh boundary layer {} has been added "
"to another mesh boundary layer.\n"
.format(elems, mr_obj.Name)
)
setting = {}
setting['hwall_n'] = Units.Quantity(mr_obj.MinimumThickness).Value
setting['ratio'] = mr_obj.GrowthRate
@@ -407,9 +448,13 @@ class GmshTools():
FreeCAD.Console.PrintError("boundary layer is only supported for 2D and 3D mesh")
self.bl_setting_list.append(setting)
else:
FreeCAD.Console.PrintError("The mesh boundary layer: " + mr_obj.Name + " is not used to create the mesh because the reference list is empty.\n")
FreeCAD.Console.PrintError(
"The mesh boundary layer: " + mr_obj.Name + " is not used to create the mesh because the reference list is empty.\n"
)
else:
FreeCAD.Console.PrintError("The mesh boundary layer: " + mr_obj.Name + " is not used to create the mesh because the min thickness is 0.0 mm.\n")
FreeCAD.Console.PrintError(
"The mesh boundary layer: " + mr_obj.Name + " is not used to create the mesh because the min thickness is 0.0 mm.\n"
)
print(' {}'.format(self.bl_setting_list))
def write_boundary_layer(self, geo):
@@ -553,7 +598,8 @@ class GmshTools():
if self.group_elements and self.group_nodes_export:
geo.write("// For each group save not only the elements but the nodes too.;\n")
geo.write("Mesh.SaveGroupsOfNodes = 1;\n")
geo.write("// Needed for Group meshing too, because for one material there is no group defined;\n") # belongs to Mesh.SaveAll but only needed if there are groups
# belongs to Mesh.SaveAll but only needed if there are groups
geo.write("// Needed for Group meshing too, because for one material there is no group defined;\n")
geo.write("// Ignore Physical definitions and save all elements;\n")
geo.write("Mesh.SaveAll = 1;\n")
geo.write('Save "' + self.temp_file_mesh + '";\n')

246
src/Mod/Fem/femmesh/meshtools.py
View File

@@ -36,9 +36,15 @@ def get_femnodes_by_femobj_with_references(femmesh, femobj):
node_set = get_femmesh_groupdata_sets_by_name(femmesh, femobj, 'Node')
# FreeCAD.Console.PrintMessage('node_set_group: {}\n'.format(node_set))
if node_set:
FreeCAD.Console.PrintMessage(" Finite element mesh nodes where retrieved from existent finite element mesh group data.\n")
FreeCAD.Console.PrintMessage(
" Finite element mesh nodes where retrieved "
"from existent finite element mesh group data.\n"
)
if not node_set:
FreeCAD.Console.PrintMessage(" Finite element mesh nodes will be retrieved by searching the appropriate nodes in the finite element mesh.\n")
FreeCAD.Console.PrintMessage(
" Finite element mesh nodes will be retrieved "
"by searching the appropriate nodes in the finite element mesh.\n"
)
node_set = get_femnodes_by_references(femmesh, femobj['Object'].References)
# FreeCAD.Console.PrintMessage('node_set_nogroup: {}\n'.format(node_set))
@@ -50,13 +56,23 @@ def get_femelements_by_references(femmesh, femelement_table, references, femnode
'''
references_femelements = []
for ref in references:
ref_femnodes = get_femnodes_by_refshape(femmesh, ref) # femnodes for the current ref
# femnodes for the current ref
ref_femnodes = get_femnodes_by_refshape(femmesh, ref)
if femnodes_ele_table:
# blind fast binary search, works for volumes only
references_femelements += get_femelements_by_femnodes_bin(femelement_table, femnodes_ele_table, ref_femnodes) # femelements for all references
# femelements for all references
references_femelements += get_femelements_by_femnodes_bin(
femelement_table,
femnodes_ele_table,
ref_femnodes
)
else:
# standard search
references_femelements += get_femelements_by_femnodes_std(femelement_table, ref_femnodes) # femelements for all references
# femelements for all references
references_femelements += get_femelements_by_femnodes_std(
femelement_table,
ref_femnodes
)
return references_femelements
@@ -189,7 +205,8 @@ def get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_set):
http://forum.freecadweb.org/viewtopic.php?f=18&p=141133&sid=013c93f496a63872951d2ce521702ffa#p141108
The bit_pattern_dict holds later an integer (bit array) for each element, which gives us
the information we are searching for:
Is this element part of the node list (searching for elements) or has this element a face we are searching for?
Is this element part of the node list (searching for elements)
or has this element a face we are searching for?
The number in the ele_dict is organized as a bit array.
The corresponding bit is set, if the node of the node_set is contained in the element.
'''
@@ -363,23 +380,32 @@ def get_femvolumeelements_by_femfacenodes(femelement_table, node_list):
if nodecount == 6 or nodecount == 8:
e.append(elementID)
else:
FreeCAD.Console.PrintError('Error in get_femvolumeelements_by_femfacenodes(): unknown volume element: ' + el_nd_ct + '\n')
FreeCAD.Console.PrintError(
'Error in get_femvolumeelements_by_femfacenodes(): unknown volume element: ' + el_nd_ct + '\n'
)
# FreeCAD.Console.PrintMessage('{}\n'.format(sorted(e)))
return e
def get_femelement_sets(femmesh, femelement_table, fem_objects, femnodes_ele_table=None): # fem_objects = FreeCAD FEM document objects
def get_femelement_sets(femmesh, femelement_table, fem_objects, femnodes_ele_table=None):
# fem_objects = FreeCAD FEM document objects
# get femelements for reference shapes of each obj.References
count_femelements = 0
referenced_femelements = []
has_remaining_femelements = None
for fem_object_i, fem_object in enumerate(fem_objects):
obj = fem_object['Object']
FreeCAD.Console.PrintMessage("Constraint: " + obj.Name + " --> " + "We're going to search in the mesh for the element ID's.\n")
FreeCAD.Console.PrintMessage(
"Constraint: " + obj.Name + " --> " + "We're going to search in the mesh for the element ID's.\n"
)
fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) # unique short identifier
if obj.References:
ref_shape_femelements = []
ref_shape_femelements = get_femelements_by_references(femmesh, femelement_table, obj.References, femnodes_ele_table)
ref_shape_femelements = get_femelements_by_references(
femmesh, femelement_table,
obj.References,
femnodes_ele_table
)
referenced_femelements += ref_shape_femelements
count_femelements += len(ref_shape_femelements)
fem_object['FEMElements'] = ref_shape_femelements
@@ -406,7 +432,8 @@ def get_femelement_sets(femmesh, femelement_table, fem_objects, femnodes_ele_tab
def get_femelement_direction1D_set(femmesh, femelement_table, beamrotation_objects, theshape=None):
'''
get for each geometry edge direction, the normal and the element ids and write all into the beamrotation_objects
get for each geometry edge direction, the normal and the element ids and
# write all into the beamrotation_objects
means no return value, we're going to write into the beamrotation_objects dictionary
FEMRotations1D is a list of dictionaries for every beamdirection of all edges
beamrot_obj['FEMRotations1D'] = [ {'ids' : [theids],
@@ -416,16 +443,19 @@ def get_femelement_direction1D_set(femmesh, femelement_table, beamrotation_objec
]
'''
if len(beamrotation_objects) == 0:
# no beamrotation document object, all beams use standard rotation of 0 degree (angle), we need theshape (the shape which was meshed)
# no beamrotation document object, all beams use standard rotation of 0 degree (angle)
# we need theshape (the shape which was meshed)
# since ccx needs to split them in sets anyway we need to take care of this too
rotations_ids = get_femelement_directions_theshape(femmesh, femelement_table, theshape)
# add normals for each direction
rotation_angle = 0
for rot in rotations_ids:
rot['normal'] = get_beam_normal(rot['direction'], rotation_angle)
beamrotation_objects.append({'FEMRotations1D': rotations_ids, 'ShortName': 'Rstd'}) # key 'Object' will be empty
# key 'Object' will be empty
beamrotation_objects.append({'FEMRotations1D': rotations_ids, 'ShortName': 'Rstd'})
elif len(beamrotation_objects) == 1:
# one beamrotation document object with no references, all beams use rotation from this object, we need theshape (the shape which was meshed)
# one beamrotation document object with no references, all beams use rotation from this object
# we need theshape (the shape which was meshed)
# since ccx needs to split them in sets anyway we need to take care of this too
rotations_ids = get_femelement_directions_theshape(femmesh, femelement_table, theshape)
# add normals for each direction
@@ -457,9 +487,12 @@ def get_femelement_directions_theshape(femmesh, femelement_table, theshape):
the_edge = {}
the_edge['direction'] = e.Vertexes[1].Point - e.Vertexes[0].Point
edge_femnodes = femmesh.getNodesByEdge(e) # femnodes for the current edge
the_edge['ids'] = get_femelements_by_femnodes_std(femelement_table, edge_femnodes) # femelements for this edge
# femelements for this edge
the_edge['ids'] = get_femelements_by_femnodes_std(femelement_table, edge_femnodes)
for rot in rotations_ids:
if rot['direction'] == the_edge['direction']: # tolerance will be managed by FreeCAD see https://forum.freecadweb.org/viewtopic.php?f=22&t=14179
# tolerance will be managed by FreeCAD
# see https://forum.freecadweb.org/viewtopic.php?f=22&t=14179
if rot['direction'] == the_edge['direction']:
rot['ids'] += the_edge['ids']
break
else:
@@ -536,7 +569,9 @@ def get_femmesh_groupdata_sets_by_name(femmesh, fem_object, group_data_type):
grp_name = femmesh.getGroupName(g)
if grp_name.startswith(obj.Name + "_"):
if femmesh.getGroupElementType(g) == group_data_type:
FreeCAD.Console.PrintMessage(" found mesh group for the IDs: " + grp_name + ', Type: ' + group_data_type + '\n')
FreeCAD.Console.PrintMessage(
" found mesh group for the IDs: " + grp_name + ', Type: ' + group_data_type + '\n'
)
return femmesh.getGroupElements(g) # == ref_shape_femelements
return () # an empty tuple is returned if no group data IDs where found
@@ -547,15 +582,20 @@ def get_femelement_sets_from_group_data(femmesh, fem_objects):
sum_group_elements = []
for fem_object_i, fem_object in enumerate(fem_objects):
obj = fem_object['Object']
FreeCAD.Console.PrintMessage("Constraint: " + obj.Name + " --> " + "We have mesh groups. We will search for appropriate group data.\n")
FreeCAD.Console.PrintMessage(
"Constraint: " + obj.Name + " --> " + "We have mesh groups. We will search for appropriate group data.\n"
)
fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) # unique short identifier
group_elements = get_femmesh_groupdata_sets_by_name(femmesh, fem_object, 'Volume') # see comments over there !
# see comments over there !
group_elements = get_femmesh_groupdata_sets_by_name(femmesh, fem_object, 'Volume')
sum_group_elements += group_elements
count_femelements += len(group_elements)
fem_object['FEMElements'] = group_elements
# check if all worked out well
if not femelements_count_ok(femmesh.VolumeCount, count_femelements):
FreeCAD.Console.PrintError('Error in get_femelement_sets_from_group_data -- > femelements_count_ok() failed!\n')
FreeCAD.Console.PrintError(
'Error in get_femelement_sets_from_group_data -- > femelements_count_ok() failed!\n'
)
return False
else:
return True
@@ -573,7 +613,9 @@ def get_elset_short_name(obj, i):
elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'Fem::FemElementGeometry2D':
return 'S' + str(i)
else:
FreeCAD.Console.PrintError('Error in creating short elset name for obj: ' + obj.Name + ' --> Proxy.Type: ' + str(obj.Proxy.Type) + '\n')
FreeCAD.Console.PrintError(
'Error in creating short elset name for obj: ' + obj.Name + ' --> Proxy.Type: ' + str(obj.Proxy.Type) + '\n'
)
# ************************************************************************************************
@@ -736,8 +778,10 @@ def get_ref_edgenodes_table(femmesh, femelement_table, refedge):
for node in femelement_table[elem]:
if node in refedge_nodes:
fe_refedge_nodes.append(node)
edge_table[elem] = fe_refedge_nodes # { volumeID : ( edgenodeID, ... , edgenodeID )} # only the refedge nodes
# FIXME: duplicate_mesh_elements: as soon as contact and springs are supported the user should decide on which edge the load is applied
# { volumeID : ( edgenodeID, ... , edgenodeID )} # only the refedge nodes
edge_table[elem] = fe_refedge_nodes
# FIXME: duplicate_mesh_elements: as soon as contact and springs are supported
# the user should decide on which edge the load is applied
edge_table = delete_duplicate_mesh_elements(edge_table)
elif is_face_femmesh(femmesh):
refedge_fem_faceelements = []
@@ -755,13 +799,16 @@ def get_ref_edgenodes_table(femmesh, femelement_table, refedge):
for node in femelement_table[elem]:
if node in refedge_nodes:
fe_refedge_nodes.append(node)
edge_table[elem] = fe_refedge_nodes # { faceID : ( edgenodeID, ... , edgenodeID )} # only the refedge nodes
# FIXME: duplicate_mesh_elements: as soon as contact and springs are supported the user should decide on which edge the load is applied
# { faceID : ( edgenodeID, ... , edgenodeID )} # only the refedge nodes
edge_table[elem] = fe_refedge_nodes
# FIXME: duplicate_mesh_elements: as soon as contact and springs are supported
# the user should decide on which edge the load is applied
edge_table = delete_duplicate_mesh_elements(edge_table)
elif is_edge_femmesh(femmesh):
refedge_fem_edgeelements = get_femelements_by_femnodes_std(femelement_table, refedge_nodes)
for elem in refedge_fem_edgeelements:
edge_table[elem] = femelement_table[elem] # { edgeID : ( nodeID, ... , nodeID )} # all nodes off this femedgeelement
# { edgeID : ( nodeID, ... , nodeID )} # all nodes off this femedgeelement
edge_table[elem] = femelement_table[elem]
return edge_table
@@ -891,32 +938,47 @@ def get_ref_facenodes_table(femmesh, femelement_table, ref_face):
face_table = {} # { meshfaceID : ( nodeID, ... , nodeID ) }
if is_solid_femmesh(femmesh):
if has_no_face_data(femmesh):
FreeCAD.Console.PrintMessage(' No face data in finite volume element mesh. FreeCAD uses getccxVolumesByFace() to retrieve the volume elements of the ref_face.\n')
FreeCAD.Console.PrintMessage(
' No face data in finite volume element mesh. '
'FreeCAD uses getccxVolumesByFace() '
'to retrieve the volume elements of the ref_face.\n'
)
# there is no face data
# the problem if we retrieve the nodes ourself is they are not sorted we just have the nodes.
# We need to sort them according the shell mesh notation of tria3, tria6, quad4, quad8
ref_face_nodes = femmesh.getNodesByFace(ref_face)
# try to use getccxVolumesByFace() to get the volume ids of element with elementfaces on the ref_face --> should work for tetra4 and tetra10
# try to use getccxVolumesByFace() to get the volume ids
# of element with elementfaces on the ref_face --> should work for tetra4 and tetra10
ref_face_volume_elements = femmesh.getccxVolumesByFace(ref_face) # list of tupels (mv, ccx_face_nr)
if ref_face_volume_elements: # mesh with tetras
FreeCAD.Console.PrintLog(' Use of getccxVolumesByFace() has returned volume elements of the ref_face.\n')
FreeCAD.Console.PrintLog(
' Use of getccxVolumesByFace() has returned volume elements of the ref_face.\n'
)
for ve in ref_face_volume_elements:
veID = ve[0]
ve_ref_face_nodes = []
for nodeID in femelement_table[veID]:
if nodeID in ref_face_nodes:
ve_ref_face_nodes.append(nodeID)
face_table[veID] = ve_ref_face_nodes # { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
# { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
face_table[veID] = ve_ref_face_nodes
else: # mesh with hexa or penta
FreeCAD.Console.PrintLog(' The use of getccxVolumesByFace() has NOT returned volume elements of the ref_face. FreeCAD tries to use get_femvolumeelements_by_femfacenodes().\n')
ref_face_volume_elements = get_femvolumeelements_by_femfacenodes(femelement_table, ref_face_nodes) # list of integer [mv]
FreeCAD.Console.PrintLog(
' The use of getccxVolumesByFace() has NOT returned '
'volume elements of the ref_face. '
'FreeCAD tries to use get_femvolumeelements_by_femfacenodes().\n'
)
# list of integer [mv]
ref_face_volume_elements = get_femvolumeelements_by_femfacenodes(femelement_table, ref_face_nodes)
for veID in ref_face_volume_elements:
ve_ref_face_nodes = []
for nodeID in femelement_table[veID]:
if nodeID in ref_face_nodes:
ve_ref_face_nodes.append(nodeID)
face_table[veID] = ve_ref_face_nodes # { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
face_table = build_mesh_faces_of_volume_elements(face_table, femelement_table) # we need to resort the nodes to make them build an element face
# { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
face_table[veID] = ve_ref_face_nodes
# we need to resort the nodes to make them build an element face
face_table = build_mesh_faces_of_volume_elements(face_table, femelement_table)
else: # the femmesh has face_data
faces = femmesh.getFacesByFace(ref_face) # (mv, mf)
for mf in faces:
@@ -1094,7 +1156,9 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table):
elif face_node_indexs == [2, 3, 4, 6, 9, 10]:
node_numbers = (2, 4, 3, 9, 10, 6)
else:
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
FreeCAD.Console.PrintError(
"Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n"
)
elif vol_node_ct == 4: # tetra4 --> tria3 face
if face_node_indexs == [1, 2, 3]: # node order of face in tetra4 volume element
node_numbers = (1, 2, 3) # node order of a tria3 face of tetra4
@@ -1105,7 +1169,9 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table):
elif face_node_indexs == [2, 3, 4]:
node_numbers = (2, 4, 3)
else:
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
FreeCAD.Console.PrintError(
"Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n"
)
elif vol_node_ct == 20: # hexa20 --> quad8 face
if face_node_indexs == [1, 2, 3, 4, 9, 10, 11, 12]: # node order of face in hexa20 volume element
node_numbers = (1, 2, 3, 4, 9, 10, 11, 12) # node order of a quad8 face of hexa20
@@ -1120,7 +1186,9 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table):
elif face_node_indexs == [2, 3, 6, 7, 10, 14, 18, 19]:
node_numbers = (2, 6, 7, 3, 18, 14, 19, 10)
else:
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
FreeCAD.Console.PrintError(
"Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n"
)
elif vol_node_ct == 8: # hexa8 --> quad4 face
face_node_indexs = sorted(face_nodenumber_table[veID])
if face_node_indexs == [1, 2, 3, 4]: # node order of face in hexa8 volume element
@@ -1136,7 +1204,9 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table):
elif face_node_indexs == [2, 3, 6, 7]:
node_numbers = (2, 6, 7, 3)
else:
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
FreeCAD.Console.PrintError(
"Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n"
)
elif vol_node_ct == 15: # penta15 --> tria6 and quad8 faces
if face_node_indexs == [1, 2, 3, 7, 8, 9]: # node order of face in penta15 volume element
node_numbers = (1, 2, 3, 7, 8, 9) # node order of a tria6 face of penta15
@@ -1149,7 +1219,9 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table):
elif face_node_indexs == [2, 3, 5, 6, 8, 11, 14, 15]:
node_numbers = (2, 5, 6, 3, 14, 11, 15, 8) # quad8
else:
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): penta15: face not found!" + str(face_node_indexs) + "\n")
FreeCAD.Console.PrintError(
"Error in build_mesh_faces_of_volume_elements(): penta15: face not found!" + str(face_node_indexs) + "\n"
)
elif vol_node_ct == 6: # penta6 --> tria3 and quad4 faces
if face_node_indexs == [1, 2, 3]: # node order of face in penta6 volume element
node_numbers = (1, 2, 3) # node order of a tria3 face of penta6
@@ -1162,9 +1234,13 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table):
elif face_node_indexs == [2, 3, 5, 6]:
node_numbers = (2, 5, 6, 3) # quad4
else:
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): pent6: face not found!" + str(face_node_indexs) + "\n")
FreeCAD.Console.PrintError(
"Error in build_mesh_faces_of_volume_elements(): pent6: face not found!" + str(face_node_indexs) + "\n"
)
else:
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): Volume not implemented: volume node count" + str(vol_node_ct) + "\n")
FreeCAD.Console.PrintError(
"Error in build_mesh_faces_of_volume_elements(): Volume not implemented: volume node count" + str(vol_node_ct) + "\n"
)
face_nodes = []
for i in node_numbers:
i -= 1 # node_number starts with 1, index starts with 0 --> index = node number - 1
@@ -1214,7 +1290,9 @@ def get_pressure_obj_faces(femmesh, femelement_table, femnodes_ele_table, femobj
# -1 if femmeshface normal opposite reference face normal direction
# easy on plane faces, but on a half sphere ... ?!?
else:
FreeCAD.Console.PrintError("Pressure on shell mesh at the moment only supported for meshes with appropriate group data.\n")
FreeCAD.Console.PrintError(
"Pressure on shell mesh at the moment only supported for meshes with appropriate group data.\n"
)
return pressure_faces
@@ -1246,7 +1324,9 @@ def get_mesh_group_elements(mesh_group_obj, aPart):
grp_ele = get_reference_group_elements(mesh_group_obj, aPart)
group_elements[grp_ele[0]] = grp_ele[1]
else:
FreeCAD.Console.PrintError(' Empty reference in mesh group object: ' + mesh_group_obj.Name + ' ' + mesh_group_obj.Label + '\n')
FreeCAD.Console.PrintError(
' Empty reference in mesh group object: ' + mesh_group_obj.Name + ' ' + mesh_group_obj.Label + '\n'
)
return group_elements
@@ -1282,7 +1362,9 @@ def get_analysis_group_elements(aAnalysis, aPart):
for g in group_elements:
# FreeCAD.Console.PrintMessage('{}\n'.format(group_elements[g]))
if len(group_elements[g]) == 0:
FreeCAD.Console.PrintError('Error: The shapes for the mesh group for the reference shapes of analysis member: ' + g + ' could not be found!\n')
FreeCAD.Console.PrintError(
'Error: The shapes for the mesh group for the reference shapes of analysis member: ' + g + ' could not be found!\n'
)
return group_elements
@@ -1320,14 +1402,22 @@ def get_reference_group_elements(obj, aPart):
if found_element is not None:
elements.append(found_element)
else:
FreeCAD.Console.PrintError('Problem: For the geometry of the following shape was no Shape found: ' + str(ref_shape) + '\n')
FreeCAD.Console.PrintError(
'Problem: For the geometry of the following shape was no Shape found: ' + str(ref_shape) + '\n'
)
FreeCAD.Console.PrintMessage(' ' + obj.Name + '\n')
FreeCAD.Console.PrintMessage(' ' + str(obj.References) + '\n')
FreeCAD.Console.PrintMessage(' ' + r[0].Name + '\n')
if parent.Name != aPart.Name:
FreeCAD.Console.PrintError('The reference Shape is not a child nor it is the shape the mesh is made of. : ' + str(ref_shape) + '\n')
FreeCAD.Console.PrintError(
'The reference Shape is not a child nor it is the shape the mesh is made of. : ' + str(ref_shape) + '\n'
)
FreeCAD.Console.PrintMessage(aPart.Name + '--> Name of the Feature we where searching in.\n')
FreeCAD.Console.PrintMessage(parent.Name + '--> Name of the parent Feature of reference Shape (Use the same as in the line before and you will have less trouble :-) !!!!!!).\n')
FreeCAD.Console.PrintMessage(
'{} --> Name of the parent Feature of reference Shape '
'(Use the same as in the line before and you will have less trouble :-) !!!!!!).\n'
.format(parent.Name)
)
# import Part
# Part.show(aShape)
# Part.show(ref_shape)
@@ -1356,7 +1446,10 @@ def get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShap
if not empty_reference_material:
empty_reference_material = m.Name
else:
FreeCAD.Console.PrintError('Problem in get_anlysis_empty_references_group_elements, we seem to have two or more materials with empty references\n')
FreeCAD.Console.PrintError(
'Problem in get_anlysis_empty_references_group_elements, '
'we seem to have two or more materials with empty references\n'
)
return {}
elif hasattr(m, "References") and m.References:
# ShapeType of the group elements, strip the number of the first group element
@@ -1365,7 +1458,9 @@ def get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShap
if not material_shape_type:
material_shape_type = group_shape_type
elif material_shape_type != group_shape_type:
FreeCAD.Console.PrintError('Problem, material shape type does not match get_anlysis_empty_references_group_elements\n')
FreeCAD.Console.PrintError(
'Problem, material shape type does not match get_anlysis_empty_references_group_elements\n'
)
for ele in group_elements[m.Name]:
material_ref_shapes.append(ele)
if material_shape_type == 'Solid':
@@ -1387,7 +1482,9 @@ def get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShap
if ele not in material_ref_shapes:
missed_material_refshapes.append(ele)
else:
FreeCAD.Console.PrintMessage(' One material with no reference shapes. No need to make a group for materials.\n')
FreeCAD.Console.PrintMessage(
' One material with no reference shapes. No need to make a group for materials.\n'
)
# make no changes group_elements
return group_elements
# FreeCAD.Console.PrintMessage('{}\n'.format(sorted(material_ref_shapes)))
@@ -1409,11 +1506,14 @@ def find_element_in_shape(aShape, anElement):
# Part.show(aShape.Solids[index])
ele = ele_st + str(index + 1)
return ele
FreeCAD.Console.PrintError('Solid ' + str(anElement) + ' not found in: ' + str(aShape) + '\n')
FreeCAD.Console.PrintError(
'Solid ' + str(anElement) + ' not found in: ' + str(aShape) + '\n'
)
if ele_st == 'Solid' and aShape.ShapeType == 'Solid':
message_part = (
'We have been searching for a Solid in a Solid and we have not found it. '
'In most cases this should be searching for a Solid inside a CompSolid. Check the ShapeType of your Part to mesh.'
'In most cases this should be searching for a Solid inside a CompSolid. '
'Check the ShapeType of your Part to mesh.'
)
FreeCAD.Console.PrintMessage(message_part + '\n')
# Part.show(anElement)
@@ -1460,7 +1560,9 @@ def get_vertexes_by_element(aShape, anElement):
ele_vertexes.append(i)
# FreeCAD.Console.PrintMessage(' ' + str(sorted(ele_vertexes)), '\n')
return ele_vertexes
FreeCAD.Console.PrintError('Error, Solid ' + str(anElement) + ' not found in: ' + str(aShape) + '\n')
FreeCAD.Console.PrintError(
'Error, Solid ' + str(anElement) + ' not found in: ' + str(aShape) + '\n'
)
elif ele_st == 'Face' or ele_st == 'Shell':
for index, face in enumerate(aShape.Faces):
if is_same_geometry(face, anElement):
@@ -1528,7 +1630,9 @@ def get_element(part, element):
if element.startswith('Solid'):
index = int(element.lstrip('Solid')) - 1
if index >= len(part.Shape.Solids):
FreeCAD.Console.PrintError('Index out of range. This Solid does not exist in the Shape!\n')
FreeCAD.Console.PrintError(
'Index out of range. This Solid does not exist in the Shape!\n'
)
return None
else:
return part.Shape.Solids[index] # Solid
@@ -1537,12 +1641,20 @@ def get_element(part, element):
def femelements_count_ok(len_femelement_table, count_femelements):
FreeCAD.Console.PrintMessage('Count finite elements as sum of constraints: {}\n'.format(count_femelements))
FreeCAD.Console.PrintMessage('Count finite elements of the finite element mesh: {}\n'.format(len_femelement_table))
FreeCAD.Console.PrintMessage(
'Count finite elements as sum of constraints: {}\n'
.format(count_femelements)
)
FreeCAD.Console.PrintMessage(
'Count finite elements of the finite element mesh: {}\n'
.format(len_femelement_table)
)
if count_femelements == len_femelement_table:
return True
else:
FreeCAD.Console.PrintMessage('ERROR: femelement_table != count_femelements\n')
FreeCAD.Console.PrintMessage(
'ERROR: femelement_table != count_femelements\n'
)
return False
@@ -1572,22 +1684,26 @@ def sortlistoflistvalues(listoflists):
def is_solid_femmesh(femmesh):
if femmesh.VolumeCount > 0: # solid femmesh
# solid femmesh
if femmesh.VolumeCount > 0:
return True
def has_no_face_data(femmesh):
if femmesh.FaceCount == 0: # femmesh has no face data, could be a edge femmesh or a solid femmesh without face data
# femmesh has no face data, could be a edge femmesh or a solid femmesh without face data
if femmesh.FaceCount == 0:
return True
def is_face_femmesh(femmesh):
if femmesh.VolumeCount == 0 and femmesh.FaceCount > 0: # face femmesh
# face femmesh
if femmesh.VolumeCount == 0 and femmesh.FaceCount > 0:
return True
def is_edge_femmesh(femmesh):
if femmesh.VolumeCount == 0 and femmesh.FaceCount == 0 and femmesh.EdgeCount > 0: # edge femmesh
# edge femmesh
if femmesh.VolumeCount == 0 and femmesh.FaceCount == 0 and femmesh.EdgeCount > 0:
return True
@@ -1771,7 +1887,9 @@ def use_correct_fluidinout_ele_def(FluidInletoutlet_ele, fileName, fluid_inout_n
node1 = int(a[j + 2])
node2 = int(a[j + 1])
node3 = int(a[j])
inout_nodes_file.write(str(node1) + ',' + str(node2) + ',' + str(node3) + ',' + FluidInletoutlet_ele[i][1] + '\n')
inout_nodes_file.write(
str(node1) + ',' + str(node2) + ',' + str(node3) + ',' + FluidInletoutlet_ele[i][1] + '\n'
)
elif j == 3:
new_line = new_line + a[j]
else:
@@ -1782,7 +1900,9 @@ def use_correct_fluidinout_ele_def(FluidInletoutlet_ele, fileName, fluid_inout_n
node1 = int(a[j - 2])
node2 = int(a[j - 1])
node3 = int(a[j])
inout_nodes_file.write(str(node1) + ',' + str(node2) + ',' + str(node3) + ',' + FluidInletoutlet_ele[i][1] + '\n')
inout_nodes_file.write(
str(node1) + ',' + str(node2) + ',' + str(node3) + ',' + FluidInletoutlet_ele[i][1] + '\n'
)
else:
new_line = new_line + a[j] + ','
if new_line == '':

200
src/Mod/Fem/femtools/ccxtools.py
View File

@@ -264,13 +264,19 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
if not self.working_dir:
message += "Working directory not set\n"
if not (os.path.isdir(self.working_dir)):
message += "Working directory \'{}\' doesn't exist.".format(self.working_dir)
message += (
"Working directory \'{}\' doesn't exist."
.format(self.working_dir)
)
# solver
if not self.solver:
message += "No solver object defined in the analysis\n"
else:
if self.solver.AnalysisType not in self.known_analysis_types:
message += "Unknown analysis type: {}\n".format(self.solver.AnalysisType)
message += (
"Unknown analysis type: {}\n"
.format(self.solver.AnalysisType)
)
if self.solver.AnalysisType == "frequency":
if not hasattr(self.solver, "EigenmodeHighLimit"):
message += "Frequency analysis: Solver has no EigenmodeHighLimit.\n"
@@ -278,22 +284,50 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
message += "Frequency analysis: Solver has no EigenmodeLowLimit.\n"
elif not hasattr(self.solver, "EigenmodesCount"):
message += "Frequency analysis: Solver has no EigenmodesCount.\n"
if hasattr(self.solver, "MaterialNonlinearity") and self.solver.MaterialNonlinearity == "nonlinear":
if hasattr(self.solver, "MaterialNonlinearity") \
and self.solver.MaterialNonlinearity == "nonlinear":
if not self.materials_nonlinear:
message += "Solver is set to nonlinear materials, but there is no nonlinear material in the analysis.\n"
if self.solver.Proxy.Type == 'Fem::FemSolverCalculixCcxTools' and self.solver.GeometricalNonlinearity != "nonlinear":
# nonlinear geometry --> should be set https://forum.freecadweb.org/viewtopic.php?f=18&t=23101&p=180489#p180489
message += "Solver CalculiX triggers nonlinear geometry for nonlinear material, thus it should to be set too.\n"
message += (
"Solver is set to nonlinear materials, "
"but there is no nonlinear material in the analysis.\n"
)
if self.solver.Proxy.Type == 'Fem::FemSolverCalculixCcxTools' \
and self.solver.GeometricalNonlinearity != "nonlinear":
# nonlinear geometry --> should be set
# https://forum.freecadweb.org/viewtopic.php?f=18&t=23101&p=180489#p180489
message += (
"Solver CalculiX triggers nonlinear geometry for nonlinear material, "
"thus it should to be set too.\n"
)
# mesh
if not self.mesh:
message += "No mesh object defined in the analysis\n"
if self.mesh:
if self.mesh.FemMesh.VolumeCount == 0 and self.mesh.FemMesh.FaceCount > 0 and not self.shell_thicknesses:
message += "FEM mesh has no volume elements, either define a shell thicknesses or provide a FEM mesh with volume elements.\n"
if self.mesh.FemMesh.VolumeCount == 0 and self.mesh.FemMesh.FaceCount == 0 and self.mesh.FemMesh.EdgeCount > 0 and not self.beam_sections and not self.fluid_sections:
message += "FEM mesh has no volume and no shell elements, either define a beam/fluid section or provide a FEM mesh with volume elements.\n"
if self.mesh.FemMesh.VolumeCount == 0 and self.mesh.FemMesh.FaceCount == 0 and self.mesh.FemMesh.EdgeCount == 0:
message += "FEM mesh has neither volume nor shell or edge elements. Provide a FEM mesh with elements!\n"
if self.mesh.FemMesh.VolumeCount == 0 \
and self.mesh.FemMesh.FaceCount > 0 \
and not self.shell_thicknesses:
message += (
"FEM mesh has no volume elements, "
"either define a shell thicknesses or "
"provide a FEM mesh with volume elements.\n"
)
if self.mesh.FemMesh.VolumeCount == 0 \
and self.mesh.FemMesh.FaceCount == 0 \
and self.mesh.FemMesh.EdgeCount > 0 \
and not self.beam_sections \
and not self.fluid_sections:
message += (
"FEM mesh has no volume and no shell elements, "
"either define a beam/fluid section or provide "
"a FEM mesh with volume elements.\n"
)
if self.mesh.FemMesh.VolumeCount == 0 \
and self.mesh.FemMesh.FaceCount == 0 \
and self.mesh.FemMesh.EdgeCount == 0:
message += (
"FEM mesh has neither volume nor shell or edge elements. "
"Provide a FEM mesh with elements!\n"
)
# material linear and nonlinear
if not self.materials_linear:
message += "No material object defined in the analysis\n"
@@ -301,7 +335,10 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
for m in self.materials_linear:
if len(m['Object'].References) == 0:
if has_no_references is True:
message += "More than one material has an empty references list (Only one empty references list is allowed!).\n"
message += (
"More than one material has an empty references list "
"(Only one empty references list is allowed!).\n"
)
has_no_references = True
mat_ref_shty = ''
for m in self.materials_linear:
@@ -309,10 +346,12 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
if not mat_ref_shty:
mat_ref_shty = ref_shty
if mat_ref_shty and ref_shty and ref_shty != mat_ref_shty:
# mat_ref_shty could be empty in one material, only the not empty ones should have the same shape type
# mat_ref_shty could be empty in one material
# only the not empty ones should have the same shape type
message += (
'Some material objects do not have the same reference shape type '
'(all material objects must have the same reference shape type, at the moment).\n'
'(all material objects must have the same reference shape type, '
'at the moment).\n'
)
for m in self.materials_linear:
mat_map = m['Object'].Material
@@ -325,7 +364,8 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
else:
message += "No YoungsModulus defined for at least one material.\n"
if 'PoissonRatio' not in mat_map:
message += "No PoissonRatio defined for at least one material.\n" # PoissonRatio is allowed to be 0.0 (in ccx), but it should be set anyway.
# PoissonRatio is allowed to be 0.0 (in ccx), but it should be set anyway.
message += "No PoissonRatio defined for at least one material.\n"
if self.solver.AnalysisType == "frequency" or self.selfweight_constraints:
if 'Density' not in mat_map:
message += "No Density defined for at least one material.\n"
@@ -334,15 +374,27 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
if not Units.Quantity(mat_map['ThermalConductivity']).Value:
message += "Value of ThermalConductivity is set to 0.0.\n"
else:
message += "Thermomechanical analysis: No ThermalConductivity defined for at least one material.\n"
message += (
"Thermomechanical analysis: No ThermalConductivity defined "
"for at least one material.\n"
)
if 'ThermalExpansionCoefficient' not in mat_map and mat_obj.Category == 'Solid':
message += "Thermomechanical analysis: No ThermalExpansionCoefficient defined for at least one material.\n" # allowed to be 0.0 (in ccx)
message += (
"Thermomechanical analysis: No ThermalExpansionCoefficient defined "
"for at least one material.\n" # allowed to be 0.0 (in ccx)
)
if 'SpecificHeat' not in mat_map:
message += "Thermomechanical analysis: No SpecificHeat defined for at least one material.\n" # allowed to be 0.0 (in ccx)
message += (
"Thermomechanical analysis: No SpecificHeat "
"defined for at least one material.\n" # allowed to be 0.0 (in ccx)
)
if len(self.materials_linear) == 1:
mobj = self.materials_linear[0]['Object']
if hasattr(mobj, 'References') and mobj.References:
FreeCAD.Console.PrintError('Only one material object, but this one has a reference shape. The reference shape will be ignored.\n')
FreeCAD.Console.PrintError(
'Only one material object, but this one has a reference shape. '
'The reference shape will be ignored.\n'
)
for m in self.materials_linear:
has_nonlinear_material = False
for nlm in self.materials_nonlinear:
@@ -355,11 +407,14 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
"Only one nonlinear material for each linear material allowed.\n"
)
# which analysis needs which constraints
# no check in the regard of loads existence (constraint force, pressure, self weight) is done
# because an analysis without loads at all is an valid analysis too
# no check in the regard of loads existence (constraint force, pressure, self weight)
# is done, because an analysis without loads at all is an valid analysis too
if self.solver.AnalysisType == "static":
if not (self.fixed_constraints or self.displacement_constraints):
message += "Static analysis: Neither constraint fixed nor constraint displacement defined.\n"
message += (
"Static analysis: Neither constraint fixed nor "
"constraint displacement defined.\n"
)
if self.solver.AnalysisType == "thermomech":
if not self.initialtemperature_constraints:
if not self.fluid_sections:
@@ -416,23 +471,38 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
if self.beam_sections:
if self.shell_thicknesses:
# this needs to be checked only once either here or in shell_thicknesses
message += "Beam sections and shell thicknesses in one analysis is not supported at the moment.\n"
message += (
"Beam sections and shell thicknesses in one analysis "
"is not supported at the moment.\n"
)
if self.fluid_sections:
# this needs to be checked only once either here or in shell_thicknesses
message += "Beam sections and fluid sections in one analysis is not supported at the moment.\n"
message += (
"Beam sections and fluid sections in one analysis "
"is not supported at the moment.\n"
)
has_no_references = False
for b in self.beam_sections:
if len(b['Object'].References) == 0:
if has_no_references is True:
message += "More than one beam section has an empty references list (Only one empty references list is allowed!).\n"
message += (
"More than one beam section has an empty references "
"list (Only one empty references list is allowed!).\n"
)
has_no_references = True
if self.mesh:
if self.mesh.FemMesh.FaceCount > 0 or self.mesh.FemMesh.VolumeCount > 0:
message += "Beam sections defined but FEM mesh has volume or shell elements.\n"
message += (
"Beam sections defined but FEM mesh has volume or shell elements.\n"
)
if self.mesh.FemMesh.EdgeCount == 0:
message += "Beam sections defined but FEM mesh has no edge elements.\n"
message += (
"Beam sections defined but FEM mesh has no edge elements.\n"
)
if len(self.beam_rotations) > 1:
message += "Multiple beam rotations in one analysis are not supported at the moment.\n"
message += (
"Multiple beam rotations in one analysis are not supported at the moment.\n"
)
# beam rotations
if self.beam_rotations and not self.beam_sections:
message += "Beam rotations in the analysis but no beam sections defined.\n"
@@ -442,7 +512,10 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
for s in self.shell_thicknesses:
if len(s['Object'].References) == 0:
if has_no_references is True:
message += "More than one shell thickness has an empty references list (Only one empty references list is allowed!).\n"
message += (
"More than one shell thickness has an empty references "
"list (Only one empty references list is allowed!).\n"
)
has_no_references = True
if self.mesh:
if self.mesh.FemMesh.VolumeCount > 0:
@@ -452,25 +525,33 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
# fluid section
if self.fluid_sections:
if not self.selfweight_constraints:
message += "A fluid network analysis requires self weight constraint to be applied"
message += (
"A fluid network analysis requires self weight constraint to be applied\n"
)
if self.solver.AnalysisType != "thermomech":
message += "A fluid network analysis can only be done in a thermomech analysis"
message += "A fluid network analysis can only be done in a thermomech analysis\n"
has_no_references = False
for f in self.fluid_sections:
if len(f['Object'].References) == 0:
if has_no_references is True:
message += "More than one fluid section has an empty references list (Only one empty references list is allowed!).\n"
message += (
"More than one fluid section has an empty references list "
"(Only one empty references list is allowed!).\n"
)
has_no_references = True
if self.mesh:
if self.mesh.FemMesh.FaceCount > 0 or self.mesh.FemMesh.VolumeCount > 0:
message += "Fluid sections defined but FEM mesh has volume or shell elements.\n"
message += (
"Fluid sections defined but FEM mesh has volume or shell elements.\n"
)
if self.mesh.FemMesh.EdgeCount == 0:
message += "Fluid sections defined but FEM mesh has no edge elements.\n"
return message
## Sets base_name
# @param self The python object self
# @param base_name base name of .inp/.frd file (without extension). It is used to construct .inp file path that is passed to CalculiX ccx
# @param base_name base name of .inp/.frd file (without extension).
# It is used to construct .inp file path that is passed to CalculiX ccx
def set_base_name(self, base_name=None):
if base_name is None:
self.base_name = ""
@@ -483,7 +564,8 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
# Normally inp file name is set set by write_inp_file
# Can be used to read mock calculations file
# @param self The python object self
# @inp_file_name .inp file name. If empty the .inp file path is constructed from working_dir, base_name and string ".inp"
# @inp_file_name .inp file name. If empty the .inp file path is constructed
# from working_dir, base_name and string ".inp"
def set_inp_file_name(self, inp_file_name=None):
if inp_file_name is not None:
self.inp_file_name = inp_file_name
@@ -504,7 +586,11 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
FreeCAD.Console.PrintMessage("Dir given as parameter \'{}\' doesn't exist.\n".format(self.working_dir))
pass
else:
FreeCAD.Console.PrintError("Dir given as parameter \'{}\' doesn't exist and create parameter is set to False.\n".format(self.working_dir))
FreeCAD.Console.PrintError(
"Dir given as parameter \'{}\' doesn't exist "
"and create parameter is set to False.\n"
.format(self.working_dir)
)
self.working_dir = femutils.get_pref_working_dir(self.solver)
FreeCAD.Console.PrintMessage("Dir \'{}\' will be used instead.\n".format(self.working_dir))
elif fem_general_prefs.GetBool("OverwriteSolverWorkingDirectory", True) is False:
@@ -725,7 +811,11 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
progress_bar.stop()
if ret_code or self.ccx_stderr:
if ret_code == 201 and self.solver.AnalysisType == 'check':
FreeCAD.Console.PrintMessage('It seams we run into NOANALYSIS problem, thus workaround for wrong exit code for *NOANALYSIS check and set ret_code to 0.\n')
FreeCAD.Console.PrintMessage(
'It seams we run into NOANALYSIS problem, '
'thus workaround for wrong exit code for *NOANALYSIS check '
'and set ret_code to 0.\n'
)
# https://forum.freecadweb.org/viewtopic.php?f=18&t=31303&start=10#p260743
ret_code = 0
else:
@@ -801,12 +891,22 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
FreeCAD.Console.PrintMessage(command_for_withoutmatnodes + '\n')
if FreeCAD.GuiUp:
import FreeCADGui
FreeCADGui.doCommand(command_for_withoutmatnodes) # with this the list without_material_elemnodes will be available for further user interaction
# with this the list without_material_elemnodes
# will be available for further user interaction
FreeCADGui.doCommand(command_for_withoutmatnodes)
FreeCAD.Console.PrintMessage('\n')
FreeCADGui.doCommand(command_to_highlight)
FreeCAD.Console.PrintMessage('\nFollowing some commands to copy which highlight the elements without materials or to reset the highlighted nodes:\n')
FreeCAD.Console.PrintMessage(
'\nFollowing some commands to copy. '
'They will highlight the elements without materials '
'or to reset the highlighted nodes:\n'
)
FreeCAD.Console.PrintMessage(command_to_highlight + '\n')
FreeCAD.Console.PrintMessage('Gui.ActiveDocument.' + self.mesh.Name + '.HighlightedNodes = []\n\n') # command to reset the Highlighted Nodes
# command to reset the Highlighted Nodes
FreeCAD.Console.PrintMessage(
'Gui.ActiveDocument.{}.HighlightedNodes = []\n\n'
.format(self.mesh.Name)
)
return True
else:
return False
@@ -836,12 +936,22 @@ class FemToolsCcx(QtCore.QRunnable, QtCore.QObject):
FreeCAD.Console.PrintMessage(command_for_nonposjacnodes + '\n')
if FreeCAD.GuiUp:
import FreeCADGui
FreeCADGui.doCommand(command_for_nonposjacnodes) # with this the list nonpositive_jacobian_elenodes will be available for further user interaction
# with this the list nonpositive_jacobian_elenodes
# will be available for further user interaction
FreeCADGui.doCommand(command_for_nonposjacnodes)
FreeCAD.Console.PrintMessage('\n')
FreeCADGui.doCommand(command_to_highlight)
FreeCAD.Console.PrintMessage('\nFollowing some commands to copy which highlight the nonpositive jacobians or to reset the highlighted nodes:\n')
FreeCAD.Console.PrintMessage(
'\nFollowing some commands to copy. '
'They highlight the nonpositive jacobians '
'or to reset the highlighted nodes:\n'
)
FreeCAD.Console.PrintMessage(command_to_highlight + '\n')
FreeCAD.Console.PrintMessage('Gui.ActiveDocument.' + self.mesh.Name + '.HighlightedNodes = []\n\n') # command to reset the Highlighted Nodes
# command to reset the Highlighted Nodes
FreeCAD.Console.PrintMessage(
'Gui.ActiveDocument.{}.HighlightedNodes = []\n\n'
.format(self.mesh.Name)
)
return True
else:
return False