Merge b0ebdcf into b8119ab

composites/core.pxd

@@ -11,25 +11,26 @@ cdef extern from "math.h":
 cdef inline double deg2rad(double thetadeg) nogil:
     return thetadeg*4*atan(1.)/180.
 
-cdef class LaminationParameters(object):
+cdef class LaminationParameters:
     cdef public double xiA1, xiA2, xiA3, xiA4
     cdef public double xiB1, xiB2, xiB3, xiB4
     cdef public double xiD1, xiD2, xiD3, xiD4
     cdef public double xiE1, xiE2, xiE3, xiE4
 
-cdef class MatLamina(object):
+cdef class MatLamina:
     cdef public double e1, e2, e3, g12, g13, g23, nu12, nu21, nu13, nu31, nu23, nu32
     cdef public double rho, a1, a2, a3, tref
     cdef public double st1, st2, sc1, sc2, ss12
     cdef public double q11, q12, q13, q21, q22, q23, q31, q32, q33, q44, q55, q66
     cdef public double c11, c12, c13, c22, c23, c33, c44, c55, c66
     cdef public double u1, u2, u3, u4, u5, u6, u7
     cpdef void rebuild(MatLamina)
+    cpdef void trace_normalize_plane_stress(MatLamina)
     cpdef cDOUBLE[:, :] get_constitutive_matrix(MatLamina)
     cpdef cDOUBLE[:, :] get_invariant_matrix(MatLamina)
 
 
-cdef class Lamina(object):
+cdef class Lamina:
     cdef public cINT plyid
     cdef public double h, thetadeg, cost, cos2t, cos4t, sint, sin2t, sin4t
     cdef public double q11L, q12L, q22L, q16L, q26L, q66L, q44L, q45L, q55L
@@ -40,13 +41,13 @@ cdef class Lamina(object):
     cpdef cDOUBLE[:, :] get_transf_matrix_stress_to_lamina(Lamina)
     cpdef cDOUBLE[:, :] get_transf_matrix_stress_to_laminate(Lamina)
 
-cdef class Laminate(object):
+cdef class Laminate:
     cdef public double A11, A12, A16, A22, A26, A66
     cdef public double B11, B12, B16, B22, B26, B66
     cdef public double D11, D12, D16, D22, D26, D66
     cdef public double E44, E45, E55
     cdef public double e1, e2, g12, nu12, nu21
-    cdef public double scf_k13, scf_k23, h, offset, rho, intrho, intrhoz, intrhoz2
+    cdef public double scf_k13, scf_k23, h, offset, intrho, intrhoz, intrhoz2
     cdef public list plies
     cdef public list stack
     cdef cDOUBLE[:, :] get_A(Laminate)
@@ -55,7 +56,6 @@ cdef class Laminate(object):
     cdef cDOUBLE[:, :] get_E(Laminate)
     cdef cDOUBLE[:, :] get_ABD(Laminate)
     cdef cDOUBLE[:, :] get_ABDE(Laminate)
-    cpdef void rebuild(Laminate)
     cpdef void calc_scf(Laminate)
     cpdef void calc_equivalent_properties(Laminate)
     cpdef void calc_constitutive_matrix(Laminate)

composites/core.pyx

@@ -15,19 +15,19 @@ import numpy as np
 INT = np.int64
 DOUBLE = np.float64
 
-cdef class LaminationParameters(object):
+cdef class LaminationParameters:
     r"""Lamination parameters
 
     Attributes
     ----------
     xiA1, xiA2, xiA3, xiA4 : float
-        Lamination parameters `\xi_{Ai}` (in-plane)
+        Lamination parameters :math:`\xi_{Ai}` (in-plane)
     xiB1, xiB2, xiB3, xiB4 : float
-        Lamination parameters `\xi_{Bi}` (in-plane coupling with bending)
+        Lamination parameters :math:`\xi_{Bi}` (in-plane coupling with bending)
     xiD1, xiD2, xiD3, xiD4 : float
-        Lamination parameters `\xi_{Di}` (bending)
+        Lamination parameters :math:`\xi_{Di}` (bending)
     xiE1, xiE2, xiE3, xiE4 : float
-        Lamination parameters `\xi_{Ei}` (transverse shear)
+        Lamination parameters :math:`\xi_{Ei}` (transverse shear)
 
     """
     def __init__(LaminationParameters self):
@@ -36,7 +36,7 @@ cdef class LaminationParameters(object):
         self.xiD1=0; self.xiD2=0; self.xiD3=0; self.xiD4=0
         self.xiE1=0; self.xiE2=0; self.xiE3=0; self.xiE4=0
 
-cdef class MatLamina(object):
+cdef class MatLamina:
     r"""
     Orthotropic material lamina
 
@@ -112,10 +112,11 @@ cdef class MatLamina(object):
 
     Notes
     -----
-    For isotropic materials when the user defines `\nu` and `E`, `G` will be
-    recaculated based on equation: `G = E/(2 \times (1+\nu))`; in a lower
-    priority if the user defines `\nu` and `G`, `E` will be recaculated based
-    on equation: `E = 2 \times (1+\nu) \times G`.
+    For isotropic materials when the user defines :math:`\nu` and :math:`E`,
+    :math:`G` will be recaculated based on equation: :math:`G = E/(2 \times
+    (1+\nu))`; in a lower priority if the user defines :math:`\nu` and
+    :math:`G`, :math:`E` will be recaculated based on equation: :math:`E = 2
+    \times (1+\nu) \times G`.
 
     """
     def __init__(MatLamina self):
@@ -187,10 +188,46 @@ cdef class MatLamina(object):
         self.u2 = (self.q11 - self.q22) / 2.
         self.u3 = (self.q11 + self.q22 - 2*self.q12 - 4*self.q44) / 8.
         self.u4 = (self.q11 + self.q22 + 6*self.q12 - 4*self.q44) / 8.
-        self.u5 = (self.u1-self.u4) / 2.
+        self.u5 = (self.u1 - self.u4) / 2.
         self.u6 = (self.q55 + self.q66) / 2.
         self.u7 = (self.q55 - self.q66) / 2.
 
+    cpdef void trace_normalize_plane_stress(MatLamina self):
+        """Trace-normalize the lamina properties for plane stress
+
+        Modify the original :class:`.MatLamina` object with a
+        trace-normalization performed after calculating the trace according to
+        Eq. 1 of reference:
+
+            Melo, J. D. D., Bi, J., and Tsai, S. W., 2017, “A Novel
+            Invariant-Based Design Approach to Carbon Fiber Reinforced
+            Laminates,” Compos. Struct., 159, pp. 44–52.
+
+        The trace calculated as :math:`tr = Q_{11} + Q_{22} + 2Q_{66}`.  The
+        universal in-plane stress stiffness components
+        :math:`Q_{11},Q_{12},Q_{22},Q_{44},Q_{55},Q_{66}` are divided by
+        :math:`tr`, and the invariants `U_1,U_2,U_3,U_4,U_5,U_6,U_7` are
+        calculated with the normalized stiffnesses, such they also become
+        trace-normalized invariants. These can be accessed using the
+        ``u1,u2,u3,u4,u5,u6,u7`` attributes.
+
+        """
+        cdef double tr
+        tr = self.q11 + self.q22 + 2*self.q66
+        self.q11 /= tr
+        self.q12 /= tr
+        self.q22 /= tr
+        self.q44 /= tr
+        self.q55 /= tr
+        self.q66 /= tr
+        self.u1 /= tr
+        self.u2 /= tr
+        self.u3 /= tr
+        self.u4 /= tr
+        self.u5 /= tr
+        self.u6 /= tr
+        self.u7 /= tr
+
     cpdef cDOUBLE[:, :] get_constitutive_matrix(MatLamina self):
         """Return the constitutive matrix
         """
@@ -217,7 +254,7 @@ cdef class MatLamina(object):
              [ 0,   0, self.u2/2.,   0, -self.u3]], dtype=DOUBLE) # q24
 
 
-cdef class Lamina(object):
+cdef class Lamina:
     """
     Attributes
     ----------
@@ -343,7 +380,7 @@ cdef class Lamina(object):
              [sincos, -sincos, 0, 0, 0, cos2-sin2]], dtype=DOUBLE)
 
 
-cdef class Laminate(object):
+cdef class Laminate:
     r"""
     Attributes
     ----------
@@ -358,25 +395,34 @@ cdef class Laminate(object):
         Equivalent laminate moduli in directions 1 and 2
     g12 : float
         Equivalent laminate shear modulus in the 12 direction
-    nu12, n21 : float
+    nu12, nu21 : float
         Equivalent laminate Poisson ratios in the 12 and 21 directions
     scf_k13, scf_k23 : float
         Shear correction factor in the 13 and 23 directions
     intrho : float
-        Integral `\int_{-h/2+offset}^{+h/2+offset} \rho(z) dz`, used in
+        Integral :math:`\int_{-h/2+offset}^{+h/2+offset} \rho(z) dz`, used in
         equivalent single layer finite element mass matrices
     intrhoz : float
-        Integral `\int_{-h/2+offset}^{+h/2+offset} \rho(z)z dz`, used in
+        Integral :math:`\int_{-h/2+offset}^{+h/2+offset} \rho(z)z dz`, used in
         equivalent single layer finite element mass matrices
     intrhoz2 : float
-        Integral `\int_{-h/2+offset}^{+h/2+offset} \rho(z)z^2 dz`, used in
-        equivalent single layer finite element mass matrices
+        Integral :math:`\int_{-h/2+offset}^{+h/2+offset} \rho(z)z^2 dz`, used
+        in equivalent single layer finite element mass matrices
 
     """
     def __init__(Laminate self):
+        self.h = 0.
+        self.e1 = 0.
+        self.e2 = 0.
+        self.g12 = 0.
+        self.nu12 = 0.
+        self.nu21 = 0.
         self.offset = 0.
         self.scf_k13 = 5/6.
         self.scf_k23 = 5/6.
+        self.intrho = 0.
+        self.intrhoz = 0.
+        self.intrhoz2 = 0.
         self.plies = []
         self.stack = []
 
@@ -410,7 +456,7 @@ cdef class Laminate(object):
                          [self.B12, self.B22, self.B26, self.D12, self.D22, self.D26, 0, 0],
                          [self.B16, self.B26, self.B66, self.D16, self.D26, self.D66, 0, 0],
                          [0, 0, 0, 0, 0, 0, self.E44, self.E45],
-                         [0, 0, 0, 0, 0, 0, self.E44, self.E55]], dtype=DOUBLE)
+                         [0, 0, 0, 0, 0, 0, self.E45, self.E55]], dtype=DOUBLE)
     @property
     def A(self):
         return np.asarray(self.get_A())
@@ -430,18 +476,6 @@ cdef class Laminate(object):
     def ABDE(self):
         return np.asarray(self.get_ABDE())
 
-    cpdef void rebuild(Laminate self):
-        """Update thickness and density"""
-        cdef double rhoh
-        self.h = 0.
-        rhoh = 0.
-        for ply in self.plies:
-            self.h += ply.h
-        for ply in self.plies:
-            ply.rebuild()
-            self.h += ply.h
-            rhoh += ply.matlamina.rho * ply.h
-        self.rho = rhoh / self.h
 
     cpdef void calc_scf(Laminate self):
         """Update shear correction factors of the :class:`.Laminate` object
@@ -583,11 +617,8 @@ cdef class Laminate(object):
     cpdef void force_orthotropic(Laminate self):
         r"""Force an orthotropic laminate
 
-        The attributes
-
-        `A_{16}`, `A_{26}`, `B_{16}`, `B_{26}`, `D_{16}`, `D_{26}`
-
-        are set to zero to force an orthotropic laminate.
+        The attributes :math:`A_{16}`, `A_{26}`, `B_{16}`, `B_{26}`, `D_{16}`,
+        `D_{26}` are set to zero to force an orthotropic laminate.
 
         """
         if self.offset != 0.:
@@ -602,7 +633,7 @@ cdef class Laminate(object):
     cpdef void force_symmetric(Laminate self):
         """Force a symmetric laminate
 
-        The `B_{ij}` terms of the constitutive matrix are set to zero.
+        The :math:`B_{ij}` terms of the constitutive matrix are set to zero.
 
         """
         if self.offset != 0.:
@@ -671,8 +702,8 @@ cdef class Laminate(object):
 cpdef LaminationParameters force_balanced_LP(LaminationParameters lp):
     r"""Force balanced lamination parameters
 
-    The lamination parameters `\xi_{A2}` and `\xi_{A4}` are set to null
-    to force a balanced laminate.
+    The lamination parameters :math:`\xi_{A2}` and :math:`\xi_{A4}` are set to
+    null to force a balanced laminate.
 
     """
     lp.xiA2 = 0
@@ -682,8 +713,8 @@ cpdef LaminationParameters force_balanced_LP(LaminationParameters lp):
 cpdef LaminationParameters force_symmetric_LP(LaminationParameters lp):
     r"""Force symmetric lamination parameters
 
-    The lamination parameters `\xi_{Bi}` are set to null
-    to force a symmetric laminate.
+    The lamination parameters :math:`\xi_{Bi}` are set to null to force a
+    symmetric laminate.
 
     """
     lp.xiB1 = 0
@@ -748,9 +779,9 @@ cpdef Laminate laminate_from_lamination_parameters2(double thickness, MatLamina
     matlamina : :class:`.MatLamina` object
         Material object
     xiA1 to xiD4 : float
-        The 16 lamination parameters `\xi_{A1} \cdots \xi_{A4}`,  `\xi_{B1}
-        \cdots \xi_{B4}`, `\xi_{C1} \cdots \xi_{C4}`,  `\xi_{D1} \cdots
-        \xi_{D4}`, `\xi_{E1} \cdots \xi_{E4}`
+        The 16 lamination parameters :math:`\xi_{A1} \cdots \xi_{A4}`,
+        :math:`\xi_{B1} \cdots \xi_{B4}`, :math:`\xi_{C1} \cdots \xi_{C4}`,
+        :math:`\xi_{D1} \cdots \xi_{D4}`, :math:`\xi_{E1} \cdots \xi_{E4}`
 
 
     Returns