quadratic trial functions for a 2d FEM calculation

Question

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?

Answer 1

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.

Here's a link to the real code: https://github.com/tjolsen/Mesh_Utilities/blob/master/linear_to_quad_tris/linear_to_quadratic.m

function [QuadTris,Points] = linear_to_quadratic(t,p)

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

Answer 2

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.

Answer 3

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