v0.10.0
fe_approximation_point_source.cpp
Go to the documentation of this file.
1 /** \file fe_approximation.cpp
2 
3  Calculates finite element (Galerkin) approximation for incident wave problem.
4 
5  Note:
6 
7  In this implementation, first acoustic potential field is approximated on
8  boundary and then finite element problem is solved. For more rigorous
9  convergence study, trace of best approximations on boundary can be calculated
10  and then finite element for domain and Neumann/mix boundary. That will give
11  exact pollution error.
12 
13  */
14 
15 /*
16  * This file is part of MoFEM.
17  * MoFEM is free software: you can redistribute it and/or modify it under
18  * the terms of the GNU Lesser General Public License as published by the
19  * Free Software Foundation, either version 3 of the License, or (at your
20  * option) any later version.
21  *
22  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
23  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
24  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
25  * License for more details.
26  *
27  * You should have received a copy of the GNU Lesser General Public
28  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
29  *
30  */
31 
32 #include <MoFEM.hpp>
33 using namespace MoFEM;
34 
35 #include <boost/numeric/ublas/vector_proxy.hpp>
36 #include <boost/numeric/ublas/symmetric.hpp>
37 #include <boost/shared_ptr.hpp>
38 #include <boost/ptr_container/ptr_map.hpp>
39 
40 #include <DirichletBC.hpp>
41 #include <PostProcOnRefMesh.hpp>
42 #include <ReynoldsStress.hpp>
43 
45 #include <petsctime.h>
46 #include <fstream>
47 #include <iostream>
48 
49 #include <stdexcept>
50 #include <cmath>
51 #include <boost/math/special_functions.hpp>
52 #include <complex>
53 
54 #include <FieldApproximation.hpp>
55 
56 using namespace std;
57 using namespace boost::math;
58 
59 #include <boost/shared_array.hpp>
60 #include <kiss_fft.h>
61 #include <kiss_fft.c>
62 
63 #include <AnalyticalSolutions.hpp>
64 #include <AnalyticalDirichlet.hpp>
65 #include <HelmholtzElement.hpp>
66 #include <TimeSeries.hpp>
67 
69 
70  Range tRis;
71 
72 };
73 
74 static char help[] = "...\n\n";
75 
76 int main(int argc, char *argv[]) {
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...
Operators and data structures for wave propagation analyze (Galerkin Element)
PetscErrorCode generateReferenceElementMesh()
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.
Post-process fields on refined mesh.
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.
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
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.
MoFEMErrorCode solveProblem(MoFEM::Interface &m_field, string problem, string fe, DirichletBC &bc, Range &tris)
solve boundary problem
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
char mesh_file_name[255]
Projection of edge entities with one mid-node on hierarchical basis.
Read impulse and apply FFTFirst signal is read from text file. Currently is assumed that this signal ...
PetscErrorCode readData(const char *str, map< double, double > &series)
Read signal from text file.
Definition: TimeSeries.hpp:105
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 forwardSpaceDft()
Apply Forward FFT and compose wave.
Definition: TimeSeries.hpp:430
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
PetscErrorCode solveForwardDFT(KSP solver, Mat A, Vec F, Vec T)
Solve problem for each wave number.
Definition: TimeSeries.hpp:590
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
PetscErrorCode createPressureSeries(Vec T)
Create vectors for each wave lengths.
Definition: TimeSeries.hpp:526
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
PetscErrorCode pressureInverseDft()
Aplly Inverse FFT to pressure degrees of freedom.
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)
MoFEMErrorCode setApproxOps(MoFEM::Interface &m_field, const std::string field_name, boost::shared_ptr< FUNEVAL > function_evaluator, const int field_number=0, const string nodals_positions="MESH_NODE_POSITIONS")
Set operators used to calculate the rhs vector and the lhs matrix.
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
PetscErrorCode destroyPressureSeries()
Destroy pressure series.
Definition: TimeSeries.hpp:550
PetscErrorCode pressureForwardDft()
Aplly Forward FFT to pressure degrees of freedom.
PetscErrorCode saveResults()
Save data on mesh.
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
int main(int argc, char *argv[])
Core (interface) class.
Definition: Core.hpp:77
Operators and data structures for the calculation of momentum flux and Reynolds Stress of the complex...
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)
MoFEMErrorCode destroyProblem()
Destroy problem.
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)