v0.14.0
Loading...
Searching...
No Matches
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...
 
PetscBool isSurfaceForceAle
 If true surface pressure is considered in ALE. More...
 
PetscBool ignoreMaterialForce
 If true surface pressure is considered in ALE. More...
 
std::string mwlsApproxFile
 Name of file with internal stresses. More...
 
std::string mwlsStressTagName
 Name of tag with internal stresses. More...
 
std::string mwlsEigenStressTagName
 Name of tag with eigen stresses. More...
 
std::string mwlsRhoTagName
 Name of tag with density. More...
 
int residualStressBlock
 Block on which residual stress is applied. More...
 
int densityMapBlock
 
std::string configFile
 
bool isStressTagSavedOnNodes
 
double partitioningWeightPower
 
PetscBool resetMWLSCoeffsEveryPropagationStep
 
PetscBool addAnalyticalInternalStressOperators
 
PetscBool solveEigenStressProblem
 Solve eigen problem. More...
 
PetscBool useEigenPositionsSimpleContact
 
int postProcLevel
 level of postprocessing (amount of output files) More...
 
int startStep
 
double crackFrontLength
 
double crackSurfaceArea
 
std::map< std::string, BitRefLevel > mapBitLevel
 

Static Public Attributes

static bool debug = false
 
static bool parallelReadAndBroadcast
 

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 > > surfaceForceAle
 assemble surface pressure (ALE) More...
 
boost::shared_ptr< NeumannForcesSurface::DataAtIntegrationPtscommonDataSurfaceForceAle
 common data at integration points (ALE) More...
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressure
 assemble surface pressure More...
 
boost::shared_ptr< boost::ptr_map< string, NeumannForcesSurface > > surfacePressureAle
 assemble surface pressure (ALE) More...
 
boost::shared_ptr< NeumannForcesSurface::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< FaceElementForcesAndSourcesCorefeLhsSurfaceForceALEMaterial
 
boost::shared_ptr< FaceElementForcesAndSourcesCorefeRhsSurfaceForceALEMaterial
 
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 declareSurfaceForceAleFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
 Declare FE for pressure BC in ALE formulation (in material domain) More...
 
MoFEMErrorCode declarePressureAleFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
 Declare FE for pressure BC in ALE formulation (in material domain) More...
 
MoFEMErrorCode declareSpringsAleFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
 Declare FE for spring BC in ALE formulation (in material domain) More...
 
MoFEMErrorCode declareSimpleContactAleFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true)
 Declare FE for pressure BC in ALE formulation (in material domain) More...
 
MoFEMErrorCode declareMaterialFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 declare material finite elements More...
 
MoFEMErrorCode declareFrontFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 
MoFEMErrorCode declareBothSidesFE (const BitRefLevel bit_spatial, const BitRefLevel bit_material, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 
MoFEMErrorCode declareSmoothingFE (const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set(), const bool proc_only=true, const bool verb=QUIET)
 declare mesh smoothing finite elements More...
 
MoFEMErrorCode declareSurfaceFE (std::string fe_name, const BitRefLevel bit, const BitRefLevel mask, const std::vector< int > &ids, const bool proc_only=true, const int verb=QUIET, const bool debug=false)
 declare surface sliding elements More...
 
MoFEMErrorCode declareEdgeFE (std::string fe_name, const BitRefLevel bit, const BitRefLevel mask, const bool proc_only=true, const int verb=QUIET, const bool debug=false)
 
MoFEMErrorCode declareCrackSurfaceFE (std::string fe_name, const BitRefLevel bit, const BitRefLevel mask, const bool proc_only=true, const int verb=QUIET, const bool debug=false)
 

Create DMs & problems

MoFEMErrorCode createCrackPropagationDM (SmartPetscObj< DM > &dm, const std::string prb_name, SmartPetscObj< DM > dm_elastic, SmartPetscObj< DM > dm_material, const BitRefLevel bit, const BitRefLevel mask, const std::vector< std::string > fe_list)
 Create DM by composition of elastic DM and material DM. More...
 
MoFEMErrorCode createElasticDM (SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
 Create elastic problem DM. More...
 
MoFEMErrorCode createEigenElasticDM (SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
 Create elastic problem DM. More...
 
MoFEMErrorCode createBcDM (SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask=BitRefLevel().set())
 Create problem to calculate boundary conditions. More...
 
MoFEMErrorCode createMaterialDM (SmartPetscObj< DM > &dm, const std::string prb_name, const BitRefLevel bit, const BitRefLevel mask, const std::vector< std::string > fe_list, const bool debug=false)
 Create DM fto calculate material problem. More...
 
MoFEMErrorCode createMaterialForcesDM (SmartPetscObj< DM > &dm, SmartPetscObj< DM > dm_material, const std::string prb_name, const int verb=QUIET)
 Create DM for calculation of material forces (sub DM of DM material) More...
 
MoFEMErrorCode createSurfaceProjectionDM (SmartPetscObj< DM > &dm, SmartPetscObj< DM > dm_material, const std::string prb_name, const std::vector< int > surface_ids, const std::vector< std::string > fe_list, const int verb=QUIET)
 create DM to calculate projection matrices (sub DM of DM material) More...
 
MoFEMErrorCode createCrackFrontAreaDM (SmartPetscObj< DM > &dm, SmartPetscObj< DM > dm_material, const std::string prb_name, const bool verb=QUIET)
 create DM to calculate Griffith energy More...
 

Create finite element instances

MoFEMErrorCode assembleElasticDM (const std::string mwls_stress_tag_name, const int verb=QUIET, const bool debug=false)
 create elastic finite element instance for spatial assembly More...
 
MoFEMErrorCode assembleElasticDM (const int verb=QUIET, const bool debug=false)
 create elastic finite element instance for spatial assembly More...
 
MoFEMErrorCode addElasticFEInstancesToSnes (DM dm, Mat m, Vec q, Vec f, boost::shared_ptr< FEMethod > arc_method=boost::shared_ptr< FEMethod >(), boost::shared_ptr< ArcLengthCtx > arc_ctx=nullptr, const int verb=QUIET, const bool debug=false)
 
MoFEMErrorCode assembleMaterialForcesDM (DM dm, const int verb=QUIET, const bool debug=false)
 create material element instance More...
 
MoFEMErrorCode assembleSmootherForcesDM (DM dm, const std::vector< int > ids, const int verb=QUIET, const bool debug=false)
 create smoothing element instance More...
 
MoFEMErrorCode assembleCouplingForcesDM (DM dm, const int verb=QUIET, const bool debug=false)
 assemble coupling element instances More...
 
MoFEMErrorCode addMaterialFEInstancesToSnes (DM dm, const bool fix_crack_front, const int verb=QUIET, const bool debug=false)
 add material elements instances to SNES More...
 
MoFEMErrorCode addSmoothingFEInstancesToSnes (DM dm, const bool fix_crack_front, const int verb=QUIET, const bool debug=false)
 add softening elements instances to SNES More...
 
MoFEMErrorCode addPropagationFEInstancesToSnes (DM dm, boost::shared_ptr< FEMethod > arc_method, boost::shared_ptr< ArcLengthCtx > arc_ctx, const std::vector< int > &surface_ids, const int verb=QUIET, const bool debug=false)
 add finite element to SNES for crack propagation problem More...
 
MoFEMErrorCode testJacobians (const BitRefLevel bit, const BitRefLevel mask, tangent_tests test)
 test LHS Jacobians More...
 
MoFEMErrorCode addMWLSStressOperators (boost::shared_ptr< CrackFrontElement > &fe_rhs, boost::shared_ptr< CrackFrontElement > &fe_lhs)
 
MoFEMErrorCode addMWLSDensityOperators (boost::shared_ptr< CrackFrontElement > &fe_rhs, boost::shared_ptr< CrackFrontElement > &fe_lhs)
 
MoFEMErrorCode updateMaterialFixedNode (const bool fix_front, const bool fix_small_g, const bool debug=false)
 Update fixed nodes. More...
 

Assemble and calculate

MoFEMErrorCode calculateElasticEnergy (DM dm, const std::string msg="")
 assemble elastic part of matrix, by running elastic finite element instance More...
 
MoFEMErrorCode calculateMaterialForcesDM (DM dm, Vec q, Vec f, const int verb=QUIET, const bool debug=false)
 assemble material forces, by running material finite element instance More...
 
MoFEMErrorCode calculateSmoothingForcesDM (DM dm, Vec q, Vec f, const int verb=QUIET, const bool debug=false)
 assemble smoothing forces, by running material finite element instance More...
 
MoFEMErrorCode calculateSmoothingForceFactor (const int verb=QUIET, const bool debug=true)
 
MoFEMErrorCode calculateSurfaceProjectionMatrix (DM dm_front, DM dm_project, const std::vector< int > &ids, const int verb=QUIET, const bool debug=false)
 assemble projection matrices More...
 
MoFEMErrorCode calculateFrontProjectionMatrix (DM dm_surface, DM dm_project, const int verb=QUIET, const bool debug=false)
 assemble crack front projection matrix (that constrains crack area growth) More...
 
MoFEMErrorCode projectMaterialForcesDM (DM dm_project, Vec f, Vec f_proj, const int verb=QUIET, const bool debug=false)
 project material forces along the crack elongation direction More...
 
MoFEMErrorCode calculateGriffithForce (DM dm, const double gc, Vec f_griffith, const int verb=QUIET, const bool debug=false)
 calculate Griffith (driving) force More...
 
MoFEMErrorCode projectGriffithForce (DM dm, Vec f_griffith, Vec f_griffith_proj, const int verb=QUIET, const bool debug=false)
 project Griffith forces More...
 
MoFEMErrorCode calculateReleaseEnergy (DM dm, Vec f_material_proj, Vec f_griffith_proj, Vec f_lambda, const double gc, const int verb=QUIET, const bool debug=true)
 calculate release energy More...
 

Solvers

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

Postprocess results

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

Set data from mesh to vectors and back

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

Detailed Description

Crack propagation data and methods

Definition at line 77 of file CrackPropagation.hpp.

Constructor & Destructor Documentation

◆ CrackPropagation()

FractureMechanics::CrackPropagation::CrackPropagation ( MoFEM::Interface m_field,
const int  approx_order = 2,
const int  geometry_order = 1 
)

Constructor of Crack Propagation, prints fracture module version and registers the three interfaces used to solve fracture problem, i.e. CrackPropagation, CPSolvers and CPMeshCut.

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

Definition at line 521 of file CrackPropagation.cpp.

526 geometryOrder(geometry_order), gC(0), betaGc(0),
527 isGcHeterogeneous(PETSC_FALSE), isPartitioned(PETSC_FALSE),
528 propagateCrack(PETSC_FALSE), refAtCrackTip(0), refOrderAtTip(0),
532 nBone(0), cpSolversPtr(new CPSolvers(*this)),
533 cpMeshCutPtr(new CPMeshCut(*this)) {
534
535 if (!LogManager::checkIfChannelExist("CPWorld")) {
536 auto core_log = logging::core::get();
537
538 core_log->add_sink(
540 core_log->add_sink(
542 core_log->add_sink(
544
545 LogManager::setLog("CPWorld");
546 LogManager::setLog("CPSync");
547 LogManager::setLog("CPSelf");
548
549 MOFEM_LOG_TAG("CPWorld", "CP");
550 MOFEM_LOG_TAG("CPSync", "CP");
551 MOFEM_LOG_TAG("CPSelf", "CP");
552 }
553
554 MOFEM_LOG("CPWorld", Sev::noisy) << "CPSolve created";
555
556#ifdef GIT_UM_SHA1_NAME
557 MOFEM_LOG_C("CPWorld", Sev::inform, "User module git commit id %s",
558 GIT_UM_SHA1_NAME);
559#endif
560
561 Version version;
562 getInterfaceVersion(version);
563 MOFEM_LOG("CPWorld", Sev::inform)
564 << "Fracture module version " << version.strVersion();
565
566#ifdef GIT_FM_SHA1_NAME
567 MOFEM_LOG_C("CPWorld", Sev::inform, "Fracture module git commit id %s",
568 GIT_FM_SHA1_NAME);
569#endif
570
571 ierr = registerInterface<CrackPropagation>();
572 CHKERRABORT(PETSC_COMM_SELF, ierr);
573 ierr = registerInterface<CPSolvers>();
574 CHKERRABORT(PETSC_COMM_SELF, ierr);
575 ierr = registerInterface<CPMeshCut>();
576 CHKERRABORT(PETSC_COMM_SELF, ierr);
577
578 mapBitLevel["mesh_cut"] = BitRefLevel().set(0);
579 mapBitLevel["spatial_domain"] = BitRefLevel().set(1);
580 mapBitLevel["material_domain"] = BitRefLevel().set(2);
581
582 if (!moabCommWorld)
583 moabCommWorld = boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
584}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
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.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
Definition: LogManager.cpp:404
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:340
std::string strVersion()

◆ ~CrackPropagation()

FractureMechanics::CrackPropagation::~CrackPropagation ( )
virtual

Definition at line 586 of file CrackPropagation.cpp.

586{}

Member Function Documentation

◆ addElasticFEInstancesToSnes()

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

Definition at line 4610 of file CrackPropagation.cpp.

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

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

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

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

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

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

◆ analyticalStrainFunction()

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

Definition at line 10030 of file CrackPropagation.cpp.

10034 {
10036 constexpr double alpha = 1.e-5;
10040 // FIXME put here formula from test
10041 double temp = 250.;
10042 double z = t_coords(2);
10043 if ((-10. < z && z < -1.) || std::abs(z + 1.) < 1e-15) {
10044 temp = 10. / 3. * (35. - 4. * z);
10045 }
10046 if ((-1. < z && z < 2.) || std::abs(z - 2.) < 1e-15) {
10047 temp = 10. / 3. * (34. - 5. * z);
10048 }
10049 if (2. < z && z < 10.) {
10050 temp = 5. / 4. * (30. + 17. * z);
10051 }
10052
10053 t_thermal_strain(i, j) = alpha * (temp - 250.) * t_kd(i, j);
10054 return t_thermal_strain;
10055}
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 6178 of file CrackPropagation.cpp.

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

◆ assembleElasticDM() [1/2]

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

create elastic finite element instance for spatial assembly

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

Definition at line 4601 of file CrackPropagation.cpp.

4602 {
4604
4606
4608}
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

Rhs

Lhs

Definition at line 3795 of file CrackPropagation.cpp.

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

Operators to assemble rhs

Definition at line 5531 of file CrackPropagation.cpp.

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

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