From eb82923869b77697314c4f635a964623bf409c0e Mon Sep 17 00:00:00 2001 From: marioalexis Date: Wed, 7 May 2025 18:47:09 -0300 Subject: [PATCH] Fem: Search elements for electrostatic constraints --- src/Mod/Fem/femmesh/meshsetsgetter.py | 73 +++++++++++++ src/Mod/Fem/femmesh/meshtools.py | 152 ++++++++++++++++++++++++++ 2 files changed, 225 insertions(+) diff --git a/src/Mod/Fem/femmesh/meshsetsgetter.py b/src/Mod/Fem/femmesh/meshsetsgetter.py index 0f9b62fa9e..1b267eb540 100644 --- a/src/Mod/Fem/femmesh/meshsetsgetter.py +++ b/src/Mod/Fem/femmesh/meshsetsgetter.py @@ -144,11 +144,14 @@ class MeshSetsGetter: self.get_constraints_sectionprint_faces() self.get_constraints_transform_nodes() self.get_constraints_temperature_nodes() + self.get_constraints_electrostatic_nodes() # constraints sets with constraint data self.get_constraints_force_nodeloads() self.get_constraints_pressure_faces() self.get_constraints_heatflux_faces() + self.get_constraints_electrostatic_faces() + self.get_constraints_electricchargedensity_faces() setstime = round((time.process_time() - time_start), 3) FreeCAD.Console.PrintMessage(f"Getting mesh data time: {setstime} seconds.\n") @@ -251,6 +254,18 @@ class MeshSetsGetter: print_obj_info(femobj["Object"]) femobj["Nodes"] = meshtools.get_femnodes_by_femobj_with_references(self.femmesh, femobj) + def get_constraints_electrostatic_nodes(self): + if not self.member.cons_electrostatic: + return + # get nodes + for femobj in self.member.cons_electrostatic: + # femobj --> dict, FreeCAD document object is femobj["Object"] + if femobj["Object"].BoundaryCondition == "Dirichlet": + print_obj_info(femobj["Object"]) + femobj["Nodes"] = meshtools.get_femnodes_by_femobj_with_references( + self.femmesh, femobj + ) + def get_constraints_force_nodeloads(self): if not self.member.cons_force: return @@ -355,6 +370,64 @@ class MeshSetsGetter: femobj["PressureFaces"] = [(some_string, pressure_faces)] FreeCAD.Console.PrintLog("{}\n".format(femobj["PressureFaces"])) + def get_constraints_electrostatic_faces(self): + if not self.member.cons_electrostatic: + return + if not self.femnodes_mesh: + self.femnodes_mesh = self.femmesh.Nodes + if not self.femelement_table: + self.femelement_table = meshtools.get_femelement_table(self.femmesh) + if not self.femnodes_ele_table: + self.femnodes_ele_table = meshtools.get_femnodes_ele_table( + self.femnodes_mesh, self.femelement_table + ) + + for femobj in self.member.cons_electrostatic: + if femobj["Object"].BoundaryCondition == "Neumann": + print_obj_info(femobj["Object"]) + + charged_faces = meshtools.get_charge_density_obj_faces( + self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj + ) + some_string = "{}: face electric flux".format(femobj["Object"].Name) + femobj["ElectricFluxFaces"] = [(some_string, charged_faces)] + FreeCAD.Console.PrintLog("{}\n".format(femobj["ElectricFluxFaces"])) + + def get_constraints_electricchargedensity_faces(self): + if not self.member.cons_electricchargedensity: + return + if not self.femnodes_mesh: + self.femnodes_mesh = self.femmesh.Nodes + if not self.femelement_table: + self.femelement_table = meshtools.get_femelement_table(self.femmesh) + if not self.femnodes_ele_table: + self.femnodes_ele_table = meshtools.get_femnodes_ele_table( + self.femnodes_mesh, self.femelement_table + ) + + for femobj in self.member.cons_electricchargedensity: + if femobj["Object"].Mode in ["Interface", "Total Interface"]: + print_obj_info(femobj["Object"]) + + charged_faces = meshtools.get_charge_density_obj_faces( + self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj + ) + some_string = "{}: face electric charge density".format(femobj["Object"].Name) + femobj["ChargeDensityFaces"] = [(some_string, charged_faces)] + FreeCAD.Console.PrintLog("{}\n".format(femobj["ChargeDensityFaces"])) + + elif femobj["Object"].Mode in ["Source", "Total Source"]: + print_obj_info(femobj["Object"]) + + charged_volumes = meshtools.get_charge_density_obj_elements( + self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj + ) + some_string = "{}: Elements with electric charge density".format( + femobj["Object"].Name + ) + femobj["ChargeDensityElements"] = (some_string, charged_volumes) + FreeCAD.Console.PrintLog("{}\n".format(femobj["ChargeDensityElements"])) + def get_constraints_contact_faces(self): if not self.member.cons_contact: return diff --git a/src/Mod/Fem/femmesh/meshtools.py b/src/Mod/Fem/femmesh/meshtools.py index 2f25bf3b43..0c7919e972 100644 --- a/src/Mod/Fem/femmesh/meshtools.py +++ b/src/Mod/Fem/femmesh/meshtools.py @@ -224,6 +224,7 @@ def get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_set): The number in the ele_dict is organized as a bit array. The corresponding bit is set, if the node of the node_set is contained in the element. """ + # print("BIT PATTERN", 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.PrintLog(f"node_set: {node_set}\n") @@ -241,6 +242,78 @@ def get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_set): # ************************************************************************************************ +def get_ccxelement_volumes_elements_from_binary_search(bit_pattern_dict): + tet10_mask = {0b1111111111: 1} + tet4_mask = {0b1111: 1} + hex8_mask = {0b11111111: 1} + hex20_mask = {0b11111111111111111111: 1} + pent6_mask = {0b111111: 1} + pent15_mask = {0b111111111111111: 1} + vol_dict = { + 4: tet4_mask, + 6: pent6_mask, + 8: hex8_mask, + 10: tet10_mask, + 15: pent15_mask, + 20: hex20_mask, + } + volumes = [] + for ele in bit_pattern_dict: + mask_dict = vol_dict[bit_pattern_dict[ele][0]] + for key in mask_dict: + if (key & bit_pattern_dict[ele][1]) == key: + volumes.append(ele) + # print("VOLUMES:", volumes) + FreeCAD.Console.PrintLog(f"found Volumes: {len(volumes)}\n") + # FreeCAD.Console.PrintMessage("faces: {}\n".format(faces)) + return volumes + + +def get_ccxelement_faces_elements_from_binary_search(bit_pattern_dict): + tria3_mask = {0b111: 1} + tria6_mask = {0b111111: 1} + quad4_mask = {0b1111: 1} + quad8_mask = {0b11111111: 1} + vol_dict = { + 3: tria3_mask, + 6: tria6_mask, + 4: quad4_mask, + 8: quad8_mask, + } + faces = [] + for ele in bit_pattern_dict: + mask_dict = vol_dict[bit_pattern_dict[ele][0]] + for key in mask_dict: + if (key & bit_pattern_dict[ele][1]) == key: + faces.append(ele) + # print("CARAS:", faces) + FreeCAD.Console.PrintMessage(f"found Edges: {len(faces)}\n") + return faces + + +def get_ccxelement_edges_from_binary_search(bit_pattern_dict): + tria3_mask = {0b011: 1, 0b110: 2, 0b101: 3} + tria6_mask = {0b001011: 1, 0b010110: 2, 0b100101: 3} + quad4_mask = {0b0011: 1, 0b0110: 2, 0b1100: 3, 0b1001: 4} + quad8_mask = {0b00010011: 1, 0b00100110: 2, 0b01001100: 3, 0b10001001: 4} + vol_dict = { + 3: tria3_mask, + 6: tria6_mask, + 4: quad4_mask, + 8: quad8_mask, + } + faces = [] + for ele in bit_pattern_dict: + mask_dict = vol_dict[bit_pattern_dict[ele][0]] + for key in mask_dict: + if (key & bit_pattern_dict[ele][1]) == key: + faces.append([ele, mask_dict[key]]) + # print("EDGES:", faces) + FreeCAD.Console.PrintMessage(f"found Edges: {len(faces)}\n") + + return faces + + def get_ccxelement_faces_from_binary_search(bit_pattern_dict): """get the CalculiX element face numbers""" # the forum topic discussion with ulrich1a and others ... Better mesh last instead of mesh first @@ -268,6 +341,7 @@ def get_ccxelement_faces_from_binary_search(bit_pattern_dict): for key in mask_dict: if (key & bit_pattern_dict[ele][1]) == key: faces.append([ele, mask_dict[key]]) + # print("FACES:", faces) FreeCAD.Console.PrintLog(f"found Faces: {len(faces)}\n") # FreeCAD.Console.PrintMessage("faces: {}\n".format(faces)) return faces @@ -1427,6 +1501,84 @@ def get_ref_shape_node_sum_geom_table(node_geom_table): # ************************************************************************************************ # ***** methods for retrieving element face sets ************************************************* +# ***** charged faces **************************************************************************** +def get_charge_density_obj_elements(femmesh, femelement_table, femnodes_ele_table, femobj): + node_set = [] + res = [] + # TODO get elements from mesh groups + # if femmesh.GroupCount: + # node_set = get_femmesh_groupdata_sets_by_name(femmesh, femobj, "Node") + # # FreeCAD.Console.PrintMessage("node_set_group: {}\n".format(node_set)) + # if node_set: + # FreeCAD.Console.PrintLog( + # " Finite element mesh nodes where retrieved " + # "from existent finite element mesh group data.\n" + # ) + if not node_set: + FreeCAD.Console.PrintLog( + " Finite element mesh nodes will be retrieved " + "by searching the appropriate nodes in the finite element mesh.\n" + ) + for feat, ref in femobj["Object"].References: + for sub_ref in ref: + sub = (feat, (sub_ref,)) + node_set = get_femnodes_by_references(femmesh, [sub]) + charged_volume_node_set = sorted(set(node_set)) + + bit_pattern_dict = get_bit_pattern_dict( + femelement_table, femnodes_ele_table, charged_volume_node_set + ) + sh = feat.getSubObject(sub_ref) + if sh.ShapeType == "Solid": + charged_elem = get_ccxelement_volumes_elements_from_binary_search( + bit_pattern_dict + ) + elif sh.ShapeType == "Face": + charged_elem = get_ccxelement_faces_elements_from_binary_search( + bit_pattern_dict + ) + res.append((sub, charged_elem)) + + return res + + +def get_charge_density_obj_faces(femmesh, femelement_table, femnodes_ele_table, femobj): + node_set = [] + res = [] + # TODO get elements from mesh groups + # if femmesh.GroupCount: + # node_set = get_femmesh_groupdata_sets_by_name(femmesh, femobj, "Node") + # # FreeCAD.Console.PrintMessage("node_set_group: {}\n".format(node_set)) + # if node_set: + # FreeCAD.Console.PrintLog( + # " Finite element mesh nodes where retrieved " + # "from existent finite element mesh group data.\n" + # ) + if not node_set: + FreeCAD.Console.PrintLog( + " Finite element mesh nodes will be retrieved " + "by searching the appropriate nodes in the finite element mesh.\n" + ) + for feat, ref in femobj["Object"].References: + for sub_ref in ref: + sub = (feat, (sub_ref,)) + node_set = get_femnodes_by_references(femmesh, [sub]) + charged_face_node_set = sorted(set(node_set)) + + bit_pattern_dict = get_bit_pattern_dict( + femelement_table, femnodes_ele_table, charged_face_node_set + ) + sh = feat.getSubObject(sub_ref) + if sh.ShapeType == "Face": + charged_faces = get_ccxelement_faces_from_binary_search(bit_pattern_dict) + elif sh.ShapeType == "Edge": + charged_faces = get_ccxelement_edges_from_binary_search(bit_pattern_dict) + + res.append((sub, charged_faces)) + + return res + + # ***** pressure faces *************************************************************************** def get_pressure_obj_faces(femmesh, femelement_table, femnodes_ele_table, femobj): # see get_ccxelement_faces_from_binary_search for more information