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_operator.cpp
19 \author E. Lunéville
20 \since 23 feb 2012
21 \date 11 may 2012
22
23 Low level tests of Operator 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
33 #include <iostream>
34 #include <fstream>
35
36 using namespace xlifepp;
37
38 namespace unit_operator {
39
40 // scalar spectral function
sinBasis(const Point & P,Parameters & pa=defaultParameters)41 Real sinBasis(const Point& P, Parameters& pa = defaultParameters)
42 {
43 Real x = P.x();
44 Real h = pa("h"); // get the parameter h (user definition)
45 Int n = pa("basis index"); // get the index of function to compute
46 return std::sqrt(2. / h) * std::sin(n * pi_ * x / h); //computation
47 }
48
49 //vector spectral function
VecBasis(const Point & P,Parameters & pa=defaultParameters)50 Vector<Real> VecBasis(const Point& P, Parameters& pa = defaultParameters)
51 {
52 Real x = P.x();
53 Real h = pa("h"); // get the parameter h (user definition)
54 Int n = pa("basis index"); // get the index of function to compute
55 Vector<Real> res(2);
56 res(1) = std::sqrt(2. / h) * std::sin(n * pi_ * x / h); // computation
57 res(2) = std::sqrt(2. / h) * std::cos(n * pi_ * x / h);
58 return res;
59 }
60
61 //vector spectral function
VecBasis3(const Point & P,Parameters & pa=defaultParameters)62 Vector<Complex> VecBasis3(const Point& P, Parameters& pa = defaultParameters)
63 {
64 Real x = P.x();
65 Real h = pa("h"); // get the parameter h (user definition)
66 Int n = pa("basis index"); // get the index of function to compute
67 Vector<Complex> res(3);
68 res(1) = std::sqrt(2. / h) * std::sin(n * pi_ * x / h);
69 res(2) = std::sqrt(2. / h) * std::cos(n * pi_ * x / h);
70 res(3) = Complex(0,1.);
71 return res;
72 }
73
scalF(const Point & P,Parameters & pa=defaultParameters)74 Real scalF(const Point& P, Parameters& pa = defaultParameters)
75 {
76 Real x = P.x();
77 return x; //computation
78 }
79
VecF(const Point & P,Parameters & pa=defaultParameters)80 Vector<Real> VecF(const Point& P, Parameters& pa = defaultParameters)
81 {
82 return Vector<Real>(1,0.);
83 }
84
85
86 // matrix function
MatF(const Point & P,Parameters & pa=defaultParameters)87 Matrix<Real> MatF(const Point& P, Parameters& pa = defaultParameters)
88 {
89 Matrix<Real> res(2, 2);
90 Real x = P.x();
91 res(1, 1) = x; res(1, 2) = 0.;
92 res(2, 1) = 0.5; res(2, 2) = 2 * x;
93 return res;
94 }
95
MatFc(const Point & P,Parameters & pa=defaultParameters)96 Matrix<Complex> MatFc(const Point& P, Parameters& pa = defaultParameters)
97 {
98 Matrix<Complex> res(3, 3);
99 Real x = P.x();
100 res(1, 1) = Complex(0,1); res(2, 2) = x; res(3, 3) = 1.;
101 return res;
102 }
103
104 // scalar kernel function
G(const Point & M,const Point & P,Parameters & pa=defaultParameters)105 Real G(const Point& M, const Point& P, Parameters& pa = defaultParameters)
106 {
107 // Real r = M.distance(P);
108 return 1.; // computation
109 }
110
unit_operator(bool check)111 String unit_operator(bool check)
112 {
113 String rootname = "unit_operator";
114 trace_p->push(rootname);
115 std::stringstream out; // string stream receiving results
116 out.precision(testPrec);
117
118 // test of DifferentialOperator
119 DifferentialOperator* do1_p = findDifferentialOperator(_id);
120 out << "difop id : name=" << do1_p->name() << " order=" << do1_p->order()
121 << " normal=" << words(do1_p->normalRequired())
122 << " extension=" << words(do1_p->extensionRequired()) << "\n";
123 findDifferentialOperator(_dt);
124 findDifferentialOperator(_dx);
125 findDifferentialOperator(_dy);
126 findDifferentialOperator(_dz);
127 findDifferentialOperator(_grad);
128 findDifferentialOperator(_div);
129 findDifferentialOperator(_curl);
130 findDifferentialOperator(_gradS);
131 findDifferentialOperator(_gradG);
132 findDifferentialOperator(_divS);
133 findDifferentialOperator(_divG);
134 findDifferentialOperator(_curlS);
135 findDifferentialOperator(_curlG);
136 findDifferentialOperator(_epsilon);
137 findDifferentialOperator(_epsilonG);
138 findDifferentialOperator(_ntimes);
139 findDifferentialOperator(_timesn);
140 findDifferentialOperator(_ncrossntimes);
141 findDifferentialOperator(_timesncrossn);
142 findDifferentialOperator(_timesn);
143 findDifferentialOperator(_ndot);
144 findDifferentialOperator(_ndotgrad);
145 findDifferentialOperator(_ndiv);
146 findDifferentialOperator(_ncross);
147 findDifferentialOperator(_ncrosscurl);
148 findDifferentialOperator(_ncrossncross);
149 verboseLevel(2);
150 printListDiffOp(out);
151
152 // test of Value
153 Value::printValueTypeRTINames(out);
154 Real s = 1.;
155 Vector<Real> vr(3); Vector<Complex> vc(3, i_);
156 Matrix<Real> mr(3, 3, 1.); Matrix<Complex> mc(3, 3, i_);
157 Value Vs(s); out << "real Value Vs(s) :" << Vs;
158 Value Vi(i_); out << "complex Value Vi(i) :" << Vi;
159 Value Vvr(vr); out << "real vector Value Vvr(vr) :" << Vvr;
160 Value Vvc(vc); out << "complex vector Value Vvc(vc) :" << Vvc;
161 Value Vmr(mr); out << "real matrix Value Vmr(mr) :" << Vmr;
162 Value Vmc(mc); out << "complex matrix Value Vmc(mc) :" << Vmc;
163 out << "get real Value s= " << Vs.value<Real>() << "\n";
164 out << "get complex Value i= " << Vi.value<Complex>() << "\n";
165 out << "get real vector Value vr= " << Vvr.value<Vector<Real> >() << "\n";
166 out << "get complex vector Value vc= " << Vvc.value<Vector<Complex> >() << "\n";
167 out << "get real matrix Value mr= " << Vmr.value<Matrix<Real> >() << "\n";
168 out << "get complex matrix Value mc= " << Vmc.value<Matrix<Complex> >() << "\n";
169
170 // test of OperatorOnUnknown
171 Mesh msh; //fake mesh
172 GeomDomain Omega(msh, "Omega", 3);
173 Parameters ps;
174 ps<<Parameter(1.,"h")<<Parameter(1,"basis index");
175 Function F(sinBasis, "sinBasis", ps);
176 Function VF(VecBasis, "(sin,cos)", ps);
177 Function MF(MatF, "mat function", ps);
178 Space V(Omega, Function(sinBasis, "sinBasis", ps), 10, "sinus basis");
179 Unknown u(V,"u");
180 out << "definition of unknown u : " << u;
181 Unknown w(V, "w", 3);
182 out << "definition of unknown w : " << w;
183
184 // test constructor and assignment
185 OperatorOnUnknown opu(u), opgu(grad(u));
186 out << "test OperatorOnUnknown opu(u) : " << opu;
187 out << "test OperatorOnUnknown opgu(grad(u)) : " << opgu;
188 opu = d0(u);
189 out << "test opu=d0(u) : " << opu;
190 out << "test opgu=opu : " << opgu;
191 // test all primary operator
192 out << "test id(u) : " << id(u);
193 out << "test d0(u) : " << d0(u);
194 out << "test dt(u) : " << dt(u);
195 out << "test d1(u) : " << d1(u);
196 out << "test dx(u) : " << dx(u);
197 out << "test d2(u) : " << d2(u);
198 out << "test dy(u) : " << dy(u);
199 out << "test d3(u) : " << d3(u);
200 out << "test dz(u) : " << dz(u);
201 out << "test grad(u) : " << grad(u);
202 out << "test nabla(u) : " << nabla(u);
203 out << "test div(w) : " << div(w);
204 out << "test curl(w) : " << curl(w);
205 out << "test rot(w) : " << rot(w);
206 // out << "test gradS(u) : " << gradS(u);
207 // out << "test nablaS(u) : " << nablaS(u);
208 // out << "test divS(w) : " << divS(w);
209 // out << "test curlS(w) : " << curlS(w);
210 // out << "test rotS(w) : " << rotS(w);
211 out << "test nx(u) : " << nx(u);
212 out << "test ndot(w) : " << ndot(w);
213 out << "test ncross(w) : " << ncross(w);
214 out << "test ncrossncross(w) : " << ncrossncross(w);
215 out << "test ndotgrad(u) : " << ndotgrad(u);
216 out << "test ndiv(w) : " << ndiv(w);
217 out << "test ncrosscurl(w) : " << ncrosscurl(w);
218 out << "test epsilon(w) : " << epsilon(w);
219 out << "test epsilonG(w) : " << epsilonG(w,1,0,0);
220 out << "test epsilonG(w) : " << epsilonG(w,1,0,1,0);
221 // out << "test gradG(u,1.) : " << gradG(u, 1.);
222 // out << "test nablaG(u,1,2) : " << nablaG(u, 1, 2);
223 // out << "test divG(w,1,2,3) : " << divG(w, 1, 2, 3);
224 // out << "test curlG(w,1,2,3) : " << curlG(w, 1, 2, 3);
225 // out << "test rotG(w,0,1) : " << rotG(w, 0, 1);
226
227 // test aliases
228 out << "test n*u : " << _n* u;
229 out << "test u*n : " << u* _n;
230 out << "test n|w : " << (_n | w);
231 out << "test w|n : " << (w | _n);
232 out << "test n^w : " << (_n ^ w);
233 out << "test n^(n^w) : " << (_n ^ (_n ^ w));
234 out << "test n|grad(u) : " << (_n | grad(u));
235 out << "test grad(u)|n : " << (grad(u) | _n);
236 out << "test n*div(w) : " << (_n * div(w));
237 out << "test div(w)*n : " << (div(w)*_n);
238 out << "test n^curl(w) : " << (_n ^ curl(w));
239 out << "test n^rot(w) : " << (_n ^ rot(w));
240
241 // test left and right operations with function
242 out << "test F*opu : " << (F * opu);
243 out << "test opu*F : " << (opu * F);
244 out << "test F*opgu*F : " << (F * opgu * F);
245 out << "test F*u : " << (F * u);
246 out << "test u*F : " << (u * F);
247 out << "test F*grad(u)*F : " << (F * grad(u)*F);
248 out << "test VF|grad(u)*F : " << (VF | grad(u)*F);
249 out << "test VF^grad(u)*VF : " << ((VF ^ grad(u)) * VF); //cross product in 2D returns a scalar !
250 out << "test MF*grad(u)*F : " << (MF * grad(u)*F);
251 out << "test (MF*grad(u))|VF : " << ((MF * grad(u)) | VF);
252
253 // test using explicit C++ functions
254 out << "test u*scalF : " << (u * scalF);
255 out << "test scalF*u : " << (scalF * u);
256 out << "test grad(u)*scalF : " << (grad(u)*scalF);
257 out << "test (adj(MatF)*grad(u))|VecF : " << ((adj(MatF)*grad(u)) | VecF);
258
259 // test left and right operations with constant
260 out << "test u*vr: " << (u * vr);
261 out << "test s*grad(u) : " << (s * grad(u));
262 out << "test i*grad(u) : " << (i_ * grad(u));
263 out << "test vr|grad(u) : " << (vr | grad(u));
264 out << "test vc|grad(u) : " << (vc | grad(u));
265 out << "test mr*grad(u) : " << (mr * grad(u));
266 out << "test mc*grad(u) : " << (mc * grad(u));
267 out << "test s*grad(u)*i : " << (s * grad(u)*i_);
268 out << "test i*grad(u)|vr : " << (i_ * grad(u) | vr);
269 out << "test (vr|grad(u))*mr : " << ((vr | grad(u))*mr);
270
271 // test conjugate/transpose and adjoint operator
272 out << "test grad(conj(u)) : " << grad(conj(u));
273 out << "test trans(mr)*grad(conj(u)) : " << (trans(mr)*grad(conj(u)));
274 out << "test grad(conj(u))|conj(vc) : " << (grad(conj(u)) | conj(vc));
275 out << "test (adj(mc)*grad(conj(u)))|conj(vc) : " << (adj(mc)*grad(conj(u)) | conj(vc));
276 out << "test (adj(mc)*grad(conj(u)))|adj(vc) : " << (adj(mc)*grad(conj(u)) | adj(vc));
277
278 // test expressions with component
279 out << "test grad(w[1]) : " << grad(w[1]);
280 out << "test grad(w[2])*s : " << (grad(w[2])*s);
281 out << "test (adj(mc)*grad(conj(w[2])))|adj(vc) : " << (adj(mc)*grad(conj(w[2])) | adj(vc));
282 out << "test adj(mc)*(grad(conj(w[2])))*adj(mc)) : " << (adj(mc) * (grad(conj(w[2]))*adj(mc)));
283
284 //test kernel form
285 Kernel ker=Helmholtz3dKernel(1.);
286 //Kernel ker=Laplace3dKernel();
287 out << "\ntest ker*u : " << ker*u;
288 out << "\ntest u*ker : " << u*ker;
289 out << "\ntest u*ker*u : " << u*ker*u;
290
291 //out << "\ntest G*u : " << G*u; //not working on visual studio, ambiguity ?
292 out << "\ntest u*G : " << u*G;
293 out << "\ntest u*G*u : " << u*G*u;
294 out << "\ntest u*(grad_x(ker)|_nx)*u : " << u*(grad_x(ker)|_nx)*u;
295 out << "\ntest u*grad_y(grad_x(ker)|_nx)*u : " << u*grad_y(grad_x(ker)|_nx)*u;
296 out << "\ntest u*grad_y(grad_x(ker)|_nx)|ny*u : " << u*(grad_y(grad_x(ker)|_nx)|_ny)*u;
297 out << "\ntest u*(_nx^(_ny^ker))*u : " << u*(_nx^(_ny^ker))*u;
298 out << "\ntest u*(grad_y(grad_x(ker))*u : " << u*(grad_y(grad_x(ker)))*u;
299 //out << "\ntest u*(_nx^(_ny^G))*u : " << u*(_nx^(_ny^G))*u; //ambiguity between ny^Function et ny^Kernel, G may be cast to both
300
301 //test operator on function
302 out << "\ntest id(F) : " << id(F);
303 out << "\ntest ntimes(F) : " << ntimes(F);
304 out << "\ntest ndot(VF) : " << ndot(VF);
305 out << "\ntest ncross(VF) : " << ncross(VF);
306 out << "\ntest ncrossncross(VF) : " << ncrossncross(VF);
307 out << "\ntest n*F : " << _n*F;
308 out << "\ntest F*n : " << F *_n;
309 out << "\ntest (n^n)*F : " << ((_n^_n)*F);
310 out << "\ntest n|VF : " << (_n|VF);
311 out << "\ntest VF|*n : " << (VF |_n);
312 out << "\ntest n^VF : " << (_n^VF);
313 out << "\ntest n^(n^VF) : " << (_n^(_n^VF));
314 out << "\ntest (n^n)^VF) : " << ((_n^_n)^VF);
315
316 //evaluate operator on function
317 Point M(1.);
318 Real val;
319 Vector<Real> v_val, n(2); n[0]=1.; n[1]=1.;
320 out << "\ntest id(F)(M) : " << id(F).eval(M,val);
321 out << "\ntest id(VF)(M) : " << id(VF).eval(M,v_val);
322 out << "\ntest (VF|_n)(M) : " << (VF|_n).eval(M,val,&n);
323 out << "\ntest (_n^VF)(M) : " << (_n^VF).eval(M,val,&n); //2D ->scalar result
324 out << "\ntest (MatF*_n)(M) : " << (MatF*_n).eval(M,v_val,&n);
325 out << "\ntest (_n*MatF)(M) : " << (_n*MatF).eval(M,v_val,&n);
326
327 Matrix<Real> A(2,2);A(1,1)=10;
328 out << "\ntest (A*_n)(M) : " << (A*_n).eval(M,v_val,&n);
329
330 Function H(VecBasis3, "(sin,cos,1)", ps);
331 n.resize(3);n[2]=0.;
332 Complex cval;
333 Vector<Complex> cv_val;
334 out << "\ntest id(H)(M) : " << id(H).eval(M,cv_val);
335 out << "\ntest (H|_n)(M) : " << (H|_n).eval(M,cval,&n);
336 out << "\ntest (_n^H)(M) : " << (_n^H).eval(M,cv_val,&n);
337 out << "\ntest (MatFc*_n)(M) : " << (MatFc*_n).eval(M,cv_val,&n);
338
339 //test operator involving TermVector
340 verboseLevel(25);
341 Mesh msq(Square(_origin=Point(0.,0.),_length=1.,_nnodes=5,_domain_name="Omega"),_triangle,1,_structured,"msq");
342 Domain omega=msq.domain("Omega");
343 Space Vh(omega,P1,"Vh"); Unknown uh(Vh,"uh");
344 TermVector x1(uh,omega,_x1,"x1");
345 out << "test x1*uh : " << x1*uh;
346 out << "test x1*uh*x1 : " << x1*uh*x1;
347 out << "test (x1^3)*uh : " << (x1^3)*uh;
348 out << "test intg(omega,(x1^2)*uh) : " << intg(omega,(x1^2)*uh);
349
350
351 //------------------------------------------------------------------------------------
352 // save results in a file or compare results with some references value in a file
353 //------------------------------------------------------------------------------------
354 trace_p->pop();
355 if (check) { return diffResFile(out, rootname); }
356 else { return saveResToFile(out, rootname); }
357 }
358
359 }
360