kindred/create
1
1
Fork 0
Code Issues 37 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity
Files
b5e2e8c007a32b1df4499b91bfb73939fce166a5
create/src/Mod/Fem/femresult/resulttools.py
FEA-eng f1995ef76c FEM: ResultsPurge - delete all results objects including pipelines and filters (#18328) 
* FEM: Update resulttools.py

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* FEM: Update resulttools.py

* FEM: Update commands.py

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* FEM: Update resulttools.py

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
2024-12-07 07:17:02 +01:00

957 lines
30 KiB
Python
Raw Blame History

# ***************************************************************************
# * Copyright (c) 2017 Bernd Hahnebach <bernd@bimstatik.org> *
# * *
# * 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 *
# * *
# ***************************************************************************
__title__ = "Fem Tools for results"
__author__ = "Bernd Hahnebach"
__url__ = "https://www.freecad.org"
## \addtogroup FEM
# @{
import numpy as np
from math import isnan
import FreeCAD
from femtools.femutils import is_of_type
def purge_results(analysis):
"""Removes all result objects and result meshes from an analysis group.
Parameters
----------
analysis : Fem::FemAnalysis
analysis group as a container for all objects needed for the analysis
"""
# if analysis type check is used, result mesh
# without result obj is created in the analysis
# we could run into trouble in one loop because
# we will delete objects and try to access them later
# result object
for m in analysis.Group:
if m.isDerivedFrom("Fem::FemResultObject"):
if m.Mesh and is_of_type(m.Mesh, "Fem::MeshResult"):
analysis.Document.removeObject(m.Mesh.Name)
analysis.Document.removeObject(m.Name)
analysis.Document.recompute()
# result mesh
for m in analysis.Group:
if is_of_type(m, "Fem::MeshResult"):
analysis.Document.removeObject(m.Name)
analysis.Document.recompute()
# text object
for m in analysis.Group:
if is_of_type(m, "App::TextDocument"):
analysis.Document.removeObject(m.Name)
analysis.Document.recompute()
# result pipeline and filter
for m in analysis.Group:
if is_of_type(m, "Fem::FemPostPipeline"):
# delete associated filters
for filter_obj in m.Filter:
analysis.Document.removeObject(filter_obj.Name)
# delete the pipeline itself
analysis.Document.removeObject(m.Name)
analysis.Document.recompute()
def reset_mesh_deformation(resultobj):
"""Resets result mesh deformation.
Parameters
----------
resultobj : Fem::ResultMechanical
FreeCAD FEM mechanical result object
"""
if FreeCAD.GuiUp:
if resultobj.Mesh:
resultobj.Mesh.ViewObject.applyDisplacement(0.0)
def reset_mesh_color(resultobj):
"""Resets result mesh color
Parameters
----------
resultobj : Fem::ResultMechanical
FreeCAD FEM mechanical result object
"""
if FreeCAD.GuiUp:
if resultobj.Mesh:
resultobj.Mesh.ViewObject.NodeColor = {}
resultobj.Mesh.ViewObject.ElementColor = {}
node_numbers = resultobj.Mesh.FemMesh.Nodes.keys()
zero_values = [0] * len(node_numbers)
resultobj.Mesh.ViewObject.setNodeColorByScalars(node_numbers, zero_values)
def show_displacement(resultobj, displacement_factor=0.0):
if FreeCAD.GuiUp:
if resultobj.Mesh.ViewObject.Visibility is False:
resultobj.Mesh.ViewObject.Visibility = True
resultobj.Mesh.ViewObject.setNodeDisplacementByVectors(
resultobj.NodeNumbers, resultobj.DisplacementVectors
)
resultobj.Mesh.ViewObject.applyDisplacement(displacement_factor)
def show_result(resultobj, result_type="Sabs", limit=None):
"""Sets mesh color using selected type of results.
Parameters
----------
resultobj : Fem::ResultMechanical
FreeCAD FEM mechanical result object
result_type : str, optional
default is Sabs
FreeCAD FEM mechanical result object
- U1, U2, U3 - deformation
- Uabs - absolute deformation
- Sabs - Von Mises stress
limit : float
limit cutoff value. All values over the limit are treated
as equal to the limit. Useful for filtering out hotspots.
"""
if result_type == "None":
reset_mesh_color(resultobj.Mesh)
return
if resultobj:
if result_type == "Sabs":
values = resultobj.vonMises
elif result_type == "Uabs":
values = resultobj.DisplacementLengths
# TODO: the result object does have more result types to show, implement them
else:
match = {"U1": 0, "U2": 1, "U3": 2}
d = zip(*resultobj.DisplacementVectors)
values = list(d[match[result_type]])
show_color_by_scalar_with_cutoff(resultobj, values, limit)
else:
FreeCAD.Console.PrintError("Error, No result object given.\n")
def show_color_by_scalar_with_cutoff(resultobj, values, limit=None):
"""Sets mesh color using list of values. Internally used by show_result function.
Parameters
----------
resultobj : Fem::ResultMechanical
FreeCAD FEM mechanical result object
values : list of floats
the values to be colored and cutoff
has to be the same length as resultobj.NodeNumbers
resultobj.NodeNumbers has to be present in the resultobj
limit : float
limit cutoff value. All values over the limit are treated
as equal to the limit. Useful for filtering out hotspots.
"""
if limit:
filtered_values = []
for v in values:
if v > limit:
filtered_values.append(limit)
else:
filtered_values.append(v)
else:
filtered_values = values
if FreeCAD.GuiUp:
if resultobj.Mesh.ViewObject.Visibility is False:
resultobj.Mesh.ViewObject.Visibility = True
resultobj.Mesh.ViewObject.setNodeColorByScalars(resultobj.NodeNumbers, filtered_values)
def get_stats(res_obj, result_type):
"""Returns minimum and maximum value for provided result type
Parameters
----------
resultobj : Fem::ResultMechanical
FreeCAD FEM mechanical result object
result_type : str
type of FEM result
allowed are: see dict keys in def get_all_stats()
None - always return (0.0, 0.0)
"""
match_table = get_all_stats(res_obj)
match_table["None"] = (0.0, 0.0)
stats = ()
if result_type in match_table:
stats = match_table[result_type]
return stats
# - U1, U2, U3 - deformation
# - Uabs - absolute deformation
# - Sabs - Von Mises stress
# - Prin1 - Principal stress 1
# - Prin2 - Principal stress 2
# - Prin3 - Principal stress 3
# - MaxSear - maximum shear stress
# - Peeq - Equivalent plastic strain
# - Temp - Temperature
# - MFlow - MassFlowRate
# - NPress - NetworkPressure
def get_all_stats(res_obj):
"""Returns all stats for provided result type.
- U1, U2, U3 - deformation
- Uabs - absolute deformation
- Sabs - Von Mises stress
- MaxPrin - Principal stress 1
- MidPrin - Principal stress 2
- MinPrin - Principal stress 3
- MaxShear - maximum shear stress
- Peeq - peeq strain
- Temp - Temperature
- MFlow - MassFlowRate
- NPress - NetworkPressure
for more information on result types and names
see in code file src/Mod/Fem/App/FemVTKTools.cpp
the methods _getFreeCADMechResultVectorProperties()
and _getFreeCADMechResultScalarProperties()
as well as forum topic
https://forum.freecad.org/viewtopic.php?f=18&t=33106&start=30#p277434
Parameters
----------
resultobj : Fem::ResultMechanical
FreeCAD FEM mechanical result object
"""
m = res_obj.Stats
stats_dict = {
"U1": (m[0], m[1]),
"U2": (m[2], m[3]),
"U3": (m[4], m[5]),
"Uabs": (m[6], m[7]),
"Sabs": (m[8], m[9]),
"MaxPrin": (m[10], m[11]),
"MidPrin": (m[12], m[13]),
"MinPrin": (m[14], m[15]),
"MaxShear": (m[16], m[17]),
"Peeq": (m[18], m[19]),
"Temp": (m[20], m[21]),
"MFlow": (m[22], m[23]),
"NPress": (m[24], m[25]),
}
return stats_dict
def fill_femresult_stats(res_obj):
"""Fills a FreeCAD FEM mechanical result object with stats data.
Parameters
----------
resultobj : Fem::ResultMechanical
FreeCAD FEM mechanical result object
"""
FreeCAD.Console.PrintLog("Calculate stats list for result obj: " + res_obj.Name + "\n")
# set stats values to 0, they may not exist in res_obj
x_min = y_min = z_min = x_max = y_max = z_max = 0
a_max = a_min = s_max = s_min = 0
p1_min = p1_max = p2_min = p2_max = p3_min = p3_max = 0
ms_min = ms_max = peeq_min = peeq_max = 0
temp_min = temp_max = 0
mflow_min = mflow_max = npress_min = npress_max = 0
if res_obj.DisplacementVectors:
x_max, y_max, z_max = map(max, zip(*res_obj.DisplacementVectors))
x_min, y_min, z_min = map(min, zip(*res_obj.DisplacementVectors))
if res_obj.DisplacementLengths:
a_min = min(res_obj.DisplacementLengths)
a_max = max(res_obj.DisplacementLengths)
if res_obj.vonMises:
s_min = min(res_obj.vonMises)
s_max = max(res_obj.vonMises)
if res_obj.PrincipalMax:
p1_min = min(res_obj.PrincipalMax)
p1_max = max(res_obj.PrincipalMax)
if res_obj.PrincipalMed:
p2_min = min(res_obj.PrincipalMed)
p2_max = max(res_obj.PrincipalMed)
if res_obj.PrincipalMin:
p3_min = min(res_obj.PrincipalMin)
p3_max = max(res_obj.PrincipalMin)
if res_obj.MaxShear:
ms_min = min(res_obj.MaxShear)
ms_max = max(res_obj.MaxShear)
if res_obj.Peeq:
peeq_min = min(res_obj.Peeq)
peeq_max = max(res_obj.Peeq)
if res_obj.Temperature:
temp_min = min(res_obj.Temperature)
temp_max = max(res_obj.Temperature)
if res_obj.MassFlowRate:
# DisplacementVectors is empty, no_of_values needs to be set
mflow_min = min(res_obj.MassFlowRate)
mflow_max = max(res_obj.MassFlowRate)
if res_obj.NetworkPressure:
npress_min = min(res_obj.NetworkPressure)
npress_max = max(res_obj.NetworkPressure)
res_obj.Stats = [
x_min,
x_max,
y_min,
y_max,
z_min,
z_max,
a_min,
a_max,
s_min,
s_max,
p1_min,
p1_max,
p2_min,
p2_max,
p3_min,
p3_max,
ms_min,
ms_max,
peeq_min,
peeq_max,
temp_min,
temp_max,
mflow_min,
mflow_max,
npress_min,
npress_max,
]
"""
stat_types = [
"U1",
"U2",
"U3",
"Uabs",
"Sabs",
"MaxPrin",
"MidPrin",
"MinPrin",
"MaxShear",
"Peeq",
"Temp",
"MFlow",
"NPress"
]
"""
# len(stat_types) == 13*3 == 39
# do not forget to adapt initialization of all Stats items in modules:
# - module femobjects/_FemResultMechanical.py
# do not forget to adapt the def get_stats in:
# - get_stats in module femresult/resulttools.py
# - module femtest/testccxtools.py
# TODO: all stats stuff should be reimplemented
# maybe a dictionary would be far more robust than a list
FreeCAD.Console.PrintLog("Stats list for result obj: " + res_obj.Name + " calculated\n")
return res_obj
def add_disp_apps(res_obj):
res_obj.DisplacementLengths = calculate_disp_abs(res_obj.DisplacementVectors)
FreeCAD.Console.PrintLog("Added DisplacementLengths.\n")
return res_obj
def add_von_mises(res_obj):
mstress = []
iterator = zip(
res_obj.NodeStressXX,
res_obj.NodeStressYY,
res_obj.NodeStressZZ,
res_obj.NodeStressXY,
res_obj.NodeStressXZ,
res_obj.NodeStressYZ,
)
for Sxx, Syy, Szz, Sxy, Sxz, Syz in iterator:
mstress.append(calculate_von_mises((Sxx, Syy, Szz, Sxy, Sxz, Syz)))
res_obj.vonMises = mstress
FreeCAD.Console.PrintLog("Added von Mises stress.\n")
return res_obj
def add_principal_stress_std(res_obj):
# saved into PrincipalMax, PrincipalMed, PrincipalMin
# TODO may be use only one container for principal stresses in result object
# https://forum.freecad.org/viewtopic.php?f=18&t=33106&p=416006#p416006
# but which one is better
prinstress1 = []
prinstress2 = []
prinstress3 = []
shearstress = []
iterator = zip(
res_obj.NodeStressXX,
res_obj.NodeStressYY,
res_obj.NodeStressZZ,
res_obj.NodeStressXY,
res_obj.NodeStressXZ,
res_obj.NodeStressYZ,
)
for Sxx, Syy, Szz, Sxy, Sxz, Syz in iterator:
prin1, prin2, prin3, shear = calculate_principal_stress_std((Sxx, Syy, Szz, Sxy, Sxz, Syz))
prinstress1.append(prin1)
prinstress2.append(prin2)
prinstress3.append(prin3)
shearstress.append(shear)
res_obj.PrincipalMax = prinstress1
res_obj.PrincipalMed = prinstress2
res_obj.PrincipalMin = prinstress3
res_obj.MaxShear = shearstress
FreeCAD.Console.PrintLog("Added standard principal stresses and max shear values.\n")
#
# Add critical strain ratio using the Stress Modified Critical Strain (SMCS) criterion
# Forum Discussion: https://forum.freecad.org/viewtopic.php?f=18&t=35893#p303392
# Background: https://www.vtt.fi/inf/julkaisut/muut/2017/VTT-R-01177-17.pdf
#
# critical strain ratio = peeq / critical_strain (>1.0 indicates ductile rupture)
# peeq = equivalent plastic strain
# critical strain = alpha * np.exp(-beta * T)
# alpha and beta are material parameters,
# where alpha can be related to unixial test data (user input) and
# beta is normally kept fixed at 1.5,
# unless available from extensive research experiments
# T = pressure / von Mises stress (stress triaxiality)
#
MatMechNon = FreeCAD.ActiveDocument.getObject("MaterialMechanicalNonlinear")
if MatMechNon:
stress_strain = MatMechNon.YieldPoints
if stress_strain:
i = -1
while stress_strain[i] == "":
i -= 1
critical_uniaxial_strain = float(stress_strain[i].split(",")[1])
# stress triaxiality T = 1/3 for uniaxial test
alpha = np.sqrt(np.e) * critical_uniaxial_strain
beta = 1.5
if res_obj.Peeq:
res_obj.CriticalStrainRatio = calculate_csr(
prinstress1, prinstress2, prinstress3, alpha, beta, res_obj
)
return res_obj
def calculate_csr(ps1, ps2, ps3, alpha, beta, res_obj):
"""Calculate critical strain ratio.
Forum Discussion: https://forum.freecad.org/viewtopic.php?f=18&t=35893#p303392
Background: https://www.vtt.fi/inf/julkaisut/muut/2017/VTT-R-01177-17.pdf
critical strain ratio = peeq / critical_strain (>1.0 indicates ductile rupture)
peeq = equivalent plastic strain
critical strain = alpha * np.exp(-beta * T)
alpha and beta are material parameters,
where alpha can be related to unixial test data (user input) and
beta is normally kept fixed at 1.5,
unless available from extensive research experiments
T = pressure / von Mises stress (stress triaxiality)
"""
csr = [] # critical strain ratio
nsr = len(ps1) # number of stress results
for i in range(nsr):
p = (ps1[i] + ps2[i] + ps3[i]) / 3.0 # pressure
svm = np.sqrt(
1.5 * (ps1[i] - p) ** 2 + 1.5 * (ps2[i] - p) ** 2 + 1.5 * (ps3[i] - p) ** 2
) # von Mises stress: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion
if svm != 0.0:
T = p / svm # stress triaxiality
else:
T = 0.0
critical_strain = alpha * np.exp(-beta * T) # critical strain
csr.append(abs(res_obj.Peeq[i]) / critical_strain) # critical strain ratio
return csr
def get_concrete_nodes(res_obj):
"""Determine concrete / non-concrete nodes."""
from femmesh.meshtools import get_femnodes_by_refshape
femmesh = res_obj.Mesh.FemMesh
nsr = femmesh.NodeCount # nsr number of stress results
# ic[iic]:
# ic = flag for material type; iic = node number
# ic = 0: NOT ASSIGNED
# ic = 1: ReinforcedMaterial
# ic = 2: NOT ReinforcedMaterial
ic = np.zeros(nsr)
for obj in res_obj.getParentGroup().Group:
if obj.isDerivedFrom("App::MaterialObjectPython") and is_of_type(
obj, "Fem::MaterialReinforced"
):
FreeCAD.Console.PrintMessage("ReinforcedMaterial\n")
if obj.References == []:
for iic in range(nsr):
if ic[iic] == 0:
ic[iic] = 1
else:
for ref in obj.References:
concrete_nodes = get_femnodes_by_refshape(femmesh, ref)
for cn in concrete_nodes:
ic[cn - 1] = 1
elif obj.isDerivedFrom("App::MaterialObjectPython") and is_of_type(
obj, "Fem::MaterialCommon"
):
FreeCAD.Console.PrintMessage("No ReinforcedMaterial\n")
if obj.References == []:
for iic in range(nsr):
if ic[iic] == 0:
ic[iic] = 2
else:
for ref in obj.References:
non_concrete_nodes = get_femnodes_by_refshape(femmesh, ref)
for ncn in non_concrete_nodes:
ic[ncn - 1] = 2
return ic
def add_principal_stress_reinforced(res_obj):
#
# determine concrete / non-concrete nodes
#
ic = get_concrete_nodes(res_obj)
#
# calculate principal and max Shear and fill them in res_obj
#
# saved into PS1Vector, PS2Vector, PS3Vector
# TODO may be use only one container for principal stresses in result object
# https://forum.freecad.org/viewtopic.php?f=18&t=33106&p=416006#p416006
# but which one is better
prinstress1 = []
prinstress2 = []
prinstress3 = []
shearstress = []
ps1v = []
ps2v = []
ps3v = []
#
# additional arrays to hold reinforcement ratios
# and mohr coulomb stress
#
rhx = []
rhy = []
rhz = []
moc = []
# material parameter
for obj in res_obj.getParentGroup().Group:
if is_of_type(obj, "Fem::MaterialReinforced"):
matrix_af = float(
FreeCAD.Units.Quantity(obj.Material["AngleOfFriction"]).getValueAs("rad")
)
matrix_cs = float(
FreeCAD.Units.Quantity(obj.Material["CompressiveStrength"]).getValueAs("MPa")
)
reinforce_yield = float(
FreeCAD.Units.Quantity(obj.Reinforcement["YieldStrength"]).getValueAs("MPa")
)
# print(matrix_af)
# print(matrix_cs)
# print(reinforce_yield)
iterator = zip(
res_obj.NodeStressXX,
res_obj.NodeStressYY,
res_obj.NodeStressZZ,
res_obj.NodeStressXY,
res_obj.NodeStressXZ,
res_obj.NodeStressYZ,
)
for isv, stress_tensor in enumerate(iterator):
rhox = 0.0
rhoy = 0.0
rhoz = 0.0
mc = 0.0
if ic[isv] == 1:
#
# for concrete scxx etc. are affected by
# reinforcement (see calculate_rho(stress_tensor)). for all other
# materials scxx etc. are the original stresses
#
rhox, rhoy, rhoz = calculate_rho(stress_tensor, reinforce_yield)
prin1, prin2, prin3, shear, psv = calculate_principal_stress_reinforced(stress_tensor)
prinstress1.append(prin1)
prinstress2.append(prin2)
prinstress3.append(prin3)
shearstress.append(shear)
ps1v.append(psv[0])
ps2v.append(psv[1])
ps3v.append(psv[2])
#
# reinforcement ratios and mohr coulomb criterion
#
rhx.append(rhox)
rhy.append(rhoy)
rhz.append(rhoz)
if ic[isv] == 1:
mc = calculate_mohr_coulomb(prin1, prin3, matrix_af, matrix_cs)
moc.append(mc)
res_obj.PrincipalMax = prinstress1
res_obj.PrincipalMed = prinstress2
res_obj.PrincipalMin = prinstress3
res_obj.MaxShear = shearstress
#
# additional concrete and principal stress plot
# results for use in _ViewProviderFemResultMechanical
#
res_obj.ReinforcementRatio_x = rhx
res_obj.ReinforcementRatio_y = rhy
res_obj.ReinforcementRatio_z = rhz
res_obj.MohrCoulomb = moc
res_obj.PS1Vector = ps1v
res_obj.PS2Vector = ps2v
res_obj.PS3Vector = ps3v
FreeCAD.Console.PrintLog(
"Added reinforcement principal stresses and max shear values as well as "
"reinforcment ratios, Mohr Coloumb values.\n"
)
return res_obj
def compact_result(res_obj):
"""
compacts result.Mesh and appropriate result.NodeNumbers
"""
# as workaround for https://www.freecad.org/tracker/view.php?id=2873
# get compact mesh data
from femmesh.meshtools import compact_mesh
compact_femmesh_data = compact_mesh(res_obj.Mesh.FemMesh)
compact_femmesh = compact_femmesh_data[0]
node_map = compact_femmesh_data[1]
# FreeCAD result obj does not support elem results ATM
# elem_map = compact_femmesh_data[2]
# set result mesh
res_obj.Mesh.FemMesh = compact_femmesh
# set result node numbers
new_node_numbers = []
for old_node_id in res_obj.NodeNumbers:
new_node_numbers.append(node_map[old_node_id])
res_obj.NodeNumbers = new_node_numbers
return res_obj
def calculate_von_mises(stress_tensor):
"""Calculate Von mises stress.
See http://en.wikipedia.org/wiki/Von_Mises_yield_criterion
Simplification: https://forum.freecad.org/viewtopic.php?f=18&t=33974&p=296542#p296542
stress_tensor ... (Sxx, Syy, Szz, Sxy, Sxz, Syz)
"""
normal = stress_tensor[:3]
shear = stress_tensor[3:]
pressure = np.average(normal)
von_mises = np.sqrt(
1.5 * np.linalg.norm(normal - pressure) ** 2 + 3.0 * np.linalg.norm(shear) ** 2
)
return von_mises
def calculate_principal_stress_std(stress_tensor):
# if NaN is inside the array, which can happen on Calculix frd result files return NaN
# https://forum.freecad.org/viewtopic.php?f=22&t=33911&start=10#p284229
# https://forum.freecad.org/viewtopic.php?f=18&t=32649#p274291
for s in stress_tensor:
if isnan(s) is True:
return (float("NaN"), float("NaN"), float("NaN"), float("NaN"))
s11 = stress_tensor[0] # Sxx
s22 = stress_tensor[1] # Syy
s33 = stress_tensor[2] # Szz
s12 = stress_tensor[3] # Sxy
s31 = stress_tensor[4] # Sxz
s23 = stress_tensor[5] # Syz
sigma = np.array(
[[s11, s12, s31], [s12, s22, s23], [s31, s23, s33]]
) # https://forum.freecad.org/viewtopic.php?f=18&t=24637&start=10#p240408
eigvals = list(np.linalg.eigvalsh(sigma))
eigvals.sort()
eigvals.reverse()
maxshear = (eigvals[0] - eigvals[2]) / 2.0
return (eigvals[0], eigvals[1], eigvals[2], maxshear)
def calculate_principal_stress_reinforced(stress_tensor):
"""Calculate principal stress vectors and values.
For total stresses use:
stress_tensor[0], stress_tensor[1], stress_tensor[2]
on the diagonal of the stress tensor
Difference with the original method:
https://forum.freecad.org/viewtopic.php?f=18&t=33106&start=90#p296539
"""
s11 = stress_tensor[0] # Sxx
s22 = stress_tensor[1] # Syy
s33 = stress_tensor[2] # Szz
s12 = stress_tensor[3] # Sxy
s31 = stress_tensor[4] # Sxz
s23 = stress_tensor[5] # Syz
sigma = np.array(
[[s11, s12, s31], [s12, s22, s23], [s31, s23, s33]]
) # https://forum.freecad.org/viewtopic.php?f=18&t=24637&start=10#p240408
eigenvalues, eigenvectors = np.linalg.eig(sigma)
#
# suppress complex eigenvalue and vectors that may occur for
# near-zero (numerical noise) stress fields
#
eigenvalues = eigenvalues.real
eigenvectors = eigenvectors.real
eigenvectors[:, 0] = eigenvalues[0] * eigenvectors[:, 0]
eigenvectors[:, 1] = eigenvalues[1] * eigenvectors[:, 1]
eigenvectors[:, 2] = eigenvalues[2] * eigenvectors[:, 2]
idx = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
maxshear = (eigenvalues[0] - eigenvalues[2]) / 2.0
return (
eigenvalues[0],
eigenvalues[1],
eigenvalues[2],
maxshear,
tuple([tuple(row) for row in eigenvectors.T]),
)
def calculate_rho(stress_tensor, fy):
"""Calculation of Reinforcement Ratios and Concrete Stresses
(in accordance with http://heronjournal.nl/53-4/3.pdf)
Parameters
----------
- fy: factored yield strength of reinforcement bars
See post:
https://forum.freecad.org/viewtopic.php?f=18&t=28821
"""
rmin = 1.0e9
eqmin = 14
sxx = stress_tensor[0]
syy = stress_tensor[1]
szz = stress_tensor[2]
sxy = stress_tensor[3]
syz = stress_tensor[5]
sxz = stress_tensor[4]
rhox = np.zeros(15)
rhoy = np.zeros(15)
rhoz = np.zeros(15)
# i1=sxx+syy+szz NOT USED
# i2=sxx*syy+syy*szz+szz*sxx-sxy**2-sxz**2-syz**2 NOT USED
i3 = sxx * syy * szz + 2 * sxy * sxz * syz - sxx * syz**2 - syy * sxz**2 - szz * sxy**2
# Solution (5)
d = sxx * syy - sxy**2
if d != 0.0:
rhoz[0] = i3 / d / fy
# Solution (6)
d = sxx * szz - sxz**2
if d != 0.0:
rhoy[1] = i3 / d / fy
# Solution (7)
d = syy * szz - syz**2
if d != 0.0:
rhox[2] = i3 / d / fy
# Solution (9)
if sxx != 0.0:
fc = sxz * sxy / sxx - syz
fxy = sxy**2 / sxx
fxz = sxz**2 / sxx
# Solution (9+)
rhoy[3] = syy - fxy + fc
rhoy[3] /= fy
rhoz[3] = szz - fxz + fc
rhoz[3] /= fy
# Solution (9-)
rhoy[4] = syy - fxy - fc
rhoy[4] /= fy
rhoz[4] = szz - fxz - fc
rhoz[4] /= fy
# Solution (10)
if syy != 0.0:
fc = syz * sxy / syy - sxz
fxy = sxy**2 / syy
fyz = syz**2 / syy
# Solution (10+)
rhox[5] = sxx - fxy + fc
rhox[5] /= fy
rhoz[5] = szz - fyz + fc
rhoz[5] /= fy
# Solution (10-)vm
rhox[6] = sxx - fxy - fc
rhox[6] /= fy
rhoz[6] = szz - fyz - fc
rhoz[6] /= fy
# Solution (11)
if szz != 0.0:
fc = sxz * syz / szz - sxy
fxz = sxz**2 / szz
fyz = syz**2 / szz
# Solution (11+)
rhox[7] = sxx - fxz + fc
rhox[7] /= fy
rhoy[7] = syy - fyz + fc
rhoy[7] /= fy
# Solution (11-)
rhox[8] = sxx - fxz - fc
rhox[8] /= fy
rhoy[8] = syy - fyz - fc
rhoy[8] /= fy
# Solution (13)
rhox[9] = (sxx + sxy + sxz) / fy
rhoy[9] = (syy + sxy + syz) / fy
rhoz[9] = (szz + sxz + syz) / fy
# Solution (14)
rhox[10] = (sxx + sxy - sxz) / fy
rhoy[10] = (syy + sxy - syz) / fy
rhoz[10] = (szz - sxz - syz) / fy
# Solution (15)
rhox[11] = (sxx - sxy - sxz) / fy
rhoy[11] = (syy - sxy + syz) / fy
rhoz[11] = (szz - sxz + syz) / fy
# Solution (16)
rhox[12] = (sxx - sxy + sxz) / fy
rhoy[12] = (syy - sxy - syz) / fy
rhoz[12] = (szz + sxz - syz) / fy
# Solution (17)
if syz != 0.0:
rhox[13] = (sxx - sxy * sxz / syz) / fy
if sxz != 0.0:
rhoy[13] = (syy - sxy * syz / sxz) / fy
if sxy != 0.0:
rhoz[13] = (szz - sxz * syz / sxy) / fy
for ir in range(0, rhox.size):
if rhox[ir] >= -1.0e-10 and rhoy[ir] >= -1.0e-10 and rhoz[ir] > -1.0e-10:
# Concrete Stresses
scxx = sxx - rhox[ir] * fy
scyy = syy - rhoy[ir] * fy
sczz = szz - rhoz[ir] * fy
ic1 = scxx + scyy + sczz
ic2 = scxx * scyy + scyy * sczz + sczz * scxx - sxy**2 - sxz**2 - syz**2
ic3 = (
scxx * scyy * sczz
+ 2 * sxy * sxz * syz
- scxx * syz**2
- scyy * sxz**2
- sczz * sxy**2
)
if ic1 <= 1.0e-6 and ic2 >= -1.0e-6 and ic3 <= 1.0e-6:
rsum = rhox[ir] + rhoy[ir] + rhoz[ir]
if rsum < rmin and rsum > 0.0:
rmin = rsum
eqmin = ir
return rhox[eqmin], rhoy[eqmin], rhoz[eqmin]
def calculate_mohr_coulomb(prin1, prin3, phi, fck):
"""Calculation of Mohr Coulomb yield criterion to judge
concrete crushing and shear failure.
Parameters
----------
- phi: angle of internal friction
- fck: factored compressive strength of the matrix material (usually concrete)
"""
coh = fck * (1 - np.sin(phi)) / 2 / np.cos(phi)
mc_stress = (prin1 - prin3) + (prin1 + prin3) * np.sin(phi) - 2.0 * coh * np.cos(phi)
if mc_stress < 0.0:
mc_stress = 0.0
return mc_stress
def calculate_disp_abs(displacements):
# see https://forum.freecad.org/viewtopic.php?f=18&t=33106&start=100#p296657
return [np.linalg.norm(nd) for nd in displacements]
## @}
Reference in New Issue View Git Blame Copy Permalink