v0.13.1
CrackPropagation.hpp
Go to the documentation of this file.
1/** \file CrackPropagation.hpp
2 \brief Main class for crack propagation
3*/
4
5/* This file is part of MoFEM.
6 * MoFEM is free software: you can redistribute it and/or modify it under
7 * the terms of the GNU Lesser General Public License as published by the
8 * Free Software Foundation, either version 3 of the License, or (at your
9 * option) any later version.
10 *
11 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 * License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18
19#ifndef __CRACK_PROPAGATION_HPP__
20#define __CRACK_PROPAGATION_HPP__
21
22namespace FractureMechanics {
23
25 boost::shared_ptr<WrapMPIComm> moab_comm_wrap);
26
28 moab::Interface &moab_tmp,
29 const int from_proc, Range &entities,
30 const bool adjacencies, const bool tags);
31
32/**
33 * @brief Tapes numbers used by ADOL-C.
34 *
35 */
48};
49
50/**
51 * @brief Names of tangent matrices tests
52 *
53 */
59};
60/**
61 * Crack propagation data and methods
62 * \nosubgrouping
63 *
64 * - Functions strarting with build, for example buildElastic,
65 * add and build fields.
66 *
67 * - Functions starting with declare, for example declareElasticFE,
68 * add and build finite elements.
69 *
70 * - Functions starting with assemble, for example assembleElasticDM,
71 * create finite element instances and operators.
72 *
73 * - Functions starting with create, for example createElasticDM, create
74 * problems.
75 *
76 */
78
79 static bool debug;
80
81 static bool parallelReadAndBroadcast; ///< That is unstable development, for
82 ///< some meshses (propably generated by
83 ///< cubit) this is not working. Error
84 ///< can be attributed to bug in MOAB.
85
86 Version fileVersion;
87
88 MoFEMErrorCode getInterfaceVersion(Version &version) const {
90 version = Version(FM_VERSION_MAJOR, FM_VERSION_MINOR, FM_VERSION_BUILD);
92 }
93
94 /**
95 * \brief Getting interface of core database
96 * @param uuid unique ID of interface that can be either for
97 * CrackPropagation, CPSolvers or CPMeshCut interface
98 * @param iface returned pointer to interface
99 * @return error code
100 */
101 MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
102 UnknownInterface **iface) const;
103
104 boost::shared_ptr<WrapMPIComm> moabCommWorld;
105
110
111 double gC; ///< Griffith energy
112 double griffithE; ///< Griffith stability parameter
113 double griffithR; ///< Griffith regularisation parameter
114
115 double betaGc; ///< heterogeneous Griffith energy exponent
116 PetscBool isGcHeterogeneous; ///< flag for heterogeneous gc
117
118 PetscBool isPartitioned;
119 PetscBool propagateCrack; ///< If true crack propagation is calculated
122 PetscBool addSingularity;
124
125 PetscBool isPressureAle; ///< If true surface pressure is considered in ALE
126 PetscBool areSpringsAle; ///< If true surface spring is considered in ALE
127
128 std::string mwlsApproxFile; ///< Name of file with internal stresses
129 std::string mwlsStressTagName; ///< Name of tag with internal stresses
130 std::string mwlsEigenStressTagName; ///< Name of tag with eigen stresses
131 std::string mwlsRhoTagName; ///< Name of tag with density
132 int residualStressBlock; ///< Block on which residual stress is applied
133 int densityMapBlock; ///< Block on which density is applied (Note: if used
134 ///< with internal stress, it has to be the same ID!)
135 std::string configFile;
137
138 double
139 partitioningWeightPower; ///< Weight controling load balance. Elements at
140 ///< the crack front are higher order, also have
141 ///< more fildes, associated with ALE
142 ///< formulation. Higher number put more weight
143 ///< on those elements while partitioning.
144
145 PetscBool resetMWLSCoeffsEveryPropagationStep; ///< If true, MWLS coefficients
146 ///< are recalulated every
147 ///< propagaton step.
148
150
152
153 );
154
155 PetscBool
156 addAnalyticalInternalStressOperators; ///< Add analytical internal stress
157 ///< operators for MWLS test
158
159 PetscBool solveEigenStressProblem; ///< Solve eigen problem
160 PetscBool useEigenPositionsSimpleContact; ///< Use eigen positions for
161 ///< matching meshes contact
162
163 int postProcLevel; ///< level of postprocessing (amount of output files)
165
168
169 /**
170 * Map associating string literals to bit refinement levels
171 * The following convention is currently used:
172 * "mesh_cut" : level 0
173 * "spatial_domain" : level 1
174 * "material_domain" : level 2
175 * */
176 std::map<std::string, BitRefLevel> mapBitLevel;
177
178 /**
179 * \brief Constructor of Crack Propagation, prints fracture module version and
180 registers the three interfaces
181 * used to solve fracture problem, i.e. CrackPropagation, CPSolvers and
182 * CPMeshCut
183 * @param m_field interface to MoAB database
184 * @param approx_order order of approximation space, default is 2
185 * @param geometry_order order of geometry, default is 1 (can be up to 2)
186 */
187 CrackPropagation(MoFEM::Interface &m_field, const int approx_order = 2,
188 const int geometry_order = 1);
189
190 virtual ~CrackPropagation();
191
192 /**
193 * \brief Get options form command line
194 * @return error code
195 */
197
198 /**
199 * \brief This is run with ctest
200 * @return error code
201 */
203
204 inline int &getNbLoadSteps() { return nbLoadSteps; }
205 inline int getNbLoadSteps() const { return nbLoadSteps; }
206 inline double &getLoadScale() { return loadScale; }
207 inline double getLoadScale() const { return loadScale; }
208 inline int &getNbCutSteps() { return nbCutSteps; }
209 inline int getNbCutSteps() const { return nbCutSteps; }
210
211 /// \brief read mesh file
213
214 /**
215 * @brief Save entities on ech processor.
216 *
217 * Usually used for debugging, to check if meshes are consistent on all
218 * processors.
219 *
220 * @param prefix file name prefix
221 * @param ents entities to save
222 * @return MoFEMErrorCode error code
223 */
224 MoFEMErrorCode saveEachPart(const std::string prefix, const Range &ents);
225
226 /** \name Resolve shared entities and partition mesh */
227
228 /**@{*/
229
230 /// \brief partotion mesh
232 int verb = QUIET, const bool debug = false);
233
234 /// \brief resolve shared entities
235 MoFEMErrorCode resolveShared(const Range &tets, Range &proc_ents,
236 const int verb = QUIET,
237 const bool debug = false);
238
239 /// \brief resole shared entities by bit level
241 const int verb = QUIET,
242 const bool debug = false);
243
244 /**@}*/
245
246 /** \name Set bit level for material problem */
247
248 /**@{*/
249
250 /**
251 * @brief Set bit ref level for entities adjacent to crack front
252 *
253 * @param from_bit
254 * @param bit
255 * @param nb_levels adjacency bridge levels
256 * @return MoFEMErrorCode
257 */
259 const int nb_levels = 2,
260 const bool debug = false);
261
262 /**@}*/
263
264 /** \name Build problem */
265
266 /**@{*/
267
268 /**
269 * \brief Build problem fields
270 * @param bit1 bit level for computational domain
271 * @param mask1 bit mask
272 * @param bit2 bit level for elements near crack front
273 * @param mask2 bit mask 2
274 * @param verb verbosity level
275 * @param debug
276 * @return error code
277 *
278 */
280 const BitRefLevel &mask1,
281 const BitRefLevel &bit2,
282 const int verb = QUIET,
283 const bool debug = false);
284
285 /**
286 * \brief Build problem finite elements
287 * @param bit1 bit level for computational domain
288 * @param bit2 bit level for elements near crack front
289 * @param surface_ids meshsets Ids on body skin
290 * @param verb verbosity level
291 * @param debug
292 * @return error code
293 *
294 */
296 std::vector<int> surface_ids,
297 const int verb = QUIET,
298 const bool debug = false);
299
300 /**
301 * \brief Construct problem data structures
302 * @param surface_ids surface_ids meshsets Ids on body skin
303 * @param verb verbosity level
304 * @param debug
305 * @return error code
306 */
307 MoFEMErrorCode createProblemDataStructures(const std::vector<int> surface_ids,
308 const int verb = QUIET,
309 const bool debug = false);
310
311 /**
312 * \brief Crate DMs for all problems and sub problems
313 * @param dm_elastic elastic DM
314 * @param dm_eiegn_elastci eigen_elastic DM to solve proble wih eigen
315 * stresses
316 * @param dm_material material DM
317 * @param dm_crack_propagation composition of elastic DM and material DM
318 * @param dm_material_forces used to calculate material forces, sub-dm of
319 * dm_crack_propagation
320 * @param dm_surface_projection dm to project material forces on surfaces,
321 * sub-dm of dm_crack_propagation
322 * @param dm_crack_srf_area dm to project material forces on crack
323 * surfaces, sub-dm of dm_crack_propagation
324 * @param surface_ids IDs of surface meshsets
325 * @param fe_surf_proj_list name of elements on surface elements
326 * @return error code
327 */
328 MoFEMErrorCode createDMs(SmartPetscObj<DM> &dm_elastic,
329 SmartPetscObj<DM> &dm_eigen_elastic,
330 SmartPetscObj<DM> &dm_material,
331 SmartPetscObj<DM> &dm_crack_propagation,
332 SmartPetscObj<DM> &dm_material_forces,
333 SmartPetscObj<DM> &dm_surface_projection,
334 SmartPetscObj<DM> &dm_crack_srf_area,
335 std::vector<int> surface_ids,
336 std::vector<std::string> fe_surf_proj_list);
337
338 MoFEMErrorCode deleteEntities(const int verb = QUIET,
339 const bool debug = false);
340
341 /**@}*/
342
343 /** \name Build fields */
344
345 /**@{*/
346
347 /**
348 * \brief Declate fields for elastic analysis
349 * @param bit bit refinement level
350 * @return error code
351 */
353 const BitRefLevel bit, const BitRefLevel mask = BitRefLevel().set(),
354 const bool proc_only = true, const bool build_fields = true,
355 const int verb = QUIET, const bool debug = false);
356
357 /**
358 * \brief Declate field for arc-length
359 * @param bit bit refinement level
360 * @return error code
361 */
363 const bool build_fields = true,
364 const int verb = QUIET);
365
366 /// \brief build fields with Lagrange multipliers to constrain surfaces
368 const bool proc_only = true,
369 const bool build_fields = true,
370 const int verb = QUIET,
371 const bool debug = false);
372
373 /// \brief build fields with Lagrange multipliers to constrain edges
375 const bool proc_only = true,
376 const bool build_fields = true,
377 const int verb = QUIET,
378 const bool debug = false);
379
380 /// \brief declare crack surface files
382 const bool proc_only = true,
383 const bool build_fields = true,
384 const int verb = QUIET,
385 const bool debug = false);
386
387 /// \brief declare crack surface files
389 const bool build_fields = true,
390 const int verb = QUIET,
391 const bool debug = false);
392
393 /**
394 * \brief Lagrange multipliers field which constrains material displacements
395 *
396 * Nodes on both sides are constrainded to have the same material
397 * displacements
398 *
399 */
401 const BitRefLevel bit_material,
402 const bool proc_only = false,
403 const bool build_fields = true,
404 const int verb = QUIET,
405 const bool debug = false);
406
407 /** \brief Zero fields with lagrange multipliers
408 */
410
411 /**@}*/
412
413 /** \name Build finite elements */
414
415 /**@{*/
416
417 /// \brief declare elastic finite elements
419 declareElasticFE(const BitRefLevel bit1, const BitRefLevel mask1,
420 const BitRefLevel bit2, const BitRefLevel mask2,
421 const bool add_forces = true, const bool proc_only = true,
422 const int verb = QUIET);
423
424 /// \brief create arc-length element entity and declare elemets
426 const int verb = QUIET);
427
430 const BitRefLevel mask = BitRefLevel().set(),
431 const bool proc_only = true);
432
433 /** \brief Declare FE for pressure BC in ALE formulation (in material domain)
434 * @param bit material domain bit ref level
435 * @param mask bit ref level mask for the problem
436 * @return error code
437 */
440 const BitRefLevel mask = BitRefLevel().set(),
441 const bool proc_only = true);
442
443 /** \brief Declare FE for spring BC in ALE formulation (in material domain)
444 * @param bit material domain bit ref level
445 * @param mask bit ref level mask for the problem
446 * @return error code
447 */
450 const BitRefLevel mask = BitRefLevel().set(),
451 const bool proc_only = true);
452
453 /** \brief Declare FE for pressure BC in ALE formulation (in material domain)
454 * @param bit material domain bit ref level
455 * @param mask bit ref level mask for the problem
456 * @return error code
457 */
460 const BitRefLevel mask = BitRefLevel().set(),
461 const bool proc_only = true);
462
463 /// \brief declare material finite elements
465 const BitRefLevel mask = BitRefLevel().set(),
466 const bool proc_only = true,
467 const bool verb = QUIET);
468
470 const BitRefLevel mask = BitRefLevel().set(),
471 const bool proc_only = true,
472 const bool verb = QUIET);
473
475 declareBothSidesFE(const BitRefLevel bit_spatial,
476 const BitRefLevel bit_material,
477 const BitRefLevel mask = BitRefLevel().set(),
478 const bool proc_only = true, const bool verb = QUIET);
479
480 /// \brief declare mesh smoothing finite elements
483 const BitRefLevel mask = BitRefLevel().set(),
484 const bool proc_only = true, const bool verb = QUIET);
485
486 /// \brief declare surface sliding elements
487 MoFEMErrorCode declareSurfaceFE(std::string fe_name, const BitRefLevel bit,
488 const BitRefLevel mask,
489 const std::vector<int> &ids,
490 const bool proc_only = true,
491 const int verb = QUIET,
492 const bool debug = false);
493
494 // \brief declare edge sliding elements
495 MoFEMErrorCode declareEdgeFE(std::string fe_name, const BitRefLevel bit,
496 const BitRefLevel mask,
497 const bool proc_only = true,
498 const int verb = QUIET,
499 const bool debug = false);
500
502 declareCrackSurfaceFE(std::string fe_name, const BitRefLevel bit,
503 const BitRefLevel mask, const bool proc_only = true,
504 const int verb = QUIET, const bool debug = false);
505
506 /**@}*/
507
508 /** \name Create DMs & problems */
509
510 /**@{*/
511
512 /** \brief Create DM by composition of elastic DM and material DM
513 * @param dm composition of elastic DM and material DM
514 * @param prb_name problem name
515 * @param dm_elastic elastic DM
516 * @param dm_material material DM
517 * @param bit domain bit ref level
518 * @param mask mask ref level for problem
519 * @return error code
520 */
522 createCrackPropagationDM(SmartPetscObj<DM> &dm, const std::string prb_name,
523 SmartPetscObj<DM> dm_elastic,
524 SmartPetscObj<DM> dm_material, const BitRefLevel bit,
525 const BitRefLevel mask,
526 const std::vector<std::string> fe_list);
527
528 /** \brief Create elastic problem DM
529 * @param dm elastic DM
530 * @param prb_name problem name
531 * @param bit domain bit ref level
532 * @param mask mask ref level for problem
533 * @return error code
534 */
535 MoFEMErrorCode createElasticDM(SmartPetscObj<DM> &dm,
536 const std::string prb_name,
537 const BitRefLevel bit,
538 const BitRefLevel mask = BitRefLevel().set());
539
540 /** \brief Create elastic problem DM
541 * @param dm elastic DM
542 * @param prb_name problem name
543 * @param bit domain bit ref level
544 * @param mask mask ref level for problem
545 * @return error code
546 */
548 createEigenElasticDM(SmartPetscObj<DM> &dm, const std::string prb_name,
549 const BitRefLevel bit,
550 const BitRefLevel mask = BitRefLevel().set());
551
552 /** \brief Create problem to calculate boundary conditions
553 * @param dm newly created boundary condition DM
554 * @param prb_name problem name
555 * @param bit auxiliary bit ref level
556 * @param mask mask ref level for problem
557 * @return error code
558 * \note BcDM is only useed for testing, to enforce boundary conditions
559staisifing anlylitical solution to which numerical solution is compared.
560 */
561 MoFEMErrorCode createBcDM(SmartPetscObj<DM> &dm, const std::string prb_name,
562 const BitRefLevel bit,
563 const BitRefLevel mask = BitRefLevel().set());
564
565 /** \brief Create DM fto calculate material problem
566 * @param dm used to calculate material forces, sub-dm of
567 * dm_crack_propagation
568 * @param prb_name problem name
569 * @param bit material domain
570 * @param mask dof mask ref level for problem
571 * @param fe_list name of elements on surface elements
572 * @param debug flag for debugging
573 * @return error code
574 */
575 MoFEMErrorCode createMaterialDM(SmartPetscObj<DM> &dm, const std::string prb_name,
576 const BitRefLevel bit, const BitRefLevel mask,
577 const std::vector<std::string> fe_list,
578 const bool debug = false);
579
580 /** \brief Create DM for calculation of material forces (sub DM of DM
581 * material)
582 * @param dm used to calculate material forces, sub-dm of
583 * dm_crack_propagation
584 * @param prb_name problem name
585 * @param verb compilation parameter determining the amount
586 * of information printed
587 * @return error code
588 */
589 MoFEMErrorCode createMaterialForcesDM(SmartPetscObj<DM> &dm,
590 SmartPetscObj<DM> dm_material,
591 const std::string prb_name,
592 const int verb = QUIET);
593
594 /** \brief create DM to calculate projection matrices (sub DM of DM
595 * material)
596 * @param dm DM to project material forces on surfaces,
597 * sub-dm of dm_crack_propagation
598 * @param dm_material material DM
599 * @param prb_name problem name
600 * @param surface_ids IDs of surface meshsets
601 * @param fe_list name of elements on surface elements
602 * @param verb compilation parameter determining the amount
603 * of information printed
604 * @return error code
605 */
607 SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
608 const std::string prb_name, const std::vector<int> surface_ids,
609 const std::vector<std::string> fe_list, const int verb = QUIET);
610
611 /** \brief create DM to calculate Griffith energy
612 * @param dm dm to project material forces on crack
613 * surfaces, sub-dm of dm_crack_propagation
614 * @param dm_material material DM
615 * @param prb_name problem name
616 * @param verb compilation parameter determining the amount
617 * of information printed
618 * @return error code
619 */
620 MoFEMErrorCode createCrackFrontAreaDM(SmartPetscObj<DM> &dm,
621 SmartPetscObj<DM> dm_material,
622 const std::string prb_name,
623 const bool verb = QUIET);
624
625 /**@}*/
626
627 /** \name Create finite element instances */
628
629 /**@{*/
630
631 /** \brief create elastic finite element instance for spatial assembly
632 * @param mwls_stress_tag_name Name of internal stress tag
633 * @param close_crack if true, crack surface is closed
634 * @param verb compilation parameter determining the amount
635 * of information printed (not in use here)
636 * @param debug flag for debugging
637 * @return error code
638 */
639 MoFEMErrorCode assembleElasticDM(const std::string mwls_stress_tag_name,
640 const int verb = QUIET,
641 const bool debug = false);
642
643 /** \brief create elastic finite element instance for spatial assembly
644 * @param verb compilation parameter determining the amount
645 * of information printed (not in use here)
646 * @param debug flag for debugging
647 * @return error code
648 */
649 MoFEMErrorCode assembleElasticDM(const int verb = QUIET,
650 const bool debug = false);
651
652 // \brief add elastic element instances to SNES
654 DM dm, Mat m, Vec q, Vec f,
655 boost::shared_ptr<FEMethod> arc_method = boost::shared_ptr<FEMethod>(),
656 boost::shared_ptr<ArcLengthCtx> arc_ctx = nullptr, const int verb = QUIET,
657 const bool debug = false);
658
659 /// \brief create material element instance
660 MoFEMErrorCode assembleMaterialForcesDM(DM dm, const int verb = QUIET,
661 const bool debug = false);
662
663 /// \brief create smoothing element instance
664 MoFEMErrorCode assembleSmootherForcesDM(DM dm, const std::vector<int> ids,
665 const int verb = QUIET,
666 const bool debug = false);
667
668 /// \brief assemble coupling element instances
669 MoFEMErrorCode assembleCouplingForcesDM(DM dm, const int verb = QUIET,
670 const bool debug = false);
671
672 /// \brief add material elements instances to SNES
673 MoFEMErrorCode addMaterialFEInstancesToSnes(DM dm, const bool fix_crack_front,
674 const int verb = QUIET,
675 const bool debug = false);
676
677 /// \brief add softening elements instances to SNES
679 const bool fix_crack_front,
680 const int verb = QUIET,
681 const bool debug = false);
682
683 /// \brief add finite element to SNES for crack propagation problem
685 addPropagationFEInstancesToSnes(DM dm, boost::shared_ptr<FEMethod> arc_method,
686 boost::shared_ptr<ArcLengthCtx> arc_ctx,
687 const std::vector<int> &surface_ids,
688 const int verb = QUIET,
689 const bool debug = false);
690
691 /** \brief test LHS Jacobians
692 * @param bit bit level
693 * @param mask bitref mask
694 * @param test_nb name of the test
695 */
697 tangent_tests test);
698
700 addMWLSStressOperators(boost::shared_ptr<CrackFrontElement> &fe_rhs,
701 boost::shared_ptr<CrackFrontElement> &fe_lhs);
702
704 addMWLSDensityOperators(boost::shared_ptr<CrackFrontElement> &fe_rhs,
705 boost::shared_ptr<CrackFrontElement> &fe_lhs);
706
707 /// \brief Update fixed nodes
708 MoFEMErrorCode updateMaterialFixedNode(const bool fix_front,
709 const bool fix_small_g,
710 const bool debug = false);
711
712 /**@}*/
713
714 /** \name Assemble and calculate */
715
716 /**@{*/
717
718 /// \brief assemble elastic part of matrix, by running elastic finite element
719 /// instance
720 MoFEMErrorCode calculateElasticEnergy(DM dm, const std::string msg = "");
721
722 /** \brief assemble material forces, by running material finite element
723 * instance
724 * @param dm dm used to calculate material forces, sub-dm
725 * of dm_crack_propagation
726 * @param q DoF sub-vector for crack propagation DM
727 * @param f Rhs sub-vector for crack propagation DM
728 * @param verb compilation parameter determining the amount
729 * of information printed
730 * @param debug flag for debugging
731 * @return error code
732 */
734 const int verb = QUIET,
735 const bool debug = false);
736
737 /** \brief assemble smoothing forces, by running material finite element
738 * instance
739 * @param dm crack propagation dm that is composition of
740 * elastic DM and material DM
741 * @param q DoF sub-vector for crack propagation DM
742 * @param f Rhs sub-vector for crack propagation DM
743 * @param verb compilation parameter determining the amount
744 * of information printed
745 * @param debug flag for debugging
746 * @return error code
747 */
749 const int verb = QUIET,
750 const bool debug = false);
751
753 const bool debug = true);
754
755 /// \brief assemble projection matrices
756 MoFEMErrorCode calculateSurfaceProjectionMatrix(DM dm_front, DM dm_project,
757 const std::vector<int> &ids,
758 const int verb = QUIET,
759 const bool debug = false);
760
761 /// \brief assemble crack front projection matrix (that constrains crack area
762 /// growth)
763 MoFEMErrorCode calculateFrontProjectionMatrix(DM dm_surface, DM dm_project,
764 const int verb = QUIET,
765 const bool debug = false);
766
767 /** \brief project material forces along the crack elongation direction
768 * @param dm_project material forces DM, sub-dm of
769 * dm_crack_propagation
770 * @param f_proj Rhs sub-vector of projected material forces
771 * for material forces DM
772 * @param f Rhs sub-vector for material forces DM
773 * @param verb compilation parameter determining the amount
774 * of information printed
775 * @param debug flag for debugging
776 * @return error code
777 */
778 MoFEMErrorCode projectMaterialForcesDM(DM dm_project, Vec f, Vec f_proj,
779 const int verb = QUIET,
780 const bool debug = false);
781
782 /// \brief calculate Griffith (driving) force
783 MoFEMErrorCode calculateGriffithForce(DM dm, const double gc, Vec f_griffith,
784 const int verb = QUIET,
785 const bool debug = false);
786
787 /// \brief project Griffith forces
788 MoFEMErrorCode projectGriffithForce(DM dm, Vec f_griffith,
789 Vec f_griffith_proj,
790 const int verb = QUIET,
791 const bool debug = false);
792
793 /// \brief calculate release energy
794 MoFEMErrorCode calculateReleaseEnergy(DM dm, Vec f_material_proj,
795 Vec f_griffith_proj, Vec f_lambda,
796 const double gc, const int verb = QUIET,
797 const bool debug = true);
798
799 /** \name Solvers */
800
801 /**@{*/
802
803 /** \brief solve elastic problem
804 * @param dm elastic DM
805 * @param snes snes solver for elastic DM
806 * @param m stiffness sub-matrix for elastic DM
807 * @param q DoF sub-vector for elastic DM
808 * @param f Rhs sub-vector for material DM
809 * @param snes_set_up boolean flag wheather to determin snes solver
810 * @param shell_m stiffness sub-matrix for elastic DM
811 * @return error code
812 */
813 MoFEMErrorCode solveElasticDM(DM dm, SNES snes, Mat m, Vec q, Vec f,
814 bool snes_set_up, Mat *shell_m);
815
816 /** \brief solve crack propagation problem
817 * @param dm crack propagation DM that is composed of
818 * elastic DM and material DM
819 * @param dm_elastic elastic DM
820 * @param snes snes solver for crack propagation DM
821 * @param m stiffness sub-matrix for crack propagation DM
822 * @param q DoF sub-vector for crack propagation DM
823 * @param f Rhs sub-vector for crack propagation DM
824 * @param snes_set_up boolean flag wheather to determin snes solver
825 * @return error code
826 */
827 MoFEMErrorCode solvePropagationDM(DM dm, DM dm_elastic, SNES snes, Mat m,
828 Vec q, Vec f);
829
830 /**@}*/
831
832 /**@}*/
833
834 /** \name Postprocess results */
835
836 /**@{*/
837
839 const int verb = QUIET,
840 const bool debug = false);
841
842 MoFEMErrorCode writeCrackFont(const BitRefLevel &bit, const int step);
843
844 /// \brief post-process results for elastic solution
845 MoFEMErrorCode postProcessDM(DM dm, const int step, const std::string fe_name,
846 const bool approx_internal_stress);
847
848 /**@}*/
849
850 /** \name Set data from mesh to vectors and back */
851
852 /**@{*/
853
854 /// \brief set field from node positions
856
857 /// \brief set material field from nodes
859
860 /// \brief set spatial field from nodes
862
863 /// \brief set coordinates from field
865 setCoordsFromField(const std::string field_name = "MESH_NODE_POSITIONS");
866
867 /// \brief set singular dofs (on edges adjacent to crack front) from geometry
869 const int verb = QUIET);
870
871 /// \brief set singular dofs of material field (on edges adjacent to crack
872 /// front) from geometry
874
875 /// \brief remove singularity
877
878 /// \brief set maetrial field order to one
880
881 /**@}*/
882
883 /** \name Auxiliary method */
884
885 /**@{*/
886
887 boost::shared_ptr<DofEntity> arcLengthDof;
888
889 /// \brief set pointer to arc-length DOF
891
892 inline boost::shared_ptr<ArcLengthCtx> &getArcCtx() { return arcCtx; }
893 inline boost::shared_ptr<ArcLengthCtx> &getEigenArcCtx() {
894 return arcEigenCtx;
895 }
896
897 inline boost::shared_ptr<FEMethod>
898 getFrontArcLengthControl(boost::shared_ptr<ArcLengthCtx> arc_ctx) {
899 boost::shared_ptr<FEMethod> arc_method(
902 griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA",
903 arc_ctx));
904 return arc_method;
905 }
906
907 inline boost::shared_ptr<AnalyticalDirichletBC::DirichletBC>
910 }
911
912 boost::shared_ptr<MWLSApprox> &getMWLSApprox() { return mwlsApprox; }
913
914 /**@}*/
915
916 /**
917 * \brief operator to post-process results on crack front nodes
918 */
920
923 std::string tagName;
924 PostProcVertexMethod(MoFEM::Interface &m_field, Vec f, std::string tag_name)
925 : mField(m_field), F(f), tagName(tag_name) {}
926
927 Tag tH;
928 double *aRray;
930
932
934 };
935
936 /**
937 * \brief Constrains material displacement on both sides
938 *
939 * Both sides of crack surface have the same material displacements.
940 *
941 */
943
944 /**
945 * @brief Construct a new Both Surface Constrains object
946 *
947 * @param m_field
948 * @param lambda_field_name lagrange multiplayers
949 * @param spatial_field_name constrained spatial field
950 */
952 MoFEM::Interface &m_field,
953 const std::string lambda_field_name = "LAMBDA_BOTH_SIDES",
954 const std::string spatial_field_name = "MESH_NODE_POSITIONS")
955 : mField(m_field), lambdaFieldName(lambda_field_name),
956 spatialFieldName(spatial_field_name), aLpha(1) {
957 ierr = getOptions();
958 CHKERRABORT(PETSC_COMM_WORLD, ierr);
959 }
960
963 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
964 "Get Booth surface constrains element scaling",
965 "none");
966 CHKERRQ(ierr);
967 CHKERR PetscOptionsScalar("-booth_side_alpha", "scaling parameter", "",
968 aLpha, &aLpha, PETSC_NULL);
969 ierr = PetscOptionsEnd();
970 CHKERRQ(ierr);
972 }
973
977
978 double aLpha; ///! Scaling parameter
979 Range masterNodes; ///! Node on which lagrange multiplier is set
980
981 private:
983 const std::string lambdaFieldName;
984 const std::string spatialFieldName;
985 };
986
987 /// \brief Determine face orientation
993 FaceOrientation(bool crack_front, Range contact_faces,
994 BitRefLevel bit_level)
995 : useProjectionFromCrackFront(crack_front), contactFaces(contact_faces),
996 defaultBitLevel(bit_level) {}
997 virtual ~FaceOrientation() {}
999 const EntityHandle face,
1000 const BitRefLevel bit);
1002 const FEMethod *fe_method_ptr) {
1004 EntityHandle face = fe_method_ptr->numeredEntFiniteElementPtr->getEnt();
1005 BitRefLevel bit = fe_method_ptr->problemPtr->getBitRefLevel();
1006 CHKERR getElementOrientation(m_field, face, bit);
1008 }
1010 const EntityHandle face) {
1014 }
1015 };
1016
1017 struct ArcLengthSnesCtx : public SnesCtx {
1018
1019 ArcLengthSnesCtx(MoFEM::Interface &m_field, const std::string problem_name,
1020 boost::shared_ptr<ArcLengthCtx> &arc_ptr)
1021 : SnesCtx(m_field, problem_name), arcPtr(arc_ptr.get()) {
1022 MOFEM_LOG("CPWorld", Sev::noisy) << "ArcLengthSnesCtx create";
1023 }
1025 MOFEM_LOG("CPWorld", Sev::noisy) << "ArcLengthSnesCtx destroyed";
1026 };
1027
1028 inline auto getArcPtr() { return arcPtr; }
1029 inline void setVecDiagM(Vec diag_m) { diagM = diag_m; }
1030 inline auto getVecDiagM() { return diagM; }
1031
1032 private:
1033
1036
1037 };
1038
1039 // private:
1040
1041 /** \name Pointers to finite element instances */
1042
1043 /**@{*/
1044
1045 bool onlyHooke; ///< True if only Hooke material is applied
1046 PetscBool onlyHookeFromOptions; ///< True if only Hooke material is applied
1047
1048 /// Integrate spatial stresses, matrix and vector (not only element but all
1049 /// data structure)
1050 boost::shared_ptr<NonlinearElasticElement> elasticFe;
1051 boost::shared_ptr<NonlinearElasticElement> materialFe;
1052 boost::shared_ptr<CrackFrontElement> feLhs; ///< Integrate elastic FE
1053 boost::shared_ptr<CrackFrontElement> feRhs; ///< Integrate elastic FE
1054 boost::shared_ptr<CrackFrontElement>
1055 feMaterialRhs; ///< Integrate material stresses, assemble vector
1056 boost::shared_ptr<CrackFrontElement>
1057 feMaterialLhs; ///< Integrate material stresses, assemble matrix
1058 boost::shared_ptr<CrackFrontElement> feEnergy; ///< Integrate energy
1059 boost::shared_ptr<ConstrainMatrixCtx>
1060 projSurfaceCtx; ///< Data structure to project on the body surface
1061 boost::shared_ptr<ConstrainMatrixCtx>
1062 projFrontCtx; ///< Data structure to project on crack front
1063 boost::shared_ptr<MWLSApprox>
1064 mwlsApprox; ///< Data struture for MWLS approximation
1065 /// Surface elment to calculate tractions in material space
1066 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
1068 boost::shared_ptr<DirichletSpatialPositionsBc>
1069 spatialDirichletBc; ///< apply Dirichlet BC to sparial positions
1070 boost::shared_ptr<AnalyticalDirichletBC::DirichletBC>
1071 analyticalDirichletBc; ///< apply analytical Dirichlet BC to sparial
1072 ///< positions
1073 boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
1074 surfaceForces; ///< assemble surface forces
1075 boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
1076 surfacePressure; ///< assemble surface pressure
1077 boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
1078 surfacePressureAle; ///< assemble surface pressure (ALE)
1079 boost::shared_ptr<NeumannForcesSurface::DataAtIntegrationPts>
1080 commonDataSurfacePressureAle; ///< common data at integration points (ALE)
1081 boost::shared_ptr<boost::ptr_map<string, EdgeForce>>
1082 edgeForces; ///< assemble edge forces
1083 boost::shared_ptr<boost::ptr_map<string, NodalForce>>
1084 nodalForces; ///< assemble nodal forces
1085 boost::shared_ptr<FEMethod> assembleFlambda; ///< assemble F_lambda vector
1086 boost::shared_ptr<FEMethod> zeroFlambda; ///< assemble F_lambda vector
1087
1088 boost::shared_ptr<CrackFrontElement>
1089 feCouplingElasticLhs; ///< FE instance to assemble coupling terms
1090 boost::shared_ptr<CrackFrontElement>
1091 feCouplingMaterialLhs; ///< FE instance to assemble coupling terms
1092
1093 boost::shared_ptr<Smoother> smootherFe;
1094 boost::shared_ptr<Smoother::MyVolumeFE>
1095 feSmootherRhs; ///< Integrate smoothing operators
1096 boost::shared_ptr<Smoother::MyVolumeFE>
1097 feSmootherLhs; ///< Integrate smoothing operators
1098 boost::shared_ptr<VolumeLengthQuality<double>> volumeLengthDouble;
1099 boost::shared_ptr<VolumeLengthQuality<adouble>> volumeLengthAdouble;
1100
1101 boost::shared_ptr<
1103 tangentConstrains; ///< Constrains crack front in tangent directiona
1104 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
1106 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
1108 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
1110 map<int, boost::shared_ptr<SurfaceSlidingConstrains>> surfaceConstrain;
1111 map<int, boost::shared_ptr<EdgeSlidingConstrains>> edgeConstrains;
1112
1113 boost::shared_ptr<BothSurfaceConstrains> bothSidesConstrains;
1114 boost::shared_ptr<BothSurfaceConstrains> closeCrackConstrains;
1115 boost::shared_ptr<BothSurfaceConstrains> bothSidesContactConstrains;
1116 boost::shared_ptr<DirichletFixFieldAtEntitiesBc> fixMaterialEnts;
1117
1118 boost::shared_ptr<GriffithForceElement> griffithForceElement;
1119 boost::shared_ptr<GriffithForceElement::MyTriangleFE> feGriffithForceRhs;
1120 boost::shared_ptr<GriffithForceElement::MyTriangleFE> feGriffithForceLhs;
1121 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrainsDelta>
1123 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>
1125 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>
1127
1128 boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
1129 analiticalSurfaceElement; ///< Finite element to integrate tractions given
1130 ///< by analytical formula
1131 boost::shared_ptr<FaceElementForcesAndSourcesCore> feSpringLhsPtr;
1132 boost::shared_ptr<FaceElementForcesAndSourcesCore> feSpringRhsPtr;
1133 /**@}*/
1134
1135 /** \name Entity ranges */
1136
1137 /**@{*/
1138
1150
1151 EntityHandle
1152 arcLengthVertex; ///< Vertex associated with dof for arc-length control
1153
1154 /**@}*/
1155
1156 /** \name Auxiliary data */
1157
1158 /**@{*/
1159
1161
1165
1166 double maxG1;
1167 double maxJ;
1169 map<EntityHandle, double> mapG1; ///< hashmap of g1 - release energy at nodes
1170 map<EntityHandle, double> mapG3; ///< hashmap of g3 - release energy at nodes
1171 map<EntityHandle, double> mapJ; ///< hashmap of J - release energy at nodes
1172 map<EntityHandle, VectorDouble3>
1173 mapMatForce; ///< hashmap of material force at nodes
1174 map<EntityHandle, VectorDouble3>
1175 mapGriffith; ///< hashmap of Griffith energy at nodes
1176
1177 map<EntityHandle, double> mapSmoothingForceFactor;
1178
1179 bool doSurfaceProjection; ///< If crack surface is immersed in the body do not
1180 ///< projection
1181
1182 double smootherAlpha; ///< Controls mesh smoothing
1183 double initialSmootherAlpha; ///< Controls mesh smoothing, set at start of
1184 ///< analysis
1185 double smootherGamma; ///< Controls mesh smoothing
1186
1187 /**@}*/
1188
1189 /** \name arc-length options */
1190
1191 /**@{*/
1192
1193 double arcAlpha;
1194 double arcBeta;
1195 double arcS;
1196
1197 boost::shared_ptr<ArcLengthCtx> arcCtx;
1198 boost::shared_ptr<ArcLengthCtx> arcEigenCtx;
1199
1200 /**@}*/
1201
1202 /** \name Contact parameters, entity ranges and data structures */
1203
1204 /**@{*/
1205
1206 int contactOrder; ///< Approximation order for spatial positions
1207 ///< used in contact elements
1208 int contactLambdaOrder; ///< Approximation order for field of Lagrange
1209 ///< multipliers used to enforce contact
1210 double rValue; ///< Parameter for regularizing absolute value
1211 ///< function in the complementarity function
1212 double cnValue; ///< Augmentation parameter in the complementarity function
1213 ///< used to enforce KKT constraints
1214
1215 PetscBool ignoreContact; ///< If set true, contact interfaces are ignored
1216 ///< and prism interfaces are not created
1217 PetscBool fixContactNodes; ///< If set true, contact nodes are fixed in
1218 ///< the material configuration
1220
1222
1226
1230
1231 boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
1237
1238 EntityHandle meshsetFaces;
1239
1240 boost::shared_ptr<FaceElementForcesAndSourcesCore> feLhsSpringALE;
1241 boost::shared_ptr<FaceElementForcesAndSourcesCore> feLhsSpringALEMaterial;
1242 boost::shared_ptr<FaceElementForcesAndSourcesCore> feRhsSpringALEMaterial;
1243 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
1245
1246 boost::shared_ptr<SimpleContactProblem> contactProblem;
1247 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1249 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1251 boost::shared_ptr<SimpleContactProblem::CommonDataSimpleContact>
1253
1254 boost::shared_ptr<MortarContactProblem> mortarContactProblemPtr;
1255
1256 boost::shared_ptr<MortarContactProblem::MortarContactElement>
1258 boost::shared_ptr<MortarContactProblem::MortarContactElement>
1260 boost::shared_ptr<MortarContactProblem::CommonDataMortarContact>
1262 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1264 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1266 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1268 boost::shared_ptr<SimpleContactProblem::CommonDataSimpleContact>
1270
1272 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1274 boost::shared_ptr<MortarContactProblem::MortarContactElement>
1278
1279 Range constrainedInterface; ///< Range of faces on the constrained interface
1280
1281 /**@}*/
1282
1283 /** \name Data for bone */
1284
1285 /**@{*/
1286
1287 double rHo0; ///< Reference density if bone is analyzed
1288 double nBone; ///< Exponent parameter in bone density
1289
1290 /**@}*/
1291
1292 /** \name Interfaces */
1293
1294 /**@{*/
1295
1296 const boost::scoped_ptr<UnknownInterface> cpSolversPtr;
1297 const boost::scoped_ptr<UnknownInterface> cpMeshCutPtr;
1298
1299 /**@}*/
1300
1301 /** \name Cut/Refine/Split options */
1302
1303 /**@{*/
1304
1305 PetscBool doCutMesh;
1308
1309 /**@}*/
1310};
1311
1312} // namespace FractureMechanics
1313
1314#endif // __CRACK_PROPAGATION_HPP__
@ QUIET
Definition: definitions.h:208
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'm', SPACE_DIM > m
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
auto bit
set bit
const FTensor::Tensor2< T, Dim, Dim > Vec
MoFEMErrorCode broadcast_entities(moab::Interface &moab, moab::Interface &moab_tmp, const int from_proc, Range &entities, const bool adjacencies, const bool tags)
MoFEMErrorCode clean_pcomms(moab::Interface &moab, boost::shared_ptr< WrapMPIComm > moab_comm_wrap)
AdolcTags
Tapes numbers used by ADOL-C.
tangent_tests
Names of tangent matrices tests.
auto f
Definition: HenckyOps.hpp:5
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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
constexpr auto field_name
static constexpr int approx_order
Store variables for ArcLength analysis.
ArcLengthSnesCtx(MoFEM::Interface &m_field, const std::string problem_name, boost::shared_ptr< ArcLengthCtx > &arc_ptr)
Constrains material displacement on both sides.
BothSurfaceConstrains(MoFEM::Interface &m_field, const std::string lambda_field_name="LAMBDA_BOTH_SIDES", const std::string spatial_field_name="MESH_NODE_POSITIONS")
Construct a new Both Surface Constrains object.
MoFEM::Interface & mField
! Node on which lagrange multiplier is set
FaceOrientation(bool crack_front, Range contact_faces, BitRefLevel bit_level)
MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const EntityHandle face)
MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const FEMethod *fe_method_ptr)
MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const EntityHandle face, const BitRefLevel bit)
operator to post-process results on crack front nodes
MoFEM::Interface & mField
Tag tH
Vec F
double * aRray
MoFEMErrorCode operator()()
function is run for every finite element
std::string tagName
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
PostProcVertexMethod(MoFEM::Interface &m_field, Vec f, std::string tag_name)
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringRhsPtr
MoFEMErrorCode zeroLambdaFields()
Zero fields with lagrange multipliers.
boost::shared_ptr< BothSurfaceConstrains > bothSidesConstrains
MoFEMErrorCode solveElasticDM(DM dm, SNES snes, Mat m, Vec q, Vec f, bool snes_set_up, Mat *shell_m)
solve elastic problem
boost::shared_ptr< AnalyticalDirichletBC::DirichletBC > analyticalDirichletBc
MoFEMErrorCode declareSpringsAleFE(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
Declare FE for spring BC in ALE formulation (in material domain)
boost::shared_ptr< CrackFrontElement > feMaterialRhs
Integrate material stresses, assemble vector.
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringLhsPtr
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressure
assemble surface pressure
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > contactOrientation
MoFEMErrorCode getArcLengthDof()
set pointer to arc-length DOF
MoFEMErrorCode declareSimpleContactAleFE(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
Declare FE for pressure BC in ALE formulation (in material domain)
MoFEMErrorCode createMaterialDM(SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask, const std::vector< std::string > fe_list, const bool debug=false)
Create DM fto calculate material problem.
MoFEMErrorCode deleteEntities(const int verb=QUIET, const bool debug=false)
MoFEMErrorCode declareSmoothingFE(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
declare mesh smoothing finite elements
MoFEMErrorCode declareBothSidesFE(const BitRefLevel bit_spatial, const BitRefLevel bit_material, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
boost::shared_ptr< DirichletSpatialPositionsBc > spatialDirichletBc
apply Dirichlet BC to sparial positions
MoFEMErrorCode declareExternalForcesFE(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
MoFEMErrorCode solvePropagationDM(DM dm, DM dm_elastic, SNES snes, Mat m, Vec q, Vec f)
solve crack propagation problem
std::string mwlsApproxFile
Name of file with internal stresses.
map< int, boost::shared_ptr< SurfaceSlidingConstrains > > surfaceConstrain
MoFEMErrorCode readMedFile()
read mesh file
boost::shared_ptr< CrackFrontElement > feRhs
Integrate elastic FE.
MoFEMErrorCode buildBothSidesFieldId(const BitRefLevel bit_spatial, const BitRefLevel bit_material, const bool proc_only=false, const bool build_fields=true, const int verb=QUIET, const bool debug=false)
Lagrange multipliers field which constrains material displacements.
PetscBool propagateCrack
If true crack propagation is calculated.
boost::shared_ptr< MortarContactProblem::MortarContactElement > fePostProcMortarContact
std::map< std::string, BitRefLevel > mapBitLevel
MoFEMErrorCode declareArcLengthFE(const BitRefLevel bits, const int verb=QUIET)
create arc-length element entity and declare elemets
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsLhs
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherLhs
Integrate smoothing operators.
double rHo0
Reference density if bone is analyzed.
double smootherGamma
Controls mesh smoothing.
MoFEMErrorCode setCrackFrontBitLevel(BitRefLevel from_bit, BitRefLevel bit, const int nb_levels=2, const bool debug=false)
Set bit ref level for entities adjacent to crack front.
MoFEMErrorCode testJacobians(const BitRefLevel bit, const BitRefLevel mask, tangent_tests test)
test LHS Jacobians
MoFEMErrorCode declareEdgeFE(std::string fe_name, const BitRefLevel bit, const BitRefLevel mask, const bool proc_only=true, const int verb=QUIET, const bool debug=false)
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContact
MoFEMErrorCode createSurfaceProjectionDM(SmartPetscObj< DM > &dm, SmartPetscObj< DM > dm_material, const std::string prb_name, const std::vector< int > surface_ids, const std::vector< std::string > fe_list, const int verb=QUIET)
create DM to calculate projection matrices (sub DM of DM material)
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > analiticalSurfaceElement
boost::shared_ptr< ConstrainMatrixCtx > projFrontCtx
Data structure to project on crack front.
PetscBool onlyHookeFromOptions
True if only Hooke material is applied.
MoFEMErrorCode getInterfaceVersion(Version &version) const
MoFEMErrorCode buildProblemFiniteElements(BitRefLevel bit1, BitRefLevel bit2, std::vector< int > surface_ids, const int verb=QUIET, const bool debug=false)
Build problem finite elements.
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsDelta > feGriffithConstrainsDelta
boost::shared_ptr< NonlinearElasticElement > elasticFe
boost::shared_ptr< MWLSApprox > & getMWLSApprox()
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 declareFrontFE(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
MoFEMErrorCode buildProblemFields(const BitRefLevel &bit1, const BitRefLevel &mask1, const BitRefLevel &bit2, const int verb=QUIET, const bool debug=false)
Build problem fields.
MoFEMErrorCode postProcessDM(DM dm, const int step, const std::string fe_name, const bool approx_internal_stress)
post-process results for elastic solution
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContactALE
boost::shared_ptr< SimpleContactProblem::CommonDataSimpleContact > commonDataSimpleContact
double betaGc
heterogeneous Griffith energy exponent
CrackPropagation(MoFEM::Interface &m_field, const int approx_order=2, const int geometry_order=1)
Constructor of Crack Propagation, prints fracture module version and registers the three interfaces u...
EntityHandle arcLengthVertex
Vertex associated with dof for arc-length control.
double griffithR
Griffith regularisation parameter.
std::string mwlsEigenStressTagName
Name of tag with eigen stresses.
MoFEMErrorCode partitionMesh(BitRefLevel bit1, BitRefLevel bit2, int verb=QUIET, const bool debug=false)
partotion mesh
PetscBool isGcHeterogeneous
flag for heterogeneous gc
MoFEMErrorCode addMWLSStressOperators(boost::shared_ptr< CrackFrontElement > &fe_rhs, boost::shared_ptr< CrackFrontElement > &fe_lhs)
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherRhs
Integrate smoothing operators.
std::string mwlsStressTagName
Name of tag with internal stresses.
map< EntityHandle, double > mapJ
hashmap of J - release energy at nodes
boost::shared_ptr< DirichletFixFieldAtEntitiesBc > fixMaterialEnts
boost::shared_ptr< FEMethod > zeroFlambda
assemble F_lambda vector
boost::shared_ptr< NonlinearElasticElement > materialFe
boost::shared_ptr< boost::ptr_map< string, EdgeForce > > edgeForces
assemble edge forces
MoFEMErrorCode buildSurfaceFields(const BitRefLevel bit, const bool proc_only=true, const bool build_fields=true, const int verb=QUIET, const bool debug=false)
build fields with Lagrange multipliers to constrain surfaces
MoFEMErrorCode setSingularDofs(const string field_name, const int verb=QUIET)
set singular dofs (on edges adjacent to crack front) from geometry
MoFEMErrorCode savePositionsOnCrackFrontDM(DM dm, Vec q, const int verb=QUIET, const bool debug=false)
boost::shared_ptr< NeumannForcesSurface::MyTriangleFE > feMaterialAnaliticalTraction
Surface elment to calculate tractions in material space.
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsRhs
MoFEMErrorCode getOptions()
Get options form command line.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
MoFEMErrorCode assembleMaterialForcesDM(DM dm, const int verb=QUIET, const bool debug=false)
create material element instance
map< EntityHandle, double > mapSmoothingForceFactor
MoFEMErrorCode calculateGriffithForce(DM dm, const double gc, Vec f_griffith, const int verb=QUIET, const bool debug=false)
calculate Griffith (driving) force
boost::shared_ptr< FEMethod > getFrontArcLengthControl(boost::shared_ptr< ArcLengthCtx > arc_ctx)
boost::shared_ptr< MortarContactProblem > mortarContactProblemPtr
boost::shared_ptr< FaceElementForcesAndSourcesCore > feLhsSpringALEMaterial
boost::shared_ptr< Smoother > smootherFe
MoFEMErrorCode setSpatialPositionFromCoords()
set spatial field from nodes
boost::shared_ptr< MortarContactProblem::CommonDataMortarContact > commonDataMortarContact
boost::shared_ptr< FaceElementForcesAndSourcesCore > feRhsSpringALEMaterial
MoFEMErrorCode updateMaterialFixedNode(const bool fix_front, const bool fix_small_g, const bool debug=false)
Update fixed nodes.
const boost::scoped_ptr< UnknownInterface > cpSolversPtr
MoFEMErrorCode createBcDM(SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
Create problem to calculate boundary conditions.
double smootherAlpha
Controls mesh smoothing.
MoFEMErrorCode createEigenElasticDM(SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
Create elastic problem DM.
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > skinOrientation
MoFEMErrorCode buildElasticFields(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool build_fields=true, const int verb=QUIET, const bool debug=false)
Declate fields for elastic analysis.
MoFEMErrorCode setSingularElementMatrialPositions(const int verb=QUIET)
set singular dofs of material field (on edges adjacent to crack front) from geometry
boost::shared_ptr< DofEntity > arcLengthDof
MoFEMErrorCode declarePressureAleFE(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
Declare FE for pressure BC in ALE formulation (in material domain)
MoFEMErrorCode setMaterialPositionFromCoords()
set material field from nodes
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)
map< EntityHandle, VectorDouble3 > mapMatForce
hashmap of material force at nodes
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressureAle
assemble surface pressure (ALE)
MoFEMErrorCode calculateElasticEnergy(DM dm, const std::string msg="")
assemble elastic part of matrix, by running elastic finite element instance
MoFEMErrorCode createMaterialForcesDM(SmartPetscObj< DM > &dm, SmartPetscObj< DM > dm_material, const std::string prb_name, const int verb=QUIET)
Create DM for calculation of material forces (sub DM of DM material)
boost::shared_ptr< VolumeLengthQuality< adouble > > volumeLengthAdouble
MoFEMErrorCode projectGriffithForce(DM dm, Vec f_griffith, Vec f_griffith_proj, const int verb=QUIET, const bool debug=false)
project Griffith forces
double griffithE
Griffith stability parameter.
MoFEMErrorCode declareSurfaceFE(std::string fe_name, const BitRefLevel bit, const BitRefLevel mask, const std::vector< int > &ids, const bool proc_only=true, const int verb=QUIET, const bool debug=false)
declare surface sliding elements
boost::shared_ptr< ArcLengthCtx > arcEigenCtx
Range constrainedInterface
Range of faces on the constrained interface.
MoFEMErrorCode resolveSharedBitRefLevel(const BitRefLevel bit, const int verb=QUIET, const bool debug=false)
resole shared entities by bit level
MoFEMErrorCode assembleSmootherForcesDM(DM dm, const std::vector< int > ids, const int verb=QUIET, const bool debug=false)
create smoothing element instance
boost::shared_ptr< boost::ptr_map< string, NodalForce > > nodalForces
assemble nodal forces
boost::shared_ptr< VolumeLengthQuality< double > > volumeLengthDouble
MoFEMErrorCode setCoordsFromField(const std::string field_name="MESH_NODE_POSITIONS")
set coordinates from field
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 buildEdgeFields(const BitRefLevel bit, const bool proc_only=true, const bool build_fields=true, const int verb=QUIET, const bool debug=false)
build fields with Lagrange multipliers to constrain edges
map< EntityHandle, double > mapG3
hashmap of g3 - release energy at nodes
PetscBool isPressureAle
If true surface pressure is considered in ALE.
MoFEMErrorCode createElasticDM(SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
Create elastic problem DM.
map< EntityHandle, double > mapG1
hashmap of g1 - release energy at nodes
MoFEMErrorCode declareCrackSurfaceFE(std::string fe_name, const BitRefLevel bit, const BitRefLevel mask, const bool proc_only=true, const int verb=QUIET, const bool debug=false)
MoFEMErrorCode createCrackPropagationDM(SmartPetscObj< DM > &dm, const std::string prb_name, SmartPetscObj< DM > dm_elastic, SmartPetscObj< DM > dm_material, const BitRefLevel bit, const BitRefLevel mask, const std::vector< std::string > fe_list)
Create DM by composition of elastic DM and material DM.
boost::shared_ptr< MortarContactProblem::MortarContactElement > feLhsMortarContact
MoFEMErrorCode tetsingReleaseEnergyCalculation()
This is run with ctest.
MoFEMErrorCode saveEachPart(const std::string prefix, const Range &ents)
Save entities on ech processor.
MoFEMErrorCode calculateSmoothingForceFactor(const int verb=QUIET, const bool debug=true)
boost::shared_ptr< ConstrainMatrixCtx > projSurfaceCtx
Data structure to project on the body surface.
static FTensor::Tensor2_symmetric< double, 3 > analyticalStrainFunction(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords)
boost::shared_ptr< ArcLengthCtx > arcCtx
MoFEMErrorCode cleanSingularElementMatrialPositions()
set maetrial field order to one
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > crackOrientation
map< int, boost::shared_ptr< EdgeSlidingConstrains > > edgeConstrains
PetscBool solveEigenStressProblem
Solve eigen problem.
MoFEMErrorCode writeCrackFont(const BitRefLevel &bit, const int step)
MoFEMErrorCode setFieldFromCoords(const std::string field_name)
set field from node positions
boost::shared_ptr< MetaSpringBC::DataAtIntegrationPtsSprings > commonDataSpringsALE
MoFEMErrorCode addMWLSDensityOperators(boost::shared_ptr< CrackFrontElement > &fe_rhs, boost::shared_ptr< CrackFrontElement > &fe_lhs)
MoFEMErrorCode createCrackFrontAreaDM(SmartPetscObj< DM > &dm, SmartPetscObj< DM > dm_material, const std::string prb_name, const bool verb=QUIET)
create DM to calculate Griffith energy
MoFEMErrorCode declareElasticFE(const BitRefLevel bit1, const BitRefLevel mask1, const BitRefLevel bit2, const BitRefLevel mask2, const bool add_forces=true, const bool proc_only=true, const int verb=QUIET)
declare elastic finite elements
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContact
boost::shared_ptr< CrackFrontElement > feCouplingMaterialLhs
FE instance to assemble coupling terms.
double nBone
Exponent parameter in bone density.
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< boost::ptr_map< string, NeumannForcesSurface > > surfaceForces
assemble surface forces
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceLhs
const boost::scoped_ptr< UnknownInterface > cpMeshCutPtr
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContactALEMaterial
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
MoFEMErrorCode addMaterialFEInstancesToSnes(DM dm, const bool fix_crack_front, const int verb=QUIET, const bool debug=false)
add material elements instances to SNES
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
bool onlyHooke
True if only Hooke material is applied.
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > fePostProcSimpleContact
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
boost::shared_ptr< CrackFrontElement > feLhs
Integrate elastic FE.
MoFEMErrorCode calculateSurfaceProjectionMatrix(DM dm_front, DM dm_project, const std::vector< int > &ids, const int verb=QUIET, const bool debug=false)
assemble projection matrices
boost::shared_ptr< BothSurfaceConstrains > bothSidesContactConstrains
MoFEMErrorCode buildArcLengthField(const BitRefLevel bit, const bool build_fields=true, const int verb=QUIET)
Declate field for arc-length.
int postProcLevel
level of postprocessing (amount of output files)
MoFEMErrorCode resolveShared(const Range &tets, Range &proc_ents, const int verb=QUIET, const bool debug=false)
resolve shared entities
boost::shared_ptr< AnalyticalDirichletBC::DirichletBC > getAnalyticalDirichletBc()
boost::shared_ptr< FaceElementForcesAndSourcesCore > feLhsSpringALE
boost::shared_ptr< FEMethod > assembleFlambda
assemble F_lambda vector
boost::shared_ptr< BothSurfaceConstrains > closeCrackConstrains
PetscBool areSpringsAle
If true surface spring is considered in ALE.
boost::shared_ptr< ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain > tangentConstrains
Constrains crack front in tangent directiona.
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
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< ContactSearchKdTree::ContactCommonData_multiIndex > contactSearchMultiIndexPtr
boost::shared_ptr< MortarContactProblem::MortarContactElement > feRhsMortarContact
map< EntityHandle, VectorDouble3 > mapGriffith
hashmap of Griffith energy at nodes
MoFEMErrorCode addSmoothingFEInstancesToSnes(DM dm, const bool fix_crack_front, const int verb=QUIET, const bool debug=false)
add softening elements instances to SNES
MoFEMErrorCode declareMaterialFE(const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
declare material finite elements
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< SimpleContactProblem::CommonDataSimpleContact > commonDataSimpleContactALE
MoFEMErrorCode buildCrackSurfaceFieldId(const BitRefLevel bit, const bool proc_only=true, const bool build_fields=true, const int verb=QUIET, const bool debug=false)
declare crack surface files
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceRhs
MoFEMErrorCode buildCrackFrontFieldId(const BitRefLevel bit, const bool build_fields=true, const int verb=QUIET, const bool debug=false)
declare crack surface files
boost::shared_ptr< CrackFrontElement > feEnergy
Integrate energy.
boost::shared_ptr< SimpleContactProblem > contactProblem
boost::shared_ptr< CrackFrontElement > feMaterialLhs
Integrate material stresses, assemble matrix.
int residualStressBlock
Block on which residual stress is applied.
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContactALEMaterial
MoFEMErrorCode unsetSingularElementMatrialPositions()
remove singularity
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)
std::string mwlsRhoTagName
Name of tag with density.
boost::shared_ptr< NeumannForcesSurface::DataAtIntegrationPts > commonDataSurfacePressureAle
common data at integration points (ALE)
boost::shared_ptr< ArcLengthCtx > & getArcCtx()
Deprecated interface functions.
Data structure to exchange data between mofem and User Loop Methods on entities.
base class for all interface classes
Class implemented by user to detect face orientation.