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 Test Laplace problem with Neumann condition
19    -delta(u)+u=f
20    du/dn=0
21 
22 \file sys_Laplace_Neumann.cpp
23 \author E. Lunéville
24 \since 16 mai 2013
25 \date 16 mai 2013
26 */
27 
28 #include "xlife++-libs.h"
29 #include "testUtils.hpp"
30 
31 using namespace std;
32 using namespace xlifepp;
33 
34 namespace sys_Laplace_Neumann {
35 //#define _TermVector(x,y) TermVector(x)((y),#x)
36 //#define _TermVector(x,y) TermVector x(y);x.name()=#x;
37 //      _TermVector(Er,U-Uex);
38 //#define _(x,y,z1) x y(z1);y.name()=#y;
39 //#define _(x,y,...) x y(__VA_ARGS__);y.name()=#y;
40 //      _(TermVector,Er,U-Uex);
41 //#define _Term(x,y,...) Term ##x y(__VA_ARGS__);y.name()=#y;
42 //       _Term(Vector,Er,U-Uex);
43 
uno(const Point & P,Parameters & pa=defaultParameters)44 Real uno(const Point& P, Parameters& pa = defaultParameters)
45 {
46   return 1.;
47 }
48 
cosx1(const Point & P,Parameters & pa=defaultParameters)49 Real cosx1(const Point& P, Parameters& pa = defaultParameters)
50 {
51   Real x=P(1);
52   return cos(pi_ * x);
53 }
54 
cosx2(const Point & P,Parameters & pa=defaultParameters)55 Real cosx2(const Point& P, Parameters& pa = defaultParameters)
56 {
57   Real x=P(1), y=P(2);
58   return cos(pi_ * x) * cos(pi_ * y);
59 }
60 
cosx3(const Point & P,Parameters & pa=defaultParameters)61 Real cosx3(const Point& P, Parameters& pa = defaultParameters)
62 {
63   Real x=P(1), y=P(2), z=P(3);
64   return cos(pi_ * x) * cos(pi_ * y) * cos(pi_ * z);
65 }
66 
fx(const Point & P,Parameters & pa=defaultParameters)67 Real fx(const Point& P, Parameters& pa=defaultParameters)
68 {
69   return P(1);
70 }
fy(const Point & P,Parameters & pa=defaultParameters)71 Real fy(const Point& P, Parameters& pa=defaultParameters)
72 {
73   return P(2);
74 }
fz(const Point & P,Parameters & pa=defaultParameters)75 Real fz(const Point& P, Parameters& pa=defaultParameters)
76 {
77   return P(3);
78 }
79 
fx2(const Point & P,Parameters & pa=defaultParameters)80 Real fx2(const Point& P, Parameters& pa=defaultParameters)
81 {
82   return P(1)*P(1);
83 }
fy2(const Point & P,Parameters & pa=defaultParameters)84 Real fy2(const Point& P, Parameters& pa=defaultParameters)
85 {
86   return P(2)*P(2);
87 }
fz2(const Point & P,Parameters & pa=defaultParameters)88 Real fz2(const Point& P, Parameters& pa=defaultParameters)
89 {
90   return P(3)*P(3);
91 }
92 
93 ofstream fout;
94 
Laplace_Neumann_Pk(const GeomDomain & omg,int k,Real h,const String & E,std::ostream & out)95 void Laplace_Neumann_Pk(const GeomDomain& omg, int k,  Real h, const String& E, std::ostream& out)
96 {
97   stringstream s;
98   trace_p->push("Laplace_Neumann_Pk");
99 
100   //create Lagrange Pind space and unknown
101   Interpolation& Pk=interpolation(Lagrange, standard, k, H1);
102   Space Vk(omg, Pk, "Vk", false);
103   Unknown u(Vk,"u");  TestFunction v(u, "v");
104 
105   //create bilinear form and linear form
106   BilinearForm muv=intg(omg, u * v);
107   BilinearForm auv = intg(omg, grad(u) | grad(v)) + intg(omg, u * v);
108   BilinearForm ruv=intg(omg, grad(u) | grad(v));
109 
110   LinearForm fv;
111   Real a;
112   TermVector Uex;
113   TermVector un(u,omg,1.,"un");
114 
115   if(omg.dim()==1)
116   {
117     fv=intg(omg, cosx1 * v);
118     Uex = TermVector(u, omg, cosx1);
119     a=1. + pi_ * pi_;
120   }
121 
122   if(omg.dim()==2)
123   {
124     fv=intg(omg, cosx2 * v);
125     Uex = TermVector(u, omg, cosx2);
126     a=1. + 2.*pi_ * pi_;
127   }
128   if(omg.dim()==3)
129   {
130     fv=intg(omg, cosx3 * v);
131     Uex=TermVector(u, omg, cosx3);
132     a = 1. + 3.*pi_ * pi_;
133   }
134   Uex /=a;
135 
136   elapsedTime();
137 
138   // create FEterm
139   TermMatrix A(auv, "a(u,v)");
140   TermMatrix R(ruv, "R");
141   TermMatrix M(muv, "M");
142   TermVector B(fv, "f(v)");
143   elapsedTime();
144   elapsedTime("compute matrices");
145 
146   // check FE computation
147   if(omg.dim()==2)
148   {
149     TermVector X(u,omg,fx,"x");
150     TermVector Y(u,omg,fy,"y");
151     TermVector X2(u,omg,fx2,"X2");
152     TermVector Y2(u,omg,fy2,"Y2");
153     out<<"-----------------------------------------------------------------------------------------"<<eol;
154     out<<"FE computation check "<< E << k << ", h=" << h << ", nb dl=" << Uex.size() << eol;
155     out<<"   norminfty(R*1) = "<<norminfty(R*un);
156     out<<", (M*1|1) = "<<(M*un|un);
157     out<<", (A*1|1) = "<<(A*un|un)<<eol;
158     out<<"  (R*x|x) = "<<(R*X|X);
159     out<<", (R*y|y) = "<<(R*Y|Y);
160     out<<", (R*x2|x2) = "<<(R*X2|X2);
161     out<<", (R*y2|y2) = "<<(R*Y2|Y2);
162     out<<", (B|1) = "<<(B|un)<<eol;
163     out<<"-----------------------------------------------------------------------------------------"<<eol;
164   }
165 
166   if(omg.dim()==3)
167   {
168     TermVector X(u,omg,fx,"x");
169     TermVector Y(u,omg,fy,"y");
170     TermVector Z(u,omg,fz,"z");
171     TermVector X2(u,omg,fx2,"X2");
172     TermVector Y2(u,omg,fy2,"Y2");
173     TermVector Z2(u,omg,fz2,"Z2");
174     out<<"-----------------------------------------------------------------------------------------"<<eol;
175     out<<"FE computation check "<< E << k << ", h=" << h << ", nb dl=" << Uex.size() << eol;
176     out<<"   norminfty(R*1) = "<<norminfty(R*un);
177     out<<", (M*1|1) = "<<(M*un|un);
178     out<<", (A*1|1) = "<<(A*un|un)<<eol;
179     out<<"  (R*x|x) = "<<(R*X|X);
180     out<<", (R*y|y) = "<<(R*Y|Y);
181     out<<", (R*z|z) = "<<(R*Z|Z);
182     out<<", (R*x2|x2) = "<<(R*X2|X2);
183     out<<", (R*y2|y2) = "<<(R*Y2|Y2);
184     out<<", (R*z2|z2) = "<<(R*Z2|Z2);
185     out<<", (B|1) = "<<(B|un)<<eol;
186     out<<"-----------------------------------------------------------------------------------------"<<eol;
187   }
188 
189   SuTermMatrix* suterm = A.begin()->second;
190 
191   LargeMatrix<Complex>* cA = suterm->entries()->cEntries_p;
192   Number cRow, cCol, rRow, rCol,nZA;
193   if (0 != cA)
194   {
195     out << "Complex matrix ";
196     out << " storage :" << words("storage type",cA->storageType());
197     cRow = cA->nbRows;
198     cCol = cA->nbCols;
199     out << " matrix size " << cRow << " x " << cCol ;
200     out << " mumber non-zero " << cA->nbNonZero() << std::endl;
201   }
202 
203   LargeMatrix<Real>* rA = suterm->entries()->rEntries_p;
204   if (0 != rA)
205   {
206     out << "Real matrix ";
207     out << " storage :" << words("storage type",rA->storageType());
208     rRow = rA->nbRows;
209     rCol = rA->nbCols;
210     out << " size " <<  rRow << " x " << rCol;
211     out << " number non-zero " << rA->nbNonZero() << std::endl;
212     nZA=rA->nbNonZero();
213   }
214   //solve linear systesm AX=F using direct method
215   TermVector U = directSolve(A, B);
216 
217   TermVector res=A*U-B;
218   out<<"residu = "<<norminfty(res)<<eol;
219 
220   //compute errors
221   TermVector Dq = U - Uex;
222   Real errL2, errH1, errC0, nL2, nH1, nC0, cpu;
223   errL2 = sqrt(abs(M * Dq | Dq)) ;
224   errH1 = sqrt(abs(A * Dq | Dq)) ;
225   errC0= norminfty(Dq);
226   nL2 = sqrt(abs(M * Uex | Uex));
227   nH1 = sqrt(abs(A * Uex | Uex));
228   nC0 = norminfty(Uex);
229   cpu = elapsedTime();
230 
231   theCout << "--------------- " <<E<< k << ", h=" << h << ", nb dl=" << U.size() << " ---------------" << eol;
232   theCout << "L2 Error = " << errL2 << ", Rel. L2 error = " << errL2 / nL2 << eol;
233   theCout << "H1 error = " << errH1 << ", Rel. H1 error = " << errH1 / nH1 << eol;
234   theCout << "C0 error = " << errC0 << ", Rel. C0 error = " << errC0 / nC0 << eol;
235   theCout << " -> cpu time = " << cpu << eol;
236 
237   out << "--------------- "<<E << k << ", h=" << h << ", nb dl=" << U.size() << " ---------------" << eol;
238   out << "L2 norm Uex = "<<nL2<<", L2 Error = " << errL2 << ", Rel. L2 error = " << errL2 / nL2 << eol;
239   out << "H1 norm Uex = "<<nH1<<", H1 error = " << errH1 << ", Rel. H1 error = " << errH1 / nH1 << eol;
240   out << "C0 norm Uex = "<<nC0<<", C0 error = " << errC0 << ", Rel. C0 error = " << errC0 / nC0 << eol;
241   out << " -> cpu time = " << cpu << eol;
242 
243   fout<<k<<" "<<h<<" "<<U.size()<<" "<<rA->nbNonZero()<<" "<<errL2<<" "<<errH1<<" "<<cpu<<eol;
244 
245   saveToFile("U",U,_vtu);
246   saveToFile("Uex",Uex,_vtu);
247 //  out<<U<<eol;
248 //  out<<Uex<<eol;
249 //  out<<Dq<<eol;
250   trace_p->pop();
251 }
252 
253 
254 //
255 // Programme principal:
sys_Laplace_Neumann(bool check)256 String sys_Laplace_Neumann(bool check)
257 {
258   String rootname = "sys_Laplace_Neumann";
259   trace_p->push(rootname);
260 
261   //numberOfThreads(1);
262 
263    bool doSegment     = false,
264         doTriangle    = false,
265         doQuadrangle  = false,
266         doTetrahedron = false,
267         doHexahedron  = false,
268         doPrism       = false,
269         doPyramid     = false,
270         doAll         = true;
271 
272    std::stringstream out;
273    out.precision(testPrec);
274    string meshInfo;
275    verboseLevel(0);
276    Number nbDivMin = 1, nbDivMax=3;
277    Number ordMin =1, ordMax=3;
278    std::vector<String> sidenames(4, "");
279    Real h;
280 
281   if (doSegment || doAll)
282   {
283     nbDivMax=1;
284     ordMin=1;
285     ordMax=3;
286     Number fac=8;
287      out <<"====================================================================================" << eol;
288     out << "=                                TEST 1D Segment                                   =" << eol;
289     out << "====================================================================================" << eol<<eol;;
290     // create a mesh and Domains SEGMENT
291     meshInfo= " mesh of [0,1] with nbStep =  ";
292 
293     fout.open("segment_errors.dat");
294     for(Number div = nbDivMin; div <= nbDivMax; div++)
295     {
296       Number n=fac*div;
297       meshInfo += tostring(n);
298       h=1./n;
299       Mesh mesh1d(Segment(_xmin=0,_xmax=1,_nnodes=n+1,_domain_name="Omega"), 1, _structured, meshInfo);
300       Domain omega = mesh1d.domain("Omega");
301       theCout << "# SEGMENT mesh:  " << n << " sudivisions, h = " <<h <<eol;
302       theCout << "===========================================" << eol;
303       out << "# SEGMENT mesh:  " << n << " sudivisions, h = " <<h <<eol;
304       out << "===========================================" << eol;
305 
306       for(Number ord=ordMin; ord<=ordMax; ord++)
307       {
308         theCout << " -> Segments P" <<ord<< "-------------------" <<eol;
309         out << " -> Segments P" <<ord<< "-------------------" <<eol;
310         Laplace_Neumann_Pk(omega, ord, h,"P", out);
311       }
312     }
313     fout.close();
314   }
315 
316   if (doTriangle || doAll)
317   {
318     nbDivMin = 2;
319     nbDivMax = 2;
320     ordMin=1;
321     ordMax=3;
322     out << "====================================================================================" << eol;
323     out << "=                             TEST 2D Triangle                                     =" << eol;
324     out << "====================================================================================" << eol<<eol;;
325     // create a mesh and Domains TRIANGLES
326     meshInfo= " TRIANGULAR mesh of [0,1]x[0,1] with nbStep =  ";
327 
328     fout.open("triangle_errors.dat");
329     for(Number div = nbDivMin; div <= nbDivMax; div++)
330     {
331       meshInfo += tostring(div);
332       h=1./(4*div);
333       Mesh mesh2d(Rectangle(_xmin=0,_xmax=1,_ymin=0,_ymax=1,_nnodes=4*div+1), _triangle, 1, _structured, meshInfo);
334       Domain omega = mesh2d.domain("Omega");
335       theCout << "# TRIANGULAR mesh:  " << div << " sudivisions, h = " <<h <<eol;
336       theCout << "===========================================" << eol;
337       out << "# TRIANGULAR mesh:  " << div << " sudivisions, h = " <<h <<eol;
338       out << "===========================================" << eol;
339       for(Number ord=ordMin; ord<=ordMax; ord++)
340       {
341         theCout << " -> Triangles P" <<ord<< "-------------------" <<eol;
342         out << " -> Triangles P" <<ord<< "-------------------" <<eol;
343         Laplace_Neumann_Pk(omega, ord, h,"P", out);
344       }
345     }
346     fout.close();
347   }
348 
349   if (doQuadrangle || doAll)
350   {
351     nbDivMin = 2;
352     nbDivMax = 2;
353     ordMin=1;
354     ordMax=3;
355     out << "====================================================================================" << eol;
356     out << "=                           TEST 2D Quadrangle                                     =" << eol;
357     out << "====================================================================================" << eol<<eol;;
358     // create a mesh and Domains QUADRANGLES sidenames[0] = "x=0";
359     meshInfo= " QUADRANGULAR mesh of [0,1]x[0,1] with nbStep =  ";
360     fout.open("quadrangle_errors.dat");
361     for(Number div = nbDivMin; div <= nbDivMax; div++)
362     {
363       meshInfo += tostring(div);
364       h=1./(4*div);
365       Mesh mesh2d(Rectangle(_xmin=0,_xmax=1,_ymin=0,_ymax=1,_nnodes=4*div+1), _quadrangle, 1, _structured, meshInfo);
366       Domain omega=mesh2d.domain("Omega");
367       theCout << " #  QUADRANGULAR mesh           " << div << " sudivisions, h = " <<h <<eol;
368       theCout << "=============================================" << eol;
369       out << " #  QUADRANGULAR mesh           " << div << " sudivisions, h = " <<h <<eol;
370       out << "==============================================" << eol;
371       for(Number ord=ordMin; ord<=ordMax; ord++)
372       {
373         theCout << " -> Quadrangles Q" <<ord<< eol;
374         out << " -> Quadrangles Q" <<ord<< eol;
375         Laplace_Neumann_Pk( omega, ord, h,"Q", out);
376       }
377     }
378     fout.close();
379   }
380 
381   // create a mesh and Domains TETRAHEDRON
382   if (doTetrahedron || doAll)
383   {
384     out << "====================================================================================" << eol;
385     out << "=                           TEST 3D Tetrahedron                                    =" << eol;
386     out << "====================================================================================" << eol<<eol;
387     meshInfo=" TETRAHEDRAL mesh of cube ";
388     nbDivMin = 2;
389     nbDivMax = 2;
390     ordMin=1;
391     ordMax=3;
392     fout.open("tetrahedron_errors.dat");
393     for(Number div = nbDivMin; div <= nbDivMax; div++)
394     {
395       Mesh mesh3d(Cube(_center=Point(0.5,0.5,0.5),_length=1.,_nnodes=4*div), _tetrahedron, 1, _subdiv, meshInfo);
396       mesh3d.renameDomain(0,"Omega");
397       Domain omega = mesh3d.domain("Omega");
398       h=1./(4*div);
399       //h = 2. / std::pow(2,static_cast<int>(div));
400       meshInfo += tostring(div);
401       theCout << "# TETRAHEDRAL mesh :  " << div << " sudivisions, h = " <<h <<eol;
402       theCout << "=========================================" << eol;
403       out << "# TETRAHEDRAL mesh :  " << div << " sudivisions, h = " <<h <<eol;
404       out << "=========================================" << eol;
405       for(Number ord = ordMin; ord <= ordMax; ord++)
406       {
407         theCout << " -> Tetrahedra P" <<ord<< eol;
408         out << " -> Tetrahedra P" <<ord<< eol;
409         Laplace_Neumann_Pk(omega, ord, h,"P", out);
410         theCout << eol;
411       }
412     }
413     fout.close();
414   }
415 
416   // create a hexaedron mesh and domains
417   if (doHexahedron || doAll)
418   {
419     nbDivMin=2;
420     nbDivMax=2;
421     ordMin=1;
422     ordMax=3;
423     out << "====================================================================================" << eol;
424     out << "=                           TEST 3D Hexahedron                                     =" << eol;
425     out << "====================================================================================" << eol<<eol;;
426     meshInfo=" HEXAHEDRAL mesh of cube [0,1]x[0,1]x[0,1] ";
427     fout.open("hexahedron_errors.dat");
428     for(Number div = nbDivMin; div <= nbDivMax; div++)
429     {
430       Mesh mesh3d(Cuboid(_xmin=0,_xmax=1,_ymin=0,_ymax=1,_zmin=0,_zmax=1,_nnodes=4*div+1), _hexahedron, 1, _structured, meshInfo);
431       Domain omega = mesh3d.domain("Omega");
432       h =1./(4*div);
433       meshInfo += tostring(div);
434       theCout << " #  HEXAHEDRAL mesh:   " << div << " sudivisions, h = " <<h <<eol;
435       theCout << "==============================================" << eol;
436       out << " #  HEXAHEDRAL mesh:   " << div << " sudivisions, h = " <<h <<eol;
437       out << "==============================================" << eol;
438       for(Number ord = ordMin; ord <= ordMax; ord++)
439       {
440         theCout << " -> Hexahedra Q" <<ord<< eol;
441         out << " -> Hexahedra Q" <<ord<< eol;
442         Laplace_Neumann_Pk(omega, ord, h,"Q", out);
443       }
444     }
445     fout.close();
446   }
447 
448   if (doPrism || doAll)
449   {
450     nbDivMin=2;
451     nbDivMax=2;
452     ordMin=1;
453     ordMax=2;
454     out << "====================================================================================" << eol;
455     out << "=                           TEST 3D Prism order 1/2                                =" << eol;
456     out << "====================================================================================" << eol<<eol;;
457     meshInfo=" mesh of cube [0,1]x[0,1]x[0,1] from 2D extrusion";
458     fout.open("prism_errors.dat");
459     for(Number div = nbDivMin; div <= nbDivMax; div++)
460     {
461       h =1./(4*div);
462       Mesh mesh2d(Rectangle(_xmin=0,_xmax=1,_ymin=0,_ymax=1,_nnodes=4*div+1), _triangle, 1, _structured, meshInfo);
463       Mesh mesh3d(mesh2d,Point(0.,0.,0.),Point(0.,0.,1.),4*div);
464       Domain omega = mesh3d.domain("Omega");
465       mesh3d.saveToFile("mesh3d_prism",_msh);
466 
467       meshInfo += tostring(div);
468       theCout << " #  PRISMATIC mesh:   " << div << " sudivisions, h = " <<h <<eol;
469       theCout << "==============================================" << eol;
470       out << " #  PRISMATIC mesh:   " << div << " sudivisions, h = " <<h <<eol;
471       out << "==============================================" << eol;
472       for(Number ord = ordMin; ord <= ordMax; ord++)
473       {
474         theCout << " -> Prism order " <<ord<< eol;
475         out << " -> Prism order " <<ord<< eol;
476         Laplace_Neumann_Pk(omega, ord, h,"Pr", out);
477       }
478     }
479     fout.close();
480   }
481 
482   if (doPyramid || doAll)
483   {
484     nbDivMin=2;
485     nbDivMax=2;
486     ordMin=1;
487     ordMax=2;
488     out << "====================================================================================" << eol;
489     out << "=                           TEST 3D Pyramid order 1/2                                =" << eol;
490     out << "====================================================================================" << eol<<eol;;
491     meshInfo=" mesh of cube [0,1]x[0,1]x[0,1]";
492     fout.open("prism_errors.dat");
493     for(Number div = nbDivMin; div <= nbDivMax; div++)
494     {
495       h =1./(4*div);
496       Mesh mesh3d(Cuboid(_xmin=0,_xmax=1,_ymin=0,_ymax=1,_zmin=0,_zmax=1,_nnodes=4*div+1), _hexahedron, 1, _structured, meshInfo);
497       Mesh mesh3dPy(mesh3d, _pyramid);
498       Domain omega = mesh3d.domain("Omega");
499       mesh3dPy.saveToFile("mesh3d_pyramid",_msh);
500 
501       meshInfo += tostring(div);
502       theCout << " #  PYRAMIDAL mesh:   " << div << " sudivisions, h = " <<h <<eol;
503       theCout << "==============================================" << eol;
504       out << " #  PYRAMIDAL mesh:   " << div << " sudivisions, h = " <<h <<eol;
505       out << "==============================================" << eol;
506       for(Number ord = ordMin; ord <= ordMax; ord++)
507       {
508         theCout << " -> Pyramid order " <<ord<< eol;
509         out << " -> Pyramid order " <<ord<< eol;
510         Laplace_Neumann_Pk(omega, ord, h,"Py", out);
511       }
512     }
513     fout.close();
514   }
515 
516   //------------------------------------------------------------------------------------
517   // save results in a file or compare results with some reference values in a file
518   //------------------------------------------------------------------------------------
519   trace_p->pop();
520   if (check) { return diffResFile(out, rootname); }
521   else { return saveResToFile(out, rootname); }
522 }
523 
524 }
525 
526