v0.9.0
HelmholtzElementObsolete.hpp
Go to the documentation of this file.
1 /** \file HelmholtzElement.hpp
2  \ingroup mofem_helmholtz_elem
3 
4  \brief Operators and data structures for wave propagation analyse (Galerkin Element)
5 
6  Implementation of Helmholtz element for wave propagation problem
7 
8  */
9 
10 /*
11  This work is part of PhD thesis by on Micro-fluids: Thomas Felix Xuan Meng
12  */
13 
14  /*
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
22  * The header file should contains as least #include as possible for speed
23  */
24 
25 #ifndef __HELMHOLTZ_ELEMENT_HPP
26 #define __HELMHOLTZ_ELEMENT_HPP
27 
28 #include<moab/Skinner.hpp>
29 using namespace std;
30 using namespace boost::math;
31 
32 
33 
34 namespace MoFEM {
35 
36 /** \brief struture grouping operators and data used for helmholtz problems
37 * \ingroup mofem_helmholtz_elem
38 *
39 * In order to assemble matrices and right hand side vectors, the loops over
40 * elements, enetities over that elememnts and finally loop over intergration
41 * points are executed.
42 *
43 * complex-real transformation are implemented for computational efficenicy
44 * See Ercegovac, Milos, and Jean-Michel Muller. "Solving Systems of Linear
45 * Equations in Complex Domain: Complex E-Method." (2007). for details.
46 *
47 *
48 * Following implementation separte those three cegories of loops and to each
49 * loop attach operator.
50 *
51 */
52 // Three dimensional homogeneous isotropic medium and time harmonic Helmholtz operator
54  /// \brief definition of volume element
57 
58  /** \brief it is used to calculate nb. of Gauss integartion points
59  *
60  * for more details pleas look
61  * Reference:
62  *
63  * Albert Nijenhuis, Herbert Wilf,
64  * Combinatorial Algorithms for Computers and Calculators,
65  * Second Edition,
66  * Academic Press, 1978,
67  * ISBN: 0-12-519260-6,
68  * LC: QA164.N54.
69  *
70  * More details about algorithm
71  * http://people.sc.fsu.edu/~jburkardt/cpp_src/gm_rule/gm_rule.html
72  **/
73  int getRule(int order) { return order; };
74  };
75 
76  MyVolumeFE feRhs; ///< cauclate right hand side for tetrahedral elements
77  MyVolumeFE& getLoopFeRhs() { return feRhs; } ///< get rhs volume element
78  MyVolumeFE feLhs; //< calculate left hand side for tetrahedral elements
79  MyVolumeFE& getLoopFeLhs() { return feLhs; } ///< get lhs volume element
80 
81  /** \brief define surface element
82  *
83  *
84  */
87  int getRule(int order) { return order; };
88  };
89 
90  MyTriFE feRhs_Tri; ///< cauclate right hand side for tetrahedral elements
91  MyTriFE& getLoopFeRhs_Tri() { return feRhs_Tri; } ///< get rhs volume element
92 
94  MyTriFE& getLoopFeFlux() { return feFlux; }
95 
96  MyTriFE feIncidentWave; //incident wave
97  MyTriFE& getLoopfeIncidentWave() { return feIncidentWave; } //incident wave
98 
100  MyTriFE& getLoopFeImpedanceLhs() { return feImpedanceLhs; }
101 
102 
103 
106  MoFEM::Interface &m_field):
107  feRhs(m_field),feRhs_Tri(m_field),feLhs(m_field),feFlux(m_field),feIncidentWave(m_field),feImpedanceLhs(m_field),mField(m_field) {}
108 
109  /** \brief data for calulation Angular Frequency and wave speed elements
110  * \infroup mofem_helmholtz_elem
111  */
112  struct BlockData {
113  double aNgularfreq; // Angular frequency
114  double sPeed; // Wave Speed
115 
116  Range tRis;
117  Range tEts; ///< constatins elements in block set
118  };
119  map<int,BlockData> setOfBlocks; ///< maps block set id with appropiate BlockData
121 
122  /** \brief data for calulation heat flux
123  * \infroup mofem_helmholtz_elem
124  */
125  struct FluxData {
126  HeatfluxCubitBcData dAta; ///< for more details look to BCMultiIndices.hpp to see details of HeatfluxCubitBcData
127  Range tRis; ///< suraface triangles where hate flux is applied
128  };
129  map<int,FluxData> setOfFluxes; ///< maps side set id with appropiate FluxData
130 
131  struct IncidentData {
132  double amplitude; ///< for more details look to BCMultiIndices.hpp to see details of HeatfluxCubitBcData
133  Range tRis; ///< suraface triangles where hate flux is applied
134  };
135  map<int,IncidentData> setOfIncident; ///< maps side set id with appropiate FluxData
136 
137  /** \brief data for Robin Boundary condition
138  * \infroup mofem_helmholtz_elem
139  */
140  struct ImpedanceData {
141  double g; /*the zero velocity term of the Robin BC*/
142  double sIgma; //admittance of the surface, if g and sIgma both zero the surface is called sound hard.
143  Range tRis; ///< those will be on body skin, except thos with contact whith other body where pressure is applied
144  };
145  map<int,ImpedanceData> setOfImpedance; //< maps block set id with appropiate data
146 
147 
148 
149  /** \brief common data used by volume elements
150  * \infroup mofem_helmholtz_elem
151  */
152  struct CommonData {
153 
154  VectorDouble pressureReAtGaussPts; //Helmholtz Pressure real
155  VectorDouble pressureImAtGaussPts; //Helmholtz Pressure imag
158 
159  inline ublas::matrix_row<MatrixDouble > getGradReAtGaussPts(const int gg) {
160  return ublas::matrix_row<MatrixDouble >(gradReAtGaussPts,gg);
161  }
162 
163  inline ublas::matrix_row<MatrixDouble > getGradImAtGaussPts(const int gg) {
164  return ublas::matrix_row<MatrixDouble >(gradImAtGaussPts,gg);
165  }
166 
167  MatrixDouble Ann,Anv,Avv;
168  ublas::vector<MatrixDouble > Aen,Afn,Aev,Afv;
169  ublas::matrix<MatrixDouble > Aee,Aef,Aff;
170 
171  };
173 
174  /// \brief operator to calculate pressure gradient at Gauss points
176 
178  OpGetGradReAtGaussPts(const string field_name,CommonData &common_data):
180  commonData(common_data) {}
181 
182  /** \brief operator calculating pressure gradients
183  *
184  * pressure gradient is calculated multiplying derivatives of shape functions by degrees of freedom.
185  */
186  PetscErrorCode doWork(
187  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
188  PetscFunctionBegin;
189  try {
190 
191  if(data.getIndices().size()==0) PetscFunctionReturn(0);
192  int nb_dofs = data.getFieldData().size();
193  int nb_gauss_pts = data.getN().size1();
194 
195  //initialize
196  commonData.gradReAtGaussPts.resize(nb_gauss_pts,3);
197  if(type == MBVERTEX) {
198  fill(commonData.gradReAtGaussPts.data().begin(),commonData.gradReAtGaussPts.data().end(),0);
199  }
200 
201  for(int gg = 0;gg<nb_gauss_pts;gg++) {
202  ublas::noalias(commonData.getGradReAtGaussPts(gg)) += prod( trans(data.getDiffN(gg,nb_dofs)), data.getFieldData() );
203  }
204 
205  } catch (const std::exception& ex) {
206  ostringstream ss;
207  ss << "throw in method: " << ex.what() << endl;
208  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
209  }
210 
211  PetscFunctionReturn(0);
212  }
213 
214  };
215 
216  /// \brief operator to calculate pressure gradient at Gauss points
218 
220  OpGetGradImAtGaussPts(const string field_name,CommonData &common_data):
222  commonData(common_data) {}
223 
224  /** \brief operator calculating pressure gradients
225  *
226  * pressure gradient is calculated multiplying derivatives of shape functions by degrees of freedom.
227  */
228  PetscErrorCode doWork(
229  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
230  PetscFunctionBegin;
231  try {
232 
233  if(data.getIndices().size()==0) PetscFunctionReturn(0);
234  int nb_dofs = data.getFieldData().size();
235  int nb_gauss_pts = data.getN().size1();
236 
237  //initialize
238  commonData.gradImAtGaussPts.resize(nb_gauss_pts,3);
239  if(type == MBVERTEX) {
240  fill(commonData.gradImAtGaussPts.data().begin(),commonData.gradImAtGaussPts.data().end(),0);
241  }
242 
243  for(int gg = 0;gg<nb_gauss_pts;gg++) {
244  ublas::noalias(commonData.getGradImAtGaussPts(gg)) += prod( trans(data.getDiffN(gg,nb_dofs)), data.getFieldData() );
245  }
246 
247  } catch (const std::exception& ex) {
248  ostringstream ss;
249  ss << "throw in method: " << ex.what() << endl;
250  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
251  }
252 
253  PetscFunctionReturn(0);
254  }
255 
256  };
257 
258  /** \brief opearator to caulate pressure and rate of pressure at Gauss points
259  * \infroup mofem_helmholtz_elem
260  */
261  template<typename OP>
263 
265  OpGetFieldAtGaussPts(const string field_name,VectorDouble &field_at_gauss_pts):
266  OP::UserDataOperator(field_name),
267  fieldAtGaussPts(field_at_gauss_pts) {}
268  /** \brief operator calculating pressure and rate of pressure
269  *
270  * pressure pressure or rate of pressure is calculated multiplyingshape functions by degrees of freedom
271  */
272  PetscErrorCode doWork(
273  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
274  PetscFunctionBegin;
275  try {
276 
277  if(data.getFieldData().size()==0) PetscFunctionReturn(0);
278  int nb_dofs = data.getFieldData().size();
279  int nb_gauss_pts = data.getN().size1();
280 
281  //initialize
282  fieldAtGaussPts.resize(nb_gauss_pts);
283  if(type == MBVERTEX) {
284  //loop over shape functions on entities allways start from
285  //vertices, so if nodal shape functions are processed, vector of
286  //field values is zeroad at initialization
287  fill(fieldAtGaussPts.begin(),fieldAtGaussPts.end(),0);
288  }
289 
290  for(int gg = 0;gg<nb_gauss_pts;gg++) {
291  fieldAtGaussPts[gg] += inner_prod(data.getN(gg,nb_dofs),data.getFieldData());
292 
293  }
294 
295  } catch (const std::exception& ex) {
296  ostringstream ss;
297  ss << "throw in method: " << ex.what() << endl;
298  SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
299  }
300 
301  PetscFunctionReturn(0);
302  }
303 
304  };
305 
306  /** \brief operator to calculate pressure at Gauss pts
307  * \infroup mofem_helmholtz_elem
308  */
309  struct OpGetTetPressureReAtGaussPts: public OpGetFieldAtGaussPts<TetElementForcesAndSourcesCore> {
310  OpGetTetPressureReAtGaussPts(const string field_name,CommonData &common_data):
311  OpGetFieldAtGaussPts<TetElementForcesAndSourcesCore>(field_name,common_data.pressureReAtGaussPts) {}
312  };
313  struct OpGetTetPressureImAtGaussPts: public OpGetFieldAtGaussPts<TetElementForcesAndSourcesCore> {
314  OpGetTetPressureImAtGaussPts(const string field_name,CommonData &common_data):
315  OpGetFieldAtGaussPts<TetElementForcesAndSourcesCore>(field_name,common_data.pressureImAtGaussPts) {}
316  };
317 
318  /** \brief operator to calculate pressure at Gauss pts
319  * \infroup mofem_helmholtz_elem
320  */
321  struct OpGetTriPressureReAtGaussPts: public OpGetFieldAtGaussPts<TriElementForcesAndSurcesCore> {
322  OpGetTriPressureReAtGaussPts(const string field_name,CommonData &common_data):
323  OpGetFieldAtGaussPts<TriElementForcesAndSurcesCore>(field_name,common_data.pressureReAtGaussPts) {}
324  };
325 
326  struct OpGetTriPressureImAtGaussPts: public OpGetFieldAtGaussPts<TriElementForcesAndSurcesCore> {
327  OpGetTriPressureImAtGaussPts(const string field_name,CommonData &common_data):
328  OpGetFieldAtGaussPts<TriElementForcesAndSurcesCore>(field_name,common_data.pressureImAtGaussPts) {}
329  };
330 
331 
333 
334  //get the higher order approximation coordinates.
336 
338  OpHoCoordTri(const string field_name,MatrixDouble &ho_coords):
340  hoCoordsTri(ho_coords) {}
341 
342  /*
343  Cartesian coordinates for gaussian points inside elements
344  X^coordinates = DOF dot* N
345  */
346  PetscErrorCode doWork(
347  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
348  PetscFunctionBegin;
349 
350  try {
351 
352  if(data.getFieldData().size()==0) PetscFunctionReturn(0);
353 
354  hoCoordsTri.resize(data.getN().size1(),3);
355  if(type == MBVERTEX) {
356  hoCoordsTri.clear();
357  }
358 
359  int nb_dofs = data.getFieldData().size();
360 
361  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
362  for(int dd = 0;dd<3;dd++) {
363  hoCoordsTri(gg,dd) += cblas_ddot(nb_dofs/3,&data.getN(gg)[0],1,&data.getFieldData()[dd],3); //calculate x,y,z at each GaussPts
364  }
365  }
366 
367  } catch (const std::exception& ex) {
368  ostringstream ss;
369  ss << "throw in method: " << ex.what() << endl;
370  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
371  }
372 
373  PetscFunctionReturn(0);
374  }
375 
376  };
377 
379 
380  //get the higher order approximation coordinates.
382 
384  OpHoCoordTet(const string field_name,MatrixDouble &ho_coords):
386  hoCoordsTet(ho_coords) {}
387 
388  /*
389  Cartesian coordinates for gaussian points inside elements
390  X^coordinates = DOF dot* N
391  */
392  PetscErrorCode doWork(
393  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
394  PetscFunctionBegin;
395 
396  try {
397 
398  if(data.getFieldData().size()==0) PetscFunctionReturn(0);
399 
400  hoCoordsTet.resize(data.getN().size1(),3);
401  if(type == MBVERTEX) {
402  hoCoordsTet.clear();
403  }
404 
405  int nb_dofs = data.getFieldData().size();
406 
407  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
408  for(int dd = 0;dd<3;dd++) {
409  hoCoordsTet(gg,dd) += cblas_ddot(nb_dofs/3,&data.getN(gg)[0],1,&data.getFieldData()[dd],3); //calculate x,y,z at each GaussPts
410  }
411  }
412 
413  } catch (const std::exception& ex) {
414  ostringstream ss;
415  ss << "throw in method: " << ex.what() << endl;
416  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
417  }
418 
419  PetscFunctionReturn(0);
420  }
421 
422  };
423 
424 
425 
426  /** \biref operator to calculate right hand side of stiffness terms
427  * \infroup mofem_helmholtz_elem
428  */
430 
433  bool useTsF;
434 
435  OpHelmholtzRhs_Re(const string re_field_name,const string im_field_name,BlockData &data,CommonData &common_data):
436  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
437  dAta(data),commonData(common_data),useTsF(true) {}
438 
439  Vec F;
440  OpHelmholtzRhs_Re(const string re_field_name,const string im_field_name,Vec _F,BlockData &data,CommonData &common_data):
441  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
442  dAta(data),commonData(common_data),useTsF(false),F(_F) { }
443 
444  /* tip:all vectors should be declared as class members instead of as local function variables. */
447 
448  /** \brief calculate helmholtz operator apply on lift.
449  *
450  * F = int diffN^T gard_T0 - N^T k T0 dOmega^2
451  *
452  */
453  PetscErrorCode doWork(
454  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
455  PetscFunctionBegin;
456 
457  if(dAta.tEts.find(getMoFEMFEPtr()->get_ent()) == dAta.tEts.end()) {
458  PetscFunctionReturn(0);
459  }
460 
461  try {
462 
463  if(data.getIndices().size()==0) PetscFunctionReturn(0);
464  if(dAta.tEts.find(getMoFEMFEPtr()->get_ent())==dAta.tEts.end()) PetscFunctionReturn(0);
465 
466  PetscErrorCode ierr;
467 
468  /* moFEM tip: use "row_data.getIndices().size()" instead of "row_data.getN().size2()" for nb_row_dofs
469  the former can cause problems with mix approximation orders */
470 
471  int nb_row_dofs = data.getIndices().size();
472  Nf.resize(nb_row_dofs);
473  bzero(&*Nf.data().begin(),data.getIndices().size()*sizeof(FieldData));
474 
475  //wave number K is the propotional to the frequency of incident wave
476  //and represents number of waves per wave length 2Pi - 2Pi/K
477  double wAvenumber = dAta.aNgularfreq/dAta.sPeed;
478  double wAvenUmber = pow(wAvenumber,2.0);
479  dNdT_Re.resize(nb_row_dofs);
480 
481  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
482 
483  double val = getVolume()*getGaussPts()(3,gg);
484 
485  if(getHoGaussPtsDetJac().size()>0) {
486  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
487  }
488 
489  dNdT_Re.clear();
490  /* k^2 N^T q */
491  double const_kT_Re = wAvenUmber*commonData.pressureReAtGaussPts[gg];
492  /* diffN^T diffN q */
493  ublas::noalias(dNdT_Re) = prod(data.getDiffN(gg,nb_row_dofs),commonData.getGradReAtGaussPts(gg));
494  /* diffN^T diffN q - k^2 N^T N q lift of helmholtz operator*/
495 
496  cblas_daxpy(nb_row_dofs, -const_kT_Re, &data.getN(gg,nb_row_dofs)[0], 1, &dNdT_Re[0], 1);
497  /* insufficient, FIX ME*/
498 
499  ublas::noalias(Nf) += val*(dNdT_Re);
500 
501  }
502 
503  if(useTsF) {
504  ierr = VecSetValues(getFEMethod()->ts_F,data.getIndices().size(),
505  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
506  } else {
507  ierr = VecSetValues(F,data.getIndices().size(),
508  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
509 
510  }
511 
512  } catch (const std::exception& ex) {
513  ostringstream ss;
514  ss << "throw in method: " << ex.what() << endl;
515  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
516  }
517 
518  PetscFunctionReturn(0);
519  }
520 
521  };
522 
523  /** \biref operator to calculate right hand side of stiffness terms
524  * \infroup mofem_helmholtz_elem
525  */
527 
530  bool useTsF;
531 
532  OpHelmholtzRhs_Im(const string re_field_name,const string im_field_name,BlockData &data,CommonData &common_data):
533  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
534  dAta(data),commonData(common_data),useTsF(true) {}
535 
536  Vec F;
537  OpHelmholtzRhs_Im(const string re_field_name,const string im_field_name,Vec _F,BlockData &data,CommonData &common_data):
538  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
539  dAta(data),commonData(common_data),useTsF(false),F(_F) { }
540 
541  /* tip:all vectors should be declared as class members instead of as local function variables. */
544 
545  /** \brief calculate helmholtz operator apply on lift.
546  *
547  * F = int diffN^T gard_T0 - N^T k T0 dOmega^2
548  *
549  */
550  PetscErrorCode doWork(
551  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
552  PetscFunctionBegin;
553 
554  if(dAta.tEts.find(getMoFEMFEPtr()->get_ent()) == dAta.tEts.end()) {
555  PetscFunctionReturn(0);
556  }
557 
558  try {
559 
560  if(data.getIndices().size()==0) PetscFunctionReturn(0);
561  if(dAta.tEts.find(getMoFEMFEPtr()->get_ent())==dAta.tEts.end()) PetscFunctionReturn(0);
562 
563  PetscErrorCode ierr;
564 
565  /* moFEM tip: use "row_data.getIndices().size()" instead of "row_data.getN().size2()" for nb_row_dofs
566  the former can cause problems with mix approximation orders */
567 
568  int nb_row_dofs = data.getIndices().size();
569  Nf.resize(nb_row_dofs);
570  bzero(&*Nf.data().begin(),data.getIndices().size()*sizeof(FieldData));
571 
572 
573  double wAvenumber = dAta.aNgularfreq/dAta.sPeed;
574  double wAvenUmber = pow(wAvenumber,2.0);
575  dNdT_Im.resize(nb_row_dofs);
576 
577  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
578 
579  double val = getVolume()*getGaussPts()(3,gg);
580 
581  if(getHoGaussPtsDetJac().size()>0) {
582  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
583  }
584 
585  dNdT_Im.clear();
586  /* k^2 N^T q */
587  double const_kT_Im = wAvenUmber*commonData.pressureImAtGaussPts[gg];
588  /* diffN^T diffN q */
589  ublas::noalias(dNdT_Im) = prod(data.getDiffN(gg,nb_row_dofs),commonData.getGradImAtGaussPts(gg));
590  /* diffN^T diffN q - k^2 N^T N q lift of helmholtz operator*/
591 
592  cblas_daxpy(nb_row_dofs, -const_kT_Im, &data.getN(gg,nb_row_dofs)[0], 1, &dNdT_Im[0], 1);
593  /* insufficient, FIX ME*/
594  ublas::noalias(Nf) += val*(dNdT_Im);
595 
596  }
597 
598  if(useTsF) {
599  ierr = VecSetValues(getFEMethod()->ts_F,data.getIndices().size(),
600  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
601  } else {
602  ierr = VecSetValues(F,data.getIndices().size(),
603  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
604 
605  }
606 
607  } catch (const std::exception& ex) {
608  ostringstream ss;
609  ss << "throw in method: " << ex.what() << endl;
610  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
611  }
612 
613  PetscFunctionReturn(0);
614  }
615 
616  };
617 
618 
619  /** \biref operator to calculate right hand side of stiffness terms
620  * \infroup mofem_helmholtz_elem
621  */
623 
626  bool useTsF;
627  bool useReal;
628 
629  OpHelmholtzRhs_impedance(const string re_field_name,const string im_field_name,BlockData &data,CommonData &common_data,bool use_real):
630  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
631  dAta(data),commonData(common_data),useReal(use_real),useTsF(true) {}
632 
633  Vec F;
634  OpHelmholtzRhs_impedance(const string re_field_name,const string im_field_name,Vec _F,BlockData &data,CommonData &common_data,bool use_real):
635  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
636  dAta(data),commonData(common_data),useReal(use_real),useTsF(false),F(_F) { }
637 
639 
640  /* tip:all vectors should be declared as class members instead of as local function variables. */
641 
642 
643  /** \brief calculate helmholtz operator apply on lift.
644  *
645  * F = int i K T N^T dS
646  *
647  */
648  PetscErrorCode doWork(
649  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
650  PetscFunctionBegin;
651 
652  try {
653 
654  if(data.getIndices().size()==0) PetscFunctionReturn(0);
655  if(dAta.tRis.find(getMoFEMFEPtr()->get_ent())==dAta.tRis.end()) PetscFunctionReturn(0);
656 
657  PetscErrorCode ierr;
658 
659  /* moFEM tip: use "row_data.getIndices().size()" instead of "row_data.getN().size2()" for nb_row_dofs
660  the former can cause problems with mix approximation orders */
661 
662  const FENumeredDofMoFEMEntity *dof_ptr;
663  ierr = getMoFEMFEPtr()->get_row_dofs_by_petsc_gloabl_dof_idx(data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
664  int rank = dof_ptr->get_max_rank();
665  int nb_row_dofs = data.getIndices().size()/rank;
666 
667  Nf.resize(nb_row_dofs);
668  bzero(&*Nf.data().begin(),data.getIndices().size()*sizeof(FieldData));
669 
670  //wave number K is the propotional to the frequency of incident wave
671  //and represents number of waves per wave length 2Pi - 2Pi/K
672  double k = dAta.aNgularfreq/dAta.sPeed;
673  bool ho_geometry = true; //FIXME
674 
675  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
676 
677  double val = getGaussPts()(2,gg);
678 
679  if(ho_geometry) {
680  double area = norm_2(getNormals_at_GaussPt(gg))*0.5;
681  val *= area;
682  } else {
683  val *= getArea();
684  }
685 
686 
687  if(useReal) {
688 
689  ublas::noalias(Nf) += val*( - k*commonData.pressureImAtGaussPts(gg)*data.getN(gg,nb_row_dofs));
690 
691  } else{
692 
693  ublas::noalias(Nf) += val*(k*commonData.pressureReAtGaussPts(gg)*data.getN(gg,nb_row_dofs));
694 
695 
696  }
697 
698 
699  }
700 
701  if(useTsF) {
702  ierr = VecSetValues(getFEMethod()->ts_F,data.getIndices().size(),
703  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
704  } else {
705  ierr = VecSetValues(F,data.getIndices().size(),
706  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
707 
708  }
709 
710  } catch (const std::exception& ex) {
711  ostringstream ss;
712  ss << "throw in method: " << ex.what() << endl;
713  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
714  }
715 
716  PetscFunctionReturn(0);
717  }
718 
719  };
720 
721 
722  /** \biref operator to calculate right hand side wave source term F
723  * \infroup mofem_helmholtz_elem
724  */
726 
730  bool useTsF;
731  bool useScalar;
732  bool useReal;
733 
734  OpHelmholtzRhs_F(const string re_field_name,const string im_field_name,BlockData &data,CommonData &common_data,MatrixDouble &ho_coords,bool usescalar,bool use_real):
735  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
736  dAta(data),commonData(common_data),useTsF(false),hoCoords(ho_coords),useScalar(usescalar),useReal(use_real) {} //here useTsF(ture) originally.
737 
738  Vec F;
739  OpHelmholtzRhs_F(const string re_field_name,const string im_field_name,Vec _F,BlockData &data,CommonData &common_data,MatrixDouble &ho_coords,bool usescalar,bool use_real):
740  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
741  dAta(data),commonData(common_data),useTsF(false),F(_F),hoCoords(ho_coords),useScalar(usescalar),useReal(use_real) { }
742 
744 
745  /** \brief calculate Helmholtz RHS source term.
746  *
747  * F = int diffN^T F(x,y,z) dOmega^2
748  *
749  */
750 
751  PetscErrorCode doWork(
752  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
753  PetscFunctionBegin;
754 
755  if(dAta.tEts.find(getMoFEMFEPtr()->get_ent()) == dAta.tEts.end()) {
756  PetscFunctionReturn(0);
757  }
758 
759 
760  try {
761 
762  PetscErrorCode ierr;
763  if(data.getIndices().size()==0) PetscFunctionReturn(0);
764  int nb_row = data.getIndices().size();
765  Nf.resize(nb_row);
766  bzero(&Nf[0],nb_row*sizeof(double));
767  double x,y,z;
768  double iNcidentwave;
769  double wAvenumber = dAta.aNgularfreq/dAta.sPeed;
770 
771  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
772 
773  double val = getVolume()*getGaussPts()(3,gg);
774 
775  if(getHoGaussPtsDetJac().size()>0) {
776  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
777  x = hoCoords(gg,0);
778  y = hoCoords(gg,1);
779  z = hoCoords(gg,2);
780 
781  } else {
782  x = getCoordsAtGaussPts()(gg,0);
783  y = getCoordsAtGaussPts()(gg,1);
784  z = getCoordsAtGaussPts()(gg,2);
785  }
786 
787  VectorDouble F_S;
788  //Scalar value source term F = a
789  if(useScalar) {
790  const double f_scalar = -1; //scalar source vector F=1, wait for modification.
791  ublas::scalar_vector<double> F (data.getN().size1(),f_scalar); //scalar source vector
792  F_S = F;
793  }
794 
795 
796  //Functional value source term F =f(x)
797 
798 
799  else {
800 
801 
802  /***** incident wave in x direction *****/
803  /*** incident wave from Finite Element Analysis of Acoustic Scattering by Frank Ihlenburg **/
804 
805  double theta = atan2(y,x)+2*M_PI; //the arctan of radians (y/x)
806 
807  const double k = wAvenumber; //Wave number
808  const double a = 0.5; //radius of the sphere,wait to modify by user
809  const double const1 = k * a;
810 
811  const complex< double > i( 0.0, 1.0 );
812 
813  //// magnitude of incident wave
814  const double phi_incident_mag = 1.0;
815 
816  const double tol = 1.0e-10;
817  double max = 0.0;
818  double min = 999999.0;
819  complex< double > result = 0.0;
820  complex< double > prev_result;
821  double error = 100.0;
822 
823 
824  /*** Cyclindrical incident wave ***/
825  //double R = sqrt(pow(x,2.0)+pow(y,2.0)); //radius
826  double R = 0.5;
827  unsigned int n = 1; //initialized the infinite series loop
828 
829  double Jn_der_zero = ( - cyl_bessel_j( 1, const1 ));
830 
831  //n=0;
832  result -= Jn_der_zero;
833 
834  while( error > tol ) //finding the acoustic potential in one single point.
835  {
836  prev_result = result;
837  //The derivative of bessel function
838  double Jn_der = (n / const1 * cyl_bessel_j( n, const1 ) - cyl_bessel_j( n + 1, const1 ));
839 
840  result -= 2.0 * pow( i, n ) * Jn_der * cos(n*theta);
841  error = abs( abs( result ) - abs( prev_result ) );
842  ++n;
843  }
844 
845  /** End **/
846  //result = i * k * cos( theta ) * exp( i * k * R * cos( theta ) ); //derivative of incident wave
847 
848  //cout << "\n rhs_F is running \n" << endl;
849 
850  if(useReal){
851  iNcidentwave = std::real(k*result);
852 
853  } else if(!useReal) {
854  iNcidentwave = std::imag(k*result);
855 
856  }
857 
858  }
859 
860  val *= iNcidentwave;
861  ublas::noalias(Nf) += val*data.getN(gg,nb_row);
862 
863  }
864 
865 
866 
867  if(useTsF) {
868  ierr = VecSetValues(getFEMethod()->ts_F,data.getIndices().size(),
869  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
870  } else {
871  ierr = VecSetValues(F,data.getIndices().size(),
872  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
873 
874  }
875 
876 
877  } catch (const std::exception& ex) {
878  ostringstream ss;
879  ss << "throw in method: " << ex.what() << endl;
880  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
881  }
882 
883  PetscFunctionReturn(0);
884  }
885 
886  };
887 
888 
889  /** \biref operator to calculate left hand side of stiffness terms
890  * \infroup mofem_helmholtz_elem
891  */
893 
896  bool cAlculate;
897 
898  bool useTsB;
899  OpHelmholtzLhs_A(const string re_field_name,const string im_field_name,BlockData &data,CommonData &common_data,bool calculate):
900  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
901  dAta(data),commonData(common_data),useTsB(true),cAlculate(calculate) { }
902 
903  Mat A;
904  OpHelmholtzLhs_A(const string re_field_name,const string im_field_name,Mat _A,BlockData &data,CommonData &common_data,bool calculate):
905  TetElementForcesAndSourcesCore::UserDataOperator(re_field_name,im_field_name),
906  dAta(data),commonData(common_data),useTsB(false),A(_A),cAlculate(calculate) {}
907 
909 
910  /** \brief calculate helmholtz stiffness matrix
911  *
912  * K = int diffN^T * diffN^T - K^2 N^T N dOmega^2
913  *
914  */
915  PetscErrorCode doWork(
916  int row_side,int col_side,
917  EntityType row_type,EntityType col_type,
920  PetscFunctionBegin;
921 
922  if(dAta.tEts.find(getMoFEMFEPtr()->get_ent()) == dAta.tEts.end()) {
923  PetscFunctionReturn(0);
924  }
925 
926  MatrixDouble *K_ptr;
927 
928  try {
929 
930  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
931  if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
932 
933 
934  switch(row_type) {
935 
936  case MBVERTEX:
937  switch(col_type) {
938  case MBVERTEX:
939  K_ptr = &commonData.Ann;
940  break;
941  case MBTET:
942  K_ptr = &commonData.Anv;
943  break;
944  default:
945  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCT,"row_type = MBVERTEX, error col_type not match");
946  break;
947  }
948  break;
949 
950  case MBEDGE:
951  switch(col_type) {
952  case MBVERTEX:
953  commonData.Aen.resize(6);
954  K_ptr = &commonData.Aen[row_side];
955  break;
956  case MBTET:
957  commonData.Aev.resize(6);
958  K_ptr = &commonData.Aev[row_side];
959  break;
960  case MBEDGE:
961  commonData.Aee.resize(6,6);
962  K_ptr = &commonData.Aee(row_side,col_side);
963  break;
964  case MBTRI:
965  commonData.Aef.resize(6,4);
966  K_ptr = &commonData.Aef(row_side,col_side);
967  break;
968 
969  default:
970  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCT,"row_type = MBEDGE, error col_type not match");
971  break;
972  }
973  break;
974 
975  case MBTRI:
976  switch(col_type) {
977  case MBVERTEX:
978  commonData.Afn.resize(4);
979  K_ptr = &commonData.Afn[row_side];
980  break;
981  case MBTET:
982  commonData.Afv.resize(4);
983  K_ptr = &commonData.Afv[row_side];
984  break;
985  case MBTRI:
986  commonData.Aff.resize(4,4);
987  K_ptr = &commonData.Aff(row_side,col_side);
988  break;
989  default:
990  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCT,"row_type = MBTRI, error col_type not match");
991  break;
992  }
993  break;
994 
995  case MBTET:
996  K_ptr = &commonData.Avv;
997  break;
998  default:
999  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCT,"row_type not match");
1000  break;
1001  //never should happen
1002  }
1003 
1004  MatrixDouble &K = *K_ptr;
1005 
1006  int nb_row = row_data.getIndices().size();
1007  int nb_col = col_data.getIndices().size();
1008 
1009  if(cAlculate) {
1010 
1011  K.resize(nb_row,nb_col);
1012  bzero(&*K.data().begin(),nb_row*nb_col*sizeof(double));
1013 
1014  for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
1015 
1016  double wAvenumber = dAta.aNgularfreq/dAta.sPeed;
1017  double k_2 = pow(wAvenumber,2.0);
1018  /* check whether if double wAvenUmber = pow(dAta.aNgularfreq/dAta.sPeed,2.0); takes
1019  less time the the calculation above */
1020 
1021  double val = getVolume()*getGaussPts()(3,gg);
1022  if(getHoGaussPtsDetJac().size()>0) {
1023  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
1024  }
1025  noalias(K) += val*(prod(row_data.getDiffN(gg,nb_row),trans(col_data.getDiffN(gg,nb_col)))-k_2*outer_prod( row_data.getN(gg,nb_row),col_data.getN(gg,nb_col) ));
1026  }
1027 
1028  }
1029  PetscErrorCode ierr;
1030  if(!useTsB) {
1031  const_cast<FEMethod*>(getFEMethod())->ts_B = A; //FEMethod does not belong to MoFEM::Interface anymore.
1032  }
1033 
1034  ierr = MatSetValues(
1035  A, //(getFEMethod()->ts_B), instead in New thermal element. wait
1036  nb_row,&row_data.getIndices()[0],
1037  nb_col,&col_data.getIndices()[0],
1038  &K(0,0),ADD_VALUES); CHKERRQ(ierr);
1039 
1040  if(row_side != col_side || row_type != col_type) {
1041  transK.resize(nb_col,nb_row);
1042  noalias(transK) = trans( K );
1043  ierr = MatSetValues(
1044  A,
1045  nb_col,&col_data.getIndices()[0],
1046  nb_row,&row_data.getIndices()[0],
1047  &transK(0,0),ADD_VALUES); CHKERRQ(ierr);
1048  }
1049 
1050  } catch (const std::exception& ex) {
1051  ostringstream ss;
1052  ss << "throw in method: " << ex.what() << endl;
1053  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1054  }
1055 
1056  PetscFunctionReturn(0);
1057  }
1058 
1059  };
1060 
1061 
1062  /** \brief operator for calculate Impedance flux and assemble to right hand side
1063  * \infroup mofem_helmholtz_elem
1064  */
1066 
1070  bool useTsF;
1071  OpHelmholtzFlux(const string re_field_name,const string im_field_name,BlockData &DATA,FluxData &data,bool _ho_geometry = false):
1072  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1073  datA(DATA),dAta(data),ho_geometry(_ho_geometry),useTsF(true) { }
1074 
1075  Vec F;
1076  OpHelmholtzFlux(const string re_field_name,const string im_field_name,Vec _F,
1077  BlockData &DATA,FluxData &data,bool _ho_geometry = false):
1078  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1079  datA(DATA),dAta(data),ho_geometry(_ho_geometry),useTsF(false),F(_F) { }
1080 
1082 
1083  /** \brief calculate Impedance flux
1084  *
1085  * F = int_S N^T * g dS
1086  *
1087  */
1088  PetscErrorCode doWork(
1089  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
1090  PetscFunctionBegin;
1091 
1092  if(data.getIndices().size()==0) PetscFunctionReturn(0);
1093  if(dAta.tRis.find(getMoFEMFEPtr()->get_ent())==dAta.tRis.end()) PetscFunctionReturn(0);
1094 
1095  PetscErrorCode ierr;
1096 
1097  const FENumeredDofMoFEMEntity *dof_ptr;
1098  ierr = getMoFEMFEPtr()->get_row_dofs_by_petsc_gloabl_dof_idx(data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
1099  int rank = dof_ptr->get_max_rank();
1100 
1101  int nb_dofs = data.getIndices().size()/rank;
1102 
1103  Nf.resize(data.getIndices().size());
1104  fill(Nf.begin(),Nf.end(),0);
1105 
1106  //cerr << getNormal() << endl;
1107  //cerr << getNormals_at_GaussPt() << endl;
1108 
1109  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
1110 
1111  double val = getGaussPts()(2,gg);
1112  double flux;
1113  cout << "\n flux occurs \n" << endl;
1114 
1115  if(ho_geometry) {
1116  double area = norm_2(getNormals_at_GaussPt(gg)); //cblas_dnrm2(3,&getNormals_at_GaussPt()(gg,0),1);
1117  flux = dAta.dAta.data.value1*area; //FluxData.HeatfluxCubitBcData.data.value1 * area
1118  } else {
1119  flux = dAta.dAta.data.value1*getArea();
1120  }
1121  //cblas_daxpy(nb_row_dofs,val*flux,&data.getN()(gg,0),1,&*Nf.data().begin(),1);
1122  ublas::noalias(Nf) += val*flux*data.getN(gg,nb_dofs);
1123 
1124  }
1125 
1126  if(useTsF) {
1127  ierr = VecSetValues(getFEMethod()->ts_F,data.getIndices().size(),
1128  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
1129  } else {
1130  ierr = VecSetValues(F,data.getIndices().size(),
1131  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
1132  }
1133 
1134  PetscFunctionReturn(0);
1135  }
1136 
1137  };
1138 
1139 
1140  /** \brief operator for calculate Impedance flux and assemble to right hand side
1141  * \infroup mofem_helmholtz_elem
1142  */
1144 
1149  bool useTsF;
1150  bool useReal;
1151 
1152  OpHelmholtzIncidentWave(const string re_field_name,const string im_field_name,BlockData &DATA,IncidentData &data,MatrixDouble &ho_coords,bool use_real,bool _ho_geometry = false):
1153  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1154  datA(DATA),dAta(data),ho_geometry(_ho_geometry),useTsF(true),hoCoords(ho_coords),useReal(use_real) { }
1155 
1156  Vec F;
1157  OpHelmholtzIncidentWave(const string re_field_name,const string im_field_name,Vec _F,
1158  BlockData &DATA,IncidentData &data,MatrixDouble &ho_coords,bool use_real,bool _ho_geometry = false):
1159  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1160  datA(DATA),dAta(data),ho_geometry(_ho_geometry),useTsF(false),F(_F),hoCoords(ho_coords),useReal(use_real) { }
1161 
1163 
1164  /** \brief calculate Incident wave on rigid scatterer
1165  *
1166  * F = int_S N^T * g gradient (exp(i*k*x dot d) ) * n dS , F = -du_inc/dr represents sound hard wave (mere reflecting no absorption)
1167  *
1168  */
1169  PetscErrorCode doWork(
1170  int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
1171  PetscFunctionBegin;
1172 
1173  if(data.getIndices().size()==0) PetscFunctionReturn(0);
1174  if(dAta.tRis.find(getMoFEMFEPtr()->get_ent())==dAta.tRis.end()) PetscFunctionReturn(0);
1175 
1176  PetscErrorCode ierr;
1177 
1178  const FENumeredDofMoFEMEntity *dof_ptr;
1179  ierr = getMoFEMFEPtr()->get_row_dofs_by_petsc_gloabl_dof_idx(data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
1180  int rank = dof_ptr->get_max_rank();
1181 
1182  int nb_dofs = data.getIndices().size();
1183 
1184  Nf.resize(data.getIndices().size());
1185  fill(Nf.begin(),Nf.end(),0);
1186 
1187  //cerr << getNormal() << endl;
1188  //cerr << getNormals_at_GaussPt() << endl;
1189  double wAvenumber = datA.aNgularfreq/datA.sPeed;
1190  unsigned int nb_gauss_pts = data.getN().size1();
1191 
1192  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
1193 
1194  double val = getGaussPts()(2,gg);
1195 
1196  double iNcidentwave;
1197  ////////////////////********************************/////////////
1198 
1199  double x,y,z;
1200 
1201 
1202  if(hoCoords.size1() == nb_gauss_pts) {
1203 
1204 
1205  double area = norm_2(getNormals_at_GaussPt(gg))*0.5;
1206  val *= area;
1207  x = hoCoords(gg,0);
1208  y = hoCoords(gg,1);
1209  z = hoCoords(gg,2);
1210  } else {
1211  val *= getArea();
1212  x = getCoordsAtGaussPts()(gg,0);
1213  y = getCoordsAtGaussPts()(gg,1);
1214  z = getCoordsAtGaussPts()(gg,2);
1215  }
1216  /***** incident wave in x direction *****/
1217  /*** incident wave from Finite Element Analysis of Acoustic Scattering by Frank Ihlenburg **/
1218  //const double pi = atan( 1.0 ) * 4.0;
1219 
1220  //double theta = atan2(y,x); //the arctan of radians (y/x)
1221  double theta = atan2(y,x)+2*M_PI;
1222  //cout << "\n theta = \n" << theta << "\n M_PI = \N" << M_PI << endl;
1223  const double k = wAvenumber; //Wave number
1224  const double a = 0.5; //radius of the sphere,wait to modify by user
1225  const double const1 = k * a;
1226 
1227  const complex< double > i( 0.0, 1.0 );
1228 
1229  //// magnitude of incident wave
1230  const double phi_incident_mag = 1.0;
1231 
1232  const double tol = 1.0e-10;
1233  double max = 0.0;
1234  double min = 999999.0;
1235  complex< double > result = 0.0;
1236  complex< double > prev_result;
1237  double error = 100.0;
1238 
1239  /**** Spherical incident wave ***/
1240 
1241  unsigned int n = 0; //initialized the infinite series loop
1242 
1243  while( error > tol ) //finding the acoustic potential
1244  {
1245  double jn_der = ( n / const1 * sph_bessel( n, const1 ) - sph_bessel( n + 1, const1 ) ); //The derivative of bessel function
1246 
1247  double Pn = legendre_p( n, cos( theta ) ); //Legendre
1248 
1249  prev_result = result;
1250 
1251  result += pow( i, n ) * ( 2.0 * n + 1.0 ) * Pn * jn_der; //edition from Papers
1252  error = abs( abs( result ) - abs( prev_result ) );
1253  ++n;
1254  }
1255  /** End **/
1256 
1257  /*** Cyclindrical incident wave ***/
1258  //double R = sqrt(pow(x,2.0)+pow(y,2.0)); //radius
1259 
1260  //unsigned int n = 1; //initialized the infinite series loop
1261 
1262  //double Jn_der_zero = ( - cyl_bessel_j( 1, const1 ));
1263 
1264  ////n=0;
1265  //result -= Jn_der_zero;
1266  //
1267  //while( error > tol ) //finding the acoustic potential in one single point.
1268  //{
1269  // prev_result = result;
1270  // //The derivative of bessel function
1271  // double Jn_der = (n / const1 * cyl_bessel_j( n, const1 ) - cyl_bessel_j( n + 1, const1 ));
1272  //
1273  // result -= 2.0 * pow( i, n ) * Jn_der * cos(n*theta);
1274  // error = abs( abs( result ) - abs( prev_result ) );
1275  // ++n;
1276  //}
1277 
1278  /** End **/
1279  //result = i * k * cos( theta ) * exp( i * k * R * cos( theta ) ); //derivative of incident wave
1280 
1281 
1282  if(useReal){
1283  iNcidentwave = std::real(result);
1284 
1285  } else if(!useReal) {
1286  iNcidentwave = std::imag(result);
1287 
1288  }
1289 
1290 
1291  ublas::noalias(Nf) += val*iNcidentwave*data.getN(gg,nb_dofs);
1292  }
1293 
1294  if(useTsF) {
1295  ierr = VecSetValues(getFEMethod()->ts_F,data.getIndices().size(),
1296  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
1297  } else {
1298  ierr = VecSetValues(F,data.getIndices().size(),
1299  &data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);
1300  }
1301 
1302  PetscFunctionReturn(0);
1303  }
1304 
1305  };
1306 
1307 
1308  /// \biref operator to calculate Impedance on body surface and assemble to imaginary lhs of equations
1310 
1315  bool useTsB;
1316 
1317  OpImpedanceLhs_reimC(const string re_field_name,const string im_field_name,
1318  ImpedanceData &data,BlockData &block_data,MatrixDouble &ho_coords,bool _ho_geometry = false):
1319  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1320  dAta(data),ho_geometry(_ho_geometry),useTsB(true),blockData(block_data),hoCoords(ho_coords) {}
1321 
1322  Mat A;
1323  OpImpedanceLhs_reimC(const string re_field_name,const string im_field_name,Mat _A,
1324  ImpedanceData &data,BlockData &block_data,MatrixDouble &ho_coords,bool _ho_geometry = false):
1325  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1326  dAta(data),ho_geometry(_ho_geometry),useTsB(false),A(_A),blockData(block_data),hoCoords(ho_coords) {}
1327 
1329  /** \brief calculate helmholtz Impedance term in the imaginary lhs of equations
1330  *
1331  * K = intS N^T i sIgma K N dS
1332  */
1333  PetscErrorCode doWork(
1334  int row_side,int col_side,
1335  EntityType row_type,EntityType col_type,
1338  PetscFunctionBegin;
1339 
1340  try {
1341 
1342  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
1343  if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
1344 
1345  int nb_row = row_data.getIndices().size();
1346  int nb_col = col_data.getIndices().size();
1347  K.resize(nb_row,nb_col);
1348  bzero(&*K.data().begin(),nb_row*nb_col*sizeof(double));
1349 
1350  for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) { /* Integrate the shape functions in one element */
1351  double waveNumber = blockData.aNgularfreq/blockData.sPeed;
1352  double val = getGaussPts()(2,gg);
1353 
1354  unsigned int nb_gauss_pts = row_data.getN().size1();
1355  double impedanceConst = dAta.sIgma*waveNumber;
1356  if(hoCoords.size1() == nb_gauss_pts) {
1357  double area = norm_2(getNormals_at_GaussPt(gg))*0.5;
1358  val *= area;
1359  } else {
1360  val *= getArea();
1361  }
1362  noalias(K) += val*impedanceConst*outer_prod( row_data.getN(gg,nb_row),col_data.getN(gg,nb_col) );
1363 
1364  }
1365 
1366  PetscErrorCode ierr;
1367  if(!useTsB) {
1368  const_cast<FEMethod*>(getFEMethod())->ts_B = A; //to change the pointer getFEMethod()'s private member ts_B value with A.
1369  } //getFEMethod() belong to class MoFEM::Interface::FEMethod
1370  ierr = MatSetValues(
1371  (getFEMethod()->ts_B),
1372  nb_row,&row_data.getIndices()[0], //MatSetValues(A,i,m,j,n,value,INSERT_VALUES,ierr), i*j the size of input matrix,
1373  nb_col,&col_data.getIndices()[0], // m and n - the position of parent matrix to allowed the input matrix.
1374  &K(0,0),ADD_VALUES); CHKERRQ(ierr);
1375  if(row_side != col_side || row_type != col_type) {
1376  transK.resize(nb_col,nb_row);
1377  noalias(transK) = trans( K );
1378  ierr = MatSetValues(
1379  (getFEMethod()->ts_B),
1380  nb_col,&col_data.getIndices()[0],
1381  nb_row,&row_data.getIndices()[0],
1382  &transK(0,0),ADD_VALUES); CHKERRQ(ierr);
1383  }
1384 
1385 
1386  } catch (const std::exception& ex) {
1387  ostringstream ss;
1388  ss << "throw in method: " << ex.what() << endl;
1389  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1390  }
1391 
1392  PetscFunctionReturn(0);
1393  }
1394 
1395  };
1396 
1397 
1398  /// \biref operator to calculate Impedance on body surface and assemble to imaginary lhs of equations
1400 
1405  bool useTsB;
1406 
1407  OpImpedanceLhs_imreC(const string re_field_name,const string im_field_name,
1408  ImpedanceData &data,BlockData &block_data,MatrixDouble &ho_coords,bool _ho_geometry = false):
1409  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1410  dAta(data),ho_geometry(_ho_geometry),useTsB(true),blockData(block_data),hoCoords(ho_coords) {}
1411 
1412  Mat A;
1413  OpImpedanceLhs_imreC(const string re_field_name,const string im_field_name,Mat _A,
1414  ImpedanceData &data,BlockData &block_data,MatrixDouble &ho_coords,bool _ho_geometry = false):
1415  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1416  dAta(data),ho_geometry(_ho_geometry),useTsB(false),A(_A),blockData(block_data),hoCoords(ho_coords) {}
1417 
1419  /** \brief calculate helmholtz Impedance term in the imaginary lhs of equations
1420  *
1421  * K = - intS N^T i sIgma K N dS
1422  */
1423  PetscErrorCode doWork(
1424  int row_side,int col_side,
1425  EntityType row_type,EntityType col_type,
1428  PetscFunctionBegin;
1429 
1430  try {
1431 
1432  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
1433  if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
1434 
1435  int nb_row = row_data.getIndices().size();
1436  int nb_col = col_data.getIndices().size();
1437  K.resize(nb_row,nb_col);
1438  bzero(&*K.data().begin(),nb_row*nb_col*sizeof(double));
1439 
1440  for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
1441 
1442  double waveNumber = blockData.aNgularfreq/blockData.sPeed;
1443  double val = getGaussPts()(2,gg);
1444 
1445  unsigned int nb_gauss_pts = row_data.getN().size1();
1446  double impedanceConst = -dAta.sIgma*waveNumber;
1447  if(hoCoords.size1() == nb_gauss_pts) {
1448  double area = norm_2(getNormals_at_GaussPt(gg))*0.5;
1449  val *= area;
1450  } else {
1451  val *= getArea();
1452  }
1453  noalias(K) += val*impedanceConst*outer_prod( row_data.getN(gg,nb_row),col_data.getN(gg,nb_col) );
1454 
1455  }
1456 
1457  PetscErrorCode ierr;
1458  if(!useTsB) {
1459  const_cast<FEMethod*>(getFEMethod())->ts_B = A; //to change the pointer getFEMethod()'s private member ts_B value with A.
1460  } //getFEMethod() belong to class MoFEM::Interface::FEMethod
1461  ierr = MatSetValues(
1462  (getFEMethod()->ts_B),
1463  nb_row,&row_data.getIndices()[0], //MatSetValues(A,i,m,j,n,value,INSERT_VALUES,ierr), i*j the size of input matrix,
1464  nb_col,&col_data.getIndices()[0], // m and n - the position of parent matrix to allowed the input matrix.
1465  &K(0,0),ADD_VALUES); CHKERRQ(ierr);
1466  if(row_side != col_side || row_type != col_type) {
1467  transK.resize(nb_col,nb_row);
1468  noalias(transK) = trans( K );
1469  ierr = MatSetValues(
1470  (getFEMethod()->ts_B),
1471  nb_col,&col_data.getIndices()[0],
1472  nb_row,&row_data.getIndices()[0],
1473  &transK(0,0),ADD_VALUES); CHKERRQ(ierr);
1474  }
1475 
1476 
1477  } catch (const std::exception& ex) {
1478  ostringstream ss;
1479  ss << "throw in method: " << ex.what() << endl;
1480  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1481  }
1482 
1483  PetscFunctionReturn(0);
1484  }
1485 
1486  };
1487 
1489 
1494  bool useTsB;
1495 
1496  OpImpedanceLhs_D(const string re_field_name,const string im_field_name,
1497  ImpedanceData &data,BlockData &block_data,MatrixDouble &ho_coords,bool _ho_geometry = false):
1498  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1499  dAta(data),ho_geometry(_ho_geometry),useTsB(true),blockData(block_data),hoCoords(ho_coords) {}
1500 
1501  Mat A;
1502  OpImpedanceLhs_D(const string re_field_name,const string im_field_name,Mat _A,
1503  ImpedanceData &data,BlockData &block_data,MatrixDouble &ho_coords,bool _ho_geometry = false):
1504  TriElementForcesAndSurcesCore::UserDataOperator(re_field_name,im_field_name),
1505  dAta(data),ho_geometry(_ho_geometry),useTsB(false),A(_A),blockData(block_data),hoCoords(ho_coords) {}
1506 
1508  /** \brief calculate helmholtz Impedance term in the real part of lhs equations
1509  *
1510  * K = intS 1/R N^T N dS
1511  * see <<Plane-wave basis finite elements and boundary elements for three-dimensional wave scattering>>
1512  * by E. Perrey-Debain, O. Laghrouche, P. Bettess and J. Trevelyan for details.
1513  */
1514  PetscErrorCode doWork(
1515  int row_side,int col_side,
1516  EntityType row_type,EntityType col_type,
1519  PetscFunctionBegin;
1520 
1521  try {
1522 
1523  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
1524  if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
1525 
1526  int nb_row = row_data.getIndices().size();
1527  int nb_col = col_data.getIndices().size();
1528  K.resize(nb_row,nb_col);
1529  bzero(&*K.data().begin(),nb_row*nb_col*sizeof(double));
1530 
1531  for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
1532  double val = getGaussPts()(2,gg);
1533 
1534  unsigned int nb_gauss_pts = row_data.getN().size1();
1535  double x,y,z;
1536  if(hoCoords.size1() == nb_gauss_pts) {
1537  double area = norm_2(getNormals_at_GaussPt(gg))*0.5;
1538  val *= area;
1539  x = hoCoords(gg,0);
1540  y = hoCoords(gg,1);
1541  z = hoCoords(gg,2);
1542  } else {
1543  val *= getArea();
1544  x = getCoordsAtGaussPts()(gg,0);
1545  y = getCoordsAtGaussPts()(gg,1);
1546  z = getCoordsAtGaussPts()(gg,2);
1547  }
1548  double r_inv = 1/(2*(sqrt(pow(x,2.0)+pow(y,2.0)))); //inverse radius
1549  noalias(K) += val*r_inv*outer_prod( row_data.getN(gg,nb_row),col_data.getN(gg,nb_col) );
1550 
1551  }
1552 
1553 
1554  PetscErrorCode ierr;
1555  if(!useTsB) {
1556  const_cast<FEMethod*>(getFEMethod())->ts_B = A;
1557  }
1558  ierr = MatSetValues(
1559  (getFEMethod()->ts_B),
1560  nb_row,&row_data.getIndices()[0],
1561  nb_col,&col_data.getIndices()[0],
1562  &K(0,0),ADD_VALUES); CHKERRQ(ierr);
1563  if(row_side != col_side || row_type != col_type) {
1564  transK.resize(nb_col,nb_row);
1565  noalias(transK) = trans( K );
1566  ierr = MatSetValues(
1567  (getFEMethod()->ts_B),
1568  nb_col,&col_data.getIndices()[0],
1569  nb_row,&row_data.getIndices()[0],
1570  &transK(0,0),ADD_VALUES); CHKERRQ(ierr);
1571  }
1572 
1573 
1574  } catch (const std::exception& ex) {
1575  ostringstream ss;
1576  ss << "throw in method: " << ex.what() << endl;
1577  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1578  }
1579 
1580  PetscFunctionReturn(0);
1581  }
1582 
1583  };
1584 
1585 
1586  /** \brief add helmholtz element on tets
1587  * \infroup mofem_helmholtz_elem
1588  *
1589  * It get data from block set and define element in moab
1590  *w
1591  * \param problem name
1592  * \param field name
1593  * \param name of mesh nodal positions (if not defined nodal coordinates are used)
1594  */
1595  PetscErrorCode addHelmholtzElements(
1596  const string problem_name,const string re_field_name,const string im_field_name,const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
1597  PetscFunctionBegin;
1598 
1599  PetscErrorCode ierr;
1600  ErrorCode rval;
1601 
1602  //to do matrix A and C; add Finite Element "HELMHOLTZ_FE"
1603  ierr = mField.add_finite_element("HELMHOLTZ_FE",MF_ZERO); CHKERRQ(ierr );
1604  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_FE",re_field_name); CHKERRQ(ierr);
1605  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_FE",re_field_name); CHKERRQ(ierr);
1606  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FE",re_field_name); CHKERRQ(ierr);
1607  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_FE",im_field_name); CHKERRQ(ierr);
1608  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_FE",im_field_name); CHKERRQ(ierr);
1609  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FE",im_field_name); CHKERRQ(ierr);
1610 
1611  if(mField.check_field(mesh_nodals_positions)) {
1612  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FE",mesh_nodals_positions); CHKERRQ(ierr);
1613  }
1614  ierr = mField.modify_problem_add_finite_element(problem_name,"HELMHOLTZ_FE"); CHKERRQ(ierr);
1615 
1616  //takes skin of block of entities
1617  //Skinner skin(&mField.get_moab());
1618  // loop over all blocksets and get data which name is FluidPressure
1619 
1621 
1622  if(it->get_Cubit_name().compare(0,13,"MAT_HELMHOLTZ") == 0) {
1623 
1624  //get block attributes
1625  vector<double> attributes;
1626  ierr = it->get_Cubit_attributes(attributes); CHKERRQ(ierr);
1627  if(attributes.size()<2) {
1628  SETERRQ1(PETSC_COMM_SELF,1,"not enough block attributes to deffine fluid pressure element, attributes.size() = %d ",attributes.size());
1629  }
1630 
1631  blockData.aNgularfreq = attributes[0];
1632  blockData.sPeed = attributes[1];
1633  setOfBlocks[it->get_msId()].aNgularfreq = attributes[0];
1634  setOfBlocks[it->get_msId()].sPeed = attributes[1];
1635 
1636  rval = mField.get_moab().get_entities_by_type(it->meshset,MBTET,setOfBlocks[it->get_msId()].tEts,true); CHKERR_PETSC(rval);
1637  ierr = mField.add_ents_to_finite_element_by_TETs(setOfBlocks[it->get_msId()].tEts,"HELMHOLTZ_FE"); CHKERRQ(ierr);
1638  }
1639  }
1640 
1641  PetscFunctionReturn(0);
1642  }
1643 
1644  /** \brief add helmholtz flux element
1645  * \infroup mofem_helmholtz_elem
1646  *
1647  * It get data from helmholtz flux set and define elemenet in moab. Alternatively
1648  * uses block set with name HELMHOLTZ_FLUX.
1649  *
1650  * \param problem name
1651  * \param field name
1652  * \param name of mesh nodal positions (if not defined nodal coordinates are used)
1653  */
1654  PetscErrorCode addHelmholtzFluxElement(
1655  const string problem_name,const string re_field_name,const string im_field_name) {
1656  PetscFunctionBegin;
1657 
1658  PetscErrorCode ierr;
1659  ErrorCode rval;
1660 
1661  ierr = mField.add_finite_element("HELMHOLTZ_FLUX_FE",MF_ZERO); CHKERRQ(ierr);
1662  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_FLUX_FE",re_field_name); CHKERRQ(ierr);
1663  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_FLUX_FE",re_field_name); CHKERRQ(ierr);
1664  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FLUX_FE",re_field_name); CHKERRQ(ierr);
1665  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_FLUX_FE",im_field_name); CHKERRQ(ierr);
1666  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_FLUX_FE",im_field_name); CHKERRQ(ierr);
1667  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FLUX_FE",im_field_name); CHKERRQ(ierr);
1668 
1669  if(mField.check_field("MESH_NODE_POSITIONS")) {
1670  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FLUX_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1671  }
1672 
1673  ierr = mField.modify_problem_add_finite_element(problem_name,"HELMHOLTZ_FLUX_FE"); CHKERRQ(ierr);
1674 
1676  ierr = it->get_cubit_bc_data_structure(setOfFluxes[it->get_msId()].dAta); CHKERRQ(ierr);
1677  rval = mField.get_moab().get_entities_by_type(it->meshset,MBTRI,setOfFluxes[it->get_msId()].tRis,true); CHKERR_PETSC(rval);
1678  ierr = mField.add_ents_to_finite_element_by_TRIs(setOfFluxes[it->get_msId()].tRis,"HELMHOLTZ_FLUX_FE"); CHKERRQ(ierr);
1679  }
1680 
1681  //this is alternative method for setting boundary conditions, to bypass bu in cubit file reader.
1682  //not elegant, but good enough
1684  if(it->get_Cubit_name().compare(0,9,"HELMHOLTZ_FLUX") == 0) {
1685  vector<double> data;
1686  ierr = it->get_Cubit_attributes(data); CHKERRQ(ierr);
1687  if(data.size()!=1) {
1688  SETERRQ(PETSC_COMM_SELF,1,"Data inconsistency");
1689  }
1690  strcpy(setOfFluxes[it->get_msId()].dAta.data.name,"HelmholtzFlux");
1691  setOfFluxes[it->get_msId()].dAta.data.flag1 = 1;
1692  setOfFluxes[it->get_msId()].dAta.data.value1 = data[0];
1693  //cerr << setOfFluxes[it->get_msId()].dAta << endl;
1694  rval = mField.get_moab().get_entities_by_type(it->meshset,MBTRI,setOfFluxes[it->get_msId()].tRis,true); CHKERR_PETSC(rval);
1695  ierr = mField.add_ents_to_finite_element_by_TRIs(setOfFluxes[it->get_msId()].tRis,"HELMHOLTZ_FLUX_FE"); CHKERRQ(ierr);
1696 
1697  }
1698  }
1699 
1700 
1701  PetscFunctionReturn(0);
1702  }
1703 
1704 
1705 
1706  /** \brief add helmholtz flux element
1707  * \infroup mofem_helmholtz_elem
1708  *
1709  * It get data from helmholtz flux set and define elemenet in moab. Alternatively
1710  * uses block set with name HELMHOLTZ_FLUX.
1711  *
1712  * \param problem name
1713  * \param field name
1714  * \param name of mesh nodal positions (if not defined nodal coordinates are used)
1715  */
1717  const string problem_name,const string re_field_name,const string im_field_name) {
1718  PetscFunctionBegin;
1719 
1720  PetscErrorCode ierr;
1721  ErrorCode rval;
1722 
1723  ierr = mField.add_finite_element("HELMHOLTZ_INCIDENT_FE",MF_ZERO); CHKERRQ(ierr);
1724  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_INCIDENT_FE",re_field_name); CHKERRQ(ierr);
1725  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_INCIDENT_FE",re_field_name); CHKERRQ(ierr);
1726  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_INCIDENT_FE",re_field_name); CHKERRQ(ierr);
1727  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_INCIDENT_FE",im_field_name); CHKERRQ(ierr);
1728  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_INCIDENT_FE",im_field_name); CHKERRQ(ierr);
1729  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_INCIDENT_FE",im_field_name); CHKERRQ(ierr);
1730 
1731  if(mField.check_field("MESH_NODE_POSITIONS")) {
1732  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_INCIDENT_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1733  }
1734 
1735  ierr = mField.modify_problem_add_finite_element(problem_name,"HELMHOLTZ_INCIDENT_FE"); CHKERRQ(ierr);
1736 
1737  //this is alternative method for setting boundary conditions, to bypass bu in cubit file reader.
1738  //not elegant, but good enough
1740  if(it->get_Cubit_name().compare(0,9,"INCIDENT") == 0) {
1741  vector<double> data;
1742  ierr = it->get_Cubit_attributes(data); CHKERRQ(ierr);
1743  if(data.size()!=1) {
1744  SETERRQ(PETSC_COMM_SELF,1,"Data inconsistency");
1745  }
1746 
1747  setOfIncident[it->get_msId()].amplitude = data[0];
1748 
1749  //cerr << setOfIncident[it->get_msId()].dAta << endl;
1750  rval = mField.get_moab().get_entities_by_type(it->meshset,MBTRI,setOfIncident[it->get_msId()].tRis,true); CHKERR_PETSC(rval);
1751  ierr = mField.add_ents_to_finite_element_by_TRIs(setOfIncident[it->get_msId()].tRis,"HELMHOLTZ_INCIDENT_FE"); CHKERRQ(ierr);
1752 
1753  }
1754  }
1755 
1756 
1757  PetscFunctionReturn(0);
1758  }
1759 
1760 
1761  /** \brief add Impedance element
1762  * \infroup mofem_helmholtz_elem
1763  *
1764  * It get data from convection set and define elemenet in moab. Alternatively
1765  * uses block set with name IMPEDANCE.
1766  *
1767  * \param problem name
1768  * \param field name
1769  * \param name of mesh nodal positions (if not defined nodal coordinates are used)
1770  */
1772  const string problem_name,const string re_field_name,const string im_field_name) {
1773  PetscFunctionBegin;
1774 
1775  PetscErrorCode ierr;
1776  ErrorCode rval;
1777 
1778  ierr = mField.add_finite_element("HELMHOLTZ_IMPEDANCE_FE",MF_ZERO); CHKERRQ(ierr);
1779  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_IMPEDANCE_FE",re_field_name); CHKERRQ(ierr);
1780  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_IMPEDANCE_FE",re_field_name); CHKERRQ(ierr);
1781  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_IMPEDANCE_FE",re_field_name); CHKERRQ(ierr);
1782  ierr = mField.modify_finite_element_add_field_row("HELMHOLTZ_IMPEDANCE_FE",im_field_name); CHKERRQ(ierr);
1783  ierr = mField.modify_finite_element_add_field_col("HELMHOLTZ_IMPEDANCE_FE",im_field_name); CHKERRQ(ierr);
1784  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_IMPEDANCE_FE",im_field_name); CHKERRQ(ierr);
1785 
1786  ierr = mField.modify_problem_add_finite_element(problem_name,"HELMHOLTZ_IMPEDANCE_FE"); CHKERRQ(ierr);
1787 
1789  if(it->get_Cubit_name().compare(0,9,"IMPEDANCE") == 0) {
1790 
1791  vector<double> data;
1792  ierr = it->get_Cubit_attributes(data); CHKERRQ(ierr);
1793  if(data.size()!=2) {
1794  SETERRQ(PETSC_COMM_SELF,1,"Data inconsistency");
1795  }
1796  setOfImpedance[it->get_msId()].g = data[0];
1797  setOfImpedance[it->get_msId()].sIgma = data[1];
1798  rval = mField.get_moab().get_entities_by_type(it->meshset,MBTRI,setOfImpedance[it->get_msId()].tRis,true); CHKERR_PETSC(rval);
1799  ierr = mField.add_ents_to_finite_element_by_TRIs(setOfImpedance[it->get_msId()].tRis,"HELMHOLTZ_IMPEDANCE_FE"); CHKERRQ(ierr);
1800 
1801  }
1802  }
1803 
1804 
1805  PetscFunctionReturn(0);
1806  }
1807 
1808 
1809  /** \brief this function is used in case of stationary problem to set elements for rhs
1810  * \infroup mofem_helmholtz_elem
1811  */
1812  PetscErrorCode setHelmholtzFiniteElementRhs_FOperators(string re_field_name,string im_field_name,Vec &F,bool useScalar,const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
1813  PetscFunctionBegin;
1814  map<int,BlockData>::iterator sit = setOfBlocks.begin();
1815 
1816  if(mField.check_field(mesh_nodals_positions)) {
1817 
1818  //feIncidentWave.get_op_to_do_Rhs().push_back(new OpHoCoordTri(mesh_nodals_positions,hoCoordsTri));
1819 
1820  feRhs.get_op_to_do_Rhs().push_back(new OpHoCoordTet(mesh_nodals_positions,hoCoordsTet));
1821 
1822  }
1823 
1824  for(;sit!=setOfBlocks.end();sit++) {
1825  //add finite element
1826  feRhs.get_op_to_do_Rhs().push_back(new OpHelmholtzRhs_F(re_field_name,re_field_name,F,sit->second,commonData,hoCoordsTet,useScalar,true));
1827  feRhs.get_op_to_do_Rhs().push_back(new OpHelmholtzRhs_F(im_field_name,im_field_name,F,sit->second,commonData,hoCoordsTet,useScalar,false));
1828  }
1829  PetscFunctionReturn(0);
1830  }
1831 
1832  PetscErrorCode setHelmholtzFiniteElementRhsOperators(string re_field_name,string im_field_name,Vec &F,bool useImpedance) {
1833  PetscFunctionBegin;
1834 
1835  map<int,BlockData>::iterator sit = setOfBlocks.begin();
1836 
1837  for(;sit!=setOfBlocks.end();sit++) {
1838 
1839  feRhs.get_op_to_do_Rhs().push_back(new OpGetGradReAtGaussPts(re_field_name,commonData));
1840  feRhs.get_op_to_do_Rhs().push_back(new OpGetTetPressureReAtGaussPts(re_field_name,commonData));
1841  feRhs.get_op_to_do_Rhs().push_back(new OpHelmholtzRhs_Re(re_field_name,re_field_name,F,sit->second,commonData));
1842 
1843  feRhs.get_op_to_do_Rhs().push_back(new OpGetGradImAtGaussPts(im_field_name,commonData));
1844  feRhs.get_op_to_do_Rhs().push_back(new OpGetTetPressureImAtGaussPts(im_field_name,commonData));
1845  feRhs.get_op_to_do_Rhs().push_back(new OpHelmholtzRhs_Im(im_field_name,im_field_name,F,sit->second,commonData));
1846 
1847  //cout << "\n useImpedance = \n" << useImpedance << std::endl;
1848  //FIXE ME, NOT WORKING
1849  //if(useImpedance) {
1850  // feIncidentWave.get_op_to_do_Rhs().push_back(new OpHelmholtzRhs_impedance(re_field_name,re_field_name,F,sit->second,commonData,true));
1851  // feIncidentWave.get_op_to_do_Rhs().push_back(new OpHelmholtzRhs_impedance(im_field_name,im_field_name,F,sit->second,commonData,false));
1852  //}
1853  //add finite element
1854  }
1855 
1856  PetscFunctionReturn(0);
1857  }
1858 
1859 
1860 
1861  /** \brief this fucntion is used in case of stationary Helmholtz problem for lhs stiffness term
1862  * \infroup mofem_helmholtz_elem
1863  */
1864  PetscErrorCode setHelmholtzFiniteElementLhsOperators(string re_field_name,string im_field_name,Mat A) {
1865  PetscFunctionBegin;
1866  map<int,BlockData>::iterator sit = setOfBlocks.begin();
1867  for(;sit!=setOfBlocks.end();sit++) {
1868  //add finite elemen
1869  feLhs.get_op_to_do_Lhs().push_back(new OpHelmholtzLhs_A(re_field_name,re_field_name,A,sit->second,commonData,true));
1870  feLhs.get_op_to_do_Lhs().push_back(new OpHelmholtzLhs_A(im_field_name,im_field_name,A,sit->second,commonData,false));
1871  }
1872  PetscFunctionReturn(0);
1873  }
1874 
1875 
1876  PetscErrorCode setHelmholtzFluxFiniteElementRhsOperators(string re_field_name,string im_field_name,Vec &F,const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
1877  PetscFunctionBegin;
1878  bool ho_geometry = false;
1879  if(mField.check_field(mesh_nodals_positions)) {
1880  ho_geometry = true;
1881  }
1882  map<int,FluxData>::iterator sit = setOfFluxes.begin();
1883  for(;sit!=setOfFluxes.end();sit++) {
1884  //add finite element
1885  feFlux.get_op_to_do_Rhs().push_back(new OpHelmholtzFlux(re_field_name,im_field_name,F,blockData,sit->second,ho_geometry));
1886  }
1887  PetscFunctionReturn(0);
1888  }
1889  /** \brief this function is used in case of statonary for helmholtz incident wave flux terms
1890 
1891  */
1892 
1893  PetscErrorCode setHelmholtzIncidentWaveFiniteElementRhsOperators(string re_field_name,string im_field_name,Vec &F,const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
1894  PetscFunctionBegin;
1895  bool ho_geometry = false;
1896  map<int,IncidentData>::iterator sit = setOfIncident.begin();
1897 
1898  if(mField.check_field(mesh_nodals_positions)) {
1899 
1900  feIncidentWave.get_op_to_do_Rhs().push_back(new OpHoCoordTri(mesh_nodals_positions,hoCoordsTri));
1901  }
1902 
1903  for(;sit!=setOfIncident.end();sit++) {
1904  //add finite element
1905  feIncidentWave.get_op_to_do_Rhs().push_back(new OpHelmholtzIncidentWave(re_field_name,re_field_name,F,blockData,sit->second,hoCoordsTri,true,ho_geometry));
1906  feIncidentWave.get_op_to_do_Rhs().push_back(new OpHelmholtzIncidentWave(im_field_name,im_field_name,F,blockData,sit->second,hoCoordsTri,false,ho_geometry));
1907 
1908  }
1909  PetscFunctionReturn(0);
1910  }
1911 
1912  PetscErrorCode setHelmholtzImpedanceFiniteElementLhsOperators(string re_field_name,string im_field_name,Mat A,const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
1913  PetscFunctionBegin;
1914  bool ho_geometry = false;
1915  //if(mField.check_field(mesh_nodals_positions)) {
1916  // ho_geometry = true;
1917  //}
1918  map<int,ImpedanceData>::iterator sit = setOfImpedance.begin();
1919 
1920  if(mField.check_field(mesh_nodals_positions)) {
1921  // ho_geometry = true;
1922  feIncidentWave.get_op_to_do_Rhs().push_back(new OpHoCoordTri(mesh_nodals_positions,hoCoordsTri));
1923  }
1924 
1925  for(;sit!=setOfImpedance.end();sit++) {
1926  //add finite element
1927  feImpedanceLhs.get_op_to_do_Lhs().push_back(new OpImpedanceLhs_D(re_field_name,re_field_name,A,sit->second,blockData,hoCoordsTri,ho_geometry));
1928  feImpedanceLhs.get_op_to_do_Lhs().push_back(new OpImpedanceLhs_D(im_field_name,im_field_name,A,sit->second,blockData,hoCoordsTri,ho_geometry));
1929  feImpedanceLhs.get_op_to_do_Lhs().push_back(new OpImpedanceLhs_reimC(re_field_name,im_field_name,A,sit->second,blockData,hoCoordsTri,ho_geometry));
1930  feImpedanceLhs.get_op_to_do_Lhs().push_back(new OpImpedanceLhs_imreC(im_field_name,re_field_name,A,sit->second,blockData,hoCoordsTri,ho_geometry));
1931  }
1932  PetscFunctionReturn(0);
1933  }
1934 
1935  };
1936 
1937  }
1938 
1939  #endif //__HELMHOLTZ_ELEMENT_HPP
1940 
1941  /***************************************************************************//**
1942  * \defgroup mofem_helmholtz_elem Helmholtz element
1943  * \ingroup user_modules
1944  ******************************************************************************/
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
calculate Incident wave on rigid scatterer
struture grouping operators and data used for helmholtz problemsIn order to assemble matrices and rig...
PetscErrorCode setHelmholtzFiniteElementRhs_FOperators(string re_field_name, string im_field_name, Vec &F, bool useScalar, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
this function is used in case of stationary problem to set elements for rhs \infroup mofem_helmholtz_...
ublas::matrix_row< MatrixDouble > getGradImAtGaussPts(const int gg)
OpGetTriPressureImAtGaussPts(const string field_name, CommonData &common_data)
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
PetscErrorCode addHelmholtzFluxElement(const string problem_name, const string re_field_name, const string im_field_name)
add helmholtz flux element \infroup mofem_helmholtz_elem
OpImpedanceLhs_D(const string re_field_name, const string im_field_name, ImpedanceData &data, BlockData &block_data, MatrixDouble &ho_coords, bool _ho_geometry=false)
MyVolumeFE & getLoopFeRhs()
get rhs volume element
map< int, BlockData > setOfBlocks
maps block set id with appropiate BlockData
operator to calculate pressure gradient at Gauss points
OpHelmholtzRhs_Re(const string re_field_name, const string im_field_name, Vec _F, BlockData &data, CommonData &common_data)
virtual moab::Interface & get_moab()=0
OpHelmholtzIncidentWave(const string re_field_name, const string im_field_name, BlockData &DATA, IncidentData &data, MatrixDouble &ho_coords, bool use_real, bool _ho_geometry=false)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSurcesCore::EntData &row_data, DataForcesAndSurcesCore::EntData &col_data)
calculate helmholtz Impedance term in the real part of lhs equations
operator for calculate Impedance flux and assemble to right hand side \infroup mofem_helmholtz_elem
map< int, FluxData > setOfFluxes
maps side set id with appropiate FluxData
OpGetFieldAtGaussPts(const string field_name, VectorDouble &field_at_gauss_pts)
PetscErrorCode setHelmholtzFiniteElementLhsOperators(string re_field_name, string im_field_name, Mat A)
this fucntion is used in case of stationary Helmholtz problem for lhs stiffness term \infroup mofem_h...
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
ublas::matrix_row< MatrixDouble > getGradReAtGaussPts(const int gg)
OpImpedanceLhs_imreC(const string re_field_name, const string im_field_name, ImpedanceData &data, BlockData &block_data, MatrixDouble &ho_coords, bool _ho_geometry=false)
OpHelmholtzLhs_A(const string re_field_name, const string im_field_name, Mat _A, BlockData &data, CommonData &common_data, bool calculate)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
operator calculating pressure gradients
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
Range tRis
those will be on body skin, except thos with contact whith other body where pressure is applied
PetscErrorCode addHelmholtzElements(const string problem_name, const string re_field_name, const string im_field_name, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
add helmholtz element on tets \infroup mofem_helmholtz_elem
virtual DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_TRIs(const Range &tris, const std::string &name)=0
add TRI entities from range to finite element database given by name
OpGetGradImAtGaussPts(const string field_name, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
calculate helmholtz operator apply on lift.
PetscErrorCode setHelmholtzFiniteElementRhsOperators(string re_field_name, string im_field_name, Vec &F, bool useImpedance)
Range tEts
constatins elements in block set
MyVolumeFE feRhs
cauclate right hand side for tetrahedral elements
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
calculate Helmholtz RHS source term.
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSurcesCore::EntData &row_data, DataForcesAndSurcesCore::EntData &col_data)
calculate helmholtz stiffness matrix
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
calculate helmholtz operator apply on lift.
OpImpedanceLhs_imreC(const string re_field_name, const string im_field_name, Mat _A, ImpedanceData &data, BlockData &block_data, MatrixDouble &ho_coords, bool _ho_geometry=false)
data for calulation Angular Frequency and wave speed elements \infroup mofem_helmholtz_elem
OpGetTetPressureReAtGaussPts(const string field_name, CommonData &common_data)
OpGetGradReAtGaussPts(const string field_name, CommonData &common_data)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
operator calculating pressure and rate of pressure
OpHelmholtzFlux(const string re_field_name, const string im_field_name, BlockData &DATA, FluxData &data, bool _ho_geometry=false)
map< int, ImpedanceData > setOfImpedance
Range tRis
suraface triangles where hate flux is applied
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
PetscErrorCode addHelmholtzImpedanceElement(const string problem_name, const string re_field_name, const string im_field_name)
add Impedance element \infroup mofem_helmholtz_elem
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
calculate helmholtz operator apply on lift.
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
operator to calculate pressure at Gauss pts \infroup mofem_helmholtz_elem
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MyTriFE feRhs_Tri
cauclate right hand side for tetrahedral elements
data for calulation heat flux \infroup mofem_helmholtz_elem
operator to calculate pressure gradient at Gauss points
common data used by volume elements \infroup mofem_helmholtz_elem
operator to calculate pressure at Gauss pts \infroup mofem_helmholtz_elem
int getRule(int order)
it is used to calculate nb. of Gauss integartion points
Range tRis
suraface triangles where hate flux is applied
OpGetTriPressureReAtGaussPts(const string field_name, CommonData &common_data)
double tol
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode setHelmholtzImpedanceFiniteElementLhsOperators(string re_field_name, string im_field_name, Mat A, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
OpImpedanceLhs_reimC(const string re_field_name, const string im_field_name, ImpedanceData &data, BlockData &block_data, MatrixDouble &ho_coords, bool _ho_geometry=false)
OpHelmholtzRhs_Im(const string re_field_name, const string im_field_name, Vec _F, BlockData &data, CommonData &common_data)
CHKERRQ(ierr)
DataForcesAndSourcesCore::EntData EntData
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
map< int, IncidentData > setOfIncident
maps side set id with appropiate FluxData
OpImpedanceLhs_reimC(const string re_field_name, const string im_field_name, Mat _A, ImpedanceData &data, BlockData &block_data, MatrixDouble &ho_coords, bool _ho_geometry=false)
OpHoCoordTri(const string field_name, MatrixDouble &ho_coords)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSurcesCore::EntData &row_data, DataForcesAndSurcesCore::EntData &col_data)
calculate helmholtz Impedance term in the imaginary lhs of equations
OpHelmholtzRhs_impedance(const string re_field_name, const string im_field_name, BlockData &data, CommonData &common_data, bool use_real)
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
OpHelmholtzRhs_F(const string re_field_name, const string im_field_name, BlockData &data, CommonData &common_data, MatrixDouble &ho_coords, bool usescalar, bool use_real)
OpHelmholtzRhs_Re(const string re_field_name, const string im_field_name, BlockData &data, CommonData &common_data)
static moab::Error error
OpHelmholtzRhs_impedance(const string re_field_name, const string im_field_name, Vec _F, BlockData &data, CommonData &common_data, bool use_real)
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpHelmholtzRhs_F(const string re_field_name, const string im_field_name, Vec _F, BlockData &data, CommonData &common_data, MatrixDouble &ho_coords, bool usescalar, bool use_real)
MyVolumeFE & getLoopFeLhs()
get lhs volume element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
PetscErrorCode setHelmholtzFluxFiniteElementRhsOperators(string re_field_name, string im_field_name, Vec &F, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
double FieldData
Field data type.
Definition: Types.hpp:36
PetscErrorCode setHelmholtzIncidentWaveFiniteElementRhsOperators(string re_field_name, string im_field_name, Vec &F, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
this function is used in case of statonary for helmholtz incident wave flux terms
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
\biref operator to calculate Impedance on body surface and assemble to imaginary lhs of equations
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MyTriFE & getLoopFeRhs_Tri()
get rhs volume element
virtual DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_TETs(const Range &tets, const std::string &name)=0
add TET entities from range to finite element database given by name
OpHelmholtzRhs_Im(const string re_field_name, const string im_field_name, BlockData &data, CommonData &common_data)
data for Robin Boundary condition \infroup mofem_helmholtz_elem
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSurcesCore::EntData &row_data, DataForcesAndSurcesCore::EntData &col_data)
calculate helmholtz Impedance term in the imaginary lhs of equations
PetscErrorCode addHelmholtzIncidentElement(const string problem_name, const string re_field_name, const string im_field_name)
add helmholtz flux element \infroup mofem_helmholtz_elem
HeatfluxCubitBcData dAta
for more details look to BCMultiIndices.hpp to see details of HeatfluxCubitBcData
virtual bool check_field(const std::string &name) const =0
check if field is in database
\biref operator to calculate Impedance on body surface and assemble to imaginary lhs of equations
OpImpedanceLhs_D(const string re_field_name, const string im_field_name, Mat _A, ImpedanceData &data, BlockData &block_data, MatrixDouble &ho_coords, bool _ho_geometry=false)
opearator to caulate pressure and rate of pressure at Gauss points \infroup mofem_helmholtz_elem
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
OpGetTetPressureImAtGaussPts(const string field_name, CommonData &common_data)
OpHelmholtzIncidentWave(const string re_field_name, const string im_field_name, Vec _F, BlockData &DATA, IncidentData &data, MatrixDouble &ho_coords, bool use_real, bool _ho_geometry=false)
operator for calculate Impedance flux and assemble to right hand side \infroup mofem_helmholtz_elem
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
HelmholtzElement(MoFEM::Interface &m_field)
OpHoCoordTet(const string field_name, MatrixDouble &ho_coords)
#define _F(n)
Definition: quad.c:25
OpHelmholtzLhs_A(const string re_field_name, const string im_field_name, BlockData &data, CommonData &common_data, bool calculate)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
calculate Impedance flux
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
operator calculating pressure gradients
double amplitude
for more details look to BCMultiIndices.hpp to see details of HeatfluxCubitBcData
OpHelmholtzFlux(const string re_field_name, const string im_field_name, Vec _F, BlockData &DATA, FluxData &data, bool _ho_geometry=false)