Merge branch 'FreeCAD:master' into FixDefaultShapeVertexColor

This commit is contained in:
Benjamin Alterauge
2021-08-04 16:26:17 +02:00
committed by GitHub
9 changed files with 403 additions and 27 deletions

View File

@@ -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

View File

@@ -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()

View File

@@ -0,0 +1,110 @@
# ***************************************************************************
# * Copyright (c) 2021 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. *
# * *
# * 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

View File

@@ -0,0 +1,119 @@
# ***************************************************************************
# * Copyright (c) 2021 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. *
# * *
# * 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

View File

@@ -0,0 +1,111 @@
# ***************************************************************************
# * Copyright (c) 2021 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. *
# * *
# * 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

View File

@@ -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

View File

@@ -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

View File

@@ -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

View File

@@ -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