FEM: move setup flow 1D thermomech example from unit tests to fem examples

This commit is contained in:
Bernd Hahnebach
2019-09-26 07:22:08 +02:00
parent d0422fa91c
commit 587612f21d
3 changed files with 263 additions and 292 deletions

View File

@@ -39,6 +39,7 @@ SET(FemExamples_SRCS
femexamples/manager.py
femexamples/multimaterial_twoboxes.py
femexamples/rc_wall_2d.py
femexamples/thermomech_flow1d.py
femexamples/thermomech_spine.py
)

View File

@@ -0,0 +1,258 @@
# ***************************************************************************
# * Copyright (c) 2019 Bernd Hahnebach <bernd@bimstatik.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. *
# * *
# * 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 Library General Public License for more details. *
# * *
# * You should have received a copy of the GNU Library General Public *
# * License along with FreeCAD; if not, write to the Free Software *
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 *
# * USA *
# * *
# ***************************************************************************
import FreeCAD
import ObjectsFem
import Fem
mesh_name = "Mesh" # needs to be Mesh to work with unit tests
def init_doc(doc=None):
if doc is None:
doc = FreeCAD.newDocument()
return doc
def setup(doc=None, solver="ccxtools"):
# setup model
if doc is None:
doc = init_doc()
p1 = FreeCAD.Vector(0, 0, 50)
p2 = FreeCAD.Vector(0, 0, -50)
p3 = FreeCAD.Vector(0, 0, -4300)
p4 = FreeCAD.Vector(4950, 0, -4300)
p5 = FreeCAD.Vector(5000, 0, -4300)
p6 = FreeCAD.Vector(8535.53, 0, -7835.53)
p7 = FreeCAD.Vector(8569.88, 0, -7870.88)
p8 = FreeCAD.Vector(12105.41, 0, -11406.41)
p9 = FreeCAD.Vector(12140.76, 0, -11441.76)
p10 = FreeCAD.Vector(13908.53, 0, -13209.53)
p11 = FreeCAD.Vector(13943.88, 0, -13244.88)
p12 = FreeCAD.Vector(15046.97, 0, -14347.97)
p13 = FreeCAD.Vector(15046.97, 0, -7947.97)
p14 = FreeCAD.Vector(15046.97, 0, -7847.97)
p15 = FreeCAD.Vector(0, 0, 0)
p16 = FreeCAD.Vector(0, 0, -2175)
p17 = FreeCAD.Vector(2475, 0, -4300)
p18 = FreeCAD.Vector(4975, 0, -4300)
p19 = FreeCAD.Vector(6767.765, 0, -6067.765)
p20 = FreeCAD.Vector(8552.705, 0, -7853.205)
p21 = FreeCAD.Vector(10337.645, 0, -9638.645)
p22 = FreeCAD.Vector(12123.085, 0, -11424.085)
p23 = FreeCAD.Vector(13024.645, 0, -12325.645)
p24 = FreeCAD.Vector(13926.205, 0, -13227.205)
p25 = FreeCAD.Vector(14495.425, 0, -13796.425)
p26 = FreeCAD.Vector(15046.97, 0, -11147.97)
p27 = FreeCAD.Vector(15046.97, 0, -7897.97)
points = [
p1, p2, p3, p4, p5, p6, p7, p8, p9, p10,
p11, p12, p13, p14, p15, p16, p17, p18, p19, p20,
p21, p22, p23, p24, p25, p26, p27
]
from Draft import makeWire
line = makeWire(
points,
closed=False,
face=False,
support=None
)
doc.recompute()
if FreeCAD.GuiUp:
import FreeCADGui
FreeCADGui.ActiveDocument.activeView().viewAxonometric()
FreeCADGui.SendMsgToActiveView("ViewFit")
# analysis
analysis = ObjectsFem.makeAnalysis(doc, "Analysis")
# solver
# TODO How to pass multiple solver for one analysis in one doc
if solver == "calculix":
solver_object = analysis.addObject(
ObjectsFem.makeSolverCalculix(doc, "SolverCalculiX")
)[0]
elif solver == "ccxtools":
solver_object = analysis.addObject(
ObjectsFem.makeSolverCalculixCcxTools(doc, "CalculiXccxTools") # CalculiX
)[0]
solver_object.WorkingDir = u""
if solver == "calculix" or solver == "ccxtools":
solver_object.AnalysisType = "thermomech"
solver_object.GeometricalNonlinearity = "linear"
solver_object.ThermoMechSteadyState = True
solver_object.MatrixSolverType = "default"
solver_object.IterationsThermoMechMaximum = 2000
solver_object.IterationsControlParameterTimeUse = False
# material
material_object = analysis.addObject(
ObjectsFem.makeMaterialFluid(doc, "FluidMaterial")
)[0]
mat = material_object.Material
mat["Name"] = "Water"
mat["Density"] = "998 kg/m^3"
mat["SpecificHeat"] = "4.182 J/kg/K"
mat["DynamicViscosity"] = "1.003e-3 kg/m/s"
mat["VolumetricThermalExpansionCoefficient"] = "2.07e-4 m/m/K"
mat["ThermalConductivity"] = "0.591 W/m/K"
material_object.Material = mat
Flow1d_inlet = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_inlet.SectionType = "Liquid"
Flow1d_inlet.LiquidSectionType = "PIPE INLET"
Flow1d_inlet.InletPressure = 0.1
Flow1d_inlet.References = [(line, "Edge1")]
Flow1d_entrance = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_entrance.SectionType = "Liquid"
Flow1d_entrance.LiquidSectionType = "PIPE ENTRANCE"
Flow1d_entrance.EntrancePipeArea = 31416.00
Flow1d_entrance.EntranceArea = 25133.00
Flow1d_entrance.References = [(line, "Edge2")]
Flow1d_manning = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_manning.SectionType = "Liquid"
Flow1d_manning.LiquidSectionType = "PIPE MANNING"
Flow1d_manning.ManningArea = 31416
Flow1d_manning.ManningRadius = 50
Flow1d_manning.ManningCoefficient = 0.002
Flow1d_manning.References = [(line, "Edge3"), (line, "Edge5")]
Flow1d_bend = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_bend.SectionType = "Liquid"
Flow1d_bend.LiquidSectionType = "PIPE BEND"
Flow1d_bend.BendPipeArea = 31416
Flow1d_bend.BendRadiusDiameter = 1.5
Flow1d_bend.BendAngle = 45
Flow1d_bend.BendLossCoefficient = 0.4
Flow1d_bend.References = [(line, "Edge4")]
Flow1d_enlargement = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_enlargement.SectionType = "Liquid"
Flow1d_enlargement.LiquidSectionType = "PIPE ENLARGEMENT"
Flow1d_enlargement.EnlargeArea1 = 31416.00
Flow1d_enlargement.EnlargeArea2 = 70686.00
Flow1d_enlargement.References = [(line, "Edge6")]
Flow1d_manning1 = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_manning1.SectionType = "Liquid"
Flow1d_manning1.LiquidSectionType = "PIPE MANNING"
Flow1d_manning1.ManningArea = 70686.00
Flow1d_manning1.ManningRadius = 75
Flow1d_manning1.ManningCoefficient = 0.002
Flow1d_manning1.References = [(line, "Edge7")]
Flow1d_contraction = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_contraction.SectionType = "Liquid"
Flow1d_contraction.LiquidSectionType = "PIPE CONTRACTION"
Flow1d_contraction.ContractArea1 = 70686
Flow1d_contraction.ContractArea2 = 17671
Flow1d_contraction.References = [(line, "Edge8")]
Flow1d_manning2 = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_manning2.SectionType = "Liquid"
Flow1d_manning2.LiquidSectionType = "PIPE MANNING"
Flow1d_manning2.ManningArea = 17671.00
Flow1d_manning2.ManningRadius = 37.5
Flow1d_manning2.ManningCoefficient = 0.002
Flow1d_manning2.References = [(line, "Edge11"), (line, "Edge9")]
Flow1d_gate_valve = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_gate_valve.SectionType = "Liquid"
Flow1d_gate_valve.LiquidSectionType = "PIPE GATE VALVE"
Flow1d_gate_valve.GateValvePipeArea = 17671
Flow1d_gate_valve.GateValveClosingCoeff = 0.5
Flow1d_gate_valve.References = [(line, "Edge10")]
Flow1d_enlargement1 = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_enlargement1.SectionType = "Liquid"
Flow1d_enlargement1.LiquidSectionType = "PIPE ENLARGEMENT"
Flow1d_enlargement1.EnlargeArea1 = 17671
Flow1d_enlargement1.EnlargeArea2 = 1e12
Flow1d_enlargement1.References = [(line, "Edge12")]
Flow1d_outlet = analysis.addObject(
ObjectsFem.makeElementFluid1D(doc, "ElementFluid1D")
)[0]
Flow1d_outlet.SectionType = "Liquid"
Flow1d_outlet.LiquidSectionType = "PIPE OUTLET"
Flow1d_outlet.OutletPressure = 0.1
Flow1d_outlet.References = [(line, "Edge13")]
Flow1d_self_weight = analysis.addObject(
ObjectsFem.makeConstraintSelfWeight(doc, "ConstraintSelfWeight")
)[0]
Flow1d_self_weight.Gravity_x = 0.0
Flow1d_self_weight.Gravity_y = 0.0
Flow1d_self_weight.Gravity_z = -1.0
# mesh
# from femexamples.meshes.mesh_thermomech_spine import create_nodes, create_elements
from femtest.data.ccx.Flow1D_mesh import create_nodes_Flow1D as create_nodes
from femtest.data.ccx.Flow1D_mesh import create_elements_Flow1D as create_elements
fem_mesh = Fem.FemMesh()
control = create_nodes(fem_mesh)
if not control:
FreeCAD.Console.PrintError("Error on creating nodes.\n")
control = create_elements(fem_mesh)
if not control:
FreeCAD.Console.PrintError("Error on creating elements.\n")
femmesh_obj = analysis.addObject(
doc.addObject("Fem::FemMeshObject", mesh_name)
)[0]
femmesh_obj.FemMesh = fem_mesh
doc.recompute()
return doc
"""
from femexamples import thermomech_flow1d as flow
flow.setup()
"""

View File

@@ -608,298 +608,10 @@ class TestCcxTools(unittest.TestCase):
):
fcc_print("--------------- Start of 1D Flow FEM tests ---------------")
p1 = FreeCAD.Vector(0, 0, 50)
p2 = FreeCAD.Vector(0, 0, -50)
p3 = FreeCAD.Vector(0, 0, -4300)
p4 = FreeCAD.Vector(4950, 0, -4300)
p5 = FreeCAD.Vector(5000, 0, -4300)
p6 = FreeCAD.Vector(8535.53, 0, -7835.53)
p7 = FreeCAD.Vector(8569.88, 0, -7870.88)
p8 = FreeCAD.Vector(12105.41, 0, -11406.41)
p9 = FreeCAD.Vector(12140.76, 0, -11441.76)
p10 = FreeCAD.Vector(13908.53, 0, -13209.53)
p11 = FreeCAD.Vector(13943.88, 0, -13244.88)
p12 = FreeCAD.Vector(15046.97, 0, -14347.97)
p13 = FreeCAD.Vector(15046.97, 0, -7947.97)
p14 = FreeCAD.Vector(15046.97, 0, -7847.97)
p15 = FreeCAD.Vector(0, 0, 0)
p16 = FreeCAD.Vector(0, 0, -2175)
p17 = FreeCAD.Vector(2475, 0, -4300)
p18 = FreeCAD.Vector(4975, 0, -4300)
p19 = FreeCAD.Vector(6767.765, 0, -6067.765)
p20 = FreeCAD.Vector(8552.705, 0, -7853.205)
p21 = FreeCAD.Vector(10337.645, 0, -9638.645)
p22 = FreeCAD.Vector(12123.085, 0, -11424.085)
p23 = FreeCAD.Vector(13024.645, 0, -12325.645)
p24 = FreeCAD.Vector(13926.205, 0, -13227.205)
p25 = FreeCAD.Vector(14495.425, 0, -13796.425)
p26 = FreeCAD.Vector(15046.97, 0, -11147.97)
p27 = FreeCAD.Vector(15046.97, 0, -7897.97)
points = [
p1, p2, p3, p4, p5, p6, p7, p8, p9, p10,
p11, p12, p13, p14, p15, p16, p17, p18, p19, p20,
p21, p22, p23, p24, p25, p26, p27
]
from Draft import makeWire
line = makeWire(
points,
closed=False,
face=False,
support=None
)
fcc_print("Checking FEM new analysis...")
analysis = ObjectsFem.makeAnalysis(
self.active_doc,
"Analysis"
)
self.assertTrue(analysis, "FemTest of new analysis failed")
fcc_print("Checking FEM new solver...")
solver_object = ObjectsFem.makeSolverCalculixCcxTools(
self.active_doc,
"CalculiX"
)
solver_object.AnalysisType = "thermomech"
solver_object.GeometricalNonlinearity = "linear"
solver_object.ThermoMechSteadyState = True
solver_object.MatrixSolverType = "default"
solver_object.IterationsThermoMechMaximum = 2000
solver_object.IterationsControlParameterTimeUse = False
self.assertTrue(
solver_object,
"FemTest of new solver failed"
)
analysis.addObject(solver_object)
fcc_print("Checking FEM new material...")
material_object = ObjectsFem.makeMaterialFluid(
self.active_doc,
"FluidMaterial"
)
mat = material_object.Material
mat["Name"] = "Water"
mat["Density"] = "998 kg/m^3"
mat["SpecificHeat"] = "4.182 J/kg/K"
mat["DynamicViscosity"] = "1.003e-3 kg/m/s"
mat["VolumetricThermalExpansionCoefficient"] = "2.07e-4 m/m/K"
mat["ThermalConductivity"] = "0.591 W/m/K"
material_object.Material = mat
self.assertTrue(
material_object,
"FemTest of new material failed"
)
analysis.addObject(material_object)
fcc_print("Checking FEM Flow1D inlet constraint...")
Flow1d_inlet = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_inlet.SectionType = "Liquid"
Flow1d_inlet.LiquidSectionType = "PIPE INLET"
Flow1d_inlet.InletPressure = 0.1
Flow1d_inlet.References = [(line, "Edge1")]
self.assertTrue(
Flow1d_inlet,
"FemTest of new Flow1D inlet constraint failed"
)
analysis.addObject(Flow1d_inlet)
fcc_print("Checking FEM new Flow1D entrance constraint...")
Flow1d_entrance = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_entrance.SectionType = "Liquid"
Flow1d_entrance.LiquidSectionType = "PIPE ENTRANCE"
Flow1d_entrance.EntrancePipeArea = 31416.00
Flow1d_entrance.EntranceArea = 25133.00
Flow1d_entrance.References = [(line, "Edge2")]
self.assertTrue(
Flow1d_entrance, "FemTest of new Flow1D entrance constraint failed"
)
analysis.addObject(Flow1d_entrance)
fcc_print("Checking FEM new Flow1D manning constraint...")
Flow1d_manning = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_manning.SectionType = "Liquid"
Flow1d_manning.LiquidSectionType = "PIPE MANNING"
Flow1d_manning.ManningArea = 31416
Flow1d_manning.ManningRadius = 50
Flow1d_manning.ManningCoefficient = 0.002
Flow1d_manning.References = [(line, "Edge3"), (line, "Edge5")]
self.assertTrue(
Flow1d_manning,
"FemTest of new Flow1D manning constraint failed"
)
analysis.addObject(Flow1d_manning)
fcc_print("Checking FEM new Flow1D bend constraint...")
Flow1d_bend = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_bend.SectionType = "Liquid"
Flow1d_bend.LiquidSectionType = "PIPE BEND"
Flow1d_bend.BendPipeArea = 31416
Flow1d_bend.BendRadiusDiameter = 1.5
Flow1d_bend.BendAngle = 45
Flow1d_bend.BendLossCoefficient = 0.4
Flow1d_bend.References = [(line, "Edge4")]
self.assertTrue(
Flow1d_bend,
"FemTest of new Flow1D bend constraint failed"
)
analysis.addObject(Flow1d_bend)
fcc_print("Checking FEM new Flow1D enlargement constraint...")
Flow1d_enlargement = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_enlargement.SectionType = "Liquid"
Flow1d_enlargement.LiquidSectionType = "PIPE ENLARGEMENT"
Flow1d_enlargement.EnlargeArea1 = 31416.00
Flow1d_enlargement.EnlargeArea2 = 70686.00
Flow1d_enlargement.References = [(line, "Edge6")]
self.assertTrue(
Flow1d_enlargement,
"FemTest of new Flow1D enlargement constraint failed"
)
analysis.addObject(Flow1d_enlargement)
fcc_print("Checking FEM new Flow1D manning constraint...")
Flow1d_manning1 = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_manning1.SectionType = "Liquid"
Flow1d_manning1.LiquidSectionType = "PIPE MANNING"
Flow1d_manning1.ManningArea = 70686.00
Flow1d_manning1.ManningRadius = 75
Flow1d_manning1.ManningCoefficient = 0.002
Flow1d_manning1.References = [(line, "Edge7")]
self.assertTrue(
Flow1d_manning1,
"FemTest of new Flow1D manning constraint failed"
)
analysis.addObject(Flow1d_manning1)
fcc_print("Checking FEM new Flow1D contraction constraint...")
Flow1d_contraction = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_contraction.SectionType = "Liquid"
Flow1d_contraction.LiquidSectionType = "PIPE CONTRACTION"
Flow1d_contraction.ContractArea1 = 70686
Flow1d_contraction.ContractArea2 = 17671
Flow1d_contraction.References = [(line, "Edge8")]
self.assertTrue(
Flow1d_contraction,
"FemTest of new Flow1D contraction constraint failed"
)
analysis.addObject(Flow1d_contraction)
fcc_print("Checking FEM new Flow1D manning constraint...")
Flow1d_manning2 = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_manning2.SectionType = "Liquid"
Flow1d_manning2.LiquidSectionType = "PIPE MANNING"
Flow1d_manning2.ManningArea = 17671.00
Flow1d_manning2.ManningRadius = 37.5
Flow1d_manning2.ManningCoefficient = 0.002
Flow1d_manning2.References = [(line, "Edge11"), (line, "Edge9")]
self.assertTrue(
Flow1d_manning2,
"FemTest of new Flow1D manning constraint failed"
)
analysis.addObject(Flow1d_manning2)
fcc_print("Checking FEM new Flow1D gate valve constraint...")
Flow1d_gate_valve = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_gate_valve.SectionType = "Liquid"
Flow1d_gate_valve.LiquidSectionType = "PIPE GATE VALVE"
Flow1d_gate_valve.GateValvePipeArea = 17671
Flow1d_gate_valve.GateValveClosingCoeff = 0.5
Flow1d_gate_valve.References = [(line, "Edge10")]
self.assertTrue(
Flow1d_gate_valve,
"FemTest of new Flow1D gate valve constraint failed"
)
analysis.addObject(Flow1d_gate_valve)
fcc_print("Checking FEM new Flow1D enlargement constraint...")
Flow1d_enlargement1 = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_enlargement1.SectionType = "Liquid"
Flow1d_enlargement1.LiquidSectionType = "PIPE ENLARGEMENT"
Flow1d_enlargement1.EnlargeArea1 = 17671
Flow1d_enlargement1.EnlargeArea2 = 1e12
Flow1d_enlargement1.References = [(line, "Edge12")]
self.assertTrue(
Flow1d_enlargement1,
"FemTest of new Flow1D enlargement constraint failed"
)
analysis.addObject(Flow1d_enlargement1)
fcc_print("Checking FEM Flow1D outlet constraint...")
Flow1d_outlet = ObjectsFem.makeElementFluid1D(
self.active_doc,
"ElementFluid1D"
)
Flow1d_outlet.SectionType = "Liquid"
Flow1d_outlet.LiquidSectionType = "PIPE OUTLET"
Flow1d_outlet.OutletPressure = 0.1
Flow1d_outlet.References = [(line, "Edge13")]
self.assertTrue(
Flow1d_outlet,
"FemTest of new Flow1D inlet constraint failed"
)
analysis.addObject(Flow1d_outlet)
fcc_print("Checking FEM self weight constraint...")
Flow1d_self_weight = ObjectsFem.makeConstraintSelfWeight(
self.active_doc,
"ConstraintSelfWeight"
)
Flow1d_self_weight.Gravity_x = 0.0
Flow1d_self_weight.Gravity_y = 0.0
Flow1d_self_weight.Gravity_z = -1.0
self.assertTrue(
Flow1d_outlet,
"FemTest of new Flow1D self weight constraint failed"
)
analysis.addObject(Flow1d_self_weight)
fcc_print("Checking FEM new mesh...")
from ..data.ccx.Flow1D_mesh import create_nodes_Flow1D
from ..data.ccx.Flow1D_mesh import create_elements_Flow1D
mesh = Fem.FemMesh()
ret = create_nodes_Flow1D(mesh)
self.assertTrue(ret, "Import of mesh nodes failed")
ret = create_elements_Flow1D(mesh)
self.assertTrue(ret, "Import of mesh volumes failed")
mesh_object = self.active_doc.addObject(
"Fem::FemMeshObject",
self.mesh_name
)
mesh_object.FemMesh = mesh
self.assertTrue(mesh, "FemTest of new mesh failed")
analysis.addObject(mesh_object)
self.active_doc.recompute()
# set up the thermomech flow1d example
from femexamples.thermomech_flow1d import setup as flow1d
flow1d(self.active_doc, "ccxtools")
analysis = self.active_doc.Analysis
Flow1D_thermomech_analysis_dir = testtools.get_unit_test_tmp_dir(
self.temp_dir,