v0.10.0
Classes | Functions | Variables
fe_approximation_point_source.cpp File Reference
#include <MoFEM.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/ptr_container/ptr_map.hpp>
#include <DirichletBC.hpp>
#include <PostProcOnRefMesh.hpp>
#include <ReynoldsStress.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <cmath>
#include <boost/math/special_functions.hpp>
#include <complex>
#include <FieldApproximation.hpp>
#include <boost/shared_array.hpp>
#include <kiss_fft.h>
#include <kiss_fft.c>
#include <AnalyticalSolutions.hpp>
#include <AnalyticalDirichlet.hpp>
#include <HelmholtzElement.hpp>
#include <TimeSeries.hpp>

Go to the source code of this file.

Classes

struct  PlaneIncidentWaveSacttrerData
 Read impulse and apply FFTFirst signal is read from text file. Currently is assumed that this signal gives profile for plane wave. Spherical or other wave types can be easily added and implemented. More...
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 76 of file fe_approximation_point_source.cpp.

76  {
77 
78 
79 
80 
81  PetscInitialize(&argc,&argv,(char *)0,help);
82 
83  //Core mb_instance;
84  moab::Core mb_instance;
85  Interface& moab = mb_instance;
86  int rank;
87  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
88 
89  PetscBool flg = PETSC_TRUE;
90  char mesh_file_name[255];
91  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
92  if(flg != PETSC_TRUE) {
93  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"*** ERROR -my_file (MESH FILE NEEDED)");
94  }
95 
96  // Create moab parallel communicator
97  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
98  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
99 
100  PetscBool is_partitioned = PETSC_FALSE;
101  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
102  if(is_partitioned) {
103  //Read mesh to MOAB
104  const char *option;
105  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
106  rval = moab.load_file(mesh_file_name, 0, option); CHKERR_PETSC(rval);
107  } else {
108  const char *option;
109  option = "";
110  rval = moab.load_file(mesh_file_name, 0, option); CHKERR_PETSC(rval);
111  }
112 
113  // Create MoFEM (cephas) database
114  MoFEM::Core core(moab);
115  MoFEM::Interface& m_field = core;
116 
117  // Get start time for analyse
118  PetscLogDouble t1,t2;
119  PetscLogDouble v1,v2;
120  ierr = PetscTime(&v1); CHKERRQ(ierr);
121  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
122 
123  //set entitities bit level
124  BitRefLevel bit_level0;
125  bit_level0.set(0);
126  EntityHandle meshset_level0;
127  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERR_PETSC(rval);
128  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
129 
130  //Fields
131  ierr = m_field.add_field("rePRES",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
132  ierr = m_field.add_field("imPRES",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
133  ierr = m_field.add_field("P",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr); // in time domain
134 
135  //meshset consisting all entities in mesh
136  EntityHandle root_set = moab.get_root_set();
137  //add entities to field
138  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"rePRES"); CHKERRQ(ierr);
139  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"imPRES"); CHKERRQ(ierr);
140  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"P"); CHKERRQ(ierr);
141 
142  //set app. order
143  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
144  PetscInt order;
145  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
146  if(flg != PETSC_TRUE) {
147  order = 2;
148  }
149 
150  ierr = m_field.set_field_order(root_set,MBTET,"rePRES",order); CHKERRQ(ierr);
151  ierr = m_field.set_field_order(root_set,MBTRI,"rePRES",order); CHKERRQ(ierr);
152  ierr = m_field.set_field_order(root_set,MBEDGE,"rePRES",order); CHKERRQ(ierr);
153  ierr = m_field.set_field_order(root_set,MBVERTEX,"rePRES",1); CHKERRQ(ierr);
154  ierr = m_field.set_field_order(root_set,MBTET,"imPRES",order); CHKERRQ(ierr);
155  ierr = m_field.set_field_order(root_set,MBTRI,"imPRES",order); CHKERRQ(ierr);
156  ierr = m_field.set_field_order(root_set,MBEDGE,"imPRES",order); CHKERRQ(ierr);
157  ierr = m_field.set_field_order(root_set,MBVERTEX,"imPRES",1); CHKERRQ(ierr);
158 
159  ierr = m_field.set_field_order(root_set,MBTET,"P",order); CHKERRQ(ierr);
160  ierr = m_field.set_field_order(root_set,MBTRI,"P",order); CHKERRQ(ierr);
161  ierr = m_field.set_field_order(root_set,MBEDGE,"P",order); CHKERRQ(ierr);
162  ierr = m_field.set_field_order(root_set,MBVERTEX,"P",1); CHKERRQ(ierr);
163 
164  if(!m_field.check_field("MESH_NODE_POSITIONS")) {
165 
166  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
167  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
168  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
169  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
170  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
171  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
172 
173  ierr = m_field.build_fields(); CHKERRQ(ierr);
174  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
175  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
176 
177  } else {
178 
179  ierr = m_field.build_fields(); CHKERRQ(ierr);
180 
181  }
182 
183  // Finite Elements
184 
185  HelmholtzElement helmholtz_element(m_field);
186  ierr = helmholtz_element.getGlobalParametersFromLineCommandOptions(); CHKERRQ(ierr);
187  ierr = helmholtz_element.addHelmholtzElements("rePRES","imPRES"); CHKERRQ(ierr);
188  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
189  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","reEX"); CHKERRQ(ierr);
190  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","imEX"); CHKERRQ(ierr);
191  }
192 
193  bool Dirichlet_bc_set = false;
194  Range bc_dirichlet_tris,analytical_bc_tris;
195  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"ANALYTICAL_BC",it)) {
196  rval = moab.get_entities_by_type(it->get_meshset(),MBTRI,analytical_bc_tris,true); CHKERR_PETSC(rval);
197  Dirichlet_bc_set = true;
198  }
199  bc_dirichlet_tris.merge(analytical_bc_tris);
200  AnalyticalDirichletBC analytical_bc_real(m_field);
201  AnalyticalDirichletBC analytical_bc_imag(m_field);
202  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",analytical_bc_tris); CHKERRQ(ierr);
203  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",analytical_bc_tris); CHKERRQ(ierr);
204 
205  PetscBool wavenumber_flg;
206  double wavenumber;
207  // set wave number from line command, that overwrite numbre form block set
208  ierr = PetscOptionsGetScalar(NULL,"-wave_number",&wavenumber,&wavenumber_flg); CHKERRQ(ierr);
209  if(!wavenumber_flg) {
210 
211  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"wave number not given, set in line command -wave_number to fix problem");
212 
213  }
214 
215  double power_of_incident_wave = 1;
216  ierr = PetscOptionsGetScalar(NULL,"-power_of_incident_wave",&power_of_incident_wave,NULL); CHKERRQ(ierr);
217 
218 
219  // This is added for a case than on some surface, defined by the user a
220  // incident plane wave is scattered.
221  map<int,PlaneIncidentWaveSacttrerData> planeWaveScatterData;
222  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"SOFT_INCIDENT_WAVE_BC",it)) {
223  rval = moab.get_entities_by_type(it->get_meshset(),MBTRI,planeWaveScatterData[it->getMeshsetId()].tRis,true); CHKERR_PETSC(rval);
224  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",planeWaveScatterData[it->getMeshsetId()].tRis); CHKERRQ(ierr);
225  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",planeWaveScatterData[it->getMeshsetId()].tRis); CHKERRQ(ierr);
226  bc_dirichlet_tris.merge(planeWaveScatterData[it->getMeshsetId()].tRis);
227 
228  Dirichlet_bc_set = true;
229 
230  }
231 
232  //ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
233  //// Build adjacencies
234  //ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
235 
236  // Problem
237  ierr = m_field.add_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
238  ierr = m_field.add_problem("BCREAL_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for real field
239  ierr = m_field.add_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for imag field
240 
241  // Set refinement level for problem
242  ierr = m_field.modify_problem_ref_level_add_bit("ACOUSTIC_PROBLEM",bit_level0); CHKERRQ(ierr);
243  ierr = m_field.modify_problem_ref_level_add_bit("BCREAL_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
244  ierr = m_field.modify_problem_ref_level_add_bit("BCIMAG_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
245 
246  // Add elements to problems
247  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
248  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_IMIM_FE"); CHKERRQ(ierr);
249  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_REIM_FE"); CHKERRQ(ierr);
250 
251  ierr = m_field.modify_problem_add_finite_element("BCREAL_PROBLEM","BCREAL_FE"); CHKERRQ(ierr);
252  ierr = m_field.modify_problem_add_finite_element("BCIMAG_PROBLEM","BCIMAG_FE"); CHKERRQ(ierr);
253 
254 
255  // Build problems
256  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
257  // Build adjacencies
258  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
259 
260  // build porblems
261  if(is_partitioned) {
262  // if mesh is partitioned
263 
264  ierr = m_field.build_partitioned_problem("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
265  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
266 
267  if(Dirichlet_bc_set) {
268  ierr = m_field.build_partitioned_problem("BCREAL_PROBLEM",true); CHKERRQ(ierr);
269  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM",true); CHKERRQ(ierr);
270 
271  ierr = m_field.build_partitioned_problem("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
272  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
273 
274  }
275 
276 
277  } else {
278  // if not partitioned mesh is load to all processes
279 
280  ierr = m_field.build_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
281  ierr = m_field.partition_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
282  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
283 
284  if(Dirichlet_bc_set) {
285  ierr = m_field.build_problem("BCREAL_PROBLEM"); CHKERRQ(ierr);
286  ierr = m_field.partition_problem("BCREAL_PROBLEM"); CHKERRQ(ierr);
287  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM"); CHKERRQ(ierr);
288 
289 
290 
291  ierr = m_field.build_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr);
292  ierr = m_field.partition_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr);
293  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM"); CHKERRQ(ierr);
294  }
295  }
296 
297  ierr = m_field.partition_ghost_dofs("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
298  if(Dirichlet_bc_set) {
299  ierr = m_field.partition_ghost_dofs("BCREAL_PROBLEM"); CHKERRQ(ierr);
300  ierr = m_field.partition_ghost_dofs("BCIMAG_PROBLEM"); CHKERRQ(ierr);
301  }
302 
303  // Get problem matrices and vectors
304  Vec F; //Right hand side vector
305  ierr = m_field.VecCreateGhost("ACOUSTIC_PROBLEM",ROW,&F); CHKERRQ(ierr);
306  Vec T; //Solution vector
307  ierr = VecDuplicate(F,&T); CHKERRQ(ierr);
308  Mat A; //Left hand side matrix
309  ierr = m_field.MatCreateMPIAIJWithArrays("ACOUSTIC_PROBLEM",&A); CHKERRQ(ierr);
310  ierr = helmholtz_element.setOperators(A,F,"rePRES","imPRES"); CHKERRQ(ierr);
311 
312  //wave direction unit vector=[x,y,z]^T
313  VectorDouble wave_direction;
314  wave_direction.resize(3);
315  wave_direction.clear();
316  wave_direction[2] = 1; // default:X direction [0,0,1]
317 
318  int nmax = 3;
319  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-wave_direction",&wave_direction[0],&nmax,NULL); CHKERRQ(ierr);
320  if(nmax > 0 && nmax != 3) {
321  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"*** ERROR -wave_direction [3*1 vector] default:X direction [0,0,1]");
322  }
323 
324  PetscInt choise_value = NO_ANALYTICAL_SOLUTION;
325  // set type of analytical solution
326  ierr = PetscOptionsGetEList(NULL,"-analytical_solution_type",analytical_solution_types,6,&choise_value,NULL); CHKERRQ(ierr);
327 
328  switch((AnalyticalSolutionTypes)choise_value) {
329 
331 
332  {
333  double scattering_sphere_radius = 0.5;
334  ierr = PetscOptionsGetScalar(NULL,"-scattering_sphere_radius",&scattering_sphere_radius,NULL); CHKERRQ(ierr);
335 
336 
337 
338  boost::shared_ptr<HardSphereScatterWave> function_evaluator = boost::shared_ptr<HardSphereScatterWave>(
339  new HardSphereScatterWave(wavenumber,scattering_sphere_radius)
340  );
341  ierr = analytical_bc_real.setApproxOps(
342  m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL
343  ); CHKERRQ(ierr);
344  ierr = analytical_bc_imag.setApproxOps(
345  m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG
346  ); CHKERRQ(ierr);
347  Dirichlet_bc_set = true;
348 
349  }
350 
351  break;
352 
354 
355  {
356 
357  double scattering_sphere_radius = 0.5;
358  ierr = PetscOptionsGetScalar(NULL,"-scattering_sphere_radius",&scattering_sphere_radius,NULL); CHKERRQ(ierr);
359 
360 
361  boost::shared_ptr<SoftSphereScatterWave> function_evaluator = boost::shared_ptr<SoftSphereScatterWave>(new SoftSphereScatterWave(wavenumber,scattering_sphere_radius));
362  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
363  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
364  Dirichlet_bc_set = true;
365 
366  }
367 
368  break;
369 
370  case PLANE_WAVE:
371 
372  {
373 
374  double angle = 0.25;
375  // set wave number from line command, that overwrite numbre form block set
376  ierr = PetscOptionsGetScalar(NULL,"-wave_guide_angle",&angle,NULL); CHKERRQ(ierr);
377 
378 
379  boost::shared_ptr<PlaneWave> function_evaluator = boost::shared_ptr<PlaneWave>(new PlaneWave(wavenumber,angle*M_PI));
380  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
381  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
382  Dirichlet_bc_set = true;
383 
384  }
385 
386  break;
387 
389 
390  {
391 
392  boost::shared_ptr<HardCylinderScatterWave> function_evaluator = boost::shared_ptr<HardCylinderScatterWave>(new HardCylinderScatterWave(wavenumber));
393  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
394  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
395  Dirichlet_bc_set = true;
396 
397 
398  }
399 
400  break;
401 
403 
404  {
405 
406  boost::shared_ptr<SoftCylinderScatterWave> function_evaluator = boost::shared_ptr<SoftCylinderScatterWave>(new SoftCylinderScatterWave(wavenumber));
407  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
408  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
409  Dirichlet_bc_set = true;
410 
411  }
412 
413  break;
414 
415  case INCIDENT_WAVE:
416 
417  {
418 
419  boost::shared_ptr<IncidentWave> function_evaluator =
420  boost::shared_ptr<IncidentWave>(new IncidentWave(wavenumber,wave_direction,power_of_incident_wave));
421  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
422  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
423  Dirichlet_bc_set = true;
424 
425  }
426 
427  break;
428 
430 
431  {
432  // Dirichlet_bc_set = false;
433  }
434 
435  break;
436 
437  }
438 
439  // Analytical boundary conditions
440  AnalyticalDirichletBC::DirichletBC analytical_ditihlet_bc_real(m_field,"rePRES",A,T,F);
441  AnalyticalDirichletBC::DirichletBC analytical_ditihlet_bc_imag(m_field,"imPRES",A,T,F);
442 
443  if(Dirichlet_bc_set) {
444 
445  {
446 
447  map<int,PlaneIncidentWaveSacttrerData>::iterator mit = planeWaveScatterData.begin();
448  for(;mit!=planeWaveScatterData.end();mit++) {
449 
450  // note negative field, scatter field should cancel incident wave
451  boost::shared_ptr<IncidentWave> function_evaluator = boost::shared_ptr<IncidentWave>(
452  new IncidentWave(wavenumber,wave_direction,-power_of_incident_wave)
453  );
454  ierr = analytical_bc_real.setApproxOps(
455  m_field,"rePRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::REAL
456  ); CHKERRQ(ierr);
457  ierr = analytical_bc_imag.setApproxOps(
458  m_field,"imPRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::IMAG
459  ); CHKERRQ(ierr);
460 
461  }
462 
463  }
464  // Solve for analytical Dirichlet bc dofs
465  ierr = analytical_bc_real.setProblem(m_field,"BCREAL_PROBLEM"); CHKERRQ(ierr);
466  ierr = analytical_bc_imag.setProblem(m_field,"BCIMAG_PROBLEM"); CHKERRQ(ierr);
467  ierr = analytical_bc_real.solveProblem(
468  m_field,"BCREAL_PROBLEM","BCREAL_FE",analytical_ditihlet_bc_real,bc_dirichlet_tris
469  ); CHKERRQ(ierr);
470  ierr = analytical_bc_imag.solveProblem(
471  m_field,"BCIMAG_PROBLEM","BCIMAG_FE",analytical_ditihlet_bc_imag,bc_dirichlet_tris
472  ); CHKERRQ(ierr);
473 
474  ierr = analytical_bc_real.destroyProblem(); CHKERRQ(ierr);
475  ierr = analytical_bc_imag.destroyProblem(); CHKERRQ(ierr);
476 
477  }
478 
479 
480  PetscBool monochromatic_wave = PETSC_TRUE;
481  ierr = PetscOptionsGetBool(PETSC_NULL,"-monochromatic_wave",&monochromatic_wave,NULL); CHKERRQ(ierr);
482 
483  PetscBool add_incident_wave = PETSC_FALSE;
484  ierr = PetscOptionsGetBool(NULL,"-add_incident_wave",&add_incident_wave,NULL); CHKERRQ(ierr);
485  if(add_incident_wave) {
486 
487  // define problem
488  ierr = m_field.add_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
489  // set finite elements for problem
490  ierr = m_field.modify_problem_add_finite_element("INCIDENT_WAVE","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
491  // set refinment level for problem
492  ierr = m_field.modify_problem_ref_level_add_bit("INCIDENT_WAVE",bit_level0); CHKERRQ(ierr);
493 
494  // build porblems
495  if(is_partitioned) {
496  // if mesh is partitioned
497  ierr = m_field.build_partitioned_problem("INCIDENT_WAVE",true); CHKERRQ(ierr);
498  ierr = m_field.partition_finite_elements("INCIDENT_WAVE",true); CHKERRQ(ierr);
499  } else {
500  // if not partitioned mesh is load to all processes
501  ierr = m_field.build_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
502  ierr = m_field.partition_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
503  ierr = m_field.partition_finite_elements("INCIDENT_WAVE"); CHKERRQ(ierr);
504  }
505  ierr = m_field.partition_ghost_dofs("INCIDENT_WAVE"); CHKERRQ(ierr);
506 
507  }
508 
509 
510  KSP solver;
511  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
512 
513  if(monochromatic_wave) {
514 
515  // Zero vectors
516  ierr = VecZeroEntries(T); CHKERRQ(ierr);
517  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
518  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
519  ierr = VecZeroEntries(F); CHKERRQ(ierr);
520  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
521  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
522  ierr = MatZeroEntries(A); CHKERRQ(ierr);
523 
524  // Assemble problem
525  if(Dirichlet_bc_set) {
526  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
527  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
528  }
529 
530  ierr = helmholtz_element.calculateA("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
531  ierr = helmholtz_element.calculateF("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
532 
533  if(Dirichlet_bc_set) {
534  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
535  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
536  }
537 
538  ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
539  ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
540  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
541  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
542  ierr = VecScale(F,-1); CHKERRQ(ierr);
543 
544  // Solve problem
545  ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
546  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
547  ierr = KSPSetUp(solver); CHKERRQ(ierr);
548 
549  ierr = KSPSolve(solver,F,T); CHKERRQ(ierr);
550 
551  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
552  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
553 
554 
555  //Save data on mesh
556  if(is_partitioned) {
557 
558  // no need for global communication
559  ierr = m_field.set_local_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
560 
561  } else {
562 
563  ierr = m_field.set_global_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
564 
565  }
566 
567  } else {
568 
569 
570  // define problem
571  ierr = m_field.add_problem("PRESSURE_IN_TIME"); CHKERRQ(ierr);
572  // set finite elements for problem
573  ierr = m_field.modify_problem_add_finite_element("PRESSURE_IN_TIME","PRESSURE_FE"); CHKERRQ(ierr);
574  // set refinment level for problem
575  ierr = m_field.modify_problem_ref_level_add_bit("PRESSURE_IN_TIME",bit_level0); CHKERRQ(ierr);
576 
577  // build porblems
578  if(is_partitioned) {
579  // if mesh is partitioned
580  ierr = m_field.build_partitioned_problem("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
581  ierr = m_field.partition_finite_elements("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
582  } else {
583  // if not partitioned mesh is load to all processes
584  ierr = m_field.build_problem("PRESSURE_IN_TIME"); CHKERRQ(ierr);
585  ierr = m_field.partition_problem("PRESSURE_IN_TIME"); CHKERRQ(ierr);
586  ierr = m_field.partition_finite_elements("PRESSURE_IN_TIME"); CHKERRQ(ierr);
587  }
588  ierr = m_field.partition_ghost_dofs("PRESSURE_IN_TIME"); CHKERRQ(ierr);
589 
590  TimeSeries time_series(m_field,helmholtz_element,
591  analytical_ditihlet_bc_real,analytical_ditihlet_bc_imag,
592  Dirichlet_bc_set);
593 
594  ierr = time_series.readData(); CHKERRQ(ierr);
595  ierr = time_series.createPressureSeries(T); CHKERRQ(ierr);
596  ierr = time_series.forwardSpaceDft(); CHKERRQ(ierr);
597  ierr = time_series.pressureForwardDft(); CHKERRQ(ierr);
598  ierr = time_series.solveForwardDFT(solver,A,F,T); CHKERRQ(ierr);
599  ierr = time_series.pressureInverseDft(); CHKERRQ(ierr);
600  ierr = time_series.generateReferenceElementMesh(); CHKERRQ(ierr);
601  ierr = time_series.saveResults(); CHKERRQ(ierr);
602  ierr = time_series.destroyPressureSeries(); CHKERRQ(ierr);
603 
604  }
605 
606 
607  //Vec P,M;
608  //ierr = m_field.VecCreateGhost("EX1_PROBLEM",COL,&M); CHKERRQ(ierr);
609  //ierr = VecDuplicate(M,&P); CHKERRQ(ierr);
610 
611  //ierr = m_field.set_local_ghost_vector("EX1_PROBLEM",COL,M,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
612  //ierr = m_field.set_other_global_ghost_vector("EX1_PROBLEM","reEX","imEX",COL,P,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
613 
614  //double nrm2_M;
615  //ierr = VecNorm(M,NORM_2,&nrm2_M); CHKERRQ(ierr);
616 
617  //Vec V;
618  //ierr = m_field.VecCreateGhost("ACOUSTIC_PROBLEM",COL,&V); CHKERRQ(ierr);
619  //ierr = m_field.set_local_ghost_vector("ACOUSTIC_PROBLEM",COL,V,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
620 
621  //VecScatter scatter_real,scatter_imag;
622 
623  //ierr = m_field.VecScatterCreate(V,"ACOUSTIC_PROBLEM","rePRES",COL,M,"EX1_PROBLEM","reEX",COL,&scatter_real); CHKERRQ(ierr);
624 
625  //ierr = m_field.VecScatterCreate(V,"ACOUSTIC_PROBLEM","imPRES",COL,P,"EX1_PROBLEM","reEX",COL,&scatter_imag); CHKERRQ(ierr);
626 
627  //VecScale(V,-1);
628 
629  //ierr = VecScatterBegin(scatter_real,V,M,ADD_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
630  //ierr = VecScatterEnd(scatter_real,V,M,ADD_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
631 
632  //double nrm2_ErrM;
633  //ierr = VecNorm(M,NORM_2,&nrm2_ErrM); CHKERRQ(ierr);
634  //PetscPrintf(PETSC_COMM_WORLD,"L2 relative error on real field of acoustic problem %6.4e\n",(nrm2_ErrM)/(nrm2_M));
635 
636 
637  //ierr = VecDestroy(&M); CHKERRQ(ierr);
638  //ierr = VecDestroy(&P); CHKERRQ(ierr);
639  //ierr = VecDestroy(&V); CHKERRQ(ierr);
640 
641  // Destroy the KSP solvers
642  ierr = MatDestroy(&A); CHKERRQ(ierr);
643  ierr = VecDestroy(&F); CHKERRQ(ierr);
644  ierr = VecDestroy(&T); CHKERRQ(ierr);
645  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
646  if(monochromatic_wave) {
647 
648  if(add_incident_wave) {
649 
650  IncidentWave function_evaluator(wavenumber,wave_direction,power_of_incident_wave);
651 
652  ierr = solve_problem(m_field,"INCIDENT_WAVE","HELMHOLTZ_RERE_FE","rePRES","imPRES",
653  ADD_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
654 
655  }
656 
657  PetscBool save_postproc_mesh = PETSC_TRUE;
658  ierr = PetscOptionsGetBool(NULL,"-save_postproc_mesh",&save_postproc_mesh,NULL); CHKERRQ(ierr);
659  if(save_postproc_mesh) {
660 
662  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
663  ierr = post_proc.addFieldValuesPostProc("rePRES"); CHKERRQ(ierr);
664  ierr = post_proc.addFieldValuesPostProc("imPRES"); CHKERRQ(ierr);
665  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
666  /* Particle Velocity The velocity of the
667  movement of these material particles is called particle velocity */
668  ierr = post_proc.addFieldValuesGradientPostProc("rePRES","reVELOCITY"); CHKERRQ(ierr);
669  ierr = post_proc.addFieldValuesGradientPostProc("imPRES","imVELOCITY"); CHKERRQ(ierr);
670 
671 
672  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
673  ierr = post_proc.addFieldValuesPostProc("reEX"); CHKERRQ(ierr);
674  ierr = post_proc.addFieldValuesPostProc("imEX"); CHKERRQ(ierr);
675  }
676 
677  PetscBool reynolds_stress = PETSC_FALSE;
678  ierr = PetscOptionsGetBool(NULL,"-reynolds_stress",&reynolds_stress,NULL); CHKERRQ(ierr);
679 
680  if(reynolds_stress) {
681 
682  post_proc.getOpPtrVector().push_back(
683  new ReynoldsStress(
684  m_field,
685  post_proc.postProcMesh,
686  post_proc.mapGaussPts,
687  "imPRES",
688  post_proc.commonData));
689 
690  post_proc.getOpPtrVector().push_back(
691  new ReynoldsStress(
692  m_field,
693  post_proc.postProcMesh,
694  post_proc.mapGaussPts,
695  "rePRES",
696  post_proc.commonData));
697 
698  }
699 
700  ierr = m_field.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE",post_proc); CHKERRQ(ierr);
701  rval = post_proc.postProcMesh.write_file("fe_solution_mesh_post_proc.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
702 
703  }
704 
705  if(is_partitioned) {
706  rval = moab.write_file("fe_solution.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
707  } else {
708 
709  if(!pcomm->rank()) {
710  rval = moab.write_file("fe_solution.h5m"); CHKERR_PETSC(rval);
711  }
712 
713  }
714 
715  }
716 
717  ierr = PetscTime(&v2);CHKERRQ(ierr);
718  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
719  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f S CPU Time = %f S \n",pcomm->rank(),v2-v1,t2-t1);
720 
721  ierr = PetscFinalize(); CHKERRQ(ierr);
722 
723  return 0;
724 
725 }
const double T
struture grouping operators and data used for helmholtz problemsIn order to assemble matrices and rig...
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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 ...
CoreTmp< 0 > Core
Definition: Core.hpp:1129
Calculate the analytical solution of impinging wave on cylinder.
auto post_proc
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
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
Structure used to enforce analytical boundary conditions.
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
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Analytical Dirichlet boundary conditions.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
Calculate the analytical solution of impinging wave on cylinder.
char mesh_file_name[255]
Projection of edge entities with one mid-node on hierarchical basis.
Calculate the analytical solution of plane wave guide propagating in direction theta.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
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[]
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DEPRECATED MoFEMErrorCode partition_finite_elements(const std::string &name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=-1)
partition finite elements
constexpr int order
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
CHKERRQ(ierr)
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
DEPRECATED MoFEMErrorCode build_problem(const std::string &name, const bool square_matrix, int verb=-1)
build problem data structures
Post processing.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
DEPRECATED MoFEMErrorCode partition_problem(const std::string &name, int verb=-1)
partition problem dofs (collective)
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...
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
AnalyticalSolutionTypes
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
static char help[]
DEPRECATED MoFEMErrorCode partition_ghost_dofs(const std::string &name, int verb=-1)
determine ghost nodes
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual bool check_field(const std::string &name) const =0
check if field is in database
Core (interface) class.
Definition: Core.hpp:77
continuous field
Definition: definitions.h:177
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
DEPRECATED MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethodThis function set data about problem, adjacencies and other multi-indices in ...
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 74 of file fe_approximation_point_source.cpp.