Implementation of bone remodeling finite element. More...

#include "users_modules/bone_remodelling/src/Remodeling.hpp"

Collaboration diagram for BoneRemodeling::Remodeling:

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.

MoFEMErrorCode	addFields ()
	Set and add entities to approximation fields.

MoFEMErrorCode	addElements ()
	Set and add finite elements.

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

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

MoFEMErrorCode	addMomentumFluxes ()
	Finite elements to calculate tractions.

MoFEMErrorCode	buildDM ()
	Set problem and DM.

MoFEMErrorCode	solveDM ()
	Solve problem set up in DM.

Public Attributes
MoFEM::Interface &	mField

CommonData &	commonData

Detailed Description

Implementation of bone remodeling finite element.

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

Examples: mofem/users_modules/bone_remodelling/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) {}

BoneRemodeling::Remodeling::mField

MoFEM::Interface & mField

Definition Remodeling.hpp:206

BoneRemodeling::Remodeling::commonData

CommonData & commonData

Definition Remodeling.hpp:207

Member Function Documentation

◆ addElements()

MoFEMErrorCode BoneRemodeling::Remodeling::addElements ( )

Set and add finite elements.

Returns: Error code

Examples: mofem/users_modules/bone_remodelling/bone_adaptation.cpp, mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 1690 of file Remodeling.cpp.

                                       {
  MoFEMFunctionBegin;
 
  CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
  CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING",
                                                    "DISPLACEMENTS");
  CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING",
                                                    "DISPLACEMENTS");
  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
                                                     "DISPLACEMENTS");
  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
  CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
                                                     "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
  CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
  CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_all, MBTET,
                                                   "FE_REMODELLING");
 
  CHKERR mField.add_finite_element("ELASTIC");
  CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
  CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
  CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
                                                     "DISPLACEMENTS");
  CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
                                                     "MESH_NODE_POSITIONS");
 
  if (commonData.nOremodellingBlock) {
    CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_block,
                                                     MBTET, "ELASTIC");
  }
  
  commonData.elasticMaterialsPtr =
      boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
  commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
      new NonlinearElasticElement(mField, 10));
  CHKERR commonData.elasticMaterialsPtr->setBlocks(
      commonData.elasticPtr->setOfBlocks);
  CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
  CHKERR commonData.elasticPtr->setOperators(
      "DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
 
  CHKERR mField.build_finite_elements();
  CHKERR mField.build_adjacencies(commonData.bitLevel);
 
  // Allocate memory for density and gradient of displacements at integration
  // points
  commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
  commonData.data.gradRhoPtr =
      boost::shared_ptr<MatrixDouble>(new MatrixDouble());
  commonData.data.gradDispPtr =
      boost::shared_ptr<MatrixDouble>(new MatrixDouble());
 
  // Add operators Rhs
  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feRhs), true, false,
                     false, false);
  auto &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
  // Get density at integration points
  list_of_operators_rhs.push_back(
      new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
  list_of_operators_rhs.push_back(
      new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
  list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
  // Get displacement gradient at integration points
  list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
      "DISPLACEMENTS", commonData.data.gradDispPtr));
  list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
  list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
  list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
 
  // Add operators Lhs
  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feLhs), true, false,
                     false, false);
  auto &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
  // Get density at integration points
  list_of_operators_lhs.push_back(
      new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
  list_of_operators_rhs.push_back(
      new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
  // Get displacement gradient at integration points
  list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
      "DISPLACEMENTS", commonData.data.gradDispPtr));
  if (commonData.with_adol_c) {
    list_of_operators_lhs.push_back(
        new OpCalculateStressTangentWithAdolc(commonData));
    list_of_operators_lhs.push_back(
        new OpAssmbleStressLhs_dF<true>(commonData));
  } else {
    list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
    list_of_operators_lhs.push_back(
        new OpAssmbleStressLhs_dF<false>(commonData));
  }
  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
  list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
  list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
 
  MoFEMFunctionReturn(0);
}

◆ addElementsTestingDensity()

MoFEMErrorCode BoneRemodeling::Remodeling::addElementsTestingDensity ( )

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

Returns: Error code

Examples: mofem/users_modules/bone_remodelling/src/Remodeling.hpp.

◆ addElementsTestingElasticty()

MoFEMErrorCode BoneRemodeling::Remodeling::addElementsTestingElasticty ( )

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

Returns: Error code

Examples: mofem/users_modules/bone_remodelling/src/Remodeling.hpp.

◆ addFields()

MoFEMErrorCode BoneRemodeling::Remodeling::addFields ( )

Set and add entities to approximation fields.

Returns: Error code

Examples: mofem/users_modules/bone_remodelling/bone_adaptation.cpp, mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 1603 of file Remodeling.cpp.

                                     {
 
  MoFEMFunctionBegin;
 
  // Seed all mesh entities to MoFEM database, those entities can be potentially
  // used as finite elements or as entities which carry some approximation
  // field.
  commonData.bitLevel.set(0);
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
      0, 3, commonData.bitLevel);
  int order = commonData.oRder;
 
  // Add displacement field
  CHKERR mField.add_field("DISPLACEMENTS", H1,
                          /*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
                          MB_TAG_SPARSE, MF_ZERO);
  // Add field representing ho-geometry
  CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
                          MB_TAG_SPARSE, MF_ZERO);
 
  // Check if density is available, if not add density field.
  bool add_rho_field = false;
  if (!mField.check_field("RHO")) {
    CHKERR mField.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1,
                            MB_TAG_SPARSE, MF_ZERO);
    // FIXME
    CHKERR mField.add_ents_to_field_by_type(commonData.tEts_all, MBTET, "RHO");
    // CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
    CHKERR mField.synchronise_field_entities("RHO");
    add_rho_field = true;
 
    CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
    CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
    CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
    CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
  }
 
  // Add entities to field
  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
 
  // Set approximation order to entities
  CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
  CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
  CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
  CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
 
  // Assumes that geometry is approximated using 2nd order polynomials.
 
  // Apply 2nd order only on skin
  {
    // Skinner skin(&mField.get_moab());
    // Range faces,tets;
    // CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
    // CHKERR skin.find_skin(0,tets,false,faces);
    // Range edges;
    // CHKERR mField.get_moab().get_adjacencies(
    // faces,1,false,edges,moab::Interface::UNION
    // );
    CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
  }
 
  CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
 
  // Build fields
  CHKERR mField.build_fields();
 
  // If order was not set from CT scans set homogenous order equal to
  // reference bone density
  if (add_rho_field) {
    CHKERR mField.getInterface<FieldBlas>()->setField(commonData.rHo_ref,
                                                      MBVERTEX, "RHO");
    // const DofEntity_multiIndex *dofs_ptr;
    // CHKERR mField.get_dofs(&dofs_ptr);
    // for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
    //   cerr << (*dit)->getFieldData() << endl;
    // }
  }
 
  // Project geometry given on 10-node tets on ho-geometry
  Projection10NodeCoordsOnField ent_method_material(mField,
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
 
  MoFEMFunctionReturn(0);
}

◆ addMomentumFluxes()

MoFEMErrorCode BoneRemodeling::Remodeling::addMomentumFluxes ( )

Finite elements to calculate tractions.

Returns: Error code

Examples: mofem/users_modules/bone_remodelling/bone_adaptation.cpp, mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 1788 of file Remodeling.cpp.

                                             {
 
  MoFEMFunctionBegin;
  // Add Neumann forces elements
  CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
  CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
  CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
 
  // forces and pressures on surface
  CHKERR MetaNeummanForces::setMomentumFluxOperators(
      mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
  // noadl forces
  CHKERR MetaNodalForces::setOperators(mField, commonData.nodalForces,
                                       PETSC_NULL, "DISPLACEMENTS");
  // edge forces
  CHKERR MetaEdgeForces::setOperators(mField, commonData.edgeForces, PETSC_NULL,
                                      "DISPLACEMENTS");
 
  for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
           commonData.neumannForces.begin();
       mit != commonData.neumannForces.end(); mit++) {
    mit->second->methodsOp.push_back(
        new TimeForceScale("-my_load_history", false));
  }
  for (boost::ptr_map<string, NodalForce>::iterator mit =
           commonData.nodalForces.begin();
       mit != commonData.nodalForces.end(); mit++) {
    mit->second->methodsOp.push_back(
        new TimeForceScale("-my_load_history", false));
  }
  for (boost::ptr_map<string, EdgeForce>::iterator mit =
           commonData.edgeForces.begin();
       mit != commonData.edgeForces.end(); mit++) {
    mit->second->methodsOp.push_back(
        new TimeForceScale("-my_load_history", false));
  }
 
  MoFEMFunctionReturn(0);
}

◆ buildDM()

MoFEMErrorCode BoneRemodeling::Remodeling::buildDM ( )

Set problem and DM.

Returns: Error code

Examples: mofem/users_modules/bone_remodelling/bone_adaptation.cpp, mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 1828 of file Remodeling.cpp.

                                   {
 
  MoFEMFunctionBegin;
  commonData.dm_name = "DM_REMODELING";
  // Register DM problem
  CHKERR DMRegister_MoFEM(commonData.dm_name);
  CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
  CHKERR DMSetType(commonData.dm, commonData.dm_name);
  CHKERR DMMoFEMSetIsPartitioned(commonData.dm, PETSC_TRUE);
  // Create DM instance
  CHKERR DMMoFEMCreateMoFEM(commonData.dm, &mField, commonData.dm_name,
                            commonData.bitLevel);
  // Configure DM form line command options (DM itself, solvers,
  // pre-conditioners, ... )
  CHKERR DMSetFromOptions(commonData.dm);
  // Add elements to dm (only one here)
  CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
  CHKERR DMMoFEMAddElement(commonData.dm, "ELASTIC");
 
  if (mField.check_finite_element("FORCE_FE")) {
    CHKERR DMMoFEMAddElement(commonData.dm, "FORCE_FE");
  }
  if (mField.check_finite_element("PRESSURE_FE")) {
    CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
  }
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
  CHKERR DMSetUp(commonData.dm);
 
  MoFEMFunctionReturn(0);
}

◆ 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: mofem/users_modules/bone_remodelling/bone_adaptation.cpp, mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 1576 of file Remodeling.cpp.

                                         {
 
  MoFEMFunctionBegin;
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
 
  commonData.oRder = 2;
  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
                         &commonData.oRder, PETSC_NULL);
 
  PetscOptionsEnd();
 
  CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    string name = it->getName();
    if (name.compare(0, 14, "NO_REMODELLING") == 0) {
      commonData.nOremodellingBlock = true;
      EntityHandle meshset = it->getMeshset();
      CHKERR this->mField.get_moab().get_entities_by_type(
          meshset, MBTET, commonData.tEts_block, true);
      commonData.tEts_all =
          subtract(commonData.tEts_all, commonData.tEts_block);
    }
  }
 
  MoFEMFunctionReturn(0);
}

◆ solveDM()

MoFEMErrorCode BoneRemodeling::Remodeling::solveDM ( )

Solve problem set up in DM.

Returns: Error code

Examples: mofem/users_modules/bone_remodelling/bone_adaptation.cpp, mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 1859 of file Remodeling.cpp.

                                   {
 
  MoFEMFunctionBegin;
 
  Vec F = commonData.F;
  Vec D = commonData.D;
  Mat A = commonData.A;
  DM dm = commonData.dm;
  TS ts = commonData.ts;
 
  CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
  CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
  double ftime = 1;
#if PETSC_VERSION_GE(3, 8, 0)
  CHKERR TSSetMaxTime(ts, ftime);
#endif
  CHKERR TSSetFromOptions(ts);
  CHKERR TSSetDM(ts, dm);
 
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
 
  // Set up FIELDSPLIT
  // Only is user set -pc_type fieldsplit
  if (is_pcfs == PETSC_TRUE) {
    IS is_disp, is_rho;
    const MoFEM::Problem *problem_ptr;
    CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
    CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
        problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
    CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
        problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
    CHKERR ISSort(is_disp);
    CHKERR ISSort(is_rho);
    CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
    CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
    CHKERR ISDestroy(&is_disp);
    CHKERR ISDestroy(&is_rho);
  }
 
  // Monitor
  MonitorPostProc post_proc(mField, commonData);
  TsCtx *ts_ctx;
  DMMoFEMGetTsCtx(dm, &ts_ctx);
  {
    ts_ctx->getPostProcessMonitor().push_back(&post_proc);
    CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
  }
 
#if PETSC_VERSION_GE(3, 7, 0)
  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
  CHKERR TSSolve(ts, D);
  CHKERR TSGetTime(ts, &ftime);
  PetscInt steps, snesfails, rejects, nonlinits, linits;
#if PETSC_VERSION_LE(3, 8, 0)
  CHKERR TSGetTimeStepNumber(ts, &steps);
#else
  CHKERR TSGetStepNumber(ts, &steps);
#endif
 
  CHKERR TSGetSNESFailures(ts, &snesfails);
  CHKERR TSGetStepRejections(ts, &rejects);
  CHKERR TSGetSNESIterations(ts, &nonlinits);
  CHKERR TSGetKSPIterations(ts, &linits);
  PetscPrintf(PETSC_COMM_WORLD,
              "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
              "linits %D\n",
              steps, rejects, snesfails, ftime, nonlinits, linits);
 
  if (commonData.is_atom_testing) {
    if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
  }
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

◆ commonData

CommonData& BoneRemodeling::Remodeling::commonData

Examples: mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 207 of file Remodeling.hpp.

◆ mField

MoFEM::Interface& BoneRemodeling::Remodeling::mField

Examples: mofem/users_modules/bone_remodelling/src/Remodeling.hpp, and mofem/users_modules/bone_remodelling/src/impl/Remodeling.cpp.

Definition at line 206 of file Remodeling.hpp.

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

users_modules/bone_remodelling/src/Remodeling.hpp
users_modules/bone_remodelling/src/impl/Remodeling.cpp

Classes

Public Member Functions

Public Attributes