v0.9.0
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 [47] 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 
)
Examples
Remodeling.hpp.

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
bone_adaptation.cpp, and Remodeling.hpp.

Definition at line 1715 of file Remodeling.cpp.

1715  {
1717 
1718  CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
1720  "DISPLACEMENTS");
1722  "DISPLACEMENTS");
1724  "DISPLACEMENTS");
1725  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
1727  "MESH_NODE_POSITIONS");
1728  CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
1729  CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
1731  "FE_REMODELLING");
1732 
1733  CHKERR mField.add_finite_element("ELASTIC");
1734  CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
1735  CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
1737  "DISPLACEMENTS");
1739  "MESH_NODE_POSITIONS");
1740 
1743  MBTET, "ELASTIC");
1744  }
1745 
1747  boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
1748  commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
1749  new NonlinearElasticElement(mField, 10));
1751  commonData.elasticPtr->setOfBlocks);
1752  CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
1753  CHKERR commonData.elasticPtr->setOperators(
1754  "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
1755 
1758 
1759  // Allocate memory for density and gradient of displacements at integration
1760  // points
1761  commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
1763  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1765  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1766 
1767  // Add operators Rhs
1768  boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
1769  &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
1770 
1771  // Get density at integration points
1772  list_of_operators_rhs.push_back(
1774  list_of_operators_rhs.push_back(
1776  list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
1777  // Get displacement gradient at integration points
1778  list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1779  "DISPLACEMENTS", commonData.data.gradDispPtr));
1780  list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
1781  list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
1782  list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
1783 
1784  // Add operators Lhs
1785  boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
1786  &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
1787  // Get density at integration points
1788  list_of_operators_lhs.push_back(
1790  list_of_operators_rhs.push_back(
1792  // Get displacement gradient at integration points
1793  list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1794  "DISPLACEMENTS", commonData.data.gradDispPtr));
1795  if (commonData.with_adol_c) {
1796  list_of_operators_lhs.push_back(
1797  new OpCalculateStressTangentWithAdolc(commonData));
1798  list_of_operators_lhs.push_back(
1799  new OpAssmbleStressLhs_dF<true>(commonData));
1800  } else {
1801  list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
1802  list_of_operators_lhs.push_back(
1803  new OpAssmbleStressLhs_dF<false>(commonData));
1804  }
1805  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
1806  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
1807  list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
1808 
1810 }
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
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
Definition: Remodeling.hpp:187
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 build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
Definition: Remodeling.hpp:188
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
Definition: Remodeling.hpp:141
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
Definition: Remodeling.hpp:185
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
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
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
Get value at integration points for scalar field.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131

◆ 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
bone_adaptation.cpp, and Remodeling.hpp.

Definition at line 1628 of file Remodeling.cpp.

1628  {
1629 
1631 
1632  // Seed all mesh entities to MoFEM database, those entities can be potentially
1633  // used as finite elements or as entities which carry some approximation
1634  // field.
1635  commonData.bitLevel.set(0);
1636  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
1637  0, 3, commonData.bitLevel);
1638  int order = commonData.oRder;
1639 
1640  // Add displacement field
1641  CHKERR mField.add_field("DISPLACEMENTS", H1,
1642  /*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
1643  MB_TAG_SPARSE, MF_ZERO);
1644  // Add field representing ho-geometry
1645  CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
1646  MB_TAG_SPARSE, MF_ZERO);
1647 
1648  // Check if density is available, if not add density field.
1649  bool add_rho_field = false;
1650  if (!mField.check_field("RHO")) {
1652  MB_TAG_SPARSE, MF_ZERO);
1653  // FIXME
1655  // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
1657  add_rho_field = true;
1658 
1659  CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
1660  CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
1661  CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
1662  CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
1663  }
1664 
1665  // Add entities to field
1666  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
1667  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
1668 
1669  // Set approximation order to entities
1670  CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
1671  CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
1672  CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
1673  CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
1674 
1675  // Assumes that geometry is approximated using 2nd order polynomials.
1676 
1677  // Apply 2nd order only on skin
1678  {
1679  // Skinner skin(&mField.get_moab());
1680  // Range faces,tets;
1681  // CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
1682  // CHKERR skin.find_skin(0,tets,false,faces);
1683  // Range edges;
1684  // CHKERR mField.get_moab().get_adjacencies(
1685  // faces,1,false,edges,moab::Interface::UNION
1686  // );
1687  CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
1688  }
1689 
1690  CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
1691 
1692  // Build fields
1694 
1695  // If order was not set from CT scans set homogenous order equal to
1696  // reference bone density
1697  if (add_rho_field) {
1699  MBVERTEX, "RHO");
1700  // const DofEntity_multiIndex *dofs_ptr;
1701  // CHKERR mField.get_dofs(&dofs_ptr);
1702  // for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
1703  // cerr << (*dit)->getFieldData() << endl;
1704  // }
1705  }
1706 
1707  // Project geometry given on 10-node tets on ho-geometry
1708  Projection10NodeCoordsOnField ent_method_material(mField,
1709  "MESH_NODE_POSITIONS");
1710  CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
1711 
1713 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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.
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.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
Basic algebra on fields.
Definition: FieldBlas.hpp:34
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.
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
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.
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
continuous field
Definition: definitions.h:171

◆ addMomentumFluxes()

MoFEMErrorCode BoneRemodeling::Remodeling::addMomentumFluxes ( )

Finite elements to calculate tractions.

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

Definition at line 1812 of file Remodeling.cpp.

1812  {
1813 
1815  // Add Neumann forces elements
1816  CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
1817  CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
1818  CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
1819 
1820  // forces and pressures on surface
1821  CHKERR MetaNeummanForces::setMomentumFluxOperators(
1822  mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
1823  // noadl forces
1825  PETSC_NULL, "DISPLACEMENTS");
1826  // edge forces
1828  "DISPLACEMENTS");
1829 
1830  for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
1831  commonData.neumannForces.begin();
1832  mit != commonData.neumannForces.end(); mit++) {
1833  mit->second->methodsOp.push_back(
1834  new TimeForceScale("-my_load_history", false));
1835  }
1836  for (boost::ptr_map<string, NodalForce>::iterator mit =
1837  commonData.nodalForces.begin();
1838  mit != commonData.nodalForces.end(); mit++) {
1839  mit->second->methodsOp.push_back(
1840  new TimeForceScale("-my_load_history", false));
1841  }
1842  for (boost::ptr_map<string, EdgeForce>::iterator mit =
1843  commonData.edgeForces.begin();
1844  mit != commonData.edgeForces.end(); mit++) {
1845  mit->second->methodsOp.push_back(
1846  new TimeForceScale("-my_load_history", false));
1847  }
1848 
1850 }
Force scale operator for reading two columns.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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:74
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:140
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
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:104
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
#define CHKERR
Inline error check.
Definition: definitions.h:596
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: EdgeForce.hpp:109
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138

◆ buildDM()

MoFEMErrorCode BoneRemodeling::Remodeling::buildDM ( )

Set problem and DM.

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

Definition at line 1852 of file Remodeling.cpp.

1852  {
1853 
1855  commonData.dm_name = "DM_REMODELING";
1856  // Register DM problem
1858  CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
1859  CHKERR DMSetType(commonData.dm, commonData.dm_name);
1861  // Create DM instance
1864  // Configure DM form line command options (DM itself, solvers,
1865  // pre-conditioners, ... )
1866  CHKERR DMSetFromOptions(commonData.dm);
1867  // Add elements to dm (only one here)
1868  CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
1869  CHKERR DMMoFEMAddElement(commonData.dm, "ELASTIC");
1870 
1871  if (mField.check_finite_element("FORCE_FE")) {
1872  CHKERR DMMoFEMAddElement(commonData.dm, "FORCE_FE");
1873  }
1874  if (mField.check_finite_element("PRESSURE_FE")) {
1875  CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
1876  }
1877  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1878  CHKERR DMSetUp(commonData.dm);
1879 
1881 }
Problem manager is used to build and partition problems.
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:414
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:53
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: DMMMoFEM.cpp:109
#define CHKERR
Inline error check.
Definition: definitions.h:596
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:971

◆ 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
bone_adaptation.cpp, and Remodeling.hpp.

Definition at line 1600 of file Remodeling.cpp.

1600  {
1601 
1603  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
1604 
1605  commonData.oRder = 2;
1606  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
1607  &commonData.oRder, PETSC_NULL);
1608 
1609  ierr = PetscOptionsEnd();
1610  CHKERRQ(ierr);
1611 
1612  CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
1614  string name = it->getName();
1615  if (name.compare(0, 14, "NO_REMODELLING") == 0) {
1617  EntityHandle meshset = it->getMeshset();
1618  CHKERR this->mField.get_moab().get_entities_by_type(
1619  meshset, MBTET, commonData.tEts_block, true);
1622  }
1623  }
1624 
1626 }
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
CHKERRQ(ierr)
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ solveDM()

MoFEMErrorCode BoneRemodeling::Remodeling::solveDM ( )

Solve problem set up in DM.

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

Definition at line 1883 of file Remodeling.cpp.

1883  {
1884 
1886 
1887  Vec F = commonData.F;
1888  Vec D = commonData.D;
1889  Mat A = commonData.A;
1890  DM dm = commonData.dm;
1891  TS ts = commonData.ts;
1892 
1893  CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
1894  CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
1895  double ftime = 1;
1896 #if PETSC_VERSION_GE(3, 8, 0)
1897  CHKERR TSSetMaxTime(ts, ftime);
1898 #endif
1899  CHKERR TSSetFromOptions(ts);
1900  CHKERR TSSetDM(ts, dm);
1901 
1902  SNES snes;
1903  CHKERR TSGetSNES(ts, &snes);
1904  KSP ksp;
1905  CHKERR SNESGetKSP(snes, &ksp);
1906  PC pc;
1907  CHKERR KSPGetPC(ksp, &pc);
1908  PetscBool is_pcfs = PETSC_FALSE;
1909  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1910 
1911  // Set up FIELDSPLIT
1912  // Only is user set -pc_type fieldsplit
1913  if (is_pcfs == PETSC_TRUE) {
1914  IS is_disp, is_rho;
1915  const MoFEM::Problem *problem_ptr;
1916  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
1917  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1918  problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
1919  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1920  problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
1921  CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
1922  CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
1923  CHKERR ISDestroy(&is_disp);
1924  CHKERR ISDestroy(&is_rho);
1925  }
1926 
1927  // Monitor
1928  MonitorPostProc post_proc(mField, commonData);
1929  TsCtx *ts_ctx;
1930  DMMoFEMGetTsCtx(dm, &ts_ctx);
1931  {
1932  ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
1933  CHKERR TSMonitorSet(ts, f_TSMonitorSet, ts_ctx, PETSC_NULL);
1934  }
1935 
1936 #if PETSC_VERSION_GE(3, 7, 0)
1937  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1938 #endif
1939  CHKERR TSSolve(ts, D);
1940  CHKERR TSGetTime(ts, &ftime);
1941  PetscInt steps, snesfails, rejects, nonlinits, linits;
1942 #if PETSC_VERSION_LE(3, 8, 0)
1943  CHKERR TSGetTimeStepNumber(ts, &steps);
1944 #else
1945  CHKERR TSGetStepNumber(ts, &steps);
1946 #endif
1947 
1948  CHKERR TSGetSNESFailures(ts, &snesfails);
1949  CHKERR TSGetStepRejections(ts, &rejects);
1950  CHKERR TSGetSNESIterations(ts, &nonlinits);
1951  CHKERR TSGetKSPIterations(ts, &linits);
1952  PetscPrintf(PETSC_COMM_WORLD,
1953  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
1954  "linits %D\n",
1955  steps, rejects, snesfails, ftime, nonlinits, linits);
1956 
1958  if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
1959  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
1960  }
1961 
1963 }
double cUrrent_psi
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:158
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
BasicMethodsSequence & get_postProcess_to_do_Monitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:187
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:27
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
keeps basic data about problemThis is low level structure with information about problem,...
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
PetscBool is_atom_testing
for atom tests
Definition: Remodeling.hpp:162
#define CHKERR
Inline error check.
Definition: definitions.h:596
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:990
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:337
DEPRECATED PetscErrorCode f_TSMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Definition: TsCtx.hpp:327
std::string getName() const

Member Data Documentation

◆ commonData

CommonData& BoneRemodeling::Remodeling::commonData
Examples
Remodeling.hpp.

Definition at line 207 of file Remodeling.hpp.

◆ mField

MoFEM::Interface& BoneRemodeling::Remodeling::mField
Examples
Remodeling.hpp.

Definition at line 206 of file Remodeling.hpp.


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