1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18
19 \file sys_FE_IM.cpp
20 \author E. Lunéville
21 \since 1 March 2014
22 \date 1 March 2014
23
24 Test of FEM-IEM method on 2D Helmholtz problem on a disk
25
26 */
27
28 #include "xlife++-libs.h"
29 #include "testUtils.hpp"
30 #include <exception>
31
32 using namespace xlifepp;
33
34 namespace sys_FE_IM
35 {
36
data_g(const Point & P,Parameters & pa=defaultParameters)37 Complex data_g(const Point& P, Parameters& pa = defaultParameters)
38 {
39 Real x=P(1), k=pa("k");
40 Vector<Complex> g(2,0.);
41 g(1) = i_*k*exp(i_*k*x);
42 return dot(g,P/norm2(P)); //dr(e^{ikx}
43 }
44
u_inc(const Point & P,Parameters & pa=defaultParameters)45 Complex u_inc(const Point& P, Parameters& pa = defaultParameters)
46 {
47 Real x=P(1), k=pa("k");
48 return exp(i_*k*x);
49 }
50
sys_FE_IM(bool check)51 String sys_FE_IM(bool check)
52 {
53 String rootname = "sys_FE_IM";
54 trace_p->push(rootname);
55 std::stringstream out;
56 out.precision(testPrec);
57 verboseLevel(100);
58
59 //numberOfThreads(8);
60
61 //parameters
62 Number nh = 20; // number of elements on Gamma
63 Real h=2.*pi_/nh; // size of mesh
64 Real re=1.+2.*h; // exterior radius
65 Number ne=Number(2*pi_*re/h); // number of elements on Sigma
66 Real l = 4.*re; // length of exterior square
67 Number nr=Number(4*l/h); // number of elements on exterior square
68 Real k= 4, k2=k*k; // wavenumber
69 Parameters pars;
70 pars << Parameter(k,"k") << Parameter(k2,"k2");
71 Kernel H=Helmholtz2dKernel(pars);
72 Function g(data_g,pars);
73 Function ui(u_inc,pars);
74 //Mesh and domains definition
75 Disk d1(_center=Point(0.,0.),_radius=1,_nnodes=nh,
76 _side_names=Strings(4,"Gamma"));
77 Disk d2(_center=Point(0.,0.),_radius=re,_nnodes=ne,_domain_name="Omega",
78 _side_names=Strings(4,"Sigma"));
79 Square rect(_center=Point(0.,0.),_length=l,_nnodes=nr,_domain_name="Omega_ext");
80 Mesh mesh(rect+(d2-d1),_triangle,1,_gmsh);
81 Domain omega=mesh.domain("Omega");
82 Domain sigma=mesh.domain("Sigma");
83 Domain gamma=mesh.domain("Gamma");
84 Domain omega_ext=mesh.domain("Omega_ext"); //for integral representation
85
86 //create P2 Lagrange interpolation
87 Space V(omega,P2,"V",true);
88 Unknown u(V,"u"); TestFunction v(u,"v");
89
90 theCout<<intg(sigma,grad_x(grad_y(H)|_ny)*v)<<eol;
91
92 // create bilinear form and linear form
93 Complex lambda=-i_*k;
94 BilinearForm auv =
95 intg(omega,grad(u)|grad(v))-k2*intg(omega,u*v)+lambda*intg(sigma,u*v)
96 +intg(sigma,gamma,u*(grad_y(grad_x(H)|_nx)|_ny)*v)
97 +lambda*intg(sigma,gamma,u*(grad_y(H)|_ny)*v);
98 BilinearForm alv =
99 intg(sigma,gamma,u*(grad_x(H)|_nx)*v)+lambda*intg(sigma,gamma,u*H*v);
100 TermMatrix ALV(alv), A(auv);
101 TermVector B(intg(gamma,g*v));
102 TermVector G(u,gamma,g);
103 B+=ALV*G;
104
105 //solve linear system AU=F
106 TermVector U=directSolve(A,B);
107 saveToFile("U.vtk",U,vtk);
108 // sigma.saveNormalsToFile("n2_sigma",_vtk);
109 // gamma.saveNormalsToFile("n2_gamma",_vtk);
110 //integral representation on omega_ext
111 Space Vext(omega_ext,P2,"Vext",false);
112 Unknown uext(Vext,"uext");
113 TermVector Uext =
114 -integralRepresentation(uext, omega_ext, intg(gamma,(grad_y(H)|_ny)*U))
115 +integralRepresentation(uext, omega_ext, intg(gamma,H*G));
116 saveToFile("Uext.vtk",Uext,vtk);
117 //total field
118 TermVector Ui(u, omega, ui), Utot=Ui+U;
119 TermVector Uiext(uext, omega_ext, ui), Utotext=Uiext+Uext;
120 saveToFile("Utot.vtk",Utot,vtk);
121 saveToFile("Utotext.vtk",Utotext,vtk);
122
123 //------------------------------------------------------------------------------------
124 // save results in a file or compare results with some references value in a file
125 //------------------------------------------------------------------------------------
126 trace_p->pop();
127 if(check) { return diffResFile(out, rootname); }
128 else { return saveResToFile(out, rootname); }
129 }
130
131 }
132