165 {
166
167
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
184 if (!
static_cast<bool>(ifstream(
param_file))) {
185 std::ofstream file(
param_file.c_str(), std::ios::ate);
186 if (file.is_open()) {
188 file.close();
189 }
190 }
191
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;
209 if (flg != PETSC_TRUE) {
210 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
211 }
212
213
214
215 PetscBool is_partitioned = PETSC_FALSE;
217 &is_partitioned, &flg);
218
219 if (is_partitioned == PETSC_TRUE) {
220
221 const char *option;
222 option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
223 "PARALLEL_PARTITION;";
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 = "";
232 }
233
236
237 Range CubitSIDESETs_meshsets;
239 SIDESET, CubitSIDESETs_meshsets);
240
241
243 bit_level0.set(0);
245 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
247 0, 3, bit_level0);
250
251
259
260 bool check_if_spatial_field_exist = m_field.
check_field(
"SPATIAL_POSITION");
267
268
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
279 CHKERR elastic.setBlocks(mat_double, mat_adouble);
280 CHKERR elastic.addElement(
"ELASTIC",
"SPATIAL_POSITION");
282 "EIGEN_VECTOR");
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
295
297 "ELASTIC");
298
300 bit_level0);
301
302
303
304 PetscInt disp_order;
306 &flg);
307 if (flg != PETSC_TRUE) {
308 disp_order = 1;
309 }
310
323
326 "SPATIAL_POSITION");
328 "SPATIAL_POSITION");
330 "SPATIAL_POSITION");
332 "MESH_NODE_POSITIONS");
334 "NEUAMNN_FE");
336 it)) {
338 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
340 "NEUAMNN_FE");
341 }
345 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
347 "NEUAMNN_FE");
348 }
349
352 "FORCE_FE");
353
354
356
357 if (!check_if_spatial_field_exist) {
359 "MESH_NODE_POSITIONS");
362 "SPATIAL_POSITION");
364
365
366
367 }
368
369
371
373
374
377 if (is_partitioned) {
379 true);
381 pcomm->size(), 1);
382 } else {
386 }
388
389
392 "ELASTIC_MECHANICS",
ROW, &F);
395 Mat Aij;
397 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
398 &Aij);
399
400
403 neumann_forces.getLoopSpatialFe();
405 it)) {
407 }
411 }
414
416 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
417 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
418 CHKERR MatZeroEntries(Aij);
419
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
426
428 my_Dirichlet_bc.snes_x =
D;
429 my_Dirichlet_bc.snes_f = F;
431 my_Dirichlet_bc);
432 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
433 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
435 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
436
437
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));
442 "SPATIAL_POSITION");
443 boost::ptr_map<std::string, NodalForce>::iterator fit =
444 nodal_forces.begin();
445 for (; fit != nodal_forces.end(); fit++) {
447 fit->second->getLoopFe());
448 }
449
454
456 elastic.getLoopFeRhs().snes_x =
D;
457 elastic.getLoopFeRhs().snes_f = F;
459 elastic.getLoopFeRhs());
460
462 my_Dirichlet_bc);
463
464
465
467 my_Dirichlet_bc.snes_B = Aij;
469 my_Dirichlet_bc);
470
471
472
473
474
475
477 elastic.getLoopFeLhs().snes_B = Aij;
479 elastic.getLoopFeLhs());
480
482 my_Dirichlet_bc);
483
484 CHKERR VecAssemblyBegin(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
493
494
495
496
497
498
499 KSP solver;
500 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
501 CHKERR KSPSetOperators(solver, Aij, Aij);
502 CHKERR KSPSetFromOptions(solver);
503
505
507 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
508 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
509
511 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
512 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
513
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
521 CHKERR MatZeroEntries(Bij);
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540 mat_adouble->doAotherwiseB = false;
541
543 my_Dirichlet_bc.snes_B = Bij;
545 my_Dirichlet_bc);
546
549 PetscBool is_conservative = PETSC_FALSE;
551 &is_conservative, &flg);
552 if (is_conservative) {
554 neumann);
555 }
556
558 elastic.getLoopFeLhs().snes_B = Bij;
560 elastic.getLoopFeLhs());
561
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
570
571
572
573
574
576 ST st;
579 PetscInt nev, maxit, its;
580
581
582
583
584 CHKERR EPSCreate(PETSC_COMM_WORLD, &
eps);
585
586
587
589
590
591
593
594
595
597
598
599
600
602 PetscPrintf(PETSC_COMM_WORLD, " Number of iterations of the method: %D\n",
603 its);
605
606
607
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);
614 PetscPrintf(PETSC_COMM_WORLD, " Stopping condition: tol=%.4g, maxit=%D\n",
616
617
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");
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";
642 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"EIGEN_VECTOR",
COL,
D,
643 INSERT_VALUES, SCATTER_REVERSE);
645 post_proc);
646 CHKERR post_proc.writeFile(o1.str().c_str());
647 }
648
649 CHKERR KSPDestroy(&solver);
655 }
657
658 SlepcFinalize();
660
661
662 return 0;
663}
const std::string default_options
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERR
Inline error check.
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
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.
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)
Set Dirichlet boundary conditions on spatial displacements.
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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.
Mat & snes_B
preconditioner of jacobian matrix
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
MoFEMErrorCode addPressure(int ms_id)
MoFEMErrorCode addForce(int ms_id)
NonLinear surface pressure element (obsolete implementation)
structure grouping operators and data used for calculation of nonlinear elastic element