v0.14.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 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>(
104 new Remodeling::FePrePostProcessRhs());
105 common_data.preProcLhs = boost::shared_ptr<Remodeling::FePrePostProcessLhs>(
106 new 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;
121 for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
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
225 MoFEM::Core::Finalize();
226 return 0;
227}
help
static char help[]
Definition: bone_adaptation.cpp:35
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
NODESET
@ NODESET
Definition: definitions.h:146
DISPLACEMENTSET
@ DISPLACEMENTSET
Definition: definitions.h:150
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
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:786
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:509
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:839
_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:23
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
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::dm
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
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::Fe
Volume finite element.
Definition: Remodeling.hpp:41
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
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:112
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:725
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce 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.9.5 and hosted at University of Glasgow