XLiFE++ website is now hosted here with up-to-date content. The current website, with potentially deprecated content, will be removed soon.

Example 17/18 - 2D Bilaplacian problem

The 2d bilaplacian problem illustrates how to use Morley or Argyris element in XLiFE++:
{ Δ2u = f  in Ω u = 0 and u.n = 0 on Ω

The variational formulation in V = {v H2(Ω),u = 0 and u.n = 0 on Ω } is:

| find u V  such that  Ω (xxuxxv + yyuyyv + 2xyuxyv ) =Ωfvv V.

The implementation in XLiFE++ using Morley elements is the following :

 
#include "xlife++.h" 
using namespace xlifepp; 
 
// data function 
Real uex(const Point& P, Parameters& pa = defaultParameters) 
{ 
  Real x=P(1), y=P(2), kp=pi_; 
  Real r=sin(kp*x)*sin(kp*y); 
  return r*r; 
} 
 
Real f(const Point& P, Parameters& pa = defaultParameters) 
{ 
  Real x=P(1), y=P(2); 
  Real dkp=2*k*pi_; 
  Real cx=cos(dkp*x), cy=cos(dkp*y); 
  return 0.25*dkp*dkp*dkp*dkp*(4*cx*cy-cx-cy); 
} 
 
int main(int argc, char** argv) 
{ 
  init(argc, argv, _lang=en); // mandatory initialization of xlifepp 
// mesh rectangle 
  SquareGeo sq(_xmin=0, _xmax=1, _ymin=0, _ymax=1., _hsteps=0.05, _domain_name="Omega", _side_names="Gamma"); 
  Mesh mesh(sq, triangle, 1, gmsh); 
  Domain omega=mesh.domain("Omega"), gamma=mesh.domain("Gamma"); 
// create space 
  Space V(_domain=omega, _FE_type=Morley, _order=1, _name="V"); 
  Unknown u(V, "u");  TestFunction v(u, "v"); 
// create problem 
  EssentialConditions ecs= (u|gamma = 0) & (ndotgrad(u)|gamma = 0); 
  TermMatrix A(intg(omega, dxx(u)*dxx(v))+intg(omega, dyy(u)*dyy(v))+2*intg(omega, dxy(u)*dxy(v)), ecs); 
  TermVector B(intg(omega, f*v), "B"); 
// solve problem 
  TermVector U=directSolve(A, B); 
// interpolate on P1 and export to vtu file 
  Space W(_domain=omega, _interpolation=P1, _name="W", _notOptimizeNumbering); 
  Unknown w(W, "w"); 
  TermVector Up1=interpolate(w, omega, U, "Up1"); 
  saveToFile("U", Up1, vtu); 
  TermVector Ue1(w, omega, uex, "Ue1"); 
  TermVector Er=Up1-Ue1; 
  saveToFile("Er", Er, vtu); 
  return 0; 
}

pict

Figure 13.1: Approximate solution and difference with exact solution of the bilaplacian 2D problem

In order to use Argyris element in place of Morley element, replace the space definition:

Space V(_domain=omega, _FE_type=Argyris, _order=1, _name="V");

The next figure shows the L2 convergence rate of the Morley approximation when solving the 2d bilaplacian problem; in agreement with the order 2 predicting by the theory.


pict

Figure 13.2: L2 convergence rate of the Morley approximation