v0.13.1
Classes | Functions | Variables
stability.cpp File Reference
#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");
295 "EIGEN_VECTOR");
296 CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "D0");
297
298 elastic.feRhs.getOpPtrVector().push_back(
300 "D0", elastic.commonData));
301 elastic.feLhs.getOpPtrVector().push_back(
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");
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 }
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
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;
410 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
411 &Aij);
412
413 // surface forces
414 NeumannForcesSurfaceComplexForLazy neumann_forces(m_field, Aij, F);
416 neumann_forces.getLoopSpatialFe();
418 it)) {
419 CHKERR neumann.addForce(it->getMeshsetId());
420 }
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 }
670
671 SlepcFinalize();
673 // PetscFinalize();
674
675 return 0;
676}
const std::string default_options
std::string param_file
static const double eps
@ COL
Definition: definitions.h:136
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
@ PRESSURESET
Definition: definitions.h:165
@ FORCESET
Definition: definitions.h:164
@ NODESET
Definition: definitions.h:159
@ SIDESET
Definition: definitions.h:160
#define CHKERR
Inline error check.
Definition: definitions.h:548
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
virtual bool check_field(const std::string &name) const =0
check if field is in database
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 problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
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
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
double tol
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
static char help[]
Definition: stability.cpp:31
Set Dirichlet boundary conditions on spatial displacements.
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
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
Managing BitRefLevels.
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
Core (interface) class.
Definition: Core.hpp:92
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Matrix manager is used to build and partition problems.
Interface for managing meshsets containing materials and boundary conditions.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Vec & snes_f
residual
Vec & snes_x
state vector
Mat & snes_B
preconditioner of jacobian matrix
SNESContext snes_ctx
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:33
NonLinear surface pressure element (obsolete implementation)
Force applied to nodes.
Definition: NodalForce.hpp:25
structure grouping operators and data used for calculation of nonlinear elastic element
Post processing.

Variable Documentation

◆ help

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

Definition at line 31 of file stability.cpp.