allow for rotational springs in kernel

Author: TchilDill

samples/usage3.ipynb

@@ -307,7 +307,8 @@ def simple_winkler_analysis(model, solver="NR", max_iter: int = 100):
             elif solver == "MNR":
                 pass
 
-            # reset displacements in case of displacement-driven analysis
+            # reset prescribed displacements to converge properly in case 
+            # of displacement-driven analysis
             U[:] = 0.0
 
         # Internal forces calculations with dim(nelem,6,6)

src/openpile/analyses.py

@@ -753,7 +753,7 @@ def create_springs() -> np.ndarray:
                 shape=(self.element_number, 2, 2, spring_dim), dtype=np.float32
             )
             mt = np.zeros(
-                shape=(self.element_number, 2, 2, spring_dim), dtype=np.float32
+                shape=(self.element_number, 2, 2, spring_dim, spring_dim), dtype=np.float32
             )
             Hb = np.zeros(shape=(1, 1, 2, spring_dim), dtype=np.float32)
             Mb = np.zeros(shape=(1, 1, 2, spring_dim), dtype=np.float32)
@@ -766,6 +766,7 @@ def create_springs() -> np.ndarray:
                     (self.soil_properties["x_top [m]"] <= layer.top)
                     & (self.soil_properties["x_bottom [m]"] >= layer.bottom)
                 ].index
+                
                 # py curve
                 if layer.lateral_model.spring_signature[
                     0

src/openpile/construct.py

@@ -15,7 +15,7 @@
 from typing_extensions import Literal
 
 from openpile.construct import Model, Pile
-
+from openpile.core.validation import UserInputError
 from openpile.core import misc
 
 
@@ -34,6 +34,34 @@ def reverse_indices(A, B):
     return np.array([x for x in A if x not in B])
 
 
+@njit(f8[:](f8[:]), cache=True)
+def double_inner_njit(in_arr):
+    # d = misc.repeat_inner(u)
+    out_arr = np.zeros(((len(in_arr) - 1) * 2), dtype=np.float64)
+    i = 0
+    while i < len(out_arr) + 1:
+        if i == 0:
+            out_arr[i] = in_arr[0]
+            i += 1
+        elif i == len(out_arr):
+            out_arr[i] = in_arr[-1]
+            i += 1
+        elif i == len(out_arr) - 1:
+            out_arr[i] = in_arr[-2]
+            i += 1
+        elif i % 2 != 0:
+            x = int((i + 1) / 2)
+            out_arr[i] = in_arr[x]
+            out_arr[i + 1] = in_arr[x]
+            i += 2
+        elif i % 2 == 0:
+            x = int(i / 2)
+            out_arr[i] = in_arr[x]
+            out_arr[i + 1] = in_arr[x]
+            i += 2
+    
+    return out_arr
+
 @njit(parallel=True, cache=True)
 def numba_ix(arr, rows, cols):
     """
@@ -279,7 +307,7 @@ def elem_py_stiffness_matrix(model, u, kind):
     L = mesh_to_element_length(model)
 
     # calculate the spring stiffness
-    ksoil = calculate_springs_stiffness(u=u[1::3], springs=model._py_springs, kind=kind)
+    ksoil = calculate_py_springs_stiffness(u=u[1::3], springs=model._py_springs, kind=kind)
 
     N = 0 * L
     A = 2 * L / 7
@@ -323,6 +351,73 @@ def elem_py_stiffness_matrix(model, u, kind):
     return ktop + kbottom
 
 
+def elem_mt_stiffness_matrix(model, u, kind):
+    """creates soil element stiffness matrix based on model info and element number.
+
+    The soil stiffness matrix assumes the soil stiffness to vary lineraly along the elements.
+
+    #TODO: proof here
+
+    Parameters
+    ----------
+    model : openpile class object `openpile.construct.Model`
+        includes information on soil/structure, elements, nodes and other model-related data.
+    u: np.ndarray
+        Global displacement vector
+    kind: str
+        "initial", "secant" or "tangent"
+
+    Returns
+    -------
+    k: numpy array (3d)
+        soil consistent stiffness matrix of all elememts related to the p-y soil springs' stiffness
+
+    Raises
+    ------
+    ValueError
+        ndof per node can be either 2 or 3
+
+    """
+
+    # calculate length vector
+    L = mesh_to_element_length(model)
+
+    # calculate the py spring stiffness
+    kspy = calculate_py_springs_stiffness(u=u[1::3], springs=model._py_springs, kind='secant')
+    upy = double_inner_njit(u[1::3]).reshape(-1,2,1,1)
+    p_mobilised = kspy * upy
+
+    # calculate the spring stiffness
+    ksoil = calculate_mt_springs_stiffness(u=u[2::3], 
+                                           mt_springs=model._mt_springs, 
+                                           py_springs=model._py_springs,
+                                           p_mobilised=p_mobilised,
+                                           kind=kind,
+                                           )
+
+    N = 0 * L
+    A = 6 / (5*L)
+    B = N + 1 / 10 
+    C = 2*L / 15
+    D = -L / 30
+
+    km = (
+        np.block(
+            [
+                [N, N, N, N, N, N],
+                [N, A, B, N,-A, B],
+                [N, B, C, N,-B, D],
+                [N, N, N, N, N, N],
+                [N,-A,-B, N, A,-B],
+                [N, B, D, N,-B, C],
+            ]
+        )
+        * (0.5*ksoil[:,0] + 0.5*ksoil[:,1])
+    )
+
+    return km
+
+
 @njit(parallel=True, cache=True)
 def jit_build(k, ndim, n_elem, node_per_element, ndof_per_node):
     # pre-allocate stiffness matrix
@@ -342,7 +437,7 @@ def jit_build(k, ndim, n_elem, node_per_element, ndof_per_node):
     return K
 
 
-def build_stiffness_matrix(model, u=None, kind=None):
+def build_stiffness_matrix(model, u=None, kind=None, p_mobilised=None):
     """Builds the stiffness matrix based on the model(element and node) properties
 
     Element stiffness matrices are first computed for each element and then loaded in the global stiffness matrix through summation.
@@ -354,6 +449,8 @@ def build_stiffness_matrix(model, u=None, kind=None):
         stores all information about nodes elements, and other model-related items, see #TODO:link to model class
     u: np.ndarray, Optional
         Global displacement vector. Must be given if soil_flag is not None
+    p_mobilised: np.ndarray, Optional
+        Mobilised p-value. Must be given if soil_flag is not None
     kind: str, Optional
         "initial", "secant" or "tangent". Must be given if soil_flag is not None
 
@@ -381,12 +478,17 @@ def build_stiffness_matrix(model, u=None, kind=None):
             UserWarning("'u' and 'kind' must be stipulated.")
         else:
             if model.distributed_lateral:
-                k += elem_py_stiffness_matrix(model, u, kind)
-            elif model.distributed_moment:
-                k += 0
-            elif model.base_shear:
+                k += elem_py_stiffness_matrix(model, u, kind) 
+            
+            if model.distributed_moment:
+                if not model.distributed_lateral:
+                    raise UserInputError('Error: Distributed moment cannot be calculated without distributed lateral springs.')
+                k += elem_mt_stiffness_matrix(model, u, kind) 
+            
+            if model.base_shear:
                 k += 0
-            elif model.base_moment:
+            
+            if model.base_moment:
                 k += 0
 
     K = jit_build(k, ndim_global, n_elem, node_per_element, ndof_per_node)
@@ -434,7 +536,7 @@ def struct_internal_force(model, u) -> np.ndarray:
         if model.distributed_lateral:
             k += elem_py_stiffness_matrix(model, u, kind)
         elif model.distributed_moment:
-            k += 0
+            k += elem_mt_stiffness_matrix(model, u, kind)
         elif model.base_shear:
             k += 0
         elif model.base_moment:
@@ -451,20 +553,19 @@ def struct_internal_force(model, u) -> np.ndarray:
 
 
 @njit(cache=True)
-def calculate_springs_stiffness(
+def calculate_py_springs_stiffness(
     u: np.ndarray, springs: np.ndarray, kind: Literal["initial", "secant", "tangent"]
 ):
-    """Calculate springs stiffness
+    """Calculate springs stiffness for py or t-z springs.
 
     Parameters
     ----------
     u : np.ndarray
         displacements to calculate stiffness.
         For dofs related to t-z curves, u = U[::3] where U is the global displacement vector.
         For dofs related to p-y curves, u = U[1::3] where U is the global displacement vector.
-        For dofs related to m-t curves, u = U[2::3] where U is the global displacement vector.
     springs : np.ndarray
-        soil-structure interaction springs array of shape (n_elem, 2, 2, spring_dim)
+        soil-structure interaction py springs array of shape (n_elem, 2, 2, spring_dim)
     kind : str
         defines whether it is initial, secant of tangent stiffness to define
 
@@ -474,30 +575,8 @@ def calculate_springs_stiffness(
         secant or tangent stiffness for all elements. Array of shape(n_elem,2,1,1)
     """
 
-    # double inner values
-    # d = misc.repeat_inner(u)
-    d = np.zeros(((len(u) - 1) * 2), dtype=np.float64)
-    i = 0
-    while i < len(d) + 1:
-        if i == 0:
-            d[i] = u[0]
-            i += 1
-        elif i == len(d):
-            d[i] = u[-1]
-            i += 1
-        elif i == len(d) - 1:
-            d[i] = u[-2]
-            i += 1
-        elif i % 2 != 0:
-            x = int((i + 1) / 2)
-            d[i] = u[x]
-            d[i + 1] = u[x]
-            i += 2
-        elif i % 2 == 0:
-            x = int(i / 2)
-            d[i] = u[x]
-            d[i + 1] = u[x]
-            i += 2
+    # double inner values for u 
+    d = double_inner_njit(u)
 
     # displacemet with same dimension as spring
     d = np.abs(d).reshape((-1, 2, 1, 1))
@@ -535,11 +614,96 @@ def calculate_springs_stiffness(
                             d[i, j, 0, 0], springs[i, j, 1], springs[i, j, 0]
                         )
 
-            k[i, j, 0, 0] = abs((p1 - p0) / dx)
+                k[i, j, 0, 0] = abs((p1 - p0) / dx)
 
     return k
 
 
+
+# @njit(cache=True)
+def calculate_mt_springs_stiffness(
+    u: np.ndarray, mt_springs: np.ndarray, py_springs:np.ndarray, p_mobilised: np.ndarray, kind: Literal["initial", "secant", "tangent"]
+):
+    """Calculate springs stiffness for rotational springs
+
+    .. note::
+        The difference with the py function is that rotational springs can depend on lateral springs
+
+    Parameters
+    ----------
+    u : np.ndarray
+        displacements to calculate stiffness.
+        For dofs related to m-t curves, u = U[2::3] where U is the global displacement vector.
+    mt_springs : np.ndarray
+        soil-structure interaction m-t springs array of shape (n_elem, 2, 2, py_spring_dim, mt_spring_dim)
+    py_springs : np.ndarray
+        soil-structure interaction p-y springs array of shape (n_elem, 2, 2, py_spring_dim)
+    p_mobilised : np.ndarray
+        current p value of p-y springs of shape (nelem, 2, 1, 1)
+    kind : str
+        defines whether it is initial, secant of tangent stiffness to define
+
+    Returns
+    -------
+    k: np.ndarray
+        secant or tangent stiffness for all elements. Array of shape(n_elem,2,1,1)
+    """
+
+    # double inner values for u 
+    d = double_inner_njit(u)
+
+    # displacemet and p_mobilised with same dimension as spring
+    d = np.abs(d).reshape((-1, 2, 1, 1))
+    p = p_mobilised.reshape((-1, 2, 1, 1))
+
+    k = np.zeros(d.shape, dtype=np.float64)
+
+    # i for each element
+    for i in range(k.shape[0]):
+        # j for top and bottom values of spring
+        for j in range(k.shape[1]):
+            if np.sum(mt_springs[i, j, 1, 0]) == 0:
+                #check if first m-vector is not a defined spring
+                #if that is the case, we do not calculate any stiffness
+                pass
+            else:
+                # get the proper m-t spring 
+                m = np.zeros(k.shape[4])
+                t = np.zeros(k.shape[4])
+
+                for ii in range((k.shape[4])):
+                    if ii == 0:
+                        pass
+                    else:
+                        m[ii] = np.interp(p[i,j,1,1], py_springs[i,j,0,:], mt_springs[i,j,0,:,ii])
+                        t[ii] = np.interp(p[i,j,1,1], py_springs[i,j,0,:], mt_springs[i,j,1,:,ii])
+                
+                if kind == "initial" or d[i, j, 0, 0] == 0.0:
+                    dt = t[1] - t[0]
+                    m0 = m[0]
+                    m1 = m[1]
+                elif kind == "secant":
+                    dt = d[i, j, 0, 0]
+                    m0 = m[0]
+                    if d[i, j, 0, 0] > t[-1]:
+                        m1 = m[-1]
+                    else:
+                        m1 = np.interp(dt, t, m)
+                elif kind == "tangent":
+                    dt = min(0.0005, d[i, j, 0, 0])
+                    if (d[i, j, 0, 0] - dt) > t[-1]:
+                        m0 = m[-1]
+                    else:
+                        m0 = np.interp(d[i, j, 0, 0] - dt, t, m)
+                    if (d[i, j, 0, 0] - dt) > t[-1]:
+                        m1 = m[-1]
+                    else:
+                        m1 = np.interp(d[i, j, 0, 0], t, m)
+
+                k[i, j, 0, 0] = abs((m1 - m0) / dt)
+
+    return k
+
 def computer():
     """This function is the solver of openpile.