1
$\begingroup$

I have output from a 3D linear elasticity finite element simulation which uses linear tetrahedral elements, such that the displacement is continuous over the nodes but the gradient is not ($C_0$ continuous). I wish to "upgrade" the displacement data so that it can be reported as a $C_1$ continuous quantity which I can directly differentiate to produce the gradient, which would yield a consistency that I find appealing. Instead of directly interpolating gradient data over the nodes (or using a higher order basis function, which of course would make this exercise moot), I wish to interpolate the primary variable, displacement, as a quadratic polynomial with splines.

For 1D quadratic splines, each element or segment takes the form: $$f_e(x)=a_e+b_ex+c_ex^2$$ $$\nabla f_e(x)=b_e+2 c_ex$$ To form the system of equations, the function value and gradient are matched at end points. And since this is 1D, each end point connects to at most one other segment. $$f_e(x_i)=f_{e+1}(x_i)$$ and $$\nabla f_e(x_i)=\nabla f_{e+1}(x_i)$$ Such is the predictable nature of this arrangement that you are guaranteed to have $3(n-1)$ variables and at least $n+2(n-1)=3n-2$ equations. This can be supplemented with an additional assumption or solved in a linear least squares fashion. Awesome.

For the 3D case, it looks similar: $$f_e(x)=a_e+x^Tb_e+x^Tc_ex$$ and $$\nabla f_e(x)=b_e+2c_ex$$where $x$ is a $3$ x $1$ vector, $a_e$ is a scalar, $b_e$ is a $3$ x $1$ vector, and $c_e$ is a $3$ x $3$ matrix. Because $x^Tc_ex$ is a quadratic form, only the symmetric part of $c_e$ is of any consequence, so it can be safely assumed $c_e$ is a symmetric matrix with 6 independent variables:

$$f_e(x)=a_e+x^Tb_e+\begin{pmatrix}x_0x_0 & x_0x_1 & x_0x_2 & x_1x_1 & x_1x_2 & x_2x_2\end{pmatrix}^T\begin{pmatrix}c_{e00} & c_{e01} & c_{e02} & c_{e11} & c_{e12} & c_{e22}\end{pmatrix}$$ $$\nabla f_e(x)=b_e+\begin{pmatrix}2x_0 & x_1 & x_2 & 0 & 0 & 0\\\ 0 & x0 & 0 & 2x_1 & x_2 & 0\\\ 0 & 0 & x_0 & 0 & x_1 & 2x_2\end{pmatrix}^T\begin{pmatrix}c_{e00} & c_{e01} & c_{e02} & c_{e11} & c_{e12} & c_{e22}\end{pmatrix}$$

The hard part is that each node is not limited to 2 neighboring elements, in some cases it can be upward of 20! In one very modest case (I usually work with much finer grids), I have approximately 2000 tetrahedra and 500 nodes, with 51,000 (!) combinations where I must enforce gradient continuity. For example: if two tetrahedra (T1, T2) share 3 nodes (A, B, C), I must enforce gradient continuity for at A, B and C using the spline variables for T1 and T2.

This yields a sparse matrix ($10^{-4}$ sparsity) with $4t+3*51,000 = 161,000$ equations and $10t=20,000$ variables, and a $161,000$ x $3$ right-hand side (each component of displacement is independently interpolated). This is obviously way more complicated than calculating the original finite element system, which seems counterintuitive to me. Nevertheless, it does work...slowly. In the case above, I solved the original finite element system in 500 milliseconds and the interpolation took 30 seconds. I didn't even have the memory needed (it blew up to 120 gigs or something silly like that) to interpolate my other grids which takes a whopping 4 seconds to solve per time step. I am using SuiteSparse's SPQR solver, which produces a least squares solution using a QR factorization of the giant sparse system.

Questions:

  1. What do you call this interpolation procedure? Splines? Super-splines? B-Splines? I can't make heads or tails of what I read in literature on anything beyond 1D splines.

  2. Am I "double counting" some conditions and making the system overly complicated? I thought I could ignore a few conditions, such as pairs of tetrahedra that touch at one point only, but ignoring even that destroys the interpolation. I am dumbfounded by how much more complicated this interpolation procedure is than the original finite element calculations.

  3. Any ideas on how to simplify this system? I am confident that the algorithms I am using are top notch -- I'm just giving them too much work.

Examples plotting pressure (trace of the gradient) without interpolation:

enter image description here

and with the procedure above:

enter image description here

Improve this question
$\endgroup$
3
  • $\begingroup$ Since it's so expensive, why not use a local "recovery" operation to obtain a piecewise linear and continuous gradient field? That is, instead of postprocessing the solution, just fine some kind of "weak" but continuous gradient of it. $\endgroup$
    – Wolfgang Bangerth
    Commented Mar 2, 2021 at 2:36
  • $\begingroup$ @WolfgangBangerth If I understand you correctly, that is like creating a interpolation of the gradient using the linear basis functions used to solve for displacement. I described my approach in this answer (which takes 10s-100s of milliseconds to compute): scicomp.stackexchange.com/questions/36852/… If creating splines for 3D data is known to be a futile process, then I am OK settling for the weak continuous gradient approach. $\endgroup$
    – Charlie S
    Commented Mar 2, 2021 at 3:38
  • $\begingroup$ S: Something like this. It's not an interpolation but a projection operation. The usual approach, maybe equivalent to what you do, is to determine the values of the "recovered" stress at each node point as the average of the gradients as taken from the cells adjacent to the vertex, weighted by the area of these cells. You can do this one vertex at a time, so it's all local operations and no global linear system needs to be solved. $\endgroup$
    – Wolfgang Bangerth
    Commented Mar 2, 2021 at 12:17

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.