Merge pull request #21100 from marioalexis84/fem-elmer_quantity_format

Fem: Use string format to round Elmer quantity values - fixes #20938
This commit is contained in:
Chris Hennes
2025-05-04 16:57:31 -05:00
committed by GitHub
5 changed files with 165 additions and 144 deletions

View File

@@ -29,6 +29,9 @@ __url__ = "https://www.freecad.org"
# \ingroup FEM
# \brief constraint electric charge density object
from FreeCAD import Units
from . import base_fempythonobject
_PropHelper = base_fempythonobject._PropHelper
@@ -85,3 +88,43 @@ class ConstraintElectricChargeDensity(base_fempythonobject.BaseFemPythonObject):
)
return prop
def get_total_source_density(self, obj):
"""
Calculate density for `Total Source` mode.
"""
size = 0
items = []
for feat, sub_elem in obj.References:
for name in sub_elem:
sub = feat.getSubObject(name)
if sub.ShapeType == "Solid":
size += sub.Volume
items.append(name)
elif sub.ShapeType == "Face":
size += sub.Area
items.append(name)
if items:
vol = Units.Quantity(f"{size} mm^3")
return obj.TotalCharge / vol
def get_total_interface_density(self, obj):
"""
Calculate density for `Total Interface` mode.
"""
size = 0
items = []
for feat, sub_elem in obj.References:
for name in sub_elem:
sub = feat.getSubObject(name)
if sub.ShapeType == "Face":
size += sub.Area
items.append(name)
elif sub.ShapeType == "Edge":
size += sub.Length
items.append(name)
if items:
area = Units.Quantity(f"{size} mm^2")
return obj.TotalCharge / area

View File

@@ -106,10 +106,9 @@ class ESwriter:
equation.PotentialDifference = 0.0
def handleElectrostaticConstants(self):
permittivity = self.write.convert(
self.write.constsdef["PermittivityOfVacuum"], "T^4*I^2/(L^3*M)"
permittivity = Units.Quantity(self.write.constsdef["PermittivityOfVacuum"]).getValueAs(
"F/m"
)
permittivity = round(permittivity, 20) # to get rid of numerical artifacts
self.write.constant("Permittivity Of Vacuum", permittivity)
def handleElectrostaticMaterial(self, bodies):
@@ -132,15 +131,11 @@ class ESwriter:
self.write.boundary(name, "! FreeCAD Name", obj.Label)
if obj.BoundaryCondition == "Dirichlet":
if obj.PotentialEnabled:
self.write.boundary(
name, "Potential", obj.Potential.getValueAs("V").Value
)
potential = obj.Potential.getValueAs("V")
self.write.boundary(name, "Potential", potential)
elif obj.BoundaryCondition == "Neumann":
self.write.boundary(
name,
"Electric Flux",
obj.ElectricFluxDensity.getValueAs("C/m^2").Value,
)
flux_density = obj.ElectricFluxDensity.getValueAs("C/m^2")
self.write.boundary(name, "Electric Flux", flux_density)
if obj.PotentialConstant:
self.write.boundary(name, "Potential Constant", True)
if obj.ElectricInfinity:
@@ -152,58 +147,38 @@ class ESwriter:
self.write.handled(obj)
for obj in self.write.getMember("Fem::ConstraintElectricChargeDensity"):
if obj.Mode not in ["Interface", "Total Interface"]:
continue
match obj.Mode:
case "Interface":
density = obj.InterfaceChargeDensity
case "Total Interface":
density = obj.Proxy.get_total_interface_density(obj)
case _:
continue
size = 0
items = []
for feat, sub_elem in obj.References:
for name in sub_elem:
sub = feat.getSubObject(name)
if sub.ShapeType == "Face":
size += sub.Area
items.append(name)
elif sub.ShapeType == "Edge":
size += sub.Length
items.append(name)
if items:
if obj.Mode == "Interface":
density = obj.InterfaceChargeDensity.getValueAs("C/m^2").Value
elif obj.Mode == "Total Interface":
area = Units.Quantity(f"{size} mm^2")
density = (obj.TotalCharge / area).getValueAs("C/m^2").Value
for name in items:
self.write.boundary(name, "! FreeCAD Name", obj.Label)
self.write.boundary(name, "Surface Charge Density", round(density, 6))
self.write.boundary(
name,
"Surface Charge Density",
density.getValueAs("C/m^2"),
)
self.write.handled(obj)
def handleElectrostaticBodyForces(self):
for obj in self.write.getMember("Fem::ConstraintElectricChargeDensity"):
if obj.Mode not in ["Source", "Total Source"]:
continue
match obj.Mode:
case "Source":
density = obj.SourceChargeDensity
case "Total Source":
density = obj.Proxy.get_total_source_density(obj)
case _:
continue
size = 0
items = []
for feat, sub_elem in obj.References:
for name in sub_elem:
sub = feat.getSubObject(name)
if sub.ShapeType == "Solid":
size += sub.Volume
items.append(name)
elif sub.ShapeType == "Face":
size += sub.Area
items.append(name)
if items:
if obj.Mode == "Source":
density = obj.SourceChargeDensity.getValueAs("C/m^3").Value
elif obj.Mode == "Total Source":
vol = Units.Quantity(f"{size} mm^3")
density = (obj.TotalCharge / vol).getValueAs("C/m^3").Value
for name in items:
self.write.bodyForce(name, "! FreeCAD Name", obj.Label)
self.write.bodyForce(name, "Charge Density", round(density, 6))
self.write.bodyForce(name, "Charge Density", density.getValueAs("C/m^3"))
self.write.handled(obj)

View File

@@ -65,7 +65,7 @@ class MgDyn2Dwriter:
s["Exec Solver"] = "Always"
s["Procedure"] = sifio.FileAttr("MagnetoDynamics/MagnetoDynamicsCalcFields")
if equation.IsHarmonic:
s["Angular Frequency"] = float(Units.Quantity(equation.AngularFrequency).Value)
s["Angular Frequency"] = equation.AngularFrequency.getValueAs("Hz")
s["Potential Variable"] = "Potential"
if equation.CalculateCurrentDensity is True:
s["Calculate Current Density"] = True
@@ -92,16 +92,15 @@ class MgDyn2Dwriter:
return s
def handleMagnetodynamic2DConstants(self):
permeability = self.write.convert(
self.write.constsdef["PermeabilityOfVacuum"], "M*L/(T^2*I^2)"
permeability = Units.Quantity(self.write.constsdef["PermeabilityOfVacuum"]).getValueAs(
"H/m"
)
# we round in the following to get rid of numerical artifacts
self.write.constant("Permeability Of Vacuum", round(permeability, 20))
self.write.constant("Permeability Of Vacuum", permeability)
permittivity = self.write.convert(
self.write.constsdef["PermittivityOfVacuum"], "T^4*I^2/(L^3*M)"
permittivity = Units.Quantity(self.write.constsdef["PermittivityOfVacuum"]).getValueAs(
"F/m"
)
self.write.constant("Permittivity Of Vacuum", round(permittivity, 20))
self.write.constant("Permittivity Of Vacuum", permittivity)
def handleMagnetodynamic2DMaterial(self, bodies):
# check that all bodies have a set material
@@ -125,9 +124,9 @@ class MgDyn2Dwriter:
"The relative permeability must be specified for all materials.\n\n"
)
self.write.material(name, "Name", m["Name"])
conductivity = self.write.convert(m["ElectricalConductivity"], "T^3*I^2/(L^3*M)")
conductivity = round(conductivity, 10) # to get rid of numerical artifacts
conductivity = Units.Quantity(m["ElectricalConductivity"]).getValueAs("S/m")
self.write.material(name, "Electric Conductivity", conductivity)
self.write.material(name, "Relative Permeability", float(m["RelativePermeability"]))
# permittivity might be necessary for the post processor
if "RelativePermittivity" in m:
@@ -137,29 +136,29 @@ class MgDyn2Dwriter:
def _outputMagnetodynamic2DBodyForce(self, obj, name, equation):
if femutils.is_derived_from(obj, "Fem::ConstraintCurrentDensity") and obj.Mode == "Normal":
currentDensity = obj.NormalCurrentDensity_re.getValueAs("A/m^2").Value
self.write.bodyForce(name, "Current Density", round(currentDensity, 6))
current_density = obj.NormalCurrentDensity_re.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density", current_density)
# imaginaries are only needed for harmonic equation
if equation.IsHarmonic:
currentDensity = obj.NormalCurrentDensity_im.getValueAs("A/m^2").Value
self.write.bodyForce(name, "Current Density Im", round(currentDensity, 6))
current_density = obj.NormalCurrentDensity_im.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density Im", current_density)
if femutils.is_derived_from(obj, "Fem::ConstraintMagnetization"):
# output only if magnetization is enabled and needed
if obj.EnableMagnetization_1:
magnetization = float(obj.Magnetization_re_1.getValueAs("A/m"))
self.write.material(name, "Magnetization 1", round(magnetization, 6))
magnetization = obj.Magnetization_re_1.getValueAs("A/m")
self.write.material(name, "Magnetization 1", magnetization)
if obj.EnableMagnetization_2:
magnetization = float(obj.Magnetization_re_2.getValueAs("A/m"))
self.write.material(name, "Magnetization 2", round(magnetization, 6))
magnetization = obj.Magnetization_re_2.getValueAs("A/m")
self.write.material(name, "Magnetization 2", magnetization)
# imaginaries are only needed for harmonic equation
if equation.IsHarmonic:
if obj.EnableMagnetization_1:
magnetization = float(obj.Magnetization_im_1.getValueAs("A/m"))
self.write.material(name, "Magnetization Im 1", round(magnetization, 6))
magnetization = obj.Magnetization_im_1.getValueAs("A/m")
self.write.material(name, "Magnetization Im 1", magnetization)
if obj.EnableMagnetization_2:
magnetization = float(obj.Magnetization_im_2.getValueAs("A/m"))
self.write.material(name, "Magnetization Im 2", round(magnetization, 6))
magnetization = obj.Magnetization_im_2.getValueAs("A/m")
self.write.material(name, "Magnetization Im 2", magnetization)
def handleMagnetodynamic2DBodyForces(self, bodies, equation):
currentDensities = self.write.getMember("Fem::ConstraintCurrentDensity")
@@ -209,8 +208,8 @@ class MgDyn2Dwriter:
self.write.boundary(name, "! FreeCAD Name", obj.Label)
if obj.PotentialEnabled:
if hasattr(obj, "Potential"):
potential = float(obj.Potential.getValueAs("V"))
self.write.boundary(name, "Potential", round(potential, 6))
potential = obj.Potential.getValueAs("V")
self.write.boundary(name, "Potential", potential)
if obj.ElectricInfinity:
self.write.boundary(name, "Infinity BC", True)
self.write.handled(obj)
@@ -221,8 +220,8 @@ class MgDyn2Dwriter:
raise general_writer.WriteError("The angular frequency must not be zero.\n\n")
self.write.equation(b, "Name", equation.Name)
if equation.IsHarmonic:
frequency = float(Units.Quantity(equation.AngularFrequency).Value)
self.write.equation(b, "Angular Frequency", round(frequency, 6))
frequency = equation.AngularFrequency.getValueAs("Hz")
self.write.equation(b, "Angular Frequency", frequency)
## @}

View File

@@ -53,9 +53,8 @@ class MgDynwriter:
s["Equation"] = "MgDynHarmonic"
s["Procedure"] = sifio.FileAttr("MagnetoDynamics/WhitneyAVHarmonicSolver")
s["Variable"] = "av[av re:1 av im:1]"
# round to get rid of numerical artifacts
frequency = float(Units.Quantity(equation.AngularFrequency).Value)
s["Angular Frequency"] = round(frequency, 6)
frequency = equation.AngularFrequency.getValueAs("Hz")
s["Angular Frequency"] = frequency
s["Exec Solver"] = "Always"
s["Optimize Bandwidth"] = True
s["Stabilize"] = equation.Stabilize
@@ -88,8 +87,8 @@ class MgDynwriter:
s["Exec Solver"] = "Before Saving"
s["Procedure"] = sifio.FileAttr("MagnetoDynamics/MagnetoDynamicsCalcFields")
if equation.IsHarmonic:
frequency = float(Units.Quantity(equation.AngularFrequency).Value)
s["Angular Frequency"] = round(frequency, 6)
frequency = equation.AngularFrequency.getValueAs("Hz")
s["Angular Frequency"] = frequency
s["Potential Variable"] = "av"
if equation.CalculateCurrentDensity is True:
s["Calculate Current Density"] = True
@@ -118,16 +117,15 @@ class MgDynwriter:
return s
def handleMagnetodynamicConstants(self):
permeability = self.write.convert(
self.write.constsdef["PermeabilityOfVacuum"], "M*L/(T^2*I^2)"
permeability = Units.Quantity(self.write.constsdef["PermeabilityOfVacuum"]).getValueAs(
"H/m"
)
# we round in the following to get rid of numerical artifacts
self.write.constant("Permeability Of Vacuum", round(permeability, 20))
self.write.constant("Permeability Of Vacuum", permeability)
permittivity = self.write.convert(
self.write.constsdef["PermittivityOfVacuum"], "T^4*I^2/(L^3*M)"
permittivity = Units.Quantity(self.write.constsdef["PermittivityOfVacuum"]).getValueAs(
"F/m"
)
self.write.constant("Permittivity Of Vacuum", round(permittivity, 20))
self.write.constant("Permittivity Of Vacuum", permittivity)
def handleMagnetodynamicMaterial(self, bodies):
# check that all bodies have a set material
@@ -151,8 +149,7 @@ class MgDynwriter:
"The relative permeability must be specified for all materials.\n\n"
)
self.write.material(name, "Name", m["Name"])
conductivity = self.write.convert(m["ElectricalConductivity"], "T^3*I^2/(L^3*M)")
conductivity = round(conductivity, 10) # to get rid of numerical artifacts
conductivity = Units.Quantity(m["ElectricalConductivity"]).getValueAs("S/m")
self.write.material(name, "Electric Conductivity", conductivity)
self.write.material(name, "Relative Permeability", float(m["RelativePermeability"]))
# permittivity might be necessary for the post processor
@@ -165,58 +162,58 @@ class MgDynwriter:
if femutils.is_derived_from(obj, "Fem::ConstraintCurrentDensity") and obj.Mode == "Custom":
# output only if current density is enabled and needed
if obj.EnableCurrentDensity_re_1:
currentDensity = float(obj.CurrentDensity_re_1.getValueAs("A/m^2"))
self.write.bodyForce(name, "Current Density 1", round(currentDensity, 6))
current_density = obj.CurrentDensity_re_1.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density 1", current_density)
if obj.EnableCurrentDensity_re_2:
currentDensity = float(obj.CurrentDensity_re_2.getValueAs("A/m^2"))
self.write.bodyForce(name, "Current Density 2", round(currentDensity, 6))
current_density = obj.CurrentDensity_re_2.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density 2", current_density)
if obj.EnableCurrentDensity_re_3:
currentDensity = float(obj.CurrentDensity_re_3.getValueAs("A/m^2"))
self.write.bodyForce(name, "Current Density 3", round(currentDensity, 6))
current_density = obj.CurrentDensity_re_3.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density 3", current_density)
# imaginaries are only needed for harmonic equation
if equation.IsHarmonic:
if obj.EnableCurrentDensity_im_1:
currentDensity = float(obj.CurrentDensity_im_1.getValueAs("A/m^2"))
self.write.bodyForce(name, "Current Density Im 1", round(currentDensity, 6))
current_density = obj.CurrentDensity_im_1.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density Im 1", current_density)
if obj.EnableCurrentDensity_im_2:
currentDensity = float(obj.CurrentDensity_im_2.getValueAs("A/m^2"))
self.write.bodyForce(name, "Current Density Im 2", round(currentDensity, 6))
current_density = obj.CurrentDensity_im_2.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density Im 2", current_density)
if obj.EnableCurrentDensity_im_3:
currentDensity = float(obj.CurrentDensity_im_3.getValueAs("A/m^2"))
self.write.bodyForce(name, "Current Density Im 3", round(currentDensity, 6))
current_density = obj.CurrentDensity_im_3.getValueAs("A/m^2")
self.write.bodyForce(name, "Current Density Im 3", current_density)
if femutils.is_derived_from(obj, "Fem::ConstraintMagnetization"):
# output only if magnetization is enabled and needed
if obj.EnableMagnetization_1:
magnetization = float(obj.Magnetization_re_1.getValueAs("A/m"))
magnetization = obj.Magnetization_re_1.getValueAs("A/m")
self.write.bodyForce(name, "Magnetization 1", magnetization)
if obj.EnableMagnetization_2:
magnetization = float(obj.Magnetization_re_2.getValueAs("A/m"))
magnetization = obj.Magnetization_re_2.getValueAs("A/m")
self.write.bodyForce(name, "Magnetization 2", magnetization)
if obj.EnableMagnetization_3:
magnetization = float(obj.Magnetization_re_3.getValueAs("A/m"))
magnetization = obj.Magnetization_re_3.getValueAs("A/m")
self.write.bodyForce(name, "Magnetization 3", magnetization)
# imaginaries are only needed for harmonic equation
if equation.IsHarmonic:
if obj.EnableMagnetization_1:
magnetization = float(obj.Magnetization_im_1.getValueAs("A/m"))
magnetization = obj.Magnetization_im_1.getValueAs("A/m")
self.write.bodyForce(name, "Magnetization Im 1", magnetization)
if obj.EnableMagnetization_2:
magnetization = float(obj.Magnetization_im_2.getValueAs("A/m"))
magnetization = obj.Magnetization_im_2.getValueAs("A/m")
self.write.bodyForce(name, "Magnetization Im 2", magnetization)
if obj.EnableMagnetization_3:
magnetization = float(obj.Magnetization_im_3.getValueAs("A/m"))
magnetization = obj.Magnetization_im_3.getValueAs("A/m")
self.write.bodyForce(name, "Magnetization Im 3", magnetization)
if femutils.is_derived_from(obj, "Fem::ConstraintElectrostaticPotential"):
if obj.PotentialEnabled:
# output only if potential is enabled and needed
potential = float(obj.Potential.getValueAs("V"))
self.write.bodyForce(name, "Electric Potential", round(potential, 6))
potential = obj.Potential.getValueAs("V")
self.write.bodyForce(name, "Electric Potential", potential)
# imaginary is only needed for harmonic equation
if equation.IsHarmonic:
potential = float(obj.AV_im.getValueAs("V"))
self.write.bodyForce(name, "Electric Potential Im", round(potential, 6))
potential = obj.AV_im.getValueAs("V")
self.write.bodyForce(name, "Electric Potential Im", potential)
def handleMagnetodynamicBodyForces(self, bodies, equation):
# the current density can either be a body force or a boundary constraint
@@ -277,46 +274,46 @@ class MgDynwriter:
def _outputMagnetodynamicBndConditions(self, obj, name, equation):
if femutils.is_derived_from(obj, "Fem::ConstraintCurrentDensity") and obj.Mode == "Normal":
currentDensity = float(obj.NormalCurrentDensity_re.getValueAs("A/m^2"))
self.write.boundary(name, "Electric Current Density", round(currentDensity, 6))
current_density = obj.NormalCurrentDensity_re.getValueAs("A/m^2")
self.write.boundary(name, "Electric Current Density", current_density)
# imaginaries are only needed for harmonic equation
if equation.IsHarmonic:
currentDensity = float(obj.NormalCurrentDensity_im.getValueAs("A/m^2"))
self.write.boundary(name, "Electric Current Density Im", round(currentDensity, 6))
current_density = obj.NormalCurrentDensity_im.getValueAs("A/m^2")
self.write.boundary(name, "Electric Current Density Im", current_density)
if femutils.is_derived_from(obj, "Fem::ConstraintElectrostaticPotential"):
if obj.EnableAV:
potential = obj.AV_re.getValueAs("V").Value
potential = obj.AV_re.getValueAs("V")
if equation.IsHarmonic:
self.write.boundary(name, "AV re", round(potential, 6))
potential = obj.AV_im.getValueAs("V").Value
self.write.boundary(name, "AV im", round(potential, 6))
self.write.boundary(name, "AV re", potential)
potential = obj.AV_im.getValueAs("V")
self.write.boundary(name, "AV im", potential)
else:
self.write.boundary(name, "AV", round(potential, 6))
self.write.boundary(name, "AV", potential)
if obj.EnableAV_1:
potential = obj.AV_re_1.getValueAs("Wb/m").Value
potential = obj.AV_re_1.getValueAs("Wb/m")
if equation.IsHarmonic:
self.write.boundary(name, "AV re {e} 1", round(potential, 6))
potential = obj.AV_im_1.getValueAs("Wb/m").Value
self.write.boundary(name, "AV im {e} 1", round(potential, 6))
self.write.boundary(name, "AV re {e} 1", potential)
potential = obj.AV_im_1.getValueAs("Wb/m")
self.write.boundary(name, "AV im {e} 1", potential)
else:
self.write.boundary(name, "AV {e} 1", round(potential, 6))
self.write.boundary(name, "AV {e} 1", potential)
if obj.EnableAV_2:
potential = obj.AV_re_2.getValueAs("Wb/m").Value
potential = obj.AV_re_2.getValueAs("Wb/m")
if equation.IsHarmonic:
self.write.boundary(name, "AV re {e} 2", round(potential, 6))
potential = obj.AV_im_2.getValueAs("Wb/m").Value
self.write.boundary(name, "AV im {e} 2", round(potential, 6))
self.write.boundary(name, "AV re {e} 2", potential)
potential = obj.AV_im_2.getValueAs("Wb/m")
self.write.boundary(name, "AV im {e} 2", potential)
else:
self.write.boundary(name, "AV {e} 2", round(potential, 6))
self.write.boundary(name, "AV {e} 2", potential)
if obj.EnableAV_3:
potential = obj.AV_re_3.getValueAs("Wb/m").Value
potential = obj.AV_re_3.getValueAs("Wb/m")
if equation.IsHarmonic:
self.write.boundary(name, "AV re {e} 3", round(potential, 6))
potential = obj.AV_im_3.getValueAs("Wb/m").Value
self.write.boundary(name, "AV im {e} 3", round(potential, 6))
self.write.boundary(name, "AV re {e} 3", potential)
potential = obj.AV_im_3.getValueAs("Wb/m")
self.write.boundary(name, "AV im {e} 3", potential)
else:
self.write.boundary(name, "AV {e} 3", round(potential, 6))
self.write.boundary(name, "AV {e} 3", potential)
def handleMagnetodynamicBndConditions(self, equation):
# the current density can either be a body force or a boundary constraint

View File

@@ -29,6 +29,7 @@ __url__ = "https://www.freecad.org"
# @{
import collections
from FreeCAD import Units
SIMULATION = "Simulation"
CONSTANTS = "Constants"
@@ -410,22 +411,28 @@ class _Writer:
def _getSifDataType(self, dataType):
if issubclass(dataType, Section):
return _TYPE_INTEGER
if issubclass(dataType, bool):
elif issubclass(dataType, bool):
return _TYPE_LOGICAL
if issubclass(dataType, int):
elif issubclass(dataType, int):
return _TYPE_INTEGER
if issubclass(dataType, float):
elif issubclass(dataType, float):
return _TYPE_REAL
if issubclass(dataType, str):
elif issubclass(dataType, str):
return _TYPE_STRING
raise ValueError("Unsupported data type: %s" % dataType)
elif issubclass(dataType, Units.Quantity):
return _TYPE_REAL
else:
raise ValueError("Unsupported data type: %s" % dataType)
def _preprocess(self, data, dataType):
if issubclass(dataType, Section):
return str(self._idMgr.getId(data))
if issubclass(dataType, str):
elif issubclass(dataType, str):
return '"%s"' % data
return str(data)
elif issubclass(dataType, Units.Quantity):
return f"{data.Value:.13g}"
else:
return str(data)
def _getAttrTypeScalar(self, data):
return self._getSifDataType(type(data))