Skip to content

Latest commit

 

History

History
125 lines (89 loc) · 4.75 KB

3.Vectorize_j_integral.md

File metadata and controls

125 lines (89 loc) · 4.75 KB

From Nested Loops to Vectorized Integration: Efficient Computation of the J-Integral

Introduction

In fracture mechanics, the J-integral is a path-independent contour integral used to characterize the intensity of the stress and strain field near the tip of a crack. In numerical implementations (e.g., finite element methods), it is often evaluated along a closed path around the crack tip using discrete integration points. This process naturally involves nested loops: over elements and over Gauss points within each element, with conditional logic applied (e.g., to skip damaged material regions or irrelevant elements).

This article walks through a basic loop-based implementation and demonstrates how the same calculation can be rewritten in a vectorized, broadcast-friendly form using NumPy. This provides a significant speedup while preserving physical correctness.

1. Problem Setup: J-Integral Formulation

The J-integral is given by:

$$ J = \int_{\Gamma} \left( W \delta_{1j} - \sigma_{ij} \frac{\partial u_i}{\partial x_1} \right) n_j , ds $$

Where:

  • $W$: strain energy density
  • $\sigma_{ij}$: stress tensor
  • $\partial u_i / \partial x_j$: displacement gradient
  • $n_j$: normal vector to the path
  • $\Gamma$: integration path around the crack tip

2. Classic Loop-Based Implementation

def compute_j_integral(elements):
    j_total = 0.0
    for elem in elements:
        if not elem.on_integration_path:
            continue
        for gp in elem.gauss_points:
            if gp.state == 'damaged':
                continue
            w = strain_energy(gp)
            sigma = stress_tensor(gp)
            grad_u = displacement_gradient(gp)
            n = normal_vector(gp)

            delta_1j = np.array([1.0, 0.0])
            term1 = w * delta_1j
            term2 = sigma @ grad_u[:, 0]
            integrand = term1 - term2
            j_local = integrand @ n

            j_total += j_local * gp.weight
    return j_total

3. Vectorization Strategy

We observe the mathematical structure of the integral and note that the contribution from each Gauss point is independent. As long as the physical quantities of all valid Gauss points are organized into arrays, the computation can be performed in batch:

  • All valid $W$: array W
  • All valid $\sigma$: tensor array sigma[N, 2, 2]
  • All valid $\nabla u$: tensor array grad_u[N, 2, 2]
  • All valid $n$: array n[N, 2]

With this structure, each term in the integrand can be efficiently computed using NumPy broadcasting or einsum.

4. Vectorized Implementation with NumPy

import numpy as np

# Simulate data: N elements with M Gauss points each
num_elements = 20
gauss_per_elem = 4
total_gps = num_elements * gauss_per_elem

# Boolean mask indicating whether each element is on the integration path
on_path_mask = np.random.rand(num_elements) > 0.5
on_path_mask = np.repeat(on_path_mask, gauss_per_elem)

# Material states: 'intact' or 'damaged'
state = np.random.choice(['intact', 'damaged'], size=total_gps, p=[0.8, 0.2])
valid_mask = (state == 'intact') & on_path_mask

# Physical quantities (simulated)
W = np.random.uniform(1.0, 2.0, size=total_gps)
sigma = np.random.rand(total_gps, 2, 2)
grad_u = np.random.rand(total_gps, 2, 2)
n = np.tile(np.array([1.0, 0.0]), (total_gps, 1))

# Compute integrand
delta_1j = np.array([1.0, 0.0])
term1 = W[:, np.newaxis] * delta_1j
term2 = np.einsum('nij,nj->ni', sigma, grad_u[:, :, 0])
integrand = term1 - term2
j_contrib = np.einsum('ni,ni->n', integrand, n)
weights = np.ones_like(W)

# Final integration
j_total = np.sum(j_contrib[valid_mask] * weights[valid_mask])
print("Vectorized J-integral =", j_total)

5. Comparison

Feature Loop-Based Vectorized
Performance Slower (pure Python loops) Faster (NumPy C-level ops)
Condition handling if/continue Boolean masking
Data layout Matches FE structure Requires flattening and reshaping
GPU/parallel ready Harder to extend Easily extensible with JAX/CUDA

6. Conclusion

Even when calculations involve complex material logic and conditional selection (e.g., skipping damaged regions), vectorization is still possible as long as each integration point is independent. With a bit of restructuring, loop-heavy finite element post-processing can be made far more efficient and ready for GPU acceleration.

This example demonstrates a clear path from numerical physics expression to optimized broadcast implementation — a valuable tool for computational mechanics workflows.