FEM: constraint plane rotation: add implementation for solver CalculiX

src/Mod/Fem/FemInputWriter.py

@@ -45,7 +45,7 @@ def __init__(self,
                  analysis_obj, solver_obj,
                  mesh_obj, mat_obj,
                  fixed_obj, displacement_obj,
-                 contact_obj,
+                 contact_obj, planerotation_obj,
                  selfweight_obj, force_obj, pressure_obj,
                  beamsection_obj, shellthickness_obj,
                  analysis_type, eigenmode_parameters,
@@ -58,6 +58,7 @@ def __init__(self,
         self.fixed_objects = fixed_obj
         self.displacement_objects = displacement_obj
         self.contact_objects = contact_obj
+        self.planerotation_objects = planerotation_obj
         self.selfweight_objects = selfweight_obj
         self.force_objects = force_obj
         self.pressure_objects = pressure_obj
@@ -80,16 +81,28 @@ def __init__(self,
         self.femmesh = self.mesh_object.FemMesh
         self.femnodes_mesh = {}
         self.femelement_table = {}
+        self.constraint_conflict_nodes = []
 
     def get_constraints_fixed_nodes(self):
         # get nodes
         for femobj in self.fixed_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
             femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References)
+            # add nodes to constraint_conflict_nodes, needed by constraint plane rotation
+            for node in femobj['Nodes']:
+                self.constraint_conflict_nodes.append(node)
 
     def get_constraints_displacement_nodes(self):
         # get nodes
         for femobj in self.displacement_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
             femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References)
+            # add nodes to constraint_conflict_nodes, needed by constraint plane rotation
+            for node in femobj['Nodes']:
+                self.constraint_conflict_nodes.append(node)
+
+    def get_constraints_planerotation_nodes(self):
+        # get nodes
+        for femobj in self.planerotation_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References)
 
     def get_constraints_force_nodeloads(self):
         # check shape type of reference shape

src/Mod/Fem/FemInputWriterCcx.py

@@ -40,7 +40,7 @@ def __init__(self,
                  analysis_obj, solver_obj,
                  mesh_obj, mat_obj,
                  fixed_obj, displacement_obj,
-                 contact_obj,
+                 contact_obj, planerotation_obj,
                  selfweight_obj, force_obj, pressure_obj,
                  beamsection_obj, shellthickness_obj,
                  analysis_type=None, eigenmode_parameters=None,
@@ -50,7 +50,7 @@ def __init__(self,
                                                analysis_obj, solver_obj,
                                                mesh_obj, mat_obj,
                                                fixed_obj, displacement_obj,
-                                               contact_obj,
+                                               contact_obj, planerotation_obj,
                                                selfweight_obj, force_obj, pressure_obj,
                                                beamsection_obj, shellthickness_obj,
                                                analysis_type, eigenmode_parameters,
@@ -71,10 +71,14 @@ def write_calculix_input_file(self):
             self.write_node_sets_constraints_fixed(inpfile)
         if self.displacement_objects:
             self.write_node_sets_constraints_displacement(inpfile)
+        if self.planerotation_objects:
+            self.write_node_sets_constraints_planerotation(inpfile)
         if self.contact_objects:
             self.write_surfaces_contraints_contact(inpfile)
         self.write_materials(inpfile)
         self.write_femelementsets(inpfile)
+        if self.planerotation_objects:
+            self.write_constraints_planerotation(inpfile)
         if self.contact_objects:
             self.write_constraints_contact(inpfile)
         self.write_step_begin(inpfile)
@@ -93,6 +97,7 @@ def write_calculix_input_file(self):
             self.write_analysis_frequency(inpfile)
         self.write_outputs_types(inpfile)
         self.write_step_end(inpfile)
+
         self.write_footer(inpfile)
         inpfile.close()
         return self.file_name
@@ -158,6 +163,43 @@ def write_node_sets_constraints_displacement(self, f):
             for n in femobj['Nodes']:
                 f.write(str(n) + ',\n')
 
+    def write_node_sets_constraints_planerotation(self, f):
+        # get nodes
+        self.get_constraints_planerotation_nodes()
+        # write nodes to file
+        if not self.femnodes_mesh:
+            self.femnodes_mesh = self.femmesh.Nodes
+        f.write('\n***********************************************************\n')
+        f.write('** Node set for plane rotation constraint\n')
+        f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
+        # info about self.constraint_conflict_nodes:
+        # is used to check if MPC and constraint fixed and constraint displacement share same nodes,
+        # because MPC's and constriants fixed an constraints displacement can't share same nodes.
+        # thus call write_node_sets_constraints_planerotation has to be after constraint fixed and constraint displacement
+        for femobj in self.planerotation_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            l_nodes = femobj['Nodes']
+            fric_obj = femobj['Object']
+            f.write('*NSET,NSET=' + fric_obj.Name + '\n')
+            # Code to extract nodes and coordinates on the PlaneRotation support face
+            nodes_coords = []
+            for node in l_nodes:
+                nodes_coords.append((node, self.femnodes_mesh[node].x, self.femnodes_mesh[node].y, self.femnodes_mesh[node].z))
+            node_planerotation = FemMeshTools.get_three_non_colinear_nodes(nodes_coords)
+            for i in range(len(l_nodes)):
+                if l_nodes[i] not in node_planerotation:
+                    node_planerotation.append(l_nodes[i])
+            MPC_nodes = []
+            for i in range(len(node_planerotation)):
+                cnt = 0
+                for j in range(len(self.constraint_conflict_nodes)):
+                    if node_planerotation[i] == self.constraint_conflict_nodes[j]:
+                        cnt = cnt + 1
+                if cnt == 0:
+                    MPC = node_planerotation[i]
+                    MPC_nodes.append(MPC)
+            for i in range(len(MPC_nodes)):
+                f.write(str(MPC_nodes[i]) + ',\n')
+
     def write_surfaces_contraints_contact(self, f):
         # get surface nodes and write them to file
         f.write('\n***********************************************************\n')
@@ -338,6 +380,16 @@ def write_constraints_contact(self, f):
                 stick = (slope / 10.0)
                 f.write(str(friction) + ', ' + str(stick) + ' \n')
 
+    def write_constraints_planerotation(self, f):
+        f.write('\n***********************************************************\n')
+        f.write('** PlaneRotation Constraints\n')
+        f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
+        for femobj in self.planerotation_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            fric_obj_name = femobj['Object'].Name
+            f.write('*MPC\n')
+            f.write('PLANE,' + fric_obj_name + '\n')
+            f.write('\n')
+
     def write_constraints_selfweight(self, f):
         f.write('\n***********************************************************\n')
         f.write('** Self weight\n')

src/Mod/Fem/FemInputWriterZ88.py

@@ -37,7 +37,7 @@ def __init__(self,
                  analysis_obj, solver_obj,
                  mesh_obj, mat_obj,
                  fixed_obj, displacement_obj,
-                 contact_obj,
+                 contact_obj, planerotation_obj,
                  selfweight_obj, force_obj, pressure_obj,
                  beamsection_obj, shellthickness_obj,
                  analysis_type=None, eigenmode_parameters=None,
@@ -47,7 +47,7 @@ def __init__(self,
                                                analysis_obj, solver_obj,
                                                mesh_obj, mat_obj,
                                                fixed_obj, displacement_obj,
-                                               contact_obj,
+                                               contact_obj, planerotation_obj,
                                                selfweight_obj, force_obj, pressure_obj,
                                                beamsection_obj, shellthickness_obj,
                                                analysis_type, eigenmode_parameters,

src/Mod/Fem/FemMeshTools.py

@@ -578,6 +578,50 @@ def is_edge_femmesh(femmesh):
         return True
 
 
+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:
+        print(len(nodes_coords))
+        print('Error: No nodes in nodes_coords')
+        return []
+    dum_max = [1, 2, 3, 4, 5, 6, 7, 8, 0]
+    for i in range(len(nodes_coords)):
+        for j in range(len(nodes_coords) - 1 - i):
+            x_1 = nodes_coords[j][1]
+            x_2 = nodes_coords[j + 1][1]
+            y_1 = nodes_coords[j][2]
+            y_2 = nodes_coords[j + 1][2]
+            z_1 = nodes_coords[j][3]
+            z_2 = nodes_coords[j + 1][3]
+            node_1 = nodes_coords[j][0]
+            node_2 = nodes_coords[j + 1][0]
+            distance = ((x_1 - x_2) ** 2 + (y_1 - y_2) ** 2 + (z_1 - z_2) ** 2) ** 0.5
+            if distance > dum_max[8]:
+                dum_max = [node_1, x_1, y_1, z_1, node_2, x_2, y_2, z_2, distance]
+    node_dis = [1, 0]
+    for i in range(len(nodes_coords)):
+        x_1 = dum_max[1]
+        x_2 = dum_max[5]
+        x_3 = nodes_coords[i][1]
+        y_1 = dum_max[2]
+        y_2 = dum_max[6]
+        y_3 = nodes_coords[i][2]
+        z_1 = dum_max[3]
+        z_2 = dum_max[7]
+        z_3 = nodes_coords[i][3]
+        node_3 = int(nodes_coords[j][0])
+        distance_1 = ((x_1 - x_3) ** 2 + (y_1 - y_3) ** 2 + (z_1 - z_3) ** 2) ** 0.5
+        distance_2 = ((x_3 - x_2) ** 2 + (y_3 - y_2) ** 2 + (z_3 - z_2) ** 2) ** 0.5
+        tot = distance_1 + distance_2
+        if tot > node_dis[1]:
+            node_dis = [node_3, tot]
+    node_1 = int(dum_max[0])
+    node_2 = int(dum_max[4])
+    print([node_1, node_2, node_3])
+    return [node_1, node_2, node_3]
+
+
 def make_femmesh(mesh_data):
     ''' makes an FreeCAD FEM Mesh object from FEM Mesh data
     '''

src/Mod/Fem/FemTools.py

@@ -190,6 +190,10 @@ def update_objects(self):
         # set of displacements for the analysis. Updated with update_objects
         # Individual displacement_constraints are Proxy.Type "FemConstraintDisplacement"
         self.displacement_constraints = []
+        ## @var planerotation_constraints
+        #  set of plane rotation constraints from the analysis. Updated with update_objects
+        #  Individual constraints are "Fem::ConstraintPlaneRotation" type
+        self.planerotation_constraints = []
         ## @var contact_constraints
         #  set of contact constraints from the analysis. Updated with update_objects
         #  Individual constraints are "Fem::ConstraintContact" type
@@ -240,6 +244,10 @@ def update_objects(self):
                 displacement_constraint_dict = {}
                 displacement_constraint_dict['Object'] = m
                 self.displacement_constraints.append(displacement_constraint_dict)
+            elif m.isDerivedFrom("Fem::ConstraintPlaneRotation"):
+                planerotation_constraint_dict = {}
+                planerotation_constraint_dict['Object'] = m
+                self.planerotation_constraints.append(planerotation_constraint_dict)
             elif m.isDerivedFrom("Fem::ConstraintContact"):
                 contact_constraint_dict = {}
                 contact_constraint_dict['Object'] = m

src/Mod/Fem/FemToolsCcx.py

@@ -92,7 +92,7 @@ def write_inp_file(self):
                                               self.analysis, self.solver,
                                               self.mesh, self.materials,
                                               self.fixed_constraints, self.displacement_constraints,
-                                              self.contact_constraints,
+                                              self.contact_constraints, self.planerotation_constraints,
                                               self.selfweight_constraints, self.force_constraints, self.pressure_constraints,
                                               self.beam_sections, self.shell_thicknesses,
                                               self.analysis_type, self.eigenmode_parameters,

src/Mod/Fem/FemToolsZ88.py

@@ -84,7 +84,7 @@ def write_inp_file(self):
                                               self.analysis, self.solver,
                                               self.mesh, self.materials,
                                               self.fixed_constraints, self.displacement_constraints,
-                                              self.contact_constraints,
+                                              self.contact_constraints, self.planerotation_constraints,
                                               self.selfweight_constraints, self.force_constraints, self.pressure_constraints,
                                               self.beam_sections, self.shell_thicknesses,
                                               self.analysis_type, None,