1 #include "xlife++.h"
2 using namespace xlifepp;
3
4 Real omg=1, eps=1, mu=1, a=pi_, ome=omg* omg* mu* eps;
5
f(const Point & P,Parameters & pa=defaultParameters)6 Vector<Real> f(const Point& P, Parameters& pa = defaultParameters)
7 {
8 Real x=P(1), y=P(2);
9 Vector<Real> res(2);
10 Real c=2*a*a-ome;
11 res(1)=-c*cos(a*x)*sin(a*y);
12 res(2)= c*sin(a*x)*cos(a*y);
13 return res;
14 }
15
solex(const Point & P,Parameters & pa=defaultParameters)16 Vector<Real> solex(const Point& P, Parameters& pa = defaultParameters)
17 {
18 Real x=P(1), y=P(2);
19 Vector<Real> res(2);
20 res(1)=-cos(a*x)*sin(a*y);
21 res(2)= sin(a*x)*cos(a*y);
22 return res;
23 }
24
main(int argc,char ** argv)25 int main(int argc, char** argv)
26 {
27 init(_lang=en);
28 //mesh square using gmsh
29 Mesh mesh2d(Rectangle(_xmin=0, _xmax=1, _ymin=0, _ymax=1, _nnodes=50, _side_names="Gamma"), triangle, 1, gmsh);
30 Domain omega=mesh2d.domain("Omega");
31 Domain gamma=mesh2d.domain("Gamma");
32 //define space and unknown
33 Space V(omega, Nedelec, 1, "V");
34 Unknown e(V, "E");
35 TestFunction q(e, "q");
36 //define forms, matrices and vectors
37 BilinearForm aev=intg(omega, curl(e)|curl(q)) - ome*intg(omega, e|q);
38 LinearForm l=intg(omega, f|q);
39 EssentialConditions ecs = (ncross(e)|gamma=0);
40 //compute
41 TermMatrix A(aev, ecs, "A");
42 TermVector b(l, "B");
43 //solve
44 TermVector E=directSolve(A, b);
45 // P1 interpolation, L2 projection on H1
46 Space W(omega, P1, "W");
47 TermVector EP1=projection(E, W, 2);
48 EP1.name("E");
49 saveToFile("E", EP1, vtu);
50 return 0;
51 }
52