41 {
42
44
45 try {
46
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",
54
56 "-elastic_only", "is used for testing, calculates only elastic problem",
57 "", PETSC_FALSE, &flg_elastic_only, PETSC_NULL);
58
59 ierr = PetscOptionsEnd();
61
62 moab::Core mb_instance;
63 moab::Interface &moab = mb_instance;
64
65
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
72 option = "PARALLEL=READ_PART;"
73 "PARALLEL_RESOLVE_SHARED_ENTS;"
74 "PARTITION=PARALLEL_PARTITION;";
76 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
77 if (pcomm == NULL)
78 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
79
82
83
87
88 PetscPrintf(PETSC_COMM_WORLD, "Read meshset. Added meshsets for bc.cfg\n");
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
103 common_data.
preProcRhs = boost::shared_ptr<Remodeling::FePrePostProcessRhs>(
105 common_data.
preProcLhs = boost::shared_ptr<Remodeling::FePrePostProcessLhs>(
107
109
111
112
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;
123 flag_cubit_disp = true;
124 }
125
126 boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
127
128
129 if (!flag_cubit_disp) {
130 dirihlet_bc_ptr =
132 m_field, "DISPLACEMENTS", "DISPLACEMENT"));
135 } else {
136 dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
140 }
141
142
143
147 dirihlet_bc_ptr.get(), NULL);
149 common_data.
feRhs.get(), NULL, NULL);
153 NULL, NULL);
154 }
155
156
157 for (boost::ptr_map<string, NeummanForcesSurface>::iterator fit =
161 &fit->second->getLoopFe(), NULL, NULL);
162 }
163
164
165 for (boost::ptr_map<std::string, EdgeForce>::iterator fit =
169 &fit->second->getLoopFe(), NULL, NULL);
170 }
171
172
173 for (boost::ptr_map<std::string, NodalForce>::iterator fit =
177 &fit->second->getLoopFe(), NULL, NULL);
178 }
179
181 dirihlet_bc_ptr.get());
184
185
186
190 dirihlet_bc_ptr.get(), NULL);
192 common_data.
feLhs.get(), NULL, NULL);
196 NULL, NULL);
197 }
199 dirihlet_bc_ptr.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
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);
220 }
221
222 }
224
226 return 0;
227}
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#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
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
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
#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.
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.
MoFEMErrorCode getParameters()
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
DM dm
Discretization manager.
boost::shared_ptr< FePrePostProcessLhs > preProcLhs
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
boost::shared_ptr< NonlinearElasticElement > elasticPtr
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.
Set Dirichlet boundary conditions on displacements.
boost::ptr_vector< MethodForForceScaling > methodsOp
Add boundary conditions form block set having 6 attributes.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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.