Search code examples
pythonnumpymatrix

Python numpy how to split a matrix into 4 not equal matrixes?


from sympy import *
import numpy as np

# Python 3.13.2
u1 = Symbol("u1")
u2 = Symbol("u2")
q3 = Symbol("q3")
u4 = Symbol("u4")
q5 = Symbol("q5")

disp_vector = np.array([u1, u2, q3, u4, q5])

stiffness_matrix = np.array([[1, 0, 0, -1, 0],
                             [0, 0.12, 0.6, 0, 0.6],
                             [0, 0.6, 4, 0, 2],
                             [-1, 0, 0, 1, 0],
                             [0, 0.6, 2, 0, 4]])

force_vector = np.array([0, 40, -26.7, 0, 0])

I am trying to code static condensation. This allows me to reduce the size of the stiffness matrix above but for that I need to be able to divide the stiffness matrix into following matrices.

krr = np.array([[1, 0, 0],
                [0, 0.12, 0.6],
                [0, 0.6, 4]])

krc = np.array([[-1, 0],
                [0, 0.6],
                [0, 2]])

kcr = np.array([[-1, 0, 0],
               [0, 0.6, 2]])

kcc = np.array([[1, 0],
               [0, 4]])

How would I do this ?


Solution

  • You need to divide your stiffness matrix into four submatrices based on the degrees of freedom you want to r and c. So to retain the first 3 rows/columns and condense the last 2

    krr: The top left block
    krc: The top right block
    kcr: The bottom left block
    kcc: The bottom right block
    

    The code

    n_retain = 3  
    n_condense = 2  
    
    
    krr = stiffness_matrix[:n_retain, :n_retain]
    krc = stiffness_matrix[:n_retain, n_retain:]
    kcr = stiffness_matrix[n_retain:, :n_retain]
    kcc = stiffness_matrix[n_retain:, n_retain:]
    
    print("krr =\n", krr)
    print("\nkrc =\n", krc)
    print("\nkcr =\n", kcr)
    print("\nkcc =\n", kcc)
    

    submatrices

    krr =
     [[1.   0.   0.  ]
     [0.   0.12 0.6 ]
     [0.   0.6  4.  ]]
    
    krc =
     [[-1.   0. ]
     [ 0.   0.6]
     [ 0.   2. ]]
    
    kcr =
     [[-1.   0.   0. ]
     [ 0.   0.6  2. ]]
    
    kcc =
     [[1. 0.]
     [0. 4.]]
    

    Animations Credits - https://ocw.mit.edu/courses/1-571-structural-analysis-and-control-spring-2004/pages/readings/

    <!DOCTYPE html>
    <html>
    <head>
      <style>
        body {
          font-family: Arial, sans-serif;
          margin: 20px;
          background-color: #f5f5f5;
        }
        .container {
          max-width: 800px;
          margin: 0 auto;
          background-color: white;
          padding: 20px;
          border-radius: 8px;
          box-shadow: 0 2px 10px rgba(0,0,0,0.1);
        }
        h1, h2 {
          color: #2c3e50;
          text-align: center;
        }
        .matrix-container {
          display: flex;
          justify-content: center;
          align-items: center;
          margin: 20px 0;
          flex-wrap: wrap;
        }
        .matrix {
          border: 2px solid #3498db;
          margin: 10px;
          padding: 5px;
          background-color: white;
          border-radius: 4px;
          transition: all 0.5s ease;
        }
        .matrix-row {
          display: flex;
          justify-content: center;
        }
        .matrix-cell {
          width: 40px;
          height: 40px;
          display: flex;
          align-items: center;
          justify-content: center;
          margin: 2px;
          font-weight: bold;
          transition: all 0.5s ease;
        }
        .retained {
          background-color: #2ecc71;
          color: white;
        }
        .condensed {
          background-color: #e74c3c;
          color: white;
        }
        .coupling {
          background-color: #f39c12;
          color: white;
        }
        .controls {
          display: flex;
          justify-content: center;
          margin: 20px 0;
        }
        button {
          background-color: #3498db;
          color: white;
          border: none;
          padding: 10px 20px;
          margin: 0 10px;
          border-radius: 4px;
          cursor: pointer;
          font-size: 16px;
          transition: background-color 0.3s;
        }
        button:hover {
          background-color: #2980b9;
        }
        button:disabled {
          background-color: #bdc3c7;
          cursor: not-allowed;
        }
        .explanation {
          margin: 20px 0;
          padding: 15px;
          background-color: #f8f9fa;
          border-left: 4px solid #3498db;
          border-radius: 4px;
        }
        .equation {
          font-family: 'Times New Roman', Times, serif;
          font-style: italic;
          text-align: center;
          margin: 15px 0;
          font-size: 18px;
        }
        .force-vector {
          display: flex;
          flex-direction: column;
          align-items: center;
          margin: 10px;
        }
        .force-cell {
          width: 40px;
          height: 30px;
          display: flex;
          align-items: center;
          justify-content: center;
          margin: 2px;
          font-weight: bold;
          border: 1px solid #3498db;
          background-color: #fff;
          transition: all 0.5s ease;
        }
        .progress {
          width: 100%;
          height: 8px;
          background-color: #ecf0f1;
          margin-bottom: 20px;
          border-radius: 4px;
          overflow: hidden;
        }
        .progress-bar {
          height: 100%;
          background-color: #3498db;
          width: 0%;
          transition: width 0.3s ease;
        }
      </style>
    </head>
    <body>
      <div class="container">
        <h1>Static Condensation Animation</h1>
        
        <div class="progress">
          <div class="progress-bar" id="progressBar"></div>
        </div>
        
        <div class="explanation" id="explanation">
          <p>Static condensation is a technique used in structural analysis to reduce the size of the stiffness matrix by eliminating degrees of freedom (DOFs) that are not of primary interest.</p>
          <p>Click "Next" to walk through the process step by step.</p>
        </div>
        
        <div class="matrix-container" id="matrixDisplay">
          <!-- Matrices will be displayed here -->
        </div>
        
        <div class="controls">
          <button id="prevBtn" disabled>Previous</button>
          <button id="nextBtn">Next</button>
          <button id="resetBtn">Reset</button>
        </div>
      </div>
    
      <script>
        // Original stiffness matrix and force vector
        const stiffnessMatrix = [
          [1, 0, 0, -1, 0],
          [0, 0.12, 0.6, 0, 0.6],
          [0, 0.6, 4, 0, 2],
          [-1, 0, 0, 1, 0],
          [0, 0.6, 2, 0, 4]
        ];
        
        const forceVector = [0, 40, -26.7, 0, 0];
        
        // Steps for the animation
        const steps = [
          {
            title: "Original System",
            explanation: "We start with the full stiffness matrix K and force vector F. Our goal is to solve Ku = F for the displacement vector u.",
            showMatrices: ["K", "F"]
          },
          {
            title: "Partition the System",
            explanation: "We partition the matrix into 4 submatrices: Krr (retained DOFs), Kcc (condensed DOFs), and Krc/Kcr (coupling terms). The force vector is similarly split into Fr and Fc.",
            showMatrices: ["K_partitioned", "F_partitioned"]
          },
          {
            title: "Extract Submatrices",
            explanation: "We extract the four submatrices explicitly: Krr (3×3), Krc (3×2), Kcr (2×3), and Kcc (2×2).",
            showMatrices: ["Krr", "Krc", "Kcr", "Kcc", "Fr", "Fc"]
          },
          {
            title: "Condensation Equations",
            explanation: "The static condensation process can be represented by these equations:\n\nCondensed stiffness matrix: K* = Krr - Krc × Kcc⁻¹ × Kcr\nCondensed force vector: F* = Fr - Krc × Kcc⁻¹ × Fc",
            showMatrices: ["equation"]
          },
          {
            title: "Condensed System",
            explanation: "We now have a reduced system K* u_r = F* that only involves the retained DOFs. This smaller system is easier to solve.",
            showMatrices: ["K_condensed", "F_condensed"]
          },
          {
            title: "Solve for Retained DOFs",
            explanation: "We solve the condensed system to find the values of the retained DOFs (u_r).",
            showMatrices: ["u_r"]
          },
          {
            title: "Back-Calculate Condensed DOFs",
            explanation: "Once we have u_r, we can back-calculate the condensed DOFs (u_c) using: u_c = Kcc⁻¹ × (Fc - Kcr × u_r)",
            showMatrices: ["u_r", "u_c"]
          },
          {
            title: "Complete Solution",
            explanation: "Finally, we have the complete solution vector u combining u_r and u_c. We've solved the full system by working with a reduced matrix!",
            showMatrices: ["u_full"]
          }
        ];
        
        // Function to create matrix display
        function createMatrixDisplay(matrix, title, className = "") {
          const matrixDiv = document.createElement("div");
          matrixDiv.className = `matrix ${className}`;
          
          const titleDiv = document.createElement("div");
          titleDiv.textContent = title;
          titleDiv.style.textAlign = "center";
          titleDiv.style.marginBottom = "5px";
          matrixDiv.appendChild(titleDiv);
          
          for (let i = 0; i < matrix.length; i++) {
            const rowDiv = document.createElement("div");
            rowDiv.className = "matrix-row";
            
            if (Array.isArray(matrix[i])) {
              // It's a 2D matrix
              for (let j = 0; j < matrix[i].length; j++) {
                const cellDiv = document.createElement("div");
                cellDiv.className = "matrix-cell";
                cellDiv.textContent = typeof matrix[i][j] === 'number' ? 
                  matrix[i][j].toFixed(2).replace(/\.00$/, "") : matrix[i][j];
                
                // Apply colors based on partitioning for the full matrix
                if (title === "K (partitioned)") {
                  if (i < 3 && j < 3) cellDiv.classList.add("retained");
                  else if (i >= 3 && j >= 3) cellDiv.classList.add("condensed");
                  else cellDiv.classList.add("coupling");
                }
                
                rowDiv.appendChild(cellDiv);
              }
            } else {
              // It's a vector
              const cellDiv = document.createElement("div");
              cellDiv.className = "force-cell";
              cellDiv.textContent = typeof matrix[i] === 'number' ? 
                matrix[i].toFixed(2).replace(/\.00$/, "") : matrix[i];
              
              // Apply colors for partitioned force vector
              if (title === "F (partitioned)") {
                if (i < 3) cellDiv.classList.add("retained");
                else cellDiv.classList.add("condensed");
              }
              
              rowDiv.appendChild(cellDiv);
            }
            
            matrixDiv.appendChild(rowDiv);
          }
          
          return matrixDiv;
        }
        
        // Function to create force vector display
        function createForceVector(vector, title, className = "") {
          const vectorDiv = document.createElement("div");
          vectorDiv.className = `force-vector ${className}`;
          
          const titleDiv = document.createElement("div");
          titleDiv.textContent = title;
          titleDiv.style.textAlign = "center";
          titleDiv.style.marginBottom = "5px";
          vectorDiv.appendChild(titleDiv);
          
          for (let i = 0; i < vector.length; i++) {
            const cellDiv = document.createElement("div");
            cellDiv.className = "force-cell";
            cellDiv.textContent = typeof vector[i] === 'number' ? 
              vector[i].toFixed(2).replace(/\.00$/, "") : vector[i];
            
            // Apply colors for partitioned force vector
            if (title === "F (partitioned)") {
              if (i < 3) cellDiv.classList.add("retained");
              else cellDiv.classList.add("condensed");
            }
            
            vectorDiv.appendChild(cellDiv);
          }
          
          return vectorDiv;
        }
        
        // Function to display an equation
        function createEquation(text) {
          const eqDiv = document.createElement("div");
          eqDiv.className = "equation";
          eqDiv.innerHTML = text;
          return eqDiv;
        }
        
        // Get DOM elements
        const matrixDisplay = document.getElementById("matrixDisplay");
        const explanation = document.getElementById("explanation");
        const prevBtn = document.getElementById("prevBtn");
        const nextBtn = document.getElementById("nextBtn");
        const resetBtn = document.getElementById("resetBtn");
        const progressBar = document.getElementById("progressBar");
        
        let currentStep = 0;
        
        // Submatrices
        const krr = [
          [1, 0, 0],
          [0, 0.12, 0.6],
          [0, 0.6, 4]
        ];
        
        const krc = [
          [-1, 0],
          [0, 0.6],
          [0, 2]
        ];
        
        const kcr = [
          [-1, 0, 0],
          [0, 0.6, 2]
        ];
        
        const kcc = [
          [1, 0],
          [0, 4]
        ];
        
        const fr = [0, 40, -26.7];
        const fc = [0, 0];
        
        // Solution (mockup values for demonstration)
        const ur = [0.78, 345.58, -9.50];
        const uc = [-0.78, 1.49];
        const uFull = [0.78, 345.58, -9.50, -0.78, 1.49];
        
        // Condensed matrices (mockup values for demonstration)
        const kCondensed = [
          [2, 0, 0],
          [0, 0.21, 0.6],
          [0, 0.6, 3.5]
        ];
        
        const fCondensed = [0, 40, -26.7];
        
        // Update display based on current step
        function updateDisplay() {
          // Update progress bar
          progressBar.style.width = `${(currentStep / (steps.length - 1)) * 100}%`;
          
          // Clear previous content
          matrixDisplay.innerHTML = "";
          
          // Set explanation
          explanation.innerHTML = `<h2>${steps[currentStep].title}</h2><p>${steps[currentStep].explanation}</p>`;
          
          // Display relevant matrices
          steps[currentStep].showMatrices.forEach(matrixName => {
            switch (matrixName) {
              case "K":
                matrixDisplay.appendChild(createMatrixDisplay(stiffnessMatrix, "K (Stiffness Matrix)"));
                break;
              case "F":
                matrixDisplay.appendChild(createForceVector(forceVector, "F (Force Vector)"));
                break;
              case "K_partitioned":
                matrixDisplay.appendChild(createMatrixDisplay(stiffnessMatrix, "K (partitioned)"));
                break;
              case "F_partitioned":
                matrixDisplay.appendChild(createForceVector(forceVector, "F (partitioned)"));
                break;
              case "Krr":
                matrixDisplay.appendChild(createMatrixDisplay(krr, "Krr (retained)", "retained"));
                break;
              case "Krc":
                matrixDisplay.appendChild(createMatrixDisplay(krc, "Krc (coupling)", "coupling"));
                break;
              case "Kcr":
                matrixDisplay.appendChild(createMatrixDisplay(kcr, "Kcr (coupling)", "coupling"));
                break;
              case "Kcc":
                matrixDisplay.appendChild(createMatrixDisplay(kcc, "Kcc (condensed)", "condensed"));
                break;
              case "Fr":
                matrixDisplay.appendChild(createForceVector(fr, "Fr (retained)", "retained"));
                break;
              case "Fc":
                matrixDisplay.appendChild(createForceVector(fc, "Fc (condensed)", "condensed"));
                break;
              case "equation":
                const eqDiv = document.createElement("div");
                eqDiv.style.width = "100%";
                eqDiv.style.textAlign = "center";
                
                const eq1 = document.createElement("div");
                eq1.className = "equation";
                eq1.innerHTML = "K* = K<sub>rr</sub> - K<sub>rc</sub> · K<sub>cc</sub><sup>-1</sup> · K<sub>cr</sub>";
                
                const eq2 = document.createElement("div");
                eq2.className = "equation";
                eq2.innerHTML = "F* = F<sub>r</sub> - K<sub>rc</sub> · K<sub>cc</sub><sup>-1</sup> · F<sub>c</sub>";
                
                eqDiv.appendChild(eq1);
                eqDiv.appendChild(eq2);
                matrixDisplay.appendChild(eqDiv);
                break;
              case "K_condensed":
                matrixDisplay.appendChild(createMatrixDisplay(kCondensed, "K* (Condensed Stiffness Matrix)"));
                break;
              case "F_condensed":
                matrixDisplay.appendChild(createForceVector(fCondensed, "F* (Condensed Force Vector)"));
                break;
              case "u_r":
                matrixDisplay.appendChild(createForceVector(ur, "ur (Retained DOFs)", "retained"));
                break;
              case "u_c":
                matrixDisplay.appendChild(createForceVector(uc, "uc (Condensed DOFs)", "condensed"));
                break;
              case "u_full":
                matrixDisplay.appendChild(createForceVector(uFull, "u (Complete Solution)"));
                break;
            }
          });
          
          // Update button states
          prevBtn.disabled = currentStep === 0;
          nextBtn.disabled = currentStep === steps.length - 1;
        }
        
        // Event listeners
        prevBtn.addEventListener("click", () => {
          if (currentStep > 0) {
            currentStep--;
            updateDisplay();
          }
        });
        
        nextBtn.addEventListener("click", () => {
          if (currentStep < steps.length - 1) {
            currentStep++;
            updateDisplay();
          }
        });
        
        resetBtn.addEventListener("click", () => {
          currentStep = 0;
          updateDisplay();
        });
        
        // Initialize display
        updateDisplay();
      </script>
    </body>
    </html>