diff --git a/src/Mod/Fem/CMakeLists.txt b/src/Mod/Fem/CMakeLists.txt index e90422449e..a0390b1035 100755 --- a/src/Mod/Fem/CMakeLists.txt +++ b/src/Mod/Fem/CMakeLists.txt @@ -248,6 +248,8 @@ SET(FemSolverCalculix_SRCS femsolver/calculix/write_constraint_centrif.py femsolver/calculix/write_constraint_contact.py femsolver/calculix/write_constraint_displacement.py + femsolver/calculix/write_constraint_electricchargedensity.py + femsolver/calculix/write_constraint_electrostatic.py femsolver/calculix/write_constraint_fixed.py femsolver/calculix/write_constraint_fluidsection.py femsolver/calculix/write_constraint_force.py diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_electricchargedensity.py b/src/Mod/Fem/femsolver/calculix/write_constraint_electricchargedensity.py new file mode 100644 index 0000000000..28903d558e --- /dev/null +++ b/src/Mod/Fem/femsolver/calculix/write_constraint_electricchargedensity.py @@ -0,0 +1,149 @@ +# SPDX-License-Identifier: LGPL-2.1-or-later + +# *************************************************************************** +# * Copyright (c) 2025 Mario Passaglia * +# * * +# * This file is part of FreeCAD. * +# * * +# * FreeCAD is free software: you can redistribute it and/or modify it * +# * under the terms of the GNU Lesser General Public License as * +# * published by the Free Software Foundation, either version 2.1 of the * +# * License, or (at your option) any later version. * +# * * +# * FreeCAD is distributed in the hope that it will be useful, but * +# * WITHOUT ANY WARRANTY; without even the implied warranty of * +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * +# * Lesser General Public License for more details. * +# * * +# * You should have received a copy of the GNU Lesser General Public * +# * License along with FreeCAD. If not, see * +# * . * +# * * +# *************************************************************************** + +__title__ = "FreeCAD FEM calculix constraint electric charge density" +__author__ = "Mario Passaglia" +__url__ = "https://www.freecad.org" + + +def get_analysis_types(): + return ["electromagnetic"] + + +def get_sets_name(): + return "constraints_electricchargedensity_node_sets" + + +def get_constraint_title(): + return "Electric charge density constraint applied" + + +def write_meshdata_constraint(f, femobj, den_obj, ccxwriter): + + # print("MESHDATA: ", femobj) + if ccxwriter.solver_obj.ElectromagneticMode != "electrostatic": + return + + if den_obj.Mode in ["Source", "Total Source"]: + msg, data = femobj["ChargeDensityElements"] + f.write(f"** {msg}\n") + f.write(f"*ELSET,ELSET={den_obj.Name}\n") + for ref, elem in data: + for e in elem: + f.write(f"{e},\n") + + +# if isinstance(femobj["FEMElements"], str): +# f.write("{}\n".format(femobj["FEMElements"])) +# else: +# for e in femobj["FEMElements"]: +# f.write(f"{e},\n") + + +# if femobj["Object"].BoundaryCondition == "Dirichlet": +# f.write(f"*NSET,NSET={den_obj.Name}\n") +# for n in femobj["Nodes"]: +# f.write(f"{n},\n") + + +def get_before_write_meshdata_constraint(): + return "" + + +def get_after_write_meshdata_constraint(): + return "" + + +def get_before_write_constraint(): + return "" + + +def get_after_write_constraint(): + return "" + + +def write_constraint(f, femobj, den_obj, ccxwriter): + + # floats read from ccx should use {:.13G}, see comment in writer module + # if den_obj.BoundaryCondition == "Dirichlet": + # f.write("*BOUNDARY\n") + # f.write("{},11,11,{:.13G}\n".format(den_obj.Name, den_obj.Potential.getValueAs("mV").Value)) + # f.write("\n") + if ccxwriter.solver_obj.ElectromagneticMode != "electrostatic": + return + + match den_obj.Mode: + case "Source": + density = den_obj.SourceChargeDensity + case "Total Source": + density = den_obj.Proxy.get_total_source_density(den_obj) + case "Interface": + density = den_obj.InterfaceChargeDensity + case "Total Interface": + density = den_obj.Proxy.get_total_interface_density(den_obj) + + if den_obj.Mode in ["Source", "Total Source"]: + f.write("*DFLUX\n") + f.write("{},BF,{:.13G}\n".format(den_obj.Name, density.getValueAs("C/mm^3").Value)) + f.write("\n") + + elif den_obj.Mode in ["Interface", "Total Interface"]: + # check internal interface + density = density.getValueAs("C/mm^2").Value + internal = _check_shared_interface(den_obj) + for feat, refs in femobj["ChargeDensityFaces"]: + f.write("** " + feat + "\n") + f.write("*DFLUX\n") + for ref in refs: + d = density + if ref[0] in internal: + d = density / 2 + for face, fno in ref[1]: + f.write("{},S{},{:.13G}\n".format(face, fno, d)) + + f.write("\n") + + +def _check_shared_interface(den_obj): + """ + Check if reference is internal shared subshape + For example, shared face in compsolid + """ + internal = [] + for o, sub in den_obj.References: + for elem in sub: + found = [] + elem_i = o.getSubObject(elem) + if elem_i.ShapeType == "Face": + for s in o.Shape.Solids: + found.append(any([q.isSame(elem_i) for q in s.Faces])) + if sum(found) > 1: + internal.append((o, (elem,))) + + if elem_i.ShapeType == "Edge": + for s in o.Shape.Faces: + found.append(any([q.isSame(elem_i) for q in s.Edges])) + if sum(found) > 1: + internal.append((o, (elem,))) + + return internal diff --git a/src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py b/src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py new file mode 100644 index 0000000000..1b5b9ef65f --- /dev/null +++ b/src/Mod/Fem/femsolver/calculix/write_constraint_electrostatic.py @@ -0,0 +1,118 @@ +# SPDX-License-Identifier: LGPL-2.1-or-later + +# *************************************************************************** +# * Copyright (c) 2025 Mario Passaglia * +# * * +# * This file is part of FreeCAD. * +# * * +# * FreeCAD is free software: you can redistribute it and/or modify it * +# * under the terms of the GNU Lesser General Public License as * +# * published by the Free Software Foundation, either version 2.1 of the * +# * License, or (at your option) any later version. * +# * * +# * FreeCAD is distributed in the hope that it will be useful, but * +# * WITHOUT ANY WARRANTY; without even the implied warranty of * +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * +# * Lesser General Public License for more details. * +# * * +# * You should have received a copy of the GNU Lesser General Public * +# * License along with FreeCAD. If not, see * +# * . * +# * * +# *************************************************************************** + +__title__ = "FreeCAD FEM calculix constraint electrostatic" +__author__ = "Mario Passaglia" +__url__ = "https://www.freecad.org" + +import FreeCAD + + +def get_analysis_types(): + return ["electromagnetic"] + + +def get_sets_name(): + return "constraints_electrostaticpotential_node_sets" + + +def get_constraint_title(): + return "Fixed electrostatic constraint applied" + + +def write_meshdata_constraint(f, femobj, pot_obj, ccxwriter): + + if ccxwriter.solver_obj.ElectromagneticMode != "electrostatic": + return + + 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(): + return "" + + +def get_after_write_meshdata_constraint(): + return "" + + +def get_before_write_constraint(): + return "" + + +def get_after_write_constraint(): + return "" + + +def write_constraint(f, femobj, pot_obj, ccxwriter): + + if ccxwriter.solver_obj.ElectromagneticMode != "electrostatic": + return + + # floats read from ccx should use {:.13G}, see comment in writer module + 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)) + elif pot_obj.BoundaryCondition == "Neumann": + density = pot_obj.ElectricFluxDensity.getValueAs("C/mm^2").Value + # check internal interface + internal = _check_shared_interface(pot_obj) + for feat, refs in femobj["ElectricFluxFaces"]: + f.write("** " + feat + "\n") + f.write("*DFLUX\n") + for ref in refs: + d = density + if ref[0] in internal: + d = density / 2 + for face, fno in ref[1]: + f.write("{},S{},{:.13G}\n".format(face, fno, d)) + + f.write("\n") + + +def _check_shared_interface(pot_obj): + """ + Check if reference is internal shared subshape + For example, shared face in compsolid + """ + internal = [] + for o, sub in pot_obj.References: + for elem in sub: + found = [] + elem_i = o.getSubObject(elem) + if elem_i.ShapeType == "Face": + for s in o.Shape.Solids: + found.append(any([q.isSame(elem_i) for q in s.Faces])) + if sum(found) > 1: + internal.append((o, (elem,))) + + if elem_i.ShapeType == "Edge": + for s in o.Shape.Faces: + found.append(any([q.isSame(elem_i) for q in s.Edges])) + if sum(found) > 1: + internal.append((o, (elem,))) + + return internal diff --git a/src/Mod/Fem/femsolver/calculix/write_femelement_material.py b/src/Mod/Fem/femsolver/calculix/write_femelement_material.py index 1b08ce23ec..33ce74188f 100644 --- a/src/Mod/Fem/femsolver/calculix/write_femelement_material.py +++ b/src/Mod/Fem/femsolver/calculix/write_femelement_material.py @@ -27,6 +27,7 @@ __url__ = "https://www.freecad.org" import FreeCAD +from femtools import constants def write_femelement_material(f, ccxwriter): @@ -92,7 +93,13 @@ def write_femelement_material(f, ccxwriter): KV = FreeCAD.Units.Quantity(mat_obj.Material["KinematicViscosity"]) KV_in_mm2s = KV.getValueAs("mm^2/s").Value DV_in_tmms = KV_in_mm2s * density_in_tonne_per_mm3 - + if ( + ccxwriter.analysis_type == "electromagnetic" + and ccxwriter.solver_obj.ElectromagneticMode == "electrostatic" + ): + rel_perm = FreeCAD.Units.Quantity(mat_obj.Material["RelativePermittivity"]).Value + vacuum_perm = FreeCAD.Units.Quantity(constants.vacuum_permittivity()) + abs_perm = vacuum_perm.getValueAs("C/(mV*mm)").Value * rel_perm # write material properties f.write(f"** FreeCAD material name: {mat_info_name}\n") f.write(f"** {mat_label}\n") @@ -114,6 +121,12 @@ def write_femelement_material(f, ccxwriter): elif mat_obj.Category == "Fluid": f.write("*FLUID CONSTANTS\n") f.write(f"{SH_in_JkgK:.13G},{DV_in_tmms:.13G}\n") + if ( + ccxwriter.analysis_type == "electromagnetic" + and ccxwriter.solver_obj.ElectromagneticMode == "electrostatic" + ): + f.write("*CONDUCTIVITY\n") + f.write(f"{abs_perm:.13G}\n") # nonlinear material properties if ccxwriter.solver_obj.MaterialNonlinearity == "nonlinear": diff --git a/src/Mod/Fem/femsolver/calculix/write_step_equation.py b/src/Mod/Fem/femsolver/calculix/write_step_equation.py index 5c57dfb27a..83a97692b5 100644 --- a/src/Mod/Fem/femsolver/calculix/write_step_equation.py +++ b/src/Mod/Fem/femsolver/calculix/write_step_equation.py @@ -78,6 +78,8 @@ def write_step_equation(f, ccxwriter): analysis_type = "*NO ANALYSIS" elif ccxwriter.analysis_type == "buckling": analysis_type = "*BUCKLE" + elif ccxwriter.analysis_type == "electromagnetic": + 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": @@ -95,9 +97,7 @@ def write_step_equation(f, ccxwriter): # analysis line --> user defined incrementations --> parameter DIRECT # --> completely switch off ccx automatic incrementation if ccxwriter.solver_obj.IterationsUserDefinedIncrementations: - if ccxwriter.analysis_type == "static": - analysis_type += ", DIRECT" - elif ccxwriter.analysis_type == "thermomech": + if ccxwriter.analysis_type in ["static", "thermomech", "electromagnetic"]: analysis_type += ", DIRECT" elif ccxwriter.analysis_type == "frequency": FreeCAD.Console.PrintMessage( diff --git a/src/Mod/Fem/femsolver/calculix/write_step_output.py b/src/Mod/Fem/femsolver/calculix/write_step_output.py index 6b62441485..9bc00687fa 100644 --- a/src/Mod/Fem/femsolver/calculix/write_step_output.py +++ b/src/Mod/Fem/femsolver/calculix/write_step_output.py @@ -47,6 +47,8 @@ def write_step_output(f, ccxwriter): f.write("U, NT\n") else: f.write("MF, PS\n") + elif ccxwriter.analysis_type == "electromagnetic": + f.write("NT\n") else: f.write("U\n") if not ccxwriter.member.geos_fluidsection: @@ -56,6 +58,8 @@ def write_step_output(f, ccxwriter): variables += ", HFL" if ccxwriter.solver_obj.MaterialNonlinearity == "nonlinear": variables += ", PEEQ" + if ccxwriter.analysis_type == "electromagnetic": + variables = "HFL" f.write(variables + "\n") diff --git a/src/Mod/Fem/femsolver/calculix/writer.py b/src/Mod/Fem/femsolver/calculix/writer.py index db2a38c9ce..8f830bf159 100644 --- a/src/Mod/Fem/femsolver/calculix/writer.py +++ b/src/Mod/Fem/femsolver/calculix/writer.py @@ -53,6 +53,8 @@ from . import write_constraint_selfweight as con_selfweight from . import write_constraint_temperature as con_temperature from . import write_constraint_tie as con_tie from . import write_constraint_transform as con_transform +from . import write_constraint_electricchargedensity as con_electricchargedensity +from . import write_constraint_electrostatic as con_electrostatic from . import write_femelement_geometry from . import write_femelement_material from . import write_femelement_matgeosets @@ -157,6 +159,10 @@ class FemInputWriterCcx(writerbase.FemInputWriter): self.write_constraints_meshsets(inpfile, self.member.cons_planerotation, con_planerotation) self.write_constraints_meshsets(inpfile, self.member.cons_transform, con_transform) self.write_constraints_meshsets(inpfile, self.member.cons_temperature, con_temperature) + self.write_constraints_meshsets( + inpfile, self.member.cons_electricchargedensity, con_electricchargedensity + ) + self.write_constraints_meshsets(inpfile, self.member.cons_electrostatic, con_electrostatic) # surface sets self.write_constraints_meshsets(inpfile, self.member.cons_contact, con_contact) @@ -194,6 +200,10 @@ class FemInputWriterCcx(writerbase.FemInputWriter): self.write_constraints_meshsets(inpfile, self.member.cons_pressure, con_pressure) self.write_constraints_propdata(inpfile, self.member.cons_temperature, con_temperature) self.write_constraints_meshsets(inpfile, self.member.cons_heatflux, con_heatflux) + self.write_constraints_propdata( + inpfile, self.member.cons_electricchargedensity, con_electricchargedensity + ) + self.write_constraints_propdata(inpfile, self.member.cons_electrostatic, con_electrostatic) con_fluidsection.write_constraints_fluidsection(inpfile, self) # output and step end diff --git a/src/Mod/Fem/femsolver/writerbase.py b/src/Mod/Fem/femsolver/writerbase.py index f480b1f78c..aea5614294 100644 --- a/src/Mod/Fem/femsolver/writerbase.py +++ b/src/Mod/Fem/femsolver/writerbase.py @@ -128,6 +128,8 @@ class FemInputWriter: self.temperature_objects = member.cons_temperature self.tie_objects = member.cons_tie self.transform_objects = member.cons_transform + self.electrostatic_objects = member.cons_electrostatic + self.electricchargedensity_objects = member.cons_electricchargedensity # meshdatagetter, for compatibility, same with all getter methods self.meshdatagetter = meshsetsgetter.MeshSetsGetter( diff --git a/src/Mod/Fem/femtools/membertools.py b/src/Mod/Fem/femtools/membertools.py index 48d97c162a..1163861ca3 100644 --- a/src/Mod/Fem/femtools/membertools.py +++ b/src/Mod/Fem/femtools/membertools.py @@ -282,6 +282,10 @@ class AnalysisMember: self.cons_temperature = self.get_several_member("Fem::ConstraintTemperature") self.cons_tie = self.get_several_member("Fem::ConstraintTie") self.cons_transform = self.get_several_member("Fem::ConstraintTransform") + self.cons_electrostatic = self.get_several_member("Fem::ConstraintElectrostaticPotential") + self.cons_electricchargedensity = self.get_several_member( + "Fem::ConstraintElectricChargeDensity" + ) def get_several_member(self, t): return get_several_member(self.analysis, t)