v0.10.0
AnalyticalSolutions.hpp
Go to the documentation of this file.
1 /** \file AnalyticalSolutions.hpp
2  \ingroup mofem_helmholtz_elem
3 
4  Analytical solutions
5 
6  */
7 
8 /* This file is part of MoFEM.
9  * MoFEM is free software: you can redistribute it and/or modify it under
10  * the terms of the GNU Lesser General Public License as published by the
11  * Free Software Foundation, either version 3 of the License, or (at your
12  * option) any later version.
13 
14  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17  * License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
21 */
22 
23 /** \brief Generic structure for analytical function
24  \ingroup mofem_helmholtz_elem
25  \bug point source example not implemented.
26 */
27 
28 
30 
32 
33  virtual vector<VectorDouble >& operator()(double x, double y, double z) = 0;
35 
36 };
37 
38 /** List of analytical solution
39  \ingroup mofem_helmholtz_elem
40 */
50 };
51 
52 /** Line command list of analytical solutions
53  \ingroup mofem_helmholtz_elem
54 */
55 const char *analytical_solution_types[] = {
56  "hard_sphere_incident_wave",
57  "soft_sphere_incident_wave",
58  "plane_wave",
59  "hard_cylinder_scatter_wave",
60  "soft_cylinder_scatter_wave",
61  "singular_scatter_wave",
62  "incident_wave",
63  "no_analytical_solution"
64 };
65 
66 // **** Analytical solutions ****
67 
68 /** Incident wave
69  \ingroup mofem_helmholtz_elem
70 
71 
72  Equation from:
73  Ihlenburg,Finite element analysis of acoustic scattering Springer Science & Business Media.
74 
75  Some details can be found here:
76  <http://ansol.us/Products/Coustyx/Validation/MultiDomain/Scattering/PlaneWave/HardSphere/Downloads/dataset_description.pdf>
77 
78  \f[
79  p_\textrm{inc} = \exp(ikd \cdot \mathbf{x})
80  \f]
81  FIXME
82  Rayleigh with attenuation of the incident wave need to be implemented.
83 
84  */
86 
87  vector<VectorDouble > rEsult;
88  double wAvenumber;
89  double aTtenuation;
90  double rAdius;
92  double fRequency;
93  double sPeed;
97  double amplitudeReal; ///< The real amplitude of the incident wave
98  double amplitudeImag;
99  double pHase;
101 
103  IncidentWave(double wavenumber,VectorDouble d,bool is_rayleigh = false,double attenuation = 0,double radius = 0, double leaky_saw = 0,double frequency = 0,double velocity = 0,double transmission_coefficient = 1,double r_amplitude = 1,double i_amplitude = 0,double phase = 0):
104  wAvenumber(wavenumber),
105  dIrection(d),
106  isRayleigh(is_rayleigh),
107  aTtenuation(attenuation),
108  rAdius(radius),
109  complexWaveNumber(leaky_saw),
110  fRequency(frequency),
111  sPeed(velocity),
112  transmissionCoefficient(transmission_coefficient),
113  amplitudeReal(r_amplitude),
114  amplitudeImag(i_amplitude),
115  pHase(phase)
116  {}
117 
119 
120  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
121 
122  const complex< double > i( 0.0, 1.0 );
123  complex< double > result = 0.0;
124  cOordinate.resize(3);
125  cOordinate[0] = x;
126  cOordinate[1] = y;
127  cOordinate[2] = z;
128  complex< double > attenuationR(0,0);
129  complex< double > waveNumber = wAvenumber + i*complexWaveNumber;
130  if(isRayleigh) {
131  /* speed of SAW over speed of longitudinal wave */
132  complex< double > vs_over_v1 = (2.0*M_PI*fRequency/waveNumber) / sPeed;
133  attenuationR = sqrt(1.0 - pow(vs_over_v1,2.0))*(y);
134  /* attenuationR of Rayleigh inside fluid */
135 
136  /* double totalViscosity = (4/3) * 0.001139 + 0.0034 = 0.004900 */
137  /* sec(Rayleigh angle) = sec(22 degree) = 1.0790 */
138  // double oMega = 2*M_PI*fRequency;
139  // double dEnsity = 998;
140  // attenuationR = (0.0049*pow(oMega,2.0)) / (dEnsity*pow(sPeed,3.0));
141  // attenuationR *= - 1.0790*abs(y);
142  }
143 
144  if(aTtenuation != 0) {
145  double xOriginal = - std::sqrt(std::pow(rAdius,2.0) - std::pow(z,2.0));
146  aTtenuation *= (x - xOriginal);
147  }
148 
149 
150  result = transmissionCoefficient*(amplitudeReal+i*amplitudeImag)*exp(i*waveNumber*inner_prod(dIrection,cOordinate)+i*pHase+waveNumber*attenuationR - aTtenuation);
151 
152 
153  rEsult.resize(2);
154  rEsult[REAL].resize(1);
155  (rEsult[REAL])[0] = std::real(result);
156  rEsult[IMAG].resize(1);
157  (rEsult[IMAG])[0] = std::imag(result);
158 
159  return rEsult;
160 
161  }
162 
163 };
164 
166 
167  vector<VectorDouble > rEsult;
168  double wAvenumber;
169  double aTtenuation;
170  double rAdius;
172  double fRequency;
173  double sPeed;
177  double amplitudeReal; ///< The real amplitude of the incident wave
180 
182  IncidentWavePointSource(double wavenumber,VectorDouble d,bool is_rayleigh = false,double attenuation = 0,double radius = 0,double leaky_saw = 0,double frequency = 0,double velocity = 0,double transmission_coefficient = 1,double r_amplitude = 1,double i_amplitude = 0):
183  wAvenumber(wavenumber),
184  aTtenuation(attenuation),
185  rAdius(radius),
186  complexWaveNumber(leaky_saw),
187  fRequency(frequency),
188  sPeed(velocity),
189  transmissionCoefficient(transmission_coefficient),
191  isRayleigh(is_rayleigh),
192  amplitudeReal(r_amplitude),
193  amplitudeImag(i_amplitude)
194  {}
195 
197 
198  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
199 
200  const complex< double > i( 0.0, 1.0 );
201  complex< double > result = 0.0;
202  pOsition.resize(3);
203 
204  pOsition[0] = (x-sourceCoordinate[0])*(x-sourceCoordinate[0]);
205  pOsition[1] = (y-sourceCoordinate[1])*(y-sourceCoordinate[1]);
206  pOsition[2] = (z-sourceCoordinate[2])*(z-sourceCoordinate[2]);
207  double absPoint = sqrt(pOsition[0]+pOsition[1]+pOsition[2]);
208  complex<double> amplitude = amplitudeReal+i*amplitudeImag;
209  complex< double > attenuationR(0,0);
210  complex< double > waveNumber = wAvenumber + i*complexWaveNumber;
211  if(isRayleigh) {
212  /* speed of SAW over speed of longitudinal wave */
213  // complex< double > vs_over_v1 = (2.0*M_PI*fRequency/waveNumber) / sPeed;
214  // attenuationR = -sqrt(1.0 - pow(vs_over_v1,2.0))*abs(y);
215  /* attenuationR of Rayleigh inside fluid */
216  double oMega = 2*M_PI*fRequency;
217  /* double totalViscosity = (4/3) * 0.001139 + 0.0034 = 0.004900 */
218  /* sec(Rayleigh angle) = sec(22 degree) = 1.0790 */
219  double dEnsity = 998;
220  attenuationR = (0.0049*pow(oMega,2.0)) / (dEnsity*pow(sPeed,3.0));
221  attenuationR *= - 1.0790*abs(y);
222  }
223 
224  if(aTtenuation != 0) {
225  double xOriginal = - std::sqrt(std::pow(rAdius,2.0) - std::pow(z,2.0));
226  aTtenuation *= (x - xOriginal);
227  }
228 
229  result = transmissionCoefficient*amplitude*(exp(-i*waveNumber*absPoint+waveNumber*attenuationR - aTtenuation)) / (4*M_PI*absPoint);
230 
231  rEsult.resize(2);
232  rEsult[REAL].resize(1);
233  (rEsult[REAL])[0] = std::real(result);
234  rEsult[IMAG].resize(1);
235  (rEsult[IMAG])[0] = std::imag(result);
236 
237  return rEsult;
238 
239  }
240 
241 };
242 
243 
244 #ifdef KISS_FFT_H
245 
246 /**
247 
248  Incident wave, i.t. is inverse Fourier transform evaluated at arbitrary
249  spatial points.
250 
251  \f[
252  \left. \left\{ \frac{1}{n} (A_{0} e^{ik \mathbf{d} \cdot \mathbf{x} + i \phi})
253  \right\}
254  \f]
255  where \f$\phi\f$ is
256  \f[
257 
258 
259  \f]
260 
261  */
262 struct IncidentWaveDFT: public GenericAnalyticalSolution {
263 
264  vector<VectorDouble > rEsult;
265  double signalLength;
266  double signalDuration;
267  double aTtenuation;
268  double rAdius;
269  double complexWaveNumber;
270  double fRequency;
271  double sPeed;
272  double transmissionCoefficient;
273  double dEnsity;
274  VectorDouble dIrection;
275  boost::shared_array<kiss_fft_cpx> complexOut;
276  int sizeFreq;
277  int sizeTime;
278  int timeStep;
279  VectorDouble cOordinate;
280  bool isRadiation;
281  bool isRayleigh;
282 
283  IncidentWaveDFT() {}
284  IncidentWaveDFT(
285  double signal_length,
286  double signal_duration,
287  VectorDouble &d,
288  double attenuation,
289  double radius,
290  double leaky_saw,
291  double frequency,
292  double velocity,
293  double transmission_coefficient,
294  double density,
295  boost::shared_array<kiss_fft_cpx> complex_out,
296  int size_frequency,
297  int size_time,
298  int time_step,
299  bool is_radiation,
300  bool is_rayleigh
301  ):
302  signalLength(signal_length),
303  signalDuration(signal_duration),
304  dIrection(d),
305  aTtenuation(attenuation),
306  rAdius(radius),
307  complexWaveNumber(leaky_saw),
308  fRequency(frequency),
309  sPeed(velocity),
310  transmissionCoefficient(transmission_coefficient),
311  dEnsity(density),
312  complexOut(complex_out),
313  sizeFreq(size_frequency),
314  sizeTime(size_time),
315  timeStep(time_step),
316  isRadiation(is_radiation),
317  isRayleigh(is_rayleigh)
318  {}
319 
320  ~IncidentWaveDFT() {}
321 
322  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
323 
324  const complex< double > i( 0.0, 1.0 );
325  cOordinate.resize(3);
326  cOordinate[0] = x;
327  cOordinate[1] = y;
328  cOordinate[2] = z;
329  complex< double > attenuationR(0,0);
330  if(isRayleigh) {
331 
332  /* speed of SAW over speed of longitudinal wave */
333  complex< double > vs_over_v1 = (signalLength*fRequency) / sPeed;
334  attenuationR = -sqrt(1.0 - pow(vs_over_v1,2.0))*abs(y);
335 
336  /* attenuationR of Rayleigh inside fluid */
337  // double oMega = 2*M_PI*fRequency;
338  /* double totalViscosity = (4/3) * 0.001139 + 0.0034; */
339  /* sec(Rayleigh angle) = sec(22 degree) = 1.0790 */
340  // attenuationR = (0.0049*pow(oMega,2.0)) / (dEnsity*pow(sPeed,3.0));
341  // attenuationR *= - 1.0790*abs(y);
342 
343  }
344 
345  if(aTtenuation != 0) {
346  double xOriginal = - std::sqrt(std::abs(std::pow(rAdius,2.0) - std::pow(z,2.0)));
347  aTtenuation *= (x - xOriginal);
348  }
349 
350  //FIXME can apply zero padding to signal length and sIze in this section when deal with large number of time steps.
351  complex< double > result = 0.0;
352  for(int f = 0;f < sizeTime && sizeFreq;f++) {
353  double speed = signalLength/signalDuration;
354  /* 2*M_PI/signalLength : 2*M_PI*f/signalLength */
355  complex< double > wave_number = (2*M_PI*f/signalLength) + i*complexWaveNumber;
356  double delta_t = signalDuration/sizeTime;
357  double distance = speed*delta_t*timeStep;
358  double phase= 2*M_PI*f*(distance/signalLength);
359 
360  // cerr << "\n exp(+ i*wave_number*inner_prod(dIrection,cOordinate)+i*phase) = \n " << exp(+ i*wave_number*inner_prod(dIrection,cOordinate)+i*phase) << endl;
361 
362  if(isRadiation) {
363  result += (complexOut[f].r+i*complexOut[f].i)*exp(+ i*wave_number*inner_prod(dIrection,cOordinate)+i*phase+wave_number*attenuationR - aTtenuation);
364  } else{
365  result += (complexOut[f].r+i*complexOut[f].i)*exp(- i*wave_number*inner_prod(dIrection,cOordinate)+i*phase+wave_number*attenuationR - aTtenuation);
366  }
367  }
368 
369  result /= sizeTime;
370  result *= transmissionCoefficient;
371 
372  // cerr << "\n result = \n " << result << endl;
373 
374  rEsult.resize(2);
375  rEsult[REAL].resize(1);
376  (rEsult[REAL])[0] = std::real(result);
377  rEsult[IMAG].resize(1);
378  (rEsult[IMAG])[0] = std::imag(result);
379 
380  return rEsult;
381 
382  }
383 
384 };
385 
386 struct IncidentWavePointSourceDFT: public GenericAnalyticalSolution {
387 
388  vector<VectorDouble > rEsult;
389  double signalLength;
390  double signalDuration;
391  double aTtenuation;
392  double rAdius;
393  double complexWaveNumber;
394  double fRequency;
395  double sPeed;
396  double transmissionCoefficient;
397  double dEnsity;
398  boost::shared_array<kiss_fft_cpx> complexOut;
399  int sizeFreq;
400  int sizeTime;
401  int timeStep;
402  bool isRadiation;
403  bool isRayleigh;
404  VectorDouble pOsition;
405  VectorDouble sourceCoordinate;
406 
407  IncidentWavePointSourceDFT() {}
408  IncidentWavePointSourceDFT(
409  double signal_length,
410  double signal_duration,
411  VectorDouble &d,
412  double attenuation,
413  double radius,
414  double leaky_saw,
415  double frequency,
416  double velocity,
417  double transmission_coefficient,
418  double density,
419  boost::shared_array<kiss_fft_cpx> complex_out,
420  int size_frequency,
421  int size_time,
422  int time_step,
423  bool is_radiation,
424  bool is_rayleigh
425  ):
426  signalLength(signal_length),
427  signalDuration(signal_duration),
428  sourceCoordinate(d),
429  aTtenuation(attenuation),
430  rAdius(radius),
431  complexWaveNumber(leaky_saw),
432  fRequency(frequency),
433  sPeed(velocity),
434  transmissionCoefficient(transmission_coefficient),
435  dEnsity(density),
436  complexOut(complex_out),
437  sizeFreq(size_frequency),
438  sizeTime(size_time),
439  timeStep(time_step),
440  isRadiation(is_radiation),
441  isRayleigh(is_rayleigh)
442  {}
443 
444  ~IncidentWavePointSourceDFT() {}
445 
446  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
447 
448  const complex< double > i( 0.0, 1.0 );
449  complex< double > result = 0.0;
450  pOsition.resize(3);
451  pOsition[0] = (x-sourceCoordinate[0])*(x-sourceCoordinate[0]);
452  pOsition[1] = (y-sourceCoordinate[1])*(y-sourceCoordinate[1]);
453  pOsition[2] = (z-sourceCoordinate[2])*(z-sourceCoordinate[2]);
454  double absPoint = sqrt(pOsition[0]+pOsition[1]+pOsition[2]);
455  complex< double > attenuationR(0,0);
456  if(isRayleigh) {
457 
458  /* speed of SAW over speed of longitudinal wave */
459  complex< double > vs_over_v1 = (signalLength*fRequency) / sPeed;
460  attenuationR = -sqrt(1.0 - pow(vs_over_v1,2.0))*abs(y);
461 
462  /* attenuationR of Rayleigh inside fluid */
463  // double oMega = 2*M_PI*fRequency;
464  /* double totalViscosity = (4/3) * 0.001139 + 0.0034; */
465  /* sec(Rayleigh angle) = sec(22 degree) = 1.0790 */
466  // attenuationR = (0.0049*pow(oMega,2.0)) / (dEnsity*pow(sPeed,3.0));
467  // attenuationR *= - 1.0790*abs(y);
468 
469  }
470 
471  if(aTtenuation != 0) {
472  double xOriginal = - std::sqrt(std::pow(rAdius,2.0) - std::pow(z,2.0));
473  aTtenuation *= (x - xOriginal);
474  }
475 
476  for(int f = 0;f < sizeTime && sizeFreq;f++) {
477  double speed = signalLength/signalDuration;
478  complex< double > wave_number = (2*M_PI*f/signalLength) + i*complexWaveNumber;
479  double delta_t = signalDuration/sizeTime;
480  double distance = speed*delta_t*timeStep;
481  double phase= 2*M_PI*f*(distance/signalLength);
482  if(isRadiation) {
483  result += (-complexOut[f].r-i*complexOut[f].i)*(exp(-i*wave_number*absPoint+i*phase+wave_number*attenuationR - aTtenuation)) / (4*M_PI*absPoint);
484  } else{
485  result += (complexOut[f].r+i*complexOut[f].i)*(exp(-i*wave_number*absPoint+i*phase+wave_number*attenuationR - aTtenuation)) / (4*M_PI*absPoint);
486  }
487  }
488 
489  result /= sizeTime;
490  result *= transmissionCoefficient;
491 
492  rEsult.resize(2);
493  rEsult[REAL].resize(1);
494  (rEsult[REAL])[0] = std::real(result);
495  rEsult[IMAG].resize(1);
496  (rEsult[IMAG])[0] = std::imag(result);
497 
498  return rEsult;
499 
500  }
501 
502 };
503 
504 /**
505 
506  Transverse Wave of SAW, i.t. is inverse Fourier transform evaluated at arbitrary
507  spatial points.
508 
509  \f[
510  \left. \left\{ \frac{1}{n} (A_{0} e^{ik \mathbf{d} \cdot \mathbf{x} + i \phi})
511  \right\}
512  \f]
513  where \f$\phi\f$ is
514  \f[
515 
516 
517  \f]
518 
519  */
520 // struct TransverseWaveDFT: public GenericAnalyticalSolution {
521 //
522 // vector<VectorDouble > rEsult;
523 // double signalLength;
524 // double signalDuration;
525 // double aTtenuation;
526 // double complexWaveNumber;
527 // double fRequency;
528 // double sPeed;
529 // VectorDouble dIrection;
530 // boost::shared_array<kiss_fft_cpx> complexOut;
531 // int sizeFreq;
532 // int sizeTime;
533 // int timeStep;
534 // VectorDouble cOordinate;
535 // bool isRadiation;
536 // bool isRayleigh;
537 // vector<complex<double> > vecAl; ///< this is to calculate constant values of series only once
538 //
539 // double a;
540 //
541 // TransverseWaveDFT() {}
542 // TransverseWaveDFT(
543 // double signal_length,
544 // double signal_duration,
545 // VectorDouble &d,
546 // double loss_tangent,
547 // double leaky_saw,
548 // double frequency,
549 // double velocity,
550 // boost::shared_array<kiss_fft_cpx> complex_out,
551 // int size_frequency,
552 // int size_time,
553 // int time_step,
554 // bool is_radiation,
555 // bool is_rayleigh,
556 // double sphere_radius = 0.5
557 // ):
558 // signalLength(signal_length),
559 // signalDuration(signal_duration),
560 // dIrection(d),
561 // aTtenuation(loss_tangent),
562 // complexWaveNumber(leaky_saw),
563 // fRequency(frequency),
564 // sPeed(velocity),
565 // complexOut(complex_out),
566 // sizeFreq(size_frequency),
567 // sizeTime(size_time),
568 // timeStep(time_step),
569 // isRadiation(is_radiation),
570 // isRayleigh(is_rayleigh),
571 // a(sphere_radius)
572 // {}
573 //
574 // ~TransverseWaveDFT() {}
575 //
576 // virtual vector<VectorDouble >& operator()(double x, double y, double z) {
577 //
578 // const complex< double > i( 0.0, 1.0 );
579 // cOordinate.resize(3);
580 // cOordinate[0] = x;
581 // cOordinate[1] = y;
582 // cOordinate[2] = z;
583 // complex< double > attenuation(0,0);
584 // double x2 = x*x,y2 = y*y;
585 // double R = sqrt(x2+y2);
586 // double theta = atan2(y,x)+2*M_PI;
587 //
588 // if(isRayleigh) {
589 // /* speed of SAW over speed of longitudinal wave */
590 // complex< double > vs_over_v1 = (signalLength*fRequency) / sPeed;
591 // attenuation = -sqrt(1.0 - pow(vs_over_v1,2.0))*abs(y);
592 // }
593 //
594 // //FIXME can apply zero padding to signal length and sIze in this section when deal with large number of time steps.
595 // complex< double > result = 0.0;
596 // for(int f = 0;f < sizeTime && sizeFreq;f++) {
597 //
598 // double speed = signalLength/signalDuration;
599 // /* 2*M_PI/signalLength : 2*M_PI*f/signalLength */
600 // double k0 = 2*M_PI*f/signalLength;
601 // complex< double > wave_number = (2*M_PI*f/signalLength)*(1.0+i*aTtenuation) + i*complexWaveNumber;
602 // double const1 = k0 * a;
603 // complex< double > const2 = wave_number * R;
604 // complex< double > const3 = wave_number * a;
605 // double delta_t = signalDuration/sizeTime;
606 // double distance = speed*delta_t*timeStep;
607 // double phase= 2*M_PI*f*(distance/signalLength);
608 // complex< double > Bm;
609 //
610 // /* */
611 // const double tol = 1.0e-10;
612 // complex< double > results = 0.0;
613 // complex< double > prev_results;
614 // double error = 100.0;
615 // unsigned int n = 1; //initialized the infinite series loop
616 //
617 // double Jn_der_zero = ( - cyl_bessel_j( 1, const1 ));
618 // complex< double > Jn_der_zero_a = ( - cyl_bessel_j( 1, const3 ));
619 // double Jn_zero = cyl_bessel_j( 0, const1 );
620 // complex< double > Jn_zero_a = cyl_bessel_j( 0, const3 );
621 // complex< double > Jn_zero_r = cyl_bessel_j( 0, const2 );
622 // complex< double > Hn_der_zero = ( - cyl_hankel_1( 1, const1 ));
623 // complex< double >Hn_zero = cyl_hankel_1( 0, const1 ); //S Hankel first kind function
624 // complex< double > Bm0 = ( (complexOut[f].r+i*complexOut[f].i) * k0 * (Hn_zero*Jn_der_zero - Hn_der_zero*Jn_zero) ) /
625 // ( wave_number*Hn_zero*Jn_der_zero_a - k0*Hn_der_zero*Jn_zero_a );
626 // //n=0;
627 // results += Bm0*Jn_zero_r;
628 //
629 // while( error > tol ) //finding the acoustic potential in one single point.
630 // {
631 // if(vecAl.size()>n) {
632 // Bm = vecAl[n-1];
633 // } else {
634 // // cylindrical Bessel function
635 // double Jn_ka = cyl_bessel_j( n, const1 );
636 // complex< double > Jn_wave_number_a = cyl_bessel_j( n, const3 );
637 // double Jn_der_ka = n / const1 * cyl_bessel_j( n, const1 ) - cyl_bessel_j( n + 1, const1 );
638 // complex< double > Jn_der_waver_number_a = double(n) / const3 * cyl_bessel_j( n, const3 ) - cyl_bessel_j( n + 1, const3 );
639 // complex< double > Hn_ka = cyl_hankel_1( n, const1 ); //S Hankel first kind function
640 // complex< double > Hn_der_ka = n / const1 * cyl_hankel_1( n, const1 ) - cyl_hankel_1( n + 1, const1 );
641 //
642 // //Constant term
643 // Bm = ( (complexOut[f].r+i*complexOut[f].i) * 2.0*pow(i,n) * k0 * (Hn_ka*Jn_der_ka - Hn_der_ka*Jn_ka) ) /
644 // ( wave_number*Hn_ka*Jn_der_waver_number_a - k0*Hn_der_ka*Jn_wave_number_a );;
645 // vecAl.push_back(Bm);
646 // }
647 //
648 // complex< double > Jn_wave_number_r = cyl_bessel_j( n, const2 );
649 // prev_results = results;
650 //
651 // results += Bm * Jn_wave_number_r * cos(n*theta) * (exp(i*phase+wave_number));
652 // error = abs( abs( results ) - abs( prev_results ) );
653 // ++n;
654 //
655 // }
656 //
657 //
658 //
659 // // if(isRadiation) {
660 // // result += (-complexOut[f].r-i*complexOut[f].i)*exp(i*wave_number*inner_prod(dIrection,cOordinate)+i*phase+wave_number*attenuation);
661 // // } else{
662 // // result += (complexOut[f].r+i*complexOut[f].i)*exp(i*wave_number*inner_prod(dIrection,cOordinate)+i*phase+wave_number*attenuation);
663 // // }
664 //
665 // }
666 //
667 // result /= sizeTime;
668 //
669 //
670 // rEsult.resize(2);
671 // rEsult[REAL].resize(1);
672 // (rEsult[REAL])[0] = std::real(result);
673 // rEsult[IMAG].resize(1);
674 // (rEsult[IMAG])[0] = std::imag(result);
675 //
676 // return rEsult;
677 //
678 // }
679 //
680 // };
681 
682 
683 #endif // KISS_FFT_H
684 
685 /** Calculate the analytical solution of impinging wave on sphere
686  \ingroup mofem_helmholtz_elem
687 
688 
689  equation from:
690  Ihlenburg,Finite element analysis of acoustic scattering Springer Science & Business Media.
691 
692  \f[
693  p_\textrm{scattered} = \sum_0^N A_l h_l(kr)P_l(\cos(\phi))
694  \f]
695 
696  where \f$h_l\f$ is the Hankel function of the first kind, \f$\phi\f$ is polar
697  angle and \f$A_l\f$ is a constant. Constant is should be calculated such that
698  it satisfies both the Helmholtz wave equation and the Sommerfeld radiation
699  condition.
700 
701  \f[
702  A_l = -(2l+1)i^l \frac{j_{l}'(ka)}{h_{l}'(ka)}
703  \f]
704  where a is scatter sphere radius and \f$j_l\f$ Spherical Bessel function.
705 
706  */
708 
709  vector<complex<double> > vecAl; ///< this is to calculate constant values of series only once
710  vector<VectorDouble > rEsult;
711  double wAvenumber;
712  double sphereRadius;
713 
714 
715  HardSphereScatterWave(double wavenumber,double sphere_radius = 0.5):
716 
717  wAvenumber(wavenumber),sphereRadius(sphere_radius) {}
719 
720  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
721 
722  const double tol = 1.0e-10;
723 
724  double x2 = x*x;
725  double y2 = y*y;
726  double z2 = z*z;
727 
728  double R = sqrt(x2+y2+z2);
729  // double cos_theta = z/R; //Incident wave in Z direction, X =>sin_theta*sin_phi, Y =>sin_theta*cos_phi
730  double cos_theta = x/R; //Incident wave in X direction.
731 
732 
733  const double k = wAvenumber; //Wave number
734  const double a = sphereRadius; //radius of the sphere,wait to modify by user
735 
736  const complex< double > i( 0.0, 1.0 );
737  complex< double > Al;
738 
739  complex< double > result = 0.0;
740  complex< double > prev_result;
741 
742  double error = 100.0;
743  unsigned int n = 0; //initialized the infinite series loop
744 
745  while( error > tol ) //finding the acoustic potential in one single point.
746  {
747 
748  if(vecAl.size()>n) {
749  Al = vecAl[n];
750  } else {
751  // spherical Bessel function
752  double const1 = k*a;
753  double jn_der = n / const1 * sph_bessel( n, const1 ) - sph_bessel( n + 1, const1 );
754 
755  // spherical Hankel function
756  complex< double > hn_der = n / const1 * sph_hankel_1( n, const1 ) - sph_hankel_1( n + 1, const1 );
757  //Constant term
758  Al = -(2.0*n+1)*pow(i,n)*jn_der/hn_der;
759  vecAl.push_back(Al);
760  }
761 
762  prev_result = result;
763 
764  // Legendre function
765  double Pn = legendre_p(n,cos_theta);
766  result += Al*sph_hankel_1(n,k*R)*Pn;
767 
768  error = abs( abs( result ) - abs( prev_result ) );
769 
770  ++n;
771 
772  }
773 
774  rEsult.resize(2);
775  rEsult[REAL].resize(1);
776  (rEsult[REAL])[0] = std::real(result);
777  rEsult[IMAG].resize(1);
778  (rEsult[IMAG])[0] = std::imag(result);
779 
780  return rEsult;
781 
782  }
783 
784 };
785 
786 /** Calculate the analytical solution of impinging wave on sphere
787  \ingroup mofem_helmholtz_elem
788 
789 
790  Equations from:
791  <http://ansol.us/Products/Coustyx/Validation/MultiDomain/Scattering/PlaneWave/SoftSphere/Downloads/dataset_description.pdf>
792 
793  \f[
794  p_\textrm{scattered} = \sum_0^N A_l h_l(kr)P_l(\cos(\phi))
795  \f]
796 
797  where \f$h_l\f$ is the Hankel function of the first kind, \f$\phi\f$ is polar
798  angle and \f$A_l\f$ is a constant. Constant is should be caculated such that
799  it satisfies both the Helmholtz wave equation and the Sommerfeld radiation
800  condition.
801 
802  \f[
803  A_l = -(2l+1)i^l \frac{j_l(ka)}{h_l(ka)}
804  \f]
805 
806  where a is scatter sphere radius and \f$j_l\f$ Spherical Bessel function.
807 
808 
809  */
811 
812  vector<complex<double> > vecAl; ///< this is to calculate constant values of series only once
813  vector<VectorDouble > rEsult;
814  double wAvenumber;
815  double sphereRadius;
816 
817 
818  SoftSphereScatterWave(double wavenumber,double sphere_radius = 0.5):
819 
820  wAvenumber(wavenumber),sphereRadius(sphere_radius) {}
822 
823  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
824 
825  const double tol = 1.0e-10;
826 
827  double x2 = x*x;
828  double y2 = y*y;
829  double z2 = z*z;
830 
831  double R = sqrt(x2+y2+z2);
832  // double cos_theta = z/R; //incident wave in Z direction
833  double cos_theta = x/R; //incident wave in X direction
834 
835  const double k = wAvenumber; //Wave number
836  const double a = sphereRadius; //radius of the sphere,wait to modify by user
837 
838  const complex< double > i( 0.0, 1.0 );
839  complex< double > Al;
840 
841  complex< double > result = 0.0;
842  complex< double > prev_result;
843 
844  double error = 100.0;
845  unsigned int n = 0; //initialized the infinite series loop
846 
847  while( error > tol ) //finding the acoustic potential in one single point.
848  {
849 
850  if(vecAl.size()>n) {
851  Al = vecAl[n];
852  } else {
853  // spherical Bessel function
854  double jn = sph_bessel(n,k*a);
855  // spherical Hankel function
856  complex<double> hn = sph_hankel_1(n,k*a);
857  //Constant term
858  Al = -(2.0*n+1)*pow(i,n)*jn/hn;
859  vecAl.push_back(Al);
860  }
861 
862  prev_result = result;
863 
864  // Legendre function
865  double Pn = legendre_p(n,cos_theta);
866  result += Al*sph_hankel_1(n,k*R)*Pn;
867 
868  error = abs( abs( result ) - abs( prev_result ) );
869 
870  ++n;
871 
872  }
873 
874  rEsult.resize(2);
875  rEsult[REAL].resize(1);
876  (rEsult[REAL])[0] = std::real(result);
877  rEsult[IMAG].resize(1);
878  (rEsult[IMAG])[0] = std::imag(result);
879 
880  return rEsult;
881 
882  }
883 
884 };
885 
886 /** \brief Calculate the analytical solution of plane wave guide propagating in direction theta
887  \ingroup mofem_helmholtz_elem
888 
889  \f[
890  p_\textrm{scattered} = exp^{ik\mathbf{x}\Theta}
891  \f]
892 
893  where:
894 
895  \f[
896  \mathbf{x} = [x,y]
897  \f]
898 
899  \f[
900  \Theta = k[\cos(\theta),\sin(\theta)]
901  \f]
902 
903  \theta is the wave propagating direction from range [0,2\pi]
904 
905 
906  Paper: Gillman, A., Djellouli, R., & Amara, M. (2007).
907  A mixed hybrid formulation based on oscillated finite element polynomials for solving Helmholtz problems.
908  Journal of computational and applied mathematics
909 
910 */
912 
913  vector<VectorDouble > rEsult;
914  double wAvenumber;
915  double tHeta;
916 
917  PlaneWave(double wavenumber,double theta):
918  wAvenumber(wavenumber),tHeta(theta) {}
919  virtual ~PlaneWave() {}
920 
921  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
922 
923  const double k = wAvenumber; //Wave number
924 
925  const complex< double > i( 0.0, 1.0 );
926  complex< double > result = 0.0;
927 
928 
929  //const complex< double > inc_field = exp( i * k * R * cos( theta ) );
930  //const complex< double > total_field = inc_field + result;
931 
932  result = exp(i*(k*cos(tHeta)*x+k*sin(tHeta)*y));
933 
934  rEsult.resize(2);
935  rEsult[REAL].resize(1);
936  (rEsult[REAL])[0] = std::real(result);
937  rEsult[IMAG].resize(1);
938  (rEsult[IMAG])[0] = std::imag(result);
939 
940  return rEsult;
941 
942  }
943 
944 };
945 
946 /** \brief Calculate the analytical solution of impinging wave on cylinder
947  \ingroup mofem_helmholtz_elem
948 
949  \f[
950  p_\textrm{scattered} = \sum_0^N A_l H_l(kr)\cos(l\theta)
951  \f]
952 
953  where \f$H_l\f$ is the cylindrical Hankel function of the first kind, \f$\theta\f$
954  is the polar angle in polar coordinate and \f$A_l\f$ is a constant. Constant is
955  should be caculated such that it satisfies both the Helmholtz wave equation and the
956  Sommerfeld radiation condition.
957 
958  \f[
959  A_l = -\epsilon_{l} i^l \frac{J_{l}'(ka)}{H_{l}'(ka)}
960  \f]
961 
962  where a is scatter sphere radius and \f$J_l\f$ Cylindrical Bessel function.
963 
964  \f[
965  \epsilon_{l} = 1 \textrm{when}l=0
966  \f]
967 
968  \f[
969  \epsilon_{l} = 2 \textrm{when}l \neq 0
970  \f]
971 
972  Paper:
973  Kechroud, R., Soulaimani, A., & Antoine, X. (2009).
974  A performance study of plane wave finite element methods with a Padé-type artificial boundary condition in acoustic scattering.
975  Advances in Engineering Software, 40(8), 738-750.
976 
977 */
979 
980  vector<complex<double> > vecAl; ///< this is to calculate constant values of series only once
981  vector<VectorDouble > rEsult;
982  double wAvenumber;
983  //double shereRadius;
984  double a;
985 
986  HardCylinderScatterWave(double wavenumber,double sphere_radius = 0.5): wAvenumber(wavenumber),a(sphere_radius) {}
988 
989  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
990 
991  const double tol = 1.0e-10;
992  double x2 = x*x,y2 = y*y;
993  double R = sqrt(x2+y2);
994  double theta = atan2(y,x)+2*M_PI;
995 
996  const double k = wAvenumber; //Wave number
997  const double const1 = k * a;
998  double const2 = k * R;
999 
1000  const complex< double > i( 0.0, 1.0 );
1001  complex< double > Al;
1002  // magnitude of incident wave
1003  //const double phi_incident_mag = 1.0;
1004 
1005  complex< double > result = 0.0;
1006  complex< double > prev_result;
1007 
1008  double error = 100.0;
1009  unsigned int n = 1; //initialized the infinite series loop
1010 
1011  double Jn_der_zero = ( - cyl_bessel_j( 1, const1 ));
1012  complex< double > Hn_der_zero = ( - cyl_hankel_1( 1, const1 ));
1013  complex< double >Hn_zero = cyl_hankel_1( 0, const2 ); //S Hankel first kind function
1014 
1015  //n=0;
1016  result -= (Jn_der_zero * Hn_zero)/Hn_der_zero;
1017 
1018  while( error > tol ) //finding the acoustic potential in one single point.
1019  {
1020  if(vecAl.size()>n) {
1021  Al = vecAl[n-1];
1022  } else {
1023  // cylindrical Bessel function
1024  double Jn_der_ka = n / const1 * cyl_bessel_j( n, const1 ) - cyl_bessel_j( n + 1, const1 );
1025  // cylindrical Hankel function
1026  complex<double> Hn_der_ka = n / const1 * cyl_hankel_1( n, const1 ) - cyl_hankel_1( n + 1, const1 );
1027  //Constant term
1028  Al = -2.0*pow(i,n)*Jn_der_ka/Hn_der_ka;
1029  vecAl.push_back(Al);
1030  }
1031 
1032  prev_result = result;
1033  complex< double >Hn_kr = cyl_hankel_1( n, const2 ); //S Hankel first kind function
1034 
1035  result += Al * Hn_kr * cos(n*theta);
1036  error = abs( abs( result ) - abs( prev_result ) );
1037  ++n;
1038 
1039  }
1040 
1041  //result *= phi_incident_mag;
1042 
1043  //const complex< double > inc_field = exp( i * k * R * cos( theta ) ); //
1044  //const complex< double > total_field = inc_field + result;
1045  //ofs << theta << "\t" << abs( result ) << "\t" << abs( inc_field ) << "\t" << abs( total_field ) << "\t" << R << endl; //write the file
1046 
1047  rEsult.resize(2);
1048  rEsult[REAL].resize(1);
1049  (rEsult[REAL])[0] = std::real(result);
1050  rEsult[IMAG].resize(1);
1051  (rEsult[IMAG])[0] = std::imag(result);
1052 
1053  return rEsult;
1054 
1055  }
1056 
1057 };
1058 
1059 
1060 /** \brief Calculate the analytical solution of impinging wave on cylinder
1061  \ingroup mofem_helmholtz_elem
1062 
1063  \f[
1064  p_\textrm{scattered} = \sum_0^N A_l H_l(kr)\cos(l\theta)
1065  \f]
1066 
1067  where \f$H_l\f$ is the cylindrical Hankel function of the first kind, \f$\theta\f$
1068  is the polar angle in polar coordinate and \f$A_l\f$ is a constant. Constant is
1069  should be caculated such that it satisfies both the Helmholtz wave equation and the
1070  Sommerfeld radiation condition.
1071 
1072  \f[
1073  A_l = -\epsilon_{l} i^l \frac{J_{l}(ka)}{H_{l}(ka)}
1074  \f]
1075 
1076  where a is scatter sphere radius and \f$J_l\f$ Cylindrical Bessel function.
1077 
1078  \f[
1079  \epsilon_{l} = 1 \textrm{when}l=0
1080  \f]
1081 
1082  \f[
1083  \epsilon_{l} = 2 \textrm{when}l \neq 0
1084  \f]
1085 
1086  Paper:
1087  Kechroud, R., Soulaimani, A., & Antoine, X. (2009).
1088  A performance study of plane wave finite element methods with a Padé-type artificial boundary condition in acoustic scattering.
1089  Advances in Engineering Software, 40(8), 738-750.
1090 
1091 
1092 */
1094 
1095  vector<complex<double> > vecAl; ///< this is to calculate constant values of series only once
1096  vector<VectorDouble > rEsult;
1097  double wAvenumber;
1098  //double shereRadius;
1099  double a;
1100 
1101  SoftCylinderScatterWave(double wavenumber,double sphere_radius = 0.5): wAvenumber(wavenumber),a(sphere_radius) {}
1103 
1104  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
1105 
1106  const double tol = 1.0e-10;
1107  double x2 = x*x,y2 = y*y;
1108  double R = sqrt(x2+y2);
1109  double theta = atan2(y,x)+2*M_PI;
1110 
1111  const double k = wAvenumber; //Wave number
1112  const double const1 = k * a;
1113  double const2 = k * R;
1114 
1115  const complex< double > i( 0.0, 1.0 );
1116  complex< double > Al;
1117  // magnitude of incident wave
1118  //const double phi_incident_mag = 1.0;
1119 
1120  complex< double > result = 0.0;
1121  complex< double > prev_result;
1122 
1123  double error = 100.0;
1124  unsigned int n = 1; //initialized the infinite series loop
1125 
1126  double Jn_zero = cyl_bessel_j( 0, const1 );
1127  complex< double > Hn_zero_kr = cyl_hankel_1( 0, const2 ); //S Hankel first kind function
1128  complex< double > Hn_zero_ka = cyl_hankel_1( 0, const1 ); //S Hankel first kind function
1129  //n=0;
1130  result -= (Jn_zero * Hn_zero_kr)/Hn_zero_ka;
1131 
1132  while( error > tol ) //finding the acoustic potential in one single point.
1133  {
1134 
1135 
1136  if(vecAl.size()>n) {
1137  Al = vecAl[n-1];
1138  } else {
1139  // cylindrical Bessel function
1140  double Jn_ka = cyl_bessel_j( n, const1 );
1141  // cylindrical Hankel function
1142  complex<double> Hn_ka = cyl_hankel_1( n, const1 );
1143  //Constant term
1144  Al = -2.0*pow(i,n)*Jn_ka/Hn_ka;
1145  vecAl.push_back(Al);
1146  }
1147 
1148 
1149  prev_result = result;
1150 
1151  complex< double >Hn_kr = cyl_hankel_1( n, const2 ); //S Hankel first kind function
1152 
1153  result += Al * Hn_kr * cos(n*theta);
1154  error = abs( abs( result ) - abs( prev_result ) );
1155  ++n;
1156  }
1157 
1158  //result *= phi_incident_mag;
1159 
1160  //const complex< double > inc_field = exp( i * k * R * cos( theta ) );
1161  //const complex< double > total_field = inc_field + result;
1162  //ofs << theta << "\t" << abs( result ) << "\t" << abs( inc_field ) << "\t" << abs( total_field ) << "\t" << R << endl; //write the file
1163 
1164  rEsult.resize(2);
1165  rEsult[REAL].resize(1);
1166  (rEsult[REAL])[0] = std::real(result);
1167  rEsult[IMAG].resize(1);
1168  (rEsult[IMAG])[0] = std::imag(result);
1169 
1170  return rEsult;
1171 
1172  }
1173 
1174 };
1175 
1176 
1177 
1178 /** Calculate the analytical solution of impinging wave on L shaped domain with singularity
1179  \ingroup mofem_helmholtz_elem
1180 
1181 
1182  Equations from:
1183  <The p and h-p Versions of the Finite Element Method, Basic Principles and Properties Read More: http://epubs.siam.org/doi/abs/10.1137/1036141>
1184 
1185  \f[
1186  p_\textrm{scattered} = j_{\frac{2}{3}}(kr)\sin(\frac{2}{3}\phi)
1187  \f]
1188  for 0 < \f$\phi\f$ < \f$\frac{2}{3}\pi\f$
1189  where r is scatter sphere radius and \f$j_l\f$ Spherical Bessel function,
1190  \f$\phi\f$ is polar angle
1191 
1192  Constant should be caculated such that
1193  it satisfies both the Helmholtz wave equation and the Robin boundary condition.
1194 
1195 
1196 
1197 
1198  */
1200 
1201 
1202  vector<VectorDouble > rEsult;
1203  double wAvenumber;
1204 
1205 
1206  SingularScatterWave(double wavenumber):
1207 
1208  wAvenumber(wavenumber) {}
1210 
1211  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
1212 
1213  const double tol = 1.0e-10;
1214 
1215  double x2 = x*x,y2 = y*y;
1216 
1217  double R = sqrt(x2+y2);
1218  double theta = atan2(y,x)+2*M_PI;
1219  // double cos_theta = z/R; //incident wave in Z direction
1220  // double cos_theta = x/R; //incident wave in X direction
1221  // double sin_theta = y/R;
1222 
1223  const double k = wAvenumber; //Wave number
1224 
1225  double const1 = k * R;
1226  double const2 = sin(2*theta/3);
1227  // double const2 = sin_theta;
1228  double Jn_two_third = cyl_bessel_j( 0.666666666666667, const1 );
1229 
1230  complex< double > result = 0.0;
1231  result = Jn_two_third * const2;
1232 
1233  // result = exp(theta);
1234 
1235  rEsult.resize(2);
1236  rEsult[REAL].resize(1);
1237  (rEsult[REAL])[0] = std::real(result);
1238  rEsult[IMAG].resize(1);
1239  (rEsult[IMAG])[0] = std::imag(result);
1240 
1241  return rEsult;
1242 
1243  }
1244 
1245 };
1246 
1247 /** Calculate the incident wave of impinging wave on L shaped domain with singularity
1248  \ingroup mofem_helmholtz_elem
1249 
1250 
1251  Equations from:
1252  <The p and h-p Versions of the Finite Element Method, Basic Principles and Properties Read More: http://epubs.siam.org/doi/abs/10.1137/1036141>
1253 
1254  \f[
1255  p_\textrm{scattered} = j_{\frac{2}{3}}(kr)\sin(\frac{2}{3}\phi)
1256  \f]
1257  for 0 < \f$\phi\f$ < \f$\frac{2}{3}\pi\f$
1258  where r is scatter sphere radius and \f$j_l\f$ Spherical Bessel function,
1259  \f$\phi\f$ is polar angle
1260 
1261  Constant should be caculated such that
1262  it satisfies both the Helmholtz wave equation and the Robin boundary condition.
1263 
1264 
1265 
1266 
1267  */
1269 
1270 
1271  vector<VectorDouble > rEsult;
1272  double wAvenumber;
1273 
1274 
1275  SingularIncidentWave(double wavenumber):
1276 
1277  wAvenumber(wavenumber) {}
1279 
1280  virtual vector<VectorDouble >& operator()(double x, double y, double z) {
1281 
1282  const double tol = 1.0e-10;
1283 
1284  double x2 = x*x,y2 = y*y;
1285 
1286  double R = sqrt(x2+y2);
1287  double theta = atan2(y,x)+2*M_PI;
1288  // double cos_theta = z/R; //incident wave in Z direction
1289  // double cos_theta = x/R; //incident wave in X direction
1290  // double sin_theta = y/R;
1291 
1292  const double k = wAvenumber; //Wave number
1293 
1294  // double const1 = k * R;
1295  // double const2 = sin(2*theta/3);
1296  // // double const2 = sin_theta;
1297  //
1298  // double Jn_two_third = cyl_bessel_j( 0.66667, const1 );
1299  // double D_Jn_two_third = (0.66667/const1) * Jn_two_third - cyl_bessel_j( 1.66667, const1 )
1300  //
1301  complex< double > result = 0.0;
1302  // result = D_Jn_two_third * const2;
1303 
1304  // result = exp(theta);
1305 
1306  // double const1 = pow(R,-0.5);
1307  // double const2 = sinh(R) * cos(theta/2)
1308  // result =
1309 
1310 
1311  rEsult.resize(2);
1312  rEsult[REAL].resize(1);
1313  (rEsult[REAL])[0] = std::real(result);
1314  rEsult[IMAG].resize(1);
1315  (rEsult[IMAG])[0] = std::imag(result);
1316 
1317  return rEsult;
1318 
1319  }
1320 
1321 };
1322 
1323 
1324 
1325 
1326 // **** Surface boundary conditions ( used to enforce BC on surface ) ****
1327 
1328 // ** Dirichlet BC **
1329 
1330 // **** Function to solve best approximation problem ****
1331 
1332 template <typename FUNVAL>
1334  MoFEM::Interface& m_field,const string &problem_name,const string &fe_name,const string &re_field,
1335  Mat A,vector<Vec> vec_F,FUNVAL &fun_evaluator
1336 ) {
1337  PetscFunctionBegin;
1338 
1339  PetscErrorCode ierr;
1340  FieldApproximationH1 field_approximation(m_field);
1341  // This increase rule for numerical integration. In case of 10 node
1342  // elements jacobian is varying linearly across element, that way to element
1343  // rule is add 1.
1344  // field_approximation.addToRule = 1;
1345  // field_approximation.feVolume = HelmholtzElement::MyVolumeFE(m_field,true,1);
1346 
1347  ierr = field_approximation.loopMatrixAndVectorVolume(
1348  problem_name,fe_name,re_field,A,vec_F,fun_evaluator
1349  ); CHKERRQ(ierr);
1350 
1351 
1352  if(A) {
1353  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1354  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1355  }
1356 
1357  PetscFunctionReturn(0);
1358 }
1359 
1360 PetscErrorCode save_data_on_mesh(
1361  MoFEM::Interface& m_field,const string &problem_name,const string& re_field,const string &im_field,
1362  Vec D,int vt,InsertMode mode,PetscBool is_partitioned
1363 ) {
1364  PetscFunctionBegin;
1365 
1366  PetscErrorCode ierr;
1367  // save data on mesh
1369  /* set data to field from solution vec */
1370  if(is_partitioned) {
1371  ierr = m_field.set_local_ghost_vector(problem_name,COL,D,mode,SCATTER_REVERSE); CHKERRQ(ierr);
1372  } else {
1373  ierr = m_field.set_global_ghost_vector(problem_name,COL,D,mode,SCATTER_REVERSE); CHKERRQ(ierr);
1374  }
1375 
1376  } else {
1377 
1378  if(is_partitioned) {
1379  ierr = m_field.set_other_local_ghost_vector(problem_name,re_field,im_field,COL,D,mode,SCATTER_REVERSE); CHKERRQ(ierr);
1380  } else {
1381  ierr = m_field.set_other_global_ghost_vector(problem_name,re_field,im_field,COL,D,mode,SCATTER_REVERSE); CHKERRQ(ierr);
1382  }
1383 
1384  }
1385 
1386  PetscFunctionReturn(0);
1387 }
1388 
1389 template <typename FUNEVAL>
1390 PetscErrorCode solve_problem(MoFEM::Interface& m_field,
1391  const string& problem_name,const string& fe_name,
1392  const string& re_field,const string &im_field,
1393  InsertMode mode,
1394  FUNEVAL &fun_evaluator,PetscBool is_partitioned) {
1395  PetscFunctionBegin;
1396 
1397  PetscErrorCode ierr;
1398 
1399  Mat A;
1400  ierr = m_field.MatCreateMPIAIJWithArrays(problem_name,&A); CHKERRQ(ierr);
1401  ierr = MatZeroEntries(A); CHKERRQ(ierr);
1402  Vec D;
1403 
1404  vector<Vec> vec_F;
1405  vec_F.resize(2);
1406 
1407  ierr = m_field.VecCreateGhost(problem_name,ROW,&vec_F[0]); CHKERRQ(ierr); /* real */
1408  ierr = m_field.VecCreateGhost(problem_name,ROW,&vec_F[1]); CHKERRQ(ierr); /* imag */
1409  ierr = m_field.VecCreateGhost(problem_name,COL,&D); CHKERRQ(ierr);
1410 
1411  for(int ss = 0;ss<2;ss++) {
1412  ierr = VecZeroEntries(vec_F[ss]); CHKERRQ(ierr);
1413  ierr = VecGhostUpdateBegin(vec_F[ss],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1414  ierr = VecGhostUpdateEnd(vec_F[ss],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1415  }
1416 
1417  ierr = calculate_matrix_and_vector(m_field,problem_name,fe_name,re_field,A,vec_F,fun_evaluator); CHKERRQ(ierr);
1418 
1419  KSP solver;
1420  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
1421  ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
1422  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
1423  ierr = KSPSetUp(solver); CHKERRQ(ierr);
1424 
1425  for(int ss = 0;ss<GenericAnalyticalSolution::LAST_VAL_TYPE;ss++) {
1426  // solve problem
1427  ierr = KSPSolve(solver,vec_F[ss],D); CHKERRQ(ierr);
1428  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1429  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1430 
1431 
1432  ierr = save_data_on_mesh(m_field,problem_name,re_field,im_field,D,ss,mode,is_partitioned); CHKERRQ(ierr);
1433 
1434  }
1435 
1436  // clean
1437  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1438  ierr = VecDestroy(&vec_F[GenericAnalyticalSolution::REAL]); CHKERRQ(ierr);
1439  ierr = VecDestroy(&vec_F[GenericAnalyticalSolution::IMAG]); CHKERRQ(ierr);
1440 
1441  ierr = VecDestroy(&D); CHKERRQ(ierr);
1442  ierr = MatDestroy(&A); CHKERRQ(ierr);
1443 
1444  PetscFunctionReturn(0);
1445 }
vector< VectorDouble > rEsult
SoftCylinderScatterWave(double wavenumber, double sphere_radius=0.5)
Deprecated interface functions.
DEPRECATED MoFEMErrorCode set_global_ghost_vector(const Problem *problem_ptr, RowColData rc, Vec V, InsertMode mode, ScatterMode scatter_mode) const
set values of vector from/to mesh database (collective)collective - need tu be run on all processors ...
virtual vector< VectorDouble > & operator()(double x, double y, double z)
Calculate the analytical solution of impinging wave on cylinder.
virtual ~PlaneWave()
vector< VectorDouble > rEsult
const double D
diffusivity
vector< VectorDouble > rEsult
virtual vector< VectorDouble > & operator()(double x, double y, double z)
virtual vector< VectorDouble > & operator()(double x, double y, double z)
SoftSphereScatterWave(double wavenumber, double sphere_radius=0.5)
DEPRECATED MoFEMErrorCode set_local_ghost_vector(const Problem *problem_ptr, RowColData rc, Vec V, InsertMode mode, ScatterMode scatter_mode) const
set values of vector from/to meshdatabase
IncidentWave(double wavenumber, VectorDouble d, bool is_rayleigh=false, double attenuation=0, double radius=0, double leaky_saw=0, double frequency=0, double velocity=0, double transmission_coefficient=1, double r_amplitude=1, double i_amplitude=0, double phase=0)
vector< complex< double > > vecAl
this is to calculate constant values of series only once
vector< complex< double > > vecAl
this is to calculate constant values of series only once
double amplitudeReal
The real amplitude of the incident wave.
Calculate the analytical solution of impinging wave on cylinder.
virtual vector< VectorDouble > & operator()(double x, double y, double z)
static Index< 'n', 3 > n
vector< complex< double > > vecAl
this is to calculate constant values of series only once
Calculate the analytical solution of plane wave guide propagating in direction theta.
Generic structure for analytical function.
DEPRECATED MoFEMErrorCode set_other_global_ghost_vector(const Problem *problem_ptr, const std::string &field_name, const std::string &cpy_field_name, RowColData rc, Vec V, InsertMode mode, ScatterMode scatter_mode, int verb=-1)
Copy vector to field which is not part of the problem (collective)collective - need tu be run on all ...
PetscErrorCode solve_problem(MoFEM::Interface &m_field, const string &problem_name, const string &fe_name, const string &re_field, const string &im_field, InsertMode mode, FUNEVAL &fun_evaluator, PetscBool is_partitioned)
const char * analytical_solution_types[]
virtual vector< VectorDouble > & operator()(double x, double y, double z)
virtual vector< VectorDouble > & operator()(double x, double y, double z)
virtual vector< VectorDouble > & operator()(double x, double y, double z)
PetscErrorCode save_data_on_mesh(MoFEM::Interface &m_field, const string &problem_name, const string &re_field, const string &im_field, Vec D, int vt, InsertMode mode, PetscBool is_partitioned)
MoFEMErrorCode loopMatrixAndVectorVolume(const std::string &problem_name, const std::string &fe_name, const std::string &field_name, Mat A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
assemble matrix and vector
vector< VectorDouble > rEsult
vector< VectorDouble > rEsult
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
VectorDouble cOordinate
vector< VectorDouble > rEsult
double tol
static Index< 'i', 3 > i
DEPRECATED MoFEMErrorCode set_other_local_ghost_vector(const Problem *problem_ptr, const std::string &fiel_name, const std::string &cpy_field_name, RowColData rc, Vec V, InsertMode mode, ScatterMode scatter_mode, int verb=-1)
Copy vector to field which is not part of the problem.
PetscErrorCode calculate_matrix_and_vector(MoFEM::Interface &m_field, const string &problem_name, const string &fe_name, const string &re_field, Mat A, vector< Vec > vec_F, FUNVAL &fun_evaluator)
virtual vector< VectorDouble > & operator()(double x, double y, double z)
CHKERRQ(ierr)
PlaneWave(double wavenumber, double theta)
HardCylinderScatterWave(double wavenumber, double sphere_radius=0.5)
double amplitudeReal
The real amplitude of the incident wave.
IncidentWavePointSource(double wavenumber, VectorDouble d, bool is_rayleigh=false, double attenuation=0, double radius=0, double leaky_saw=0, double frequency=0, double velocity=0, double transmission_coefficient=1, double r_amplitude=1, double i_amplitude=0)
vector< VectorDouble > rEsult
virtual vector< VectorDouble > & operator()(double x, double y, double z)=0
DEPRECATED MoFEMErrorCode VecCreateGhost(const std::string &name, RowColData rc, Vec *V) const
create ghost vector for problem (collective)collective - need to be run on all processors in communic...
VectorDouble dIrection
const double R
SingularIncidentWave(double wavenumber)
AnalyticalSolutionTypes
vector< VectorDouble > rEsult
virtual vector< VectorDouble > & operator()(double x, double y, double z)
static Index< 'k', 3 > k
Finite element for approximating analytical filed on the mesh.
double transmissionCoefficient
vector< complex< double > > vecAl
this is to calculate constant values of series only once
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
DEPRECATED MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)
HardSphereScatterWave(double wavenumber, double sphere_radius=0.5)
vector< VectorDouble > rEsult
SingularScatterWave(double wavenumber)