v0.15.0
Loading...
Searching...
No Matches
Functions | Variables
bone_adaptation.cpp File Reference
#include <boost/program_options.hpp>
#include <BasicFiniteElements.hpp>
#include <ElasticMaterials.hpp>
#include <Remodeling.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 41 of file bone_adaptation.cpp.

41 {
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 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 PetscOptionsEnd();
60
61 moab::Core mb_instance;
62 moab::Interface &moab = mb_instance;
63
64 // Read parameters from line command
65 if (flg_file != PETSC_TRUE) {
66 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
67 }
68 const char *option;
69
70 // option = "";
71 option = "PARALLEL=READ_PART;"
72 "PARALLEL_RESOLVE_SHARED_ENTS;"
73 "PARTITION=PARALLEL_PARTITION;";
74 CHKERR moab.load_file(mesh_file_name, 0, option);
75 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
76 if (pcomm == NULL)
77 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
78
79 MoFEM::Core core(moab);
80 MoFEM::Interface &m_field = core;
81
82 // Add meshsets with material and boundary conditions
83 MeshsetsManager *meshsets_manager_ptr;
84 CHKERR m_field.getInterface<MeshsetsManager>(meshsets_manager_ptr);
85 CHKERR meshsets_manager_ptr->setMeshsetFromFile();
86
87 PetscPrintf(PETSC_COMM_WORLD, "Read meshset. Added meshsets for bc.cfg\n");
88 for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
89 PetscPrintf(PETSC_COMM_WORLD, "%s",
90 static_cast<std::ostringstream &>(
91 std::ostringstream().seekp(0) << *it << endl)
92 .str()
93 .c_str());
94 cerr << *it << endl;
95 }
96
97 Remodeling::CommonData common_data;
98 common_data.feRhs =
99 boost::shared_ptr<Remodeling::Fe>(new Remodeling::Fe(m_field));
100 common_data.feLhs =
101 boost::shared_ptr<Remodeling::Fe>(new Remodeling::Fe(m_field));
102 common_data.preProcRhs = boost::shared_ptr<Remodeling::FePrePostProcessRhs>(
103 new Remodeling::FePrePostProcessRhs());
104 common_data.preProcLhs = boost::shared_ptr<Remodeling::FePrePostProcessLhs>(
105 new Remodeling::FePrePostProcessLhs());
106
107 CHKERR common_data.getParameters();
108
109 Remodeling remodelling(m_field, common_data);
110
111 // Remodeling
112 if (!flg_elastic_only) {
113 remodelling.getParameters();
114 remodelling.addFields();
115 remodelling.addMomentumFluxes();
116 remodelling.addElements();
117 remodelling.buildDM();
118
119 bool flag_cubit_disp = false;
120 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
121 m_field, NODESET | DISPLACEMENTSET, it)) {
122 flag_cubit_disp = true;
123 }
124
125 boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
126 // if normally defined boundary conditions are not found, try to use
127 // DISPLACEMENT blockset
128 if (!flag_cubit_disp) {
129 dirihlet_bc_ptr =
130 boost::shared_ptr<FEMethod>(new DirichletSetFieldFromBlockWithFlags(
131 m_field, "DISPLACEMENTS", "DISPLACEMENT"));
132 dynamic_cast<DirichletSetFieldFromBlockWithFlags &>(*dirihlet_bc_ptr)
133 .methodsOp.push_back(new TimeForceScale("-my_load_history", false));
134 } else {
135 dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
136 new DirichletDisplacementBc(m_field, "DISPLACEMENTS"));
137 dynamic_cast<DirichletDisplacementBc &>(*dirihlet_bc_ptr)
138 .methodsOp.push_back(new TimeForceScale("-my_load_history", false));
139 }
140
141 // Add elements Time Solver via DM interface.
142 // This part is to calculate residual vector.
143 CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL,
144 common_data.preProcRhs.get(), NULL);
145 CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL,
146 dirihlet_bc_ptr.get(), NULL);
147 CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING",
148 common_data.feRhs.get(), NULL, NULL);
149 if (common_data.nOremodellingBlock) {
150 CHKERR DMMoFEMTSSetIFunction(common_data.dm, "ELASTIC",
151 &(common_data.elasticPtr->getLoopFeRhs()),
152 NULL, NULL);
153 }
154
155 // Add surface forces
156 for (boost::ptr_map<string, NeummanForcesSurface>::iterator fit =
157 common_data.neumannForces.begin();
158 fit != common_data.neumannForces.end(); fit++) {
159 CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
160 &fit->second->getLoopFe(), NULL, NULL);
161 }
162
163 // Add edge forces
164 for (boost::ptr_map<std::string, EdgeForce>::iterator fit =
165 common_data.edgeForces.begin();
166 fit != common_data.edgeForces.end(); fit++) {
167 CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
168 &fit->second->getLoopFe(), NULL, NULL);
169 }
170
171 // Add nodal forces
172 for (boost::ptr_map<std::string, NodalForce>::iterator fit =
173 common_data.nodalForces.begin();
174 fit != common_data.nodalForces.end(); fit++) {
175 CHKERR DMMoFEMTSSetIFunction(common_data.dm, fit->first.c_str(),
176 &fit->second->getLoopFe(), NULL, NULL);
177 }
178
179 CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL, NULL,
180 dirihlet_bc_ptr.get());
181 CHKERR DMMoFEMTSSetIFunction(common_data.dm, "FE_REMODELLING", NULL, NULL,
182 common_data.preProcRhs.get());
183
184 // Add elements Time Solver via DM interface.
185 // This part is to calculate tangent matrix.
186 CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL,
187 common_data.preProcLhs.get(), NULL);
188 CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL,
189 dirihlet_bc_ptr.get(), NULL);
190 CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING",
191 common_data.feLhs.get(), NULL, NULL);
192 if (common_data.nOremodellingBlock) {
193 CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "ELASTIC",
194 &common_data.elasticPtr->getLoopFeLhs(),
195 NULL, NULL);
196 }
197 CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL, NULL,
198 dirihlet_bc_ptr.get());
199 CHKERR DMMoFEMTSSetIJacobian(common_data.dm, "FE_REMODELLING", NULL, NULL,
200 common_data.preProcLhs.get());
201
202 CHKERR DMCreateMatrix(common_data.dm, &common_data.A);
203 CHKERR DMCreateGlobalVector(common_data.dm, &common_data.D);
204 CHKERR DMCreateGlobalVector(common_data.dm, &common_data.F);
205
206 // Initialize vectors with initial data from mesh
207 CHKERR DMoFEMMeshToLocalVector(common_data.dm, common_data.D,
208 INSERT_VALUES, SCATTER_FORWARD);
209
210 CHKERR TSCreate(PETSC_COMM_WORLD, &common_data.ts);
211 CHKERR TSSetType(common_data.ts, TSBEULER);
212
213 CHKERR remodelling.solveDM();
214
215 CHKERR MatDestroy(&common_data.A);
216 CHKERR VecDestroy(&common_data.D);
217 CHKERR VecDestroy(&common_data.F);
218 CHKERR TSDestroy(&common_data.ts);
219 }
220
221 }
222 CATCH_ERRORS;
223
224 MoFEM::Core::Finalize();
225 return 0;
226}
help
static char help[]
Definition bone_adaptation.cpp:35
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
NODESET
@ NODESET
Definition definitions.h:159
DISPLACEMENTSET
@ DISPLACEMENTSET
Definition definitions.h:163
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
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:790
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:514
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:843
_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
_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
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
BoneRemodeling::Remodeling::CommonData::preProcRhs
boost::shared_ptr< FePrePostProcessRhs > preProcRhs
Definition Remodeling.hpp:132
BoneRemodeling::Remodeling::CommonData::feLhs
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition Remodeling.hpp:130
BoneRemodeling::Remodeling::CommonData::feRhs
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition Remodeling.hpp:131
BoneRemodeling::Remodeling::CommonData::neumannForces
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition Remodeling.hpp:136
BoneRemodeling::Remodeling::CommonData::edgeForces
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition Remodeling.hpp:138
BoneRemodeling::Remodeling::CommonData::preProcLhs
boost::shared_ptr< FePrePostProcessLhs > preProcLhs
Definition Remodeling.hpp:133
BoneRemodeling::Remodeling::CommonData::nodalForces
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition Remodeling.hpp:137
BoneRemodeling::Remodeling::CommonData::elasticPtr
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition Remodeling.hpp:140
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::FePrePostProcessRhs
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
Definition Remodeling.hpp:53
BoneRemodeling::Remodeling::Fe
Volume finite element.
Definition Remodeling.hpp:41
BoneRemodeling::Remodeling
Implementation of bone remodeling finite element.
Definition Remodeling.hpp:36
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition DirichletBC.hpp:57
DirichletDisplacementBc::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition DirichletBC.hpp:76
DirichletSetFieldFromBlockWithFlags
Add boundary conditions form block set having 6 attributes.
Definition DirichletBC.hpp:356
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:104
MoFEM::MeshsetsManager::setMeshsetFromFile
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
Definition MeshsetsManager.cpp:802
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
TimeForceScale
Force scale operator for reading two columns.
Definition TimeForceScale.hpp:18

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 35 of file bone_adaptation.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow