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 
22 namespace 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::Wunsym_x::clone
virtual Wunsym_x< dim > * clone() const
concepts::Wsym_y::get_lambda
virtual Cmplx get_lambda() const
concepts::Wunsym_x::get_lambda
virtual Cmplx get_lambda() const
concepts::Wunsym_x::ells
const Cmplx ells
Definition: DtNmap2D_visc.hh:176
concepts::Wunsym_y::info
virtual std::ostream & info(std::ostream &os) const
concepts::Cos_n_x
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:132
concepts::Wsym_x::get_lambda
virtual Cmplx get_lambda() const
concepts::visc_ell_fast
Cmplx visc_ell_fast(const Cmplx lambda, const Real nu, const Real omega, const Real rho0)
Definition: DtNmap2D_visc.hh:44
concepts::Wunsym_y
Definition: DtNmap2D_visc.hh:181
concepts::Wunsym_x::omega_
const Real omega_
Definition: DtNmap2D_visc.hh:168
concepts::Wunsym_y::rho0_
const Real rho0_
Definition: DtNmap2D_visc.hh:200
concepts::Wunsym_y::lambda
const Cmplx lambda
Definition: DtNmap2D_visc.hh:204
concepts::Wsym_x::lambda
const Cmplx lambda
Definition: DtNmap2D_visc.hh:113
concepts::Wsym_x::L_
const Real L_
Definition: DtNmap2D_visc.hh:111
concepts::Wsym_y::ellf
const Cmplx ellf
Definition: DtNmap2D_visc.hh:144
concepts::Formula
Interface for a formula.
Definition: lform.hh:18
concepts::sgn
int sgn(const Real d)
Definition: DtNmap2D_visc.hh:37
concepts::Wsym_x::n_
const int n_
Definition: DtNmap2D_visc.hh:106
conceptsAssert
#define conceptsAssert(cond, exc)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:394
concepts::Wsym_y::L_
const Real L_
Definition: DtNmap2D_visc.hh:141
concepts::Wunsym_x::ellf
const Cmplx ellf
Definition: DtNmap2D_visc.hh:175
alglib::conj
alglib::complex conj(const alglib::complex &z)
concepts::SparseMatrix::addInto
void addInto(Matrix< H > &dest, const I fact, const uint rowoffset=0, const uint coloffset=0) const
This matrix is added as block to the given matrix dest.
Definition: sparseMatrix.hh:538
concepts::addExactDtN_X_2Dcos_wp
void addExactDtN_X_2Dcos_wp(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, Vector< Cmplx > &rhs, const G &frm, const Sequence< F > DtNCoeffRhs, const Real L=1.0)
Definition: DtNmap2D_visc.hh:214
concepts::lambda_limit
Cmplx lambda_limit(const Real omega, const Real c0, const int n, const Real L)
Definition: DtNmap2D_visc.hh:82
concepts::cmplx_i
const Cmplx cmplx_i(0, 1)
M_PI
#define M_PI
Definition: typedefs.hh:13
concepts::Wsym_y::nu_
const Real nu_
Definition: DtNmap2D_visc.hh:140
hp1D::Riesz
Linear form on edges in nD.
Definition: linearForm.hh:67
concepts::Wunsym_x::rho0_
const Real rho0_
Definition: DtNmap2D_visc.hh:170
concepts::Wunsym_x::L_
const Real L_
Definition: DtNmap2D_visc.hh:172
concepts::Wsym_y::info
virtual std::ostream & info(std::ostream &os) const
concepts::Wsym_x::ellf
const Cmplx ellf
Definition: DtNmap2D_visc.hh:114
concepts::Wsym_x
Definition: DtNmap2D_visc.hh:90
concepts::Wunsym_x::n_
const int n_
Definition: DtNmap2D_visc.hh:167
concepts::Assertion
Exception class for assertions.
Definition: exceptions.hh:258
concepts::Wsym_y::omega_
const Real omega_
Definition: DtNmap2D_visc.hh:137
concepts::Wsym_y::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Application operator.
concepts::Wsym_x::rho0_
const Real rho0_
Definition: DtNmap2D_visc.hh:109
concepts::Vector::add
virtual Function< F > & add(const Function< F > &fnc, F sc)
concepts::Wunsym_y::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Application operator.
concepts::Wunsym_y::c0_
const Real c0_
Definition: DtNmap2D_visc.hh:199
concepts::Wunsym_x::info
virtual std::ostream & info(std::ostream &os) const
concepts::Matrix
Abstract class for an operator.
Definition: compositions.hh:31
concepts::Wsym_y::c0_
const Real c0_
Definition: DtNmap2D_visc.hh:138
concepts::Wsym_x::c0_
const Real c0_
Definition: DtNmap2D_visc.hh:108
concepts::Wunsym_y::cf
const Cmplx cf
Definition: DtNmap2D_visc.hh:207
concepts::Wsym_y::clone
virtual Wsym_y< dim > * clone() const
concepts::addExactDtN_X_2Dsincos_wp
void addExactDtN_X_2Dsincos_wp(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, Vector< Cmplx > &rhs, const G &frm, const Sequence< F > DtNCoeffRhs, const Real L=1.0)
Definition: DtNmap2D_visc.hh:274
concepts::Wsym_x::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Application operator.
concepts::Sequence
Sequence with operations, output operator, and method of the particular element types.
Definition: sequence.hh:39
concepts::Cmplx
std::complex< Real > Cmplx
Type for a complex number. It also depends on the setting of Real.
Definition: typedefs.hh:39
concepts::Wsym_x::cf
const Cmplx cf
Definition: DtNmap2D_visc.hh:116
concepts::Wunsym_x::c0_
const Real c0_
Definition: DtNmap2D_visc.hh:169
concepts::Wsym_x::ells
const Cmplx ells
Definition: DtNmap2D_visc.hh:115
concepts::Wsym_y::rho0_
const Real rho0_
Definition: DtNmap2D_visc.hh:139
concepts::Sin_n_x
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:162
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::addExactDtN_X_2Dsin_wp
void addExactDtN_X_2Dsin_wp(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, Vector< Cmplx > &rhs, const G &frm, const Sequence< F > DtNCoeffRhs, const Real L=1.0)
Definition: DtNmap2D_visc.hh:233
concepts::Wsym_x::nu_
const Real nu_
Definition: DtNmap2D_visc.hh:110
concepts::visc_ell_slow
Cmplx visc_ell_slow(const Cmplx lambda, const Real nu, const Real omega, const Real rho0, const Real c0)
Definition: DtNmap2D_visc.hh:63
concepts::Wsym_y::lambda
const Cmplx lambda
Definition: DtNmap2D_visc.hh:143
concepts::Wunsym_x::nu_
const Real nu_
Definition: DtNmap2D_visc.hh:171
concepts::Wunsym_y::ellf
const Cmplx ellf
Definition: DtNmap2D_visc.hh:205
concepts::Wunsym_y::omega_
const Real omega_
Definition: DtNmap2D_visc.hh:198
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::SparseMatrix::addIntoT
void addIntoT(Matrix< H > &dest, const I fact, const uint rowoffset=0, const uint coloffset=0) const
The transposed of this matrix is added as block to the given matrix.
Definition: sparseMatrix.hh:567
concepts::Mumps
MUMPS : MUltifrontal Massively Parallel sparse direct Solver.
Definition: mumps.hh:72
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::addExactDtN_X_2Dcossin_wp
void addExactDtN_X_2Dcossin_wp(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, Vector< Cmplx > &rhs, const G &frm, const Sequence< F > DtNCoeffRhs, const Real L=1.0)
Definition: DtNmap2D_visc.hh:251
concepts::Wunsym_x
Definition: DtNmap2D_visc.hh:151
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_x::cf
const Cmplx cf
Definition: DtNmap2D_visc.hh:177
concepts::Wsym_y
Definition: DtNmap2D_visc.hh:120
concepts::Wsym_x::omega_
const Real omega_
Definition: DtNmap2D_visc.hh:107
concepts::sqr
Real sqr(const Real x)
Definition: DtNmap2D_visc.hh:27
frmE_product.hh
concepts::Wsym_y::cf
const Cmplx cf
Definition: DtNmap2D_visc.hh:146
concepts::Wunsym_x::lambda
const Cmplx lambda
Definition: DtNmap2D_visc.hh:174
concepts::ReadEigenValuesFromFile
concepts::Sequence< Cmplx > ReadEigenValuesFromFile(const int NDtN_given, const std::string Filename)
Definition: DtNmap2D_visc.hh:762
concepts::Wsym_y::ells
const Cmplx ells
Definition: DtNmap2D_visc.hh:145
concepts::Wunsym_y::L_
const Real L_
Definition: DtNmap2D_visc.hh:202
concepts::Wsym_x::info
virtual std::ostream & info(std::ostream &os) const
concepts::Wunsym_y::get_lambda
virtual Cmplx get_lambda() const
concepts::Wunsym_y::nu_
const Real nu_
Definition: DtNmap2D_visc.hh:201
concepts::Wsym_x::clone
virtual Wsym_x< dim > * clone() const
concepts::Wunsym_y::n_
const int n_
Definition: DtNmap2D_visc.hh:197
concepts::Wsym_y::n_
const int n_
Definition: DtNmap2D_visc.hh:136
concepts::Wunsym_y::clone
virtual Wunsym_y< dim > * clone() const
concepts::Real
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
concepts::Wunsym_y::ells
const Cmplx ells
Definition: DtNmap2D_visc.hh:206
concepts::Wunsym_x::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Application operator.
concepts
Basic namespace for Concepts-2.
Definition: pml_formula.h:16
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich