diff --git a/src/Mod/Fem/CMakeLists.txt b/src/Mod/Fem/CMakeLists.txt index b1049e25ec..bad64aa37f 100755 --- a/src/Mod/Fem/CMakeLists.txt +++ b/src/Mod/Fem/CMakeLists.txt @@ -51,6 +51,9 @@ SET(FemExamples_SRCS femexamples/ccx_cantilever_base_edge.py femexamples/ccx_cantilever_base_face.py femexamples/ccx_cantilever_base_solid.py + femexamples/ccx_cantilever_beam_circle.py + femexamples/ccx_cantilever_beam_pipe.py + femexamples/ccx_cantilever_beam_rect.py femexamples/ccx_cantilever_ele_hexa20.py femexamples/ccx_cantilever_ele_quad4.py femexamples/ccx_cantilever_ele_quad8.py diff --git a/src/Mod/Fem/femexamples/ccx_cantilever_base_edge.py b/src/Mod/Fem/femexamples/ccx_cantilever_base_edge.py index f1bb01361c..0731a1dbba 100644 --- a/src/Mod/Fem/femexamples/ccx_cantilever_base_edge.py +++ b/src/Mod/Fem/femexamples/ccx_cantilever_base_edge.py @@ -59,6 +59,7 @@ def setup_cantilever_base_edge(doc=None, solvertype="ccxtools"): doc.recompute() if FreeCAD.GuiUp: + load_line.ViewObject.Visibility = False geom_obj.ViewObject.Document.activeView().viewAxonometric() geom_obj.ViewObject.Document.activeView().fitAll() diff --git a/src/Mod/Fem/femexamples/ccx_cantilever_beam_circle.py b/src/Mod/Fem/femexamples/ccx_cantilever_beam_circle.py new file mode 100644 index 0000000000..f4e26bdc68 --- /dev/null +++ b/src/Mod/Fem/femexamples/ccx_cantilever_beam_circle.py @@ -0,0 +1,110 @@ +# *************************************************************************** +# * Copyright (c) 2021 Bernd Hahnebach * +# * * +# * 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 * +# * * +# *************************************************************************** + +import FreeCAD + +from . import manager +from .ccx_cantilever_base_edge import setup_cantilever_base_edge +from .manager import init_doc + + +def get_information(): + return { + "name": "CCX cantilever beam circle", + "meshtype": "edge", + "meshelement": "Seg3", + "constraints": ["fixed", "force"], + "solvers": ["calculix"], + "material": "solid", + "equation": "mechanical" + } + + +def get_explanation(header=""): + return header + """ + +To run the example from Python console use: +from femexamples.ccx_cantilever_beam_circle import setup +setup() + + +See forum topic post: +https://forum.freecadweb.org/viewtopic.php?f=18&t=16044 + +CalculiX cantilever: +- modeled with seg3 beam elements + +CrossSection: +- circle +- diameter 1000 mm + +max deflection: +doc = App.ActiveDocument +len = doc.CantileverLine.Shape.Length +iyy = doc.CrossSectionFace.Shape.MatrixOfInertia.A22 # Iyy +force = doc.ConstraintForce.Force +from FreeCAD import Units +ym = Units.Quantity(doc.MechanicalMaterial.Material["YoungsModulus"]).getValueAs("MPa") +w = force * len**3 / (3 * ym * iyy) +w # should print 149.0 mm + +CalculiX FEM max deflection: +- 146.7 mm +- Delta ca. 1.5 % + +""" + + +def setup(doc=None, solvertype="ccxtools"): + + # init FreeCAD document + if doc is None: + doc = init_doc() + + # explanation object + # just keep the following line and change text string in get_explanation method + manager.add_explanation_obj(doc, get_explanation(manager.get_header(get_information()))) + + diameter = 1000 + cs_wire = doc.addObject("Part::Circle", "WireOfCrossSection") + cs_wire.Radius = diameter / 2 + cs_wire.Placement = FreeCAD.Placement( + FreeCAD.Vector(0, 500, 500), + FreeCAD.Rotation(0, 90, 0), + FreeCAD.Vector(0, 0, 0), + ) + doc.recompute() + cs_face = doc.addObject("Part::Face", "CrossSectionFace") + cs_face.Sources = cs_wire + doc.recompute() + + # setup CalculiX cantilever + doc = setup_cantilever_base_edge(doc, solvertype) + beamsection_obj = doc.getObject("BeamCrossSection") + + # change cross section to a circular section + beamsection_obj.SectionType = "Circular" + beamsection_obj.CircDiameter = diameter + + doc.recompute() + return doc diff --git a/src/Mod/Fem/femexamples/ccx_cantilever_beam_pipe.py b/src/Mod/Fem/femexamples/ccx_cantilever_beam_pipe.py new file mode 100644 index 0000000000..0fb9b85e78 --- /dev/null +++ b/src/Mod/Fem/femexamples/ccx_cantilever_beam_pipe.py @@ -0,0 +1,119 @@ +# *************************************************************************** +# * Copyright (c) 2021 Bernd Hahnebach * +# * * +# * 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 * +# * * +# *************************************************************************** + +import FreeCAD + +from . import manager +from .ccx_cantilever_base_edge import setup_cantilever_base_edge +from .manager import init_doc + + +def get_information(): + return { + "name": "CCX cantilever beam pipe", + "meshtype": "edge", + "meshelement": "Seg3", + "constraints": ["fixed", "force"], + "solvers": ["calculix"], + "material": "solid", + "equation": "mechanical" + } + + +def get_explanation(header=""): + return header + """ + +To run the example from Python console use: +from femexamples.ccx_cantilever_beam_pipe import setup +setup() + + +See forum topic post: +https://forum.freecadweb.org/viewtopic.php?f=18&t=16044 + +CalculiX cantilever: +- modeled with seg3 beam elements + +CrossSection: +- pipe +- outer diameter 1000 mm +- thickness 100 mm +- inner diameter = 1000 - 2x100 = 800 +- inner radius = 1000/2 - 100 = 400 + +max deflection: +doc = App.ActiveDocument +len = doc.CantileverLine.Shape.Length +iyy = doc.CrossSectionFace.Shape.MatrixOfInertia.A22 # Iyy +force = doc.ConstraintForce.Force +from FreeCAD import Units +ym = Units.Quantity(doc.MechanicalMaterial.Material["YoungsModulus"]).getValueAs("MPa") +w = force * len**3 / (3 * ym * iyy) +w # should print 252.4 mm + +CalculiX FEM max deflection: +- edit the solver input: element type B32R instead B32R +- 249.8 mm +- Delta ca. 1.0 % + +""" + + +def setup(doc=None, solvertype="ccxtools"): + + # init FreeCAD document + if doc is None: + doc = init_doc() + + # explanation object + # just keep the following line and change text string in get_explanation method + manager.add_explanation_obj(doc, get_explanation(manager.get_header(get_information()))) + + diameter = 1000 + thickness = 100 + cs_wire_outer = doc.addObject("Part::Circle", "OuterWireOfCrossSection") + cs_wire_outer.Radius = diameter / 2 + cs_wire_outer.Placement = FreeCAD.Placement( + FreeCAD.Vector(0, 500, 500), + FreeCAD.Rotation(0, 90, 0), + FreeCAD.Vector(0, 0, 0), + ) + cs_wire_inner = doc.addObject("Part::Circle", "InnerWireOfCrossSection") + cs_wire_inner.Radius = (diameter / 2) - thickness + cs_wire_inner.Placement = cs_wire_outer.Placement + doc.recompute() + cs_face = doc.addObject("Part::Face", "CrossSectionFace") + cs_face.Sources = [cs_wire_outer, cs_wire_inner] + doc.recompute() + + # setup CalculiX cantilever + doc = setup_cantilever_base_edge(doc, solvertype) + beamsection_obj = doc.getObject("BeamCrossSection") + + # change cross section to a pipe section + beamsection_obj.SectionType = "Pipe" + beamsection_obj.PipeDiameter = diameter + beamsection_obj.PipeThickness = thickness + + doc.recompute() + return doc diff --git a/src/Mod/Fem/femexamples/ccx_cantilever_beam_rect.py b/src/Mod/Fem/femexamples/ccx_cantilever_beam_rect.py new file mode 100644 index 0000000000..6ac31dd33a --- /dev/null +++ b/src/Mod/Fem/femexamples/ccx_cantilever_beam_rect.py @@ -0,0 +1,111 @@ +# *************************************************************************** +# * Copyright (c) 2021 Bernd Hahnebach * +# * * +# * 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 * +# * * +# *************************************************************************** + +import FreeCAD + +from . import manager +from .ccx_cantilever_base_edge import setup_cantilever_base_edge +from .manager import init_doc + + +def get_information(): + return { + "name": "CCX cantilever beam rectangle", + "meshtype": "edge", + "meshelement": "Seg3", + "constraints": ["fixed", "force"], + "solvers": ["calculix"], + "material": "solid", + "equation": "mechanical" + } + + +def get_explanation(header=""): + return header + """ + +To run the example from Python console use: +from femexamples.ccx_cantilever_beam_rect import setup +setup() + + +See forum topic post: +https://forum.freecadweb.org/viewtopic.php?f=18&t=16044 + +CalculiX cantilever: +- modeled with seg3 beam elements + +CrossSection: +- rectangle +- width = 400 mm +- height = 1250 mm + +max deflection: +doc = App.ActiveDocument +len = doc.CantileverLine.Shape.Length +iyy = doc.CrossSectionFace.Shape.MatrixOfInertia.A22 # Iyy +force = doc.ConstraintForce.Force +from FreeCAD import Units +ym = Units.Quantity(doc.MechanicalMaterial.Material["YoungsModulus"]).getValueAs("MPa") +w = force * len**3 / (3 * ym * iyy) +w # should print 112.3 mm + +CalculiX FEM max deflection: +- 112.2 mm +- but the rotation seams 90 degree rotated (FIXME) + +""" + + +def setup(doc=None, solvertype="ccxtools"): + + # init FreeCAD document + if doc is None: + doc = init_doc() + + # explanation object + # just keep the following line and change text string in get_explanation method + manager.add_explanation_obj(doc, get_explanation(manager.get_header(get_information()))) + + width = 400 + height = 1250 + cs_face = doc.addObject("Part::Plane", "CrossSectionFace") + cs_face.Width = width + cs_face.Length = height + cs_face.Placement = FreeCAD.Placement( + FreeCAD.Vector(0, 500 - 0.5 * width, 500 + 0.5 * height), + FreeCAD.Rotation(0, 90, 0), + FreeCAD.Vector(0, 0, 0), + ) + doc.recompute() + + # setup CalculiX cantilever + doc = setup_cantilever_base_edge(doc, solvertype) + beamsection_obj = doc.getObject("BeamCrossSection") + + # change cross section to a circular section + beamsection_obj.SectionType = "Rectangular" + beamsection_obj.RectWidth = width + beamsection_obj.RectHeight = height + + doc.recompute() + return doc diff --git a/src/Mod/Fem/femmesh/meshsetsgetter.py b/src/Mod/Fem/femmesh/meshsetsgetter.py index 7dd017e71c..7bc5ab334a 100644 --- a/src/Mod/Fem/femmesh/meshsetsgetter.py +++ b/src/Mod/Fem/femmesh/meshsetsgetter.py @@ -729,7 +729,7 @@ class MeshSetsGetter(): # "beamsection_obj" : "beamsection_obj" if exists # "fluidsection_obj" : "fluidsection_obj" if exists # "shellthickness_obj" : shellthickness_obj" if exists - # "beam_normal" : normal vector for beams only + # "beam_axis_m" : main local beam axis for beams only # }, # {}, ... , {} ] @@ -756,8 +756,8 @@ class MeshSetsGetter(): matgeoset["mat_obj_name"] = mat_obj.Name matgeoset["ccx_mat_name"] = mat_obj.Material["Name"] matgeoset["beamsection_obj"] = beamsec_obj - # normal for this direction - matgeoset["beam_normal"] = beamdirection["normal"] + # beam_axis_m for this direction + matgeoset["beam_axis_m"] = beamdirection["beam_axis_m"] self.mat_geo_sets.append(matgeoset) def get_mat_geo_sets_single_mat_multiple_beam(self): @@ -783,8 +783,8 @@ class MeshSetsGetter(): matgeoset["mat_obj_name"] = mat_obj.Name matgeoset["ccx_mat_name"] = mat_obj.Material["Name"] matgeoset["beamsection_obj"] = beamsec_obj - # normal for this direction - matgeoset["beam_normal"] = beamdirection["normal"] + # beam_axis_m for this direction + matgeoset["beam_axis_m"] = beamdirection["beam_axis_m"] self.mat_geo_sets.append(matgeoset) def get_mat_geo_sets_multiple_mat_single_beam(self): @@ -809,8 +809,8 @@ class MeshSetsGetter(): matgeoset["mat_obj_name"] = mat_obj.Name matgeoset["ccx_mat_name"] = mat_obj.Material["Name"] matgeoset["beamsection_obj"] = beamsec_obj - # normal for this direction - matgeoset["beam_normal"] = beamdirection["normal"] + # beam_axis_m for this direction + matgeoset["beam_axis_m"] = beamdirection["beam_axis_m"] self.mat_geo_sets.append(matgeoset) def get_mat_geo_sets_multiple_mat_multiple_beam(self): @@ -840,8 +840,8 @@ class MeshSetsGetter(): matgeoset["mat_obj_name"] = mat_obj.Name matgeoset["ccx_mat_name"] = mat_obj.Material["Name"] matgeoset["beamsection_obj"] = beamsec_obj - # normal for this direction - matgeoset["beam_normal"] = beamdirection["normal"] + # beam_axis_m for this direction + matgeoset["beam_axis_m"] = beamdirection["beam_axis_m"] self.mat_geo_sets.append(matgeoset) # fluid diff --git a/src/Mod/Fem/femmesh/meshtools.py b/src/Mod/Fem/femmesh/meshtools.py index ee0a614da5..ac4a7aca22 100644 --- a/src/Mod/Fem/femmesh/meshtools.py +++ b/src/Mod/Fem/femmesh/meshtools.py @@ -537,14 +537,14 @@ def get_femelement_direction1D_set( theshape=None ): """ - get for each geometry edge direction, the normal and the element ids and + get for each geometry edge direction, the local direction m and the element ids and # write all into the beamrotation_objects means no return value, we're going to write into the beamrotation_objects dictionary FEMRotations1D is a list of dictionaries for every beamdirection of all edges beamrot_obj["FEMRotations1D"] = [{ "ids" : [theids], "direction" : direction, - "normal" : normal + "beam_axis_m" : beam_axis_m }, ... ] """ if len(beamrotation_objects) == 0: @@ -555,7 +555,7 @@ def get_femelement_direction1D_set( # add normals for each direction rotation_angle = 0 for rot in rotations_ids: - rot["normal"] = get_beam_normal(rot["direction"], rotation_angle) + rot["beam_axis_m"] = get_beam_main_axis_m(rot["direction"], rotation_angle) # key "Object" will be empty beamrotation_objects.append({"FEMRotations1D": rotations_ids, "ShortName": "Rstd"}) elif len(beamrotation_objects) == 1: @@ -567,7 +567,7 @@ def get_femelement_direction1D_set( # add normals for each direction rotation_angle = beamrotation_objects[0]["Object"].Rotation for rot in rotations_ids: - rot["normal"] = get_beam_normal(rot["direction"], rotation_angle) + rot["beam_axis_m"] = get_beam_main_axis_m(rot["direction"], rotation_angle) beamrotation_objects[0]["FEMRotations1D"] = rotations_ids beamrotation_objects[0]["ShortName"] = "R0" elif len(beamrotation_objects) > 1: @@ -611,7 +611,30 @@ def get_femelement_directions_theshape(femmesh, femelement_table, theshape): # ************************************************************************************************ -def get_beam_normal(beam_direction, defined_angle): +def get_beam_main_axis_m(beam_direction, defined_angle): + + # former name was get_beam_normal + # see forum topic https://forum.freecadweb.org/viewtopic.php?f=18&t=24878 + # beam_direction ... FreeCAD vector + # defined_angle ... degree + # base for the rotation: + # a beam_direction = (1, 0, 0) and angle = 0, returns (-0, 1, 0) + # https://forum.freecadweb.org/viewtopic.php?f=18&t=24878&start=30#p195567 + # https://forum.freecadweb.org/viewtopic.php?f=13&t=59239&start=140#p521999 + # changing the angle, changes the normal accordingly, 360 would again return (0,1,0) + + # CalxuliX uses negative z axis as base, if nothing is given in input file + # see the standard direction of 1-direction in CalxuliX manual + + # here the local main axis is called beam_axis_m the minor axis will be beam_axis_n + + # eventually a better name might be get_beam_rotation + + # FIXME: since we fix the base ange we would get this information out of the mesh edge too + print("beam_axis_m is retrieved from the geometry but we could get if from mesh edge too") + # print("beam_direction: {}".format(beam_direction)) + # print("defined_angle: {}".format(defined_angle)) + import math vector_a = beam_direction angle_rad = (math.pi / 180) * defined_angle @@ -667,6 +690,7 @@ def get_beam_normal(beam_direction, defined_angle): # dummy usage of the axis dot to get flake8 quiet del dot_x, dot_y, dot_z, dot, dot_nt + # print("normal_n: {}".format(normal_n)) return normal_n diff --git a/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py b/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py index 75991bc31d..c6b23c36f3 100644 --- a/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py +++ b/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py @@ -39,12 +39,25 @@ def write_femelement_geometry(f, ccxwriter): if "beamsection_obj"in matgeoset: # beam mesh beamsec_obj = matgeoset["beamsection_obj"] - normal = matgeoset["beam_normal"] + beam_axis_m = matgeoset["beam_axis_m"] + # in CalxuliX called the 1direction + # see meshtools.get_beam_main_axis_m(beam_direction, defined_angle) + section_nor = "{:.13G}, {:.13G}, {:.13G}\n".format( + beam_axis_m[0], + beam_axis_m[1], + beam_axis_m[2] + ) + print(section_nor) if beamsec_obj.SectionType == "Rectangular": - height = beamsec_obj.RectHeight.getValueAs("mm").Value - width = beamsec_obj.RectWidth.getValueAs("mm").Value + # see meshtools.get_beam_main_axis_m(beam_direction, defined_angle) + # the method get_beam_main_axis_m() which calculates the beam_axis_m vector + # returns for a line in x direction and angle 0 degree + # the y-axis as local main direction (beam_axis_m vector) + # in users head and in beam section object this is the Width + len_beam_axis_m = beamsec_obj.RectWidth.getValueAs("mm").Value + len_beam_axis_n = beamsec_obj.RectHeight.getValueAs("mm").Value section_type = ", SECTION=RECT" - section_geo = "{:.13G},{:.13G}\n".format(height, width) + section_geo = "{:.13G},{:.13G}\n".format(len_beam_axis_m, len_beam_axis_n) section_def = "*BEAM SECTION, {}{}{}\n".format( elsetdef, material, @@ -69,17 +82,10 @@ def write_femelement_geometry(f, ccxwriter): material, section_type ) - # see forum topic for output formatting of rotation - # https://forum.freecadweb.org/viewtopic.php?f=18&t=46133&p=405142#p405142 - section_nor = "{:.13G}, {:.13G}, {:.13G}\n".format( - normal[0], - normal[1], - normal[2] - ) f.write(section_def) f.write(section_geo) f.write(section_nor) - elif "fluidsection_obj"in matgeoset: # fluid mesh + elif "fluidsection_obj" in matgeoset: # fluid mesh fluidsec_obj = matgeoset["fluidsection_obj"] if fluidsec_obj.SectionType == "Liquid": section_type = fluidsec_obj.LiquidSectionType @@ -101,7 +107,7 @@ def write_femelement_geometry(f, ccxwriter): """ f.write(section_def) f.write(section_geo) - elif "shellthickness_obj"in matgeoset: # shell mesh + elif "shellthickness_obj" in matgeoset: # shell mesh shellth_obj = matgeoset["shellthickness_obj"] section_def = "*SHELL SECTION, {}{}\n".format(elsetdef, material) thickness = shellth_obj.Thickness.getValueAs("mm").Value diff --git a/src/Mod/Fem/femsolver/calculix/write_step_output.py b/src/Mod/Fem/femsolver/calculix/write_step_output.py index cf5d8b5613..117cf15f44 100644 --- a/src/Mod/Fem/femsolver/calculix/write_step_output.py +++ b/src/Mod/Fem/femsolver/calculix/write_step_output.py @@ -58,6 +58,8 @@ def write_step_output(f, ccxwriter): # dat file # reaction forces: freecadweb.org/tracker/view.php?id=2934 + # some hint can be found in this topic: + # https://forum.freecadweb.org/viewtopic.php?f=18&t=20664&start=10#p520642 if ccxwriter.member.cons_fixed: f.write("** outputs --> dat file\n") # reaction forces for all Constraint fixed