v0.14.0
bone_adaptation.cpp
Go to the documentation of this file.
1 /**
2  * \file bone_adaptation.cpp
3  * \example bone_adaptation.cpp
4  * \ingroup bone_remodelling
5  *
6  * \brief Main program for running bone adaptation
7  *
8  */
9 
10 /* This file is part of MoFEM.
11 
12 * MoFEM is free software: you can redistribute it and/or modify it under
13 * the terms of the GNU Lesser General Public License as published by the
14 * Free Software Foundation, either version 3 of the License, or (at your
15 * option) any later version.
16 *
17 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
18 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20 * License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
24 
25 #include <boost/program_options.hpp>
26 using namespace std;
27 namespace po = boost::program_options;
28 
29 #include <BasicFiniteElements.hpp>
30 #include <ElasticMaterials.hpp>
31 #include <Remodeling.hpp>
32 
33 using namespace BoneRemodeling;
34 
35 static char help[] = "\n";
36 
37 // #ifndef WITH_METAIO
38 // #error "You need compile users modules with MetaIO -DWITH_METAIO=1"
39 // #endif
40 
41 int main(int argc, char *argv[]) {
42 
43  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
44 
45  try {
46 
47  char mesh_file_name[255];
48  PetscBool flg_file = PETSC_FALSE;
49  PetscBool flg_elastic_only = PETSC_FALSE;
50 
51  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Main set up", "none");
52  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
53  mesh_file_name, 255, &flg_file);
54 
55  CHKERR PetscOptionsBool(
56  "-elastic_only", "is used for testing, calculates only elastic problem",
57  "", PETSC_FALSE, &flg_elastic_only, PETSC_NULL);
58 
59  ierr = PetscOptionsEnd();
60  CHKERRQ(ierr);
61 
62  moab::Core mb_instance;
63  moab::Interface &moab = mb_instance;
64 
65  // Read parameters from line command
66  if (flg_file != PETSC_TRUE) {
67  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
68  }
69  const char *option;
70 
71  // option = "";
72  option = "PARALLEL=READ_PART;"
73  "PARALLEL_RESOLVE_SHARED_ENTS;"
74  "PARTITION=PARALLEL_PARTITION;";
75  CHKERR moab.load_file(mesh_file_name, 0, option);
76  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
77  if (pcomm == NULL)
78  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
79 
80  MoFEM::Core core(moab);
81  MoFEM::Interface &m_field = core;
82 
83  // Add meshsets with material and boundary conditions
84  MeshsetsManager *meshsets_manager_ptr;
85  CHKERR m_field.getInterface<MeshsetsManager>(meshsets_manager_ptr);
86  CHKERR meshsets_manager_ptr->setMeshsetFromFile();
87 
88  PetscPrintf(PETSC_COMM_WORLD, "Read meshset. Added meshsets for bc.cfg\n");
89  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
90  PetscPrintf(PETSC_COMM_WORLD, "%s",
91  static_cast<std::ostringstream &>(
92  std::ostringstream().seekp(0) << *it << endl)
93  .str()
94  .c_str());
95  cerr << *it << endl;
96  }
97 
98  Remodeling::CommonData common_data;
99  common_data.feRhs =
100  boost::shared_ptr<Remodeling::Fe>(new Remodeling::Fe(m_field));
101  common_data.feLhs =
102  boost::shared_ptr<Remodeling::Fe>(new Remodeling::Fe(m_field));
103  common_data.preProcRhs = boost::shared_ptr<Remodeling::FePrePostProcessRhs>(
105  common_data.preProcLhs = boost::shared_ptr<Remodeling::FePrePostProcessLhs>(
107 
108  CHKERR common_data.getParameters();
109 
110  Remodeling remodelling(m_field, common_data);
111 
112  // Remodeling
113  if (!flg_elastic_only) {
114  remodelling.getParameters();
115  remodelling.addFields();
116  remodelling.addMomentumFluxes();
117  remodelling.addElements();
118  remodelling.buildDM();
119 
120  bool flag_cubit_disp = false;
122  m_field, NODESET | DISPLACEMENTSET, it)) {
123  flag_cubit_disp = true;
124  }
125 
126  boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
127  // if normally defined boundary conditions are not found, try to use
128  // DISPLACEMENT blockset
129  if (!flag_cubit_disp) {
130  dirihlet_bc_ptr =
131  boost::shared_ptr<FEMethod>(new DirichletSetFieldFromBlockWithFlags(
132  m_field, "DISPLACEMENTS", "DISPLACEMENT"));
133  dynamic_cast<DirichletSetFieldFromBlockWithFlags &>(*dirihlet_bc_ptr)
134  .methodsOp.push_back(new TimeForceScale("-my_load_history", false));
135  } else {
136  dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
137  new DirichletDisplacementBc(m_field, "DISPLACEMENTS"));
138  dynamic_cast<DirichletDisplacementBc &>(*dirihlet_bc_ptr)
139  .methodsOp.push_back(new TimeForceScale("-my_load_history", false));
140  }
141 
142  // Add elements Time Solver via DM interface.
143  // This part is to calculate residual vector.
144  CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL,
145  common_data.preProcRhs.get(), NULL);
146  CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL,
147  dirihlet_bc_ptr.get(), NULL);
148  CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING",
149  common_data.feRhs.get(), NULL, NULL);
150  if (common_data.nOremodellingBlock) {
151  CHKERR DMMoFEMTSSetIFunction(common_data.dm, "ELASTIC",
152  &(common_data.elasticPtr->getLoopFeRhs()),
153  NULL, NULL);
154  }
155 
156  // Add surface forces
157  for (boost::ptr_map<string, NeummanForcesSurface>::iterator fit =
158  common_data.neumannForces.begin();
159  fit != common_data.neumannForces.end(); fit++) {
160  CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
161  &fit->second->getLoopFe(), NULL, NULL);
162  }
163 
164  // Add edge forces
165  for (boost::ptr_map<std::string, EdgeForce>::iterator fit =
166  common_data.edgeForces.begin();
167  fit != common_data.edgeForces.end(); fit++) {
168  CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
169  &fit->second->getLoopFe(), NULL, NULL);
170  }
171 
172  // Add nodal forces
173  for (boost::ptr_map<std::string, NodalForce>::iterator fit =
174  common_data.nodalForces.begin();
175  fit != common_data.nodalForces.end(); fit++) {
176  CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
177  &fit->second->getLoopFe(), NULL, NULL);
178  }
179 
180  CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL, NULL,
181  dirihlet_bc_ptr.get());
182  CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL, NULL,
183  common_data.preProcRhs.get());
184 
185  // Add elements Time Solver via DM interface.
186  // This part is to calculate tangent matrix.
187  CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL,
188  common_data.preProcLhs.get(), NULL);
189  CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL,
190  dirihlet_bc_ptr.get(), NULL);
191  CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING",
192  common_data.feLhs.get(), NULL, NULL);
193  if (common_data.nOremodellingBlock) {
194  CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "ELASTIC",
195  &common_data.elasticPtr->getLoopFeLhs(),
196  NULL, NULL);
197  }
198  CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL, NULL,
199  dirihlet_bc_ptr.get());
200  CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL, NULL,
201  common_data.preProcLhs.get());
202 
203  CHKERR DMCreateMatrix(common_data.dm, &common_data.A);
204  CHKERR DMCreateGlobalVector(common_data.dm, &common_data.D);
205  CHKERR DMCreateGlobalVector(common_data.dm, &common_data.F);
206 
207  // Initialize vectors with initial data from mesh
208  CHKERR DMoFEMMeshToLocalVector(common_data.dm, common_data.D,
209  INSERT_VALUES, SCATTER_FORWARD);
210 
211  CHKERR TSCreate(PETSC_COMM_WORLD, &common_data.ts);
212  CHKERR TSSetType(common_data.ts, TSBEULER);
213 
214  CHKERR remodelling.solveDM();
215 
216  CHKERR MatDestroy(&common_data.A);
217  CHKERR VecDestroy(&common_data.D);
218  CHKERR VecDestroy(&common_data.F);
219  CHKERR TSDestroy(&common_data.ts);
220  }
221 
222  }
223  CATCH_ERRORS;
224 
226  return 0;
227 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
BoneRemodeling::Remodeling::CommonData::nOremodellingBlock
bool nOremodellingBlock
Definition: Remodeling.hpp:164
BoneRemodeling::Remodeling::buildDM
MoFEMErrorCode buildDM()
Set problem and DM.
Definition: Remodeling.cpp:1832
BoneRemodeling::Remodeling::addFields
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
Definition: Remodeling.cpp:1607
BoneRemodeling::Remodeling::CommonData::feRhs
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
BoneRemodeling::Remodeling::CommonData::ts
TS ts
Time solver.
Definition: Remodeling.hpp:128
BoneRemodeling::Remodeling::CommonData::neumannForces
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
BoneRemodeling::Remodeling::CommonData::nodalForces
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
BoneRemodeling::Remodeling
Implementation of bone remodeling finite element.
Definition: Remodeling.hpp:36
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
DirichletSetFieldFromBlockWithFlags
Add boundary conditions form block set having 6 attributes.
Definition: DirichletBC.hpp:356
Remodeling.hpp
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
NODESET
@ NODESET
Definition: definitions.h:146
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::DMMoFEMTSSetIJacobian
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:853
BoneRemodeling::Remodeling::getParameters
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
Definition: Remodeling.cpp:1579
BoneRemodeling::Remodeling::CommonData::D
Vec D
Definition: Remodeling.hpp:120
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
BoneRemodeling::Remodeling::CommonData::elasticPtr
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
MoFEM::MeshsetsManager::setMeshsetFromFile
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
Definition: MeshsetsManager.cpp:788
DISPLACEMENTSET
@ DISPLACEMENTSET
Definition: definitions.h:150
BoneRemodeling::Remodeling::CommonData::edgeForces
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
help
static char help[]
Definition: bone_adaptation.cpp:35
BoneRemodeling
Definition: DensityMaps.hpp:27
BoneRemodeling::Remodeling::CommonData::F
Vec F
Definition: Remodeling.hpp:120
BoneRemodeling::Remodeling::Fe
Volume finite element.
Definition: Remodeling.hpp:41
BoneRemodeling::Remodeling::FePrePostProcessRhs
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Definition: Remodeling.hpp:53
BoneRemodeling::Remodeling::addMomentumFluxes
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
Definition: Remodeling.cpp:1792
ElasticMaterials.hpp
Elastic materials.
main
int main(int argc, char *argv[])
Definition: bone_adaptation.cpp:41
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
BoneRemodeling::Remodeling::solveDM
MoFEMErrorCode solveDM()
Solve problem set up in DM.
Definition: Remodeling.cpp:1863
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
BoneRemodeling::Remodeling::CommonData::feLhs
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
BoneRemodeling::Remodeling::CommonData
Definition: Remodeling.hpp:117
BoneRemodeling::Remodeling::addElements
MoFEMErrorCode addElements()
Set and add finite elements.
Definition: Remodeling.cpp:1694
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
std
Definition: enable_if.hpp:5
BoneRemodeling::Remodeling::CommonData::getParameters
MoFEMErrorCode getParameters()
Definition: Remodeling.cpp:62
BoneRemodeling::Remodeling::CommonData::dm
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::DMMoFEMTSSetIFunction
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:800
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
DirichletDisplacementBc::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: DirichletBC.hpp:76
BoneRemodeling::Remodeling::CommonData::A
Mat A
Definition: Remodeling.hpp:119
_IT_CUBITMESHSETS_FOR_LOOP_
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
Definition: MeshsetsManager.hpp:34
BoneRemodeling::Remodeling::FePrePostProcessLhs
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Definition: Remodeling.hpp:79
BoneRemodeling::Remodeling::CommonData::preProcLhs
boost::shared_ptr< FePrePostProcessLhs > preProcLhs
Definition: Remodeling.hpp:133
BoneRemodeling::Remodeling::CommonData::preProcRhs
boost::shared_ptr< FePrePostProcessRhs > preProcRhs
Definition: Remodeling.hpp:132