[FEM] Elmer: add support for 3D magnetodynamics

- adds the corresponding Elmer equation (it is now possible to do Elmer's tutorial example no. 14)

src/Mod/Fem/femsolver/elmer/equations/magnetodynamic.py

@@ -0,0 +1,206 @@
+# ***************************************************************************
+# *   Copyright (c) 2023 Uwe Stöhr <uwestoehr@lyx.org>                      *
+# *                                                                         *
+# *   This file is part of the FreeCAD CAx development system.              *
+# *                                                                         *
+# *   This program is free software; you can redistribute it and/or modify  *
+# *   it under the terms of the GNU Lesser General Public License (LGPL)    *
+# *   as published by the Free Software Foundation; either version 2 of     *
+# *   the License, or (at your option) any later version.                   *
+# *   for detail see the LICENCE text file.                                 *
+# *                                                                         *
+# *   This program 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 Library General Public License for more details.                  *
+# *                                                                         *
+# *   You should have received a copy of the GNU Library General Public     *
+# *   License along with this program; if not, write to the Free Software   *
+# *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  *
+# *   USA                                                                   *
+# *                                                                         *
+# ***************************************************************************
+
+__title__ = "FreeCAD FEM solver Elmer equation object Magnetodynamic"
+__author__ = "Uwe Stöhr"
+__url__ = "https://www.freecadweb.org"
+#
+## \addtogroup FEM
+#  @{
+
+from femtools import femutils
+from . import nonlinear
+from ... import equationbase
+
+def create(doc, name="Magnetodynamic"):
+    return femutils.createObject(
+        doc, name, Proxy, ViewProxy)
+
+
+class Proxy(nonlinear.Proxy, equationbase.MagnetodynamicProxy):
+
+    Type = "Fem::EquationElmerMagnetodynamic"
+
+    def __init__(self, obj):
+        super(Proxy, self).__init__(obj)
+
+        obj.addProperty(
+            "App::PropertyBool",
+            "IsHarmonic",
+            "Magnetodynamic",
+            "If the magnetic source is harmonically driven"
+        )
+        obj.addProperty(
+            "App::PropertyFrequency",
+            "AngularFrequency",
+            "Magnetodynamic",
+            "Frequency of the driving current"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "UsePiolaTransform",
+            "Magnetodynamic",
+            "Must be True if basis functions for edge element interpolation\n"
+            "are selected to be members of optimal edge element family\n"
+            "or if second-order approximation is used."
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "QuadraticApproximation",
+            "Magnetodynamic",
+            "Enables second-order approximation of driving current"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "StaticConductivity",
+            "Magnetodynamic",
+            "See Elmer models manual for info"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "FixInputCurrentDensity",
+            "Magnetodynamic",
+            "Ensures divergence-freeness of current density"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "AutomatedSourceProjectionBCs",
+            "Magnetodynamic",
+            "See Elmer models manual for info"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "UseLagrangeGauge",
+            "Magnetodynamic",
+            "See Elmer models manual for info"
+        )
+        obj.addProperty(
+            "App::PropertyFloat",
+            "LagrangeGaugePenalizationCoefficient",
+            "Magnetodynamic",
+            "See Elmer models manual for info"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "UseTreeGauge",
+            "Magnetodynamic",
+            "See Elmer models manual for info\n"
+            "Will be ignored if 'UsePiolaTransform' is True"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "LinearSystemRefactorize",
+            "Linear System",
+            ""
+        )
+
+        obj.IsHarmonic = False
+        obj.AngularFrequency = 0
+        obj.Priority = 10
+
+        # the post processor options
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateCurrentDensity",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateElectricField",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateElementalFields",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateHarmonicLoss",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateJouleHeating",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateMagneticFieldStrength",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateMaxwellStress",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateNodalFields",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateNodalForces",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "CalculateNodalHeating",
+            "Results",
+            ""
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "DiscontinuousBodies",
+            "Results",
+            ""
+        )
+        obj.CalculateCurrentDensity = False
+        obj.CalculateElectricField = False
+        # FIXME: at the moment FreeCAD's post processor cannot display elementary field
+        # results, therefore disable despite this is by default on in Elmer
+        obj.CalculateElementalFields = False
+        obj.CalculateHarmonicLoss = False
+        obj.CalculateJouleHeating = False
+        obj.CalculateMagneticFieldStrength = False
+        obj.CalculateMaxwellStress = False
+        obj.CalculateNodalFields = True
+        obj.CalculateNodalForces = False
+        obj.CalculateNodalHeating = False
+        obj.DiscontinuousBodies = False
+        
+
+class ViewProxy(nonlinear.ViewProxy, equationbase.MagnetodynamicViewProxy):
+    pass
+
+##  @}

src/Mod/Fem/femsolver/elmer/equations/magnetodynamic_writer.py

@@ -0,0 +1,375 @@
+# ***************************************************************************
+# *   Copyright (c) 2023 Uwe Stöhr <uwestoehr@lyx.org>                      *
+# *                                                                         *
+# *   This file is part of the FreeCAD CAx development system.              *
+# *                                                                         *
+# *   This program is free software; you can redistribute it and/or modify  *
+# *   it under the terms of the GNU Lesser General Public License (LGPL)    *
+# *   as published by the Free Software Foundation; either version 2 of     *
+# *   the License, or (at your option) any later version.                   *
+# *   for detail see the LICENCE text file.                                 *
+# *                                                                         *
+# *   This program 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 Library General Public License for more details.                  *
+# *                                                                         *
+# *   You should have received a copy of the GNU Library General Public     *
+# *   License along with this program; if not, write to the Free Software   *
+# *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  *
+# *   USA                                                                   *
+# *                                                                         *
+# ***************************************************************************
+
+__title__ = "FreeCAD FEM Magnetodynamics Elmer writer"
+__author__ = "Uwe Stöhr"
+__url__ = "https://www.freecad.org"
+
+## \addtogroup FEM
+#  @{
+
+from FreeCAD import Console
+from FreeCAD import Units
+
+from .. import sifio
+from .. import writer as general_writer
+
+class MgDynwriter:
+
+    def __init__(self, writer, solver):
+        self.write = writer
+        self.solver = solver
+
+    def getMagnetodynamicSolver(self, equation):
+        # output the equation parameters
+        s = self.write.createNonlinearSolver(equation)
+        if not equation.IsHarmonic:
+            s["Equation"] = "MgDyn"
+            s["Procedure"] = sifio.FileAttr("MagnetoDynamics/WhitneyAVSolver")
+            s["Variable"] = "av"
+        else:
+            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)
+        s["Exec Solver"] = "Always"
+        s["Optimize Bandwidth"] = True
+        s["Stabilize"] = equation.Stabilize
+        if equation.LinearSystemRefactorize is True:
+            s["Linear System Refactorize"] = True
+        if equation.UsePiolaTransform is True:
+            s["Use Piola Transform"] = True
+        if equation.QuadraticApproximation is True:
+            s["Quadratic Approximation"] = True
+        if equation.StaticConductivity is True:
+            s["Static Conductivity"] = True
+        if equation.FixInputCurrentDensity is True:
+            s["Fix Input Current Density"] = True
+        if equation.AutomatedSourceProjectionBCs is True:
+            s["Automated Source Projection BCs"] = True
+        if equation.UseLagrangeGauge is True:
+            s["Use Lagrange Gauge"] = True
+        if equation.LagrangeGaugePenalizationCoefficient != 0.0:
+            s["Lagrange Gauge Penalization Coefficient"] = \
+                equation.LagrangeGaugePenalizationCoefficient
+        if equation.UseTreeGauge is True:
+            s["Use Tree Gauge"] = True
+        return s
+
+    def getMagnetodynamicSolverPost(self, equation):
+        # output the equation parameters
+        s = self.write.createNonlinearSolver(equation)
+        s["Equation"] = "MgDynPost"
+        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)
+        s["Potential Variable"] = "av"
+        if equation.CalculateCurrentDensity is True:
+            s["Calculate Current Density"] = True
+        if equation.CalculateElectricField is True:
+            s["Calculate Electric Field"] = True
+        if equation.CalculateElementalFields is False:
+            s["Calculate Elemental Fields"] = False
+        if equation.CalculateHarmonicLoss is True:
+            s["Calculate Harmonic Loss"] = True
+        if equation.CalculateJouleHeating is True:
+            s["Calculate Joule Heating"] = True
+        if equation.CalculateMagneticFieldStrength is True:
+            s["Calculate Magnetic Field Strength"] = True
+        if equation.CalculateMaxwellStress is True:
+            s["Calculate Maxwell Stress"] = True
+        if equation.CalculateNodalFields is False:
+            s["Calculate Nodal Fields"] = False
+        if equation.CalculateNodalForces is True:
+            s["Calculate Nodal Forces"] = True
+        if equation.CalculateNodalHeating is True:
+            s["Calculate Nodal Heating"] = True
+        if equation.DiscontinuousBodies is True:
+            s["Discontinuous Bodies"] = True
+        s["Optimize Bandwidth"] = True
+        s["Stabilize"] = equation.Stabilize
+        return s
+
+    def handleMagnetodynamicConstants(self):
+        permeability = self.write.convert(
+            self.write.constsdef["PermeabilityOfVacuum"],
+            "M*L/(T^2*I^2)"
+        )
+        # we round in the following to get rid of numerical artifacts
+        self.write.constant("Permeability Of Vacuum", round(permeability, 20))
+
+        permittivity = self.write.convert(
+            self.write.constsdef["PermittivityOfVacuum"],
+            "T^4*I^2/(L^3*M)"
+        )
+        self.write.constant("Permittivity Of Vacuum", round(permittivity, 20))
+
+    def handleMagnetodynamicMaterial(self, bodies):
+        # check that all bodies have a set material
+        for name in bodies:
+            if self.write.getBodyMaterial(name) == None:
+                raise general_writer.WriteError(
+                    "The body {} is not referenced in any material.\n\n".format(name)
+                )
+        for obj in self.write.getMember("App::MaterialObject"):
+            m = obj.Material
+            refs = (
+                obj.References[0][1]
+                if obj.References
+                else self.write.getAllBodies())
+            for name in (n for n in refs if n in bodies):
+                if "ElectricalConductivity" not in m:
+                    Console.PrintMessage("m: {}\n".format(m))
+                    raise general_writer.WriteError(
+                        "The electrical conductivity must be specified for all materials.\n\n"
+                    )
+                if "RelativePermeability" not in m:
+                    Console.PrintMessage("m: {}\n".format(m))
+                    raise general_writer.WriteError(
+                        "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
+                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:
+                    self.write.material(
+                        name, "Relative Permittivity",
+                        float(m["RelativePermittivity"])
+                    )
+
+    def _outputMagnetodynamicBodyForce(self, obj, name, equation):
+        if hasattr(obj, "CurrentDensity_re_1"):
+            # output only if current density is enabled and needed
+            if not obj.CurrentDensity_re_1_Disabled:
+                currentDensity = float(obj.CurrentDensity_re_1.getValueAs("A/m^2"))
+                self.write.bodyForce(name, "Current Density 1", round(currentDensity, 6))
+            if not obj.CurrentDensity_re_2_Disabled:
+                currentDensity = float(obj.CurrentDensity_re_2.getValueAs("A/m^2"))
+                self.write.bodyForce(name, "Current Density 2", round(currentDensity, 6))
+            if not obj.CurrentDensity_re_3_Disabled:
+                currentDensity = float(obj.CurrentDensity_re_3.getValueAs("A/m^2"))
+                self.write.bodyForce(name, "Current Density 3", round(currentDensity, 6))
+            # imaginaries are only needed for harmonic equation
+            if equation.IsHarmonic:
+                if not obj.CurrentDensity_im_1_Disabled:
+                    currentDensity = float(obj.CurrentDensity_im_1.getValueAs("A/m^2"))
+                    self.write.bodyForce(name, "Current Density Im 1", round(currentDensity, 6))
+                if not obj.CurrentDensity_im_2_Disabled:
+                    currentDensity = float(obj.CurrentDensity_im_2.getValueAs("A/m^2"))
+                    self.write.bodyForce(name, "Current Density Im 2", round(currentDensity, 6))
+                if not obj.CurrentDensity_im_3_Disabled:
+                    currentDensity = float(obj.CurrentDensity_im_3.getValueAs("A/m^2"))
+                    self.write.bodyForce(name, "Current Density Im 3", round(currentDensity, 6))
+
+        if hasattr(obj, "Magnetization_im_1"):
+            # output only if magnetization is enabled and needed
+            if not obj.Magnetization_re_1_Disabled:
+                magnetization = float(obj.Magnetization_re_1.getValueAs("A/m"))
+                self.write.bodyForce(name, "Magnetization 1", magnetization)
+            if not obj.Magnetization_re_2_Disabled:
+                magnetization = float(obj.Magnetization_re_2.getValueAs("A/m"))
+                self.write.bodyForce(name, "Magnetization 2", magnetization)
+            if not obj.Magnetization_re_3_Disabled:
+                magnetization = float(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 not obj.Magnetization_im_1_Disabled:
+                    magnetization = float(obj.Magnetization_im_1.getValueAs("A/m"))
+                    self.write.bodyForce(name, "Magnetization Im 1", magnetization)
+                if not obj.Magnetization_im_2_Disabled:
+                    magnetization = float(obj.Magnetization_im_2.getValueAs("A/m"))
+                    self.write.bodyForce(name, "Magnetization Im 2", magnetization)
+                if not obj.Magnetization_im_3_Disabled:
+                    magnetization = float(obj.Magnetization_im_3.getValueAs("A/m"))
+                    self.write.bodyForce(name, "Magnetization Im 3", magnetization)
+
+        if hasattr(obj, "PotentialEnabled"):
+            # check for PotentialEnabled not Potential since PotentialEnabled was
+            # added later and only with this the imaginary property is available
+            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))
+            # imaginary is only needed for harmonic equation
+            if equation.IsHarmonic:
+                if not obj.AV_im_Disabled:
+                    potential = float(obj.AV_im.getValueAs("V"))
+                    self.write.bodyForce(name, "Electric Potential Im", round(potential, 6))
+
+    def handleMagnetodynamicBodyForces(self, bodies, equation):
+        # the current density can either be a body force or a boundary constraint
+        # therefore only output here if a solid is referenced
+        currentDensities = self.write.getMember("Fem::ConstraintCurrentDensity")
+        for obj in currentDensities:
+            if obj.References:
+                firstName = obj.References[0][1][0]
+                firstName = firstName.rstrip('0123456789')
+                if firstName == "Solid":
+                    for name in obj.References[0][1]:
+                        self._outputMagnetodynamicBodyForce(obj, name, equation)
+                    self.write.handled(obj)
+            else:
+                # if there is only one current density without a reference,
+                # add it to all bodies
+                if len(currentDensities) == 1:
+                    for name in bodies:
+                        self._outputMagnetodynamicBodyForce(obj, name, equation)
+                else:
+                    raise general_writer.WriteError(
+                        "Several current density constraints found without reference to a body.\n"
+                        "Please set a body for each current density constraint."
+                    )
+            self.write.handled(obj)
+
+        magnetizations = self.write.getMember("Fem::ConstraintMagnetization")
+        for obj in magnetizations:
+            if obj.References:
+                for name in obj.References[0][1]:
+                    self._outputMagnetodynamicBodyForce(obj, name, equation)
+                self.write.handled(obj)
+            else:
+                # if there is only one magnetization without a reference,
+                # add it to all bodies
+                if len(magnetizations) == 1:
+                    for name in bodies:
+                        self._outputMagnetodynamicBodyForce(obj, name, equation)
+                else:
+                    raise general_writer.WriteError(
+                        "Several magnetization constraints found without reference to a body.\n"
+                        "Please set a body for each current density constraint."
+                    )
+            self.write.handled(obj)
+
+        # the potential can either be a body force or a boundary constraint
+        # therefore only output here if a solid is referenced
+        potentials = self.write.getMember("Fem::ConstraintElectrostaticPotential")
+        for obj in potentials:
+            if obj.References:
+                firstName = obj.References[0][1][0]
+                firstName = firstName.rstrip('0123456789')
+                if firstName == "Solid":
+                    for name in obj.References[0][1]:
+                        # output only if potentiual is enabled and needed
+                        self._outputMagnetodynamicBodyForce(obj, name, equation)
+                    self.write.handled(obj)
+
+    def _outputMagnetodynamicBndConditions(self, obj, name, equation):
+        if hasattr(obj, "CurrentDensity_re_1"):
+            # output only if current density is enabled and needed
+            if not obj.CurrentDensity_re_1_Disabled:
+                currentDensity = float(obj.CurrentDensity_re_1.getValueAs("A/m^2"))
+                self.write.boundary(name, "Current Density 1", round(currentDensity, 6))
+            # imaginaries are only needed for harmonic equation
+            if equation.IsHarmonic:
+                if not obj.CurrentDensity_im_1_Disabled:
+                    currentDensity = float(obj.CurrentDensity_im_1.getValueAs("A/m^2"))
+                    self.write.boundary(name, "Current Density Im 1", round(currentDensity, 6))
+
+        if hasattr(obj, "PotentialEnabled"):
+            # check for PotentialEnabled not Potential since PotentialEnabled was
+            # added later and only with this the vectorial properties are available
+            if obj.PotentialEnabled:
+                potential = float(obj.Potential.getValueAs("V"))
+                if equation.IsHarmonic:
+                    self.write.boundary(name, "AV re", round(potential, 6))
+                else:
+                    self.write.boundary(name, "AV", round(potential, 6))
+            if not obj.AV_re_1_Disabled:
+                potential = float(obj.AV_re_1.getValueAs("V"))
+                if equation.IsHarmonic:
+                    self.write.boundary(name, "AV re {e} 1", round(potential, 6))
+                else:
+                    self.write.boundary(name, "AV {e} 1", round(potential, 6))
+            if not obj.AV_re_2_Disabled:
+                potential = float(obj.AV_re_2.getValueAs("V"))
+                if equation.IsHarmonic:
+                    self.write.boundary(name, "AV re {e} 2", round(potential, 6))
+                else:
+                    self.write.boundary(name, "AV {e} 2", round(potential, 6))
+            if not obj.AV_re_3_Disabled:
+                potential = float(obj.AV_re_3.getValueAs("V"))
+                if equation.IsHarmonic:
+                    self.write.boundary(name, "AV re {e} 3", round(potential, 6))
+                else:
+                    self.write.boundary(name, "AV {e} 3", round(potential, 6))
+            # imaginaries are only needed for harmonic equation
+            if equation.IsHarmonic:
+                if not obj.AV_im_Disabled:
+                    potential = float(obj.AV_im.getValueAs("V"))
+                    self.write.boundary(name, "AV im", round(potential, 6))
+                if not obj.AV_im_1_Disabled:
+                    potential = float(obj.AV_im_1.getValueAs("V"))
+                    self.write.boundary(name, "AV im {e} 1", round(potential, 6))
+                if not obj.AV_im_2_Disabled:
+                    potential = float(obj.AV_im_2.getValueAs("V"))
+                    self.write.boundary(name, "AV im {e} 2", round(potential, 6))
+                if not obj.AV_im_3_Disabled:
+                    potential = float(obj.AV_im_3.getValueAs("V"))
+                    self.write.boundary(name, "AV im {e} 3", round(potential, 6))
+
+    def handleMagnetodynamicBndConditions(self, equation):
+        # the current density can either be a body force or a boundary constraint
+        # therefore only output here if a face is referenced
+        currentDensities = self.write.getMember("Fem::ConstraintCurrentDensity")
+        for obj in currentDensities:
+            if obj.References:
+                firstName = obj.References[0][1][0]
+                firstName = firstName.rstrip('0123456789')
+                if firstName == "Face":
+                    for name in obj.References[0][1]:
+                        self._outputMagnetodynamicBndConditions(obj, name, equation)
+                    self.write.handled(obj)
+        
+        # the potential can either be a body force or a boundary constraint
+        # therefore only output here if a face is referenced
+        potentials = self.write.getMember("Fem::ConstraintElectrostaticPotential")
+        if len(potentials) == 0:
+            raise general_writer.WriteError(
+                "The Magnetodynamic equation needs at least one ElectrostaticPotential"
+                "constraint."
+            )
+        for obj in potentials:
+            if obj.References:
+                firstName = obj.References[0][1][0]
+                firstName = firstName.rstrip('0123456789')
+                if firstName == "Face":
+                    for name in obj.References[0][1]:
+                        # output the FreeCAD label as comment
+                        if obj.Label:
+                            self.write.boundary(name, "! FreeCAD Name", obj.Label)
+                        # output only if potentiual is enabled and needed
+                        self._outputMagnetodynamicBndConditions(obj, name, equation)
+                    self.write.handled(obj)
+
+##  @}

src/Mod/Fem/femsolver/elmer/writer.py

@@ -54,6 +54,7 @@ from .equations import electrostatic_writer as ES_writer
 from .equations import flow_writer
 from .equations import flux_writer
 from .equations import heat_writer
+from .equations import magnetodynamic_writer as MgDyn_writer
 from .equations import magnetodynamic2D_writer as MgDyn2D_writer
 
 
@@ -99,6 +100,7 @@ class Writer(object):
         self._handleHeat()
         self._handleFlow()
         self._handleFlux()
+        self._handleMagnetodynamic()
         self._handleMagnetodynamic2D()
         self._addOutputSolver()
 
@@ -534,6 +536,30 @@ class Writer(object):
             HeatW.handleHeatBodyForces(activeIn)
             HeatW.handleHeatMaterial(activeIn)
 
+    #-------------------------------------------------------------------------------------------
+    # Magnetodynamic
+
+    def _handleMagnetodynamic(self):
+        MgDyn = MgDyn_writer.MgDynwriter(self, self.solver)
+        activeIn = []
+        for equation in self.solver.Group:
+            if femutils.is_of_type(equation, "Fem::EquationElmerMagnetodynamic"):
+                if equation.References:
+                    activeIn = equation.References[0][1]
+                else:
+                    activeIn = self.getAllBodies()
+                    
+                solverSection = MgDyn.getMagnetodynamicSolver(equation)
+                solverPostSection = MgDyn.getMagnetodynamicSolverPost(equation)
+                for body in activeIn:
+                    self._addSolver(body, solverSection)
+                    self._addSolver(body, solverPostSection)
+        if activeIn:
+            MgDyn.handleMagnetodynamicConstants()
+            MgDyn.handleMagnetodynamicBndConditions(equation)
+            MgDyn.handleMagnetodynamicBodyForces(activeIn, equation)
+            MgDyn.handleMagnetodynamicMaterial(activeIn)
+
     #-------------------------------------------------------------------------------------------
     # Magnetodynamic2D
 

src/Mod/Fem/femsolver/equationbase.py

@@ -129,6 +129,16 @@ class HeatViewProxy(BaseViewProxy):
         return ":/icons/FEM_EquationHeat.svg"
 
 
+class MagnetodynamicProxy(BaseProxy):
+    pass
+
+
+class MagnetodynamicViewProxy(BaseViewProxy):
+
+    def getIcon(self):
+        return ":/icons/FEM_EquationMagnetodynamic.svg"
+
+
 class Magnetodynamic2DProxy(BaseProxy):
     pass