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

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

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

Classes

struct  ArcLengthSnesCtx
 
struct  BothSurfaceConstrains
 Constrains material displacement on both sides. More...
 
struct  FaceOrientation
 Determine face orientation. More...
 
struct  PostProcVertexMethod
 operator to post-process results on crack front nodes More...
 

Public Member Functions

MoFEMErrorCode getInterfaceVersion (Version &version) const
 
MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 Getting interface of core database. More...
 
 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 used to solve fracture problem, i.e. CrackPropagation, CPSolvers and CPMeshCut. More...
 
virtual ~CrackPropagation ()
 
MoFEMErrorCode getOptions ()
 Get options form command line. More...
 
MoFEMErrorCode tetsingReleaseEnergyCalculation ()
 This is run with ctest. More...
 
int & getNbLoadSteps ()
 
int getNbLoadSteps () const
 
double & getLoadScale ()
 
double getLoadScale () const
 
int & getNbCutSteps ()
 
int getNbCutSteps () const
 
MoFEMErrorCode readMedFile ()
 read mesh file More...
 
MoFEMErrorCode saveEachPart (const std::string prefix, const Range &ents)
 Save entities on ech processor. 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 reference 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
 

Static Public Member Functions

static FTensor::Tensor2_symmetric< double, 3 > analyticalStrainFunction (FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords)
 
- 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...
 

Public Attributes

Version fileVersion
 
boost::shared_ptr< WrapMPIComm > moabCommWorld
 
MoFEM::Interface & mField
 
bool setSingularCoordinates
 
int approxOrder
 
int geometryOrder
 
double gC
 Griffith energy. More...
 
double griffithE
 Griffith stability parameter. More...
 
double griffithR
 Griffith regularisation parameter. More...
 
double betaGc
 heterogeneous Griffith energy exponent More...
 
PetscBool isGcHeterogeneous
 flag for heterogeneous gc More...
 
PetscBool isPartitioned
 
PetscBool propagateCrack
 If true crack propagation is calculated. More...
 
int refAtCrackTip
 
int refOrderAtTip
 
PetscBool addSingularity
 
PetscBool otherSideConstrains
 
PetscBool isPressureAle
 If true surface pressure is considered in ALE. More...
 
PetscBool areSpringsAle
 If true surface spring is considered in ALE. More...
 
PetscBool isSurfaceForceAle
 If true surface pressure is considered in ALE. More...
 
PetscBool ignoreMaterialForce
 If true surface pressure is considered in ALE. More...
 
std::string mwlsApproxFile
 Name of file with internal stresses. More...
 
std::string mwlsStressTagName
 Name of tag with internal stresses. More...
 
std::string mwlsEigenStressTagName
 Name of tag with eigen stresses. More...
 
std::string mwlsRhoTagName
 Name of tag with density. More...
 
int residualStressBlock
 Block on which residual stress is applied. More...
 
int densityMapBlock
 
std::string configFile
 
bool isStressTagSavedOnNodes
 
double partitioningWeightPower
 
PetscBool resetMWLSCoeffsEveryPropagationStep
 
PetscBool addAnalyticalInternalStressOperators
 
PetscBool solveEigenStressProblem
 Solve eigen problem. More...
 
PetscBool useEigenPositionsSimpleContact
 
int postProcLevel
 level of postprocessing (amount of output files) More...
 
int startStep
 
double crackFrontLength
 
double crackSurfaceArea
 
std::map< std::string, BitRefLevel > mapBitLevel
 

Static Public Attributes

static bool debug = false
 
static bool parallelReadAndBroadcast
 

Resolve shared entities and partition mesh

MoFEMErrorCode partitionMesh (BitRefLevel bit1, BitRefLevel bit2, int verb=QUIET, const bool debug=false)
 partotion mesh More...
 
MoFEMErrorCode resolveShared (const Range &tets, Range &proc_ents, const int verb=QUIET, const bool debug=false)
 resolve shared entities More...
 
MoFEMErrorCode resolveSharedBitRefLevel (const BitRefLevel bit, const int verb=QUIET, const bool debug=false)
 resole shared entities by bit level More...
 

Set bit level for material problem

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. More...
 

Build problem

MoFEMErrorCode buildProblemFields (const BitRefLevel &bit1, const BitRefLevel &mask1, const BitRefLevel &bit2, const int verb=QUIET, const bool debug=false)
 Build problem fields. More...
 
MoFEMErrorCode buildProblemFiniteElements (BitRefLevel bit1, BitRefLevel bit2, std::vector< int > surface_ids, const int verb=QUIET, const bool debug=false)
 Build problem finite elements. More...
 
MoFEMErrorCode createProblemDataStructures (const std::vector< int > surface_ids, const int verb=QUIET, const bool debug=false)
 Construct problem data structures. More...
 
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. More...
 
MoFEMErrorCode deleteEntities (const int verb=QUIET, const bool debug=false)
 

Build fields

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. More...
 
MoFEMErrorCode buildArcLengthField (const BitRefLevel bit, const bool build_fields=true, const int verb=QUIET)
 Declate field for arc-length. More...
 
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 More...
 
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 More...
 
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 More...
 
MoFEMErrorCode buildCrackFrontFieldId (const BitRefLevel bit, const bool build_fields=true, const int verb=QUIET, const bool debug=false)
 declare crack surface files More...
 
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. More...
 
MoFEMErrorCode zeroLambdaFields ()
 Zero fields with lagrange multipliers. More...
 

Build finite elements

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 More...
 
MoFEMErrorCode declareArcLengthFE (const BitRefLevel bits, const int verb=QUIET)
 create arc-length element entity and declare elemets More...
 
MoFEMErrorCode declareExternalForcesFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
 
MoFEMErrorCode declareSurfaceForceAleFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
 Declare FE for pressure BC in ALE formulation (in material domain) More...
 
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) More...
 
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) More...
 
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) More...
 
MoFEMErrorCode declareMaterialFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 declare material finite elements More...
 
MoFEMErrorCode declareFrontFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 
MoFEMErrorCode declareBothSidesFE (const BitRefLevel bit_spatial, const BitRefLevel bit_material, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 
MoFEMErrorCode declareSmoothingFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 declare mesh smoothing finite elements More...
 
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 More...
 
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)
 
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)
 

Create DMs & problems

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. More...
 
MoFEMErrorCode createElasticDM (SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
 Create elastic problem DM. More...
 
MoFEMErrorCode createEigenElasticDM (SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
 Create elastic problem DM. More...
 
MoFEMErrorCode createBcDM (SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
 Create problem to calculate boundary conditions. More...
 
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. More...
 
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) More...
 
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) More...
 
MoFEMErrorCode createCrackFrontAreaDM (SmartPetscObj< DM > &dm, SmartPetscObj< DM > dm_material, const std::string prb_name, const bool verb=QUIET)
 create DM to calculate Griffith energy More...
 

Create finite element instances

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 More...
 
MoFEMErrorCode assembleElasticDM (const int verb=QUIET, const bool debug=false)
 create elastic finite element instance for spatial assembly More...
 
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)
 
MoFEMErrorCode assembleMaterialForcesDM (DM dm, const int verb=QUIET, const bool debug=false)
 create material element instance More...
 
MoFEMErrorCode assembleSmootherForcesDM (DM dm, const std::vector< int > ids, const int verb=QUIET, const bool debug=false)
 create smoothing element instance More...
 
MoFEMErrorCode assembleCouplingForcesDM (DM dm, const int verb=QUIET, const bool debug=false)
 assemble coupling element instances More...
 
MoFEMErrorCode addMaterialFEInstancesToSnes (DM dm, const bool fix_crack_front, const int verb=QUIET, const bool debug=false)
 add material elements instances to SNES More...
 
MoFEMErrorCode addSmoothingFEInstancesToSnes (DM dm, const bool fix_crack_front, const int verb=QUIET, const bool debug=false)
 add softening elements instances to SNES More...
 
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 More...
 
MoFEMErrorCode testJacobians (const BitRefLevel bit, const BitRefLevel mask, tangent_tests test)
 test LHS Jacobians More...
 
MoFEMErrorCode addMWLSStressOperators (boost::shared_ptr< CrackFrontElement > &fe_rhs, boost::shared_ptr< CrackFrontElement > &fe_lhs)
 
MoFEMErrorCode addMWLSDensityOperators (boost::shared_ptr< CrackFrontElement > &fe_rhs, boost::shared_ptr< CrackFrontElement > &fe_lhs)
 
MoFEMErrorCode updateMaterialFixedNode (const bool fix_front, const bool fix_small_g, const bool debug=false)
 Update fixed nodes. More...
 

Assemble and calculate

MoFEMErrorCode calculateElasticEnergy (DM dm, const std::string msg="")
 assemble elastic part of matrix, by running elastic finite element instance More...
 
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 More...
 
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 More...
 
MoFEMErrorCode calculateSmoothingForceFactor (const int verb=QUIET, const bool debug=true)
 
MoFEMErrorCode calculateSurfaceProjectionMatrix (DM dm_front, DM dm_project, const std::vector< int > &ids, const int verb=QUIET, const bool debug=false)
 assemble projection matrices More...
 
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) More...
 
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 More...
 
MoFEMErrorCode calculateGriffithForce (DM dm, const double gc, Vec f_griffith, const int verb=QUIET, const bool debug=false)
 calculate Griffith (driving) force More...
 
MoFEMErrorCode projectGriffithForce (DM dm, Vec f_griffith, Vec f_griffith_proj, const int verb=QUIET, const bool debug=false)
 project Griffith forces More...
 
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 More...
 

Solvers

MoFEMErrorCode solveElasticDM (DM dm, SNES snes, Mat m, Vec q, Vec f, bool snes_set_up, Mat *shell_m)
 solve elastic problem More...
 
MoFEMErrorCode solvePropagationDM (DM dm, DM dm_elastic, SNES snes, Mat m, Vec q, Vec f)
 solve crack propagation problem More...
 

Postprocess results

MoFEMErrorCode savePositionsOnCrackFrontDM (DM dm, Vec q, const int verb=QUIET, const bool debug=false)
 
MoFEMErrorCode writeCrackFont (const BitRefLevel &bit, const int step)
 
MoFEMErrorCode postProcessDM (DM dm, const int step, const std::string fe_name, const bool approx_internal_stress)
 post-process results for elastic solution More...
 

Set data from mesh to vectors and back

MoFEMErrorCode setFieldFromCoords (const std::string field_name)
 set field from node positions More...
 
MoFEMErrorCode setMaterialPositionFromCoords ()
 set material field from nodes More...
 
MoFEMErrorCode setSpatialPositionFromCoords ()
 set spatial field from nodes More...
 
MoFEMErrorCode setCoordsFromField (const std::string field_name="MESH_NODE_POSITIONS")
 set coordinates from field More...
 
MoFEMErrorCode setSingularDofs (const string field_name, const int verb=QUIET)
 set singular dofs (on edges adjacent to crack front) from geometry More...
 
MoFEMErrorCode setSingularElementMatrialPositions (const int verb=QUIET)
 set singular dofs of material field (on edges adjacent to crack front) from geometry More...
 
MoFEMErrorCode unsetSingularElementMatrialPositions ()
 remove singularity More...
 
MoFEMErrorCode cleanSingularElementMatrialPositions ()
 set maetrial field order to one More...
 

Auxiliary method

boost::shared_ptr< DofEntity > arcLengthDof
 
MoFEMErrorCode getArcLengthDof ()
 set pointer to arc-length DOF More...
 
boost::shared_ptr< ArcLengthCtx > & getArcCtx ()
 
boost::shared_ptr< ArcLengthCtx > & getEigenArcCtx ()
 
boost::shared_ptr< FEMethod > getFrontArcLengthControl (boost::shared_ptr< ArcLengthCtx > arc_ctx)
 
boost::shared_ptr< AnalyticalDirichletBC::DirichletBC > getAnalyticalDirichletBc ()
 
boost::shared_ptr< MWLSApprox > & getMWLSApprox ()
 

Pointers to finite element instances

bool onlyHooke
 True if only Hooke material is applied. More...
 
PetscBool onlyHookeFromOptions
 True if only Hooke material is applied. More...
 
boost::shared_ptr< NonlinearElasticElement > elasticFe
 
boost::shared_ptr< NonlinearElasticElement > materialFe
 
boost::shared_ptr< CrackFrontElement > feLhs
 Integrate elastic FE. More...
 
boost::shared_ptr< CrackFrontElement > feRhs
 Integrate elastic FE. More...
 
boost::shared_ptr< CrackFrontElement > feMaterialRhs
 Integrate material stresses, assemble vector. More...
 
boost::shared_ptr< CrackFrontElement > feMaterialLhs
 Integrate material stresses, assemble matrix. More...
 
boost::shared_ptr< CrackFrontElement > feEnergy
 Integrate energy. More...
 
boost::shared_ptr< ConstrainMatrixCtx > projSurfaceCtx
 Data structure to project on the body surface. More...
 
boost::shared_ptr< ConstrainMatrixCtx > projFrontCtx
 Data structure to project on crack front. More...
 
boost::shared_ptr< MWLSApprox > mwlsApprox
 
boost::shared_ptr< NeumannForcesSurface::MyTriangleFE > feMaterialAnaliticalTraction
 Surface elment to calculate tractions in material space. More...
 
boost::shared_ptr< DirichletSpatialPositionsBc > spatialDirichletBc
 apply Dirichlet BC to sparial positions More...
 
boost::shared_ptr< AnalyticalDirichletBC::DirichletBC > analyticalDirichletBc
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfaceForces
 assemble surface forces More...
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfaceForceAle
 assemble surface pressure (ALE) More...
 
boost::shared_ptr< NeumannForcesSurface::DataAtIntegrationPts > commonDataSurfaceForceAle
 common data at integration points (ALE) More...
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressure
 assemble surface pressure More...
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressureAle
 assemble surface pressure (ALE) More...
 
boost::shared_ptr< NeumannForcesSurface::DataAtIntegrationPts > commonDataSurfacePressureAle
 common data at integration points (ALE) More...
 
boost::shared_ptr< boost::ptr_map< string, EdgeForce > > edgeForces
 assemble edge forces More...
 
boost::shared_ptr< boost::ptr_map< string, NodalForce > > nodalForces
 assemble nodal forces More...
 
boost::shared_ptr< FEMethod > assembleFlambda
 assemble F_lambda vector More...
 
boost::shared_ptr< FEMethod > zeroFlambda
 assemble F_lambda vector More...
 
boost::shared_ptr< CrackFrontElement > feCouplingElasticLhs
 FE instance to assemble coupling terms. More...
 
boost::shared_ptr< CrackFrontElement > feCouplingMaterialLhs
 FE instance to assemble coupling terms. More...
 
boost::shared_ptr< Smoother > smootherFe
 
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherRhs
 Integrate smoothing operators. More...
 
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherLhs
 Integrate smoothing operators. More...
 
boost::shared_ptr< VolumeLengthQuality< double > > volumeLengthDouble
 
boost::shared_ptr< VolumeLengthQuality< adouble > > volumeLengthAdouble
 
boost::shared_ptr< ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain > tangentConstrains
 Constrains crack front in tangent directiona. More...
 
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > skinOrientation
 
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > crackOrientation
 
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > contactOrientation
 
map< int, boost::shared_ptr< SurfaceSlidingConstrains > > surfaceConstrain
 
map< int, boost::shared_ptr< EdgeSlidingConstrains > > edgeConstrains
 
boost::shared_ptr< BothSurfaceConstrains > bothSidesConstrains
 
boost::shared_ptr< BothSurfaceConstrains > closeCrackConstrains
 
boost::shared_ptr< BothSurfaceConstrains > bothSidesContactConstrains
 
boost::shared_ptr< DirichletFixFieldAtEntitiesBc > fixMaterialEnts
 
boost::shared_ptr< GriffithForceElement > griffithForceElement
 
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceRhs
 
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceLhs
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsDelta > feGriffithConstrainsDelta
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsRhs
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsLhs
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > analiticalSurfaceElement
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringLhsPtr
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringRhsPtr
 

Entity ranges

Range bitEnts
 
Range bitProcEnts
 
Range bodySkin
 
Range crackFaces
 
Range oneSideCrackFaces
 
Range otherSideCrackFaces
 
Range crackFront
 
Range crackFrontNodes
 
Range crackFrontNodesEdges
 
Range crackFrontNodesTris
 
Range crackFrontElements
 
EntityHandle arcLengthVertex
 Vertex associated with dof for arc-length control. More...
 

Auxiliary data

int defaultMaterial
 
int nbLoadSteps
 
int nbCutSteps
 
double loadScale
 
double maxG1
 
double maxJ
 
double fractionOfFixedNodes
 
map< EntityHandle, double > mapG1
 hashmap of g1 - release energy at nodes More...
 
map< EntityHandle, double > mapG3
 hashmap of g3 - release energy at nodes More...
 
map< EntityHandle, double > mapJ
 hashmap of J - release energy at nodes More...
 
map< EntityHandle, VectorDouble3 > mapMatForce
 hashmap of material force at nodes More...
 
map< EntityHandle, VectorDouble3 > mapGriffith
 hashmap of Griffith energy at nodes More...
 
map< EntityHandle, double > mapSmoothingForceFactor
 
bool doSurfaceProjection
 
double smootherAlpha
 Controls mesh smoothing. More...
 
double initialSmootherAlpha
 
double smootherGamma
 Controls mesh smoothing. More...
 

arc-length options

double arcAlpha
 
double arcBeta
 
double arcS
 
boost::shared_ptr< ArcLengthCtx > arcCtx
 
boost::shared_ptr< ArcLengthCtx > arcEigenCtx
 

Contact parameters, entity ranges and data structures

int contactOrder
 
int contactLambdaOrder
 
double rValue
 
double cnValue
 
PetscBool ignoreContact
 
PetscBool fixContactNodes
 
PetscBool printContactState
 
Range contactMeshsetFaces
 
Range contactElements
 
Range contactMasterFaces
 
Range contactSlaveFaces
 
Range mortarContactElements
 
Range mortarContactMasterFaces
 
Range mortarContactSlaveFaces
 
boost::shared_ptr< ContactSearchKdTree::ContactCommonData_multiIndex > contactSearchMultiIndexPtr
 
Range contactBothSidesMasterFaces
 
Range contactBothSidesSlaveFaces
 
Range contactBothSidesMasterNodes
 
Range contactTets
 
EntityHandle meshsetFaces
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > feLhsSpringALE
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > feLhsSpringALEMaterial
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > feRhsSpringALEMaterial
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > feLhsSurfaceForceALEMaterial
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > feRhsSurfaceForceALEMaterial
 
boost::shared_ptr< MetaSpringBC::DataAtIntegrationPtsSprings > commonDataSpringsALE
 
boost::shared_ptr< SimpleContactProblem > contactProblem
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContact
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContact
 
boost::shared_ptr< SimpleContactProblem::CommonDataSimpleContact > commonDataSimpleContact
 
boost::shared_ptr< MortarContactProblem > mortarContactProblemPtr
 
boost::shared_ptr< MortarContactProblem::MortarContactElement > feRhsMortarContact
 
boost::shared_ptr< MortarContactProblem::MortarContactElement > feLhsMortarContact
 
boost::shared_ptr< MortarContactProblem::CommonDataMortarContact > commonDataMortarContact
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContactALEMaterial
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContactALEMaterial
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContactALE
 
boost::shared_ptr< SimpleContactProblem::CommonDataSimpleContact > commonDataSimpleContactALE
 
PetscBool contactOutputIntegPts
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > fePostProcSimpleContact
 
boost::shared_ptr< MortarContactProblem::MortarContactElement > fePostProcMortarContact
 
moab::Core contactPostProcCore
 
moab::Interface & contactPostProcMoab
 
Range constrainedInterface
 Range of faces on the constrained interface. More...
 

Data for bone

double rHo0
 Reference density if bone is analyzed. More...
 
double nBone
 Exponent parameter in bone density. More...
 

Interfaces

const boost::scoped_ptr< UnknownInterface > cpSolversPtr
 
const boost::scoped_ptr< UnknownInterface > cpMeshCutPtr
 

Cut/Refine/Split options

PetscBool doCutMesh
 
double crackAccelerationFactor
 
PetscBool doElasticWithoutCrack
 

Detailed Description

Crack propagation data and methods

Definition at line 77 of file CrackPropagation.hpp.

Constructor & Destructor Documentation

◆ CrackPropagation()

FractureMechanics::CrackPropagation::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 used to solve fracture problem, i.e. CrackPropagation, CPSolvers and CPMeshCut.

Parameters
m_fieldinterface to MoAB database
approx_orderorder of approximation space, default is 2
geometry_orderorder of geometry, default is 1 (can be up to 2)

Definition at line 521 of file CrackPropagation.cpp.

524  : mField(m_field), contactPostProcMoab(contactPostProcCore),
525  setSingularCoordinates(false), approxOrder(approx_order),
526  geometryOrder(geometry_order), gC(0), betaGc(0),
527  isGcHeterogeneous(PETSC_FALSE), isPartitioned(PETSC_FALSE),
528  propagateCrack(PETSC_FALSE), refAtCrackTip(0), refOrderAtTip(0),
529  residualStressBlock(-1), densityMapBlock(-1),
530  addAnalyticalInternalStressOperators(PETSC_FALSE), postProcLevel(0),
531  startStep(0), crackFrontLength(-1), crackSurfaceArea(-1), rHo0(1),
532  nBone(0), cpSolversPtr(new CPSolvers(*this)),
533  cpMeshCutPtr(new CPMeshCut(*this)) {
534 
535  if (!LogManager::checkIfChannelExist("CPWorld")) {
536  auto core_log = logging::core::get();
537 
538  core_log->add_sink(
539  LogManager::createSink(LogManager::getStrmWorld(), "CPWorld"));
540  core_log->add_sink(
541  LogManager::createSink(LogManager::getStrmSync(), "CPSync"));
542  core_log->add_sink(
543  LogManager::createSink(LogManager::getStrmSelf(), "CPSelf"));
544 
545  LogManager::setLog("CPWorld");
546  LogManager::setLog("CPSync");
547  LogManager::setLog("CPSelf");
548 
549  MOFEM_LOG_TAG("CPWorld", "CP");
550  MOFEM_LOG_TAG("CPSync", "CP");
551  MOFEM_LOG_TAG("CPSelf", "CP");
552  }
553 
554  MOFEM_LOG("CPWorld", Sev::noisy) << "CPSolve created";
555 
556 #ifdef GIT_UM_SHA1_NAME
557  MOFEM_LOG_C("CPWorld", Sev::inform, "User module git commit id %s",
558  GIT_UM_SHA1_NAME);
559 #endif
560 
561  Version version;
562  getInterfaceVersion(version);
563  MOFEM_LOG("CPWorld", Sev::inform)
564  << "Fracture module version " << version.strVersion();
565 
566 #ifdef GIT_FM_SHA1_NAME
567  MOFEM_LOG_C("CPWorld", Sev::inform, "Fracture module git commit id %s",
568  GIT_FM_SHA1_NAME);
569 #endif
570 
571  ierr = registerInterface<CrackPropagation>();
572  CHKERRABORT(PETSC_COMM_SELF, ierr);
573  ierr = registerInterface<CPSolvers>();
574  CHKERRABORT(PETSC_COMM_SELF, ierr);
575  ierr = registerInterface<CPMeshCut>();
576  CHKERRABORT(PETSC_COMM_SELF, ierr);
577 
578  mapBitLevel["mesh_cut"] = BitRefLevel().set(0);
579  mapBitLevel["spatial_domain"] = BitRefLevel().set(1);
580  mapBitLevel["material_domain"] = BitRefLevel().set(2);
581 
582  if (!moabCommWorld)
583  moabCommWorld = boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
584 }

◆ ~CrackPropagation()

FractureMechanics::CrackPropagation::~CrackPropagation ( )
virtual

Definition at line 586 of file CrackPropagation.cpp.

586 {}

Member Function Documentation

◆ addElasticFEInstancesToSnes()

MoFEMErrorCode FractureMechanics::CrackPropagation::addElasticFEInstancesToSnes ( DM  dm,
Mat  m,
Vec  q,
Vec  f,
boost::shared_ptr< FEMethod >  arc_method = boost::shared_ptr<FEMethod>(),
boost::shared_ptr< ArcLengthCtx >  arc_ctx = nullptr,
const int  verb = QUIET,
const bool  debug = false 
)

Definition at line 4610 of file CrackPropagation.cpp.

4612  {
4613  boost::shared_ptr<FEMethod> null;
4614  MoFEMFunctionBegin;
4615 
4616  // Assemble F_lambda force
4617  assembleFlambda = boost::shared_ptr<FEMethod>(
4618  new AssembleFlambda(arc_ctx, spatialDirichletBc));
4619  zeroFlambda = boost::shared_ptr<FEMethod>(new ZeroFLmabda(arc_ctx));
4620 
4621  if (mwlsApprox) {
4622  mwlsApprox->F_lambda = arc_ctx->F_lambda;
4623  CHKERR getArcLengthDof();
4624  mwlsApprox->arcLengthDof = arcLengthDof;
4625  }
4626  // Set up DM specific vectors and data
4627  spatialDirichletBc->mapZeroRows.clear();
4628 
4629  VecSetOption(arc_ctx->F_lambda, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
4630 
4631  // Set-up F_lambda to elements & operator evalulating forces
4632  // Surface force
4633  for (auto fit = surfaceForces->begin(); fit != surfaceForces->end(); fit++) {
4634  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4635  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4636  for (; oit != hi_oit; oit++) {
4637  if (boost::typeindex::type_id_runtime(*oit) ==
4638  boost::typeindex::type_id<NeumannForcesSurface::OpNeumannForce>()) {
4639  dynamic_cast<NeumannForcesSurface::OpNeumannForce &>(*oit).F =
4640  arc_ctx->F_lambda;
4641  }
4642  }
4643  }
4644  for (auto fit = surfacePressure->begin(); fit != surfacePressure->end();
4645  fit++) {
4646  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4647  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4648  for (; oit != hi_oit; oit++) {
4649  if (boost::typeindex::type_id_runtime(*oit) ==
4650  boost::typeindex::type_id<
4651  NeumannForcesSurface::OpNeumannPressure>()) {
4652  dynamic_cast<NeumannForcesSurface::OpNeumannPressure &>(*oit).F =
4653  arc_ctx->F_lambda;
4654  }
4655  if (boost::typeindex::type_id_runtime(*oit) ==
4656  boost::typeindex::type_id<
4657  NeumannForcesSurface::OpNeumannForceAnalytical>()) {
4658  dynamic_cast<NeumannForcesSurface::OpNeumannForceAnalytical &>(*oit).F =
4659  arc_ctx->F_lambda;
4660  }
4661  }
4662  }
4663  for (auto fit = edgeForces->begin(); fit != edgeForces->end(); fit++) {
4664  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4665  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4666  for (; oit != hi_oit; oit++) {
4667  if (boost::typeindex::type_id_runtime(*oit) ==
4668  boost::typeindex::type_id<EdgeForce::OpEdgeForce>()) {
4669  dynamic_cast<EdgeForce::OpEdgeForce &>(*oit).F = arc_ctx->F_lambda;
4670  }
4671  }
4672  }
4673  for (auto fit = nodalForces->begin(); fit != nodalForces->end(); fit++) {
4674  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4675  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4676  for (; oit != hi_oit; oit++) {
4677  if (boost::typeindex::type_id_runtime(*oit) ==
4678  boost::typeindex::type_id<NodalForce::OpNodalForce>()) {
4679  dynamic_cast<NodalForce::OpNodalForce &>(*oit).F = arc_ctx->F_lambda;
4680  }
4681  }
4682  }
4683 
4684  // Lhs
4685  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, spatialDirichletBc,
4686  null);
4687 #ifdef __ANALITICAL_DISPLACEMENT__
4688  if (mField.check_problem("BC_PROBLEM")) {
4689  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null,
4690  analyticalDirichletBc, null);
4691  }
4692 #endif //__ANALITICAL_DISPLACEMENT__
4693  CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", feLhs, null, null);
4694  CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", feSpringLhsPtr, null, null);
4695  CHKERR DMMoFEMSNESSetJacobian(dm, "CONTACT", feLhsSimpleContact, null, null);
4696  CHKERR DMMoFEMSNESSetJacobian(dm, "MORTAR_CONTACT", feLhsMortarContact, null,
4697  null);
4698  {
4699  const MoFEM::Problem *problem_ptr;
4700  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
4701  if (problem_ptr->getName() == "EIGEN_ELASTIC") {
4702  CHKERR DMMoFEMSNESSetJacobian(dm, "CLOSE_CRACK", closeCrackConstrains,
4703  null, null);
4704  }
4705  }
4706 
4707 #ifdef __ANALITICAL_DISPLACEMENT__
4708  if (mField.check_problem("BC_PROBLEM")) {
4709  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, null,
4710  analyticalDirichletBc);
4711  }
4712 #endif //__ANALITICAL_DISPLACEMENT__
4713  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, null,
4714  spatialDirichletBc);
4715  CHKERR DMMoFEMSNESSetJacobian(dm, "ARC_LENGTH", arc_method, null, null);
4716 
4717  // Rhs
4718  auto fe_set_option = boost::make_shared<FEMethod>();
4719  fe_set_option->preProcessHook = [fe_set_option]() {
4720  return VecSetOption(fe_set_option->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
4721  PETSC_TRUE);
4722  };
4723  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, fe_set_option, null);
4724  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, spatialDirichletBc,
4725  null);
4726  CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", null, zeroFlambda, null);
4727 #ifdef __ANALITICAL_DISPLACEMENT__
4728  if (mField.check_problem("BC_PROBLEM")) {
4729  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null,
4730  analyticalDirichletBc, null);
4731  }
4732 #endif //__ANALITICAL_DISPLACEMENT__
4733  CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", feRhs, null, null);
4734  CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", feSpringRhsPtr, null, null);
4735  CHKERR DMMoFEMSNESSetFunction(dm, "CONTACT", feRhsSimpleContact, null, null);
4736  CHKERR DMMoFEMSNESSetFunction(dm, "MORTAR_CONTACT", feRhsMortarContact, null,
4737  null);
4738  {
4739  const MoFEM::Problem *problem_ptr;
4740  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
4741  if (problem_ptr->getName() == "EIGEN_ELASTIC") {
4742  CHKERR DMMoFEMSNESSetFunction(dm, "CLOSE_CRACK", closeCrackConstrains,
4743  null, null);
4744  }
4745  }
4746 
4747  for (auto fit = surfaceForces->begin(); fit != surfaceForces->end(); fit++) {
4748  CHKERR DMMoFEMSNESSetFunction(
4749  dm, fit->first.c_str(),
4750  boost::shared_ptr<FEMethod>(surfaceForces, &fit->second->getLoopFe()),
4751  null, null);
4752  }
4753  for (auto fit = surfacePressure->begin(); fit != surfacePressure->end();
4754  fit++) {
4755  CHKERR DMMoFEMSNESSetFunction(
4756  dm, fit->first.c_str(),
4757  boost::shared_ptr<FEMethod>(surfacePressure, &fit->second->getLoopFe()),
4758  null, null);
4759  }
4760  for (auto fit = edgeForces->begin(); fit != edgeForces->end(); fit++) {
4761  CHKERR DMMoFEMSNESSetFunction(
4762  dm, fit->first.c_str(),
4763  boost::shared_ptr<FEMethod>(edgeForces, &fit->second->getLoopFe()),
4764  null, null);
4765  }
4766  for (auto fit = nodalForces->begin(); fit != nodalForces->end(); fit++) {
4767  CHKERR DMMoFEMSNESSetFunction(
4768  dm, fit->first.c_str(),
4769  boost::shared_ptr<FEMethod>(nodalForces, &fit->second->getLoopFe()),
4770  null, null);
4771  }
4772 #ifdef __ANALITICAL_TRACTION__
4773  if (analiticalSurfaceElement) {
4774  for (auto fit = analiticalSurfaceElement->begin();
4775  fit != analiticalSurfaceElement->end(); fit++) {
4776  CHKERR DMMoFEMSNESSetFunction(
4777  dm, fit->first.c_str(),
4778  boost::shared_ptr<FEMethod>(analiticalSurfaceElement,
4779  &fit->second->fe),
4780  null, null);
4781  }
4782  }
4783 #endif //__ANALITICAL_TRACTION__
4784 #ifdef __ANALITICAL_DISPLACEMENT__
4785  if (mField.check_problem("BC_PROBLEM")) {
4786  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, null,
4787  analyticalDirichletBc);
4788  }
4789 #endif //__ANALITICAL_DISPLACEMENT__
4790  CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", assembleFlambda, null, null);
4791  CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", arc_method, null, null);
4792  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, null,
4793  spatialDirichletBc);
4794 
4795  if (debug) {
4796  if (m == PETSC_NULL || q == PETSC_NULL || f == PETSC_NULL) {
4797  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
4798  "problem matrix or vectors are set to null");
4799  }
4800  SNES snes;
4801  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
4802  SnesCtx *snes_ctx;
4803  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
4804  CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
4805  CHKERR SNESSetJacobian(snes, m, m, SnesMat, snes_ctx);
4806  CHKERR SnesMat(snes, q, m, m, snes_ctx);
4807  CHKERR SnesRhs(snes, q, f, snes_ctx);
4808  CHKERR SNESDestroy(&snes);
4809  MatView(m, PETSC_VIEWER_DRAW_WORLD);
4810  std::string wait;
4811  std::cin >> wait;
4812  }
4813 
4814  MoFEMFunctionReturn(0);
4815 }

◆ addMaterialFEInstancesToSnes()

MoFEMErrorCode FractureMechanics::CrackPropagation::addMaterialFEInstancesToSnes ( DM  dm,
const bool  fix_crack_front,
const int  verb = QUIET,
const bool  debug = false 
)

add material elements instances to SNES

Definition at line 6631 of file CrackPropagation.cpp.

6632  {
6633  boost::shared_ptr<FEMethod> null;
6634  MoFEMFunctionBegin;
6635 
6636  // Set up DM specific vectors and data
6637  Vec front_f, tangent_front_f;
6638  CHKERR DMCreateGlobalVector(dm, &front_f);
6639  CHKERR VecDuplicate(front_f, &tangent_front_f);
6640  CHKERR VecSetOption(front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6641  CHKERR VecSetOption(tangent_front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6642  if (smootherFe->smootherData.ownVectors) {
6643  if (smootherFe->smootherData.frontF != PETSC_NULL)
6644  CHKERR VecDestroy(&smootherFe->smootherData.frontF);
6645  if (smootherFe->smootherData.tangentFrontF != PETSC_NULL)
6646  CHKERR VecDestroy(&smootherFe->smootherData.tangentFrontF);
6647  }
6648  smootherFe->smootherData.frontF = front_f;
6649  smootherFe->smootherData.tangentFrontF = tangent_front_f;
6650  smootherFe->smootherData.ownVectors = true;
6651  CHKERR PetscObjectReference((PetscObject)front_f);
6652  CHKERR PetscObjectReference((PetscObject)tangent_front_f);
6653  if (tangentConstrains->ownVectors) {
6654  if (tangentConstrains->frontF != PETSC_NULL)
6655  VecDestroy(&tangentConstrains->frontF);
6656  if (tangentConstrains->tangentFrontF != PETSC_NULL)
6657  VecDestroy(&tangentConstrains->tangentFrontF);
6658  }
6659  tangentConstrains->frontF = front_f;
6660  tangentConstrains->tangentFrontF = tangent_front_f;
6661  tangentConstrains->ownVectors = true;
6662  CHKERR PetscObjectReference((PetscObject)front_f);
6663  CHKERR PetscObjectReference((PetscObject)tangent_front_f);
6664  CHKERR VecDestroy(&front_f);
6665  CHKERR VecDestroy(&tangent_front_f);
6666 
6667  const MoFEM::Problem *problem_ptr;
6668  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
6669  CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
6670  problem_ptr->getName(), COL, feGriffithConstrainsDelta->deltaVec);
6671  CHKERR VecSetOption(feGriffithConstrainsDelta->deltaVec,
6672  VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6673  feGriffithConstrainsRhs->deltaVec = feGriffithConstrainsDelta->deltaVec;
6674  feGriffithConstrainsLhs->deltaVec = feGriffithConstrainsDelta->deltaVec;
6675 
6676  // Fix nodes
6677  CHKERR updateMaterialFixedNode(true, false, debug);
6678 
6679  if (debug && !mField.get_comm_rank()) {
6680  const FieldEntity_multiIndex *field_ents;
6681  CHKERR mField.get_field_ents(&field_ents);
6682  for (Range::iterator vit = crackFrontNodes.begin();
6683  vit != crackFrontNodes.end(); ++vit) {
6684  cerr << "Node " << endl;
6685  FieldEntity_multiIndex::index<Ent_mi_tag>::type::iterator eit, hi_eit;
6686  eit = field_ents->get<Ent_mi_tag>().lower_bound(*vit);
6687  hi_eit = field_ents->get<Ent_mi_tag>().upper_bound(*vit);
6688  for (; eit != hi_eit; ++eit) {
6689  cerr << **eit << endl;
6690  }
6691  }
6692  }
6693 
6694  SnesCtx *snes_ctx;
6695  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6696  CHKERR snes_ctx->clearLoops();
6697 
6698  // Add to SNES Jacobian
6699  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, fixMaterialEnts, null);
6700  CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feSmootherLhs, null, null);
6701  CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feMaterialLhs, null, null);
6702  CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6703  null, null);
6704  CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CONTACT",
6705  bothSidesContactConstrains, null, null);
6706  for (auto &m : surfaceConstrain) {
6707  if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6708  CHKERR DMMoFEMSNESSetJacobian(dm, "SURFACE", m.second->feLhsPtr, null,
6709  null);
6710  } else {
6711  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACK_SURFACE", m.second->feLhsPtr,
6712  null, null);
6713  }
6714  }
6715  for (auto &m : edgeConstrains) {
6716  CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE", m.second->feLhsPtr, null, null);
6717  }
6718  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
6719  feGriffithConstrainsDelta, null, null);
6720  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
6721  feGriffithConstrainsLhs, null, null);
6722  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
6723  tangentConstrains, null, null);
6724  CHKERR DMMoFEMSNESSetJacobian(dm, "GRIFFITH_FORCE_ELEMENT",
6725  feGriffithForceLhs, null, null);
6726  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, null, fixMaterialEnts);
6727 
6728  // Add to SNES residuals
6729  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, fixMaterialEnts, null);
6730  CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feSmootherRhs, null, null);
6731  CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feMaterialRhs, null, null);
6732  CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6733  null, null);
6734  CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CONTACT",
6735  bothSidesContactConstrains, null, null);
6736  for (auto &m : surfaceConstrain) {
6737  if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6738  CHKERR DMMoFEMSNESSetFunction(dm, "SURFACE", m.second->feRhsPtr, null,
6739  null);
6740  } else {
6741  CHKERR DMMoFEMSNESSetFunction(dm, "CRACK_SURFACE", m.second->feRhsPtr,
6742  null, null);
6743  }
6744  }
6745  for (auto &m : edgeConstrains) {
6746  CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", m.second->feRhsPtr, null, null);
6747  }
6748  CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
6749  feGriffithConstrainsDelta, null, null);
6750  CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
6751  feGriffithConstrainsRhs, null, null);
6752  CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
6753  tangentConstrains, null, null);
6754  CHKERR DMMoFEMSNESSetFunction(dm, "GRIFFITH_FORCE_ELEMENT",
6755  feGriffithForceRhs, null, null);
6756  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, null, fixMaterialEnts);
6757 
6758  if (debug) {
6759  Vec q, f;
6760  CHKERR DMCreateGlobalVector(dm, &q);
6761  CHKERR VecDuplicate(q, &f);
6762  CHKERR DMoFEMMeshToLocalVector(dm, q, INSERT_VALUES, SCATTER_FORWARD);
6763  CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6764  Mat m;
6765  CHKERR DMCreateMatrix(dm, &m);
6766  SNES snes;
6767  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
6768  SnesCtx *snes_ctx;
6769  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6770  CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
6771  CHKERR SNESSetJacobian(snes, m, m, SnesMat, snes_ctx);
6772  CHKERR SnesMat(snes, q, m, m, snes_ctx);
6773  CHKERR SnesRhs(snes, q, f, snes_ctx);
6774  CHKERR SNESDestroy(&snes);
6775  // MatView(m,PETSC_VIEWER_DRAW_WORLD);
6776  MatView(m, PETSC_VIEWER_STDOUT_WORLD);
6777  std::string wait;
6778  std::cin >> wait;
6779  CHKERR VecDestroy(&q);
6780  CHKERR VecDestroy(&f);
6781  CHKERR MatDestroy(&m);
6782  }
6783  MoFEMFunctionReturn(0);
6784 }

◆ addMWLSDensityOperators()

MoFEMErrorCode FractureMechanics::CrackPropagation::addMWLSDensityOperators ( boost::shared_ptr< CrackFrontElement > &  fe_rhs,
boost::shared_ptr< CrackFrontElement > &  fe_lhs 
)

Definition at line 5027 of file CrackPropagation.cpp.

5029  {
5030  MoFEMFunctionBegin;
5031 
5032  if (densityMapBlock == -1)
5033  MoFEMFunctionReturnHot(0);
5034 
5035  if (!elasticFe) {
5036  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
5037  "Elastic element instance not created");
5038  }
5039 
5040  // Create elastic element finite element instance for residuals. Note that
5041  // this element approx. singularity at crack front.
5042  fe_rhs = boost::make_shared<CrackFrontElement>(
5043  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
5044  crackFrontElements, addSingularity);
5045  // Create finite element instance for assembly of tangent matrix
5046  fe_lhs = boost::make_shared<CrackFrontElement>(
5047  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
5048  crackFrontElements, addSingularity);
5049 
5050  bool test_with_mwls = false; // switch for debugging
5051 
5052  // Arbitrary lagrangian formulation, mesh node positions are taken into
5053  // account by operators.
5054  fe_rhs->meshPositionsFieldName = "NONE";
5055  fe_lhs->meshPositionsFieldName = "NONE";
5056  fe_rhs->addToRule = 0;
5057  fe_lhs->addToRule = 0;
5058 
5059  boost::shared_ptr<HookeElement::DataAtIntegrationPts>
5060  data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
5061  boost::shared_ptr<map<int, NonlinearElasticElement::BlockData>>
5062  block_sets_ptr(elasticFe, &elasticFe->setOfBlocks);
5063 
5064  if (!mwlsApprox) {
5065  // Create instance for moving least square approximation
5066  mwlsApprox = boost::make_shared<MWLSApprox>(mField);
5067  // Load mesh with stresses
5068  CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
5069  MeshsetsManager *block_manager_ptr;
5070  CHKERR mField.getInterface(block_manager_ptr);
5071  CHKERR block_manager_ptr->getEntitiesByDimension(
5072  densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
5073  Tag th_rho;
5074  if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(), th_rho) !=
5075  MB_SUCCESS) {
5076  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5077  "Density tag name %s cannot be found. Please "
5078  "provide the correct name.",
5079  mwlsRhoTagName.c_str());
5080  }
5081  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(), th_rho);
5082  CHKERR mwlsApprox->getValuesToNodes(th_rho);
5083  } else {
5084  mwlsApprox->tetsInBlock.clear();
5085  CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
5086  densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
5087  }
5088 
5089  auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
5090  fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5091  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5092  fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5093  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5094 
5095  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_MWLS_ptr(
5096  mwlsApprox, &mwlsApprox->mwlsMaterialPositions);
5097  fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5098  "MESH_NODE_POSITIONS", mat_pos_at_pts_MWLS_ptr));
5099  fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5100  "MESH_NODE_POSITIONS", mat_pos_at_pts_MWLS_ptr));
5101  auto mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
5102  fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
5103  "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
5104  fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
5105  "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
5106  fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
5107  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5108  fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
5109  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5110  fe_rhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
5111  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5112  fe_rhs->singularElement, fe_rhs->invSJac));
5113  fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
5114  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5115  fe_lhs->singularElement, fe_lhs->invSJac));
5116  auto space_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
5117  fe_rhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
5118  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_rhs->singularElement,
5119  fe_rhs->invSJac));
5120  fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
5121  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
5122  fe_lhs->invSJac));
5123 
5124  fe_rhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
5125  fe_rhs->singularElement, fe_rhs->detS, fe_rhs->invSJac));
5126 
5127  // data_hooke_element_at_pts->HMat = mat_grad_pos_at_pts_ptr;
5128  // data_hooke_element_at_pts->hMat = space_grad_pos_at_pts_ptr;
5129 
5130  // RHS //
5131  if (test_with_mwls) {
5132  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5133  fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
5134  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
5135  mwlsRhoTagName, true, false));
5136  else {
5137  fe_rhs->getOpPtrVector().push_back(
5138  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
5139  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox));
5140  fe_rhs->getOpPtrVector().push_back(
5141  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5142  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
5143  mwlsRhoTagName, true, false));
5144  }
5145  } else {
5146 
5147  fe_rhs->getOpPtrVector().push_back(new OpGetDensityFieldForTesting(
5148  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr, mwlsApprox->rhoAtGaussPts,
5149  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5150  fe_rhs, mwlsApprox->singularInitialDisplacement,
5151  mat_grad_pos_at_pts_ptr));
5152  }
5153 
5154  fe_rhs->getOpPtrVector().push_back(
5155  new HookeElement::OpCalculateStiffnessScaledByDensityField(
5156  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5157  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone, rHo0));
5158  fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStrainAle(
5159  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5160  fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStress<1>(
5161  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5162  fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEnergy(
5163  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5164 
5165  fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEshelbyStress(
5166  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5167 
5168  fe_rhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dX(
5169  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5170  fe_rhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dx(
5171  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
5172 
5173  fe_rhs->getOpPtrVector().push_back(new OpRhsBoneExplicitDerivariveWithHooke(
5174  *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5175  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->tetsInBlock, rHo0, nBone,
5176  &crackFrontNodes));
5177 
5178  // LHS //
5179  fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
5180  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
5181  fe_lhs->invSJac));
5182  fe_lhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
5183  fe_lhs->singularElement, fe_lhs->detS, fe_lhs->invSJac));
5184 
5185  if (test_with_mwls) {
5186  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5187  fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
5188  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
5189  mwlsRhoTagName, true, true));
5190  else {
5191  fe_lhs->getOpPtrVector().push_back(
5192  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
5193  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox));
5194  fe_lhs->getOpPtrVector().push_back(
5195  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5196  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
5197  mwlsRhoTagName, true, true));
5198  }
5199  } else {
5200  fe_lhs->getOpPtrVector().push_back(new OpGetDensityFieldForTesting(
5201  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr, mwlsApprox->rhoAtGaussPts,
5202  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5203  fe_lhs, mwlsApprox->singularInitialDisplacement,
5204  mat_grad_pos_at_pts_ptr));
5205  }
5206 
5207  fe_lhs->getOpPtrVector().push_back(
5208  new HookeElement::OpCalculateStiffnessScaledByDensityField(
5209  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5210  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone, rHo0));
5211 
5212  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStrainAle(
5213  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5214  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStress<1>(
5215  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5216  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEnergy(
5217  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5218  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEshelbyStress(
5219  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5220 
5221  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dX_dX<1>(
5222  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5223  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhsPre_dX_dx<1>(
5224  "MESH_NODE_POSITIONS", "SPATIAL_POSITION", data_hooke_element_at_pts));
5225  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dX_dx(
5226  "MESH_NODE_POSITIONS", "SPATIAL_POSITION", data_hooke_element_at_pts));
5227  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<1>(
5228  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
5229  fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dX<1>(
5230  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5231 
5232  fe_lhs->getOpPtrVector().push_back(
5233  new HookeElement::OpAleLhsWithDensity_dX_dX(
5234  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5235  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5236  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
5237 
5238  fe_lhs->getOpPtrVector().push_back(
5239  new HookeElement::OpAleLhsWithDensity_dx_dX(
5240  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", data_hooke_element_at_pts,
5241  mwlsApprox->rhoAtGaussPts, mwlsApprox->diffRhoAtGaussPts, nBone,
5242  rHo0));
5243 
5244  boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
5245  if (addSingularity) {
5246  mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
5247  mwlsApprox, &mwlsApprox->singularInitialDisplacement);
5248 
5249  fe_lhs->getOpPtrVector().push_back(
5250  new OpAleLhsWithDensitySingularElement_dX_dX(
5251  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5252  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5253  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0, mat_singular_disp_ptr));
5254 
5255  fe_lhs->getOpPtrVector().push_back(
5256  new OpAleLhsWithDensitySingularElement_dx_dX(
5257  "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
5258  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5259  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0, mat_singular_disp_ptr));
5260  }
5261 
5262  fe_lhs->getOpPtrVector().push_back(
5263  new OpLhsBoneExplicitDerivariveWithHooke_dX(
5264  *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5265  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5266  mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
5267  rHo0, nBone, &crackFrontNodes));
5268 
5269  fe_lhs->getOpPtrVector().push_back(
5270  new OpLhsBoneExplicitDerivariveWithHooke_dx(
5271  *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5272  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5273  mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
5274  rHo0, nBone, &crackFrontNodes));
5275 
5276  MoFEMFunctionReturn(0);
5277 }

◆ addMWLSStressOperators()

MoFEMErrorCode FractureMechanics::CrackPropagation::addMWLSStressOperators ( boost::shared_ptr< CrackFrontElement > &  fe_rhs,
boost::shared_ptr< CrackFrontElement > &  fe_lhs 
)

Definition at line 4902 of file CrackPropagation.cpp.

4904  {
4905  MoFEMFunctionBegin;
4906 
4907  if (residualStressBlock == -1)
4908  MoFEMFunctionReturnHot(0);
4909 
4910  if (!elasticFe) {
4911  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
4912  "Elastic element instance not created");
4913  }
4914 
4915  // Create elastic element finite element instance for residuals. Note that
4916  // this element approx. singularity at crack front.
4917  fe_rhs = boost::make_shared<CrackFrontElement>(
4918  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
4919  crackFrontElements, addSingularity);
4920  // Create finite element instance for assembly of tangent matrix
4921  fe_lhs = boost::make_shared<CrackFrontElement>(
4922  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
4923  crackFrontElements, addSingularity);
4924 
4925  // Arbitrary lagrangian formulation, mesh node positions are taken into
4926  // account by operators.
4927  fe_rhs->meshPositionsFieldName = "NONE";
4928  fe_lhs->meshPositionsFieldName = "NONE";
4929  fe_rhs->addToRule = 0;
4930  fe_lhs->addToRule = 0;
4931 
4932  if (!mwlsApprox) {
4933  // Create instance for moving least square approximation
4934  mwlsApprox = boost::make_shared<MWLSApprox>(mField);
4935  // Load mesh with stresses
4936  CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
4937  MeshsetsManager *block_manager_ptr;
4938  CHKERR mField.getInterface(block_manager_ptr);
4939  CHKERR block_manager_ptr->getEntitiesByDimension(
4940  residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
4941  Tag th;
4942  if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsStressTagName.c_str(), th) !=
4943  MB_SUCCESS) {
4944  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4945  "Internal stress tag name %s cannot be found. Please "
4946  "provide the correct name.",
4947  mwlsStressTagName.substr(4).c_str());
4948  }
4949  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsStressTagName.c_str(), th);
4950  CHKERR mwlsApprox->getValuesToNodes(th);
4951  } else {
4952  mwlsApprox->tetsInBlock.clear();
4953  CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
4954  residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
4955  }
4956 
4957  auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4958  fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4959  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4960  fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4961  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4962  boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr(new MatrixDouble());
4963  fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4964  "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4965  fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4966  "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4967 
4968  boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr(new MatrixDouble());
4969  fe_rhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4970  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_rhs->singularElement,
4971  fe_rhs->invSJac));
4972  fe_rhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4973  fe_rhs->singularElement, fe_rhs->detS, fe_rhs->invSJac));
4974 
4975  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4976  fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSStressAtGaussPts(
4977  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
4978  mwlsStressTagName, false, false));
4979  else {
4980  fe_rhs->getOpPtrVector().push_back(
4981  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
4982  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox));
4983  fe_rhs->getOpPtrVector().push_back(
4984  new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
4985  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
4986  mwlsStressTagName, false, false));
4987  }
4988 
4989  fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressRhs(
4990  mat_grad_pos_at_pts_ptr, mwlsApprox, false));
4991  fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSMaterialStressRhs(
4992  space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
4993  &crackFrontNodes));
4994 
4995  fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4996  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
4997  fe_lhs->invSJac));
4998  fe_lhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4999  fe_lhs->singularElement, fe_lhs->detS, fe_lhs->invSJac));
5000 
5001  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5002  fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSStressAtGaussPts(
5003  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
5004  mwlsStressTagName, true, false));
5005  else {
5006  fe_lhs->getOpPtrVector().push_back(
5007  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
5008  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox));
5009  fe_lhs->getOpPtrVector().push_back(
5010  new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
5011  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
5012  mwlsStressTagName, true, false));
5013  }
5014 
5015  fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressLhs_DX(
5016  mat_grad_pos_at_pts_ptr, mwlsApprox));
5017  fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSMaterialStressLhs_DX(
5018  space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
5019  &crackFrontNodes));
5020  fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSMaterialStressLhs_Dx(
5021  space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
5022  &crackFrontNodes));
5023 
5024  MoFEMFunctionReturn(0);
5025 }

◆ addPropagationFEInstancesToSnes()

MoFEMErrorCode FractureMechanics::CrackPropagation::addPropagationFEInstancesToSnes ( DM  dm,
boost::shared_ptr< FEMethod >  arc_method,
boost::shared_ptr< ArcLengthCtx >  arc_ctx,
const std::vector< int > &  surface_ids,
const int  verb = QUIET,
const bool  debug = false 
)

add finite element to SNES for crack propagation problem

Definition at line 6921 of file CrackPropagation.cpp.

6924  {
6925  boost::shared_ptr<FEMethod> null;
6926  MoFEMFunctionBegin;
6927 
6928  // Assemble F_lambda force
6929  assembleFlambda = boost::shared_ptr<FEMethod>(
6930  new AssembleFlambda(arc_ctx, spatialDirichletBc));
6931  zeroFlambda = boost::shared_ptr<FEMethod>(new ZeroFLmabda(arc_ctx));
6932 
6933  if (mwlsApprox) {
6934  mwlsApprox->F_lambda = arc_ctx->F_lambda;
6935  CHKERR getArcLengthDof();
6936  mwlsApprox->arcLengthDof = arcLengthDof;
6937  }
6938 
6939  // Set up DM specific vectors and data
6940  spatialDirichletBc->mapZeroRows.clear();
6941 
6942  // Set-up F_lambda to elements & operator evalulating forces
6943  // Surface fores
6944  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6945  surfaceForces->begin();
6946  fit != surfaceForces->end(); fit++) {
6947  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
6948  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6949  for (; oit != hi_oit; oit++) {
6950  if (boost::typeindex::type_id_runtime(*oit) ==
6951  boost::typeindex::type_id<NeumannForcesSurface::OpNeumannForce>()) {
6952  dynamic_cast<NeumannForcesSurface::OpNeumannForce &>(*oit).F =
6953  arc_ctx->F_lambda;
6954  }
6955  }
6956  }
6957  // Surface pressure (ALE)
6958  if (isSurfaceForceAle) {
6959  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6960  surfaceForceAle->begin();
6961  fit != surfaceForceAle->end(); fit++) {
6962 
6963  CHKERR getArcLengthDof();
6964  commonDataSurfaceForceAle->arcLengthDof = arcLengthDof;
6965 
6966  auto oit = fit->second->getLoopFeMatRhs().getOpPtrVector().begin();
6967  auto hi_oit = fit->second->getLoopFeMatRhs().getOpPtrVector().end();
6968  for (; oit != hi_oit; oit++) {
6969  if (boost::typeindex::type_id_runtime(*oit) ==
6970  boost::typeindex::type_id<
6971  NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX>()) {
6972  dynamic_cast<
6973  NeumannForcesSurface::OpNeumannSurfaceForceMaterialRhs_dX &>(*oit)
6974  .F = arc_ctx->F_lambda;
6975  }
6976  }
6977  }
6978  }
6979 
6980  // Surface pressure
6981  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6982  surfacePressure->begin();
6983  fit != surfacePressure->end(); fit++) {
6984 
6985  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
6986  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6987  for (; oit != hi_oit; oit++) {
6988  if (boost::typeindex::type_id_runtime(*oit) ==
6989  boost::typeindex::type_id<
6990  NeumannForcesSurface::OpNeumannPressure>()) {
6991  dynamic_cast<NeumannForcesSurface::OpNeumannPressure &>(*oit).F =
6992  arc_ctx->F_lambda;
6993  }
6994  if (boost::typeindex::type_id_runtime(*oit) ==
6995  boost::typeindex::type_id<
6996  NeumannForcesSurface::OpNeumannForceAnalytical>()) {
6997  dynamic_cast<NeumannForcesSurface::OpNeumannForceAnalytical &>(*oit).F =
6998  arc_ctx->F_lambda;
6999  }
7000  }
7001  }
7002  // Surface pressure (ALE)
7003  if (isPressureAle) {
7004  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7005  surfacePressureAle->begin();
7006  fit != surfacePressureAle->end(); fit++) {
7007 
7008  CHKERR getArcLengthDof();
7009  commonDataSurfacePressureAle->arcLengthDof = arcLengthDof;
7010 
7011  auto oit = fit->second->getLoopFeMatRhs().getOpPtrVector().begin();
7012  auto hi_oit = fit->second->getLoopFeMatRhs().getOpPtrVector().end();
7013  for (; oit != hi_oit; oit++) {
7014  if (boost::typeindex::type_id_runtime(*oit) ==
7015  boost::typeindex::type_id<
7016  NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX>()) {
7017  dynamic_cast<NeumannForcesSurface::OpNeumannPressureMaterialRhs_dX &>(
7018  *oit)
7019  .F = arc_ctx->F_lambda;
7020  }
7021  }
7022  }
7023  }
7024  // Surface edge
7025  for (boost::ptr_map<string, EdgeForce>::iterator fit = edgeForces->begin();
7026  fit != edgeForces->end(); fit++) {
7027  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
7028  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
7029  for (; oit != hi_oit; oit++) {
7030  if (boost::typeindex::type_id_runtime(*oit) ==
7031  boost::typeindex::type_id<EdgeForce::OpEdgeForce>()) {
7032  dynamic_cast<EdgeForce::OpEdgeForce &>(*oit).F = arc_ctx->F_lambda;
7033  }
7034  }
7035  }
7036  // Nodal force
7037  for (boost::ptr_map<string, NodalForce>::iterator fit = nodalForces->begin();
7038  fit != nodalForces->end(); fit++) {
7039  auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
7040  auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
7041  for (; oit != hi_oit; oit++) {
7042  if (boost::typeindex::type_id_runtime(*oit) ==
7043  boost::typeindex::type_id<NodalForce::OpNodalForce>()) {
7044  dynamic_cast<NodalForce::OpNodalForce &>(*oit).F = arc_ctx->F_lambda;
7045  }
7046  }
7047  }
7048  // Set up DM specific vectors and data
7049  Vec front_f, tangent_front_f;
7050  CHKERR DMCreateGlobalVector(dm, &front_f);
7051  CHKERR VecDuplicate(front_f, &tangent_front_f);
7052  CHKERR VecSetOption(front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7053  CHKERR VecSetOption(tangent_front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7054 
7055  {
7056  if (smootherFe->smootherData.ownVectors) {
7057  if (smootherFe->smootherData.frontF != PETSC_NULL)
7058  CHKERR VecDestroy(&smootherFe->smootherData.frontF);
7059  if (smootherFe->smootherData.tangentFrontF != PETSC_NULL)
7060  CHKERR VecDestroy(&smootherFe->smootherData.tangentFrontF);
7061  }
7062  smootherFe->smootherData.frontF = front_f;
7063  smootherFe->smootherData.tangentFrontF = tangent_front_f;
7064  smootherFe->smootherData.ownVectors = true;
7065 
7066  if (tangentConstrains->ownVectors) {
7067  if (tangentConstrains->frontF != PETSC_NULL)
7068  CHKERR VecDestroy(&tangentConstrains->frontF);
7069  if (tangentConstrains->tangentFrontF != PETSC_NULL)
7070  CHKERR VecDestroy(&tangentConstrains->tangentFrontF);
7071  }
7072  tangentConstrains->frontF = front_f;
7073  tangentConstrains->tangentFrontF = tangent_front_f;
7074  tangentConstrains->ownVectors = false;
7075  }
7076 
7077  const MoFEM::Problem *problem_ptr;
7078  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
7079  CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
7080  problem_ptr->getName(), COL, feGriffithConstrainsDelta->deltaVec);
7081  CHKERR VecSetOption(feGriffithConstrainsDelta->deltaVec,
7082  VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7083  feGriffithConstrainsRhs->deltaVec = feGriffithConstrainsDelta->deltaVec;
7084  feGriffithConstrainsLhs->deltaVec = feGriffithConstrainsDelta->deltaVec;
7085 
7086  // Fix nodes
7087  CHKERR updateMaterialFixedNode(false, true, debug);
7088  dynamic_cast<AssembleFlambda *>(assembleFlambda.get())
7089  ->pushDirichletBC(fixMaterialEnts);
7090 
7091  // Add to SNES Jacobian
7092  // Lhs
7093  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, fixMaterialEnts, null);
7094  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, spatialDirichletBc,
7095  null);
7096  CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC_NOT_ALE", feLhs, null, null);
7097  CHKERR DMMoFEMSNESSetJacobian(dm, "MATERIAL", feCouplingElasticLhs, null,
7098  null);
7099  CHKERR DMMoFEMSNESSetJacobian(dm, "MATERIAL", feCouplingMaterialLhs, null,
7100  null);
7101  CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", feSpringLhsPtr, null, null);
7102  CHKERR DMMoFEMSNESSetJacobian(dm, "CONTACT", feLhsSimpleContact, null, null);
7103  CHKERR DMMoFEMSNESSetJacobian(dm, "MORTAR_CONTACT", feLhsMortarContact, null,
7104  null);
7105  CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
7106  null, null);
7107  CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CONTACT",
7108  bothSidesContactConstrains, null, null);
7109  for (auto &m : surfaceConstrain) {
7110  if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
7111  CHKERR DMMoFEMSNESSetJacobian(dm, "SURFACE", m.second->feLhsPtr, null,
7112  null);
7113  } else {
7114  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACK_SURFACE", m.second->feLhsPtr,
7115  null, null);
7116  }
7117  }
7118  for (auto &m : edgeConstrains)
7119  CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE", m.second->feLhsPtr, null, null);
7120 
7121  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
7122  feGriffithConstrainsDelta, null, null);
7123  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
7124  feGriffithConstrainsLhs, null, null);
7125  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
7126  tangentConstrains, null, null);
7127  CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feSmootherLhs, null, null);
7128  CHKERR DMMoFEMSNESSetJacobian(dm, "GRIFFITH_FORCE_ELEMENT",
7129  feGriffithForceLhs, null, null);
7130  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, null,
7131  spatialDirichletBc);
7132  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, null, fixMaterialEnts);
7133  CHKERR DMMoFEMSNESSetJacobian(dm, "ARC_LENGTH", arc_method, null, null);
7134 
7135  // Rhs
7136  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, fixMaterialEnts, null);
7137  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, spatialDirichletBc,
7138  null);
7139  CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", null, zeroFlambda, null);
7140  CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", feRhs, null, null);
7141  CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", feSpringRhsPtr, null, null);
7142  CHKERR DMMoFEMSNESSetFunction(dm, "CONTACT", feRhsSimpleContact, null, null);
7143  CHKERR DMMoFEMSNESSetFunction(dm, "MORTAR_CONTACT", feRhsMortarContact, null,
7144  null);
7145  CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
7146  null, null);
7147  CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CONTACT",
7148  bothSidesContactConstrains, null, null);
7149  for (auto &m : surfaceConstrain) {
7150  if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
7151  CHKERR DMMoFEMSNESSetFunction(dm, "SURFACE", m.second->feRhsPtr, null,
7152  null);
7153  } else {
7154  CHKERR DMMoFEMSNESSetFunction(dm, "CRACK_SURFACE", m.second->feRhsPtr,
7155  null, null);
7156  }
7157  }
7158  for (auto &m : edgeConstrains) {
7159  CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", m.second->feRhsPtr, null, null);
7160  }
7161  CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feSmootherRhs, null, null);
7162  CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
7163  feGriffithConstrainsDelta, null, null);
7164  CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
7165  feGriffithConstrainsRhs, null, null);
7166  CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
7167  tangentConstrains, null, null);
7168  CHKERR DMMoFEMSNESSetFunction(dm, "MATERIAL", feMaterialRhs, null, null);
7169  CHKERR DMMoFEMSNESSetFunction(dm, "GRIFFITH_FORCE_ELEMENT",
7170  feGriffithForceRhs, null, null);
7171  // Add force elements
7172  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7173  surfaceForces->begin();
7174  fit != surfaceForces->end(); fit++) {
7175  CHKERR DMMoFEMSNESSetFunction(
7176  dm, fit->first.c_str(),
7177  boost::shared_ptr<FEMethod>(surfaceForces, &fit->second->getLoopFe()),
7178  null, null);
7179  }
7180 
7181  // Add surface force elements (ALE)
7182  if (isSurfaceForceAle) {
7183  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7184  surfaceForceAle->begin();
7185  fit != surfaceForceAle->end(); fit++) {
7186  CHKERR DMMoFEMSNESSetJacobian(
7187  dm, fit->first.c_str(),
7188  boost::shared_ptr<FEMethod>(surfaceForceAle,
7189  &fit->second->getLoopFeLhs()),
7190  null, null);
7191  if (!ignoreMaterialForce) {
7192  CHKERR DMMoFEMSNESSetFunction(
7193  dm, fit->first.c_str(),
7194  boost::shared_ptr<FEMethod>(surfaceForceAle,
7195  &fit->second->getLoopFeMatRhs()),
7196  null, null);
7197  CHKERR DMMoFEMSNESSetJacobian(
7198  dm, fit->first.c_str(),
7199  boost::shared_ptr<FEMethod>(surfaceForceAle,
7200  &fit->second->getLoopFeMatLhs()),
7201  null, null);
7202  }
7203  }
7204  }
7205 
7206  // Add pressure elements
7207  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7208  surfacePressure->begin();
7209  fit != surfacePressure->end(); fit++) {
7210  CHKERR DMMoFEMSNESSetFunction(
7211  dm, fit->first.c_str(),
7212  boost::shared_ptr<FEMethod>(surfacePressure, &fit->second->getLoopFe()),
7213  null, null);
7214  }
7215 
7216  if (isPressureAle) {
7217  for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7218  surfacePressureAle->begin();
7219  fit != surfacePressureAle->end(); fit++) {
7220  CHKERR DMMoFEMSNESSetJacobian(
7221  dm, fit->first.c_str(),
7222  boost::shared_ptr<FEMethod>(surfacePressureAle,
7223  &fit->second->getLoopFeLhs()),
7224  null, null);
7225  CHKERR DMMoFEMSNESSetFunction(
7226  dm, fit->first.c_str(),
7227  boost::shared_ptr<FEMethod>(surfacePressureAle,
7228  &fit->second->getLoopFeMatRhs()),
7229  null, null);
7230  CHKERR DMMoFEMSNESSetJacobian(
7231  dm, fit->first.c_str(),
7232  boost::shared_ptr<FEMethod>(surfacePressureAle,
7233  &fit->second->getLoopFeMatLhs()),
7234  null, null);
7235  }
7236  }
7237 
7238  if (areSpringsAle) {
7239  CHKERR DMMoFEMSNESSetFunction(
7240  dm, "SPRING_ALE", feRhsSpringALEMaterial.get(), PETSC_NULL, PETSC_NULL);
7241  CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING_ALE", feLhsSpringALE.get(), NULL,
7242  NULL);
7243  CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING_ALE",
7244  feLhsSpringALEMaterial.get(), NULL, NULL);
7245  }
7246 
7247  if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
7248  CHKERR DMMoFEMSNESSetFunction(dm, "SIMPLE_CONTACT_ALE",
7249  feRhsSimpleContactALEMaterial.get(),
7250  PETSC_NULL, PETSC_NULL);
7251  CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ALE",
7252  feLhsSimpleContactALEMaterial.get(), NULL,
7253  NULL);
7254  CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ALE",
7255  feLhsSimpleContactALE.get(), NULL, NULL);
7256  }
7257  for (boost::ptr_map<string, EdgeForce>::iterator fit = edgeForces->begin();
7258  fit != edgeForces->end(); fit++) {
7259  CHKERR DMMoFEMSNESSetFunction(
7260  dm, fit->first.c_str(),
7261  boost::shared_ptr<FEMethod>(edgeForces, &fit->second->getLoopFe()),
7262  null, null);
7263  }
7264  for (boost::ptr_map<string, NodalForce>::iterator fit = nodalForces->begin();
7265  fit != nodalForces->end(); fit++) {
7266  CHKERR DMMoFEMSNESSetFunction(
7267  dm, fit->first.c_str(),
7268  boost::shared_ptr<FEMethod>(nodalForces, &fit->second->getLoopFe()),
7269  null, null);
7270  }
7271  CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", assembleFlambda, null, null);
7272  CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", arc_method, null, null);
7273  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, null,
7274  spatialDirichletBc);
7275  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, null, fixMaterialEnts);
7276 
7277  MoFEMFunctionReturn(0);
7278 }

◆ addSmoothingFEInstancesToSnes()

MoFEMErrorCode FractureMechanics::CrackPropagation::addSmoothingFEInstancesToSnes ( DM  dm,
const bool  fix_crack_front,
const int  verb = QUIET,
const bool  debug = false 
)

add softening elements instances to SNES

Definition at line 6786 of file CrackPropagation.cpp.

6787  {
6788  boost::shared_ptr<FEMethod> null;
6789  MoFEMFunctionBegin;
6790 
6791  // Set up DM specific vectors and data
6792  if (smootherFe->smootherData.ownVectors) {
6793  if (smootherFe->smootherData.frontF != PETSC_NULL)
6794  CHKERR VecDestroy(&smootherFe->smootherData.frontF);
6795  if (smootherFe->smootherData.tangentFrontF != PETSC_NULL)
6796  CHKERR VecDestroy(&smootherFe->smootherData.tangentFrontF);
6797  }
6798  smootherFe->smootherData.frontF = PETSC_NULL;
6799  smootherFe->smootherData.tangentFrontF = PETSC_NULL;
6800  smootherFe->smootherData.ownVectors = false;
6801 
6802  // Fix nodes
6803  CHKERR updateMaterialFixedNode(true, false, debug);
6804 
6805  if (debug && !mField.get_comm_rank()) {
6806  const FieldEntity_multiIndex *field_ents;
6807  CHKERR mField.get_field_ents(&field_ents);
6808  for (Range::iterator vit = crackFrontNodes.begin();
6809  vit != crackFrontNodes.end(); ++vit) {
6810  cerr << "Node " << endl;
6811  FieldEntity_multiIndex::index<Ent_mi_tag>::type::iterator eit, hi_eit;
6812  eit = field_ents->get<Ent_mi_tag>().lower_bound(*vit);
6813  hi_eit = field_ents->get<Ent_mi_tag>().upper_bound(*vit);
6814  for (; eit != hi_eit; ++eit) {
6815  cerr << **eit << endl;
6816  }
6817  }
6818  }
6819 
6820  SnesCtx *snes_ctx;
6821  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6822  CHKERR snes_ctx->clearLoops();
6823 
6824  // Add to SNES Jacobian
6825  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, fixMaterialEnts, null);
6826  CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6827  null, null);
6828  CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CONTACT",
6829  bothSidesContactConstrains, null, null);
6830  for (auto &m : surfaceConstrain) {
6831  if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6832  CHKERR DMMoFEMSNESSetJacobian(dm, "SURFACE", m.second->feLhsPtr, null,
6833  null);
6834  } else {
6835  CHKERR DMMoFEMSNESSetJacobian(dm, "CRACK_SURFACE", m.second->feLhsPtr,
6836  null, null);
6837  }
6838  }
6839  for (auto &m : edgeConstrains) {
6840  CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE", m.second->feLhsPtr, null, null);
6841  }
6842  CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feSmootherLhs, null, null);
6843  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, null, fixMaterialEnts);
6844 
6845  // Add to SNES residuals
6846  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, fixMaterialEnts, null);
6847  CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6848  null, null);
6849  CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CONTACT",
6850  bothSidesContactConstrains, null, null);
6851  for (auto &m : surfaceConstrain) {
6852  if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6853  CHKERR DMMoFEMSNESSetFunction(dm, "SURFACE", m.second->feRhsPtr, null,
6854  null);
6855  } else {
6856  CHKERR DMMoFEMSNESSetFunction(dm, "CRACK_SURFACE", m.second->feRhsPtr,
6857  null, null);
6858  }
6859  }
6860  for (auto &m : edgeConstrains) {
6861  CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", m.second->feRhsPtr, null, null);
6862  }
6863 
6864  CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feSmootherRhs, null, null);
6865  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, null, fixMaterialEnts);
6866  for (auto &m : edgeConstrains) {
6867  CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", null, null, m.second->feRhsPtr);
6868  }
6869 
6870  if (debug) {
6871  Vec q, f;
6872  CHKERR DMCreateGlobalVector(dm, &q);
6873  CHKERR VecDuplicate(q, &f);
6874  CHKERR DMoFEMMeshToLocalVector(dm, q, INSERT_VALUES, SCATTER_FORWARD);
6875  CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6876  Mat m;
6877  CHKERR DMCreateMatrix(dm, &m);
6878  SNES snes;
6879  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
6880  SnesCtx *snes_ctx;
6881  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6882  CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
6883  CHKERR SNESSetJacobian(snes, m, m, SnesMat, snes_ctx);
6884  CHKERR SnesMat(snes, q, m, m, snes_ctx);
6885  CHKERR SnesRhs(snes, q, f, snes_ctx);
6886  CHKERR SNESDestroy(&snes);
6887  // MatView(m,PETSC_VIEWER_DRAW_WORLD);
6888  MatView(m, PETSC_VIEWER_STDOUT_WORLD);
6889  std::string wait;
6890  std::cin >> wait;
6891  CHKERR VecDestroy(&q);
6892  CHKERR VecDestroy(&f);
6893  CHKERR MatDestroy(&m);
6894  }
6895  MoFEMFunctionReturn(0);
6896 }

◆ analyticalStrainFunction()

FTensor::Tensor2_symmetric< double, 3 > FractureMechanics::CrackPropagation::analyticalStrainFunction ( FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &  t_coords)
static

Definition at line 10030 of file CrackPropagation.cpp.

10034  {
10035  FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
10036  constexpr double alpha = 1.e-5;
10037  FTensor::Index<'i', 3> i;
10038  FTensor::Index<'j', 3> j;
10039  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
10040  // FIXME put here formula from test
10041  double temp = 250.;
10042  double z = t_coords(2);
10043  if ((-10. < z && z < -1.) || std::abs(z + 1.) < 1e-15) {
10044  temp = 10. / 3. * (35. - 4. * z);
10045  }
10046  if ((-1. < z && z < 2.) || std::abs(z - 2.) < 1e-15) {
10047  temp = 10. / 3. * (34. - 5. * z);
10048  }
10049  if (2. < z && z < 10.) {
10050  temp = 5. / 4. * (30. + 17. * z);
10051  }
10052 
10053  t_thermal_strain(i, j) = alpha * (temp - 250.) * t_kd(i, j);
10054  return t_thermal_strain;
10055 }

◆ assembleCouplingForcesDM()

MoFEMErrorCode FractureMechanics::CrackPropagation::assembleCouplingForcesDM ( DM  dm,
const int  verb = QUIET,
const bool  debug = false 
)

assemble coupling element instances

Definition at line 6178 of file CrackPropagation.cpp.

6179  {
6180  MoFEMFunctionBeginHot;
6181 
6182  if (!elasticFe) {
6183  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
6184  "Elastic element instance not created");
6185  }
6186  feCouplingElasticLhs = boost::make_shared<CrackFrontElement>(
6187  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
6188  crackFrontElements, addSingularity);
6189  feCouplingElasticLhs->meshPositionsFieldName = "NONE";
6190  feCouplingElasticLhs->addToRule = 0;
6191 
6192  using BlockData = NonlinearElasticElement::BlockData;
6193  boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
6194  elasticFe, &elasticFe->setOfBlocks);
6195 
6196  // Calculate spatial positions and gradients (of deformation) at integration
6197  // pts.
6198  // This operator takes into account singularity at crack front
6199  if (!onlyHooke) {
6200  feCouplingElasticLhs->getOpPtrVector().push_back(
6201  new OpGetCrackFrontCommonDataAtGaussPts(
6202  "SPATIAL_POSITION", elasticFe->commonData,
6203  feCouplingElasticLhs->singularElement,
6204  feCouplingElasticLhs->invSJac));
6205  // Calculate material positions and gradients at integration pts.
6206  feCouplingElasticLhs->getOpPtrVector().push_back(
6207  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
6208  "MESH_NODE_POSITIONS", elasticFe->commonData));
6209  // Transform base functions to get singular base functions at crack front
6210  feCouplingElasticLhs->getOpPtrVector().push_back(
6211  new OpTransfromSingularBaseFunctions(
6212  feCouplingElasticLhs->singularElement, feCouplingElasticLhs->detS,
6213  feCouplingElasticLhs->invSJac));
6214  // Iterate over material blocks
6215  map<int, NonlinearElasticElement::BlockData>::iterator sit =
6216  elasticFe->setOfBlocks.begin();
6217  for (; sit != elasticFe->setOfBlocks.end(); sit++) {
6218  // Evaluate stress at integration pts
6219  feCouplingElasticLhs->getOpPtrVector().push_back(
6220  new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
6221  "SPATIAL_POSITION", sit->second, elasticFe->commonData,
6222  ELASTIC_TAG, true, true, false));
6223  // Assemble tangent matrix, derivative of spatial positions
6224  feCouplingElasticLhs->getOpPtrVector().push_back(
6225  new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(
6226  "SPATIAL_POSITION", "SPATIAL_POSITION", sit->second,
6227  elasticFe->commonData));
6228  // Assemble tangent matrix, derivative of material positions
6229  feCouplingElasticLhs->getOpPtrVector().push_back(
6230  new NonlinearElasticElement::OpLhsPiolaKirchhoff_dX(
6231  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", sit->second,
6232  elasticFe->commonData));
6233  }
6234  } else {
6235  boost::shared_ptr<HookeElement::DataAtIntegrationPts>
6236  data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
6237  feCouplingElasticLhs->getOpPtrVector().push_back(
6238  new OpCalculateVectorFieldGradient<3, 3>(
6239  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
6240 
6241  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(new MatrixDouble());
6242  feCouplingElasticLhs->getOpPtrVector().push_back(
6243  new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
6244  mat_pos_at_pts_ptr));
6245 
6246  if (defaultMaterial == BONEHOOKE) {
6247  if (!mwlsApprox) {
6248  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
6249  "mwlsApprox not allocated");
6250  }
6251  // calculate position at integration points
6252  boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
6253  mat_grad_pos_at_pts_ptr =
6254  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
6255  feCouplingElasticLhs->getOpPtrVector().push_back(
6256  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
6257  mat_grad_pos_at_pts_ptr));
6258  feCouplingElasticLhs->getOpPtrVector().push_back(
6259  new OpGetCrackFrontDataGradientAtGaussPts(
6260  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6261  feCouplingElasticLhs->singularElement,
6262  feCouplingElasticLhs->invSJac));
6263 
6264  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6265  feCouplingElasticLhs->getOpPtrVector().push_back(
6266  new MWLSApprox::OpMWLSRhoAtGaussPts(
6267  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6268  feCouplingElasticLhs, mwlsApprox, mwlsRhoTagName, true, false));
6269  else {
6270  feCouplingElasticLhs->getOpPtrVector().push_back(
6271  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
6272  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6273  feCouplingElasticLhs, mwlsApprox));
6274  feCouplingElasticLhs->getOpPtrVector().push_back(
6275  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
6276  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6277  feCouplingElasticLhs, mwlsApprox, mwlsRhoTagName, true, false));
6278  }
6279  feCouplingElasticLhs->getOpPtrVector().push_back(
6280  new HookeElement::OpCalculateStiffnessScaledByDensityField(
6281  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6282  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
6283  rHo0));
6284  // Calculate spatial positions and gradients (of deformation) at
6285  // integration pts. This operator takes into account singularity at crack
6286  // front
6287  feCouplingElasticLhs->getOpPtrVector().push_back(
6288  new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
6289  "SPATIAL_POSITION",
6290  data_hooke_element_at_pts));
6291  feCouplingElasticLhs->getOpPtrVector().push_back(
6292  new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
6293  "SPATIAL_POSITION",
6294  data_hooke_element_at_pts));
6295  // Transform base functions to get singular base functions at crack front
6296  feCouplingElasticLhs->getOpPtrVector().push_back(
6297  new OpTransfromSingularBaseFunctions(
6298  feCouplingElasticLhs->singularElement, feCouplingElasticLhs->detS,
6299  feCouplingElasticLhs->invSJac));
6300  feCouplingElasticLhs->getOpPtrVector().push_back(
6301  new HookeElement::OpAleLhs_dx_dx<1>("SPATIAL_POSITION",
6302  "SPATIAL_POSITION",
6303  data_hooke_element_at_pts));
6304  feCouplingElasticLhs->getOpPtrVector().push_back(
6305  new HookeElement::OpAleLhs_dx_dX<1>("SPATIAL_POSITION",
6306  "MESH_NODE_POSITIONS",
6307  data_hooke_element_at_pts));
6308 
6309  feCouplingElasticLhs->getOpPtrVector().push_back(
6310  new HookeElement::OpAleLhsWithDensity_dx_dX(
6311  "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
6312  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6313  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
6314 
6315  boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
6316  if (addSingularity) {
6317 
6318  mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
6319  mwlsApprox, &mwlsApprox->singularInitialDisplacement);
6320 
6321  feCouplingElasticLhs->getOpPtrVector().push_back(
6322  new OpAleLhsWithDensitySingularElement_dx_dX(
6323  "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
6324  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6325  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0,
6326  mat_singular_disp_ptr));
6327  }
6328 
6329  } else {
6330 
6331  feCouplingElasticLhs->getOpPtrVector().push_back(
6332  new HookeElement::OpCalculateHomogeneousStiffness<0>(
6333  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6334  data_hooke_element_at_pts));
6335  // Calculate spatial positions and gradients (of deformation) at
6336  // integration pts. This operator takes into account singularity at crack
6337  // front
6338  feCouplingElasticLhs->getOpPtrVector().push_back(
6339  new OpGetCrackFrontDataGradientAtGaussPts(
6340  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6341  feCouplingElasticLhs->singularElement,
6342  feCouplingElasticLhs->invSJac));
6343  feCouplingElasticLhs->getOpPtrVector().push_back(
6344  new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
6345  "SPATIAL_POSITION",
6346  data_hooke_element_at_pts));
6347  feCouplingElasticLhs->getOpPtrVector().push_back(
6348  new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
6349  "SPATIAL_POSITION",
6350  data_hooke_element_at_pts));
6351  // Transform base functions to get singular base functions at crack front
6352  feCouplingElasticLhs->getOpPtrVector().push_back(
6353  new OpTransfromSingularBaseFunctions(
6354  feCouplingElasticLhs->singularElement, feCouplingElasticLhs->detS,
6355  feCouplingElasticLhs->invSJac));
6356  feCouplingElasticLhs->getOpPtrVector().push_back(
6357  new HookeElement::OpAleLhs_dx_dx<0>("SPATIAL_POSITION",
6358  "SPATIAL_POSITION",
6359  data_hooke_element_at_pts));
6360  feCouplingElasticLhs->getOpPtrVector().push_back(
6361  new HookeElement::OpAleLhs_dx_dX<0>("SPATIAL_POSITION",
6362  "MESH_NODE_POSITIONS",
6363  data_hooke_element_at_pts));
6364  }
6365  }
6366 
6367  if (!materialFe) {
6368  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
6369  "Material element instance not created");
6370  }
6371  feCouplingMaterialLhs = boost::make_shared<CrackFrontElement>(
6372  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
6373  crackFrontElements, addSingularity);
6374  feCouplingMaterialLhs->meshPositionsFieldName = "NONE";
6375  feCouplingMaterialLhs->addToRule = 0;
6376 
6377  // calculate positions at integration points
6378  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(new MatrixDouble());
6379  feCouplingMaterialLhs->getOpPtrVector().push_back(
6380  new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
6381  mat_pos_at_pts_ptr));
6382  boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
6383  if (!onlyHooke && residualStressBlock != -1) {
6384  mat_grad_pos_at_pts_ptr =
6385  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
6386  feCouplingMaterialLhs->getOpPtrVector().push_back(
6387  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
6388  mat_grad_pos_at_pts_ptr));
6389  }
6390  boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr;
6391  if (!onlyHooke && residualStressBlock != -1) {
6392  space_grad_pos_at_pts_ptr =
6393  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
6394  feCouplingMaterialLhs->getOpPtrVector().push_back(
6395  new OpGetCrackFrontDataGradientAtGaussPts(
6396  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
6397  feCouplingMaterialLhs->singularElement,
6398  feCouplingMaterialLhs->invSJac));
6399  }
6400 
6401  // assemble tangent matrix
6402  if (!onlyHooke) {
6403  feCouplingMaterialLhs->getOpPtrVector().push_back(
6404  new OpGetCrackFrontCommonDataAtGaussPts(
6405  "SPATIAL_POSITION", materialFe->commonData,
6406  feCouplingMaterialLhs->singularElement,
6407  feCouplingMaterialLhs->invSJac));
6408  feCouplingMaterialLhs->getOpPtrVector().push_back(
6409  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
6410  "MESH_NODE_POSITIONS", materialFe->commonData));
6411  feCouplingMaterialLhs->getOpPtrVector().push_back(
6412  new OpTransfromSingularBaseFunctions(
6413  feCouplingMaterialLhs->singularElement, feCouplingMaterialLhs->detS,
6414  feCouplingMaterialLhs->invSJac));
6415  for (map<int, NonlinearElasticElement::BlockData>::iterator sit =
6416  materialFe->setOfBlocks.begin();
6417  sit != materialFe->setOfBlocks.end(); sit++) {
6418  feCouplingMaterialLhs->getOpPtrVector().push_back(
6419  new NonlinearElasticElement::OpJacobianEshelbyStress(
6420  "SPATIAL_POSITION", sit->second, materialFe->commonData,
6421  MATERIAL_TAG, true, true));
6422  feCouplingMaterialLhs->getOpPtrVector().push_back(
6423  new NonlinearElasticElement::OpLhsEshelby_dX(
6424  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", sit->second,
6425  materialFe->commonData));
6426  feCouplingMaterialLhs->getOpPtrVector().push_back(
6427  new NonlinearElasticElement::OpLhsEshelby_dx(
6428  "MESH_NODE_POSITIONS", "SPATIAL_POSITION", sit->second,
6429  materialFe->commonData));
6430  }
6431  } else {
6432 
6433  boost::shared_ptr<HookeElement::DataAtIntegrationPts>
6434  data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
6435  data_hooke_element_at_pts->forcesOnlyOnEntitiesRow = crackFrontNodes;
6436 
6437  feCouplingMaterialLhs->getOpPtrVector().push_back(
6438  new OpCalculateVectorFieldGradient<3, 3>(
6439  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
6440  mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
6441 
6442  if (defaultMaterial == BONEHOOKE) {
6443  if (!mwlsApprox) {
6444  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
6445  "mwlsApprox not allocated");
6446  }
6447 
6448  feCouplingMaterialLhs->getOpPtrVector().push_back(
6449  new OpGetCrackFrontDataGradientAtGaussPts(
6450  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6451  feCouplingMaterialLhs->singularElement,
6452  feCouplingMaterialLhs->invSJac));
6453 
6454  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6455  feCouplingMaterialLhs->getOpPtrVector().push_back(
6456  new MWLSApprox::OpMWLSRhoAtGaussPts(
6457  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6458  feCouplingMaterialLhs, mwlsApprox, mwlsRhoTagName, true, true));
6459  else {
6460  feCouplingMaterialLhs->getOpPtrVector().push_back(
6461  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
6462  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6463  feCouplingMaterialLhs, mwlsApprox));
6464  feCouplingMaterialLhs->getOpPtrVector().push_back(
6465  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
6466  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6467  feCouplingMaterialLhs, mwlsApprox, mwlsRhoTagName, true, true));
6468  }
6469 
6470  feCouplingMaterialLhs->getOpPtrVector().push_back(
6471  new HookeElement::OpCalculateStiffnessScaledByDensityField(
6472  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6473  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
6474  rHo0));
6475  space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
6476  feCouplingMaterialLhs->getOpPtrVector().push_back(
6477  new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
6478  "MESH_NODE_POSITIONS",
6479  data_hooke_element_at_pts));
6480  feCouplingMaterialLhs->getOpPtrVector().push_back(
6481  new HookeElement::OpCalculateStress<1>("MESH_NODE_POSITIONS",
6482  "MESH_NODE_POSITIONS",
6483  data_hooke_element_at_pts));
6484  feCouplingMaterialLhs->getOpPtrVector().push_back(
6485  new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
6486  "MESH_NODE_POSITIONS",
6487  data_hooke_element_at_pts));
6488  feCouplingMaterialLhs->getOpPtrVector().push_back(
6489  new HookeElement::OpCalculateEshelbyStress(
6490  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6491  data_hooke_element_at_pts));
6492  feCouplingMaterialLhs->getOpPtrVector().push_back(
6493  new OpTransfromSingularBaseFunctions(
6494  feCouplingMaterialLhs->singularElement,
6495  feCouplingMaterialLhs->detS, feCouplingMaterialLhs->invSJac));
6496  feCouplingMaterialLhs->getOpPtrVector().push_back(
6497  new HookeElement::OpAleLhs_dX_dX<1>("MESH_NODE_POSITIONS",
6498  "MESH_NODE_POSITIONS",
6499  data_hooke_element_at_pts));
6500  feCouplingMaterialLhs->getOpPtrVector().push_back(
6501  new HookeElement::OpAleLhsPre_dX_dx<1>("MESH_NODE_POSITIONS",
6502  "SPATIAL_POSITION",
6503  data_hooke_element_at_pts));
6504  feCouplingMaterialLhs->getOpPtrVector().push_back(
6505  new HookeElement::OpAleLhs_dX_dx("MESH_NODE_POSITIONS",
6506  "SPATIAL_POSITION",
6507  data_hooke_element_at_pts));
6508 
6509  feCouplingMaterialLhs->getOpPtrVector().push_back(
6510  new HookeElement::OpAleLhsWithDensity_dX_dX(
6511  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6512  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6513  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
6514 
6515  boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
6516  if (addSingularity) {
6517 
6518  mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
6519  mwlsApprox, &mwlsApprox->singularInitialDisplacement);
6520 
6521  feCouplingMaterialLhs->getOpPtrVector().push_back(
6522  new OpAleLhsWithDensitySingularElement_dX_dX(
6523  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6524  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6525  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0,
6526  mat_singular_disp_ptr));
6527  }
6528 
6529  feCouplingMaterialLhs->getOpPtrVector().push_back(
6530  new OpLhsBoneExplicitDerivariveWithHooke_dX(
6531  *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6532  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
6533  mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
6534  rHo0, nBone, &crackFrontNodes));
6535 
6536  feCouplingMaterialLhs->getOpPtrVector().push_back(
6537  new OpLhsBoneExplicitDerivariveWithHooke_dx(
6538  *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6539  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
6540  mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
6541  rHo0, nBone, &crackFrontNodes));
6542 
6543  } else {
6544 
6545  feCouplingMaterialLhs->getOpPtrVector().push_back(
6546  new HookeElement::OpCalculateHomogeneousStiffness<0>(
6547  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6548  data_hooke_element_at_pts));
6549  // Calculate spatial positions and gradients (of deformation) at
6550  // integration pts. This operator takes into account singularity at crack
6551  // front
6552  feCouplingMaterialLhs->getOpPtrVector().push_back(
6553  new OpGetCrackFrontDataGradientAtGaussPts(
6554  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6555  feCouplingMaterialLhs->singularElement,
6556  feCouplingMaterialLhs->invSJac));
6557  space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
6558  feCouplingMaterialLhs->getOpPtrVector().push_back(
6559  new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
6560  "MESH_NODE_POSITIONS",
6561  data_hooke_element_at_pts));
6562  feCouplingMaterialLhs->getOpPtrVector().push_back(
6563  new HookeElement::OpCalculateStress<0>("MESH_NODE_POSITIONS",
6564  "MESH_NODE_POSITIONS",
6565  data_hooke_element_at_pts));
6566  feCouplingMaterialLhs->getOpPtrVector().push_back(
6567  new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
6568  "MESH_NODE_POSITIONS",
6569  data_hooke_element_at_pts));
6570  feCouplingMaterialLhs->getOpPtrVector().push_back(
6571  new HookeElement::OpCalculateEshelbyStress(
6572  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6573  data_hooke_element_at_pts));
6574  feCouplingMaterialLhs->getOpPtrVector().push_back(
6575  new OpTransfromSingularBaseFunctions(
6576  feCouplingMaterialLhs->singularElement,
6577  feCouplingMaterialLhs->detS, feCouplingMaterialLhs->invSJac));
6578  feCouplingMaterialLhs->getOpPtrVector().push_back(
6579  new HookeElement::OpAleLhs_dX_dX<0>("MESH_NODE_POSITIONS",
6580  "MESH_NODE_POSITIONS",
6581  data_hooke_element_at_pts));
6582  feCouplingMaterialLhs->getOpPtrVector().push_back(
6583  new HookeElement::OpAleLhsPre_dX_dx<0>("MESH_NODE_POSITIONS",
6584  "SPATIAL_POSITION",
6585  data_hooke_element_at_pts));
6586  feCouplingMaterialLhs->getOpPtrVector().push_back(
6587  new HookeElement::OpAleLhs_dX_dx("MESH_NODE_POSITIONS",
6588  "SPATIAL_POSITION",
6589  data_hooke_element_at_pts));
6590  }
6591  }
6592 
6593  if (residualStressBlock != -1) {
6594  if (!mwlsApprox) {
6595  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
6596  "mwlsApprox not allocated");
6597  }
6598 
6599  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6600  feCouplingMaterialLhs->getOpPtrVector().push_back(
6601  new MWLSApprox::OpMWLSStressAtGaussPts(
6602  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6603  feCouplingMaterialLhs, mwlsApprox, mwlsStressTagName, true));
6604  else {
6605  feCouplingMaterialLhs->getOpPtrVector().push_back(
6606  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
6607  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6608  feCouplingMaterialLhs, mwlsApprox));
6609  feCouplingMaterialLhs->getOpPtrVector().push_back(
6610  new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
6611  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6612  feCouplingMaterialLhs, mwlsApprox, mwlsStressTagName, true));
6613  }
6614 
6615  feCouplingMaterialLhs->getOpPtrVector().push_back(
6616  new MWLSApprox::OpMWLSSpatialStressLhs_DX(mat_grad_pos_at_pts_ptr,
6617  mwlsApprox));
6618  feCouplingMaterialLhs->getOpPtrVector().push_back(
6619  new MWLSApprox::OpMWLSMaterialStressLhs_DX(
6620  space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
6621  &crackFrontNodes));
6622  feCouplingMaterialLhs->getOpPtrVector().push_back(
6623  new MWLSApprox::OpMWLSMaterialStressLhs_Dx(
6624  space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
6625  &crackFrontNodes));
6626  }
6627 
6628  MoFEMFunctionReturnHot(0);
6629 }

◆ assembleElasticDM() [1/2]

MoFEMErrorCode FractureMechanics::CrackPropagation::assembleElasticDM ( const int  verb = QUIET,
const bool  debug = false 
)

create elastic finite element instance for spatial assembly

Parameters
verbcompilation parameter determining the amount of information printed (not in use here)
debugflag for debugging
Returns
error code

Definition at line 4601 of file CrackPropagation.cpp.

4602  {
4603  MoFEMFunctionBegin;
4604 
4605  CHKERR assembleElasticDM(mwlsStressTagName, verb, debug);
4606 
4607  MoFEMFunctionReturn(0);
4608 }

◆ assembleElasticDM() [2/2]

MoFEMErrorCode FractureMechanics::CrackPropagation::assembleElasticDM ( const std::string  mwls_stress_tag_name,
const int  verb = QUIET,
const bool  debug = false 
)

create elastic finite element instance for spatial assembly

Parameters
mwls_stress_tag_nameName of internal stress tag
close_crackif true, crack surface is closed
verbcompilation parameter determining the amount of information printed (not in use here)
debugflag for debugging
Returns
error code

Rhs

Lhs

Definition at line 3795 of file CrackPropagation.cpp.

3796  {
3797  MoFEMFunctionBegin;
3798 
3799  // Create elastic element data structure
3800  elasticFe = boost::make_shared<NonlinearElasticElement>(mField, ELASTIC_TAG);
3801 
3802  // Create elastic element finite element instance for residuals. Note that
3803  // this element approx. singularity at crack front.
3804  feRhs = boost::make_shared<CrackFrontElement>(
3805  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
3806  crackFrontElements, addSingularity);
3807 
3808  // Create finite element instance for assembly of tangent matrix
3809  feLhs = boost::make_shared<CrackFrontElement>(
3810  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
3811  crackFrontElements, addSingularity);
3812 
3813  // Arbitrary lagrangian formulation, mesh node positions are taken into
3814  // account by operators.
3815  feRhs->meshPositionsFieldName = "NONE";
3816  feLhs->meshPositionsFieldName = "NONE";
3817  feRhs->addToRule = 0;
3818  feLhs->addToRule = 0;
3819 
3820  // Create operators to calculate finite element matrices for elastic element
3821  onlyHooke = true;
3822  using BlockData = NonlinearElasticElement::BlockData;
3823  boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
3824  elasticFe, &elasticFe->setOfBlocks);
3825  {
3826  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
3827  mField, BLOCKSET | MAT_ELASTICSET, it)) {
3828 
3829  // Get data from block
3830  Mat_Elastic mydata;
3831  CHKERR it->getAttributeDataStructure(mydata);
3832  int id = it->getMeshsetId();
3833  EntityHandle meshset = it->getMeshset();
3834  CHKERR mField.get_moab().get_entities_by_type(
3835  meshset, MBTET, elasticFe->setOfBlocks[id].tEts, true);
3836  elasticFe->setOfBlocks[id].iD = id;
3837  elasticFe->setOfBlocks[id].E = mydata.data.Young;
3838  elasticFe->setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
3839 
3840  // Print material properties of blocks after being read
3841  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\nMaterial block %d \n", id);
3842  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung's modulus %6.4e \n",
3843  mydata.data.Young);
3844  CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPoisson's ratio %6.4e \n\n",
3845  mydata.data.Poisson);
3846 
3847  // Set default material to elastic blocks. Create material instances
3848  switch (defaultMaterial) {
3849  case HOOKE:
3850  elasticFe->setOfBlocks[id].materialDoublePtr =
3851  boost::make_shared<Hooke<double>>();
3852  elasticFe->setOfBlocks[id].materialAdoublePtr =
3853  boost::make_shared<Hooke<adouble>>();
3854  break;
3855  case KIRCHHOFF:
3856  elasticFe->setOfBlocks[id].materialDoublePtr = boost::make_shared<
3857  NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<
3858  double>>();
3859  elasticFe->setOfBlocks[id].materialAdoublePtr = boost::make_shared<
3860  NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<
3861  adouble>>();
3862  onlyHooke = false;
3863  break;
3864  case NEOHOOKEAN:
3865  elasticFe->setOfBlocks[id].materialDoublePtr =
3866  boost::make_shared<NeoHookean<double>>();
3867  elasticFe->setOfBlocks[id].materialAdoublePtr =
3868  boost::make_shared<NeoHookean<adouble>>();
3869  onlyHooke = false;
3870  break;
3871  case BONEHOOKE:
3872  elasticFe->setOfBlocks[id].materialDoublePtr =
3873  boost::make_shared<Hooke<double>>();
3874  elasticFe->setOfBlocks[id].materialAdoublePtr =
3875  boost::make_shared<Hooke<adouble>>();
3876  break;
3877  default:
3878  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
3879  "Material type not yet implemented");
3880  break;
3881  }
3882  }
3883  if (onlyHookeFromOptions == PETSC_FALSE)
3884  onlyHooke = false;
3885 
3886  elasticFe->commonData.spatialPositions = "SPATIAL_POSITION";
3887  elasticFe->commonData.meshPositions = "MESH_NODE_POSITIONS";
3888 
3889  auto data_hooke_element_at_pts =
3890  boost::make_shared<HookeElement::DataAtIntegrationPts>();
3891 
3892  // calculate material positions at integration points
3893  auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3894  feRhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3895  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
3896  boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
3897  if (residualStressBlock != -1 || densityMapBlock != -1) {
3898  mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3899  feRhs->getOpPtrVector().push_back(
3900  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
3901  mat_grad_pos_at_pts_ptr));
3902  }
3903 
3904  /// Rhs
3905  if (!onlyHooke) {
3906  // Calculate spatial positions and gradients (of deformation) at
3907  // integration pts. This operator takes into account singularity at crack
3908  // front
3909  feRhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
3910  "SPATIAL_POSITION", elasticFe->commonData, feRhs->singularElement,
3911  feRhs->invSJac));
3912  // Calculate material positions and gradients at integration pts.
3913  feRhs->getOpPtrVector().push_back(
3914  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
3915  "MESH_NODE_POSITIONS", elasticFe->commonData));
3916  // Transfrom base functions to create singularity at crack front
3917  feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
3918  feRhs->singularElement, feRhs->detS, feRhs->invSJac));
3919  // Iterate over material blocks
3920  map<int, NonlinearElasticElement::BlockData>::iterator sit =
3921  elasticFe->setOfBlocks.begin();
3922  for (; sit != elasticFe->setOfBlocks.end(); sit++) {
3923  // Evaluate stress at integration pts.
3924  feRhs->getOpPtrVector().push_back(
3925  new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
3926  "SPATIAL_POSITION", sit->second, elasticFe->commonData,
3927  ELASTIC_TAG, false, true, false));
3928  // Assemble internal forces for right hand side
3929  feRhs->getOpPtrVector().push_back(
3930  new NonlinearElasticElement::OpRhsPiolaKirchhoff(
3931  "SPATIAL_POSITION", sit->second, elasticFe->commonData));
3932  }
3933  } else {
3934  feRhs->getOpPtrVector().push_back(
3935  new OpCalculateVectorFieldGradient<3, 3>(
3936  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
3937  mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
3938 
3939  // Evaluate density at integration points fro material like bone
3940  if (defaultMaterial == BONEHOOKE) {
3941  if (!mwlsApprox) {
3942  // Create instance for moving least square approximation
3943  mwlsApprox = boost::make_shared<MWLSApprox>(mField);
3944  // Load mesh with stresses
3945  CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
3946  MeshsetsManager *block_manager_ptr;
3947  CHKERR mField.getInterface(block_manager_ptr);
3948  CHKERR block_manager_ptr->getEntitiesByDimension(
3949  densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3950  Tag th_rho;
3951  if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
3952  th_rho) != MB_SUCCESS) {
3953  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3954  "Density tag name %s cannot be found. Please "
3955  "provide the correct name.",
3956  mwlsRhoTagName.c_str());
3957  }
3958  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
3959  th_rho);
3960  CHKERR mwlsApprox->getValuesToNodes(th_rho);
3961  } else {
3962  mwlsApprox->tetsInBlock.clear();
3963  CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
3964  densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3965  }
3966  // Calculate spatial positions and gradients (of deformation) at
3967  // integration pts. This operator takes into account singularity at
3968  // crack front
3969  feRhs->getOpPtrVector().push_back(
3970  new OpGetCrackFrontDataGradientAtGaussPts(
3971  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
3972  feRhs->singularElement, feRhs->invSJac));
3973  // Evaluate density at integration pts.
3974  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
3975  feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
3976  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
3977  mwlsRhoTagName, true, true));
3978  else {
3979  feRhs->getOpPtrVector().push_back(
3980  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
3981  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
3982  mwlsApprox));
3983  feRhs->getOpPtrVector().push_back(
3984  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
3985  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
3986  mwlsApprox, mwlsRhoTagName, true, false));
3987  }
3988 
3989  feRhs->getOpPtrVector().push_back(
3990  new HookeElement::OpCalculateStiffnessScaledByDensityField(
3991  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
3992  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
3993  rHo0));
3994  feRhs->getOpPtrVector().push_back(
3995  new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
3996  "SPATIAL_POSITION",
3997  data_hooke_element_at_pts));
3998  feRhs->getOpPtrVector().push_back(
3999  new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
4000  "SPATIAL_POSITION",
4001  data_hooke_element_at_pts));
4002  } else {
4003  feRhs->getOpPtrVector().push_back(
4004  new HookeElement::OpCalculateHomogeneousStiffness<0>(
4005  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
4006  data_hooke_element_at_pts));
4007  // Calculate spatial positions and gradients (of deformation) at
4008  // integration pts. This operator takes into account singularity at
4009  // crack front
4010  feRhs->getOpPtrVector().push_back(
4011  new OpGetCrackFrontDataGradientAtGaussPts(
4012  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4013  feRhs->singularElement, feRhs->invSJac));
4014  feRhs->getOpPtrVector().push_back(
4015  new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
4016  "SPATIAL_POSITION",
4017  data_hooke_element_at_pts));
4018  feRhs->getOpPtrVector().push_back(
4019  new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
4020  "SPATIAL_POSITION",
4021  data_hooke_element_at_pts));
4022  }
4023 
4024  feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4025  feRhs->singularElement, feRhs->detS, feRhs->invSJac));
4026  feRhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dx(
4027  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
4028  }
4029 
4030  // Add operators if internal stresses are present
4031  if (residualStressBlock != -1) {
4032  if (!mwlsApprox) {
4033  // Create instance for moving least square approximation
4034  mwlsApprox = boost::make_shared<MWLSApprox>(mField);
4035  // Load mesh with stresses
4036  CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
4037  MeshsetsManager *block_manager_ptr;
4038  CHKERR mField.getInterface(block_manager_ptr);
4039  CHKERR block_manager_ptr->getEntitiesByDimension(
4040  residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
4041  Tag th;
4042  if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
4043  th) != MB_SUCCESS) {
4044  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4045  "Internal stress tag name %s cannot be found. Please "
4046  "provide the correct name.",
4047  mwls_stress_tag_name.substr(4).c_str());
4048  }
4049  }
4050 
4051  Tag th;
4052  if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
4053  th) != MB_SUCCESS) {
4054  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4055  "Internal stress tag name %s cannot be found. Please "
4056  "provide the correct name.",
4057  mwls_stress_tag_name.substr(4).c_str());
4058  }
4059  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
4060  th);
4061  CHKERR mwlsApprox->getValuesToNodes(th);
4062 
4063  // Evaluate internal stresses at integration pts.
4064  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4065  feRhs->getOpPtrVector().push_back(
4066  new MWLSApprox::OpMWLSStressAtGaussPts(
4067  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
4068  mwls_stress_tag_name, false));
4069  else {
4070  feRhs->getOpPtrVector().push_back(
4071  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
4072  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
4073  mwlsApprox));
4074  feRhs->getOpPtrVector().push_back(
4075  new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
4076  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
4077  mwls_stress_tag_name, false));
4078  }
4079 
4080  // Assemble internall stresses forces
4081  feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressRhs(
4082  mat_grad_pos_at_pts_ptr, mwlsApprox));
4083  if (addAnalyticalInternalStressOperators) {
4084  feRhs->getOpPtrVector().push_back(
4085  new HookeElement::OpAnalyticalInternalAleStrain_dx<0>(
4086  "SPATIAL_POSITION", data_hooke_element_at_pts,
4087  analyticalStrainFunction,
4088  boost::shared_ptr<MatrixDouble>(
4089  mwlsApprox, &mwlsApprox->mwlsMaterialPositions)));
4090  }
4091  }
4092 
4093  /// Lhs
4094  if (!onlyHooke) {
4095  // Calculate spatial positions and gradients (of deformation) at
4096  // integration pts. This operator takes into account singularity at crack
4097  // front
4098  feLhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
4099  "SPATIAL_POSITION", elasticFe->commonData, feLhs->singularElement,
4100  feLhs->invSJac));
4101  // Calculate material positions and gradients at integration pts.
4102  feLhs->getOpPtrVector().push_back(
4103  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
4104  "MESH_NODE_POSITIONS", elasticFe->commonData));
4105  // Transform base functions to get singular base functions at crack front
4106  feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4107  feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4108  // Iterate over material blocks
4109  map<int, NonlinearElasticElement::BlockData>::iterator sit =
4110  elasticFe->setOfBlocks.begin();
4111  for (; sit != elasticFe->setOfBlocks.end(); sit++) {
4112  // Evaluate stress at integration pts
4113  feLhs->getOpPtrVector().push_back(
4114  new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
4115  "SPATIAL_POSITION", sit->second, elasticFe->commonData,
4116  ELASTIC_TAG, true, true, false));
4117  // Assemble tanget matrix, derivative of spatial poisotions
4118  feLhs->getOpPtrVector().push_back(
4119  new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(
4120  "SPATIAL_POSITION", "SPATIAL_POSITION", sit->second,
4121  elasticFe->commonData));
4122  }
4123  } else {
4124  // calculate material positions at integration points
4125  feLhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4126  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4127  boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
4128  if (residualStressBlock != -1 || densityMapBlock != -1) {
4129  mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4130  feLhs->getOpPtrVector().push_back(
4131  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
4132  mat_grad_pos_at_pts_ptr));
4133  }
4134 
4135  feLhs->getOpPtrVector().push_back(
4136  new OpCalculateVectorFieldGradient<3, 3>(
4137  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
4138  mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
4139 
4140  if (defaultMaterial == BONEHOOKE) {
4141  if (!mwlsApprox) {
4142  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4143  "mwlsApprox not allocated");
4144  } else {
4145  mwlsApprox->tetsInBlock.clear();
4146  CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
4147  densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock,
4148  true); // FIXME: do we still need this?
4149  }
4150 
4151  // Calculate spatial positions and gradients (of deformation) at
4152  // integration pts. This operator takes into account singularity at
4153  // crack front
4154  feLhs->getOpPtrVector().push_back(
4155  new OpGetCrackFrontDataGradientAtGaussPts(
4156  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4157  feLhs->singularElement, feLhs->invSJac));
4158 
4159  // Evaluate density at integration pts.
4160  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4161  feLhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
4162  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs, mwlsApprox,
4163  mwlsRhoTagName, true, true));
4164  else {
4165  feLhs->getOpPtrVector().push_back(
4166  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
4167  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
4168  mwlsApprox));
4169  feLhs->getOpPtrVector().push_back(
4170  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
4171  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
4172  mwlsApprox, mwlsRhoTagName, true, false));
4173  }
4174 
4175  feLhs->getOpPtrVector().push_back(
4176  new HookeElement::OpCalculateStiffnessScaledByDensityField(
4177  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
4178  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
4179  rHo0));
4180  feLhs->getOpPtrVector().push_back(
4181  new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
4182  "SPATIAL_POSITION",
4183  data_hooke_element_at_pts));
4184  feLhs->getOpPtrVector().push_back(
4185  new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
4186  "SPATIAL_POSITION",
4187  data_hooke_element_at_pts));
4188 
4189  feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4190  feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4191  feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<1>(
4192  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
4193 
4194  } else {
4195  feLhs->getOpPtrVector().push_back(
4196  new HookeElement::OpCalculateHomogeneousStiffness<0>(
4197  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
4198  data_hooke_element_at_pts));
4199  // Calculate spatial positions and gradients (of deformation) at
4200  // integration pts. This operator takes into account singularity at
4201  // crack front
4202  feLhs->getOpPtrVector().push_back(
4203  new OpGetCrackFrontDataGradientAtGaussPts(
4204  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4205  feLhs->singularElement, feLhs->invSJac));
4206  feLhs->getOpPtrVector().push_back(
4207  new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
4208  "SPATIAL_POSITION",
4209  data_hooke_element_at_pts));
4210  feLhs->getOpPtrVector().push_back(
4211  new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
4212  "SPATIAL_POSITION",
4213  data_hooke_element_at_pts));
4214  // Transform base functions to get singular base functions at crack
4215  // front
4216  feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4217  feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4218  feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<0>(
4219  "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
4220  }
4221  }
4222  }
4223 
4224  // Assembly forces
4225 
4226  surfaceForces =
4227  boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4228  // Surface forces
4229  {
4230  string fe_name_str = "FORCE_FE";
4231  surfaceForces->insert(fe_name_str, new NeumannForcesSurface(mField));
4232  auto &nf = surfaceForces->at(fe_name_str);
4233  auto &fe = nf.getLoopFe();
4234 
4235  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe, false, false);
4236 
4237  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
4238  it)) {
4239  CHKERR surfaceForces->at(fe_name_str)
4240  .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), true);
4241  }
4242 
4243  const string block_set_force_name("FORCE");
4244  // search for block named FORCE and add its attributes to FORCE_FE element
4245  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
4246  if (it->getName().compare(0, block_set_force_name.length(),
4247  block_set_force_name) == 0) {
4248  CHKERR surfaceForces->at(fe_name_str)
4249  .addForce("SPATIAL_POSITION", PETSC_NULL, (it->getMeshsetId()),
4250  true, true);
4251  }
4252  }
4253  }
4254 
4255  if (isSurfaceForceAle) {
4256  surfaceForceAle =
4257  boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4258  commonDataSurfaceForceAle =
4259  boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
4260  commonDataSurfaceForceAle->forcesOnlyOnEntitiesRow = crackFrontNodes;
4261 
4262  string fe_name_str = "FORCE_FE_ALE";
4263  surfaceForceAle->insert(fe_name_str, new NeumannForcesSurface(mField));
4264 
4265  auto &nf = surfaceForceAle->at(fe_name_str);
4266  auto &fe_lhs = nf.getLoopFeLhs();
4267  auto &fe_mat_rhs = nf.getLoopFeMatRhs();
4268  auto &fe_mat_lhs = nf.getLoopFeMatLhs();
4269 
4270  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_lhs, false, false);
4271  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_rhs, false, false);
4272  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_lhs, false, false);
4273 
4274  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
4275  bit)) {
4276  CHKERR surfaceForceAle->at(fe_name_str)
4277  .addForceAle("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
4278  commonDataSurfaceForceAle, "MATERIAL", PETSC_NULL,
4279  PETSC_NULL, bit->getMeshsetId(), true, false,
4280  ignoreMaterialForce);
4281  }
4282 
4283  const string block_set_force_name("FORCE");
4284  // search for block named FORCE and add its attributes to FORCE_FE element
4285  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
4286  if (it->getName().compare(0, block_set_force_name.length(),
4287  block_set_force_name) == 0) {
4288  CHKERR surfaceForceAle->at(fe_name_str)
4289  .addForceAle("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
4290  commonDataSurfaceForceAle, "MATERIAL", PETSC_NULL,
4291  PETSC_NULL, it->getMeshsetId(), true, true,
4292  ignoreMaterialForce);
4293  }
4294  }
4295  }
4296 
4297  // assemble surface pressure
4298  surfacePressure =
4299  boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4300  {
4301  string fe_name_str = "PRESSURE_FE";
4302  surfacePressure->insert(fe_name_str, new NeumannForcesSurface(mField));
4303 
4304  auto &nf = surfacePressure->at(fe_name_str);
4305  auto &fe = nf.getLoopFe();
4306 
4307  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe, false, false);
4308 
4309  // iterate over sidestep where pressure is applied
4310  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
4311  mField, SIDESET | PRESSURESET, it)) {
4312  CHKERR surfacePressure->at(fe_name_str)
4313  .addPressure("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
4314  true);
4315  }
4316 
4317  const string block_set_pressure_name("PRESSURE");
4318  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
4319  if (it->getName().compare(0, block_set_pressure_name.length(),
4320  block_set_pressure_name) == 0) {
4321  CHKERR surfacePressure->at(fe_name_str)
4322  .addPressure("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
4323  true, true);
4324  }
4325  }
4326 
4327  const string block_set_linear_pressure_name("LINEAR_PRESSURE");
4328  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
4329  if (it->getName().compare(0, block_set_linear_pressure_name.length(),
4330  block_set_linear_pressure_name) == 0) {
4331  CHKERR surfacePressure->at(fe_name_str)
4332  .addLinearPressure("SPATIAL_POSITION", PETSC_NULL,
4333  it->getMeshsetId(), true);
4334  }
4335  }
4336  }
4337  // assemble surface pressure (ALE)
4338  if (isPressureAle) {
4339  surfacePressureAle =
4340  boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4341  commonDataSurfacePressureAle =
4342  boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
4343  commonDataSurfacePressureAle->forcesOnlyOnEntitiesRow = crackFrontNodes;
4344  string fe_name_str = "PRESSURE_ALE";
4345  surfacePressureAle->insert(fe_name_str, new NeumannForcesSurface(mField));
4346 
4347  auto &nf = surfacePressureAle->at(fe_name_str);
4348  auto &fe_lhs = nf.getLoopFeLhs();
4349  auto &fe_mat_rhs = nf.getLoopFeMatRhs();
4350  auto &fe_mat_lhs = nf.getLoopFeMatLhs();
4351 
4352  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_lhs, false, false);
4353  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_rhs, false, false);
4354  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_lhs, false, false);
4355 
4356  // iterate over sidesets where pressure is applied
4357  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
4358  mField, SIDESET | PRESSURESET, it)) {
4359  CHKERR surfacePressureAle->at(fe_name_str)
4360  .addPressureAle("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
4361  commonDataSurfacePressureAle, "MATERIAL", PETSC_NULL,
4362  PETSC_NULL, it->getMeshsetId(), true);
4363  }
4364  }
4365 
4366  if (areSpringsAle) {
4367  feRhsSpringALEMaterial =
4368  boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
4369 
4370  feLhsSpringALE =
4371  boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
4372  feLhsSpringALEMaterial =
4373  boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
4374  commonDataSpringsALE =
4375  boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(mField);
4376  commonDataSpringsALE->forcesOnlyOnEntitiesRow = crackFrontNodes;
4377 
4378  CHKERR MetaSpringBC::setSpringOperatorsMaterial(
4379  mField, feLhsSpringALE, feLhsSpringALEMaterial, feRhsSpringALEMaterial,
4380  commonDataSpringsALE, "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
4381  "MATERIAL");
4382  }
4383 
4384  // Implementation of spring element
4385  // Create new instances of face elements for springs
4386  feSpringLhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
4387  new FaceElementForcesAndSourcesCore(mField));
4388  feSpringRhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
4389  new FaceElementForcesAndSourcesCore(mField));
4390 
4391  CHKERR MetaSpringBC::setSpringOperators(mField, feSpringLhsPtr,
4392  feSpringRhsPtr, "SPATIAL_POSITION",
4393  "MESH_NODE_POSITIONS");
4394 
4395  bool use_eigen_pos =
4396  solveEigenStressProblem && (mwls_stress_tag_name == mwlsStressTagName);
4397 
4398  bool use_eigen_pos_simple_contact =
4399  use_eigen_pos && useEigenPositionsSimpleContact;
4400 
4401  // Set contact operators
4402  feRhsSimpleContact =
4403  boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4404  feLhsSimpleContact =
4405  boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4406  commonDataSimpleContact =
4407  boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(mField);
4408 
4409  if (!ignoreContact) {
4410  if (printContactState) {
4411  feRhsSimpleContact->contactStateVec =
4412  commonDataSimpleContact->gaussPtsStateVec;
4413  }
4414  contactProblem->setContactOperatorsRhs(
4415  feRhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4416  "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4417  "EIGEN_SPATIAL_POSITIONS");
4418  contactProblem->setMasterForceOperatorsRhs(
4419  feRhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4420  "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4421  "EIGEN_SPATIAL_POSITIONS");
4422 
4423  contactProblem->setContactOperatorsLhs(
4424  feLhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4425  "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4426  "EIGEN_SPATIAL_POSITIONS");
4427 
4428  contactProblem->setMasterForceOperatorsLhs(
4429  feLhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4430  "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4431  "EIGEN_SPATIAL_POSITIONS");
4432  }
4433 
4434  // Close crack constrains
4435  if (solveEigenStressProblem) {
4436  closeCrackConstrains =
4437  boost::shared_ptr<BothSurfaceConstrains>(new BothSurfaceConstrains(
4438  mField, "LAMBDA_CLOSE_CRACK", "SPATIAL_POSITION"));
4439  closeCrackConstrains->aLpha = 1;
4440  }
4441 
4442  // Set contact operators
4443 
4444  if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
4445  // Set contact operators
4446  feRhsSimpleContactALEMaterial =
4447  boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4448  feLhsSimpleContactALEMaterial =
4449  boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4450  feLhsSimpleContactALE =
4451  boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4452  commonDataSimpleContactALE =
4453  boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
4454  mField);
4455 
4456  commonDataSimpleContactALE->forcesOnlyOnEntitiesRow = crackFrontNodes;
4457 
4458  contactProblem->setContactOperatorsRhsALEMaterial(
4459  feRhsSimpleContactALEMaterial, commonDataSimpleContactALE,
4460  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAMBDA_CONTACT",
4461  "MAT_CONTACT");
4462 
4463  contactProblem->setContactOperatorsLhsALEMaterial(
4464  feLhsSimpleContactALEMaterial, commonDataSimpleContactALE,
4465  "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAMBDA_CONTACT",
4466  "MAT_CONTACT");
4467 
4468  contactProblem->setContactOperatorsLhsALE(
4469  feLhsSimpleContactALE, commonDataSimpleContactALE, "SPATIAL_POSITION",
4470  "MESH_NODE_POSITIONS", "LAMBDA_CONTACT", use_eigen_pos_simple_contact,
4471  "EIGEN_SPATIAL_POSITIONS");
4472  }
4473 
4474  feRhsMortarContact =
4475  boost::make_shared<MortarContactProblem::MortarContactElement>(
4476  mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
4477  "MESH_NODE_POSITIONS");
4478  feLhsMortarContact =
4479  boost::make_shared<MortarContactProblem::MortarContactElement>(
4480  mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
4481  "MESH_NODE_POSITIONS");
4482  commonDataMortarContact =
4483  boost::make_shared<MortarContactProblem::CommonDataMortarContact>(mField);
4484 
4485  if (!ignoreContact) {
4486  if (printContactState) {
4487  feRhsMortarContact->contactStateVec =
4488  commonDataMortarContact->gaussPtsStateVec;
4489  }
4490  mortarContactProblemPtr->setContactOperatorsRhs(
4491  feRhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4492  "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4493 
4494  mortarContactProblemPtr->setMasterForceOperatorsRhs(
4495  feRhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4496  "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4497 
4498  mortarContactProblemPtr->setContactOperatorsLhs(
4499  feLhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4500  "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4501 
4502  mortarContactProblemPtr->setMasterForceOperatorsLhs(
4503  feLhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4504  "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4505  }
4506 
4507  // set contact post proc operators
4508 
4509  fePostProcSimpleContact =
4510  boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4511 
4512  fePostProcMortarContact =
4513  boost::make_shared<MortarContactProblem::MortarContactElement>(
4514  mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
4515  "MESH_NODE_POSITIONS");
4516 
4517  if (!ignoreContact) {
4518  CHKERR contactProblem->setContactOperatorsForPostProc(
4519  fePostProcSimpleContact, commonDataSimpleContact, mField,
4520  "SPATIAL_POSITION", "LAMBDA_CONTACT", contactPostProcMoab, false,
4521  use_eigen_pos_simple_contact, "EIGEN_SPATIAL_POSITIONS");
4522 
4523  CHKERR mortarContactProblemPtr->setContactOperatorsForPostProc(
4524  fePostProcMortarContact, commonDataMortarContact, mField,
4525  "SPATIAL_POSITION", "LAMBDA_CONTACT", contactPostProcMoab, false,
4526  use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4527  }
4528 
4529  // assemble nodal forces
4530  nodalForces = boost::make_shared<boost::ptr_map<string, NodalForce>>();
4531  {
4532  string fe_name_str = "FORCE_FE";
4533  nodalForces->insert(fe_name_str, new NodalForce(mField));
4534  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
4535  it)) {
4536  CHKERR nodalForces->at(fe_name_str)
4537  .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), false);
4538  }
4539  }
4540  // assemble edge forces
4541  edgeForces = boost::make_shared<boost::ptr_map<string, EdgeForce>>();
4542  {
4543  string fe_name_str = "FORCE_FE";
4544  edgeForces->insert(fe_name_str, new EdgeForce(mField));
4545  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
4546  it)) {
4547  CHKERR edgeForces->at(fe_name_str)
4548  .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), false);
4549  }
4550  }
4551 
4552  // Kinematic boundary conditions
4553  spatialDirichletBc = boost::shared_ptr<DirichletSpatialPositionsBc>(
4554  new DirichletSpatialPositionsBc(mField, "SPATIAL_POSITION", PETSC_NULL,
4555  PETSC_NULL, PETSC_NULL));
4556 
4557 // Add boundary conditions for displacements given by analytical function
4558 #ifdef __ANALITICAL_DISPLACEMENT__
4559  analyticalDirichletBc = boost::shared_ptr<AnalyticalDirichletBC::DirichletBC>(
4560  new AnalyticalDirichletBC::DirichletBC(
4561  mField, "SPATIAL_POSITION", PETSC_NULL, PETSC_NULL, PETSC_NULL));
4562 #endif //__ANALITICAL_DISPLACEMENT__
4563 
4564 // Analytical forces
4565 #ifdef __ANALITICAL_TRACTION__
4566  MeshsetsManager *meshset_manager_ptr;
4567  CHKERR mField.getInterface(meshset_manager_ptr);
4568  if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
4569  if (!analiticalSurfaceElement) {
4570  analiticalSurfaceElement =
4571  boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4572  {
4573  string fe_name_str = "ANALITICAL_TRACTION";
4574  analiticalSurfaceElement->insert(fe_name_str,
4575  new NeumannForcesSurface(mField));
4576  analiticalSurfaceElement->at(fe_name_str).fe.addToRule = 1;
4577  const CubitMeshSets *cubit_meshset_ptr;
4578  meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
4579  &cubit_meshset_ptr);
4580  Range faces;
4581  CHKERR meshset_manager_ptr->getEntitiesByDimension(
4582  cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, faces, true);
4583  analiticalSurfaceElement->at(fe_name_str)
4584  .analyticalForceOp.push_back(new AnalyticalForces());
4585  analiticalSurfaceElement->at(fe_name_str)
4586  .fe.getOpPtrVector()
4587  .push_back(new OpAnalyticalSpatialTraction(
4588  mField, setSingularCoordinates, crackFrontNodes,
4589  crackFrontNodesEdges, crackFrontElements, addSingularity,
4590  "SPATIAL_POSITION", faces,
4591  analiticalSurfaceElement->at(fe_name_str).methodsOp,
4592  analiticalSurfaceElement->at(fe_name_str).analyticalForceOp));
4593  }
4594  }
4595  }
4596 #endif //__ANALITICAL_TRACTION__
4597 
4598  MoFEMFunctionReturn(0);
4599 }

◆ assembleMaterialForcesDM()

MoFEMErrorCode FractureMechanics::CrackPropagation::assembleMaterialForcesDM ( DM  dm,
const int  verb = QUIET,
const bool  debug = false 
)

create material element instance

Operators to assemble rhs

Definition at line 5531 of file CrackPropagation.cpp.

5532  {
5533  MoFEMFunctionBegin;
5534  if (!elasticFe)
5535  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
5536  "Elastic element instance not created");
5537 
5538  // Create element
5539  materialFe = boost::shared_ptr<NonlinearElasticElement>(
5540  new NonlinearElasticElement(mField, MATERIAL_TAG));
5541 
5542  // Create material data blocks
5543  for (std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
5544  elasticFe->setOfBlocks.begin();
5545  sit != elasticFe->setOfBlocks.end(); sit++) {
5546  materialFe->setOfBlocks[sit->first] = sit->second;
5547  materialFe->setOfBlocks[sit->first].forcesOnlyOnEntitiesRow =
5548  crackFrontNodes;
5549  }
5550  materialFe->commonData.spatialPositions = "SPATIAL_POSITION";
5551  materialFe->commonData.meshPositions = "MESH_NODE_POSITIONS";
5552 
5553  // create finite element instances for the right hand side and left hand side
5554  feMaterialRhs = boost::make_shared<CrackFrontElement>(
5555  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
5556  crackFrontElements, addSingularity);
5557  feMaterialLhs = boost::make_shared<CrackFrontElement>(
5558  mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
5559  crackFrontElements, addSingularity);
5560  feMaterialRhs->meshPositionsFieldName = "NONE";
5561  feMaterialLhs->meshPositionsFieldName = "NONE";
5562  feMaterialRhs->addToRule = 0;
5563  feMaterialLhs->addToRule = 0;
5564 
5565  using BlockData = NonlinearElasticElement::BlockData;
5566  boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
5567  elasticFe, &elasticFe->setOfBlocks);
5568  boost::shared_ptr<HookeElement::DataAtIntegrationPts>
5569  data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
5570  data_hooke_element_at_pts->forcesOnlyOnEntitiesRow = crackFrontNodes;
5571 
5572  // calculate position at integration points
5573  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(new MatrixDouble());
5574  feMaterialRhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5575  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5576  feMaterialLhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5577  "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5578  boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
5579  if (!onlyHooke && (residualStressBlock != -1 || densityMapBlock != -1)) {
5580  mat_grad_pos_at_pts_ptr =
5581  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
5582  feMaterialRhs->getOpPtrVector().push_back(
5583  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
5584  mat_grad_pos_at_pts_ptr));
5585  feMaterialLhs->getOpPtrVector().push_back(
5586  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
5587  mat_grad_pos_at_pts_ptr));
5588  }
5589  boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr;
5590  if (!onlyHooke && (residualStressBlock != -1 || densityMapBlock != -1)) {
5591  space_grad_pos_at_pts_ptr =
5592  boost::shared_ptr<MatrixDouble>(new MatrixDouble());
5593  feMaterialRhs->getOpPtrVector().push_back(
5594  new OpGetCrackFrontDataGradientAtGaussPts(
5595  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
5596  feMaterialRhs->singularElement, feMaterialRhs->invSJac));
5597  feMaterialLhs->getOpPtrVector().push_back(
5598  new OpGetCrackFrontDataGradientAtGaussPts(
5599  "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
5600  feMaterialLhs->singularElement, feMaterialLhs->invSJac));
5601  }
5602 
5603  /// Operators to assemble rhs
5604  if (!onlyHooke) {
5605  feMaterialRhs->getOpPtrVector().push_back(
5606  new OpGetCrackFrontCommonDataAtGaussPts(
5607  "SPATIAL_POSITION", materialFe->commonData,
5608  feMaterialRhs->singularElement, feMaterialRhs->invSJac));
5609  feMaterialRhs->getOpPtrVector().push_back(
5610  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
5611  "MESH_NODE_POSITIONS", materialFe->commonData));
5612 
5613  feMaterialRhs->getOpPtrVector().push_back(
5614  new OpTransfromSingularBaseFunctions(feMaterialRhs->singularElement,
5615  feMaterialRhs->detS,
5616  feMaterialRhs->invSJac));
5617  for (map<int, NonlinearElasticElement::BlockData>::iterator sit =
5618  materialFe->setOfBlocks.begin();
5619  sit != materialFe->setOfBlocks.end(); sit++) {
5620  feMaterialRhs->getOpPtrVector().push_back(
5621  new NonlinearElasticElement::OpJacobianEshelbyStress(
5622  "SPATIAL_POSITION", sit->second, materialFe->commonData,
5623  MATERIAL_TAG, false, true));
5624  feMaterialRhs->getOpPtrVector().push_back(
5625  new NonlinearElasticElement::OpRhsEshelbyStress(
5626  "MESH_NODE_POSITIONS", sit->second, materialFe->commonData));
5627  }
5628 
5629  } else { // HOOKE
5630 
5631  feMaterialRhs->getOpPtrVector().push_back(
5632  new OpCalculateVectorFieldGradient<3, 3>(
5633  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5634  mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
5635 
5636  feMaterialRhs->getOpPtrVector().push_back(
5637  new OpGetCrackFrontDataGradientAtGaussPts(
5638  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5639  feMaterialRhs->singularElement, feMaterialRhs->invSJac));
5640  space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
5641 
5642  if (defaultMaterial == BONEHOOKE) {
5643  if (!mwlsApprox) {
5644  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5645  "mwlsApprox not allocated");
5646  }
5647 
5648  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5649  feMaterialRhs->getOpPtrVector().push_back(
5650  new MWLSApprox::OpMWLSRhoAtGaussPts(
5651  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5652  mwlsApprox, mwlsRhoTagName, true, true));
5653  else {
5654  feMaterialRhs->getOpPtrVector().push_back(
5655  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
5656  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5657  mwlsApprox));
5658  feMaterialRhs->getOpPtrVector().push_back(
5659  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5660  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5661  mwlsApprox, mwlsRhoTagName, true, false));
5662  }
5663 
5664  feMaterialRhs->getOpPtrVector().push_back(
5665  new HookeElement::OpCalculateStiffnessScaledByDensityField(
5666  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5667  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
5668  rHo0));
5669  feMaterialRhs->getOpPtrVector().push_back(
5670  new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5671  "MESH_NODE_POSITIONS",
5672  data_hooke_element_at_pts));
5673  feMaterialRhs->getOpPtrVector().push_back(
5674  new HookeElement::OpCalculateStress<1>("MESH_NODE_POSITIONS",
5675  "MESH_NODE_POSITIONS",
5676  data_hooke_element_at_pts));
5677 
5678  } else {
5679  feMaterialRhs->getOpPtrVector().push_back(
5680  new HookeElement::OpCalculateHomogeneousStiffness<0>(
5681  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5682  data_hooke_element_at_pts));
5683  feMaterialRhs->getOpPtrVector().push_back(
5684  new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5685  "MESH_NODE_POSITIONS",
5686  data_hooke_element_at_pts));
5687  feMaterialRhs->getOpPtrVector().push_back(
5688  new HookeElement::OpCalculateStress<0>("MESH_NODE_POSITIONS",
5689  "MESH_NODE_POSITIONS",
5690  data_hooke_element_at_pts));
5691  }
5692  feMaterialRhs->getOpPtrVector().push_back(
5693  new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
5694  "MESH_NODE_POSITIONS",
5695  data_hooke_element_at_pts));
5696 
5697  feMaterialRhs->getOpPtrVector().push_back(
5698  new HookeElement::OpCalculateEshelbyStress("MESH_NODE_POSITIONS",
5699  "MESH_NODE_POSITIONS",
5700  data_hooke_element_at_pts));
5701  feMaterialRhs->getOpPtrVector().push_back(
5702  new OpTransfromSingularBaseFunctions(feMaterialRhs->singularElement,
5703  feMaterialRhs->detS,
5704  feMaterialRhs->invSJac));
5705 
5706  feMaterialRhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dX(
5707  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5708  data_hooke_element_at_pts));
5709 
5710  if (defaultMaterial == BONEHOOKE) {
5711  feMaterialRhs->getOpPtrVector().push_back(
5712  new OpRhsBoneExplicitDerivariveWithHooke(
5713  *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5714  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->tetsInBlock, rHo0,
5715  nBone, &crackFrontNodes));
5716  }
5717  }
5718 
5719  if (residualStressBlock != -1) {
5720  if (!mwlsApprox) {
5721  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5722  "mwlsApprox not allocated");
5723  } else {
5724  mwlsApprox->tetsInBlock.clear();
5725  CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
5726  residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
5727  }
5728 
5729  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5730  feMaterialRhs->getOpPtrVector().push_back(
5731  new MWLSApprox::OpMWLSStressAtGaussPts(
5732  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5733  mwlsApprox, mwlsStressTagName, false));
5734  else {
5735  feMaterialRhs->getOpPtrVector().push_back(
5736  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
5737  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5738  mwlsApprox));
5739  feMaterialRhs->getOpPtrVector().push_back(
5740  new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
5741  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5742  mwlsApprox, mwlsStressTagName, false));
5743  }
5744 
5745  feMaterialRhs->getOpPtrVector().push_back(
5746  new MWLSApprox::OpMWLSMaterialStressRhs(space_grad_pos_at_pts_ptr,
5747  mat_grad_pos_at_pts_ptr,
5748  mwlsApprox, &crackFrontNodes));
5749  if (addAnalyticalInternalStressOperators) {
5750  feMaterialRhs->getOpPtrVector().push_back(
5751  new HookeElement::OpAnalyticalInternalAleStrain_dX<0>(
5752  "MESH_NODE_POSITIONS", data_hooke_element_at_pts,
5753  analyticalStrainFunction,
5754  boost::shared_ptr<MatrixDouble>(
5755  mwlsApprox, &mwlsApprox->mwlsMaterialPositions)));
5756  }
5757  }
5758 
5759 // Analytical forces
5760 #ifdef __ANALITICAL_TRACTION__
5761  {
5762  if (analiticalSurfaceElement) {
5763  feMaterialAnaliticalTraction =
5764  boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>(
5765  // analiticalSurfaceElement,
5766  new NeumannForcesSurface::MyTriangleFE(mField));
5767  feMaterialAnaliticalTraction->addToRule = 1;
5768  MeshsetsManager *meshset_manager_ptr;
5769  CHKERR mField.getInterface(meshset_manager_ptr);
5770  if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
5771  const CubitMeshSets *cubit_meshset_ptr;
5772  meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
5773  &cubit_meshset_ptr);
5774  Range faces;
5775  CHKERR meshset_manager_ptr->getEntitiesByDimension(
5776  cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, faces, true);
5777  feMaterialAnaliticalTraction->getOpPtrVector().push_back(
5778  new OpAnalyticalMaterialTraction(
5779  mField, setSingularCoordinates, crackFrontNodes,
5780  crackFrontNodesEdges, crackFrontElements, addSingularity,
5781  "MESH_NODE_POSITIONS", faces,
5782  analiticalSurfaceElement->at("ANALITICAL_TRACTION").methodsOp,
5783  analiticalSurfaceElement->at("ANALITICAL_TRACTION")
5784  .analyticalForceOp,
5785  &crackFrontNodes));
5786  }
5787  }
5788  }
5789 #endif //__ANALITICAL_TRACTION__
5790 
5791  // assemble tangent matrix
5792  if (!onlyHooke) {
5793  feMaterialLhs->getOpPtrVector().push_back(
5794  new OpGetCrackFrontCommonDataAtGaussPts(
5795  "SPATIAL_POSITION", materialFe->commonData,
5796  feMaterialLhs->singularElement, feMaterialLhs->invSJac));
5797  feMaterialLhs->getOpPtrVector().push_back(
5798  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
5799  "MESH_NODE_POSITIONS", materialFe->commonData));
5800  feMaterialLhs->getOpPtrVector().push_back(
5801  new OpTransfromSingularBaseFunctions(
5802  feMaterialLhs->singularElement, feMaterialLhs->detS,
5803  feMaterialLhs
5804  ->invSJac //,AINSWORTH_LOBATTO_BASE//,AINSWORTH_LEGENDRE_BASE
5805  ));
5806  for (auto sit = materialFe->setOfBlocks.begin();
5807  sit != materialFe->setOfBlocks.end(); sit++) {
5808  feMaterialLhs->getOpPtrVector().push_back(
5809  new NonlinearElasticElement::OpJacobianEshelbyStress(
5810  "SPATIAL_POSITION", sit->second, materialFe->commonData,
5811  MATERIAL_TAG, true, true));
5812  feMaterialLhs->getOpPtrVector().push_back(
5813  new NonlinearElasticElement::OpLhsEshelby_dX(
5814  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", sit->second,
5815  materialFe->commonData));
5816  }
5817  } else { // HOOKE
5818  feMaterialLhs->getOpPtrVector().push_back(
5819  new OpCalculateVectorFieldGradient<3, 3>(
5820  "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5821  mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
5822  feMaterialLhs->getOpPtrVector().push_back(
5823  new OpGetCrackFrontDataGradientAtGaussPts(
5824  "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5825  feMaterialLhs->singularElement, feMaterialLhs->invSJac));
5826  space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
5827 
5828  if (defaultMaterial == BONEHOOKE) {
5829  if (!mwlsApprox) {
5830  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5831  "mwlsApprox not allocated");
5832  }
5833 
5834  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5835  feMaterialLhs->getOpPtrVector().push_back(
5836  new MWLSApprox::OpMWLSRhoAtGaussPts(
5837  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5838  mwlsApprox, mwlsRhoTagName, true, true));
5839  else {
5840  feMaterialLhs->getOpPtrVector().push_back(
5841  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
5842  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5843  mwlsApprox));
5844  feMaterialLhs->getOpPtrVector().push_back(
5845  new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5846  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5847  mwlsApprox, mwlsRhoTagName, true, true));
5848  }
5849 
5850  feMaterialLhs->getOpPtrVector().push_back(
5851  new HookeElement::OpCalculateStiffnessScaledByDensityField(
5852  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5853  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
5854  rHo0));
5855 
5856  feMaterialLhs->getOpPtrVector().push_back(
5857  new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5858  "MESH_NODE_POSITIONS",
5859  data_hooke_element_at_pts));
5860 
5861  feMaterialLhs->getOpPtrVector().push_back(
5862  new HookeElement::OpCalculateStress<1>("MESH_NODE_POSITIONS",
5863  "MESH_NODE_POSITIONS",
5864  data_hooke_element_at_pts));
5865 
5866  feMaterialLhs->getOpPtrVector().push_back(
5867  new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
5868  "MESH_NODE_POSITIONS",
5869  data_hooke_element_at_pts));
5870  feMaterialLhs->getOpPtrVector().push_back(
5871  new HookeElement::OpCalculateEshelbyStress(
5872  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5873  data_hooke_element_at_pts));
5874  feMaterialLhs->getOpPtrVector().push_back(
5875  new OpTransfromSingularBaseFunctions(feMaterialLhs->singularElement,
5876  feMaterialLhs->detS,
5877  feMaterialLhs->invSJac));
5878  feMaterialLhs->getOpPtrVector().push_back(
5879  new HookeElement::OpAleLhs_dX_dX<1>("MESH_NODE_POSITIONS",
5880  "MESH_NODE_POSITIONS",
5881  data_hooke_element_at_pts));
5882 
5883  feMaterialLhs->getOpPtrVector().push_back(
5884  new HookeElement::OpAleLhsWithDensity_dX_dX(
5885  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5886  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5887  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
5888 
5889  boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
5890  if (addSingularity) {
5891 
5892  mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
5893  mwlsApprox, &mwlsApprox->singularInitialDisplacement);
5894 
5895  feMaterialLhs->getOpPtrVector().push_back(
5896  new OpAleLhsWithDensitySingularElement_dX_dX(
5897  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5898  data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5899  mwlsApprox->diffRhoAtGaussPts, nBone, rHo0,
5900  mat_singular_disp_ptr));
5901  }
5902 
5903  feMaterialLhs->getOpPtrVector().push_back(
5904  new OpLhsBoneExplicitDerivariveWithHooke_dX(
5905  *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5906  mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5907  mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
5908  rHo0, nBone, &crackFrontNodes));
5909 
5910  } else {
5911  feMaterialLhs->getOpPtrVector().push_back(
5912  new HookeElement::OpCalculateHomogeneousStiffness<0>(
5913  "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5914  data_hooke_element_at_pts));
5915 
5916  feMaterialLhs->getOpPtrVector().push_back(
5917  new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5918  "MESH_NODE_POSITIONS",
5919  data_hooke_element_at_pts));
5920  feMaterialLhs->getOpPtrVector().push_back(
5921  new HookeElement::OpCalculateStress<0>("MESH_NODE_POSITIONS",
5922  "MESH_NODE_POSITIONS",
5923  data_hooke_element_at_pts));
5924 
5925  feMaterialLhs->getOpPtrVector().push_back(
5926  new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
5927  "MESH_NODE_POSITIONS",
5928  data_hooke_element_at_pts));
5929  feMaterialLhs->getOpPtrVector().push_back(
5930  new HookeElement::OpCalculateEshelbyStress(
5931  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5932  data_hooke_element_at_pts));
5933  feMaterialLhs->getOpPtrVector().push_back(
5934  new OpTransfromSingularBaseFunctions(feMaterialLhs->singularElement,
5935  feMaterialLhs->detS,
5936  feMaterialLhs->invSJac));
5937  feMaterialLhs->getOpPtrVector().push_back(
5938  new HookeElement::OpAleLhs_dX_dX<0>("MESH_NODE_POSITIONS",
5939  "MESH_NODE_POSITIONS",
5940  data_hooke_element_at_pts));
5941  }
5942  }
5943  if (residualStressBlock != -1) {
5944  if (!mwlsApprox) {
5945  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5946  "mwlsApprox not allocated");
5947  }
5948 
5949  if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5950  feMaterialLhs->getOpPtrVector().push_back(
5951  new MWLSApprox::OpMWLSStressAtGaussPts(
5952  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5953  mwlsApprox, mwlsStressTagName, true));
5954  else {
5955  feMaterialLhs->getOpPtrVector().push_back(
5956  new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
5957  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5958  mwlsApprox));
5959  feMaterialLhs->getOpPtrVector().push_back(
5960  new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
5961  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5962  mwlsApprox, mwlsStressTagName, false));
5963  }
5964 
5965  feMaterialLhs->getOpPtrVector().push_back(
5966  new MWLSApprox::OpMWLSMaterialStressLhs_DX(
5967  space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
5968  &crackFrontNodes));
5969  }
5970 
5971  MoFEMFunctionReturn(0);
5972 }

◆ assembleSmootherForcesDM()

MoFEMErrorCode FractureMechanics::CrackPropagation::assembleSmootherForcesDM ( DM  dm,
const std::vector< int >  ids,
const int  verb = QUIET,
const bool  debug = false 
)

create smoothing element instance

Definition at line 5975 of file CrackPropagation.cpp.

5976  {
5977  MoFEMFunctionBeginHot;
5978  if (!elasticFe)
5979  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
5980  "Elastic element instance not created");
5981 
5982  smootherFe = boost::make_shared<Smoother>(mField);
5983  volumeLengthAdouble = boost::shared_ptr<VolumeLengthQuality<adouble>>(
5984  new VolumeLengthQuality<adouble>(
5985  BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME, smootherAlpha,
5986  smootherGamma));
5987  volumeLengthDouble = boost::shared_ptr<VolumeLengthQuality<double>>(
5988  new VolumeLengthQuality<double>(
5989  BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME, smootherAlpha,
5990  smootherGamma));
5991 
5992  // set block data
5993  for (std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
5994  elasticFe->setOfBlocks.begin();
5995  sit != elasticFe->setOfBlocks.end(); sit++) {
5996  // E = fmax(E,sit->second.E);
5997  smootherFe->setOfBlocks[0].tEts.merge(sit->second.tEts);
5998  }
5999  smootherFe->setOfBlocks[0].materialDoublePtr = volumeLengthDouble;
6000  smootherFe->setOfBlocks[0].materialAdoublePtr = volumeLengthAdouble;
6001  smootherFe->setOfBlocks[0].forcesOnlyOnEntitiesRow = crackFrontNodes;
6002 
6003  // set element data
6004  smootherFe->commonData.spatialPositions = "MESH_NODE_POSITIONS";
6005  smootherFe->commonData.meshPositions = "NONE";
6006 
6007  // mesh field name
6008  smootherFe->feRhs.meshPositionsFieldName = "NONE";
6009  smootherFe->feLhs.meshPositionsFieldName = "NONE";
6010  smootherFe->feRhs.addToRule = 0;
6011  smootherFe->feLhs.addToRule = 0;
6012 
6013  // Create finite element instances for the right hand side and left hand side
6014  feSmootherRhs = smootherFe->feRhsPtr;
6015  feSmootherLhs = smootherFe->feLhsPtr;
6016 
6017  // Smoother right hand side
6018  feSmootherRhs->getOpPtrVector().push_back(
6019  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
6020  "MESH_NODE_POSITIONS", smootherFe->commonData));
6021  map<int, NonlinearElasticElement::BlockData>::iterator sit =
6022  smootherFe->setOfBlocks.begin();
6023  feSmootherRhs->getOpPtrVector().push_back(new Smoother::OpJacobianSmoother(
6024  "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
6025  smootherFe->commonData, SMOOTHING_TAG, false));
6026  feSmootherRhs->getOpPtrVector().push_back(new Smoother::OpRhsSmoother(
6027  "MESH_NODE_POSITIONS", sit->second, smootherFe->commonData,
6028  smootherFe->smootherData));
6029 
6030  // Smoother left hand side
6031  feSmootherLhs->getOpPtrVector().push_back(
6032  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
6033  "MESH_NODE_POSITIONS", smootherFe->commonData));
6034  feSmootherLhs->getOpPtrVector().push_back(new Smoother::OpJacobianSmoother(
6035  "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
6036  smootherFe->commonData, SMOOTHING_TAG, true));
6037  feSmootherLhs->getOpPtrVector().push_back(new Smoother::OpLhsSmoother(
6038  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6039  smootherFe->setOfBlocks.at(0), smootherFe->commonData,
6040  smootherFe->smootherData, "LAMBDA_CRACKFRONT_AREA_TANGENT"));
6041 
6042  // Constrains at crack front
6043  tangentConstrains = boost::shared_ptr<
6044  ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain>(
6045  new ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain(
6046  mField, "LAMBDA_CRACKFRONT_AREA_TANGENT"));
6047 
6048  Range contact_faces;
6049  contact_faces.merge(contactSlaveFaces);
6050  contact_faces.merge(contactMasterFaces);
6051 
6052  // Surface sliding constrains
6053  surfaceConstrain.clear();
6054  skinOrientation =
6055  boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
6056  new FaceOrientation(false, contact_faces, mapBitLevel["mesh_cut"]));
6057  for (auto id : ids) {
6058  if (id != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6059  surfaceConstrain[id] = boost::make_shared<SurfaceSlidingConstrains>(
6060  mField, *skinOrientation);
6061  surfaceConstrain[id]->setOperators(
6062  SURFACE_SLIDING_TAG,
6063  "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(id),
6064  "MESH_NODE_POSITIONS", &smootherAlpha);
6065  }
6066  }
6067  // Add crack surface sliding constrains
6068  crackOrientation =
6069  boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
6070  new FaceOrientation(true, contact_faces, mapBitLevel["mesh_cut"]));
6071  surfaceConstrain[getInterface<CPMeshCut>()->getCrackSurfaceId()] =
6072  boost::make_shared<SurfaceSlidingConstrains>(mField, *crackOrientation);
6073  surfaceConstrain[getInterface<CPMeshCut>()->getCrackSurfaceId()]
6074  ->setOperators(SURFACE_SLIDING_TAG,
6075  "LAMBDA_SURFACE" +
6076  boost::lexical_cast<std::string>(
6077  getInterface<CPMeshCut>()->getCrackSurfaceId()),
6078  "MESH_NODE_POSITIONS", &smootherAlpha);
6079 
6080  edgeConstrains.clear();
6081  auto get_edges_block_set = [this]() {
6082  return getInterface<CPMeshCut>()->getEdgesBlockSet();
6083  };
6084  if (get_edges_block_set()) {
6085  edgeConstrains[get_edges_block_set()] =
6086  boost::shared_ptr<EdgeSlidingConstrains>(
6087  new EdgeSlidingConstrains(mField));
6088  edgeConstrains[get_edges_block_set()]->setOperators(
6089  EDGE_SLIDING_TAG, "LAMBDA_EDGE", "MESH_NODE_POSITIONS", &smootherAlpha);
6090  }
6091 
6092  bothSidesConstrains = boost::shared_ptr<BothSurfaceConstrains>(
6093  new BothSurfaceConstrains(mField));
6094  bothSidesConstrains->aLpha = smootherAlpha;
6095 
6096  bothSidesContactConstrains = boost::shared_ptr<BothSurfaceConstrains>(
6097  new BothSurfaceConstrains(mField, "LAMBDA_BOTH_SIDES_CONTACT"));
6098  bothSidesContactConstrains->aLpha = smootherAlpha;
6099  bothSidesContactConstrains->masterNodes = contactBothSidesMasterNodes;
6100 
6101  Range fix_material_nodes;
6102  fixMaterialEnts = boost::make_shared<DirichletFixFieldAtEntitiesBc>(
6103  mField, "MESH_NODE_POSITIONS", fix_material_nodes);
6104  const Field_multiIndex *fields_ptr;
6105  CHKERR mField.get_fields(&fields_ptr);
6106  for (auto f : *fields_ptr) {
6107  if (f->getName().compare(0, 6, "LAMBDA") == 0 &&
6108  f->getName() != "LAMBDA_ARC_LENGTH" &&
6109  f->getName() != "LAMBDA_CONTACT" &&
6110  f->getName() != "LAMBDA_CLOSE_CRACK") {
6111  fixMaterialEnts->fieldNames.push_back(f->getName());
6112  }
6113  }
6114 
6115  // Create griffith force finite element instances and operators
6116  griffithForceElement = boost::make_shared<GriffithForceElement>(mField);
6117  griffithForceElement->blockData[0].gc = gC;
6118  griffithForceElement->blockData[0].E = gC;
6119  griffithForceElement->blockData[0].r = 1.0;
6120  griffithForceElement->blockData[0].frontNodes = crackFrontNodes;
6121  CHKERR GriffithForceElement::getElementOptions(
6122  griffithForceElement->blockData[0]);
6123  gC = griffithForceElement->blockData[0].gc;
6124  griffithE = griffithForceElement->blockData[0].E;
6125  griffithR = griffithForceElement->blockData[0].r;
6126 
6127  feGriffithForceRhs = griffithForceElement->feRhsPtr;
6128  feGriffithForceLhs = griffithForceElement->feLhsPtr;
6129 
6130  feGriffithForceRhs->getOpPtrVector().push_back(
6131  new GriffithForceElement::OpCalculateGriffithForce(
6132  GRIFFITH_FORCE_TAG, griffithForceElement->blockData[0],
6133  griffithForceElement->commonData));
6134  feGriffithForceRhs->getOpPtrVector().push_back(
6135  new GriffithForceElement::OpRhs(GRIFFITH_FORCE_TAG,
6136  griffithForceElement->blockData[0],
6137  griffithForceElement->commonData));
6138  feGriffithForceLhs->getOpPtrVector().push_back(
6139  new GriffithForceElement::OpLhs(GRIFFITH_FORCE_TAG,
6140  griffithForceElement->blockData[0],
6141  griffithForceElement->commonData));
6142 
6143  // Creating Griffith constrains finite element instances and operators
6144  feGriffithConstrainsDelta =
6145  boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrainsDelta>(
6146  new GriffithForceElement::MyTriangleFEConstrainsDelta(
6147  mField, "LAMBDA_CRACKFRONT_AREA"));
6148  feGriffithConstrainsDelta->getOpPtrVector().push_back(
6149  new GriffithForceElement::OpConstrainsDelta(
6150  mField, GRIFFITH_CONSTRAINS_TAG, griffithForceElement->blockData[0],
6151  griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA",
6152  feGriffithConstrainsDelta->deltaVec));
6153  feGriffithConstrainsRhs =
6154  boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>(
6155  new GriffithForceElement::MyTriangleFEConstrains(
6156  mField, "LAMBDA_CRACKFRONT_AREA",
6157  griffithForceElement->blockData[0],
6158  feGriffithConstrainsDelta->deltaVec));
6159  feGriffithConstrainsLhs =
6160  boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>(
6161  new GriffithForceElement::MyTriangleFEConstrains(
6162  mField, "LAMBDA_CRACKFRONT_AREA",
6163  griffithForceElement->blockData[0],
6164  feGriffithConstrainsDelta->deltaVec));
6165  feGriffithConstrainsRhs->getOpPtrVector().push_back(
6166  new GriffithForceElement::OpConstrainsRhs(
6167  mField, GRIFFITH_CONSTRAINS_TAG, griffithForceElement->blockData[0],
6168  griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA"));
6169  feGriffithConstrainsLhs->getOpPtrVector().push_back(
6170  new GriffithForceElement::OpConstrainsLhs(
6171  mField, GRIFFITH_CONSTRAINS_TAG, griffithForceElement->blockData[0],
6172  griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA",
6173  feGriffithConstrainsLhs->deltaVec));
6174 
6175  MoFEMFunctionReturnHot(0);
6176 }

◆ buildArcLengthField()

MoFEMErrorCode FractureMechanics::CrackPropagation::buildArcLengthField ( const BitRefLevel  bit,
const bool  build_fields = true,
const int  verb = QUIET 
)

Declate field for arc-length.

Parameters
bitbit refinement level
Returns
error code

Definition at line 1726 of file CrackPropagation.cpp.

1728  {
1729  const RefEntity_multiIndex *refined_ents_ptr;
1730  const FieldEntity_multiIndex *field_ents_ptr;
1731  MoFEMFunctionBegin;
1732  CHKERR mField.get_ref_ents(&refined_ents_ptr);
1733  CHKERR mField.get_field_ents(&field_ents_ptr);
1734  if (mField.check_field("LAMBDA_ARC_LENGTH")) {
1735  auto vit = refined_ents_ptr->get<Ent_mi_tag>().find(arcLengthVertex);
1736  if (vit == refined_ents_ptr->get<Ent_mi_tag>().end()) {
1737  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1738  "Vertex of arc-length field not found");
1739  }
1740  *(const_cast<RefEntity *>(vit->get())->getBitRefLevelPtr()) |= bit;
1741  auto fit = field_ents_ptr->get<Unique_mi_tag>().find(
1742  FieldEntity::getLocalUniqueIdCalculate(
1743  mField.get_field_bit_number("LAMBDA_ARC_LENGTH"), arcLengthVertex));
1744  if (fit == field_ents_ptr->get<Unique_mi_tag>().end())
1745  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1746  "Arc length field entity not found");
1747 
1748  MOFEM_LOG("CPWorld", Sev::noisy) << **fit;
1749  MOFEM_LOG_C("CPWorld", Sev::inform, "Arc-Length field lambda = %3.4e",
1750  (*fit)->getEntFieldData()[0]);
1751  MoFEMFunctionReturnHot(0);
1752  }
1753 
1754  CHKERR mField.add_field("LAMBDA_ARC_LENGTH", NOFIELD, NOBASE, 1,
1755  MB_TAG_SPARSE, MF_ZERO, verb);
1756  EntityHandle lambda_meshset = mField.get_field_meshset("LAMBDA_ARC_LENGTH");
1757  CHKERR mField.get_moab().add_entities(lambda_meshset, &arcLengthVertex, 1);
1758  // Add vertex to MoFEM database
1759  std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
1760  const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
1761  ->insert(boost::make_shared<RefEntity>(
1762  mField.get_basic_entity_data_ptr(), arcLengthVertex));
1763  *(const_cast<RefEntity *>(p_ent.first->get())->getBitRefLevelPtr()) |=
1764  (bit | BitRefLevel().set(BITREFLEVEL_SIZE - 2));
1765  if (build_fields) {
1766  CHKERR mField.build_fields(verb);
1767  }
1768  MoFEMFunctionReturn(0);
1769 }

◆ buildBothSidesFieldId()

MoFEMErrorCode FractureMechanics::CrackPropagation::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.

Nodes on both sides are constrainded to have the same material displacements

Definition at line 1973 of file CrackPropagation.cpp.

1976  {
1977  MoFEMFunctionBegin;
1978 
1979  auto add_both_sides_field = [&](const std::string lambda_field_name,
1980  Range master_nodes, Range nodes_to_remove,
1981  bool add_ho) {
1982  MoFEMFunctionBegin;
1983  // This field set constrains that both side of crack surface have the same
1984  // node positions
1985  CHKERR mField.add_field(lambda_field_name, H1, AINSWORTH_LEGENDRE_BASE, 3,
1986  MB_TAG_SPARSE, MF_ZERO);
1987 
1988  if (debug) {
1989  std::ostringstream file;
1990  file << lambda_field_name << ".vtk";
1991  EntityHandle meshset;
1992  CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
1993  CHKERR mField.get_moab().add_entities(meshset, master_nodes);
1994  CHKERR mField.get_moab().write_file(file.str().c_str(), "VTK", "",
1995  &meshset, 1);
1996  CHKERR mField.get_moab().delete_entities(&meshset, 1);
1997  }
1998 
1999  // remove old ents from field which are no longer part of the field
2000  Range ents_to_remove;
2001  CHKERR mField.get_field_entities_by_handle(lambda_field_name,
2002  ents_to_remove);
2003  ents_to_remove = subtract(ents_to_remove, master_nodes);
2004  CHKERR mField.remove_ents_from_field(lambda_field_name, ents_to_remove,
2005  verb);
2006 
2007  CHKERR mField.add_ents_to_field_by_type(master_nodes, MBVERTEX,
2008  lambda_field_name);
2009  if (approxOrder > 1 && add_ho) {
2010  CHKERR mField.add_ents_to_field_by_type(master_nodes, MBEDGE,
2011  lambda_field_name);
2012  CHKERR mField.add_ents_to_field_by_type(master_nodes, MBTRI,
2013  lambda_field_name);
2014  }
2015  CHKERR mField.set_field_order(master_nodes.subset_by_dimension(0),
2016  lambda_field_name, 1);
2017  if (approxOrder > 1 && add_ho) {
2018  CHKERR mField.set_field_order(master_nodes.subset_by_dimension(1),
2019  lambda_field_name, approxOrder);
2020  CHKERR mField.set_field_order(master_nodes.subset_by_dimension(2),
2021  lambda_field_name, approxOrder);
2022  }
2023 
2024  // Since we state constrains that mesh nodes, in material space, have to
2025  // move the same on the top and bottom side of crack surface, we have to
2026  // remove constraint from the all surface constrains fields, except
2027  // "LAMBDA_SURFACE", otherwise system will be over-constrained.
2028  const Field_multiIndex *fields_ptr;
2029  CHKERR mField.get_fields(&fields_ptr);
2030  for (auto const &field : (*fields_ptr)) {
2031  const std::string field_name = field->getName();
2032  if (field_name.compare(0, 14, "LAMBDA_SURFACE") == 0) {
2033  CHKERR mField.remove_ents_from_field(field_name, nodes_to_remove, verb);
2034  }
2035  if (field_name.compare(0, 11, "LAMBDA_EDGE") == 0) {
2036  CHKERR mField.remove_ents_from_field(field_name, nodes_to_remove, verb);
2037  }
2038  }
2039  MoFEMFunctionReturn(0);
2040  };
2041 
2042  auto write_file = [&](const std::string file_name, Range ents) {
2043  MoFEMFunctionBegin;
2044  EntityHandle meshset;
2045  CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
2046  CHKERR mField.get_moab().add_entities(meshset, ents);
2047  if (mField.get_comm_rank() == 0) {
2048  CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "",
2049  &meshset, 1);
2050  }
2051  CHKERR mField.get_moab().delete_entities(&meshset, 1);
2052  MoFEMFunctionReturn(0);
2053  };
2054 
2055  auto get_prisms_level_nodes = [&](auto bit) {
2056  Range prisms_level;
2057  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2058  bit, BitRefLevel().set(), MBPRISM, prisms_level);
2059  if (proc_only) {
2060  prisms_level = intersect(bitProcEnts, prisms_level);
2061  }
2062  Range prisms_level_nodes;
2063  CHKERR mField.get_moab().get_connectivity(prisms_level, prisms_level_nodes,
2064  true);
2065  return prisms_level_nodes;
2066  };
2067 
2068  auto get_other_side_crack_faces_nodes = [&](auto prisms_level_nodes) {
2069  Range other_side_crack_faces_nodes;
2070  CHKERR mField.get_moab().get_connectivity(
2071  otherSideCrackFaces, other_side_crack_faces_nodes, true);
2072  other_side_crack_faces_nodes =
2073  intersect(other_side_crack_faces_nodes, prisms_level_nodes);
2074  other_side_crack_faces_nodes =
2075  subtract(other_side_crack_faces_nodes, crackFrontNodes);
2076 
2077  return other_side_crack_faces_nodes;
2078  };
2079 
2080  auto prisms_level_nodes = get_prisms_level_nodes(bit_material);
2081  auto other_side_crack_faces_nodes =
2082  get_other_side_crack_faces_nodes(prisms_level_nodes);
2083 
2084  CHKERR add_both_sides_field("LAMBDA_BOTH_SIDES", other_side_crack_faces_nodes,
2085  other_side_crack_faces_nodes, false);
2086 
2087  if (solveEigenStressProblem) {
2088  Range other_side_crack_faces_edges;
2089  if (approxOrder > 1) {
2090  CHKERR mField.get_moab().get_adjacencies(otherSideCrackFaces, 1, false,
2091  other_side_crack_faces_edges,
2092  moab::Interface::UNION);
2093  other_side_crack_faces_edges =
2094  subtract(other_side_crack_faces_edges, crackFront);
2095  other_side_crack_faces_edges.merge(otherSideCrackFaces);
2096  }
2097  other_side_crack_faces_edges.merge(
2098  get_other_side_crack_faces_nodes(get_prisms_level_nodes(bit_spatial)));
2099  CHKERR add_both_sides_field("LAMBDA_CLOSE_CRACK",
2100  other_side_crack_faces_edges, Range(), true);
2101  }
2102 
2103  if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
2104  contactBothSidesMasterNodes.clear();
2105  CHKERR mField.get_moab().get_connectivity(
2106  contactBothSidesMasterFaces, contactBothSidesMasterNodes, true);
2107  contactBothSidesMasterNodes =
2108  intersect(contactBothSidesMasterNodes, prisms_level_nodes);
2109  contactBothSidesMasterNodes =
2110  subtract(contactBothSidesMasterNodes, other_side_crack_faces_nodes);
2111  Range contact_both_sides_slave_nodes;
2112  for (Range::iterator nit = contactBothSidesMasterNodes.begin();
2113  nit != contactBothSidesMasterNodes.end(); nit++) {
2114  Range contact_prisms;
2115  CHKERR mField.get_moab().get_adjacencies(
2116  &*nit, 1, 3, false, contact_prisms, moab::Interface::UNION);
2117  contact_prisms = intersect(contact_prisms, contactElements);
2118  EntityHandle prism = contact_prisms.front();
2119  const EntityHandle *conn;
2120  int side_number, other_side_number, sense, offset, number_nodes = 0;
2121  CHKERR mField.get_moab().get_connectivity(prism, conn, number_nodes);
2122  CHKERR mField.get_moab().side_number(prism, *nit, side_number, sense,
2123  offset);
2124  if (side_number < 3)
2125  contact_both_sides_slave_nodes.insert(conn[side_number + 3]);
2126  else
2127  contact_both_sides_slave_nodes.insert(conn[side_number - 3]);
2128  }
2129  contact_both_sides_slave_nodes =
2130  intersect(contact_both_sides_slave_nodes, prisms_level_nodes);
2131 
2132  if (debug) {
2133  CHKERR write_file("contact_both_sides_slave_nodes.vtk",
2134  contact_both_sides_slave_nodes);
2135  }
2136 
2137  CHKERR add_both_sides_field("LAMBDA_BOTH_SIDES_CONTACT",
2138  contactBothSidesMasterNodes,
2139  contact_both_sides_slave_nodes, false);
2140  }
2141 
2142  MoFEMFunctionReturn(0);
2143 }

◆ buildCrackFrontFieldId()

MoFEMErrorCode FractureMechanics::CrackPropagation::buildCrackFrontFieldId ( const BitRefLevel  bit,
const bool  build_fields = true,
const int  verb = QUIET,
const bool  debug = false 
)

declare crack surface files

Definition at line 2145 of file CrackPropagation.cpp.

2148  {
2149  moab::Interface &moab = mField.get_moab();
2150  MoFEMFunctionBegin;
2151 
2152  Range level_edges;
2153  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2154  bit, BitRefLevel().set(), MBEDGE, level_edges);
2155  Range crack_front_edges = intersect(crackFront, level_edges);
2156  Range crack_front_nodes;
2157  CHKERR moab.get_connectivity(crack_front_edges, crack_front_nodes, true);
2158 
2159  CHKERR mField.add_field("LAMBDA_CRACKFRONT_AREA", H1, AINSWORTH_LEGENDRE_BASE,
2160  1, MB_TAG_SPARSE, MF_ZERO, verb);
2161  // Add dofs to LAMBDA_CRACKFRONT_AREA field
2162  {
2163  Range ents_to_remove;
2164  CHKERR mField.get_field_entities_by_handle("LAMBDA_CRACKFRONT_AREA",
2165  ents_to_remove);
2166  ents_to_remove = subtract(ents_to_remove, crack_front_nodes);
2167  CHKERR mField.remove_ents_from_field("LAMBDA_CRACKFRONT_AREA",
2168  ents_to_remove, verb);
2169  CHKERR mField.add_ents_to_field_by_type(crack_front_nodes, MBVERTEX,
2170  "LAMBDA_CRACKFRONT_AREA", verb);
2171  CHKERR mField.set_field_order(crack_front_nodes, "LAMBDA_CRACKFRONT_AREA",
2172  1, verb);
2173  }
2174 
2175  CHKERR mField.add_field("LAMBDA_CRACKFRONT_AREA_TANGENT", H1,
2176  AINSWORTH_LEGENDRE_BASE, 1, MB_TAG_SPARSE, MF_ZERO,
2177  verb);
2178  // Add dofs to LAMBDA_CRACKFRONT_AREA_TANGENT field
2179  {
2180  // Remove dofs on crack front ends
2181  Skinner skin(&mField.get_moab());
2182  Range skin_crack_front;
2183  rval = skin.find_skin(0, crack_front_edges, false, skin_crack_front);
2184  Range field_ents = subtract(crack_front_nodes, skin_crack_front);
2185  // remove old ents from field which are no longer part of the field
2186  Range ents_to_remove;
2187  CHKERR mField.get_field_entities_by_handle("LAMBDA_CRACKFRONT_AREA_TANGENT",
2188  ents_to_remove);
2189  ents_to_remove = subtract(ents_to_remove, field_ents);
2190  CHKERR mField.remove_ents_from_field("LAMBDA_CRACKFRONT_AREA_TANGENT",
2191  ents_to_remove, verb);
2192  CHKERR mField.add_ents_to_field_by_type(
2193  field_ents, MBVERTEX, "LAMBDA_CRACKFRONT_AREA_TANGENT", verb);
2194  CHKERR mField.set_field_order(field_ents, "LAMBDA_CRACKFRONT_AREA_TANGENT",
2195  1, verb);
2196  }
2197 
2198  if (build_fields) {
2199  CHKERR mField.