From 6bdff4769e77bdaf07873dbd5452e2dc9106e051 Mon Sep 17 00:00:00 2001 From: lyphrowny <79705170+lyphrowny@users.noreply.github.com> Date: Mon, 18 Mar 2024 17:16:10 +0000 Subject: [PATCH] FEA: fix 1d beam cross section orientation (#12833) * FEA: fix 1d beam cross section orientation * FEM: revert 1,2-directions to FreeCAD way; fix normal direction for z<0 * FEM: change 1-direction to +y axis Previous commit changes the 1-direction to -y axis, however originally it was directed at +y axis. This commit makes the 1-direction face +y axis * FEM: update cantilever tests * FEM: relax math.isclose by adding abs tolerance --- src/Mod/Fem/femmesh/meshtools.py | 128 +++++++++--------- .../calculix/write_femelement_geometry.py | 12 +- .../calculix/ccx_cantilever_beam_circle.inp | 2 +- .../calculix/ccx_cantilever_beam_pipe.inp | 2 +- .../calculix/ccx_cantilever_beam_rect.inp | 2 +- .../data/calculix/ccx_cantilever_ele_seg2.inp | 2 +- .../data/calculix/ccx_cantilever_ele_seg3.inp | 2 +- 7 files changed, 79 insertions(+), 71 deletions(-) diff --git a/src/Mod/Fem/femmesh/meshtools.py b/src/Mod/Fem/femmesh/meshtools.py index da4d614763..35af3b8c93 100644 --- a/src/Mod/Fem/femmesh/meshtools.py +++ b/src/Mod/Fem/femmesh/meshtools.py @@ -573,7 +573,7 @@ def get_femelement_direction1D_set( # since ccx needs to split them in sets anyway we need to take care of this too rotations_ids = get_femelement_directions_theshape(femmesh, femelement_table, theshape) # add normals for each direction - rotation_angle = beamrotation_objects[0]["Object"].Rotation + rotation_angle = beamrotation_objects[0]["Object"].Rotation.getValueAs("deg").Value for rot in rotations_ids: rot["beam_axis_m"] = get_beam_main_axis_m(rot["direction"], rotation_angle) beamrotation_objects[0]["FEMRotations1D"] = rotations_ids @@ -619,87 +619,91 @@ def get_femelement_directions_theshape(femmesh, femelement_table, theshape): # ************************************************************************************************ -def get_beam_main_axis_m(beam_direction, defined_angle): +def get_beam_main_axis_m(beam_direction: FreeCAD.Vector, defined_angle: int) -> FreeCAD.Vector: # former name was get_beam_normal # see forum topic https://forum.freecad.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.freecad.org/viewtopic.php?f=18&t=24878&start=30#p195567 - # https://forum.freecad.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 + # a beam_direction = (1, 0, 0) and angle = 0, returns (0, -1, 0) + # changing the angle, changes the normal accordingly, 360 would again return (0, -1, 0) + # + # original thread discussing this (kept for history): + # https://forum.freecad.org/viewtopic.php?f=18&t=24878&start=30#p195567 + # https://forum.freecad.org/viewtopic.php?f=13&t=59239&start=140#p521999 + # thread with updated analysis, that fixes a found bug: + # https://forum.freecad.org/viewtopic.php?t=85951 + # + # CalculiX 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 + # FIXME: since we fix the base angle we would get this information out of the mesh edge too + # upd: angle is not fixed anymore (?) (except for the vertical beam) 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 - nx = abs(math.cos(angle_rad)) - ny = abs(math.sin(angle_rad)) - if nx < 0.0000001: - nx = 0 - if ny < 0.0000001: - ny = 0 - # vector_n = [nx, ny] # not used ATM - if vector_a[0] != 0: - temp_valx = -(vector_a[1] + vector_a[2]) / vector_a[0] - else: - temp_valx = 0 - if vector_a[1] != 0: - temp_valy = -(vector_a[0] + vector_a[2]) / vector_a[1] - else: - temp_valy = 0 - if vector_a[2] != 0: - temp_valz = -(vector_a[0] + vector_a[1]) / vector_a[2] - else: - temp_valz = 0 + def normalize(vec: FreeCAD.Vector) -> FreeCAD.Vector: + return vec / vec.Length - # Dot_product_check - dot_x = None - dot_y = None - dot_z = None - dot_nt = None - if vector_a[0] != 0 and vector_a[1] == 0 and vector_a[2] == 0: - normal_n = [temp_valx, nx, ny] - dot_x = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - elif vector_a[0] == 0 and vector_a[1] != 0 and vector_a[2] == 0: - normal_n = [nx, temp_valy, ny] - dot_y = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - elif vector_a[0] == 0 and vector_a[1] == 0 and vector_a[2] != 0: - normal_n = [nx, ny, temp_valz] - dot_z = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - elif vector_a[0] == 0 and vector_a[1] != 0 and vector_a[2] != 0: - normal_n = [nx, temp_valy, ny] - dot_y = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - elif vector_a[0] != 0 and vector_a[1] == 0 and vector_a[2] != 0: - normal_n = [nx, ny, temp_valz] - dot_z = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - else: - normal_n = [temp_valx, nx, ny] - dot_nt = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] + def get_normal(vec: FreeCAD.Vector) -> FreeCAD.Vector: + x, y, z = vec - dot = vector_a[0] * normal_n[0] + vector_a[1] * normal_n[1] + vector_a[2] * normal_n[2] - FreeCAD.Console.PrintLog("{}\n".format(dot)) - FreeCAD.Console.PrintLog("{}\n".format(normal_n)) + if z == 0: # in xy plane + n = (0, 0, -1) # default in CalculiX + elif z != 0 and x == y == 0: # vertical beam + n = (-1, 0, 0) # or (0, -1, 0) + else: + # . (x, y, z) + # / \ + # vec / \ n + # / \ + # (0,0,0) .-------. (x, y, ?) + # + # since `n` should be perpendicular to `vec` + # calculate the dot product to find `?` + n = (x, y, -(x**2 + y**2) / z) + # we want the normal to point downwards for `vec` with + # z < 0 as well + n = normalize((1 | -(z > 0)) * FreeCAD.Vector(n)) - # dummy usage of the axis dot to get flake8 quiet - del dot_x, dot_y, dot_z, dot, dot_nt + # sanity check + dot = vec.dot(n) + if not math.isclose(dot, 0, abs_tol=1e-13): + FreeCAD.Console.PrintError( + f"Wrong calculation of normal vector for {vec = }! {n = }, {dot = }\n" + "Please consider submitting an issue\n" + ) - # print("normal_n: {}".format(normal_n)) - return normal_n + return n + + def rotate_around_vector(axis: FreeCAD.Vector, angle: int) -> FreeCAD.Vector: + """ + * - is the axis vector, which is looking at us + the rotation is done counter-clockwise around this vector + + More info is available at: + https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula + """ + angle = math.radians(angle) + axis = normalize(axis) + n = get_normal(axis) + # the mentioned formula rotates both the perpendicular and not parts of the vector + # but since we're rotating a normal (which is already perpendicular), one can + # simplify the formula a bit + rot = n * math.cos(angle) + axis.cross(n) * math.sin(angle) + # replace all the close-to-zero elements with 0 itself + # I was getting values of 1e-17 order, hence the threshold is of 1e-15 order + return FreeCAD.Vector([coord * (abs(coord) > 1e-15) for coord in rot]) + + # `+90` here is because the normal is calculated towards -z axis + # and we use 1,2-directions, that are rotated by `+90` degrees + return rotate_around_vector(beam_direction, defined_angle + 90) # ************************************************************************************************ diff --git a/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py b/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py index bbf368185a..e2c0176317 100644 --- a/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py +++ b/src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py @@ -51,11 +51,15 @@ def write_femelement_geometry(f, ccxwriter): if beamsec_obj.SectionType == "Rectangular": # 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 + # unless rotated, this vector points towards +y axis + # doesn't follow 1,2-direction order of CalculiX + # ^ (n, 2-direction) + # | + # | + # .----> (m, 1-direction) + # len_beam_axis_n = beamsec_obj.RectHeight.getValueAs("mm").Value + len_beam_axis_m = beamsec_obj.RectWidth.getValueAs("mm").Value section_type = ", SECTION=RECT" section_geo = "{:.13G},{:.13G}\n".format(len_beam_axis_m, len_beam_axis_n) section_def = "*BEAM SECTION, {}{}{}\n".format( diff --git a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_circle.inp b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_circle.inp index 8a491a68be..31bcd23b49 100644 --- a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_circle.inp +++ b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_circle.inp @@ -58,7 +58,7 @@ Eedges ** Sections *BEAM SECTION, ELSET=M0B0RstdD0, MATERIAL=MechanicalMaterial, SECTION=CIRC 1000 --0, 1, 0 +0, 1, -0 *********************************************************** ** At least one step is needed to run an CalculiX analysis of FreeCAD diff --git a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_pipe.inp b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_pipe.inp index 3e6f83f32b..3b00eee774 100644 --- a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_pipe.inp +++ b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_pipe.inp @@ -58,7 +58,7 @@ Eedges ** Sections *BEAM SECTION, ELSET=M0B0RstdD0, MATERIAL=MechanicalMaterial, SECTION=PIPE 500,100 --0, 1, 0 +0, 1, -0 *********************************************************** ** At least one step is needed to run an CalculiX analysis of FreeCAD diff --git a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_rect.inp b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_rect.inp index c7232b1c9c..6d3eb29682 100644 --- a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_rect.inp +++ b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_rect.inp @@ -58,7 +58,7 @@ Eedges ** Sections *BEAM SECTION, ELSET=M0B0RstdD0, MATERIAL=MechanicalMaterial, SECTION=RECT 400,1250 --0, 1, 0 +0, 1, -0 *********************************************************** ** At least one step is needed to run an CalculiX analysis of FreeCAD diff --git a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg2.inp b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg2.inp index 48ef5ce230..73b099e0ef 100644 --- a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg2.inp +++ b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg2.inp @@ -200,7 +200,7 @@ Eedges ** Sections *BEAM SECTION, ELSET=M0B0RstdD0, MATERIAL=MechanicalMaterial, SECTION=RECT 1000,1000 --0, 1, 0 +0, 1, -0 *********************************************************** ** At least one step is needed to run an CalculiX analysis of FreeCAD diff --git a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg3.inp b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg3.inp index 20c0708829..49fb62c6b9 100644 --- a/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg3.inp +++ b/src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg3.inp @@ -58,7 +58,7 @@ Eedges ** Sections *BEAM SECTION, ELSET=M0B0RstdD0, MATERIAL=MechanicalMaterial, SECTION=RECT 1000,1000 --0, 1, 0 +0, 1, -0 *********************************************************** ** At least one step is needed to run an CalculiX analysis of FreeCAD