kindred/create
1
1
Fork 0
Code Issues 30 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity

FEA: fix 1d beam cross section orientation (#12833)

Browse Source 
* 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
This commit is contained in:
lyphrowny
2024-03-18 17:16:10 +00:00
committed by GitHub
parent 555d06e875
commit 6bdff4769e
7 changed files with 79 additions and 71 deletions
Download Patch File Download Diff File Expand all files Collapse all files

128
src/Mod/Fem/femmesh/meshtools.py
View File

@@ -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)
# ************************************************************************************************

12
src/Mod/Fem/femsolver/calculix/write_femelement_geometry.py
View File

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

2
src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_circle.inp
View File

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

2
src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_pipe.inp
View File

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

2
src/Mod/Fem/femtest/data/calculix/ccx_cantilever_beam_rect.inp
View File

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

2
src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg2.inp
View File

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

2
src/Mod/Fem/femtest/data/calculix/ccx_cantilever_ele_seg3.inp
View File

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