Skip to content

Feature/corotational shell #37

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pruliere wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Feature/corotational shell #37

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
5c037e6
Add corotational formulation for shell elements
pruliere Apr 30, 2026
c59814e
Improve localframe management.
pruliere May 3, 2026
8f51504
Correct the drill constraint and add option to derive a consistant me…
pruliere May 6, 2026
95fda0a
Add selective reduced integration elements for plate
pruliere May 11, 2026
7bafb54
make the navier stokes example work
pruliere May 11, 2026
342d1cb
A few updates
pruliere May 11, 2026
f13a458
add some plate formulations
pruliere May 12, 2026
21f98e2
typo
pruliere May 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
415 changes: 246 additions & 169 deletions examples/user_equation/navier_stokes.py
View file
Open in desktop

Large diffs are not rendered by default.

137 changes: 44 additions & 93 deletions fedoo/constitutivelaw/shell.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,32 +10,16 @@
class ShellBase(ConstitutiveLaw):
# base model class that should derive any other shell constitutive laws
def __init__(self, thickness, k=1, name=""):
# assert get_Dimension() == '3D', "No 2D model for a shell kinematic. Choose '3D' problem dimension."

ConstitutiveLaw.__init__(self, name) # heritage

self.__thickness = thickness
self.__k = k
self.__GeneralizedStrain = None
self.__GeneralizedStress = None

def GetThickness(self):
return self.__thickness

def Get_k(self):
return self.__k

def GetShellRigidityMatrix(self):
raise NameError('"GetShellRigidityMatrix" not implemented, contact developer.')

def GetShellRigidityMatrix_RI(self):
raise NameError(
'"GetShellRigidityMatrix_RI" not implemented, contact developer.'
)
self.thickness = thickness
"""Shell thickness."""
self.k = k
"""Shear shape factor."""

def GetShellRigidityMatrix_FI(self):
def get_shell_stiffness_matrix(self):
raise NameError(
'"GetShellRigidityMatrix_FI" not implemented, contact developer.'
'"get_shell_stiffness_matrix" not implemented, contact developer.'
)

def update(self, assembly, pb):
Expand All @@ -46,13 +30,13 @@ def update(self, assembly, pb):
assembly.sv["GeneralizedStrain"] = 0
assembly.sv["GeneralizedStress"] = 0
else:
GeneralizedStrainOp = assembly.weakform.GetGeneralizedStrainOperator()
GeneralizedStrainOp = assembly.weakform.generalized_strain_operator()
GeneralizedStrain = [
op if np.isscalar(op) else assembly.get_gp_results(op, U)
for op in GeneralizedStrainOp
]

H = self.GetShellRigidityMatrix()
H = self.get_shell_stiffness_matrix()

# all terms are computed. Perhaps could be optimized by computed only the termes related to the associated weak form (eg shear for selective integration)
assembly.sv["GeneralizedStress"] = [
Expand All @@ -67,11 +51,10 @@ def update(self, assembly, pb):
assembly.sv["GeneralizedStrain"] = GeneralizedStrain

def get_strain(self, assembly, **kargs):
"""
Return the last computed strain associated to the given assembly
"""Return the last computed strain associated to the given assembly.

Parameters
------------------------
----------
assembly: Assembly

position : float (optional)
Expand All @@ -87,7 +70,7 @@ def get_strain(self, assembly, **kargs):
StrainTensorList object containing the strain at integration point
"""
position = kargs.get("position", 1)
z = position * self.GetThickness() / 2
z = position * self.thickness / 2

Strain = StrainTensorList([0 for i in range(6)])
GeneralizedStrain = assembly.sv["GeneralizedStrain"]
Expand All @@ -107,69 +90,38 @@ def get_stress(self, **kargs):


class ShellHomogeneous(ShellBase):
def __init__(self, MatConstitutiveLaw, thickness, k=1, name=""):
# assert get_Dimension() == '3D', "No 2D model for a shell kinematic. Choose '3D' problem dimension."
def __init__(self, material, thickness, k=1, name=""):
# k: shear shape factor

if isinstance(MatConstitutiveLaw, str):
MatConstitutiveLaw = ConstitutiveLaw.get_all()[MatConstitutiveLaw]
if isinstance(material, str):
material = ConstitutiveLaw.get_all()[material]

ShellBase.__init__(self, thickness, k, name) # heritage

self.__material = MatConstitutiveLaw

def GetMaterial(self):
return self.__material
self.material = material

def GetShellRigidityMatrix(self):
Hplane = self.__material.get_elastic_matrix(
def get_shell_stiffness_matrix(self):
Hplane = self.material.get_elastic_matrix(
"2Dstress"
) # membrane rigidity matrix with plane stress assumption
Hplane = np.array(
[[Hplane[i][j] for j in [0, 1, 3]] for i in [0, 1, 3]], dtype="object"
)
Hshear = self.__material.get_elastic_matrix()
Hshear = self.material.get_elastic_matrix()
Hshear = np.array(
[[Hshear[i][j] for j in [4, 5]] for i in [4, 5]], dtype="object"
)

H = np.zeros((8, 8), dtype="object")
H[:3, :3] = self.GetThickness() * Hplane # Membrane
H[3:6, 3:6] = (
self.GetThickness() ** 3 / 12
) * Hplane # Flexual rigidity matrix
H[6:8, 6:8] = (self.Get_k() * self.GetThickness()) * Hshear

return H

def GetShellRigidityMatrix_RI(self):
# only shear component are given for reduce integration part
Hshear = self.__material.get_elastic_matrix()
Hshear = np.array(
[[Hshear[i][j] for j in [4, 5]] for i in [4, 5]], dtype="object"
)

return (self.Get_k() * self.GetThickness()) * Hshear

def GetShellRigidityMatrix_FI(self):
# membrane and flexural component are given for full integration part
Hplane = self.__material.get_elastic_matrix(
"2Dstress"
) # membrane rigidity matrix with plane stress assumption
Hplane = np.array(
[[Hplane[i][j] for j in [0, 1, 3]] for i in [0, 1, 3]], dtype="object"
)

H = np.zeros((6, 6), dtype="object")
H[:3, :3] = self.GetThickness() * Hplane # Membrane
H[3:6, 3:6] = (
self.GetThickness() ** 3 / 12
) * Hplane # Flexual rigidity matrix
H[:3, :3] = self.thickness * Hplane # Membrane
H[3:6, 3:6] = (self.thickness**3 / 12) * Hplane # Flexual rigidity matrix
H[6:8, 6:8] = (self.k * self.thickness) * Hshear

return H

def get_stress(self, assembly, **kargs):
Strain = self.get_strain(assembly, **kargs)
Hplane = self.__material.get_elastic_matrix(
Hplane = self.material.get_elastic_matrix(
"2Dstress"
) # membrane rigidity matrix with plane stress assumption
Stress = [
Expand All @@ -185,7 +137,7 @@ def get_stress(self, assembly, **kargs):
)
for i in range(4)
] # SXX, SYY, SXY (SZZ should be = 0)
Hshear = self.__material.get_elastic_matrix()
Hshear = self.material.get_elastic_matrix()
Stress += [
sum(
[
Expand All @@ -203,7 +155,7 @@ def get_stress(self, assembly, **kargs):
return StressTensorList(Stress)

def GetStressDistribution(self, assembly, pg, resolution=100):
h = self.GetThickness()
h = self.thickness
z = np.arange(-h / 2, h / 2, h / resolution)

Strain = StrainTensorList([0 for i in range(6)])
Expand All @@ -218,7 +170,7 @@ def GetStressDistribution(self, assembly, pg, resolution=100):
Strain[4] = GeneralizedStrain[6][pg] * np.ones_like(z) # 2epsXZ -> shear
Strain[5] = GeneralizedStrain[6][pg] * np.ones_like(z) # 2epsYZ -> shear

Hplane = self.__material.get_elastic_matrix(
Hplane = self.material.get_elastic_matrix(
"2Dstress"
) # membrane rigidity matrix with plane stress assumption
Stress = [
Expand All @@ -234,7 +186,7 @@ def GetStressDistribution(self, assembly, pg, resolution=100):
)
for i in range(4)
] # SXX, SYY, SXY (SZZ should be = 0)
Hshear = self.__material.get_elastic_matrix()
Hshear = self.material.get_elastic_matrix()
Stress += [
sum(
[
Expand All @@ -253,25 +205,25 @@ def GetStressDistribution(self, assembly, pg, resolution=100):


class ShellLaminate(ShellBase):
def __init__(self, listMat, listThickness, k=1, name=""):
def __init__(self, listMat, list_thickness, k=1, name=""):
# assert get_Dimension() == '3D', "No 2D model for a shell kinematic. Choose '3D' problem dimension."

self.__listMat = [
ConstitutiveLaw.get_all()[mat] if isinstance(mat, str) else mat
for mat in listMat
]
thickness = sum(listThickness) # total thickness
thickness = sum(list_thickness) # total thickness

self.__layer = (
np.hstack((0, np.cumsum(listThickness))) - np.sum(listThickness) / 2
np.hstack((0, np.cumsum(list_thickness))) - np.sum(list_thickness) / 2
) # z coord of layers interfaces
self.__listThickness = listThickness
self.list_thickness = list_thickness

ShellBase.__init__(self, thickness, k, name) # heritage

def GetShellRigidityMatrix(self):
def get_shell_stiffness_matrix(self):
H = np.zeros((8, 8), dtype="object")
for i in range(len(self.__listThickness)):
for i in range(len(self.list_thickness)):
Hplane = self.__listMat[i].get_elastic_matrix(
"2Dstress"
) # membrane rigidity matrix with plane stress assumption
Expand All @@ -283,7 +235,7 @@ def GetShellRigidityMatrix(self):
[[Hshear[i][j] for j in [4, 5]] for i in [4, 5]], dtype="object"
)

H[0:3, 0:3] += self.__listThickness[i] * Hplane # Membrane
H[0:3, 0:3] += self.list_thickness[i] * Hplane # Membrane
H[0:3, 3:6] += (
0.5 * (self.__layer[i + 1] ** 2 - self.__layer[i] ** 2) * Hplane
)
Expand All @@ -293,34 +245,34 @@ def GetShellRigidityMatrix(self):
H[3:6, 3:6] += (
(1 / 3) * (self.__layer[i + 1] ** 3 - self.__layer[i] ** 3) * Hplane
) # Flexual rigidity matrix
H[6:8, 6:8] += (self.Get_k() * self.__listThickness[i]) * Hshear
H[6:8, 6:8] += (self.k * self.list_thickness[i]) * Hshear

return H

def GetShellRigidityMatrix_RI(self):
def get_shell_stiffness_matrix_RI(self):
# only shear component are given for reduce integration part
H = np.zeros((2, 2), dtype="object")
for i in range(len(self.__listThickness)):
for i in range(len(self.list_thickness)):
Hshear = self.__listMat[i].get_elastic_matrix()
Hshear = np.array(
[[Hshear[i][j] for j in [4, 5]] for i in [4, 5]], dtype="object"
)
H += (self.Get_k() * self.__listThickness[i]) * Hshear
H += (self.k * self.list_thickness[i]) * Hshear

return H

def GetShellRigidityMatrix_FI(self):
def get_shell_stiffness_matrix_FI(self):
# membrane and flexural component are given for full integration part
H = np.zeros((6, 6), dtype="object")
for i in range(len(self.__listThickness)):
for i in range(len(self.list_thickness)):
Hplane = self.__listMat[i].get_elastic_matrix(
"2Dstress"
) # membrane rigidity matrix with plane stress assumption
Hplane = np.array(
[[Hplane[i][j] for j in [0, 1, 3]] for i in [0, 1, 3]], dtype="object"
)

H[0:3, 0:3] += self.__listThickness[i] * Hplane # Membrane
H[0:3, 0:3] += self.list_thickness[i] * Hplane # Membrane
H[0:3, 3:6] += (
0.5 * (self.__layer[i + 1] ** 2 - self.__layer[i] ** 2) * Hplane
)
Expand Down Expand Up @@ -374,7 +326,7 @@ def get_stress(self, assembly, **kargs):
return StressTensorList(Stress)

def GetStressDistribution(self, assembly, pg, resolution=100):
h = self.GetThickness()
h = self.thickness
z = np.linspace(-h / 2, h / 2, resolution)

Strain = StrainTensorList([0 for i in range(6)])
Expand Down Expand Up @@ -449,8 +401,7 @@ def GetStressDistribution(self, assembly, pg, resolution=100):
return z, Stress

def find_layer(self, position=1):
"""
Returns the num of layer at a given position in the thickness
"""Return the id of layer at a given position in the thickness.

Parameters
----------
Expand All @@ -471,5 +422,5 @@ def find_layer(self, position=1):
), "position should be a float with value in [-1,1]"
if position == -1:
return 0 # 1st layer = bottom layer
z = position * self.GetThickness() / 2
z = position * self.thickness / 2
return list((z - self.__layer) <= 0).index(True) - 1
Loading
Loading