import sympy as sp
from IPython.display import Math, display
from symbolic_fem_workbench import (
ReferenceTriangleP1, AffineTriangleMap2D,
local_trial_expansion, local_test_expansion,
grad_2d, pullback_gradient_2d,
extract_coefficient_matrix, extract_coefficient_vector,
build_poisson_triangle_p1_local_problem,
)
from symbolic_fem_workbench.quadrature import integrate_reference_triangle_exact7 2D Poisson on a P1 Triangle: Gradient Pullback and Local Tensors
This notebook extends finite element analysis to 2D by deriving local element matrices for a Poisson equation on a P1 triangle. The key new concepts are the Jacobian-based gradient transformation (pullback) and handling 2D integrands over distorted elements.
7.1 Step 1: Reference triangle and shape functions
The reference triangle is defined on the unit triangle in the \(\xi\)-\(\eta\) parameter space: \[\hat{\Omega} = \{(\xi, \eta) : 0 \le \xi, \eta, \xi + \eta \le 1\}\]
Linear shape functions (P1) are barycentric coordinates.
xi, eta = sp.symbols("xi eta")
ref = ReferenceTriangleP1(xi, eta)
print("Reference nodes:", ref.nodes)
for i, N in enumerate(ref.shape_functions):
display(Math(f"N_{{{i+1}}} = " + sp.latex(N)))
print("\nReference gradients (constant for P1):")
for i, grad in enumerate(ref.shape_gradients_reference):
display(Math(r"\nabla_{\hat{}} N_" + str(i+1) + " = " + sp.latex(grad)))Reference nodes: ((0, 0), (1, 0), (0, 1))\(\displaystyle N_{1} = - \eta - \xi + 1\)
\(\displaystyle N_{2} = \xi\)
\(\displaystyle N_{3} = \eta\)
Reference gradients (constant for P1):\(\displaystyle \nabla_{\hat{}} N_1 = \left[\begin{matrix}-1\\-1\end{matrix}\right]\)
\(\displaystyle \nabla_{\hat{}} N_2 = \left[\begin{matrix}1\\0\end{matrix}\right]\)
\(\displaystyle \nabla_{\hat{}} N_3 = \left[\begin{matrix}0\\1\end{matrix}\right]\)
7.2 Step 2: Affine mapping to physical triangle
The affine mapping transforms the reference triangle to a physical triangle with vertices \((x_1, y_1)\), \((x_2, y_2)\), \((x_3, y_3)\): \[\mathbf{x}(\xi, \eta) = \mathbf{x}_1 + (\mathbf{x}_2 - \mathbf{x}_1)\xi + (\mathbf{x}_3 - \mathbf{x}_1)\eta\]
The Jacobian matrix is constant and relates derivatives in the two spaces.
x1, y1, x2, y2, x3, y3 = sp.symbols("x1 y1 x2 y2 x3 y3", real=True)
geom = AffineTriangleMap2D(xi=xi, eta=eta, x1=x1, y1=y1, x2=x2, y2=y2, x3=x3, y3=y3)
print("Jacobian matrix J:")
display(Math(r"J = " + sp.latex(geom.jacobian)))
print("\nDeterminant of J:")
display(Math(r"|J| = " + sp.latex(geom.det_jacobian)))
print("\nTriangle area = |J|/2 =", geom.area)
print("\nInverse transpose J^{-T} (for gradient pullback):")
display(Math(r"J^{-T} = " + sp.latex(geom.inverse_transpose_jacobian)))Jacobian matrix J:\(\displaystyle J = \left[\begin{matrix}- x_{1} + x_{2} & - x_{1} + x_{3}\\- y_{1} + y_{2} & - y_{1} + y_{3}\end{matrix}\right]\)
Determinant of J:\(\displaystyle |J| = x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\)
Triangle area = |J|/2 = x1*y2/2 - x1*y3/2 - x2*y1/2 + x2*y3/2 + x3*y1/2 - x3*y2/2
Inverse transpose J^{-T} (for gradient pullback):\(\displaystyle J^{-T} = \left[\begin{matrix}\frac{- y_{1} + y_{3}}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}} & \frac{y_{1} - y_{2}}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}}\\\frac{x_{1} - x_{3}}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}} & \frac{- x_{1} + x_{2}}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}}\end{matrix}\right]\)
7.3 Step 3: Gradient pullback
The key transformation: gradients in physical space relate to reference gradients via the inverse Jacobian transpose: \[\nabla u = J^{-T} \nabla_{\hat{}} \hat{u}\]
This ensures the chain rule is satisfied correctly for distorted elements.
d0, d1, d2 = sp.symbols("d0 d1 d2", real=True)
w0, w1, w2 = sp.symbols("w0 w1 w2", real=True)
uh = local_trial_expansion(ref.shape_functions, [d0, d1, d2])
vh = local_test_expansion(ref.shape_functions, [w0, w1, w2])
grad_uh_ref = grad_2d(uh, xi, eta)
grad_uh_phys = pullback_gradient_2d(grad_uh_ref, geom.jacobian)
print("Gradient of u_h in reference coords:")
display(Math(sp.latex(grad_uh_ref)))
print("\nGradient of u_h in physical coords (after J^{-T} pullback):")
display(Math(sp.latex(grad_uh_phys)))Gradient of u_h in reference coords:\(\displaystyle \left[\begin{matrix}- d_{0} + d_{1}\\- d_{0} + d_{2}\end{matrix}\right]\)
Gradient of u_h in physical coords (after J^{-T} pullback):\(\displaystyle \left[\begin{matrix}\frac{\left(d_{0} - d_{1}\right) \left(y_{1} - y_{3}\right) - \left(d_{0} - d_{2}\right) \left(y_{1} - y_{2}\right)}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}}\\\frac{- \left(d_{0} - d_{1}\right) \left(x_{1} - x_{3}\right) + \left(d_{0} - d_{2}\right) \left(x_{1} - x_{2}\right)}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}}\end{matrix}\right]\)
7.4 Step 4: Form the integrand and extract local tensors
The weak form of the Poisson equation is: \[a(u, v) = \int_\Omega \nabla u \cdot \nabla v \, dx\] \[F(v) = \int_\Omega f v \, dx\]
We evaluate the integrals on the reference element and extract coefficients.
grad_vh_ref = grad_2d(vh, xi, eta)
grad_vh_phys = pullback_gradient_2d(grad_vh_ref, geom.jacobian)
bilinear_integrand = sp.simplify((grad_uh_phys.T * grad_vh_phys)[0] * geom.det_jacobian)
f = sp.Symbol("f", real=True)
linear_integrand = sp.simplify(f * vh * geom.det_jacobian)
bilinear_val = integrate_reference_triangle_exact(bilinear_integrand, xi, eta)
linear_val = integrate_reference_triangle_exact(linear_integrand, xi, eta)
Ke = extract_coefficient_matrix(bilinear_val, [d0, d1, d2], [w0, w1, w2])
fe = extract_coefficient_vector(linear_val, [w0, w1, w2])
print("Element stiffness matrix Ke (generic triangle):")
display(Math(r"K_e = " + sp.latex(sp.simplify(Ke))))
print("\nElement load vector fe (constant source f):")
display(Math(r"f_e = " + sp.latex(sp.simplify(fe))))Element stiffness matrix Ke (generic triangle):\(\displaystyle K_e = \left[\begin{matrix}\frac{\frac{x_{2}^{2}}{2} - x_{2} x_{3} + \frac{x_{3}^{2}}{2} + \frac{y_{2}^{2}}{2} - y_{2} y_{3} + \frac{y_{3}^{2}}{2}}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}} & \frac{- x_{1} x_{2} + x_{1} x_{3} + x_{2} x_{3} - x_{3}^{2} - y_{1} y_{2} + y_{1} y_{3} + y_{2} y_{3} - y_{3}^{2}}{2 \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)} & \frac{x_{1} x_{2} - x_{1} x_{3} - x_{2}^{2} + x_{2} x_{3} + y_{1} y_{2} - y_{1} y_{3} - y_{2}^{2} + y_{2} y_{3}}{2 \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)}\\\frac{- x_{1} x_{2} + x_{1} x_{3} + x_{2} x_{3} - x_{3}^{2} - y_{1} y_{2} + y_{1} y_{3} + y_{2} y_{3} - y_{3}^{2}}{2 \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)} & \frac{\frac{x_{1}^{2}}{2} - x_{1} x_{3} + \frac{x_{3}^{2}}{2} + \frac{y_{1}^{2}}{2} - y_{1} y_{3} + \frac{y_{3}^{2}}{2}}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}} & \frac{- x_{1}^{2} + x_{1} x_{2} + x_{1} x_{3} - x_{2} x_{3} - y_{1}^{2} + y_{1} y_{2} + y_{1} y_{3} - y_{2} y_{3}}{2 \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)}\\\frac{x_{1} x_{2} - x_{1} x_{3} - x_{2}^{2} + x_{2} x_{3} + y_{1} y_{2} - y_{1} y_{3} - y_{2}^{2} + y_{2} y_{3}}{2 \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)} & \frac{- x_{1}^{2} + x_{1} x_{2} + x_{1} x_{3} - x_{2} x_{3} - y_{1}^{2} + y_{1} y_{2} + y_{1} y_{3} - y_{2} y_{3}}{2 \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)} & \frac{\frac{x_{1}^{2}}{2} - x_{1} x_{2} + \frac{x_{2}^{2}}{2} + \frac{y_{1}^{2}}{2} - y_{1} y_{2} + \frac{y_{2}^{2}}{2}}{x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}}\end{matrix}\right]\)
Element load vector fe (constant source f):\(\displaystyle f_e = \left[\begin{matrix}\frac{f \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)}{6}\\\frac{f \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)}{6}\\\frac{f \left(x_{1} y_{2} - x_{1} y_{3} - x_{2} y_{1} + x_{2} y_{3} + x_{3} y_{1} - x_{3} y_{2}\right)}{6}\end{matrix}\right]\)
7.5 Step 5: Specialize to the unit right triangle
For the special case of a unit right triangle with vertices at \((0,0)\), \((1,0)\), \((0,1)\), the Jacobian is the identity, and formulas simplify.
unit_subs = {x1: 0, y1: 0, x2: 1, y2: 0, x3: 0, y3: 1}
print("Unit right triangle Ke:")
display(Math(r"K_e = " + sp.latex(sp.simplify(Ke.subs(unit_subs)))))
print("\nUnit right triangle fe:")
display(Math(r"f_e = " + sp.latex(sp.simplify(fe.subs(unit_subs)))))Unit right triangle Ke:\(\displaystyle K_e = \left[\begin{matrix}1 & - \frac{1}{2} & - \frac{1}{2}\\- \frac{1}{2} & \frac{1}{2} & 0\\- \frac{1}{2} & 0 & \frac{1}{2}\end{matrix}\right]\)
Unit right triangle fe:\(\displaystyle f_e = \left[\begin{matrix}\frac{f}{6}\\\frac{f}{6}\\\frac{f}{6}\end{matrix}\right]\)
7.6 Verify with pre-built workflow
data = build_poisson_triangle_p1_local_problem()
print("Pre-built Ke (unit right triangle):")
display(Math(sp.latex(data["Ke_unit_right_triangle"])))
print("Match:", sp.simplify(Ke.subs(unit_subs) - data["Ke_unit_right_triangle"]) == sp.zeros(3))Pre-built Ke (unit right triangle):\(\displaystyle \left[\begin{matrix}1 & - \frac{1}{2} & - \frac{1}{2}\\- \frac{1}{2} & \frac{1}{2} & 0\\- \frac{1}{2} & 0 & \frac{1}{2}\end{matrix}\right]\)
Match: True7.7 Exercise
Question: Evaluate \(K_e\) for an equilateral triangle with vertices at \((0,0)\), \((1,0)\), \((1/2, \sqrt{3}/2)\). Is the matrix still symmetric? Compare the diagonal entries with the unit right triangle case. Why do they differ?
# Exercise: Evaluate Ke for an equilateral triangle with vertices
# (0,0), (1,0), (0.5, sqrt(3)/2). Is the matrix still symmetric?
# Compare the diagonal entries with the unit right triangle case.