v0.14.0
Public Member Functions | Public Attributes | List of all members
FractureMechanics::CPSolvers Struct Reference

#include <users_modules/fracture_mechanics/src/CPSolvers.hpp>

Inheritance diagram for FractureMechanics::CPSolvers:
[legend]
Collaboration diagram for FractureMechanics::CPSolvers:
[legend]

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
 Getting interface of core database. More...
 
 CPSolvers (CrackPropagation &cp)
 
virtual ~CPSolvers ()
 
MoFEMErrorCode getOptions ()
 
MoFEMErrorCode solveElastic (DM dm_elastic, Vec init, const double load_factor)
 Solve elastic problem. More...
 
MoFEMErrorCode calculateGc (DM dm_material_forces, DM dm_surface_projection, DM dm_crack_srf_area, const std::vector< int > &surface_ids)
 Solve for energy realease rate. More...
 
MoFEMErrorCode solvePropagation (DM dm_crack_propagation, DM dm_elastic, DM dm_material, DM dm_material_forces, DM dm_surface_projection, DM dm_crack_srf_area, const std::vector< int > &surface_ids, int cut_mesh_it=0, const bool set_cut_surface=false)
 Solve crack propagation problem. More...
 
MoFEMErrorCode solveCutMesh (SmartPetscObj< DM > &dm_crack_propagation, SmartPetscObj< DM > &dm_elastic, SmartPetscObj< DM > &dm_eigen_elastic, SmartPetscObj< DM > &dm_material, SmartPetscObj< DM > &dm_material_forces, SmartPetscObj< DM > &dm_surface_projection, SmartPetscObj< DM > &dm_crack_srf_area, std::vector< int > &surface_ids, const bool debug=false)
 Solve cutting mesh problem. More...
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Public Attributes

CrackPropagation & cP
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

Crack propagation solvers

Definition at line 28 of file CPSolvers.hpp.

Constructor & Destructor Documentation

◆ CPSolvers()

FractureMechanics::CPSolvers::CPSolvers ( CrackPropagation &  cp)

Definition at line 56 of file CPSolvers.cpp.

56  : cP(cp) {
57 
58  if (!LogManager::checkIfChannelExist("CPSolverWorld")) {
59  auto core_log = logging::core::get();
60 
61  core_log->add_sink(
62  LogManager::createSink(LogManager::getStrmWorld(), "CPSolverWorld"));
63  core_log->add_sink(
64  LogManager::createSink(LogManager::getStrmSync(), "CPSolverSync"));
65  core_log->add_sink(
66  LogManager::createSink(LogManager::getStrmSelf(), "CPSolverSelf"));
67 
68  LogManager::setLog("CPSolverWorld");
69  LogManager::setLog("CPSolverSync");
70  LogManager::setLog("CPSolverSelf");
71 
72  MOFEM_LOG_TAG("CPSolverWorld", "CPSolve");
73  MOFEM_LOG_TAG("CPSolverSync", "CPSolve");
74  MOFEM_LOG_TAG("CPSolverSelf", "CPSolve");
75  }
76 
77  MOFEM_LOG("CPSolverWorld", Sev::noisy) << "CPSolve created";
78 }

◆ ~CPSolvers()

virtual FractureMechanics::CPSolvers::~CPSolvers ( )
inlinevirtual

Definition at line 41 of file CPSolvers.hpp.

41 {}

Member Function Documentation

◆ calculateGc()

MoFEMErrorCode FractureMechanics::CPSolvers::calculateGc ( DM  dm_material_forces,
DM  dm_surface_projection,
DM  dm_crack_srf_area,
const std::vector< int > &  surface_ids 
)

Solve for energy realease rate.

Parameters
dm_material_forces
dm_surface_projection
dm_crack_srf_area
surface_ids
Returns
MoFEMErrorCode

Definition at line 186 of file CPSolvers.cpp.

189  {
190  MoFEM::Interface &m_field = cP.mField;
191  MoFEMFunctionBegin;
192 
193  // Create vector of material forces and set operators
194  Vec q_material, f_material;
195  CHKERR DMCreateGlobalVector(dm_material_forces, &q_material);
196  CHKERR VecDuplicate(q_material, &f_material);
197  Vec f_material_proj, f_griffith, f_griffith_proj;
198  CHKERR VecDuplicate(f_material, &f_material_proj);
199  CHKERR VecDuplicate(f_material, &f_griffith);
200  CHKERR VecDuplicate(f_griffith, &f_griffith_proj);
201  Vec f_lambda;
202  {
203  const MoFEM::Problem *prb_ptr;
204  CHKERR DMMoFEMGetProblemPtr(dm_crack_srf_area, &prb_ptr);
205  // Need to use directly MoFEM interface to create vector, matrix for non
206  // square matrix
207  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
208  prb_ptr->getName(), ROW, &f_lambda);
209  }
210 
211  // Calculate projection matrices
212  CHKERR cP.calculateSurfaceProjectionMatrix(
213  dm_surface_projection, dm_material_forces, surface_ids, VERBOSE, false);
214  CHKERR cP.calculateFrontProjectionMatrix(dm_crack_srf_area,
215  dm_material_forces, VERBOSE, false);
216  // Calculate griffith force vector (material resistance)
217  CHKERR cP.calculateGriffithForce(dm_material_forces, cP.gC ? cP.gC : 1,
218  f_griffith, VERBOSE, false);
219  CHKERR cP.projectGriffithForce(dm_material_forces, f_griffith,
220  f_griffith_proj, VERBOSE, false);
221  // Calculate material forces
222  CHKERR DMoFEMMeshToLocalVector(dm_material_forces, q_material, INSERT_VALUES,
223  SCATTER_FORWARD);
224 
225  CHKERR cP.calculateMaterialForcesDM(dm_material_forces, q_material,
226  f_material, VERBOSE, false);
227 
228  // project material forces
229  CHKERR cP.projectMaterialForcesDM(dm_material_forces, f_material,
230  f_material_proj, VERBOSE, false);
231  // calculate griffith energy
232  CHKERR cP.calculateReleaseEnergy(dm_crack_srf_area, f_material_proj,
233  f_griffith_proj, f_lambda, cP.gC, VERBOSE,
234  false);
235 
236  CHKERR VecDestroy(&q_material);
237  CHKERR VecDestroy(&f_material);
238  CHKERR VecDestroy(&f_material_proj);
239  CHKERR VecDestroy(&f_griffith);
240  CHKERR VecDestroy(&f_griffith_proj);
241  CHKERR VecDestroy(&f_lambda);
242 
243  MoFEMFunctionReturn(0);
244 }

◆ getOptions()

MoFEMErrorCode FractureMechanics::CPSolvers::getOptions ( )

Definition at line 80 of file CPSolvers.cpp.

80  {
81  MoFEMFunctionBegin;
82 
83  // CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Solvers options", "none");
84  // ierr = PetscOptionsEnd();
85  // CHKERRG(ierr);
86 
87  MoFEMFunctionReturn(0);
88 }

◆ query_interface()

MoFEMErrorCode FractureMechanics::CPSolvers::query_interface ( boost::typeindex::type_index  type_index,
MoFEM::UnknownInterface **  iface 
) const
virtual

Getting interface of core database.

Parameters
uuidunique ID of interface of CPSolvers solely
ifacereturned pointer to interface
Returns
error code

Implements MoFEM::UnknownInterface.

Definition at line 50 of file CPSolvers.cpp.

51  {
52  *iface = const_cast<CPSolvers *>(this);
53  return 0;
54 }

◆ solveCutMesh()

MoFEMErrorCode FractureMechanics::CPSolvers::solveCutMesh ( SmartPetscObj< DM > &  dm_crack_propagation,
SmartPetscObj< DM > &  dm_elastic,
SmartPetscObj< DM > &  dm_eigen_elastic,
SmartPetscObj< DM > &  dm_material,
SmartPetscObj< DM > &  dm_material_forces,
SmartPetscObj< DM > &  dm_surface_projection,
SmartPetscObj< DM > &  dm_crack_srf_area,
std::vector< int > &  surface_ids,
const bool  debug = false 
)

Solve cutting mesh problem.

Parameters
dm_crack_propagationcomposition of elastic DM and material DM
dm_elasticelastic DM
dm_eigen_elasticeigen elastic DM
dm_materialmaterial DM
dm_material_forcesused to calculate material forces, sub-dm of dm_crack_propagation
dm_surface_projectiondm to project material forces on surfaces, sub-dm of dm_crack_propagation
dm_crack_srf_areadm to project material forces on crack surfaces, sub-dm of dm_crack_propagation
surface_idsIDs of surface meshsets
debugflag for debugging
Returns
MoFEMErrorCode

Definition at line 850 of file CPSolvers.cpp.

858  {
859  SolveFunctionBegin;
860 
861  auto calculate_load_factor = [this]() {
862  double a = fabs(cP.maxG1) / pow(cP.getArcCtx()->getFieldData(), 2);
863  double load_factor =
864  copysign(sqrt(cP.gC / a), cP.getArcCtx()->getFieldData());
865  return load_factor;
866  };
867  cP.getArcCtx()->getFieldData() = calculate_load_factor();
868  MOFEM_LOG_C("CPSolverWorld", Sev::inform, "Calculated load factor %g",
869  cP.getArcCtx()->getFieldData());
870  MOFEM_LOG_C("CPSolverSync", Sev::noisy, "Calculated load factor %g",
871  cP.getArcCtx()->getFieldData());
872  MOFEM_LOG_SYNCHRONISE(cP.mField.get_comm());
873 
874  auto save_position_on_mesh_tags_from_coords =
875  [this](const std::string tag_name) {
876  Tag th;
877  MoFEMFunctionBegin;
878  rval = cP.mField.get_moab().tag_get_handle(tag_name.c_str(), th);
879  if (rval != MB_SUCCESS) {
880  double def[] = {0, 0, 0};
881  CHKERR cP.mField.get_moab().tag_get_handle(
882  tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
883  MB_TAG_CREAT | MB_TAG_SPARSE, def);
884  }
885  Range verts;
886  CHKERR cP.mField.getInterface<BitRefManager>()
887  ->getEntitiesByTypeAndRefLevel(cP.mapBitLevel["spatial_domain"],
888  BitRefLevel().set(), MBVERTEX,
889  verts);
890  std::vector<double> coords(verts.size() * 3);
891  CHKERR cP.mField.get_moab().get_coords(verts, &*coords.begin());
892  CHKERR cP.mField.get_moab().tag_set_data(th, verts, &*coords.begin());
893  MoFEMFunctionReturn(0);
894  };
895 
896  int nb_cuts = 0;
897  for (;;) {
898 
899  CHKERR save_position_on_mesh_tags_from_coords("MESH_NODE_POSITIONS");
900  CHKERR solvePropagation(dm_crack_propagation, dm_elastic, dm_material,
901  dm_material_forces, dm_surface_projection,
902  dm_crack_srf_area, surface_ids, nb_cuts, true);
903 
904  // Destroy DMs
905  auto reset_dms = [&]() {
906  dm_crack_propagation.reset();
907  dm_elastic.reset();
908  dm_eigen_elastic.reset();
909  dm_material.reset();
910  dm_material_forces.reset();
911  dm_surface_projection.reset();
912  dm_crack_srf_area.reset();
913  };
914 
915  auto reset_fes = [&]() {
916  std::vector<std::pair<boost::weak_ptr<FEMethod>, std::string>> vec_fe;
917  std::vector<
918  std::pair<boost::weak_ptr<NonlinearElasticElement>, std::string>>
919  v_elastic_ele_str;
920 
921  auto push_fes = [&]() {
922  v_elastic_ele_str.emplace_back(
923  std::make_pair(cP.elasticFe, "elasticFe"));
924  v_elastic_ele_str.emplace_back(
925  std::make_pair(cP.materialFe, "materialFe"));
926  vec_fe.emplace_back(std::make_pair(cP.feLhs, "feLhs"));
927  vec_fe.emplace_back(std::make_pair(cP.feRhs, "feRhs"));
928  vec_fe.emplace_back(std::make_pair(cP.feMaterialRhs, "feMaterialRhs"));
929  vec_fe.emplace_back(std::make_pair(cP.feMaterialLhs, "feMaterialLhs"));
930  vec_fe.emplace_back(std::make_pair(cP.feEnergy, "feEnergy"));
931  vec_fe.emplace_back(
932  std::make_pair(cP.feCouplingElasticLhs, "feCouplingElasticLhs"));
933  vec_fe.emplace_back(
934  std::make_pair(cP.feCouplingMaterialLhs, "feCouplingMaterialLhs"));
935  vec_fe.emplace_back(std::make_pair(cP.feSmootherRhs, "feSmootherRhs"));
936  vec_fe.emplace_back(std::make_pair(cP.feSmootherLhs, "feSmootherLhs"));
937  vec_fe.emplace_back(
938  std::make_pair(cP.feGriffithForceRhs, "feGriffithForceRhs"));
939  vec_fe.emplace_back(std::make_pair(cP.feGriffithConstrainsDelta,
940  "feGriffithConstrainsDelta"));
941  vec_fe.emplace_back(std::make_pair(cP.feGriffithConstrainsRhs,
942  "feGriffithConstrainsRhs"));
943  vec_fe.emplace_back(std::make_pair(cP.feGriffithConstrainsLhs,
944  "feGriffithConstrainsLhs"));
945  vec_fe.emplace_back(
946  std::make_pair(cP.feSpringLhsPtr, "feSpringLhsPtr"));
947  vec_fe.emplace_back(
948  std::make_pair(cP.feSpringRhsPtr, "feSpringRhsPtr"));
949  vec_fe.emplace_back(
950  std::make_pair(cP.feRhsSimpleContact, "feRhsSimpleContact"));
951  vec_fe.emplace_back(
952  std::make_pair(cP.feLhsSimpleContact, "feLhsSimpleContact"));
953  vec_fe.emplace_back(std::make_pair(cP.feMaterialAnaliticalTraction,
954  "feMaterialAnaliticalTraction"));
955 
956  vec_fe.emplace_back(
957  std::make_pair(cP.bothSidesConstrains, "bothSidesConstrains"));
958  vec_fe.emplace_back(
959  std::make_pair(cP.closeCrackConstrains, "closeCrackConstrains"));
960  vec_fe.emplace_back(std::make_pair(cP.bothSidesContactConstrains,
961  "bothSidesContactConstrains"));
962  vec_fe.emplace_back(
963  std::make_pair(cP.fixMaterialEnts, "fixMaterialEnts"));
964  vec_fe.emplace_back(
965  std::make_pair(cP.feSpringLhsPtr, "feSpringLhsPtr"));
966  vec_fe.emplace_back(
967  std::make_pair(cP.feSpringRhsPtr, "feSpringRhsPtr"));
968  vec_fe.emplace_back(
969  std::make_pair(cP.assembleFlambda, "assembleFlambda"));
970  vec_fe.emplace_back(std::make_pair(cP.zeroFlambda, "zeroFlambda"));
971  };
972 
973  auto reset = [&]() {
974  cP.smootherFe.reset();
975  cP.feSmootherRhs.reset();
976  cP.feSmootherLhs.reset();
977  cP.volumeLengthDouble.reset();
978  cP.volumeLengthAdouble.reset();
979  cP.griffithForceElement.reset();
980 
981  cP.tangentConstrains.reset();
982  cP.skinOrientation.reset();
983  cP.crackOrientation.reset();
984  cP.contactOrientation.reset();
985  for (auto &m : cP.surfaceConstrain)
986  m.second.reset();
987  for (auto &m : cP.edgeConstrains)
988  m.second.reset();
989 
990  cP.projSurfaceCtx.reset();
991  cP.projFrontCtx.reset();
992 
993  cP.elasticFe.reset();
994  cP.materialFe.reset();
995  cP.feLhs.reset();
996  cP.feRhs.reset();
997  cP.feMaterialRhs.reset();
998  cP.feMaterialLhs.reset();
999  cP.feEnergy.reset();
1000  cP.feCouplingElasticLhs.reset();
1001  cP.feCouplingMaterialLhs.reset();
1002  cP.feSmootherRhs.reset();
1003  cP.feSmootherLhs.reset();
1004  cP.feGriffithForceRhs.reset();
1005  cP.feGriffithConstrainsDelta.reset();
1006  cP.feGriffithConstrainsRhs.reset();
1007  cP.feGriffithConstrainsLhs.reset();
1008  cP.feSpringLhsPtr.reset();
1009  cP.feSpringRhsPtr.reset();
1010  cP.feRhsSimpleContact.reset();
1011  cP.feLhsSimpleContact.reset();
1012  cP.feMaterialAnaliticalTraction.reset();
1013  cP.bothSidesConstrains.reset();
1014  cP.bothSidesContactConstrains.reset();
1015  cP.fixMaterialEnts.reset();
1016  cP.assembleFlambda.reset();
1017  cP.zeroFlambda.reset();
1018  cP.closeCrackConstrains.reset();
1019  };
1020 
1021  auto check = [&]() {
1022  for (auto v : vec_fe)
1023  if (auto a = v.first.lock()) {
1024  MOFEM_LOG("CPSolverWorld", Sev::error)
1025  << "fe " << v.second << " not destroyed " << a.use_count();
1026  }
1027 
1028  for (auto v : v_elastic_ele_str)
1029  if (auto a = v.first.lock()) {
1030  MOFEM_LOG("CPSolverWorld", Sev::error)
1031  << "fe elastic " << v.second << " not destroyed "
1032  << a.use_count();
1033  }
1034 
1035  if (dm_crack_propagation.use_count())
1036  MOFEM_LOG("CPSolverWorld", Sev::error)
1037  << "dm_crack_propagation not destroyed "
1038  << dm_crack_propagation.use_count();
1039  if (dm_elastic.use_count())
1040  MOFEM_LOG("CPSolverWorld", Sev::error)
1041  << "dm_elastic not destroyed" << dm_elastic.use_count();
1042  if (dm_material_forces.use_count())
1043  MOFEM_LOG("CPSolverWorld", Sev::error)
1044  << "dm_material_forces not destroyed "
1045  << dm_material_forces.use_count();
1046  if (dm_surface_projection.use_count())
1047  MOFEM_LOG("CPSolverWorld", Sev::error)
1048  << "dm_surface_projection not destroyed "
1049  << dm_surface_projection.use_count();
1050  if (dm_crack_srf_area.use_count())
1051  MOFEM_LOG("CPSolverWorld", Sev::error)
1052  << "dm_crack_srf_area not destroyed "
1053  << dm_crack_srf_area.use_count();
1054 
1055  if (cP.tangentConstrains.use_count())
1056  MOFEM_LOG("CPSolverWorld", Sev::error)
1057  << "tangentConstrains not destroyed";
1058  if (cP.skinOrientation.use_count())
1059  MOFEM_LOG("CPSolverWorld", Sev::error)
1060  << "skinOrientation not destroyed";
1061  if (cP.crackOrientation.use_count())
1062  MOFEM_LOG("CPSolverWorld", Sev::error)
1063  << "crackOrientation not destroyed";
1064  if (cP.contactOrientation.use_count())
1065  MOFEM_LOG("CPSolverWorld", Sev::error)
1066  << "contactOrientation not destroyed";
1067  for (auto &m : cP.surfaceConstrain)
1068  if (m.second.use_count())
1069  MOFEM_LOG("CPSolverWorld", Sev::error)
1070  << "surfaceConstrain not destroyed : " << m.first;
1071  for (auto &m : cP.edgeConstrains)
1072  if (m.second.use_count())
1073  MOFEM_LOG("CPSolverWorld", Sev::error)
1074  << "edgeConstrains not destroyed : " << m.first;
1075  cP.surfaceConstrain.clear();
1076  cP.edgeConstrains.clear();
1077 
1078  if (cP.projSurfaceCtx.use_count())
1079  MOFEM_LOG("CPSolverWorld", Sev::error)
1080  << "projSurfaceCtx not destroyed";
1081  if (cP.projFrontCtx.use_count())
1082  MOFEM_LOG("CPSolverWorld", Sev::error)
1083  << "projFrontCtx not destroyed";
1084  };
1085 
1086  push_fes();
1087  reset();
1088  check();
1089  };
1090 
1091  reset_dms();
1092  reset_fes();
1093 
1094  if (cP.mwlsApprox) {
1095  MOFEM_LOG("MWLSWorld", Sev::inform)
1096  << "Resest MWLS approximation coefficients. Coefficients will be "
1097  "recalulated for current material positions.";
1098  cP.mwlsApprox->invABMap.clear();
1099  cP.mwlsApprox->influenceNodesMap.clear();
1100  cP.mwlsApprox->dmNodesMap.clear();
1101  }
1102 
1103  ++nb_cuts;
1104  if (nb_cuts <= cP.getNbCutSteps()) {
1105 
1106  // clear database
1107  CHKERR cP.mField.clear_database(QUIET);
1108  CHKERR cP.getInterface<CPMeshCut>()->clearData();
1109 
1110  auto tag_gid = cP.mField.get_moab().globalId_tag();
1111  std::vector<Tag> tag_handles;
1112  CHKERR cP.mField.get_moab().tag_get_tags(tag_handles);
1113  for (auto th : tag_handles)
1114  if (th != tag_gid)
1115  CHKERR cP.mField.get_moab().tag_delete(th);
1116  CHKERR cP.mField.get_moab().delete_mesh();
1117 
1118  // reaload mesh
1119  CHKERR PetscBarrier(PETSC_NULL);
1120  const std::string file_name =
1121  "restart_" + boost::lexical_cast<std::string>(cP.startStep - 1) +
1122  ".h5m";
1123 
1124  auto add_bits_tets = [&]() {
1125  MoFEMFunctionBegin;
1126  CHKERR cP.mField.getInterface<BitRefManager>()
1127  ->addToDatabaseBitRefLevelByType(MBTET, BitRefLevel().set(),
1128  BitRefLevel().set());
1129  MoFEMFunctionReturn(0);
1130  };
1131 
1132  std::vector<int> surface_ids;
1133  surface_ids.push_back(cP.getInterface<CPMeshCut>()->getSkinOfTheBodyId());
1134  std::vector<std::string> fe_surf_proj_list;
1135  fe_surf_proj_list.push_back("SURFACE");
1136 
1137  ParallelComm *pcomm =
1138  ParallelComm::get_pcomm(&cP.mField.get_moab(), MYPCOMM_INDEX);
1139  if (pcomm == NULL)
1140  pcomm = new ParallelComm(&cP.mField.get_moab(),
1141  cP.moabCommWorld->get_comm());
1142 
1143  if (cP.mField.get_comm_size() == 1 ||
1144  !CrackPropagation::parallelReadAndBroadcast) {
1145  const char *option = "";
1146  CHKERR cP.mField.get_moab().load_file(file_name.c_str(), 0, option);
1147  } else {
1148  if (cP.mField.get_comm_rank() == 0) {
1149  // Read mesh file
1150  moab::Core mb_instance_read;
1151  moab::Interface &moab_read = mb_instance_read;
1152  ParallelComm *pcomm_read =
1153  ParallelComm::get_pcomm(&moab_read, MYPCOMM_INDEX);
1154  if (pcomm_read == NULL)
1155  pcomm_read =
1156  new ParallelComm(&moab_read, cP.moabCommWorld->get_comm());
1157 
1158  const char *option = "";
1159  CHKERR moab_read.load_file(file_name.c_str(), 0, option);
1160  Range ents;
1161  CHKERR moab_read.get_entities_by_handle(0, ents, false);
1162  CHKERR moab_read.get_entities_by_type(0, MBENTITYSET, ents, false);
1163  CHKERR FractureMechanics::broadcast_entities(
1164  cP.mField.get_moab(), moab_read, 0, ents, false, true);
1165  } else {
1166  Range ents;
1167  CHKERR FractureMechanics::broadcast_entities(
1168  cP.mField.get_moab(), cP.mField.get_moab(), 0, ents, false, true);
1169  }
1170  }
1171 
1172  CHKERR MoFEM::Interface::setFileVersion(cP.mField.get_moab(),
1173  cP.fileVersion);
1174 
1175  CHKERR cP.mField.set_moab_interface(cP.mField.get_moab(), QUIET);
1176  CHKERR cP.getInterface<CPMeshCut>()->getInterfacesPtr();
1177  CHKERR add_bits_tets();
1178 
1179  Range ents;
1180  CHKERR cP.mField.getInterface<BitRefManager>()
1181  ->getAllEntitiesNotInDatabase(ents);
1182  Range meshsets;
1183  CHKERR cP.mField.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
1184  false);
1185  for (auto m : meshsets)
1186  CHKERR cP.mField.get_moab().remove_entities(m, ents);
1187  CHKERR cP.mField.get_moab().delete_entities(ents);
1188 
1189  string cutting_surface_name =
1190  "cutting_surface_" + boost::lexical_cast<std::string>(cP.startStep) +
1191  ".vtk";
1192 
1193  if (cP.mField.get_comm_rank() == 0)
1194  CHKERR cP.getInterface<CPMeshCut>()->rebuildCrackSurface(
1195  cP.crackAccelerationFactor, cutting_surface_name, QUIET, debug);
1196  else
1197  CHKERR cP.getInterface<CPMeshCut>()->rebuildCrackSurface(
1198  cP.crackAccelerationFactor, "", QUIET, debug);
1199 
1200  CHKERR cP.getInterface<CPMeshCut>()->refineMesh(nullptr, false, QUIET,
1201  debug);
1202  CHKERR cP.getInterface<CPMeshCut>()->cutRefineAndSplit(QUIET, debug);
1203 
1204  // Create crack propagation fields and elements
1205  CHKERR cP.mField.build_field("LAMBDA_ARC_LENGTH");
1206  CHKERR cP.mField.build_finite_elements("ARC_LENGTH");
1207 
1208  CHKERR cP.createProblemDataStructures(surface_ids, QUIET, debug);
1209 
1210  // set inital coordinates
1211  CHKERR cP.setMaterialPositionFromCoords();
1212  CHKERR cP.setSpatialPositionFromCoords();
1213  if (cP.addSingularity == PETSC_TRUE)
1214  CHKERR cP.setSingularDofs("SPATIAL_POSITION");
1215 
1216  CHKERR cP.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
1217  dm_crack_propagation, dm_material_forces,
1218  dm_surface_projection, dm_crack_srf_area, surface_ids,
1219  fe_surf_proj_list);
1220 
1221  auto set_up_arc_length = [&]() {
1222  MoFEMFunctionBegin;
1223  const MoFEM::Problem *problem_ptr;
1224  CHKERR DMMoFEMGetProblemPtr(dm_elastic, &problem_ptr);
1225  cP.arcCtx = boost::shared_ptr<ArcLengthCtx>(new ArcLengthCtx(
1226  cP.mField, problem_ptr->getName(), "LAMBDA_ARC_LENGTH"));
1227  CHKERR cP.arcCtx->setAlphaBeta(0, 1);
1228  auto arc_snes_ctx =
1229  boost::make_shared<CrackPropagation::ArcLengthSnesCtx>(
1230  cP.mField, problem_ptr->getName(), cP.arcCtx);
1231  CHKERR DMMoFEMSetSnesCtx(dm_elastic, arc_snes_ctx);
1232  MoFEMFunctionReturn(0);
1233  };
1234 
1235  CHKERR set_up_arc_length();
1236 
1237  CHKERR cP.setSpatialPositionFromCoords();
1238  if (cP.addSingularity == PETSC_TRUE)
1239  CHKERR cP.setSingularDofs("SPATIAL_POSITION");
1240 
1241  const auto load_factor = std::abs(cP.getArcCtx()->getFieldData());
1242  MOFEM_LOG("CPSolverSync", Sev::noisy) << "Lambda factor " << load_factor;
1243  MOFEM_LOG_SYNCHRONISE(cP.mField.get_comm());
1244  if (cP.solveEigenStressProblem) {
1245  MOFEM_LOG("CPSolverWorld", Sev::verbose)
1246  << "Solve Eigen ELASTIC problem";
1247  CHKERR cP.getInterface<CPSolvers>()->solveElastic(
1248  dm_eigen_elastic, cP.getEigenArcCtx()->x0, 1);
1249  }
1250  cP.getLoadScale() = load_factor;
1251  MOFEM_LOG("CPSolverSync", Sev::noisy)
1252  << "Reset lambda factor " << cP.getLoadScale();
1253  MOFEM_LOG("CPSolverWorld", Sev::noisy)
1254  << "Solve ELASTIC problem and calculate g";
1255  CHKERR cP.getInterface<CPSolvers>()->solveElastic(
1256  dm_elastic, cP.getArcCtx()->x0, load_factor);
1257 
1258  // set finite element instances and user data operators on instances
1259  CHKERR cP.assembleElasticDM(VERBOSE, false);
1260  CHKERR cP.assembleMaterialForcesDM(PETSC_NULL);
1261  CHKERR cP.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE,
1262  false);
1263  CHKERR cP.assembleCouplingForcesDM(PETSC_NULL, VERBOSE, false);
1264 
1265  CHKERR cP.calculateElasticEnergy(dm_elastic);
1266  CHKERR cP.getInterface<CPSolvers>()->calculateGc(
1267  dm_material_forces, dm_surface_projection, dm_crack_srf_area,
1268  surface_ids);
1269  cP.getArcCtx()->getFieldData() = calculate_load_factor();
1270  MOFEM_LOG_C("CPSolverWorld", Sev::inform, "Calculated load factor %g",
1271  cP.getArcCtx()->getFieldData());
1272 
1273  } else
1274  break;
1275  }
1276 
1277  MoFEMFunctionReturn(0);
1278 }

◆ solveElastic()

MoFEMErrorCode FractureMechanics::CPSolvers::solveElastic ( DM  dm_elastic,
Vec  init,
const double  load_factor 
)

Solve elastic problem.

Parameters
dm_elastic
Returns
MoFEMErrorCode

Clear finite elements added SNES solver

Definition at line 90 of file CPSolvers.cpp.

91  {
92  MoFEM::Interface &m_field = cP.mField;
93  MoFEMFunctionBegin;
94 
95  // Create matrix
96  auto m_elastic = smartCreateDMMatrix(dm_elastic);
97  // Create vector
98  auto q_elastic = smartCreateDMVector(dm_elastic);
99  auto f_elastic = smartVectorDuplicate(q_elastic);
100 
101  const MoFEM::Problem *problem_ptr;
102  CHKERR DMMoFEMGetProblemPtr(dm_elastic, &problem_ptr);
103 
104  auto solve_elastic_problem = [&](auto init_x) {
105  MoFEMFunctionBegin;
106 
107  /// Clear finite elements added SNES solver
108  smartGetDMSnesCtx(dm_elastic)->clearLoops();
109 
110  boost::shared_ptr<ArcLengthCtx> arc_ctx;
111  if (problem_ptr->getName() == "ELASTIC")
112  arc_ctx = cP.getArcCtx();
113  if (problem_ptr->getName() == "EIGEN_ELASTIC")
114  arc_ctx = cP.getEigenArcCtx();
115 
116  // Set operators for elasticity analysis
117  boost::shared_ptr<SimpleArcLengthControl> arc_method(
118  new SimpleArcLengthControl(arc_ctx));
119 
120  // Add operators to DM SNES
121  CHKERR cP.addElasticFEInstancesToSnes(dm_elastic, PETSC_NULL, PETSC_NULL,
122  PETSC_NULL, arc_method, arc_ctx,
123  VERBOSE, false);
124 
125  CHKERR arc_ctx->setAlphaBeta(0, 1);
126  CHKERR arc_ctx->setS(0);
127  arc_ctx->getFieldData() = load_factor;
128 
129  CHKERR cP.getArcLengthDof();
130 
131  if (m_field.check_problem("BC_PROBLEM")) {
132  // solve problem to find dofs for analytical displacements
133  cP.getAnalyticalDirichletBc()->snes_B = m_elastic;
134  cP.getAnalyticalDirichletBc()->snes_x = q_elastic;
135  cP.getAnalyticalDirichletBc()->snes_f = f_elastic;
136  // solve for Dirichlet bc dofs
137  AnalyticalDirichletBC analytical_bc(m_field);
138  CHKERR analytical_bc.setUpProblem(m_field, "BC_PROBLEM");
139  boost::shared_ptr<AnalyticalDisp> testing_function =
140  boost::shared_ptr<AnalyticalDisp>(new AnalyticalDisp());
141  // EntityHandle fe_meshset =
142  // m_field.get_finite_element_meshset("BC_FE");
143  CHKERR analytical_bc.setApproxOps(m_field, "SPATIAL_POSITION",
144  testing_function, 0,
145  "MESH_NODE_POSITIONS");
146 
147  analytical_bc.approxField.getLoopFeApprox().addToRule = 1;
148  CHKERR analytical_bc.solveProblem(m_field, "BC_PROBLEM", "BC_FE",
149  *cP.getAnalyticalDirichletBc());
150 
151  CHKERR analytical_bc.destroyProblem();
152  }
153 
154  // set vector values from field data
155  CHKERR DMoFEMMeshToLocalVector(dm_elastic, q_elastic, INSERT_VALUES,
156  SCATTER_FORWARD);
157 
158  auto snes = createSNES(m_field.get_comm());
159  Mat shell_m;
160 
161  // for(int ll = 0;ll!=cP.getNbLoadSteps();ll++) {
162  CHKERR VecCopy(q_elastic, init_x);
163  CHKERR cP.solveElasticDM(dm_elastic, snes, m_elastic, q_elastic, f_elastic,
164  true, &shell_m);
165  CHKERR DMoFEMMeshToLocalVector(dm_elastic, q_elastic, INSERT_VALUES,
166  SCATTER_REVERSE);
167 
168  CHKERR MatDestroy(&shell_m);
169 
170  MoFEMFunctionReturn(0);
171  };
172 
173  // If some flag is set, we solve eigen problem, i.e. calculate eigen
174  // displacements
175  if (problem_ptr->getName() == "EIGEN_ELASTIC")
176  CHKERR cP.assembleElasticDM(cP.mwlsEigenStressTagName);
177  if (problem_ptr->getName() == "ELASTIC")
178  CHKERR cP.assembleElasticDM(cP.mwlsStressTagName);
179 
180  // + and ops which will close crack
181  CHKERR solve_elastic_problem(init);
182 
183  MoFEMFunctionReturn(0);
184 }

◆ solvePropagation()

MoFEMErrorCode FractureMechanics::CPSolvers::solvePropagation ( DM  dm_crack_propagation,
DM  dm_elastic,
DM  dm_material,
DM  dm_material_forces,
DM  dm_surface_projection,
DM  dm_crack_srf_area,
const std::vector< int > &  surface_ids,
int  cut_mesh_it = 0,
const bool  set_cut_surface = false 
)

Solve crack propagation problem.

Parameters
dm_crack_propagationcomposition of elastic DM and material DM
dm_elasticelastic DM
dm_materialmaterial DM
dm_material_forcesused to calculate material forces, sub-dm of dm_crack_propagation
dm_surface_projectiondm to project material forces on surfaces, sub-dm of dm_crack_propagation
dm_crack_srf_areadm to project material forces on crack surfaces, sub-dm of dm_crack_propagation
surface_idsIDs of surface meshsets
cut_mesh_itnumber of catting steps
set_cut_surfaceflag to cut surfaces
Returns
MoFEMErrorCode

Definition at line 247 of file CPSolvers.cpp.

251  {
252  MoFEM::Interface &m_field = cP.mField;
253  MoFEMFunctionBegin;
254 
255  const MoFEM::Problem *problem_elastic_ptr;
256  CHKERR DMMoFEMGetProblemPtr(dm_elastic, &problem_elastic_ptr);
257  cP.arcCtx = boost::make_shared<ArcLengthCtx>(
258  cP.mField, problem_elastic_ptr->getName(), "LAMBDA_ARC_LENGTH");
259  auto arc_snes_ctx = boost::make_shared<CrackPropagation::ArcLengthSnesCtx>(
260  cP.mField, problem_elastic_ptr->getName(), cP.arcCtx);
261  CHKERR DMMoFEMSetSnesCtx(dm_elastic, arc_snes_ctx);
262 
263  // set finite element instances and user data operators on instances
264  CHKERR cP.assembleElasticDM(VERBOSE, false);
265  MOFEM_LOG("CPSolverWorld", Sev::noisy)
266  << "Solve ELASTIC problem in solvePropagation";
267  CHKERR solveElastic(dm_elastic, cP.getArcCtx()->x0,
268  cP.getArcCtx()->getFieldData());
269 
270  CHKERR cP.assembleMaterialForcesDM(PETSC_NULL);
271  CHKERR cP.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE, false);
272  CHKERR cP.assembleCouplingForcesDM(PETSC_NULL, VERBOSE, false);
273 
274  struct StepAdaptivity {
275 
276  CrackPropagation &cP;
277 
278  int itsD;
279  double gAmma;
280  double minStep;
281  double maxStep;
282  PetscBool minStepFlg;
283  PetscBool maxStepFlg;
284 
285  StepAdaptivity(CrackPropagation &cp, int its_d, double gamma)
286  : cP(cp), itsD(its_d), gAmma(gamma), minStep(0), maxStep(0) {
287  ierr = getOptions();
288  CHKERRABORT(PETSC_COMM_WORLD, ierr);
289  }
290 
291  MoFEMErrorCode getOptions() {
292  SolveFunctionBegin;
293  ierr =
294  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Fracture options", "none");
295  CHKERRQ(ierr);
296  {
297  CHKERR PetscOptionsInt("-adapt_step_its_d",
298  "number of desired iterations", "", itsD, &itsD,
299  PETSC_NULL);
300  CHKERR PetscOptionsScalar("-adapt_step_min_s", "minimal setp", "",
301  minStep, &minStep, &minStepFlg);
302  CHKERR PetscOptionsScalar("-adapt_step_max_s", "maximal setp", "",
303  maxStep, &maxStep, &maxStepFlg);
304  }
305  ierr = PetscOptionsEnd();
306  CHKERRQ(ierr);
307 
308  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
309  "### Input parameter: -adapt_step_its_d %d", itsD);
310  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
311  "### Input parameter: -adapt_step_min_s %6.4e", minStep);
312  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
313  "### Input parameter: -adapt_step_max_s %6.4e", maxStep);
314 
315  if (minStepFlg == PETSC_TRUE && maxStepFlg == PETSC_TRUE)
316  if (minStep > maxStep) {
317  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
318  "Minimal step size cannot be bigger than maximal step size");
319  }
320  MoFEMFunctionReturn(0);
321  }
322 
323  MoFEMErrorCode operator()(DM dm, SNES snes,
324  boost::shared_ptr<ArcLengthCtx> arc_ctx,
325  double &s) const {
326  SolveFunctionBegin;
327 
328  auto min_max = [this](const double s) {
329  double new_s = s;
330  if (minStepFlg == PETSC_TRUE) {
331  new_s = std::max(new_s, minStep);
332  }
333  if (maxStepFlg == PETSC_TRUE) {
334  new_s = std::min(new_s, maxStep);
335  }
336  return new_s;
337  };
338 
339  SNESConvergedReason reason;
340  CHKERR SNESGetConvergedReason(snes, &reason);
341  if (reason < 0) {
342  s = min_max(0.25 * s);
343  MOFEM_LOG_C("CPSolverWorld", Sev::warning,
344  "*** Diverged ***"
345  "Setting reduced load step arc_s = %3.4g",
346  s);
347  CHKERR DMoFEMMeshToLocalVector(dm, arc_ctx->x0, INSERT_VALUES,
348  SCATTER_REVERSE);
349 
350  } else {
351  int its;
352  CHKERR SNESGetIterationNumber(snes, &its);
353  s = min_max(s * pow((double)itsD / (double)(its + 1), gAmma));
354  }
355  MoFEMFunctionReturn(0);
356  }
357  };
358 
359  struct SmoothingAlphaAdaptivity {
360 
361  CrackPropagation &cP;
362  double gAmma;
363  double dAlpha;
364  double minAlpha;
365  double maxAlpha;
366  double incrementAlpha;
367 
368  SmoothingAlphaAdaptivity(CrackPropagation &cp, const double gamma,
369  const double d_alpha, const double m_alpha)
370  : cP(cp), gAmma(gamma), dAlpha(d_alpha), minAlpha(m_alpha) {
371  ierr = getOptions();
372  CHKERRABORT(PETSC_COMM_WORLD, ierr);
373  }
374 
375  MoFEMErrorCode getOptions() {
376  MoFEMFunctionBegin;
377  ierr =
378  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Fracture options", "none");
379  CHKERRQ(ierr);
380  CHKERR PetscOptionsScalar("-adapt_min_smoother_alpha",
381  "minimal smoother alpha", "", minAlpha,
382  &minAlpha, PETSC_NULL);
383  maxAlpha = 1e3;
384  CHKERR PetscOptionsScalar("-adapt_max_smoother_alpha",
385  "minimal smoother alpha", "", maxAlpha,
386  &maxAlpha, PETSC_NULL);
387  incrementAlpha = 10;
388  CHKERR PetscOptionsScalar("-adapt_incitement_smoother_alpha",
389  "minimal smoother alpha", "", incrementAlpha,
390  &incrementAlpha, PETSC_NULL);
391 
392  CHKERR PetscOptionsScalar(
393  "-adapt_desired_alpha",
394  "Desired alpha (allowed error for gc calculation)", "", dAlpha,
395  &dAlpha, PETSC_NULL);
396 
397  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
398  "### Input parameter: -adapt_min_smoother_alpha %6.4e",
399  minAlpha);
400  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
401  "### Input parameter: -adapt_max_smoother_alpha %6.4e",
402  maxAlpha);
403  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
404  "### Input parameter: -adapt_desired_alpha %6.4e", dAlpha);
405 
406  ierr = PetscOptionsEnd();
407 
408  MoFEMFunctionReturn(0);
409  }
410 
411  MoFEMErrorCode operator()(DM dm, SNES snes, Vec q) const {
412  MoFEMFunctionBegin;
413 
414  if (cP.smootherFe->smootherData.sTabilised) {
415  SNESConvergedReason reason;
416  CHKERR SNESGetConvergedReason(snes, &reason);
417  if (reason < 0) {
418  const double old_alpha = cP.smootherAlpha;
419  cP.smootherAlpha =
420  std::min(incrementAlpha * cP.smootherAlpha, maxAlpha);
421  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
422  "Increase smoothing alpha = %6.5g (old %6.5e)",
423  cP.smootherAlpha, old_alpha);
424  } else {
425  auto f_smoothing = smartVectorDuplicate(q);
426  auto f_griffith = smartVectorDuplicate(q);
427 
428  CHKERR cP.calculateGriffithForce(dm, cP.gC ? cP.gC : 1, f_griffith,
429  VERBOSE, false);
430  // Calculate material forces
431  CHKERR DMoFEMMeshToLocalVector(dm, q, INSERT_VALUES, SCATTER_FORWARD);
432  CHKERR cP.calculateSmoothingForcesDM(dm, q, f_smoothing, VERBOSE,
433  false);
434  CHKERR cP.calculateSmoothingForceFactor(VERBOSE, false);
435 
436  double max_alpha = 0;
437  for (const auto &kv : cP.mapSmoothingForceFactor)
438  max_alpha = std::max(max_alpha, kv.second);
439 
440  const double old_alpha = cP.smootherAlpha;
441  cP.smootherAlpha *= pow(dAlpha / max_alpha, gAmma);
442  cP.smootherAlpha = std::max(cP.smootherAlpha, minAlpha);
443  MOFEM_LOG_C("CPSolverWorld", Sev::inform,
444  "Set smoothing alpha = %6.5g (old %6.5e)",
445  cP.smootherAlpha, old_alpha);
446  }
447 
448  cP.volumeLengthAdouble->aLpha = cP.smootherAlpha;
449  cP.volumeLengthDouble->aLpha = cP.smootherAlpha;
450  for (auto &kv : cP.surfaceConstrain)
451  kv.second->aLpha = cP.smootherAlpha;
452 
453  for (auto &kv : cP.edgeConstrains)
454  kv.second->aLpha = cP.smootherAlpha;
455 
456  cP.bothSidesConstrains->aLpha = cP.smootherAlpha;
457  cP.bothSidesContactConstrains->aLpha = cP.smootherAlpha;
458  }
459 
460  MoFEMFunctionReturn(0);
461  }
462  };
463 
464  const MoFEM::Problem *problem_cp_ptr;
465  CHKERR DMMoFEMGetProblemPtr(dm_crack_propagation, &problem_cp_ptr);
466 
467  auto cp_arc_ctx = boost::make_shared<ArcLengthCtx>(
468  m_field, problem_cp_ptr->getName(), "LAMBDA_ARC_LENGTH");
469 
470  CHKERR cP.mField.getInterface<VecManager>()->vecCreateGhost(
471  problem_cp_ptr->getName(), COL, cp_arc_ctx->x0);
472  cp_arc_ctx->dx = smartVectorDuplicate(cp_arc_ctx->x0);
473  arc_snes_ctx = boost::make_shared<CrackPropagation::ArcLengthSnesCtx>(
474  m_field, problem_cp_ptr->getName(), cp_arc_ctx);
475  CHKERR DMMoFEMSetSnesCtx(dm_crack_propagation, arc_snes_ctx);
476 
477  boost::shared_ptr<FEMethod> arc_method =
478  cP.getFrontArcLengthControl(cp_arc_ctx);
479 
480  // Add operators to DM SNES
481 
482  CHKERR cP.addPropagationFEInstancesToSnes(dm_crack_propagation, arc_method,
483  cp_arc_ctx, surface_ids, 1, false);
484 
485  {
486  // create vectors and matrix
487  auto m = smartCreateDMMatrix(dm_crack_propagation);
488  auto f = smartCreateDMVector(dm_crack_propagation);
489  auto q = smartVectorDuplicate(f);
490 
491  // Create snes
492  auto snes_crack_propagation = createSNES(m_field.get_comm());
493  CHKERR SNESAppendOptionsPrefix(snes_crack_propagation, "propagation_");
494 
495  for (int ll = 0; ll != cP.getNbLoadSteps(); ll++) {
496 
497  int run_step = cP.getNbLoadSteps() * cut_mesh_it + ll;
498  int step = cP.startStep++;
499 
500  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\n\nLoad step %d ( %d ) [ %d ]\n\n",
501  ll, run_step, step);
502 
503  // set vector values from mesh data
504  CHKERR cP.zeroLambdaFields();
505  CHKERR DMoFEMMeshToLocalVector(dm_crack_propagation, q, INSERT_VALUES,
506  SCATTER_FORWARD);
507 
508  CHKERR VecCopy(q, cp_arc_ctx->x0);
509  CHKERR cp_arc_ctx->setAlphaBeta(cP.arcAlpha, cP.arcBeta);
510 
511  auto solve_problem = [&]() {
512  MoFEMFunctionBegin;
513 
514  const double arc_s = cP.arcS;
515 
516  cP.arcS *= cP.fractionOfFixedNodes;
517  CHKERR cp_arc_ctx->setS(cP.arcS);
518 
519  CHKERR cP.solvePropagationDM(dm_crack_propagation, dm_elastic,
520  snes_crack_propagation, m, q, f);
521  cP.arcS = arc_s;
522  CHKERR cp_arc_ctx->setS(cP.arcS);
523 
524  // Adapt step
525  CHKERR StepAdaptivity(cP, 6, 0.5)(
526  dm_crack_propagation, snes_crack_propagation, cp_arc_ctx, cP.arcS);
527  CHKERR SmoothingAlphaAdaptivity(cP, 0.25, 1e-2,
528  cP.initialSmootherAlpha * 1e-3)(
529  dm_crack_propagation, snes_crack_propagation, q);
530  MoFEMFunctionReturn(0);
531  };
532 
533  CHKERR solve_problem();
534 
535  auto if_not_converged = [&](Vec q) {
536  MoFEMFunctionBegin;
537  SNESConvergedReason reason;
538  CHKERR SNESGetConvergedReason(snes_crack_propagation, &reason);
539  if (reason < 0) {
540  CHKERR VecCopy(cp_arc_ctx->x0, q);
541  CHKERR DMoFEMMeshToLocalVector(dm_crack_propagation, q, INSERT_VALUES,
542  SCATTER_FORWARD);
543  }
544  MoFEMFunctionReturn(0);
545  };
546  CHKERR if_not_converged(q);
547 
548  auto scatter_forward_vec = [](Vec q) {
549  MoFEMFunctionBegin;
550  CHKERR VecGhostUpdateBegin(q, INSERT_VALUES, SCATTER_FORWARD);
551  CHKERR VecGhostUpdateEnd(q, INSERT_VALUES, SCATTER_FORWARD);
552  MoFEMFunctionReturn(0);
553  };
554  CHKERR scatter_forward_vec(q);
555 
556  // Save material displacements on the mesh
557  CHKERR VecAYPX(cp_arc_ctx->x0, -1, q);
558  CHKERR scatter_forward_vec(cP.getArcCtx()->x0);
559 
560  CrackPropagation::PostProcVertexMethod ent_front_disp(
561  cP.mField, cp_arc_ctx->x0, "W");
562  CHKERR cP.mField.loop_dofs(problem_cp_ptr->getName(),
563  "MESH_NODE_POSITIONS", ROW, ent_front_disp, 0,
564  cP.mField.get_comm_rank());
565 
566  // Displace mesh
567  CHKERR cP.setCoordsFromField("MESH_NODE_POSITIONS");
568 
569  // calculate energy
570  SNESConvergedReason reason;
571  CHKERR SNESGetConvergedReason(snes_crack_propagation, &reason);
572  std::string step_msg;
573  if (reason < 0) {
574  step_msg = "Not Converged Propagation step " +
575  boost::lexical_cast<std::string>(step) + " ( " +
576  boost::lexical_cast<std::string>(ll) + " )";
577  } else {
578  step_msg = "Propagation step " +
579  boost::lexical_cast<std::string>(step) + " ( " +
580  boost::lexical_cast<std::string>(ll) + " )";
581  }
582  CHKERR cP.calculateElasticEnergy(dm_elastic, step_msg);
583  // calculate griffith energy
584  CHKERR calculateGc(dm_material_forces, dm_surface_projection,
585  dm_crack_srf_area, surface_ids);
586  CHKERR cP.updateMaterialFixedNode(false, true, false);
587 
588  // post-processing
589  CHKERR cP.postProcessDM(dm_elastic, step, "ELASTIC",
590  false || run_step == 0);
591  // save crack front
592  CHKERR cP.savePositionsOnCrackFrontDM(dm_elastic, q, QUIET, false);
593 
594  // Save positions
595  auto save_positions_on_mesh_tags_from_field =
596  [this, problem_cp_ptr, &q](const std::string field_name) {
597  MoFEMFunctionBegin;
598  CrackPropagation::PostProcVertexMethod ent_post_proc(cP.mField, q,
599  field_name);
600  CHKERR cP.mField.loop_dofs(problem_cp_ptr->getName(), field_name,
601  ROW, ent_post_proc, 0,
602  cP.mField.get_comm_size());
603  MoFEMFunctionReturn(0);
604  };
605  CHKERR save_positions_on_mesh_tags_from_field("MESH_NODE_POSITIONS");
606  CHKERR save_positions_on_mesh_tags_from_field("SPATIAL_POSITION");
607 
608  // write state
609  CHKERR cP.writeCrackFont(cP.mapBitLevel["spatial_domain"], step);
610 
611  auto save_restart_file = [&]() {
612  auto bit_interface = cP.mField.getInterface<BitRefManager>();
613 
614  MoFEMFunctionBegin;
615 
616  auto bit = BitRefLevel().set(BITREFLEVEL_SIZE - 1);
617  bit.set(BITREFLEVEL_SIZE - 2);
618  auto meshset_ptr = get_temp_meshset_ptr(m_field.get_moab());
619 
620  Range ents_on_last_bits;
621  CHKERR bit_interface->getEntitiesByRefLevel(bit, BitRefLevel().set(),
622  ents_on_last_bits);
623  CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
624  ents_on_last_bits);
625 
626  Range ents_spatial_domain;
627  CHKERR bit_interface->getEntitiesByRefLevel(
628  cP.mapBitLevel.at("spatial_domain"), BitRefLevel().set(),
629  ents_spatial_domain);
630 
631  Range entities_on_skin;
632  entities_on_skin.merge(cP.bodySkin);
633  entities_on_skin.merge(cP.crackFaces);
634  entities_on_skin.merge(cP.crackFront);
635 
636  entities_on_skin = intersect(ents_spatial_domain, entities_on_skin);
637  CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
638  entities_on_skin);
639  CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
640  ents_spatial_domain);
641 
642  Range entities_bit;
643  entities_bit.merge(ents_on_last_bits);
644  entities_bit.merge(ents_spatial_domain);
645 
646  auto add_field = [&](auto field_name, auto &tags_list) {
647  MoFEMFunctionBegin;
648  auto field_ptr = cP.mField.get_field_structure(field_name);
649  auto meshset = field_ptr->getMeshset();
650  CHKERR cP.mField.get_moab().add_entities(*meshset_ptr, &meshset, 1);
651 
652  Tag th_field_id;
653  CHKERR cP.mField.get_moab().tag_get_handle("_FieldId", th_field_id);
654  Tag th_field_base;
655  CHKERR cP.mField.get_moab().tag_get_handle("_FieldBase",
656  th_field_base);
657  Tag th_field_space;
658  CHKERR cP.mField.get_moab().tag_get_handle("_FieldSpace",
659  th_field_space);
660  Tag th_field_name;
661  CHKERR cP.mField.get_moab().tag_get_handle("_FieldName",
662  th_field_name);
663  Tag th_field_name_data_name_prefix;
664  CHKERR cP.mField.get_moab().tag_get_handle(
665  "_FieldName_DataNamePrefix", th_field_name_data_name_prefix);
666 
667  tags_list.push_back(th_field_id);
668  tags_list.push_back(th_field_base);
669  tags_list.push_back(th_field_space);
670  tags_list.push_back(th_field_name);
671  tags_list.push_back(th_field_name_data_name_prefix);
672  tags_list.push_back(field_ptr->th_FieldDataVerts);
673  tags_list.push_back(field_ptr->th_FieldData);
674  tags_list.push_back(field_ptr->th_AppOrder);
675  tags_list.push_back(field_ptr->th_FieldRank);
676  MoFEMFunctionReturn(0);
677  };
678 
679  auto add_fe = [&](auto fe_name, auto &tags_list) {
680  MoFEMFunctionBegin;
681  auto arc_length_finite_element =
682  cP.mField.get_finite_element_meshset(fe_name);
683  CHKERR cP.mField.get_moab().add_entities(
684  *meshset_ptr, &arc_length_finite_element, 1);
685  Tag th_fe_id;
686  CHKERR cP.mField.get_moab().tag_get_handle("_FEId", th_fe_id);
687  Tag th_fe_name;
688  CHKERR cP.mField.get_moab().tag_get_handle("_FEName", th_fe_name);
689  Tag th_fe_id_col, th_fe_id_row, th_fe_id_data;
690  CHKERR cP.mField.get_moab().tag_get_handle("_FEIdCol", th_fe_id_col);
691  CHKERR cP.mField.get_moab().tag_get_handle("_FEIdRow", th_fe_id_row);
692  CHKERR cP.mField.get_moab().tag_get_handle("_FEIdData",
693  th_fe_id_data);
694  tags_list.push_back(th_fe_id);
695  tags_list.push_back(th_fe_name);
696  tags_list.push_back(th_fe_id_col);
697  tags_list.push_back(th_fe_id_row);
698  tags_list.push_back(th_fe_id_data);
699  MoFEMFunctionReturn(0);
700  };
701 
702  auto arc_field_meshset =
703  cP.mField.get_field_meshset("LAMBDA_ARC_LENGTH");
704  auto arc_length_finite_element =
705  cP.mField.get_finite_element_meshset("ARC_LENGTH");
706  CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
707  &arc_field_meshset, 1);
708  CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
709  &arc_length_finite_element, 1);
710 
711  auto get_tags_list = [&]() {
712  Tag th_griffith_force;
713  CHKERR cP.mField.get_moab().tag_get_handle("GRIFFITH_FORCE",
714  th_griffith_force);
715  Tag th_griffith_force_projected;
716  CHKERR cP.mField.get_moab().tag_get_handle(
717  "GRIFFITH_FORCE_PROJECTED", th_griffith_force_projected);
718  Tag th_w;
719  CHKERR cP.mField.get_moab().tag_get_handle("W", th_w);
720  Tag th_material_force;
721  CHKERR cP.mField.get_moab().tag_get_handle("MATERIAL_FORCE",
722  th_material_force);
723  Tag th_interface_side;
724  CHKERR cP.mField.get_moab().tag_get_handle("INTERFACE_SIDE",
725  th_interface_side);
726  Tag th_org_pos;
727  CHKERR cP.mField.get_moab().tag_get_handle("ORG_POSITION",
728  th_org_pos);
729  Tag th_version;
730  CHKERR cP.mField.get_moab().tag_get_handle("MOFEM_VERSION",
731  th_version);
732 
733  std::vector<Tag> tags_list;
734  tags_list.push_back(bit_interface->get_th_RefBitLevel());
735 
736  tags_list.push_back(th_version);
737  tags_list.push_back(th_griffith_force);
738  tags_list.push_back(th_griffith_force_projected);
739  tags_list.push_back(th_w);
740  tags_list.push_back(th_material_force);
741  tags_list.push_back(th_interface_side);
742  tags_list.push_back(th_org_pos);
743  return tags_list;
744  };
745 
746  auto add_meshsets = [&](Range bit, auto &tags_list) {
747  auto meshsets_mng = cP.mField.getInterface<MeshsetsManager>();
748  std::vector<boost::shared_ptr<TempMeshset>> meshses_tmp_list;
749  auto &list = meshsets_mng->getMeshsetsMultindex();
750  for (auto &m : list) {
751  meshses_tmp_list.push_back(
752  get_temp_meshset_ptr(cP.mField.get_moab()));
753  EntityHandle new_meshset = *meshses_tmp_list.back();
754  auto meshset = m.getMeshset();
755  std::vector<Tag> tmp_tags_list;
756  CHKERR cP.mField.get_moab().tag_get_tags_on_entity(meshset,
757  tmp_tags_list);
758  Range ents;
759  CHKERR cP.mField.get_moab().get_entities_by_handle(meshset, ents,
760  true);
761  ents = intersect(ents, bit);
762  CHKERR cP.mField.get_moab().add_entities(new_meshset, ents);
763  for (auto t : tmp_tags_list) {
764  void *tag_vals[1];
765  int tag_size[1];
766  CHKERR cP.mField.get_moab().tag_get_by_ptr(
767  t, &meshset, 1, (const void **)tag_vals, tag_size);
768  CHKERR cP.mField.get_moab().tag_set_by_ptr(t, &new_meshset, 1,
769  tag_vals, tag_size);
770  }
771 
772  std::vector<std::string> remove_tags;
773  remove_tags.push_back("CATEGORY");
774  remove_tags.push_back("OBB");
775  remove_tags.push_back("ROOTSETSURF");
776  remove_tags.push_back("__PARALLEL_");
777 
778  for (auto t : tmp_tags_list) {
779 
780  std::string tag_name;
781  CHKERR cP.mField.get_moab().tag_get_name(t, tag_name);
782  bool add = true;
783 
784  for (auto &p : remove_tags) {
785  if (tag_name.compare(0, p.size(), p) == 0) {
786  add = false;
787  break;
788  }
789  }
790 
791  if (add)
792  tags_list.push_back(t);
793  }
794  }
795  return meshses_tmp_list;
796  };
797 
798  auto tags_list = get_tags_list();
799  CHKERR add_field("LAMBDA_ARC_LENGTH", tags_list);
800  CHKERR add_fe("ARC_LENGTH", tags_list);
801 
802  auto meshses_tmp_list = add_meshsets(entities_bit, tags_list);
803  for (auto &m_ptr : meshses_tmp_list) {
804  EntityHandle m = *m_ptr;
805  CHKERR cP.mField.get_moab().add_entities(*meshset_ptr, &m, 1);
806  }
807 
808  {
809  Tag gid_tag = cP.mField.get_moab().globalId_tag();
810  int negative = -1;
811  CHKERR cP.mField.get_moab().tag_clear_data(gid_tag, entities_bit,
812  &negative);
813  Range verts_ents = entities_bit.subset_by_type(MBVERTEX);
814  int gid = 1; // moab indexing from 1
815  for (auto v : verts_ents) {
816  CHKERR cP.mField.get_moab().tag_set_data(gid_tag, &v, 1, &gid);
817  gid++;
818  }
819  }
820 
821  std::sort(tags_list.begin(), tags_list.end());
822  auto new_end = std::unique(tags_list.begin(), tags_list.end());
823  tags_list.resize(std::distance(tags_list.begin(), new_end));
824 
825  for (auto t : tags_list) {
826  std::string tag_name;
827  CHKERR cP.mField.get_moab().tag_get_name(t, tag_name);
828  MOFEM_LOG("CPSolverWorld", Sev::noisy)
829  << "Add tag to restart file: " << tag_name;
830  }
831 
832  EntityHandle save_meshset = *meshset_ptr;
833  const std::string file_name =
834  "restart_" + boost::lexical_cast<std::string>(step) + ".h5m";
835  CHKERR cP.mField.get_moab().write_file(file_name.c_str(), "MOAB", "",
836  &save_meshset, 1, &tags_list[0],
837  tags_list.size());
838 
839  MoFEMFunctionReturn(0);
840  };
841 
842  if (cP.mField.get_comm_rank() == 0)
843  CHKERR save_restart_file();
844  }
845  }
846 
847  MoFEMFunctionReturn(0);
848 }

Member Data Documentation

◆ cP

CrackPropagation& FractureMechanics::CPSolvers::cP

Definition at line 39 of file CPSolvers.hpp.

The documentation for this struct was generated from the following files:
FractureMechanics::CrackPropagation::calculateFrontProjectionMatrix
MoFEMErrorCode calculateFrontProjectionMatrix(DM dm_surface, DM dm_project, const int verb=QUIET, const bool debug=false)
assemble crack front projection matrix (that constrains crack area growth)
Definition: CrackPropagation.cpp:7637
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
FractureMechanics::CrackPropagation::feGriffithConstrainsLhs
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsLhs
Definition: CrackPropagation.hpp:1144
FractureMechanics::CrackPropagation::writeCrackFont
MoFEMErrorCode writeCrackFont(const BitRefLevel &bit, const int step)
Definition: CrackPropagation.cpp:8093
FractureMechanics::CPSolvers::CPSolvers
CPSolvers(CrackPropagation &cp)
Definition: CPSolvers.cpp:56
MoFEM::DMMoFEMSetSnesCtx
PetscErrorCode DMMoFEMSetSnesCtx(DM dm, boost::shared_ptr< MoFEM::SnesCtx > snes_ctx)
Set MoFEM::SnesCtx data structure.
Definition: DMMoFEM.cpp:1115
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreInterface::get_finite_element_meshset
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
MoFEM::smartGetDMSnesCtx
DEPRECATED auto smartGetDMSnesCtx(DM dm)
Definition: DMMoFEM.hpp:1054
FractureMechanics::CrackPropagation::parallelReadAndBroadcast
static bool parallelReadAndBroadcast
Definition: CrackPropagation.hpp:81
SimpleArcLengthControl
Definition: ArcLengthTools.hpp:334
FractureMechanics::CrackPropagation::calculateSmoothingForcesDM
MoFEMErrorCode calculateSmoothingForcesDM(DM dm, Vec q, Vec f, const int verb=QUIET, const bool debug=false)
assemble smoothing forces, by running material finite element instance
Definition: CrackPropagation.cpp:7476
FractureMechanics::CrackPropagation::calculateSurfaceProjectionMatrix
MoFEMErrorCode calculateSurfaceProjectionMatrix(DM dm_front, DM dm_project, const std::vector< int > &ids, const int verb=QUIET, const bool debug=false)
assemble projection matrices
Definition: CrackPropagation.cpp:7563
EntityHandle
FractureMechanics::CrackPropagation::feLhs
boost::shared_ptr< CrackFrontElement > feLhs
Integrate elastic FE.
Definition: CrackPropagation.hpp:1065
FractureMechanics::CrackPropagation::feMaterialAnaliticalTraction
boost::shared_ptr< NeumannForcesSurface::MyTriangleFE > feMaterialAnaliticalTraction
Surface elment to calculate tractions in material space.
Definition: CrackPropagation.hpp:1080
FractureMechanics::CrackPropagation::crackFaces
Range crackFaces
Definition: CrackPropagation.hpp:1160
FractureMechanics::CrackPropagation::zeroLambdaFields
MoFEMErrorCode zeroLambdaFields()
Zero fields with lagrange multipliers.
Definition: CrackPropagation.cpp:9553
MoFEM::CoreInterface::set_moab_interface
virtual MoFEMErrorCode set_moab_interface(moab::Interface &new_moab, int verb=VERBOSE)=0
Set the moab interface object.
FractureMechanics::CrackPropagation::maxG1
double maxG1
Definition: CrackPropagation.hpp:1184
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
FractureMechanics::CrackPropagation::mapBitLevel
std::map< std::string, BitRefLevel > mapBitLevel
Definition: CrackPropagation.hpp:180
ArcLengthCtx
Store variables for ArcLength analysis.
Definition: ArcLengthTools.hpp:65
FractureMechanics::CrackPropagation::feGriffithConstrainsRhs
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsRhs
Definition: CrackPropagation.hpp:1142
FractureMechanics::CrackPropagation::surfaceConstrain
map< int, boost::shared_ptr< SurfaceSlidingConstrains > > surfaceConstrain
Definition: CrackPropagation.hpp:1128
MoFEM::createSNES
auto createSNES(MPI_Comm comm)
Definition: PetscSmartObj.hpp:251
MoFEM::Field::getMeshset
EntityHandle getMeshset() const
Get field meshset.
Definition: FieldMultiIndices.hpp:122
FractureMechanics::CrackPropagation::feMaterialLhs
boost::shared_ptr< CrackFrontElement > feMaterialLhs
Integrate material stresses, assemble matrix.
Definition: CrackPropagation.hpp:1070
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
FractureMechanics::CPSolvers::solveElastic
MoFEMErrorCode solveElastic(DM dm_elastic, Vec init, const double load_factor)
Solve elastic problem.
Definition: CPSolvers.cpp:90
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FractureMechanics::CrackPropagation::bothSidesConstrains
boost::shared_ptr< BothSurfaceConstrains > bothSidesConstrains
Definition: CrackPropagation.hpp:1131
FractureMechanics::CrackPropagation::setSingularDofs
MoFEMErrorCode setSingularDofs(const string field_name, const int verb=QUIET)
set singular dofs (on edges adjacent to crack front) from geometry
Definition: CrackPropagation.cpp:9334
FractureMechanics::CrackPropagation::solveEigenStressProblem
PetscBool solveEigenStressProblem
Solve eigen problem.
Definition: CrackPropagation.hpp:163
FractureMechanics::CrackPropagation::getArcCtx
boost::shared_ptr< ArcLengthCtx > & getArcCtx()
Definition: CrackPropagation.hpp:907
FractureMechanics::broadcast_entities
MoFEMErrorCode broadcast_entities(moab::Interface &moab, moab::Interface &moab_tmp, const int from_proc, Range &entities, const bool adjacencies, const bool tags)
Definition: CrackPropagation.cpp:149
FractureMechanics::CrackPropagation::assembleCouplingForcesDM
MoFEMErrorCode assembleCouplingForcesDM(DM dm, const int verb=QUIET, const bool debug=false)
assemble coupling element instances
Definition: CrackPropagation.cpp:6178
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
FractureMechanics::CrackPropagation::feGriffithForceRhs
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceRhs
Definition: CrackPropagation.hpp:1137
SolveFunctionBegin
#define SolveFunctionBegin
Definition: CPSolvers.cpp:43
FractureMechanics::CrackPropagation::arcAlpha
double arcAlpha
Definition: CrackPropagation.hpp:1211
FractureMechanics::CrackPropagation::initialSmootherAlpha
double initialSmootherAlpha
Definition: CrackPropagation.hpp:1201
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:527
FractureMechanics::CrackPropagation::addPropagationFEInstancesToSnes
MoFEMErrorCode addPropagationFEInstancesToSnes(DM dm, boost::shared_ptr< FEMethod > arc_method, boost::shared_ptr< ArcLengthCtx > arc_ctx, const std::vector< int > &surface_ids, const int verb=QUIET, const bool debug=false)
add finite element to SNES for crack propagation problem
Definition: CrackPropagation.cpp:6921
FractureMechanics::CrackPropagation::getEigenArcCtx
boost::shared_ptr< ArcLengthCtx > & getEigenArcCtx()
Definition: CrackPropagation.hpp:908
FractureMechanics::CrackPropagation::setSpatialPositionFromCoords
MoFEMErrorCode setSpatialPositionFromCoords()
set spatial field from nodes
Definition: CrackPropagation.cpp:9330
FractureMechanics::CPSolvers::calculateGc
MoFEMErrorCode calculateGc(DM dm_material_forces, DM dm_surface_projection, DM dm_crack_srf_area, const std::vector< int > &surface_ids)
Solve for energy realease rate.
Definition: CPSolvers.cpp:186
MoFEM::CoreInterface::get_field_meshset
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
FractureMechanics::CrackPropagation::getNbLoadSteps
int & getNbLoadSteps()
Definition: CrackPropagation.hpp:208
FractureMechanics::CrackPropagation::feEnergy
boost::shared_ptr< CrackFrontElement > feEnergy
Integrate energy.
Definition: CrackPropagation.hpp:1071
FractureMechanics::CrackPropagation::calculateGriffithForce
MoFEMErrorCode calculateGriffithForce(DM dm, const double gc, Vec f_griffith, const int verb=QUIET, const bool debug=false)
calculate Griffith (driving) force
Definition: CrackPropagation.cpp:7962
FractureMechanics::CrackPropagation::volumeLengthAdouble
boost::shared_ptr< VolumeLengthQuality< adouble > > volumeLengthAdouble
Definition: CrackPropagation.hpp:1117
MoFEM::CoreInterface::check_problem
virtual bool check_problem(const std::string name)=0
check if problem exist
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
FractureMechanics::CrackPropagation::mField
MoFEM::Interface & mField
Definition: CrackPropagation.hpp:106
FractureMechanics::CrackPropagation::solvePropagationDM
MoFEMErrorCode solvePropagationDM(DM dm, DM dm_elastic, SNES snes, Mat m, Vec q, Vec f)
solve crack propagation problem
Definition: CrackPropagation.cpp:8494
FractureMechanics::CrackPropagation::setMaterialPositionFromCoords
MoFEMErrorCode setMaterialPositionFromCoords()
set material field from nodes
Definition: CrackPropagation.cpp:9326
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
VERBOSE
@ VERBOSE
Definition: definitions.h:209
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FractureMechanics::CrackPropagation::calculateMaterialForcesDM
MoFEMErrorCode calculateMaterialForcesDM(DM dm, Vec q, Vec f, const int verb=QUIET, const bool debug=false)
assemble material forces, by running material finite element instance
Definition: CrackPropagation.cpp:7416
FractureMechanics::CrackPropagation::updateMaterialFixedNode
MoFEMErrorCode updateMaterialFixedNode(const bool fix_front, const bool fix_small_g, const bool debug=false)
Update fixed nodes.
Definition: CrackPropagation.cpp:6898
FractureMechanics::CrackPropagation::fixMaterialEnts
boost::shared_ptr< DirichletFixFieldAtEntitiesBc > fixMaterialEnts
Definition: CrackPropagation.hpp:1134
MoFEM::UnknownInterface::setFileVersion
static MoFEMErrorCode setFileVersion(moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
Get database major version.
Definition: UnknownInterface.cpp:41
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
FractureMechanics::CrackPropagation::getNbCutSteps
int & getNbCutSteps()
Definition: CrackPropagation.hpp:212
MoFEM::SmartPetscObj::use_count
int use_count() const
Definition: PetscSmartObj.hpp:105
FractureMechanics::CrackPropagation::arcS
double arcS
Definition: CrackPropagation.hpp:1213
a
constexpr double a
Definition: approx_sphere.cpp:30
FractureMechanics::CrackPropagation::calculateSmoothingForceFactor
MoFEMErrorCode calculateSmoothingForceFactor(const int verb=QUIET, const bool debug=true)
Definition: CrackPropagation.cpp:7524
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1857
FractureMechanics::CrackPropagation::bothSidesContactConstrains
boost::shared_ptr< BothSurfaceConstrains > bothSidesContactConstrains
Definition: CrackPropagation.hpp:1133
FractureMechanics::CrackPropagation::closeCrackConstrains
boost::shared_ptr< BothSurfaceConstrains > closeCrackConstrains
Definition: CrackPropagation.hpp:1132
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
FractureMechanics::CrackPropagation::assembleSmootherForcesDM
MoFEMErrorCode assembleSmootherForcesDM(DM dm, const std::vector< int > ids, const int verb=QUIET, const bool debug=false)
create smoothing element instance
Definition: CrackPropagation.cpp:5975
FractureMechanics::CrackPropagation::smootherFe
boost::shared_ptr< Smoother > smootherFe
Definition: CrackPropagation.hpp:1111
COL
@ COL
Definition: definitions.h:123
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::smartCreateDMMatrix
DEPRECATED auto smartCreateDMMatrix(DM dm)
Definition: DMMoFEM.hpp:1011
FractureMechanics::CrackPropagation::getFrontArcLengthControl
boost::shared_ptr< FEMethod > getFrontArcLengthControl(boost::shared_ptr< ArcLengthCtx > arc_ctx)
Definition: CrackPropagation.hpp:913
FractureMechanics::CrackPropagation::crackFront
Range crackFront
Definition: CrackPropagation.hpp:1163
MoFEM::smartVectorDuplicate
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
Definition: PetscSmartObj.hpp:226
FractureMechanics::CrackPropagation::feSpringLhsPtr
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringLhsPtr
Definition: CrackPropagation.hpp:1149
FractureMechanics::CrackPropagation::fileVersion
Version fileVersion
Definition: CrackPropagation.hpp:86
FractureMechanics::CrackPropagation::griffithForceElement
boost::shared_ptr< GriffithForceElement > griffithForceElement
Definition: CrackPropagation.hpp:1136
MoFEM::smartCreateDMVector
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1026
FractureMechanics::CrackPropagation::arcBeta
double arcBeta
Definition: CrackPropagation.hpp:1212
FractureMechanics::CPSolvers::getOptions
MoFEMErrorCode getOptions()
Definition: CPSolvers.cpp:80
FractureMechanics::CrackPropagation::feCouplingElasticLhs
boost::shared_ptr< CrackFrontElement > feCouplingElasticLhs
FE instance to assemble coupling terms.
Definition: CrackPropagation.hpp:1107
FractureMechanics::CrackPropagation::gC
double gC
Griffith energy.
Definition: CrackPropagation.hpp:111
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FractureMechanics::CrackPropagation::assembleElasticDM
MoFEMErrorCode assembleElasticDM(const std::string mwls_stress_tag_name, const int verb=QUIET, const bool debug=false)
create elastic finite element instance for spatial assembly
Definition: CrackPropagation.cpp:3795
FractureMechanics::CrackPropagation::calculateElasticEnergy
MoFEMErrorCode calculateElasticEnergy(DM dm, const std::string msg="")
assemble elastic part of matrix, by running elastic finite element instance
Definition: CrackPropagation.cpp:7280
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
FractureMechanics::CrackPropagation::crackAccelerationFactor
double crackAccelerationFactor
Definition: CrackPropagation.hpp:1328
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
FractureMechanics::CrackPropagation::postProcessDM
MoFEMErrorCode postProcessDM(DM dm, const int step, const std::string fe_name, const bool approx_internal_stress)
post-process results for elastic solution
Definition: CrackPropagation.cpp:9010
FractureMechanics::CrackPropagation::crackOrientation
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > crackOrientation
Definition: CrackPropagation.hpp:1125
FractureMechanics::CrackPropagation::feRhs
boost::shared_ptr< CrackFrontElement > feRhs
Integrate elastic FE.
Definition: CrackPropagation.hpp:1066
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
FractureMechanics::CPSolvers::solvePropagation
MoFEMErrorCode solvePropagation(DM dm_crack_propagation, DM dm_elastic, DM dm_material, DM dm_material_forces, DM dm_surface_projection, DM dm_crack_srf_area, const std::vector< int > &surface_ids, int cut_mesh_it=0, const bool set_cut_surface=false)
Solve crack propagation problem.
Definition: CPSolvers.cpp:247
FractureMechanics::CrackPropagation::feSmootherLhs
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherLhs
Integrate smoothing operators.
Definition: CrackPropagation.hpp:1115
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FractureMechanics::CrackPropagation::calculateReleaseEnergy
MoFEMErrorCode calculateReleaseEnergy(DM dm, Vec f_material_proj, Vec f_griffith_proj, Vec f_lambda, const double gc, const int verb=QUIET, const bool debug=true)
calculate release energy
Definition: CrackPropagation.cpp:7731
MoFEM::CoreInterface::clear_database
virtual MoFEMErrorCode clear_database(int verb=DEFAULT_VERBOSITY)=0
Clear database.
FractureMechanics::CrackPropagation::bodySkin
Range bodySkin
Definition: CrackPropagation.hpp:1159
FractureMechanics::CrackPropagation::mwlsEigenStressTagName
std::string mwlsEigenStressTagName
Name of tag with eigen stresses.
Definition: CrackPropagation.hpp:134
FractureMechanics::CrackPropagation::solveElasticDM
MoFEMErrorCode solveElasticDM(DM dm, SNES snes, Mat m, Vec q, Vec f, bool snes_set_up, Mat *shell_m)
solve elastic problem
Definition: CrackPropagation.cpp:8224
FractureMechanics::CrackPropagation::savePositionsOnCrackFrontDM
MoFEMErrorCode savePositionsOnCrackFrontDM(DM dm, Vec q, const int verb=QUIET, const bool debug=false)
Definition: CrackPropagation.cpp:8064
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
FractureMechanics::CrackPropagation::projSurfaceCtx
boost::shared_ptr< ConstrainMatrixCtx > projSurfaceCtx
Data structure to project on the body surface.
Definition: CrackPropagation.hpp:1073
Range
MoFEM::MeshsetsManager::getMeshsetsMultindex
CubitMeshSet_multiIndex & getMeshsetsMultindex()
Definition: MeshsetsManager.hpp:229
FractureMechanics::CrackPropagation::feGriffithConstrainsDelta
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsDelta > feGriffithConstrainsDelta
Definition: CrackPropagation.hpp:1140
FractureMechanics::CrackPropagation::feSpringRhsPtr
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringRhsPtr
Definition: CrackPropagation.hpp:1150
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FractureMechanics::CrackPropagation::moabCommWorld
boost::shared_ptr< WrapMPIComm > moabCommWorld
Definition: CrackPropagation.hpp:104
FractureMechanics::CrackPropagation::volumeLengthDouble
boost::shared_ptr< VolumeLengthQuality< double > > volumeLengthDouble
Definition: CrackPropagation.hpp:1116
FractureMechanics::CrackPropagation::projectMaterialForcesDM
MoFEMErrorCode projectMaterialForcesDM(DM dm_project, Vec f, Vec f_proj, const int verb=QUIET, const bool debug=false)
project material forces along the crack elongation direction
Definition: CrackPropagation.cpp:7908
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
FractureMechanics::CrackPropagation::elasticFe
boost::shared_ptr< NonlinearElasticElement > elasticFe
Definition: CrackPropagation.hpp:1063
FractureMechanics::CrackPropagation::contactOrientation
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > contactOrientation
Definition: CrackPropagation.hpp:1127
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
FractureMechanics::CrackPropagation::feMaterialRhs
boost::shared_ptr< CrackFrontElement > feMaterialRhs
Integrate material stresses, assemble vector.
Definition: CrackPropagation.hpp:1068
BITREFLEVEL_SIZE
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
FractureMechanics::CrackPropagation::getLoadScale
double & getLoadScale()
Definition: CrackPropagation.hpp:210
FractureMechanics::CrackPropagation::mapSmoothingForceFactor
map< EntityHandle, double > mapSmoothingForceFactor
Definition: CrackPropagation.hpp:1195
FractureMechanics::CrackPropagation::mwlsStressTagName
std::string mwlsStressTagName
Name of tag with internal stresses.
Definition: CrackPropagation.hpp:133
FractureMechanics::CrackPropagation::assembleMaterialForcesDM
MoFEMErrorCode assembleMaterialForcesDM(DM dm, const int verb=QUIET, const bool debug=false)
create material element instance
Definition: CrackPropagation.cpp:5531
FractureMechanics::CrackPropagation::smootherAlpha
double smootherAlpha
Controls mesh smoothing.
Definition: CrackPropagation.hpp:1200
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FractureMechanics::CrackPropagation::skinOrientation
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > skinOrientation
Definition: CrackPropagation.hpp:1123
FractureMechanics::CrackPropagation::assembleFlambda
boost::shared_ptr< FEMethod > assembleFlambda
assemble F_lambda vector
Definition: CrackPropagation.hpp:1103
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:430
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
FractureMechanics::CrackPropagation::zeroFlambda
boost::shared_ptr< FEMethod > zeroFlambda
assemble F_lambda vector
Definition: CrackPropagation.hpp:1104
FractureMechanics::CrackPropagation::feRhsSimpleContact
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContact
Definition: CrackPropagation.hpp:1270
FractureMechanics::CrackPropagation::getArcLengthDof
MoFEMErrorCode getArcLengthDof()
set pointer to arc-length DOF
Definition: CrackPropagation.cpp:9662
FractureMechanics::CrackPropagation::feCouplingMaterialLhs
boost::shared_ptr< CrackFrontElement > feCouplingMaterialLhs
FE instance to assemble coupling terms.
Definition: CrackPropagation.hpp:1109
AnalyticalDirichletBC
Analytical Dirichlet boundary conditions.
Definition: AnalyticalDirichlet.hpp:18
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_field
virtual MoFEMErrorCode build_field(const std::string field_name, int verb=DEFAULT_VERBOSITY)=0
build field by name
FractureMechanics::CrackPropagation::arcCtx
boost::shared_ptr< ArcLengthCtx > arcCtx
Definition: CrackPropagation.hpp:1215
QUIET
@ QUIET
Definition: definitions.h:208
FractureMechanics::CrackPropagation::projectGriffithForce
MoFEMErrorCode projectGriffithForce(DM dm, Vec f_griffith, Vec f_griffith_proj, const int verb=QUIET, const bool debug=false)
project Griffith forces
Definition: CrackPropagation.cpp:8011
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
FractureMechanics::CrackPropagation::tangentConstrains
boost::shared_ptr< ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain > tangentConstrains
Constrains crack front in tangent directiona.
Definition: CrackPropagation.hpp:1121
FractureMechanics::CrackPropagation::fractionOfFixedNodes
double fractionOfFixedNodes
Definition: CrackPropagation.hpp:1186
FractureMechanics::CrackPropagation::feLhsSimpleContact
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContact
Definition: CrackPropagation.hpp:1272
FractureMechanics::CrackPropagation::createDMs
MoFEMErrorCode createDMs(SmartPetscObj< DM > &dm_elastic, SmartPetscObj< DM > &dm_eigen_elastic, SmartPetscObj< DM > &dm_material, SmartPetscObj< DM > &dm_crack_propagation, SmartPetscObj< DM > &dm_material_forces, SmartPetscObj< DM > &dm_surface_projection, SmartPetscObj< DM > &dm_crack_srf_area, std::vector< int > surface_ids, std::vector< std::string > fe_surf_proj_list)
Crate DMs for all problems and sub problems.
Definition: CrackPropagation.cpp:9571
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
FractureMechanics::CrackPropagation::setCoordsFromField
MoFEMErrorCode setCoordsFromField(const std::string field_name="MESH_NODE_POSITIONS")
set coordinates from field
Definition: CrackPropagation.cpp:9304
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FractureMechanics::CrackPropagation::materialFe
boost::shared_ptr< NonlinearElasticElement > materialFe
Definition: CrackPropagation.hpp:1064
FractureMechanics::CrackPropagation::edgeConstrains
map< int, boost::shared_ptr< EdgeSlidingConstrains > > edgeConstrains
Definition: CrackPropagation.hpp:1129
FractureMechanics::CrackPropagation::addElasticFEInstancesToSnes
MoFEMErrorCode addElasticFEInstancesToSnes(DM dm, Mat m, Vec q, Vec f, boost::shared_ptr< FEMethod > arc_method=boost::shared_ptr< FEMethod >(), boost::shared_ptr< ArcLengthCtx > arc_ctx=nullptr, const int verb=QUIET, const bool debug=false)
Definition: CrackPropagation.cpp:4610
FractureMechanics::CrackPropagation::mwlsApprox
boost::shared_ptr< MWLSApprox > mwlsApprox
Definition: CrackPropagation.hpp:1077
FractureMechanics::CrackPropagation::createProblemDataStructures
MoFEMErrorCode createProblemDataStructures(const std::vector< int > surface_ids, const int verb=QUIET, const bool debug=false)
Construct problem data structures.
Definition: CrackPropagation.cpp:9525
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
FractureMechanics::CrackPropagation::projFrontCtx
boost::shared_ptr< ConstrainMatrixCtx > projFrontCtx
Data structure to project on crack front.
Definition: CrackPropagation.hpp:1075
FractureMechanics::CPSolvers::cP
CrackPropagation & cP
Definition: CPSolvers.hpp:39
convert.init
def init(l)
Definition: convert.py:55
FractureMechanics::CrackPropagation::addSingularity
PetscBool addSingularity
Definition: CrackPropagation.hpp:122
FractureMechanics::CrackPropagation::startStep
int startStep
Definition: CrackPropagation.hpp:168
FractureMechanics::CrackPropagation::getAnalyticalDirichletBc
boost::shared_ptr< AnalyticalDirichletBC::DirichletBC > getAnalyticalDirichletBc()
Definition: CrackPropagation.hpp:923
FractureMechanics::CrackPropagation::feSmootherRhs
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherRhs
Integrate smoothing operators.
Definition: CrackPropagation.hpp:1113
Generated by Doxygen 1.8.17 and hosted at University of Glasgow