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 elasticity problem in 2D and 3D
19
20 \file unit_epsilon.cpp
21 \author E. Lunéville
22 \since 16 mai 2013
23 \date 16 mai 2013
24 */
25
26 #include "xlife++-libs.h"
27 #include "testUtils.hpp"
28
29 #include <iostream>
30 #include <fstream>
31 #include <vector>
32
33
34 using namespace std;
35 using namespace xlifepp;
36
37 namespace unit_epsilon {
38
Epsilon_Pk(const GeomDomain & omg,int k,Real h,std::ostream & out)39 void Epsilon_Pk(const GeomDomain& omg, int k, Real h, std::ostream& out)
40 {
41 stringstream s;
42 trace_p->push("Epsilon_Pk");
43
44 //create Lagrange Pk space and unknown
45 Interpolation& Pk=interpolation(Lagrange, standard, k, H1);
46 Space Vk(omg, Pk, "Vk", true);
47 Unknown u(Vk, "u", 2); TestFunction v(u, "v");
48 TermVector UN(u,omg,1.,"un");
49
50 // Create bilinear form
51 TermMatrix M(intg(omg, u | v), "M");
52 TermMatrix S(intg( omg, epsilon(u)%epsilon(v) ), "S");
53 TermMatrix D(intg( omg, div(u)*div(v) ), "D");
54 BilinearForm auv=intg( omg, epsilon(u)%epsilon(v) ) + intg(omg, div(u) * div(v) ) + intg(omg, u | v);
55 TermMatrix A(auv,"a(u,v)");
56 out<<"M*UN|UN = "<<real((M*UN|UN))<<" S*UN|UN = "<<real((S*UN|UN))<<" D*UN|UN = "<<real((D*UN|UN))<<" A*UN|UN = "<<real((A*UN|UN))<<eol;
57
58 Number NbP = A.begin()->second->entries()->rmEntries_p->nbRows;
59 Number NbPs = A.begin()->second->entries()->rmEntries_p->nbRowsSub;
60 out << "Matrice de " << NbP << "x" << NbP << " blocs: " << NbPs << "x" << NbPs << std::endl;
61 TermVector B = A*UN;
62 TermVector E=directSolve(A,B);
63 E-=UN;
64 out << "C0 error " << norminfty(E) << std::endl;
65 trace_p->pop();
66 }
67
68 // Programme principal:
unit_epsilon(bool check)69 String unit_epsilon(bool check)
70 {
71 String rootname = "unit_epsilon";
72 trace_p->push(rootname);
73
74 std::stringstream out;
75 out.precision(testPrec);
76 string meshInfo;
77 verboseLevel(100);
78 Number nbDivMin , nbDivMax;
79 Number ordMin, ordMax;
80 //Real h;
81
82 //numberOfThreads(1);
83
84 //create a mesh and Domains TRIANGLES
85 out << " =====================================================================================" << eol;
86 out << " 2D test on square" << eol;
87 out << " =====================================================================================" << eol;
88 meshInfo= " Triangle mesh of square";
89 nbDivMin =5, nbDivMax=5;
90 ordMin = 1, ordMax=1;
91 for(Number div = nbDivMin; div <= nbDivMax; div+=5)
92 {
93 Real h=1./div;
94 Mesh mesh2d(Rectangle(_xmin=0, _xmax=1, _ymin=0, _ymax=1, _nnodes=div+1, _domain_name = "Omega"), _triangle, 1, _structured, meshInfo);
95 Domain omega = mesh2d.domain("Omega");
96 out << "# Triangle mesh: " << div << " sudivisions, h = " << h << eol;
97 for(Number ord=ordMin; ord<=ordMax; ord++)
98 {
99 out << " -> Triangles P" << ord << " -------------------" << eol;
100 Epsilon_Pk(omega, ord, h, out);
101 }
102 }
103
104 // out << " =====================================================================================" << eol;
105 // out << " 3D test on cube" << eol;
106 // out << " =====================================================================================" << eol;
107 // meshInfo=" Tetrahedra mesh of cube ";
108 // nbDivMin =5, nbDivMax=15;
109 // ordMin = 1, ordMax=1;
110 // for(Number div = nbDivMin; div <= nbDivMax; div+=5)
111 // {
112 // Mesh mesh3d(Cube(_origin=Point(0.,0.,0.),_length=1.,_nnodes=div+1,_domain_name="Omega"), _tetrahedron, 1, _gmsh, meshInfo);
113 // Domain omega = mesh3d.domain("Omega");
114 // h = 1. /div;
115 // out << "# Tetrahedra mesh : " << div << " sudivisions, h = " << h << eol;
116 // for(Number ord = ordMin; ord <= ordMax; ord++)
117 // {
118 // out << " -> Tetrahedra P" << ord << " -------------------" << eol;
119 // Epsilon_Pk(omega, ord, h, out);
120 // }
121 // }
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 }
133