v0.13.2
Loading...
Searching...
No Matches
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}
default_options
const std::string default_options
Definition: DefaultParams.hpp:1
param_file
std::string param_file
Definition: DefaultParams.hpp:113
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
COL
@ COL
Definition: definitions.h:123
ROW
@ ROW
Definition: definitions.h:123
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
H1
@ H1
continuous field
Definition: definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
PRESSURESET
@ PRESSURESET
Definition: definitions.h:152
FORCESET
@ FORCESET
Definition: definitions.h:151
NODESET
@ NODESET
Definition: definitions.h:146
SIDESET
@ SIDESET
Definition: definitions.h:147
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::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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::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_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::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:49
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::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2341
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::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1685
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:2169
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
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
D
double D
Definition: initial_diffusion.cpp:44
tol
double tol
Definition: mesh_smoothing.cpp:27
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
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::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
help
static char help[]
Definition: stability.cpp:18
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition: DirichletBC.hpp:211
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
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
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
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:82
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
MoFEM::SnesMethod::snes_x
Vec & snes_x
state vector
Definition: LoopMethods.hpp:121
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition: LoopMethods.hpp:124
MoFEM::SnesMethod::snes_ctx
SNESContext snes_ctx
Definition: LoopMethods.hpp:118
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
NeumannForcesSurfaceComplexForLazy
NonLinear surface pressure element (obsolete implementation)
Definition: SurfacePressureComplexForLazy.hpp:16
NodalForce
Force applied to nodes.
Definition: NodalForce.hpp:13
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:957

Variable Documentation

◆ help

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

Definition at line 18 of file stability.cpp.

Generated by Doxygen 1.9.5 and hosted at University of Glasgow