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

Go to the source code of this file.

Macros

#define BOOST_UBLAS_NDEBUG   1
 

Functions

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

Variables

static char help []
 

Macro Definition Documentation

◆ BOOST_UBLAS_NDEBUG

#define BOOST_UBLAS_NDEBUG   1

Definition at line 31 of file memcheck_leak.cpp.

Function Documentation

◆ main()

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

Definition at line 57 of file memcheck_leak.cpp.

57 {
58
59 const char param_file[] = "param_file.petsc";
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;
73 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-debug", &flg, PETSC_NULL);
74 if (flg == PETSC_TRUE)
76
77 char mesh_file_name[255];
78 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
79 255, &flg);
80
81 moab::Core mb_instance;
82 moab::Interface &moab = mb_instance;
83
84 // Read mesh file
85 if (flg) {
86 const char *option;
87 option = "";
88 CHKERR moab.load_file(mesh_file_name, 0, option);
89 }
90
91 // Create mofem database
93 MoFEM::Core core(moab, PETSC_COMM_WORLD, VERBOSE);
94 MoFEM::Interface &m_field = core;
95
96 auto core_log = logging::core::get();
97 core_log->add_sink(
99 LogManager::setLog("MEMCHECK");
100 MOFEM_LOG_TAG("MEMCHECK", "memcheck");
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 // Create data structure for crack propagation
111 CrackPropagation cp(m_field, 3, 2);
112
113 CHKERR MoFEM::Interface::getFileVersion(moab, cp.fileVersion);
114 MOFEM_LOG("MEMCHECK", Sev::inform)
115 << "File version " << cp.fileVersion.strVersion();
116
117 // Configure meshes
118 CHKERR cp.getOptions();
119 // Register MOFEM discrete data managers
120 {
121 DMType dm_name = "MOFEM";
122 CHKERR DMRegister_MoFEM(dm_name);
123 }
124
125 auto add_bits = [&]() {
128 ->addToDatabaseBitRefLevelByType(MBTET, BitRefLevel().set(),
129 BitRefLevel().set());
131 };
132 CHKERR add_bits();
133
134 std::vector<int> surface_ids;
135 surface_ids.push_back(cp.getInterface<CPMeshCut>()->getSkinOfTheBodyId());
136 std::vector<std::string> fe_surf_proj_list;
137 fe_surf_proj_list.push_back("SURFACE");
138
139 // create bit levels, fields and finite elements structures in database
140 if (cp.doCutMesh) {
141
142 if (string(mesh_file_name).compare(0, 8, "restart_") == 0) {
143
144 // restart code from file
145 CHKERR m_field.build_field("LAMBDA_ARC_LENGTH");
146 CHKERR m_field.build_finite_elements("ARC_LENGTH");
147
148 Range ents;
150 ->getAllEntitiesNotInDatabase(ents);
151 Range meshsets;
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 =
166 first_number_string(std::string(mesh_file_name));
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
175 if (m_field.get_comm_rank() == 0)
176 CHKERR cp.getInterface<CPMeshCut>()->rebuildCrackSurface(
177 cp.crackAccelerationFactor, cutting_surface_name, VERBOSE,
179 else
180 CHKERR cp.getInterface<CPMeshCut>()->rebuildCrackSurface(
181 cp.crackAccelerationFactor, "", NOISY, CrackPropagation::debug);
182
183 CHKERR cp.getInterface<CPMeshCut>()->refineMesh(
184 nullptr, false, QUIET, CrackPropagation::debug);
185
186 } else {
187
188 string cutting_surface_name =
189 "cutting_surface_" +
190 boost::lexical_cast<std::string>(cp.startStep) + ".vtk";
191
192 if (!cp.getInterface<CPMeshCut>()->getCutSurfMeshName().empty())
193 CHKERR cp.getInterface<CPMeshCut>()->setCutSurfaceFromFile();
194 CHKERR cp.getInterface<CPMeshCut>()->copySurface(cutting_surface_name);
195 Range vol;
196 CHKERR moab.get_entities_by_type(0, MBTET, vol, false);
197 CHKERR cp.getInterface<CPMeshCut>()->refineMesh(
198 &vol, true, QUIET, CrackPropagation::debug);
199 }
200
201 CHKERR cp.getInterface<CPMeshCut>()->cutRefineAndSplit(
203
204 } else {
205
206 Range vol;
207 CHKERR moab.get_entities_by_type(0, MBTET, vol, false);
208 BitRefLevel bit0 = cp.mapBitLevel["mesh_cut"];
209 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
210 vol, bit0.set(BITREFLEVEL_SIZE - 1), true, CrackPropagation::debug);
211 CHKERR cp.getInterface<CPMeshCut>()->refineAndSplit(
213 }
214
215 CHKERR cp.createProblemDataStructures(surface_ids, QUIET,
217
218 if (cp.propagateCrack == PETSC_TRUE) {
219
220 // set inital coordinates
221 CHKERR cp.setMaterialPositionFromCoords();
222 CHKERR cp.setSpatialPositionFromCoords();
223 if (cp.addSingularity == PETSC_TRUE)
224 CHKERR cp.setSingularDofs("SPATIAL_POSITION");
225
226 SmartPetscObj<DM> dm_elastic;
227 SmartPetscObj<DM> dm_eigen_elastic;
228 SmartPetscObj<DM> dm_material;
229 SmartPetscObj<DM> dm_crack_propagation;
230 SmartPetscObj<DM> dm_material_forces;
231 SmartPetscObj<DM> dm_surface_projection;
232 SmartPetscObj<DM> dm_crack_srf_area;
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 }
244 if (string(mesh_file_name).compare(0, 8, "restart_") == 0)
245 cp.getLoadScale() = cp.getArcCtx()->getFieldData();
246
247 if (cp.solveEigenStressProblem) {
248 MOFEM_LOG("MEMCHECK", Sev::verbose) << "Solve Eigen ELASTIC problem";
249 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
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";
255 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
256 dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
257 CHKERR cp.postProcessDM(dm_elastic, -1, "ELASTIC", true);
258
259 // set finite element instances and user data operators on instances
260 CHKERR cp.assembleElasticDM(VERBOSE, false);
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);
267 CHKERR cp.getInterface<CPSolvers>()->calculateGc(
268 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
269 surface_ids);
270
271 if (cp.doCutMesh == PETSC_TRUE) {
272 CHKERR cp.getInterface<CPSolvers>()->solveCutMesh(
273 dm_crack_propagation, dm_elastic, dm_eigen_elastic, dm_material,
274 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
275 surface_ids, CrackPropagation::debug);
276 } else {
277 CHKERR cp.getInterface<CPSolvers>()->solvePropagation(
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 // set inital coordinates
285 CHKERR cp.setMaterialPositionFromCoords();
286 CHKERR cp.setSpatialPositionFromCoords();
287 if (cp.addSingularity == PETSC_TRUE)
288 CHKERR cp.setSingularDofs("SPATIAL_POSITION");
289
290 SmartPetscObj<DM> dm_elastic;
291 SmartPetscObj<DM> dm_eigen_elastic;
292 SmartPetscObj<DM> dm_material;
293 SmartPetscObj<DM> dm_crack_propagation;
294 SmartPetscObj<DM> dm_material_forces;
295 SmartPetscObj<DM> dm_surface_projection;
296 SmartPetscObj<DM> dm_crack_srf_area;
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
302 SmartPetscObj<DM> dm_bc;
304 "ANALITICAL_DISP"))
305 CHKERR cp.createBcDM(dm_bc, "BC_PROBLEM",
306 cp.mapBitLevel["spatial_domain"],
307 BitRefLevel().set());
308
309 if (cp.solveEigenStressProblem) {
310 MOFEM_LOG("MEMCHECK", Sev::verbose) << "Solve Eigen ELASTIC problem";
311 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
312 dm_eigen_elastic, cp.getEigenArcCtx()->x0, 1);
313 CHKERR cp.postProcessDM(dm_eigen_elastic, -2, "ELASTIC", true);
314 }
315
316 // set finite element instances and user data operators on instances
317 MOFEM_LOG("MEMCHECK", Sev::verbose) << "Solve ELASTIC problem";
318 CHKERR cp.getInterface<CPSolvers>()->solveElastic(
319 dm_elastic, cp.getArcCtx()->x0, cp.getLoadScale());
320
321 // set finite element instances and user data operators on instances
322 CHKERR cp.assembleElasticDM(VERBOSE, false);
323 CHKERR cp.assembleMaterialForcesDM(PETSC_NULL);
324 if (cp.doCutMesh)
325 CHKERR cp.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE,
326 false);
327
328 CHKERR cp.getInterface<CPSolvers>()->calculateGc(
329 dm_material_forces, dm_surface_projection, dm_crack_srf_area,
330 surface_ids);
331 CHKERR cp.calculateElasticEnergy(dm_elastic);
332
333 // post-processing
334 CHKERR cp.postProcessDM(dm_elastic, 0, "ELASTIC", true);
335 // set vector values from field data
336 CHKERR cp.savePositionsOnCrackFrontDM(dm_elastic, PETSC_NULL, 2, false);
337 CHKERR cp.writeCrackFont(cp.mapBitLevel["spatial_domain"], 0);
338
339 // Run tests
340 CHKERR cp.tetsingReleaseEnergyCalculation();
341
342 auto refined_ents_ptr = m_field.get_ref_ents();
343 v_ref_entity.reserve(refined_ents_ptr->size());
344 for (auto e : *refined_ents_ptr)
345 v_ref_entity.emplace_back(e);
346
347 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
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
352 auto fields_ptr = m_field.get_fields();
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()) {
384 MOFEM_LOG("MEMCHECK", Sev::error)
385 << "cp.arcLengthDof.use_count() " << a.use_count();
386 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
387 }
388
389 for (auto v : vec_fe)
390 if (auto a = v.lock()) {
391 MOFEM_LOG("MEMCHECK", Sev::error) << "fe " << a.use_count();
392 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
393 }
394
395 for (auto v : v_elastic_ele_str)
396 if (auto a = v.lock()) {
397 MOFEM_LOG("MEMCHECK", Sev::error)
398 << "elastic structure " << a.use_count();
399 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
400 }
401
402 if (auto a = mwls_approx.lock()) {
403 MOFEM_LOG("MEMCHECK", Sev::error)
404 << "cp.mwlsApprox.use_count() " << a.use_count();
405 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
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();
411 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
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();
418 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
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();
425 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "cyclic shared pointer");
426 }
427 }
428
430 return 0;
431}
std::string param_file
if(!static_cast< bool >(ifstream(param_file)))
#define MOFEM_LOG_C(channel, severity, format,...)
constexpr double a
@ QUIET
@ NOISY
@ VERBOSE
#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 ...
@ MOFEM_INVALID_DATA
Definition definitions.h:36
#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.
Definition DMMoFEM.cpp:47
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)
static char help[]
char mesh_file_name[255]
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
const std::string & getCutSurfMeshName() const
Managing BitRefLevels.
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
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Interface for managing meshsets containing materials and boundary conditions.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.

Variable Documentation

◆ help

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

Definition at line 44 of file memcheck_leak.cpp.