178 {
179
180
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
197 if (!
static_cast<bool>(ifstream(
param_file))) {
198 std::ofstream file(
param_file.c_str(), std::ios::ate);
199 if (file.is_open()) {
201 file.close();
202 }
203 }
204
207
208 try {
209
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;
222 if (flg != PETSC_TRUE) {
223 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
224 }
225
226
227
228 PetscBool is_partitioned = PETSC_FALSE;
230 &is_partitioned, &flg);
231
232 if (is_partitioned == PETSC_TRUE) {
233
234 const char *option;
235 option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
236 "PARALLEL_PARTITION;";
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 = "";
245 }
246
249
250 Range CubitSIDESETs_meshsets;
252 SIDESET, CubitSIDESETs_meshsets);
253
254
256 bit_level0.set(0);
257 EntityHandle meshset_level0;
258 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
260 0, 3, bit_level0);
263
264
272
273 bool check_if_spatial_field_exist = m_field.
check_field(
"SPATIAL_POSITION");
280
281
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
292 CHKERR elastic.setBlocks(mat_double, mat_adouble);
293 CHKERR elastic.addElement(
"ELASTIC",
"SPATIAL_POSITION");
295 "EIGEN_VECTOR");
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
308
310 "ELASTIC");
311
313 bit_level0);
314
315
316
317 PetscInt disp_order;
319 &flg);
320 if (flg != PETSC_TRUE) {
321 disp_order = 1;
322 }
323
336
339 "SPATIAL_POSITION");
341 "SPATIAL_POSITION");
343 "SPATIAL_POSITION");
345 "MESH_NODE_POSITIONS");
347 "NEUAMNN_FE");
349 it)) {
350 Range tris;
351 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
353 "NEUAMNN_FE");
354 }
357 Range tris;
358 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
360 "NEUAMNN_FE");
361 }
362
365 "FORCE_FE");
366
367
369
370 if (!check_if_spatial_field_exist) {
372 "MESH_NODE_POSITIONS");
375 "SPATIAL_POSITION");
377
378
379
380 }
381
382
384
386
387
390 if (is_partitioned) {
392 true);
394 pcomm->size(), 1);
395 } else {
399 }
401
402
405 "ELASTIC_MECHANICS",
ROW, &F);
408 Mat Aij;
410 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
411 &Aij);
412
413
416 neumann_forces.getLoopSpatialFe();
418 it)) {
420 }
424 }
427
429 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
430 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
431 CHKERR MatZeroEntries(Aij);
432
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
439
440 my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
441 my_Dirichlet_bc.snes_x =
D;
442 my_Dirichlet_bc.snes_f = F;
444 my_Dirichlet_bc);
445 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
446 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
448 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
449
450
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));
455 "SPATIAL_POSITION");
456 boost::ptr_map<std::string, NodalForce>::iterator fit =
457 nodal_forces.begin();
458 for (; fit != nodal_forces.end(); fit++) {
460 fit->second->getLoopFe());
461 }
462
463 neumann.
snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
467
468 elastic.getLoopFeRhs().snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
469 elastic.getLoopFeRhs().snes_x =
D;
470 elastic.getLoopFeRhs().snes_f = F;
472 elastic.getLoopFeRhs());
473
475 my_Dirichlet_bc);
476
477
478
479 my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
480 my_Dirichlet_bc.snes_B = Aij;
482 my_Dirichlet_bc);
483
484
485
486
487
488
489 elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
490 elastic.getLoopFeLhs().snes_B = Aij;
492 elastic.getLoopFeLhs());
493
495 my_Dirichlet_bc);
496
497 CHKERR VecAssemblyBegin(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
506
507
508
509
510
511
512 KSP solver;
513 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
514 CHKERR KSPSetOperators(solver, Aij, Aij);
515 CHKERR KSPSetFromOptions(solver);
516
518
520 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
521 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
522
524 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
525 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
526
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
534 CHKERR MatZeroEntries(Bij);
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553 mat_adouble->doAotherwiseB = false;
554
555 my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
556 my_Dirichlet_bc.snes_B = Bij;
558 my_Dirichlet_bc);
559
560 neumann.
snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
562 PetscBool is_conservative = PETSC_FALSE;
564 &is_conservative, &flg);
565 if (is_conservative) {
567 neumann);
568 }
569
570 elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
571 elastic.getLoopFeLhs().snes_B = Bij;
573 elastic.getLoopFeLhs());
574
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
583
584
585
586
587
589 ST st;
592 PetscInt nev, maxit, its;
593
594
595
596
597 CHKERR EPSCreate(PETSC_COMM_WORLD, &
eps);
598
599
600
602
603
604
606
607
608
610
611
612
613
615 PetscPrintf(PETSC_COMM_WORLD, " Number of iterations of the method: %D\n",
616 its);
618
619
620
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);
627 PetscPrintf(PETSC_COMM_WORLD, " Stopping condition: tol=%.4g, maxit=%D\n",
629
630
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");
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";
655 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"EIGEN_VECTOR",
COL,
D,
656 INSERT_VALUES, SCATTER_REVERSE);
658 post_proc);
659 CHKERR post_proc.writeFile(o1.str().c_str());
660 }
661
662 CHKERR KSPDestroy(&solver);
668 }
670
671 SlepcFinalize();
673
674
675 return 0;
676}
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 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
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)
DeprecatedCoreInterface Interface
const double D
diffusivity
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