v0.10.0
spherical_acoustic.cpp
Go to the documentation of this file.
1 /* Copyright (C) 2013, Lukasz Kaczmarczyk (likask AT wp.pl)
2  * --------------------------------------------------------------
3  * FIXME: DESCRIPTION
4  */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *PhD student Thomas Xuan Meng/Users/xuanmeng/Documents/mofem-cephas/mofem_v0.2/users_modules/helmholtz/README
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
19  * abs = sqrt(rePRES^2+imPRES^2) */
20 
21 
22 #include <MoFEM.hpp>
23 
24 
25 #include <DirichletBC.hpp>
26 #include <PotsProcOnRefMesh.hpp>
27 #include <HelmholtzElement.hpp>
28 
31 #include <boost/iostreams/tee.hpp>
32 #include <boost/iostreams/stream.hpp>
33 #include <petsctime.h>
34 #include <fstream>
35 #include <iostream>
36 
37 #include <stdexcept>
38 #include <cmath>
39 #include <boost/math/special_functions.hpp>
40 #include <complex>
41 
42 
43 using namespace std;
44 using namespace boost::math;
45 
46 namespace bio = boost::iostreams;
47 using bio::tee_device;
48 using bio::stream;
49 
50 using namespace MoFEM;
51 
52 static char help[] = "...\n\n";
53 //argc = argument counts, argv = argument vectors
54 int main(int argc, char *argv[]) {
55 
56 
57 
58 
59  PetscInitialize(&argc,&argv,(char *)0,help);
60 
61  //Core mb_instance;
62  moab::Core mb_instance;
63  Interface& moab = mb_instance;
64  int rank;
65  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
66 
67  bool useImpedance;
68  PetscBool flg = PETSC_TRUE;
69  char mesh_file_name[255];
70  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
71  if(flg != PETSC_TRUE) {
72  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
73  }
74 
75  char impedance[255];
76  ierr = PetscOptionsGetString(PETSC_NULL,"-use_impedance",impedance,255,&flg); CHKERRQ(ierr);
77  if(flg != PETSC_TRUE) {
78  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -use_impedance (true of false needed)");
79  }
80  if (strcmp ("true",impedance ) == 0) {useImpedance = true;}
81  else if(strcmp ("false",impedance ) == 0) {useImpedance = false;}
82 
83  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
84  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
85 
86  const char *option;
87  option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
88  BARRIER_RANK_START(pcomm)
89  rval = moab.load_file(mesh_file_name, 0, option); CHKERR_PETSC(rval);
90  BARRIER_RANK_END(pcomm)
91 
92  //Create MoFEM (cephas) database
93  MoFEM::Core core(moab);
94  MoFEM::Interface& mField = core;
95 
96  //count the comsumption of time by single run
97  PetscLogDouble t1,t2;
98  PetscLogDouble v1,v2;
99  ierr = PetscTime(&v1); CHKERRQ(ierr);
100  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
101 
102  //set entitities bit level
103  BitRefLevel bit_level0;
104  bit_level0.set(0);
105  EntityHandle meshset_level0;
106  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERR_PETSC(rval);
107  ierr = mField.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
108 
109  //Fields
110  ierr = mField.add_field("rePRES",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr); //field order distinguish the scalar field and vector field.
111  ierr = mField.add_field("imPRES",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
112 
113  //Problem
114  ierr = mField.add_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
115  ierr = mField.add_problem("BC1_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for real field
116  ierr = mField.add_problem("BC2_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for imag field
117 
118  //set refinment level for problem
119  ierr = mField.modify_problem_ref_level_add_bit("ACOUSTIC_PROBLEM",bit_level0); CHKERRQ(ierr);
120  ierr = mField.modify_problem_ref_level_add_bit("BC1_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
121  ierr = mField.modify_problem_ref_level_add_bit("BC2_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
122 
123  //meshset consisting all entities in mesh
124  EntityHandle root_set = moab.get_root_set();
125  //add entities to field
126  ierr = mField.add_ents_to_field_by_type(root_set,MBTET,"rePRES"); CHKERRQ(ierr);
127  ierr = mField.add_ents_to_field_by_type(root_set,MBTET,"imPRES"); CHKERRQ(ierr);
128 
129  //set app. order
130  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
131  PetscInt order;
132  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
133  if(flg != PETSC_TRUE) {
134  order = 2;
135  }
136 
137  ierr = mField.set_field_order(root_set,MBTET,"rePRES",order); CHKERRQ(ierr);
138  ierr = mField.set_field_order(root_set,MBTRI,"rePRES",order); CHKERRQ(ierr);
139  ierr = mField.set_field_order(root_set,MBEDGE,"rePRES",order); CHKERRQ(ierr);
140  ierr = mField.set_field_order(root_set,MBVERTEX,"rePRES",1); CHKERRQ(ierr);
141 
142  if(!mField.check_field("MESH_NODE_POSITIONS")) {
143  ierr = mField.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
144  ierr = mField.add_ents_to_field_by_type(root_set,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
145  ierr = mField.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
146  ierr = mField.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
147  ierr = mField.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
148  ierr = mField.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
149  }
150  ierr = mField.set_field_order(root_set,MBTET,"imPRES",order); CHKERRQ(ierr);
151  ierr = mField.set_field_order(root_set,MBTRI,"imPRES",order); CHKERRQ(ierr);
152  ierr = mField.set_field_order(root_set,MBEDGE,"imPRES",order); CHKERRQ(ierr);
153  ierr = mField.set_field_order(root_set,MBVERTEX,"imPRES",1); CHKERRQ(ierr);
154 
155  HelmholtzElement helmholtz_elements(mField); //Create the HelmholtzElement class in the header-file
156 
157  ierr = helmholtz_elements.addHelmholtzElements("ACOUSTIC_PROBLEM","rePRES","imPRES"); CHKERRQ(ierr);
158  ierr = helmholtz_elements.addHelmholtzFluxElement("ACOUSTIC_PROBLEM","rePRES","imPRES"); CHKERRQ(ierr);
159  ierr = helmholtz_elements.addHelmholtzIncidentElement("ACOUSTIC_PROBLEM","rePRES","imPRES"); CHKERRQ(ierr);
160  if(useImpedance) {
161  ierr = helmholtz_elements.addHelmholtzImpedanceElement("ACOUSTIC_PROBLEM","rePRES","imPRES"); CHKERRQ(ierr);
162  }
163 
164 
165  /*Set up the analytical Dirichlet BC */
166  Range bc_tris;
167  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(mField,"ANALYTICAL_BC",it)) {
168  rval = moab.get_entities_by_type(it->get_meshset(),MBTRI,bc_tris,true); CHKERR_PETSC(rval);
169  }
170 
171  AnalyticalDirihletBC analytical_bc1(mField,bc_tris);
172  AnalyticalDirihletBC analytical_bc2(mField,bc_tris);
173  ierr = analytical_bc1.initializeProblem(mField,"BC1_PROBLEM","BC1_FE","rePRES"); CHKERRQ(ierr);
174  ierr = analytical_bc2.initializeProblem(mField,"BC2_PROBLEM","BC2_FE","imPRES"); CHKERRQ(ierr);
175  ///*End of Dirichlet set up */
176 
177  /*** add exact solution data in finite element */
178  if(mField.check_field("reEX") && mField.check_field("imEX")) {
179  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FE","reEX"); CHKERRQ(ierr);
180  ierr = mField.modify_finite_element_add_field_data("HELMHOLTZ_FE","imEX"); CHKERRQ(ierr);
181  }
182 
183 
184  /****/
185  //build database
186  //build field
187  ierr = mField.build_fields(); CHKERRQ(ierr);
188  //build finite elemnts
189  ierr = mField.build_finite_elements(); CHKERRQ(ierr);
190  //build adjacencies
191  ierr = mField.build_adjacencies(bit_level0); CHKERRQ(ierr);
192  //build problem
193  ierr = mField.build_problems(); CHKERRQ(ierr);
194 
195  Projection10NodeCoordsOnField ent_method_material(mField,"MESH_NODE_POSITIONS");
196  ierr = mField.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
197 
198  /****/
199  //mesh partitioning
200  //partition
201  ierr = mField.partition_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
202  ierr = mField.partition_finite_elements("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
203  //what are ghost nodes, see Petsc Manual
204  ierr = mField.partition_ghost_dofs("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
205 
206  //////mesh partitinoning for analytical Dirichlet
207  ierr = mField.simple_partition_problem("BC1_PROBLEM"); CHKERRQ(ierr);
208  ierr = mField.partition_finite_elements("BC1_PROBLEM"); CHKERRQ(ierr);
209  ierr = mField.simple_partition_problem("BC2_PROBLEM"); CHKERRQ(ierr);
210  ierr = mField.partition_finite_elements("BC2_PROBLEM"); CHKERRQ(ierr);
211  ////what are ghost nodes, see Petsc Manual
212  ierr = mField.partition_ghost_dofs("BC1_PROBLEM"); CHKERRQ(ierr);
213  ierr = mField.partition_ghost_dofs("BC2_PROBLEM"); CHKERRQ(ierr);
214 
215  Vec F; //Right hand side vector
216  ierr = mField.VecCreateGhost("ACOUSTIC_PROBLEM",ROW,&F); CHKERRQ(ierr);
217  Vec T; //Solution vector
218  ierr = VecDuplicate(F,&T); CHKERRQ(ierr);
219  Mat A; //Left hand side matrix
220  ierr = mField.MatCreateMPIAIJWithArrays("ACOUSTIC_PROBLEM",&A); CHKERRQ(ierr);
221 
222  //bool useScalar = true;
223  //ierr = helmholtz_elements.setHelmholtzFiniteElementRhs_FOperators("rePRES","rePRES",F,useScalar); CHKERRQ(ierr); //The analytical F source vector
224  ierr = helmholtz_elements.setHelmholtzFiniteElementRhsOperators("rePRES","imPRES",F,useImpedance); CHKERRQ(ierr); //the Rhs of Dirichlet BC
225  ierr = helmholtz_elements.setHelmholtzFiniteElementLhsOperators("rePRES","imPRES",(A)); CHKERRQ(ierr);//Stiffness Matrix
226  //ierr = helmholtz_elements.setHelmholtzFluxFiniteElementRhsOperators("rePRES","rePRES",F); CHKERRQ(ierr); //real Neumann BC
227  //ierr = helmholtz_elements.setHelmholtzFluxFiniteElementRhsOperators("imPRES","imPRES",F); CHKERRQ(ierr); //Imag Neumann BC
228  ierr = helmholtz_elements.setHelmholtzIncidentWaveFiniteElementRhsOperators("rePRES","imPRES",F); CHKERRQ(ierr); // Incident wave flux.
229  if(useImpedance) {
230  //The boundary Impedance BC.
231  ierr = helmholtz_elements.setHelmholtzImpedanceFiniteElementLhsOperators("rePRES","imPRES",(A)); CHKERRQ(ierr);
232  }
233 
234 
235  ierr = VecZeroEntries(T); CHKERRQ(ierr);
236  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
237  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
238  ierr = VecZeroEntries(F); CHKERRQ(ierr);
239  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
240  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
241  ierr = MatZeroEntries(A); CHKERRQ(ierr);
242 
243  ////analytical dirichlet bc
244  AnalyticalDirihletBC::DirichletBC analytical_ditihlet_bc1(mField,"rePRES",A,T,F);
245  AnalyticalDirihletBC::DirichletBC analytical_ditihlet_bc2(mField,"imPRES",A,T,F);
246 
247  ////solve for analytical dirichlet bc dofs
248  ierr = analytical_bc1.setProblem(mField,"BC1_PROBLEM"); CHKERRQ(ierr);
249  ierr = analytical_bc2.setProblem(mField,"BC2_PROBLEM"); CHKERRQ(ierr);
250 
251  static double aNgularfreq;
252  static double sPeed; //Without static. got error:use of local variable with automatic storage from containing function
253 
255  {
256  cout << endl << *it << endl;
257 
258 
259  if(it->get_Cubit_name().compare(0,13,"MAT_HELMHOLTZ") == 0) {
260 
261  //get block attributes
262  vector<double> attributes;
263  ierr = it->get_Cubit_attributes(attributes); CHKERRQ(ierr);
264  if(attributes.size()<2) {
265  SETERRQ1(PETSC_COMM_SELF,1,"not enough block attributes to deffine fluid pressure element, attributes.size() = %d ",attributes.size());
266  }
267 
268  aNgularfreq = attributes[0];
269  sPeed = attributes[1];
270 
271  }
272  }
273 
274 
275 
276  /* this function compute the scattered field of helmholtz operator */
277  struct AnaliticalFunction {
278  static double fUN(double x,double y,double z,bool use_real) {
279 
280 
281  const double pi = atan( 1.0 ) * 4.0;
282  double R = sqrt(pow(x,2.0)+pow(y,2.0)+pow(z,2.0)); //radius
283  double theta = atan2(y,x)+2*pi; //the arctan of radians (y/x)
284 
285  const double wAvenumber = aNgularfreq/sPeed;
286 
287  const double k = wAvenumber; //Wave number
288  const double a = 0.5; //radius of the sphere,wait to modify by user
289  const double const1 = k * a;
290  double const2 = k * R;
291  const complex< double > i( 0.0, 1.0 );
292 
293  // magnitude of incident wave
294  const double phi_incident_mag = 1.0;
295 
296  const double tol = 1.0e-10;
297  double max = 0.0;
298  double min = 999999.0;
299 
300  complex< double > result = 0.0;
301  //complex< double > prev_result;
302 
303  //double error = 100.0;
304  //unsigned int n = 0; //initialized the infinite series loop
305 
306  //while( error > tol ) //finding the acoustic potential in one single point.
307  //{
308  // double jn_der = n / const1 * sph_bessel( n, const1 ) - sph_bessel( n + 1, const1 ); //The derivative of bessel function
309 
310  //
311  // complex< double > hn_der = n / const1 * sph_hankel_1( n, const1 ) - sph_hankel_1( n + 1, const1 );
312  // //complex< double > hn_der = 0.5 * ( sph_hankel_1( n - 1, const1 ) -
313  // //( sph_hankel_1( n, const1 ) + const1 * sph_hankel_1( n + 1, const1 ) ) / const1 );
314  // double Pn = legendre_p( n, cos( theta ) );
315 
316  // complex< double >hn = sph_hankel_1( n, const2 ); //S Hankel first kind function
317  //
318  // prev_result = result;
319  // result -= pow( i, n ) * ( 2.0 * n + 1.0 ) * jn_der / hn_der * Pn * hn;
320  // error = abs( abs( result ) - abs( prev_result ) );
321  // ++n;
322  //}
323 
324 
325  const complex< double > inc_field = -exp( i * k * R * cos( theta ) ); //incident wave
326  //const complex< double > total_field = inc_field + result;
327 
328  result = inc_field;
329 
330  if(use_real) {
331  return std::real(result);
332  } else {
333  return std::imag(result);
334  }
335 
336 
337  }
338 
339  };
340 
341 
342 
343  ierr = analytical_bc1.setApproxOps(mField,"rePRES",AnaliticalFunction::fUN,true); CHKERRQ(ierr); //Triangles
344  ierr = analytical_bc2.setApproxOps(mField,"imPRES",AnaliticalFunction::fUN,false); CHKERRQ(ierr);
345 
346  ierr = analytical_bc1.solveProblem(mField,"BC1_PROBLEM","BC1_FE",analytical_ditihlet_bc1); CHKERRQ(ierr);
347  ierr = analytical_bc2.solveProblem(mField,"BC2_PROBLEM","BC2_FE",analytical_ditihlet_bc2); CHKERRQ(ierr);
348 
349 
350  ierr = analytical_bc1.destroyProblem(); CHKERRQ(ierr);
351  ierr = analytical_bc2.destroyProblem(); CHKERRQ(ierr);
352 
353  /* preproc */
354  /* Preprocess the analytical Dirichlet BC */
355  ierr = mField.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc1); CHKERRQ(ierr);
356  ierr = mField.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc2); CHKERRQ(ierr);
357 
358 
359 
360  ierr = mField.set_global_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
361 
362  ierr = mField.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_FE",helmholtz_elements.getLoopFeRhs()); CHKERRQ(ierr);
363  ierr = mField.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_FE",helmholtz_elements.getLoopFeLhs()); CHKERRQ(ierr);
364  //ierr = mField.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_FLUX_FE",helmholtz_elements.getLoopFeFlux()); CHKERRQ(ierr); //scalar flux
365  ierr = mField.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_INCIDENT_FE",helmholtz_elements.getLoopfeIncidentWave()); CHKERRQ(ierr); //Incident wave flux
366 
367 
368  if(useImpedance) {
369  ierr = mField.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_IMPEDANCE_FE",helmholtz_elements.getLoopFeImpedanceLhs()); CHKERRQ(ierr);
370  }
371  /*above terms related to operators in HelmholtzElement.hpp */
372 
373 
374  /* postproc */
375  /* Postprocess the Analytical Dirichlet BC */
376 
377  ierr = mField.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc1); CHKERRQ(ierr);
378  ierr = mField.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc2); CHKERRQ(ierr);
379 
380  ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
381  ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
382  ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
383  ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
384  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
385  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
386  ierr = VecScale(F,-1); CHKERRQ(ierr);
387 
388  //Solver
389  KSP solver;
390  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
391  ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
392  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
393  ierr = KSPSetUp(solver); CHKERRQ(ierr);
394 
395  ierr = KSPSolve(solver,F,T); CHKERRQ(ierr);
396  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
397  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
398 
399  ierr = mField.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc1); CHKERRQ(ierr);
400  ierr = mField.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc2); CHKERRQ(ierr);
401 
402  //Save data on mesh
403  ierr = mField.set_global_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
404 
405  ////Wait to putput the data in format
406  //PetscViewer viewer;
407  //ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"Acoustic_Impining_Sphere.txt",&viewer); CHKERRQ(ierr);
408  //VecView(T,viewer);
409  //ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
410 
411 
412 
413  //if(pcomm->rank()==0) {
414  rval = moab.write_file("impinging_finite_element.h5m"); CHKERR_PETSC(rval);
415  //}
416 
417  //destroy the KSP solvers
418  ierr = MatDestroy(&A); CHKERRQ(ierr);
419  ierr = VecDestroy(&F); CHKERRQ(ierr);
420  ierr = VecDestroy(&T); CHKERRQ(ierr);
421  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
422 
423 
424  PostProcVolumeOnRefinedMesh post_proc1(mField);
425 
426  ierr = post_proc1.generateRefereneElemenMesh(); CHKERRQ(ierr);
427  ierr = post_proc1.addFieldValuesPostProc("rePRES"); CHKERRQ(ierr);
428  ierr = post_proc1.addFieldValuesPostProc("imPRES"); CHKERRQ(ierr);
429 
430  if(mField.check_field("reEX") && mField.check_field("imEX")) {
431 
432  ierr = post_proc1.addFieldValuesPostProc("reEX"); CHKERRQ(ierr);
433  ierr = post_proc1.addFieldValuesPostProc("imEX"); CHKERRQ(ierr);
434 
435  ierr = post_proc1.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
436  ierr = mField.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_FE",post_proc1); CHKERRQ(ierr);
437 
438  rval = post_proc1.postProcMesh.write_file("four_fields_finite_element.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
439 
440  ////output the results from Docker
441  //char command1[] = "mbconvert ./four_fields_finite_element.h5m ./four_fields_finite_element.vtk && cp ./four_fields_finite_element.vtk ../../../../../mnt/home/Desktop/U_pan/helmholtz_results/";
442  //int todo1 = system( command1 );
443 
444  } else {
445 
446  ierr = post_proc1.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
447  ierr = mField.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_FE",post_proc1); CHKERRQ(ierr);
448  rval = post_proc1.postProcMesh.write_file("acoustic_finite_element_out.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
449 
450  }
451 
452  /** get the time interval **/
453  ierr = PetscTime(&v2);CHKERRQ(ierr);
454  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
455  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f S CPU Time = %f S \n",pcomm->rank(),v2-v1,t2-t1);
456 
457 
458 
459 
460 
461  //wait to calculate the the magnitude of acoustic potential-abs(phi) using real and imag components.
462  //ierr = VecView(T,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
463  //std::cin >> wait;
464 
465  ////Calculate the L^2 Norm manually below
466  //ierr = VecPow(C,2.0); CHKERRQ(ierr);
467  //ierr = VecAbs(C); CHKERRQ(ierr);
468  //Vec H;
469  //ierr = VecDuplicate(T,&H); CHKERRQ(ierr);
470  //ierr = VecCopy(T,H); CHKERRQ(ierr);
471  ////ierr = VecPow(H,2.0); CHKERRQ(ierr);
472  //ierr = VecAbs(H); CHKERRQ(ierr);
473  ////std::cout << "\n \n \n \n \n \n \n " << std::endl;
474  ////ierr = VecView(H,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
475  ////std::cin >> wait;
476  //ierr = VecAXPBY(C,-1.0,1.0,H); CHKERRQ(ierr);
477  //PetscReal l2norm,pointwisenorm,sum;
478  //PetscScalar pow = 2;
479  ////Compute the Global L^2 norm error and point wise error
480  //Vec G;
481  //ierr = VecDuplicate(C,&G); CHKERRQ(ierr);
482  //ierr = VecCopy(C,G); CHKERRQ(ierr);
483  //ierr = VecPow(G,pow); CHKERRQ(ierr);
484  //ierr = VecSum(G,&sum); CHKERRQ(ierr);
485  //l2norm = sqrt(sum);
486  //ierr = VecDestroy(&G); CHKERRQ(ierr);
487 
488  ////ierr = VecNorm(C,NORM_FROBENIUS,&l2norm);;
489  ////ierr = VecNorm(C,NORM_MAX,&pointwisenorm);
490  //ierr = VecMax(C,NULL,&pointwisenorm);
491  ////std::cout << "\n The Global L2 Norm of error for this problem is : --\n" << l2norm << std::endl;
492  //std::cout << "\n The Global Pointwise Norm of error for this problem is : --\n" << pointwisenorm << std::endl;
493  //std::cin >> wait;
494  ////ierr = VecSqrtAbs(C); CHKERRQ(ierr);
495  //ierr = VecView(C,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
496  //std::cin >> wait;
497  //std::cout << "\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n " << std::endl;
498  //give the right value of Im(Scalar Solution)
499 
500 
501  ierr = PetscFinalize(); CHKERRQ(ierr);
502 
503  return 0;
504 
505 }
506 
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 addHelmholtzFluxElement(const string problem_name, const string re_field_name, const string im_field_name)
add helmholtz flux element \infroup mofem_helmholtz_elem
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MyVolumeFE & getLoopFeRhs()
get rhs volume element
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
Analytical Dirichlet boundary conditions.
PetscErrorCode setHelmholtzFiniteElementLhsOperators(string re_field_name, string im_field_name, Mat A)
this fucntion is used in case of stationary Helmholtz problem for lhs stiffness term \infroup mofem_h...
#define BARRIER_RANK_END(PCMB)
Definition: definitions.h:351
static char help[]
moab::Interface & postProcMesh
virtual DEPRECATED MoFEMErrorCode build_problems(int verb=-1)=0
build problem data structures
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
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
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")
PetscErrorCode setHelmholtzFiniteElementRhsOperators(string re_field_name, string im_field_name, Vec &F, bool useImpedance)
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.
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.
PetscErrorCode addHelmholtzImpedanceElement(const string problem_name, const string re_field_name, const string im_field_name)
add Impedance element \infroup mofem_helmholtz_elem
int main(int argc, char *argv[])
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
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
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
constexpr int order
double tol
static Index< 'i', 3 > i
#define BARRIER_RANK_START(PCMB)
Definition: definitions.h:328
PetscErrorCode setHelmholtzImpedanceFiniteElementLhsOperators(string re_field_name, string im_field_name, Mat A, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
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
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.
MyVolumeFE & getLoopFeLhs()
get lhs volume element
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
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
PetscErrorCode setHelmholtzIncidentWaveFiniteElementRhsOperators(string re_field_name, string im_field_name, Vec &F, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
this function is used in case of statonary for helmholtz incident wave flux terms
const double R
#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.
static Index< 'k', 3 > k
PetscErrorCode addHelmholtzIncidentElement(const string problem_name, const string re_field_name, const string im_field_name)
add helmholtz flux element \infroup mofem_helmholtz_elem
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
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...
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...
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)