diff --git a/src/Mod/Fem/femmesh/gmshtools.py b/src/Mod/Fem/femmesh/gmshtools.py index d4162b093d..b7455dc50b 100644 --- a/src/Mod/Fem/femmesh/gmshtools.py +++ b/src/Mod/Fem/femmesh/gmshtools.py @@ -29,7 +29,7 @@ __url__ = "http://www.freecadweb.org" import FreeCAD import Fem -from . import meshtools as FemMeshTools +from . import meshtools from FreeCAD import Units import subprocess import tempfile @@ -75,7 +75,6 @@ class GmshTools(): self.dimension = self.mesh_obj.ElementDimension # Algorithm2D - # known_mesh_algorithm_2D = ['Automatic', 'MeshAdapt', 'Delaunay', 'Frontal', 'BAMG', 'DelQuad'] algo2D = self.mesh_obj.Algorithm2D if algo2D == 'Automatic': self.algorithm2D = '2' @@ -93,7 +92,6 @@ class GmshTools(): self.algorithm2D = '2' # Algorithm3D - # known_mesh_algorithm_3D = ['Automatic', 'Delaunay', 'New Delaunay', 'Frontal', 'Frontal Delaunay', 'Frontal Hex', 'MMG3D', 'R-tree'] algo3D = self.mesh_obj.Algorithm3D if algo3D == 'Automatic': self.algorithm3D = '1' @@ -180,11 +178,16 @@ class GmshTools(): self.dimension = '0' elif shty == 'Compound': # print(' Found a ' + shty) - FreeCAD.Console.PrintLog(" Found a Compound. Since it could contain any kind of shape dimension 3 is used.\n") + FreeCAD.Console.PrintLog( + " Found a Compound. Since it could contain" + "any kind of shape dimension 3 is used.\n" + ) self.dimension = '3' # dimension 3 works for 2D and 1d shapes as well else: self.dimension = '0' - FreeCAD.Console.PrintError('Could not retrieve Dimension from shape type. Please choose dimension.') + FreeCAD.Console.PrintError( + 'Could not retrieve Dimension from shape type. Please choose dimension.' + ) elif self.dimension == '3D': self.dimension = '3' elif self.dimension == '2D': @@ -215,11 +218,15 @@ class GmshTools(): print(' ' + self.temp_file_geo) def get_gmsh_command(self): - gmsh_std_location = FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Fem/Gmsh").GetBool("UseStandardGmshLocation") + gmsh_std_location = FreeCAD.ParamGet( + "User parameter:BaseApp/Preferences/Mod/Fem/Gmsh" + ).GetBool("UseStandardGmshLocation") if gmsh_std_location: if system() == "Windows": gmsh_path = FreeCAD.getHomePath() + "bin/gmsh.exe" - FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Fem/Gmsh").SetString("gmshBinaryPath", gmsh_path) + FreeCAD.ParamGet( + "User parameter:BaseApp/Preferences/Mod/Fem/Gmsh" + ).SetString("gmshBinaryPath", gmsh_path) self.gmsh_bin = gmsh_path elif system() == "Linux": p1 = subprocess.Popen(['which', 'gmsh'], stdout=subprocess.PIPE) @@ -231,7 +238,8 @@ class GmshTools(): 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" + "Please install Gmsh or set path to binary " + "in FEM preferences tab Gmsh.\n" ) FreeCAD.Console.PrintError(error_message) raise Exception(error_message) @@ -245,7 +253,9 @@ class GmshTools(): raise Exception(error_message) else: if not self.gmsh_bin: - self.gmsh_bin = FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Fem/Gmsh").GetString("gmshBinaryPath", "") + self.gmsh_bin = FreeCAD.ParamGet( + "User parameter:BaseApp/Preferences/Mod/Fem/Gmsh" + ).GetString("gmshBinaryPath", "") if not self.gmsh_bin: # in prefs not set, we will try to use something reasonable if system() == "Linux": self.gmsh_bin = "gmsh" @@ -266,7 +276,7 @@ class GmshTools(): else: print(' Mesh group objects, we need to get the elements.') for mg in self.mesh_obj.MeshGroupList: - new_group_elements = FemMeshTools.get_mesh_group_elements(mg, self.part_obj) + new_group_elements = meshtools.get_mesh_group_elements(mg, self.part_obj) for ge in new_group_elements: if ge not in self.group_elements: self.group_elements[ge] = new_group_elements[ge] @@ -274,11 +284,16 @@ class GmshTools(): FreeCAD.Console.PrintError(" A group with this name exists already.\n") # group meshing for analysis - analysis_group_meshing = FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Fem/General").GetBool("AnalysisGroupMeshing", False) + analysis_group_meshing = FreeCAD.ParamGet( + "User parameter:BaseApp/Preferences/Mod/Fem/General" + ).GetBool("AnalysisGroupMeshing", False) if self.analysis and analysis_group_meshing: print(' Group meshing for analysis.') self.group_nodes_export = True - new_group_elements = FemMeshTools.get_analysis_group_elements(self.analysis, self.part_obj) + new_group_elements = meshtools.get_analysis_group_elements( + self.analysis, + self.part_obj + ) for ge in new_group_elements: if ge not in self.group_elements: self.group_elements[ge] = new_group_elements[ge] @@ -305,7 +320,9 @@ class GmshTools(): if self.mesh_obj.MeshRegionList: # 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"): + 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. ' @@ -325,20 +342,30 @@ 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, + # 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.") - # print(" But we are going to try to find it in the Shape to mesh :-)") + FreeCAD.Console.PrintLog( + " One element of the meshregion {} is " + "not an element of the Part to mesh.\n" + "But we are going to try to find it in " + "the Shape to mesh :-)\n" + .format(mr_obj.Name) + ) search_ele_in_shape_to_mesh = True for elems in sub[1]: # 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 - # 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) + # we're going to try to find the element in the + # Shape to mesh and use the found element as elems + # the method getElement(element) + # does not return Solid elements + ele_shape = meshtools.get_element(sub[0], elems) + found_element = meshtools.find_element_in_shape( + self.part_obj.Shape, ele_shape + ) if found_element: elems = found_element else: @@ -349,30 +376,41 @@ class GmshTools(): ) # print(elems) # element if elems not in self.ele_length_map: - self.ele_length_map[elems] = Units.Quantity(mr_obj.CharacteristicLength).Value + 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" + "The element {} of the meshregion {} has " + "been added to another mesh region.\n" + .format(elems, mr_obj.Name) ) else: FreeCAD.Console.PrintError( - "The meshregion: " + mr_obj.Name + " is not used to create the mesh because the reference list is empty.\n" + 'The meshregion: {} is not used to create the mesh ' + 'because the reference list is empty.\n' + .format(mr_obj.Name) ) else: FreeCAD.Console.PrintError( - "The meshregion: " + mr_obj.Name + " is not used to create the mesh because the CharacteristicLength is 0.0 mm.\n" + 'The meshregion: {} is not used to create the ' + 'mesh because the CharacteristicLength is 0.0 mm.\n' + .format(mr_obj.Name) ) 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) + # the method getElement(element) does not return Solid elements + ele_shape = meshtools.get_element(self.part_obj, eleml) + ele_vertexes = meshtools.get_vertexes_by_element(self.part_obj.Shape, ele_shape) self.ele_node_map[eleml] = ele_vertexes print(' {}'.format(self.ele_length_map)) print(' {}'.format(self.ele_node_map)) def get_boundary_layer_data(self): # mesh boundary layer - # currently only one boundary layer setting object is allowed, but multiple boundary can be selected - # Mesh.CharacteristicLengthMin, must be zero, or a value less than first inflation layer height + # currently only one boundary layer setting object is allowed + # but multiple boundary can be selected + # Mesh.CharacteristicLengthMin, must be zero + # or a value less than first inflation layer height if not self.mesh_obj.MeshBoundaryLayerList: # print(' No mesh boundary layer setting document object.') pass @@ -398,46 +436,66 @@ class GmshTools(): # 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.") - # print(" But we're going to find it in the Shape to mesh :-)") + FreeCAD.Console.PrintLog( + " One element of the mesh boundary layer {} is " + "not an element of the Part to mesh.\n" + "But we are going to try to find it in " + "the Shape to mesh :-)\n" + .format(mr_obj.Name) + ) search_ele_in_shape_to_mesh = True for elems in sub[1]: # 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 + # we try to find the element it in the Shape to mesh + # and use the found element as elems # 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) + ele_shape = meshtools.get_element(sub[0], elems) + found_element = meshtools.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 {} could not be found " - "in the Part to mesh. It will be ignored.\n" + "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 + # 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 {} of the mesh boundary layer {} has been added " + "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 - setting['thickness'] = sum([setting['hwall_n'] * setting['ratio'] ** i for i in range(mr_obj.NumberOfLayers)]) - setting['hwall_t'] = setting['thickness'] # setting['hwall_n'] * 5 # tangential cell dimension + setting['thickness'] = sum([ + setting['hwall_n'] * setting['ratio'] ** i for i in range( + mr_obj.NumberOfLayers + ) + ]) + # setting['hwall_n'] * 5 # tangential cell dimension + setting['hwall_t'] = setting['thickness'] - # hfar: cell dimension outside boundary should be set later if some character length is set - if self.clmax > setting['thickness'] * 0.8 and self.clmax < setting['thickness'] * 1.6: + # hfar: cell dimension outside boundary + # should be set later if some character length is set + if self.clmax > setting['thickness'] * 0.8 \ + and self.clmax < setting['thickness'] * 1.6: setting['hfar'] = self.clmax else: - setting['hfar'] = setting['thickness'] # set a value for safety, it may works as background mesh cell size + # set a value for safety, it may works as background mesh cell size + setting['hfar'] = setting['thickness'] # from face name -> face id is done in geo file write up # TODO: fan angle setup is not implemented yet if self.dimension == '2': @@ -445,15 +503,21 @@ class GmshTools(): elif self.dimension == '3': setting['FacesList'] = belem_list else: - FreeCAD.Console.PrintError("boundary layer is only supported for 2D and 3D mesh") + 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" + 'The mesh boundary layer: {} is not used to create ' + 'the mesh because the reference list is empty.\n' + .format(mr_obj.Name) ) 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" + 'The mesh boundary layer: {} is not used to create ' + 'the mesh because the min thickness is 0.0 mm.\n' + .format(mr_obj.Name) ) print(' {}'.format(self.bl_setting_list)) @@ -469,7 +533,8 @@ class GmshTools(): for k in item: v = item[k] if k in set(['EdgesList', 'FacesList']): - # the element name of FreeCAD which starts with 1 (example: 'Face1'), same as Gmsh + # the element name of FreeCAD which starts + # with 1 (example: 'Face1'), same as Gmsh # el_id = int(el[4:]) # FIXME: strip `face` or `edge` prefix ele_nodes = (''.join((str(el[4:]) + ', ') for el in v)).rstrip(', ') line = prefix + '.' + str(k) + ' = {' + ele_nodes + ' };\n' @@ -501,8 +566,10 @@ class GmshTools(): if self.group_elements: # print(' We are going to have to find elements to make mesh groups for.') geo.write("// group data\n") - # we use the element name of FreeCAD which starts with 1 (example: 'Face1'), same as Gmsh - for group in sorted(self.group_elements.keys()): # for unit test we need them to have a fixed order + # we use the element name of FreeCAD which starts + # with 1 (example: 'Face1'), same as Gmsh + # for unit test we need them to have a fixed order + for group in sorted(self.group_elements.keys()): gdata = self.group_elements[group] # print(gdata) # geo.write("// " + group + "\n") @@ -526,25 +593,35 @@ class GmshTools(): if ele_nr: ele_nr = ele_nr.rstrip(', ') # print(ele_nr) - geo.write('Physical ' + physical_type + '("' + group + '") = {' + ele_nr + '};\n') + geo.write( + 'Physical ' + physical_type + '("' + group + '") = {' + ele_nr + '};\n' + ) geo.write("\n") geo.write("// Characteristic Length\n") if self.ele_length_map: - # we use the index FreeCAD which starts with 0, we need to add 1 for the index in Gmsh + # we use the index FreeCAD which starts with 0 + # we need to add 1 for the index in Gmsh geo.write("// Characteristic Length according CharacteristicLengthMap\n") for e in self.ele_length_map: - ele_nodes = (''.join((str(n + 1) + ', ') for n in self.ele_node_map[e])).rstrip(', ') + ele_nodes = ( + ''.join((str(n + 1) + ', ') for n in self.ele_node_map[e]) + ).rstrip(', ') geo.write("// " + e + "\n") - geo.write("Characteristic Length { " + ele_nodes + " } = " + str(self.ele_length_map[e]) + ";\n") + geo.write( + "Characteristic Length { {} } = {};\n" + .format(ele_nodes, self.ele_length_map[e]) + ) geo.write("\n") - # boundary layer generation may need special setup of Gmsh properties, set them in Gmsh TaskPanel + # boundary layer generation may need special setup + # of Gmsh properties, set them in Gmsh TaskPanel self.write_boundary_layer(geo) geo.write("// min, max Characteristic Length\n") geo.write("Mesh.CharacteristicLengthMax = " + str(self.clmax) + ";\n") if len(self.bl_setting_list): - # if minLength must smaller than first layer of boundary_layer, it is safer to set it as zero (default value) to avoid error + # if minLength must smaller than first layer of boundary_layer + # it is safer to set it as zero (default value) to avoid error geo.write("Mesh.CharacteristicLengthMin = " + str(0) + ";\n") else: geo.write("Mesh.CharacteristicLengthMin = " + str(self.clmin) + ";\n") @@ -566,28 +643,48 @@ class GmshTools(): geo.write("Mesh.OptimizeNetgen = 0;\n") # higher order mesh optimizing if hasattr(self.mesh_obj, 'HighOrderOptimize') and self.mesh_obj.HighOrderOptimize is True: - geo.write("Mesh.HighOrderOptimize = 1; // for more HighOrderOptimize parameter check http://gmsh.info/doc/texinfo/gmsh.html\n") + geo.write( + "Mesh.HighOrderOptimize = 1; // for more HighOrderOptimize " + "parameter check http://gmsh.info/doc/texinfo/gmsh.html\n" + ) else: - geo.write("Mesh.HighOrderOptimize = 0; // for more HighOrderOptimize parameter check http://gmsh.info/doc/texinfo/gmsh.html\n") + geo.write( + "Mesh.HighOrderOptimize = 0; // for more HighOrderOptimize " + "parameter check http://gmsh.info/doc/texinfo/gmsh.html\n" + ) geo.write("\n") geo.write("// mesh order\n") geo.write("Mesh.ElementOrder = " + self.order + ";\n") geo.write("\n") - geo.write("// mesh algorithm, only a few algorithms are usable with 3D boundary layer generation\n") - geo.write("// 2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=BAMG, 8=DelQuad)\n") + geo.write( + "// mesh algorithm, only a few algorithms are " + "usable with 3D boundary layer generation\n" + ) + geo.write( + "// 2D mesh algorithm (1=MeshAdapt, 2=Automatic, " + "5=Delaunay, 6=Frontal, 7=BAMG, 8=DelQuad)\n" + ) if len(self.bl_setting_list) and self.dimension == 3: geo.write("Mesh.Algorithm = " + 'DelQuad' + ";\n") # Frontal/DelQuad are tested else: geo.write("Mesh.Algorithm = " + self.algorithm2D + ";\n") - geo.write("// 3D mesh algorithm (1=Delaunay, 2=New Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)\n") + geo.write( + "// 3D mesh algorithm (1=Delaunay, 2=New Delaunay, 4=Frontal, " + "5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)\n" + ) geo.write("Mesh.Algorithm3D = " + self.algorithm3D + ";\n") geo.write("\n") geo.write("// meshing\n") - # remove duplicate vertices, see https://forum.freecadweb.org/viewtopic.php?f=18&t=21571&start=20#p179443 + # remove duplicate vertices + # see https://forum.freecadweb.org/viewtopic.php?f=18&t=21571&start=20#p179443 if hasattr(self.mesh_obj, 'CoherenceMesh') and self.mesh_obj.CoherenceMesh is True: - geo.write("Geometry.Tolerance = " + str(self.geotol) + "; // set geometrical tolerance (also used for merging nodes)\n") + geo.write( + "Geometry.Tolerance = {}; // set geometrical " + "tolerance (also used for merging nodes)\n" + .format(self.geotol) + ) geo.write("Mesh " + self.dimension + ";\n") geo.write("Coherence Mesh; // Remove duplicate vertices\n") else: @@ -599,7 +696,9 @@ class GmshTools(): geo.write("// For each group save not only the elements but the nodes too.;\n") geo.write("Mesh.SaveGroupsOfNodes = 1;\n") # 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( + "// 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') @@ -608,7 +707,10 @@ class GmshTools(): geo.write("// Gmsh documentation:\n") geo.write("// http://gmsh.info/doc/texinfo/gmsh.html#Mesh\n") geo.write("//\n") - geo.write("// We do not check if something went wrong, like negative jacobians etc. You can run Gmsh manually yourself: \n") + geo.write( + "// We do not check if something went wrong, like negative " + "jacobians etc. You can run Gmsh manually yourself: \n" + ) geo.write("//\n") geo.write("// to see full Gmsh log, run in bash:\n") geo.write("// " + self.gmsh_bin + " - " + self.temp_file_geo + "\n") @@ -621,12 +723,19 @@ class GmshTools(): comandlist = [self.gmsh_bin, '-', self.temp_file_geo] # print(comandlist) try: - p = subprocess.Popen(comandlist, shell=False, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + p = subprocess.Popen( + comandlist, + shell=False, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE + ) output, error = p.communicate() if sys.version_info.major >= 3: output = output.decode('utf-8') error = error.decode('utf-8') - # print(output) # stdout is still cut at some point but the warnings are in stderr and thus printed :-) + # stdout is still cut at some point + # but the warnings are in stderr and thus printed :-) + # print(output) # print(error) except: error = 'Error executing: {}\n'.format(" ".join(comandlist)) @@ -638,9 +747,9 @@ class GmshTools(): if not self.error: fem_mesh = Fem.read(self.temp_file_mesh) self.mesh_obj.FemMesh = fem_mesh - print(' The Part should have a pretty new FEM mesh!') + FreeCAD.Console.PrintError(' The Part should have a pretty new FEM mesh!\n') else: - print('No mesh was created.') + FreeCAD.Console.PrintError('No mesh was created.\n') del self.temp_file_geometry del self.temp_file_mesh diff --git a/src/Mod/Fem/femmesh/meshtools.py b/src/Mod/Fem/femmesh/meshtools.py index 63c285303d..aa906e352f 100644 --- a/src/Mod/Fem/femmesh/meshtools.py +++ b/src/Mod/Fem/femmesh/meshtools.py @@ -30,7 +30,11 @@ __url__ = "http://www.freecadweb.org" import FreeCAD -def get_femnodes_by_femobj_with_references(femmesh, femobj): +# ************************************************************************************************ +def get_femnodes_by_femobj_with_references( + femmesh, + femobj +): node_set = [] if femmesh.GroupCount: node_set = get_femmesh_groupdata_sets_by_name(femmesh, femobj, 'Node') @@ -51,7 +55,13 @@ def get_femnodes_by_femobj_with_references(femmesh, femobj): return node_set -def get_femelements_by_references(femmesh, femelement_table, references, femnodes_ele_table=None): +# ************************************************************************************************ +def get_femelements_by_references( + femmesh, + femelement_table, + references, + femnodes_ele_table=None +): '''get the femelements for a list of references ''' references_femelements = [] @@ -76,7 +86,11 @@ def get_femelements_by_references(femmesh, femelement_table, references, femnode return references_femelements -def get_femnodes_by_references(femmesh, references): +# ************************************************************************************************ +def get_femnodes_by_references( + femmesh, + references +): '''get the femnodes for a list of references ''' references_femnodes = [] @@ -89,7 +103,10 @@ def get_femnodes_by_references(femmesh, references): return list(set(references_femnodes)) # removes duplicate nodes, sorts node order -def get_femnodes_by_refshape(femmesh, ref): +def get_femnodes_by_refshape( + femmesh, + ref +): nodes = [] for refelement in ref[1]: # the following method getElement(element) does not return Solid elements @@ -118,7 +135,10 @@ def get_femnodes_by_refshape(femmesh, ref): return nodes -def get_femelement_table(femmesh): +# ************************************************************************************************ +def get_femelement_table( + femmesh +): """ get_femelement_table(femmesh): { elementid : [ nodeid, nodeid, ... , nodeid ] }""" femelement_table = {} if is_solid_femmesh(femmesh): @@ -135,7 +155,10 @@ def get_femelement_table(femmesh): return femelement_table -def get_femelement_volumes_table(femmesh): +# ************************************************************************************************ +def get_femelement_volumes_table( + femmesh +): """ get_femelement_volumes_table(femmesh): { elementid : [ nodeid, nodeid, ... , nodeid ] }""" table = {} for i in femmesh.Volumes: @@ -143,7 +166,11 @@ def get_femelement_volumes_table(femmesh): return table -def get_femelement_faces_table(femmesh, faces_only=None): +# ************************************************************************************************ +def get_femelement_faces_table( + femmesh, + faces_only=None +): """ get_femelement_faces_table(femmesh): { elementid : [ nodeid, nodeid, ... , nodeid ] }""" table = {} if not faces_only: @@ -153,7 +180,11 @@ def get_femelement_faces_table(femmesh, faces_only=None): return table -def get_femelement_edges_table(femmesh, edges_only=None): +# ************************************************************************************************ +def get_femelement_edges_table( + femmesh, + edges_only=None +): """ get_femelement_edges_table(femmesh): { elementid : [ nodeid, nodeid, ... , nodeid ] }""" table = {} if not edges_only: @@ -163,16 +194,25 @@ def get_femelement_edges_table(femmesh, edges_only=None): return table -def get_femnodes_ele_table(femnodes_mesh, femelement_table): +# ************************************************************************************************ +def get_femnodes_ele_table( + femnodes_mesh, + femelement_table +): '''the femnodes_ele_table contains for each node its membership in elements {nodeID : [[eleID, NodePosition], [], ...], nodeID : [[], [], ...], ...} stored informatation are: - element number, the number of nodes per element, the position of the node in the element. - The position of the node in the element is coded as a set bit at that position in a bit array (integer) - Fixme: the number of nodes per element should be replaced by the type of the element + element number, the number of nodes per element + the position of the node in the element. + The position of the node in the element is coded + as a set bit at that position in a bit array (integer) + FIXME: the number of nodes per element should be + replaced by the type of the element but I did not know, how to get this from the mesh. - Since the femelement_table contains either volume or face or edgemesh the femnodes_ele_table only - has either volume or face or edge elements, see get_femelement_table() + Since the femelement_table contains either + volume or face or edgemesh the femnodes_ele_table only + has either volume or face or edge elements + see get_femelement_table() ''' femnodes_ele_table = {} # node_dict in ulrichs class for n in femnodes_mesh: # initialize it with sorted node keys and empty lists @@ -184,12 +224,18 @@ def get_femnodes_ele_table(femnodes_mesh, femelement_table): for ele_node in ele_list: femnodes_ele_table[ele_node].append([ele, pos]) pos = pos << 1 - FreeCAD.Console.PrintLog('len femnodes_ele_table: ' + str(len(femnodes_ele_table)) + '\n') - # FreeCAD.Console.PrintMessage('femnodes_ele_table: {}\n'.format(femnodes_ele_table)) + FreeCAD.Console.PrintLog( + 'len femnodes_ele_table: {}\n' + .format(len(femnodes_ele_table)) + ) + FreeCAD.Console.PrintLog('femnodes_ele_table: {}\n'.format(femnodes_ele_table)) return femnodes_ele_table -def get_copy_of_empty_femelement_table(femelement_table): +# ************************************************************************************************ +def get_copy_of_empty_femelement_table( + femelement_table +): '''{eleID : 0, eleID : 0, ...} ''' empty_femelement_table = {} @@ -198,11 +244,16 @@ def get_copy_of_empty_femelement_table(femelement_table): return empty_femelement_table.copy() -def get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_set): +# ************************************************************************************************ +def get_bit_pattern_dict( + femelement_table, + femnodes_ele_table, + node_set +): '''Now we are looking for nodes inside of the Faces = filling the bit_pattern_dict {eleID : [lenEleNodes, binary_position]} see forum post for a very good explanation of what's really happening - http://forum.freecadweb.org/viewtopic.php?f=18&p=141133&sid=013c93f496a63872951d2ce521702ffa#p141108 + https://forum.freecadweb.org/viewtopic.php?f=18&t=17318&start=50#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) @@ -212,7 +263,7 @@ def get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_set): ''' FreeCAD.Console.PrintLog('len femnodes_ele_table: ' + str(len(femnodes_ele_table)) + '\n') FreeCAD.Console.PrintLog('len node_set: ' + str(len(node_set)) + '\n') - # FreeCAD.Console.PrintMessage('node_set: {}\n'.format(node_set)) + FreeCAD.Console.PrintLog('node_set: {}\n'.format(node_set)) bit_pattern_dict = get_copy_of_empty_femelement_table(femelement_table) # # initializing the bit_pattern_dict for ele in femelement_table: @@ -226,7 +277,10 @@ def get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_set): return bit_pattern_dict -def get_ccxelement_faces_from_binary_search(bit_pattern_dict): +# ************************************************************************************************ +def get_ccxelement_faces_from_binary_search( + bit_pattern_dict +): '''get the CalculiX element face numbers ''' tet10_mask = { @@ -283,7 +337,12 @@ def get_ccxelement_faces_from_binary_search(bit_pattern_dict): return faces -def get_femelements_by_femnodes_bin(femelement_table, femnodes_ele_table, node_list): +# ************************************************************************************************ +def get_femelements_by_femnodes_bin( + femelement_table, + femnodes_ele_table, + node_list +): '''for every femelement of femelement_table if all nodes of the femelement are in node_list, the femelement is added to the list which is returned @@ -298,12 +357,17 @@ def get_femelements_by_femnodes_bin(femelement_table, femnodes_ele_table, node_l 15: 32767, 20: 1048575} # Now we are looking for nodes inside of the Volumes = filling the bit_pattern_dict - FreeCAD.Console.PrintMessage('len femnodes_ele_table: ' + str(len(femnodes_ele_table)) + '\n') + FreeCAD.Console.PrintMessage( + 'len femnodes_ele_table: {}\n' + .format(len(femnodes_ele_table)) + ) bit_pattern_dict = get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_list) # search ele_list = [] # The ele_list contains the result of the search. for ele in bit_pattern_dict: - # FreeCAD.Console.PrintMessage('bit_pattern_dict[ele][0]: {}\n'.format(bit_pattern_dict[ele][0])) + FreeCAD.Console.PrintLog( + 'bit_pattern_dict[ele][0]: {}\n'.format(bit_pattern_dict[ele][0]) + ) if bit_pattern_dict[ele][1] == vol_masks[bit_pattern_dict[ele][0]]: ele_list.append(ele) FreeCAD.Console.PrintMessage('found Volumes: {}\n'.format(len(ele_list))) @@ -311,7 +375,11 @@ def get_femelements_by_femnodes_bin(femelement_table, femnodes_ele_table, node_l return ele_list -def get_femelements_by_femnodes_std(femelement_table, node_list): +# ************************************************************************************************ +def get_femelements_by_femnodes_std( + femelement_table, + node_list +): '''for every femelement of femelement_table if all nodes of the femelement are in node_list, the femelement is added to the list which is returned @@ -324,19 +392,27 @@ def get_femelements_by_femnodes_std(femelement_table, node_list): for nodeID in femelement_table[elementID]: if nodeID in node_list: nodecount = nodecount + 1 - if nodecount == len(femelement_table[elementID]): # all nodes of the element are in the node_list! + # all nodes of the element are in the node_list! + if nodecount == len(femelement_table[elementID]): e.append(elementID) return e -def get_femvolumeelements_by_femfacenodes(femelement_table, node_list): +def get_femvolumeelements_by_femfacenodes( + femelement_table, + node_list +): '''assume femelement_table only has volume elements for every femvolumeelement of femelement_table for tetra4 and tetra10 the C++ methods could be used --> test again to be sure - if hexa8 volume element --> if exact 4 element nodes are in node_list --> add femelement - if hexa20 volume element --> if exact 8 element nodes are in node_list --> add femelement - if penta6 volume element --> if exact 3 or 6 element nodes are in node_list --> add femelement - if penta15 volume element --> if exact 6 or 8 element nodes are in node_list --> add femelement + if hexa8 volume element + --> if exact 4 element nodes are in node_list --> add femelement + if hexa20 volume element + --> if exact 8 element nodes are in node_list --> add femelement + if penta6 volume element + --> if exact 3 or 6 element nodes are in node_list --> add femelement + if penta15 volume element + --> if exact 6 or 8 element nodes are in node_list --> add femelement e: elementlist nodes: nodelist ''' e = [] # elementlist @@ -381,13 +457,21 @@ def get_femvolumeelements_by_femfacenodes(femelement_table, node_list): e.append(elementID) else: FreeCAD.Console.PrintError( - 'Error in get_femvolumeelements_by_femfacenodes(): unknown volume element: ' + el_nd_ct + '\n' + 'Error in get_femvolumeelements_by_femfacenodes(): ' + 'unknown volume element: {}\n' + .format(el_nd_ct) ) # FreeCAD.Console.PrintMessage('{}\n'.format(sorted(e))) return e -def get_femelement_sets(femmesh, femelement_table, fem_objects, femnodes_ele_table=None): +# ************************************************************************************************ +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 @@ -396,9 +480,12 @@ def get_femelement_sets(femmesh, femelement_table, fem_objects, femnodes_ele_tab 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" + "Constraint: {} --> We're going to search " + "in the mesh for the element ID's.\n" + .format(obj.Name) ) - fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) # unique short identifier + # unique short identifier + fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) if obj.References: ref_shape_femelements = [] ref_shape_femelements = get_femelements_by_references( @@ -426,21 +513,29 @@ def get_femelement_sets(femmesh, femelement_table, fem_objects, femnodes_ele_tab if femelements_count_ok(len(femelement_table), count_femelements): return True else: - FreeCAD.Console.PrintError('Error in get_femelement_sets -- > femelements_count_ok() failed!\n') + FreeCAD.Console.PrintError( + 'Error in get_femelement_sets -- > femelements_count_ok() failed!\n' + ) return False -def get_femelement_direction1D_set(femmesh, femelement_table, beamrotation_objects, theshape=None): +# ************************************************************************************************ +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 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], - 'direction' : direction, - 'normal' : normal}, - ... - ] + beamrot_obj['FEMRotations1D'] = [{ + 'ids' : [theids], + 'direction' : direction, + 'normal' : normal + }, ... ] ''' if len(beamrotation_objects) == 0: # no beamrotation document object, all beams use standard rotation of 0 degree (angle) @@ -454,7 +549,8 @@ def get_femelement_direction1D_set(femmesh, femelement_table, beamrotation_objec # 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 + # 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) @@ -465,20 +561,24 @@ def get_femelement_direction1D_set(femmesh, femelement_table, beamrotation_objec beamrotation_objects[0]['FEMRotations1D'] = rotations_ids beamrotation_objects[0]['ShortName'] = 'R0' elif len(beamrotation_objects) > 1: - # multiple beam rotation document objects, rotations defined by reference shapes, TODO implement this + # multiple beam rotation document objects + # rotations defined by reference shapes, TODO implement this # do not forget all the corner cases: # one beam rotation object, but not all edges are ref shapes # more than one beam rotation object, but not all edges are in the ref shapes # for the both cases above, all other edges get standard rotation. - # more than one beam rotation objects and on has no ref shapes, all edges no in an rotation object use this rotation + # more than one beam rotation objects and on has no ref shapes + # all edges no in an rotation object use this rotation # one edge is in more than one beam rotation object, error # pre check, only one beam rotation with empty ref shapes is allowed - # we need theshape for multiple rotations too, because of the corner cases mentioned above + # we need theshape for multiple rotations too + # because of the corner cases mentioned above FreeCAD.Console.PrintError('Multiple Rotations not yet supported!\n') for rot_object in beamrotation_objects: # debug output FreeCAD.Console.PrintMessage('{}\n'.format(rot_object['FEMRotations1D'])) +# ************************************************************************************************ def get_femelement_directions_theshape(femmesh, femelement_table, theshape): # see get_femelement_direction1D_set rotations_ids = [] @@ -500,6 +600,7 @@ def get_femelement_directions_theshape(femmesh, femelement_table, theshape): return rotations_ids +# ************************************************************************************************ def get_beam_normal(beam_direction, defined_angle): import math vector_a = beam_direction @@ -525,44 +626,52 @@ def get_beam_normal(beam_direction, defined_angle): else: temp_valz = 0 - Dot_product_check_x = None - Dot_product_check_y = None - Dot_product_check_z = None - Dot_product_check_nt = None + # Dot_product_check + dot_x = None + dot_y = None + dot_z = None + dot_nt = None if vector_a[0] != 0 and vector_a[1] == 0 and vector_a[2] == 0: normal_n = [temp_valx, nx, ny] - Dot_product_check_x = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + dot_x = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] elif vector_a[0] == 0 and vector_a[1] != 0 and vector_a[2] == 0: normal_n = [nx, temp_valy, ny] - Dot_product_check_y = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + dot_y = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] elif vector_a[0] == 0 and vector_a[1] == 0 and vector_a[2] != 0: normal_n = [nx, ny, temp_valz] - Dot_product_check_z = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + dot_z = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] elif vector_a[0] == 0 and vector_a[1] != 0 and vector_a[2] != 0: normal_n = [nx, temp_valy, ny] - Dot_product_check_y = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + dot_y = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] elif vector_a[0] != 0 and vector_a[1] == 0 and vector_a[2] != 0: normal_n = [nx, ny, temp_valz] - Dot_product_check_z = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + dot_z = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] else: normal_n = [temp_valx, nx, ny] - Dot_product_check_nt = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + dot_nt = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - Dot_product_check = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - # FreeCAD.Console.PrintMessage('{}\n'.format(Dot_product_check)) - # FreeCAD.Console.PrintMessage('{}\n'.format(normal_n)) + dot = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + FreeCAD.Console.PrintLog('{}\n'.format(dot)) + FreeCAD.Console.PrintLog('{}\n'.format(normal_n)) - # dummy usage of the axis Dot_product_check to get flake8 quiet - del Dot_product_check_x, Dot_product_check_y, Dot_product_check_z, Dot_product_check, Dot_product_check_nt + # dummy usage of the axis dot to get flake8 quiet + del dot_x, dot_y, dot_z, dot, dot_nt return normal_n -def get_femmesh_groupdata_sets_by_name(femmesh, fem_object, group_data_type): +# ************************************************************************************************ +def get_femmesh_groupdata_sets_by_name( + femmesh, + fem_object, + group_data_type +): # get ids from femmesh groupdata for reference shapes of each obj.References - # we assume the mesh group data fits with the reference shapes, no check is done in this regard !!! + # we assume the mesh group data fits with the reference shapes + # no check is done in this regard !!! # we just check for the group name and the group data type - # what happens if a reference shape was changed, but the mesh and the mesh groups were not created new !?! + # what happens if a reference shape was changed + # but the mesh and the mesh groups were not created new !?! obj = fem_object['Object'] if femmesh.GroupCount: for g in femmesh.Groups: @@ -570,22 +679,30 @@ def get_femmesh_groupdata_sets_by_name(femmesh, fem_object, group_data_type): 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' + ' found mesh group for the IDs: {}, Type: {}\n' + .format(grp_name, group_data_type) ) return femmesh.getGroupElements(g) # == ref_shape_femelements return () # an empty tuple is returned if no group data IDs where found -def get_femelement_sets_from_group_data(femmesh, fem_objects): +# ************************************************************************************************ +def get_femelement_sets_from_group_data( + femmesh, + fem_objects +): # get femelements from femmesh groupdata for reference shapes of each obj.References count_femelements = 0 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" + "Constraint: {} --> We have mesh groups. " + "We will search for appropriate group data.\n" + .format(obj.Name) ) - fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) # unique short identifier + # unique short identifier + fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) # see comments over there ! group_elements = get_femmesh_groupdata_sets_by_name(femmesh, fem_object, 'Volume') sum_group_elements += group_elements @@ -594,14 +711,19 @@ def get_femelement_sets_from_group_data(femmesh, fem_objects): # 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' + 'Error in get_femelement_sets_from_group_data -- > ' + 'femelements_count_ok() failed!\n' ) return False else: return True -def get_elset_short_name(obj, i): +# ************************************************************************************************ +def get_elset_short_name( + obj, + i +): if hasattr(obj, "Proxy") and obj.Proxy.Type == 'Fem::Material': return 'M' + str(i) elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'Fem::FemElementGeometry1D': @@ -614,16 +736,25 @@ def get_elset_short_name(obj, i): 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' + 'Error in creating short elset name ' + 'for obj: {} --> Proxy.Type: {}\n' + .format(obj.Name, obj.Proxy.Type) ) # ************************************************************************************************ # ***** methods for retrieving nodes and node load values for constraint force ******************* # ***** Vertex loads ***************************************************************************** -def get_force_obj_vertex_nodeload_table(femmesh, frc_obj): +def get_force_obj_vertex_nodeload_table( + femmesh, + frc_obj +): # force_obj_node_load_table: - # [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)] + # [ + # ('refshape_name.elemname', node_load_table), + # ..., + # ('refshape_name.elemname', node_load_table) + # ] force_obj_node_load_table = [] node_load = frc_obj.Force / len(frc_obj.References) for o, elem_tup in frc_obj.References: @@ -640,7 +771,9 @@ def get_force_obj_vertex_nodeload_table(femmesh, frc_obj): ) node = femmesh.getNodesByVertex(ref_node) elem_info_string = 'node load on shape: ' + o.Name + ':' + elem - force_obj_node_load_table.append((elem_info_string, {node[0]: node_load / node_count})) + force_obj_node_load_table.append( + (elem_info_string, {node[0]: node_load / node_count}) + ) return force_obj_node_load_table @@ -648,9 +781,18 @@ def get_force_obj_vertex_nodeload_table(femmesh, frc_obj): # get_force_obj_edge_nodeload_table # get_ref_edgenodes_table # get_ref_edgenodes_lengths -def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, frc_obj): +def get_force_obj_edge_nodeload_table( + femmesh, + femelement_table, + femnodes_mesh, + frc_obj +): # force_obj_node_load_table: - # [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)] + # [ + # ('refshape_name.elemname', node_load_table), + # ..., + # ('refshape_name.elemname', node_load_table) + # ] force_obj_node_load_table = [] sum_ref_edge_length = 0 sum_ref_edge_node_length = 0 # for debugging @@ -678,15 +820,18 @@ def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, edge_table = get_ref_edgenodes_table(femmesh, femelement_table, ref_edge) # node_length_table: - # [ (nodeID, length), ... , (nodeID, length) ] some nodes will have more than one entry + # [ (nodeID, length), ... , (nodeID, length) ] + # some nodes will have more than one entry node_length_table = get_ref_edgenodes_lengths(femnodes_mesh, edge_table) # node_sum_length_table: - # { nodeID : Length, ... , nodeID : Length } LengthSum for each node, one entry for each node + # { nodeID : Length, ... , nodeID : Length } + # LengthSum for each node, one entry for each node node_sum_length_table = get_ref_shape_node_sum_geom_table(node_length_table) # node_load_table: - # { nodeID : NodeLoad, ... , nodeID : NodeLoad } NodeLoad for each node, one entry for each node + # { nodeID : NodeLoad, ... , nodeID : NodeLoad } + # NodeLoad for each node, one entry for each node node_load_table = {} sum_node_lengths = 0 # for debugging for node in node_sum_length_table: @@ -694,7 +839,9 @@ def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, node_load_table[node] = node_sum_length_table[node] * force_per_sum_ref_edge_length ratio_refedge_lengths = sum_node_lengths / ref_edge.Length if ratio_refedge_lengths < 0.99 or ratio_refedge_lengths > 1.01: - FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n') + FreeCAD.Console.PrintError( + 'Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n' + ) FreeCAD.Console.PrintMessage(' sum_node_lengths: {}\n'.format(sum_node_lengths)) FreeCAD.Console.PrintMessage(' refedge_length: {}\n'.format(ref_edge.Length)) bad_refedge = ref_edge @@ -709,13 +856,34 @@ def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, ratio = sum_node_load / frc_obj.Force if ratio < 0.99 or ratio > 1.01: - FreeCAD.Console.PrintMessage('Deviation sum_node_load to frc_obj.Force is more than 1% : {}\n'.format(ratio)) - FreeCAD.Console.PrintMessage(' sum_ref_edge_node_length: {}\n'.format(sum_ref_edge_node_length)) - FreeCAD.Console.PrintMessage(' sum_ref_edge_length: {}\n'.format(sum_ref_edge_length)) - FreeCAD.Console.PrintMessage(' sum_node_load: {}\n'.format(sum_node_load)) - FreeCAD.Console.PrintMessage(' frc_obj.Force: {}\n'.format(frc_obj.Force)) - FreeCAD.Console.PrintMessage(' the reason could be simply a circle length --> see method get_ref_edge_node_lengths\n') - FreeCAD.Console.PrintMessage(' the reason could also be a problem in retrieving the ref_edge_node_length\n') + FreeCAD.Console.PrintMessage( + 'Deviation sum_node_load to frc_obj.Force is more than 1% : {}\n' + .format(ratio) + ) + FreeCAD.Console.PrintMessage( + ' sum_ref_edge_node_length: {}\n' + .format(sum_ref_edge_node_length) + ) + FreeCAD.Console.PrintMessage( + ' sum_ref_edge_length: {}\n' + .format(sum_ref_edge_length) + ) + FreeCAD.Console.PrintMessage( + ' sum_node_load: {}\n' + .format(sum_node_load) + ) + FreeCAD.Console.PrintMessage( + ' frc_obj.Force: {}\n' + .format(frc_obj.Force) + ) + FreeCAD.Console.PrintMessage( + ' the reason could be simply a circle length --> ' + 'see method get_ref_edge_node_lengths\n' + ) + FreeCAD.Console.PrintMessage( + ' the reason could also be a problem in ' + 'retrieving the ref_edge_node_length\n' + ) # try debugging of the last bad refedge FreeCAD.Console.PrintMessage('DEBUGGING\n') @@ -740,11 +908,13 @@ def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, if node not in bad_edge_table_nodes: bad_edge_table_nodes.append(node) FreeCAD.Console.PrintMessage('sorted(bad_edge_table_nodes)\n') - FreeCAD.Console.PrintMessage('{}\n'.format(sorted(bad_edge_table_nodes))) # should be == bad_refedge_nodes + # should be == bad_refedge_nodes + FreeCAD.Console.PrintMessage('{}\n'.format(sorted(bad_edge_table_nodes))) # import FreeCADGui # FreeCADGui.ActiveDocument.Compound_Mesh.HighlightedNodes = bad_edge_table_nodes # bad_node_length_table: - # [ (nodeID, length), ... , (nodeID, length) ] some nodes will have more than one entry + # [ (nodeID, length), ... , (nodeID, length) ] + # some nodes will have more than one entry FreeCAD.Console.PrintMessage('good_edge_table\n') good_edge_table = delete_duplicate_mesh_elements(bad_edge_table) @@ -759,12 +929,18 @@ def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, return force_obj_node_load_table -def get_ref_edgenodes_table(femmesh, femelement_table, refedge): +# ************************************************************************************************ +def get_ref_edgenodes_table( + femmesh, + femelement_table, + refedge +): edge_table = {} # { meshedgeID : ( nodeID, ... , nodeID ) } refedge_nodes = femmesh.getNodesByEdge(refedge) if is_solid_femmesh(femmesh): refedge_fem_volumeelements = [] - # if at least two nodes of a femvolumeelement are in refedge_nodes the volume is added to refedge_fem_volumeelements + # if at least two nodes of a femvolumeelement are in + # refedge_nodes the volume is added to refedge_fem_volumeelements for elem in femelement_table: nodecount = 0 for node in femelement_table[elem]: @@ -772,7 +948,8 @@ def get_ref_edgenodes_table(femmesh, femelement_table, refedge): nodecount += 1 if nodecount > 1: refedge_fem_volumeelements.append(elem) - # for every refedge_fem_volumeelement look which of its nodes is in refedge_nodes --> add all these nodes to edge_table + # for every refedge_fem_volumeelement look which of its nodes + # is in refedge_nodes --> add all these nodes to edge_table for elem in refedge_fem_volumeelements: fe_refedge_nodes = [] for node in femelement_table[elem]: @@ -785,7 +962,8 @@ def get_ref_edgenodes_table(femmesh, femelement_table, refedge): edge_table = delete_duplicate_mesh_elements(edge_table) elif is_face_femmesh(femmesh): refedge_fem_faceelements = [] - # if at least two nodes of a femfaceelement are in refedge_nodes the volume is added to refedge_fem_volumeelements + # if at least two nodes of a femfaceelement are in + # refedge_nodes the volume is added to refedge_fem_volumeelements for elem in femelement_table: nodecount = 0 for node in femelement_table[elem]: @@ -793,7 +971,8 @@ def get_ref_edgenodes_table(femmesh, femelement_table, refedge): nodecount += 1 if nodecount > 1: refedge_fem_faceelements.append(elem) - # for every refedge_fem_faceelement look which of his nodes is in refedge_nodes --> add all these nodes to edge_table + # for every refedge_fem_faceelement look which of his nodes is in + # refedge_nodes --> add all these nodes to edge_table for elem in refedge_fem_faceelements: fe_refedge_nodes = [] for node in femelement_table[elem]: @@ -805,20 +984,31 @@ def get_ref_edgenodes_table(femmesh, femelement_table, refedge): # 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) + refedge_fem_edgeelements = get_femelements_by_femnodes_std( + femelement_table, + refedge_nodes + ) for elem in refedge_fem_edgeelements: # { edgeID : ( nodeID, ... , nodeID )} # all nodes off this femedgeelement edge_table[elem] = femelement_table[elem] return edge_table -def get_ref_edgenodes_lengths(femnodes_mesh, edge_table): +# ************************************************************************************************ +def get_ref_edgenodes_lengths( + femnodes_mesh, + edge_table +): # calculate the appropriate node_length for every node of every mesh edge (me) # G. Lakshmi Narasaiah, Finite Element Analysis, p206ff - # [ (nodeID, length), ... , (nodeID, length) ] some nodes will have more than one entry + # [ (nodeID, length), ... , (nodeID, length) ] + # some nodes will have more than one entry if (not femnodes_mesh) or (not edge_table): - FreeCAD.Console.PrintError("Error in get_ref_edgenodes_lengths(): Empty femnodes_mesh or edge_table!\n") + FreeCAD.Console.PrintError( + "Error in get_ref_edgenodes_lengths(): " + "Empty femnodes_mesh or edge_table!\n" + ) return [] node_length_table = [] mesh_edge_length = 0 @@ -863,9 +1053,18 @@ def get_ref_edgenodes_lengths(femnodes_mesh, edge_table): # get_ref_facenodes_table # get_ref_facenodes_areas # build_mesh_faces_of_volume_elements -def get_force_obj_face_nodeload_table(femmesh, femelement_table, femnodes_mesh, frc_obj): +def get_force_obj_face_nodeload_table( + femmesh, + femelement_table, + femnodes_mesh, + frc_obj +): # force_obj_node_load_table: - # [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)] + # [ + # ('refshape_name.elemname',node_load_table), + # ..., + # ('refshape_name.elemname',node_load_table) + # ] force_obj_node_load_table = [] sum_ref_face_area = 0 sum_ref_face_node_area = 0 # for debugging @@ -893,15 +1092,18 @@ def get_force_obj_face_nodeload_table(femmesh, femelement_table, femnodes_mesh, face_table = get_ref_facenodes_table(femmesh, femelement_table, ref_face) # node_area_table: - # [ (nodeID, Area), ... , (nodeID, Area) ] some nodes will have more than one entry + # [ (nodeID, Area), ... , (nodeID, Area) ] + # some nodes will have more than one entry node_area_table = get_ref_facenodes_areas(femnodes_mesh, face_table) # node_sum_area_table: - # { nodeID : Area, ... , nodeID : Area } AreaSum for each node, one entry for each node + # { nodeID : Area, ... , nodeID : Area } + # AreaSum for each node, one entry for each node node_sum_area_table = get_ref_shape_node_sum_geom_table(node_area_table) # node_load_table: - # { nodeID : NodeLoad, ... , nodeID : NodeLoad } NodeLoad for each node, one entry for each node + # { nodeID : NodeLoad, ... , nodeID : NodeLoad } + # NodeLoad for each node, one entry for each node node_load_table = {} sum_node_areas = 0 # for debugging for node in node_sum_area_table: @@ -909,7 +1111,9 @@ def get_force_obj_face_nodeload_table(femmesh, femelement_table, femnodes_mesh, node_load_table[node] = node_sum_area_table[node] * force_per_sum_ref_face_area ratio_refface_areas = sum_node_areas / ref_face.Area if ratio_refface_areas < 0.99 or ratio_refface_areas > 1.01: - FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n') + FreeCAD.Console.PrintError( + 'Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n' + ) FreeCAD.Console.PrintMessage(' sum_node_areas: {}\n'.format(sum_node_areas)) FreeCAD.Console.PrintMessage(' ref_face_area: {}\n'.format(ref_face.Area)) sum_ref_face_node_area += sum_node_areas @@ -923,18 +1127,44 @@ def get_force_obj_face_nodeload_table(femmesh, femelement_table, femnodes_mesh, ratio = sum_node_load / frc_obj.Force if ratio < 0.99 or ratio > 1.01: - FreeCAD.Console.PrintMessage('Deviation sum_node_load to frc_obj.Force is more than 1% : {}\n'.format(ratio)) - FreeCAD.Console.PrintMessage(' sum_ref_face_node_area: {}\n'.format(sum_ref_face_node_area)) - FreeCAD.Console.PrintMessage(' sum_ref_face_area: {}\n'.format(sum_ref_face_area)) - FreeCAD.Console.PrintMessage(' sum_node_load: {}\n'.format(sum_node_load)) - FreeCAD.Console.PrintMessage(' frc_obj.Force: {}\n'.format(frc_obj.Force)) - FreeCAD.Console.PrintMessage(' the reason could be simply a circle area --> see method get_ref_face_node_areas\n') - FreeCAD.Console.PrintMessage(' the reason could also be a problem in retrieving the ref_face_node_area\n') + FreeCAD.Console.PrintMessage( + 'Deviation sum_node_load to frc_obj.Force is more than 1% : {}\n' + .format(ratio) + ) + FreeCAD.Console.PrintMessage( + ' sum_ref_face_node_area: {}\n' + .format(sum_ref_face_node_area) + ) + FreeCAD.Console.PrintMessage( + ' sum_ref_face_area: {}\n' + .format(sum_ref_face_area) + ) + FreeCAD.Console.PrintMessage( + ' sum_node_load: {}\n' + .format(sum_node_load) + ) + FreeCAD.Console.PrintMessage( + ' frc_obj.Force: {}\n' + .format(frc_obj.Force) + ) + FreeCAD.Console.PrintMessage( + ' the reason could be simply a circle area --> ' + 'see method get_ref_face_node_areas\n' + ) + FreeCAD.Console.PrintMessage( + ' the reason could also be a problem in ' + 'retrieving the ref_face_node_area\n' + ) return force_obj_node_load_table -def get_ref_facenodes_table(femmesh, femelement_table, ref_face): +# ************************************************************************************************ +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): @@ -944,15 +1174,20 @@ def get_ref_facenodes_table(femmesh, femelement_table, ref_face): '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 + # if we retrieve the nodes ourself we gone have a problem: + # 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 - ref_face_volume_elements = femmesh.getccxVolumesByFace(ref_face) # list of tupels (mv, ccx_face_nr) + # of element with elementfaces on the ref_face + # --> should work for tetra4 and tetra10 + # list of tupels (mv, ccx_face_nr) + ref_face_volume_elements = femmesh.getccxVolumesByFace(ref_face) if ref_face_volume_elements: # mesh with tetras FreeCAD.Console.PrintLog( - ' Use of getccxVolumesByFace() has returned volume elements of the ref_face.\n' + ' Use of getccxVolumesByFace() has ' + 'returned volume elements of the ref_face.\n' ) for ve in ref_face_volume_elements: veID = ve[0] @@ -969,7 +1204,10 @@ def get_ref_facenodes_table(femmesh, femelement_table, 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) + 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]: @@ -992,13 +1230,18 @@ def get_ref_facenodes_table(femmesh, femelement_table, ref_face): return face_table -def get_ref_facenodes_areas(femnodes_mesh, face_table): +# ************************************************************************************************ +def get_ref_facenodes_areas( + femnodes_mesh, + face_table +): # calculate the appropriate node_areas for every node of every mesh face (mf) # G. Lakshmi Narasaiah, Finite Element Analysis, p206ff # FIXME: only gives exact results in case of a real triangle. If for S6 or C3D10 elements # the midnodes are not on the line between the end nodes the area will not be a triangle # see http://forum.freecadweb.org/viewtopic.php?f=18&t=10939&start=40#p91355 and ff - # same applies for the quads, results are exact only if mid nodes are on the line between corner nodes + # same applies for the quads + # results are exact only if mid nodes are on the line between corner nodes # [ (nodeID,Area), ... , (nodeID,Area) ] some nodes will have more than one entry if (not femnodes_mesh) or (not face_table): @@ -1075,7 +1318,12 @@ def get_ref_facenodes_areas(femnodes_mesh, face_table): mesh_face_t2_area = get_triangle_area(P2, P5, P4) mesh_face_t3_area = get_triangle_area(P3, P6, P5) mesh_face_t4_area = get_triangle_area(P4, P5, P6) - mesh_face_area = mesh_face_t1_area + mesh_face_t2_area + mesh_face_t3_area + mesh_face_t4_area + mesh_face_area = ( + mesh_face_t1_area + + mesh_face_t2_area + + mesh_face_t3_area + + mesh_face_t4_area + ) middle_node_area = mesh_face_area / 3.0 node_area_table.append((face_table[mf][0], 0)) @@ -1114,7 +1362,14 @@ def get_ref_facenodes_areas(femnodes_mesh, face_table): mesh_face_t4_area = get_triangle_area(P7, P4, P8) mesh_face_t5_area = get_triangle_area(P5, P7, P8) mesh_face_t6_area = get_triangle_area(P5, P6, P7) - mesh_face_area = mesh_face_t1_area + mesh_face_t2_area + mesh_face_t3_area + mesh_face_t4_area + mesh_face_t5_area + mesh_face_t6_area + mesh_face_area = ( + mesh_face_t1_area + + mesh_face_t2_area + + mesh_face_t3_area + + mesh_face_t4_area + + mesh_face_t5_area + + mesh_face_t6_area + ) corner_node_area = -mesh_face_area / 12.0 middle_node_area = mesh_face_area / 3.0 @@ -1129,26 +1384,35 @@ def get_ref_facenodes_areas(femnodes_mesh, face_table): return node_area_table -def build_mesh_faces_of_volume_elements(face_table, femelement_table): +# ************************************************************************************************ +def build_mesh_faces_of_volume_elements( + face_table, + femelement_table +): # node index of facenodes in femelementtable volume element - # if we know the position of the node we can build the element face out of the unsorted face nodes + # if we know the position of the node + # we can build the element face out of the unsorted face nodes face_nodenumber_table = {} # { volumeID : ( index, ... , index ) } for veID in face_table: face_nodenumber_table[veID] = [] for n in face_table[veID]: index = femelement_table[veID].index(n) # FreeCAD.Console.PrintMessage('{}\n'.format(index)) - face_nodenumber_table[veID].append(index + 1) # local node number = index + 1 - # FreeCAD.Console.PrintMessage('VolElement: {}\n'.format(veID)) - # FreeCAD.Console.PrintMessage(' --> {}\n'.format(femelement_table[veID])) - # FreeCAD.Console.PrintMessage(' --> {}\n'.format(face_table[veID])) - # FreeCAD.Console.PrintMessage(' --> {}\n'.format(face_nodenumber_table[veID])) + # local node number = index + 1 + face_nodenumber_table[veID].append(index + 1) + FreeCAD.Console.PrintLog('VolElement: {}\n'.format(veID)) + FreeCAD.Console.PrintLog(' --> {}\n'.format(femelement_table[veID])) + FreeCAD.Console.PrintLog(' --> {}\n'.format(face_table[veID])) + FreeCAD.Console.PrintLog(' --> {}\n'.format(face_nodenumber_table[veID])) for veID in face_nodenumber_table: vol_node_ct = len(femelement_table[veID]) face_node_indexs = sorted(face_nodenumber_table[veID]) - if vol_node_ct == 10: # tetra10 --> tria6 face - if face_node_indexs == [1, 2, 3, 5, 6, 7]: # node order of face in tetra10 volume element - node_numbers = (1, 2, 3, 5, 6, 7) # node order of a tria6 face of tetra10 + # tetra10 --> tria6 face + if vol_node_ct == 10: + # node order of face in tetra10 volume element + if face_node_indexs == [1, 2, 3, 5, 6, 7]: + # node order of a tria6 face of tetra10 + node_numbers = (1, 2, 3, 5, 6, 7) elif face_node_indexs == [1, 2, 4, 5, 8, 9]: node_numbers = (1, 4, 2, 8, 9, 5) elif face_node_indexs == [1, 3, 4, 7, 8, 10]: @@ -1157,11 +1421,16 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table): 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" + "Error in build_mesh_faces_of_volume_elements(): " + "hexa20: face not found! {}\n" + .format(face_node_indexs) ) - 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 + # tetra4 --> tria3 face + elif vol_node_ct == 4: + # node order of face in tetra4 volume element + if face_node_indexs == [1, 2, 3]: + # node order of a tria3 face of tetra4 + node_numbers = (1, 2, 3) elif face_node_indexs == [1, 2, 4]: node_numbers = (1, 4, 2, 8) elif face_node_indexs == [1, 3, 4]: @@ -1170,11 +1439,15 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table): 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" + "Error in build_mesh_faces_of_volume_elements(): " + "hexa20: face not found! {}\n" + .format(face_node_indexs) ) 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 + # node order of face in hexa20 volume element + if face_node_indexs == [1, 2, 3, 4, 9, 10, 11, 12]: + # node order of a quad8 face of hexa20 + node_numbers = (1, 2, 3, 4, 9, 10, 11, 12) elif face_node_indexs == [5, 6, 7, 8, 13, 14, 15, 16]: node_numbers = (5, 8, 7, 6, 16, 15, 14, 13) elif face_node_indexs == [1, 2, 5, 6, 9, 13, 17, 18]: @@ -1187,12 +1460,16 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table): 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" + "Error in build_mesh_faces_of_volume_elements(): " + "hexa20: face not found! {}\n" + .format(face_node_indexs) ) 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 - node_numbers = (1, 2, 3, 4) # node order of a quad8 face of hexa8 + # node order of face in hexa8 volume element + if face_node_indexs == [1, 2, 3, 4]: + # node order of a quad8 face of hexa8 + node_numbers = (1, 2, 3, 4) elif face_node_indexs == [5, 6, 7, 8]: node_numbers = (5, 8, 7, 6) elif face_node_indexs == [1, 2, 5, 6]: @@ -1205,11 +1482,16 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table): 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" + "Error in build_mesh_faces_of_volume_elements(): " + "hexa20: face not found! {}\n" + .format(face_node_indexs) ) - 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 + # penta15 --> tria6 and quad8 faces + elif vol_node_ct == 15: + # node order of face in penta15 volume element + if face_node_indexs == [1, 2, 3, 7, 8, 9]: + # node order of a tria6 face of penta15 + node_numbers = (1, 2, 3, 7, 8, 9) elif face_node_indexs == [4, 5, 6, 10, 11, 12]: node_numbers = (4, 6, 5, 12, 11, 10) # tria6 elif face_node_indexs == [1, 2, 4, 5, 7, 10, 13, 14]: @@ -1220,11 +1502,16 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table): 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" + "Error in build_mesh_faces_of_volume_elements(): " + "penta15: face not found! {}\n" + .format(face_node_indexs) ) - 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 + # penta6 --> tria3 and quad4 faces + elif vol_node_ct == 6: + # node order of face in penta6 volume element + if face_node_indexs == [1, 2, 3]: + # node order of a tria3 face of penta6 + node_numbers = (1, 2, 3) elif face_node_indexs == [4, 5, 6]: node_numbers = (4, 6, 5) # tria3 elif face_node_indexs == [1, 2, 4, 5]: @@ -1235,15 +1522,21 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table): 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" + "Error in build_mesh_faces_of_volume_elements(): " + "pent6: face not found! {}\n" + .format(face_node_indexs) ) else: FreeCAD.Console.PrintError( - "Error in build_mesh_faces_of_volume_elements(): Volume not implemented: volume node count" + str(vol_node_ct) + "\n" + "Error in build_mesh_faces_of_volume_elements(): " + "Volume not implemented: volume node count {}\n" + .format(vol_node_ct) ) face_nodes = [] for i in node_numbers: - i -= 1 # node_number starts with 1, index starts with 0 --> index = node number - 1 + # node_number starts with 1 + # index starts with 0 --> + # index = node number - 1i -= 1 face_nodes.append(femelement_table[veID][i]) face_table[veID] = face_nodes # reset the entry in face_table # FreeCAD.Console.PrintMessage(' --> {}\n'.format(face_table[veID])) @@ -1251,7 +1544,9 @@ def build_mesh_faces_of_volume_elements(face_table, femelement_table): # ***** helper for Face and Edge loads *********************************************************** -def get_ref_shape_node_sum_geom_table(node_geom_table): +def get_ref_shape_node_sum_geom_table( + node_geom_table +): # shape could be Edge or Face, geom could be length or area # sum of length or area for each node of the ref_shape node_sum_geom_table = {} @@ -1266,22 +1561,35 @@ def get_ref_shape_node_sum_geom_table(node_geom_table): # ************************************************************************************************ # ***** methods for retrieving faces for constraint pressure ************************************* -def get_pressure_obj_faces(femmesh, femelement_table, femnodes_ele_table, femobj): +def get_pressure_obj_faces( + femmesh, + femelement_table, + femnodes_ele_table, + femobj +): if is_solid_femmesh(femmesh): # get the nodes - prs_face_node_set = get_femnodes_by_femobj_with_references(femmesh, femobj) # sorted and duplicates removed + # sorted and duplicates removed + prs_face_node_set = get_femnodes_by_femobj_with_references(femmesh, femobj) # FreeCAD.Console.PrintMessage('prs_face_node_set: {}\n'.format(prs_face_node_set)) # fill the bit_pattern_dict and search for the faces - bit_pattern_dict = get_bit_pattern_dict(femelement_table, femnodes_ele_table, prs_face_node_set) + bit_pattern_dict = get_bit_pattern_dict( + femelement_table, + femnodes_ele_table, + prs_face_node_set + ) pressure_faces = get_ccxelement_faces_from_binary_search(bit_pattern_dict) elif is_face_femmesh(femmesh): pressure_faces = [] - # normally we should call get_femelements_by_references and the group check should be integrated there + # normally we should call get_femelements_by_references and + # the group check should be integrated there if femmesh.GroupCount: meshfaces = get_femmesh_groupdata_sets_by_name(femmesh, femobj, 'Face') # FreeCAD.Console.PrintMessage('{}\n'.format(meshfaces)) if not meshfaces: - FreeCAD.Console.PrintError("Error: Something went wrong in getting the group element faces.\n") + FreeCAD.Console.PrintError( + "Error: Something went wrong in getting the group element faces.\n" + ) else: for mf in meshfaces: # pressure_faces.append([mf, 0]) @@ -1291,16 +1599,21 @@ def get_pressure_obj_faces(femmesh, femelement_table, femnodes_ele_table, femobj # 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" + "Pressure on shell mesh at the moment only " + "supported for meshes with appropriate group data.\n" ) return pressure_faces +# ************************************************************************************************ # depreciated method for pressure faces for constraint pressure and finite solid element mesh # why did we switch to the get_ccxelement_faces_from_binary_search? # just performance? # TODO: Find the forum topic discussion with ulrich1a and post a link -def get_pressure_obj_faces_depreciated(femmesh, femobj): +def get_pressure_obj_faces_depreciated( + femmesh, + femobj +): pressure_faces = [] for o, elem_tup in femobj['Object'].References: for elem in elem_tup: @@ -1308,13 +1621,18 @@ def get_pressure_obj_faces_depreciated(femmesh, femobj): elem_info_string = 'face load on shape: ' + o.Name + ':' + elem FreeCAD.Console.PrintMessage('{}\n'.format(elem_info_string)) if ref_shape.ShapeType == 'Face': - pressure_faces.append((elem_info_string, femmesh.getccxVolumesByFace(ref_shape))) + pressure_faces.append( + (elem_info_string, femmesh.getccxVolumesByFace(ref_shape)) + ) return pressure_faces # ************************************************************************************************ # ***** groups *********************************************************************************** -def get_mesh_group_elements(mesh_group_obj, aPart): +def get_mesh_group_elements( + mesh_group_obj, + aPart +): '''the Reference shapes of the mesh_group_object are searched in the Shape of aPart. If found in shape they are added to a dict {MeshGroupIdentifier : ['ShapeType of the Elements'], [ElementID, ElementID, ...], ...} @@ -1325,12 +1643,17 @@ def get_mesh_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' + ' Empty reference in mesh group object: {} {}\n' + .format(mesh_group_obj.Name, mesh_group_obj.Label) ) return group_elements -def get_analysis_group_elements(aAnalysis, aPart): +# ************************************************************************************************ +def get_analysis_group_elements( + aAnalysis, + aPart +): ''' all Reference shapes of all Analysis member are searched in the Shape of aPart. If found in shape they are added to a dict {ConstraintName : ['ShapeType of the Elements'], [ElementID, ElementID, ...], ...} @@ -1339,7 +1662,8 @@ def get_analysis_group_elements(aAnalysis, aPart): empty_references = [] for m in aAnalysis.Group: if hasattr(m, "References") and "ReadOnly" not in m.getEditorMode("References"): - # some C++ Constraints have a not used References Property, it is set to Hidden in ReadOnly and PropertyEditor + # some C++ Constraints have a not used References Property + # it is set to Hidden in ReadOnly and PropertyEditor if m.References: grp_ele = get_reference_group_elements(m, aPart) group_elements[grp_ele[0]] = grp_ele[1] @@ -1348,40 +1672,67 @@ def get_analysis_group_elements(aAnalysis, aPart): empty_references.append(m) if empty_references: if len(empty_references) == 1: - group_elements = get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aPart.Shape) + group_elements = get_anlysis_empty_references_group_elements( + group_elements, + aAnalysis, + aPart.Shape + ) else: - FreeCAD.Console.PrintError('Problem: more than one object with empty references.\n') - FreeCAD.Console.PrintMessage('We are going to try to get the empty material references anyway.\n') - # FemElementGeometry2D, ElementGeometry1D and FemElementFluid1D could have empty references, + FreeCAD.Console.PrintError( + 'Problem: more than one object with empty references.\n' + ) + FreeCAD.Console.PrintMessage( + 'We are going to try to get the empty material references anyway.\n' + ) + # FemElementGeometry2D, ElementGeometry1D and + # FemElementFluid1D could have empty references, # but on solid meshes only materials should have empty references for er in empty_references: FreeCAD.Console.PrintMessage(er.Name + '\n') - group_elements = get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aPart.Shape) + group_elements = get_anlysis_empty_references_group_elements( + group_elements, + aAnalysis, + aPart.Shape + ) # check if all groups have at least one element, # it doesn't mean ALL reference shapes for a group have been found 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' + 'Error: The shapes for the mesh group for the reference ' + 'shapes of analysis member: {} could not be found!\n' + .format(g) ) return group_elements -def get_reference_group_elements(obj, aPart): - ''' obj is an FEM object which has reference shapes like the group object, the material, most of the constraints - aPart is geometry feature normally CompSolid, the method searches all reference shapes of obj inside aPart even if +# ************************************************************************************************ +def get_reference_group_elements( + obj, + aPart +): + ''' obj is an FEM object which has reference shapes like the group object + the material, most of the constraints + aPart is geometry feature normally CompSolid + the method searches all reference shapes of obj inside aPart even if the reference shapes are a totally different geometry feature. a tuple is returned ('Name or Label of the FEMobject', ['Element1', 'Element2', ...]) - The names in the list are the Elements of the geometry aPart whereas 'Solid1' == aPart.Shape.Solids[0] - !!! It is strongly recommended to use as reference shapes the Solids of a CompSolid an not the Solids the CompSolid is made of !!! - see https://forum.freecadweb.org/viewtopic.php?f=18&t=12212&p=175777#p175777 and following posts - Occt might change the Solids a CompSolid is made of during creation of the CompSolid by adding Edges and vertices + The names in the list are the Elements of the geometry aPart + whereas 'Solid1' == aPart.Shape.Solids[0] + !!! It is strongly recommended to use as reference shapes the Solids of a CompSolid + and not the Solids the CompSolid is made of !!! + see https://forum.freecadweb.org/viewtopic.php?f=18&t=12212&p=175777#p175777 + and following posts + Occt might change the Solids a CompSolid is made of during + creation of the CompSolid by adding Edges and vertices Thus the Elements do not have the same geometry anymore ''' aShape = aPart.Shape if hasattr(obj, "UseLabel") and obj.UseLabel: - key = obj.Label # TODO: check the character of the Label, only allow underline and standard english character + # TODO: check the character of the Label + # only allow underline and standard english character + key = obj.Label else: key = obj.Name elements = [] @@ -1392,30 +1743,41 @@ def get_reference_group_elements(obj, aPart): # FreeCAD.Console.PrintMessage('{}\n'.format(parent)) # FreeCAD.Console.PrintMessage('{}\n'.format(childs)) for child in childs: - ref_shape = get_element(parent, child) # the method getElement(element) does not return Solid elements + # the method getElement(element) does not return Solid elements + ref_shape = get_element(parent, child) if not stype: stype = ref_shape.ShapeType elif stype != ref_shape.ShapeType: - FreeCAD.Console.PrintError('Error, two refshapes in References with different ShapeTypes.\n') - # FreeCAD.Console.PrintMessage('\n'.format(ref_shape)) + FreeCAD.Console.PrintError( + 'Error, two refshapes in References with different ShapeTypes.\n' + ) + FreeCAD.Console.PrintLog('\n'.format(ref_shape)) found_element = find_element_in_shape(aShape, ref_shape) 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' + 'Problem: For the geometry of the ' + 'following shape was no Shape found: {}\n' + .format(ref_shape) ) 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' + 'The reference Shape is not a child ' + 'nor it is the shape the mesh is made of. : {}\n' + .format(ref_shape) + ) + FreeCAD.Console.PrintMessage( + '{}--> Name of the Feature we where searching in.\n' + .format(aPart.Name) ) - FreeCAD.Console.PrintMessage(aPart.Name + '--> Name of the Feature we where searching in.\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' + '(Use the same as in the line before and you ' + 'will have less trouble :-) !!!!!!).\n' .format(parent.Name) ) # import Part @@ -1423,17 +1785,25 @@ def get_reference_group_elements(obj, aPart): # Part.show(ref_shape) else: FreeCAD.Console.PrintError('This should not happen, please debug!\n') - # in this case we would not have needed to use the is_same_geometry() inside find_element_in_shape() + # in this case we would not have needed to use the + # is_same_geometry() inside find_element_in_shape() # AFAIK we could have used the Part methods isPartner() or even isSame() # We're going to find out when we need to debug this :-)! return (key, sorted(elements)) -def get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShape): +# ************************************************************************************************ +def get_anlysis_empty_references_group_elements( + group_elements, + aAnalysis, + aShape +): '''get the elementIDs if the Reference shape is empty see get_analysis_group_elements() for more informatations - on solid meshes only material objects could have an empty reference without being something wrong! - face meshes could have empty ShellThickness and edge meshes could have empty BeamSection/FluidSection + on solid meshes only material objects could have an + empty reference without being something wrong! + face meshes could have empty ShellThickness and + edge meshes could have empty BeamSection/FluidSection ''' # FreeCAD.Console.PrintMessage('{}\n'.format(group_elements)) material_ref_shapes = [] @@ -1459,7 +1829,8 @@ def get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShap 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' + '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) @@ -1495,7 +1866,11 @@ def get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShap return group_elements -def find_element_in_shape(aShape, anElement): +# ************************************************************************************************ +def find_element_in_shape( + aShape, + anElement +): # import Part ele_st = anElement.ShapeType if ele_st == 'Solid' or ele_st == 'CompSolid': @@ -1546,7 +1921,11 @@ def find_element_in_shape(aShape, anElement): FreeCAD.Console.PrintError('Compound is not supported.\n') -def get_vertexes_by_element(aShape, anElement): +# ************************************************************************************************ +def get_vertexes_by_element( + aShape, + anElement +): # we're going to extend the method find_element_in_shape and return the vertexes # import Part ele_vertexes = [] @@ -1591,7 +1970,11 @@ def get_vertexes_by_element(aShape, anElement): FreeCAD.Console.PrintError('Compound is not supported.\n') -def is_same_geometry(shape1, shape2): +# ************************************************************************************************ +def is_same_geometry( + shape1, + shape2 +): # the vertexes and the CenterOfMass are compared # it is a hack, but I do not know any better ! # check of Volume and Area before starting with the vertices could be added @@ -1626,7 +2009,11 @@ def is_same_geometry(shape1, shape2): return False -def get_element(part, element): +# ************************************************************************************************ +def get_element( + part, + element +): if element.startswith('Solid'): index = int(element.lstrip('Solid')) - 1 if index >= len(part.Shape.Solids): @@ -1640,7 +2027,11 @@ def get_element(part, element): return part.Shape.getElement(element) # Face, Edge, Vertex -def femelements_count_ok(len_femelement_table, count_femelements): +# ************************************************************************************************ +def femelements_count_ok( + len_femelement_table, + count_femelements +): FreeCAD.Console.PrintMessage( 'Count finite elements as sum of constraints: {}\n' .format(count_femelements) @@ -1658,7 +2049,10 @@ def femelements_count_ok(len_femelement_table, count_femelements): return False -def delete_duplicate_mesh_elements(refelement_table): +# ************************************************************************************************ +def delete_duplicate_mesh_elements( + refelement_table +): new_refelement_table = {} # duplicates deleted for elem, nodes in refelement_table.items(): if sorted(nodes) not in sortlistoflistvalues(new_refelement_table.values()): @@ -1666,7 +2060,12 @@ def delete_duplicate_mesh_elements(refelement_table): return new_refelement_table -def get_triangle_area(P1, P2, P3): +# ************************************************************************************************ +def get_triangle_area( + P1, + P2, + P3 +): # import Part # W = Part.Wire([Part.makeLine(P1,P2), Part.makeLine(P2,P3), Part.makeLine(P3,P1)]) # Part.show(Part.Face(W)) @@ -1676,38 +2075,56 @@ def get_triangle_area(P1, P2, P3): return 0.5 * vec3.Length -def sortlistoflistvalues(listoflists): +# ************************************************************************************************ +def sortlistoflistvalues( + listoflists +): new_list = [] for l in listoflists: new_list.append(sorted(l)) return new_list -def is_solid_femmesh(femmesh): +# ************************************************************************************************ +def is_solid_femmesh( + femmesh +): # solid femmesh if femmesh.VolumeCount > 0: return True -def has_no_face_data(femmesh): +# ************************************************************************************************ +def has_no_face_data( + femmesh +): # 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): +# ************************************************************************************************ +def is_face_femmesh( + femmesh +): # face femmesh if femmesh.VolumeCount == 0 and femmesh.FaceCount > 0: return True -def is_edge_femmesh(femmesh): +# ************************************************************************************************ +def is_edge_femmesh( + femmesh +): # edge femmesh if femmesh.VolumeCount == 0 and femmesh.FaceCount == 0 and femmesh.EdgeCount > 0: return True -def is_zplane_2D_mesh(femmesh): +# ************************************************************************************************ +def is_zplane_2D_mesh( + femmesh +): # used in oofem writer to distinguish between 3D and 2D plane stress if is_face_femmesh(femmesh) is True: tol = 0.0001 @@ -1720,7 +2137,10 @@ def is_zplane_2D_mesh(femmesh): return False -def get_three_non_colinear_nodes(nodes_coords): +# ************************************************************************************************ +def get_three_non_colinear_nodes( + nodes_coords +): # Code to obtain three non-colinear nodes on the PlaneRotation support face # nodes_coords --> [(nodenumber, x, y, z), (nodenumber, x, y, z), ...] if not nodes_coords: @@ -1764,7 +2184,10 @@ def get_three_non_colinear_nodes(nodes_coords): return [node_1, node_2, node_3] -def get_rectangular_coords(obj): +# ************************************************************************************************ +def get_rectangular_coords( + obj +): from math import cos, sin, radians A = [1, 0, 0] B = [0, 1, 0] @@ -1811,7 +2234,10 @@ def get_rectangular_coords(obj): return coords -def get_cylindrical_coords(obj): +# ************************************************************************************************ +def get_cylindrical_coords( + obj +): vec = obj.Axis base = obj.BasePoint Ax = base[0] + 10 * vec[0] @@ -1828,7 +2254,10 @@ def get_cylindrical_coords(obj): return coords -def write_D_network_element_to_inputfile(fileName): +# ************************************************************************************************ +def write_D_network_element_to_inputfile( + fileName +): # replace B32 elements with D elements for fluid section f = open(fileName, 'r+') lines = f.readlines() @@ -1843,7 +2272,12 @@ def write_D_network_element_to_inputfile(fileName): f.close() -def use_correct_fluidinout_ele_def(FluidInletoutlet_ele, fileName, fluid_inout_nodes_file): +# ************************************************************************************************ +def use_correct_fluidinout_ele_def( + FluidInletoutlet_ele, + fileName, + fluid_inout_nodes_file +): f = open(fileName, 'r') cnt = 0 line = f.readline() @@ -1872,7 +2306,10 @@ def use_correct_fluidinout_ele_def(FluidInletoutlet_ele, fileName, fluid_inout_n f.seek(0) cnt = 0 elem_counter = 0 - FreeCAD.Console.PrintMessage('1DFlow inout nodes file: ' + fluid_inout_nodes_file + '\n') + FreeCAD.Console.PrintMessage( + '1DFlow inout nodes file: {}\n' + .format(fluid_inout_nodes_file) + ) inout_nodes_file = open(fluid_inout_nodes_file, "w") for line in lines: new_line = '' @@ -1887,9 +2324,12 @@ 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('{},{},{},{}\n'.format( + node1, + node2, + node3, + FluidInletoutlet_ele[i][1] + )) elif j == 3: new_line = new_line + a[j] else: @@ -1900,9 +2340,12 @@ 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('{},{},{},{}\n'.format( + node1, + node2, + node3, + FluidInletoutlet_ele[i][1] + )) else: new_line = new_line + a[j] + ',' if new_line == '': @@ -1915,7 +2358,10 @@ def use_correct_fluidinout_ele_def(FluidInletoutlet_ele, fileName, fluid_inout_n inout_nodes_file.close() -def compact_mesh(old_femmesh): +# ************************************************************************************************ +def compact_mesh( + old_femmesh +): ''' removes all gaps in node and element ids, start ids with 1 returns a tuple (FemMesh, node_assignment_map, element_assignment_map)