quadratic trial functions for a 2d FEM calculation

I want to solve \begin{align} \nabla^4\psi+\alpha\nabla^2\psi+\beta\psi=F(x,y);\quad \nabla \psi\cdot \hat{n}=\nabla^3\psi\cdot \hat{n}=0\quad \text{on boundaries} \end{align} with a 2d FEM scheme. With the biharmonic operator, linear test functions are ruled out. I can integrate the first term by parts twice, so that the weak version of the PDE includes a term \begin{align} \int_A\nabla^2\phi_i \nabla^2\phi_j\,da \end{align} Quadratic test functions will allow me to evaluate the integral. I understand that quadratic test functions have nodes at the vertices and at the mid-point of each side of each triangle . My question relates to how to deal with this at assembly. If the coordinates of the vertices of each triangle are in p, and the connectivity matrix is t, then what is the procedure to extend these for the additional nodes located on the sides of each triangle?

This wonderful reference http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch23.d/AFEM.Ch23.pdf '23.2.7. Six Node Quadratic Interpolation' describes the quadratic elements perfectly, but I can't figure out how to integrate the midway points into my node list?

• Are you trying to figure out how to inject new nodes at the midpoints of the edges of a linear mesh? One way is to explicitly build up a list of mesh edges (with each edge holding the ID of the triangle on each side of it, at a minimum). Then, create new nodes in the middle of each of these edges, updating the element connectivity appropriately. If you're using Matlab (or similar), all of this can probably be done in just a few lines of code reasonably efficiently by making use of the "sort" and "unique" functions. – Tyler Olsen Sep 10 '17 at 4:53
• @TylerOlsen That is exactly what I am trying to do, I will try your suggestion. Thanks – Clinton Winant Sep 10 '17 at 16:10

Here's the version that I had in mind. It's a bit terse, but there is a commented version up on github (use that one, it actually defines the inputs/outputs and explains the logic for the various intermediate calculations). Based on some initial benchmarking, it beats out your loop/branch-heavy code by about a factor of 300 (Matlab 2017a on a large-ish mesh), since it spends all of its time in vectorized calls.

nt = size(t,1);
t_idx = (1:nt)';

all_edges = [t(:,1) t(:,2) t_idx t_idx 3*ones(nt,1);
t(:,2) t(:,3) t_idx t_idx 1*ones(nt,1);
t(:,3) t(:,1) t_idx t_idx 2*ones(nt,1)];

forward_edges = all_edges(:,2) > all_edges(:,1);
all_edges(forward_edges,4) = 0;
all_edges(~forward_edges,3) = 0;
all_edges(~forward_edges,[1,2]) = all_edges(~forward_edges,[2,1]);

[unique_edges,~,ic] = unique(all_edges(:,[1,2]), 'rows');
Edges = [unique_edges,
accumarray(ic, all_edges(:,3),[],@sum),
accumarray(ic, all_edges(:,4),[],@sum)];
t2f = accumarray([all_edges(:,3)+all_edges(:,4),all_edges(:,5)],
sign(all_edges(:,3) - all_edges(:,4)).*ic, [], @sum);

newPoints = [(p(Edges(:,1),1)+p(Edges(:,2),1))/2,
(p(Edges(:,1),2)+p(Edges(:,2),2))/2];
Points = [p; newPoints];

np = size(p,1);
QuadTris = [t, np+(abs(t2f))];

end
• I would have voted this up, but don't have enough brownie points to do that – Clinton Winant Sep 11 '17 at 16:23
• Glad it helped! Stick around, answer some questions, and those brownie points will come :) – Tyler Olsen Sep 11 '17 at 17:23

I suspect that you will be disappointed to find out that, once you are able to implement the quadratic shape functions, the solution you will get is not correct. That's because for the bilaplacian operator ($\nabla^4$) it is not enough to have shape functions that are quadratic polynomials on each cell, but you also need to have special kind of (Hermite) shape functions that have no "kinks" across cell interfaces but are continuously differentiable at faces. These are substantially more difficult to implement than the quadratic shape functions you are looking for, but necessary to produce the correct solution of the equation you are trying to solve.

• Good to know. Out of curiosity, do you have a reference that expands on this? I could imagine constructing for a tensor-product element out of the cubics typically used for 1D beams, but I wouldn't know how to approach the problem for simplex elements. – Tyler Olsen Sep 11 '17 at 19:45
• @WolfgangBangerth Following on Tyler's comment, Any reference to solving the equation I wrote in the OP would be very much appreciated. – Clinton Winant Sep 11 '17 at 19:58
• In structural mechanics, the pde with the bilaplacian operator represents bending of a thin plate. The set of notes you reference above show one reasonable way of implementing a finite element for this pde. The difficulty that @Wolfgangbangerth alludes to is one of the reasons this has been an active research area since the 60s. – Bill Greene Sep 11 '17 at 21:10
• Indeed. @ClintonWinant -- look through the typical finite element books and look for the term "biharmonic equation". You will find considerations that come down to the fact that after integration by parts, you still have to take second derivatives of shape functions -- but this doesn't work with the usual functions because they have a kink: their first derivative exists (but is discontinuous) and their second derivatives does not exist any more at cell interfaces. – Wolfgang Bangerth Sep 12 '17 at 1:58
• @WolfgangBangerth This seems consistent with my experiece with the oned version of the equation in my OP. With 1d cubic elements that make both the solution and the first derivative continuous, the FEM solution reproduces exact solutions to whatever precision one needs. Perhaps I would be better off dealing with the system of two equations in two unknowns that give rise to my fouth order problem in 2d? – Clinton Winant Sep 12 '17 at 2:09

Following @TylerOlson 's suggestion here is an Octave/Matlab script that generates a mesh with six points in each triangle: the three vertices and the three mid-points. If anyone knows a way to simplify this inelegant code, I would appreciate it.

clear
%% Persson's routine to make simple triangles in a rectangular grid
nx=4;ny=3;
[x,y]=ndgrid(linspace(-1,1,nx),linspace(-1,1,ny)); % forms x and y lists
p=[x(:),y(:)];
t=[1,nx+2,nx+1;1,2,nx+2];
t=kron(t,ones(nx-1,1))+kron(ones(size(t)),(0:nx-2)');
t=kron(t,ones(ny-1,1))+kron(ones(size(t)),(0:ny-2)'*nx);
np=length(p(:,1));nt=length(t(:,1));
%%make a list of all sides from t;  There will be redundant entries
s=zeros(3*nt,6);
for kt=1:nt
for ks=1:3
if ks==1
s(3*(kt-1)+ks,1)=min(t(kt,1:2));
s(3*(kt-1)+ks,2)=max(t(kt,1:2));
s(3*(kt-1)+ks,3)=kt;
elseif ks==2
s(3*(kt-1)+ks,1)=min(t(kt,2:3));
s(3*(kt-1)+ks,2)=max(t(kt,2:3));
s(3*(kt-1)+ks,3)=kt;
else
s(3*(kt-1)+ks,1)=min(t(kt,[1,3]));
s(3*(kt-1)+ks,2)=max(t(kt,[1,3]));
s(3*(kt-1)+ks,3)=kt;
endif
endfor
endfor
%% once this has been done for all triangles, work on s to eliminate dupes
for kk=1:max(s(:,1))
order1=find(s(:,1)==kk);%s(order1,:)
[junk,order2]=sort(s(order1,2));
dupe=find(junk(1:end-1)==junk(2:end));
n=s(order1(order2),:);
for kl=1:length(dupe)
n(dupe(end+1-kl),4)=n(dupe(end+1-kl)+1,3);% add on triangle no of dupe
endfor
[junk,order3,mojunk]=unique(n(:,2),"first");
n=n(order3,:);
if kk==1
snodupe=n;
else
snodupe=[snodupe;n];
endif
endfor
%% add to position table the position of new nodes
for kn=1:length(snodupe)
p(np+kn,1)=0.5*(p(snodupe(kn,1),1)+p(snodupe(kn,2),1));
p(np+kn,2)=0.5*(p(snodupe(kn,1),2)+p(snodupe(kn,2),2));
endfor
%% build a new connectivity table, t6
t6=zeros(nt,6);
for kt=1:nt
%t(kt,:)
for ks=1:3
if ks==1
ss=t(kt,1);
se=t(kt,2);
elseif ks==2
ss=t(kt,2);
se=t(kt,3);
else
ss=t(kt,3);
se=t(kt,1);
endif
%% where is this side in snodupe table?
%% to lookup in snodupe need ends to be in order
s1=min([ss se]);s2=max([ss se]);
ind1=find(snodupe(:,1)==s1);
ind2=find(snodupe(ind1,2)==s2);
sideno=ind1(1)-1+ind2;
t6(kt,2*ks-1)=ss;
t6(kt,2*ks)=np+sideno;
endfor
endfor