v0.9.0
Public Member Functions | Public Attributes | List of all members
MoFEM::HelmholtzElement::OpHelmholtzIncidentWave Struct Reference

operator for calculate Impedance flux and assemble to right hand side \infroup mofem_helmholtz_elem More...

#include <users_modules/helmholtz/src/HelmholtzElementObsolete.hpp>

Inherits UserDataOperator.

Collaboration diagram for MoFEM::HelmholtzElement::OpHelmholtzIncidentWave:
[legend]

Public Member Functions

 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)
 
 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)
 
PetscErrorCode doWork (int side, EntityType type, DataForcesAndSurcesCore::EntData &data)
 calculate Incident wave on rigid scatterer More...
 

Public Attributes

MatrixDoublehoCoords
 
BlockDatadatA
 
IncidentDatadAta
 
bool ho_geometry
 
bool useTsF
 
bool useReal
 
Vec F
 
VectorDouble Nf
 

Detailed Description

operator for calculate Impedance flux and assemble to right hand side \infroup mofem_helmholtz_elem

Definition at line 1143 of file HelmholtzElementObsolete.hpp.

Constructor & Destructor Documentation

◆ OpHelmholtzIncidentWave() [1/2]

MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::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 
)

◆ OpHelmholtzIncidentWave() [2/2]

MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::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 
)

Member Function Documentation

◆ doWork()

PetscErrorCode MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::doWork ( int  side,
EntityType  type,
DataForcesAndSurcesCore::EntData data 
)

calculate Incident wave on rigid scatterer

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)

Definition at line 1169 of file HelmholtzElementObsolete.hpp.

1170  {
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  }
Range tRis
suraface triangles where hate flux is applied
double tol
CHKERRQ(ierr)
static moab::Error error
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.

Member Data Documentation

◆ datA

BlockData& MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::datA

Definition at line 1146 of file HelmholtzElementObsolete.hpp.

◆ dAta

IncidentData& MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::dAta

Definition at line 1147 of file HelmholtzElementObsolete.hpp.

◆ F

Vec MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::F

Definition at line 1156 of file HelmholtzElementObsolete.hpp.

◆ ho_geometry

bool MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::ho_geometry

Definition at line 1148 of file HelmholtzElementObsolete.hpp.

◆ hoCoords

MatrixDouble& MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::hoCoords

Definition at line 1145 of file HelmholtzElementObsolete.hpp.

◆ Nf

VectorDouble MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::Nf

Definition at line 1162 of file HelmholtzElementObsolete.hpp.

◆ useReal

bool MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::useReal

Definition at line 1150 of file HelmholtzElementObsolete.hpp.

◆ useTsF

bool MoFEM::HelmholtzElement::OpHelmholtzIncidentWave::useTsF

Definition at line 1149 of file HelmholtzElementObsolete.hpp.


The documentation for this struct was generated from the following file: