Reputation: 23

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 ?

Upvotes: 2

Answers (1)

Bhargav

Reputation: 4251

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>

Upvotes: 2

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

Answers (1)

Related Questions