v0.14.0
Loading...
Searching...
No Matches
crack_propagation.cpp File Reference
#include <tetgen.h>
#include <MoFEM.hpp>
#include <BasicFiniteElements.hpp>
#include <Mortar.hpp>
#include <NeoHookean.hpp>
#include <Hooke.hpp>
#include <CrackFrontElement.hpp>
#include <ComplexConstArea.hpp>
#include <MWLS.hpp>
#include <GriffithForceElement.hpp>
#include <VolumeLengthQuality.hpp>
#include <CrackPropagation.hpp>
#include <CPSolvers.hpp>
#include <CPMeshCut.hpp>
#include <AnalyticalFun.hpp>

Go to the source code of this file.

Macros

#define BOOST_UBLAS_NDEBUG   1
 

Functions

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

Variables

static char help []
 

Macro Definition Documentation

◆ BOOST_UBLAS_NDEBUG

#define BOOST_UBLAS_NDEBUG   1

Definition at line 31 of file crack_propagation.cpp.

Function Documentation

◆ main()

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

Definition at line 56 of file crack_propagation.cpp.

56 {
57
58 const char param_file[] = "param_file.petsc";
59 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
60
61 try {
62
63 {
64 PetscBool flg = PETSC_FALSE;
65 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-debug", &flg, PETSC_NULL);
66 if (flg == PETSC_TRUE)
68 }
69
70 {
71 PetscBool flg = PETSC_FALSE;
72 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-read_and_brodcast", &flg,
73 PETSC_NULL);
74 if (flg == PETSC_TRUE)
76 }
77
78 PetscBool flg = PETSC_FALSE;
79 char mesh_file_name[255];
80 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
81 255, &flg);
82
83 const char *tan_list[] = {"tangent_internal_stress", "density_tangent",
84 "griffith_tangent"};
85 int choice_test = MWLS_STRESS_TAN;
86 PetscBool isTangentTest;
87 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-testing_mwls", tan_list,
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
98 if (pcomm->size() == 1 || !CrackPropagation::parallelReadAndBroadcast) {
99 if (flg) {
100 const char *option = "";
101 CHKERR moab.load_file(mesh_file_name, 0, option);
102 }
103 } else {
104 if (pcomm->rank() == 0) {
105 // Read mesh file
106 if (flg) {
107 moab::Core mb_instance_read;
108 moab::Interface &moab_read = mb_instance_read;
109 ParallelComm *pcomm =
110 ParallelComm::get_pcomm(&moab_read, MYPCOMM_INDEX);
111 if (pcomm == NULL)
112 pcomm = new ParallelComm(&moab_read, moab_comm_wrap->get_comm());
113
114 const char *option = "";
115 CHKERR moab_read.load_file(mesh_file_name, 0, option);
116 Range ents;
117 CHKERR moab_read.get_entities_by_handle(0, ents, false);
118 CHKERR moab_read.get_entities_by_type(0, MBENTITYSET, ents, false);
119 CHKERR FractureMechanics::broadcast_entities(moab, moab_read, 0, ents,
120 false, true);
121 }
122 } else {
123
124 Range ents;
125 CHKERR FractureMechanics::broadcast_entities(moab, moab, 0, ents, false,
126 true);
127 }
128 CHKERR FractureMechanics::clean_pcomms(moab, moab_comm_wrap);
129 }
130
131 // Create mofem database
133 MoFEM::Core core(moab, PETSC_COMM_WORLD, VERBOSE);
134 MoFEM::Interface &m_field = core;
135
136 auto core_log = logging::core::get();
137 core_log->add_sink(
139 LogManager::setLog("CP");
140 MOFEM_LOG_TAG("CP", "cp");
141 MOFEM_LOG_C("CP", Sev::inform, "### Input parameter: -my_file %s",
143
144 // Create data structure for crack propagation
145 CrackPropagation cp(m_field, 3, 2);
146 cp.moabCommWorld = moab_comm_wrap;
147
148 CHKERR MoFEM::Interface::getFileVersion(moab, cp.fileVersion);
149 MOFEM_LOG("CP", Sev::inform)
150 << "File version " << cp.fileVersion.strVersion();
151
152 // Configure meshes
153 CHKERR cp.getOptions();
154
155 // Register MOFEM discrete data managers
156 {
157 DMType dm_name = "MOFEM";
158 CHKERR DMRegister_MoFEM(dm_name);
159 }
160
161 auto add_bits_tets = [&]() {
164 ->addToDatabaseBitRefLevelByType(MBTET, BitRefLevel().set(),
165 BitRefLevel().set());
167 };
168
169 CHKERR add_bits_tets();
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
178 if (string(mesh_file_name).compare(0, 8, "restart_") == 0) {
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 // create bit levels, fields and finite elements structures in database
194 if (cp.doCutMesh) {
195
196 if (string(mesh_file_name).compare(0, 8, "restart_") == 0) {
197
198 Range ents;
200 ->getAllEntitiesNotInDatabase(ents);
201 Range meshsets;
202 CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
203 false);
204 for (auto m : meshsets)
205 CHKERR m_field.get_moab().remove_entities(m, ents);
206 CHKERR m_field.get_moab().delete_entities(ents);
207
208 string cutting_surface_name =
209 "cutting_surface_" +
210 boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
211
212 if (m_field.get_comm_rank() == 0)
213 CHKERR cp.getInterface<CPMeshCut>()->rebuildCrackSurface(
214 cp.crackAccelerationFactor, cutting_surface_name, QUIET,
216 else
217 CHKERR cp.getInterface<CPMeshCut>()->rebuildCrackSurface(
218 cp.crackAccelerationFactor, "", QUIET, CrackPropagation::debug);
219
220 CHKERR cp.getInterface<CPMeshCut>()->refineMesh(
221 nullptr, false, QUIET, CrackPropagation::debug);
222
223 } else {
224
225 string cutting_surface_name =
226 "cutting_surface_" +
227 boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
228
229 if (!cp.getInterface<CPMeshCut>()->getCutSurfMeshName().empty())
230 CHKERR cp.getInterface<CPMeshCut>()->setCutSurfaceFromFile();
231 CHKERR cp.getInterface<CPMeshCut>()->copySurface(cutting_surface_name);
232 Range vol;
233 CHKERR moab.get_entities_by_type(0, MBTET, vol, false);
234 CHKERR cp.getInterface<CPMeshCut>()->refineMesh(
235 &vol, true, QUIET, CrackPropagation::debug);
236 }
237
238 CHKERR cp.getInterface<CPMeshCut>()->cutRefineAndSplit(
240
241 } else {
242
243 Range vol;
244 CHKERR moab.get_entities_by_type(0, MBTET, vol, false);
245 BitRefLevel bit0 = cp.mapBitLevel["mesh_cut"];
246 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
247 vol, bit0.set(BITREFLEVEL_SIZE - 1), true, CrackPropagation::debug);
248
249 // intact mesh will be saved in "spatial_domain" mapBitLevel to solve
250 // the elastic problem in the absence of the crack (restart files won't
251 // work)
252 if (!cp.doCutMesh && !cp.propagateCrack && cp.doElasticWithoutCrack) {
253 BitRefLevel bit1 = cp.mapBitLevel["spatial_domain"];
254 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
255 vol, bit1.set(BITREFLEVEL_SIZE - 1), true, CrackPropagation::debug);
256 } else {
257 CHKERR cp.getInterface<CPMeshCut>()->refineAndSplit(
259 }
260 }
261
262 if (string(mesh_file_name).compare(0, 8, "restart_") == 0) {
263 // restart code from file
264 CHKERR m_field.build_field("LAMBDA_ARC_LENGTH");
265 CHKERR m_field.build_finite_elements("ARC_LENGTH");
266 }
267
268 std::vector<int> surface_ids;
269 surface_ids.push_back(cp.getInterface<CPMeshCut>()->getSkinOfTheBodyId());
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 // set inital coordinates
281 CHKERR cp.setMaterialPositionFromCoords();
282 CHKERR cp.setSpatialPositionFromCoords();
283 if (cp.addSingularity == PETSC_TRUE) {
284 CHKERR cp.setSingularDofs("SPATIAL_POSITION");
285 }
286
287 SmartPetscObj<DM> dm_elastic;
288 SmartPetscObj<DM> dm_eigen_elastic;
289 SmartPetscObj<DM> dm_material;
290 SmartPetscObj<DM> dm_crack_propagation;
291 SmartPetscObj<DM> dm_material_forces;
292 SmartPetscObj<DM> dm_surface_projection;
293 SmartPetscObj<DM> dm_crack_srf_area;
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 // set finite element instances and user data operators on instances
300 CHKERR cp.assembleElasticDM(VERBOSE, false);
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 // set inital coordinates
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) {
313 case MWLS_STRESS_TAN:
314 CHKERR cp.testJacobians(cp.mapBitLevel["material_domain"],
316 break;
317 case MWLS_DENSITY_TAN:
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 // set inital coordinates
334 CHKERR cp.setMaterialPositionFromCoords();
335 CHKERR cp.setSpatialPositionFromCoords();
336 if (cp.addSingularity == PETSC_TRUE)
337 CHKERR cp.setSingularDofs("SPATIAL_POSITION");
338
339 SmartPetscObj<DM> dm_elastic;
340 SmartPetscObj<DM> dm_eigen_elastic;
341 SmartPetscObj<DM> dm_material;
342 SmartPetscObj<DM> dm_crack_propagation;
343 SmartPetscObj<DM> dm_material_forces;
344 SmartPetscObj<DM> dm_surface_projection;
345 SmartPetscObj<DM> dm_crack_srf_area;
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 }
355 if (string(mesh_file_name).compare(0, 8, "restart_") == 0)
356 cp.getLoadScale() = cp.getArcCtx()->getFieldData();
357
358 // Solve eigen problem
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";
364 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
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;
370 MOFEM_LOG("CPSolverSync", Sev::noisy)
371 << "Reset lambda factor " << cp.getLoadScale();
373 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
374 dm_elastic, cp.getArcCtx()->x0, load_factor);
375 CHKERR cp.postProcessDM(dm_elastic, -1, "ELASTIC", true);
376
377 // set finite element instances and user data operators on instances
378 CHKERR cp.assembleElasticDM(VERBOSE, false);
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);
385 CHKERR cp.getInterface<CPSolvers>()->calculateGc(
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";
391 CHKERR cp.getInterface<CPSolvers>()->solveCutMesh(
392 dm_crack_propagation, dm_elastic, dm_eigen_elastic, dm_material,
393 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
394 surface_ids, CrackPropagation::debug);
395 } else {
396 MOFEM_LOG("CP", Sev::verbose) << "Propagate crack (not cut mesh)";
397 CHKERR cp.getInterface<CPSolvers>()->solvePropagation(
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 // set inital coordinates
407 CHKERR cp.setMaterialPositionFromCoords();
408 CHKERR cp.setSpatialPositionFromCoords();
409 if (cp.addSingularity == PETSC_TRUE)
410 CHKERR cp.setSingularDofs("SPATIAL_POSITION");
411
412 SmartPetscObj<DM> dm_elastic;
413 SmartPetscObj<DM> dm_eigen_elastic;
414 SmartPetscObj<DM> dm_material;
415 SmartPetscObj<DM> dm_crack_propagation;
416 SmartPetscObj<DM> dm_material_forces;
417 SmartPetscObj<DM> dm_surface_projection;
418 SmartPetscObj<DM> dm_crack_srf_area;
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
424 SmartPetscObj<DM> dm_bc;
426 "ANALITICAL_DISP"))
427 CHKERR cp.createBcDM(dm_bc, "BC_PROBLEM",
428 cp.mapBitLevel["spatial_domain"],
429 BitRefLevel().set());
430
431 // Solve eigen problem
432 if (cp.solveEigenStressProblem) {
433 MOFEM_LOG("CP", Sev::verbose) << "Solve Eigen ELASTIC problem";
434 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
435 dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
436 CHKERR cp.postProcessDM(dm_eigen_elastic, -2, "ELASTIC", true);
437 }
438 // set finite element instances and user data operators on instances
439 MOFEM_LOG("CP", Sev::noisy) << "Solve ELASTIC problem";
440 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
441 dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
442
443 // set finite element instances and user data operators on instances
444 CHKERR cp.assembleElasticDM(VERBOSE, false);
445 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
446 if (cp.doCutMesh)
447 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE,
448 false);
449
450 CHKERR cp.getInterface<CPSolvers>()->calculateGc(
451 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
452 surface_ids);
453 CHKERR cp.calculateElasticEnergy(dm_elastic);
454
455 // post-processing
456 CHKERR cp.postProcessDM(dm_elastic, 0, "ELASTIC", true);
457 // set vector values from field data
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 // Run tests
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,...)
static char help[]
@ QUIET
@ VERBOSE
#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.
Definition DMMoFEM.cpp:43
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
char mesh_file_name[255]
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.
Definition Types.hpp:40
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
const std::string & getCutSurfMeshName() const
Managing BitRefLevels.
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
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
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.
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.

Variable Documentation

◆ help

char help[]
static
Initial value:
= "Calculate crack release energy and crack propagation"
"\n\n"

Definition at line 43 of file crack_propagation.cpp.