You are here: Concepts>Concepts Web>Class documentation

Class documentation of Concepts

Loading...
Searching...
No Matches
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
45namespace 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
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::ComposeFormulaMatVec::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:334
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::clone
virtual ComposeFormulaMatVec * clone() const
Virtual copy constructor.
Definition sources.hh:310
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::ComposeFormulaVecEntry::clone
virtual ComposeFormulaVecEntry * clone() const
Virtual copy constructor.
Definition sources.hh:356
concepts::ComposeFormulaVecEntry::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:378
concepts::Edge2d::jacobian
Real2d jacobian(Real x) const
concepts::ElementWithCell::cell
virtual const Cell & cell() const =0
Returns the cell on which the element is built.
concepts::FormulaExpImag1D::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Definition pml_formula.h:31
concepts::FormulaExpImag1D::clone
virtual FormulaExpImag1D * clone() const
Virtual copy constructor.
Definition sources.hh:56
concepts::FormulaExpImag1D::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:71
concepts::FormulaExpImag2DGrad::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:197
concepts::FormulaExpImag2DGrad::operator()
virtual Cmplx2d operator()(const Real p, const Real t=0.0) const
Definition pml_formula.h:147
concepts::FormulaExpImag2DGrad::clone
virtual FormulaExpImag2DGrad * clone() const
Virtual copy constructor.
Definition sources.hh:170
concepts::FormulaExpImag2DRadialDer::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Definition pml_formula.h:101
concepts::FormulaExpImag2DRadialDer::clone
virtual FormulaExpImag2DRadialDer * clone() const
Virtual copy constructor.
Definition sources.hh:126
concepts::FormulaExpImag2DRadialDer::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:150
concepts::FormulaExpImag2D::clone
virtual FormulaExpImag2D * clone() const
Virtual copy constructor.
Definition sources.hh:88
concepts::FormulaExpImag2D::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:108
concepts::FormulaExpImag2D::operator()
virtual Cmplx operator()(const Real p, const Real t=0.0) const
Definition pml_formula.h:63
concepts::FormulaIncPlaneWaveSource::clone
virtual FormulaIncPlaneWaveSource * clone() const
Virtual copy constructor.
Definition sources.hh:402
concepts::FormulaIncPlaneWaveSource::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:441
concepts::FormulaLayerPlaneWaveLayerGrad::a_
Cmplx a_
Complex coefficient associated to the plane wave .
Definition sources.hh:1025
concepts::FormulaLayerPlaneWaveLayerGrad::clone
FormulaLayerPlaneWaveLayerGrad * clone() const
Virtual copy constructor.
Definition sources.hh:1044
concepts::FormulaLayerPlaneWaveLayerGrad::b_
Cmplx b_
Complex coefficient associated to the plane wave .
Definition sources.hh:1027
concepts::FormulaLayerPlaneWaveLayerGrad::info
std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:1068
concepts::FormulaLayerPlaneWaveLayerGrad::kx_
Real kx_
Wave number in the -direction.
Definition sources.hh:1029
concepts::FormulaLayerPlaneWaveLayerGrad::ky_
Real ky_
Wave number in the -direction.
Definition sources.hh:1031
concepts::FormulaLayerPlaneWaveLayer::kx_
Real kx_
Wave number in the -direction.
Definition sources.hh:966
concepts::FormulaLayerPlaneWaveLayer::clone
FormulaLayerPlaneWaveLayer * clone() const
Virtual copy constructor.
Definition sources.hh:981
concepts::FormulaLayerPlaneWaveLayer::info
std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:1003
concepts::FormulaLayerPlaneWaveLayer::ky_
Real ky_
Wave number in the -direction.
Definition sources.hh:968
concepts::FormulaLayerPlaneWaveLayer::b_
Cmplx b_
Complex coefficient associated to the plane wave .
Definition sources.hh:964
concepts::FormulaLayerPlaneWaveLayer::a_
Cmplx a_
Complex coefficient associated to the plane wave .
Definition sources.hh:962
concepts::FormulaLayerPlaneWaveSourceGrad::clone
virtual FormulaLayerPlaneWaveSourceGrad * clone() const
Virtual copy constructor.
Definition sources.hh:889
concepts::FormulaLayerPlaneWaveSourceGrad::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:932
concepts::FormulaLayerPlaneWaveSource::clone
virtual FormulaLayerPlaneWaveSource * clone() const
Virtual copy constructor.
Definition sources.hh:675
concepts::FormulaLayerPlaneWaveSource::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:716
concepts::FormulaLayerPlaneWaveTotalGrad::FormulaLayerPlaneWaveTotalGrad
FormulaLayerPlaneWaveTotalGrad(const Vector< Real > &epsilon, const Vector< Real > &d, const Real kx, const Real omega, const std::map< int, int > phystolayer)
Definition sources.hh:1215
concepts::FormulaLayerPlaneWaveTotal::FormulaLayerPlaneWaveTotal
FormulaLayerPlaneWaveTotal(const Vector< Real > &epsilon, const Vector< Real > &d, const Real kx, const Real omega, const std::map< int, int > phystolayer)
Definition sources.hh:1102
concepts::FormulaNormalOuterSP2D::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition sources.hh:283
concepts::FormulaNormalOuterSP2D::clone
virtual FormulaNormalOuterSP2D * clone() const
Virtual copy constructor.
Definition sources.hh:237
concepts::PiecewiseFormula< Cmplx >::set
void set(const uint attrib, const Formula< Cmplx > &formula)
concepts::Point::ortho
Point< F, dim > & ortho(const Point< F, dim > &a)
concepts::Point::dot
F dot(const Point< F, dim > &b) const
Inner product, i.e. for complex arguments: this * conjugate(b)
concepts::Point::l2
Real l2() const
Returns the Euclidian norm of the vector.
Definition vectorsMatrices.hh:182
conceptsAssert
#define conceptsAssert(cond, exc)
Definition exceptions.hh:394
DEBUGL
#define DEBUGL(doit, msg)
Definition debug.hh:40
conceptsAssert3
#define conceptsAssert3(cond, exc, msg)
Definition exceptions.hh:442
conceptsThrowSimpleE
#define conceptsThrowSimpleE(msg)
Definition exceptions.hh:411
concepts::typeOf
std::string typeOf(const T &t)
Definition output.hh:43
concepts::Real
double Real
Definition typedefs.hh:39
concepts::makeSet
Set< F > makeSet(uint n, const F &first,...)
Definition set.hh:320
concepts::Cmplx
std::complex< Real > Cmplx
Type for a complex number. It also depends on the setting of Real.
Definition typedefs.hh:42
Generated on Wed Sep 13 2023 21:06:25 for Concepts by doxygen 1.9.8