v0.13.1
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
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
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

CrackPropagationcP
 

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}
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
CrackPropagation & cP
Definition: CPSolvers.hpp:39

◆ ~CPSolvers()

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

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;
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
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
244}
@ VERBOSE
Definition: definitions.h:209
@ ROW
Definition: definitions.h:123
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
const FTensor::Tensor2< T, Dim, Dim > Vec
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
MoFEMErrorCode calculateGriffithForce(DM dm, const double gc, Vec f_griffith, const int verb=QUIET, const bool debug=false)
calculate Griffith (driving) force
MoFEMErrorCode projectGriffithForce(DM dm, Vec f_griffith, Vec f_griffith_proj, const int verb=QUIET, const bool debug=false)
project Griffith forces
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
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
MoFEMErrorCode calculateSurfaceProjectionMatrix(DM dm_front, DM dm_project, const std::vector< int > &ids, const int verb=QUIET, const bool debug=false)
assemble projection matrices
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)
Deprecated interface functions.
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23

◆ getOptions()

MoFEMErrorCode FractureMechanics::CPSolvers::getOptions ( )

Definition at line 80 of file CPSolvers.cpp.

80 {
82
83 // CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Solvers options", "none");
84 // ierr = PetscOptionsEnd();
85 // CHKERRG(ierr);
86
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}
CPSolvers(CrackPropagation &cp)
Definition: CPSolvers.cpp:56

◆ 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 {
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());
873
874 auto save_position_on_mesh_tags_from_coords =
875 [this](const std::string tag_name) {
876 Tag th;
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;
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());
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();
1008 cP.feSpringLhsPtr.reset();
1009 cP.feSpringRhsPtr.reset();
1010 cP.feRhsSimpleContact.reset();
1011 cP.feLhsSimpleContact.reset();
1013 cP.bothSidesConstrains.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
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 = [&]() {
1127 ->addToDatabaseBitRefLevelByType(MBTET, BitRefLevel().set(),
1128 BitRefLevel().set());
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 ||
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);
1164 cP.mField.get_moab(), moab_read, 0, ents, false, true);
1165 } else {
1166 Range ents;
1168 cP.mField.get_moab(), cP.mField.get_moab(), 0, ents, false, true);
1169 }
1170 }
1171
1173 cP.fileVersion);
1174
1176 CHKERR cP.getInterface<CPMeshCut>()->getInterfacesPtr();
1177 CHKERR add_bits_tets();
1178
1179 Range ents;
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(
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
1209
1210 // set inital coordinates
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 = [&]() {
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);
1233 };
1234
1235 CHKERR set_up_arc_length();
1236
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;
1245 MOFEM_LOG("CPSolverWorld", Sev::verbose)
1246 << "Solve Eigen ELASTIC problem";
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";
1256 dm_elastic, cP.getArcCtx()->x0, load_factor);
1257
1258 // set finite element instances and user data operators on instances
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);
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
1278}
#define SolveFunctionBegin
Definition: CPSolvers.cpp:43
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:338
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
constexpr double a
@ QUIET
Definition: definitions.h:208
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
static const bool debug
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSetSnesCtx(DM dm, boost::shared_ptr< MoFEM::SnesCtx > snes_ctx)
Set MoFEM::SnesCtx data structure.
Definition: DMMMoFEM.cpp:1058
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
const double v
phase velocity of light in medium (cm/ns)
MoFEMErrorCode broadcast_entities(moab::Interface &moab, moab::Interface &moab_tmp, const int from_proc, Range &entities, const bool adjacencies, const bool tags)
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
CoreTmp< 0 > Core
Definition: Core.hpp:1086
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
Store variables for ArcLength analysis.
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
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
MoFEMErrorCode solveElastic(DM dm_elastic, Vec init, const double load_factor)
Solve elastic problem.
Definition: CPSolvers.cpp:90
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringRhsPtr
boost::shared_ptr< BothSurfaceConstrains > bothSidesConstrains
boost::shared_ptr< CrackFrontElement > feMaterialRhs
Integrate material stresses, assemble vector.
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringLhsPtr
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > contactOrientation
map< int, boost::shared_ptr< SurfaceSlidingConstrains > > surfaceConstrain
boost::shared_ptr< CrackFrontElement > feRhs
Integrate elastic FE.
std::map< std::string, BitRefLevel > mapBitLevel
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsLhs
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherLhs
Integrate smoothing operators.
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContact
boost::shared_ptr< ConstrainMatrixCtx > projFrontCtx
Data structure to project on crack front.
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsDelta > feGriffithConstrainsDelta
boost::shared_ptr< NonlinearElasticElement > elasticFe
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherRhs
Integrate smoothing operators.
boost::shared_ptr< DirichletFixFieldAtEntitiesBc > fixMaterialEnts
boost::shared_ptr< FEMethod > zeroFlambda
assemble F_lambda vector
boost::shared_ptr< NonlinearElasticElement > materialFe
MoFEMErrorCode setSingularDofs(const string field_name, const int verb=QUIET)
set singular dofs (on edges adjacent to crack front) from geometry
boost::shared_ptr< NeumannForcesSurface::MyTriangleFE > feMaterialAnaliticalTraction
Surface elment to calculate tractions in material space.
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsRhs
MoFEMErrorCode assembleMaterialForcesDM(DM dm, const int verb=QUIET, const bool debug=false)
create material element instance
boost::shared_ptr< Smoother > smootherFe
MoFEMErrorCode setSpatialPositionFromCoords()
set spatial field from nodes
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > skinOrientation
MoFEMErrorCode setMaterialPositionFromCoords()
set material field from nodes
MoFEMErrorCode calculateElasticEnergy(DM dm, const std::string msg="")
assemble elastic part of matrix, by running elastic finite element instance
boost::shared_ptr< VolumeLengthQuality< adouble > > volumeLengthAdouble
MoFEMErrorCode assembleSmootherForcesDM(DM dm, const std::vector< int > ids, const int verb=QUIET, const bool debug=false)
create smoothing element instance
boost::shared_ptr< VolumeLengthQuality< double > > volumeLengthDouble
boost::shared_ptr< ConstrainMatrixCtx > projSurfaceCtx
Data structure to project on the body surface.
boost::shared_ptr< ArcLengthCtx > arcCtx
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > crackOrientation
map< int, boost::shared_ptr< EdgeSlidingConstrains > > edgeConstrains
PetscBool solveEigenStressProblem
Solve eigen problem.
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContact
boost::shared_ptr< CrackFrontElement > feCouplingMaterialLhs
FE instance to assemble coupling terms.
boost::shared_ptr< CrackFrontElement > feCouplingElasticLhs
FE instance to assemble coupling terms.
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
boost::shared_ptr< GriffithForceElement > griffithForceElement
MoFEMErrorCode assembleCouplingForcesDM(DM dm, const int verb=QUIET, const bool debug=false)
assemble coupling element instances
boost::shared_ptr< MWLSApprox > mwlsApprox
boost::shared_ptr< CrackFrontElement > feLhs
Integrate elastic FE.
boost::shared_ptr< BothSurfaceConstrains > bothSidesContactConstrains
boost::shared_ptr< FEMethod > assembleFlambda
assemble F_lambda vector
boost::shared_ptr< BothSurfaceConstrains > closeCrackConstrains
boost::shared_ptr< ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain > tangentConstrains
Constrains crack front in tangent directiona.
boost::shared_ptr< WrapMPIComm > moabCommWorld
MoFEMErrorCode createProblemDataStructures(const std::vector< int > surface_ids, const int verb=QUIET, const bool debug=false)
Construct problem data structures.
boost::shared_ptr< ArcLengthCtx > & getEigenArcCtx()
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.
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceRhs
boost::shared_ptr< CrackFrontElement > feEnergy
Integrate energy.
boost::shared_ptr< CrackFrontElement > feMaterialLhs
Integrate material stresses, assemble matrix.
boost::shared_ptr< ArcLengthCtx > & getArcCtx()
Managing BitRefLevels.
virtual int get_comm_size() const =0
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 MoFEMErrorCode set_moab_interface(moab::Interface &new_moab, int verb=VERBOSE)=0
Set the moab interface object.
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode clear_database(int verb=DEFAULT_VERBOSITY)=0
Clear database.
virtual int get_comm_rank() const =0
static MoFEMErrorCode setFileVersion(moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
Get database major version.

◆ solveElastic()

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

Solve elastic problem.

Parameters
dm_elastic
Returns
MoFEMErrorCode

Definition at line 90 of file CPSolvers.cpp.

91 {
92 MoFEM::Interface &m_field = cP.mField;
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) {
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
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",
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
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")
177 if (problem_ptr->getName() == "ELASTIC")
179
180 // + and ops which will close crack
181 CHKERR solve_elastic_problem(init);
182
184}
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:953
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
virtual bool check_problem(const std::string name)=0
check if problem exist
auto createSNES(MPI_Comm comm)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:987
def init(l)
Definition: convert.py:55
Analytical Dirichlet boundary conditions.
MoFEMErrorCode solveElasticDM(DM dm, SNES snes, Mat m, Vec q, Vec f, bool snes_set_up, Mat *shell_m)
solve elastic problem
MoFEMErrorCode getArcLengthDof()
set pointer to arc-length DOF
std::string mwlsEigenStressTagName
Name of tag with eigen stresses.
std::string mwlsStressTagName
Name of tag with internal stresses.
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)
boost::shared_ptr< AnalyticalDirichletBC::DirichletBC > getAnalyticalDirichletBc()

◆ 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;
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
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
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
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 }
321 }
322
323 MoFEMErrorCode operator()(DM dm, SNES snes,
324 boost::shared_ptr<ArcLengthCtx> arc_ctx,
325 double &s) const {
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 }
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
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
409 }
410
411 MoFEMErrorCode operator()(DM dm, SNES snes, Vec q) const {
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;
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);
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
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
458 }
459
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
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 = [&]() {
513
514 const double arc_s = cP.arcS;
515
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);
531 };
532
533 CHKERR solve_problem();
534
535 auto if_not_converged = [&](Vec q) {
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 }
545 };
546 CHKERR if_not_converged(q);
547
548 auto scatter_forward_vec = [](Vec q) {
550 CHKERR VecGhostUpdateBegin(q, INSERT_VALUES, SCATTER_FORWARD);
551 CHKERR VecGhostUpdateEnd(q, INSERT_VALUES, SCATTER_FORWARD);
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,
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) {
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,
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
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) {
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);
677 };
678
679 auto add_fe = [&](auto fe_name, auto &tags_list) {
681 auto arc_length_finite_element =
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);
700 };
701
702 auto arc_field_meshset =
703 cP.mField.get_field_meshset("LAMBDA_ARC_LENGTH");
704 auto arc_length_finite_element =
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(
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
840 };
841
842 if (cP.mField.get_comm_rank() == 0)
843 CHKERR save_restart_file();
844 }
845 }
846
848}
static Index< 'p', 3 > p
@ COL
Definition: definitions.h:123
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
virtual EntityHandle get_finite_element_meshset(const std::string &name) const =0
virtual Field * get_field_structure(const std::string &name)=0
get field structure
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.
auto bit
set bit
auto f
Definition: HenckyOps.hpp:5
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary meshset.
Definition: Templates.hpp:1465
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
MoFEMErrorCode getOptions()
Definition: CPSolvers.cpp:80
MoFEMErrorCode zeroLambdaFields()
Zero fields with lagrange multipliers.
MoFEMErrorCode solvePropagationDM(DM dm, DM dm_elastic, SNES snes, Mat m, Vec q, Vec f)
solve crack propagation problem
MoFEMErrorCode postProcessDM(DM dm, const int step, const std::string fe_name, const bool approx_internal_stress)
post-process results for elastic solution
MoFEMErrorCode savePositionsOnCrackFrontDM(DM dm, Vec q, const int verb=QUIET, const bool debug=false)
map< EntityHandle, double > mapSmoothingForceFactor
boost::shared_ptr< FEMethod > getFrontArcLengthControl(boost::shared_ptr< ArcLengthCtx > arc_ctx)
MoFEMErrorCode updateMaterialFixedNode(const bool fix_front, const bool fix_small_g, const bool debug=false)
Update fixed nodes.
double smootherAlpha
Controls mesh smoothing.
MoFEMErrorCode setCoordsFromField(const std::string field_name="MESH_NODE_POSITIONS")
set coordinates from field
MoFEMErrorCode calculateSmoothingForceFactor(const int verb=QUIET, const bool debug=true)
MoFEMErrorCode writeCrackFont(const BitRefLevel &bit, const int step)
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
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
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
EntityHandle getMeshset() const
Get field meshset.
Interface for managing meshsets containing materials and boundary conditions.
CubitMeshSet_multiIndex & getMeshsetsMultindex()

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: