v0.10.0
Functions | Variables
fe_approximation.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 <TimeSeries.hpp>

Go to the source code of this file.

Functions

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

Variables

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

Detailed Description

Calculates finite element (Galerkin) approximation for incident wave problem.

Note:

In this implementation, first acoustic potential field is approximated on boundary and then finite element problem is solved. For more rigorous convergence study, trace of best approximations on boundary can be calculated and then finite element for domain and Neumann/mix boundary. That will give exact pollution error.

Definition in file fe_approximation.cpp.

Function Documentation

◆ main()

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

Definition at line 78 of file fe_approximation.cpp.

78  {
79 
80  ErrorCode rval;
81  PetscErrorCode ierr;
82 
83  PetscInitialize(&argc,&argv,(char *)0,help);
84 
85  //Core mb_instance;
86  moab::Core mb_instance;
87  moab::Interface& moab = mb_instance;
88  int rank;
89  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
90 
91  PetscBool flg = PETSC_TRUE;
92  char mesh_file_name[255];
93  ierr = PetscOptionsGetString(PETSC_NULL,PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
94  if(flg != PETSC_TRUE) {
95  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"*** ERROR -my_file (MESH FILE NEEDED)");
96  }
97 
98  // Create moab parallel communicator
99  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
100  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
101 
102  PetscBool is_partitioned = PETSC_FALSE;
103  ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
104  if(is_partitioned) {
105  //Read mesh to MOAB PARALLEL=READ_DELETE PARALLEL=READ_PART
106  const char *option;
107  // option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
108  option = "PARALLEL=READ_PART;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
109  // rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
110  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
111  } else {
112  const char *option;
113  option = "";
114  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
115  }
116 
117 
118  PetscBool is_labatto = PETSC_TRUE;
119  ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-lobatto",&is_labatto,NULL); CHKERRQ(ierr);
120 
121 
122  // Create MoFEM (cephas) database
123  MoFEM::Core core(moab);
124  MoFEM::Interface& m_field = core;
125 
126 
127  //set entitities bit level
128  BitRefLevel bit_level0;
129  bit_level0.set(0);
130  EntityHandle meshset_level0;
131  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
132  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
133 
134  //Fields
135  if(is_labatto) {
136  ierr = m_field.add_field("rePRES",H1,AINSWORTH_LOBATTO_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
137  ierr = m_field.add_field("imPRES",H1,AINSWORTH_LOBATTO_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
138  ierr = m_field.add_field("P",H1,AINSWORTH_LOBATTO_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr); // in time domain
139  } else {
140  ierr = m_field.add_field("rePRES",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
141  ierr = m_field.add_field("imPRES",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
142  ierr = m_field.add_field("P",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr); // in time domain
143  }
144 
145  //meshset consisting all entities in mesh
146  EntityHandle root_set = moab.get_root_set();
147  //add entities to field
148  ierr = m_field.add_ents_to_field_by_TETs(root_set,"rePRES"); CHKERRQ(ierr);
149  ierr = m_field.add_ents_to_field_by_TETs(root_set,"imPRES"); CHKERRQ(ierr);
150  ierr = m_field.add_ents_to_field_by_TETs(root_set,"P"); CHKERRQ(ierr);
151  /* FIXME add wave number A priori error estimator in here as a field, which
152  store the error on the elements of computational domain. */
153 
154  //set app. order
155  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
156  PetscInt order;
157  ierr = PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
158  if(flg != PETSC_TRUE) {
159  order = 2;
160  }
161 
162  PetscBool wavenumber_flg;
163  double wavenumber;
164  // set wave number from line command, that overwrite numbre form block set
165  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-wave_number",&wavenumber,&wavenumber_flg); CHKERRQ(ierr);
166  if(!wavenumber_flg) {
167 
168  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"wave number not given, set in line command -wave_number to fix problem");
169 
170  }
171 
172 
173  PetscBool aDaptivity = PETSC_FALSE;
174  ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-adaptivity",&aDaptivity,NULL); CHKERRQ(ierr);
175 
176 
177  if(aDaptivity) {
178  /* A priori error indicator */
179  /* Introduced by BĂ©riot from Efficient implementation of high-order finite
180  elements for Helmholtz problems */
181  PetscInt threeErrorlv = 1;
182  ierr = PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-error_level",&threeErrorlv,NULL); CHKERRQ(ierr);
183  if(threeErrorlv > 2) {
184  threeErrorlv = 2;
185  }
186 
187  /* 1 => l2, 2 => h1, 3 => anisworth */
188  PetscInt l2 = 1;
189  ierr = PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-error_type",&l2,NULL); CHKERRQ(ierr);
190  if(l2 > 3) {
191  l2 = 3;
192  }
193 
194  /* 1 => maximum edge length, 2 => average edge length, 3 => minimum edge length */
195  pair<PetscInt,PetscBool> isCriterion;
196  isCriterion.first = 1;
197  ierr = PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-error_criterion",&isCriterion.first,&isCriterion.second); CHKERRQ(ierr);
198  if(isCriterion.first > 3) {
199  isCriterion.first = 3;
200  }
201 
202  Range tets;
203  map<int,Range> orders_map;
204 
205  rval = m_field.get_moab().get_entities_by_type(0,MBTET,tets,false); CHKERRQ_MOAB(rval);
206  // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,"MAT_HELMHOLTZ",it)) {
207  // rval = mField.get_moab().get_entities_by_type(it->meshset,MBTET,tets,true); CHKERRQ_MOAB(rval);
208  // }
209 
210  for(Range::iterator tit = tets.begin();tit!=tets.end();tit++) {
211 
212  double sum = 0;
213  double eDges [6];
214  // std::vector<double> valarryEdge(6);
215  for(int ee = 0;ee<6;ee++) {
216 
217  EntityHandle edge;
218  rval = m_field.get_moab().side_element(*tit,1,ee,edge); CHKERRQ_MOAB(rval);
219  const EntityHandle* conn;
220  int num_nodes;
221  rval = m_field.get_moab().get_connectivity(edge,conn,num_nodes,true); CHKERRQ_MOAB(rval);
222 
223  double coords[num_nodes*3];
224  rval = m_field.get_moab().get_coords(conn,num_nodes,coords); CHKERRQ_MOAB(rval);
225  cblas_daxpy(3,-1,&coords[0],1,&coords[3],1);
226  // for(int ii = 0;ii<6;ii++) {
227  // cout << "\n number " << ii << endl;
228  // cout << " coords after cblas_daxpy = \n" << coords[ii] << endl;
229  // }
230  eDges[ee] = cblas_dnrm2(3,&coords[3],1);
231  sum += eDges[ee];
232  // valarryEdge[ee] = cblas_dnrm2(3,&coords[3],1);
233 
234  }
235 
236  /* more robust to select average length of edges of elements than
237  max or min */
238 
239  double averageEdgeK = (sum/6.) * wavenumber;
240  double maxEdge = *std::max_element(eDges,eDges+6);
241  double minEdge = *std::min_element(eDges,eDges+6);
242 
243  /* assign element size definition based on mesh quality */
244  // double qualityIndicator = 3.*minEdge/maxEdge;
245 
246  // if(qualityIndicator <= 1.0) {
247  // cerr << "\n qualityIndicator = \n" << qualityIndicator << endl;
248  // }
249  // if( (qualityIndicator < 1.5) || (qualityIndicator > 0.) ) {
250  // averageEdgeK = maxEdge * wavenumber;
251  // }
252 
253  if(isCriterion.first == 1) {
254  averageEdgeK = maxEdge * wavenumber;
255  } else if(isCriterion.first == 2) {
256  averageEdgeK = minEdge * wavenumber;;
257  }
258 
259  // cerr << "\n averageEdgeK = \n" << averageEdgeK << endl;
260 
261 
262  // Caulate order from table of 1D error test.
263  /* FIXME the order of p needs to be extend to more than 10 in future implementation */
264 
265  std::vector<double> errorLevel (10);
266 
267  double error_level [10];
268  if(threeErrorlv == 0) {
269  if(l2 == 1) {
270  /* target l2 error 15% */
271  /* error estimator in paper */
272  // 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);
273  /* error estimator 1 */
274  const double a[] = {0.5,2.900,4.8,7.5,7.97,12.33,12.54,14.6,15.2,19.43};
275  memcpy(error_level,a , sizeof error_level);
276  /* error estimator 3 */
277  // 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);
278  } else if(l2 == 2) {
279  /* error estimator h1 1 */
280  const double a[] = {0.2,2.3,3.88,6.8,7.1,11,11.17,14.54,15.28,19.023};
281  memcpy(error_level,a , sizeof error_level);
282  /* error estimator h1 3 */
283  // 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);
284  } else {
285  /* Anisworth error estimator */
286  // error_level [ ] = { };
287  }
288  errorLevel.insert (errorLevel.begin(),error_level,error_level+10);
289  // for(int ii = 0;ii<9;ii++) {
290  // cerr << " \n error_level[ii] = \n " << error_level[ii] << endl;
291  // }
292  } else if(threeErrorlv == 1) {
293 
294  if(l2 == 1) {
295  /* target l2 error 5% */
296  // 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);
297  const double a[] = {0.1,1.85,3.67,6.6,6.76,8.96,10.3,13.43,13.7,17};
298  memcpy(error_level,a , sizeof error_level);
299  // 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);
300  } else if(l2 == 2) {
301  const double a[] = {0.05,1.85,3.25,5.5,5.6,8.4,9.67,12.74,12.75,15.04};
302  memcpy(error_level,a , sizeof error_level);
303  // 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);
304  } else {
305  // memcpy(error_level, (double[]) { };
306  }
307  errorLevel.insert (errorLevel.begin(),error_level,error_level+10);
308 
309  } else if(threeErrorlv == 2) {
310 
311  if(l2 == 1) {
312  /* target l2 error 0.5% */
313  // 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);
314  const double a[] = {0.05,1,1.33,3.6,4.6,7.05,7.35,9.4,10.76,13.625};
315  memcpy(error_level,a , sizeof error_level);
316  // 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);
317  } else if(l2 == 2) {
318  const double a[] = {0.02,0.41,0.83,2.35,3.64,6.25,6.28,8.46,9.83,12.62};
319  memcpy(error_level,a , sizeof error_level);
320  // 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);
321  } else {
322  // memcpy(error_level, (double[]) { };
323  }
324  errorLevel.insert (errorLevel.begin(),error_level,error_level+10);
325 
326  }
327 
328  int orderP = 1;
329  for(int ii = 1;ii<8;ii++) {
330  if(errorLevel[ii] < averageEdgeK) {
331  orderP = ii + 1;
332  } else {
333  break;
334  }
335 
336  }
337 
338  // input a prior error indicator by ascending order with Range.
339  orders_map[orderP].insert(*tit);
340 
341  }
342 
343  map<int,Range>::iterator mit = orders_map.begin();
344  int max_value = mit->first;
345  int min_value = 10;
346  int avg_value = 0;
347 
348  for(;mit!=orders_map.end();mit++) {
349 
350  if(mit->first > max_value) {
351  max_value = mit->first;
352  }
353  if(min_value > mit->first) {
354  min_value = mit->first;
355  }
356  avg_value += mit->first;
357 
358 
359  Range tris;
360  rval = moab.get_adjacencies(mit->second,2,false,tris,moab::Interface::UNION);
361  Range edges;
362  rval = moab.get_adjacencies(mit->second,1,false,edges,moab::Interface::UNION);
363 
364  /* WARNING : the order of triangle and edges are consistent with tets in current implementation */
365  // ierr = m_field.set_field_order(mit->second,"rePRES",mit->first,2); CHKERRQ(ierr);
366  ierr = m_field.set_field_order(mit->second,"rePRES",mit->first); CHKERRQ(ierr);
367  ierr = m_field.set_field_order(tris,"rePRES",mit->first); CHKERRQ(ierr);
368  ierr = m_field.set_field_order(edges,"rePRES",mit->first); CHKERRQ(ierr);
369  //
370  ierr = m_field.set_field_order(mit->second,"imPRES",mit->first); CHKERRQ(ierr);
371  ierr = m_field.set_field_order(tris,"imPRES",mit->first); CHKERRQ(ierr);
372  ierr = m_field.set_field_order(edges,"imPRES",mit->first); CHKERRQ(ierr);
373 
374  ierr = m_field.set_field_order(mit->second,"P",mit->first); CHKERRQ(ierr);
375  ierr = m_field.set_field_order(tris,"P",mit->first); CHKERRQ(ierr);
376  ierr = m_field.set_field_order(edges,"P",mit->first); CHKERRQ(ierr);
377 
378  }
379 
380  if(!pcomm->rank()) {
381  printf("max_value_order_p %u\n",max_value);
382  printf("min_value_order_p %u\n",min_value);
383  printf("avg_value_order_p %u\n",int (avg_value/orders_map.size()));
384  }
385 
386  } else {
387 
388  ierr = m_field.set_field_order(root_set,MBTET,"rePRES",order); CHKERRQ(ierr);
389  ierr = m_field.set_field_order(root_set,MBTRI,"rePRES",order); CHKERRQ(ierr);
390  ierr = m_field.set_field_order(root_set,MBEDGE,"rePRES",order); CHKERRQ(ierr);
391 
392  ierr = m_field.set_field_order(root_set,MBTET,"imPRES",order); CHKERRQ(ierr);
393  ierr = m_field.set_field_order(root_set,MBTRI,"imPRES",order); CHKERRQ(ierr);
394  ierr = m_field.set_field_order(root_set,MBEDGE,"imPRES",order); CHKERRQ(ierr);
395 
396  ierr = m_field.set_field_order(root_set,MBTET,"P",order); CHKERRQ(ierr);
397  ierr = m_field.set_field_order(root_set,MBTRI,"P",order); CHKERRQ(ierr);
398  ierr = m_field.set_field_order(root_set,MBEDGE,"P",order); CHKERRQ(ierr);
399 
400  }
401 
402 
403  ierr = m_field.set_field_order(root_set,MBVERTEX,"rePRES",1); CHKERRQ(ierr);
404  ierr = m_field.set_field_order(root_set,MBVERTEX,"imPRES",1); CHKERRQ(ierr);
405  ierr = m_field.set_field_order(root_set,MBVERTEX,"P",1); CHKERRQ(ierr);
406 
407  if(!m_field.check_field("MESH_NODE_POSITIONS")) {
408 
409  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
410  ierr = m_field.add_ents_to_field_by_TETs(root_set,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
411  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
412  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
413  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
414  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
415 
416  ierr = m_field.build_fields(); CHKERRQ(ierr);
417  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
418  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
419 
420  } else {
421 
422  ierr = m_field.build_fields(); CHKERRQ(ierr);
423 
424  }
425 
426  // Finite Elements
427 
428  HelmholtzElement helmholtz_element(m_field);
429  ierr = helmholtz_element.getGlobalParametersFromLineCommandOptions(); CHKERRQ(ierr);
430  ierr = helmholtz_element.addHelmholtzElements("rePRES","imPRES"); CHKERRQ(ierr);
431  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
432  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","reEX"); CHKERRQ(ierr);
433  ierr = m_field.modify_finite_element_add_field_data("HELMHOLTZ_RERE_FE","imEX"); CHKERRQ(ierr);
434  }
435 
436  bool Dirichlet_bc_set = false;
437  Range bc_dirichlet_tris,analytical_bc_tris;
438  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"ANALYTICAL_BC",it)) {
439  rval = moab.get_entities_by_type(it->getMeshset(),MBTRI,analytical_bc_tris,true); CHKERRQ_MOAB(rval);
440  Dirichlet_bc_set = true;
441  }
442  bc_dirichlet_tris.merge(analytical_bc_tris);
443  AnalyticalDirichletHelmholtzBC analytical_bc_real(m_field);
444  AnalyticalDirichletHelmholtzBC analytical_bc_imag(m_field);
445  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",analytical_bc_tris); CHKERRQ(ierr);
446  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",analytical_bc_tris); CHKERRQ(ierr);
447 
448 
449  double aTtenuation = 0;
450  // set loss tangent from line command, the complex wave number to include attenuation.
451  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-attenuation",&aTtenuation,NULL); CHKERRQ(ierr);
452 
453  ierr = PetscOptionsBegin(m_field.get_comm(),NULL,"Helmholtz problem options","none"); CHKERRQ(ierr);
454  PetscBool isRayleigh = PETSC_FALSE;
455  ierr = PetscOptionsBool(
456  "-rayleigh_wave",
457  "If true the incident is Rayleigh wave, otherwise, it is an non-attenuated wave","",
458  PETSC_FALSE,
459  &isRayleigh,
460  NULL
461  ); CHKERRQ(ierr);
462 
463  double power_of_incident_wave = 1;
464  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-amplitude_of_incident_wave",&power_of_incident_wave,NULL); CHKERRQ(ierr);
465 
466 
467  // This is added for a case than on some surface, defined by the user a
468  // incident plane wave is scattered.
469  map<int,PlaneIncidentWaveSacttrerData> planeWaveScatterData;
470  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"SOFT_INCIDENT_WAVE_BC",it)) {
471  rval = moab.get_entities_by_type(it->getMeshset(),MBTRI,planeWaveScatterData[it->getMeshsetId()].tRis,true); CHKERRQ_MOAB(rval);
472  ierr = analytical_bc_real.initializeProblem(m_field,"BCREAL_FE","rePRES",planeWaveScatterData[it->getMeshsetId()].tRis); CHKERRQ(ierr);
473  ierr = analytical_bc_imag.initializeProblem(m_field,"BCIMAG_FE","imPRES",planeWaveScatterData[it->getMeshsetId()].tRis); CHKERRQ(ierr);
474  bc_dirichlet_tris.merge(planeWaveScatterData[it->getMeshsetId()].tRis);
475 
476  Dirichlet_bc_set = true;
477 
478  }
479 
480  //ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
481  //// Build adjacencies
482  //ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
483 
484  // Problem
485  ierr = m_field.add_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
486  ierr = m_field.add_problem("BCREAL_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for real field
487  ierr = m_field.add_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr); //analytical Dirichlet for imag field
488 
489  // Set refinement level for problem
490  ierr = m_field.modify_problem_ref_level_add_bit("ACOUSTIC_PROBLEM",bit_level0); CHKERRQ(ierr);
491  ierr = m_field.modify_problem_ref_level_add_bit("BCREAL_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
492  ierr = m_field.modify_problem_ref_level_add_bit("BCIMAG_PROBLEM",bit_level0); CHKERRQ(ierr); //analytical Dirichlet
493 
494  // Add elements to problems
495  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
496  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_IMIM_FE"); CHKERRQ(ierr);
497  ierr = m_field.modify_problem_add_finite_element("ACOUSTIC_PROBLEM","HELMHOLTZ_REIM_FE"); CHKERRQ(ierr);
498 
499  ierr = m_field.modify_problem_add_finite_element("BCREAL_PROBLEM","BCREAL_FE"); CHKERRQ(ierr);
500  ierr = m_field.modify_problem_add_finite_element("BCIMAG_PROBLEM","BCIMAG_FE"); CHKERRQ(ierr);
501 
502  // Get start time for analyse
503  PetscLogDouble t1,t2;
504  PetscLogDouble v1,v2;
505  ierr = PetscTime(&v1); CHKERRQ(ierr);
506  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
507  PetscLogDouble m1;
508  PetscLogDouble m2;
509  ierr = PetscMemoryGetCurrentUsage(&m1); CHKERRQ(ierr);
510 
511  // Build problems
512  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
513  // Build adjacencies
514  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
515 
516  // build porblems
517  if(is_partitioned) {
518  // if mesh is partitioned
519 
520  ierr = m_field.build_problem_on_distributed_mesh("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
521  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
522 
523  if(Dirichlet_bc_set) {
524 
525  ierr = m_field.build_problem_on_distributed_mesh("BCREAL_PROBLEM",true); CHKERRQ(ierr);
526  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM",true); CHKERRQ(ierr);
527 
528  ierr = m_field.build_problem_on_distributed_mesh("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
529  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
530 
531  }
532 
533 
534  } else {
535  // if not partitioned mesh is load to all processes
536 
537  ierr = m_field.build_problem("ACOUSTIC_PROBLEM",true); CHKERRQ(ierr);
538  ierr = m_field.partition_problem("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
539  ierr = m_field.partition_finite_elements("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
540 
541  if(Dirichlet_bc_set) {
542  ierr = m_field.build_problem("BCREAL_PROBLEM",true); CHKERRQ(ierr);
543  ierr = m_field.partition_problem("BCREAL_PROBLEM"); CHKERRQ(ierr);
544  ierr = m_field.partition_finite_elements("BCREAL_PROBLEM"); CHKERRQ(ierr);
545 
546 
547 
548  ierr = m_field.build_problem("BCIMAG_PROBLEM",true); CHKERRQ(ierr);
549  ierr = m_field.partition_problem("BCIMAG_PROBLEM"); CHKERRQ(ierr);
550  ierr = m_field.partition_finite_elements("BCIMAG_PROBLEM"); CHKERRQ(ierr);
551  }
552  }
553 
554  ierr = m_field.partition_ghost_dofs("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
555  if(Dirichlet_bc_set) {
556  ierr = m_field.partition_ghost_dofs("BCREAL_PROBLEM"); CHKERRQ(ierr);
557  ierr = m_field.partition_ghost_dofs("BCIMAG_PROBLEM"); CHKERRQ(ierr);
558  }
559 
560  // Get problem matrices and vectors
561  Vec F; //Right hand side vector
562  ierr = m_field.VecCreateGhost("ACOUSTIC_PROBLEM",ROW,&F); CHKERRQ(ierr);
563  Vec T; //Solution vector
564  ierr = VecDuplicate(F,&T); CHKERRQ(ierr);
565  Mat A; //Left hand side matrix
566  ierr = m_field.MatCreateMPIAIJWithArrays("ACOUSTIC_PROBLEM",&A); CHKERRQ(ierr);
567  /* change matrix type to suit ilu */
568  // Mat B;
569  // ierr = MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&A);
570  // ierr = MatDestroy(&B); CHKERRQ(ierr);
571  // ierr = MatSetType(A,MATBAIJ);
572  // ierr = m_field.MatCreateMPIBAIJWithArrays("ACOUSTIC_PROBLEM",&A); CHKERRQ(ierr);
573  /* change matrix type to suit ilu end */
574  ierr = helmholtz_element.setOperators(A,F,"rePRES","imPRES"); CHKERRQ(ierr);
575 
576 
577  //wave direction unit vector=[x,y,z]^T
578  VectorDouble wave_direction;
579  wave_direction.resize(3);
580  wave_direction.clear();
581  wave_direction[2] = 1; // default:X direction [0,0,1]
582 
583  int nmax = 3;
584  ierr = PetscOptionsGetRealArray(PETSC_NULL,PETSC_NULL,"-wave_direction",&wave_direction[0],&nmax,NULL); CHKERRQ(ierr);
585  if(nmax > 0 && nmax != 3) {
586  SETERRQ(PETSC_COMM_SELF,MOFEM_INVALID_DATA,"*** ERROR -wave_direction [3*1 vector] default:X direction [0,0,1]");
587  }
588 
589  PetscInt choise_value = NO_ANALYTICAL_SOLUTION;
590  // set type of analytical solution
591  ierr = PetscOptionsGetEList(PETSC_NULL,NULL,"-analytical_solution_type",analytical_solution_types,6,&choise_value,NULL); CHKERRQ(ierr);
592 
593  switch((AnalyticalSolutionTypes)choise_value) {
594 
596 
597  {
598  double scattering_sphere_radius = 0.5;
599  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-scattering_sphere_radius",&scattering_sphere_radius,NULL); CHKERRQ(ierr);
600 
601 
602 
603  boost::shared_ptr<HardSphereScatterWave> function_evaluator = boost::shared_ptr<HardSphereScatterWave>(
604  new HardSphereScatterWave(wavenumber,scattering_sphere_radius)
605  );
606  ierr = analytical_bc_real.setApproxOps(
607  m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL
608  ); CHKERRQ(ierr);
609  ierr = analytical_bc_imag.setApproxOps(
610  m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG
611  ); CHKERRQ(ierr);
612  Dirichlet_bc_set = true;
613 
614  }
615 
616  break;
617 
619 
620  {
621 
622  double scattering_sphere_radius = 0.5;
623  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-scattering_sphere_radius",&scattering_sphere_radius,NULL); CHKERRQ(ierr);
624 
625 
626  boost::shared_ptr<SoftSphereScatterWave> function_evaluator = boost::shared_ptr<SoftSphereScatterWave>(new SoftSphereScatterWave(wavenumber,scattering_sphere_radius));
627  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
628  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
629  Dirichlet_bc_set = true;
630 
631  }
632 
633  break;
634 
635  case PLANE_WAVE:
636 
637  {
638 
639  double angle = 0.25;
640  // set wave number from line command, that overwrite numbre form block set
641  ierr = PetscOptionsGetScalar(PETSC_NULL,NULL,"-wave_guide_angle",&angle,NULL); CHKERRQ(ierr);
642 
643 
644  boost::shared_ptr<PlaneWave> function_evaluator = boost::shared_ptr<PlaneWave>(new PlaneWave(wavenumber,angle*M_PI));
645  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
646  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
647  Dirichlet_bc_set = true;
648 
649  }
650 
651  break;
652 
654 
655  {
656 
657  boost::shared_ptr<HardCylinderScatterWave> function_evaluator = boost::shared_ptr<HardCylinderScatterWave>(new HardCylinderScatterWave(wavenumber));
658  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
659  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
660  Dirichlet_bc_set = true;
661 
662 
663  }
664 
665  break;
666 
668 
669  {
670 
671  boost::shared_ptr<SoftCylinderScatterWave> function_evaluator = boost::shared_ptr<SoftCylinderScatterWave>(new SoftCylinderScatterWave(wavenumber));
672  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
673  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
674  Dirichlet_bc_set = true;
675 
676  }
677 
678  break;
679 
680  case INCIDENT_WAVE:
681 
682  {
683 
684  boost::shared_ptr<IncidentWave> function_evaluator =
685  boost::shared_ptr<IncidentWave>(new IncidentWave(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
686  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,power_of_incident_wave));
687  ierr = analytical_bc_real.setApproxOps(m_field,"rePRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::REAL); CHKERRQ(ierr);
688  ierr = analytical_bc_imag.setApproxOps(m_field,"imPRES",analytical_bc_tris,function_evaluator,GenericAnalyticalSolution::IMAG); CHKERRQ(ierr);
689  Dirichlet_bc_set = true;
690 
691  }
692 
693  break;
694 
696 
697  {
698  // Dirichlet_bc_set = false;
699  }
700 
701  break;
702 
703  }
704 
705  // Analytical boundary conditions
706  AnalyticalDirichletHelmholtzBC::DirichletBC analytical_ditihlet_bc_real(m_field,"rePRES",A,T,F);
707  AnalyticalDirichletHelmholtzBC::DirichletBC analytical_ditihlet_bc_imag(m_field,"imPRES",A,T,F);
708 
709 
710  PetscBool monochromatic_wave = PETSC_TRUE;
711  ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-monochromatic_wave",&monochromatic_wave,NULL); CHKERRQ(ierr);
712 
713  PetscBool add_incident_wave = PETSC_FALSE;
714  ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-add_incident_wave",&add_incident_wave,NULL); CHKERRQ(ierr);
715  if(add_incident_wave) {
716 
717  // define problem
718  ierr = m_field.add_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
719  // set finite elements for problem
720  ierr = m_field.modify_problem_add_finite_element("INCIDENT_WAVE","HELMHOLTZ_RERE_FE"); CHKERRQ(ierr);
721  // set refinment level for problem
722  ierr = m_field.modify_problem_ref_level_add_bit("INCIDENT_WAVE",bit_level0); CHKERRQ(ierr);
723 
724  // build porblems
725  if(is_partitioned) {
726  // if mesh is partitioned
727  ierr = m_field.build_problem_on_distributed_mesh("INCIDENT_WAVE",true); CHKERRQ(ierr);
728  ierr = m_field.partition_finite_elements("INCIDENT_WAVE",true); CHKERRQ(ierr);
729  } else {
730  // if not partitioned mesh is load to all processes
731  ierr = m_field.build_problem("INCIDENT_WAVE",true); CHKERRQ(ierr);
732  ierr = m_field.partition_problem("INCIDENT_WAVE"); CHKERRQ(ierr);
733  ierr = m_field.partition_finite_elements("INCIDENT_WAVE"); CHKERRQ(ierr);
734  }
735  ierr = m_field.partition_ghost_dofs("INCIDENT_WAVE"); CHKERRQ(ierr);
736 
737  }
738 
739 
740  KSP solver;
741  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
742 
743  if(monochromatic_wave) {
744 
745  if(Dirichlet_bc_set) {
746 
747  {
748 
749  map<int,PlaneIncidentWaveSacttrerData>::iterator mit = planeWaveScatterData.begin();
750 
751  boost::shared_ptr<IncidentWave> function_evaluator;
752  boost::shared_ptr<IncidentWavePointSource> functionEvaluator;
753  for(;mit!=planeWaveScatterData.end();mit++) {
754 
755  if(!helmholtz_element.globalParameters.sourceCoordinate.second) {
756  if(helmholtz_element.globalParameters.isRadiation.first) {
757  // note negative field, scatter field should cancel incident wave
758  function_evaluator = boost::shared_ptr<IncidentWave>(
759  new IncidentWave(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
760  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,-helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
761  );
762  } else {
763  function_evaluator = boost::shared_ptr<IncidentWave>(
764  new IncidentWave(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
765  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
766  );
767  }
768  ierr = analytical_bc_real.setApproxOps(
769  m_field,"rePRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::REAL
770  ); CHKERRQ(ierr);
771  ierr = analytical_bc_imag.setApproxOps(
772  m_field,"imPRES",mit->second.tRis,function_evaluator,GenericAnalyticalSolution::IMAG
773  ); CHKERRQ(ierr);
774  } else {
775  if(helmholtz_element.globalParameters.isRadiation.first) {
776  // note negative field, scatter field should cancel incident wave
777  functionEvaluator = boost::shared_ptr<IncidentWavePointSource>(
778  new IncidentWavePointSource(wavenumber,helmholtz_element.globalParameters.sourceCoordinate.first,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
779  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,-helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
780  );
781  } else {
782  functionEvaluator = boost::shared_ptr<IncidentWavePointSource>(
783  new IncidentWavePointSource(wavenumber,helmholtz_element.globalParameters.sourceCoordinate.first,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
784  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first)
785  );
786  }
787  ierr = analytical_bc_real.setApproxOps(
788  m_field,"rePRES",mit->second.tRis,functionEvaluator,GenericAnalyticalSolution::REAL
789  ); CHKERRQ(ierr);
790  ierr = analytical_bc_imag.setApproxOps(
791  m_field,"imPRES",mit->second.tRis,functionEvaluator,GenericAnalyticalSolution::IMAG
792  ); CHKERRQ(ierr);
793  }
794 
795  }
796 
797  }
798  // Solve for analytical Dirichlet bc dofs
799  ierr = analytical_bc_real.setProblem(m_field,"BCREAL_PROBLEM"); CHKERRQ(ierr);
800  ierr = analytical_bc_imag.setProblem(m_field,"BCIMAG_PROBLEM"); CHKERRQ(ierr);
801  ierr = analytical_bc_real.solveProblem(
802  m_field,"BCREAL_PROBLEM","BCREAL_FE",analytical_ditihlet_bc_real,bc_dirichlet_tris
803  ); CHKERRQ(ierr);
804  ierr = analytical_bc_imag.solveProblem(
805  m_field,"BCIMAG_PROBLEM","BCIMAG_FE",analytical_ditihlet_bc_imag,bc_dirichlet_tris
806  ); CHKERRQ(ierr);
807 
808  ierr = analytical_bc_real.destroyProblem(); CHKERRQ(ierr);
809  ierr = analytical_bc_imag.destroyProblem(); CHKERRQ(ierr);
810 
811  }
812 
813  /* calculate the interior problem to set the fulid velocity */
814  if(helmholtz_element.globalParameters.fRequency.second) {
815  // const complex< double > i( 0.0, 1.0 );
816  complex< double > wave_number = (2*M_PI*helmholtz_element.globalParameters.fRequency.first / helmholtz_element.globalParameters.vElocity.first);
817  // FIXME Why it is not working with first method ?
818  // map<int,helmholtz_element::VolumeData>::iterator vit = helmholtz_element.volumeData.begin();
819  //
820  // for(;vit != helmholtz_element.volumeData.end();vit++) {
821  // vit->second.waveNumber = wave_number;
822  // }
823 
825  if(it->getName().compare(0,13,"MAT_HELMHOLTZ") == 0) {
826  helmholtz_element.volumeData[it->getMeshsetId()].waveNumber = wave_number;
827  }
828  }
829 
830  }
831 
832  // Zero vectors
833  ierr = VecZeroEntries(T); CHKERRQ(ierr);
834  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
835  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
836  ierr = VecZeroEntries(F); CHKERRQ(ierr);
837  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
838  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
839  ierr = MatZeroEntries(A); CHKERRQ(ierr);
840 
841  // Assemble problem
842  if(Dirichlet_bc_set) {
843  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
844  ierr = m_field.problem_basic_method_preProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
845  }
846 
847  ierr = helmholtz_element.calculateA("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
848  ierr = helmholtz_element.calculateF("ACOUSTIC_PROBLEM"); CHKERRQ(ierr);
849  // ierr = MatView(A,PETSC_VIEWER_DRAW_WORLD);
850  // std::string wait;
851  // std::cin >> wait;
852 
853  if(Dirichlet_bc_set) {
854  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_real); CHKERRQ(ierr);
855  ierr = m_field.problem_basic_method_postProcess("ACOUSTIC_PROBLEM",analytical_ditihlet_bc_imag); CHKERRQ(ierr);
856  }
857 
858  // FIXME Point Source Right hand sides for vertex source for Total Acoustic potential
859  // \Phi = \Phi_{S} + \Phi_{I}
860  // const Problem *problem_ptr;
861  // ierr = m_field.get_problem("ACOUSTIC_PROBLEM",&problem_ptr); CHKERRQ(ierr);
862  // for(_IT_NUMEREDDOF_ROW_BY_NAME_ENT_PART_FOR_LOOP_(problem_ptr,"rePRES",node,m_field.getCommRank(),dit)) {
863  // int global_idx = dit-> getPetscGlobalDofIdx();
864  // cerr << "\n global_idx = \n\n " << global_idx << endl;
865  // ierr = VecSetValue(F,global_idx,helmholtz_element.globalParameters.amplitudeOfIncidentWaveReal.first,INSERT_VALUES); CHKERRQ(ierr);
866  // }
867 
868 
869  ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
870  ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
871  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
872  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
873 
874  // if(monochromatic_wave) {
875  /* get the sparsity information of matrix */
876  MatInfo info;
877  // double nonZeroes;
878  ierr = MatGetInfo(A,MAT_GLOBAL_SUM,&info);
879  // nonZeroes = info.nz_allocated;
880  // MatGetSize(A,&m,&n);
881  if(!pcomm->rank()) {
882  printf("info.nz_used %g\n",info.nz_used);
883  printf("info.memory %g\n",info.memory/2073741824);
884  // printf("info.factor_mallocs %g\n",info.factor_mallocs);
885  // printf("info.assemblies %g\n",info.assemblies);
886  }
887  // }
888  // ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
889 
890  ierr = VecScale(F,-1); CHKERRQ(ierr);
891  // Solve problem
892  ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
893  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
894  ierr = KSPSetUp(solver); CHKERRQ(ierr);
895 
896  // /* flops count */
897  // PetscLogEvent USER_EVENT;
898  // ierr = PetscLogEventRegister("User event",0,&USER_EVENT); CHKERRQ(ierr);
899  // ierr = PetscLogEventBegin(USER_EVENT,solver,0,0,0); CHKERRQ(ierr);
900  // /* flops count */
901 
902  ierr = KSPSolve(solver,F,T); CHKERRQ(ierr);
903  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
904  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
905 
906  // /* flops count */
907  // PetscScalar flops,user_flops;
908  // // ierr = PetscGetFlops(&flops); CHKERRQ(ierr);
909  // ierr = PetscLogFlops(user_flops); CHKERRQ(ierr);
910  // ierr = PetscLogEventEnd(USER_EVENT,solver,0,0,0); CHKERRQ(ierr);
911  // // PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total flops %d \n",pcomm->rank(),flops);
912  // PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Log flops %d \n",pcomm->rank(),user_flops);
913  // /* flops count */
914 
915 
916 
917  //Save data on mesh
918  if(is_partitioned) {
919 
920  // no need for global communication
921  ierr = m_field.set_local_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
922 
923  } else {
924 
925  ierr = m_field.set_global_ghost_vector("ACOUSTIC_PROBLEM",ROW,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
926 
927  }
928 
929  } else {
930 
931  // define problem
932  ierr = m_field.add_problem("PRESSURE_IN_TIME"); CHKERRQ(ierr);
933  // set finite elements for problem
934  ierr = m_field.modify_problem_add_finite_element("PRESSURE_IN_TIME","PRESSURE_FE"); CHKERRQ(ierr);
935  // set refinment level for problem
936  ierr = m_field.modify_problem_ref_level_add_bit("PRESSURE_IN_TIME",bit_level0); CHKERRQ(ierr);
937 
938  // build porblems
939  if(is_partitioned) {
940  // if mesh is partitioned
941  ierr = m_field.build_problem_on_distributed_mesh("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
942  ierr = m_field.partition_finite_elements("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
943  } else {
944  // if not partitioned mesh is load to all processes
945  ierr = m_field.build_problem("PRESSURE_IN_TIME",true); CHKERRQ(ierr);
946  ierr = m_field.partition_problem("PRESSURE_IN_TIME"); CHKERRQ(ierr);
947  ierr = m_field.partition_finite_elements("PRESSURE_IN_TIME"); CHKERRQ(ierr);
948  }
949  ierr = m_field.partition_ghost_dofs("PRESSURE_IN_TIME"); CHKERRQ(ierr);
950 
951  TimeSeries time_series(m_field,helmholtz_element,
952  analytical_bc_real,analytical_bc_imag,
953  analytical_ditihlet_bc_real,analytical_ditihlet_bc_imag,
954  Dirichlet_bc_set,bc_dirichlet_tris,planeWaveScatterData);
955 
956  ierr = time_series.readData(); CHKERRQ(ierr);
957  ierr = time_series.createPressureSeries(T); CHKERRQ(ierr);
958  ierr = time_series.forwardSpaceDft(); CHKERRQ(ierr);
959  // ierr = time_series.saveFrequencyData(); // save the frequency data to file
960  ierr = time_series.solveForwardDFT(solver,A,F,T); CHKERRQ(ierr);
961  ierr = time_series.pressureForwardDft(); CHKERRQ(ierr);
962  ierr = time_series.pressureInverseDft(); CHKERRQ(ierr);
963  ierr = time_series.generateReferenceElementMesh(); CHKERRQ(ierr);
964  ierr = time_series.saveResults(); CHKERRQ(ierr);
965  ierr = time_series.destroyPressureSeries(); CHKERRQ(ierr);
966 
967  }
968  // PetscInt its;
969  // ierr = KSPGetIterationNumber(solver,&its);
970  ierr = PetscTime(&v2);CHKERRQ(ierr);
971  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
972  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f S CPU Time = %f S \n",pcomm->rank(),v2-v1,t2-t1);
973  ierr = PetscMemoryGetCurrentUsage(&m2); CHKERRQ(ierr);
974  /* PetscMemoryGetMaximumUsage */
975  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Memory usage = %f GBytes \n",pcomm->rank(),(m2-m1)/1073741824);
976 
977 
978 
979  //Vec P,M;
980  //ierr = m_field.VecCreateGhost("EX1_PROBLEM",COL,&M); CHKERRQ(ierr);
981  //ierr = VecDuplicate(M,&P); CHKERRQ(ierr);
982 
983  //ierr = m_field.set_local_ghost_vector("EX1_PROBLEM",COL,M,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
984  //ierr = m_field.set_other_global_ghost_vector("EX1_PROBLEM","reEX","imEX",COL,P,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
985 
986  //double nrm2_M;
987  //ierr = VecNorm(M,NORM_2,&nrm2_M); CHKERRQ(ierr);
988 
989  //Vec V;
990  //ierr = m_field.VecCreateGhost("ACOUSTIC_PROBLEM",COL,&V); CHKERRQ(ierr);
991  //ierr = m_field.set_local_ghost_vector("ACOUSTIC_PROBLEM",COL,V,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
992 
993  //VecScatter scatter_real,scatter_imag;
994 
995  //ierr = m_field.VecScatterCreate(V,"ACOUSTIC_PROBLEM","rePRES",COL,M,"EX1_PROBLEM","reEX",COL,&scatter_real); CHKERRQ(ierr);
996 
997  //ierr = m_field.VecScatterCreate(V,"ACOUSTIC_PROBLEM","imPRES",COL,P,"EX1_PROBLEM","reEX",COL,&scatter_imag); CHKERRQ(ierr);
998 
999  //VecScale(V,-1);
1000 
1001  //ierr = VecScatterBegin(scatter_real,V,M,ADD_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1002  //ierr = VecScatterEnd(scatter_real,V,M,ADD_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1003 
1004  //double nrm2_ErrM;
1005  //ierr = VecNorm(M,NORM_2,&nrm2_ErrM); CHKERRQ(ierr);
1006  //PetscPrintf(PETSC_COMM_WORLD,"L2 relative error on real field of acoustic problem %6.4e\n",(nrm2_ErrM)/(nrm2_M));
1007 
1008 
1009  //ierr = VecDestroy(&M); CHKERRQ(ierr);
1010  //ierr = VecDestroy(&P); CHKERRQ(ierr);
1011  //ierr = VecDestroy(&V); CHKERRQ(ierr);
1012 
1013 
1014 
1015  if(monochromatic_wave) {
1016 
1017  if(add_incident_wave) {
1018  /* FIXME problem with HELMHOLTZ_RERE_FE since this element contains 2 fields,and
1019  create big matrix did not consistent with field_approximation.hpp calculation method. */
1020 
1021  if(helmholtz_element.globalParameters.sourceCoordinate.second) {
1022  IncidentWavePointSource function_evaluator(wavenumber,helmholtz_element.globalParameters.sourceCoordinate.first,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
1023  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,power_of_incident_wave);
1024  ierr = solve_problem(m_field,"INCIDENT_WAVE","HELMHOLTZ_RERE_FE","rePRES","imPRES",
1025  ADD_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
1026  } else {
1027  if(helmholtz_element.globalParameters.complexWaveNumber.second) {
1028  wave_direction = helmholtz_element.globalParameters.waveOscilationDirection.first;
1029  }
1030  IncidentWave function_evaluator(wavenumber,wave_direction,isRayleigh,aTtenuation,helmholtz_element.globalParameters.rAdius.first,helmholtz_element.globalParameters.complexWaveNumber.first,
1031  helmholtz_element.globalParameters.fRequency.first,helmholtz_element.globalParameters.vElocity.first,helmholtz_element.globalParameters.transmissionCoefficient.first,power_of_incident_wave);
1032  ierr = solve_problem(m_field,"INCIDENT_WAVE","HELMHOLTZ_RERE_FE","rePRES","imPRES",
1033  ADD_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
1034  }
1035 
1036  }
1037 
1038  if(is_partitioned) {
1039 
1040  rval = moab.write_file("fe_solution.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERRQ_MOAB(rval);
1041 
1042  } else {
1043 
1044  if(!pcomm->rank()) {
1045  rval = moab.write_file("fe_solution.h5m"); CHKERRQ_MOAB(rval);
1046  }
1047 
1048  }
1049 
1050 
1051 
1052  PetscBool save_postproc_mesh = PETSC_TRUE;
1053  ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-save_postproc_mesh",&save_postproc_mesh,NULL); CHKERRQ(ierr);
1054  if(save_postproc_mesh) {
1055 
1057  // PetscBool pressure_field = PETSC_FALSE;
1058  /* FIXME post-process the pressure field from acoustic potential */
1059  // ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-pressure_field",&pressure_field,NULL); CHKERRQ(ierr);
1060  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
1061  // PetscScalar scaling_pressure;
1062  // if(pressure_field) {
1063  // scaling_pressure = - 2.0 * M_PI * helmholtzElement.globalParameters.fRequency.first * helmholtzElement.globalParameters.dEnsity.first;
1064  // ierr = VecScale(p_scatter_wave_real,scaling_pressure);
1065  // }
1066  ierr = post_proc.addFieldValuesPostProc("rePRES"); CHKERRQ(ierr);
1067  // if(pressure_field) {
1068  // scaling_pressure *= - 1.0;
1069  // ierr = VecScale(p_scatter_wave_imag,scaling_pressure);
1070  // }
1071  ierr = post_proc.addFieldValuesPostProc("imPRES"); CHKERRQ(ierr);
1072  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1073  /* Particle Velocity The velocity of the
1074  movement of these material particles is called particle velocity */
1075  ierr = post_proc.addFieldValuesGradientPostProc("rePRES","reVELOCITY"); CHKERRQ(ierr);
1076  ierr = post_proc.addFieldValuesGradientPostProc("imPRES","imVELOCITY"); CHKERRQ(ierr);
1077 
1078  /* the potential and velocity of analytical solution fields */
1079  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
1080  ierr = post_proc.addFieldValuesPostProc("reEX"); CHKERRQ(ierr);
1081  ierr = post_proc.addFieldValuesPostProc("imEX"); CHKERRQ(ierr);
1082  // ierr = post_proc.addFieldValuesGradientPostProc("reEX","reVELOCITYex"); CHKERRQ(ierr);
1083  // ierr = post_proc.addFieldValuesGradientPostProc("imEX","imVELOCITYex"); CHKERRQ(ierr);
1084  }
1085 
1086  PetscBool reynolds_stress = PETSC_FALSE;
1087  ierr = PetscOptionsGetBool(PETSC_NULL,NULL,"-reynolds_stress",&reynolds_stress,NULL); CHKERRQ(ierr);
1088 
1089  if(reynolds_stress) {
1090 
1091  post_proc.getOpPtrVector().push_back(
1092  new ReynoldsStress(
1093  m_field,
1094  helmholtz_element,
1095  post_proc.postProcMesh,
1096  post_proc.mapGaussPts,
1097  "imPRES",
1098  "imPRES",
1099  post_proc.commonData));
1100 
1101  post_proc.getOpPtrVector().push_back(
1102  new ReynoldsStress(
1103  m_field,
1104  helmholtz_element,
1105  post_proc.postProcMesh,
1106  post_proc.mapGaussPts,
1107  "rePRES",
1108  "rePRES",
1109  post_proc.commonData));
1110 
1111  }
1112 
1113  struct OpPostProcOrder: public VolumeElementForcesAndSourcesCore::UserDataOperator {
1114 
1115  string FieldName;
1116  moab::Interface &postProcMesh;
1117  vector<EntityHandle> &mapGaussPts;
1118 
1119  OpPostProcOrder(
1120  const string& FieldName,moab::Interface &post_proc_mesh,vector<EntityHandle> &map_gauss_pts):VolumeElementForcesAndSourcesCore::UserDataOperator(FieldName,ForcesAndSurcesCore::UserDataOperator::OPROW),
1121  postProcMesh(post_proc_mesh),
1122  mapGaussPts(map_gauss_pts) { }
1123 
1124  PetscErrorCode doWork(
1125  int side,
1126  EntityType type,
1128 
1129  PetscFunctionBegin;
1131  int nb_gauss_pts = data.getN().size1();
1132 
1133  int order = getVolumeFE()->getMaxDataOrder();
1134  // cerr << "\n order = \n" << order << endl;
1135  Tag th_post_proc_order;
1136  int tag_length = 1;
1137  // double def_VAL[tag_length];
1138  // bzero(def_VAL,tag_length*sizeof(double));
1139  double def_VAL = 1;
1140  rval = postProcMesh.tag_get_handle(
1141  "ORDER_ADAPT",tag_length,MB_TYPE_INTEGER,th_post_proc_order,MB_TAG_CREAT|MB_TAG_SPARSE,&def_VAL
1142  ); CHKERRQ_MOAB(rval);
1143 
1144  for(int gg = 0;gg<nb_gauss_pts;gg++) {
1145  rval = postProcMesh.tag_set_data(
1146  th_post_proc_order,&mapGaussPts[gg],1,&order
1147  ); CHKERRQ_MOAB(rval);
1148  }
1149 
1150  PetscFunctionReturn(0);
1151  }
1152 
1153  };
1154 
1155  if(aDaptivity) {
1156  post_proc.getOpPtrVector().push_back(
1157  new OpPostProcOrder("rePRES",
1158  post_proc.postProcMesh,
1159  post_proc.mapGaussPts)
1160  );
1161  }
1162  // ierr = m_field.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_IMIM_FE",post_proc); CHKERRQ(ierr);
1163  ierr = m_field.loop_finite_elements("ACOUSTIC_PROBLEM","HELMHOLTZ_RERE_FE",post_proc); CHKERRQ(ierr);
1164 
1165  // if(!pcomm->rank()) {
1166  ierr = post_proc.writeFile("fe_solution_mesh_post_proc.h5m"); CHKERRQ(ierr);
1167  // }
1168 
1169  }
1170 
1171 
1172  }
1173 
1174  // Destroy the KSP solvers
1175  ierr = MatDestroy(&A); CHKERRQ(ierr);
1176  ierr = VecDestroy(&F); CHKERRQ(ierr);
1177  ierr = VecDestroy(&T); CHKERRQ(ierr);
1178  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1179 
1180  ierr = PetscOptionsEnd(); CHKERRQ(ierr);
1181  ierr = PetscFinalize(); CHKERRQ(ierr);
1182 
1183  return 0;
1184 
1185 }
#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.
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.
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
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
CHKERRQ(ierr)
DataForcesAndSourcesCore::EntData EntData
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 char help[]
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

◆ help

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

Definition at line 76 of file fe_approximation.cpp.