v0.14.0
Loading...
Searching...
No Matches
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>
26using namespace std;
27namespace po = boost::program_options;
28
30#include <ElasticMaterials.hpp>
31#include <Remodeling.hpp>
32
33using namespace BoneRemodeling;
34
35static char help[] = "\n";
36
37// #ifndef WITH_METAIO
38// #error "You need compile users modules with MetaIO -DWITH_METAIO=1"
39// #endif
40
41int 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 }
224
226 return 0;
227}
Elastic materials.
int main()
static char help[]
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ NODESET
@ DISPLACEMENTSET
#define CHKERR
Inline error check.
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:786
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:509
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:839
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
boost::shared_ptr< FePrePostProcessRhs > preProcRhs
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
boost::shared_ptr< FePrePostProcessLhs > preProcLhs
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Volume finite element.
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Implementation of bone remodeling finite element.
MoFEMErrorCode addElements()
Set and add finite elements.
MoFEMErrorCode solveDM()
Solve problem set up in DM.
MoFEMErrorCode buildDM()
Set problem and DM.
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
Set Dirichlet boundary conditions on displacements.
boost::ptr_vector< MethodForForceScaling > methodsOp
Add boundary conditions form block set having 6 attributes.
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.