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++:The variational formulation in is:
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;
}
In order to use Argyris element in place of Morley element, replace the space definition:
The next figure shows the convergence rate of the Morley approximation when solving the 2d bilaplacian problem; in agreement with the order 2 predicting by the theory.