57 {
58
61
62 std::vector<boost::weak_ptr<RefEntity>> v_ref_entity;
63 std::vector<boost::weak_ptr<RefElement>> v_ref_element;
64 std::vector<boost::weak_ptr<Field>> v_field;
65 boost::weak_ptr<DofEntity> arc_dof;
66 boost::weak_ptr<MWLSApprox> mwls_approx;
67 std::vector<boost::weak_ptr<FEMethod>> vec_fe;
68 std::vector<boost::weak_ptr<NonlinearElasticElement>> v_elastic_ele_str;
69
70 try {
71
72 PetscBool flg = PETSC_FALSE;
74 if (flg == PETSC_TRUE)
76
79 255, &flg);
80
81 moab::Core mb_instance;
82 moab::Interface &moab = mb_instance;
83
84
85 if (flg) {
86 const char *option;
87 option = "";
89 }
90
91
95
96 auto core_log = logging::core::get();
97 core_log->add_sink(
101
103 auto set_log_attr = [&](
auto c) {
105 };
106 set_log_attr("MEMCHECK");
107 MOFEM_LOG_C(
"MEMCHECK", Sev::inform,
"### Input parameter: -my_file %s",
109
110
112
115 << "File version " << cp.fileVersion.strVersion();
116
117
119
120 {
121 DMType dm_name = "MOFEM";
123 }
124
125 auto add_bits = [&]() {
128 ->addToDatabaseBitRefLevelByType(MBTET,
BitRefLevel().set(),
131 };
133
134 std::vector<int> surface_ids;
136 std::vector<std::string> fe_surf_proj_list;
137 fe_surf_proj_list.push_back("SURFACE");
138
139
140 if (cp.doCutMesh) {
141
143
144
147
150 ->getAllEntitiesNotInDatabase(ents);
152 CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets,
false);
153 for (
auto m : meshsets)
154 CHKERR moab.remove_entities(
m, ents);
155 CHKERR moab.delete_entities(ents);
156
157 auto first_number_string = [](std::string const &str) {
158 std::size_t
const n = str.find_first_of(
"0123456789");
159 if (
n != std::string::npos) {
160 std::size_t
const m = str.find_first_not_of(
"0123456789",
n);
161 return str.substr(
n,
m != std::string::npos ?
m -
n :
m);
162 }
163 return std::string();
164 };
165 std::string str_number =
167
168 if (!str_number.empty())
169 cp.startStep = atoi(str_number.c_str()) + 1;
170
171 string cutting_surface_name =
172 "cutting_surface_" +
173 boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
174
177 cp.crackAccelerationFactor, cutting_surface_name,
VERBOSE,
179 else
182
185
186 } else {
187
188 string cutting_surface_name =
189 "cutting_surface_" +
190 boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
191
194 CHKERR cp.getInterface<
CPMeshCut>()->copySurface(cutting_surface_name);
196 CHKERR moab.get_entities_by_type(0, MBTET, vol,
false);
199 }
200
203
204 } else {
205
207 CHKERR moab.get_entities_by_type(0, MBTET, vol,
false);
213 }
214
215 CHKERR cp.createProblemDataStructures(surface_ids,
QUIET,
217
218 if (cp.propagateCrack == PETSC_TRUE) {
219
220
221 CHKERR cp.setMaterialPositionFromCoords();
222 CHKERR cp.setSpatialPositionFromCoords();
223 if (cp.addSingularity == PETSC_TRUE)
224 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
225
233 CHKERR cp.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
234 dm_crack_propagation, dm_material_forces,
235 dm_surface_projection, dm_crack_srf_area, surface_ids,
236 fe_surf_proj_list);
239
240 CHKERR cp.setSpatialPositionFromCoords();
241 if (cp.addSingularity == PETSC_TRUE) {
242 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
243 }
245 cp.getLoadScale() = cp.getArcCtx()->getFieldData();
246
247 if (cp.solveEigenStressProblem) {
248 MOFEM_LOG(
"MEMCHECK", Sev::verbose) <<
"Solve Eigen ELASTIC problem";
250 dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
251 CHKERR cp.postProcessDM(dm_eigen_elastic, -2,
"ELASTIC",
true);
252 }
253
254 MOFEM_LOG(
"MEMCHECK", Sev::verbose) <<
"Solve ELASTIC problem";
256 dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
257 CHKERR cp.postProcessDM(dm_elastic, -1,
"ELASTIC",
true);
258
259
261 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
262 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
263 false);
264 CHKERR cp.assembleCouplingForcesDM(PETSC_NULL,
VERBOSE,
false);
265
266 CHKERR cp.calculateElasticEnergy(dm_elastic);
268 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
269 surface_ids);
270
271 if (cp.doCutMesh == PETSC_TRUE) {
273 dm_crack_propagation, dm_elastic, dm_eigen_elastic, dm_material,
274 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
276 } else {
278 dm_crack_propagation, dm_elastic, dm_material, dm_material_forces,
279 dm_surface_projection, dm_crack_srf_area, surface_ids);
280 }
281
282 } else {
283
284
285 CHKERR cp.setMaterialPositionFromCoords();
286 CHKERR cp.setSpatialPositionFromCoords();
287 if (cp.addSingularity == PETSC_TRUE)
288 CHKERR cp.setSingularDofs(
"SPATIAL_POSITION");
289
297 CHKERR cp.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
298 dm_crack_propagation, dm_material_forces,
299 dm_surface_projection, dm_crack_srf_area, surface_ids,
300 fe_surf_proj_list);
301
304 "ANALITICAL_DISP"))
305 CHKERR cp.createBcDM(dm_bc,
"BC_PROBLEM",
306 cp.mapBitLevel["spatial_domain"],
308
309 if (cp.solveEigenStressProblem) {
310 MOFEM_LOG(
"MEMCHECK", Sev::verbose) <<
"Solve Eigen ELASTIC problem";
312 dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
313 CHKERR cp.postProcessDM(dm_eigen_elastic, -2,
"ELASTIC",
true);
314 }
315
316
317 MOFEM_LOG(
"MEMCHECK", Sev::verbose) <<
"Solve ELASTIC problem";
319 dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
320
321
323 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
324 if (cp.doCutMesh)
325 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids,
VERBOSE,
326 false);
327
329 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
330 surface_ids);
331 CHKERR cp.calculateElasticEnergy(dm_elastic);
332
333
334 CHKERR cp.postProcessDM(dm_elastic, 0,
"ELASTIC",
true);
335
336 CHKERR cp.savePositionsOnCrackFrontDM(dm_elastic, PETSC_NULL, 2,
false);
337 CHKERR cp.writeCrackFont(cp.mapBitLevel[
"spatial_domain"], 0);
338
339
340 CHKERR cp.tetsingReleaseEnergyCalculation();
341
343 v_ref_entity.reserve(refined_ents_ptr->size());
344 for (auto e : *refined_ents_ptr)
345 v_ref_entity.emplace_back(e);
346
348 v_ref_element.reserve(refined_finite_elements_ptr->size());
349 for (auto e : *refined_finite_elements_ptr)
350 v_ref_element.emplace_back(e);
351
353 v_field.reserve(fields_ptr->size());
354 for (auto e : *fields_ptr)
355 v_field.emplace_back(e);
356
357 arc_dof = cp.arcLengthDof;
358 mwls_approx = cp.mwlsApprox;
359 v_elastic_ele_str.emplace_back(cp.elasticFe);
360 v_elastic_ele_str.emplace_back(cp.materialFe);
361 vec_fe.emplace_back(cp.feLhs);
362 vec_fe.emplace_back(cp.feRhs);
363 vec_fe.emplace_back(cp.feMaterialRhs);
364 vec_fe.emplace_back(cp.feMaterialLhs);
365 vec_fe.emplace_back(cp.feEnergy);
366 vec_fe.emplace_back(cp.feCouplingElasticLhs);
367 vec_fe.emplace_back(cp.feCouplingMaterialLhs);
368 vec_fe.emplace_back(cp.feSmootherRhs);
369 vec_fe.emplace_back(cp.feSmootherLhs);
370 vec_fe.emplace_back(cp.feGriffithForceRhs);
371 vec_fe.emplace_back(cp.feGriffithConstrainsDelta);
372 vec_fe.emplace_back(cp.feGriffithConstrainsRhs);
373 vec_fe.emplace_back(cp.feGriffithConstrainsLhs);
374 vec_fe.emplace_back(cp.feSpringLhsPtr);
375 vec_fe.emplace_back(cp.feSpringRhsPtr);
376 vec_fe.emplace_back(cp.feRhsSimpleContact);
377 vec_fe.emplace_back(cp.feLhsSimpleContact);
378 vec_fe.emplace_back(cp.feMaterialAnaliticalTraction);
379 }
380 }
382
383 if (
auto a = arc_dof.lock()) {
385 <<
"cp.arcLengthDof.use_count() " <<
a.use_count();
387 }
388
389 for (
auto v : vec_fe)
390 if (auto
a =
v.lock()) {
391 MOFEM_LOG(
"MEMCHECK", Sev::error) <<
"fe " <<
a.use_count();
393 }
394
395 for (
auto v : v_elastic_ele_str)
396 if (auto
a =
v.lock()) {
398 <<
"elastic structure " <<
a.use_count();
400 }
401
402 if (
auto a = mwls_approx.lock()) {
404 <<
"cp.mwlsApprox.use_count() " <<
a.use_count();
406 }
407
408 for (auto e : v_ref_entity) {
409 if (
auto a = e.lock()) {
410 MOFEM_LOG(
"MEMCHECK", Sev::error) <<
"v_ref_entity " <<
a.use_count();
412 }
413 }
414
415 for (auto e : v_ref_element) {
416 if (
auto a = e.lock()) {
417 MOFEM_LOG(
"MEMCHECK", Sev::error) <<
"v_ref_element " <<
a.use_count();
419 }
420 }
421
422 for (auto e : v_field) {
423 if (
auto a = e.lock()) {
424 MOFEM_LOG(
"MEMCHECK", Sev::error) <<
"v_field " <<
a.use_count();
426 }
427 }
428
430 return 0;
431}
if(!static_cast< bool >(ifstream(param_file)))
#define MOFEM_LOG_C(channel, severity, format,...)
#define CATCH_ERRORS
Catch errors.
#define BITREFLEVEL_SIZE
max number of refinements
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
virtual const RefElement_multiIndex * get_ref_finite_elements() const =0
Get the ref finite elements object.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
#define MOFEM_LOG_FUNCTION()
Set scope.
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
int getSkinOfTheBodyId() const
const std::string & getCutSurfMeshName() const
virtual MoFEMErrorCode build_field(const std::string field_name, int verb=DEFAULT_VERBOSITY)=0
build field by name
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Interface for managing meshsets containing materials and boundary conditions.
static bool brodcastMeshsets
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.