sources.hh

Go to the documentation of this file.
1 
10 #ifndef waveprop_sources_hh
11 #define waveprop_sources_hh
12 
13 #include <iostream>
14 #include <memory>
15 #include "stdio.h"
16 #include "matfile.hh"
17 #include "basics.hh"
18 #include "formula/elementFormulaRCP.hh"
19 #include "formula/formula.hh"
20 #include "formula/piecewiseFormula.hh"
21 
22 #include "geometry/formula.hh"
23 #include "operator/compositions.hh"
24 #include "operator/sparseMatrix.hh"
25 #ifdef HAS_MUMPS
26 #include "operator/mumps.hh"
27 #else
28 #ifdef HAS_SuperLU
29 #include "operator/superLU.hh"
30 #else
31 #include "operator/gmres.hh"
32 #endif
33 #endif
34 #include "hp2D.hh"
35 
36 //#include "app-juliette/PWScatLay.hh"
37 
38 #ifdef HAS_MATLABINTERFACE
39 #include "matfile/matlabMatfile.hh"
40 #endif
41 
42 #define FRMLAYPW_D 0
43 #define FRMLAYTOT_D 0
44 
45 namespace concepts {
46 
47  // ******************************************************** FormulaExpImag1D **
48 
51  class FormulaExpImag1D : public Formula<Cmplx> {
52  public:
53  FormulaExpImag1D(const Real k, const Cmplx u = 1.0, Real x0 = 0.0)
54  : k_(k), x0_(x0), u_(u) {}
55 
56  virtual FormulaExpImag1D* clone() const {
57  return new FormulaExpImag1D(k_, u_, x0_);
58  }
59 
60  virtual Cmplx operator() (const Real p, const Real t = 0.0) const {
61  const Real arg = k_*(p-x0_);
62  return Cmplx(cos(arg), sin(arg)) * u_;
63  }
64  virtual Cmplx operator() (const Real2d& p, const Real t = 0.0) const {
65  return (*this)(p[0], t);
66  }
67  virtual Cmplx operator() (const Real3d& p, const Real t = 0.0) const {
68  return (*this)(p[0], t);
69  }
70  protected:
71  virtual std::ostream& info(std::ostream& os) const {
72  return os << concepts::typeOf(*this)<<"(" << u_ << " * exp(i " << k_ << "(x - " << x0_ << ")))";
73  }
74  private:
75  const Real k_, x0_;
76  const Cmplx u_;
77  };
78 
79  // ******************************************************** FormulaExpImag2D **
80 
83  class FormulaExpImag2D : public Formula<Cmplx> {
84  public:
85  FormulaExpImag2D(const Real2d k, const Cmplx u = 1.0, Real2d x0 = 0.0)
86  : k(k), x0(x0), u(u) {}
87 
88  virtual FormulaExpImag2D* clone() const {
89  return new FormulaExpImag2D(k, u, x0);
90  }
91 
92  virtual Cmplx operator() (const Real p, const Real t = 0.0) const {
93  conceptsThrowSimpleE("FormulaExpImag2D currently only supports 2D!"); assert(false);
94  }
95 
96  virtual Cmplx operator() (const Real2d& p, const Real t = 0.0) const
97  {
98  const Real arg = k[0]*(p[0]-x0[0])
99  + k[1]*(p[1]-x0[1]) ;
100  //cout << "e^{i k x}: " << p << ", " << Cmplx(cos(arg), sin(arg)) * u << endl;
101  return Cmplx(cos(arg), sin(arg)) * u;
102  }
103 
104  virtual Cmplx operator() (const Real3d& p, const Real t = 0.0) const {
105  return (*this)( Real2d(p[0], p[1]), t);
106  }
107  protected:
108  virtual std::ostream& info(std::ostream& os) const {
109  return os << concepts::typeOf(*this)<<"(" << u << " * exp(i " << k << "(x - " << x0 << ")))";
110  }
111  public:
112  const Real2d k, x0;
113  const Cmplx u;
114  };
115 
116  // ********************************************* FormulaExpImag2DRadialDer **
117 
120  class FormulaExpImag2DRadialDer : public Formula<Cmplx> {
121  public:
122  //FormulaExpImag2DRadialDer(const Real2d k, const Cmplx u = 1.0, Real2d x0 = Real2d(0,0))
123  FormulaExpImag2DRadialDer(const Real2d k, const Cmplx u = 1.0)
124  : k(k), x0(0,0), u(u) {}
125 
126  virtual FormulaExpImag2DRadialDer* clone() const {
127  return new FormulaExpImag2DRadialDer(*this);
128  }
129 
130  virtual Cmplx operator() (const Real p, const Real t = 0.0) const {
131  conceptsThrowSimpleE("FormulaExpImag2DRadialDer currently only supports 2D!"); assert(false);
132  }
133 
134  virtual Cmplx operator() (const Real2d& p, const Real t = 0.0) const
135  {
136  const Real arg = k[0]*(p[0]-x0[0])
137  + k[1]*(p[1]-x0[1]) ;
138  Cmplx val(cos(arg), sin(arg));
139  val *= u;
140 
141  const Real len = sqrt( p[0]*p[0] + p[1]*p[1] );
142 
143  return Cmplx(0, (k[0]*p[0] + k[1]*p[1])/len) * val;
144  }
145 
146  virtual Cmplx operator() (const Real3d& p, const Real t = 0.0) const {
147  return (*this)( Real2d(p[0], p[1]), t);
148  }
149  protected:
150  virtual std::ostream& info(std::ostream& os) const {
151  return os << concepts::typeOf(*this)<<"(" << u << " * exp(i " << k << "(x - " << x0 << ")))";
152  }
153  public:
154  const Real2d k, x0;
155  const Cmplx u;
156  };
157 
158  // ************************************************** FormulaExpImag2DGrad **
159 
163  class FormulaExpImag2DGrad : public Formula<Point<Cmplx, 2> >
164  {
165  public:
166  FormulaExpImag2DGrad(const Real2d k, const Cmplx u = 1.0, Real2d x0 = Real2d(0,0))
167  : k(k), x0(x0), u(u)
168  { }
169 
170  virtual FormulaExpImag2DGrad* clone() const {
171  return new FormulaExpImag2DGrad(k, u, x0);
172  }
173 
174  virtual Cmplx2d operator() (const Real p, const Real t = 0.0) const {
175  conceptsThrowSimpleE("FormulaExpImag2DGrad currently only supports 2D!"); assert(false);
176  }
177 
178  virtual Cmplx2d operator() (const Real2d& p, const Real t = 0.0) const
179  {
180  const Cmplx imag(0,1);
181 
182  const Real arg = k[0]*(p[0]-x0[0])
183  + k[1]*(p[1]-x0[1]) ;
184  Cmplx val(cos(arg), sin(arg));
185  val *= u;
186 
187  //const Real len = sqrt( p[0]*p[0] + p[1]*p[1] );
188 
189  return Cmplx2d( Cmplx(0, k[0])*val, Cmplx(0, k[1])*val );
190  }
191 
192  virtual Cmplx2d operator() (const Real3d& p, const Real t = 0.0) const {
193  //conceptsThrowSimpleE("FormulaExpImag2DGrad currently only supports 2D!"); assert(false);
194  return (*this)( Real2d(p[0], p[1]), t);
195  }
196  protected:
197  virtual std::ostream& info(std::ostream& os) const {
198  return os << concepts::typeOf(*this)<<"(" << u << " * exp(i " << k << "(x - " << x0 << ")))";
199  }
200  public:
201  const Real2d k, x0;
202  const Cmplx u;
203  };
204 
205 
206  // ************************************************** FormulaNormalOuterSP2D **
207 
222  template< class F >
223  class FormulaNormalOuterSP2D : public ElementFormula< F >
224  {
225  public:
226  static const bool OUTER = true;
227  static const bool INNER = false;
228 
229  FormulaNormalOuterSP2D(
230  const ElementFormulaContainer< Point<F, 2> > vf,
231  Real2d center = Real2d(0,0), bool direction = OUTER)
232  : vf(vf)
233  , center(center)
234  , direction(direction)
235  { }
236 
237  virtual FormulaNormalOuterSP2D* clone() const {
238  return new FormulaNormalOuterSP2D(*this);
239  }
240 
241  virtual F operator() (const ElementWithCell<Real>& elm,
242  const Real p, const Real t = 0.0) const
243  {
244  const Cell& c = elm.cell();
245 
246  const Edge2d* cc = dynamic_cast<const Edge2d*> (&c);
247  assert(cc);
248 
249  Real2d normal = cc->jacobian(p);
250  normal.ortho();
251  normal /= normal.l2();
252 
253  Real2d pp = elm.elemMap(p); //pp /= pp.l2();
254  pp -= center;
255 
256  if ( ( normal.dot(pp) < 0 && direction == OUTER ) ||
257  ( normal.dot(pp) > 0 && direction == INNER )
258  )
259  {
260  normal *= -1;
261  }
262 
263  const Point<F, 2> u = vf.operator()(elm, p, t);
264 
265  //std::cout << "normal " << elm.elemMap(p) << ": " << normal
266  // << "SP: " << normal[0] * u[0] + normal[1] * u[1] << "center: " << center << std::endl;
267 
268  return normal[0] * u[0] + normal[1] * u[1];
269  }
270 
271  virtual F operator() (const ElementWithCell<Real>& elm,
272  const Real2d& p, const Real t = 0.0) const
273  {
274  conceptsThrowSimpleE("FormulaNormalOuterSP2D currently only supports 1D!"); assert(false);
275  }
276 
277  virtual F operator() (const ElementWithCell<Real>& elm,
278  const Real3d& p, const Real t = 0.0) const {
279  conceptsThrowSimpleE("FormulaNormalOuterSP2D currently only supports 1D!"); assert(false);
280  //return (*this)( Real2d(p[0], p[1]), t);
281  }
282  protected:
283  virtual std::ostream& info(std::ostream& os) const {
284  return os << concepts::typeOf(*this)<<"(" << vf << ", " << center << ", " << direction << ")";
285  }
286  public:
287  const ElementFormulaContainer< Point<F, 2> > vf;
288  Real2d center;
289  bool direction;
290  };
291 
292  // ************************************************** ComposeFormulaMatVec **
293 
298  template< class F, uint DIM, typename G = typename Realtype<F>::type>
299  class ComposeFormulaMatVec : public ElementFormula< Point<F, DIM>, G >
300  {
301  public:
302  ComposeFormulaMatVec(
303  RCP<const ElementFormula< Mapping<F, DIM>, G> > A,
304  RCP<const ElementFormula< Point<F, DIM> , G> > vf
305  )
306  : A(A)
307  , vf(vf)
308  { }
309 
310  virtual ComposeFormulaMatVec* clone() const {
311  return new ComposeFormulaMatVec(*this);
312  }
313 
314  virtual Point<F, DIM> operator() (const ElementWithCell<G>& elm,
315  const Real p, const Real t = 0.0) const
316  {
317  return A->operator()(elm, p, t) * vf->operator()(elm, p, t);
318  }
319 
320  virtual Point<F, DIM> operator() (const ElementWithCell<G>& elm,
321  const Real2d& p, const Real t = 0.0) const
322  {
323  //cout << "p: " << p << A->operator()(elm, p, t) << vf->operator()(elm, p, t)
324  // << " = " << A->operator()(elm, p, t) * vf->operator()(elm, p, t) << endl;
325  return A->operator()(elm, p, t) * vf->operator()(elm, p, t);
326  }
327 
328  virtual Point<F, DIM> operator() (const ElementWithCell<G>& elm,
329  const Real3d& p, const Real t = 0.0) const
330  {
331  return A->operator()(elm, p, t) * vf->operator()(elm, p, t);
332  }
333  protected:
334  virtual std::ostream& info(std::ostream& os) const {
335  return os << concepts::typeOf(*this)<<"(A=" << *A << ", vf="<< *vf << ")";
336  }
337  public:
338  RCP<const ElementFormula< Mapping<F, DIM>, G> > A;
339  RCP<const ElementFormula< Point<F, DIM>, G> > vf;
340  };
341 
342  // *************************************************** ComposeFormulaVecEntry **
343  // * where it's needed?
344  template< class F, uint DIM, typename G = typename Realtype<F>::type>
345  class ComposeFormulaVecEntry : public ElementFormula< F, G >
346  {
347  public:
348  ComposeFormulaVecEntry(
349  RCP<const ElementFormula< Point<F, DIM> , G> > vf,
350  int index
351  )
352  : vf(vf)
353  , index(index)
354  { }
355 
356  virtual ComposeFormulaVecEntry* clone() const {
357  return new ComposeFormulaVecEntry(*this);
358  }
359 
360  virtual F operator() (const ElementWithCell<G>& elm,
361  const Real p, const Real t = 0.0) const
362  {
363  return vf->operator()(elm, p, t)[index];
364  }
365 
366  virtual F operator() (const ElementWithCell<G>& elm,
367  const Real2d& p, const Real t = 0.0) const
368  {
369  return vf->operator()(elm, p, t)[index];
370  }
371 
372  virtual F operator() (const ElementWithCell<G>& elm,
373  const Real3d& p, const Real t = 0.0) const
374  {
375  return vf->operator()(elm, p, t)[index];
376  }
377  protected:
378  virtual std::ostream& info(std::ostream& os) const {
379  return os << concepts::typeOf(*this)<<"(i=" << index << ", vf="<< *vf << ")";
380  }
381  public:
382  RCP<const ElementFormula< Point<F, DIM>, G> > vf;
383  int index;
384  };
385 
386  // *********************************************** FormulaIncPlaneWaveSource **
390  class FormulaIncPlaneWaveSource : public ElementFormula<Cmplx> {
391  public:
392  typedef std::map<int, double> MID;
393 
394  FormulaIncPlaneWaveSource(ElementFormulaContainer<Cmplx> coeff_a
395  , ElementFormulaContainer<Cmplx> coeff_b
396  , RCP< FormulaExpImag2D > u_inc )
397  : coeff_a(coeff_a)
398  , coeff_b(coeff_b)
399  , u_inc(u_inc)
400  { }
401 
402  virtual FormulaIncPlaneWaveSource* clone() const {
403  return new FormulaIncPlaneWaveSource(*this);
404  }
405 
406  virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real3d &p, const Real t=0.0) const
407  {
408  Real2d px = elm.elemMap(p);
409  return operator()(elm, p, px, t);
410  }
411 
412  virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real2d &p, const Real t=0.0) const
413  {
414  Real2d px = elm.elemMap(p);
415  return operator()(elm, p, px, t);
416  }
417 
418  virtual Cmplx operator() (const ElementWithCell< Real > &elm, const Real p, const Real t=0.0) const
419  {
420  Real2d px = elm.elemMap(p);
421  return operator()(elm, p, px, t);
422  }
423 
424  template <class RealNd>
425  Cmplx operator() (const ElementWithCell< Real > &elm, const RealNd &p, Real2d px,
426  const Real t=0.0) const
427  {
428  Real2d k = u_inc->k;
429  Cmplx a = coeff_a(elm, p);
430  Cmplx b = coeff_b(elm, p);
431  return (b - a) * (k[0]*k[0] + k[1]*k[1]) * u_inc->operator()(elm.elemMap(p), t);
432  }
433 
434  static PiecewiseConstFormula<Cmplx> genTMCoeff(
435  const MID& attToEps);
436 
437  static PiecewiseConstFormula<Cmplx> genTECoeff(
438  const MID& attToEps);
439 
440  protected:
441  virtual std::ostream& info(std::ostream& os) const {
442  return os << concepts::typeOf(*this)<<"(" << coeff_a << ", " << coeff_b << ")";
443  }
444  private:
445  ElementFormulaContainer<Cmplx> coeff_a, coeff_b;
446  RCP< FormulaExpImag2D > u_inc;
447  };
448 
449  PiecewiseConstFormula<Cmplx> FormulaIncPlaneWaveSource::genTMCoeff(
450  const MID& attrToEps)
451  {
452  PiecewiseConstFormula<Cmplx> coeffs(1.);
453  for(MID::const_iterator it = attrToEps.begin();
454  it != attrToEps.end(); ++it)
455  {
456  coeffs[it->first] = it->second;
457 
458  if(it->second < 0) {
459  double abs_fact = 1e-2;
460  printf("Adding absorbing material (%f %f) for attr %d\n",
461  -it->second, abs_fact, it->first);
462  coeffs[it->first] = -it->second;
463 #ifdef __cplusplus >= 201103L
464  coeffs[it->first].imag(abs_fact);
465 #else
466  imag(coeffs[it->first]) = abs_fact;
467 #endif
468 
469  } else {
470  printf("Adding non-absorbing material (%f %f) for attr %d\n",
471  it->second, 0., it->first);
472  }
473  }
474  return coeffs;
475  }
476 
477  PiecewiseConstFormula<Cmplx> FormulaIncPlaneWaveSource::genTECoeff(
478  const MID& attrToEps)
479  {
480  PiecewiseConstFormula<Cmplx> coeffs(1.);
481  for(MID::const_iterator it = attrToEps.begin();
482  it != attrToEps.end(); ++it)
483  {
484  coeffs[it->first] = 1. / it->second;
485 
486  if(it->second < 0) {
487  coeffs[it->first] = -it->second;
488 #ifdef __cplusplus >= 201103L
489  coeffs[it->first].imag(1e-3);
490 #else
491  imag(coeffs[it->first]) = 1e-3;
492 #endif
493  }
494  }
495  return coeffs;
496  }
497 
498 
499 
500 
501  // ********************************************* FormulaLayerPlaneWaveSource **
502 
505  class FormulaLayerPlaneWaveSource : public Formula<Cmplx> {
506  public:
507 
508 #ifdef HAS_MATLABINTERFACE
509 
510  FormulaLayerPlaneWaveSource(const std::string& filename)
511  :filename_(filename)
512  {
513  MatlabMatfile matfile(filename_);
514  d_ .reset(matfile.getVector<Real> ("d") );
515  ky_.reset(matfile.getVector<Real> ("ky"));
516  A_ .reset(matfile.getVector<Cmplx>("A") );
517  B_ .reset(matfile.getVector<Cmplx>("B") );
518  std::unique_ptr<Vector<Real> >
519  kx(matfile.getVector<Real>("kx"));
520  kx_ = **kx;
521  N_ = d_->size();
522  }
523 #endif
524 
525  FormulaLayerPlaneWaveSource(Vector<Real>& epsilon, Vector<Real>& d,
526  const Real& kx, const Real omega)
527  {
528  N_ = epsilon.size();
529  d_.reset( new Vector<Real>(N_-1) ) ;
530  epsilon_.reset( new Vector<Real>(N_) ) ;
531  kx_ = kx;
532  ky_.reset( new Vector<Real>(N_) ) ;
533  A_.reset( new Vector<Cmplx>(N_) ) ;
534  B_.reset( new Vector<Cmplx>(N_) ) ;
535  rho_.reset( new Vector<Cmplx>(N_) ) ;
536  omega_ = omega;
537 
538  for (uint i = 0; i < N_; i++){
539  (*epsilon_)(i)=epsilon(i);
540  }
541  for (uint i = 0; i < N_-1; i++){
542  (*d_)(i)=d(i);
543  }
544  //cout << "epsilon = " << *epsilon_ << endl;
545  Construct();
546  }
547 
548 
549  void Construct(){
550  Cmplx IMAG(0.0 , 1.0);
551  SparseMatrix<Cmplx> M(2*N_);
552  Vector<Real> a(N_);
553  for (uint i = 0 ; i < N_; i++)
554  {
555  // bi = 1.0
556  // ai = 1.0/epsilon(i)
557  // ky(i) = sqrt(omega^2 * bi/ai - kx^2) = -sqrt(epsilon(i) - kx^2)
558  Real w = omega_*omega_*(*epsilon_)(i) - kx_*kx_;
559 
560  conceptsAssert3(w >= 0, Assertion(),
561  "problem in the " << i << " th layer, ky is complex");
562 
563  (*ky_)(i) = -sqrt(w);
564  a(i) = 1.0/(*epsilon_)(i);
565  }
566  for (uint i=0; i < N_-1; i++)
567  {
568  M(i,i) = exp( IMAG * (*ky_)(i) *(*d_)(i) );
569  M(i,i+1) = - exp( IMAG * (*ky_)(i+1) *(*d_)(i) );
570  M(i,N_+i) = exp(-1.0 *IMAG * (*ky_)(i) *(*d_)(i) );
571  M(i,N_+i+1) = - exp(-1.0 *IMAG * (*ky_)(i+1) *(*d_)(i) );
572 
573  M(N_+i,i) = a(i) * (*ky_)(i) * exp( IMAG* (*ky_)(i) *(*d_)(i) );
574  M(N_+i,i+1) = -1.0 * a(i+1) * (*ky_)(i+1) * exp( IMAG* (*ky_)(i+1) *(*d_)(i) );
575  M(N_+i,N_+i) = -1.0 * a(i) * (*ky_)(i) * exp(-1.0*IMAG* (*ky_)(i) *(*d_)(i) );
576  M(N_+i,N_+i+1) = a(i+1) * (*ky_)(i+1) * exp(-1.0*IMAG* (*ky_)(i+1) *(*d_)(i) );
577  }
578  for (uint i=0; i<2* N_; i++)
579  {
580  M(N_-1, i ) = 0.0 ;
581  M(2*N_-1, i) = 0.0 ;
582  }
583  M(N_-1,0) = 1.0 ;
584  M(2*N_-1,2*N_-1) = 1.0 ;
585 
586  // cout << "M = " << M << endl;
587  // MatfileOutput mymat("matlabM");
588  // mymat.add(M,"M");
589 
590  // the matrix is ready to be inverted
591  // we prepare the rhs
592  Vector<Cmplx> X (2*N_);
593  for (uint i=0; i<2*N_;i++)
594  {
595  X(i) = 0.0 ;
596  }
597  X(N_-1) = 1.0 ;
598  // the right hand side is ready
599  // We solve Y = M^{-1} X
600  std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
601 #ifdef HAS_MUMPS
602  solver.reset(new concepts::Mumps<Cmplx>(M));
603 #else
604 #ifdef HAS_SuperLU
605  solver.reset(new concepts::SuperLU<Cmplx>(M));
606 #else
607  solver.reset(new concepts::GMRes<Cmplx>(M, 1e-12));
608 #endif
609 #endif
610 
611  // ** solve **
612  concepts::Vector<Cmplx> Y(2*N_);
613  (*solver)(X, Y);
614 
615  // We copy the result into A and B
616  for (uint i=0; i<N_; i++)
617  {
618  (*A_)(i) = Y(i);
619  (*B_)(i) = Y(N_+i);
620  }
621  }
622 
623  // Formula from Kong's book : does it work ?
624  void ConstructFromBook()
625  {
626 
627  Cmplx IMAG(0.0 , 1.0);
628  // We know ky(i) from Descartes Law
629  Real normk2 = (kx_*kx_+(*ky_)(0)*(*ky_)(0))/(*epsilon_)(0);
630  for (uint i = 0; i < N_; i++){
631  (*ky_)(i) = sqrt( normk2*(*epsilon_)(i)-kx_*kx_);
632  }
633  // We know the amplitudes from the continuity of u and \partial_n u
634  (*B_)(0) = 1.0;
635  for ( int i = N_-2; i>=0; i--)
636  {
637  Cmplx ei = (*epsilon_)(i);
638  Cmplx eip=(*epsilon_)(i+1);
639  Cmplx kyi = (*ky_)(i);
640  Cmplx kyip=(*ky_)(i+1);
641  Cmplx p= ei*kyip /( eip*kyi);
642  Cmplx R = (1.0-p)/(1.0+p);
643  Cmplx R_sq = R*R;
644  Cmplx atmp = 2.0*IMAG*(*ky_)(i)*(*d_)(i);
645  Cmplx btmp=2.0*IMAG* ((*ky_)(i)+(*ky_)(i+1))*(*d_)(i);
646  Cmplx ctmp = 2.0*IMAG*(*ky_)(i+1)*(*d_)(i);
647  Cmplx dtmp=exp(atmp)/R ;
648  Cmplx etmp =(1.0-1.0/R_sq) * exp( btmp)/(exp(ctmp)/R+(*rho_)(i+1));
649  (*rho_)(i) = dtmp+etmp;
650  }
651  for ( uint i = 1; i < N_; i++)
652  {
653  Cmplx atmp = -IMAG*(*ky_)(i-1)*(*d_)(i-1);
654  Cmplx btmp = IMAG*(*ky_)(i-1)*(*d_)(i-1);
655  Cmplx ctmp = -IMAG*(*ky_)(i)*(*d_)(i-1);
656  Cmplx dtmp = IMAG*(*ky_)(i)*(*d_)(i-1);
657  Cmplx u =(*rho_)(i-1)*exp(atmp)+exp(btmp);
658  Cmplx v = (*rho_)(i)*exp(ctmp)+exp(dtmp);
659  (*B_)(i) = (*B_)(i-1)*u/v;
660  (*A_)(i) = (*B_)(i)* (*rho_)(i);
661  }
662  (*A_)(0) = (*B_)(0)*(*rho_)(0);
663  }
664 
665  void Display(){
666  std::cout << " epsilon = " << *epsilon_ << std::endl;
667  std::cout << " d = " << *d_ << std::endl;
668  std::cout << " kx = " << kx_ << std::endl;
669  std::cout << " omega = " << omega_ << std::endl;
670  std::cout << " ky = " << *ky_ << std::endl;
671  std::cout << " A = " << *A_ << std::endl;
672  std::cout << " B = " << *B_ << std::endl;
673  }
674 
675  virtual FormulaLayerPlaneWaveSource* clone() const {
676 #ifdef HAS_MATLABINTERFACE
677  return new FormulaLayerPlaneWaveSource(filename_);
678 #endif
679  return new FormulaLayerPlaneWaveSource(*epsilon_, *d_, kx_, omega_);
680  }
681 
682 
683  virtual Cmplx operator() (const Real3d &p, const Real t=0.0) const
684  {
685  return (*this)(Real2d(p[0], p[1]),t);
686  }
687 
688  virtual Cmplx operator() (const Real2d &p, const Real t=0.0) const
689  {
690  Cmplx coeff_a, coeff_b;
691  Real KY;
692  // Determine the layer
693  uint index = 0, i;
694  for(i = 0; i < N_-1; index = ++i)
695  if (p[1] >= (*d_)[i])
696  break;
697 
698  DEBUGL(FRMLAYPW_D, "The point " << p << " is in the " << index << "th layer.");
699  // Assign the value and calculate
700  coeff_a = (*A_)[index];
701  coeff_b = (*B_)[index];
702  KY = ((*ky_)[index]);
703  DEBUGL(FRMLAYPW_D, "In this layer, A = " << coeff_a << ", B = "
704  << coeff_b << ", kx = " << kx_ << ", ky = " << KY);
705  return coeff_a * exp(Cmplx(0, kx_*p[0]+KY*p[1])) +
706  coeff_b * exp(Cmplx(0, kx_*p[0]-KY*p[1]));
707  }
708 
709  virtual Cmplx operator() (const Real p, const Real t=0.0) const
710  {
711  conceptsThrowSimpleE("FormulaLayerPlaneWaveSource currently only supports 2D!");
712  assert(false);
713  }
714 
715  protected:
716  virtual std::ostream& info(std::ostream& os) const {
717  return os << concepts::typeOf(*this)<<"(" << "layers: " << d_ << ", ky: " << ky_
718  << ", A_:" << A_ << ", B_:" << B_ <<")";
719  }
720  private:
721  const std::string filename_;
722  std::unique_ptr<Vector<Real> > d_;
723  std::unique_ptr<Vector<Real> > epsilon_;
724  std::unique_ptr<Vector<Cmplx> > rho_;
725  std::unique_ptr<Vector<Real> > ky_;
726  std::unique_ptr<Vector<Cmplx> > A_;
727  std::unique_ptr<Vector<Cmplx> > B_;
728  Real kx_;
729  Real omega_;
730  uint N_;
731  };
732 
733  // *************************************** FormulaLayerPlaneWaveSourceGrad **
734 
737  class FormulaLayerPlaneWaveSourceGrad : public Formula<Cmplx2d> {
738  public:
739 
740 #ifdef HAS_MATLABINTERFACE
741  FormulaLayerPlaneWaveSourceGrad(const std::string& filename)
742  :filename_(filename)
743  {
744  MatlabMatfile matfile(filename_);
745  d_ .reset(matfile.getVector<Real> ("d") );
746  ky_.reset(matfile.getVector<Real> ("ky"));
747  A_ .reset(matfile.getVector<Cmplx>("A") );
748  B_ .reset(matfile.getVector<Cmplx>("B") );
749  std::unique_ptr<Vector<Real> > kx(matfile.getVector<Real>("kx"));
750  kx_ = **kx;
751  N_ = d_->size();
752  }
753 #endif
754 
755  FormulaLayerPlaneWaveSourceGrad(Vector<Real>& epsilon, Vector<Real>& d,
756  const Real& kx, const Real omega)
757  {
758  N_ = epsilon.size();
759  d_.reset( new Vector<Real>(N_-1) ) ;
760  epsilon_.reset( new Vector<Real>(N_) ) ;
761  kx_ = kx;
762  omega_ = omega;
763  ky_.reset( new Vector<Real>(N_) ) ;
764  A_.reset( new Vector<Cmplx>(N_) ) ;
765  B_.reset( new Vector<Cmplx>(N_) ) ;
766  rho_.reset( new Vector<Cmplx>(N_) ) ;
767 
768  for (uint i = 0; i < N_; i++){
769  (*epsilon_)(i)=epsilon(i);
770  }
771  for (uint i = 0; i < N_-1; i++){
772  (*d_)(i)=d(i);
773  }
774  Construct();
775  }
776 
777 
778 
779  void Construct(){
780  Cmplx IMAG(0.0 , 1.0);
781  SparseMatrix<Cmplx> M(2*N_);
782  Vector<Real> a(N_);
783  for (uint i = 0 ; i < N_; i++)
784  { // bi = 1.0
785  // ai = 1.0/epsilon(i)
786  // ky(i) = sqrt(omega^2 * bi/ai - kx^2) = -sqrt(epsilon(i) - kx^2)
787  Real w = omega_*omega_*(*epsilon_)(i) - kx_*kx_;
788 
789  conceptsAssert3(w >= 0, Assertion(),
790  "problem in the " << i << " th layer, ky is complex");
791 
792  (*ky_)(i) = -sqrt(w);
793  a(i) = 1.0/(*epsilon_)(i);
794  }
795 
796  // cout << "M = " << M << endl;
797  // MatfileOutput mymat("matlabM");
798  // mymat.add(M,"M");
799 
800  // the matrix is ready to be inverted
801  // we prepare the rhs
802  Vector<Cmplx> X (2*N_);
803  for (uint i=0; i<2*N_;i++)
804  {
805  X(i) = 0.0 ;
806  }
807  X(N_-1) = 1.0 ;
808  // the right hand side is ready
809  // We solve Y = M^{-1} X
810  std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
811 #ifdef HAS_MUMPS
812  solver.reset(new concepts::Mumps<Cmplx>(M));
813 #else
814 #ifdef HAS_SuperLU
815  solver.reset(new concepts::SuperLU<Cmplx>(M));
816 #else
817  solver.reset(new concepts::GMRes<Cmplx>(M, 1e-12));
818 #endif
819 #endif
820 
821  // ** solve **
822  concepts::Vector<Cmplx> Y(2*N_);
823  (*solver)(X, Y);
824 
825  // We copy the result into A and B
826  for (uint i=0; i<N_; i++)
827  {
828  (*A_)(i) = Y(i);
829  (*B_)(i) = Y(N_+i);
830  }
831  }
832 
833 
834  void ConstructFromBook(){
835 
836  Cmplx IMAG(0.0 , 1.0);
837  // We know ky(i) from Descartes Law
838  Real normk2 = (kx_*kx_+(*ky_)(0)*(*ky_)(0))/(*epsilon_)(0);
839  for (uint i = 0; i < N_; i++){
840  (*ky_)(i) = sqrt( normk2*(*epsilon_)(i)-kx_*kx_);
841  }
842  // We know the amplitudes from the continuity of u and \partial_n u
843  (*B_)(0) = 1.0;
844  for ( int i = N_-2; i>=0; i--)
845  {
846  Cmplx ei = (*epsilon_)(i);
847  Cmplx eip=(*epsilon_)(i+1);
848  Cmplx kyi = (*ky_)(i);
849  Cmplx kyip=(*ky_)(i+1);
850  Cmplx p= ei*kyip /( eip*kyi);
851  Cmplx R = (1.0-p)/(1.0+p);
852  Cmplx R_sq = R*R;
853  Cmplx atmp = 2.0*IMAG*(*ky_)(i)*(*d_)(i);
854  Cmplx btmp=2.0*IMAG* ((*ky_)(i)+(*ky_)(i+1))*(*d_)(i);
855  Cmplx ctmp = 2.0*IMAG*(*ky_)(i+1)*(*d_)(i);
856  Cmplx dtmp=exp(atmp)/R ;
857  Cmplx etmp =(1.0-1.0/R_sq) * exp( btmp)/(exp(ctmp)/R+(*rho_)(i+1));
858  (*rho_)(i) = dtmp+etmp;
859  }
860  for ( uint i = 1; i < N_; i++)
861  {
862  Cmplx atmp = -IMAG*(*ky_)(i-1)*(*d_)(i-1);
863  Cmplx btmp = IMAG*(*ky_)(i-1)*(*d_)(i-1);
864  Cmplx ctmp = -IMAG*(*ky_)(i)*(*d_)(i-1);
865  Cmplx dtmp = IMAG*(*ky_)(i)*(*d_)(i-1);
866  Cmplx u =(*rho_)(i-1)*exp(atmp)+exp(btmp);
867  Cmplx v = (*rho_)(i)*exp(ctmp)+exp(dtmp);
868  (*B_)(i) = (*B_)(i-1)*u/v;
869  (*A_)(i) = (*B_)(i)* (*rho_)(i);
870  }
871  (*A_)(0) = (*B_)(0)*(*rho_)(0);
872 
873  }
874 
875  void Display(){
876  std::cout << " epsilon = " << *epsilon_ << std::endl;
877  std::cout << " d = " << *d_ << std::endl;
878  std::cout << " kx = " << kx_ << std::endl;
879  std::cout << " omega = " << omega_ << std::endl;
880  std::cout << " ky = " << *ky_ << std::endl;
881  std::cout << " A = " << *A_ << std::endl;
882  std::cout << " B = " << *B_ << std::endl;
883  }
884 
885 
886 
887 
888 
889  virtual FormulaLayerPlaneWaveSourceGrad* clone() const {
890 #ifdef HAS_MATLABINTERFACE
891  return new FormulaLayerPlaneWaveSourceGrad(filename_);
892 #endif
893  return new FormulaLayerPlaneWaveSourceGrad(*epsilon_, *d_, kx_, omega_);
894  }
895 
896 
897  virtual Cmplx2d operator() (const Real3d &p, const Real t=0.0) const
898  {
899  return (*this)(Real2d(p[0], p[1]),t);
900  }
901 
902  virtual Cmplx2d operator() (const Real2d &p, const Real t=0.0) const
903  {
904  Cmplx coeff_a, coeff_b;
905  Real KY;
906  // Determine the layer
907  uint index = 0, i;
908  for(i = 0; i < N_-1 ; index = ++i)
909  if (p[1] >= (*d_)[i])
910  break;
911 
912  DEBUGL(FRMLAYPW_D, "The point " << p << " is in the " << index << "th layer.");
913  // Assign the value and calculate
914  coeff_a = (*A_)[index];
915  coeff_b = (*B_)[index];
916  KY = (*ky_)[index];
917  DEBUGL(FRMLAYPW_D, "In this layer, A = " << coeff_a << ", B = "
918  << coeff_b << ", kx = " << kx_ << ", ky = " << KY);
919  return Cmplx2d(Cmplx(0.,1.)*kx_*coeff_a*exp(Cmplx(0.,kx_*p[0]+KY*p[1])),
920  Cmplx(0.,1.)*KY *coeff_a*exp(Cmplx(0.,kx_*p[0]+KY*p[1])))
921  +Cmplx2d(Cmplx(0., 1.)*kx_*coeff_b*exp(Cmplx(0.,kx_*p[0]-KY*p[1])),
922  Cmplx(0.,-1.)*KY *coeff_b*exp(Cmplx(0.,kx_*p[0]-KY*p[1])));
923  }
924 
925  virtual Cmplx2d operator() (const Real p, const Real t=0.0) const
926  {
927  conceptsThrowSimpleE("FormulaLayerPlaneWaveSourceGrad currently only supports 2D!");
928  assert(false);
929  }
930 
931  protected:
932  virtual std::ostream& info(std::ostream& os) const {
933  return os << concepts::typeOf(*this)<<"(" << "layers: " << d_ << ", ky: " << ky_
934  << ", A_:" << A_ << ", B_:" << B_ <<")";
935  }
936  private:
937  const std::string filename_;
938  std::unique_ptr<Vector<Real> > d_;
939  std::unique_ptr<Vector<Real> > epsilon_;
940  std::unique_ptr<Vector<Cmplx> > rho_;
941  std::unique_ptr<Vector<Real> > ky_;
942  std::unique_ptr<Vector<Cmplx> > A_;
943  std::unique_ptr<Vector<Cmplx> > B_;
944  Real kx_;
945  Real omega_;
946  uint N_;
947  };
948 
949  // ********************************************** FormulaLayerPlaneWaveLayer **
950 
958  class FormulaLayerPlaneWaveLayer : public Formula<Cmplx> {
959 
960  protected:
962  Cmplx a_;
964  Cmplx b_;
966  Real kx_;
968  Real ky_;
970  int layer_label_;
971  public:
972 
973  FormulaLayerPlaneWaveLayer(const Cmplx a, const Cmplx b, const Real kx, const Real ky, const int layer_label=0) :
974  a_(a), b_(b), kx_(kx), ky_(ky), layer_label_(layer_label) {}
975 
976  FormulaLayerPlaneWaveLayer(const FormulaLayerPlaneWaveLayer & other) :
977  a_(other.a_), b_(other.b_), kx_(other.kx_), ky_(other.ky_), layer_label_(other.layer_label_) {}
978 
979  ~FormulaLayerPlaneWaveLayer() {};
980 
981  FormulaLayerPlaneWaveLayer* clone() const {
982  return new FormulaLayerPlaneWaveLayer(*this);
983  }
984 
985  Cmplx operator() (const Real3d &p, const Real t=0.0) const
986  {
987  DEBUGL(FRMLAYTOT_D,p << " for label " << layer_label_);
988  return (*this)(Real2d(p[0], p[1]),t);
989  }
990 
991  Cmplx operator() (const Real2d &p, const Real t=0.0) const
992  {
993  return a_ * exp(Cmplx(0, kx_*p[0]+ky_*p[1])) +
994  b_ * exp(Cmplx(0, kx_*p[0]-ky_*p[1]));
995  }
996 
997  Cmplx operator() (const Real p, const Real t=0.0) const
998  {
999  conceptsThrowSimpleE("FormulaLayerPlaneWaveLayer only supports 2D!");
1000  assert(false);
1001  }
1002 
1003  std::ostream& info(std::ostream& os) const {
1004  return os << concepts::typeOf(*this)<<"("
1005  << "a_: " << a_
1006  << ", b_: " << b_
1007  << ", kx_: " << kx_
1008  << ", ky_: " << ky_ << ")";
1009  }
1010  };
1011 
1012  // ****************************************** FormulaLayerPlaneWaveLayerGrad **
1013 
1021  class FormulaLayerPlaneWaveLayerGrad : public Formula<Cmplx2d> {
1022 
1023  protected:
1025  Cmplx a_;
1027  Cmplx b_;
1029  Real kx_;
1031  Real ky_;
1033  int layer_label_;
1034  public:
1035 
1036  FormulaLayerPlaneWaveLayerGrad(const Cmplx a, const Cmplx b, const Real kx, const Real ky, const int layer_label=0) :
1037  a_(a), b_(b), kx_(kx), ky_(ky), layer_label_(layer_label) {}
1038 
1039  FormulaLayerPlaneWaveLayerGrad(const FormulaLayerPlaneWaveLayerGrad & other) :
1040  a_(other.a_), b_(other.b_), kx_(other.kx_), ky_(other.ky_), layer_label_(other.layer_label_) {}
1041 
1042  ~FormulaLayerPlaneWaveLayerGrad() {};
1043 
1044  FormulaLayerPlaneWaveLayerGrad* clone() const {
1045  return new FormulaLayerPlaneWaveLayerGrad(*this);
1046  }
1047 
1048  Cmplx2d operator() (const Real3d &p, const Real t=0.0) const
1049  {
1050  DEBUGL(FRMLAYTOT_D,p << " for label " << layer_label_);
1051  return (*this)(Real2d(p[0], p[1]),t);
1052  }
1053 
1054  Cmplx2d operator() (const Real2d &p, const Real t=0.0) const
1055  {
1056  return Cmplx2d(Cmplx(0.,1.)*kx_*a_*exp(Cmplx(0.,kx_*p[0]+ky_*p[1])),
1057  Cmplx(0.,1.)*ky_*a_*exp(Cmplx(0.,kx_*p[0]+ky_*p[1])))
1058  +Cmplx2d(Cmplx(0., 1.)*kx_*a_*exp(Cmplx(0.,kx_*p[0]-ky_*p[1])),
1059  Cmplx(0.,-1.)*ky_*b_*exp(Cmplx(0.,kx_*p[0]-ky_*p[1])));
1060  }
1061 
1062  Cmplx2d operator() (const Real p, const Real t=0.0) const
1063  {
1064  conceptsThrowSimpleE("FormulaLayerPlaneWaveLayerGrad only supports 2D!");
1065  assert(false);
1066  }
1067 
1068  std::ostream& info(std::ostream& os) const {
1069  return os << concepts::typeOf(*this)<<"("
1070  << "a_: " << a_
1071  << ", b_: " << b_
1072  << ", kx_: " << kx_
1073  << ", ky_: " << ky_ << ")";
1074  }
1075  };
1076 
1077  // ********************************************** FormulaLayerPlaneWaveTotal **
1078 
1085  class FormulaLayerPlaneWaveTotal : public PiecewiseFormula<Cmplx> {
1086  protected:
1087  public:
1102  FormulaLayerPlaneWaveTotal(const Vector<Real>& epsilon,
1103  const Vector<Real>& d,
1104  const Real kx, const Real omega,
1105  const std::map<int, int> phystolayer)
1106  : PiecewiseFormula<Cmplx>(Cmplx(0.,0.))
1107  {
1108  uint N = epsilon.size();
1109  conceptsAssert(N-1==d.size(), Assertion() );
1110  // Complex imaginary number
1111  Cmplx I(0.0 , 1.0);
1112  // Vectors that will contain the y-direction of our planewave
1113  Vector<Real> ky(N);
1114  // Vector that will contain the coefficients to build the system matrix
1115  Vector<Real> a(N);
1116  // System matrix, RHS and solution
1117  SparseMatrix<Cmplx> M(2*N);
1118  Vector<Cmplx> rhs(2*N);
1119  Vector<Cmplx> sol(2*N);
1120 
1121  // Prepare the coefficients
1122  for (uint i=0 ; i < N; i++)
1123  {
1124  Real w = omega*omega*epsilon(i)-kx*kx;
1125  conceptsAssert3(w >= 0, Assertion(),
1126  "Problem in layer " << i << ", ky is complex");
1127  ky(i) = -sqrt(w);
1128  a(i) = 1.0/epsilon(i);
1129  }
1130 
1131  // Prepare the system matrix
1132  for (uint i=0; i < N-1; i++)
1133  {
1134  M(i,i) = exp( I * ky(i) *d(i) );
1135  M(i,i+1) = - exp( I * ky(i+1) *d(i) );
1136  M(i,N+i) = exp(-1.0 *I * ky(i) *d(i) );
1137  M(i,N+i+1) = - exp(-1.0 *I * ky(i+1) *d(i) );
1138 
1139  M(N+i,i) = a(i) * ky(i) * exp( I* ky(i) *d(i) );
1140  M(N+i,i+1) = -1.0 * a(i+1) * ky(i+1) * exp( I* ky(i+1) *d(i) );
1141  M(N+i,N+i) = -1.0 * a(i) * ky(i) * exp(-1.0*I* ky(i) *d(i) );
1142  M(N+i,N+i+1) = a(i+1) * ky(i+1) * exp(-1.0*I* ky(i+1) *d(i) );
1143  }
1144  for (uint i=0; i<2*N; i++)
1145  {
1146  M(N-1, i ) = 0.0 ;
1147  M(2*N-1, i) = 0.0 ;
1148  }
1149  M(N-1,0) = 1.0 ;
1150  M(2*N-1,2*N-1) = 1.0 ;
1151  DEBUGL(FRMLAYTOT_D, "M=" << M);
1152 
1153  // Solve the system matrix
1154  rhs = Cmplx(0.0);
1155  rhs(N-1) = 1.0;
1156  std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
1157 #ifdef HAS_MUMPS
1158  solver.reset(new concepts::Mumps<Cmplx>(M));
1159 #else
1160 #ifdef HAS_SuperLU
1161  solver.reset(new concepts::SuperLU<Cmplx>(M));
1162 #else
1163  solver.reset(new concepts::GMRes<Cmplx>(M, 1e-12));
1164 #endif
1165 #endif
1166  (*solver)(rhs,sol);
1167  DEBUGL(FRMLAYTOT_D, "sol=" << sol);
1168  // For each layer i, the coeff a(i) is stored in sol(i) and the coeff b(i)
1169  // is stored in b(i)
1170  solver.reset(nullptr);
1171 
1172  // Post-post-processing
1173  std::map<int, int>::const_iterator it = phystolayer.begin();
1174  for( ; it != phystolayer.end() ; it++)
1175  {
1176  uint phys_label = it->first;
1177  uint layer_label = it->second;
1178  DEBUGL(FRMLAYTOT_D, "Adding formula of layer " << layer_label
1179  << " to domain of physical label " << phys_label);
1180  FormulaLayerPlaneWaveLayer formula(sol(layer_label),
1181  sol(N+layer_label),
1182  kx,ky(layer_label),layer_label);
1183  this->set(phys_label, formula);
1184  }
1185  // Set, is used then as a its parent class
1186  }
1187 
1188  ~FormulaLayerPlaneWaveTotal() {};
1189 
1190  };
1191 
1198  class FormulaLayerPlaneWaveTotalGrad : public PiecewiseFormula<Cmplx2d> {
1199  protected:
1200  public:
1215  FormulaLayerPlaneWaveTotalGrad(const Vector<Real>& epsilon,
1216  const Vector<Real>& d,
1217  const Real kx, const Real omega,
1218  const std::map<int, int> phystolayer)
1219  : PiecewiseFormula<Cmplx2d>(Cmplx2d(Cmplx(0.,0.),Cmplx(0.,0.)))
1220  {
1221  uint N = epsilon.size();
1222  conceptsAssert(N-1==d.size(), Assertion() );
1223  // Complex imaginary number
1224  Cmplx I(0.0 , 1.0);
1225  // Vectors that will contain the y-direction of our planewave
1226  Vector<Real> ky(N);
1227  // Vector that will contain the coefficients to build the system matrix
1228  Vector<Real> a(N);
1229  // System matrix, RHS and solution
1230  SparseMatrix<Cmplx> M(2*N);
1231  Vector<Cmplx> rhs(2*N);
1232  Vector<Cmplx> sol(2*N);
1233 
1234  // Prepare the coefficients
1235  for (uint i=0 ; i < N; i++)
1236  {
1237  Real w = omega*omega*epsilon(i)-kx*kx;
1238  conceptsAssert3(w >= 0, Assertion(),
1239  "Problem in layer " << i << ", ky is complex");
1240  ky(i) = -sqrt(w);
1241  a(i) = 1.0/epsilon(i);
1242  }
1243 
1244  // Prepare the system matrix
1245  for (uint i=0; i < N-1; i++)
1246  {
1247  M(i,i) = exp( I * ky(i) *d(i) );
1248  M(i,i+1) = - exp( I * ky(i+1) *d(i) );
1249  M(i,N+i) = exp(-1.0 *I * ky(i) *d(i) );
1250  M(i,N+i+1) = - exp(-1.0 *I * ky(i+1) *d(i) );
1251 
1252  M(N+i,i) = a(i) * ky(i) * exp( I* ky(i) *d(i) );
1253  M(N+i,i+1) = -1.0 * a(i+1) * ky(i+1) * exp( I* ky(i+1) *d(i) );
1254  M(N+i,N+i) = -1.0 * a(i) * ky(i) * exp(-1.0*I* ky(i) *d(i) );
1255  M(N+i,N+i+1) = a(i+1) * ky(i+1) * exp(-1.0*I* ky(i+1) *d(i) );
1256  }
1257  for (uint i=0; i<2*N; i++)
1258  {
1259  M(N-1, i ) = 0.0 ;
1260  M(2*N-1, i) = 0.0 ;
1261  }
1262  M(N-1,0) = 1.0 ;
1263  M(2*N-1,2*N-1) = 1.0 ;
1264  DEBUGL(FRMLAYTOT_D, "M=" << M);
1265 
1266  // Solve the system matrix
1267  rhs = Cmplx(0.0);
1268  rhs(N-1) = 1.0;
1269  std::unique_ptr<concepts::Operator<Cmplx> > solver(nullptr);
1270 #ifdef HAS_MUMPS
1271  solver.reset(new concepts::Mumps<Cmplx>(M));
1272 #else
1273 #ifdef HAS_SuperLU
1274  solver.reset(new concepts::SuperLU<Cmplx>(M));
1275 #else
1276  solver.reset(new concepts::GMRes<Cmplx>(M, 1e-12));
1277 #endif
1278 #endif
1279  (*solver)(rhs,sol);
1280  DEBUGL(FRMLAYTOT_D, "sol=" << sol);
1281  // For each layer i, the coeff a(i) is stored in sol(i) and the coeff b(i)
1282  // is stored in b(i)
1283  solver.reset(0);
1284 
1285  // Post-post-processing
1286  std::map<int, int>::const_iterator it = phystolayer.begin();
1287  for( ; it != phystolayer.end() ; it++)
1288  {
1289  uint phys_label = it->first;
1290  uint layer_label = it->second;
1291  DEBUGL(FRMLAYTOT_D, "Adding formula of layer " << layer_label
1292  << " to domain of physical label " << phys_label);
1293  FormulaLayerPlaneWaveLayerGrad formula(sol(layer_label),
1294  sol(N+layer_label),
1295  kx,ky(layer_label),layer_label);
1296  this->set(phys_label, formula);
1297  }
1298  // Set, is used then as a its parent class
1299  }
1300 
1301  ~FormulaLayerPlaneWaveTotalGrad() {};
1302 
1303  };
1304 
1305 
1306 } // namespace concepts
1307 
1308 #endif // waveprop_sources_hh
concepts::FormulaIncPlaneWaveSource::coeff_b
ElementFormulaContainer< Cmplx > coeff_b
Definition: sources.hh:445
concepts::FormulaIncPlaneWaveSource::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition: sources.hh:441
concepts::Mapping
Basic class for a 2D or 3D map.
Definition: vectorsMatrices.hh:313
concepts::ComposeFormulaMatVec::clone
virtual ComposeFormulaMatVec * clone() const
Virtual constructor.
Definition: sources.hh:310
concepts::FormulaExpImag1D::info
virtual std::ostream & info(std::ostream &os) const
Definition: sources.hh:71
concepts::FormulaLayerPlaneWaveSourceGrad
Formula for plane wave source in layered structure.
Definition: sources.hh:737
concepts::FormulaExpImag2DGrad::FormulaExpImag2DGrad
FormulaExpImag2DGrad(const Real2d k, const Cmplx u=1.0, Real2d x0=Real2d(0, 0))
Definition: sources.hh:166
concepts::FormulaNormalOuterSP2D::operator()
virtual F operator()(const concepts::ElementWithCell< Real > &elm, const Real p, const Real t=0.0) const
Definition: pml_formula.h:212
concepts::FormulaLayerPlaneWaveLayer::FormulaLayerPlaneWaveLayer
FormulaLayerPlaneWaveLayer(const Cmplx a, const Cmplx b, const Real kx, const Real ky, const int layer_label=0)
Definition: sources.hh:973
concepts::FormulaIncPlaneWaveSource::u_inc
RCP< FormulaExpImag2D > u_inc
Definition: sources.hh:446
concepts::FormulaExpImag2DGrad::clone
virtual FormulaExpImag2DGrad * clone() const
Definition: sources.hh:170
conceptsThrowSimpleE
#define conceptsThrowSimpleE(msg)
Throws a equivalent to conceptsAssert(false, exc)
Definition: exceptions.hh:411
concepts::FormulaLayerPlaneWaveSource::info
virtual std::ostream & info(std::ostream &os) const
Definition: sources.hh:716
concepts::FormulaLayerPlaneWaveSource::epsilon_
std::unique_ptr< Vector< Real > > epsilon_
Definition: sources.hh:723
concepts::PiecewiseConstFormula
Piecewise constant function defined by the attribute of a cell.
Definition: formula.hh:84
concepts::FormulaIncPlaneWaveSource::u_inc
RCP< concepts::FormulaExpImag2D > u_inc
Definition: pml_formula.h:431
concepts::FormulaNormalOuterSP2D::FormulaNormalOuterSP2D
FormulaNormalOuterSP2D(const ElementFormulaContainer< Point< F, 2 > > vf, Real2d center=Real2d(0, 0), bool direction=OUTER)
Definition: sources.hh:229
concepts::FormulaExpImag2DRadialDer
Formula for radial derivative of .
Definition: pml_formula.h:91
concepts::FormulaIncPlaneWaveSource::genTMCoeff
static PiecewiseConstFormula< Cmplx > genTMCoeff(const MID &attToEps)
concepts::SuperLU
Direct sparse solver for unsymmetric matrices.
Definition: superLU.hh:70
concepts::ComposeFormulaVecEntry::operator()
virtual F operator()(const concepts::ElementWithCell< G > &elm, const Real p, const Real t=0.0) const
Definition: pml_formula.h:330
concepts::FormulaExpImag2DRadialDer::FormulaExpImag2DRadialDer
FormulaExpImag2DRadialDer(const Real2d k, const Cmplx u=1.0)
Definition: sources.hh:123
concepts::FormulaLayerPlaneWaveLayer::operator()
Cmplx operator()(const Real3d &p, const Real t=0.0) const
Application operator.
Definition: sources.hh:985
concepts::FormulaLayerPlaneWaveLayer
Formula for plane wave source in layered structure.
Definition: sources.hh:958
concepts::FormulaNormalOuterSP2D
Computes the scalar product <n, vf> of the normal n with a vector valued formula vf,...
Definition: pml_formula.h:195
concepts::FormulaLayerPlaneWaveSource::rho_
std::unique_ptr< Vector< Cmplx > > rho_
Definition: sources.hh:724
concepts::FormulaLayerPlaneWaveLayer::ky_
Real ky_
Wave number in the -direction.
Definition: sources.hh:968
concepts::ComposeFormulaMatVec::ComposeFormulaMatVec
ComposeFormulaMatVec(RCP< const ElementFormula< Mapping< F, DIM >, G > > A, RCP< const ElementFormula< Point< F, DIM >, G > > vf)
Definition: sources.hh:302
concepts::FormulaIncPlaneWaveSource::genTECoeff
static PiecewiseConstFormula< Cmplx > genTECoeff(const MID &attToEps)
concepts::Real2d
Point< Real, 2 > Real2d
Definition: vectorsMatricesForward.hh:18
concepts::FormulaLayerPlaneWaveSource::d_
std::unique_ptr< Vector< Real > > d_
Definition: sources.hh:722
std::arg
concepts::Real arg(const concepts::Point< concepts::Real, 2 > &p)
Returns the phase angle of a real 2D vector.
concepts::Cell
A cell in a mesh consist of topological information (neighbours, connectivity, orientation) and geome...
Definition: cell.hh:39
concepts::FormulaLayerPlaneWaveLayerGrad::FormulaLayerPlaneWaveLayerGrad
FormulaLayerPlaneWaveLayerGrad(const Cmplx a, const Cmplx b, const Real kx, const Real ky, const int layer_label=0)
Definition: sources.hh:1036
concepts::FormulaLayerPlaneWaveSourceGrad::B_
std::unique_ptr< Vector< Cmplx > > B_
Definition: sources.hh:943
concepts::ElementWithCell::elemMap
Real3d elemMap(const Real coord_local) const
Definition: element.hh:86
concepts::FormulaLayerPlaneWaveSourceGrad::clone
virtual FormulaLayerPlaneWaveSourceGrad * clone() const
Definition: sources.hh:889
concepts::FormulaLayerPlaneWaveSource::operator()
virtual Cmplx operator()(const Real3d &p, const Real t=0.0) const
Application operator.
Definition: sources.hh:683
concepts::FormulaLayerPlaneWaveSourceGrad::FormulaLayerPlaneWaveSourceGrad
FormulaLayerPlaneWaveSourceGrad(Vector< Real > &epsilon, Vector< Real > &d, const Real &kx, const Real omega)
Definition: sources.hh:755
concepts::FormulaIncPlaneWaveSource::coeff_a
concepts::ElementFormulaContainer< Cmplx > coeff_a
Definition: pml_formula.h:428
concepts::Formula
Interface for a formula.
Definition: lform.hh:18
concepts::FormulaExpImag1D::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Application operator.
Definition: pml_formula.h:31
conceptsAssert
#define conceptsAssert(cond, exc)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:394
concepts::Edge2d
A 1D cell: edge in 2D.
Definition: cell1D.hh:189
concepts::Edge2d::jacobian
Real2d jacobian(Real x) const
Returns the Jacobian.
FRMLAYTOT_D
#define FRMLAYTOT_D
Definition: sources.hh:43
concepts::FormulaExpImag2D::clone
virtual FormulaExpImag2D * clone() const
Definition: sources.hh:88
FRMLAYPW_D
#define FRMLAYPW_D
Definition: sources.hh:42
concepts::ComposeFormulaMatVec::operator()
virtual concepts::Point< F, DIM > operator()(const concepts::ElementWithCell< G > &elm, const Real p, const Real t=0.0) const
Definition: pml_formula.h:284
concepts::ComposeFormulaMatVec::ComposeFormulaMatVec
ComposeFormulaMatVec(RCP< const ElementFormula< concepts::Mapping< F, DIM >, G > > A, RCP< const ElementFormula< concepts::Point< F, DIM >, G > > vf)
Definition: pml_formula.h:272
concepts::Vector::size
uint size() const
Definition: vector.hh:147
concepts::FormulaIncPlaneWaveSource::FormulaIncPlaneWaveSource
FormulaIncPlaneWaveSource(ElementFormulaContainer< Cmplx > coeff_a, ElementFormulaContainer< Cmplx > coeff_b, RCP< concepts::FormulaExpImag2D > u_inc)
Definition: pml_formula.h:362
concepts::FormulaIncPlaneWaveSource::MID
std::map< int, double > MID
Definition: sources.hh:392
concepts::FormulaLayerPlaneWaveTotal
Formula for plane wave source in layered structure.
Definition: sources.hh:1085
concepts::FormulaExpImag2DRadialDer::info
virtual std::ostream & info(std::ostream &os) const
Definition: sources.hh:150
DEBUGL
#define DEBUGL(doit, msg)
concepts::FormulaNormalOuterSP2D::clone
virtual FormulaNormalOuterSP2D * clone() const
Virtual constructor.
Definition: sources.hh:237
concepts::FormulaLayerPlaneWaveSourceGrad::A_
std::unique_ptr< Vector< Cmplx > > A_
Definition: sources.hh:942
concepts::FormulaLayerPlaneWaveLayerGrad
Formula for gradient of plane wave source in layered structure.
Definition: sources.hh:1021
concepts::FormulaLayerPlaneWaveLayerGrad::b_
Cmplx b_
Complex coefficient associated to the plane wave .
Definition: sources.hh:1027
concepts::FormulaLayerPlaneWaveLayerGrad::ky_
Real ky_
Wave number in the -direction.
Definition: sources.hh:1031
concepts::FormulaLayerPlaneWaveSource::B_
std::unique_ptr< Vector< Cmplx > > B_
Definition: sources.hh:727
concepts::FormulaLayerPlaneWaveSource
Formula for plane wave source in layered structure.
Definition: sources.hh:505
concepts::PiecewiseFormula< Cmplx >::set
void set(const uint attrib, const Formula< Cmplx > &formula)
Set the constant in an particular cell.
concepts::FormulaLayerPlaneWaveSource::clone
virtual FormulaLayerPlaneWaveSource * clone() const
Definition: sources.hh:675
concepts::FormulaExpImag2DGrad::operator()
virtual Cmplx2d operator()(const Real p, const Real t=0.0) const
Application operator.
Definition: pml_formula.h:147
concepts::FormulaExpImag2D::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Application operator.
Definition: pml_formula.h:63
concepts::Assertion
Exception class for assertions.
Definition: exceptions.hh:258
concepts::FormulaLayerPlaneWaveTotalGrad::FormulaLayerPlaneWaveTotalGrad
FormulaLayerPlaneWaveTotalGrad(const Vector< Real > &epsilon, const Vector< Real > &d, const Real kx, const Real omega, const std::map< int, int > phystolayer)
Constructor.
Definition: sources.hh:1215
concepts::FormulaLayerPlaneWaveLayer::info
std::ostream & info(std::ostream &os) const
Definition: sources.hh:1003
concepts::FormulaLayerPlaneWaveLayerGrad::info
std::ostream & info(std::ostream &os) const
Definition: sources.hh:1068
concepts::Point::l2
Real l2() const
Returns the Euclidian norm of the vector.
Definition: vectorsMatrices.hh:182
concepts::ComposeFormulaVecEntry::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition: sources.hh:378
concepts::FormulaExpImag2DGrad::info
virtual std::ostream & info(std::ostream &os) const
Definition: sources.hh:197
concepts::ElementFormula
Interface for a formula defined element by element.
Definition: elementFormula.hh:35
concepts::FormulaIncPlaneWaveSource::coeff_b
concepts::ElementFormulaContainer< Cmplx > coeff_b
Definition: pml_formula.h:428
concepts::FormulaLayerPlaneWaveSourceGrad::rho_
std::unique_ptr< Vector< Cmplx > > rho_
Definition: sources.hh:940
concepts::FormulaLayerPlaneWaveSourceGrad::d_
std::unique_ptr< Vector< Real > > d_
Definition: sources.hh:938
concepts::FormulaExpImag1D::FormulaExpImag1D
FormulaExpImag1D(const Real k, const Cmplx u=1.0, Real x0=0.0)
Definition: sources.hh:53
concepts::PiecewiseFormula
Piecewise defined function defined by attributes.
Definition: formula.hh:218
concepts::ComposeFormulaVecEntry::ComposeFormulaVecEntry
ComposeFormulaVecEntry(RCP< const ElementFormula< Point< F, DIM >, G > > vf, int index)
Definition: sources.hh:348
concepts::Point::dot
F dot(const Point< F, dim > &b) const
Inner product, i.e. for complex arguments: this * conjugate(b)
concepts::FormulaIncPlaneWaveSource::genTMCoeff
static PiecewiseConstFormula< Cmplx > genTMCoeff(const MID &attToEps)
Definition: pml_formula.h:434
concepts::FormulaIncPlaneWaveSource::clone
virtual FormulaIncPlaneWaveSource * clone() const
Virtual constructor.
Definition: sources.hh:402
concepts::ComposeFormulaMatVec::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition: sources.hh:334
concepts::ComposeFormulaVecEntry::ComposeFormulaVecEntry
ComposeFormulaVecEntry(RCP< const ElementFormula< concepts::Point< F, DIM >, G > > vf, int index)
Definition: pml_formula.h:318
concepts::ComposeFormulaVecEntry::clone
virtual ComposeFormulaVecEntry * clone() const
Virtual constructor.
Definition: sources.hh:356
concepts::FormulaNormalOuterSP2D::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition: sources.hh:283
concepts::FormulaNormalOuterSP2D::FormulaNormalOuterSP2D
FormulaNormalOuterSP2D(const ElementFormulaContainer< concepts::Point< F, 2 > > vf, Real2d center=Real2d(0, 0), bool direction=OUTER)
Definition: pml_formula.h:200
concepts::FormulaNormalOuterSP2D::vf
const ElementFormulaContainer< concepts::Point< F, 2 > > vf
Definition: pml_formula.h:258
concepts::FormulaLayerPlaneWaveLayer::a_
Cmplx a_
Complex coefficient associated to the plane wave .
Definition: sources.hh:962
concepts::FormulaLayerPlaneWaveTotalGrad
Formula for gradient of plane wave source in layered structure.
Definition: sources.hh:1198
concepts::Cmplx
std::complex< Real > Cmplx
Type for a complex number. It also depends on the setting of Real.
Definition: typedefs.hh:39
concepts::ComposeFormulaVecEntry::vf
RCP< const ElementFormula< concepts::Point< F, DIM >, G > > vf
Definition: pml_formula.h:352
concepts::FormulaLayerPlaneWaveLayerGrad::operator()
Cmplx2d operator()(const Real3d &p, const Real t=0.0) const
Application operator.
Definition: sources.hh:1048
concepts::FormulaLayerPlaneWaveLayerGrad::kx_
Real kx_
Wave number in the -direction.
Definition: sources.hh:1029
concepts::FormulaLayerPlaneWaveSourceGrad::info
virtual std::ostream & info(std::ostream &os) const
Definition: sources.hh:932
concepts::FormulaLayerPlaneWaveLayer::clone
FormulaLayerPlaneWaveLayer * clone() const
Definition: sources.hh:981
concepts::FormulaExpImag2D::FormulaExpImag2D
FormulaExpImag2D(const Real2d k, const Cmplx u=1.0, Real2d x0=0.0)
Definition: sources.hh:85
concepts::FormulaLayerPlaneWaveLayerGrad::clone
FormulaLayerPlaneWaveLayerGrad * clone() const
Definition: sources.hh:1044
concepts::Point::ortho
Point< F, dim > & ortho(const Point< F, dim > &a)
Change vector to the by 90 degrees (clockwise) rotated vector a (only 2D)
concepts::ComposeFormulaVecEntry::vf
RCP< const ElementFormula< Point< F, DIM >, G > > vf
Definition: sources.hh:382
concepts::RCP
Reference-counting pointer.
Definition: bf_iddiv.hh:15
concepts::FormulaExpImag2DRadialDer::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Application operator.
Definition: pml_formula.h:101
concepts::ComposeFormulaMatVec::A
RCP< const ElementFormula< Mapping< F, DIM >, G > > A
Definition: sources.hh:338
concepts::ComposeFormulaMatVec
Computes the Matrix-vector product A * vf, where A is a matrix valued formula and vf a vector valued ...
Definition: pml_formula.h:270
concepts::Cmplx2d
Point< Cmplx, 2 > Cmplx2d
Definition: vectorsMatricesForward.hh:21
concepts::FormulaNormalOuterSP2D::vf
const ElementFormulaContainer< Point< F, 2 > > vf
Definition: sources.hh:287
concepts::FormulaLayerPlaneWaveLayerGrad::FormulaLayerPlaneWaveLayerGrad
FormulaLayerPlaneWaveLayerGrad(const FormulaLayerPlaneWaveLayerGrad &other)
Definition: sources.hh:1039
concepts::ElementWithCell::cell
virtual const Cell & cell() const =0
Returns the cell on which the element is built.
concepts::FormulaLayerPlaneWaveLayer::b_
Cmplx b_
Complex coefficient associated to the plane wave .
Definition: sources.hh:964
concepts::Mumps
MUMPS : MUltifrontal Massively Parallel sparse direct Solver.
Definition: mumps.hh:72
concepts::FormulaIncPlaneWaveSource::genTECoeff
static PiecewiseConstFormula< Cmplx > genTECoeff(const MID &attToEps)
Definition: pml_formula.h:463
concepts::FormulaIncPlaneWaveSource::FormulaIncPlaneWaveSource
FormulaIncPlaneWaveSource(ElementFormulaContainer< Cmplx > coeff_a, ElementFormulaContainer< Cmplx > coeff_b, RCP< FormulaExpImag2D > u_inc)
Definition: sources.hh:394
concepts::FormulaExpImag2DGrad
Formula for gradient of a plane wave.
Definition: pml_formula.h:137
concepts::FormulaExpImag1D::clone
virtual FormulaExpImag1D * clone() const
Definition: sources.hh:56
concepts::FormulaLayerPlaneWaveSource::ky_
std::unique_ptr< Vector< Real > > ky_
Definition: sources.hh:725
concepts::GMRes
Solves a system of linear equations with general minimal residuals (GMRes).
Definition: gmres.hh:24
concepts::FormulaLayerPlaneWaveSource::FormulaLayerPlaneWaveSource
FormulaLayerPlaneWaveSource(Vector< Real > &epsilon, Vector< Real > &d, const Real &kx, const Real omega)
Definition: sources.hh:525
concepts::FormulaExpImag2DRadialDer::clone
virtual FormulaExpImag2DRadialDer * clone() const
Definition: sources.hh:126
concepts::FormulaExpImag2D::info
virtual std::ostream & info(std::ostream &os) const
Definition: sources.hh:108
conceptsAssert3
#define conceptsAssert3(cond, exc, msg)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:442
concepts::typeOf
std::string typeOf(const T &t)
Return the typeid name of a class object.
Definition: output.hh:43
concepts::FormulaLayerPlaneWaveLayerGrad::a_
Cmplx a_
Complex coefficient associated to the plane wave .
Definition: sources.hh:1025
concepts::ComposeFormulaMatVec::vf
RCP< const ElementFormula< concepts::Point< F, DIM >, G > > vf
Definition: pml_formula.h:309
concepts::ComposeFormulaMatVec::A
RCP< const ElementFormula< concepts::Mapping< F, DIM >, G > > A
Definition: pml_formula.h:308
concepts::FormulaLayerPlaneWaveLayer::kx_
Real kx_
Wave number in the -direction.
Definition: sources.hh:966
concepts::FormulaLayerPlaneWaveSourceGrad::epsilon_
std::unique_ptr< Vector< Real > > epsilon_
Definition: sources.hh:939
concepts::ComposeFormulaMatVec::vf
RCP< const ElementFormula< Point< F, DIM >, G > > vf
Definition: sources.hh:339
concepts::FormulaLayerPlaneWaveSourceGrad::operator()
virtual Cmplx2d operator()(const Real3d &p, const Real t=0.0) const
Application operator.
Definition: sources.hh:897
concepts::FormulaLayerPlaneWaveSource::A_
std::unique_ptr< Vector< Cmplx > > A_
Definition: sources.hh:726
concepts::Real
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
concepts::FormulaLayerPlaneWaveSourceGrad::ky_
std::unique_ptr< Vector< Real > > ky_
Definition: sources.hh:941
concepts::FormulaIncPlaneWaveSource::operator()
virtual Cmplx operator()(const ElementWithCell< Real > &elm, const Real3d &p, const Real t=0.0) const
Definition: pml_formula.h:376
concepts
Basic namespace for Concepts-2.
Definition: pml_formula.h:16
concepts::FormulaLayerPlaneWaveLayer::FormulaLayerPlaneWaveLayer
FormulaLayerPlaneWaveLayer(const FormulaLayerPlaneWaveLayer &other)
Definition: sources.hh:976
concepts::FormulaLayerPlaneWaveTotal::FormulaLayerPlaneWaveTotal
FormulaLayerPlaneWaveTotal(const Vector< Real > &epsilon, const Vector< Real > &d, const Real kx, const Real omega, const std::map< int, int > phystolayer)
Constructor.
Definition: sources.hh:1102
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich