v0.10.0
fe_approximation_obsolete.cpp
Go to the documentation of this file.
1 /** \file fe_approximation.cpp
2  \ingroup mofem_helmholtz_elem
3 
4 
5  Calculates finite element (Galerkin) approximation for incident wave problem.
6 
7  Note:
8 
9  In this implementation, first acoustic potential field is approximated on
10  boundary and then finite element problem is solved.
11 
12  For more rigorous convergence study, trace of best approximations on boundary
13  can be calculated and then finite element for domain and Neumann/mix boundary.
14  That will give exact pollution error.
15 
16  */
17 
18 /*
19  * This file is part of MoFEM.
20  * MoFEM is free software: you can redistribute it and/or modify it under
21  * the terms of the GNU Lesser General Public License as published by the
22  * Free Software Foundation, either version 3 of the License, or (at your
23  * option) any later version.
24  *
25  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
26  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
27  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
28  * License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
32  *
33  */
34 
35 #include <MoFEM.hpp>
36 
37 #include <boost/numeric/ublas/vector_proxy.hpp>
38 #include <boost/shared_ptr.hpp>
39 #include <boost/ptr_container/ptr_map.hpp>
40 
41 #include <DirichletBC.hpp>
42 #include <PotsProcOnRefMesh.hpp>
43 
45 #include <boost/iostreams/tee.hpp>
46 #include <boost/iostreams/stream.hpp>
47 #include <petsctime.h>
48 #include <fstream>
49 #include <iostream>
50 
51 #include <stdexcept>
52 #include <cmath>
53 #include <boost/math/special_functions.hpp>
54 #include <complex>
55 
56 #include <FieldApproximation.hpp>
57 
58 using namespace std;
59 using namespace boost::math;
60 
61 namespace bio = boost::iostreams;
62 using bio::tee_device;
63 using bio::stream;
64 
65 using namespace MoFEM;
66 
67 #include <AnalyticalSolutions.hpp>
68 #include <AnalyticalDirihlet.hpp>
69 
70 #include <HelmholtzElement.hpp>
71 
73 
74  double wAveNumber;
75  double pOwer;
76 
77  Range tRis;
78 
79 };
80 
81 static char help[] = "...\n\n";
82 //argc = argument counts, argv = argument vectors
83 int main(int argc, char *argv[]) {
84 
85  ErrorCode rval;
86  PetscErrorCode ierr;
87 
88  PetscInitialize(&argc,&argv,(char *)0,help);
89 
90  //Core mb_instance;
91  moab::Core mb_instance;
92  Interface& moab = mb_instance;
93  int rank;
94  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
95 
96  PetscBool flg = PETSC_TRUE;
97  char mesh_file_name[255];
98  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
99  if(flg != PETSC_TRUE) {
100  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
101  }
102 
103  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
104  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
105  PetscBool is_partitioned = PETSC_FALSE;
106  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
107  if(is_partitioned) {
108  //Read mesh to MOAB
109  const char *option;
110  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
111  rval = moab.load_file(mesh_file_name, 0, option); CHKERR_PETSC(rval);
112  } else {
113  const char *option;
114  option = "";
115  rval = moab.load_file(mesh_file_name, 0, option); CHKERR_PETSC(rval);
116  }
117 
118  // Create MoFEM (cephas) database
119  MoFEM::Core core(moab);
120  MoFEM::Interface& m_field = core;
121 
122  // Get start time for analyse
123  PetscLogDouble t1,t2;
124  PetscLogDouble v1,v2;
125  ierr = PetscTime(&v1); CHKERRQ(ierr);
126  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
127 
128  //set entitities bit level
129  BitRefLevel bit_level0;
130  bit_level0.set(0);
131  EntityHandle meshset_level0;
132  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERR_PETSC(rval);
133  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
134 
135  //Fields
136  ierr = m_field.add_field("rePRES",H1,1); CHKERRQ(ierr);
137  ierr = m_field.add_field("imPRES",H1,1); CHKERRQ(ierr);
138 
139  //meshset consisting all entities in mesh
140  EntityHandle root_set = moab.get_root_set();
141  //add entities to field
142  ierr = m_field.add_ents_to_field_by_TETs(root_set,"rePRES"); CHKERRQ(ierr);
143  ierr = m_field.add_ents_to_field_by_TETs(root_set,"imPRES"); CHKERRQ(ierr);
144 
145  //set app. order
146  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
147  PetscInt order;
148  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
149  if(flg != PETSC_TRUE) {
150  order = 2;
151  }
152 
153  ierr = m_field.set_field_order(root_set,MBTET,"rePRES",order); CHKERRQ(ierr);
154  ierr = m_field.set_field_order(root_set,MBTRI,"rePRES",order); CHKERRQ(ierr);
155  ierr = m_field.set_field_order(root_set,MBEDGE,"rePRES",order); CHKERRQ(ierr);
156  ierr = m_field.set_field_order(root_set,MBVERTEX,"rePRES",1); CHKERRQ(ierr);
157  ierr = m_field.set_field_order(root_set,MBTET,"imPRES",order); CHKERRQ(ierr);
158  ierr = m_field.set_field_order(root_set,MBTRI,"imPRES",order); CHKERRQ(ierr);
159  ierr = m_field.set_field_order(root_set,MBEDGE,"imPRES",order); CHKERRQ(ierr);
160  ierr = m_field.set_field_order(root_set,MBVERTEX,"imPRES",1); CHKERRQ(ierr);
161 
162  if(!m_field.check_field("MESH_NODE_POSITIONS")) {
163  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,3); CHKERRQ(ierr);
164  ierr = m_field.add_ents_to_field_by_TETs(root_set,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
165  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
166  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
167  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
168  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
169 
170  ierr = m_field.build_fields(); CHKERRQ(ierr);
171  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
172  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
173 
174  } else {
175 
176  ierr = m_field.build_fields(); CHKERRQ(ierr);
177 
178  }
179 
180  // Finite Elements
181 
182  HelmholtzElement helmholtz_element(m_field);
183  ierr = helmholtz_element.addHelmholtzElements("rePRES","imPRES"); CHKERRQ(ierr);
184 
185  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
186 
187  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","reEX"); CHKERRQ(ierr);
188  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","imEX"); CHKERRQ(ierr);
189 
190  }
191 
192  Range bc_dirichlet_tris,analytical_bc_tris;
193  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"ANALYTICAL_BC",it)) {
194  rval = moab.get_entities_by_type(it->get_meshset(),MBTRI,analytical_bc_tris,true); CHKERR_PETSC(rval);
195  }
196  bc_dirichlet_tris.merge(analytical_bc_tris);
197  AnalyticalDirihletBC analytical_bc_real(m_field);
198  AnalyticalDirihletBC analytical_bc_imag(m_field);
199  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",analytical_bc_tris); CHKERRQ(ierr);
200  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",analytical_bc_tris); CHKERRQ(ierr);
201 
202  // This is added for a case than on some surface, defined by the user a
203  // incident plane wave (or other wave) is scattered.
204  map<int,PlaneIncidentWaveSacttrerData> planeWaveScatterData;
205 
206  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"SOFT_INCIDENT_WAVE_BC",it)) {
207 
208  //get block attributes
209  vector<double> attributes;
210  ierr = it->get_Cubit_attributes(attributes); CHKERRQ(ierr);
211  if(attributes.size()<2) {
212  SETERRQ1(PETSC_COMM_SELF,1,"not enough block attributes to define SOFT_INCIDENT_WAVE_BC, attributes.size() = %d ",attributes.size());
213  }
214 
215  planeWaveScatterData[it->get_msId()].wAveNumber = attributes[0];
216  planeWaveScatterData[it->get_msId()].pOwer = attributes[1];
217  ierr = PetscOptionsGetScalar(NULL,"-wave_number",&(planeWaveScatterData[it->get_msId()].wAveNumber),NULL); CHKERRQ(ierr);
218 
219  rval = moab.get_entities_by_type(it->get_meshset(),MBTRI,planeWaveScatterData[it->get_msId()].tRis,true); CHKERR_PETSC(rval);
220  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",planeWaveScatterData[it->get_msId()].tRis); CHKERRQ(ierr);
221  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",planeWaveScatterData[it->get_msId()].tRis); CHKERRQ(ierr);
222  bc_dirichlet_tris.merge(planeWaveScatterData[it->get_msId()].tRis);
223 
224  }
225 
226  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
227  // Build adjacencies
228  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
229 
230  // Problem
231  ierr = m_field.add_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
232  ierr = m_field.add_problem("BCREAL_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for real field
233  ierr = m_field.add_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for imag field
234 
235  // Set refinement level for problem
236  ierr = m_field.modify_problem_ref_level_add_bit("ACOUSTIC_PROBLEM",bit_level0); CHKERRQ(ierr);
237  ierr = m_field.modify_problem_ref_level_add_bit("BCREAL_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
238  ierr = m_field.modify_problem_ref_level_add_bit("BCIMAG_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
239 
240  // Add elements to problems
241  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
242  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_IMIM_FE"); CHKERRQ(ierr);
243  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_REIM_FE"); CHKERRQ(ierr);
244  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_IMRE_FE"); CHKERRQ(ierr);
245 
246  ierr = m_field.modify_problem_add_finite_element("BCREAL_PROBLEM","BCREAL_FE"); CHKERRQ(ierr);
247  ierr = m_field.modify_problem_add_finite_element("BCIMAG_PROBLEM","BCIMAG_FE"); CHKERRQ(ierr);
248 
249  // Build problems
250 
251  // build porblems
252  if(is_partitioned) {
253  // if mesh is partitioned
254 
255  ierr = m_field.build_partitioned_problem("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
256  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
257 
258  ierr = m_field.build_partitioned_problem("BCREAL_PROBLEM",true); CHKERRQ(ierr);
259  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM",true); CHKERRQ(ierr);
260 
261  ierr = m_field.build_partitioned_problem("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
262  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
263 
264  } else {
265  // if not partitioned mesh is load to all processes
266 
267  ierr = m_field.build_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
268  ierr = m_field.partition_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
269  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
270 
271  ierr = m_field.build_problem("BCREAL_PROBLEM"); CHKERRQ(ierr);
272  ierr = m_field.partition_problem("BCREAL_PROBLEM"); CHKERRQ(ierr);
273  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM"); CHKERRQ(ierr);
274 
275  ierr = m_field.build_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr);
276  ierr = m_field.partition_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr);
277  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM"); CHKERRQ(ierr);
278 
279  }
280 
281  ierr = m_field.partition_ghost_dofs("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
282  ierr = m_field.partition_ghost_dofs("BCREAL_PROBLEM"); CHKERRQ(ierr);
283  ierr = m_field.partition_ghost_dofs("BCIMAG_PROBLEM"); CHKERRQ(ierr);
284 
285  // Get problem matrices and vectors
286 
287  Vec F; //Right hand side vector
288  ierr = m_field.VecCreateGhost("ACOUSTIC_PROBLEM",ROW,&F); CHKERRQ(ierr);
289  Vec T; //Solution vector
290  ierr = VecDuplicate(F,&T); CHKERRQ(ierr);
291  Mat A; //Left hand side matrix
292  ierr = m_field.MatCreateMPIAIJWithArrays("ACOUSTIC_PROBLEM",&A); CHKERRQ(ierr);
293 
294  // Solve for analytical Dirichlet bc dofs
295  ierr = analytical_bc_real.setProblem(m_field,"BCREAL_PROBLEM"); CHKERRQ(ierr);
296  ierr = analytical_bc_imag.setProblem(m_field,"BCIMAG_PROBLEM"); CHKERRQ(ierr);
297 
298  double angularfreq = 1;
299  double speed = 1;
300 
301  /// this works only for one block
302  int nb_of_blocks = 0;
303  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"MAT_HELMHOLTZ",it)) {
304 
305  // get block attributes
306  vector<double> attributes;
307  ierr = it->get_Cubit_attributes(attributes); CHKERRQ(ierr);
308  if(attributes.size()<2) {
309  SETERRQ1(PETSC_COMM_SELF,MOFEM_INVALID_DATA,
310  "not enough block attributes, expected 2 attributes ( angular freq., speed) , attributes.size() = %d ",attributes.size());
311  }
312  angularfreq = attributes[0];
313  speed = attributes[1];
314  nb_of_blocks++;
315 
316  }
317 
318  if(nb_of_blocks!=1) {
319  PetscPrintf(PETSC_COMM_SELF,"Warning: wave number is set to all blocks based on last evaluated block");
320  }
321  double wavenumber = angularfreq/speed;
322 
323  // set wave number from line command, that overwrite numbre form block set
324  ierr = PetscOptionsGetScalar(NULL,"-wave_number",&wavenumber,NULL); CHKERRQ(ierr);
325 
326  //wave direction unit vector=[x,y,z]^T
327  double waveDirection[3];
328  int nmax=3;
329  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-wave_direction",waveDirection,&nmax,&flg); CHKERRQ(ierr);
330 
331  VectorDouble wave_direction;
332  wave_direction.resize(3);
333  cblas_dcopy(3, &waveDirection[0], 1, &wave_direction(0), 1);
334  if(flg != PETSC_TRUE) {
335  wave_direction[0] = 1;
336  wave_direction[1] = 0;
337  wave_direction[2] = 0;
338  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -wave_direction [3*1 vector] default:X direction [1,0,0]");
339  }
340 
341 
342  PetscInt choise_value = 0;
343  // set type of analytical solution
344  ierr = PetscOptionsGetEList(NULL,"-analytical_solution_type",analytical_solution_types,6,&choise_value,NULL); CHKERRQ(ierr);
345 
346 
347  switch((AnalyticalSolutionTypes)choise_value) {
348 
350 
351  {
352  boost::shared_ptr<HardSphereScatterWave> function_evaluator = boost::shared_ptr<HardSphereScatterWave>(new HardSphereScatterWave(wavenumber));
353  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
354  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
355  }
356 
357  break;
358 
360 
361  {
362 
363  boost::shared_ptr<SoftSphereScatterWave> function_evaluator = boost::shared_ptr<SoftSphereScatterWave>(new SoftSphereScatterWave(wavenumber));
364  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
365  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
366 
367  }
368 
369  break;
370 
371  case PLANE_WAVE:
372 
373  {
374 
375  boost::shared_ptr<PlaneWave> function_evaluator = boost::shared_ptr<PlaneWave>(new PlaneWave(wavenumber,0.25*M_PI));
376  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
377  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
378 
379  }
380 
381  break;
382 
384 
385  {
386  boost::shared_ptr<HardCylinderScatterWave> function_evaluator = boost::shared_ptr<HardCylinderScatterWave>(new HardCylinderScatterWave(wavenumber));
387  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
388  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
389 
390  }
391 
392  break;
393 
395 
396  {
397  boost::shared_ptr<SoftCylinderScatterWave> function_evaluator = boost::shared_ptr<SoftCylinderScatterWave>(new SoftCylinderScatterWave(wavenumber));
398  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
399  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
400  }
401 
402  break;
403 
404  case INCIDENT_WAVE:
405 
406  {
407 
408  boost::shared_ptr<IncidentWave> function_evaluator = boost::shared_ptr<IncidentWave>(new IncidentWave(wavenumber,wave_direction,1));
409  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
410  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
411 
412  }
413 
414  break;
415 
416 
417  }
418 
419  // Analytical boundary conditions
420  AnalyticalDirihletBC::DirichletBC analytical_ditihlet_bc_real(m_field,"rePRES",A,T,F);
421  AnalyticalDirihletBC::DirichletBC analytical_ditihlet_bc_imag(m_field,"imPRES",A,T,F);
422 
423  {
424  map<int,PlaneIncidentWaveSacttrerData>::iterator mit = planeWaveScatterData.begin();
425  for(;mit!=planeWaveScatterData.end();mit++) {
426 
427  // note negative field, scatter field should cancel incident wave
428  boost::shared_ptr<IncidentWave> function_evaluator = boost::shared_ptr<IncidentWave>(new IncidentWave(wavenumber,-1));
429  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
430  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
431 
432  }
433  }
434 
435  ierr = analytical_bc_real.solveProblem(m_field,"BCREAL_PROBLEM","BCREAL_FE",analytical_ditihlet_bc_real,bc_dirichlet_tris); CHKERRQ(ierr);
436  ierr = analytical_bc_imag.solveProblem(m_field,"BCIMAG_PROBLEM","BCIMAG_FE",analytical_ditihlet_bc_imag,bc_dirichlet_tris); CHKERRQ(ierr);
437 
438  ierr = analytical_bc_real.destroyProblem(); CHKERRQ(ierr);
439  ierr = analytical_bc_imag.destroyProblem(); CHKERRQ(ierr);
440 
441  // Zero vectors
442  ierr = VecZeroEntries(T); CHKERRQ(ierr);
443  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
444  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
445  ierr = VecZeroEntries(F); CHKERRQ(ierr);
446  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
447  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
448  ierr = MatZeroEntries(A); CHKERRQ(ierr);
449 
450  // Assemble problem
451  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
452  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
453 
454  ierr = helmholtz_element.setOperators(A,F,"rePRES","imPRES"); CHKERRQ(ierr);
455  ierr = helmholtz_element.calculateA("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
456  ierr = helmholtz_element.calculateF("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
457 
458  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
459  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
460 
461  ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
462  ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
463  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
464  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
465  ierr = VecScale(F,-1); CHKERRQ(ierr);
466 
467  // Solve problem
468  KSP solver;
469  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
470  ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
471  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
472  ierr = KSPSetUp(solver); CHKERRQ(ierr);
473 
474  ierr = KSPSolve(solver,F,T); CHKERRQ(ierr);
475 
476  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
477  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
478 
479  //Save data on mesh
480  if(is_partitioned) {
481 
482  // no need for global communication
483  ierr = m_field.set_local_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
484 
485  } else {
486 
487  ierr = m_field.set_global_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
488 
489  }
490 
491  // Destroy the KSP solvers
492  ierr = MatDestroy(&A); CHKERRQ(ierr);
493  ierr = VecDestroy(&F); CHKERRQ(ierr);
494  ierr = VecDestroy(&T); CHKERRQ(ierr);
495  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
496 
497  PetscBool add_incident_wave = PETSC_TRUE;
498  ierr = PetscOptionsGetBool(NULL,"-add_incident_wave",&add_incident_wave,NULL); CHKERRQ(ierr);
499 
500  if(add_incident_wave) {
501 
502  // define problem
503  ierr = m_field.add_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
504  // set finite elements for problem
505  ierr = m_field.modify_problem_add_finite_element("INCIDENT_WAVE","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
506  // set refinment level for problem
507  ierr = m_field.modify_problem_ref_level_add_bit("INCIDENT_WAVE",bit_level0); CHKERRQ(ierr);
508 
509  // build porblems
510  if(is_partitioned) {
511  // if mesh is partitioned
512  ierr = m_field.build_partitioned_problem("INCIDENT_WAVE",true); CHKERRQ(ierr);
513  ierr = m_field.partition_finite_elements("INCIDENT_WAVE",true); CHKERRQ(ierr);
514  } else {
515  // if not partitioned mesh is load to all processes
516  ierr = m_field.build_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
517  ierr = m_field.partition_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
518  ierr = m_field.partition_finite_elements("INCIDENT_WAVE"); CHKERRQ(ierr);
519  }
520  ierr = m_field.partition_ghost_dofs("INCIDENT_WAVE"); CHKERRQ(ierr);
521 
522  IncidentWave function_evaluator(wavenumber,wave_direction,1);
523  ierr = solve_problem(m_field,"INCIDENT_WAVE","HELMHOLTZ_RERE_FE","rePRES","imPRES",
524  ADD_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
525 
526  }
527 
528  if(is_partitioned) {
529  rval = moab.write_file("fe_solution.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
530  } else {
531  if(!pcomm->rank()) {
532  rval = moab.write_file("fe_solution.h5m"); CHKERR_PETSC(rval);
533  }
534  }
535 
536  PetscBool save_postproc_mesh = PETSC_TRUE;
537  ierr = PetscOptionsGetBool(NULL,"-save_postproc_mesh",&save_postproc_mesh,NULL); CHKERRQ(ierr);
538  if(save_postproc_mesh) {
539 
541  ierr = post_proc.generateRefereneElemenMesh(); CHKERRQ(ierr);
542  ierr = post_proc.addFieldValuesPostProc("rePRES"); CHKERRQ(ierr);
543  ierr = post_proc.addFieldValuesPostProc("imPRES"); CHKERRQ(ierr);
544 
545  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
546  ierr = post_proc.addFieldValuesPostProc("reEX"); CHKERRQ(ierr);
547  ierr = post_proc.addFieldValuesPostProc("imEX"); CHKERRQ(ierr);
548  }
549 
550  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
551  ierr = m_field.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE",post_proc); CHKERRQ(ierr);
552  rval = post_proc.postProcMesh.write_file("fe_solution_mesh_post_proc.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
553 
554  }
555 
556  ierr = PetscTime(&v2);CHKERRQ(ierr);
557  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
558  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f S CPU Time = %f S \n",pcomm->rank(),v2-v1,t2-t1);
559 
560  ierr = PetscFinalize(); CHKERRQ(ierr);
561 
562  return 0;
563 
564 }
565 
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 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
PetscErrorCode solveProblem(MoFEM::Interface &m_field, string problem, string fe, DirichletBC &bc)
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.
Analytical Dirichlet boundary conditions.
auto post_proc
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
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
PetscErrorCode setProblem(MoFEM::Interface &m_field, string problem)
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
PetscErrorCode setApproxOps(MoFEM::Interface &m_field, string re_field_name, double(*fun)(double x, double y, double z, bool use_real), bool useReal, string nodals_positions="MESH_NODE_POSITIONS")
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.
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.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
char mesh_file_name[255]
PetscErrorCode initializeProblem(MoFEM::Interface &m_field, string problem, string fe, string re_field_name, string nodals_positions="MESH_NODE_POSITIONS")
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 ...
Calculate the analytical solution of plane wave guide propagating in direction theta.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
static char help[]
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)
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...
int main(int argc, char *argv[])
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
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...
DEPRECATED MoFEMErrorCode add_ents_to_field_by_TETs(const EntityHandle meshset, const std::string &name, int verb=-1)
set field entities from adjacencies of tetrahedronThe lower dimension entities are added depending on...
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)