From 8d0d38911ffb9850690e42c4cec63a89e41bfa12 Mon Sep 17 00:00:00 2001 From: Uwe Date: Tue, 7 Feb 2023 02:54:23 +0100 Subject: [PATCH] [FEM] Elmer writer: sort out flow equation - sort out flow equation - next step to refactor writer.py --- src/Mod/Fem/CMakeLists.txt | 1 + .../femsolver/elmer/equations/flow_writer.py | 264 ++++++++++++++++++ src/Mod/Fem/femsolver/elmer/writer.py | 254 ++--------------- 3 files changed, 281 insertions(+), 238 deletions(-) create mode 100644 src/Mod/Fem/femsolver/elmer/equations/flow_writer.py diff --git a/src/Mod/Fem/CMakeLists.txt b/src/Mod/Fem/CMakeLists.txt index b82d968c81..735ff3da1c 100755 --- a/src/Mod/Fem/CMakeLists.txt +++ b/src/Mod/Fem/CMakeLists.txt @@ -257,6 +257,7 @@ SET(FemSolverElmerEquations_SRCS femsolver/elmer/equations/electrostatic_writer.py femsolver/elmer/equations/equation.py femsolver/elmer/equations/flow.py + femsolver/elmer/equations/flow_writer.py femsolver/elmer/equations/flux.py femsolver/elmer/equations/flux_writer.py femsolver/elmer/equations/heat.py diff --git a/src/Mod/Fem/femsolver/elmer/equations/flow_writer.py b/src/Mod/Fem/femsolver/elmer/equations/flow_writer.py new file mode 100644 index 0000000000..786e970af9 --- /dev/null +++ b/src/Mod/Fem/femsolver/elmer/equations/flow_writer.py @@ -0,0 +1,264 @@ +# *************************************************************************** +# * Copyright (c) 2023 Uwe Stöhr * +# * * +# * 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 Flow Elmer writer" +__author__ = "Uwe Stöhr" +__url__ = "https://www.freecad.org" + +## \addtogroup FEM +# @{ + +from .. import sifio +from .. import writer as general_writer +from . import flow + +class Flowwriter: + + def __init__(self, writer, solver): + self.write = writer + self.solver = solver + + def getFlowSolver(self, equation): + # check if we need to update the equation + self._updateFlowSolver(equation) + # output the equation parameters + s = self.write.createNonlinearSolver(equation) + s["Equation"] = "Navier-Stokes" + s["Procedure"] = sifio.FileAttr("FlowSolve/FlowSolver") + if equation.DivDiscretization is True: + s["Div Discretization"] = equation.DivDiscretization + s["Exec Solver"] = "Always" + if equation.FlowModel != "Full": + s["Flow Model"] = equation.FlowModel + if equation.GradpDiscretization is True: + s["Gradp Discretization"] = equation.GradpDiscretization + s["Stabilize"] = equation.Stabilize + s["Optimize Bandwidth"] = True + if equation.Variable != "Flow Solution[Velocity:3 Pressure:1]": + s["Variable"] = equation.Variable + return s + + def handleFlowConstants(self): + gravity = self.write.convert(self.write.constsdef["Gravity"], "L/T^2") + self.write.constant("Gravity", (0.0, -1.0, 0.0, gravity)) + + def _updateFlowSolver(self, equation): + # updates older Flow equations + if not hasattr(equation, "Convection"): + equation.addProperty( + "App::PropertyEnumeration", + "Convection", + "Equation", + "Type of convection to be used" + ) + equation.Convection = flow.CONVECTION_TYPE + equation.Convection = "Computed" + if not hasattr(equation, "DivDiscretization"): + equation.addProperty( + "App::PropertyBool", + "DivDiscretization", + "Flow", + ( + "Set to true for incompressible flow for more stable\n" + "discretization when Reynolds number increases" + ) + ) + if not hasattr(equation, "FlowModel"): + equation.addProperty( + "App::PropertyEnumeration", + "FlowModel", + "Flow", + "Flow model to be used" + ) + equation.FlowModel = flow.FLOW_MODEL + equation.FlowModel = "Full" + if not hasattr(equation, "GradpDiscretization"): + equation.addProperty( + "App::PropertyBool", + "GradpDiscretization", + "Flow", + ( + "If true pressure Dirichlet boundary conditions can be used.\n" + "Also mass flux is available as a natural boundary condition." + ) + ) + if not hasattr(equation, "MagneticInduction"): + equation.addProperty( + "App::PropertyBool", + "MagneticInduction", + "Equation", + ( + "Magnetic induction equation will be solved\n" + "along with the Navier-Stokes equations" + ) + ) + if not hasattr(equation, "Variable"): + equation.addProperty( + "App::PropertyString", + "Variable", + "Flow", + "Only for a 2D model change the '3' to '2'" + ) + equation.Variable = "Flow Solution[Velocity:3 Pressure:1]" + + def handleFlowMaterial(self, bodies): + tempObj = self.write.getSingleMember("Fem::ConstraintInitialTemperature") + if tempObj is not None: + refTemp = self.write.getFromUi(tempObj.initialTemperature, "K", "O") + for name in bodies: + self.write.material(name, "Reference Temperature", refTemp) + 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): + self.write.material(name, "Name", m["Name"]) + if "Density" in m: + self.write.material( + name, "Density", + self.write.getDensity(m) + ) + if "ThermalConductivity" in m: + self.write.material( + name, "Heat Conductivity", + self.write.convert(m["ThermalConductivity"], "M*L/(T^3*O)") + ) + if "KinematicViscosity" in m: + density = self.write.getDensity(m) + kViscosity = self.write.convert(m["KinematicViscosity"], "L^2/T") + self.write.material( + name, "Viscosity", kViscosity * density) + if "ThermalExpansionCoefficient" in m: + value = self.write.convert(m["ThermalExpansionCoefficient"], "O^-1") + if value > 0: + self.write.material( + name, "Heat expansion Coefficient", value) + if "ReferencePressure" in m: + pressure = self.write.convert(m["ReferencePressure"], "M/(L*T^2)") + self.write.material(name, "Reference Pressure", pressure) + if "SpecificHeatRatio" in m: + self.write.material( + name, "Specific Heat Ratio", + float(m["SpecificHeatRatio"]) + ) + if "CompressibilityModel" in m: + self.write.material( + name, "Compressibility Model", + m["CompressibilityModel"]) + + def _outputInitialPressure(self, obj, name): + # initial pressure only makes sense for fluid material + if self.write.isBodyMaterialFluid(name): + pressure = float(obj.Pressure.getValueAs("Pa")) + self.write.initial(name, "Pressure", pressure) + + def handleFlowInitialPressure(self, bodies): + initialPressures = self.write.getMember("Fem::ConstraintInitialPressure") + for obj in initialPressures: + if obj.References: + for name in obj.References[0][1]: + self._outputInitialPressure(obj, name) + self.write.handled(obj) + else: + # if there is only one initial pressure without a reference + # add it to all fluid bodies + if len(initialPressures) == 1: + for name in bodies: + self._outputInitialPressure(obj, name) + else: + raise general_writer.WriteError( + "Several initial pressures found without reference to a body.\n" + "Please set a body for each initial pressure." + ) + self.write.handled(obj) + + def _outputInitialVelocity(self, obj, name): + # flow only makes sense for fluid material + if self.write.isBodyMaterialFluid(name): + if obj.VelocityXEnabled: + velocity = self.write.getFromUi(obj.VelocityX, "m/s", "L/T") + self.write.initial(name, "Velocity 1", velocity) + if obj.VelocityYEnabled: + velocity = self.write.getFromUi(obj.VelocityY, "m/s", "L/T") + self.write.initial(name, "Velocity 2", velocity) + if obj.VelocityZEnabled: + velocity = self.write.getFromUi(obj.VelocityZ, "m/s", "L/T") + self.write.initial(name, "Velocity 3", velocity) + + def handleFlowInitialVelocity(self, bodies): + initialVelocities = self.write.getMember("Fem::ConstraintInitialFlowVelocity") + for obj in initialVelocities: + if obj.References: + for name in obj.References[0][1]: + self._outputInitialVelocity(obj, name) + self.write.handled(obj) + else: + # if there is only one initial velocity without a reference + # add it to all fluid bodies + if len(initialVelocities) == 1: + for name in bodies: + self._outputInitialVelocity(obj, name) + else: + raise general_writer.WriteError( + "Several initial velocities found without reference to a body.\n" + "Please set a body for each initial velocity." + ) + self.write.handled(obj) + + def handleFlowBndConditions(self): + for obj in self.write.getMember("Fem::ConstraintFlowVelocity"): + if obj.References: + for name in obj.References[0][1]: + if obj.VelocityXEnabled: + velocity = self.write.getFromUi(obj.VelocityX, "m/s", "L/T") + self.write.boundary(name, "Velocity 1", velocity) + if obj.VelocityYEnabled: + velocity = self.write.getFromUi(obj.VelocityY, "m/s", "L/T") + self.write.boundary(name, "Velocity 2", velocity) + if obj.VelocityZEnabled: + velocity = self.write.getFromUi(obj.VelocityZ, "m/s", "L/T") + self.write.boundary(name, "Velocity 3", velocity) + if obj.NormalToBoundary: + self.write.boundary(name, "Normal-Tangential Velocity", True) + self.write.handled(obj) + for obj in self.write.getMember("Fem::ConstraintPressure"): + if obj.References: + for name in obj.References[0][1]: + pressure = self.write.getFromUi(obj.Pressure, "MPa", "M/(L*T^2)") + if obj.Reversed: + pressure *= -1 + self.write.boundary(name, "External Pressure", pressure) + self.write.handled(obj) + + def handleFlowEquation(self, bodies, equation): + for b in bodies: + # not for bodies with solid material + if self.write.isBodyMaterialFluid(b): + if equation.Convection != "None": + self.write.equation(b, "Convection", equation.Convection) + if equation.MagneticInduction is True: + self.write.equation(b, "Magnetic Induction", equation.MagneticInduction) + +## @} diff --git a/src/Mod/Fem/femsolver/elmer/writer.py b/src/Mod/Fem/femsolver/elmer/writer.py index 06108f15bf..8592bf615c 100644 --- a/src/Mod/Fem/femsolver/elmer/writer.py +++ b/src/Mod/Fem/femsolver/elmer/writer.py @@ -51,7 +51,7 @@ from femtools import membertools from .equations import elasticity from .equations import electricforce from .equations import electrostatic_writer as ES_writer -from .equations import flow +from .equations import flow_writer from .equations import flux_writer from .equations import heat_writer @@ -416,7 +416,7 @@ class Writer(object): activeIn = self.getAllBodies() solverSection = self._getElasticitySolver(equation) for body in activeIn: - if not self._isBodyMaterialFluid(body): + if not self.isBodyMaterialFluid(body): self._addSolver(body, solverSection) self._handleElasticityEquation(activeIn, equation) if activeIn: @@ -479,9 +479,9 @@ class Writer(object): def _handleElasticityEquation(self, bodies, equation): for b in bodies: # not for bodies with fluid material - if not self._isBodyMaterialFluid(b): + if not self.isBodyMaterialFluid(b): if equation.PlaneStress: - self._equation(b, "Plane Stress", equation.PlaneStress) + self.equation(b, "Plane Stress", equation.PlaneStress) def _updateElasticitySolver(self, equation): # updates older Elasticity equations @@ -793,7 +793,7 @@ class Writer(object): ) for name in (n for n in refs if n in bodies): # don't evaluate fluid material - if self._isBodyMaterialFluid(name): + if self.isBodyMaterialFluid(name): break if "YoungsModulus" not in m: Console.PrintMessage("m: {}\n".format(m)) @@ -895,6 +895,7 @@ class Writer(object): # Flow def _handleFlow(self): + FlowW = flow_writer.Flowwriter(self, self.solver) activeIn = [] for equation in self.solver.Group: if femutils.is_of_type(equation, "Fem::EquationElmerFlow"): @@ -906,240 +907,17 @@ class Writer(object): activeIn = equation.References[0][1] else: activeIn = self.getAllBodies() - solverSection = self._getFlowSolver(equation) + solverSection = FlowW.getFlowSolver(equation) for body in activeIn: - if self._isBodyMaterialFluid(body): + if self.isBodyMaterialFluid(body): self._addSolver(body, solverSection) - self._handleFlowEquation(activeIn, equation) + FlowW.handleFlowEquation(activeIn, equation) if activeIn: - self._handleFlowConstants() - self._handleFlowBndConditions() - self._handleFlowInitialPressure(activeIn) - self._handleFlowInitialVelocity(activeIn) - self._handleFlowMaterial(activeIn) - - def _getFlowSolver(self, equation): - # check if we need to update the equation - self._updateFlowSolver(equation) - # output the equation parameters - s = self.createNonlinearSolver(equation) - s["Equation"] = "Navier-Stokes" - s["Procedure"] = sifio.FileAttr("FlowSolve/FlowSolver") - if equation.DivDiscretization is True: - s["Div Discretization"] = equation.DivDiscretization - s["Exec Solver"] = "Always" - if equation.FlowModel != "Full": - s["Flow Model"] = equation.FlowModel - if equation.GradpDiscretization is True: - s["Gradp Discretization"] = equation.GradpDiscretization - s["Stabilize"] = equation.Stabilize - s["Optimize Bandwidth"] = True - if equation.Variable != "Flow Solution[Velocity:3 Pressure:1]": - s["Variable"] = equation.Variable - return s - - def _handleFlowConstants(self): - gravity = self.convert(self.constsdef["Gravity"], "L/T^2") - self.constant("Gravity", (0.0, -1.0, 0.0, gravity)) - - def _updateFlowSolver(self, equation): - # updates older Flow equations - if not hasattr(equation, "Convection"): - equation.addProperty( - "App::PropertyEnumeration", - "Convection", - "Equation", - "Type of convection to be used" - ) - equation.Convection = flow.CONVECTION_TYPE - equation.Convection = "Computed" - if not hasattr(equation, "DivDiscretization"): - equation.addProperty( - "App::PropertyBool", - "DivDiscretization", - "Flow", - ( - "Set to true for incompressible flow for more stable\n" - "discretization when Reynolds number increases" - ) - ) - if not hasattr(equation, "FlowModel"): - equation.addProperty( - "App::PropertyEnumeration", - "FlowModel", - "Flow", - "Flow model to be used" - ) - equation.FlowModel = flow.FLOW_MODEL - equation.FlowModel = "Full" - if not hasattr(equation, "GradpDiscretization"): - equation.addProperty( - "App::PropertyBool", - "GradpDiscretization", - "Flow", - ( - "If true pressure Dirichlet boundary conditions can be used.\n" - "Also mass flux is available as a natural boundary condition." - ) - ) - if not hasattr(equation, "MagneticInduction"): - equation.addProperty( - "App::PropertyBool", - "MagneticInduction", - "Equation", - ( - "Magnetic induction equation will be solved\n" - "along with the Navier-Stokes equations" - ) - ) - if not hasattr(equation, "Variable"): - equation.addProperty( - "App::PropertyString", - "Variable", - "Flow", - "Only for a 2D model change the '3' to '2'" - ) - equation.Variable = "Flow Solution[Velocity:3 Pressure:1]" - - def _handleFlowMaterial(self, bodies): - tempObj = self.getSingleMember("Fem::ConstraintInitialTemperature") - if tempObj is not None: - refTemp = self.getFromUi(tempObj.initialTemperature, "K", "O") - for name in bodies: - self.material(name, "Reference Temperature", refTemp) - for obj in self.getMember("App::MaterialObject"): - m = obj.Material - refs = ( - obj.References[0][1] - if obj.References - else self.getAllBodies()) - for name in (n for n in refs if n in bodies): - self.material(name, "Name", m["Name"]) - if "Density" in m: - self.material( - name, "Density", - self.getDensity(m) - ) - if "ThermalConductivity" in m: - self.material( - name, "Heat Conductivity", - self.convert(m["ThermalConductivity"], "M*L/(T^3*O)") - ) - if "KinematicViscosity" in m: - density = self.getDensity(m) - kViscosity = self.convert(m["KinematicViscosity"], "L^2/T") - self.material( - name, "Viscosity", kViscosity * density) - if "ThermalExpansionCoefficient" in m: - value = self.convert(m["ThermalExpansionCoefficient"], "O^-1") - if value > 0: - self.material( - name, "Heat expansion Coefficient", value) - if "ReferencePressure" in m: - pressure = self.convert(m["ReferencePressure"], "M/(L*T^2)") - self.material(name, "Reference Pressure", pressure) - if "SpecificHeatRatio" in m: - self.material( - name, "Specific Heat Ratio", - float(m["SpecificHeatRatio"]) - ) - if "CompressibilityModel" in m: - self.material( - name, "Compressibility Model", - m["CompressibilityModel"]) - - def _outputInitialPressure(self, obj, name): - # initial pressure only makes sense for fluid material - if self._isBodyMaterialFluid(name): - pressure = float(obj.Pressure.getValueAs("Pa")) - self.initial(name, "Pressure", pressure) - - def _handleFlowInitialPressure(self, bodies): - initialPressures = self.getMember("Fem::ConstraintInitialPressure") - for obj in initialPressures: - if obj.References: - for name in obj.References[0][1]: - self._outputInitialPressure(obj, name) - self.handled(obj) - else: - # if there is only one initial pressure without a reference - # add it to all fluid bodies - if len(initialPressures) == 1: - for name in bodies: - self._outputInitialPressure(obj, name) - else: - raise WriteError( - "Several initial pressures found without reference to a body.\n" - "Please set a body for each initial pressure." - ) - self.handled(obj) - - def _outputInitialVelocity(self, obj, name): - # flow only makes sense for fluid material - if self._isBodyMaterialFluid(name): - if obj.VelocityXEnabled: - velocity = self.getFromUi(obj.VelocityX, "m/s", "L/T") - self.initial(name, "Velocity 1", velocity) - if obj.VelocityYEnabled: - velocity = self.getFromUi(obj.VelocityY, "m/s", "L/T") - self.initial(name, "Velocity 2", velocity) - if obj.VelocityZEnabled: - velocity = self.getFromUi(obj.VelocityZ, "m/s", "L/T") - self.initial(name, "Velocity 3", velocity) - - def _handleFlowInitialVelocity(self, bodies): - initialVelocities = self.getMember("Fem::ConstraintInitialFlowVelocity") - for obj in initialVelocities: - if obj.References: - for name in obj.References[0][1]: - self._outputInitialVelocity(obj, name) - self.handled(obj) - else: - # if there is only one initial velocity without a reference - # add it to all fluid bodies - if len(initialVelocities) == 1: - for name in bodies: - self._outputInitialVelocity(obj, name) - else: - raise WriteError( - "Several initial velocities found without reference to a body.\n" - "Please set a body for each initial velocity." - ) - self.handled(obj) - - def _handleFlowBndConditions(self): - for obj in self.getMember("Fem::ConstraintFlowVelocity"): - if obj.References: - for name in obj.References[0][1]: - if obj.VelocityXEnabled: - velocity = self.getFromUi(obj.VelocityX, "m/s", "L/T") - self.boundary(name, "Velocity 1", velocity) - if obj.VelocityYEnabled: - velocity = self.getFromUi(obj.VelocityY, "m/s", "L/T") - self.boundary(name, "Velocity 2", velocity) - if obj.VelocityZEnabled: - velocity = self.getFromUi(obj.VelocityZ, "m/s", "L/T") - self.boundary(name, "Velocity 3", velocity) - if obj.NormalToBoundary: - self.boundary(name, "Normal-Tangential Velocity", True) - self.handled(obj) - for obj in self.getMember("Fem::ConstraintPressure"): - if obj.References: - for name in obj.References[0][1]: - pressure = self.getFromUi(obj.Pressure, "MPa", "M/(L*T^2)") - if obj.Reversed: - pressure *= -1 - self.boundary(name, "External Pressure", pressure) - self.handled(obj) - - def _handleFlowEquation(self, bodies, equation): - for b in bodies: - # not for bodies with solid material - if self._isBodyMaterialFluid(b): - if equation.Convection != "None": - self._equation(b, "Convection", equation.Convection) - if equation.MagneticInduction is True: - self._equation(b, "Magnetic Induction", equation.MagneticInduction) + FlowW.handleFlowConstants() + FlowW.handleFlowBndConditions() + FlowW.handleFlowInitialPressure(activeIn) + FlowW.handleFlowInitialVelocity(activeIn) + FlowW.handleFlowMaterial(activeIn) #------------------------------------------------------------------------------------------- # Flux @@ -1294,7 +1072,7 @@ class Writer(object): return True return False - def _isBodyMaterialFluid(self, body): + def isBodyMaterialFluid(self, body): # we can have the case that a body has no assigned material # then assume it is a solid if self._getBodyMaterial(body) is not None: @@ -1390,7 +1168,7 @@ class Writer(object): def material(self, body, key, attr): self._builder.material(body, key, attr) - def _equation(self, body, key, attr): + def equation(self, body, key, attr): self._builder.equation(body, key, attr) def _bodyForce(self, body, key, attr):