Fem: usar DFLUX en electrostatic

Author: marioalexis84

src/Mod/Fem/femmesh/meshsetsgetter.py

@@ -150,6 +150,7 @@ def get_mesh_sets(self):
         self.get_constraints_force_nodeloads()
         self.get_constraints_pressure_faces()
         self.get_constraints_heatflux_faces()
+        self.get_constraints_electrostatic_faces()
 
         setstime = round((time.process_time() - time_start), 3)
         FreeCAD.Console.PrintMessage(f"Getting mesh data time: {setstime} seconds.\n")
@@ -258,8 +259,11 @@ def get_constraints_electrostatic_nodes(self):
         # get nodes
         for femobj in self.member.cons_electrostatic:
             # femobj --> dict, FreeCAD document object is femobj["Object"]
-            print_obj_info(femobj["Object"])
-            femobj["Nodes"] = meshtools.get_femnodes_by_femobj_with_references(self.femmesh, femobj)
+            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:
@@ -365,6 +369,34 @@ def get_constraints_pressure_faces(self):
             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:
+            # femobj --> dict, FreeCAD document object is femobj["Object"]
+            if femobj["Object"].BoundaryCondition == "Neumann":
+                print_obj_info(femobj["Object"])
+
+                pressure_faces = meshtools.get_pressure_obj_faces(
+                    self.femmesh, self.femelement_table, self.femnodes_ele_table, femobj
+                )
+            # the data model is for compatibility reason with deprecated version
+            # get_pressure_obj_faces_depreciated returns the face ids in a tuple per ref_shape
+            # some_string was the reference_shape_element_string in deprecated method
+            # [(some_string, [ele_id, ele_face_id], [ele_id, ele_face_id], ...])]
+            some_string = "{}: face electric flux".format(femobj["Object"].Name)
+            femobj["ElectricFluxFaces"] = [(some_string, pressure_faces)]
+            FreeCAD.Console.PrintLog("{}\n".format(femobj["ElectricFluxFaces"]))
+
     def get_constraints_contact_faces(self):
         if not self.member.cons_contact:
             return

src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py

@@ -41,9 +41,11 @@ def get_constraint_title():
 
 
 def write_meshdata_constraint(f, femobj, pot_obj, ccxwriter):
-    f.write(f"*NSET,NSET={pot_obj.Name}\n")
-    for n in femobj["Nodes"]:
-        f.write(f"{n},\n")
+
+    if femobj["Object"].BoundaryCondition == "Dirichlet":
+        f.write(f"*NSET,NSET={pot_obj.Name}\n")
+        for n in femobj["Nodes"]:
+            f.write(f"{n},\n")
 
 
 def get_before_write_meshdata_constraint():
@@ -66,18 +68,37 @@ def write_constraint(f, femobj, pot_obj, ccxwriter):
 
     # floats read from ccx should use {:.13G}, see comment in writer module
 
-    NumberOfNodes = len(femobj["Nodes"])
     if pot_obj.BoundaryCondition == "Dirichlet":
         f.write("*BOUNDARY\n")
         f.write("{},11,11,{:.13G}\n".format(pot_obj.Name, pot_obj.Potential.getValueAs("mV").Value))
         f.write("\n")
     elif pot_obj.BoundaryCondition == "Neumann":
-        f.write("*CFLUX\n")
-        # CFLUX has to be specified in mW
-        f.write(
-            "{},11,{:.13G}\n".format(
-                pot_obj.Name,
-                pot_obj.SurfaceChargeDensity.getValueAs("A*s/mm^2").Value / NumberOfNodes,
-            )
-        )
+        charge_density = pot_obj.SurfaceChargeDensity.getValueAs("A*s/mm^2").Value
+
+        f.write("*DFLUX\n")
+        for ref_shape in femobj["ElectricFluxFaces"]:
+            # the loop is needed for compatibility reason
+            # in deprecated method get_pressure_obj_faces_depreciated
+            # the face ids where per ref_shape
+            f.write("** " + ref_shape[0] + "\n")
+            for face, fno in ref_shape[1]:
+                if fno > 0:  # solid mesh face
+                    f.write(f"{face},S{fno},{charge_density}\n")
+        #                # on shell mesh face: fno == 0
+        #                # normal of element face == face normal
+        #                elif fno == 0:
+        #                    f.write(f"{face},S,{press_rev}\n")
+        #                # on shell mesh face: fno == -1
+        #                # normal of element face opposite direction face normal
+        #                elif fno == -1:
+        #                    f.write(f"{face},P,{-1 * press_rev}\n")
+
+        #        f.write("*CFLUX\n")
+        #        # CFLUX has to be specified in mW
+        #        f.write(
+        #            "{},11,{:.13G}\n".format(
+        #                pot_obj.Name,
+        #                pot_obj.SurfaceChargeDensity.getValueAs("A*s/mm^2").Value,
+        #            )
+        #        )
         f.write("\n")

src/Mod/Fem/femsolver/calculix/write_femelement_material.py

@@ -93,9 +93,9 @@ def is_density_needed():
                 KV_in_mm2s = KV.getValueAs("mm^2/s").Value
                 DV_in_tmms = KV_in_mm2s * density_in_tonne_per_mm3
         if ccxwriter.analysis_type == "electrostatic":
-            TC = FreeCAD.Units.Quantity(mat_obj.Material["ThermalConductivity"])
-            TC_in_WmK = TC.getValueAs("W/m/K").Value
-
+            rel_perm = FreeCAD.Units.Quantity(mat_obj.Material["RelativePermittivity"]).Value
+            vacuum_perm = FreeCAD.Units.Quantity("8.85419e-12 F/m").getValueAs("C/(mV*mm)").Value
+            abs_perm = vacuum_perm * rel_perm
         # write material properties
         f.write(f"** FreeCAD material name: {mat_info_name}\n")
         f.write(f"** {mat_label}\n")
@@ -119,7 +119,7 @@ def is_density_needed():
                 f.write(f"{SH_in_JkgK:.13G},{DV_in_tmms:.13G}\n")
         if ccxwriter.analysis_type == "electrostatic":
             f.write("*CONDUCTIVITY\n")
-            f.write(f"{TC_in_WmK:.13G}\n")
+            f.write(f"{abs_perm:.13G}\n")
 
         # nonlinear material properties
         if ccxwriter.solver_obj.MaterialNonlinearity == "nonlinear":

src/Mod/Fem/femsolver/calculix/write_step_equation.py

@@ -79,7 +79,7 @@ def write_step_equation(f, ccxwriter):
     elif ccxwriter.analysis_type == "buckling":
         analysis_type = "*BUCKLE"
     elif ccxwriter.analysis_type == "electrostatic":
-        analysis_type = "*HEAT TRANSFER, STEADY STATE, MATRIXSTORAGE"
+        analysis_type = "*HEAT TRANSFER, STEADY STATE"
     # analysis line --> solver type
     # https://forum.freecad.org/viewtopic.php?f=18&t=43178
     if ccxwriter.solver_obj.MatrixSolverType == "default":