56 {
57
58 const char param_file[] = "param_file.petsc";
60
61 try {
62
63 {
64 PetscBool flg = PETSC_FALSE;
66 if (flg == PETSC_TRUE)
68 }
69
70 {
71 PetscBool flg = PETSC_FALSE;
73 PETSC_NULL);
74 if (flg == PETSC_TRUE)
76 }
77
78 PetscBool flg = PETSC_FALSE;
81 255, &flg);
82
83 const char *tan_list[] = {"tangent_internal_stress", "density_tangent",
84 "griffith_tangent"};
86 PetscBool isTangentTest;
88 LASTTAN, &choice_test, &isTangentTest);
89 moab::Core mb_instance;
90 moab::Interface &moab = mb_instance;
91
92 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
93 auto moab_comm_wrap =
94 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
95 if (pcomm == NULL)
96 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
97
99 if (flg) {
100 const char *option = "";
102 }
103 } else {
104 if (pcomm->rank() == 0) {
105
106 if (flg) {
107 moab::Core mb_instance_read;
108 moab::Interface &moab_read = mb_instance_read;
109 ParallelComm *pcomm =
111 if (pcomm == NULL)
112 pcomm = new ParallelComm(&moab_read, moab_comm_wrap->get_comm());
113
114 const char *option = "";
117 CHKERR moab_read.get_entities_by_handle(0, ents,
false);
118 CHKERR moab_read.get_entities_by_type(0, MBENTITYSET, ents,
false);
120 false, true);
121 }
122 } else {
123
126 true);
127 }
129 }
130
131
135
136 auto core_log = logging::core::get();
137 core_log->add_sink(
141 MOFEM_LOG_C(
"CP", Sev::inform,
"### Input parameter: -my_file %s",
143
144
146 cp.moabCommWorld = moab_comm_wrap;
147
150 << "File version " << cp.fileVersion.strVersion();
151
152
154
155
156 {
157 DMType dm_name = "MOFEM";
159 }
160
161 auto add_bits_tets = [&]() {
164 ->addToDatabaseBitRefLevelByType(MBTET,
BitRefLevel().set(),
167 };
168
170
172 std::string fn =
173 "mesh_after_cut_and_broadcast_" +
174 boost::lexical_cast<std::string>(m_field.
get_comm_rank()) +
".vtk";
175 CHKERR moab.write_file(fn.c_str(),
"VTK",
"");
176 }
177
179 auto first_number_string = [](std::string const &str) {
180 std::size_t
const n = str.find_first_of(
"0123456789");
181 if (
n != std::string::npos) {
182 std::size_t
const m = str.find_first_not_of(
"0123456789",
n);
183 return str.substr(
n,
m != std::string::npos ?
m -
n :
m);
184 }
185 return std::string();
186 };
187 std::string str_number = first_number_string(std::string(
mesh_file_name));
188
189 if (!str_number.empty())
190 cp.startStep = atoi(str_number.c_str()) + 1;
191 }
192
193
194 if (cp.doCutMesh) {
195
197
200 ->getAllEntitiesNotInDatabase(ents);
202 CHKERR m_field.
get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
203 false);
204 for (
auto m : meshsets)
207
208 string cutting_surface_name =
209 "cutting_surface_" +
210 boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
211
214 cp.crackAccelerationFactor, cutting_surface_name,
QUIET,
216 else
219
222
223 } else {
224
225 string cutting_surface_name =
226 "cutting_surface_" +
227 boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
228
231 CHKERR cp.getInterface<
CPMeshCut>()->copySurface(cutting_surface_name);
233 CHKERR moab.get_entities_by_type(0, MBTET, vol,
false);
236 }
237
240
241 } else {
242
244 CHKERR moab.get_entities_by_type(0, MBTET, vol,
false);
248
249
250
251
252 if (!cp.doCutMesh && !cp.propagateCrack && cp.doElasticWithoutCrack) {
253 BitRefLevel bit1 = cp.mapBitLevel[
"spatial_domain"];
256 } else {
259 }
260 }
261
263
266 }
267
268 std::vector<int> surface_ids;
270 std::vector<std::string> fe_surf_proj_list;
271 fe_surf_proj_list.push_back("SURFACE");
272
273 CHKERR cp.createProblemDataStructures(surface_ids,
QUIET,
275
276 if (isTangentTest) {
277
278 MOFEM_LOG(
"CP", Sev::verbose) <<
"Testing tangent variant";
279
280
281 CHKERR cp.setMaterialPositionFromCoords();
282 CHKERR cp.setSpatialPositionFromCoords();
283 if (cp.addSingularity == PETSC_TRUE) {
284 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
285 }
286
294 CHKERR cp.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
295 dm_crack_propagation, dm_material_forces,
296 dm_surface_projection, dm_crack_srf_area, surface_ids,
297 fe_surf_proj_list);
298
299
301 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
302 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
303 false);
304 CHKERR cp.assembleCouplingForcesDM(PETSC_NULL,
VERBOSE,
false);
305
306
307 CHKERR cp.setMaterialPositionFromCoords();
308 CHKERR cp.setSpatialPositionFromCoords();
309 if (cp.addSingularity == PETSC_TRUE)
310 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
311
312 switch (choice_test) {
314 CHKERR cp.testJacobians(cp.mapBitLevel[
"material_domain"],
316 break;
318 CHKERR cp.testJacobians(cp.mapBitLevel[
"material_domain"],
320 break;
322 CHKERR cp.testJacobians(cp.mapBitLevel[
"material_domain"],
324 break;
325 default:
326 break;
327 }
328
329 } else if (cp.propagateCrack == PETSC_TRUE) {
330
331 MOFEM_LOG(
"CP", Sev::verbose) <<
"Crack propagation variant is used";
332
333
334 CHKERR cp.setMaterialPositionFromCoords();
335 CHKERR cp.setSpatialPositionFromCoords();
336 if (cp.addSingularity == PETSC_TRUE)
337 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
338
346 CHKERR cp.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
347 dm_crack_propagation, dm_material_forces,
348 dm_surface_projection, dm_crack_srf_area, surface_ids,
349 fe_surf_proj_list);
350
351 CHKERR cp.setSpatialPositionFromCoords();
352 if (cp.addSingularity == PETSC_TRUE) {
353 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
354 }
356 cp.getLoadScale() = cp.getArcCtx()->getFieldData();
357
358
359 const auto load_factor = std::abs(cp.getLoadScale());
360 MOFEM_LOG(
"CPSolverSync", Sev::noisy) <<
"Lambda factor " << load_factor;
362 if (cp.solveEigenStressProblem) {
363 MOFEM_LOG(
"CP", Sev::verbose) <<
"Solve Eigen ELASTIC problem";
365 dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
366 CHKERR cp.postProcessDM(dm_eigen_elastic, -2,
"ELASTIC",
true);
367 }
368 MOFEM_LOG(
"CP", Sev::noisy) <<
"Solve ELASTIC problem";
369 cp.getLoadScale() = load_factor;
371 << "Reset lambda factor " << cp.getLoadScale();
374 dm_elastic, cp.getArcCtx()->x0, load_factor);
375 CHKERR cp.postProcessDM(dm_elastic, -1,
"ELASTIC",
true);
376
377
379 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
380 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
381 false);
382 CHKERR cp.assembleCouplingForcesDM(PETSC_NULL,
VERBOSE,
false);
383
384 CHKERR cp.calculateElasticEnergy(dm_elastic);
386 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
387 surface_ids);
388
389 if (cp.doCutMesh == PETSC_TRUE) {
390 MOFEM_LOG(
"CP", Sev::verbose) <<
"Propagate crack and cut mesh";
392 dm_crack_propagation, dm_elastic, dm_eigen_elastic, dm_material,
393 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
395 } else {
396 MOFEM_LOG(
"CP", Sev::verbose) <<
"Propagate crack (not cut mesh)";
398 dm_crack_propagation, dm_elastic, dm_material, dm_material_forces,
399 dm_surface_projection, dm_crack_srf_area, surface_ids);
400 }
401
402 } else {
403
404 MOFEM_LOG(
"CP", Sev::verbose) <<
"Do not propagate crack";
405
406
407 CHKERR cp.setMaterialPositionFromCoords();
408 CHKERR cp.setSpatialPositionFromCoords();
409 if (cp.addSingularity == PETSC_TRUE)
410 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
411
419 CHKERR cp.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
420 dm_crack_propagation, dm_material_forces,
421 dm_surface_projection, dm_crack_srf_area, surface_ids,
422 fe_surf_proj_list);
423
426 "ANALITICAL_DISP"))
427 CHKERR cp.createBcDM(dm_bc,
"BC_PROBLEM",
428 cp.mapBitLevel["spatial_domain"],
430
431
432 if (cp.solveEigenStressProblem) {
433 MOFEM_LOG(
"CP", Sev::verbose) <<
"Solve Eigen ELASTIC problem";
435 dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
436 CHKERR cp.postProcessDM(dm_eigen_elastic, -2,
"ELASTIC",
true);
437 }
438
439 MOFEM_LOG(
"CP", Sev::noisy) <<
"Solve ELASTIC problem";
441 dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
442
443
445 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
446 if (cp.doCutMesh)
447 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
448 false);
449
451 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
452 surface_ids);
453 CHKERR cp.calculateElasticEnergy(dm_elastic);
454
455
456 CHKERR cp.postProcessDM(dm_elastic, 0,
"ELASTIC",
true);
457
458 CHKERR cp.savePositionsOnCrackFrontDM(dm_elastic, PETSC_NULL, 2,
false);
459 if (!cp.doElasticWithoutCrack)
460 CHKERR cp.writeCrackFont(cp.mapBitLevel[
"spatial_domain"], 0);
461
462
463 CHKERR cp.tetsingReleaseEnergyCalculation();
464 }
465 }
467
469 return 0;
470}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
#define CATCH_ERRORS
Catch errors.
#define BITREFLEVEL_SIZE
max number of refinements
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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.
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
const double n
refractive index of diffusive medium
MoFEMErrorCode broadcast_entities(moab::Interface &moab, moab::Interface &moab_tmp, const int from_proc, Range &entities, const bool adjacencies, const bool tags)
MoFEMErrorCode clean_pcomms(moab::Interface &moab, boost::shared_ptr< WrapMPIComm > moab_comm_wrap)
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 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)
FTensor::Index< 'm', 3 > m
int getSkinOfTheBodyId() const
const std::string & getCutSurfMeshName() const
static bool parallelReadAndBroadcast
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_field(const std::string field_name, int verb=DEFAULT_VERBOSITY)=0
build field by name
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() 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 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.
Interface for managing meshsets containing materials and boundary conditions.
static bool brodcastMeshsets
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.