v0.13.1
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
 
doublegetLoadScale ()
 
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 refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

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::InterfacemField
 
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...
 
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
 

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< FEMethodgetFrontArcLengthControl (boost::shared_ptr< ArcLengthCtx > arc_ctx)
 
boost::shared_ptr< AnalyticalDirichletBC::DirichletBCgetAnalyticalDirichletBc ()
 
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< NonlinearElasticElementelasticFe
 
boost::shared_ptr< NonlinearElasticElementmaterialFe
 
boost::shared_ptr< CrackFrontElementfeLhs
 Integrate elastic FE. More...
 
boost::shared_ptr< CrackFrontElementfeRhs
 Integrate elastic FE. More...
 
boost::shared_ptr< CrackFrontElementfeMaterialRhs
 Integrate material stresses, assemble vector. More...
 
boost::shared_ptr< CrackFrontElementfeMaterialLhs
 Integrate material stresses, assemble matrix. More...
 
boost::shared_ptr< CrackFrontElementfeEnergy
 Integrate energy. More...
 
boost::shared_ptr< ConstrainMatrixCtxprojSurfaceCtx
 Data structure to project on the body surface. More...
 
boost::shared_ptr< ConstrainMatrixCtxprojFrontCtx
 Data structure to project on crack front. More...
 
boost::shared_ptr< MWLSApproxmwlsApprox
 
boost::shared_ptr< NeumannForcesSurface::MyTriangleFEfeMaterialAnaliticalTraction
 Surface elment to calculate tractions in material space. More...
 
boost::shared_ptr< DirichletSpatialPositionsBcspatialDirichletBc
 apply Dirichlet BC to sparial positions More...
 
boost::shared_ptr< AnalyticalDirichletBC::DirichletBCanalyticalDirichletBc
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfaceForces
 assemble surface forces 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::DataAtIntegrationPtscommonDataSurfacePressureAle
 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< FEMethodassembleFlambda
 assemble F_lambda vector More...
 
boost::shared_ptr< FEMethodzeroFlambda
 assemble F_lambda vector More...
 
boost::shared_ptr< CrackFrontElementfeCouplingElasticLhs
 FE instance to assemble coupling terms. More...
 
boost::shared_ptr< CrackFrontElementfeCouplingMaterialLhs
 FE instance to assemble coupling terms. More...
 
boost::shared_ptr< SmoothersmootherFe
 
boost::shared_ptr< Smoother::MyVolumeFEfeSmootherRhs
 Integrate smoothing operators. More...
 
boost::shared_ptr< Smoother::MyVolumeFEfeSmootherLhs
 Integrate smoothing operators. More...
 
boost::shared_ptr< VolumeLengthQuality< double > > volumeLengthDouble
 
boost::shared_ptr< VolumeLengthQuality< adouble > > volumeLengthAdouble
 
boost::shared_ptr< ObosleteUsersModules::TangentWithMeshSmoothingFrontConstraintangentConstrains
 Constrains crack front in tangent directiona. More...
 
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientationskinOrientation
 
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientationcrackOrientation
 
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientationcontactOrientation
 
map< int, boost::shared_ptr< SurfaceSlidingConstrains > > surfaceConstrain
 
map< int, boost::shared_ptr< EdgeSlidingConstrains > > edgeConstrains
 
boost::shared_ptr< BothSurfaceConstrainsbothSidesConstrains
 
boost::shared_ptr< BothSurfaceConstrainscloseCrackConstrains
 
boost::shared_ptr< BothSurfaceConstrainsbothSidesContactConstrains
 
boost::shared_ptr< DirichletFixFieldAtEntitiesBcfixMaterialEnts
 
boost::shared_ptr< GriffithForceElementgriffithForceElement
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEfeGriffithForceRhs
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEfeGriffithForceLhs
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsDeltafeGriffithConstrainsDelta
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsfeGriffithConstrainsRhs
 
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsfeGriffithConstrainsLhs
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > analiticalSurfaceElement
 
boost::shared_ptr< FaceElementForcesAndSourcesCorefeSpringLhsPtr
 
boost::shared_ptr< FaceElementForcesAndSourcesCorefeSpringRhsPtr
 

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, doublemapG1
 hashmap of g1 - release energy at nodes More...
 
map< EntityHandle, doublemapG3
 hashmap of g3 - release energy at nodes More...
 
map< EntityHandle, doublemapJ
 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, doublemapSmoothingForceFactor
 
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< ArcLengthCtxarcCtx
 
boost::shared_ptr< ArcLengthCtxarcEigenCtx
 

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_multiIndexcontactSearchMultiIndexPtr
 
Range contactBothSidesMasterFaces
 
Range contactBothSidesSlaveFaces
 
Range contactBothSidesMasterNodes
 
Range contactTets
 
EntityHandle meshsetFaces
 
boost::shared_ptr< FaceElementForcesAndSourcesCorefeLhsSpringALE
 
boost::shared_ptr< FaceElementForcesAndSourcesCorefeLhsSpringALEMaterial
 
boost::shared_ptr< FaceElementForcesAndSourcesCorefeRhsSpringALEMaterial
 
boost::shared_ptr< MetaSpringBC::DataAtIntegrationPtsSpringscommonDataSpringsALE
 
boost::shared_ptr< SimpleContactProblemcontactProblem
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElementfeRhsSimpleContact
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElementfeLhsSimpleContact
 
boost::shared_ptr< SimpleContactProblem::CommonDataSimpleContactcommonDataSimpleContact
 
boost::shared_ptr< MortarContactProblemmortarContactProblemPtr
 
boost::shared_ptr< MortarContactProblem::MortarContactElementfeRhsMortarContact
 
boost::shared_ptr< MortarContactProblem::MortarContactElementfeLhsMortarContact
 
boost::shared_ptr< MortarContactProblem::CommonDataMortarContactcommonDataMortarContact
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElementfeRhsSimpleContactALEMaterial
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElementfeLhsSimpleContactALEMaterial
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElementfeLhsSimpleContactALE
 
boost::shared_ptr< SimpleContactProblem::CommonDataSimpleContactcommonDataSimpleContactALE
 
PetscBool contactOutputIntegPts
 
boost::shared_ptr< SimpleContactProblem::SimpleContactElementfePostProcSimpleContact
 
boost::shared_ptr< MortarContactProblem::MortarContactElementfePostProcMortarContact
 
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
 

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

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 515 of file CrackPropagation.cpp.

520 geometryOrder(geometry_order), gC(0), betaGc(0),
521 isGcHeterogeneous(PETSC_FALSE), isPartitioned(PETSC_FALSE),
522 propagateCrack(PETSC_FALSE), refAtCrackTip(0), refOrderAtTip(0),
526 nBone(0), cpSolversPtr(new CPSolvers(*this)),
527 cpMeshCutPtr(new CPMeshCut(*this)) {
528
529 if (!LogManager::checkIfChannelExist("CPWorld")) {
530 auto core_log = logging::core::get();
531
532 core_log->add_sink(
533 LogManager::createSink(LogManager::getStrmWorld(), "CPWorld"));
534 core_log->add_sink(
535 LogManager::createSink(LogManager::getStrmSync(), "CPSync"));
536 core_log->add_sink(
537 LogManager::createSink(LogManager::getStrmSelf(), "CPSelf"));
538
539 LogManager::setLog("CPWorld");
540 LogManager::setLog("CPSync");
541 LogManager::setLog("CPSelf");
542
543 MOFEM_LOG_TAG("CPWorld", "CP");
544 MOFEM_LOG_TAG("CPSync", "CP");
545 MOFEM_LOG_TAG("CPSelf", "CP");
546 }
547
548 MOFEM_LOG("CPWorld", Sev::noisy) << "CPSolve created";
549
550#ifdef GIT_UM_SHA1_NAME
551 MOFEM_LOG_C("CPWorld", Sev::inform, "User module git commit id %s",
552 GIT_UM_SHA1_NAME);
553#endif
554
555 Version version;
556 getInterfaceVersion(version);
557 MOFEM_LOG("CPWorld", Sev::inform)
558 << "Fracture module version " << version.strVersion();
559
560#ifdef GIT_FM_SHA1_NAME
561 MOFEM_LOG_C("CPWorld", Sev::inform, "Fracture module git commit id %s",
562 GIT_FM_SHA1_NAME);
563#endif
564
565 ierr = registerInterface<CrackPropagation>();
566 CHKERRABORT(PETSC_COMM_SELF, ierr);
567 ierr = registerInterface<CPSolvers>();
568 CHKERRABORT(PETSC_COMM_SELF, ierr);
569 ierr = registerInterface<CPMeshCut>();
570 CHKERRABORT(PETSC_COMM_SELF, ierr);
571
572 mapBitLevel["mesh_cut"] = BitRefLevel().set(0);
573 mapBitLevel["spatial_domain"] = BitRefLevel().set(1);
574 mapBitLevel["material_domain"] = BitRefLevel().set(2);
575
576 if (!moabCommWorld)
577 moabCommWorld = boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
578}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
static constexpr int approx_order
PetscBool propagateCrack
If true crack propagation is calculated.
std::map< std::string, BitRefLevel > mapBitLevel
double rHo0
Reference density if bone is analyzed.
MoFEMErrorCode getInterfaceVersion(Version &version) const
double betaGc
heterogeneous Griffith energy exponent
PetscBool isGcHeterogeneous
flag for heterogeneous gc
const boost::scoped_ptr< UnknownInterface > cpSolversPtr
double nBone
Exponent parameter in bone density.
const boost::scoped_ptr< UnknownInterface > cpMeshCutPtr
int postProcLevel
level of postprocessing (amount of output files)
boost::shared_ptr< WrapMPIComm > moabCommWorld
int residualStressBlock
Block on which residual stress is applied.
std::string strVersion()

◆ ~CrackPropagation()

FractureMechanics::CrackPropagation::~CrackPropagation ( )
virtual

Definition at line 580 of file CrackPropagation.cpp.

580{}

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 4487 of file CrackPropagation.cpp.

4489 {
4490 boost::shared_ptr<FEMethod> null;
4492
4493 // Assemble F_lambda force
4494 assembleFlambda = boost::shared_ptr<FEMethod>(
4495 new AssembleFlambda(arc_ctx, spatialDirichletBc));
4496 zeroFlambda = boost::shared_ptr<FEMethod>(new ZeroFLmabda(arc_ctx));
4497
4498 if (mwlsApprox) {
4499 mwlsApprox->F_lambda = arc_ctx->F_lambda;
4501 mwlsApprox->arcLengthDof = arcLengthDof;
4502 }
4503 // Set up DM specific vectors and data
4504 spatialDirichletBc->mapZeroRows.clear();
4505
4506 VecSetOption(arc_ctx->F_lambda, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
4507
4508 // Set-up F_lambda to elements & operator evalulating forces
4509 // Surface force
4510 for (auto fit = surfaceForces->begin(); fit != surfaceForces->end(); fit++) {
4511 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4512 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4513 for (; oit != hi_oit; oit++) {
4514 if (boost::typeindex::type_id_runtime(*oit) ==
4515 boost::typeindex::type_id<NeumannForcesSurface::OpNeumannForce>()) {
4516 dynamic_cast<NeumannForcesSurface::OpNeumannForce &>(*oit).F =
4517 arc_ctx->F_lambda;
4518 }
4519 }
4520 }
4521 for (auto fit = surfacePressure->begin(); fit != surfacePressure->end();
4522 fit++) {
4523 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4524 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4525 for (; oit != hi_oit; oit++) {
4526 if (boost::typeindex::type_id_runtime(*oit) ==
4527 boost::typeindex::type_id<
4529 dynamic_cast<NeumannForcesSurface::OpNeumannPressure &>(*oit).F =
4530 arc_ctx->F_lambda;
4531 }
4532 if (boost::typeindex::type_id_runtime(*oit) ==
4533 boost::typeindex::type_id<
4536 arc_ctx->F_lambda;
4537 }
4538 }
4539 }
4540 for (auto fit = edgeForces->begin(); fit != edgeForces->end(); fit++) {
4541 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4542 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4543 for (; oit != hi_oit; oit++) {
4544 if (boost::typeindex::type_id_runtime(*oit) ==
4545 boost::typeindex::type_id<EdgeForce::OpEdgeForce>()) {
4546 dynamic_cast<EdgeForce::OpEdgeForce &>(*oit).F = arc_ctx->F_lambda;
4547 }
4548 }
4549 }
4550 for (auto fit = nodalForces->begin(); fit != nodalForces->end(); fit++) {
4551 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4552 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4553 for (; oit != hi_oit; oit++) {
4554 if (boost::typeindex::type_id_runtime(*oit) ==
4555 boost::typeindex::type_id<NodalForce::OpNodalForce>()) {
4556 dynamic_cast<NodalForce::OpNodalForce &>(*oit).F = arc_ctx->F_lambda;
4557 }
4558 }
4559 }
4560
4561 // Lhs
4563 null);
4564#ifdef __ANALITICAL_DISPLACEMENT__
4565 if (mField.check_problem("BC_PROBLEM")) {
4567 analyticalDirichletBc, null);
4568 }
4569#endif //__ANALITICAL_DISPLACEMENT__
4570 CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", feLhs, null, null);
4571 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", feSpringLhsPtr, null, null);
4572 CHKERR DMMoFEMSNESSetJacobian(dm, "CONTACT", feLhsSimpleContact, null, null);
4573 CHKERR DMMoFEMSNESSetJacobian(dm, "MORTAR_CONTACT", feLhsMortarContact, null,
4574 null);
4575 {
4576 const MoFEM::Problem *problem_ptr;
4577 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
4578 if (problem_ptr->getName() == "EIGEN_ELASTIC") {
4580 null, null);
4581 }
4582 }
4583
4584#ifdef __ANALITICAL_DISPLACEMENT__
4585 if (mField.check_problem("BC_PROBLEM")) {
4588 }
4589#endif //__ANALITICAL_DISPLACEMENT__
4592 CHKERR DMMoFEMSNESSetJacobian(dm, "ARC_LENGTH", arc_method, null, null);
4593
4594 // Rhs
4595 auto fe_set_option = boost::make_shared<FEMethod>();
4596 fe_set_option->preProcessHook = [fe_set_option]() {
4597 return VecSetOption(fe_set_option->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
4598 PETSC_TRUE);
4599 };
4600 CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, fe_set_option, null);
4602 null);
4603 CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", null, zeroFlambda, null);
4604#ifdef __ANALITICAL_DISPLACEMENT__
4605 if (mField.check_problem("BC_PROBLEM")) {
4607 analyticalDirichletBc, null);
4608 }
4609#endif //__ANALITICAL_DISPLACEMENT__
4610 CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", feRhs, null, null);
4611 CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", feSpringRhsPtr, null, null);
4612 CHKERR DMMoFEMSNESSetFunction(dm, "CONTACT", feRhsSimpleContact, null, null);
4613 CHKERR DMMoFEMSNESSetFunction(dm, "MORTAR_CONTACT", feRhsMortarContact, null,
4614 null);
4615 {
4616 const MoFEM::Problem *problem_ptr;
4617 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
4618 if (problem_ptr->getName() == "EIGEN_ELASTIC") {
4620 null, null);
4621 }
4622 }
4623
4624 for (auto fit = surfaceForces->begin(); fit != surfaceForces->end(); fit++) {
4626 dm, fit->first.c_str(),
4627 boost::shared_ptr<FEMethod>(surfaceForces, &fit->second->getLoopFe()),
4628 null, null);
4629 }
4630 for (auto fit = surfacePressure->begin(); fit != surfacePressure->end();
4631 fit++) {
4633 dm, fit->first.c_str(),
4634 boost::shared_ptr<FEMethod>(surfacePressure, &fit->second->getLoopFe()),
4635 null, null);
4636 }
4637 for (auto fit = edgeForces->begin(); fit != edgeForces->end(); fit++) {
4639 dm, fit->first.c_str(),
4640 boost::shared_ptr<FEMethod>(edgeForces, &fit->second->getLoopFe()),
4641 null, null);
4642 }
4643 for (auto fit = nodalForces->begin(); fit != nodalForces->end(); fit++) {
4645 dm, fit->first.c_str(),
4646 boost::shared_ptr<FEMethod>(nodalForces, &fit->second->getLoopFe()),
4647 null, null);
4648 }
4649#ifdef __ANALITICAL_TRACTION__
4651 for (auto fit = analiticalSurfaceElement->begin();
4652 fit != analiticalSurfaceElement->end(); fit++) {
4654 dm, fit->first.c_str(),
4655 boost::shared_ptr<FEMethod>(analiticalSurfaceElement,
4656 &fit->second->fe),
4657 null, null);
4658 }
4659 }
4660#endif //__ANALITICAL_TRACTION__
4661#ifdef __ANALITICAL_DISPLACEMENT__
4662 if (mField.check_problem("BC_PROBLEM")) {
4665 }
4666#endif //__ANALITICAL_DISPLACEMENT__
4667 CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", assembleFlambda, null, null);
4668 CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", arc_method, null, null);
4671
4672 if (debug) {
4673 if (m == PETSC_NULL || q == PETSC_NULL || f == PETSC_NULL) {
4674 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
4675 "problem matrix or vectors are set to null");
4676 }
4677 SNES snes;
4678 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
4679 SnesCtx *snes_ctx;
4680 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
4681 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
4682 CHKERR SNESSetJacobian(snes, m, m, SnesMat, snes_ctx);
4683 CHKERR SnesMat(snes, q, m, m, snes_ctx);
4684 CHKERR SnesRhs(snes, q, f, snes_ctx);
4685 CHKERR SNESDestroy(&snes);
4686 MatView(m, PETSC_VIEWER_DRAW_WORLD);
4687 std::string wait;
4688 std::cin >> wait;
4689 }
4690
4692}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:665
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:706
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1041
virtual bool check_problem(const std::string name)=0
check if problem exist
auto f
Definition: HenckyOps.hpp:5
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:136
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringRhsPtr
boost::shared_ptr< AnalyticalDirichletBC::DirichletBC > analyticalDirichletBc
boost::shared_ptr< FaceElementForcesAndSourcesCore > feSpringLhsPtr
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressure
assemble surface pressure
MoFEMErrorCode getArcLengthDof()
set pointer to arc-length DOF
boost::shared_ptr< DirichletSpatialPositionsBc > spatialDirichletBc
apply Dirichlet BC to sparial positions
boost::shared_ptr< CrackFrontElement > feRhs
Integrate elastic FE.
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContact
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > analiticalSurfaceElement
boost::shared_ptr< FEMethod > zeroFlambda
assemble F_lambda vector
boost::shared_ptr< boost::ptr_map< string, EdgeForce > > edgeForces
assemble edge forces
boost::shared_ptr< DofEntity > arcLengthDof
boost::shared_ptr< boost::ptr_map< string, NodalForce > > nodalForces
assemble nodal forces
boost::shared_ptr< MortarContactProblem::MortarContactElement > feLhsMortarContact
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContact
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfaceForces
assemble surface forces
boost::shared_ptr< MWLSApprox > mwlsApprox
boost::shared_ptr< CrackFrontElement > feLhs
Integrate elastic FE.
boost::shared_ptr< FEMethod > assembleFlambda
assemble F_lambda vector
boost::shared_ptr< BothSurfaceConstrains > closeCrackConstrains
boost::shared_ptr< MortarContactProblem::MortarContactElement > feRhsMortarContact
keeps basic data about problem
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:15
RHS-operator for pressure element (spatial configuration)
Operator to assemble nodal force into right hand side vector.
Definition: NodalForce.hpp:35
Zero F_lambda.

◆ 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 6508 of file CrackPropagation.cpp.

6509 {
6510 boost::shared_ptr<FEMethod> null;
6512
6513 // Set up DM specific vectors and data
6514 Vec front_f, tangent_front_f;
6515 CHKERR DMCreateGlobalVector(dm, &front_f);
6516 CHKERR VecDuplicate(front_f, &tangent_front_f);
6517 CHKERR VecSetOption(front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6518 CHKERR VecSetOption(tangent_front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6519 if (smootherFe->smootherData.ownVectors) {
6520 if (smootherFe->smootherData.frontF != PETSC_NULL)
6521 CHKERR VecDestroy(&smootherFe->smootherData.frontF);
6522 if (smootherFe->smootherData.tangentFrontF != PETSC_NULL)
6523 CHKERR VecDestroy(&smootherFe->smootherData.tangentFrontF);
6524 }
6525 smootherFe->smootherData.frontF = front_f;
6526 smootherFe->smootherData.tangentFrontF = tangent_front_f;
6527 smootherFe->smootherData.ownVectors = true;
6528 CHKERR PetscObjectReference((PetscObject)front_f);
6529 CHKERR PetscObjectReference((PetscObject)tangent_front_f);
6530 if (tangentConstrains->ownVectors) {
6531 if (tangentConstrains->frontF != PETSC_NULL)
6532 VecDestroy(&tangentConstrains->frontF);
6533 if (tangentConstrains->tangentFrontF != PETSC_NULL)
6534 VecDestroy(&tangentConstrains->tangentFrontF);
6535 }
6536 tangentConstrains->frontF = front_f;
6537 tangentConstrains->tangentFrontF = tangent_front_f;
6538 tangentConstrains->ownVectors = true;
6539 CHKERR PetscObjectReference((PetscObject)front_f);
6540 CHKERR PetscObjectReference((PetscObject)tangent_front_f);
6541 CHKERR VecDestroy(&front_f);
6542 CHKERR VecDestroy(&tangent_front_f);
6543
6544 const MoFEM::Problem *problem_ptr;
6545 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
6546 CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
6547 problem_ptr->getName(), COL, feGriffithConstrainsDelta->deltaVec);
6548 CHKERR VecSetOption(feGriffithConstrainsDelta->deltaVec,
6549 VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6552
6553 // Fix nodes
6554 CHKERR updateMaterialFixedNode(true, false, debug);
6555
6556 if (debug && !mField.get_comm_rank()) {
6557 const FieldEntity_multiIndex *field_ents;
6558 CHKERR mField.get_field_ents(&field_ents);
6559 for (Range::iterator vit = crackFrontNodes.begin();
6560 vit != crackFrontNodes.end(); ++vit) {
6561 cerr << "Node " << endl;
6562 FieldEntity_multiIndex::index<Ent_mi_tag>::type::iterator eit, hi_eit;
6563 eit = field_ents->get<Ent_mi_tag>().lower_bound(*vit);
6564 hi_eit = field_ents->get<Ent_mi_tag>().upper_bound(*vit);
6565 for (; eit != hi_eit; ++eit) {
6566 cerr << **eit << endl;
6567 }
6568 }
6569 }
6570
6571 SnesCtx *snes_ctx;
6572 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6573 CHKERR snes_ctx->clearLoops();
6574
6575 // Add to SNES Jacobian
6577 CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feSmootherLhs, null, null);
6578 CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feMaterialLhs, null, null);
6579 CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6580 null, null);
6581 CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CONTACT",
6582 bothSidesContactConstrains, null, null);
6583 for (auto &m : surfaceConstrain) {
6584 if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6585 CHKERR DMMoFEMSNESSetJacobian(dm, "SURFACE", m.second->feLhsPtr, null,
6586 null);
6587 } else {
6588 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACK_SURFACE", m.second->feLhsPtr,
6589 null, null);
6590 }
6591 }
6592 for (auto &m : edgeConstrains) {
6593 CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE", m.second->feLhsPtr, null, null);
6594 }
6595 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
6596 feGriffithConstrainsDelta, null, null);
6597 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
6598 feGriffithConstrainsLhs, null, null);
6599 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
6600 tangentConstrains, null, null);
6601 CHKERR DMMoFEMSNESSetJacobian(dm, "GRIFFITH_FORCE_ELEMENT",
6602 feGriffithForceLhs, null, null);
6604
6605 // Add to SNES residuals
6607 CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feSmootherRhs, null, null);
6608 CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feMaterialRhs, null, null);
6609 CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6610 null, null);
6611 CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CONTACT",
6612 bothSidesContactConstrains, null, null);
6613 for (auto &m : surfaceConstrain) {
6614 if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6615 CHKERR DMMoFEMSNESSetFunction(dm, "SURFACE", m.second->feRhsPtr, null,
6616 null);
6617 } else {
6618 CHKERR DMMoFEMSNESSetFunction(dm, "CRACK_SURFACE", m.second->feRhsPtr,
6619 null, null);
6620 }
6621 }
6622 for (auto &m : edgeConstrains) {
6623 CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", m.second->feRhsPtr, null, null);
6624 }
6625 CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
6626 feGriffithConstrainsDelta, null, null);
6627 CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
6628 feGriffithConstrainsRhs, null, null);
6629 CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
6630 tangentConstrains, null, null);
6631 CHKERR DMMoFEMSNESSetFunction(dm, "GRIFFITH_FORCE_ELEMENT",
6632 feGriffithForceRhs, null, null);
6634
6635 if (debug) {
6636 Vec q, f;
6637 CHKERR DMCreateGlobalVector(dm, &q);
6638 CHKERR VecDuplicate(q, &f);
6639 CHKERR DMoFEMMeshToLocalVector(dm, q, INSERT_VALUES, SCATTER_FORWARD);
6640 CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6641 Mat m;
6642 CHKERR DMCreateMatrix(dm, &m);
6643 SNES snes;
6644 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
6645 SnesCtx *snes_ctx;
6646 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6647 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
6648 CHKERR SNESSetJacobian(snes, m, m, SnesMat, snes_ctx);
6649 CHKERR SnesMat(snes, q, m, m, snes_ctx);
6650 CHKERR SnesRhs(snes, q, f, snes_ctx);
6651 CHKERR SNESDestroy(&snes);
6652 // MatView(m,PETSC_VIEWER_DRAW_WORLD);
6653 MatView(m, PETSC_VIEWER_STDOUT_WORLD);
6654 std::string wait;
6655 std::cin >> wait;
6656 CHKERR VecDestroy(&q);
6657 CHKERR VecDestroy(&f);
6658 CHKERR MatDestroy(&m);
6659 }
6661}
@ COL
Definition: definitions.h:123
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< FieldEntity, UId, &FieldEntity::localUId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FieldEntity::interface_type_RefEntity, EntityHandle, &FieldEntity::getEnt > > > > FieldEntity_multiIndex
MultiIndex container keeps FieldEntity.
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
const FTensor::Tensor2< T, Dim, Dim > Vec
boost::shared_ptr< BothSurfaceConstrains > bothSidesConstrains
boost::shared_ptr< CrackFrontElement > feMaterialRhs
Integrate material stresses, assemble vector.
map< int, boost::shared_ptr< SurfaceSlidingConstrains > > surfaceConstrain
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsLhs
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherLhs
Integrate smoothing operators.
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrainsDelta > feGriffithConstrainsDelta
boost::shared_ptr< Smoother::MyVolumeFE > feSmootherRhs
Integrate smoothing operators.
boost::shared_ptr< DirichletFixFieldAtEntitiesBc > fixMaterialEnts
boost::shared_ptr< GriffithForceElement::MyTriangleFEConstrains > feGriffithConstrainsRhs
boost::shared_ptr< Smoother > smootherFe
MoFEMErrorCode updateMaterialFixedNode(const bool fix_front, const bool fix_small_g, const bool debug=false)
Update fixed nodes.
map< int, boost::shared_ptr< EdgeSlidingConstrains > > edgeConstrains
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceLhs
boost::shared_ptr< BothSurfaceConstrains > bothSidesContactConstrains
boost::shared_ptr< ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain > tangentConstrains
Constrains crack front in tangent directiona.
boost::shared_ptr< GriffithForceElement::MyTriangleFE > feGriffithForceRhs
boost::shared_ptr< CrackFrontElement > feMaterialLhs
Integrate material stresses, assemble matrix.
virtual int get_comm_rank() const =0
MoFEMErrorCode clearLoops()
Clear loops.
Definition: SnesCtx.cpp:16
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23

◆ addMWLSDensityOperators()

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

Definition at line 4904 of file CrackPropagation.cpp.

4906 {
4908
4909 if (densityMapBlock == -1)
4911
4912 if (!elasticFe) {
4913 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
4914 "Elastic element instance not created");
4915 }
4916
4917 // Create elastic element finite element instance for residuals. Note that
4918 // this element approx. singularity at crack front.
4919 fe_rhs = boost::make_shared<CrackFrontElement>(
4922 // Create finite element instance for assembly of tangent matrix
4923 fe_lhs = boost::make_shared<CrackFrontElement>(
4926
4927 bool test_with_mwls = false; // switch for debugging
4928
4929 // Arbitrary lagrangian formulation, mesh node positions are taken into
4930 // account by operators.
4931 fe_rhs->meshPositionsFieldName = "NONE";
4932 fe_lhs->meshPositionsFieldName = "NONE";
4933 fe_rhs->addToRule = 0;
4934 fe_lhs->addToRule = 0;
4935
4936 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
4937 data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
4938 boost::shared_ptr<map<int, NonlinearElasticElement::BlockData>>
4939 block_sets_ptr(elasticFe, &elasticFe->setOfBlocks);
4940
4941 if (!mwlsApprox) {
4942 // Create instance for moving least square approximation
4943 mwlsApprox = boost::make_shared<MWLSApprox>(mField);
4944 // Load mesh with stresses
4945 CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
4946 MeshsetsManager *block_manager_ptr;
4947 CHKERR mField.getInterface(block_manager_ptr);
4948 CHKERR block_manager_ptr->getEntitiesByDimension(
4949 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
4950 Tag th_rho;
4951 if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(), th_rho) !=
4952 MB_SUCCESS) {
4953 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4954 "Density tag name %s cannot be found. Please "
4955 "provide the correct name.",
4956 mwlsRhoTagName.c_str());
4957 }
4958 CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(), th_rho);
4959 CHKERR mwlsApprox->getValuesToNodes(th_rho);
4960 } else {
4961 mwlsApprox->tetsInBlock.clear();
4962 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
4963 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
4964 }
4965
4966 auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4967 fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4968 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4969 fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4970 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4971
4972 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_MWLS_ptr(
4973 mwlsApprox, &mwlsApprox->mwlsMaterialPositions);
4974 fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4975 "MESH_NODE_POSITIONS", mat_pos_at_pts_MWLS_ptr));
4976 fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4977 "MESH_NODE_POSITIONS", mat_pos_at_pts_MWLS_ptr));
4978 auto mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4979 fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4980 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4981 fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4982 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4983 fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4984 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
4985 fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4986 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
4987 fe_rhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4988 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4989 fe_rhs->singularElement, fe_rhs->invSJac));
4990 fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4991 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4992 fe_lhs->singularElement, fe_lhs->invSJac));
4993 auto space_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4994 fe_rhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4995 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_rhs->singularElement,
4996 fe_rhs->invSJac));
4997 fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4998 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
4999 fe_lhs->invSJac));
5000
5001 fe_rhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
5002 fe_rhs->singularElement, fe_rhs->detS, fe_rhs->invSJac));
5003
5004 // data_hooke_element_at_pts->HMat = mat_grad_pos_at_pts_ptr;
5005 // data_hooke_element_at_pts->hMat = space_grad_pos_at_pts_ptr;
5006
5007 // RHS //
5008 if (test_with_mwls) {
5009 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5010 fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
5011 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
5012 mwlsRhoTagName, true, false));
5013 else {
5014 fe_rhs->getOpPtrVector().push_back(
5016 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox));
5017 fe_rhs->getOpPtrVector().push_back(
5018 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5019 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
5020 mwlsRhoTagName, true, false));
5021 }
5022 } else {
5023
5024 fe_rhs->getOpPtrVector().push_back(new OpGetDensityFieldForTesting(
5025 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr, mwlsApprox->rhoAtGaussPts,
5026 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5027 fe_rhs, mwlsApprox->singularInitialDisplacement,
5028 mat_grad_pos_at_pts_ptr));
5029 }
5030
5031 fe_rhs->getOpPtrVector().push_back(
5032 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5033 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5034 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone, rHo0));
5035 fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStrainAle(
5036 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5037 fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStress<1>(
5038 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5039 fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEnergy(
5040 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5041
5042 fe_rhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEshelbyStress(
5043 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5044
5045 fe_rhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dX(
5046 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5047 fe_rhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dx(
5048 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
5049
5050 fe_rhs->getOpPtrVector().push_back(new OpRhsBoneExplicitDerivariveWithHooke(
5051 *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5052 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->tetsInBlock, rHo0, nBone,
5053 &crackFrontNodes));
5054
5055 // LHS //
5056 fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
5057 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
5058 fe_lhs->invSJac));
5059 fe_lhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
5060 fe_lhs->singularElement, fe_lhs->detS, fe_lhs->invSJac));
5061
5062 if (test_with_mwls) {
5063 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5064 fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
5065 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
5066 mwlsRhoTagName, true, true));
5067 else {
5068 fe_lhs->getOpPtrVector().push_back(
5070 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox));
5071 fe_lhs->getOpPtrVector().push_back(
5072 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5073 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
5074 mwlsRhoTagName, true, true));
5075 }
5076 } else {
5077 fe_lhs->getOpPtrVector().push_back(new OpGetDensityFieldForTesting(
5078 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr, mwlsApprox->rhoAtGaussPts,
5079 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5080 fe_lhs, mwlsApprox->singularInitialDisplacement,
5081 mat_grad_pos_at_pts_ptr));
5082 }
5083
5084 fe_lhs->getOpPtrVector().push_back(
5085 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5086 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5087 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone, rHo0));
5088
5089 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStrainAle(
5090 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5091 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateStress<1>(
5092 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5093 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEnergy(
5094 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5095 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpCalculateEshelbyStress(
5096 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5097
5098 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dX_dX<1>(
5099 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5100 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhsPre_dX_dx<1>(
5101 "MESH_NODE_POSITIONS", "SPATIAL_POSITION", data_hooke_element_at_pts));
5102 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dX_dx(
5103 "MESH_NODE_POSITIONS", "SPATIAL_POSITION", data_hooke_element_at_pts));
5104 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<1>(
5105 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
5106 fe_lhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dX<1>(
5107 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5108
5109 fe_lhs->getOpPtrVector().push_back(
5110 new HookeElement::OpAleLhsWithDensity_dX_dX(
5111 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5112 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5113 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
5114
5115 fe_lhs->getOpPtrVector().push_back(
5116 new HookeElement::OpAleLhsWithDensity_dx_dX(
5117 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", data_hooke_element_at_pts,
5118 mwlsApprox->rhoAtGaussPts, mwlsApprox->diffRhoAtGaussPts, nBone,
5119 rHo0));
5120
5121 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
5122 if (addSingularity) {
5123 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
5124 mwlsApprox, &mwlsApprox->singularInitialDisplacement);
5125
5126 fe_lhs->getOpPtrVector().push_back(
5127 new OpAleLhsWithDensitySingularElement_dX_dX(
5128 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5129 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5130 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0, mat_singular_disp_ptr));
5131
5132 fe_lhs->getOpPtrVector().push_back(
5133 new OpAleLhsWithDensitySingularElement_dx_dX(
5134 "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
5135 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5136 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0, mat_singular_disp_ptr));
5137 }
5138
5139 fe_lhs->getOpPtrVector().push_back(
5140 new OpLhsBoneExplicitDerivariveWithHooke_dX(
5141 *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5142 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5143 mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
5145
5146 fe_lhs->getOpPtrVector().push_back(
5147 new OpLhsBoneExplicitDerivariveWithHooke_dx(
5148 *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5149 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5150 mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
5152
5154}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ BLOCKSET
Definition: definitions.h:148
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimension
std::string mwlsApproxFile
Name of file with internal stresses.
boost::shared_ptr< NonlinearElasticElement > elasticFe
std::string mwlsRhoTagName
Name of tag with density.
OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl< VolumeElementForcesAndSourcesCore > OpMWLSCalculateBaseCoeffcientsAtGaussPts
Definition: MWLS.hpp:296
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.

◆ addMWLSStressOperators()

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

Definition at line 4779 of file CrackPropagation.cpp.

4781 {
4783
4784 if (residualStressBlock == -1)
4786
4787 if (!elasticFe) {
4788 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
4789 "Elastic element instance not created");
4790 }
4791
4792 // Create elastic element finite element instance for residuals. Note that
4793 // this element approx. singularity at crack front.
4794 fe_rhs = boost::make_shared<CrackFrontElement>(
4797 // Create finite element instance for assembly of tangent matrix
4798 fe_lhs = boost::make_shared<CrackFrontElement>(
4801
4802 // Arbitrary lagrangian formulation, mesh node positions are taken into
4803 // account by operators.
4804 fe_rhs->meshPositionsFieldName = "NONE";
4805 fe_lhs->meshPositionsFieldName = "NONE";
4806 fe_rhs->addToRule = 0;
4807 fe_lhs->addToRule = 0;
4808
4809 if (!mwlsApprox) {
4810 // Create instance for moving least square approximation
4811 mwlsApprox = boost::make_shared<MWLSApprox>(mField);
4812 // Load mesh with stresses
4813 CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
4814 MeshsetsManager *block_manager_ptr;
4815 CHKERR mField.getInterface(block_manager_ptr);
4816 CHKERR block_manager_ptr->getEntitiesByDimension(
4817 residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
4818 Tag th;
4819 if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsStressTagName.c_str(), th) !=
4820 MB_SUCCESS) {
4821 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4822 "Internal stress tag name %s cannot be found. Please "
4823 "provide the correct name.",
4824 mwlsStressTagName.substr(4).c_str());
4825 }
4826 CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsStressTagName.c_str(), th);
4827 CHKERR mwlsApprox->getValuesToNodes(th);
4828 } else {
4829 mwlsApprox->tetsInBlock.clear();
4830 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
4831 residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
4832 }
4833
4834 auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4835 fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4836 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4837 fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4838 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4839 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr(new MatrixDouble());
4840 fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4841 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4842 fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
4843 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4844
4845 boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr(new MatrixDouble());
4846 fe_rhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4847 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_rhs->singularElement,
4848 fe_rhs->invSJac));
4849 fe_rhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4850 fe_rhs->singularElement, fe_rhs->detS, fe_rhs->invSJac));
4851
4852 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4853 fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSStressAtGaussPts(
4854 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
4855 mwlsStressTagName, false, false));
4856 else {
4857 fe_rhs->getOpPtrVector().push_back(
4859 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox));
4860 fe_rhs->getOpPtrVector().push_back(
4861 new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
4862 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs, mwlsApprox,
4863 mwlsStressTagName, false, false));
4864 }
4865
4866 fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressRhs(
4867 mat_grad_pos_at_pts_ptr, mwlsApprox, false));
4868 fe_rhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSMaterialStressRhs(
4869 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
4870 &crackFrontNodes));
4871
4872 fe_lhs->getOpPtrVector().push_back(new OpGetCrackFrontDataGradientAtGaussPts(
4873 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
4874 fe_lhs->invSJac));
4875 fe_lhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4876 fe_lhs->singularElement, fe_lhs->detS, fe_lhs->invSJac));
4877
4878 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4879 fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSStressAtGaussPts(
4880 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
4881 mwlsStressTagName, true, false));
4882 else {
4883 fe_lhs->getOpPtrVector().push_back(
4885 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox));
4886 fe_lhs->getOpPtrVector().push_back(
4887 new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
4888 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs, mwlsApprox,
4889 mwlsStressTagName, true, false));
4890 }
4891
4892 fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressLhs_DX(
4893 mat_grad_pos_at_pts_ptr, mwlsApprox));
4894 fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSMaterialStressLhs_DX(
4895 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
4896 &crackFrontNodes));
4897 fe_lhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSMaterialStressLhs_Dx(
4898 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
4899 &crackFrontNodes));
4900
4902}
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::string mwlsStressTagName
Name of tag with internal stresses.

◆ 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 6798 of file CrackPropagation.cpp.

6801 {
6802 boost::shared_ptr<FEMethod> null;
6804
6805 // Assemble F_lambda force
6806 assembleFlambda = boost::shared_ptr<FEMethod>(
6807 new AssembleFlambda(arc_ctx, spatialDirichletBc));
6808 zeroFlambda = boost::shared_ptr<FEMethod>(new ZeroFLmabda(arc_ctx));
6809
6810 if (mwlsApprox) {
6811 mwlsApprox->F_lambda = arc_ctx->F_lambda;
6813 mwlsApprox->arcLengthDof = arcLengthDof;
6814 }
6815
6816 // Set up DM specific vectors and data
6817 spatialDirichletBc->mapZeroRows.clear();
6818
6819 // Set-up F_lambda to elements & operator evalulating forces
6820 // Surface fores
6821 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6822 surfaceForces->begin();
6823 fit != surfaceForces->end(); fit++) {
6824 boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator>::iterator oit,
6825 hi_oit;
6826 oit = fit->second->getLoopFe().getOpPtrVector().begin();
6827 hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6828 for (; oit != hi_oit; oit++) {
6829 if (boost::typeindex::type_id_runtime(*oit) ==
6830 boost::typeindex::type_id<NeumannForcesSurface::OpNeumannForce>()) {
6831 dynamic_cast<NeumannForcesSurface::OpNeumannForce &>(*oit).F =
6832 arc_ctx->F_lambda;
6833 }
6834 }
6835 }
6836 // Surface pressure
6837 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6838 surfacePressure->begin();
6839 fit != surfacePressure->end(); fit++) {
6840
6841 boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator>::iterator oit,
6842 hi_oit;
6843 oit = fit->second->getLoopFe().getOpPtrVector().begin();
6844 hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6845 for (; oit != hi_oit; oit++) {
6846 if (boost::typeindex::type_id_runtime(*oit) ==
6847 boost::typeindex::type_id<
6849 dynamic_cast<NeumannForcesSurface::OpNeumannPressure &>(*oit).F =
6850 arc_ctx->F_lambda;
6851 }
6852 if (boost::typeindex::type_id_runtime(*oit) ==
6853 boost::typeindex::type_id<
6856 arc_ctx->F_lambda;
6857 }
6858 }
6859 }
6860 // Surface pressure (ALE)
6861 if (isPressureAle) {
6862 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6863 surfacePressureAle->begin();
6864 fit != surfacePressureAle->end(); fit++) {
6865
6868
6869 boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator>::iterator oit,
6870 hi_oit;
6871 oit = fit->second->getLoopFeMatRhs().getOpPtrVector().begin();
6872 hi_oit = fit->second->getLoopFeMatRhs().getOpPtrVector().end();
6873 for (; oit != hi_oit; oit++) {
6874 if (boost::typeindex::type_id_runtime(*oit) ==
6875 boost::typeindex::type_id<
6878 *oit)
6879 .F = arc_ctx->F_lambda;
6880 }
6881 }
6882 }
6883 }
6884 // Surface edge
6885 for (boost::ptr_map<string, EdgeForce>::iterator fit = edgeForces->begin();
6886 fit != edgeForces->end(); fit++) {
6887 boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator>::iterator oit,
6888 hi_oit;
6889 oit = fit->second->getLoopFe().getOpPtrVector().begin();
6890 hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6891 for (; oit != hi_oit; oit++) {
6892 if (boost::typeindex::type_id_runtime(*oit) ==
6893 boost::typeindex::type_id<EdgeForce::OpEdgeForce>()) {
6894 dynamic_cast<EdgeForce::OpEdgeForce &>(*oit).F = arc_ctx->F_lambda;
6895 }
6896 }
6897 }
6898 // Nodal force
6899 for (boost::ptr_map<string, NodalForce>::iterator fit = nodalForces->begin();
6900 fit != nodalForces->end(); fit++) {
6901 boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator>::iterator oit,
6902 hi_oit;
6903 oit = fit->second->getLoopFe().getOpPtrVector().begin();
6904 hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6905 for (; oit != hi_oit; oit++) {
6906 if (boost::typeindex::type_id_runtime(*oit) ==
6907 boost::typeindex::type_id<NodalForce::OpNodalForce>()) {
6908 dynamic_cast<NodalForce::OpNodalForce &>(*oit).F = arc_ctx->F_lambda;
6909 }
6910 }
6911 }
6912 // Set up DM specific vectors and data
6913 Vec front_f, tangent_front_f;
6914 CHKERR DMCreateGlobalVector(dm, &front_f);
6915 CHKERR VecDuplicate(front_f, &tangent_front_f);
6916 CHKERR VecSetOption(front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6917 CHKERR VecSetOption(tangent_front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6918
6919 {
6920 if (smootherFe->smootherData.ownVectors) {
6921 if (smootherFe->smootherData.frontF != PETSC_NULL)
6922 CHKERR VecDestroy(&smootherFe->smootherData.frontF);
6923 if (smootherFe->smootherData.tangentFrontF != PETSC_NULL)
6924 CHKERR VecDestroy(&smootherFe->smootherData.tangentFrontF);
6925 }
6926 smootherFe->smootherData.frontF = front_f;
6927 smootherFe->smootherData.tangentFrontF = tangent_front_f;
6928 smootherFe->smootherData.ownVectors = true;
6929
6930 if (tangentConstrains->ownVectors) {
6931 if (tangentConstrains->frontF != PETSC_NULL)
6932 CHKERR VecDestroy(&tangentConstrains->frontF);
6933 if (tangentConstrains->tangentFrontF != PETSC_NULL)
6934 CHKERR VecDestroy(&tangentConstrains->tangentFrontF);
6935 }
6936 tangentConstrains->frontF = front_f;
6937 tangentConstrains->tangentFrontF = tangent_front_f;
6938 tangentConstrains->ownVectors = false;
6939 }
6940
6941 const MoFEM::Problem *problem_ptr;
6942 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
6943 CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
6944 problem_ptr->getName(), COL, feGriffithConstrainsDelta->deltaVec);
6945 CHKERR VecSetOption(feGriffithConstrainsDelta->deltaVec,
6946 VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6949
6950 // Fix nodes
6951 CHKERR updateMaterialFixedNode(false, true, debug);
6952 dynamic_cast<AssembleFlambda *>(assembleFlambda.get())
6953 ->pushDirichletBC(fixMaterialEnts);
6954
6955 // Add to SNES Jacobian
6956 // Lhs
6959 null);
6960 CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC_NOT_ALE", feLhs, null, null);
6962 null);
6964 null);
6965 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", feSpringLhsPtr, null, null);
6966 CHKERR DMMoFEMSNESSetJacobian(dm, "CONTACT", feLhsSimpleContact, null, null);
6967 CHKERR DMMoFEMSNESSetJacobian(dm, "MORTAR_CONTACT", feLhsMortarContact, null,
6968 null);
6969 CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6970 null, null);
6971 CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CONTACT",
6972 bothSidesContactConstrains, null, null);
6973 for (auto &m : surfaceConstrain) {
6974 if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6975 CHKERR DMMoFEMSNESSetJacobian(dm, "SURFACE", m.second->feLhsPtr, null,
6976 null);
6977 } else {
6978 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACK_SURFACE", m.second->feLhsPtr,
6979 null, null);
6980 }
6981 }
6982 for (auto &m : edgeConstrains)
6983 CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE", m.second->feLhsPtr, null, null);
6984
6985 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
6986 feGriffithConstrainsDelta, null, null);
6987 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_ELEM",
6988 feGriffithConstrainsLhs, null, null);
6989 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
6990 tangentConstrains, null, null);
6991 CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feSmootherLhs, null, null);
6992 CHKERR DMMoFEMSNESSetJacobian(dm, "GRIFFITH_FORCE_ELEMENT",
6993 feGriffithForceLhs, null, null);
6997 CHKERR DMMoFEMSNESSetJacobian(dm, "ARC_LENGTH", arc_method, null, null);
6998
6999 // Rhs
7002 null);
7003 CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", null, zeroFlambda, null);
7004 CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", feRhs, null, null);
7005 CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", feSpringRhsPtr, null, null);
7006 CHKERR DMMoFEMSNESSetFunction(dm, "CONTACT", feRhsSimpleContact, null, null);
7007 CHKERR DMMoFEMSNESSetFunction(dm, "MORTAR_CONTACT", feRhsMortarContact, null,
7008 null);
7009 CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
7010 null, null);
7011 CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CONTACT",
7012 bothSidesContactConstrains, null, null);
7013 for (auto &m : surfaceConstrain) {
7014 if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
7015 CHKERR DMMoFEMSNESSetFunction(dm, "SURFACE", m.second->feRhsPtr, null,
7016 null);
7017 } else {
7018 CHKERR DMMoFEMSNESSetFunction(dm, "CRACK_SURFACE", m.second->feRhsPtr,
7019 null, null);
7020 }
7021 }
7022 for (auto &m : edgeConstrains) {
7023 CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", m.second->feRhsPtr, null, null);
7024 }
7025 CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feSmootherRhs, null, null);
7026 CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
7027 feGriffithConstrainsDelta, null, null);
7028 CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_ELEM",
7029 feGriffithConstrainsRhs, null, null);
7030 CHKERR DMMoFEMSNESSetFunction(dm, "CRACKFRONT_AREA_TANGENT_ELEM",
7031 tangentConstrains, null, null);
7032 CHKERR DMMoFEMSNESSetFunction(dm, "MATERIAL", feMaterialRhs, null, null);
7033 CHKERR DMMoFEMSNESSetFunction(dm, "GRIFFITH_FORCE_ELEMENT",
7034 feGriffithForceRhs, null, null);
7035 // Add force elements
7036 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7037 surfaceForces->begin();
7038 fit != surfaceForces->end(); fit++) {
7040 dm, fit->first.c_str(),
7041 boost::shared_ptr<FEMethod>(surfaceForces, &fit->second->getLoopFe()),
7042 null, null);
7043 }
7044 // Add pressure elements
7045 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7046 surfacePressure->begin();
7047 fit != surfacePressure->end(); fit++) {
7049 dm, fit->first.c_str(),
7050 boost::shared_ptr<FEMethod>(surfacePressure, &fit->second->getLoopFe()),
7051 null, null);
7052 }
7053 // Add pressure elements (ALE)
7054 if (isPressureAle) {
7055 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7056 surfacePressureAle->begin();
7057 fit != surfacePressureAle->end(); fit++) {
7059 dm, fit->first.c_str(),
7060 boost::shared_ptr<FEMethod>(surfacePressureAle,
7061 &fit->second->getLoopFeLhs()),
7062 null, null);
7064 dm, fit->first.c_str(),
7065 boost::shared_ptr<FEMethod>(surfacePressureAle,
7066 &fit->second->getLoopFeMatRhs()),
7067 null, null);
7069 dm, fit->first.c_str(),
7070 boost::shared_ptr<FEMethod>(surfacePressureAle,
7071 &fit->second->getLoopFeMatLhs()),
7072 null, null);
7073 }
7074 }
7075
7076 if (areSpringsAle) {
7078 dm, "SPRING_ALE", feRhsSpringALEMaterial.get(), PETSC_NULL, PETSC_NULL);
7079 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING_ALE", feLhsSpringALE.get(), NULL,
7080 NULL);
7081 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING_ALE",
7082 feLhsSpringALEMaterial.get(), NULL, NULL);
7083 }
7084
7085 if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
7086 CHKERR DMMoFEMSNESSetFunction(dm, "SIMPLE_CONTACT_ALE",
7088 PETSC_NULL, PETSC_NULL);
7089 CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ALE",
7091 NULL);
7092 CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ALE",
7093 feLhsSimpleContactALE.get(), NULL, NULL);
7094 }
7095 for (boost::ptr_map<string, EdgeForce>::iterator fit = edgeForces->begin();
7096 fit != edgeForces->end(); fit++) {
7098 dm, fit->first.c_str(),
7099 boost::shared_ptr<FEMethod>(edgeForces, &fit->second->getLoopFe()),
7100 null, null);
7101 }
7102 for (boost::ptr_map<string, NodalForce>::iterator fit = nodalForces->begin();
7103 fit != nodalForces->end(); fit++) {
7105 dm, fit->first.c_str(),
7106 boost::shared_ptr<FEMethod>(nodalForces, &fit->second->getLoopFe()),
7107 null, null);
7108 }
7109 CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", assembleFlambda, null, null);
7110 CHKERR DMMoFEMSNESSetFunction(dm, "ARC_LENGTH", arc_method, null, null);
7114
7116}
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContactALE
boost::shared_ptr< FaceElementForcesAndSourcesCore > feLhsSpringALEMaterial
boost::shared_ptr< FaceElementForcesAndSourcesCore > feRhsSpringALEMaterial
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressureAle
assemble surface pressure (ALE)
PetscBool isPressureAle
If true surface pressure is considered in ALE.
boost::shared_ptr< CrackFrontElement > feCouplingMaterialLhs
FE instance to assemble coupling terms.
boost::shared_ptr< CrackFrontElement > feCouplingElasticLhs
FE instance to assemble coupling terms.
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feLhsSimpleContactALEMaterial
boost::shared_ptr< FaceElementForcesAndSourcesCore > feLhsSpringALE
PetscBool areSpringsAle
If true surface spring is considered in ALE.
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > feRhsSimpleContactALEMaterial
boost::shared_ptr< NeumannForcesSurface::DataAtIntegrationPts > commonDataSurfacePressureAle
common data at integration points (ALE)
RHS-operator for the pressure element (material configuration)

◆ 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 6663 of file CrackPropagation.cpp.

6664 {
6665 boost::shared_ptr<FEMethod> null;
6667
6668 // Set up DM specific vectors and data
6669 if (smootherFe->smootherData.ownVectors) {
6670 if (smootherFe->smootherData.frontF != PETSC_NULL)
6671 CHKERR VecDestroy(&smootherFe->smootherData.frontF);
6672 if (smootherFe->smootherData.tangentFrontF != PETSC_NULL)
6673 CHKERR VecDestroy(&smootherFe->smootherData.tangentFrontF);
6674 }
6675 smootherFe->smootherData.frontF = PETSC_NULL;
6676 smootherFe->smootherData.tangentFrontF = PETSC_NULL;
6677 smootherFe->smootherData.ownVectors = false;
6678
6679 // Fix nodes
6680 CHKERR updateMaterialFixedNode(true, false, debug);
6681
6682 if (debug && !mField.get_comm_rank()) {
6683 const FieldEntity_multiIndex *field_ents;
6684 CHKERR mField.get_field_ents(&field_ents);
6685 for (Range::iterator vit = crackFrontNodes.begin();
6686 vit != crackFrontNodes.end(); ++vit) {
6687 cerr << "Node " << endl;
6688 FieldEntity_multiIndex::index<Ent_mi_tag>::type::iterator eit, hi_eit;
6689 eit = field_ents->get<Ent_mi_tag>().lower_bound(*vit);
6690 hi_eit = field_ents->get<Ent_mi_tag>().upper_bound(*vit);
6691 for (; eit != hi_eit; ++eit) {
6692 cerr << **eit << endl;
6693 }
6694 }
6695 }
6696
6697 SnesCtx *snes_ctx;
6698 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6699 CHKERR snes_ctx->clearLoops();
6700
6701 // Add to SNES Jacobian
6703 CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6704 null, null);
6705 CHKERR DMMoFEMSNESSetJacobian(dm, "BOTH_SIDES_OF_CONTACT",
6706 bothSidesContactConstrains, null, null);
6707 for (auto &m : surfaceConstrain) {
6708 if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6709 CHKERR DMMoFEMSNESSetJacobian(dm, "SURFACE", m.second->feLhsPtr, null,
6710 null);
6711 } else {
6712 CHKERR DMMoFEMSNESSetJacobian(dm, "CRACK_SURFACE", m.second->feLhsPtr,
6713 null, null);
6714 }
6715 }
6716 for (auto &m : edgeConstrains) {
6717 CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE", m.second->feLhsPtr, null, null);
6718 }
6719 CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feSmootherLhs, null, null);
6721
6722 // Add to SNES residuals
6724 CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CRACK", bothSidesConstrains,
6725 null, null);
6726 CHKERR DMMoFEMSNESSetFunction(dm, "BOTH_SIDES_OF_CONTACT",
6727 bothSidesContactConstrains, null, null);
6728 for (auto &m : surfaceConstrain) {
6729 if (m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6730 CHKERR DMMoFEMSNESSetFunction(dm, "SURFACE", m.second->feRhsPtr, null,
6731 null);
6732 } else {
6733 CHKERR DMMoFEMSNESSetFunction(dm, "CRACK_SURFACE", m.second->feRhsPtr,
6734 null, null);
6735 }
6736 }
6737 for (auto &m : edgeConstrains) {
6738 CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", m.second->feRhsPtr, null, null);
6739 }
6740
6741 CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feSmootherRhs, null, null);
6743 for (auto &m : edgeConstrains) {
6744 CHKERR DMMoFEMSNESSetFunction(dm, "EDGE", null, null, m.second->feRhsPtr);
6745 }
6746
6747 if (debug) {
6748 Vec q, f;
6749 CHKERR DMCreateGlobalVector(dm, &q);
6750 CHKERR VecDuplicate(q, &f);
6751 CHKERR DMoFEMMeshToLocalVector(dm, q, INSERT_VALUES, SCATTER_FORWARD);
6752 CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6753 Mat m;
6754 CHKERR DMCreateMatrix(dm, &m);
6755 SNES snes;
6756 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
6757 SnesCtx *snes_ctx;
6758 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
6759 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
6760 CHKERR SNESSetJacobian(snes, m, m, SnesMat, snes_ctx);
6761 CHKERR SnesMat(snes, q, m, m, snes_ctx);
6762 CHKERR SnesRhs(snes, q, f, snes_ctx);
6763 CHKERR SNESDestroy(&snes);
6764 // MatView(m,PETSC_VIEWER_DRAW_WORLD);
6765 MatView(m, PETSC_VIEWER_STDOUT_WORLD);
6766 std::string wait;
6767 std::cin >> wait;
6768 CHKERR VecDestroy(&q);
6769 CHKERR VecDestroy(&f);
6770 CHKERR MatDestroy(&m);
6771 }
6773}

◆ analyticalStrainFunction()

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

Definition at line 9857 of file CrackPropagation.cpp.

9861 {
9863 constexpr double alpha = 1.e-5;
9867 // FIXME put here formula from test
9868 double temp = 250.;
9869 double z = t_coords(2);
9870 if ((-10. < z && z < -1.) || std::abs(z + 1.) < 1e-15) {
9871 temp = 10. / 3. * (35. - 4. * z);
9872 }
9873 if ((-1. < z && z < 2.) || std::abs(z - 2.) < 1e-15) {
9874 temp = 10. / 3. * (34. - 5. * z);
9875 }
9876 if (2. < z && z < 10.) {
9877 temp = 5. / 4. * (30. + 17. * z);
9878 }
9879
9880 t_thermal_strain(i, j) = alpha * (temp - 250.) * t_kd(i, j);
9881 return t_thermal_strain;
9882}
Kronecker Delta class symmetric.
constexpr auto t_kd
void temp(int x, int y=10)
Definition: simple.cpp:4
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j

◆ assembleCouplingForcesDM()

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

assemble coupling element instances

Definition at line 6055 of file CrackPropagation.cpp.

6056 {
6058
6059 if (!elasticFe) {
6060 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
6061 "Elastic element instance not created");
6062 }
6063 feCouplingElasticLhs = boost::make_shared<CrackFrontElement>(
6066 feCouplingElasticLhs->meshPositionsFieldName = "NONE";
6067 feCouplingElasticLhs->addToRule = 0;
6068
6070 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
6071 elasticFe, &elasticFe->setOfBlocks);
6072
6073 // Calculate spatial positions and gradients (of deformation) at integration
6074 // pts.
6075 // This operator takes into account singularity at crack front
6076 if (!onlyHooke) {
6077 feCouplingElasticLhs->getOpPtrVector().push_back(
6078 new OpGetCrackFrontCommonDataAtGaussPts(
6079 "SPATIAL_POSITION", elasticFe->commonData,
6080 feCouplingElasticLhs->singularElement,
6081 feCouplingElasticLhs->invSJac));
6082 // Calculate material positions and gradients at integration pts.
6083 feCouplingElasticLhs->getOpPtrVector().push_back(
6085 "MESH_NODE_POSITIONS", elasticFe->commonData));
6086 // Transform base functions to get singular base functions at crack front
6087 feCouplingElasticLhs->getOpPtrVector().push_back(
6088 new OpTransfromSingularBaseFunctions(
6089 feCouplingElasticLhs->singularElement, feCouplingElasticLhs->detS,
6090 feCouplingElasticLhs->invSJac));
6091 // Iterate over material blocks
6092 map<int, NonlinearElasticElement::BlockData>::iterator sit =
6093 elasticFe->setOfBlocks.begin();
6094 for (; sit != elasticFe->setOfBlocks.end(); sit++) {
6095 // Evaluate stress at integration pts
6096 feCouplingElasticLhs->getOpPtrVector().push_back(
6098 "SPATIAL_POSITION", sit->second, elasticFe->commonData,
6099 ELASTIC_TAG, true, true, false));
6100 // Assemble tangent matrix, derivative of spatial positions
6101 feCouplingElasticLhs->getOpPtrVector().push_back(
6103 "SPATIAL_POSITION", "SPATIAL_POSITION", sit->second,
6104 elasticFe->commonData));
6105 // Assemble tangent matrix, derivative of material positions
6106 feCouplingElasticLhs->getOpPtrVector().push_back(
6108 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", sit->second,
6109 elasticFe->commonData));
6110 }
6111 } else {
6112 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
6113 data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
6114 feCouplingElasticLhs->getOpPtrVector().push_back(
6116 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
6117
6118 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(new MatrixDouble());
6119 feCouplingElasticLhs->getOpPtrVector().push_back(
6120 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
6121 mat_pos_at_pts_ptr));
6122
6123 if (defaultMaterial == BONEHOOKE) {
6124 if (!mwlsApprox) {
6125 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
6126 "mwlsApprox not allocated");
6127 }
6128 // calculate position at integration points
6129 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
6130 mat_grad_pos_at_pts_ptr =
6131 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
6132 feCouplingElasticLhs->getOpPtrVector().push_back(
6133 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
6134 mat_grad_pos_at_pts_ptr));
6135 feCouplingElasticLhs->getOpPtrVector().push_back(
6136 new OpGetCrackFrontDataGradientAtGaussPts(
6137 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6138 feCouplingElasticLhs->singularElement,
6139 feCouplingElasticLhs->invSJac));
6140
6141 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6142 feCouplingElasticLhs->getOpPtrVector().push_back(
6143 new MWLSApprox::OpMWLSRhoAtGaussPts(
6144 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6146 else {
6147 feCouplingElasticLhs->getOpPtrVector().push_back(
6149 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6151 feCouplingElasticLhs->getOpPtrVector().push_back(
6152 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
6153 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6155 }
6156 feCouplingElasticLhs->getOpPtrVector().push_back(
6157 new HookeElement::OpCalculateStiffnessScaledByDensityField(
6158 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6159 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
6160 rHo0));
6161 // Calculate spatial positions and gradients (of deformation) at
6162 // integration pts. This operator takes into account singularity at crack
6163 // front
6164 feCouplingElasticLhs->getOpPtrVector().push_back(
6165 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
6166 "SPATIAL_POSITION",
6167 data_hooke_element_at_pts));
6168 feCouplingElasticLhs->getOpPtrVector().push_back(
6169 new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
6170 "SPATIAL_POSITION",
6171 data_hooke_element_at_pts));
6172 // Transform base functions to get singular base functions at crack front
6173 feCouplingElasticLhs->getOpPtrVector().push_back(
6174 new OpTransfromSingularBaseFunctions(
6175 feCouplingElasticLhs->singularElement, feCouplingElasticLhs->detS,
6176 feCouplingElasticLhs->invSJac));
6177 feCouplingElasticLhs->getOpPtrVector().push_back(
6178 new HookeElement::OpAleLhs_dx_dx<1>("SPATIAL_POSITION",
6179 "SPATIAL_POSITION",
6180 data_hooke_element_at_pts));
6181 feCouplingElasticLhs->getOpPtrVector().push_back(
6182 new HookeElement::OpAleLhs_dx_dX<1>("SPATIAL_POSITION",
6183 "MESH_NODE_POSITIONS",
6184 data_hooke_element_at_pts));
6185
6186 feCouplingElasticLhs->getOpPtrVector().push_back(
6187 new HookeElement::OpAleLhsWithDensity_dx_dX(
6188 "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
6189 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6190 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
6191
6192 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
6193 if (addSingularity) {
6194
6195 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
6196 mwlsApprox, &mwlsApprox->singularInitialDisplacement);
6197
6198 feCouplingElasticLhs->getOpPtrVector().push_back(
6199 new OpAleLhsWithDensitySingularElement_dx_dX(
6200 "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
6201 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6202 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0,
6203 mat_singular_disp_ptr));
6204 }
6205
6206 } else {
6207
6208 feCouplingElasticLhs->getOpPtrVector().push_back(
6209 new HookeElement::OpCalculateHomogeneousStiffness<0>(
6210 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6211 data_hooke_element_at_pts));
6212 // Calculate spatial positions and gradients (of deformation) at
6213 // integration pts. This operator takes into account singularity at crack
6214 // front
6215 feCouplingElasticLhs->getOpPtrVector().push_back(
6216 new OpGetCrackFrontDataGradientAtGaussPts(
6217 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6218 feCouplingElasticLhs->singularElement,
6219 feCouplingElasticLhs->invSJac));
6220 feCouplingElasticLhs->getOpPtrVector().push_back(
6221 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
6222 "SPATIAL_POSITION",
6223 data_hooke_element_at_pts));
6224 feCouplingElasticLhs->getOpPtrVector().push_back(
6225 new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
6226 "SPATIAL_POSITION",
6227 data_hooke_element_at_pts));
6228 // Transform base functions to get singular base functions at crack front
6229 feCouplingElasticLhs->getOpPtrVector().push_back(
6230 new OpTransfromSingularBaseFunctions(
6231 feCouplingElasticLhs->singularElement, feCouplingElasticLhs->detS,
6232 feCouplingElasticLhs->invSJac));
6233 feCouplingElasticLhs->getOpPtrVector().push_back(
6234 new HookeElement::OpAleLhs_dx_dx<0>("SPATIAL_POSITION",
6235 "SPATIAL_POSITION",
6236 data_hooke_element_at_pts));
6237 feCouplingElasticLhs->getOpPtrVector().push_back(
6238 new HookeElement::OpAleLhs_dx_dX<0>("SPATIAL_POSITION",
6239 "MESH_NODE_POSITIONS",
6240 data_hooke_element_at_pts));
6241 }
6242 }
6243
6244 if (!materialFe) {
6245 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
6246 "Material element instance not created");
6247 }
6248 feCouplingMaterialLhs = boost::make_shared<CrackFrontElement>(
6251 feCouplingMaterialLhs->meshPositionsFieldName = "NONE";
6252 feCouplingMaterialLhs->addToRule = 0;
6253
6254 // calculate positions at integration points
6255 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(new MatrixDouble());
6256 feCouplingMaterialLhs->getOpPtrVector().push_back(
6257 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS",
6258 mat_pos_at_pts_ptr));
6259 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
6260 if (!onlyHooke && residualStressBlock != -1) {
6261 mat_grad_pos_at_pts_ptr =
6262 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
6263 feCouplingMaterialLhs->getOpPtrVector().push_back(
6264 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
6265 mat_grad_pos_at_pts_ptr));
6266 }
6267 boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr;
6268 if (!onlyHooke && residualStressBlock != -1) {
6269 space_grad_pos_at_pts_ptr =
6270 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
6271 feCouplingMaterialLhs->getOpPtrVector().push_back(
6272 new OpGetCrackFrontDataGradientAtGaussPts(
6273 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
6274 feCouplingMaterialLhs->singularElement,
6275 feCouplingMaterialLhs->invSJac));
6276 }
6277
6278 // assemble tangent matrix
6279 if (!onlyHooke) {
6280 feCouplingMaterialLhs->getOpPtrVector().push_back(
6281 new OpGetCrackFrontCommonDataAtGaussPts(
6282 "SPATIAL_POSITION", materialFe->commonData,
6283 feCouplingMaterialLhs->singularElement,
6284 feCouplingMaterialLhs->invSJac));
6285 feCouplingMaterialLhs->getOpPtrVector().push_back(
6287 "MESH_NODE_POSITIONS", materialFe->commonData));
6288 feCouplingMaterialLhs->getOpPtrVector().push_back(
6289 new OpTransfromSingularBaseFunctions(
6290 feCouplingMaterialLhs->singularElement, feCouplingMaterialLhs->detS,
6291 feCouplingMaterialLhs->invSJac));
6292 for (map<int, NonlinearElasticElement::BlockData>::iterator sit =
6293 materialFe->setOfBlocks.begin();
6294 sit != materialFe->setOfBlocks.end(); sit++) {
6295 feCouplingMaterialLhs->getOpPtrVector().push_back(
6297 "SPATIAL_POSITION", sit->second, materialFe->commonData,
6298 MATERIAL_TAG, true, true));
6299 feCouplingMaterialLhs->getOpPtrVector().push_back(
6301 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", sit->second,
6302 materialFe->commonData));
6303 feCouplingMaterialLhs->getOpPtrVector().push_back(
6305 "MESH_NODE_POSITIONS", "SPATIAL_POSITION", sit->second,
6306 materialFe->commonData));
6307 }
6308 } else {
6309
6310 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
6311 data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
6312 data_hooke_element_at_pts->forcesOnlyOnEntitiesRow = crackFrontNodes;
6313
6314 feCouplingMaterialLhs->getOpPtrVector().push_back(
6316 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
6317 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
6318
6319 if (defaultMaterial == BONEHOOKE) {
6320 if (!mwlsApprox) {
6321 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
6322 "mwlsApprox not allocated");
6323 }
6324
6325 feCouplingMaterialLhs->getOpPtrVector().push_back(
6326 new OpGetCrackFrontDataGradientAtGaussPts(
6327 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6328 feCouplingMaterialLhs->singularElement,
6329 feCouplingMaterialLhs->invSJac));
6330
6331 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6332 feCouplingMaterialLhs->getOpPtrVector().push_back(
6333 new MWLSApprox::OpMWLSRhoAtGaussPts(
6334 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6336 else {
6337 feCouplingMaterialLhs->getOpPtrVector().push_back(
6339 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6341 feCouplingMaterialLhs->getOpPtrVector().push_back(
6342 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
6343 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6345 }
6346
6347 feCouplingMaterialLhs->getOpPtrVector().push_back(
6348 new HookeElement::OpCalculateStiffnessScaledByDensityField(
6349 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6350 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
6351 rHo0));
6352 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
6353 feCouplingMaterialLhs->getOpPtrVector().push_back(
6354 new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
6355 "MESH_NODE_POSITIONS",
6356 data_hooke_element_at_pts));
6357 feCouplingMaterialLhs->getOpPtrVector().push_back(
6358 new HookeElement::OpCalculateStress<1>("MESH_NODE_POSITIONS",
6359 "MESH_NODE_POSITIONS",
6360 data_hooke_element_at_pts));
6361 feCouplingMaterialLhs->getOpPtrVector().push_back(
6362 new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
6363 "MESH_NODE_POSITIONS",
6364 data_hooke_element_at_pts));
6365 feCouplingMaterialLhs->getOpPtrVector().push_back(
6366 new HookeElement::OpCalculateEshelbyStress(
6367 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6368 data_hooke_element_at_pts));
6369 feCouplingMaterialLhs->getOpPtrVector().push_back(
6370 new OpTransfromSingularBaseFunctions(
6371 feCouplingMaterialLhs->singularElement,
6373 feCouplingMaterialLhs->getOpPtrVector().push_back(
6374 new HookeElement::OpAleLhs_dX_dX<1>("MESH_NODE_POSITIONS",
6375 "MESH_NODE_POSITIONS",
6376 data_hooke_element_at_pts));
6377 feCouplingMaterialLhs->getOpPtrVector().push_back(
6378 new HookeElement::OpAleLhsPre_dX_dx<1>("MESH_NODE_POSITIONS",
6379 "SPATIAL_POSITION",
6380 data_hooke_element_at_pts));
6381 feCouplingMaterialLhs->getOpPtrVector().push_back(
6382 new HookeElement::OpAleLhs_dX_dx("MESH_NODE_POSITIONS",
6383 "SPATIAL_POSITION",
6384 data_hooke_element_at_pts));
6385
6386 feCouplingMaterialLhs->getOpPtrVector().push_back(
6387 new HookeElement::OpAleLhsWithDensity_dX_dX(
6388 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6389 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6390 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
6391
6392 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
6393 if (addSingularity) {
6394
6395 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
6396 mwlsApprox, &mwlsApprox->singularInitialDisplacement);
6397
6398 feCouplingMaterialLhs->getOpPtrVector().push_back(
6399 new OpAleLhsWithDensitySingularElement_dX_dX(
6400 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6401 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6402 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0,
6403 mat_singular_disp_ptr));
6404 }
6405
6406 feCouplingMaterialLhs->getOpPtrVector().push_back(
6407 new OpLhsBoneExplicitDerivariveWithHooke_dX(
6408 *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6409 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
6410 mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
6412
6413 feCouplingMaterialLhs->getOpPtrVector().push_back(
6414 new OpLhsBoneExplicitDerivariveWithHooke_dx(
6415 *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
6416 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
6417 mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
6419
6420 } else {
6421
6422 feCouplingMaterialLhs->getOpPtrVector().push_back(
6423 new HookeElement::OpCalculateHomogeneousStiffness<0>(
6424 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
6425 data_hooke_element_at_pts));
6426 // Calculate spatial positions and gradients (of deformation) at
6427 // integration pts. This operator takes into account singularity at crack
6428 // front
6429 feCouplingMaterialLhs->getOpPtrVector().push_back(
6430 new OpGetCrackFrontDataGradientAtGaussPts(
6431 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6432 feCouplingMaterialLhs->singularElement,
6433 feCouplingMaterialLhs->invSJac));
6434 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
6435 feCouplingMaterialLhs->getOpPtrVector().push_back(
6436 new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
6437 "MESH_NODE_POSITIONS",
6438 data_hooke_element_at_pts));
6439 feCouplingMaterialLhs->getOpPtrVector().push_back(
6440 new HookeElement::OpCalculateStress<0>("MESH_NODE_POSITIONS",
6441 "MESH_NODE_POSITIONS",
6442 data_hooke_element_at_pts));
6443 feCouplingMaterialLhs->getOpPtrVector().push_back(
6444 new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
6445 "MESH_NODE_POSITIONS",
6446 data_hooke_element_at_pts));
6447 feCouplingMaterialLhs->getOpPtrVector().push_back(
6448 new HookeElement::OpCalculateEshelbyStress(
6449 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
6450 data_hooke_element_at_pts));
6451 feCouplingMaterialLhs->getOpPtrVector().push_back(
6452 new OpTransfromSingularBaseFunctions(
6453 feCouplingMaterialLhs->singularElement,
6455 feCouplingMaterialLhs->getOpPtrVector().push_back(
6456 new HookeElement::OpAleLhs_dX_dX<0>("MESH_NODE_POSITIONS",
6457 "MESH_NODE_POSITIONS",
6458 data_hooke_element_at_pts));
6459 feCouplingMaterialLhs->getOpPtrVector().push_back(
6460 new HookeElement::OpAleLhsPre_dX_dx<0>("MESH_NODE_POSITIONS",
6461 "SPATIAL_POSITION",
6462 data_hooke_element_at_pts));
6463 feCouplingMaterialLhs->getOpPtrVector().push_back(
6464 new HookeElement::OpAleLhs_dX_dx("MESH_NODE_POSITIONS",
6465 "SPATIAL_POSITION",
6466 data_hooke_element_at_pts));
6467 }
6468 }
6469
6470 if (residualStressBlock != -1) {
6471 if (!mwlsApprox) {
6472 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
6473 "mwlsApprox not allocated");
6474 }
6475
6476 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6477 feCouplingMaterialLhs->getOpPtrVector().push_back(
6478 new MWLSApprox::OpMWLSStressAtGaussPts(
6479 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6481 else {
6482 feCouplingMaterialLhs->getOpPtrVector().push_back(
6484 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6486 feCouplingMaterialLhs->getOpPtrVector().push_back(
6487 new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
6488 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6490 }
6491
6492 feCouplingMaterialLhs->getOpPtrVector().push_back(
6493 new MWLSApprox::OpMWLSSpatialStressLhs_DX(mat_grad_pos_at_pts_ptr,
6494 mwlsApprox));
6495 feCouplingMaterialLhs->getOpPtrVector().push_back(
6496 new MWLSApprox::OpMWLSMaterialStressLhs_DX(
6497 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
6498 &crackFrontNodes));
6499 feCouplingMaterialLhs->getOpPtrVector().push_back(
6500 new MWLSApprox::OpMWLSMaterialStressLhs_Dx(
6501 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
6502 &crackFrontNodes));
6503 }
6504
6506}
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
NonlinearElasticElement::BlockData BlockData
Definition: elasticity.cpp:74
boost::shared_ptr< NonlinearElasticElement > materialFe
bool onlyHooke
True if only Hooke material is applied.

◆ 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 4478 of file CrackPropagation.cpp.

4479 {
4481
4483
4485}
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

◆ 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

Definition at line 3736 of file CrackPropagation.cpp.

3737 {
3739
3740 // Create elastic element data structure
3741 elasticFe = boost::make_shared<NonlinearElasticElement>(mField, ELASTIC_TAG);
3742
3743 // Create elastic element finite element instance for residuals. Note that
3744 // this element approx. singularity at crack front.
3745 feRhs = boost::make_shared<CrackFrontElement>(
3748
3749 // Create finite element instance for assembly of tangent matrix
3750 feLhs = boost::make_shared<CrackFrontElement>(
3753
3754 // Arbitrary lagrangian formulation, mesh node positions are taken into
3755 // account by operators.
3756 feRhs->meshPositionsFieldName = "NONE";
3757 feLhs->meshPositionsFieldName = "NONE";
3758 feRhs->addToRule = 0;
3759 feLhs->addToRule = 0;
3760
3761 // Create operators to calculate finite element matrices for elastic element
3762 onlyHooke = true;
3764 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
3765 elasticFe, &elasticFe->setOfBlocks);
3766 {
3768 mField, BLOCKSET | MAT_ELASTICSET, it)) {
3769
3770 // Get data from block
3771 Mat_Elastic mydata;
3772 CHKERR it->getAttributeDataStructure(mydata);
3773 int id = it->getMeshsetId();
3774 EntityHandle meshset = it->getMeshset();
3775 CHKERR mField.get_moab().get_entities_by_type(
3776 meshset, MBTET, elasticFe->setOfBlocks[id].tEts, true);
3777 elasticFe->setOfBlocks[id].iD = id;
3778 elasticFe->setOfBlocks[id].E = mydata.data.Young;
3779 elasticFe->setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
3780
3781 // Print material properties of blocks after being read
3782 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\nMaterial block %d \n", id);
3783 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung's modulus %6.4e \n",
3784 mydata.data.Young);
3785 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPoisson's ratio %6.4e \n\n",
3786 mydata.data.Poisson);
3787
3788 // Set default material to elastic blocks. Create material instances
3789 switch (defaultMaterial) {
3790 case HOOKE:
3791 elasticFe->setOfBlocks[id].materialDoublePtr =
3792 boost::make_shared<Hooke<double>>();
3793 elasticFe->setOfBlocks[id].materialAdoublePtr =
3794 boost::make_shared<Hooke<adouble>>();
3795 break;
3796 case KIRCHHOFF:
3797 elasticFe->setOfBlocks[id].materialDoublePtr = boost::make_shared<
3799 double>>();
3800 elasticFe->setOfBlocks[id].materialAdoublePtr = boost::make_shared<
3802 adouble>>();
3803 onlyHooke = false;
3804 break;
3805 case NEOHOOKEAN:
3806 elasticFe->setOfBlocks[id].materialDoublePtr =
3807 boost::make_shared<NeoHookean<double>>();
3808 elasticFe->setOfBlocks[id].materialAdoublePtr =
3809 boost::make_shared<NeoHookean<adouble>>();
3810 onlyHooke = false;
3811 break;
3812 case BONEHOOKE:
3813 elasticFe->setOfBlocks[id].materialDoublePtr =
3814 boost::make_shared<Hooke<double>>();
3815 elasticFe->setOfBlocks[id].materialAdoublePtr =
3816 boost::make_shared<Hooke<adouble>>();
3817 break;
3818 default:
3819 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
3820 "Material type not yet implemented");
3821 break;
3822 }
3823 }
3824 if (onlyHookeFromOptions == PETSC_FALSE)
3825 onlyHooke = false;
3826
3827 elasticFe->commonData.spatialPositions = "SPATIAL_POSITION";
3828 elasticFe->commonData.meshPositions = "MESH_NODE_POSITIONS";
3829
3830 auto data_hooke_element_at_pts =
3831 boost::make_shared<HookeElement::DataAtIntegrationPts>();
3832
3833 // calculate material positions at integration points
3834 auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3835 feRhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3836 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
3837 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
3838 if (residualStressBlock != -1 || densityMapBlock != -1) {
3839 mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3840 feRhs->getOpPtrVector().push_back(
3841 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
3842 mat_grad_pos_at_pts_ptr));
3843 }
3844
3845 /// Rhs
3846 if (!onlyHooke) {
3847 // Calculate spatial positions and gradients (of deformation) at
3848 // integration pts. This operator takes into account singularity at crack
3849 // front
3850 feRhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
3851 "SPATIAL_POSITION", elasticFe->commonData, feRhs->singularElement,
3852 feRhs->invSJac));
3853 // Calculate material positions and gradients at integration pts.
3854 feRhs->getOpPtrVector().push_back(
3856 "MESH_NODE_POSITIONS", elasticFe->commonData));
3857 // Transfrom base functions to create singularity at crack front
3858 feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
3859 feRhs->singularElement, feRhs->detS, feRhs->invSJac));
3860 // Iterate over material blocks
3861 map<int, NonlinearElasticElement::BlockData>::iterator sit =
3862 elasticFe->setOfBlocks.begin();
3863 for (; sit != elasticFe->setOfBlocks.end(); sit++) {
3864 // Evaluate stress at integration pts.
3865 feRhs->getOpPtrVector().push_back(
3867 "SPATIAL_POSITION", sit->second, elasticFe->commonData,
3868 ELASTIC_TAG, false, true, false));
3869 // Assemble internal forces for right hand side
3870 feRhs->getOpPtrVector().push_back(
3872 "SPATIAL_POSITION", sit->second, elasticFe->commonData));
3873 }
3874 } else {
3875 feRhs->getOpPtrVector().push_back(
3877 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
3878 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
3879
3880 // Evaluate density at integration points fro material like bone
3881 if (defaultMaterial == BONEHOOKE) {
3882 if (!mwlsApprox) {
3883 // Create instance for moving least square approximation
3884 mwlsApprox = boost::make_shared<MWLSApprox>(mField);
3885 // Load mesh with stresses
3886 CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
3887 MeshsetsManager *block_manager_ptr;
3888 CHKERR mField.getInterface(block_manager_ptr);
3889 CHKERR block_manager_ptr->getEntitiesByDimension(
3890 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3891 Tag th_rho;
3892 if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
3893 th_rho) != MB_SUCCESS) {
3894 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3895 "Density tag name %s cannot be found. Please "
3896 "provide the correct name.",
3897 mwlsRhoTagName.c_str());
3898 }
3899 CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
3900 th_rho);
3901 CHKERR mwlsApprox->getValuesToNodes(th_rho);
3902 } else {
3903 mwlsApprox->tetsInBlock.clear();
3904 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
3905 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3906 }
3907 // Calculate spatial positions and gradients (of deformation) at
3908 // integration pts. This operator takes into account singularity at
3909 // crack front
3910 feRhs->getOpPtrVector().push_back(
3911 new OpGetCrackFrontDataGradientAtGaussPts(
3912 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
3913 feRhs->singularElement, feRhs->invSJac));
3914 // Evaluate density at integration pts.
3915 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
3916 feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
3917 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
3918 mwlsRhoTagName, true, true));
3919 else {
3920 feRhs->getOpPtrVector().push_back(
3922 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
3923 mwlsApprox));
3924 feRhs->getOpPtrVector().push_back(
3925 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
3926 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
3927 mwlsApprox, mwlsRhoTagName, true, false));
3928 }
3929
3930 feRhs->getOpPtrVector().push_back(
3931 new HookeElement::OpCalculateStiffnessScaledByDensityField(
3932 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
3933 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
3934 rHo0));
3935 feRhs->getOpPtrVector().push_back(
3936 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
3937 "SPATIAL_POSITION",
3938 data_hooke_element_at_pts));
3939 feRhs->getOpPtrVector().push_back(
3940 new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
3941 "SPATIAL_POSITION",
3942 data_hooke_element_at_pts));
3943 } else {
3944 feRhs->getOpPtrVector().push_back(
3945 new HookeElement::OpCalculateHomogeneousStiffness<0>(
3946 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
3947 data_hooke_element_at_pts));
3948 // Calculate spatial positions and gradients (of deformation) at
3949 // integration pts. This operator takes into account singularity at
3950 // crack front
3951 feRhs->getOpPtrVector().push_back(
3952 new OpGetCrackFrontDataGradientAtGaussPts(
3953 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
3954 feRhs->singularElement, feRhs->invSJac));
3955 feRhs->getOpPtrVector().push_back(
3956 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
3957 "SPATIAL_POSITION",
3958 data_hooke_element_at_pts));
3959 feRhs->getOpPtrVector().push_back(
3960 new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
3961 "SPATIAL_POSITION",
3962 data_hooke_element_at_pts));
3963 }
3964
3965 feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
3966 feRhs->singularElement, feRhs->detS, feRhs->invSJac));
3967 feRhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dx(
3968 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
3969 }
3970
3971 // Add operators if internal stresses are present
3972 if (residualStressBlock != -1) {
3973 if (!mwlsApprox) {
3974 // Create instance for moving least square approximation
3975 mwlsApprox = boost::make_shared<MWLSApprox>(mField);
3976 // Load mesh with stresses
3977 CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
3978 MeshsetsManager *block_manager_ptr;
3979 CHKERR mField.getInterface(block_manager_ptr);
3980 CHKERR block_manager_ptr->getEntitiesByDimension(
3981 residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3982 Tag th;
3983 if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
3984 th) != MB_SUCCESS) {
3985 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3986 "Internal stress tag name %s cannot be found. Please "
3987 "provide the correct name.",
3988 mwls_stress_tag_name.substr(4).c_str());
3989 }
3990 }
3991
3992 Tag th;
3993 if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
3994 th) != MB_SUCCESS) {
3995 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3996 "Internal stress tag name %s cannot be found. Please "
3997 "provide the correct name.",
3998 mwls_stress_tag_name.substr(4).c_str());
3999 }
4000 CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
4001 th);
4002 CHKERR mwlsApprox->getValuesToNodes(th);
4003
4004 // Evaluate internal stresses at integration pts.
4005 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4006 feRhs->getOpPtrVector().push_back(
4007 new MWLSApprox::OpMWLSStressAtGaussPts(
4008 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
4009 mwls_stress_tag_name, false));
4010 else {
4011 feRhs->getOpPtrVector().push_back(
4013 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
4014 mwlsApprox));
4015 feRhs->getOpPtrVector().push_back(
4016 new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
4017 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
4018 mwls_stress_tag_name, false));
4019 }
4020
4021 // Assemble internall stresses forces
4022 feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressRhs(
4023 mat_grad_pos_at_pts_ptr, mwlsApprox));
4025 feRhs->getOpPtrVector().push_back(
4026 new HookeElement::OpAnalyticalInternalAleStrain_dx<0>(
4027 "SPATIAL_POSITION", data_hooke_element_at_pts,
4029 boost::shared_ptr<MatrixDouble>(
4030 mwlsApprox, &mwlsApprox->mwlsMaterialPositions)));
4031 }
4032 }
4033
4034 /// Lhs
4035 if (!onlyHooke) {
4036 // Calculate spatial positions and gradients (of deformation) at
4037 // integration pts. This operator takes into account singularity at crack
4038 // front
4039 feLhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
4040 "SPATIAL_POSITION", elasticFe->commonData, feLhs->singularElement,
4041 feLhs->invSJac));
4042 // Calculate material positions and gradients at integration pts.
4043 feLhs->getOpPtrVector().push_back(
4045 "MESH_NODE_POSITIONS", elasticFe->commonData));
4046 // Transform base functions to get singular base functions at crack front
4047 feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4048 feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4049 // Iterate over material blocks
4050 map<int, NonlinearElasticElement::BlockData>::iterator sit =
4051 elasticFe->setOfBlocks.begin();
4052 for (; sit != elasticFe->setOfBlocks.end(); sit++) {
4053 // Evaluate stress at integration pts
4054 feLhs->getOpPtrVector().push_back(
4056 "SPATIAL_POSITION", sit->second, elasticFe->commonData,
4057 ELASTIC_TAG, true, true, false));
4058 // Assemble tanget matrix, derivative of spatial poisotions
4059 feLhs->getOpPtrVector().push_back(
4061 "SPATIAL_POSITION", "SPATIAL_POSITION", sit->second,
4062 elasticFe->commonData));
4063 }
4064 } else {
4065 // calculate material positions at integration points
4066 feLhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4067 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4068 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
4069 if (residualStressBlock != -1 || densityMapBlock != -1) {
4070 mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4071 feLhs->getOpPtrVector().push_back(
4072 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
4073 mat_grad_pos_at_pts_ptr));
4074 }
4075
4076 feLhs->getOpPtrVector().push_back(
4078 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
4079 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
4080
4081 if (defaultMaterial == BONEHOOKE) {
4082 if (!mwlsApprox) {
4083 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4084 "mwlsApprox not allocated");
4085 } else {
4086 mwlsApprox->tetsInBlock.clear();
4087 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
4088 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock,
4089 true); // FIXME: do we still need this?
4090 }
4091
4092 // Calculate spatial positions and gradients (of deformation) at
4093 // integration pts. This operator takes into account singularity at
4094 // crack front
4095 feLhs->getOpPtrVector().push_back(
4096 new OpGetCrackFrontDataGradientAtGaussPts(
4097 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4098 feLhs->singularElement, feLhs->invSJac));
4099
4100 // Evaluate density at integration pts.
4101 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4102 feLhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
4103 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs, mwlsApprox,
4104 mwlsRhoTagName, true, true));
4105 else {
4106 feLhs->getOpPtrVector().push_back(
4108 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
4109 mwlsApprox));
4110 feLhs->getOpPtrVector().push_back(
4111 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
4112 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
4113 mwlsApprox, mwlsRhoTagName, true, false));
4114 }
4115
4116 feLhs->getOpPtrVector().push_back(
4117 new HookeElement::OpCalculateStiffnessScaledByDensityField(
4118 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
4119 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
4120 rHo0));
4121 feLhs->getOpPtrVector().push_back(
4122 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
4123 "SPATIAL_POSITION",
4124 data_hooke_element_at_pts));
4125 feLhs->getOpPtrVector().push_back(
4126 new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
4127 "SPATIAL_POSITION",
4128 data_hooke_element_at_pts));
4129
4130 feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4131 feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4132 feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<1>(
4133 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
4134
4135 } else {
4136 feLhs->getOpPtrVector().push_back(
4137 new HookeElement::OpCalculateHomogeneousStiffness<0>(
4138 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
4139 data_hooke_element_at_pts));
4140 // Calculate spatial positions and gradients (of deformation) at
4141 // integration pts. This operator takes into account singularity at
4142 // crack front
4143 feLhs->getOpPtrVector().push_back(
4144 new OpGetCrackFrontDataGradientAtGaussPts(
4145 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4146 feLhs->singularElement, feLhs->invSJac));
4147 feLhs->getOpPtrVector().push_back(
4148 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
4149 "SPATIAL_POSITION",
4150 data_hooke_element_at_pts));
4151 feLhs->getOpPtrVector().push_back(
4152 new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
4153 "SPATIAL_POSITION",
4154 data_hooke_element_at_pts));
4155 // Transform base functions to get singular base functions at crack
4156 // front
4157 feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4158 feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4159 feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<0>(
4160 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
4161 }
4162 }
4163 }
4164
4165 // Assembly forces
4166
4168 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4169 // Surface forces
4170 {
4171 string fe_name_str = "FORCE_FE";
4172 surfaceForces->insert(fe_name_str, new NeumannForcesSurface(mField));
4174 it)) {
4175 CHKERR surfaceForces->at(fe_name_str)
4176 .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), true);
4177 }
4178
4179 const string block_set_force_name("FORCE");
4180 // search for block named FORCE and add its attributes to FORCE_FE element
4182 if (it->getName().compare(0, block_set_force_name.length(),
4183 block_set_force_name) == 0) {
4184 CHKERR surfaceForces->at(fe_name_str)
4185 .addForce("SPATIAL_POSITION", PETSC_NULL, (it->getMeshsetId()),
4186 true, true);
4187 }
4188 }
4189 }
4190 // assemble surface pressure
4192 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4193 {
4194 string fe_name_str = "PRESSURE_FE";
4195 surfacePressure->insert(fe_name_str, new NeumannForcesSurface(mField));
4196 // iterate over sidestep where pressure is applied
4198 mField, SIDESET | PRESSURESET, it)) {
4199 CHKERR surfacePressure->at(fe_name_str)
4200 .addPressure("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
4201 true);
4202 }
4203
4204 const string block_set_pressure_name("PRESSURE");
4206 if (it->getName().compare(0, block_set_pressure_name.length(),
4207 block_set_pressure_name) == 0) {
4208 CHKERR surfacePressure->at(fe_name_str)
4209 .addPressure("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
4210 true, true);
4211 }
4212 }
4213
4214 const string block_set_linear_pressure_name("LINEAR_PRESSURE");
4216 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
4217 block_set_linear_pressure_name) == 0) {
4218 CHKERR surfacePressure->at(fe_name_str)
4219 .addLinearPressure("SPATIAL_POSITION", PETSC_NULL,
4220 it->getMeshsetId(), true);
4221 }
4222 }
4223 }
4224 // assemble surface pressure (ALE)
4225 if (isPressureAle) {
4227 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4229 boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
4230 commonDataSurfacePressureAle->forcesOnlyOnEntitiesRow = crackFrontNodes;
4231 string fe_name_str = "PRESSURE_ALE";
4232 surfacePressureAle->insert(fe_name_str, new NeumannForcesSurface(mField));
4233 // iterate over sidesets where pressure is applied
4235 mField, SIDESET | PRESSURESET, it)) {
4236 CHKERR surfacePressureAle->at(fe_name_str)
4237 .addPressureAle("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
4238 commonDataSurfacePressureAle, "MATERIAL", PETSC_NULL,
4239 PETSC_NULL, it->getMeshsetId(), true);
4240 }
4241 }
4242
4243 if (areSpringsAle) {
4245 boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
4246
4248 boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
4250 boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
4252 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(mField);
4253 commonDataSpringsALE->forcesOnlyOnEntitiesRow = crackFrontNodes;
4254
4257 commonDataSpringsALE, "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
4258 "MATERIAL");
4259 }
4260
4261 // Implementation of spring element
4262 // Create new instances of face elements for springs
4263 feSpringLhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
4265 feSpringRhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
4267
4269 feSpringRhsPtr, "SPATIAL_POSITION",
4270 "MESH_NODE_POSITIONS");
4271
4272 bool use_eigen_pos =
4273 solveEigenStressProblem && (mwls_stress_tag_name == mwlsStressTagName);
4274
4275 bool use_eigen_pos_simple_contact =
4276 use_eigen_pos && useEigenPositionsSimpleContact;
4277
4278 // Set contact operators
4280 boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4282 boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4284 boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(mField);
4285
4286 if (!ignoreContact) {
4287 if (printContactState) {
4288 feRhsSimpleContact->contactStateVec =
4289 commonDataSimpleContact->gaussPtsStateVec;
4290 }
4291 contactProblem->setContactOperatorsRhs(
4292 feRhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4293 "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4294 "EIGEN_SPATIAL_POSITIONS");
4295 contactProblem->setMasterForceOperatorsRhs(
4296 feRhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4297 "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4298 "EIGEN_SPATIAL_POSITIONS");
4299
4300 contactProblem->setContactOperatorsLhs(
4301 feLhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4302 "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4303 "EIGEN_SPATIAL_POSITIONS");
4304
4305 contactProblem->setMasterForceOperatorsLhs(
4306 feLhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
4307 "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
4308 "EIGEN_SPATIAL_POSITIONS");
4309 }
4310
4311 // Close crack constrains
4314 boost::shared_ptr<BothSurfaceConstrains>(new BothSurfaceConstrains(
4315 mField, "LAMBDA_CLOSE_CRACK", "SPATIAL_POSITION"));
4316 closeCrackConstrains->aLpha = 1;
4317 }
4318
4319 // Set contact operators
4320
4321 if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
4322 // Set contact operators
4324 boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4326 boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4328 boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4330 boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
4331 mField);
4332
4333 commonDataSimpleContactALE->forcesOnlyOnEntitiesRow = crackFrontNodes;
4334
4335 contactProblem->setContactOperatorsRhsALEMaterial(
4337 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAMBDA_CONTACT",
4338 "MAT_CONTACT");
4339
4340 contactProblem->setContactOperatorsLhsALEMaterial(
4342 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAMBDA_CONTACT",
4343 "MAT_CONTACT");
4344
4345 contactProblem->setContactOperatorsLhsALE(
4347 "MESH_NODE_POSITIONS", "LAMBDA_CONTACT", use_eigen_pos_simple_contact,
4348 "EIGEN_SPATIAL_POSITIONS");
4349 }
4350
4352 boost::make_shared<MortarContactProblem::MortarContactElement>(
4353 mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
4354 "MESH_NODE_POSITIONS");
4356 boost::make_shared<MortarContactProblem::MortarContactElement>(
4357 mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
4358 "MESH_NODE_POSITIONS");
4360 boost::make_shared<MortarContactProblem::CommonDataMortarContact>(mField);
4361
4362 if (!ignoreContact) {
4363 if (printContactState) {
4364 feRhsMortarContact->contactStateVec =
4365 commonDataMortarContact->gaussPtsStateVec;
4366 }
4367 mortarContactProblemPtr->setContactOperatorsRhs(
4368 feRhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4369 "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4370
4371 mortarContactProblemPtr->setMasterForceOperatorsRhs(
4372 feRhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4373 "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4374
4375 mortarContactProblemPtr->setContactOperatorsLhs(
4376 feLhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4377 "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4378
4379 mortarContactProblemPtr->setMasterForceOperatorsLhs(
4380 feLhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
4381 "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4382 }
4383
4384 // set contact post proc operators
4385
4387 boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
4388
4390 boost::make_shared<MortarContactProblem::MortarContactElement>(
4391 mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
4392 "MESH_NODE_POSITIONS");
4393
4394 if (!ignoreContact) {
4395 CHKERR contactProblem->setContactOperatorsForPostProc(
4397 "SPATIAL_POSITION", "LAMBDA_CONTACT", contactPostProcMoab, false,
4398 use_eigen_pos_simple_contact, "EIGEN_SPATIAL_POSITIONS");
4399
4400 CHKERR mortarContactProblemPtr->setContactOperatorsForPostProc(
4402 "SPATIAL_POSITION", "LAMBDA_CONTACT", contactPostProcMoab, false,
4403 use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
4404 }
4405
4406 // assemble nodal forces
4407 nodalForces = boost::make_shared<boost::ptr_map<string, NodalForce>>();
4408 {
4409 string fe_name_str = "FORCE_FE";
4410 nodalForces->insert(fe_name_str, new NodalForce(mField));
4412 it)) {
4413 CHKERR nodalForces->at(fe_name_str)
4414 .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), false);
4415 }
4416 }
4417 // assemble edge forces
4418 edgeForces = boost::make_shared<boost::ptr_map<string, EdgeForce>>();
4419 {
4420 string fe_name_str = "FORCE_FE";
4421 edgeForces->insert(fe_name_str, new EdgeForce(mField));
4423 it)) {
4424 CHKERR edgeForces->at(fe_name_str)
4425 .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), false);
4426 }
4427 }
4428
4429 // Kinematic boundary conditions
4430 spatialDirichletBc = boost::shared_ptr<DirichletSpatialPositionsBc>(
4431 new DirichletSpatialPositionsBc(mField, "SPATIAL_POSITION", PETSC_NULL,
4432 PETSC_NULL, PETSC_NULL));
4433
4434// Add boundary conditions for displacements given by analytical function
4435#ifdef __ANALITICAL_DISPLACEMENT__
4436 analyticalDirichletBc = boost::shared_ptr<AnalyticalDirichletBC::DirichletBC>(
4438 mField, "SPATIAL_POSITION", PETSC_NULL, PETSC_NULL, PETSC_NULL));
4439#endif //__ANALITICAL_DISPLACEMENT__
4440
4441// Analytical forces
4442#ifdef __ANALITICAL_TRACTION__
4443 MeshsetsManager *meshset_manager_ptr;
4444 CHKERR mField.getInterface(meshset_manager_ptr);
4445 if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
4448 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4449 {
4450 string fe_name_str = "ANALITICAL_TRACTION";
4451 analiticalSurfaceElement->insert(fe_name_str,
4453 analiticalSurfaceElement->at(fe_name_str).fe.addToRule = 1;
4454 const CubitMeshSets *cubit_meshset_ptr;
4455 meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
4456 &cubit_meshset_ptr);
4457 Range faces;
4458 CHKERR meshset_manager_ptr->getEntitiesByDimension(
4459 cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, faces, true);
4460 analiticalSurfaceElement->at(fe_name_str)
4461 .analyticalForceOp.push_back(new AnalyticalForces());
4462 analiticalSurfaceElement->at(fe_name_str)
4463 .fe.getOpPtrVector()
4464 .push_back(new OpAnalyticalSpatialTraction(
4467 "SPATIAL_POSITION", faces,
4468 analiticalSurfaceElement->at(fe_name_str).methodsOp,
4469 analiticalSurfaceElement->at(fe_name_str).analyticalForceOp));
4470 }
4471 }
4472 }
4473#endif //__ANALITICAL_TRACTION__
4474
4476}
@ PRESSURESET
Definition: definitions.h:152
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ FORCESET
Definition: definitions.h:151
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Structure used to enforce analytical boundary conditions.
Set Dirichlet boundary conditions on spatial displacements.
Force on edges and lines.
Definition: EdgeForce.hpp:13
boost::shared_ptr< MortarContactProblem::MortarContactElement > fePostProcMortarContact
PetscBool onlyHookeFromOptions
True if only Hooke material is applied.
boost::shared_ptr< SimpleContactProblem::CommonDataSimpleContact > commonDataSimpleContact
boost::shared_ptr< MortarContactProblem > mortarContactProblemPtr
boost::shared_ptr< MortarContactProblem::CommonDataMortarContact > commonDataMortarContact
static FTensor::Tensor2_symmetric< double, 3 > analyticalStrainFunction(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords)
PetscBool solveEigenStressProblem
Solve eigen problem.
boost::shared_ptr< MetaSpringBC::DataAtIntegrationPtsSprings > commonDataSpringsALE
boost::shared_ptr< SimpleContactProblem::SimpleContactElement > fePostProcSimpleContact
boost::shared_ptr< ContactSearchKdTree::ContactCommonData_multiIndex > contactSearchMultiIndexPtr
boost::shared_ptr< SimpleContactProblem::CommonDataSimpleContact > commonDataSimpleContactALE
boost::shared_ptr< SimpleContactProblem > contactProblem
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode setSpringOperatorsMaterial(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr_dx, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr_dX, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, boost::shared_ptr< DataAtIntegrationPtsSprings > data_at_integration_pts, const std::string field_name, const std::string mesh_nodals_positions, std::string side_fe_name)
Implementation of spring element. Set operators to calculate LHS and RHS.
virtual moab::Interface & get_moab()=0
this struct keeps basic methods for moab meshset about material and boundary conditions
int getMeshsetId() const
get meshset id as it set in preprocessing software
Elastic material data structure.
Finite element and operators to apply force/pressures applied to surfaces.
Force applied to nodes.
Definition: NodalForce.hpp:13
Implementation of elastic (non-linear) St. Kirchhoff equation.

◆ assembleMaterialForcesDM()

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

create material element instance

Definition at line 5408 of file CrackPropagation.cpp.

5409 {
5411 if (!elasticFe)
5412 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
5413 "Elastic element instance not created");
5414
5415 // Create element
5416 materialFe = boost::shared_ptr<NonlinearElasticElement>(
5418
5419 // Create material data blocks
5420 for (std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
5421 elasticFe->setOfBlocks.begin();
5422 sit != elasticFe->setOfBlocks.end(); sit++) {
5423 materialFe->setOfBlocks[sit->first] = sit->second;
5424 materialFe->setOfBlocks[sit->first].forcesOnlyOnEntitiesRow =
5426 }
5427 materialFe->commonData.spatialPositions = "SPATIAL_POSITION";
5428 materialFe->commonData.meshPositions = "MESH_NODE_POSITIONS";
5429
5430 // create finite element instances for the right hand side and left hand side
5431 feMaterialRhs = boost::make_shared<CrackFrontElement>(
5434 feMaterialLhs = boost::make_shared<CrackFrontElement>(
5437 feMaterialRhs->meshPositionsFieldName = "NONE";
5438 feMaterialLhs->meshPositionsFieldName = "NONE";
5439 feMaterialRhs->addToRule = 0;
5440 feMaterialLhs->addToRule = 0;
5441
5443 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
5444 elasticFe, &elasticFe->setOfBlocks);
5445 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
5446 data_hooke_element_at_pts(new HookeElement::DataAtIntegrationPts());
5447 data_hooke_element_at_pts->forcesOnlyOnEntitiesRow = crackFrontNodes;
5448
5449 // calculate position at integration points
5450 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(new MatrixDouble());
5451 feMaterialRhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5452 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5453 feMaterialLhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
5454 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5455 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
5456 if (!onlyHooke && (residualStressBlock != -1 || densityMapBlock != -1)) {
5457 mat_grad_pos_at_pts_ptr =
5458 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
5459 feMaterialRhs->getOpPtrVector().push_back(
5460 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
5461 mat_grad_pos_at_pts_ptr));
5462 feMaterialLhs->getOpPtrVector().push_back(
5463 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
5464 mat_grad_pos_at_pts_ptr));
5465 }
5466 boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr;
5467 if (!onlyHooke && (residualStressBlock != -1 || densityMapBlock != -1)) {
5468 space_grad_pos_at_pts_ptr =
5469 boost::shared_ptr<MatrixDouble>(new MatrixDouble());
5470 feMaterialRhs->getOpPtrVector().push_back(
5471 new OpGetCrackFrontDataGradientAtGaussPts(
5472 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
5473 feMaterialRhs->singularElement, feMaterialRhs->invSJac));
5474 feMaterialLhs->getOpPtrVector().push_back(
5475 new OpGetCrackFrontDataGradientAtGaussPts(
5476 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
5477 feMaterialLhs->singularElement, feMaterialLhs->invSJac));
5478 }
5479
5480 /// Operators to assemble rhs
5481 if (!onlyHooke) {
5482 feMaterialRhs->getOpPtrVector().push_back(
5483 new OpGetCrackFrontCommonDataAtGaussPts(
5484 "SPATIAL_POSITION", materialFe->commonData,
5485 feMaterialRhs->singularElement, feMaterialRhs->invSJac));
5486 feMaterialRhs->getOpPtrVector().push_back(
5488 "MESH_NODE_POSITIONS", materialFe->commonData));
5489
5490 feMaterialRhs->getOpPtrVector().push_back(
5491 new OpTransfromSingularBaseFunctions(feMaterialRhs->singularElement,
5492 feMaterialRhs->detS,
5493 feMaterialRhs->invSJac));
5494 for (map<int, NonlinearElasticElement::BlockData>::iterator sit =
5495 materialFe->setOfBlocks.begin();
5496 sit != materialFe->setOfBlocks.end(); sit++) {
5497 feMaterialRhs->getOpPtrVector().push_back(
5499 "SPATIAL_POSITION", sit->second, materialFe->commonData,
5500 MATERIAL_TAG, false, true));
5501 feMaterialRhs->getOpPtrVector().push_back(
5503 "MESH_NODE_POSITIONS", sit->second, materialFe->commonData));
5504 }
5505
5506 } else { // HOOKE
5507
5508 feMaterialRhs->getOpPtrVector().push_back(
5510 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5511 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
5512
5513 feMaterialRhs->getOpPtrVector().push_back(
5514 new OpGetCrackFrontDataGradientAtGaussPts(
5515 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5516 feMaterialRhs->singularElement, feMaterialRhs->invSJac));
5517 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
5518
5519 if (defaultMaterial == BONEHOOKE) {
5520 if (!mwlsApprox) {
5521 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5522 "mwlsApprox not allocated");
5523 }
5524
5525 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5526 feMaterialRhs->getOpPtrVector().push_back(
5527 new MWLSApprox::OpMWLSRhoAtGaussPts(
5528 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5529 mwlsApprox, mwlsRhoTagName, true, true));
5530 else {
5531 feMaterialRhs->getOpPtrVector().push_back(
5533 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5534 mwlsApprox));
5535 feMaterialRhs->getOpPtrVector().push_back(
5536 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5537 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5538 mwlsApprox, mwlsRhoTagName, true, false));
5539 }
5540
5541 feMaterialRhs->getOpPtrVector().push_back(
5542 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5543 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5544 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
5545 rHo0));
5546 feMaterialRhs->getOpPtrVector().push_back(
5547 new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5548 "MESH_NODE_POSITIONS",
5549 data_hooke_element_at_pts));
5550 feMaterialRhs->getOpPtrVector().push_back(
5551 new HookeElement::OpCalculateStress<1>("MESH_NODE_POSITIONS",
5552 "MESH_NODE_POSITIONS",
5553 data_hooke_element_at_pts));
5554
5555 } else {
5556 feMaterialRhs->getOpPtrVector().push_back(
5557 new HookeElement::OpCalculateHomogeneousStiffness<0>(
5558 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5559 data_hooke_element_at_pts));
5560 feMaterialRhs->getOpPtrVector().push_back(
5561 new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5562 "MESH_NODE_POSITIONS",
5563 data_hooke_element_at_pts));
5564 feMaterialRhs->getOpPtrVector().push_back(
5565 new HookeElement::OpCalculateStress<0>("MESH_NODE_POSITIONS",
5566 "MESH_NODE_POSITIONS",
5567 data_hooke_element_at_pts));
5568 }
5569 feMaterialRhs->getOpPtrVector().push_back(
5570 new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
5571 "MESH_NODE_POSITIONS",
5572 data_hooke_element_at_pts));
5573
5574 feMaterialRhs->getOpPtrVector().push_back(
5575 new HookeElement::OpCalculateEshelbyStress("MESH_NODE_POSITIONS",
5576 "MESH_NODE_POSITIONS",
5577 data_hooke_element_at_pts));
5578 feMaterialRhs->getOpPtrVector().push_back(
5579 new OpTransfromSingularBaseFunctions(feMaterialRhs->singularElement,
5580 feMaterialRhs->detS,
5581 feMaterialRhs->invSJac));
5582
5583 feMaterialRhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dX(
5584 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5585 data_hooke_element_at_pts));
5586
5587 if (defaultMaterial == BONEHOOKE) {
5588 feMaterialRhs->getOpPtrVector().push_back(
5589 new OpRhsBoneExplicitDerivariveWithHooke(
5590 *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5591 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->tetsInBlock, rHo0,
5593 }
5594 }
5595
5596 if (residualStressBlock != -1) {
5597 if (!mwlsApprox) {
5598 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5599 "mwlsApprox not allocated");
5600 } else {
5601 mwlsApprox->tetsInBlock.clear();
5602 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
5603 residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
5604 }
5605
5606 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5607 feMaterialRhs->getOpPtrVector().push_back(
5608 new MWLSApprox::OpMWLSStressAtGaussPts(
5609 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5610 mwlsApprox, mwlsStressTagName, false));
5611 else {
5612 feMaterialRhs->getOpPtrVector().push_back(
5614 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5615 mwlsApprox));
5616 feMaterialRhs->getOpPtrVector().push_back(
5617 new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
5618 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialRhs,
5619 mwlsApprox, mwlsStressTagName, false));
5620 }
5621
5622 feMaterialRhs->getOpPtrVector().push_back(
5623 new MWLSApprox::OpMWLSMaterialStressRhs(space_grad_pos_at_pts_ptr,
5624 mat_grad_pos_at_pts_ptr,
5627 feMaterialRhs->getOpPtrVector().push_back(
5628 new HookeElement::OpAnalyticalInternalAleStrain_dX<0>(
5629 "MESH_NODE_POSITIONS", data_hooke_element_at_pts,
5631 boost::shared_ptr<MatrixDouble>(
5632 mwlsApprox, &mwlsApprox->mwlsMaterialPositions)));
5633 }
5634 }
5635
5636// Analytical forces
5637#ifdef __ANALITICAL_TRACTION__
5638 {
5641 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>(
5642 // analiticalSurfaceElement,
5644 feMaterialAnaliticalTraction->addToRule = 1;
5645 MeshsetsManager *meshset_manager_ptr;
5646 CHKERR mField.getInterface(meshset_manager_ptr);
5647 if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
5648 const CubitMeshSets *cubit_meshset_ptr;
5649 meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
5650 &cubit_meshset_ptr);
5651 Range faces;
5652 CHKERR meshset_manager_ptr->getEntitiesByDimension(
5653 cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, faces, true);
5654 feMaterialAnaliticalTraction->getOpPtrVector().push_back(
5655 new OpAnalyticalMaterialTraction(
5658 "MESH_NODE_POSITIONS", faces,
5659 analiticalSurfaceElement->at("ANALITICAL_TRACTION").methodsOp,
5660 analiticalSurfaceElement->at("ANALITICAL_TRACTION")
5661 .analyticalForceOp,
5662 &crackFrontNodes));
5663 }
5664 }
5665 }
5666#endif //__ANALITICAL_TRACTION__
5667
5668 // assemble tangent matrix
5669 if (!onlyHooke) {
5670 feMaterialLhs->getOpPtrVector().push_back(
5671 new OpGetCrackFrontCommonDataAtGaussPts(
5672 "SPATIAL_POSITION", materialFe->commonData,
5673 feMaterialLhs->singularElement, feMaterialLhs->invSJac));
5674 feMaterialLhs->getOpPtrVector().push_back(
5676 "MESH_NODE_POSITIONS", materialFe->commonData));
5677 feMaterialLhs->getOpPtrVector().push_back(
5678 new OpTransfromSingularBaseFunctions(
5679 feMaterialLhs->singularElement, feMaterialLhs->detS,
5681 ->invSJac //,AINSWORTH_LOBATTO_BASE//,AINSWORTH_LEGENDRE_BASE
5682 ));
5683 for (auto sit = materialFe->setOfBlocks.begin();
5684 sit != materialFe->setOfBlocks.end(); sit++) {
5685 feMaterialLhs->getOpPtrVector().push_back(
5687 "SPATIAL_POSITION", sit->second, materialFe->commonData,
5688 MATERIAL_TAG, true, true));
5689 feMaterialLhs->getOpPtrVector().push_back(
5691 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS", sit->second,
5692 materialFe->commonData));
5693 }
5694 } else { // HOOKE
5695 feMaterialLhs->getOpPtrVector().push_back(
5697 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5698 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
5699 feMaterialLhs->getOpPtrVector().push_back(
5700 new OpGetCrackFrontDataGradientAtGaussPts(
5701 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5702 feMaterialLhs->singularElement, feMaterialLhs->invSJac));
5703 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
5704
5705 if (defaultMaterial == BONEHOOKE) {
5706 if (!mwlsApprox) {
5707 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5708 "mwlsApprox not allocated");
5709 }
5710
5711 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5712 feMaterialLhs->getOpPtrVector().push_back(
5713 new MWLSApprox::OpMWLSRhoAtGaussPts(
5714 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5715 mwlsApprox, mwlsRhoTagName, true, true));
5716 else {
5717 feMaterialLhs->getOpPtrVector().push_back(
5719 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5720 mwlsApprox));
5721 feMaterialLhs->getOpPtrVector().push_back(
5722 new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
5723 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5724 mwlsApprox, mwlsRhoTagName, true, true));
5725 }
5726
5727 feMaterialLhs->getOpPtrVector().push_back(
5728 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5729 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5730 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
5731 rHo0));
5732
5733 feMaterialLhs->getOpPtrVector().push_back(
5734 new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5735 "MESH_NODE_POSITIONS",
5736 data_hooke_element_at_pts));
5737
5738 feMaterialLhs->getOpPtrVector().push_back(
5739 new HookeElement::OpCalculateStress<1>("MESH_NODE_POSITIONS",
5740 "MESH_NODE_POSITIONS",
5741 data_hooke_element_at_pts));
5742
5743 feMaterialLhs->getOpPtrVector().push_back(
5744 new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
5745 "MESH_NODE_POSITIONS",
5746 data_hooke_element_at_pts));
5747 feMaterialLhs->getOpPtrVector().push_back(
5748 new HookeElement::OpCalculateEshelbyStress(
5749 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5750 data_hooke_element_at_pts));
5751 feMaterialLhs->getOpPtrVector().push_back(
5752 new OpTransfromSingularBaseFunctions(feMaterialLhs->singularElement,
5753 feMaterialLhs->detS,
5754 feMaterialLhs->invSJac));
5755 feMaterialLhs->getOpPtrVector().push_back(
5756 new HookeElement::OpAleLhs_dX_dX<1>("MESH_NODE_POSITIONS",
5757 "MESH_NODE_POSITIONS",
5758 data_hooke_element_at_pts));
5759
5760 feMaterialLhs->getOpPtrVector().push_back(
5761 new HookeElement::OpAleLhsWithDensity_dX_dX(
5762 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5763 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5764 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0));
5765
5766 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr = nullptr;
5767 if (addSingularity) {
5768
5769 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
5770 mwlsApprox, &mwlsApprox->singularInitialDisplacement);
5771
5772 feMaterialLhs->getOpPtrVector().push_back(
5773 new OpAleLhsWithDensitySingularElement_dX_dX(
5774 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5775 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5776 mwlsApprox->diffRhoAtGaussPts, nBone, rHo0,
5777 mat_singular_disp_ptr));
5778 }
5779
5780 feMaterialLhs->getOpPtrVector().push_back(
5781 new OpLhsBoneExplicitDerivariveWithHooke_dX(
5782 *data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts,
5783 mwlsApprox->diffRhoAtGaussPts, mwlsApprox->diffDiffRhoAtGaussPts,
5784 mwlsApprox->singularInitialDisplacement, mwlsApprox->tetsInBlock,
5786
5787 } else {
5788 feMaterialLhs->getOpPtrVector().push_back(
5789 new HookeElement::OpCalculateHomogeneousStiffness<0>(
5790 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
5791 data_hooke_element_at_pts));
5792
5793 feMaterialLhs->getOpPtrVector().push_back(
5794 new HookeElement::OpCalculateStrainAle("MESH_NODE_POSITIONS",
5795 "MESH_NODE_POSITIONS",
5796 data_hooke_element_at_pts));
5797 feMaterialLhs->getOpPtrVector().push_back(
5798 new HookeElement::OpCalculateStress<0>("MESH_NODE_POSITIONS",
5799 "MESH_NODE_POSITIONS",
5800 data_hooke_element_at_pts));
5801
5802 feMaterialLhs->getOpPtrVector().push_back(
5803 new HookeElement::OpCalculateEnergy("MESH_NODE_POSITIONS",
5804 "MESH_NODE_POSITIONS",
5805 data_hooke_element_at_pts));
5806 feMaterialLhs->getOpPtrVector().push_back(
5807 new HookeElement::OpCalculateEshelbyStress(
5808 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5809 data_hooke_element_at_pts));
5810 feMaterialLhs->getOpPtrVector().push_back(
5811 new OpTransfromSingularBaseFunctions(feMaterialLhs->singularElement,
5812 feMaterialLhs->detS,
5813 feMaterialLhs->invSJac));
5814 feMaterialLhs->getOpPtrVector().push_back(
5815 new HookeElement::OpAleLhs_dX_dX<0>("MESH_NODE_POSITIONS",
5816 "MESH_NODE_POSITIONS",
5817 data_hooke_element_at_pts));
5818 }
5819 }
5820 if (residualStressBlock != -1) {
5821 if (!mwlsApprox) {
5822 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
5823 "mwlsApprox not allocated");
5824 }
5825
5826 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5827 feMaterialLhs->getOpPtrVector().push_back(
5828 new MWLSApprox::OpMWLSStressAtGaussPts(
5829 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5831 else {
5832 feMaterialLhs->getOpPtrVector().push_back(
5834 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5835 mwlsApprox));
5836 feMaterialLhs->getOpPtrVector().push_back(
5837 new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
5838 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feMaterialLhs,
5839 mwlsApprox, mwlsStressTagName, false));
5840 }
5841
5842 feMaterialLhs->getOpPtrVector().push_back(
5843 new MWLSApprox::OpMWLSMaterialStressLhs_DX(
5844 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, mwlsApprox,
5845 &crackFrontNodes));
5846 }
5847
5849}
boost::shared_ptr< NeumannForcesSurface::MyTriangleFE > feMaterialAnaliticalTraction
Surface elment to calculate tractions in material space.
structure grouping operators and data used for calculation of nonlinear elastic element

◆ 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 5852 of file CrackPropagation.cpp.

5853 {
5855 if (!elasticFe)
5856 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
5857 "Elastic element instance not created");
5858
5859 smootherFe = boost::make_shared<Smoother>(mField);
5860 volumeLengthAdouble = boost::shared_ptr<VolumeLengthQuality<adouble>>(
5863 smootherGamma));
5864 volumeLengthDouble = boost::shared_ptr<VolumeLengthQuality<double>>(
5867 smootherGamma));
5868
5869 // set block data
5870 for (std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
5871 elasticFe->setOfBlocks.begin();
5872 sit != elasticFe->setOfBlocks.end(); sit++) {
5873 // E = fmax(E,sit->second.E);
5874 smootherFe->setOfBlocks[0].tEts.merge(sit->second.tEts);
5875 }
5876 smootherFe->setOfBlocks[0].materialDoublePtr = volumeLengthDouble;
5877 smootherFe->setOfBlocks[0].materialAdoublePtr = volumeLengthAdouble;
5878 smootherFe->setOfBlocks[0].forcesOnlyOnEntitiesRow = crackFrontNodes;
5879
5880 // set element data
5881 smootherFe->commonData.spatialPositions = "MESH_NODE_POSITIONS";
5882 smootherFe->commonData.meshPositions = "NONE";
5883
5884 // mesh field name
5885 smootherFe->feRhs.meshPositionsFieldName = "NONE";
5886 smootherFe->feLhs.meshPositionsFieldName = "NONE";
5887 smootherFe->feRhs.addToRule = 0;
5888 smootherFe->feLhs.addToRule = 0;
5889
5890 // Create finite element instances for the right hand side and left hand side
5891 feSmootherRhs = smootherFe->feRhsPtr;
5892 feSmootherLhs = smootherFe->feLhsPtr;
5893
5894 // Smoother right hand side
5895 feSmootherRhs->getOpPtrVector().push_back(
5897 "MESH_NODE_POSITIONS", smootherFe->commonData));
5898 map<int, NonlinearElasticElement::BlockData>::iterator sit =
5899 smootherFe->setOfBlocks.begin();
5900 feSmootherRhs->getOpPtrVector().push_back(new Smoother::OpJacobianSmoother(
5901 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
5902 smootherFe->commonData, SMOOTHING_TAG, false));
5903 feSmootherRhs->getOpPtrVector().push_back(new Smoother::OpRhsSmoother(
5904 "MESH_NODE_POSITIONS", sit->second, smootherFe->commonData,
5905 smootherFe->smootherData));
5906
5907 // Smoother left hand side
5908 feSmootherLhs->getOpPtrVector().push_back(
5910 "MESH_NODE_POSITIONS", smootherFe->commonData));
5911 feSmootherLhs->getOpPtrVector().push_back(new Smoother::OpJacobianSmoother(
5912 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
5913 smootherFe->commonData, SMOOTHING_TAG, true));
5914 feSmootherLhs->getOpPtrVector().push_back(new Smoother::OpLhsSmoother(
5915 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
5916 smootherFe->setOfBlocks.at(0), smootherFe->commonData,
5917 smootherFe->smootherData, "LAMBDA_CRACKFRONT_AREA_TANGENT"));
5918
5919 // Constrains at crack front
5920 tangentConstrains = boost::shared_ptr<
5923 mField, "LAMBDA_CRACKFRONT_AREA_TANGENT"));
5924
5925 Range contact_faces;
5926 contact_faces.merge(contactSlaveFaces);
5927 contact_faces.merge(contactMasterFaces);
5928
5929 // Surface sliding constrains
5930 surfaceConstrain.clear();
5932 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
5933 new FaceOrientation(false, contact_faces, mapBitLevel["mesh_cut"]));
5934 for (auto id : ids) {
5935 if (id != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
5936 surfaceConstrain[id] = boost::make_shared<SurfaceSlidingConstrains>(
5938 surfaceConstrain[id]->setOperators(
5940 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(id),
5941 "MESH_NODE_POSITIONS", &smootherAlpha);
5942 }
5943 }
5944 // Add crack surface sliding constrains
5946 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
5947 new FaceOrientation(true, contact_faces, mapBitLevel["mesh_cut"]));
5948 surfaceConstrain[getInterface<CPMeshCut>()->getCrackSurfaceId()] =
5949 boost::make_shared<SurfaceSlidingConstrains>(mField, *crackOrientation);
5950 surfaceConstrain[getInterface<CPMeshCut>()->getCrackSurfaceId()]
5951 ->setOperators(SURFACE_SLIDING_TAG,
5952 "LAMBDA_SURFACE" +
5953 boost::lexical_cast<std::string>(
5954 getInterface<CPMeshCut>()->getCrackSurfaceId()),
5955 "MESH_NODE_POSITIONS", &smootherAlpha);
5956
5957 edgeConstrains.clear();
5958 auto get_edges_block_set = [this]() {
5959 return getInterface<CPMeshCut>()->getEdgesBlockSet();
5960 };
5961 if (get_edges_block_set()) {
5962 edgeConstrains[get_edges_block_set()] =
5963 boost::shared_ptr<EdgeSlidingConstrains>(
5965 edgeConstrains[get_edges_block_set()]->setOperators(
5966 EDGE_SLIDING_TAG, "LAMBDA_EDGE", "MESH_NODE_POSITIONS", &smootherAlpha);
5967 }
5968
5969 bothSidesConstrains = boost::shared_ptr<BothSurfaceConstrains>(
5970 new BothSurfaceConstrains(mField));
5972
5973 bothSidesContactConstrains = boost::shared_ptr<BothSurfaceConstrains>(
5974 new BothSurfaceConstrains(mField, "LAMBDA_BOTH_SIDES_CONTACT"));
5977
5978 Range fix_material_nodes;
5979 fixMaterialEnts = boost::make_shared<DirichletFixFieldAtEntitiesBc>(
5980 mField, "MESH_NODE_POSITIONS", fix_material_nodes);
5981 const Field_multiIndex *fields_ptr;
5982 CHKERR mField.get_fields(&fields_ptr);
5983 for (auto f : *fields_ptr) {
5984 if (f->getName().compare(0, 6, "LAMBDA") == 0 &&
5985 f->getName() != "LAMBDA_ARC_LENGTH" &&
5986 f->getName() != "LAMBDA_CONTACT" &&
5987 f->getName() != "LAMBDA_CLOSE_CRACK") {
5988 fixMaterialEnts->fieldNames.push_back(f->getName());
5989 }
5990 }
5991
5992 // Create griffith force finite element instances and operators
5993 griffithForceElement = boost::make_shared<GriffithForceElement>(mField);
5994 griffithForceElement->blockData[0].gc = gC;
5995 griffithForceElement->blockData[0].E = gC;
5996 griffithForceElement->blockData[0].r = 1.0;
5997 griffithForceElement->blockData[0].frontNodes = crackFrontNodes;
5999 griffithForceElement->blockData[0]);
6000 gC = griffithForceElement->blockData[0].gc;
6001 griffithE = griffithForceElement->blockData[0].E;
6002 griffithR = griffithForceElement->blockData[0].r;
6003
6006
6007 feGriffithForceRhs->getOpPtrVector().push_back(
6008 new GriffithForceElement::OpCalculateGriffithForce(
6010 griffithForceElement->commonData));
6011 feGriffithForceRhs->getOpPtrVector().push_back(
6012 new GriffithForceElement::OpRhs(GRIFFITH_FORCE_TAG,
6013 griffithForceElement->blockData[0],
6014 griffithForceElement->commonData));
6015 feGriffithForceLhs->getOpPtrVector().push_back(
6016 new GriffithForceElement::OpLhs(GRIFFITH_FORCE_TAG,
6017 griffithForceElement->blockData[0],
6018 griffithForceElement->commonData));
6019
6020 // Creating Griffith constrains finite element instances and operators
6022 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrainsDelta>(
6023 new GriffithForceElement::MyTriangleFEConstrainsDelta(
6024 mField, "LAMBDA_CRACKFRONT_AREA"));
6025 feGriffithConstrainsDelta->getOpPtrVector().push_back(
6026 new GriffithForceElement::OpConstrainsDelta(
6028 griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA",
6029 feGriffithConstrainsDelta->deltaVec));
6031 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>(
6032 new GriffithForceElement::MyTriangleFEConstrains(
6033 mField, "LAMBDA_CRACKFRONT_AREA",
6034 griffithForceElement->blockData[0],
6035 feGriffithConstrainsDelta->deltaVec));
6037 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>(
6038 new GriffithForceElement::MyTriangleFEConstrains(
6039 mField, "LAMBDA_CRACKFRONT_AREA",
6040 griffithForceElement->blockData[0],
6041 feGriffithConstrainsDelta->deltaVec));
6042 feGriffithConstrainsRhs->getOpPtrVector().push_back(
6043 new GriffithForceElement::OpConstrainsRhs(
6045 griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA"));
6046 feGriffithConstrainsLhs->getOpPtrVector().push_back(
6047 new GriffithForceElement::OpConstrainsLhs(
6049 griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA",
6050 feGriffithConstrainsLhs->deltaVec));
6051
6053}
multi_index_container< boost::shared_ptr< Field >, indexed_by< hashed_unique< tag< BitFieldId_mi_tag >, const_mem_fun< Field, const BitFieldId &, &Field::getId >, HashBit< BitFieldId >, EqBit< BitFieldId > >, ordered_unique< tag< Meshset_mi_tag >, member< Field, EntityHandle, &Field::meshSet > >, ordered_unique< tag< FieldName_mi_tag >, const_mem_fun< Field, boost::string_ref, &Field::getNameRef > >, ordered_non_unique< tag< BitFieldId_space_mi_tag >, const_mem_fun< Field, FieldSpace, &Field::getSpace > > > > Field_multiIndex
Field_multiIndex for Field.
@ BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
double smootherGamma
Controls mesh smoothing.
double griffithR
Griffith regularisation parameter.
double smootherAlpha
Controls mesh smoothing.
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > skinOrientation
boost::shared_ptr< VolumeLengthQuality< adouble > > volumeLengthAdouble
double griffithE
Griffith stability parameter.
boost::shared_ptr< VolumeLengthQuality< double > > volumeLengthDouble
boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > crackOrientation
boost::shared_ptr< GriffithForceElement > griffithForceElement
static MoFEMErrorCode getElementOptions(BlockData &block_data)
Volume Length Quality.

◆ 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 1741 of file CrackPropagation.cpp.

1743 {
1744 const RefEntity_multiIndex *refined_ents_ptr;
1745 const FieldEntity_multiIndex *field_ents_ptr;
1747 CHKERR mField.get_ref_ents(&refined_ents_ptr);
1748 CHKERR mField.get_field_ents(&field_ents_ptr);
1749 if (mField.check_field("LAMBDA_ARC_LENGTH")) {
1750 auto vit = refined_ents_ptr->get<Ent_mi_tag>().find(arcLengthVertex);
1751 if (vit == refined_ents_ptr->get<Ent_mi_tag>().end()) {
1752 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1753 "Vertex of arc-length field not found");
1754 }
1755 *(const_cast<RefEntity *>(vit->get())->getBitRefLevelPtr()) |= bit;
1756 auto fit = field_ents_ptr->get<Unique_mi_tag>().find(
1757 FieldEntity::getLocalUniqueIdCalculate(
1758 mField.get_field_bit_number("LAMBDA_ARC_LENGTH"), arcLengthVertex));
1759 if (fit == field_ents_ptr->get<