v0.13.2
Loading...
Searching...
No Matches
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 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");
282 "EIGEN_VECTOR");
283 CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "D0");
284
285 elastic.feRhs.getOpPtrVector().push_back(
287 "D0", elastic.commonData));
288 elastic.feLhs.getOpPtrVector().push_back(
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");
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 }
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
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;
397 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
398 &Aij);
399
400 // surface forces
401 NeumannForcesSurfaceComplexForLazy neumann_forces(m_field, Aij, F);
403 neumann_forces.getLoopSpatialFe();
405 it)) {
406 CHKERR neumann.addForce(it->getMeshsetId());
407 }
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
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
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 }
657
658 SlepcFinalize();
660 // PetscFinalize();
661
662 return 0;
663}
const std::string default_options
std::string param_file
static const double eps
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ PRESSURESET
Definition: definitions.h:152
@ FORCESET
Definition: definitions.h:151
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
#define CHKERR
Inline error check.
Definition: definitions.h:535
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 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 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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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 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
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
double D
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:40
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)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
static char help[]
Definition: stability.cpp:18
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:128
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
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:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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:23
NonLinear surface pressure element (obsolete implementation)
Force applied to nodes.
Definition: NodalForce.hpp:13
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 18 of file stability.cpp.