2D FEM Solver with Periodic Boundary Conditions and Rigid Body Constraints

A complete Python implementation of a 2D Finite Element Method (FEM) solver for linear elasticity problems on square meshes with periodic boundary conditions and rigid body movement constraints.

Features

• Structured Quadrilateral Mesh Generation: Automatic generation of Q4 (4-node quadrilateral) meshes for square domains
• 2D Linear Elasticity: Plane stress formulation with full stiffness matrix assembly
• Periodic Boundary Conditions: Enforces periodicity on opposite boundaries of the domain
• Rigid Body Constraints: Eliminates rigid body modes (translation and rotation) to ensure unique solutions
• Multiple Constraint Methods: Supports both penalty method and Lagrange multipliers
• Comprehensive Visualization: Built-in tools for plotting meshes, deformations, displacements, and stresses
• Flexible Loading: Supports both uniform and spatially-varying body forces
• 🆕 Multi-Material Support: Element-wise material properties for composite analysis
• 🆕 RVE Homogenization: Computational homogenization for fiber-reinforced composites with circular inclusions
• 🆕 Effective Properties: Automatic computation of effective stiffness matrix and engineering constants

Mathematical Formulation

Governing Equations

The solver implements the weak form of 2D linear elasticity:

Find u such that: ∫_Ω B^T D B u dΩ = ∫_Ω N^T f dΩ

where:

• B is the strain-displacement matrix
• D is the material constitutive matrix (plane stress)
• N are the shape functions
• f is the body force vector

Periodic Boundary Conditions

Periodic BC enforces:

u(left) = u(right)
u(bottom) = u(top)

This is implemented using either:

1. Penalty method: Adds penalty terms to enforce constraints approximately
2. Lagrange multipliers: Adds constraint equations to the system exactly

Rigid Body Constraints

To prevent singular systems, rigid body modes are eliminated by fixing:

• Translation in x and y directions
• Rotation about a reference point

This is achieved by constraining specific degrees of freedom at corner nodes.

Installation

Requirements

• Python 3.7+
• NumPy
• SciPy
• Matplotlib

Install dependencies:

pip install -r requirements.txt

Project Structure

test_fem2d/
├── mesh_generator.py                # Mesh generation for square domains
├── fem_assembly.py                  # FEM stiffness matrix and load vector assembly
├── fem_assembly_multimaterial.py    # Multi-material FEM assembly
├── boundary_conditions.py           # Periodic BC and rigid body constraints
├── fem_solver.py                    # Main solver class
├── rve_generator.py                 # RVE generation with circular fiber
├── homogenization.py                # Computational homogenization module
├── visualize.py                     # Visualization tools (matplotlib)
├── svg_export.py                    # SVG export functionality
├── validation.py                    # Validation tests with analytical solutions
├── example_simple.py                # Simple example with uniform load
├── example_advanced.py              # Advanced example with spatially-varying loads
├── example_rve_homogenization.py    # RVE homogenization example
├── test_both_methods.py             # Test penalty vs Lagrange multiplier methods
├── requirements.txt                 # Python dependencies
└── README.md                        # This file

Usage

Quick Start

Run a simple example:

python example_simple.py

This will:

1. Create an 8×8 element mesh
2. Apply periodic boundary conditions
3. Add rigid body constraints
4. Solve with uniform body force
5. Generate visualization plots

Basic Usage

from fem_solver import FEMSolver
from visualize import FEMVisualizer

# Create solver
solver = FEMSolver(nx=10, ny=10, lx=1.0, ly=1.0, E=1000.0, nu=0.3)

# Setup boundary conditions
solver.setup_periodic_bc()
solver.setup_rigid_body_constraints()

# Solve with uniform body force
body_force = [0.0, -10.0]  # Force in -y direction
displacement = solver.solve(body_force=body_force, method='lagrange')

# Compute stresses
stress, strain = solver.compute_stress_strain()

# Visualize
vis = FEMVisualizer(solver)
vis.plot_deformed(scale=50.0)
vis.plot_stress_field(component='von_mises')

Spatially-Varying Loads

Define a function for spatially-varying body forces:

def my_load(x, y):
    """Custom body force function"""
    fx = 10.0 * np.sin(2 * np.pi * x)
    fy = -20.0 * np.cos(2 * np.pi * y)
    return np.array([fx, fy])

# Solve with custom load
displacement = solver.solve(body_force=my_load, method='lagrange')

Examples

Example 1: Simple Uniform Load

python example_simple.py

Creates a comprehensive visualization showing:

• Original mesh
• Periodic boundary condition pairs
• Deformed mesh
• Displacement fields (x, y components)
• von Mises stress distribution

Example 2: Advanced Spatially-Varying Loads

python example_advanced.py

Demonstrates two types of complex loading:

1. Sinusoidal load: Wave-like force distribution
2. Vortex load: Circular force pattern

Shows comparative results for both cases.

Example 3: RVE Homogenization for Composites

python example_rve_homogenization.py

Advanced example demonstrating computational homogenization for fiber-reinforced composites:

• Creates dense RVE mesh (40×40 or finer) with circular fiber inclusion
• Uses voxel-based material assignment (elements entirely within fiber get fiber properties)
• Applies periodic boundary conditions on all edges
• Computes effective stiffness matrix through three independent load cases
• Extracts effective engineering constants (E_x, E_y, ν_xy, G_xy)
• Compares results with theoretical bounds (Voigt and Reuss)
• Includes convergence study and volume fraction parametric study

Key capabilities:

• Material Distribution: Visualizes matrix and fiber phases
• Effective Properties: Computes homogenized material constants
• Validation: Compares with Voigt (upper) and Reuss (lower) bounds
• Parametric Studies: Mesh convergence and volume fraction effects

Module Documentation

mesh_generator.py

SquareMesh: Generates structured quadrilateral meshes

mesh = SquareMesh(nx=10, ny=10, lx=1.0, ly=1.0)

Key attributes:

• nodes: Node coordinates (n_nodes × 2)
• elements: Element connectivity (n_elements × 4)
• left_nodes, right_nodes, top_nodes, bottom_nodes: Boundary node sets
• get_periodic_pairs(): Returns master-slave node pairs for periodic BC

fem_assembly.py

FEMAssembler: Assembles global stiffness matrix and load vector

assembler = FEMAssembler(mesh, E=1000.0, nu=0.3, thickness=1.0)
K = assembler.assemble_stiffness_matrix()
F = assembler.assemble_load_vector(body_force=[0.0, -10.0])

Features:

• Q4 element formulation with 2×2 Gauss quadrature
• Plane stress material model
• Efficient sparse matrix assembly

boundary_conditions.py

BoundaryConditions: Applies periodic BC and rigid body constraints

bc = BoundaryConditions(mesh, n_dofs)
bc.add_periodic_bc()
bc.add_rigid_body_constraints()

# Method 1: Penalty method
K_mod, F_mod = bc.apply_constraints_penalty(K, F, penalty=1e10)

# Method 2: Lagrange multipliers
K_aug, F_aug, n_const = bc.apply_constraints_lagrange(K, F)

fem_solver.py

FEMSolver: Main solver combining all components

solver = FEMSolver(nx=10, ny=10, lx=1.0, ly=1.0, E=1000.0, nu=0.3)
solver.setup_periodic_bc()
solver.setup_rigid_body_constraints()
displacement = solver.solve(body_force=[0.0, -10.0], method='lagrange')
stress, strain = solver.compute_stress_strain()

visualize.py

FEMVisualizer: Visualization tools (matplotlib-based)

vis = FEMVisualizer(solver)
vis.plot_mesh()                           # Plot mesh
vis.plot_periodic_bc()                    # Show periodic pairs
vis.plot_deformed(scale=50.0)            # Plot deformed mesh
vis.plot_displacement_field('magnitude')  # Displacement contours
vis.plot_stress_field('von_mises')       # Stress contours

svg_export.py

SVGExporter: Export mesh and results to SVG format

from svg_export import SVGExporter

exporter = SVGExporter(solver)

# Export undeformed mesh
exporter.export_mesh('mesh_undeformed.svg', show_nodes=True)

# Export deformed mesh
coords_def = solver.get_deformed_coordinates(scale=50.0)
exporter.export_mesh('mesh_deformed.svg', coords=coords_def)

# Export boundary conditions with color coding
exporter.export_boundary_conditions('boundary_conditions.svg')

# Export side-by-side comparison
exporter.export_deformed_comparison('comparison.svg', scale=50.0)

validation.py

ValidationTests: Test suite with analytical solutions

from validation import ValidationTests, run_all_tests

# Run all validation tests
run_all_tests()

# Or run individual tests
tests = ValidationTests()
tests.test_constant_stress_field()      # Patch test
tests.test_mesh_convergence()           # Convergence study
tests.test_periodic_bc_enforcement()    # Verify periodic BC
tests.test_rigid_body_modes()           # Check RBM suppression
tests.test_manufactured_solution()      # Manufactured solution test

rve_generator.py

CompositeRVE: Generate RVE with circular fiber in matrix

from rve_generator import CompositeRVE

# Create RVE with circular fiber
rve = CompositeRVE(nx=40, ny=40, lx=1.0, ly=1.0,
                  fiber_radius=0.3,
                  E_matrix=3.5, nu_matrix=0.35,    # Epoxy matrix
                  E_fiber=230.0, nu_fiber=0.20)    # Carbon fiber

# Get element properties
E_elements, nu_elements = rve.get_all_element_properties()

# Visualize material distribution
rve.plot_material_distribution()

homogenization.py

PeriodicHomogenization: Compute effective properties via homogenization

from homogenization import PeriodicHomogenization

# Create homogenization solver
homog = PeriodicHomogenization(rve, thickness=1.0)

# Compute effective stiffness matrix
C_eff = homog.compute_effective_stiffness()

# Get effective properties
props = homog.effective_properties
print(f"Effective E_x: {props['E_x']:.2f} GPa")
print(f"Effective E_y: {props['E_y']:.2f} GPa")
print(f"Effective ν_xy: {props['nu_xy']:.4f}")
print(f"Effective G_xy: {props['G_xy']:.2f} GPa")

# Compare with theoretical bounds
homog.compare_with_theory()

Theory

Element Formulation

The Q4 element uses bilinear shape functions in natural coordinates (ξ, η):

N_i(ξ, η) = 1/4 (1 ± ξ)(1 ± η)

The strain-displacement relation is:

ε = B u_e

where B is derived from shape function derivatives.

Material Model

Plane stress constitutive matrix:

D = E/(1-ν²) [ 1      ν      0        ]
              [ ν      1      0        ]
              [ 0      0   (1-ν)/2    ]

Constraint Enforcement

Lagrange Multipliers: The augmented system is:

[ K   C^T ] [ u ]   [ F ]
[ C    0  ] [ λ ] = [ 0 ]

where C represents constraint equations and λ are Lagrange multipliers.

Penalty Method: Adds large penalty terms to enforce constraints:

K_ii += α  (for constrained DOF i)

where α is a large penalty parameter (typically 10^10).

Computational Homogenization

For periodic RVE, effective properties are computed by:

1. Applying prescribed strains: Three independent unit strain states:

   • ε̄ = [1, 0, 0]ᵀ (uniaxial x-strain)
   • ε̄ = [0, 1, 0]ᵀ (uniaxial y-strain)
   • ε̄ = [0, 0, 1]ᵀ (shear strain)

2. Solving FEM problem: For each prescribed strain, solve with periodic BC

3. Computing average stress: Volume average over RVE

   ⟨σ⟩ = (1/V) ∫_V σ dV
   

4. Building effective stiffness:

   C_eff relates ⟨σ⟩ = C_eff : ε̄
   

The effective stiffness matrix columns are the average stresses from the three load cases.

Theoretical Bounds:

• Voigt (upper bound): Assumes uniform strain throughout RVE

  E_Voigt = V_f E_f + V_m E_m
  

• Reuss (lower bound): Assumes uniform stress throughout RVE

  E_Reuss = 1 / (V_f/E_f + V_m/E_m)
  

FEM results should lie between these bounds, typically closer to Reuss for transverse loading of fiber composites.

Performance

Simple FEM (10×10 mesh):

• Assembly time: ~0.01 seconds
• Solution time: ~0.02 seconds
• Total DOFs: 242 (with constraints)

RVE Homogenization (40×40 mesh, composite):

• Assembly time: ~0.15 seconds
• Solution time per load case: ~0.05 seconds
• Total DOFs: ~3600
• Three load cases: ~0.15 seconds total solution time

Scales well for meshes up to 100×100 elements on standard hardware. For RVE problems, 40×40 to 60×60 meshes typically provide good accuracy.

Limitations

• Currently limited to square domains with structured meshes
• Plane stress formulation only (can be extended to plane strain)
• Linear elastic materials only
• Small deformation assumption

Future Extensions

Potential enhancements:

• Unstructured mesh support
• Nonlinear materials (plasticity, hyperelasticity)
• Geometric nonlinearity (large deformations)
• 3D extension
• Adaptive mesh refinement
• Parallel assembly and solution

References

Finite Element Method:

1. Zienkiewicz, O.C. and Taylor, R.L., "The Finite Element Method", Volume 1 & 2, 5th Edition
2. Hughes, T.J.R., "The Finite Element Method: Linear Static and Dynamic Finite Element Analysis"
3. Bathe, K.J., "Finite Element Procedures"

Computational Homogenization: 4. Suquet, P.M., "Elements of Homogenization for Inelastic Solid Mechanics", Homogenization Techniques for Composite Media, 1987 5. Geers, M.G.D. et al., "Multi-scale computational homogenization: Trends and challenges", Journal of Computational and Applied Mathematics, 2010 6. Nguyen, V.P. et al., "Homogenization-based multiscale crack modelling: from micro-diffusive damage to macro-cracks", Computer Methods in Applied Mechanics and Engineering, 2011

License

MIT License - feel free to use and modify for your purposes.

Author

Created as a demonstration of FEM implementation with periodic boundary conditions and rigid body constraints.

Contributing

Contributions are welcome! Please feel free to submit pull requests or open issues for bugs and feature requests.