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 \file unit_TermMatrix.cpp
19 \author E. Lunéville
20 \since 5 october 2012
21 \date 11 september 2013
22
23 Low level tests of TestMatrix class and related classes.
24 Almost functionalities are checked.
25 This function may either creates a reference file storing the results (check=false)
26 or compares results to those stored in the reference file (check=true)
27 It returns reporting informations in a string
28 */
29
30 #include "xlife++-libs.h"
31 #include "testUtils.hpp"
32 #include "../../src/term/SymbolicTermMatrix.hpp"
33
34 using namespace xlifepp;
35
36 namespace unit_TermMatrix {
37
38 // scalar 2D function
un(const Point & P,Parameters & pa=defaultParameters)39 Real un(const Point& P, Parameters& pa = defaultParameters)
40 {
41 return 1.;
42 }
43
xfy(const Point & P,Parameters & pa=defaultParameters)44 Real xfy(const Point& P, Parameters& pa = defaultParameters)
45 {
46 return P(1)*P(2);
47 }
48
G(const Point & P,const Point & Q,Parameters & pa=defaultParameters)49 Real G(const Point& P, const Point& Q,Parameters& pa = defaultParameters)
50 {
51 return 1.;
52 }
53
unit_TermMatrix(bool check)54 String unit_TermMatrix(bool check)
55 {
56 String rootname = "unit_TermMatrix";
57 trace_p->push(rootname);
58 std::stringstream out; // string stream receiving results
59 out.precision(testPrec);
60 //numberOfThreads(1);
61
62 verboseLevel(30);
63
64 // create a mesh and Domains
65 std::vector<String> sidenames(4,"");
66 sidenames[0]="Gamma_1";
67 sidenames[3]="Gamma_2";
68 verboseLevel(10);
69 Mesh mesh2d(Rectangle(_xmin=0,_xmax=1,_ymin=0,_ymax=1,_nnodes=3,_side_names=sidenames),_triangle,1,_structured,"P1 mesh of [0,1]x[0,1]");
70 Domain omega=mesh2d.domain("Omega");
71 Domain gamma1=mesh2d.domain("Gamma_1");
72 // Domain gamma2=mesh2d.domain("Gamma_2");
73
74 // create Lagrange P1 space and unknown
75 Interpolation *LagrangeP1=findInterpolation(Lagrange,standard,1,H1);
76 Space V1(omega,*LagrangeP1,"V1",false);
77 Unknown u(V1,"u");
78 TestFunction v(u,"v");
79 Function f1(un,"one");
80 Function fxy(xfy,"xy");
81 cpuTime("time for initialisation, mesh and space definition");
82
83 // useful TermVector
84 TermVector un(u,omega,1.,"un");
85 TermVector unv(v,omega,1.,"unv");
86
87 // test some volumic mass bilinear forms
88 BilinearForm uv=intg(omega,u*v);
89 TermMatrix M(uv,"UV");
90 cpuTime("time for computation of M");
91 M.viewStorage(out);
92 out<<M;
93 out<<"----> M*un|un="<<realPart(M*un|un)<<"\n\n";
94
95 TermMatrix Md(uv,_csDualStorage,"UV");
96 cpuTime("time for computation of Md");
97 Md.viewStorage(out);
98 out<<Md;
99 out<<"----> Md*un|un="<<realPart(Md*un|un)<<"\n\n";
100 clear(Md);
101
102 BilinearForm uv1=intg(omega,f1*u*v);
103 TermMatrix M1(uv1,"f1*UV");
104 cpuTime("time for computation of M1");
105 out<<M1;
106 //clear(M1);
107
108 BilinearForm xyuv=intg(omega,xfy*u*v);
109 TermMatrix Mxy(xyuv,"xy*UV");
110 cpuTime("time for computation of Mxy");
111 out<<Mxy;
112 out<<"----> Mxy*un|un="<<realPart(Mxy*un|un)<<"\n\n";
113 clear(Mxy);
114
115 BilinearForm uxyv=intg(omega,u*xfy*v);
116 TermMatrix Mxy2(uxyv,"U*xy*V");
117 cpuTime("time for computation of Mxy2");
118 out<<Mxy2;
119 clear(Mxy2);
120
121 BilinearForm uvxy=intg(omega,u*(v*xfy));
122 TermMatrix Mxy3(uvxy,"UV*xy");
123 cpuTime("time for computation of Mxy3");
124 out<<Mxy3;
125 clear(Mxy3);
126
127 // test some volumic rigidity bilinear forms
128 BilinearForm gugv=intg(omega,grad(u)|grad(v));
129 TermMatrix K(gugv,"K");
130 cpuTime("time for computation of K");
131 out<<K;
132 out<<"----> K*un="<<K*un<<"\n\n";
133
134 BilinearForm gugv1=intg(omega,f1*grad(u)|grad(v));
135 TermMatrix K1(gugv1,"K1");
136 cpuTime("time for computation of K1");
137 out<<K1;
138 clear(K1);
139
140 BilinearForm gugvxy=intg(omega,xfy*grad(u)|grad(v));
141 TermMatrix Kxy(gugvxy,"Kxy");
142 cpuTime("time for computation of Kxy");
143 out<<Kxy;
144 clear(Kxy);
145
146 BilinearForm dxuv=intg(omega,dx(u)*v);
147 TermMatrix Dxuv(dxuv,"DxUV");
148 cpuTime("time for computation of Dxuv");
149 out<<Dxuv;
150 out<<"----> DxUV*un="<<Dxuv*un<<"\n\n";
151 clear(Dxuv);
152
153 // ============ problem of deallocation (clear) with this guy ================
154 TermMatrix Dxuvd(dxuv,"_skylineDualStorage, DxUVd");
155 cpuTime("time for computation of Dxuvd");
156 out<<Dxuvd;
157 out<<"----> DxUVd*un="<<Dxuvd*un<<"\n\n";
158 clear(Dxuvd);
159
160 BilinearForm udxv=intg(omega,u*dx(v));
161 TermMatrix Dxvu(udxv,"DxVU");
162 cpuTime("time for computation of Dxvu");
163 out<<Dxvu;
164 out<<"----> DxVU*un="<<Dxvu*un<<"\n\n";
165
166 // test some surfacic mass bilinear forms
167 BilinearForm mg1=intg(gamma1,u*v);
168 TermMatrix Mg1(mg1,"Mg1");
169 cpuTime("time for computation of Mg1");
170 out<<Mg1;
171 TermVector ung(u,gamma1,1.,"un_gamma1");
172 out<<"----> Mg1*ung|ung="<<realPart(Mg1*ung|ung)<<"\n\n";
173
174 TermMatrix Mg1d(mg1,_skylineSymStorage,"Mg1d");
175 cpuTime("time for computation of Mg1d");
176 out<<Mg1d;
177 out<<"----> Mg1d*ung|ung="<<realPart(Mg1d*ung|ung)<<"\n\n";
178
179 // test some surfacic bilinear forms requiring domain extension
180 TermMatrix Dxuv1(intg(gamma1,dx(u)*v),"Dxuv1");
181 cpuTime("time for computation of Dxuv1");
182 out<<Dxuv1;
183 TermMatrix Dxvu1(intg(gamma1,u*dx(v)),"Dxvu1");
184 cpuTime("time for computation of Dxvu1");
185 out<<Dxvu1;
186 TermMatrix Dyuv1(intg(gamma1,dy(u)*v),"Dyuv1");
187 cpuTime("time for computation of Dyuv1");
188 out<<Dyuv1;
189 TermMatrix Dnuv1(intg(gamma1,ndotgrad(u)*v),"Dnuv1");
190 cpuTime("time for computation of Dnuv1");
191 out<<Dnuv1;
192 TermMatrix Dtudtv(intg(gamma1,gradS(u)|gradS(v)),"Dtudtv");
193 cpuTime("time for computation of Dtudtv");
194 out<<Dtudtv;
195 TermMatrix DtuGdtv(intg(gamma1,gamma1,gradS(u)|(G*gradS(v))),"DtuGdtv");
196 cpuTime("time for computation of DtuGdtv");
197 out<<DtuGdtv;
198
199 //test some combination of bilinear forms
200 Real k=1,alpha=1.;
201 BilinearForm h=intg(omega,grad(u)|grad(v))-k*k*intg(omega,u*v)+alpha*intg(gamma1,u*v);
202 TermMatrix H(h,"H");
203 cpuTime("time for computation of H");
204 out<<H;
205
206 // test storage conversion
207 H.setStorage(_skyline,_sym);
208 out<<H;
209 H.setStorage(_skyline,_dual);
210 out<<H;
211 H.setStorage(_dense,_row);
212 out<<H;
213 H.setStorage(_cs,_dual);
214 out<<H;
215
216 //test of linear combination of TermMatrix's
217 TermMatrix Hlc=K-(k*k)*M+alpha*Mg1;
218 //Hlc.name()="Hlc";
219 cpuTime("time for computation of Hlc");
220 out<<Hlc;
221
222 clear(Hlc);
223
224 //test some combination of bilinear forms
225 BilinearForm h2=intg(omega,grad(u)|grad(v))+intg(omega,u*dx(v))+alpha*intg(gamma1,u*v);
226 TermMatrix H2(h2,"H2");
227 cpuTime("time for computation of H2");
228 out<<H2;
229
230 //test of linear combination of TermMatrix's
231 TermMatrix H2lc=K+Dxvu+alpha*Mg1;
232 //H2lc.name()="H2lc";
233 cpuTime("time for computation of H2lc");
234 out<<H2lc;
235 //H2lc=H2lc-H2;
236 H2lc-=H2;
237 out<<H2lc;
238 clear(H2lc);
239 clear(H2);
240
241 // test product and inverse
242 TermMatrix M12=M1*M1;
243 thePrintStream<<"M12=M1*M1="<<M12;
244 TermMatrix invM1M1f=directSolve(M1,M1,_keep);
245 out<<"invM1M1f=directSolve(M1,M1)="<<invM1M1f;
246 TermMatrix Id(M1,_idMatrix,"Id");
247 out<<Id;
248 out<<"directSolve(M1,Id,_keep)="<<directSolve(M1,Id,_keep);
249 TermMatrix invM1=inverse(M1);
250 out<<"invM1=inverse(M1)="<<invM1;
251 out<<"invM1*M1="<<invM1*M1;
252 out<<"invM1*M12="<<invM1*M12;
253 TermMatrix MT=invM1*M12-M1;
254 out<<"norm(invM1*M12-M1)="<<norm2(MT)<<eol;
255 out<<"norm(invM1*M12-M1)="<<norm2(invM1*M12-M1)<<eol;
256 out<<"out<<invM1*M12-M1="<<invM1*M12-M1<<eol;
257
258 thePrintStream<<"===== avant Mpq ====="<<eol;
259 printMatrixStorages(thePrintStream);
260
261 //test SymbolicTermMatrix
262 TermMatrix Mf;
263 factorize(M,Mf);
264 SymbolicTermMatrix S=~M+2*inv(M)*M*(~M+~M);
265 theCout<<"S="<<S<<eol;
266 TermVector X(u,omega,1.);
267 TermVector SX=S*X;
268 TermVector S5=5*M*X;
269 theCout<<"|S5-SX|="<<norm2(S5-SX)<<eol;
270
271 // //other interpolations
272 // Interpolation LagrangeP2(Lagrange,standard,2,H1);
273 // Space V2(omega,LagrangeP2,"V2",false);
274 // Unknown p(V2,"p");
275 // TestFunction q(p,"q");
276 // BilinearForm pq=intg(omega,p*q);
277 // TermMatrix Mpq(pq,"pq");
278 // cpuTime("time for computation of Mpq");
279 // out<<Mpq;
280 // TermVector pun(p,omega,1.,"pun");
281 // TermVector qun(q,omega,1.,"qun");
282 // out<<"\nMpq*un|un="<<realPart(Mpq*pun|qun)<<"\n\n";
283 // clear(Mpq);
284 //
285 // thePrintStream<<"===== avant Kpq ====="<<eol;
286 // printMatrixStorages(thePrintStream);
287 //
288 // BilinearForm gpgq=intg(omega,grad(p)|grad(q));
289 // TermMatrix Kpq(gpgq,"Kpq");
290 // cpuTime("time for computation of Kpq");
291 // out<<Kpq;
292 // out<<"\n"<<Kpq*pun<<"\n\n";
293 // clear(Kpq);
294 //
295 // thePrintStream<<"===== après Kpq ====="<<eol;
296 // printMatrixStorages(thePrintStream);thePrintStream.flush();
297 //
298 // Interpolation LagrangeP3(Lagrange,standard,3,H1);
299 // Space V3(omega,LagrangeP3,"V3",false);
300 // Unknown s(V3,"s");
301 // TestFunction t(s,"t");
302 // BilinearForm st=intg(omega,s*t);
303 // TermMatrix Mst(st,"st");
304 // cpuTime("time for computation of Mst");
305 // out<<Mst;
306 // TermVector sun(s,omega,1.,"pun");
307 // TermVector tun(t,omega,1.,"qun");
308 // out<<"\nMst*un|un="<<realPart(Mst*sun|tun)<<"\n\n";
309 // clear(Mst);
310 //
311 // BilinearForm gsgt=intg(omega,grad(s)|grad(t));
312 // TermMatrix Kst(gsgt,"Kst");
313 // cpuTime("time for computation of Kst");
314 // out<<Kst;
315 // out<<"\n"<<Kst*sun<<"\n";
316 // clear(Kst);
317 //
318 // Interpolation LagrangeP4(Lagrange,standard,4,H1);
319 // Space V4(omega,LagrangeP4,"V4",false);
320 // Unknown e(V4,"e");
321 // TestFunction f(e,"f");
322 // BilinearForm ef=intg(omega,e*f);
323 // TermMatrix Mef(ef,"ef");
324 // cpuTime("time for computation of Mef");
325 // out<<Mef;
326 // TermVector eun(e,omega,1.,"eun");
327 // TermVector fun(f,omega,1.,"fun");
328 // out<<"\nMef*un|un="<<realPart(Mef*eun|fun)<<"\n\n";
329 // clear(Mef);
330 //
331 // BilinearForm gegf=intg(omega,grad(e)|grad(f));
332 // TermMatrix Kef(gegf,"Kef");
333 // cpuTime("time for computation of Kef");
334 // out<<Kef;
335 // out<<"\n"<<Kef*eun<<"\n";
336 // clear(Kef);
337
338 //------------------------------------------------------------------------------------
339 // save results in a file or compare results with some references value in a file
340 //------------------------------------------------------------------------------------
341 trace_p->pop();
342 if (check) { return diffResFile(out, rootname); }
343 else { return saveResToFile(out, rootname); }
344
345 }
346
347 }
348