v0.14.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Public Attributes | List of all members
BoneRemodeling::Remodeling Struct Reference

Implementation of bone remodeling finite element. More...

#include <users_modules/bone_remodelling/src/Remodeling.hpp>

Collaboration diagram for BoneRemodeling::Remodeling:
[legend]

Classes

struct  CommonData
 
struct  Fe
 Volume finite element. More...
 
struct  FePrePostProcessLhs
 Not used at this stage. Could be used to do some calculations, before assembly of local elements. More...
 
struct  FePrePostProcessRhs
 Not used at this stage. Could be used to do some calculations, before assembly of local elements. More...
 

Public Member Functions

 Remodeling (MoFEM::Interface &m_field, CommonData &common_data)
 
MoFEMErrorCode getParameters ()
 Get parameters form line command or config file. More...
 
MoFEMErrorCode addFields ()
 Set and add entities to approximation fields. More...
 
MoFEMErrorCode addElements ()
 Set and add finite elements. More...
 
MoFEMErrorCode addElementsTestingDensity ()
 (Testing only) Set finite element to run mass transport problem only More...
 
MoFEMErrorCode addElementsTestingElasticty ()
 (Testing only) Set finite element to run elastic problem only More...
 
MoFEMErrorCode addMomentumFluxes ()
 Finite elements to calculate tractions. More...
 
MoFEMErrorCode buildDM ()
 Set problem and DM. More...
 
MoFEMErrorCode solveDM ()
 Solve problem set up in DM. More...
 

Public Attributes

MoFEM::InterfacemField
 
CommonDatacommonData
 

Detailed Description

Implementation of bone remodeling finite element.

Implementation base on paper [77] http://biomechanics.stanford.edu/paper/IJSS12.pdf

Examples
bone_adaptation.cpp.

Definition at line 36 of file Remodeling.hpp.

Constructor & Destructor Documentation

◆ Remodeling()

BoneRemodeling::Remodeling::Remodeling ( MoFEM::Interface m_field,
CommonData common_data 
)
inline

Definition at line 209 of file Remodeling.hpp.

210 : mField(m_field), commonData(common_data) {}
MoFEM::Interface & mField
Definition: Remodeling.hpp:206

Member Function Documentation

◆ addElements()

MoFEMErrorCode BoneRemodeling::Remodeling::addElements ( )

Set and add finite elements.

Returns
Error code
Examples
Remodeling.cpp, Remodeling.hpp, and bone_adaptation.cpp.

Definition at line 1694 of file Remodeling.cpp.

1694 {
1696
1697 CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
1699 "DISPLACEMENTS");
1701 "DISPLACEMENTS");
1703 "DISPLACEMENTS");
1704 CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
1706 "MESH_NODE_POSITIONS");
1707 CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
1708 CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
1710 "FE_REMODELLING");
1711
1712 CHKERR mField.add_finite_element("ELASTIC");
1713 CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
1714 CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
1716 "DISPLACEMENTS");
1718 "MESH_NODE_POSITIONS");
1719
1722 MBTET, "ELASTIC");
1723 }
1724
1726 boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
1727 commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
1730 commonData.elasticPtr->setOfBlocks);
1731 CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
1732 CHKERR commonData.elasticPtr->setOperators(
1733 "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
1734
1737
1738 // Allocate memory for density and gradient of displacements at integration
1739 // points
1740 commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
1742 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1744 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1745
1746 // Add operators Rhs
1747 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feRhs), true, false,
1748 false, false);
1749 auto &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
1750 // Get density at integration points
1751 list_of_operators_rhs.push_back(
1753 list_of_operators_rhs.push_back(
1755 list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
1756 // Get displacement gradient at integration points
1757 list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1758 "DISPLACEMENTS", commonData.data.gradDispPtr));
1759 list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
1760 list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
1761 list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
1762
1763 // Add operators Lhs
1764 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feLhs), true, false,
1765 false, false);
1766 auto &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
1767 // Get density at integration points
1768 list_of_operators_lhs.push_back(
1770 list_of_operators_rhs.push_back(
1772 // Get displacement gradient at integration points
1773 list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1774 "DISPLACEMENTS", commonData.data.gradDispPtr));
1775 if (commonData.with_adol_c) {
1776 list_of_operators_lhs.push_back(
1777 new OpCalculateStressTangentWithAdolc(commonData));
1778 list_of_operators_lhs.push_back(
1779 new OpAssmbleStressLhs_dF<true>(commonData));
1780 } else {
1781 list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
1782 list_of_operators_lhs.push_back(
1783 new OpAssmbleStressLhs_dF<false>(commonData));
1784 }
1785 list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
1786 list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
1787 list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
1788
1790}
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
Definition: Remodeling.hpp:185
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
Definition: Remodeling.hpp:187
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
Definition: Remodeling.hpp:188
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
Definition: Remodeling.hpp:141
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
structure grouping operators and data used for calculation of nonlinear elastic element

◆ addElementsTestingDensity()

MoFEMErrorCode BoneRemodeling::Remodeling::addElementsTestingDensity ( )

(Testing only) Set finite element to run mass transport problem only

Returns
Error code
Examples
Remodeling.hpp.

◆ addElementsTestingElasticty()

MoFEMErrorCode BoneRemodeling::Remodeling::addElementsTestingElasticty ( )

(Testing only) Set finite element to run elastic problem only

Returns
Error code
Examples
Remodeling.hpp.

◆ addFields()

MoFEMErrorCode BoneRemodeling::Remodeling::addFields ( )

Set and add entities to approximation fields.

Returns
Error code
Examples
Remodeling.cpp, Remodeling.hpp, and bone_adaptation.cpp.

Definition at line 1607 of file Remodeling.cpp.

1607 {
1608
1610
1611 // Seed all mesh entities to MoFEM database, those entities can be potentially
1612 // used as finite elements or as entities which carry some approximation
1613 // field.
1614 commonData.bitLevel.set(0);
1615 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
1616 0, 3, commonData.bitLevel);
1617 int order = commonData.oRder;
1618
1619 // Add displacement field
1620 CHKERR mField.add_field("DISPLACEMENTS", H1,
1621 /*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
1622 MB_TAG_SPARSE, MF_ZERO);
1623 // Add field representing ho-geometry
1624 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
1625 MB_TAG_SPARSE, MF_ZERO);
1626
1627 // Check if density is available, if not add density field.
1628 bool add_rho_field = false;
1629 if (!mField.check_field("RHO")) {
1631 MB_TAG_SPARSE, MF_ZERO);
1632 // FIXME
1634 // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
1636 add_rho_field = true;
1637
1638 CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
1639 CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
1640 CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
1641 CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
1642 }
1643
1644 // Add entities to field
1645 CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
1646 CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
1647
1648 // Set approximation order to entities
1649 CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
1650 CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
1651 CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
1652 CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
1653
1654 // Assumes that geometry is approximated using 2nd order polynomials.
1655
1656 // Apply 2nd order only on skin
1657 {
1658 // Skinner skin(&mField.get_moab());
1659 // Range faces,tets;
1660 // CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
1661 // CHKERR skin.find_skin(0,tets,false,faces);
1662 // Range edges;
1663 // CHKERR mField.get_moab().get_adjacencies(
1664 // faces,1,false,edges,moab::Interface::UNION
1665 // );
1666 CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
1667 }
1668
1669 CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
1670
1671 // Build fields
1673
1674 // If order was not set from CT scans set homogenous order equal to
1675 // reference bone density
1676 if (add_rho_field) {
1678 MBVERTEX, "RHO");
1679 // const DofEntity_multiIndex *dofs_ptr;
1680 // CHKERR mField.get_dofs(&dofs_ptr);
1681 // for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
1682 // cerr << (*dit)->getFieldData() << endl;
1683 // }
1684 }
1685
1686 // Project geometry given on 10-node tets on ho-geometry
1687 Projection10NodeCoordsOnField ent_method_material(mField,
1688 "MESH_NODE_POSITIONS");
1689 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
1690
1692}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
constexpr int order
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
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.
Managing BitRefLevels.
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.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ addMomentumFluxes()

MoFEMErrorCode BoneRemodeling::Remodeling::addMomentumFluxes ( )

Finite elements to calculate tractions.

Returns
Error code
Examples
Remodeling.cpp, Remodeling.hpp, and bone_adaptation.cpp.

Definition at line 1792 of file Remodeling.cpp.

1792 {
1793
1795 // Add Neumann forces elements
1796 CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
1797 CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
1798 CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
1799
1800 // forces and pressures on surface
1801 CHKERR MetaNeummanForces::setMomentumFluxOperators(
1802 mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
1803 // noadl forces
1805 PETSC_NULL, "DISPLACEMENTS");
1806 // edge forces
1808 "DISPLACEMENTS");
1809
1810 for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1811 commonData.neumannForces.begin();
1812 mit != commonData.neumannForces.end(); mit++) {
1813 mit->second->methodsOp.push_back(
1814 new TimeForceScale("-my_load_history", false));
1815 }
1816 for (boost::ptr_map<string, NodalForce>::iterator mit =
1817 commonData.nodalForces.begin();
1818 mit != commonData.nodalForces.end(); mit++) {
1819 mit->second->methodsOp.push_back(
1820 new TimeForceScale("-my_load_history", false));
1821 }
1822 for (boost::ptr_map<string, EdgeForce>::iterator mit =
1823 commonData.edgeForces.begin();
1824 mit != commonData.edgeForces.end(); mit++) {
1825 mit->second->methodsOp.push_back(
1826 new TimeForceScale("-my_load_history", false));
1827 }
1828
1830}
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
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
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
Force scale operator for reading two columns.

◆ buildDM()

MoFEMErrorCode BoneRemodeling::Remodeling::buildDM ( )

Set problem and DM.

Returns
Error code
Examples
Remodeling.cpp, Remodeling.hpp, and bone_adaptation.cpp.

Definition at line 1832 of file Remodeling.cpp.

1832 {
1833
1835 commonData.dm_name = "DM_REMODELING";
1836 // Register DM problem
1838 CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
1841 // Create DM instance
1844 // Configure DM form line command options (DM itself, solvers,
1845 // pre-conditioners, ... )
1846 CHKERR DMSetFromOptions(commonData.dm);
1847 // Add elements to dm (only one here)
1848 CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
1850
1851 if (mField.check_finite_element("FORCE_FE")) {
1853 }
1854 if (mField.check_finite_element("PRESSURE_FE")) {
1855 CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
1856 }
1857 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1858 CHKERR DMSetUp(commonData.dm);
1859
1861}
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
Problem manager is used to build and partition problems.

◆ getParameters()

MoFEMErrorCode BoneRemodeling::Remodeling::getParameters ( )

Get parameters form line command or config file.

Read command line and config file to setup material and model parameters.

Returns
Error code
Examples
Remodeling.cpp, Remodeling.hpp, and bone_adaptation.cpp.

Definition at line 1579 of file Remodeling.cpp.

1579 {
1580
1582 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
1583
1584 commonData.oRder = 2;
1585 CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
1586 &commonData.oRder, PETSC_NULL);
1587
1588 ierr = PetscOptionsEnd();
1589 CHKERRQ(ierr);
1590
1591 CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
1593 string name = it->getName();
1594 if (name.compare(0, 14, "NO_REMODELLING") == 0) {
1596 EntityHandle meshset = it->getMeshset();
1597 CHKERR this->mField.get_moab().get_entities_by_type(
1598 meshset, MBTET, commonData.tEts_block, true);
1601 }
1602 }
1603
1605}
@ BLOCKSET
Definition: definitions.h:148
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
virtual moab::Interface & get_moab()=0

◆ solveDM()

MoFEMErrorCode BoneRemodeling::Remodeling::solveDM ( )

Solve problem set up in DM.

Returns
Error code
Examples
Remodeling.cpp, Remodeling.hpp, and bone_adaptation.cpp.

Definition at line 1863 of file Remodeling.cpp.

1863 {
1864
1866
1867 Vec F = commonData.F;
1868 Vec D = commonData.D;
1869 Mat A = commonData.A;
1870 DM dm = commonData.dm;
1871 TS ts = commonData.ts;
1872
1873 CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
1874 CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
1875 double ftime = 1;
1876#if PETSC_VERSION_GE(3, 8, 0)
1877 CHKERR TSSetMaxTime(ts, ftime);
1878#endif
1879 CHKERR TSSetFromOptions(ts);
1880 CHKERR TSSetDM(ts, dm);
1881
1882 SNES snes;
1883 CHKERR TSGetSNES(ts, &snes);
1884 KSP ksp;
1885 CHKERR SNESGetKSP(snes, &ksp);
1886 PC pc;
1887 CHKERR KSPGetPC(ksp, &pc);
1888 PetscBool is_pcfs = PETSC_FALSE;
1889 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1890
1891 // Set up FIELDSPLIT
1892 // Only is user set -pc_type fieldsplit
1893 if (is_pcfs == PETSC_TRUE) {
1894 IS is_disp, is_rho;
1895 const MoFEM::Problem *problem_ptr;
1896 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
1897 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1898 problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
1899 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1900 problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
1901 CHKERR ISSort(is_disp);
1902 CHKERR ISSort(is_rho);
1903 CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1904 CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1905 CHKERR ISDestroy(&is_disp);
1906 CHKERR ISDestroy(&is_rho);
1907 }
1908
1909 // Monitor
1910 MonitorPostProc post_proc(mField, commonData);
1911 TsCtx *ts_ctx;
1912 DMMoFEMGetTsCtx(dm, &ts_ctx);
1913 {
1914 ts_ctx->getPostProcessMonitor().push_back(&post_proc);
1915 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1916 }
1917
1918#if PETSC_VERSION_GE(3, 7, 0)
1919 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1920#endif
1921 CHKERR TSSolve(ts, D);
1922 CHKERR TSGetTime(ts, &ftime);
1923 PetscInt steps, snesfails, rejects, nonlinits, linits;
1924#if PETSC_VERSION_LE(3, 8, 0)
1925 CHKERR TSGetTimeStepNumber(ts, &steps);
1926#else
1927 CHKERR TSGetStepNumber(ts, &steps);
1928#endif
1929
1930 CHKERR TSGetSNESFailures(ts, &snesfails);
1931 CHKERR TSGetStepRejections(ts, &rejects);
1932 CHKERR TSGetSNESIterations(ts, &nonlinits);
1933 CHKERR TSGetKSPIterations(ts, &linits);
1934 PetscPrintf(PETSC_COMM_WORLD,
1935 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1936 "linits %D\n",
1937 steps, rejects, snesfails, ftime, nonlinits, linits);
1938
1940 if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
1941 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
1942 }
1943
1945}
@ ROW
Definition: definitions.h:123
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ F
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1128
double D
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
constexpr AssemblyType A
PetscBool is_atom_testing
for atom tests
Definition: Remodeling.hpp:162
double cUrrent_psi
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:158
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
keeps basic data about problem
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:144

Member Data Documentation

◆ commonData

CommonData& BoneRemodeling::Remodeling::commonData
Examples
Remodeling.cpp, and Remodeling.hpp.

Definition at line 207 of file Remodeling.hpp.

◆ mField

MoFEM::Interface& BoneRemodeling::Remodeling::mField
Examples
Remodeling.cpp, and Remodeling.hpp.

Definition at line 206 of file Remodeling.hpp.


The documentation for this struct was generated from the following files: