52 {
53
54
55 const char param_file[] = "param_file.petsc";
57
58
59 auto core_log = logging::core::get();
63 core_log->add_sink(
67
68 #ifdef ENABLE_PYTHON_BINDING
69 Py_Initialize();
70 np::initialize();
71 MOFEM_LOG(
"EP", Sev::inform) <<
"Python initialised";
72#else
73 MOFEM_LOG(
"EP", Sev::inform) <<
"Python NOT initialised";
74#endif
75
76 core_log->add_sink(
80
81 try {
82
84
85
86 PetscBool flg = PETSC_FALSE;
87 char mesh_file_name[255] = "";
89 255, &flg);
90 if (!flg) {
92 mesh_file_name, 255, &flg);
93 }
94 if (!flg) {
96 "Missing required mesh option: use -file_name or -my_file");
97 }
102 double time = 0;
104
105
106 MOFEM_LOG(
"EP", Sev::inform) <<
"Mesh file: " << mesh_file_name;
108 MOFEM_LOG(
"EP", Sev::inform) <<
"Time: " << time;
109
110
111 DMType dm_name = "DMMOFEM";
113 DMType dm_name_mg = "DMMOFEM_MG";
115
116
117 moab::Core moab_core;
118 moab::Interface &moab = moab_core;
119
120 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
121 if (pcomm == NULL)
122 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
123
124 PetscBool fully_distributed = PETSC_FALSE;
126 &fully_distributed, PETSC_NULLPTR);
127 if (fully_distributed) {
128 const char *option;
129 if (pcomm->proc_config().proc_size() == 1)
130 option = "";
131 else
132 option = "PARALLEL=READ_PART;"
133 "PARALLEL_RESOLVE_SHARED_ENTS;"
134 "PARTITION=PARALLEL_PARTITION";
135 CHKERR moab.load_file(mesh_file_name, 0, option);
136 } else {
139 }
140
141
142 MOFEM_LOG(
"EP", Sev::inform) <<
"Initialise MoFEM database";
145 MOFEM_LOG(
"EP", Sev::inform) <<
"Initialise MoFEM database <- done";
148
149
151
154 0, 3, bit_level0);
155
156
158
162 auto get_adj = [&](
Range ents,
int dim) {
164 CHKERR moab.get_adjacencies(ents, dim,
false, adj,
165 moab::Interface::UNION);
166 return adj;
167 };
169 *meshset_ptr,
172 2));
174 *meshset_ptr,
177 1));
179 *meshset_ptr,
182 0));
183
184
185 CHKERR ep.createCrackSurfaceMeshset();
186
187 CHKERR ep.getSpatialDispBc();
188 CHKERR ep.getSpatialRotationBc();
189 CHKERR ep.getSpatialTractionBc();
190 CHKERR ep.getSpatialTractionFreeBc();
191 CHKERR ep.getExternalStrain();
192 CHKERR ep.setBlockTagsOnSkin();
193
194 CHKERR ep.createExchangeVectors(Sev::inform);
195
196 CHKERR ep.projectMaterialTags();
197 CHKERR ep.addFields(*meshset_ptr);
198 CHKERR ep.projectGeometry(*meshset_ptr, 0.0);
199 CHKERR ep.addVolumeFiniteElement(*meshset_ptr);
200 CHKERR ep.addBoundaryFiniteElement(*meshset_ptr);
202
204 "stvenant_kirchhoff", "mooney_rivlin", "hencky",
205 "neo_hookean", "storakers", "meta"};
208 list_materials,
210 &choice_material, PETSC_NULLPTR);
213
214 auto material_model_names_dictionary = std::map<std::string, std::string>{
215 {"grad", "stretchH1AtPts"}, {"P", "PAtPts"}, {"P_dF", "PAtPts_du"}};
217 true;
218
219 switch (choice_material) {
221 MOFEM_LOG(
"EP", Sev::inform) <<
"St. Venant-Kirchhoff material model";
222 MOFEM_LOG(
"EP", Sev::inform) <<
"Meta Material model";
223 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
227 material_model_names_dictionary),
229 CHKERR physical_equations_ptr->getOptions(&m_field);
230 CHKERR physical_equations_ptr->recordTape();
231 CHKERR ep.addMaterial_Core(physical_equations_ptr);
232 };
234 MOFEM_LOG(
"EP", Sev::inform) <<
"Mooney-Rivlin material model";
235 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
239 material_model_names_dictionary),
241 CHKERR physical_equations_ptr->getOptions(&m_field);
242 CHKERR physical_equations_ptr->recordTape();
243 CHKERR ep.addMaterial_Core(physical_equations_ptr);
244 }; break;
246 MOFEM_LOG(
"EP", Sev::inform) <<
"Hencky material model";
247 CHKERR ep.addMaterial_Hencky(5., 0.25);
248 break;
250 MOFEM_LOG(
"EP", Sev::inform) <<
"Neo-Hookean material model";
251
252
253
254 CHKERR ep.addMaterial_HMHNeohookean(1.0, 2.0);
255 break;
257 MOFEM_LOG(
"EP", Sev::inform) <<
"Storakers material model";
258 CHKERR ep.addMaterial_HMHStorakers(0.5, 1.0, 1.0);
259 break;
261 MOFEM_LOG(
"EP", Sev::inform) <<
"Meta Material model";
262 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
263 auto physical_equations_ptr =
267 material_model_names_dictionary),
269 CHKERR physical_equations_ptr->getOptions(&m_field);
270 CHKERR physical_equations_ptr->recordTape();
271 CHKERR ep.addMaterial_Core(physical_equations_ptr);
272 }; break;
273 default:
275 break;
276 }
277
279 CHKERR ep.setElasticElementToTs(ep.dmElastic);
280
283 CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
284
285 PetscInt start_step = 0;
287 PetscViewer viewer;
289 FILE_MODE_READ, &viewer);
290 CHKERR VecLoad(x_elastic, viewer);
291 CHKERR PetscViewerDestroy(&viewer);
292 CHKERR VecGhostUpdateBegin(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
293 CHKERR VecGhostUpdateEnd(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
294
296 std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
297 std::smatch match;
298 if (std::regex_search(restart_file_str, match, restart_pattern)) {
299 start_step = std::stoi(match[1]);
300 } else {
302 "Restart file name must be in the format restart_##.dat");
303 }
304 }
305
306 auto ts_elastic =
createTS(PETSC_COMM_WORLD);
307 CHKERR TSSetType(ts_elastic, TSBEULER);
309 TSAdapt adapt;
310 CHKERR TSGetAdapt(ts_elastic, &adapt);
311 CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
312
315 MOFEM_LOG(
"EP", Sev::inform) <<
"Solver type: TS";
316 CHKERR TSSetTime(ts_elastic, time);
317 CHKERR TSSetStepNumber(ts_elastic, start_step);
318 CHKERR ep.solveElastic(ts_elastic, x_elastic);
319 break;
321 MOFEM_LOG(
"EP", Sev::inform) <<
"Solver type: Dynamic Relaxation";
323 SCATTER_REVERSE);
324 CHKERR ep.solveDynamicRelaxation(ts_elastic, x_elastic, start_step, time);
325 break;
327 MOFEM_LOG(
"EP", Sev::inform) <<
"Solver type: Cohesive";
329 SCATTER_REVERSE);
330 CHKERR ep.solveCohesiveCrackGrowth(ts_elastic, x_elastic, start_step,
331 time);
332 break;
334 MOFEM_LOG(
"EP", Sev::inform) <<
"Solver type: Load Factor";
336 SCATTER_REVERSE);
337 CHKERR ep.solveLoadFactor(ts_elastic, x_elastic, start_step, time);
338 break;
340 MOFEM_LOG(
"EP", Sev::inform) <<
"Solver type: Shape Optimization";
342 SCATTER_REVERSE);
343 CHKERR ep.solveSchapeOptimisation(ts_elastic, x_elastic, start_step,
344 time);
345 break;
348 << "Solver type: Test Topological Derivative";
350 SCATTER_REVERSE);
351 CHKERR ep.solveTestTopologicalDerivative(ts_elastic, x_elastic,
352 start_step, time);
353 break;
354 default:
356 "Unknown solver type");
357 break;
358 }
359
360 }
362
363
365
366#ifdef ENABLE_PYTHON_BINDING
367 if (Py_FinalizeEx() < 0) {
368 exit(120);
369 }
370#endif
371}
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
boost::shared_ptr< PhysicalEquations > createMatOpsPhysicalEquationsPtr(boost::shared_ptr< MatOpsData > mat_ops_data_ptr, int tag)
boost::shared_ptr< MatOpsData > createMatOpsDataPtr()
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createTS(MPI_Comm comm)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
static enum SolverType solverType
@ TestTopologicalDerivative
static enum MaterialModel materialModel
static bool useDeformationGradient
static Range getPartEntities(moab::Interface &moab, int part)
static MoFEMErrorCode loadFileRootProcAllRestDistributed(moab::Interface &moab, const char *file_name, int dim, LoadFileFun proc_skin_fun=defaultProcSkinFun, const char *options="PARALLEL=BCAST;PARTITION=")
Root proc has whole mesh, other procs only part of it.
virtual MPI_Comm & get_comm() const =0
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.
static MoFEMErrorCode setMeshFileFromJson()
Set -file_name from JSON before Core is available.
static MoFEMErrorCode addCanonicalAttributeNames(const std::vector< std::string > &names)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.