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

Go to the source code of this file.

Functions

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

Variables

static int debug = 0
 
static char help [] = "...\n\n"
 

Function Documentation

◆ main()

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

Definition at line 79 of file fe_field_split.cpp.

79  {
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...
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
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)
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
VolumeElementForcesAndSourcesCoreBase * getVolumeFE() const
return pointer to Generic Volume Finite Element object
#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
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.
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
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.
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)
DEPRECATED MoFEMErrorCode partition_finite_elements(const std::string &name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=-1)
partition finite elements
constexpr int order
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
static Index< 'j', 3 > j
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.
#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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual bool check_field(const std::string &name) const =0
check if field is in database
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
static int debug
Core (interface) class.
Definition: Core.hpp:77
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:177
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
DEPRECATED MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)
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)
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
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...
int getMaxDataOrder() const
Get max order of approximation for data fields.
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)
DEPRECATED typedef ForcesAndSourcesCore ForcesAndSurcesCore

Variable Documentation

◆ debug

int debug = 0
static

Definition at line 76 of file fe_field_split.cpp.

◆ help

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

Definition at line 77 of file fe_field_split.cpp.