Coverage for pyxfem/core/nbkalgebra2d.py: 0%
69 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*
3# pylint: disable=invalid-name
4r"""
5Created march 2026
6@author: nchevaug
82D Algebraic Kernel for High-Performance Numba JIT Compilation.
10Philosophy & Architecture:
11This module provides fundamental 2D vector and symmetric tensor operations
12specifically tailored for Numba's `nopython` mode (`@njit`).
14To achieve maximum performance and avoid the massive overhead of heap
15allocations inside hot loops (e.g., cell or vertex iterations), this module
16STRICTLY avoids NumPy arrays for small local variables. Instead, it relies
17entirely on Python `Tuple` of floats.
19Why Tuples?
20When compiled with Numba and LLVM, homogenous tuples of primitive types
21(like `Tuple[float, float]`) are not treated as standard Python objects.
22They are lowered directly into C-style structs passed by value, which the
23LLVM backend maps directly into the processor's hardware registers.
24This guarantees zero-allocation overhead, completely bypasses the Python
25garbage collector, and allows aggressive compiler optimizations like loop
26unrolling and SIMD vectorization (enhanced by `fastmath=True`).
28Storage Convention:
29- Vectors (Vec2) are stored as 2-tuples: `(x, y)`
30- Symmetric 2x2 Matrices (Sy22) are stored as 3-tuples: `(A00, A01, A11)`.
31 This exploits the symmetry ($A_{01} = A_{10}$) to save register space
32 and reduce floating-point operations.
33"""
35from typing import Tuple
36from pyxfem.core.config import USE_NUMBA
38# ==========================================
39# STATIC ANALYSIS TYPE ALIASES (Type Hinting)
40# ==========================================
41# These types guide the IDE (VS Code, Pylance, MyPy) for autocompletion
42# and static type checking. They are ignored by Numba at runtime.
43# Vec2 = Tuple[float, float]
44# Sy22 = Tuple[float, float, float]
45#
47if USE_NUMBA:
48 from math import sqrt
49 import numba
50 import numba.core.types as nct
51 # pylint: disable=ungrouped-imports
52 from pyxfem.core.config import FASTMATH, nbfloat, fzero,npfloat
53 # ==========================================
54 # EXPLICIT NUMBA SIGNATURE TYPES
55 # ==========================================
56 sy22_t = nct.UniTuple(nbfloat, 3)
57 vec2_t = nct.UniTuple(nbfloat, 2)
58 # Vec2 = vec2_t
59 # Sy22 = sy22_t
60 # Vec2 = Tuple[nbfloat, nbfloat]
61 # Sy22 = Tuple[nbfloat, nbfloat, nbfloat]
62 Vec2 = Tuple[npfloat, npfloat]
63 Sy22 = Tuple[npfloat, npfloat, npfloat]
64 Ge22 = Tuple[npfloat, npfloat, npfloat, npfloat]
65 # ==========================================
66 # constants
67 # ==========================================
68 vec2_zero: Vec2 = (fzero, fzero)
69 sy22_zero: Sy22 = (fzero, fzero, fzero)
71 # ==========================================
72 # 2D VECTOR OPERATIONS
73 # ==========================================
74 @numba.njit(inline='always', fastmath = FASTMATH)
75 def set_vec2(index, array, v):
76 array[index, 0], array[index, 1] = v[0], v[1]
78 @numba.njit(cache=True, fastmath = FASTMATH)
79 def vec2_add(a: Vec2, b: Vec2) -> Vec2:
80 r"""
81 Adds two 2D vectors.
82 Math: $a + b$
83 """
84 return (a[0] + b[0], a[1] + b[1])
86 @numba.njit(cache=True, fastmath = FASTMATH)
87 def vec2_add3(a: Vec2, b: Vec2, c: Vec2) -> Vec2:
88 r"""
89 Adds three 2D vectors.
90 Math: $a + b + c$
91 """
92 return (a[0] + b[0]+c[0], a[1] + b[1]+c[1])
94 @numba.njit(cache=True, fastmath = FASTMATH)
95 def vec2_scal(alpha: npfloat, x: Vec2) -> Vec2:
96 r"""
97 Multiplies a 2D vector by a scalar.
98 Math: $\alpha x$
99 """
100 return (alpha*x[0], alpha*x[1])
101 @numba.njit(cache=True, fastmath = FASTMATH)
103 def vec2_sub(a: Vec2, b: Vec2) -> Vec2:
104 r"""
105 Subtracts one 2D vector from another.
106 Math: $a - b$
107 """
108 return (a[0] - b[0], a[1] - b[1])
110 @numba.njit(cache=True, fastmath = FASTMATH)
111 def vec2_norm(a: Vec2) -> npfloat:
112 r"""
113 Computes the Euclidean (L2) norm of a 2D vector.
114 Math: $\|a\| = \sqrt{a_0^2 + a_1^2}$
115 """
116 return npfloat(sqrt((a[0]**2 +a[1]**2)))
118 @numba.njit(cache=True, fastmath = FASTMATH)
119 def vec2_dot(a: Vec2, b: Vec2) -> npfloat:
120 r"""
121 Computes the dot product of two 2D vectors.
122 Math: $a \cdot b = a_0 b_0 + a_1 b_1$
123 """
124 return a[0]*b[0] +a[1]*b[1]
126 @numba.njit(cache=True, fastmath = FASTMATH)
127 def vec2_square(a: Vec2) -> npfloat:
128 r"""
129 Computes the squared norm of a 2D vector.
130 Math: $\|a\|^2 = a_0^2 + a_1^2$
131 """
132 return a[0]**2 +a[1]**2
134 # ==========================================
135 # 2X2 SYMMETRIC TENSOR OPERATIONS
136 # ==========================================
138 @numba.njit(cache=True, fastmath = FASTMATH)
139 def sy22_mv(A: Sy22, x: Vec2) -> Vec2:
140 r"""
141 Performs a matrix-vector multiplication for a symmetric 2x2 matrix.
142 Math: $A x$
143 Note: The matrix is implicitly reconstructed from its components $(A_{00}, A_{01}, A_{11})$.
144 """
145 A00, A01, A11 = A[0], A[1], A[2]
146 x0, x1 = x[0], x[1]
147 return (A00*x0 + A01*x1, A01*x0+ A11*x1)
149 @numba.njit(cache=True, fastmath = FASTMATH)
150 def sy22_add(A: Sy22, B: Sy22) -> Sy22:
151 r"""
152 Adds two symmetric 2x2 matrices.
153 Math: $A + B$
154 """
155 return (A[0]+B[0], A[1]+B[1], A[2]+B[2])
157 @numba.njit(cache=True, fastmath = FASTMATH)
158 def sy22_r(alpha: npfloat, x: Vec2) -> Sy22:
159 r"""
160 Computes the symmetric rank-1 update (dyadic product).
161 Math: $\alpha (x \otimes x)$
162 """
163 A00 = alpha * x[0] * x[0]
164 A01 = alpha * x[0] * x[1]
165 A11 = alpha * x[1] * x[1]
166 return (A00, A01, A11)
168 @numba.njit(cache=True, fastmath = FASTMATH)
169 def sy22_r2(alpha: npfloat, x: Vec2, y: Vec2) -> Sy22:
170 r"""
171 Computes the scalar-penalized symmetric rank-2 update.
172 Math: $\alpha (x \otimes y + y \otimes x)$
173 """
174 A00 = alpha*2*x[0] * y[0]
175 A01 = alpha*(x[0] * y[1] + y[0] * x[1])
176 A11 = alpha*2*(x[1] * y[1])
177 return (A00, A01, A11)
178 @numba.njit(cache=True, fastmath = FASTMATH)
180 def sy22_trace(A: Sy22) -> npfloat:
181 r"""
182 Computes the trace of a symmetric 2x2 matrix.
183 Math: $\text{Tr}(A) = A_{00} + A_{11}$
184 """
185 return A[0] + A[2]
187 @numba.njit(cache=True, fastmath = FASTMATH)
188 def sy22_smv(A: Sy22, b: Vec2) -> Vec2:
189 r"""
190 Solves the linear system $A x = b$ for a symmetric 2x2 matrix
191 using explicit Cramer's rule inversion.
192 Math: $x = A^{-1} b$
193 """
194 A00, A01, A11 = A[0], A[1], A[2]
195 b0, b1 = b[0], b[1]
196 detA = A00*A11 - A01*A01
197 x0 = (A11*b0 - A01*b1)/detA
198 x1 = (A00*b1 - A01*b0)/detA
199 return (x0, x1)
201 # # The two following function where not used yet :need to b checked
202 # # for efficiency
203 # @numba.njit(cache=True)
204 # def gemm22_numba(alpha, A, B, beta, C):
205 # C[0,0] = beta*C[0,0] + alpha*(A[0,0]*B[0,0] + A[0,1]*B[1,0])
206 # C[0,1] = beta*C[0,1] + alpha*(A[0,0]*B[0,1] + A[0,1]*B[1,1])
207 # C[1,0] = beta*C[1,0] + alpha*(A[1,0]*B[0,0] + A[1,1]*B[1,0])
208 # C[1,1] = beta*C[1,1] + alpha*(A[1,0]*B[0,1] + A[1,1]*B[1,1])
209 # @numba.njit(cache=True)
210 # def gemv22_numba(A,x,y, beta =fzero):
211 # x0, x1= x[0], x[1]
212 # y[0] = beta*y[0] + A[0,0]*x0 + A[0,1]*x1
213 # y[1] = beta*y[1] + A[1,0]*x0 + A[1,1]*x1