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_Maxwell3D_BEM.cpp
20 \author E. Luneville
21 \since 18 November 2017
22 \date 18 November 2017
23 */
24
25 #include "xlife++-libs.h"
26 #include "testUtils.hpp"
27
28 using namespace std;
29 using namespace xlifepp;
30
31 namespace sys_Maxwell3D_BEM
32 {
Ei(const Point & P,Parameters & pars)33 Vector<complex_t> Ei(const Point& P, Parameters& pars)
34 {
35 Vector<real_t> incPol(3,0.); incPol(1)=1.;
36 Point incDir(0.,0.,1.) ;
37 Real k = pars("k");
38 return incPol*exp(i_*k * dot(P,incDir));
39 }
40
Eixn(const Point & P,Parameters & pars)41 Vector<Complex> Eixn(const Point& P, Parameters& pars)
42 {
43 Vector<Real> n=P;
44 Real n2=norm2(n);
45 if(n2>theEpsilon) n/=n2;
46 return crossProduct(Ei(P,pars),n);
47 }
48
gc(const Point & P,const Point & Q,Parameters & pars)49 Vector<Real> gc(const Point& P, const Point& Q,Parameters& pars)
50 {
51 return Vector<Real>(3,1.);
52 }
53
Eixnxn(const Point & P,Parameters & pars)54 Vector<Complex> Eixnxn(const Point& P, Parameters& pars)
55 {
56 Vector<Real> n=P;
57 Real n2=norm2(n);
58 if(n2>theEpsilon) n/=n2;
59 return crossProduct(crossProduct(Ei(P,pars),n),n);
60 }
61
sys_Maxwell3D_BEM(bool check)62 String sys_Maxwell3D_BEM(bool check)
63 {
64 String rootname = "sys_Maxwell3D_BEM";
65 trace_p->push(rootname);
66 std::stringstream out;
67 out.precision(testPrec);
68 //numberOfThreads(1);
69 verboseLevel(20);
70 out.precision(5);
71 out.setf(ios::scientific);
72
73 //define parameters and functions
74 Real k= 1, radius=1.;
75 Parameters pars;
76 pars << Parameter(k, "k") << Parameter(radius, "radius");
77 Kernel H = Helmholtz3dKernel(pars); // load Helmholtz3D kernel
78 Function fi(Ei, pars); // define right hand side function
79 Function fixn(Eixn, pars); // define right hand side function
80 Function fixnxn(Eixnxn, pars); // define right hand side function
81 Function Eex(scatteredFieldMaxwellExn,pars);
82
83 // stuff to do integral representation and compute error on segment [(0,0,1.5),(0,0,20)]
84 Mesh mseg(Segment(_v1=Point(0,0,5),_v2=Point(0,0,20),_nnodes=100,_domain_name="S"),triangle);
85 Domain S = mseg.domain("S");
86 Space VS(S,P1,"VS");
87 Unknown WS(VS,"WS",3);
88 TermVector Eexa(WS, S, Eex);
89 TermMatrix M(intg(S,WS|WS));
90 Real nex=sqrt(real(M*Eexa|Eexa));
91
92 // meshing the unit sphere
93 Number npa=10; //nb of points by diameter/4 of sphere
94 out<<"npa="<<npa<<eol;
95 Sphere sphere(_center=Point(0,0,0), _radius=radius, _nnodes=npa, _domain_name="Gamma");
96 Mesh meshSh(sphere, triangle, 1, gmsh);
97 Domain Gamma = meshSh.domain("Gamma");
98
99 //define FE-RT1 space and unknown
100 Space RTh(Gamma, RT_1, "RTh");
101 Unknown U(RTh,"U"); TestFunction V(U,"V");
102 Space Vh(Gamma, P1, "Vh"); Unknown W(Vh,"W",3);
103
104 IntegrationMethods ims(_SauterSchwabIM,4,0.,_defaultRule,5,2.,_defaultRule,3);
105 IntegrationMethods im(_defaultRule,10,1.,_defaultRule, 5);
106 TermVector Fx(W,Gamma,fixn); TermVector FxRT=projection(Fx,RTh);FxRT.name("FxRT");
107
108 //compute indirect EFIE system and solve it
109 BilinearForm as((1./k)*intg(Gamma,Gamma, div(U)*H*div(V),ims)-k*intg(Gamma,Gamma, U*H|V,ims));
110 TermMatrix SS(as,"S");
111 TermVector BIE(intg(Gamma,fi|V),"BIE");
112 TermVector LambdaIE = directSolve(SS,-BIE,true);
113 TermVector EIE=-(1./k)*integralRepresentation(WS, S, intg(Gamma,grad_x(H)*div(U),im), LambdaIE)
114 -k*integralRepresentation(WS, S, intg(Gamma,H*U,im), LambdaIE);
115 TermVector erIE=EIE-Eexa;
116 Real nerIE=sqrt(real(M*erIE|erIE));
117 out<<"indirect EFIE L2 error = "<<nerIE<<" relative error = "<<nerIE/nex<<eol;
118
119 //compute direct EFIE system and solve it
120 TermMatrix C(intg(Gamma,Gamma,U|(grad_x(H)^V),ims),"C");
121 TermVector BDE=C*FxRT-0.5*BIE; //intg(Gamma,Gamma,fixn|(grad_x(H)^V),ims))-0.5*intg(Gamma,fixnxn|V))
122 TermVector LambdaDE = directSolve(SS,BDE,true);
123 TermVector EDE= integralRepresentation(WS, S, intg(Gamma,grad_x(H)^U,im),FxRT)
124 -(1./k)*integralRepresentation(WS, S, intg(Gamma,grad_x(H)*div(U),im), LambdaDE)
125 -k*integralRepresentation(WS, S, intg(Gamma,H*U,im), LambdaDE);
126 TermVector erDE=EDE-Eexa;
127 Real nerDE=sqrt(real(M*erDE|erDE));
128 out<<" direct EFIE L2 error = "<<nerDE<<" relative error = "<<nerDE/nex<<eol;
129
130 //compute indirect MFIE system and solve it
131 TermMatrix N(intg(Gamma,(U^_n)|V),"N");
132 TermMatrix AIM=0.5*N-C;
133 TermVector LambdaIM = directSolve(AIM,BIE,true);
134 TermVector EIM= -integralRepresentation(WS, S, intg(Gamma,grad_x(H)^U,im),LambdaIM);
135 TermVector erIM=EIM-Eexa;
136 Real nerIM=sqrt(real(M*erIM|erIM));
137 out<<"indirect MFIE L2 error = "<<nerIM<<" relative error = "<<nerIM/nex<<eol;
138
139 //compute direct MFIE system and solve it
140 TermMatrix ADM=0.5*N+C;
141 TermVector LambdaDM = directSolve(ADM,SS*FxRT,true);
142 TermVector EDM= integralRepresentation(WS, S, intg(Gamma,grad_x(H)^U,im),FxRT)
143 -(1./k)*integralRepresentation(WS, S, intg(Gamma,grad_x(H)*div(U),im), LambdaDM)
144 -k*integralRepresentation(WS, S, intg(Gamma,H*U,im), LambdaDM);
145 TermVector erDM=EDM-Eexa;
146 Real nerDM=sqrt(real(M*erDM|erDM));
147 out<<" direct MFIE L2 error = "<<nerDM<<" relative error = "<<nerDM/nex<<eol;
148
149 //compute CFIE system and solve it
150 Unknown theta(RTh,"theta"); TestFunction xi(theta,"xi");
151 Real eta=0.25; Complex ieta=i_*eta;
152 BilinearForm ac= (ieta/k)*intg(Gamma,Gamma, div(U)*H*div(V),ims)-(ieta*k)*intg(Gamma,Gamma, U*H|V,ims)
153 + 0.5*intg(Gamma,(theta^_n)|V)-intg(Gamma,Gamma,theta|(grad_x(H)^V),ims)
154 + intg(Gamma,theta|xi)+intg(Gamma,div(theta)*div(xi))+intg(Gamma,(U^_n)|xi);
155 TermMatrix AC(ac,"AC");
156 TermVector BC(intg(Gamma,fi|V));
157 TermVector UT=directSolve(AC,BC);
158 TermVector LambdaC=UT(U), Theta=UT(theta);
159 TermVector EC =(ieta/k)*integralRepresentation(WS, S, intg(Gamma,grad_x(H)*div(U),im), LambdaC)
160 +(ieta*k)*integralRepresentation(WS, S, intg(Gamma,H*U,im), LambdaC)
161 - integralRepresentation(WS, S, intg(Gamma,grad_x(H)^theta,im),Theta);
162 TermVector erC=EC-Eexa;
163 Real nerC=sqrt(real(M*erC|erC));
164 out<<" CFIE L2 error = "<<nerC<<" relative error = "<<nerC/nex<<eol;
165
166 //------------------------------------------------------------------------------------
167 // save results in a file or compare results with some references value in a file
168 //------------------------------------------------------------------------------------
169 trace_p->pop();
170 if(check) { return diffResFile(out, rootname); }
171 else { return saveResToFile(out, rootname); }
172 }
173
174
175 }
176