v0.10.0
fe_field_split.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 
41 #include <DirichletBC.hpp>
42 #include <PostProcOnRefMesh.hpp>
43 
45 #include <petsctime.h>
46 #include <fstream>
47 #include <iostream>
48 #include <valarray>
49 
50 #include <stdexcept>
51 #include <cmath>
52 #include <boost/math/special_functions.hpp>
53 #include <complex>
54 
55 #include <FieldApproximation.hpp>
56 
57 using namespace std;
58 using namespace boost::math;
59 
60 #include <boost/shared_array.hpp>
61 #include <kiss_fft.h>
62 #include <kiss_fft.c>
63 
65 #include <AnalyticalSolutions.hpp>
66 #include <HelmholtzElement.hpp>
67 #include <ReynoldsStress.hpp>
68 #include <SecondOrderStress.hpp>
69 #include <TimeSeriesFieldSplit.hpp>
70 
71 // struct PlaneIncidentWaveSacttrerData {
72 //
73 // Range tRis;
74 //
75 // };
76 static int debug = 0;
77 static char help[] = "...\n\n";
78 
79 int main(int argc, char *argv[]) {
80 
81  ErrorCode rval;
82  PetscErrorCode ierr;
83 
84  PetscInitialize(&argc,&argv,(char *)0,help);
85 
86  //Core mb_instance;
87  moab::Core mb_instance;
88  moab::Interface& moab = mb_instance;
89  int rank;
90  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
91 
92  PetscBool flg = PETSC_TRUE;
93  char mesh_file_name[255];
94  ierr = PetscOptionsGetString(PETSC_NULL,PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
95  if(flg != PETSC_TRUE) {
96  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"*** ERROR -my_file (MESH FILE NEEDED)");
97  }
98 
99  // Create moab parallel communicator
100  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
101  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
102 
103  PetscBool is_partitioned = PETSC_FALSE;
104  ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
105  if(is_partitioned) {
106  //Read mesh to MOAB PARALLEL=READ_DELETE PARALLEL=READ_PART
107  const char *option;
108  // option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
109  option = "PARALLEL=READ_PART;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
110  // rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
111  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
112  } else {
113  const char *option;
114  option = "";
115  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
116  }
117 
118 
119  PetscBool is_labatto = PETSC_TRUE;
120  ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-lobatto",&is_labatto,NULL); CHKERRQ(ierr);
121 
122 
123  // Create MoFEM (cephas) database
124  MoFEM::Core core(moab);
125  MoFEM::Interface& m_field = core;
126 
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); CHKERRQ_MOAB(rval);
133  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
134 
135  //Fields
136  if(is_labatto) {
137  ierr = m_field.add_field("rePRES",H1,AINSWORTH_LOBATTO_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
138  ierr = m_field.add_field("imPRES",H1,AINSWORTH_LOBATTO_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
139  ierr = m_field.add_field("P",H1,AINSWORTH_LOBATTO_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr); // in time domain
140  } else {
141  ierr = m_field.add_field("rePRES",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
142  ierr = m_field.add_field("imPRES",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
143  ierr = m_field.add_field("P",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr); // in time domain
144  }
145 
146  //meshset consisting all entities in mesh
147  EntityHandle root_set = moab.get_root_set();
148  //add entities to field
149  ierr = m_field.add_ents_to_field_by_TETs(root_set,"rePRES"); CHKERRQ(ierr);
150  ierr = m_field.add_ents_to_field_by_TETs(root_set,"imPRES"); CHKERRQ(ierr);
151  ierr = m_field.add_ents_to_field_by_TETs(root_set,"P"); CHKERRQ(ierr);
152  /* FIXME add wave number A priori error estimator in here as a field, which
153  store the error on the elements of computational domain. */
154 
155  //set app. order
156  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
157  PetscInt order;
158  ierr = PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
159  if(flg != PETSC_TRUE) {
160  order = 2;
161  }
162 
163  PetscBool wavenumber_flg;
164  double wavenumber;
165  // set wave number from line command, that overwrite numbre form block set
166  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-wave_number",&wavenumber,&wavenumber_flg); CHKERRQ(ierr);
167  if(!wavenumber_flg) {
168 
169  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"wave number not given, set in line command -wave_number to fix problem");
170 
171  }
172 
173 
174  PetscBool aDaptivity = PETSC_FALSE;
175  ierr = PetscOptionsGetBool(PETSC_NULL,"-adaptivity",&aDaptivity,NULL); CHKERRQ(ierr);
176 
177 
178  if(aDaptivity) {
179  /* A priori error indicator */
180  /* Introduced by BĂ©riot from Efficient implementation of high-order finite
181  elements for Helmholtz problems */
182  PetscInt threeErrorlv = 1;
183  ierr = PetscOptionsGetInt(PETSC_NULL,"-error_level",&threeErrorlv,NULL); CHKERRQ(ierr);
184  if(threeErrorlv > 2) {
185  threeErrorlv = 2;
186  }
187 
188  /* 1 => l2, 2 => h1, 3 => anisworth */
189  PetscInt l2 = 1;
190  ierr = PetscOptionsGetInt(PETSC_NULL,"-error_type",&l2,NULL); CHKERRQ(ierr);
191  if(l2 > 3) {
192  l2 = 3;
193  }
194 
195  /* 1 => maximum edge length, 2 => average edge length, 3 => minimum edge length */
196  pair<PetscInt,PetscBool> isCriterion;
197  isCriterion.first = 1;
198  ierr = PetscOptionsGetInt(PETSC_NULL,"-error_criterion",&isCriterion.first,&isCriterion.second); CHKERRQ(ierr);
199  if(isCriterion.first > 3) {
200  isCriterion.first = 3;
201  }
202 
203  Range tets;
204  map<int,Range> orders_map;
205 
206  rval = m_field.get_moab().get_entities_by_type(0,MBTET,tets,false); CHKERRQ_MOAB(rval);
207  // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,"MAT_HELMHOLTZ",it)) {
208  // rval = mField.get_moab().get_entities_by_type(it->meshset,MBTET,tets,true); CHKERRQ_MOAB(rval);
209  // }
210 
211  for(Range::iterator tit = tets.begin();tit!=tets.end();tit++) {
212 
213  double sum = 0;
214  double eDges [6];
215  // std::vector<double> valarryEdge(6);
216  for(int ee = 0;ee<6;ee++) {
217 
218  EntityHandle edge;
219  rval = m_field.get_moab().side_element(*tit,1,ee,edge); CHKERRQ_MOAB(rval);
220  const EntityHandle* conn;
221  int num_nodes;
222  rval = m_field.get_moab().get_connectivity(edge,conn,num_nodes,true); CHKERRQ_MOAB(rval);
223 
224  double coords[num_nodes*3];
225  rval = m_field.get_moab().get_coords(conn,num_nodes,coords); CHKERRQ_MOAB(rval);
226  cblas_daxpy(3,-1,&coords[0],1,&coords[3],1);
227  // for(int ii = 0;ii<6;ii++) {
228  // cout << "\n number " << ii << endl;
229  // cout << " coords after cblas_daxpy = \n" << coords[ii] << endl;
230  // }
231  eDges[ee] = cblas_dnrm2(3,&coords[3],1);
232  sum += eDges[ee];
233  // valarryEdge[ee] = cblas_dnrm2(3,&coords[3],1);
234 
235  }
236 
237  /* more robust to select average length of edges of elements than
238  max or min */
239 
240  double averageEdgeK = (sum/6.) * wavenumber;
241  double maxEdge = *std::max_element(eDges,eDges+6);
242  double minEdge = *std::min_element(eDges,eDges+6);
243 
244  /* assign element size definition based on mesh quality */
245  // double qualityIndicator = 3.*minEdge/maxEdge;
246  // cerr << "\n averageEdgeK = \n" << averageEdgeK << endl;
247  // if(qualityIndicator <= 1.0) {
248  // cerr << "\n qualityIndicator = \n" << qualityIndicator << endl;
249  // }
250  // if( (qualityIndicator < 1.5) || (qualityIndicator > 0.) ) {
251  // averageEdgeK = maxEdge * wavenumber;
252  // }
253 
254  if(isCriterion.first == 1) {
255  averageEdgeK = maxEdge * wavenumber;
256  } else if(isCriterion.first == 2) {
257  averageEdgeK = minEdge * wavenumber;;
258  }
259 
260  // Caulate order from table of 1D error test.
261  /* FIXME the order of p needs to be extend to more than 10 in future implementation */
262 
263  std::vector<double> errorLevel (10);
264 
265  double error_level [10];
266  if(threeErrorlv == 0) {
267  if(l2 == 1) {
268  /* target l2 error 15% */
269  /* error estimator in paper */
270  // memcpy(error_level, (double[]) {1.5,2.9,4.6,6.4,8.1,10.1,11.8,13.7,15.5,17.4}, sizeof error_level);
271  /* error estimator 1 */
272  const double a[] = {0.5,2.900,4.8,7.5,7.97,12.33,12.54,14.6,15.2,19.43};
273  memcpy(error_level, a , sizeof error_level);
274  /* error estimator 3 */
275  // memcpy(error_level, (double[]){1.62,2.900,4.52,6.52,8.275,9.995,11.885,13.645,15.515,17.6125}, sizeof error_level);
276  } else if(l2 == 2) {
277  /* error estimator h1 1 */
278  const double a[] = {0.2,2.3,3.88,6.8,7.1,11,11.17,14.54,15.28,19.023};
279  memcpy(error_level, a, sizeof error_level);
280  /* error estimator h1 3 */
281  // memcpy(error_level, (double[]) {0.675,1.956,4.255,6.49,8.315,9.085,11.445,13.650,15.765,17.9125}, sizeof error_level);
282  } else {
283  /* Anisworth error estimator */
284  // error_level [ ] = { };
285  }
286  errorLevel.insert (errorLevel.begin(),error_level,error_level+10);
287  // for(int ii = 0;ii<9;ii++) {
288  // cerr << " \n error_level[ii] = \n " << error_level[ii] << endl;
289  // }
290  } else if(threeErrorlv == 1) {
291 
292  if(l2 == 1) {
293  /* target l2 error 5% */
294  // memcpy(error_level, (double[]) {0.8,2.0,3.4,5.0,6.6,8.4,10.1,11.9,13.7,15.4}, sizeof error_level);a
295  const double a[] = {0.1,1.85,3.67,6.6,6.76,8.96,10.3,13.43,13.7,17};
296  memcpy(error_level,a, sizeof error_level);
297  // memcpy(error_level, (double[]) {0.776,1.98,3.755,5.25,6.585,8.385,10.395,12.170,13.765,15.3125}, sizeof error_level);
298  } else if(l2 == 2) {
299  const double a[] = {0.05,1.85,3.25,5.5,5.6,8.4,9.67,12.74,12.75,15.04};
300  memcpy(error_level,a , sizeof error_level);
301  // memcpy(error_level, (double[]) {0.049535,1.31,3.455,3.78,5.585,7.685,10.039,12.230,12.6965,14.6325}, sizeof error_level);
302  } else {
303  // memcpy(error_level, (double[]) { };
304  }
305  errorLevel.insert (errorLevel.begin(),error_level,error_level+10);
306 
307  } else if(threeErrorlv == 2) {
308 
309  if(l2 == 1) {
310  /* target l2 error 0.5% */
311  // memcpy(error_level, (double[]) {0.2,0.9,1.9,3.1,4.5,5.9,7.4,8.9,10.6,12.2}, sizeof error_level);
312  const double a[] = {0.05,1,1.33,3.6,4.6,7.05,7.35,9.4,10.76,13.625};
313  memcpy(error_level,a ,sizeof error_level);
314  // memcpy(error_level, (double[]) {0.249,1.085,2.175,3.11,4.625,6.745,7.439,8.815,10.7965,13.0825}, sizeof error_level);
315  } else if(l2 == 2) {
316  const double a[] = {0.02,0.41,0.83,2.35,3.64,6.25,6.28,8.46,9.83,12.62};
317  memcpy(error_level,a , sizeof error_level);
318  // memcpy(error_level, (double[]) {0.120,0.685,1.155,2.175,3.905,6.145,6.149,7.915,10.1865,12.4965}, sizeof error_level);
319  } else {
320  // memcpy(error_level, (double[]) { };
321  }
322  errorLevel.insert (errorLevel.begin(),error_level,error_level+10);
323 
324  }
325 
326  int orderP = 1;
327  for(int ii = 1;ii<8;ii++) {
328 
329  if(errorLevel[ii] < averageEdgeK) {
330  orderP = ii + 1;
331  } else {
332  break;
333  }
334 
335  }
336  // cerr << "\n orderP = " << orderP << endl;
337  // double max = 1.7*pow(10,308) + 1.0;
338  // cerr << "\n exceed double maximum = \n" << max << endl;
339  // std::cout << "Minimum value for int: " << std::numeric_limits<int>::min() << '\n';
340  // std::cout << "Maximum value for int: " << std::numeric_limits<int>::max() << '\n';
341 
342  // input a prior error indicator by ascending order with Range.
343  orders_map[orderP].insert(*tit);
344 
345  }
346 
347  map<int,Range>::iterator mit = orders_map.begin();
348  int max_value = mit->first;
349  int min_value = 10;
350  int avg_value = 0;
351 
352  for(;mit!=orders_map.end();mit++) {
353 
354  if(mit->first > max_value) {
355  max_value = mit->first;
356  }
357  if(min_value > mit->first) {
358  min_value = mit->first;
359  }
360  avg_value += mit->first;
361  // cout << "\n mit->first = \n " << mit->first << endl;
362 
363  Range tris;
364  rval = moab.get_adjacencies(mit->second,2,false,tris,moab::Interface::UNION);
365  Range edges;
366  rval = moab.get_adjacencies(mit->second,1,false,edges,moab::Interface::UNION);
367 
368  /* WARNING : the order of triangle and edges are consistent with tets in current implementation */
369  // ierr = m_field.set_field_order(mit->second,"rePRES",mit->first,2); CHKERRQ(ierr);
370  ierr = m_field.set_field_order(mit->second,"rePRES",mit->first); CHKERRQ(ierr);
371  ierr = m_field.set_field_order(tris,"rePRES",mit->first); CHKERRQ(ierr);
372  ierr = m_field.set_field_order(edges,"rePRES",mit->first); CHKERRQ(ierr);
373  //
374  ierr = m_field.set_field_order(mit->second,"imPRES",mit->first); CHKERRQ(ierr);
375  ierr = m_field.set_field_order(tris,"imPRES",mit->first); CHKERRQ(ierr);
376  ierr = m_field.set_field_order(edges,"imPRES",mit->first); CHKERRQ(ierr);
377 
378  ierr = m_field.set_field_order(mit->second,"P",mit->first); CHKERRQ(ierr);
379  ierr = m_field.set_field_order(tris,"P",mit->first); CHKERRQ(ierr);
380  ierr = m_field.set_field_order(edges,"P",mit->first); CHKERRQ(ierr);
381 
382  }
383 
384  if(!pcomm->rank()) {
385  printf("max_value_order_p %u\n",max_value);
386  printf("min_value_order_p %u\n",min_value);
387  printf("avg_value_order_p %u\n",int (avg_value/orders_map.size()));
388  }
389 
390  } else {
391 
392  ierr = m_field.set_field_order(root_set,MBTET,"rePRES",order); CHKERRQ(ierr);
393  ierr = m_field.set_field_order(root_set,MBTRI,"rePRES",order); CHKERRQ(ierr);
394  ierr = m_field.set_field_order(root_set,MBEDGE,"rePRES",order); CHKERRQ(ierr);
395 
396  ierr = m_field.set_field_order(root_set,MBTET,"imPRES",order); CHKERRQ(ierr);
397  ierr = m_field.set_field_order(root_set,MBTRI,"imPRES",order); CHKERRQ(ierr);
398  ierr = m_field.set_field_order(root_set,MBEDGE,"imPRES",order); CHKERRQ(ierr);
399 
400  ierr = m_field.set_field_order(root_set,MBTET,"P",order); CHKERRQ(ierr);
401  ierr = m_field.set_field_order(root_set,MBTRI,"P",order); CHKERRQ(ierr);
402  ierr = m_field.set_field_order(root_set,MBEDGE,"P",order); CHKERRQ(ierr);
403 
404  }
405 
406 
407  ierr = m_field.set_field_order(root_set,MBVERTEX,"rePRES",1); CHKERRQ(ierr);
408  ierr = m_field.set_field_order(root_set,MBVERTEX,"imPRES",1); CHKERRQ(ierr);
409  ierr = m_field.set_field_order(root_set,MBVERTEX,"P",1); CHKERRQ(ierr);
410 
411  if(!m_field.check_field("MESH_NODE_POSITIONS")) {
412 
413  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
414  ierr = m_field.add_ents_to_field_by_TETs(root_set,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
415  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
416  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
417  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
418  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
419 
420  ierr = m_field.build_fields(); CHKERRQ(ierr);
421  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
422  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
423 
424  } else {
425 
426  ierr = m_field.build_fields(); CHKERRQ(ierr);
427 
428  }
429 
430  // Finite Elements
431 
432  HelmholtzElement helmholtz_element(m_field);
433  ierr = helmholtz_element.getGlobalParametersFromLineCommandOptions(); CHKERRQ(ierr);
434  ierr = helmholtz_element.addHelmholtzElements("rePRES","imPRES"); CHKERRQ(ierr);
435  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
436  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","reEX"); CHKERRQ(ierr);
437  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","imEX"); CHKERRQ(ierr);
438  // ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_IMIM_FE","reEX"); CHKERRQ(ierr);
439  // ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_IMIM_FE","imEX"); CHKERRQ(ierr);
440  }
441 
442  bool Dirichlet_bc_set = false;
443  Range bc_dirichlet_tris,analytical_bc_tris;
444  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"ANALYTICAL_BC",it)) {
445  rval = moab.get_entities_by_type(it->getMeshset(),MBTRI,analytical_bc_tris,true); CHKERRQ_MOAB(rval);
446  Dirichlet_bc_set = true;
447  }
448  bc_dirichlet_tris.merge(analytical_bc_tris);
449  AnalyticalDirichletHelmholtzBC analytical_bc_real(m_field);
450  AnalyticalDirichletHelmholtzBC analytical_bc_imag(m_field);
451  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",analytical_bc_tris); CHKERRQ(ierr);
452  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",analytical_bc_tris); CHKERRQ(ierr);
453 
454 
455  double aTtenuation = 0;
456  // set loss tangent from line command, the complex wave number to include attenuation.
457  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-attenuation",&aTtenuation,NULL); CHKERRQ(ierr);
458 
459  ierr = PetscOptionsBegin(m_field.get_comm(),NULL,"Helmholtz problem options","none"); CHKERRQ(ierr);
460  PetscBool isRayleigh = PETSC_FALSE;
461  ierr = PetscOptionsBool(
462  "-rayleigh_wave",
463  "If true the incident is Rayleigh wave, otherwise, it is an non-attenuated wave","",
464  PETSC_FALSE,
465  &isRayleigh,
466  NULL
467  ); CHKERRQ(ierr);
468 
469  double power_of_incident_wave = 1;
470  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-amplitude_of_incident_wave",&power_of_incident_wave,NULL); CHKERRQ(ierr);
471 
472 
473  // This is added for a case than on some surface, defined by the user a
474  // incident plane wave is scattered.
475  map<int,PlaneIncidentWaveSacttrerData> planeWaveScatterData;
476  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"SOFT_INCIDENT_WAVE_BC",it)) {
477  rval = moab.get_entities_by_type(it->getMeshset(),MBTRI,planeWaveScatterData[it->getMeshsetId()].tRis,true); CHKERRQ_MOAB(rval);
478  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",planeWaveScatterData[it->getMeshsetId()].tRis); CHKERRQ(ierr);
479  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",planeWaveScatterData[it->getMeshsetId()].tRis); CHKERRQ(ierr);
480  bc_dirichlet_tris.merge(planeWaveScatterData[it->getMeshsetId()].tRis);
481 
482  Dirichlet_bc_set = true;
483 
484  }
485 
486  //ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
487  //// Build adjacencies
488  //ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
489 
490  // Problem
491  ierr = m_field.add_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
492  ierr = m_field.add_problem("BCREAL_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for real field
493  ierr = m_field.add_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for imag field
494 
495  // Set refinement level for problem
496  ierr = m_field.modify_problem_ref_level_add_bit("ACOUSTIC_PROBLEM",bit_level0); CHKERRQ(ierr);
497  ierr = m_field.modify_problem_ref_level_add_bit("BCREAL_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
498  ierr = m_field.modify_problem_ref_level_add_bit("BCIMAG_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
499 
500  // Add elements to problems
501  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
502  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_IMIM_FE"); CHKERRQ(ierr);
503  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_REIM_FE"); CHKERRQ(ierr);
504 
505  ierr = m_field.modify_problem_add_finite_element("BCREAL_PROBLEM","BCREAL_FE"); CHKERRQ(ierr);
506  ierr = m_field.modify_problem_add_finite_element("BCIMAG_PROBLEM","BCIMAG_FE"); CHKERRQ(ierr);
507 
508  // Get start time for analyse
509  PetscLogDouble t1,t2;
510  PetscLogDouble v1,v2;
511  ierr = PetscTime(&v1); CHKERRQ(ierr);
512  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
513  PetscLogDouble m1;
514  PetscLogDouble m2;
515  ierr = PetscMemoryGetCurrentUsage(&m1); CHKERRQ(ierr);
516 
517  // Build problems
518  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
519  // Build adjacencies
520  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
521 
522  // field split
523  ublas::matrix<Mat> nested_matrices(2,2);
524  ublas::vector<IS> nested_is_rows(2);
525 
526  // ublas::vector<IS> nested_is_cols(2);
527  // IS nested_is_rows_test = PETSC_NULL;
528  // IS nested_is_cols_test = PETSC_NULL;
529  for(int i = 0;i!=2;i++) {
530  nested_is_rows[i] = PETSC_NULL;
531 
532  // nested_is_cols[i] = PETSC_NULL;
533  for(int j = 0;j!=2;j++) {
534  nested_matrices(i,j) = PETSC_NULL;
535  }
536  }
537  //FIXME
538  // Mat A re re position
539  // const Problem *prb_ptr;
540  // ierr = m_field.get_problem("ACOUSTIC_PROBLEM",&prb_ptr); CHKERRQ(ierr);
541  // boost::shared_ptr<Problem::SubProblemData> sub_data = prb_ptr->getSubData();
542  // ierr = sub_data->getRowIs(&nested_is_rows[0]); CHKERRQ(ierr);
543  // ierr = sub_data->getColIs(&nested_is_cols[0]); CHKERRQ(ierr);
544  // // Mat A in im im position
545  // boost::shared_ptr<Problem::SubProblemData> sub_data = problem_ptr->getSubData();
546  // ierr = sub_data->getRowIs(&nested_is_rows[1]); CHKERRQ(ierr);
547  // ierr = sub_data->getColIs(&nested_is_cols[1]); CHKERRQ(ierr);
548  // // Mat -C in re im position
549  // boost::shared_ptr<Problem::SubProblemData> sub_data = problem_ptr->getSubData();
550  // ierr = sub_data->getRowIs(&nested_is_rows[0]); CHKERRQ(ierr);
551  // ierr = sub_data->getColIs(&nested_is_cols[1]); CHKERRQ(ierr);
552  // // Mat C in im re position
553  // boost::shared_ptr<Problem::SubProblemData> sub_data = problem_ptr->getSubData();
554  // ierr = sub_data->getRowIs(&nested_is_rows[1]); CHKERRQ(ierr);
555  // ierr = sub_data->getColIs(&nested_is_cols[0]); CHKERRQ(ierr);
556  // That will be filled at the end
557  //field split
558 
559  // build porblems
560  if(is_partitioned) {
561  // if mesh is partitioned
562 
563  ierr = m_field.build_problem_on_distributed_mesh("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
564  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
565 
566  if(Dirichlet_bc_set) {
567 
568  ierr = m_field.build_problem_on_distributed_mesh("BCREAL_PROBLEM",true); CHKERRQ(ierr);
569  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM",true); CHKERRQ(ierr);
570 
571  ierr = m_field.build_problem_on_distributed_mesh("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
572  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
573 
574  }
575 
576 
577  } else {
578  // if not partitioned mesh is load to all processes
579 
580  ierr = m_field.build_problem("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
581  ierr = m_field.partition_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
582  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
583 
584  if(Dirichlet_bc_set) {
585  ierr = m_field.build_problem("BCREAL_PROBLEM",true); CHKERRQ(ierr);
586  ierr = m_field.partition_problem("BCREAL_PROBLEM"); CHKERRQ(ierr);
587  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM"); CHKERRQ(ierr);
588 
589 
590 
591  ierr = m_field.build_problem("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
592  ierr = m_field.partition_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr);
593  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM"); CHKERRQ(ierr);
594  }
595  }
596 
597  ierr = m_field.partition_ghost_dofs("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
598  if(Dirichlet_bc_set) {
599  ierr = m_field.partition_ghost_dofs("BCREAL_PROBLEM"); CHKERRQ(ierr);
600  ierr = m_field.partition_ghost_dofs("BCIMAG_PROBLEM"); CHKERRQ(ierr);
601  }
602 
603 
604 
605 
606  //field split
607  // Mat SubA;
608  // ierr = MatCreateNest(
609  // PETSC_COMM_WORLD,
610  // 2,&sub_nested_is_rows[0],
611  // 2,&sub_nested_is_cols[0],
612  // &sub_nested_matrices(0,0),
613  // &SubA
614  // ); CHKERRQ(ierr);
615  // nested_matrices(0,0) = SubA;
616  //
617  // ierr = MatAssemblyBegin(SubA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
618  // ierr = MatAssemblyEnd(SubA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
619  //
620  // if(debug) {
621  // cerr << "Nested SubA" << endl;
622  // MatView(SubA,PETSC_VIEWER_STDOUT_WORLD);
623  // }
624  //field split
625 
626  // Get problem matrices and vectors
627  Vec F; //Right hand side vector
628  ierr = m_field.VecCreateGhost("ACOUSTIC_PROBLEM",ROW,&F); CHKERRQ(ierr);
629  Vec T; //Solution vector
630  ierr = VecDuplicate(F,&T); CHKERRQ(ierr);
631  Mat A; //Left hand side matrix
632 
633  ierr = m_field.MatCreateMPIAIJWithArrays("ACOUSTIC_PROBLEM",&A); CHKERRQ(ierr);
634  /* change matrix type to suit ilu */
635  // Mat B;
636  // ierr = MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&A);
637  // ierr = MatDestroy(&B); CHKERRQ(ierr);
638  // ierr = MatSetType(A,MATBAIJ);
639  // ierr = m_field.MatCreateMPIBAIJWithArrays("ACOUSTIC_PROBLEM",&A); CHKERRQ(ierr);
640  /* change matrix type to suit ilu end */
641  ierr = helmholtz_element.setOperators(A,F,"rePRES","imPRES"); CHKERRQ(ierr);
642 
643 
644  //wave direction unit vector=[x,y,z]^T
645  VectorDouble wave_direction;
646  wave_direction.resize(3);
647  wave_direction.clear();
648  wave_direction[2] = 1; // default:X direction [0,0,1]
649 
650  int nmax = 3;
651  ierr = PetscOptionsGetRealArray(PETSC_NULL,PETSC_NULL,"-wave_direction",&wave_direction[0],&nmax,NULL); CHKERRQ(ierr);
652  if(nmax > 0 && nmax != 3) {
653  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"*** ERROR -wave_direction [3*1 vector] default:X direction [0,0,1]");
654  }
655 
656  PetscInt choise_value = NO_ANALYTICAL_SOLUTION;
657  // set type of analytical solution
658  ierr = PetscOptionsGetEList(PETSC_NULL,NULL,"-analytical_solution_type",analytical_solution_types,7,&choise_value,NULL); CHKERRQ(ierr);
659 
660  switch((AnalyticalSolutionTypes)choise_value) {
661 
663 
664  {
665  double scattering_sphere_radius = 0.5;
666  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-scattering_sphere_radius",&scattering_sphere_radius,NULL); CHKERRQ(ierr);
667 
668 
669 
670  boost::shared_ptr<HardSphereScatterWave> function_evaluator = boost::shared_ptr<HardSphereScatterWave>(
671  new HardSphereScatterWave(wavenumber,scattering_sphere_radius)
672  );
673  ierr = analytical_bc_real.setApproxOps(
674  m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL
675  ); CHKERRQ(ierr);
676  ierr = analytical_bc_imag.setApproxOps(
677  m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG
678  ); CHKERRQ(ierr);
679  Dirichlet_bc_set = true;
680 
681  }
682 
683  break;
684 
686 
687  {
688 
689  double scattering_sphere_radius = 0.5;
690  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-scattering_sphere_radius",&scattering_sphere_radius,NULL); CHKERRQ(ierr);
691 
692 
693  boost::shared_ptr<SoftSphereScatterWave> function_evaluator = boost::shared_ptr<SoftSphereScatterWave>(new SoftSphereScatterWave(wavenumber,scattering_sphere_radius));
694  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
695  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
696  Dirichlet_bc_set = true;
697 
698  }
699 
700  break;
701 
702  case PLANE_WAVE:
703 
704  {
705 
706  double angle = 0.25;
707  // set wave number from line command, that overwrite numbre form block set
708  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-wave_guide_angle",&angle,NULL); CHKERRQ(ierr);
709 
710 
711  boost::shared_ptr<PlaneWave> function_evaluator = boost::shared_ptr<PlaneWave>(new PlaneWave(wavenumber,angle*M_PI));
712  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
713  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
714  Dirichlet_bc_set = true;
715 
716  }
717 
718  break;
719 
721 
722  {
723 
724  boost::shared_ptr<HardCylinderScatterWave> function_evaluator = boost::shared_ptr<HardCylinderScatterWave>(new HardCylinderScatterWave(wavenumber));
725  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
726  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
727  Dirichlet_bc_set = true;
728 
729 
730  }
731 
732  break;
733 
735 
736  {
737 
738  boost::shared_ptr<SoftCylinderScatterWave> function_evaluator = boost::shared_ptr<SoftCylinderScatterWave>(new SoftCylinderScatterWave(wavenumber));
739  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
740  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
741  Dirichlet_bc_set = true;
742 
743  }
744 
745  break;
746 
748 
749  {
750 
751  boost::shared_ptr<SingularScatterWave> function_evaluator = boost::shared_ptr<SingularScatterWave>(new SingularScatterWave(wavenumber));
752  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
753  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
754  Dirichlet_bc_set = true;
755 
756  }
757 
758  break;
759 
760  case INCIDENT_WAVE:
761 
762  {
763 
764  boost::shared_ptr<IncidentWave> function_evaluator =
765  boost::shared_ptr<IncidentWave>(new IncidentWave(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
766  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,power_of_incident_wave));
767  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
768  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
769  Dirichlet_bc_set = true;
770 
771  }
772 
773  break;
774 
776 
777  {
778  // Dirichlet_bc_set = false;
779  }
780 
781  break;
782 
783  }
784 
785  // Analytical boundary conditions
786  AnalyticalDirichletHelmholtzBC::DirichletBC analytical_ditihlet_bc_real(m_field,"rePRES",A,T,F);
787  AnalyticalDirichletHelmholtzBC::DirichletBC analytical_ditihlet_bc_imag(m_field,"imPRES",A,T,F);
788 
789 
790  PetscBool monochromatic_wave = PETSC_TRUE;
791  ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-monochromatic_wave",&monochromatic_wave,NULL); CHKERRQ(ierr);
792 
793  PetscBool add_incident_wave = PETSC_FALSE;
794  ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-add_incident_wave",&add_incident_wave,NULL); CHKERRQ(ierr);
795  if(add_incident_wave) {
796 
797  // define problem
798  ierr = m_field.add_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
799  // set finite elements for problem
800  ierr = m_field.modify_problem_add_finite_element("INCIDENT_WAVE","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
801  // set refinment level for problem
802  ierr = m_field.modify_problem_ref_level_add_bit("INCIDENT_WAVE",bit_level0); CHKERRQ(ierr);
803 
804  // build porblems
805  if(is_partitioned) {
806  // if mesh is partitioned
807  ierr = m_field.build_problem_on_distributed_mesh("INCIDENT_WAVE",true); CHKERRQ(ierr);
808  ierr = m_field.partition_finite_elements("INCIDENT_WAVE",true); CHKERRQ(ierr);
809  } else {
810  // if not partitioned mesh is load to all processes
811  ierr = m_field.build_problem("INCIDENT_WAVE",true); CHKERRQ(ierr);
812  ierr = m_field.partition_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
813  ierr = m_field.partition_finite_elements("INCIDENT_WAVE"); CHKERRQ(ierr);
814  }
815  ierr = m_field.partition_ghost_dofs("INCIDENT_WAVE"); CHKERRQ(ierr);
816 
817  }
818 
819 
820  KSP solver;
821  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
822 
823  if(monochromatic_wave) {
824 
825  if(Dirichlet_bc_set) {
826 
827  {
828 
829  map<int,PlaneIncidentWaveSacttrerData>::iterator mit = planeWaveScatterData.begin();
830 
831  boost::shared_ptr<IncidentWave> function_evaluator;
832  boost::shared_ptr<IncidentWavePointSource> functionEvaluator;
833  for(;mit!=planeWaveScatterData.end();mit++) {
834 
835  if(!helmholtz_element.globalParameters.sourceCoordinate.second) {
836  if(helmholtz_element.globalParameters.isRadiation.first) {
837  // note negative field, scatter field should cancel incident wave
838  function_evaluator = boost::shared_ptr<IncidentWave>(
839  new IncidentWave(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
840  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,-helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
841  );
842  } else {
843  function_evaluator = boost::shared_ptr<IncidentWave>(
844  new IncidentWave(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
845  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
846  );
847  }
848  ierr = analytical_bc_real.setApproxOps(
849  m_field,"rePRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::REAL
850  ); CHKERRQ(ierr);
851  ierr = analytical_bc_imag.setApproxOps(
852  m_field,"imPRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::IMAG
853  ); CHKERRQ(ierr);
854  } else {
855  if(helmholtz_element.globalParameters.isRadiation.first) {
856  // note negative field, scatter field should cancel incident wave
857  functionEvaluator = boost::shared_ptr<IncidentWavePointSource>(
858  new IncidentWavePointSource(wavenumber,helmholtz_element.globalParameters.sourceCoordinate.first,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
859  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,-helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
860  );
861  } else {
862  functionEvaluator = boost::shared_ptr<IncidentWavePointSource>(
863  new IncidentWavePointSource(wavenumber,helmholtz_element.globalParameters.sourceCoordinate.first,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
864  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
865  );
866  }
867  ierr = analytical_bc_real.setApproxOps(
868  m_field,"rePRES",mit->second.tRis,functionEvaluator,GenericAnalyticalSolution::REAL
869  ); CHKERRQ(ierr);
870  ierr = analytical_bc_imag.setApproxOps(
871  m_field,"imPRES",mit->second.tRis,functionEvaluator,GenericAnalyticalSolution::IMAG
872  ); CHKERRQ(ierr);
873  }
874 
875  }
876 
877  }
878  // Solve for analytical Dirichlet bc dofs
879  ierr = analytical_bc_real.setProblem(m_field,"BCREAL_PROBLEM"); CHKERRQ(ierr);
880  ierr = analytical_bc_imag.setProblem(m_field,"BCIMAG_PROBLEM"); CHKERRQ(ierr);
881  ierr = analytical_bc_real.solveProblem(
882  m_field,"BCREAL_PROBLEM","BCREAL_FE",analytical_ditihlet_bc_real,bc_dirichlet_tris
883  ); CHKERRQ(ierr);
884  ierr = analytical_bc_imag.solveProblem(
885  m_field,"BCIMAG_PROBLEM","BCIMAG_FE",analytical_ditihlet_bc_imag,bc_dirichlet_tris
886  ); CHKERRQ(ierr);
887 
888  ierr = analytical_bc_real.destroyProblem(); CHKERRQ(ierr);
889  ierr = analytical_bc_imag.destroyProblem(); CHKERRQ(ierr);
890 
891  }
892 
893  /* calculate the interior problem to set the fulid velocity */
894  if(helmholtz_element.globalParameters.fRequency.second) {
895  const complex< double > i( 0.0, 1.0 );
896  complex< double > wave_number = (2*M_PI*helmholtz_element.globalParameters.fRequency.first / helmholtz_element.globalParameters.vElocity.first);
897  // FIXME Why it is not working with first method ?
898  // map<int,helmholtz_element::VolumeData>::iterator vit = helmholtz_element.volumeData.begin();
899  //
900  // for(;vit != helmholtz_element.volumeData.end();vit++) {
901  // vit->second.waveNumber = wave_number;
902  // }
903 
905  if(it->getName().compare(0,13,"MAT_HELMHOLTZ") == 0) {
906  helmholtz_element.volumeData[it->getMeshsetId()].waveNumber = wave_number;
907  }
908  }
909 
910  }
911 
912 
913  // Zero vectors
914  ierr = VecZeroEntries(T); CHKERRQ(ierr);
915  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
916  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
917  ierr = VecZeroEntries(F); CHKERRQ(ierr);
918  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
919  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
920  ierr = MatZeroEntries(A); CHKERRQ(ierr);
921 
922  // Assemble problem
923  if(Dirichlet_bc_set) {
924  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
925  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
926  }
927 
928  ierr = helmholtz_element.calculateA("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
929  ierr = helmholtz_element.calculateF("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
930  // ierr = MatView(A,PETSC_VIEWER_DRAW_WORLD);
931  // std::string wait;
932  // std::cin >> wait;
933 
934  if(Dirichlet_bc_set) {
935  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
936  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
937  }
938 
939  // FIXME Point Source Right hand sides for vertex source for Total Acoustic potential
940  // \Phi = \Phi_{S} + \Phi_{I}
941  // const Problem *problem_ptr;
942  // ierr = m_field.get_problem("ACOUSTIC_PROBLEM",&problem_ptr); CHKERRQ(ierr);
943  // for(_IT_NUMEREDDOF_ROW_BY_NAME_ENT_PART_FOR_LOOP_(problem_ptr,"rePRES",node,m_field.getCommRank(),dit)) {
944  // int global_idx = dit-> getPetscGlobalDofIdx();
945  // cerr << "\n global_idx = \n\n " << global_idx << endl;
946  // ierr = VecSetValue(F,global_idx,helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first,INSERT_VALUES); CHKERRQ(ierr);
947  // }
948 
949 
950  ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
951  ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
952  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
953  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
954 
955 
956  // get ROW and COL indices
958  "ACOUSTIC_PROBLEM",
959  ROW,"rePRES",0,1,
960  &nested_is_rows[0]
961  ); CHKERRQ(ierr);
962  // ierr = m_field.ISCreateProblemFieldAndRank(
963  // "ACOUSTIC_PROBLEM",
964  // COL,"rePRES",0,1,
965  // &nested_is_cols[0]
966  // ); CHKERRQ(ierr);
968  "ACOUSTIC_PROBLEM",
969  ROW,"imPRES",0,1,
970  &nested_is_rows[1]
971  ); CHKERRQ(ierr);
972  // ierr = m_field.ISCreateProblemFieldAndRank(
973  // "ACOUSTIC_PROBLEM",
974  // COL,"imPRES",0,1,
975  // &nested_is_cols[1]
976  // ); CHKERRQ(ierr);
977 
978 
979  // field split
980  Mat SubA,SubC,SubConjC;
981 
982  ierr = MatGetSubMatrix(A,nested_is_rows[0],nested_is_rows[0],MAT_INITIAL_MATRIX,&SubA); CHKERRQ(ierr);
983  ierr = MatGetSubMatrix(A,nested_is_rows[1],nested_is_rows[0],MAT_INITIAL_MATRIX,&SubC); CHKERRQ(ierr);
984  ierr = MatGetSubMatrix(A,nested_is_rows[0],nested_is_rows[1],MAT_INITIAL_MATRIX,&SubConjC); CHKERRQ(ierr);
985 
986  ierr = MatAssemblyBegin(SubA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
987  ierr = MatAssemblyEnd(SubA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
988  ierr = MatAssemblyBegin(SubC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
989  ierr = MatAssemblyEnd(SubC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
990  ierr = MatAssemblyBegin(SubConjC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
991  ierr = MatAssemblyEnd(SubConjC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
992 
993  nested_matrices(0,0) = SubA;
994  nested_matrices(1,1) = SubA;
995  nested_matrices(1,0) = SubC;
996  nested_matrices(0,1) = SubConjC;
997 
998 
999  if(debug) {
1000  cerr << "Nested SubA" << endl;
1001  // MatView(SubA,PETSC_VIEWER_STDOUT_WORLD);
1002  ierr = MatView(SubA,PETSC_VIEWER_DRAW_WORLD);
1003  std::string wait;
1004  std::cin >> wait;
1005  cerr << "Nested SubC" << endl;
1006  // MatView(SubA,PETSC_VIEWER_STDOUT_WORLD);
1007  ierr = MatView(SubC,PETSC_VIEWER_DRAW_WORLD);
1008  std::cin >> wait;
1009  }
1010 
1011 
1012  Mat B;
1013  ierr = MatCreateNest(
1014  PETSC_COMM_WORLD,
1015  2,&nested_is_rows[0],
1016  2,&nested_is_rows[0],
1017  &nested_matrices(0,0),
1018  &B
1019  ); CHKERRQ(ierr);
1020 
1021  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1022  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1023 
1024  if(debug) {
1025  // PetscInt nestM,nestN,nestLocalIS1,nestGlobalIS1,nestLocalIS2,nestGlobalIS2;
1026  // ierr = MatNestGetISs(B,&nested_is_rows_global[0],NULL); CHKERRQ(ierr);
1027  // ierr = ISGetSize(nested_is_rows_global[0],&nestGlobalIS1); CHKERRQ(ierr);
1028  // ierr = ISGetLocalSize(nested_is_rows_global[0],&nestLocalIS1); CHKERRQ(ierr);
1029  // ierr = ISGetSize(nested_is_rows_global[1],&nestGlobalIS2); CHKERRQ(ierr);
1030  // ierr = ISGetLocalSize(nested_is_rows_global[1],&nestLocalIS2); CHKERRQ(ierr);
1031  // cerr << "\n nest matrix LocalIS1 = \n " << nestLocalIS1 << "\n nest matrix LocalIS2 = \n " << nestLocalIS2 << endl;
1032  // cerr << "\n nest matrix GlobalIS1 = \n " << nestGlobalIS1 << "\n nest matrix GlobalIS2 = \n " << nestGlobalIS2 << endl;
1033  std::string wait;
1034  // std::cin >> wait;
1035  ierr = MatView(A,PETSC_VIEWER_DRAW_WORLD);
1036  std::cin >> wait;
1037  }
1038 
1039 
1040 
1041  /* get the sparsity information of matrix */
1042  MatInfo info;
1043  // double nonZeroes;
1044  ierr = MatGetInfo(A,MAT_GLOBAL_SUM,&info);
1045  // nonZeroes = info.nz_allocated;
1046  // MatGetSize(A,&m,&n);
1047  if(!pcomm->rank()) {
1048  printf("info.nz_used %g\n",info.nz_used);
1049  printf("info.memory %g\n",info.memory/2073741824);
1050  // printf("info.factor_mallocs %g\n",info.factor_mallocs);
1051  // printf("info.assemblies %g\n",info.assemblies);
1052  }
1053  // ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1054 
1055  ierr = MatDestroy(&A); CHKERRQ(ierr);
1056  ierr = VecScale(F,-1); CHKERRQ(ierr);
1057 
1058  // double nrm2_F;
1059  // ierr = VecNorm(F,NORM_2,&nrm2_F); CHKERRQ(ierr);
1060  // cerr << "\n nomr of vector: \n" << nrm2_F << endl;
1061 
1062  // Solve problem
1063  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
1064  ierr = KSPSetOperators(solver,B,B); CHKERRQ(ierr);
1065 
1066  // field split
1067  PC pc;
1068  ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
1069  ierr = PCSetType(pc,PCFIELDSPLIT); CHKERRQ(ierr);
1070  PetscBool is_pcfs = PETSC_FALSE;
1071  PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);
1072  if(is_pcfs) {
1073 
1074  ierr = PCSetOperators(pc,B,B); CHKERRQ(ierr);
1075  ierr = PCFieldSplitSetIS(pc,NULL,nested_is_rows[0]); CHKERRQ(ierr);
1076  ierr = PCFieldSplitSetIS(pc,NULL,nested_is_rows[1]); CHKERRQ(ierr);
1077  // ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR); CHKERRQ(ierr);
1078  // ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_MULTIPLICATIVE); CHKERRQ(ierr);
1079  ierr = PCSetUp(pc); CHKERRQ(ierr);
1080  // ierr = PCFieldSplitSetSchurFactType(sub_pc_0,PC_FIELDSPLIT_SCHUR_FACT_LOWER); CHKERRQ(ierr);
1081  // ierr = PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_SCHUR); CHKERRQ(ierr);
1082 
1083  // ierr = PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_A11,SubA); CHKERRQ(ierr);
1084  // KSP *sub_ksp;
1085  // PetscInt n;
1086  // ierr = PCFieldSplitGetSubKSP(pc,&n,&sub_ksp); CHKERRQ(ierr);
1087  // ierr = KSPSetOperators(sub_ksp[1],SubA,SubA);
1088  // ierr = PetscFree(sub_ksp); CHKERRQ(ierr);
1089 
1090  } else {
1091  SETERRQ(
1092  PETSC_COMM_WORLD,
1094  "Only works with pre-conditioner PCFIELDSPLIT"
1095  );
1096  }
1097  // field split
1098 
1099  ierr = KSPSetUp(solver); CHKERRQ(ierr);
1100  ierr = KSPSolve(solver,F,T); CHKERRQ(ierr);
1101  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1102  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1103 
1104  if(is_partitioned) {
1105 
1106  // no need for global communication
1107  ierr = m_field.set_local_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
1108 
1109  } else {
1110 
1111  ierr = m_field.set_global_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
1112 
1113  }
1114 
1115  // ierr = PCDestroy(&pc); CHKERRQ(ierr);
1116  ierr = MatDestroy(&SubA); CHKERRQ(ierr);
1117  ierr = MatDestroy(&SubC); CHKERRQ(ierr);
1118  ierr = MatDestroy(&SubConjC); CHKERRQ(ierr);
1119  ierr = MatDestroy(&B); CHKERRQ(ierr);
1120  // ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1121 
1122 
1123  } else {
1124 
1125  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1126  // define problem
1127  ierr = m_field.add_problem("PRESSURE_IN_TIME"); CHKERRQ(ierr);
1128  // set finite elements for problem
1129  ierr = m_field.modify_problem_add_finite_element("PRESSURE_IN_TIME","PRESSURE_FE"); CHKERRQ(ierr);
1130  // set refinment level for problem
1131  ierr = m_field.modify_problem_ref_level_add_bit("PRESSURE_IN_TIME",bit_level0); CHKERRQ(ierr);
1132 
1133  // build porblems
1134  if(is_partitioned) {
1135  // if mesh is partitioned
1136  ierr = m_field.build_problem_on_distributed_mesh("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
1137  ierr = m_field.partition_finite_elements("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
1138  } else {
1139  // if not partitioned mesh is load to all processes
1140  ierr = m_field.build_problem("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
1141  ierr = m_field.partition_problem("PRESSURE_IN_TIME"); CHKERRQ(ierr);
1142  ierr = m_field.partition_finite_elements("PRESSURE_IN_TIME"); CHKERRQ(ierr);
1143  }
1144  ierr = m_field.partition_ghost_dofs("PRESSURE_IN_TIME"); CHKERRQ(ierr);
1145 
1146  TimeSeriesFieldSplit time_series(m_field,helmholtz_element,
1147  analytical_bc_real,analytical_bc_imag,
1148  analytical_ditihlet_bc_real,analytical_ditihlet_bc_imag,
1149  Dirichlet_bc_set,bc_dirichlet_tris,planeWaveScatterData);
1150 
1151  ierr = time_series.readData(); CHKERRQ(ierr);
1152  ierr = time_series.createPressureSeries(T); CHKERRQ(ierr);
1153  ierr = time_series.forwardSpaceDft(); CHKERRQ(ierr);
1154  // ierr = time_series.saveFrequencyData(); // save the frequency data to file
1155  ierr = time_series.solveForwardDFT(A,F,T); CHKERRQ(ierr);
1156  ierr = MatDestroy(&A); CHKERRQ(ierr);
1157  ierr = time_series.pressureForwardDft(); CHKERRQ(ierr);
1158  ierr = time_series.pressureInverseDft(); CHKERRQ(ierr);
1159  ierr = time_series.generateReferenceElementMesh(); CHKERRQ(ierr);
1160  ierr = time_series.saveResults(); CHKERRQ(ierr);
1161  ierr = time_series.destroyPressureSeries(); CHKERRQ(ierr);
1162 
1163  }
1164  // PetscInt its;
1165  // ierr = KSPGetIterationNumber(solver,&its);
1166  PetscScalar flops;
1167  ierr = PetscTime(&v2);CHKERRQ(ierr);
1168  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
1169  ierr = PetscGetFlops(&flops); CHKERRQ(ierr);
1170  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f S CPU Time = %f S \n",pcomm->rank(),v2-v1,t2-t1);
1171  ierr = PetscMemoryGetCurrentUsage(&m2); CHKERRQ(ierr);
1172  /* PetscMemoryGetMaximumUsage */
1173  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Memory usage = %f GBytes \n",pcomm->rank(),(m2-m1)/1073741824);
1174  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total flops %d \n",pcomm->rank(),flops);
1175 
1176 
1177  //Vec P,M;
1178  //ierr = m_field.VecCreateGhost("EX1_PROBLEM",COL,&M); CHKERRQ(ierr);
1179  //ierr = VecDuplicate(M,&P); CHKERRQ(ierr);
1180 
1181  //ierr = m_field.set_local_ghost_vector("EX1_PROBLEM",COL,M,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1182  //ierr = m_field.set_other_global_ghost_vector("EX1_PROBLEM","reEX","imEX",COL,P,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1183 
1184  //double nrm2_M;
1185  //ierr = VecNorm(M,NORM_2,&nrm2_M); CHKERRQ(ierr);
1186 
1187  //Vec V;
1188  //ierr = m_field.VecCreateGhost("ACOUSTIC_PROBLEM",COL,&V); CHKERRQ(ierr);
1189  //ierr = m_field.set_local_ghost_vector("ACOUSTIC_PROBLEM",COL,V,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1190 
1191  //VecScatter scatter_real,scatter_imag;
1192 
1193  //ierr = m_field.VecScatterCreate(V,"ACOUSTIC_PROBLEM","rePRES",COL,M,"EX1_PROBLEM","reEX",COL,&scatter_real); CHKERRQ(ierr);
1194 
1195  //ierr = m_field.VecScatterCreate(V,"ACOUSTIC_PROBLEM","imPRES",COL,P,"EX1_PROBLEM","reEX",COL,&scatter_imag); CHKERRQ(ierr);
1196 
1197  //VecScale(V,-1);
1198 
1199  //ierr = VecScatterBegin(scatter_real,V,M,ADD_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1200  //ierr = VecScatterEnd(scatter_real,V,M,ADD_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1201 
1202  //double nrm2_ErrM;
1203  //ierr = VecNorm(M,NORM_2,&nrm2_ErrM); CHKERRQ(ierr);
1204  //PetscPrintf(PETSC_COMM_WORLD,"L2 relative error on real field of acoustic problem %6.4e\n",(nrm2_ErrM)/(nrm2_M));
1205 
1206 
1207  //ierr = VecDestroy(&M); CHKERRQ(ierr);
1208  //ierr = VecDestroy(&P); CHKERRQ(ierr);
1209  //ierr = VecDestroy(&V); CHKERRQ(ierr);
1210 
1211 
1212 
1213  if(monochromatic_wave) {
1214 
1215  if(add_incident_wave) {
1216  /* FIXME problem with HELMHOLTZ_RERE_FE since this element contains 2 fields,and
1217  create big matrix did not consistent with field_approximation.hpp calculation method. */
1218 
1219  if(helmholtz_element.globalParameters.sourceCoordinate.second) {
1220  IncidentWavePointSource function_evaluator(wavenumber,helmholtz_element.globalParameters.sourceCoordinate.first,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
1221  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,power_of_incident_wave);
1222  ierr = solve_problem(m_field,"INCIDENT_WAVE","HELMHOLTZ_RERE_FE","rePRES","imPRES",
1223  ADD_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
1224  } else {
1225  if(helmholtz_element.globalParameters.complexWaveNumber.second) {
1226  wave_direction = helmholtz_element.globalParameters.waveOscilationDirection.first;
1227  }
1228  IncidentWave function_evaluator(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
1229  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,power_of_incident_wave);
1230  ierr = solve_problem(m_field,"INCIDENT_WAVE","HELMHOLTZ_RERE_FE","rePRES","imPRES",
1231  ADD_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
1232  }
1233 
1234  }
1235 
1236  if(is_partitioned) {
1237 
1238  rval = moab.write_file("fe_solution.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERRQ_MOAB(rval);
1239 
1240  } else {
1241 
1242  if(!pcomm->rank()) {
1243  rval = moab.write_file("fe_solution.h5m"); CHKERRQ_MOAB(rval);
1244  }
1245 
1246  }
1247 
1248 
1249 
1250  PetscBool save_postproc_mesh = PETSC_TRUE;
1251  ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-save_postproc_mesh",&save_postproc_mesh,NULL); CHKERRQ(ierr);
1252  if(save_postproc_mesh) {
1253 
1255  // PetscBool pressure_field = PETSC_FALSE;
1256  /* FIXME post-process the pressure field from acoustic potential */
1257  // ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-pressure_field",&pressure_field,NULL); CHKERRQ(ierr);
1258  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
1259  // PetscScalar scaling_pressure;
1260  // if(pressure_field) {
1261  // scaling_pressure = - 2.0 * M_PI * helmholtzElement.globalParameters.fRequency.first * helmholtzElement.globalParameters.dEnsity.first;
1262  // ierr = VecScale(p_scatter_wave_real,scaling_pressure);
1263  // }
1264  ierr = post_proc.addFieldValuesPostProc("rePRES"); CHKERRQ(ierr);
1265  // if(pressure_field) {
1266  // scaling_pressure *= - 1.0;
1267  // ierr = VecScale(p_scatter_wave_imag,scaling_pressure);
1268  // }
1269  ierr = post_proc.addFieldValuesPostProc("imPRES"); CHKERRQ(ierr);
1270  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1271  /* Particle Velocity The velocity of the
1272  movement of these material particles is called particle velocity */
1273  // ierr = post_proc.addFieldValuesGradientPostProc("rePRES","reVELOCITY"); CHKERRQ(ierr);
1274  // ierr = post_proc.addFieldValuesGradientPostProc("imPRES","imVELOCITY"); CHKERRQ(ierr);
1275 
1276  /* the potential and velocity of analytical solution fields */
1277  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
1278  ierr = post_proc.addFieldValuesPostProc("reEX"); CHKERRQ(ierr);
1279  ierr = post_proc.addFieldValuesPostProc("imEX"); CHKERRQ(ierr);
1280  // ierr = post_proc.addFieldValuesGradientPostProc("reEX","reVELOCITYex"); CHKERRQ(ierr);
1281  // ierr = post_proc.addFieldValuesGradientPostProc("imEX","imVELOCITYex"); CHKERRQ(ierr);
1282  }
1283 
1284  PetscBool reynolds_stress = PETSC_FALSE;
1285  ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-reynolds_stress",&reynolds_stress,NULL); CHKERRQ(ierr);
1286 
1287  if(reynolds_stress) {
1288 
1289  post_proc.getOpPtrVector().push_back(
1290  new ReynoldsStress(
1291  m_field,
1292  helmholtz_element,
1293  post_proc.postProcMesh,
1294  post_proc.mapGaussPts,
1295  "imPRES",
1296  "imPRES",
1297  post_proc.commonData));
1298 
1299  post_proc.getOpPtrVector().push_back(
1300  new ReynoldsStress(
1301  m_field,
1302  helmholtz_element,
1303  post_proc.postProcMesh,
1304  post_proc.mapGaussPts,
1305  "rePRES",
1306  "rePRES",
1307  post_proc.commonData));
1308 
1309  }
1310 
1311  struct OpPostProcOrder: public VolumeElementForcesAndSourcesCore::UserDataOperator {
1312 
1313  string FieldName;
1314  moab::Interface &postProcMesh;
1315  vector<EntityHandle> &mapGaussPts;
1316 
1317  OpPostProcOrder(
1318  const string& FieldName,moab::Interface &post_proc_mesh,vector<EntityHandle> &map_gauss_pts):VolumeElementForcesAndSourcesCore::UserDataOperator(FieldName,ForcesAndSurcesCore::UserDataOperator::OPROW),
1319  postProcMesh(post_proc_mesh),
1320  mapGaussPts(map_gauss_pts) { }
1321 
1322  PetscErrorCode doWork(
1323  int side,
1324  EntityType type,
1326 
1327  PetscFunctionBegin;
1329  int nb_gauss_pts = data.getN().size1();
1330 
1331  int order = getVolumeFE()->getMaxDataOrder();
1332  // cerr << "\n order = \n" << order << endl;
1333  Tag th_post_proc_order;
1334  int tag_length = 1;
1335  // double def_VAL[tag_length];
1336  // bzero(def_VAL,tag_length*sizeof(double));
1337  double def_VAL = 1;
1338  rval = postProcMesh.tag_get_handle(
1339  "ORDER_ADAPT",tag_length,MB_TYPE_INTEGER,th_post_proc_order,MB_TAG_CREAT|MB_TAG_SPARSE,&def_VAL
1340  ); CHKERRQ_MOAB(rval);
1341 
1342  for(int gg = 0;gg<nb_gauss_pts;gg++) {
1343  rval = postProcMesh.tag_set_data(
1344  th_post_proc_order,&mapGaussPts[gg],1,&order
1345  ); CHKERRQ_MOAB(rval);
1346  }
1347 
1348  PetscFunctionReturn(0);
1349  }
1350 
1351  };
1352 
1353  if(aDaptivity) {
1354  post_proc.getOpPtrVector().push_back(
1355  new OpPostProcOrder("rePRES",
1356  post_proc.postProcMesh,
1357  post_proc.mapGaussPts)
1358  );
1359  }
1360 
1361  ierr = m_field.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE",post_proc); CHKERRQ(ierr);
1362 
1363  // if(!pcomm->rank()) {
1364  ierr = post_proc.writeFile("fe_solution_mesh_post_proc.h5m"); CHKERRQ(ierr);
1365  // }
1366 
1367  }
1368 
1369 
1370  }
1371 
1372  // field split delete entities
1373  for(int i = 0;i!=2;i++) {
1374  if(nested_is_rows[i]) {
1375  ierr = ISDestroy(&nested_is_rows[i]); CHKERRQ(ierr);
1376  }
1377  for(int j = 0;j!=2;j++) {
1378  if(nested_matrices(i,j)) {
1379  ierr = MatDestroy(&nested_matrices(i,j)); CHKERRQ(ierr);
1380  }
1381  }
1382  }
1383  // field split delete entities
1384 
1385  // Destroy the KSP solvers
1386  ierr = VecDestroy(&F); CHKERRQ(ierr);
1387  ierr = VecDestroy(&T); CHKERRQ(ierr);
1388  /* FIXME KSP solver can not be destroyed when using field split method */
1389  // ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1390  ierr = PetscOptionsEnd(); CHKERRQ(ierr);
1391  ierr = PetscFinalize(); CHKERRQ(ierr);
1392 
1393  return 0;
1394 
1395 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:523
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 initializeProblem(MoFEM::Interface &m_field, string fe, string field, Range &tris, string nodals_positions="MESH_NODE_POSITIONS")
PetscErrorCode pressureInverseDft()
Aplly Inverse FFT to pressure degrees of freedom.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
moab::ErrorCode MoABErrorCode
MoAB error code.
Definition: Exceptions.hpp:66
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
PetscErrorCode generateReferenceElementMesh()
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 ...
virtual moab::Interface & get_moab()=0
CoreTmp< 0 > Core
Definition: Core.hpp:1129
Calculate the analytical solution of impinging wave on cylinder.
DEPRECATED MoFEMErrorCode ISCreateProblemFieldAndRank(const std::string &problem, RowColData rc, const std::string &field, int min_coeff_idx, int max_coeff_idx, IS *is, int verb=-1) const
create IS for given problem, field and rank range (collective)
Post-process fields on refined mesh.
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
auto post_proc
PetscErrorCode createPressureSeries(Vec T)
Create vectors for each wave lengths.
#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
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.
static char help[]
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.
PetscErrorCode pressureForwardDft()
Aplly Forward FFT to pressure degrees of freedom.
PetscErrorCode readData(const char *str, map< double, double > &series)
Read signal from text file.
PetscErrorCode destroyPressureSeries()
Destroy pressure series.
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.
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
DEPRECATED MoFEMErrorCode build_problem_on_distributed_mesh(const std::string &name, const bool square_matrix=true, int verb=-1)
build problem data structures, assuming that mesh is distributed (collective)
Calculate the analytical solution of plane wave guide propagating in direction theta.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
PetscErrorCode solve_problem(MoFEM::Interface &m_field, const string &problem_name, const string &fe_name, const string &re_field, const string &im_field, InsertMode mode, FUNEVAL &fun_evaluator, PetscBool is_partitioned)
const char * analytical_solution_types[]
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
int main(int argc, char *argv[])
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
Operators and data structures for the calculation of second order flux term employed to retrieve the ...
constexpr int order
static Index< 'i', 3 > i
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
CHKERRQ(ierr)
DataForcesAndSourcesCore::EntData EntData
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
static Index< 'j', 3 > j
PetscErrorCode solveProblem(MoFEM::Interface &m_field, string problem, string fe, DirichletBC &bc, Range &tris)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
DEPRECATED MoFEMErrorCode build_problem(const std::string &name, const bool square_matrix, int verb=-1)
build problem data structures
Post processing.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
DEPRECATED MoFEMErrorCode partition_problem(const std::string &name, int verb=-1)
partition problem dofs (collective)
DEPRECATED MoFEMErrorCode VecCreateGhost(const std::string &name, RowColData rc, Vec *V) const
create ghost vector for problem (collective)collective - need to be run on all processors in communic...
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
PetscErrorCode setApproxOps(MoFEM::Interface &m_field, string field_name, Range &tris, boost::shared_ptr< FUNEVAL > funtcion_evaluator, int field_number=0, string nodals_positions="MESH_NODE_POSITIONS")
#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
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
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.
PetscErrorCode forwardSpaceDft()
Apply Forward FFT and compose wave.
virtual bool check_field(const std::string &name) const =0
check if field is in database
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
static int debug
Core (interface) class.
Definition: Core.hpp:77
Operators and data structures for the calculation of momentum flux and Reynolds Stress of the complex...
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:177
PetscErrorCode setProblem(MoFEM::Interface &m_field, string problem)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
DEPRECATED MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)
PetscErrorCode solveForwardDFT(Mat A, Vec F, Vec T)
Solve problem for each wave number.
Analytical Dirichlet boundary conditions.
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)