v0.14.0
Classes | Functions | Variables
stability.cpp File Reference
User modules » NonLinear Elastic Element
#include <BasicFiniteElements.hpp>
#include <Hooke.hpp>
#include <slepceps.h>
#include <SurfacePressureComplexForLazy.hpp>

Go to the source code of this file.

Classes

struct  MyMat_double< TYPE >
 
struct  MyMat< TYPE >
 

Functions

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

Variables

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

Detailed Description

Solves stability problem. Currently uses 3d tetrahedral elements.

Definition in file stability.cpp.

Function Documentation

◆ main()

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

Definition at line 165 of file stability.cpp.

165  {
166 
167  // PetscInitialize(&argc,&argv,(char *)0,help);
168  const string default_options = "-ksp_type fgmres \n"
169  "-pc_type lu \n"
170  "-pc_factor_mat_solver_type mumps \n"
171  "-mat_mumps_icntl_20 0 \n"
172  "-ksp_atol 1e-10 \n"
173  "-ksp_rtol 1e-10 \n"
174  "-snes_monitor \n"
175  "-snes_type newtonls \n"
176  "-snes_linesearch_type basic \n"
177  "-snes_max_it 100 \n"
178  "-snes_atol 1e-7 \n"
179  "-snes_rtol 1e-7 \n"
180  "-ts_monitor \n"
181  "-ts_type alpha \n";
182 
183  string param_file = "param_file.petsc";
184  if (!static_cast<bool>(ifstream(param_file))) {
185  std::ofstream file(param_file.c_str(), std::ios::ate);
186  if (file.is_open()) {
187  file << default_options;
188  file.close();
189  }
190  }
191 
192  SlepcInitialize(&argc, &argv, param_file.c_str(), help);
193  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
194 
195  try {
196 
197  moab::Core mb_instance;
198  moab::Interface &moab = mb_instance;
199  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
200  auto moab_comm_wrap =
201  boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
202  if (pcomm == NULL)
203  pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
204 
205  PetscBool flg = PETSC_TRUE;
206  char mesh_file_name[255];
207  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
208  mesh_file_name, 255, &flg);
209  if (flg != PETSC_TRUE) {
210  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
211  }
212 
213  // use this if your mesh is partotioned and you run code on parts,
214  // you can solve very big problems
215  PetscBool is_partitioned = PETSC_FALSE;
216  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_partitioned",
217  &is_partitioned, &flg);
218 
219  if (is_partitioned == PETSC_TRUE) {
220  // Read mesh to MOAB
221  const char *option;
222  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
223  "PARALLEL_PARTITION;";
224  CHKERR moab.load_file(mesh_file_name, 0, option);
225  CHKERR pcomm->resolve_shared_ents(0, 3, 0);
226  CHKERR pcomm->resolve_shared_ents(0, 3, 1);
227  CHKERR pcomm->resolve_shared_ents(0, 3, 2);
228  } else {
229  const char *option;
230  option = "";
231  CHKERR moab.load_file(mesh_file_name, 0, option);
232  }
233 
234  MoFEM::Core core(moab);
235  MoFEM::Interface &m_field = core;
236 
237  Range CubitSIDESETs_meshsets;
238  CHKERR m_field.getInterface<MeshsetsManager>()->getMeshsetsByType(
239  SIDESET, CubitSIDESETs_meshsets);
240 
241  // ref meshset ref level 0
242  BitRefLevel bit_level0;
243  bit_level0.set(0);
244  EntityHandle meshset_level0;
245  CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
246  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
247  0, 3, bit_level0);
248  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
249  bit_level0, BitRefLevel().set(), meshset_level0);
250 
251  // Fields
252  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
253  3, MB_TAG_SPARSE, MF_ZERO);
254  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
255  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
256  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
257  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
258  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
259 
260  bool check_if_spatial_field_exist = m_field.check_field("SPATIAL_POSITION");
261  CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
262  MB_TAG_SPARSE, MF_ZERO);
263  CHKERR m_field.add_field("EIGEN_VECTOR", H1, AINSWORTH_LEGENDRE_BASE, 3,
264  MB_TAG_SPARSE, MF_ZERO);
265  CHKERR m_field.add_field("D0", H1, AINSWORTH_LEGENDRE_BASE, 3,
266  MB_TAG_SPARSE, MF_ZERO);
267 
268  // add entitities (by tets) to the field
269  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
270  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_VECTOR");
271  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "D0");
272 
273  boost::shared_ptr<Hooke<double>> mat_double =
274  boost::make_shared<Hooke<double>>();
275  boost::shared_ptr<MyMat<adouble>> mat_adouble =
276  boost::make_shared<MyMat<adouble>>();
277 
278  NonlinearElasticElement elastic(m_field, 2);
279  CHKERR elastic.setBlocks(mat_double, mat_adouble);
280  CHKERR elastic.addElement("ELASTIC", "SPATIAL_POSITION");
281  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC",
282  "EIGEN_VECTOR");
283  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "D0");
284 
285  elastic.feRhs.getOpPtrVector().push_back(
286  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
287  "D0", elastic.commonData));
288  elastic.feLhs.getOpPtrVector().push_back(
289  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
290  "D0", elastic.commonData));
291  CHKERR elastic.setOperators("SPATIAL_POSITION");
292 
293  // define problems
294  CHKERR m_field.add_problem("ELASTIC_MECHANICS", MF_ZERO);
295  // set finite elements for problems
296  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
297  "ELASTIC");
298  // set refinement level for problem
299  CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
300  bit_level0);
301 
302  // set app. order
303 
304  PetscInt disp_order;
305  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &disp_order,
306  &flg);
307  if (flg != PETSC_TRUE) {
308  disp_order = 1;
309  }
310 
311  CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", disp_order);
312  CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", disp_order);
313  CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", disp_order);
314  CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
315  CHKERR m_field.set_field_order(0, MBTET, "EIGEN_VECTOR", disp_order);
316  CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_VECTOR", disp_order);
317  CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_VECTOR", disp_order);
318  CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_VECTOR", 1);
319  CHKERR m_field.set_field_order(0, MBTET, "D0", disp_order);
320  CHKERR m_field.set_field_order(0, MBTRI, "D0", disp_order);
321  CHKERR m_field.set_field_order(0, MBEDGE, "D0", disp_order);
322  CHKERR m_field.set_field_order(0, MBVERTEX, "D0", 1);
323 
324  CHKERR m_field.add_finite_element("NEUAMNN_FE", MF_ZERO);
325  CHKERR m_field.modify_finite_element_add_field_row("NEUAMNN_FE",
326  "SPATIAL_POSITION");
327  CHKERR m_field.modify_finite_element_add_field_col("NEUAMNN_FE",
328  "SPATIAL_POSITION");
329  CHKERR m_field.modify_finite_element_add_field_data("NEUAMNN_FE",
330  "SPATIAL_POSITION");
331  CHKERR m_field.modify_finite_element_add_field_data("NEUAMNN_FE",
332  "MESH_NODE_POSITIONS");
333  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
334  "NEUAMNN_FE");
335  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
336  it)) {
337  Range tris;
338  CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, true);
339  CHKERR m_field.add_ents_to_finite_element_by_type(tris, MBTRI,
340  "NEUAMNN_FE");
341  }
342  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
343  m_field, SIDESET | PRESSURESET, it)) {
344  Range tris;
345  CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, true);
346  CHKERR m_field.add_ents_to_finite_element_by_type(tris, MBTRI,
347  "NEUAMNN_FE");
348  }
349  // add nodal force element
350  CHKERR MetaNodalForces::addElement(m_field, "SPATIAL_POSITION");
351  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
352  "FORCE_FE");
353 
354  // build field
355  CHKERR m_field.build_fields();
356  // 10 node tets
357  if (!check_if_spatial_field_exist) {
358  Projection10NodeCoordsOnField ent_method_material(m_field,
359  "MESH_NODE_POSITIONS");
360  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
361  Projection10NodeCoordsOnField ent_method_spatial(m_field,
362  "SPATIAL_POSITION");
363  CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method_spatial);
364  // CHKERR m_field.set_field(0,MBTRI,"SPATIAL_POSITION");
365  // CHKERR m_field.set_field(0,MBTET,"SPATIAL_POSITION");
366  // CHKERR m_field.field_axpy(1,"SPATIAL_POSITION","D0",true);
367  }
368 
369  // build finite elemnts
370  CHKERR m_field.build_finite_elements();
371  // build adjacencies
372  CHKERR m_field.build_adjacencies(bit_level0);
373 
374  // build database
375  ProblemsManager *prb_mng_ptr;
376  CHKERR m_field.getInterface(prb_mng_ptr);
377  if (is_partitioned) {
378  CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("ELASTIC_MECHANICS",
379  true);
380  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS", true, 0,
381  pcomm->size(), 1);
382  } else {
383  CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
384  CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
385  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS");
386  }
387  CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
388 
389  // create matrices
390  Vec F;
391  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
392  "ELASTIC_MECHANICS", ROW, &F);
393  Vec D;
394  CHKERR VecDuplicate(F, &D);
395  Mat Aij;
396  CHKERR m_field.getInterface<MatrixManager>()
397  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
398  &Aij);
399 
400  // surface forces
401  NeumannForcesSurfaceComplexForLazy neumann_forces(m_field, Aij, F);
402  NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE &neumann =
403  neumann_forces.getLoopSpatialFe();
404  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
405  it)) {
406  CHKERR neumann.addForce(it->getMeshsetId());
407  }
408  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
409  m_field, SIDESET | PRESSURESET, it)) {
410  CHKERR neumann.addPressure(it->getMeshsetId());
411  }
412  DirichletSpatialPositionsBc my_Dirichlet_bc(m_field, "SPATIAL_POSITION",
413  Aij, D, F);
414 
415  CHKERR VecZeroEntries(F);
416  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
417  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
418  CHKERR MatZeroEntries(Aij);
419 
420  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
421  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
422  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
423  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
424 
425  // F Vector
426  // preproc
427  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
428  my_Dirichlet_bc.snes_x = D;
429  my_Dirichlet_bc.snes_f = F;
430  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
431  my_Dirichlet_bc);
432  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
433  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
434  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
435  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
436  // elem loops
437  // noadl forces
438  boost::ptr_map<std::string, NodalForce> nodal_forces;
439  string fe_name_str = "FORCE_FE";
440  nodal_forces.insert(fe_name_str, new NodalForce(m_field));
441  CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
442  "SPATIAL_POSITION");
443  boost::ptr_map<std::string, NodalForce>::iterator fit =
444  nodal_forces.begin();
445  for (; fit != nodal_forces.end(); fit++) {
446  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", fit->first,
447  fit->second->getLoopFe());
448  }
449  // surface forces
450  neumann.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
451  neumann.snes_x = D;
452  neumann.snes_f = F;
453  m_field.loop_finite_elements("ELASTIC_MECHANICS", "NEUAMNN_FE", neumann);
454  // stiffnes
455  elastic.getLoopFeRhs().snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
456  elastic.getLoopFeRhs().snes_x = D;
457  elastic.getLoopFeRhs().snes_f = F;
458  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
459  elastic.getLoopFeRhs());
460  // postproc
461  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
462  my_Dirichlet_bc);
463 
464  // Aij Matrix
465  // preproc
466  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
467  my_Dirichlet_bc.snes_B = Aij;
468  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
469  my_Dirichlet_bc);
470  // surface forces
471  // neumann.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
472  // neumann.snes_B = Aij;
473  // CHKERR
474  // m_field.loop_finite_elements("ELASTIC_MECHANICS","NEUAMNN_FE",neumann);
475  // stiffnes
476  elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
477  elastic.getLoopFeLhs().snes_B = Aij;
478  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
479  elastic.getLoopFeLhs());
480  // postproc
481  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
482  my_Dirichlet_bc);
483 
484  CHKERR VecAssemblyBegin(F);
485  CHKERR VecAssemblyEnd(F);
486  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
487  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
488 
489  CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
490  CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
491 
492  // Matrix View
493  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
494  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
495  // std::string wait;
496  // std::cin >> wait;
497 
498  // Solver
499  KSP solver;
500  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
501  CHKERR KSPSetOperators(solver, Aij, Aij);
502  CHKERR KSPSetFromOptions(solver);
503 
504  CHKERR KSPSetUp(solver);
505 
506  CHKERR VecZeroEntries(D);
507  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
508  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
509 
510  CHKERR KSPSolve(solver, F, D);
511  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
512  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
513 
514  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
515  "ELASTIC_MECHANICS", "SPATIAL_POSITION", "D0", COL, D, INSERT_VALUES,
516  SCATTER_REVERSE);
517 
518  Mat Bij;
519  CHKERR MatDuplicate(Aij, MAT_SHARE_NONZERO_PATTERN, &Bij);
520  // CHKERR MatZeroEntries(Aij);
521  CHKERR MatZeroEntries(Bij);
522 
523  /*//Aij Matrix
524  //preproc
525  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
526  my_Dirichlet_bc.snes_B = Aij;
527  CHKERR
528  m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",my_Dirichlet_bc);
529  //stiffnes
530  elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
531  elastic.getLoopFeLhs().snes_B = Aij;
532  CHKERR
533  m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs());
534  //postproc
535  CHKERR
536  m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",my_Dirichlet_bc);
537  */
538 
539  // Bij Matrix
540  mat_adouble->doAotherwiseB = false;
541  // preproc
542  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
543  my_Dirichlet_bc.snes_B = Bij;
544  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
545  my_Dirichlet_bc);
546  // surface forces
547  neumann.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
548  neumann.snes_B = Bij;
549  PetscBool is_conservative = PETSC_FALSE;
550  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_conservative",
551  &is_conservative, &flg);
552  if (is_conservative) {
553  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "NEUAMNN_FE",
554  neumann);
555  }
556  // stiffnes
557  elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
558  elastic.getLoopFeLhs().snes_B = Bij;
559  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
560  elastic.getLoopFeLhs());
561  // postproc
562  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
563  my_Dirichlet_bc);
564 
565  CHKERR MatSetOption(Bij, MAT_SPD, PETSC_TRUE);
566  CHKERR MatAssemblyBegin(Bij, MAT_FINAL_ASSEMBLY);
567  CHKERR MatAssemblyEnd(Bij, MAT_FINAL_ASSEMBLY);
568 
569  // Matrix View
570  // MatView(Bij,PETSC_VIEWER_STDOUT_WORLD);
571  // MatView(Bij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
572  // std::string wait;
573  // std::cin >> wait;
574 
575  EPS eps;
576  ST st;
577  EPSType type;
578  PetscReal tol;
579  PetscInt nev, maxit, its;
580 
581  /*
582  Create eigensolver context
583  */
584  CHKERR EPSCreate(PETSC_COMM_WORLD, &eps);
585  /*
586  Set operators. In this case, it is a generalized eigenvalue problem
587  */
588  CHKERR EPSSetOperators(eps, Bij, Aij);
589  /*
590  Set solver parameters at runtime
591  */
592  CHKERR EPSSetFromOptions(eps);
593  /*
594  Optional: Get some information from the solver and display it
595  */
596  CHKERR EPSSolve(eps);
597 
598  /*
599  Optional: Get some information from the solver and display it
600  */
601  CHKERR EPSGetIterationNumber(eps, &its);
602  PetscPrintf(PETSC_COMM_WORLD, " Number of iterations of the method: %D\n",
603  its);
604  CHKERR EPSGetST(eps, &st);
605  // CHKERR STGetOperationCounters(st,NULL,&lits);
606  // PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method:
607  // %D\n",lits);
608  CHKERR EPSGetType(eps, &type);
609  PetscPrintf(PETSC_COMM_WORLD, " Solution method: %s\n", type);
610  CHKERR EPSGetDimensions(eps, &nev, NULL, NULL);
611  PetscPrintf(PETSC_COMM_WORLD, " Number of requested eigenvalues: %D\n",
612  nev);
613  CHKERR EPSGetTolerances(eps, &tol, &maxit);
614  PetscPrintf(PETSC_COMM_WORLD, " Stopping condition: tol=%.4g, maxit=%D\n",
615  (double)tol, maxit);
616 
617  // get solutions
618  PostProcVolumeOnRefinedMesh post_proc(m_field);
619  CHKERR post_proc.generateReferenceElementMesh();
620  CHKERR post_proc.addFieldValuesGradientPostProc("SPATIAL_POSITION");
621  CHKERR post_proc.addFieldValuesPostProc("SPATIAL_POSITION");
622  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
623  CHKERR post_proc.addFieldValuesPostProc("EIGEN_VECTOR");
624  CHKERR post_proc.addFieldValuesGradientPostProc("EIGEN_VECTOR");
625  CHKERR post_proc.addFieldValuesPostProc("D0");
626  CHKERR post_proc.addFieldValuesGradientPostProc("D0");
627  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
628  post_proc);
629  CHKERR post_proc.writeFile("out.h5m");
630 
631  PetscScalar eigr, eigi, nrm2r;
632  for (int nn = 0; nn < nev; nn++) {
633  CHKERR EPSGetEigenpair(eps, nn, &eigr, &eigi, D, PETSC_NULL);
634  CHKERR VecNorm(D, NORM_2, &nrm2r);
635  PetscPrintf(
636  PETSC_COMM_WORLD,
637  " ncov = %D eigr = %.4g eigi = %.4g (inv eigr = %.4g) nrm2r = %.4g\n",
638  nn, eigr, eigi, 1. / eigr, nrm2r);
639  std::ostringstream o1;
640  o1 << "eig_" << nn << ".h5m";
641  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
642  "ELASTIC_MECHANICS", "SPATIAL_POSITION", "EIGEN_VECTOR", COL, D,
643  INSERT_VALUES, SCATTER_REVERSE);
644  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
645  post_proc);
646  CHKERR post_proc.writeFile(o1.str().c_str());
647  }
648 
649  CHKERR KSPDestroy(&solver);
650  CHKERR VecDestroy(&F);
651  CHKERR VecDestroy(&D);
652  CHKERR MatDestroy(&Aij);
653  CHKERR MatDestroy(&Bij);
654  CHKERR EPSDestroy(&eps);
655  }
656  CATCH_ERRORS;
657 
658  SlepcFinalize();
659  MoFEM::Core::Finalize();
660  // PetscFinalize();
661 
662  return 0;
663 }

Variable Documentation

◆ help

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

Definition at line 18 of file stability.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
SIDESET
@ SIDESET
Definition: definitions.h:147
MoFEM::CoreInterface::problem_basic_method_postProcess
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:87
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NeumannForcesSurfaceComplexForLazy
NonLinear surface pressure element (obsolete implementation)
Definition: SurfacePressureComplexForLazy.hpp:16
MetaNodalForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
PRESSURESET
@ PRESSURESET
Definition: definitions.h:152
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::SnesMethod::snes_x
Vec & snes_x
state vector
Definition: LoopMethods.hpp:121
MoFEM::SnesMethod::snes_ctx
SNESContext snes_ctx
Definition: LoopMethods.hpp:118
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition: LoopMethods.hpp:124
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addForce
MoFEMErrorCode addForce(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:516
NodalForce
Force applied to nodes.
Definition: NodalForce.hpp:13
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MetaNodalForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
NODESET
@ NODESET
Definition: definitions.h:146
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::CoreInterface::problem_basic_method_preProcess
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
FORCESET
@ FORCESET
Definition: definitions.h:151
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::add_field
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.
convert.type
type
Definition: convert.py:64
COL
@ COL
Definition: definitions.h:123
help
static char help[]
Definition: stability.cpp:18
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE
Definition: SurfacePressureComplexForLazy.hpp:39
MoFEM::ProblemsManager::partitionFiniteElements
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
Definition: ProblemsManager.cpp:2167
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
Range
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::ProblemsManager::buildProblemOnDistributedMesh
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
Definition: ProblemsManager.cpp:290
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addPressure
MoFEMErrorCode addPressure(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:531
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CoreInterface::modify_problem_add_finite_element
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
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::set_field_order
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.
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2339
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1683
tol
double tol
Definition: mesh_smoothing.cpp:27
F
@ F
Definition: free_surface.cpp:394
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition: DirichletBC.hpp:211
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow