v0.13.1
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 178 of file stability.cpp.

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

Variable Documentation

◆ help

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

Definition at line 31 of file stability.cpp.

Generated on Tue Jul 5 2022 22:52:24 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow