1404 {
1406
1409 PetscBool analysis_mesh_flg = PETSC_FALSE;
1411
1412 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling post-process",
1413 "none");
1414
1416 "-mass_postproc",
1417 "is used for testing, calculates mass and energy at each time step", "",
1419
1420 CHKERR PetscOptionsBool(
"-analysis_mesh",
1421 "saves analysis mesh at each time step", "",
1422 analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
1423
1424 CHKERR PetscOptionsScalar(
1425 "-equilibrium_stop_rate",
1426 "is used to stop calculations when equilibium state is achieved", "",
1428 ierr = PetscOptionsEnd();
1430
1432
1433 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
1434 "Bone remodeling post-process", "none");
1435
1437 CHKERR PetscOptionsInt(
"-my_output_prt",
1438 "frequncy how often results are dumped on hard disk",
1440 ierr = PetscOptionsEnd();
1442
1448 }
1449
1450
1451
1452 {
1453 auto &list_of_operators =
postProc.getOpPtrVector();
1454
1455 list_of_operators.push_back(
1459
1465 }
1466
1471 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
1477 }
1478
1480 }
1481
1485 int mass_vec_ghost[] = {0};
1487 1, 1, mass_vec_ghost, &mass_vec);
1488
1489 CHKERR VecZeroEntries(mass_vec);
1490 CHKERR VecDuplicate(mass_vec, &energ_vec);
1491
1492 Remodeling::Fe energyCalc(
mField);
1494 true);
1495 energyCalc.getOpPtrVector().push_back(
1499 energyCalc.getOpPtrVector().push_back(
1503 energyCalc.getOpPtrVector().push_back(
1504 new OpMassAndEnergyCalculation(
"RHO",
commonData, energ_vec, mass_vec));
1506 &energyCalc);
1507
1508 CHKERR VecAssemblyBegin(mass_vec);
1509 CHKERR VecAssemblyEnd(mass_vec);
1510 double mass_sum;
1511 CHKERR VecSum(mass_vec, &mass_sum);
1512
1513 CHKERR VecAssemblyBegin(energ_vec);
1514 CHKERR VecAssemblyEnd(energ_vec);
1515 double energ_sum;
1516 CHKERR VecSum(energ_vec, &energ_sum);
1517
1518 CHKERR PetscPrintf(PETSC_COMM_WORLD,
1519 "Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
1520 mass_sum, energ_sum);
1521 CHKERR VecDestroy(&mass_vec);
1522 CHKERR VecDestroy(&energ_vec);
1523
1525 double equilibrium_rate =
1527 double equilibrium_mass_rate =
1529 if (equilibrium_rate <
rate) {
1531 PETSC_COMM_WORLD,
1532 "Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
1533 equilibrium_rate * 100);
1534 }
1535 if (equilibrium_mass_rate <
rate) {
1537 PETSC_COMM_WORLD,
1538 "Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
1539 equilibrium_mass_rate * 100);
1540 }
1543 }
1544 }
1545
1546 int step;
1547#if PETSC_VERSION_LE(3, 8, 0)
1548 CHKERR TSGetTimeStepNumber(ts, &step);
1549#else
1550 CHKERR TSGetStepNumber(ts, &step);
1551#endif
1552
1553 if ((step) %
pRT == 0) {
1554
1555 ostringstream sss;
1556 sss << "out_" << step << ".h5m";
1559 "PARALLEL=WRITE_PART");
1560 if (analysis_mesh_flg) {
1561 ostringstream ttt;
1562 ttt << "analysis_mesh_" << step << ".h5m";
1564 "PARALLEL=WRITE_PART");
1565 }
1566
1570 ostringstream ss;
1571 ss << "out_elastic_" << step << ".h5m";
1573 "PARALLEL=WRITE_PART");
1574 }
1575 }
1577}
#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 addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
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 filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
