DtNmap2D.hh

Go to the documentation of this file.
1 
9 #ifndef dtn_map2d_hh
10 #define dtn_map2d_hh
11 
12 #include "basics/exceptions.hh"
13 #include "toolbox/sequence.hh"
14 #include "space/integral.hh"
15 #include "formula/elementFormulaContainer.hh"
16 #include "formula/formulas2D.hh"
17 #include "formula/bessel.hh"
18 #include "formula/frmE_product.hh"
19 #include "function/vector.hh"
20 #include "operator/sparseMatrix.hh"
21 #include "hp1D/linearForm.hh"
22 
24 #define LaplaceDtNCoeff_D 0
25 #define LaplaceAddDtn_D 0
27 namespace concepts {
28 
29  // ******************************************* getHelmholtzDtNCoeff_Circle2D **
30 
39  Sequence<Cmplx> getHelmholtzDtNCoeff_Circle2D(const Real omega, const Real R,
40  uint N = 0)
41  {
42  // sequence n=0,1,...,N
43  Sequence<int> nSeq = makeRangeSequence(N+1);
44  // (-omega/(2*pi*R) * H'_n(omega*R,n)/H_n(omega*R,n) for n=0,1,...,N
45  return
46  hankel_1_deriv_n(omega*R, nSeq) /
47  hankel_1_n (omega*R, nSeq) * (-omega/(2*M_PI*R));
48  }
49 
50  // ********************************************* getLaplaceDtNCoeff_Circle2D **
51 
61  Sequence<Real> getLaplaceDtNCoeff_Circle2D(const Real R, uint N = 0)
62  {
63  DEBUGL(LaplaceDtNCoeff_D, "R = " << R << ", N = " << N);
64  // sequence n=0,1,...,N
65  Sequence<int> nSeq = makeRangeSequence(N+1);
66  // (R^(-n))'/R^(-n) = (-n R^(-(n+1)) / R^(-n) = -n / R
67  Sequence<Real> DtNCoeff = nSeq / (-R);
68  // log(R)' / log(R) = 1/R/log(R)
69  DtNCoeff[0] = 1./R/log(R);
70  // Divide by factor coming from Fourier transform and dphi = 1/R ds.
71  // Sign change as it is for -Delta where after integration by parts
72  // we have -int_partialOmega dnu v ds.
73  DtNCoeff = DtNCoeff / (-2*M_PI*R);
74 
75  DEBUGL(LaplaceDtNCoeff_D, "DtNCoeff = " << DtNCoeff);
76 
77  return DtNCoeff;
78  }
79 
80  // ******************************************* getLaplaceDtNCoeff_Straight2D **
81 
129  Sequence<Real> getLaplaceDtNCoeff_Straight2D(const Real L, uint N = 0)
130  {
131  DEBUGL(LaplaceDtNCoeff_D, "L = " << L << ", N = " << N);
132  // sequence n=0,1,...,N
133  Sequence<int> nSeq = makeRangeSequence(N+1);
134  // DtN-coefficients: 2 * pi * n / L^2
135  Sequence<Real> DtNCoeff = nSeq * (2. * M_PI / (L*L));
136 
137  DEBUGL(LaplaceDtNCoeff_D, "DtNCoeff = " << DtNCoeff);
138 
139  return DtNCoeff;
140  }
141 
142  // ****************************************** getHelmholtzDtNCoeff_Straight2D **
143 
190  Sequence<Cmplx> getHelmholtzDtNCoeff_Straight2D(const Real omega,
191  const Real L, uint N = 0)
192  {
193  DEBUGL(LaplaceDtNCoeff_D, "omega = " << omega << ", L = " << L
194  << ", N = " << N);
195  // sequence n=0,1,...,N
196  Sequence<Cmplx> DtNCoeff(N+1);
197  Cmplx I(0,1);
198  // DtN-coefficients: 2 * pi * n / L^2
199 
200  for (uint n=0; n <= N; n++)
201  {
202  if (n*M_PI < omega*L)
203  DtNCoeff[n] = I*sqrt(omega*omega-(n*M_PI/L)*(n*M_PI/L));
204  else
205  DtNCoeff[n] = -sqrt((n*M_PI/L)*(n*M_PI/L)-omega*omega);
206  }
207 
208  DEBUGL(LaplaceDtNCoeff_D, "DtNCoeff = " << DtNCoeff);
209 
210  return DtNCoeff;
211  }
212 
213  // **************************************************** addExactDtN_Circle2D **
214 
215 
223  void addExactDtN_Circle2D(Matrix<Real>& dest, const SpaceOnCells<Real>& spc,
224  const Sequence<Real> DtNCoeff)
225  {
226  // n = 0
227  {
228  ConstFormula<Real> frm(1.0);
229  hp1D::Riesz<Real> vLform(frm);
230  Vector<Real> vVector(spc, vLform);
231  // rank-1-product
232  SparseMatrix<Real> DtNmatrix(spc,vVector,vVector);
233  // add it to system matrix
234  DtNmatrix.addInto(dest, DtNCoeff[0]);
235  }
236  // n >= 1
237  for (uint n = 1; n < DtNCoeff.size(); ++n) {
238  // part with cos(n phi)
239  Cos_n_phi cos_nphi(n);
240  hp1D::Riesz<Real> vLform(cos_nphi);
241  Vector<Real> vVector(spc, vLform);
242  // rank-1-product
243  SparseMatrix<Real> DtNmatrix(spc,vVector,vVector);
244  // add it to system matrix
245  DtNmatrix.addInto(dest, DtNCoeff[n]*2.0);
246  // part with sin(n phi)
247  Sin_n_phi sin_nphi(n);
248  hp1D::Riesz<Real> uLform(sin_nphi);
249  Vector<Real> uVector(spc, uLform);
250  // rank-1-product
251  DtNmatrix = SparseMatrix<Real>(spc,uVector,uVector);
252  // add it to system matrix
253  DtNmatrix.addInto(dest, DtNCoeff[n]*2.0);
254  }
255  }
256 
263  template<class F>
264  void addExactDtN_Circle2D(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
265  const Sequence<F> DtNCoeff)
266  {
267  for (uint n = 0; n < DtNCoeff.size(); ++n) {
268  // v-integral
269  Exp_i_n_phi exp_p_inphi(+n);
270  hp1D::Riesz<Cmplx> vLform(exp_p_inphi);
271  Vector<Cmplx> vVector(spc, vLform);
272  // u-integral
273  Exp_i_n_phi exp_m_inphi(-n);
274  hp1D::Riesz<Cmplx> uLform(exp_m_inphi);
275  Vector<Cmplx> uVector(spc, uLform);
276  // rank-1-product of v- and u-integral
277  SparseMatrix<Cmplx> DtNmatrix(spc,vVector,uVector);
278  // add it to system matrix
279  DtNmatrix.addInto(dest, DtNCoeff[n]);
280  // add contribution for -n (test and trial functions interchanged)
281  if (n > 0)
282  DtNmatrix.addIntoT(dest, DtNCoeff[n]);
283  }
284  }
285 
286  void addExactDtN_Circle2D(Matrix<Real>& dest, const SpaceOnCells<Real>& spc,
287  const Sequence<Real> DtNCoeff, Vector<Real>& rhs,
288  const ElementFormula<Real>& frm,
289  const Sequence<Real> DtNCoeffRhs)
290  {
291  conceptsAssert(DtNCoeff.size() == DtNCoeffRhs.size(), Assertion());
292  // n = 0
293  {
294  ConstFormula<Real> one(1.0);
295  hp1D::Riesz<Real> vLform(one);
296  Vector<Real> vVector(spc, vLform);
297  // rank-1-product
298  SparseMatrix<Real> DtNmatrix(spc,vVector,vVector);
299  // add it to system matrix
300  DtNmatrix.addInto(dest, DtNCoeff[0]);
301  // integral of frm along the circular boundary
302  Real integral = integrate(spc,frm);
303  // contribution for vector of r.h.s.
304  rhs = rhs + vVector*(DtNCoeffRhs[0]*integral);
305  }
306  for (uint n = 1; n < DtNCoeff.size(); ++n) {
307  // part with cos(n phi)
308  Cos_n_phi cos_nphi(n);
309  hp1D::Riesz<Real> vLform(cos_nphi);
310  Vector<Real> vVector(spc, vLform);
311  // rank-1-product
312  SparseMatrix<Real> DtNmatrix(spc,vVector,vVector);
313  // add it to system matrix
314  DtNmatrix.addInto(dest, DtNCoeff[n]*2.0);
315  // integral of frm * cos(n phi) along the circular boundary
316  Real integral = integrate(spc,frm*cos_nphi);
317  // contribution for vector of r.h.s.
318  rhs = rhs + vVector*(DtNCoeffRhs[n]*2.0*integral);
319 
320  // part with sin(n phi)
321  Sin_n_phi sin_nphi(n);
322  hp1D::Riesz<Real> uLform(sin_nphi);
323  Vector<Real> uVector(spc, uLform);
324  // rank-1-product
325  DtNmatrix = SparseMatrix<Real>(spc,uVector,uVector);
326  // add it to system matrix
327  DtNmatrix.addInto(dest, DtNCoeff[n]*2.0);
328  // integral of frm * cos(n phi) along the circular boundary
329  integral = integrate(spc,frm*sin_nphi);
330  // contribution for vector of r.h.s.
331  rhs = rhs + vVector*(DtNCoeffRhs[n]*2.0*integral);
332  }
333  }
334 
335  template<class F, class G>
336  void addExactDtN_Circle2D(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
337  const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
338  const ElementFormula<G>& frm,
339  const Sequence<F> DtNCoeffRhs)
340  {
341  conceptsAssert(DtNCoeff.size() == DtNCoeffRhs.size(), Assertion());
342  for (uint n = 0; n < DtNCoeff.size(); ++n) {
343  // v-integral
344  Exp_i_n_phi exp_p_inphi(+n);
345  hp1D::Riesz<Cmplx> vLform(exp_p_inphi);
346  Vector<Cmplx> vVector(spc, vLform);
347  // u-integral
348  Exp_i_n_phi exp_m_inphi(-n);
349  hp1D::Riesz<Cmplx> uLform(exp_m_inphi);
350  Vector<Cmplx> uVector(spc, uLform);
351  // rank-1-product of v- and u-integral
352  SparseMatrix<Cmplx> DtNmatrix(spc,vVector,uVector);
353  // add it to system matrix
354  DtNmatrix.addInto(dest, DtNCoeff[n]);
355  // add contribution for -n (test and trial functions interchanged)
356  if (n > 0)
357  DtNmatrix.addIntoT(dest, DtNCoeff[n]);
358 
359  // integral of frm * exp_m_inphi along the circular boundary
360  Cmplx integral = integrate(spc,frm*exp_m_inphi);
361  // contribution for vector of r.h.s.
362  rhs = rhs + vVector*(DtNCoeffRhs[n]*integral);
363  // if (n > 0) also add contribution of -n
364  if (n > 0) {
365  Cmplx integral2 = integrate(spc,frm*exp_p_inphi);
366  rhs = rhs + uVector*(DtNCoeffRhs[n]*integral2);
367  }
368  }
369  }
370 
371  // ******************************************************** addExactDtN_X_2D **
372 
373  void addExactDtN_X_2D(Matrix<Real>& dest, const SpaceOnCells<Real>& spc,
374  const Sequence<Real> DtNCoeff, const Real L = 1.0)
375  {
376  DEBUGL(LaplaceAddDtn_D, "L = " << L);
377  conceptsAssert(L > 0, Assertion());
378  // concepts::MatfileOutput mo("matrices.mat");
379  // n = 0
380  {
381  ConstFormula<Real> frm(1.0);
382  hp1D::Riesz<Real> vLform(frm);
383  Vector<Real> vVector(spc, vLform);
384  // rank-1-product
385  SparseMatrix<Real> DtNmatrix(spc,vVector,vVector);
386  // mo.add(DtNmatrix, "L0");
387  // add it to system matrix
388  DtNmatrix.addInto(dest, DtNCoeff[0]);
389  }
390  // n >= 1
391  for (uint n = 1; n < DtNCoeff.size(); ++n) {
392  // part with cos(2 * pi * n * x / L)
393  Cos_n_x cos_nx(n, L);
394  DEBUGL(LaplaceAddDtn_D, cos_nx);
395  hp1D::Riesz<Real> vLform(cos_nx);
396  Vector<Real> vVector(spc, vLform);
397  // rank-1-product
398  SparseMatrix<Real> DtNmatrix(spc, vVector, vVector);
399  std::stringstream s;
400  s << "LC" << n;
401  // mo.add(DtNmatrix, s.str());
402  s.str(""); s << "lc" << n;
403  // mo.add(vVector, s.str());
404  // add it to system matrix
405  DtNmatrix.addInto(dest, DtNCoeff[n]*2.0);
406  // part with sin(2 * pi * n * x / L)
407  Sin_n_x sin_nx(n, L);
408  DEBUGL(LaplaceAddDtn_D, sin_nx);
409  hp1D::Riesz<Real> uLform(sin_nx);
410  Vector<Real> uVector(spc, uLform);
411  // rank-1-product
412  DtNmatrix = SparseMatrix<Real>(spc, uVector, uVector);
413  s.str(""); s << "LS" << n;
414  // mo.add(DtNmatrix, s.str());
415  s.str(""); s << "ls" << n;
416  // mo.add(uVector, s.str());
417  // add it to system matrix
418  DtNmatrix.addInto(dest, DtNCoeff[n]*2.0);
419  }
420  }
421 
422  template<class F>
423  void addExactDtN_X_2D(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
424  const Sequence<F> DtNCoeff, const Real L = 1.0)
425  {
426  conceptsAssert(L > 0, Assertion());
427  for (uint n = 0; n < DtNCoeff.size(); ++n) {
428  // v-integral
429  Exp_i_n_x exp_p_inx(+n, L);
430  hp1D::Riesz<Cmplx> vLform(exp_p_inx);
431  Vector<Cmplx> vVector(spc, vLform);
432  // u-integral
433  Exp_i_n_x exp_m_inx(-n, L);
434  hp1D::Riesz<Cmplx> uLform(exp_m_inx);
435  Vector<Cmplx> uVector(spc, uLform);
436  // rank-1-product of v- and u-integral
437  SparseMatrix<Cmplx> DtNmatrix(spc,vVector,uVector);
438  // add it to system matrix
439  DtNmatrix.addInto(dest, DtNCoeff[n]);
440  // add contribution for -n (test and trial functions interchanged)
441  if (n > 0)
442  DtNmatrix.addIntoT(dest, DtNCoeff[n]);
443  }
444  }
445 
446  template<class F, class G>
447  void addExactDtN_X_2D(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
448  const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
449  const ElementFormula<G>& frm, const Sequence<F> DtNCoeffRhs,
450  const Real L = 1.0)
451  {
452  conceptsAssert(L > 0, Assertion());
453  conceptsAssert(DtNCoeff.size() == DtNCoeffRhs.size(), Assertion());
454  for (uint n = 0; n < DtNCoeff.size(); ++n) {
455  // v-integral
456  Exp_i_n_x exp_p_inx(+n, L);
457  hp1D::Riesz<Cmplx> vLform(exp_p_inx);
458  Vector<Cmplx> vVector(spc, vLform);
459  // u-integral
460  Exp_i_n_x exp_m_inx(-n, L);
461  hp1D::Riesz<Cmplx> uLform(exp_m_inx);
462  Vector<Cmplx> uVector(spc, uLform);
463  // rank-1-product of v- and u-integral
464  SparseMatrix<Cmplx> DtNmatrix(spc,vVector,uVector);
465  // add it to system matrix
466  DtNmatrix.addInto(dest, DtNCoeff[n]);
467  // add contribution for -n (test and trial functions interchanged)
468  if (n > 0)
469  DtNmatrix.addIntoT(dest, DtNCoeff[n]);
470 
471  // Integral of frm * exp_m_inx along the boundary
472  Cmplx integral = integrate(spc, frm*exp_m_inx);
473  // Contribution for vector of r.h.s.
474  rhs = rhs + vVector*(DtNCoeffRhs[n]*integral);
475  // add contribution for -n
476  if (n > 0)
477  {
478  Cmplx integral2 = integrate(spc, frm*exp_p_inx);
479  rhs = rhs + uVector*(DtNCoeffRhs[n]*integral2);
480  }
481 
482  }
483  }
484 
485  // ******************************************************** addExactDtN_Y_2D **
486 
487  template<class F>
488  void addExactDtN_Y_2D(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
489  const Sequence<F> DtNCoeff, const Real L = 1.0)
490  {
491  conceptsAssert(L > 0, Assertion());
492  for (uint n = 0; n < DtNCoeff.size(); ++n) {
493  // v-integral
494  Exp_i_n_y exp_p_iny(+n, L);
495  hp1D::Riesz<Cmplx> vLform(exp_p_iny);
496  Vector<Cmplx> vVector(spc, vLform);
497  // u-integral
498  Exp_i_n_y exp_m_iny(-n, L);
499  hp1D::Riesz<Cmplx> uLform(exp_m_iny);
500  Vector<Cmplx> uVector(spc, uLform);
501  // rank-1-product of v- and u-integral
502  SparseMatrix<Cmplx> DtNmatrix(spc,vVector,uVector);
503  // add it to system matrix
504  DtNmatrix.addInto(dest, DtNCoeff[n]);
505  // add contribution for -n (test and trial functions interchanged)
506  if (n > 0)
507  DtNmatrix.addIntoT(dest, DtNCoeff[n]);
508  }
509  }
510 
511  template<class F, class G>
512  void addExactDtN_Y_2D(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
513  const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
514  const ElementFormula<G>& frm, const Sequence<F> DtNCoeffRhs,
515  const Real L = 1.0)
516  {
517  conceptsAssert(L > 0, Assertion());
518  conceptsAssert(DtNCoeff.size() == DtNCoeffRhs.size(), Assertion());
519  for (uint n = 0; n < DtNCoeff.size(); ++n) {
520  // v-integral
521  Exp_i_n_y exp_p_iny(+n, L);
522  hp1D::Riesz<Cmplx> vLform(exp_p_iny);
523  Vector<Cmplx> vVector(spc, vLform);
524  // u-integral
525  Exp_i_n_y exp_m_iny(-n, L);
526  hp1D::Riesz<Cmplx> uLform(exp_m_iny);
527  Vector<Cmplx> uVector(spc, uLform);
528  // rank-1-product of v- and u-integral
529  SparseMatrix<Cmplx> DtNmatrix(spc,vVector,uVector);
530  // add it to system matrix
531  DtNmatrix.addInto(dest, DtNCoeff[n]);
532  // add contribution for -n (test and trial functions interchanged)
533  if (n > 0)
534  DtNmatrix.addIntoT(dest, DtNCoeff[n]);
535 
536  // Integral of frm * exp_m_iny along the boundary
537  Cmplx integral = integrate(spc, frm*exp_m_iny);
538  // Contribution for vector of r.h.s.
539  rhs = rhs + vVector*(DtNCoeffRhs[n]*integral);
540  // add contribution for -n
541  if (n > 0)
542  {
543  Cmplx integral2 = integrate(spc, frm*exp_p_iny);
544  rhs = rhs + uVector*(DtNCoeffRhs[n]*integral2);
545  }
546 
547  }
548  }
549 
550  // ****************************************************** addExactDtN_Y_2Dsym **
551 
552  template<class F>
553  void addExactDtN_Y_2Dsym(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
554  const Sequence<F> DtNCoeff, const Real L = 1.0)
555  {
556  conceptsAssert(L > 0, Assertion());
557  for (uint n = 0; n < DtNCoeff.size(); ++n) {
558  Cos_n_y cos_ny(+n, L);
559  hp1D::Riesz<Real> cLform(cos_ny);
560  Vector<Real> cVector(spc, cLform);
561  SparseMatrix<Real> cDtNmatrix(spc,cVector,cVector);
562  Cmplx factor = 2.0 / L;
563  if (n == 0) factor = 1.0 / L;
564  cDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
565  if (n > 0)
566  {
567  Sin_n_y sin_ny(+n, L);
568  hp1D::Riesz<Real> sLform(sin_ny);
569  Vector<Real> sVector(spc, sLform);
570  SparseMatrix<Real> sDtNmatrix(spc,sVector,sVector);
571  }
572  }
573  }
574 
575  // ****************************************************** addExactDtN_Y_2Dcos **
576 
577  template<class F>
578  void addExactDtN_Y_2Dcos(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
579  const Sequence<F> DtNCoeff, const Real L = 1.0)
580  {
581  conceptsAssert(L > 0, Assertion());
582  for (uint n = 0; n < DtNCoeff.size(); ++n) {
583  Cos_n_y cos_ny(+n, L);
584  hp1D::Riesz<Real> cLform(cos_ny);
585  Vector<Real> cVector(spc, cLform);
586  SparseMatrix<Real> cDtNmatrix(spc,cVector,cVector);
587  Cmplx factor = 2.0 / L;
588  if (n == 0) factor = 1.0 / L;
589  cDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
590  }
591  }
592 
593  // ****************************************************** addExactDtN_Y_2Dsin **
594 
595  template<class F>
596  void addExactDtN_Y_2Dsin(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
597  const Sequence<F> DtNCoeff, const Real L = 1.0)
598  {
599  conceptsAssert(L > 0, Assertion());
600  for (uint n = 1; n < DtNCoeff.size(); ++n) {
601  Sin_n_y sin_ny(+n, L);
602  hp1D::Riesz<Real> sLform(sin_ny);
603  Vector<Real> sVector(spc, sLform);
604  SparseMatrix<Real> sDtNmatrix(spc,sVector,sVector);
605  Cmplx factor = 2.0 / L;
606  sDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
607  }
608  }
609 
610  // **************************************************** addExactDtN_Y_2Dunsym **
611 
612  template<class F>
613  void addExactDtN_Y_2Dunsym(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
614  const Sequence<F> DtNCoeff, const Real L = 1.0)
615  {
616  conceptsAssert(L > 0, Assertion());
617  for (uint n = 0; n < DtNCoeff.size(); ++n) {
618  Cos_n_y cos_ny(+n+1, L);
619  Sin_n_y sin_ny(+n+1, L);
620  hp1D::Riesz<Real> cLform(cos_ny);
621  Vector<Real> cVector(spc, cLform);
622  hp1D::Riesz<Real> sLform(sin_ny);
623  Vector<Real> sVector(spc, sLform);
624  // cos part over test fuction, sin part over unknown function
625  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
626  Cmplx factor = 2.0 / L;
627  csDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
628  csDtNmatrix.addIntoT(dest, -DtNCoeff[n]*factor);
629  }
630  }
631 
632  // *************************************************** addExactDtN_Y_2Dcossin **
633 
634  template<class F>
635  void addExactDtN_Y_2Dcossin(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
636  const Sequence<F> DtNCoeff, const Real L = 1.0)
637  {
638  conceptsAssert(L > 0, Assertion());
639  for (uint n = 0; n < DtNCoeff.size(); ++n) {
640  Cos_n_y cos_ny(+n+1, L);
641  Sin_n_y sin_ny(+n+1, L);
642  hp1D::Riesz<Real> cLform(cos_ny);
643  Vector<Real> cVector(spc, cLform);
644  hp1D::Riesz<Real> sLform(sin_ny);
645  Vector<Real> sVector(spc, sLform);
646  // cos part over test fuction, sin part over unknown function
647  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
648  Cmplx factor = 2.0 / L;
649  csDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
650  }
651  }
652 
653  // *************************************************** addExactDtN_Y_2Dsincos **
654 
655  template<class F>
656  void addExactDtN_Y_2Dsincos(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
657  const Sequence<F> DtNCoeff, const Real L = 1.0)
658  {
659  conceptsAssert(L > 0, Assertion());
660  for (uint n = 0; n < DtNCoeff.size(); ++n) {
661  Cos_n_y cos_ny(+n+1, L);
662  Sin_n_y sin_ny(+n+1, L);
663  hp1D::Riesz<Real> cLform(cos_ny);
664  Vector<Real> cVector(spc, cLform);
665  hp1D::Riesz<Real> sLform(sin_ny);
666  Vector<Real> sVector(spc, sLform);
667  // cos part over test fuction, sin part over unknown function
668  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
669  Cmplx factor = 2.0 / L;
670  csDtNmatrix.addIntoT(dest, DtNCoeff[n]*factor);
671  }
672  }
673 
674  template<class F>
675  void addExactDtN_X_2Dsym(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
676  const Sequence<F> DtNCoeff, const Real L = 1.0)
677  {
678  conceptsAssert(L > 0, Assertion());
679  for (uint n = 0; n < DtNCoeff.size(); ++n) {
680  Cos_n_x cos_nx(+n, L);
681  hp1D::Riesz<Real> cLform(cos_nx);
682  Vector<Real> cVector(spc, cLform);
683  SparseMatrix<Real> cDtNmatrix(spc,cVector,cVector);
684  Cmplx factor = 2.0 / L;
685  if (n == 0) factor = 1.0 / L;
686  cDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
687  if (n > 0)
688  {
689  Sin_n_x sin_nx(+n, L);
690  hp1D::Riesz<Real> sLform(sin_nx);
691  Vector<Real> sVector(spc, sLform);
692  SparseMatrix<Real> sDtNmatrix(spc,sVector,sVector);
693  }
694  }
695  }
696 
697  // ****************************************************** addExactDtN_X_2Dcos **
698 
699  template<class F>
700  void addExactDtN_X_2Dcos(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
701  const Sequence<F> DtNCoeff, const Real L = 1.0)
702  {
703  conceptsAssert(L > 0, Assertion());
704  for (uint n = 0; n < DtNCoeff.size(); ++n) {
705  Cos_n_x cos_nx(+n, L);
706  hp1D::Riesz<Real> cLform(cos_nx);
707  Vector<Real> cVector(spc, cLform);
708  SparseMatrix<Real> cDtNmatrix(spc,cVector,cVector);
709  Cmplx factor = 2.0 / L;
710  if (n == 0) factor = 1.0 / L;
711  cDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
712  }
713  }
714 
715  template<class F, class G>
716  void addExactDtN_X_2Dcos(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
717  const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
718  const ElementFormula<G>& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
719  {
720  conceptsAssert(L > 0, Assertion());
721  for (uint n = 0; n < DtNCoeff.size(); ++n) {
722  Cos_n_x cos_nx(+n, L);
723  hp1D::Riesz<Real> cLform(cos_nx);
724  Vector<Real> cVector(spc, cLform);
725  SparseMatrix<Real> cDtNmatrix(spc,cVector,cVector);
726  Cmplx factor = 2.0 / L;
727  if (n == 0) factor = 1.0 / L;
728  cDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
729  Cmplx integral = integrate(spc, frm*cos_nx);
730  rhs = rhs + cVector*(DtNCoeffRhs[n]*integral);
731  }
732  }
733 
734  // ****************************************************** addExactDtN_X_2Dsin **
735 
736  template<class F>
737  void addExactDtN_X_2Dsin(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
738  const Sequence<F> DtNCoeff, const Real L = 1.0)
739  {
740  conceptsAssert(L > 0, Assertion());
741  for (uint n = 1; n < DtNCoeff.size(); ++n) {
742  Sin_n_x sin_nx(+n, L);
743  hp1D::Riesz<Real> sLform(sin_nx);
744  Vector<Real> sVector(spc, sLform);
745  SparseMatrix<Real> sDtNmatrix(spc,sVector,sVector);
746  Cmplx factor = 2.0 / L;
747  sDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
748  }
749  }
750 
751  template<class F, class G>
752  void addExactDtN_X_2Dsin(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
753  const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
754  const ElementFormula<G>& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
755  {
756  conceptsAssert(L > 0, Assertion());
757  for (uint n = 1; n < DtNCoeff.size(); ++n) {
758  Sin_n_x sin_nx(+n, L);
759  hp1D::Riesz<Real> sLform(sin_nx);
760  Vector<Real> sVector(spc, sLform);
761  SparseMatrix<Real> sDtNmatrix(spc,sVector,sVector);
762  Cmplx factor = 2.0 / L;
763  sDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
764  Cmplx integral = integrate(spc, frm*sin_nx);
765  rhs = rhs + sVector*(DtNCoeffRhs[n]*integral);
766  }
767  }
768 
769  // **************************************************** addExactDtN_Y_2Dunsym **
770 
771  template<class F>
772  void addExactDtN_X_2Dunsym(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
773  const Sequence<F> DtNCoeff, const Real L = 1.0)
774  {
775  conceptsAssert(L > 0, Assertion());
776  for (uint n = 0; n < DtNCoeff.size(); ++n) {
777  Cos_n_x cos_nx(+n+1, L);
778  Sin_n_x sin_nx(+n+1, L);
779  hp1D::Riesz<Real> cLform(cos_nx);
780  Vector<Real> cVector(spc, cLform);
781  hp1D::Riesz<Real> sLform(sin_nx);
782  Vector<Real> sVector(spc, sLform);
783  // cos part over test fuction, sin part over unknown function
784  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
785  Cmplx factor = 2.0 / L;
786  csDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
787  csDtNmatrix.addIntoT(dest, -DtNCoeff[n]*factor);
788  }
789  }
790 
791  // *************************************************** addExactDtN_Y_2Dcossin **
792 
793  template<class F>
794  void addExactDtN_X_2Dcossin(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
795  const Sequence<F> DtNCoeff, const Real L = 1.0)
796  {
797  conceptsAssert(L > 0, Assertion());
798  for (uint n = 0; n < DtNCoeff.size(); ++n) {
799  Cos_n_x cos_nx(+n+1, L);
800  Sin_n_x sin_nx(+n+1, L);
801  hp1D::Riesz<Real> cLform(cos_nx);
802  Vector<Real> cVector(spc, cLform);
803  hp1D::Riesz<Real> sLform(sin_nx);
804  Vector<Real> sVector(spc, sLform);
805  // cos part over test fuction, sin part over unknown function
806  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
807  Cmplx factor = 2.0 / L;
808  csDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
809  }
810  }
811 
812  template<class F, class G>
813  void addExactDtN_X_2Dcossin(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
814  const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
815  const ElementFormula<G>& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
816  {
817  conceptsAssert(L > 0, Assertion());
818  for (uint n = 0; n < DtNCoeff.size(); ++n) {
819  Cos_n_x cos_nx(+n+1, L);
820  Sin_n_x sin_nx(+n+1, L);
821  hp1D::Riesz<Real> cLform(cos_nx);
822  Vector<Real> cVector(spc, cLform);
823  hp1D::Riesz<Real> sLform(sin_nx);
824  Vector<Real> sVector(spc, sLform);
825  // cos part over test fuction, sin part over unknown function
826  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
827  Cmplx factor = 2.0 / L;
828  csDtNmatrix.addInto(dest, DtNCoeff[n]*factor);
829  Cmplx integral = integrate(spc, frm*cos_nx);
830  rhs = rhs + sVector*(DtNCoeffRhs[n]*integral);
831  }
832  }
833 
834  // *************************************************** addExactDtN_Y_2Dsincos **
835 
836  template<class F>
837  void addExactDtN_X_2Dsincos(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
838  const Sequence<F> DtNCoeff, const Real L = 1.0)
839  {
840  conceptsAssert(L > 0, Assertion());
841  for (uint n = 0; n < DtNCoeff.size(); ++n) {
842  Cos_n_x cos_nx(+n+1, L);
843  Sin_n_x sin_nx(+n+1, L);
844  hp1D::Riesz<Real> cLform(cos_nx);
845  Vector<Real> cVector(spc, cLform);
846  hp1D::Riesz<Real> sLform(sin_nx);
847  Vector<Real> sVector(spc, sLform);
848  // cos part over test fuction, sin part over unknown function
849  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
850  Cmplx factor = 2.0 / L;
851  csDtNmatrix.addIntoT(dest, DtNCoeff[n]*factor);
852  }
853  }
854 
855 
856  template<class F, class G>
857  void addExactDtN_X_2Dsincos(Matrix<Cmplx>& dest, const SpaceOnCells<Real>& spc,
858  const Sequence<F> DtNCoeff, Vector<Cmplx>& rhs,
859  const ElementFormula<G>& frm, const Sequence<F> DtNCoeffRhs, const Real L = 1.0)
860  {
861  conceptsAssert(L > 0, Assertion());
862  for (uint n = 0; n < DtNCoeff.size(); ++n) {
863  Cos_n_x cos_nx(+n+1, L);
864  Sin_n_x sin_nx(+n+1, L);
865  hp1D::Riesz<Real> cLform(cos_nx);
866  Vector<Real> cVector(spc, cLform);
867  hp1D::Riesz<Real> sLform(sin_nx);
868  Vector<Real> sVector(spc, sLform);
869  // cos part over test fuction, sin part over unknown function
870  SparseMatrix<Real> csDtNmatrix(spc,cVector,sVector);
871  Cmplx factor = 2.0 / L;
872  csDtNmatrix.addIntoT(dest, DtNCoeff[n]*factor);
873  Cmplx integral = integrate(spc, frm*sin_nx);
874  rhs = rhs + cVector*(DtNCoeffRhs[n]*integral);
875  }
876  }
877 
878 } // namespace concepts
879 
880 #endif // dtn_map2d_hh
concepts::integrate
Real integrate(const Element< G > &elm)
Returns the area of the cell belonging to the element elm.
Definition: integral.hh:43
concepts::Cos_n_phi
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:49
concepts::hankel_1_n
Cmplx hankel_1_n(const Real x, const int n)
Evaluates the Hankel function .
concepts::Cos_n_x
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:132
concepts::Sin_n_phi
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:76
concepts::addExactDtN_X_2Dsincos
void addExactDtN_X_2Dsincos(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:837
concepts::addExactDtN_Y_2Dsincos
void addExactDtN_Y_2Dsincos(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:656
LaplaceDtNCoeff_D
#define LaplaceDtNCoeff_D
Debug flag for construction of the DtN coefficients.
Definition: DtNmap2D.hh:24
concepts::ConstFormula
Class for a constant formula.
Definition: constFormula.hh:26
concepts::addExactDtN_Y_2D
void addExactDtN_Y_2D(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:488
concepts::getLaplaceDtNCoeff_Circle2D
Sequence< Real > getLaplaceDtNCoeff_Circle2D(const Real R, uint N=0)
Returns the coefficients for a non-local DtN map for the Laplace operator for a circular boundary of ...
Definition: DtNmap2D.hh:61
concepts::getLaplaceDtNCoeff_Straight2D
Sequence< Real > getLaplaceDtNCoeff_Straight2D(const Real L, uint N=0)
Returns the coefficients for a non-local DtN map for the Laplace operator for a straight boundary of ...
Definition: DtNmap2D.hh:129
concepts::addExactDtN_X_2D
void addExactDtN_X_2D(Matrix< Real > &dest, const SpaceOnCells< Real > &spc, const Sequence< Real > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:373
conceptsAssert
#define conceptsAssert(cond, exc)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:394
concepts::addExactDtN_X_2Dcos
void addExactDtN_X_2Dcos(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:700
concepts::Exp_i_n_phi
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:22
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::Exp_i_n_x
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:103
DEBUGL
#define DEBUGL(doit, msg)
M_PI
#define M_PI
Definition: typedefs.hh:13
hp1D::Riesz
Linear form on edges in nD.
Definition: linearForm.hh:67
concepts::Sin_n_y
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:250
concepts::Cos_n_y
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:220
concepts::Assertion
Exception class for assertions.
Definition: exceptions.hh:258
concepts::addExactDtN_X_2Dcossin
void addExactDtN_X_2Dcossin(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:794
concepts::addExactDtN_Circle2D
void addExactDtN_Circle2D(Matrix< Real > &dest, const SpaceOnCells< Real > &spc, const Sequence< Real > DtNCoeff)
Add DtN contribution for a circular boundary.
Definition: DtNmap2D.hh:223
concepts::hankel_1_deriv_n
Cmplx hankel_1_deriv_n(const Real x, const int n)
Evaluates the derivative of the Hankel function .
concepts::addExactDtN_X_2Dunsym
void addExactDtN_X_2Dunsym(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:772
concepts::addExactDtN_Y_2Dsym
void addExactDtN_Y_2Dsym(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:553
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::Sin_n_x
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:162
concepts::addExactDtN_Y_2Dcossin
void addExactDtN_Y_2Dcossin(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:635
concepts::Exp_i_n_y
Class for evaluating for points in cartesian coordinates.
Definition: formulas2D.hh:191
LaplaceAddDtn_D
#define LaplaceAddDtn_D
Debug flag for construction of the DtN operators.
Definition: DtNmap2D.hh:26
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::getHelmholtzDtNCoeff_Circle2D
Sequence< Cmplx > getHelmholtzDtNCoeff_Circle2D(const Real omega, const Real R, uint N=0)
Returns the coefficients for a non-local DtN map for the Helmholtz operator with frequency omega for ...
Definition: DtNmap2D.hh:39
concepts::addExactDtN_Y_2Dunsym
void addExactDtN_Y_2Dunsym(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:613
concepts::addExactDtN_Y_2Dsin
void addExactDtN_Y_2Dsin(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:596
concepts::getHelmholtzDtNCoeff_Straight2D
Sequence< Cmplx > getHelmholtzDtNCoeff_Straight2D(const Real omega, const Real L, uint N=0)
Returns the coefficients for a non-local DtN map for the Helmholtz operator with frequency omega for ...
Definition: DtNmap2D.hh:190
concepts::makeRangeSequence
Sequence< int > makeRangeSequence(uint n)
Returns the sequence 0,1,...,n-1.
frmE_product.hh
concepts::addExactDtN_X_2Dsin
void addExactDtN_X_2Dsin(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:737
concepts::addExactDtN_Y_2Dcos
void addExactDtN_Y_2Dcos(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:578
concepts::Real
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
concepts::addExactDtN_X_2Dsym
void addExactDtN_X_2Dsym(Matrix< Cmplx > &dest, const SpaceOnCells< Real > &spc, const Sequence< F > DtNCoeff, const Real L=1.0)
Definition: DtNmap2D.hh:675
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