58 const char param_file[] =
"param_file.petsc";
64 PetscBool flg = PETSC_FALSE;
66 if (flg == PETSC_TRUE)
71 PetscBool flg = PETSC_FALSE;
74 if (flg == PETSC_TRUE)
75 CrackPropagation::parallelReadAndBroadcast =
true;
78 PetscBool flg = PETSC_FALSE;
83 const char *tan_list[] = {
"tangent_internal_stress",
"density_tangent",
86 PetscBool isTangentTest;
88 LASTTAN, &choice_test, &isTangentTest);
92 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
94 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD,
false);
96 pcomm =
new ParallelComm(&moab, moab_comm_wrap->get_comm());
98 if (pcomm->size() == 1 || !CrackPropagation::parallelReadAndBroadcast) {
100 const char *option =
"";
104 if (pcomm->rank() == 0) {
109 ParallelComm *pcomm =
112 pcomm =
new ParallelComm(&moab_read, moab_comm_wrap->get_comm());
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);
132 MeshsetsManager::brodcastMeshsets =
false;
136 auto core_log = logging::core::get();
138 LogManager::createSink(LogManager::getStrmWorld(),
"CP"));
139 LogManager::setLog(
"CP");
141 MOFEM_LOG_C(
"CP", Sev::inform,
"### Input parameter: -my_file %s",
146 cp.moabCommWorld = moab_comm_wrap;
150 <<
"File version " << cp.fileVersion.strVersion();
157 DMType dm_name =
"MOFEM";
161 auto add_bits_tets = [&]() {
164 ->addToDatabaseBitRefLevelByType(MBTET,
BitRefLevel().set(),
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",
"");
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);
185 return std::string();
187 std::string str_number = first_number_string(std::string(
mesh_file_name));
189 if (!str_number.empty())
190 cp.startStep = atoi(str_number.c_str()) + 1;
200 ->getAllEntitiesNotInDatabase(ents);
202 CHKERR m_field.
get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
204 for (
auto m : meshsets)
208 string cutting_surface_name =
210 boost::lexical_cast<std::string>(cp.startStep) +
".vtk";
214 cp.crackAccelerationFactor, cutting_surface_name,
QUIET,
225 string cutting_surface_name =
227 boost::lexical_cast<std::string>(cp.startStep) +
".vtk";
231 CHKERR cp.getInterface<
CPMeshCut>()->copySurface(cutting_surface_name);
233 CHKERR moab.get_entities_by_type(0, MBTET, vol,
false);
244 CHKERR moab.get_entities_by_type(0, MBTET, vol,
false);
252 if (!cp.doCutMesh && !cp.propagateCrack && cp.doElasticWithoutCrack) {
253 BitRefLevel bit1 = cp.mapBitLevel[
"spatial_domain"];
268 std::vector<int> surface_ids;
270 std::vector<std::string> fe_surf_proj_list;
271 fe_surf_proj_list.push_back(
"SURFACE");
273 CHKERR cp.createProblemDataStructures(surface_ids,
QUIET,
278 MOFEM_LOG(
"CP", Sev::verbose) <<
"Testing tangent variant";
281 CHKERR cp.setMaterialPositionFromCoords();
282 CHKERR cp.setSpatialPositionFromCoords();
283 if (cp.addSingularity == PETSC_TRUE) {
284 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
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,
301 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
302 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
304 CHKERR cp.assembleCouplingForcesDM(PETSC_NULL,
VERBOSE,
false);
307 CHKERR cp.setMaterialPositionFromCoords();
308 CHKERR cp.setSpatialPositionFromCoords();
309 if (cp.addSingularity == PETSC_TRUE)
310 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
312 switch (choice_test) {
314 CHKERR cp.testJacobians(cp.mapBitLevel[
"material_domain"],
318 CHKERR cp.testJacobians(cp.mapBitLevel[
"material_domain"],
322 CHKERR cp.testJacobians(cp.mapBitLevel[
"material_domain"],
329 }
else if (cp.propagateCrack == PETSC_TRUE) {
331 MOFEM_LOG(
"CP", Sev::verbose) <<
"Crack propagation variant is used";
334 CHKERR cp.setMaterialPositionFromCoords();
335 CHKERR cp.setSpatialPositionFromCoords();
336 if (cp.addSingularity == PETSC_TRUE)
337 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
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,
351 CHKERR cp.setSpatialPositionFromCoords();
352 if (cp.addSingularity == PETSC_TRUE) {
353 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
356 cp.getLoadScale() = cp.getArcCtx()->getFieldData();
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);
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);
379 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
380 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
382 CHKERR cp.assembleCouplingForcesDM(PETSC_NULL,
VERBOSE,
false);
384 CHKERR cp.calculateElasticEnergy(dm_elastic);
386 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
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,
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);
404 MOFEM_LOG(
"CP", Sev::verbose) <<
"Do not propagate crack";
407 CHKERR cp.setMaterialPositionFromCoords();
408 CHKERR cp.setSpatialPositionFromCoords();
409 if (cp.addSingularity == PETSC_TRUE)
410 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
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,
427 CHKERR cp.createBcDM(dm_bc,
"BC_PROBLEM",
428 cp.mapBitLevel[
"spatial_domain"],
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);
439 MOFEM_LOG(
"CP", Sev::noisy) <<
"Solve ELASTIC problem";
441 dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
445 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
447 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
451 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
453 CHKERR cp.calculateElasticEnergy(dm_elastic);
456 CHKERR cp.postProcessDM(dm_elastic, 0,
"ELASTIC",
true);
458 CHKERR cp.savePositionsOnCrackFrontDM(dm_elastic, PETSC_NULL, 2,
false);
459 if (!cp.doElasticWithoutCrack)
460 CHKERR cp.writeCrackFont(cp.mapBitLevel[
"spatial_domain"], 0);
463 CHKERR cp.tetsingReleaseEnergyCalculation();