You are here: Concepts>Concepts Web>Class documentation

Class documentation of Concepts

Loading...
Searching...
No Matches
DtNmap2D_visc.hh
Go to the documentation of this file.
1
7#ifndef __DtNmap2D_visc__hh__
8#define __DtNmap2D_visc__hh__
9
10#include "formula.hh"
11#include "basics/exceptions.hh"
12#include "toolbox/sequence.hh"
13#include "space/integral.hh"
14#include "formula/elementFormulaContainer.hh"
15#include "formula/formulas2D.hh"
16#include "formula/bessel.hh"
17#include "formula/frmE_product.hh"
18#include "function/vector.hh"
19#include "operator/sparseMatrix.hh"
20#include "hp1D/linearForm.hh"
21
22namespace concepts
23{
24
25 const Cmplx cmplx_i(0,1);
26
27 Real sqr(const Real x)
28 {
29 return x*x;
30 }
31
32 enum dimproj
33 {
34 dimX, dimY, dimZ, dimdiv
35 };
36
37 int sgn(const Real d)
38 {
39 if (d<0) return -1;
40 if (d>0) return 1;
41 return 0;
42 }
43
44 Cmplx visc_ell_fast(const Cmplx lambda,
45 const Real nu,
46 const Real omega,
47 const Real rho0)
48 {
49 const Cmplx ellsquare = lambda*lambda + cmplx_i * omega * rho0 / nu;
50 const Real a = ellsquare.real(); const Real b = ellsquare.imag();
51 Cmplx ell;
52#if __cplusplus >= 201103L
53 ell.real(1.0/sqrt(2.0)*sqrt(sqrt(a*a+b*b)+a));
54 ell.imag(sgn(b)/sqrt(2.0)*sqrt(sqrt(a*a+b*b)-a));
55#else
56 ell.real() = 1.0/sqrt(2.0)*sqrt(sqrt(a*a+b*b)+a);
57 ell.imag() = sgn(b)/sqrt(2.0)*sqrt(sqrt(a*a+b*b)-a);
58#endif
59 return (cmplx_i*ell);
60 }
61
62
63 Cmplx visc_ell_slow(const Cmplx lambda,
64 const Real nu,
65 const Real omega,
66 const Real rho0,
67 const Real c0)
68 {
69 const Cmplx ellsquare = lambda*lambda + omega*omega / (c0*c0 - cmplx_i*omega*nu/rho0);
70 const Real a = ellsquare.real(); const Real b = ellsquare.imag();
71 Cmplx ell;
72#if __cplusplus >= 201103L
73 ell.real(1.0/sqrt(2.0)*sqrt(sqrt(a*a+b*b)+a));
74 ell.imag(sgn(b)/sqrt(2.0)*sqrt(sqrt(a*a+b*b)-a));
75#else
76 ell.real() = 1.0/sqrt(2.0)*sqrt(sqrt(a*a+b*b)+a);
77 ell.imag() = sgn(b)/sqrt(2.0)*sqrt(sqrt(a*a+b*b)-a);
78#endif
79 return (-cmplx_i*ell);
80 }
81
82 Cmplx lambda_limit(const Real omega, const Real c0, const int n, const Real L)
83 {
84 if (n*c0*M_PI<omega*L)
85 return cmplx_i*sqrt(sqr(omega/c0)-sqr(n*M_PI/L));
86 return sqrt(sqr(n*M_PI/L)-sqr(omega/c0));
87 }
88
89 template<uint dim>
90 class Wsym_x : public Formula<Cmplx> {
91 public:
93 Wsym_x(const int n=1,
94 const Real omega=1.0,
95 const Real c0=1.0,
96 const Real rho0=1.0,
97 const Real nu=0.1,
98 const Real L=1.0);
99 virtual Wsym_x<dim>* clone() const;
100 virtual Cmplx operator() (const Real p, const Real t=0.0) const;
101 virtual Cmplx operator() (const Real2d& p, const Real t=0.0) const;
102 virtual Cmplx operator() (const Real3d& p, const Real t=0.0) const;
103 virtual std::ostream& info(std::ostream& os) const;
104 virtual Cmplx get_lambda() const;
105 protected:
106 const int n_;
107 const Real omega_;
108 const Real c0_;
109 const Real rho0_;
110 const Real nu_;
111 const Real L_;
112
113 const Cmplx lambda;
114 const Cmplx ellf;
115 const Cmplx ells;
116 const Cmplx cf;
117 };
118
119 template<uint dim>
120 class Wsym_y : public Formula<Cmplx> {
121 public:
123 Wsym_y(const int n=1,
124 const Real omega=1.0,
125 const Real c0=1.0,
126 const Real rho0=1.0,
127 const Real nu=0.1,
128 const Real L=1.0);
129 virtual Wsym_y<dim>* clone() const;
130 virtual Cmplx operator() (const Real p, const Real t=0.0) const;
131 virtual Cmplx operator() (const Real2d& p, const Real t=0.0) const;
132 virtual Cmplx operator() (const Real3d& p, const Real t=0.0) const;
133 virtual std::ostream& info(std::ostream& os) const;
134 virtual Cmplx get_lambda() const;
135 protected:
136 const int n_;
137 const Real omega_;
138 const Real c0_;
139 const Real rho0_;
140 const Real nu_;
141 const Real L_;
142
143 const Cmplx lambda;
144 const Cmplx ellf;
145 const Cmplx ells;
146 const Cmplx cf;
147 };
148
149
150 template<uint dim>
151 class Wunsym_x : public Formula<Cmplx> {
152 public:
154 Wunsym_x(const int n=1,
155 const Real omega=1.0,
156 const Real c0=1.0,
157 const Real rho0=1.0,
158 const Real nu=0.1,
159 const Real L=1.0);
160 virtual Wunsym_x<dim>* clone() const;
161 virtual Cmplx operator() (const Real p, const Real t=0.0) const;
162 virtual Cmplx operator() (const Real2d& p, const Real t=0.0) const;
163 virtual Cmplx operator() (const Real3d& p, const Real t=0.0) const;
164 virtual std::ostream& info(std::ostream& os) const;
165 virtual Cmplx get_lambda() const;
166 protected:
167 const int n_;
168 const Real omega_;
169 const Real c0_;
170 const Real rho0_;
171 const Real nu_;
172 const Real L_;
173
174 const Cmplx lambda;
175 const Cmplx ellf;
176 const Cmplx ells;
177 const Cmplx cf;
178 };
179
180 template<uint dim>
181 class Wunsym_y : public Formula<Cmplx> {
182 public:
184 Wunsym_y(const int n=1,
185 const Real omega=1.0,
186 const Real c0=1.0,
187 const Real rho0=1.0,
188 const Real nu=0.1,
189 const Real L=1.0);
190 virtual Wunsym_y<dim>* clone() const;
191 virtual Cmplx operator() (const Real p, const Real t=0.0) const;
192 virtual Cmplx operator() (const Real2d& p, const Real t=0.0) const;
193 virtual Cmplx operator() (const Real3d& p, const Real t=0.0) const;
194 virtual std::ostream& info(std::ostream& os) const;
195 virtual Cmplx get_lambda() const;
196 protected:
197 const int n_;
198 const Real omega_;
199 const Real c0_;
200 const Real rho0_;
201 const Real nu_;
202 const Real L_;
203
204 const Cmplx lambda;
205 const Cmplx ellf;
206 const Cmplx ells;
207 const Cmplx cf;
208 };
209
210
211
212
213 template<class F, class G>
214 void addExactDtN_X_2Dcos_wp(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
215 const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
216 const G& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
217 {
218 conceptsAssert(L > 0, Assertion());
219 for (uint n = 0; n < DtNCoeff.size(); ++n) {
220 Cos_n_x cos_nx(+n, L);
221 hp1D::Riesz<Real> cLform(cos_nx);
222 Vector<Real> cVector(spc, cLform);
223 SparseMatrix<Real> cDtNmatrix(spc,cVector,cVector);
224 Cmplx factor = 2.0 / L;
225 if (n == 0) factor = 1.0 / L;
226 cDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
227 Cmplx integral = frm.projection(n);
228 rhs.add(cVector, DtNCoeffRhs[n]*integral);
229 }
230 }
231
232 template<class F, class G>
233 void addExactDtN_X_2Dsin_wp(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
234 const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
235 const G& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
236 {
237 conceptsAssert(L > 0, Assertion());
238 for (uint n = 1; n < DtNCoeff.size(); ++n) {
239 Sin_n_x sin_nx(+n, L);
240 hp1D::Riesz<Real> sLform(sin_nx);
241 Vector<Real> sVector(spc, sLform);
242 SparseMatrix<Real> sDtNmatrix(spc,sVector,sVector);
243 Cmplx factor = 2.0 / L;
244 sDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
245 Cmplx integral = frm.projection(-n);
246 rhs.add(sVector, DtNCoeffRhs[n]*integral);
247 }
248 }
249
250 template<class F, class G>
251 void addExactDtN_X_2Dcossin_wp(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
252 const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
253 const G& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
254 {
255 conceptsAssert(L > 0, Assertion());
256 for (uint n = 0; n < DtNCoeff.size(); ++n) {
257 Cos_n_x cos_nx(+n+1, L);
258 Sin_n_x sin_nx(+n+1, L);
259 hp1D::Riesz<Real> cLform(cos_nx);
260 Vector<Real> cVector(spc, cLform);
261 hp1D::Riesz<Real> sLform(sin_nx);
262 Vector<Real> sVector(spc, sLform);
263 // cos part over test fuction, sin part over unknown function
264 SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
265 Cmplx factor = 2.0 / L;
266 csDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
267 Cmplx integral = frm.projection(n+1);
268 rhs.add(sVector, DtNCoeffRhs[n]*integral);
269 }
270 }
271
272
273 template<class F, class G>
274 void addExactDtN_X_2Dsincos_wp(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
275 const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
276 const G& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
277 {
278 conceptsAssert(L > 0, Assertion());
279 for (uint n = 0; n < DtNCoeff.size(); ++n) {
280 Cos_n_x cos_nx(+n+1, L);
281 Sin_n_x sin_nx(+n+1, L);
282 hp1D::Riesz<Real> cLform(cos_nx);
283 Vector<Real> cVector(spc, cLform);
284 hp1D::Riesz<Real> sLform(sin_nx);
285 Vector<Real> sVector(spc, sLform);
286 // cos part over test fuction, sin part over unknown function
287 SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
288 Cmplx factor = 2.0 / L;
289 csDtNmatrix.addIntoT(dest, DtNCoeff[n]*factor);
290 Cmplx integral = frm.projection(-n-1);
291 rhs.add(cVector, DtNCoeffRhs[n]*integral);
292 }
293 }
294
295#ifdef HAS_MUMPS
296
298 template<uint dim>
299 void Wsym_x_GrammSchmidt(SparseMatrix<Cmplx>& Beta, const int n, const Real omega, const Real c0, const Real rho0, const Real nu, const Real L, const SpaceOnCells<Real>& spc)
300 {
301 conceptsAssert(Beta.dimX() == n, Assertion());
302 conceptsAssert(Beta.dimY() == n, Assertion());
303 for (int m=0; m<n;m ++)
304 {
305 Beta(m,m)=1.0;
306 if (m > 0)
307 {
308 // Create projection matrix with coefficients
309 SparseMatrix<Cmplx> CoeffsProj(m);
310 for (int l=0; l<m; l++)
311 for (int j=0; j<m; j++)
312 {
313 Wsym_x<dim> wsl(l, omega, c0, rho0, nu, L);
314 Wsym_x<dim> wsj(j, omega, c0, rho0, nu, L);
315 hp1D::Riesz<Cmplx> wslform(wsl);
316 hp1D::Riesz<Cmplx> wsjform(wsj);
317 Vector<Cmplx> Vwsl(spc, wslform);
318 Vector<Cmplx> Vwsj(spc, wsjform);
319 Cmplx ScalProjLJ = Vwsj * Vwsl;
320 for(int i=l; i<m; i++)
321 {
322 CoeffsProj(i,j) = CoeffsProj(i,j) + conj(Beta(i,l)) * ScalProjLJ;
323 }
324 }
325
326 // Create Right hand side
327 Vector<Cmplx> Rhs(m);
328 Vector<Cmplx> Sol(m);
329 Wsym_x<dim> wsm(m, omega, c0, rho0, nu, L);
330 hp1D::Riesz<Cmplx> wsmform(wsm);
331 Vector<Cmplx> Vwsm(spc, wsmform);
332 for (int l=0; l<m; l++)
333 {
334 Wsym_x<dim> wsl(l, omega, c0, rho0, nu, L);
335 hp1D::Riesz<Cmplx> wslform(wsl);
336 Vector<Cmplx> Vwsl(spc, wslform);
337 Cmplx ScalProjLM = Vwsm * Vwsl;
338 for (int j=l; j<m; j++)
339 {
340 Rhs(j) = Rhs(j) + conj(Beta(j,l)) * ScalProjLM;
341 }
342 }
343
344 // Solve system
345 if (m==1)
346 {
347 // Explicit solve of 1x1 system
348 Beta(1,0) = Rhs(0) / CoeffsProj(0,0);
349 }
350 else
351 {
352 std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
353 solver.reset(new concepts::Mumps<Cmplx>(CoeffsProj));
354 (*solver)(Rhs, Sol);
355 for (int l=0; l<m; l++)
356 Beta(m,l) = Sol(l);
357 }
358
359 }
360
361 // Renormalization of Betas
362 Wsym_x<dim> wsm(m, omega, c0, rho0, nu, L);
363 hp1D::Riesz<Cmplx> wsmform(wsm);
364 Vector<Cmplx> Vwsm(spc, wsmform);
365 for (int l=0; l<m; l++)
366 {
367 Wsym_x<dim> wsl(l, omega, c0, rho0, nu, L);
368 hp1D::Riesz<Cmplx> wslform(wsl);
369 Vector<Cmplx> Vwsl(spc, wslform);
370 Vwsl *= Beta(m,l);
371 Vwsm = Vwsm + Vwsl;
372 }
373 Real Norme = Vwsm.l2();
374 for(int l=0; l<=m; l++)
375 Beta(m,l) = Beta(m,l) / Norme;
376 }
377 }
379 template<uint dim>
380 void Wsym_y_GrammSchmidt(SparseMatrix<Cmplx>& Beta, const int n, const Real omega, const Real c0, const Real rho0, const Real nu, const Real L, const SpaceOnCells<Real>& spc)
381 {
382 conceptsAssert(Beta.dimX() == n, Assertion());
383 conceptsAssert(Beta.dimY() == n, Assertion());
384 for (int m=0; m<n;m ++)
385 {
386 Beta(m,m)=1.0;
387 if (m > 0)
388 {
389 // Create projection matrix with coefficients
390 SparseMatrix<Cmplx> CoeffsProj(m);
391 for (int l=0; l<m; l++)
392 for (int j=0; j<m; j++)
393 {
394 Wsym_y<dim> wsl(l, omega, c0, rho0, nu, L);
395 Wsym_y<dim> wsj(j, omega, c0, rho0, nu, L);
396 hp1D::Riesz<Cmplx> wslform(wsl);
397 hp1D::Riesz<Cmplx> wsjform(wsj);
398 Vector<Cmplx> Vwsl(spc, wslform);
399 Vector<Cmplx> Vwsj(spc, wsjform);
400 Cmplx ScalProjLJ = Vwsj * Vwsl;
401 for(int i=l; i<m; i++)
402 {
403 CoeffsProj(i,j) = CoeffsProj(i,j) + conj(Beta(i,l)) * ScalProjLJ;
404 }
405 }
406
407 // Create Right hand side
408 Vector<Cmplx> Rhs(m);
409 Vector<Cmplx> Sol(m);
410 Wsym_y<dim> wsm(m, omega, c0, rho0, nu, L);
411 hp1D::Riesz<Cmplx> wsmform(wsm);
412 Vector<Cmplx> Vwsm(spc, wsmform);
413 for (int l=0; l<m; l++)
414 {
415 Wsym_y<dim> wsl(l, omega, c0, rho0, nu, L);
416 hp1D::Riesz<Cmplx> wslform(wsl);
417 Vector<Cmplx> Vwsl(spc, wslform);
418 Cmplx ScalProjLM = Vwsm * Vwsl;
419 for (int j=l; j<m; j++)
420 {
421 Rhs(j) = Rhs(j) + conj(Beta(j,l)) * ScalProjLM;
422 }
423 }
424
425 // Solve system
426 if (m==1)
427 {
428 // Explicit solve of 1x1 system
429 Beta(1,0) = Rhs(0) / CoeffsProj(0,0);
430 }
431 else
432 {
433 std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
434 solver.reset(new concepts::Mumps<Cmplx>(CoeffsProj));
435 (*solver)(Rhs, Sol);
436 for (int l=0; l<m; l++)
437 Beta(m,l) = Sol(l);
438 }
439
440 }
441
442 // Renormalization of Betas
443 Wsym_y<dim> wsm(m, omega, c0, rho0, nu, L);
444 hp1D::Riesz<Cmplx> wsmform(wsm);
445 Vector<Cmplx> Vwsm(spc, wsmform);
446 for (int l=0; l<m; l++)
447 {
448 Wsym_y<dim> wsl(l, omega, c0, rho0, nu, L);
449 hp1D::Riesz<Cmplx> wslform(wsl);
450 Vector<Cmplx> Vwsl(spc, wslform);
451 Vwsl *= Beta(m,l);
452 Vwsm = Vwsm + Vwsl;
453 }
454 Real Norme = Vwsm.l2();
455 for(int l=0; l<=m; l++)
456 Beta(m,l) = Beta(m,l) / Norme;
457 }
458 }
459
461 template<uint dim>
462 void Wunsym_x_GrammSchmidt(SparseMatrix<Cmplx>& Beta, const int n, const Real omega, const Real c0, const Real rho0, const Real nu, const Real L, const SpaceOnCells<Real>& spc)
463 {
464 conceptsAssert(Beta.dimX() == n, Assertion());
465 conceptsAssert(Beta.dimY() == n, Assertion());
466 for (int m=0; m<n;m ++)
467 {
468 Beta(m,m)=1.0;
469 if (m > 0)
470 {
471 // Create projection matrix with coefficients
472 SparseMatrix<Cmplx> CoeffsProj(m);
473 for (int l=0; l<m; l++)
474 for (int j=0; j<m; j++)
475 {
476 Wunsym_x<dim> wsl(l, omega, c0, rho0, nu, L);
477 Wunsym_x<dim> wsj(j, omega, c0, rho0, nu, L);
478 hp1D::Riesz<Cmplx> wslform(wsl);
479 hp1D::Riesz<Cmplx> wsjform(wsj);
480 Vector<Cmplx> Vwsl(spc, wslform);
481 Vector<Cmplx> Vwsj(spc, wsjform);
482 Cmplx ScalProjLJ = Vwsj * Vwsl;
483 for(int i=l; i<m; i++)
484 {
485 CoeffsProj(i,j) = CoeffsProj(i,j) + conj(Beta(i,l)) * ScalProjLJ;
486 }
487 }
488
489 // Create Right hand side
490 Vector<Cmplx> Rhs(m);
491 Vector<Cmplx> Sol(m);
492 Wunsym_x<dim> wsm(m, omega, c0, rho0, nu, L);
493 hp1D::Riesz<Cmplx> wsmform(wsm);
494 Vector<Cmplx> Vwsm(spc, wsmform);
495 for (int l=0; l<m; l++)
496 {
497 Wunsym_x<dim> wsl(l, omega, c0, rho0, nu, L);
498 hp1D::Riesz<Cmplx> wslform(wsl);
499 Vector<Cmplx> Vwsl(spc, wslform);
500 Cmplx ScalProjLM = Vwsm * Vwsl;
501 for (int j=l; j<m; j++)
502 {
503 Rhs(j) = Rhs(j) + conj(Beta(j,l)) * ScalProjLM;
504 }
505 }
506
507 // Solve system
508 if (m==1)
509 {
510 // Explicit solve of 1x1 system
511 Beta(1,0) = Rhs(0) / CoeffsProj(0,0);
512 }
513 else
514 {
515 std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
516 solver.reset(new concepts::Mumps<Cmplx>(CoeffsProj));
517 (*solver)(Rhs, Sol);
518 for (int l=0; l<m; l++)
519 Beta(m,l) = Sol(l);
520 }
521
522 }
523
524 // Renormalization of Betas
525 Wunsym_x<dim> wsm(m, omega, c0, rho0, nu, L);
526 hp1D::Riesz<Cmplx> wsmform(wsm);
527 Vector<Cmplx> Vwsm(spc, wsmform);
528 for (int l=0; l<m; l++)
529 {
530 Wunsym_x<dim> wsl(l, omega, c0, rho0, nu, L);
531 hp1D::Riesz<Cmplx> wslform(wsl);
532 Vector<Cmplx> Vwsl(spc, wslform);
533 Vwsl *= Beta(m,l);
534 Vwsm = Vwsm + Vwsl;
535 }
536 Real Norme = Vwsm.l2();
537 for(int l=0; l<=m; l++)
538 Beta(m,l) = Beta(m,l) / Norme;
539 }
540 }
541
543 template<uint dim>
544 void Wunsym_y_GrammSchmidt(SparseMatrix<Cmplx>& Beta, const int n, const Real omega, const Real c0, const Real rho0, const Real nu, const Real L, const SpaceOnCells<Real>& spc)
545 {
546 conceptsAssert(Beta.dimX() == n, Assertion());
547 conceptsAssert(Beta.dimY() == n, Assertion());
548 for (int m=0; m<n;m ++)
549 {
550 Beta(m,m)=1.0;
551 if (m > 0)
552 {
553 // Create projection matrix with coefficients
554 SparseMatrix<Cmplx> CoeffsProj(m);
555 for (int l=0; l<m; l++)
556 for (int j=0; j<m; j++)
557 {
558 Wunsym_y<dim> wsl(l, omega, c0, rho0, nu, L);
559 Wunsym_y<dim> wsj(j, omega, c0, rho0, nu, L);
560 hp1D::Riesz<Cmplx> wslform(wsl);
561 hp1D::Riesz<Cmplx> wsjform(wsj);
562 Vector<Cmplx> Vwsl(spc, wslform);
563 Vector<Cmplx> Vwsj(spc, wsjform);
564 Cmplx ScalProjLJ = Vwsj * Vwsl;
565 for(int i=l; i<m; i++)
566 {
567 CoeffsProj(i,j) = CoeffsProj(i,j) + conj(Beta(i,l)) * ScalProjLJ;
568 }
569 }
570
571 // Create Right hand side
572 Vector<Cmplx> Rhs(m);
573 Vector<Cmplx> Sol(m);
574 Wunsym_y<dim> wsm(m, omega, c0, rho0, nu, L);
575 hp1D::Riesz<Cmplx> wsmform(wsm);
576 Vector<Cmplx> Vwsm(spc, wsmform);
577 for (int l=0; l<m; l++)
578 {
579 Wunsym_y<dim> wsl(l, omega, c0, rho0, nu, L);
580 hp1D::Riesz<Cmplx> wslform(wsl);
581 Vector<Cmplx> Vwsl(spc, wslform);
582 Cmplx ScalProjLM = Vwsm * Vwsl;
583 for (int j=l; j<m; j++)
584 {
585 Rhs(j) = Rhs(j) + conj(Beta(j,l)) * ScalProjLM;
586 }
587 }
588
589 // Solve system
590 if (m==1)
591 {
592 // Explicit solve of 1x1 system
593 Beta(1,0) = Rhs(0) / CoeffsProj(0,0);
594 }
595 else
596 {
597 std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
598 solver.reset(new concepts::Mumps<Cmplx>(CoeffsProj));
599 (*solver)(Rhs, Sol);
600 for (int l=0; l<m; l++)
601 Beta(m,l) = Sol(l);
602 }
603
604 }
605
606 // Renormalization of Betas
607 Wunsym_y<dim> wsm(m, omega, c0, rho0, nu, L);
608 hp1D::Riesz<Cmplx> wsmform(wsm);
609 Vector<Cmplx> Vwsm(spc, wsmform);
610 for (int l=0; l<m; l++)
611 {
612 Wunsym_y<dim> wsl(l, omega, c0, rho0, nu, L);
613 hp1D::Riesz<Cmplx> wslform(wsl);
614 Vector<Cmplx> Vwsl(spc, wslform);
615 Vwsl *= Beta(m,l);
616 Vwsm = Vwsm + Vwsl;
617 }
618 Real Norme = Vwsm.l2();
619 for(int l=0; l<=m; l++)
620 Beta(m,l) = Beta(m,l) / Norme;
621 }
622 }
623
624 template<class F>
625 void GrammSchmidt(SparseMatrix<F>& Beta, const Sequence<Vector<F> > VBasis)
626 {
627 const int n = VBasis.size();
628 conceptsAssert(Beta.dimX() == n, Assertion());
629 conceptsAssert(Beta.dimY() == n, Assertion());
630 for (int m=0; m<n;m ++)
631 {
632 Beta(m,m)=1.0;
633 if (m > 0)
634 {
635 // Create projection matrix with coefficients
636 SparseMatrix<F> CoeffsProj(m);
637 for (int l=0; l<m; l++)
638 for (int j=0; j<m; j++)
639 {
640 F ScalProjLJ = VBasis[j] * VBasis[l];
641 for(int i=l; i<m; i++)
642 {
643 CoeffsProj(i,j) = CoeffsProj(i,j) + conj(Beta(i,l)) * ScalProjLJ;
644 }
645 }
646
647 // Create Right hand side
648 Vector<F> Rhs(m);
649 Vector<F> Sol(m);
650 for (int l=0; l<m; l++)
651 {
652 F ScalProjLM = VBasis[m] * VBasis[l];
653 for (int j=l; j<m; j++)
654 {
655 Rhs(j) = Rhs(j) + conj(Beta(j,l)) * ScalProjLM;
656 }
657 }
658
659 // Solve system
660 if (m==1)
661 {
662 // Explicit solve of 1x1 system
663 Beta(1,0) = Rhs(0) / CoeffsProj(0,0);
664 }
665 else
666 {
667 std::unique_ptr<concepts::Operator<F> > solver(nullptr);
668 solver.reset(new concepts::Mumps<F>(CoeffsProj));
669 (*solver)(Rhs, Sol);
670 for (int l=0; l<m; l++)
671 Beta(m,l) = Sol(l);
672 }
673
674 }
675
676 // Renormalization of Betas
677 Vector<F> Vwsm = VBasis[m];
678 for (int l=0; l<m; l++)
679 {
680 Vwsm = Vwsm + VBasis[l] * Beta(m,l);
681 }
682 Real Norme = Vwsm.l2();
683 for(int l=0; l<=m; l++)
684 Beta(m,l) = Beta(m,l) / Norme;
685 }
686 }
687
688 template<class F>
689 void GrammSchmidt(SparseMatrix<F>& Beta, const Sequence<Sequence< Vector<F> > > VBasis)
690 {
691 const int n = VBasis.size();
692 const int dim = VBasis[0].size();
693 conceptsAssert(Beta.dimX() == n, Assertion());
694 conceptsAssert(Beta.dimY() == n, Assertion());
695 for (int m=0; m<n;m ++)
696 {
697 Beta(m,m)=1.0;
698 if (m > 0)
699 {
700 // Create projection matrix with coefficients
701 SparseMatrix<F> CoeffsProj(m);
702 for (int l=0; l<m; l++)
703 for (int j=0; j<m; j++)
704 {
705 F ScalProjLJ = 0;
706 for(int d=0; d<dim; d++)
707 ScalProjLJ += VBasis[j][d] * VBasis[l][d];
708 for(int i=l; i<m; i++)
709 {
710 CoeffsProj(i,j) = CoeffsProj(i,j) + conj(Beta(i,l)) * ScalProjLJ;
711 }
712 }
713
714 // Create Right hand side
715 Vector<F> Rhs(m);
716 Vector<F> Sol(m);
717 for (int l=0; l<m; l++)
718 {
719 F ScalProjLM = 0;
720 for(int d=0; d<dim; d++)
721 ScalProjLM += VBasis[m][d] * VBasis[l][d];
722 for (int j=l; j<m; j++)
723 {
724 Rhs(j) = Rhs(j) + conj(Beta(j,l)) * ScalProjLM;
725 }
726 }
727
728 // Solve system
729 if (m==1)
730 {
731 // Explicit solve of 1x1 system
732 Beta(1,0) = Rhs(0) / CoeffsProj(0,0);
733 }
734 else
735 {
736 std::unique_ptr<concepts::Operator<F> > solver(nullptr);
737 solver.reset(new concepts::Mumps<F>(CoeffsProj));
738 (*solver)(Rhs, Sol);
739 for (int l=0; l<m; l++)
740 Beta(m,l) = Sol(l);
741 }
742
743 }
744
745 // Renormalization of Betas
746 Vector<F> Vwsm = VBasis[m][0];
747 for (int d=1; d<dim; d++)
748 Vwsm = Vwsm + VBasis[m][d];
749 for (int l=0; l<m; l++)
750 for (int d=0; d<dim; d++)
751 {
752 Vwsm = Vwsm + VBasis[l][d] * Beta(m,l);
753 }
754 Real Norme = Vwsm.l2();
755 for(int l=0; l<=m; l++)
756 Beta(m,l) = Beta(m,l) / Norme;
757 }
758 }
759
760#endif
761
762 concepts::Sequence<Cmplx> ReadEigenValuesFromFile(const int NDtN_given, const std::string Filename)
763 {
764 ifstream EigenFile;
765 EigenFile.open(Filename.c_str());
766 int NDtN;
767 EigenFile >> NDtN;
768
769 conceptsAssert(NDtN == NDtN_given, Assertion());
770 concepts::Sequence<Cmplx> Eigenvalues(NDtN);
771
772 Real EigenReal, EigenImag;
773 for (int idx = 0; idx < NDtN ; idx++)
774 {
775 EigenFile >> EigenReal;
776 EigenFile >> EigenImag;
777#if __cplusplus >= 201103L
778 Eigenvalues[idx].real(EigenReal);
779 Eigenvalues[idx].imag(EigenImag);
780#else
781 Eigenvalues[idx].real() = EigenReal;
782 Eigenvalues[idx].imag() = EigenImag;
783#endif
784 // std::cout << "Eigenvalue number " << idx << " is " << Eigenvalues[idx] << std::endl;
785 }
786 EigenFile.close();
787
788 return Eigenvalues;
789 }
790
791}
792
793#endif // __DtNmap2D_visc__hh__
concepts::Cos_n_x
Definition formulas2D.hh:132
concepts::Wsym_x
Definition DtNmap2D_visc.hh:90
concepts::Wsym_x::clone
virtual Wsym_x< dim > * clone() const
Virtual copy constructor.
concepts::Wsym_x::Wsym_x
Wsym_x(const int n=1, const Real omega=1.0, const Real c0=1.0, const Real rho0=1.0, const Real nu=0.1, const Real L=1.0)
Constructor.
concepts::Wsym_x::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
concepts::Wsym_x::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::Wsym_y
Definition DtNmap2D_visc.hh:120
concepts::Wsym_y::clone
virtual Wsym_y< dim > * clone() const
Virtual copy constructor.
concepts::Wsym_y::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
concepts::Wsym_y::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::Wsym_y::Wsym_y
Wsym_y(const int n=1, const Real omega=1.0, const Real c0=1.0, const Real rho0=1.0, const Real nu=0.1, const Real L=1.0)
Constructor.
concepts::Wunsym_x
Definition DtNmap2D_visc.hh:151
concepts::Wunsym_x::Wunsym_x
Wunsym_x(const int n=1, const Real omega=1.0, const Real c0=1.0, const Real rho0=1.0, const Real nu=0.1, const Real L=1.0)
Constructor.
concepts::Wunsym_x::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
concepts::Wunsym_x::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::Wunsym_x::clone
virtual Wunsym_x< dim > * clone() const
Virtual copy constructor.
concepts::Wunsym_y
Definition DtNmap2D_visc.hh:181
concepts::Wunsym_y::Wunsym_y
Wunsym_y(const int n=1, const Real omega=1.0, const Real c0=1.0, const Real rho0=1.0, const Real nu=0.1, const Real L=1.0)
Constructor.
concepts::Wunsym_y::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
concepts::Wunsym_y::clone
virtual Wunsym_y< dim > * clone() const
Virtual copy constructor.
concepts::Wunsym_y::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
frmE_product.hh
conceptsAssert
#define conceptsAssert(cond, exc)
Definition exceptions.hh:394
concepts::Real
double Real
Definition typedefs.hh:39
concepts::makeSet
Set< F > makeSet(uint n, const F &first,...)
Definition set.hh:320
concepts::Cmplx
std::complex< Real > Cmplx
Type for a complex number. It also depends on the setting of Real.
Definition typedefs.hh:42
Generated on Wed Sep 13 2023 21:06:22 for Concepts by doxygen 1.9.8