1403 {
1405
1408 PetscBool analysis_mesh_flg = PETSC_FALSE;
1410
1411 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling post-process",
1412 "none");
1413
1415 "-mass_postproc",
1416 "is used for testing, calculates mass and energy at each time step", "",
1418
1419 CHKERR PetscOptionsBool(
"-analysis_mesh",
1420 "saves analysis mesh at each time step", "",
1421 analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1422
1423 CHKERR PetscOptionsScalar(
1424 "-equilibrium_stop_rate",
1425 "is used to stop calculations when equilibium state is achieved", "",
1427 PetscOptionsEnd();
1428
1430
1431 PetscOptionsBegin(PETSC_COMM_WORLD, "",
1432 "Bone remodeling post-process", "none");
1433
1435 CHKERR PetscOptionsInt(
"-my_output_prt",
1436 "frequncy how often results are dumped on hard disk",
1438 PetscOptionsEnd();
1439
1445 }
1446
1447
1448
1449 {
1450 auto &list_of_operators =
postProc.getOpPtrVector();
1451
1452 list_of_operators.push_back(
1456
1462 }
1463
1468 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1474 }
1475
1477 }
1478
1482 int mass_vec_ghost[] = {0};
1484 1, 1, mass_vec_ghost, &mass_vec);
1485
1486 CHKERR VecZeroEntries(mass_vec);
1487 CHKERR VecDuplicate(mass_vec, &energ_vec);
1488
1489 Remodeling::Fe energyCalc(
mField);
1491 true);
1492 energyCalc.getOpPtrVector().push_back(
1496 energyCalc.getOpPtrVector().push_back(
1500 energyCalc.getOpPtrVector().push_back(
1501 new OpMassAndEnergyCalculation(
"RHO",
commonData, energ_vec, mass_vec));
1503 &energyCalc);
1504
1505 CHKERR VecAssemblyBegin(mass_vec);
1506 CHKERR VecAssemblyEnd(mass_vec);
1507 double mass_sum;
1508 CHKERR VecSum(mass_vec, &mass_sum);
1509
1510 CHKERR VecAssemblyBegin(energ_vec);
1511 CHKERR VecAssemblyEnd(energ_vec);
1512 double energ_sum;
1513 CHKERR VecSum(energ_vec, &energ_sum);
1514
1515 CHKERR PetscPrintf(PETSC_COMM_WORLD,
1516 "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1517 mass_sum, energ_sum);
1518 CHKERR VecDestroy(&mass_vec);
1519 CHKERR VecDestroy(&energ_vec);
1520
1522 double equilibrium_rate =
1524 double equilibrium_mass_rate =
1526 if (equilibrium_rate <
rate) {
1528 PETSC_COMM_WORLD,
1529 "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1530 equilibrium_rate * 100);
1531 }
1532 if (equilibrium_mass_rate <
rate) {
1534 PETSC_COMM_WORLD,
1535 "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1536 equilibrium_mass_rate * 100);
1537 }
1540 }
1541 }
1542
1543 int step;
1544#if PETSC_VERSION_LE(3, 8, 0)
1545 CHKERR TSGetTimeStepNumber(ts, &step);
1546#else
1547 CHKERR TSGetStepNumber(ts, &step);
1548#endif
1549
1550 if ((step) %
pRT == 0) {
1551
1552 ostringstream sss;
1553 sss << "out_" << step << ".h5m";
1556 "PARALLEL=WRITE_PART");
1557 if (analysis_mesh_flg) {
1558 ostringstream ttt;
1559 ttt << "analysis_mesh_" << step << ".h5m";
1561 "PARALLEL=WRITE_PART");
1562 }
1563
1567 ostringstream ss;
1568 ss << "out_elastic_" << step << ".h5m";
1570 "PARALLEL=WRITE_PART");
1571 }
1572 }
1574}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field gradient.
const FTensor::Tensor2< T, Dim, Dim > Vec
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
PetscBool equilibrium_flg
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
PetscBool less_post_proc
reduce file size
DM dm
Discretization manager.
double cUrrent_psi
current free energy for evaluating equilibrium state
boost::shared_ptr< NonlinearElasticElement > elasticPtr
double cUrrent_mass
current free energy for evaluating equilibrium state
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision scalar field values calculation.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
