diff --git a/framework/include/interfaces/XFEMInterface.h b/framework/include/interfaces/XFEMInterface.h index d94bb0a0a3aa..1998211a1f54 100644 --- a/framework/include/interfaces/XFEMInterface.h +++ b/framework/include/interfaces/XFEMInterface.h @@ -93,6 +93,7 @@ class XFEMInterface : public ConsoleStreamInterface const std::vector> & nl, AuxiliarySystem & aux) = 0; + virtual bool didNearTipEnrichmentChange() = 0; /** * Initialize the solution on newly created nodes */ diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index f31d6b51d760..8c057f1e52f7 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -8284,10 +8284,20 @@ FEProblemBase::updateMeshXFEM() restoreSolutions(); _console << "\nXFEM update complete: Mesh modified" << std::endl; } + } + // Changing near-tip enrichment requires repeating the solution, even if the mesh is not updated + bool crack_front_advanced = false; + if (!updated) + { + crack_front_advanced = _xfem->didNearTipEnrichmentChange(); + if (crack_front_advanced) + _console << "\nXFEM update complete: Mesh not modified but advancing crack changed near-tip " + "enrichment" + << std::endl; else _console << "\nXFEM update complete: Mesh not modified" << std::endl; } - return updated; + return (updated || crack_front_advanced); } void diff --git a/modules/solid_mechanics/include/materials/ComputeIncrementalStrainBase.h b/modules/solid_mechanics/include/materials/ComputeIncrementalStrainBase.h index 51d0ec9d1f6f..7c7686260e84 100644 --- a/modules/solid_mechanics/include/materials/ComputeIncrementalStrainBase.h +++ b/modules/solid_mechanics/include/materials/ComputeIncrementalStrainBase.h @@ -40,4 +40,6 @@ class ComputeIncrementalStrainBase : public ComputeStrainBase const MaterialProperty & _total_strain_old; std::vector *> _eigenstrains_old; + /// Rate of change of the displacement gradient + MaterialProperty & _grad_disp_rate; }; diff --git a/modules/solid_mechanics/include/userobjects/CrackFrontDefinition.h b/modules/solid_mechanics/include/userobjects/CrackFrontDefinition.h index 6ccd56088be8..cb0619e97cc9 100644 --- a/modules/solid_mechanics/include/userobjects/CrackFrontDefinition.h +++ b/modules/solid_mechanics/include/userobjects/CrackFrontDefinition.h @@ -41,6 +41,13 @@ class CrackFrontDefinition : public GeneralUserObject, public BoundaryRestrictab virtual void finalize() override; virtual void execute() override; + /** + * Obtain the updated set of crack points from the CrackFrontPointsProvider and + * perform other needed initialization for that set of crack front points. This is only + * to be called if a CrackFrontPointsProvider is used. + */ + void updateCrackFrontPoints(); + /// used by Actions to add CrackFrontDefinitionParams static void includeCrackFrontDefinitionParams(InputParameters & params); diff --git a/modules/solid_mechanics/src/materials/ComputeIncrementalStrainBase.C b/modules/solid_mechanics/src/materials/ComputeIncrementalStrainBase.C index 5e45e4529f2e..947a2362459f 100644 --- a/modules/solid_mechanics/src/materials/ComputeIncrementalStrainBase.C +++ b/modules/solid_mechanics/src/materials/ComputeIncrementalStrainBase.C @@ -26,7 +26,8 @@ ComputeIncrementalStrainBase::ComputeIncrementalStrainBase(const InputParameters _deformation_gradient(declareProperty(_base_name + "deformation_gradient")), _mechanical_strain_old(getMaterialPropertyOld(_base_name + "mechanical_strain")), _total_strain_old(getMaterialPropertyOld(_base_name + "total_strain")), - _eigenstrains_old(_eigenstrain_names.size()) + _eigenstrains_old(_eigenstrain_names.size()), + _grad_disp_rate(declareProperty(_base_name + "grad_disp_rate")) { for (unsigned int i = 0; i < _eigenstrains_old.size(); ++i) _eigenstrains_old[i] = &getMaterialPropertyOld(_eigenstrain_names[i]); diff --git a/modules/solid_mechanics/src/userobjects/CrackFrontDefinition.C b/modules/solid_mechanics/src/userobjects/CrackFrontDefinition.C index 75a6a54c6ba7..9973b0dc881f 100644 --- a/modules/solid_mechanics/src/userobjects/CrackFrontDefinition.C +++ b/modules/solid_mechanics/src/userobjects/CrackFrontDefinition.C @@ -357,25 +357,27 @@ CrackFrontDefinition::initialSetup() } } +void +CrackFrontDefinition::updateCrackFrontPoints() +{ + mooseAssert(_crack_front_points_provider, "_crack_front_points_provider is NULL"); + _crack_front_points = + _crack_front_points_provider->getCrackFrontPoints(_num_points_from_provider); + updateCrackFrontGeometry(); + std::size_t num_crack_front_points = getNumCrackFrontPoints(); + if (_q_function_type == "GEOMETRY") + for (std::size_t i = 0; i < num_crack_front_points; ++i) + { + bool is_point_on_intersecting_boundary = isPointWithIndexOnIntersectingBoundary(i); + _is_point_on_intersecting_boundary.push_back(is_point_on_intersecting_boundary); + } +} + void CrackFrontDefinition::initialize() { - // Update the crack front for fracture integral calculations - // This is only useful for growing cracks which are currently described by the mesh - // cutter if (_use_mesh_cutter && _is_cutter_modified) - { - _crack_front_points = - _crack_front_points_provider->getCrackFrontPoints(_num_points_from_provider); - updateCrackFrontGeometry(); - std::size_t num_crack_front_points = getNumCrackFrontPoints(); - if (_q_function_type == "GEOMETRY") - for (std::size_t i = 0; i < num_crack_front_points; ++i) - { - bool is_point_on_intersecting_boundary = isPointWithIndexOnIntersectingBoundary(i); - _is_point_on_intersecting_boundary.push_back(is_point_on_intersecting_boundary); - } - } + updateCrackFrontPoints(); } void diff --git a/modules/xfem/doc/content/source/materials/ComputeCrackTipEnrichmentIncrementalStrain.md b/modules/xfem/doc/content/source/materials/ComputeCrackTipEnrichmentIncrementalStrain.md new file mode 100644 index 000000000000..9519f49106af --- /dev/null +++ b/modules/xfem/doc/content/source/materials/ComputeCrackTipEnrichmentIncrementalStrain.md @@ -0,0 +1,13 @@ +# ComputeCrackTipEnrichmentIncrementalStrain + +!syntax description /Materials/ComputeCrackTipEnrichmentIncrementalStrain + +# Description + +With XFEM, the displacement field contains both the standard and optionally, the near-tip enrichment solution. This Material calculates the strains for near-tip enrichment from both the standard and enrichment solutions. This model is incremental, which allows for the use of nonlinear material models, but is applicable only to small-strain problems. This is a drop-in replacement for the standard strain calculators. + +!syntax parameters /Materials/ComputeCrackTipEnrichmentIncrementalStrain + +!syntax inputs /Materials/ComputeCrackTipEnrichmentIncrementalStrain + +!syntax children /Materials/ComputeCrackTipEnrichmentIncrementalStrain diff --git a/modules/xfem/include/base/XFEM.h b/modules/xfem/include/base/XFEM.h index 66f2f8e4e0fd..ce2da85052a3 100644 --- a/modules/xfem/include/base/XFEM.h +++ b/modules/xfem/include/base/XFEM.h @@ -323,6 +323,13 @@ class XFEM : public XFEMInterface const Elem * cut_elem, const Elem * parent_elem = nullptr) const; + /** + * Check whether the near-tip enrichment changed during the most recent call to + * XFEM::update(). + * @return true if it changed + */ + virtual bool didNearTipEnrichmentChange() override; + private: bool _has_secondary_cut; diff --git a/modules/xfem/include/materials/ComputeCrackTipEnrichmentIncrementalStrain.h b/modules/xfem/include/materials/ComputeCrackTipEnrichmentIncrementalStrain.h new file mode 100644 index 000000000000..9eb9b8ce8959 --- /dev/null +++ b/modules/xfem/include/materials/ComputeCrackTipEnrichmentIncrementalStrain.h @@ -0,0 +1,70 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once +#include "ComputeIncrementalStrainBase.h" +#include "Material.h" +#include "RankTwoTensor.h" +#include "RankFourTensor.h" +#include "RotationTensor.h" +#include "Assembly.h" +#include "CrackFrontDefinition.h" +#include "EnrichmentFunctionCalculation.h" + +/** + * ComputeIncrementalStrain defines a strain increment and rotation increment (=1), for small + * strains. + */ +class ComputeCrackTipEnrichmentIncrementalStrain : public ComputeIncrementalStrainBase, + public EnrichmentFunctionCalculation +{ +public: + static InputParameters validParams(); + + ComputeCrackTipEnrichmentIncrementalStrain(const InputParameters & parameters); + virtual ~ComputeCrackTipEnrichmentIncrementalStrain() {} + virtual void computeProperties() override; + +protected: + /// enrichment displacement + std::vector _enrich_disp; + + /// gradient of enrichment displacement + std::vector _grad_enrich_disp; + std::vector _grad_enrich_disp_old; + + /// enrichment displacement variables + std::vector> _enrich_variable; + + /// the current shape functions + const VariablePhiValue & _phi; + + /// gradient of the shape function + const VariablePhiGradient & _grad_phi; + + /// gradient of the enriched solution + MaterialProperty & _grad_enrich_disp_tensor; + const MaterialProperty & _grad_enrich_disp_tensor_old; + +private: + /// enrichment function value + std::vector _B; + /// derivatives of enrichment function respect to global cooridnate + std::vector _dBX; + /// derivatives of enrichment function respect to crack front cooridnate + std::vector _dBx; + /// enrichment function at node I + std::vector> _BI; + /// shape function + const std::vector> * _fe_phi; + /// gradient of shape function + const std::vector> * _fe_dphi; + NonlinearSystem * _nl; + const NumericVector * _sln; +}; diff --git a/modules/xfem/include/userobjects/GeometricCutUserObject.h b/modules/xfem/include/userobjects/GeometricCutUserObject.h index 412411912d79..e3995f81360a 100644 --- a/modules/xfem/include/userobjects/GeometricCutUserObject.h +++ b/modules/xfem/include/userobjects/GeometricCutUserObject.h @@ -195,6 +195,8 @@ class GeometricCutUserObject : public CrackFrontPointsProvider */ CutSubdomainID getCutSubdomainID(const Elem * elem) const; + virtual bool isCutterMeshChanged() const; + protected: /// Pointer to the XFEM controller object std::shared_ptr _xfem; diff --git a/modules/xfem/include/userobjects/MeshCut2DFractureUserObject.h b/modules/xfem/include/userobjects/MeshCut2DFractureUserObject.h index 80d2b995d098..9ee6a163b472 100644 --- a/modules/xfem/include/userobjects/MeshCut2DFractureUserObject.h +++ b/modules/xfem/include/userobjects/MeshCut2DFractureUserObject.h @@ -27,6 +27,8 @@ class MeshCut2DFractureUserObject : public MeshCut2DUserObjectBase virtual void initialize() override; + virtual bool isCutterMeshChanged() const override; + protected: virtual void findActiveBoundaryGrowth() override; diff --git a/modules/xfem/include/userobjects/MeshCut2DUserObjectBase.h b/modules/xfem/include/userobjects/MeshCut2DUserObjectBase.h index 55e58aa92e0f..a404584aff18 100644 --- a/modules/xfem/include/userobjects/MeshCut2DUserObjectBase.h +++ b/modules/xfem/include/userobjects/MeshCut2DUserObjectBase.h @@ -49,6 +49,8 @@ class MeshCut2DUserObjectBase : public MeshCutUserObjectBase virtual const std::vector getCrackPlaneNormals(unsigned int num_crack_front_points) const override; + bool isPointOnEdgeBoundary(const Point & point, Real tolerance = 1e-10); + protected: /// The FE solution mesh MooseMesh & _mesh; diff --git a/modules/xfem/src/base/XFEM.C b/modules/xfem/src/base/XFEM.C index e9c3612c29c1..ad7ac3f70a7d 100644 --- a/modules/xfem/src/base/XFEM.C +++ b/modules/xfem/src/base/XFEM.C @@ -183,6 +183,19 @@ XFEM::clearGeomMarkedElems() _geom_marked_elems_3d.clear(); } +bool +XFEM::didNearTipEnrichmentChange() +{ + bool cutter_mesh_changed = false; + for (unsigned int i = 0; i < _geometric_cuts.size(); ++i) + { + cutter_mesh_changed = _geometric_cuts[i]->isCutterMeshChanged(); + if (cutter_mesh_changed) + break; + } + return cutter_mesh_changed; +} + void XFEM::storeCrackTipOriginAndDirection() { @@ -614,7 +627,6 @@ XFEM::markCutEdgesByState(Real time) // crack tip origin coordinates and direction Point crack_tip_origin(0, 0, 0); Point crack_tip_direction(0, 0, 0); - if (isElemAtCrackTip(elem)) // crack tip element's crack intiation { orig_cut_side_id = CEMElem->getTipEdgeID(); diff --git a/modules/xfem/src/materials/ComputeCrackTipEnrichmentIncrementalStrain.C b/modules/xfem/src/materials/ComputeCrackTipEnrichmentIncrementalStrain.C new file mode 100644 index 000000000000..f0c34e6372ce --- /dev/null +++ b/modules/xfem/src/materials/ComputeCrackTipEnrichmentIncrementalStrain.C @@ -0,0 +1,172 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "ComputeCrackTipEnrichmentIncrementalStrain.h" +#include "MooseMesh.h" +#include "libmesh/fe_interface.h" +#include "libmesh/string_to_enum.h" + +#include "libmesh/quadrature.h" + +registerMooseObject("XFEMApp", ComputeCrackTipEnrichmentIncrementalStrain); + +InputParameters +ComputeCrackTipEnrichmentIncrementalStrain::validParams() +{ + InputParameters params = ComputeIncrementalStrainBase::validParams(); + params.addClassDescription( + "Computes the crack tip enrichment strain for an incremental small-strain formulation."); + params.addRequiredParam>("enrichment_displacements", + "The enrichment displacement"); + + params.addRequiredParam("crack_front_definition", + "The CrackFrontDefinition user object name"); + + return params; +} + +ComputeCrackTipEnrichmentIncrementalStrain::ComputeCrackTipEnrichmentIncrementalStrain( + const InputParameters & parameters) + : ComputeIncrementalStrainBase(parameters), + EnrichmentFunctionCalculation(&getUserObject("crack_front_definition")), + _enrich_disp(3), + _grad_enrich_disp(3), + _grad_enrich_disp_old(3), + _enrich_variable(4), + _phi(_assembly.phi()), + _grad_phi(_assembly.gradPhi()), + _grad_enrich_disp_tensor( + declareProperty(_base_name + "grad_enrich_disp_tensor")), + _grad_enrich_disp_tensor_old( + getMaterialPropertyOld(_base_name + "grad_enrich_disp_tensor")), + _B(4), + _dBX(4), + _dBx(4) +{ + for (unsigned int i = 0; i < _enrich_variable.size(); ++i) + _enrich_variable[i].resize(_ndisp); + + const std::vector & nl_vnames = + getParam>("enrichment_displacements"); + + if (_ndisp == 2 && nl_vnames.size() != 8) + mooseError("The number of enrichment displacements should be total 8 for 2D."); + else if (_ndisp == 3 && nl_vnames.size() != 12) + mooseError("The number of enrichment displacements should be total 12 for 3D."); + + _nl = &(_fe_problem.getNonlinearSystem(/*nl_sys_num=*/0)); + + for (unsigned int j = 0; j < _ndisp; ++j) + for (unsigned int i = 0; i < 4; ++i) + _enrich_variable[i][j] = &(_nl->getVariable(0, nl_vnames[j * 4 + i])); + + if (_ndisp == 2) + _BI.resize(4); // QUAD4 + else if (_ndisp == 3) + _BI.resize(8); // HEX8 + + for (unsigned int i = 0; i < _BI.size(); ++i) + _BI[i].resize(4); +} + +void +ComputeCrackTipEnrichmentIncrementalStrain::computeProperties() +{ + FEType fe_type(Utility::string_to_enum("first"), + Utility::string_to_enum("lagrange")); + const unsigned int dim = _current_elem->dim(); + std::unique_ptr fe(FEBase::build(dim, fe_type)); + fe->attach_quadrature_rule(const_cast(_qrule)); + _fe_phi = &(fe->get_phi()); + _fe_dphi = &(fe->get_dphi()); + + if (isBoundaryMaterial()) + fe->reinit(_current_elem, _current_side); + else + fe->reinit(_current_elem); + + _sln = _nl->currentSolution(); + + for (unsigned int i = 0; i < _BI.size(); ++i) + crackTipEnrichementFunctionAtPoint(*(_current_elem->node_ptr(i)), _BI[i]); + + for (_qp = 0; _qp < _qrule->n_points(); ++_qp) + { + // enrichment function + crackTipEnrichementFunctionAtPoint(_q_point[_qp], _B); + unsigned int crack_front_point_index = + crackTipEnrichementFunctionDerivativeAtPoint(_q_point[_qp], _dBx); + + for (unsigned int i = 0; i < 4; ++i) + rotateFromCrackFrontCoordsToGlobal(_dBx[i], _dBX[i], crack_front_point_index); + + for (unsigned int m = 0; m < _ndisp; ++m) + { + _enrich_disp[m] = 0.0; + _grad_enrich_disp[m].zero(); + _grad_enrich_disp_old[m].zero(); + for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i) + { + const Node * node_i = _current_elem->node_ptr(i); + for (unsigned int j = 0; j < 4; ++j) + { + dof_id_type dof = node_i->dof_number(_nl->number(), _enrich_variable[j][m]->number(), 0); + Real soln = (*_sln)(dof); + _enrich_disp[m] += (*_fe_phi)[i][_qp] * (_B[j] - _BI[i][j]) * soln; + RealVectorValue grad_B(_dBX[j]); + + _grad_enrich_disp[m] += + ((*_fe_dphi)[i][_qp] * (_B[j] - _BI[i][j]) + (*_fe_phi)[i][_qp] * grad_B) * soln; + } + } + } + + _grad_enrich_disp_tensor[_qp] = RankTwoTensor::initializeFromRows( + _grad_enrich_disp[0], _grad_enrich_disp[1], _grad_enrich_disp[2]); + + auto grad_disp_tensor = RankTwoTensor::initializeFromRows( + (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]); + + auto grad_disp_tensor_old = RankTwoTensor::initializeFromRows( + (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]); + + _deformation_gradient[_qp] = grad_disp_tensor + _grad_enrich_disp_tensor[_qp]; + _deformation_gradient[_qp].addIa(1.0); + + _strain_increment[_qp] = + 0.5 * ((grad_disp_tensor + _grad_enrich_disp_tensor[_qp]) + + (grad_disp_tensor + _grad_enrich_disp_tensor[_qp]).transpose()) - + 0.5 * ((grad_disp_tensor_old + _grad_enrich_disp_tensor_old[_qp]) + + (grad_disp_tensor_old + _grad_enrich_disp_tensor_old[_qp]).transpose()); + + // total strain + _total_strain[_qp] = _total_strain_old[_qp] + _strain_increment[_qp]; + + // Remove the Eigen strain increment + subtractEigenstrainIncrementFromStrain(_strain_increment[_qp]); + + // strain rate + if (_dt > 0) + { + _strain_rate[_qp] = _strain_increment[_qp] / _dt; + _grad_disp_rate[_qp] = (grad_disp_tensor - grad_disp_tensor_old) / _dt; + } + else + { + _strain_rate[_qp].zero(); + _grad_disp_rate[_qp].zero(); + } + + // Update strain in intermediate configuration: rotations are not needed + _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp]; + + // incremental small strain does not include rotation + _rotation_increment[_qp].setToIdentity(); + } +} diff --git a/modules/xfem/src/userobjects/GeometricCutUserObject.C b/modules/xfem/src/userobjects/GeometricCutUserObject.C index 7175915c133d..7500783f2395 100644 --- a/modules/xfem/src/userobjects/GeometricCutUserObject.C +++ b/modules/xfem/src/userobjects/GeometricCutUserObject.C @@ -126,6 +126,12 @@ GeometricCutUserObject::execute() } } +bool +GeometricCutUserObject::isCutterMeshChanged() const +{ + return false; +} + void GeometricCutUserObject::threadJoin(const UserObject & y) { diff --git a/modules/xfem/src/userobjects/MeshCut2DFractureUserObject.C b/modules/xfem/src/userobjects/MeshCut2DFractureUserObject.C index 32097fb99fe0..66467f667382 100644 --- a/modules/xfem/src/userobjects/MeshCut2DFractureUserObject.C +++ b/modules/xfem/src/userobjects/MeshCut2DFractureUserObject.C @@ -100,6 +100,14 @@ MeshCut2DFractureUserObject::initialize() _crack_front_definition->updateNumberOfCrackFrontPoints( _original_and_current_front_node_ids.size()); _crack_front_definition->isCutterModified(_is_mesh_modified); + if (_is_mesh_modified) + _crack_front_definition->updateCrackFrontPoints(); +} + +bool +MeshCut2DFractureUserObject::isCutterMeshChanged() const +{ + return _is_mesh_modified; } void diff --git a/modules/xfem/src/userobjects/MeshCut2DUserObjectBase.C b/modules/xfem/src/userobjects/MeshCut2DUserObjectBase.C index 9a57763017dd..fef6c0181d31 100644 --- a/modules/xfem/src/userobjects/MeshCut2DUserObjectBase.C +++ b/modules/xfem/src/userobjects/MeshCut2DUserObjectBase.C @@ -279,14 +279,61 @@ MeshCut2DUserObjectBase::findOriginalCrackFrontNodes() mooseAssert(this_node, "Node is NULL"); Point & this_point = *this_node; + bool is_on_bc = isPointOnEdgeBoundary(this_point); + const Elem * elem = (*pl)(this_point); - if (elem != NULL) + if (elem != NULL && !is_on_bc) + { _original_and_current_front_node_ids.push_back(std::make_pair(node, node)); + } } std::sort(_original_and_current_front_node_ids.begin(), _original_and_current_front_node_ids.end()); } +bool +MeshCut2DUserObjectBase::isPointOnEdgeBoundary(const Point & point, Real tolerance) +{ + for (const auto & bnd_elem : as_range(_mesh.bndElemsBegin(), _mesh.bndElemsEnd())) + { + // if (bnd_elem->_elem->processor_id() == processor_id()) + { + auto side = bnd_elem->_elem->side_ptr(bnd_elem->_side); + + if (side->n_nodes() == 2) + { + // Create vectors for the endpoints of the edge + Point p1 = side->point(0); + Point p2 = side->point(1); + + // Compute the vector from p1 to p2 (edge vector) and from p1 to the point + Point edge_vector = p2 - p1; + Point point_vector = point - p1; + + // Compute the cross product magnitude to find if the point is collinear with the edge + Real cross_product_magnitude = std::sqrt( + std::pow(edge_vector(1) * point_vector(2) - edge_vector(2) * point_vector(1), 2) + + std::pow(edge_vector(2) * point_vector(0) - edge_vector(0) * point_vector(2), 2) + + std::pow(edge_vector(0) * point_vector(1) - edge_vector(1) * point_vector(0), 2)); + + // Check if the point lies on the edge within a tolerance + if (cross_product_magnitude < tolerance) + { + // Check if the point lies within the edge segment bounds + Real dot_product = point_vector * edge_vector; + Real edge_length_squared = edge_vector * edge_vector; + + if (dot_product >= 0 && dot_product <= edge_length_squared) + return true; + } + } + else + mooseError("MeshCut2DUserObjectBase currently only supports linear elements.\n"); + } + } + return false; +} + void MeshCut2DUserObjectBase::growFront() { diff --git a/modules/xfem/test/tests/crack_tip_test/elastic_test.i b/modules/xfem/test/tests/crack_tip_test/elastic_test.i new file mode 100644 index 000000000000..5f5fca60aa90 --- /dev/null +++ b/modules/xfem/test/tests/crack_tip_test/elastic_test.i @@ -0,0 +1,207 @@ +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = false +[] + +[XFEM] + geometric_cut_userobjects = 'cut_mesh2' + qrule = moment_fitting + output_cut_plane = true + use_crack_tip_enrichment = true + crack_front_definition = crackFrontDefinition + enrichment_displacements = 'enrich1_x enrich2_x enrich3_x enrich4_x enrich1_y enrich2_y enrich3_y enrich4_y' + displacements = 'disp_x disp_y' + cut_off_boundary = all + cut_off_radius = 0.3726 +[] + +[Mesh] + file = square3_3.e +[] + +[DomainIntegral] + integrals = 'JIntegral InteractionIntegralKI InteractionIntegralKII' + displacements = 'disp_x disp_y' + crack_front_points_provider = cut_mesh2 + 2d = true + number_points_from_provider = 1 + crack_direction_method = CurvedCrackFront + radius_inner = '0.3' + radius_outer = '0.5' + youngs_modulus = 200000 + poissons_ratio = 0.3 + output_q = true + output_vpp = false + incremental = true +[] + +[UserObjects] + [cut_mesh2] + type = MeshCut2DFractureUserObject + mesh_file = line0p45.e + growth_increment = 0.165 + ki_vectorpostprocessor = "II_KI_1" + k_critical = 0 + [] +[] + +[Variables] + [disp_x] + order = FIRST + family = LAGRANGE + [] + [disp_y] + order = FIRST + family = LAGRANGE + [] +[] + +[AuxVariables] + [saved_x] + [] + [saved_y] + [] + [stress_xx] + order = CONSTANT + family = MONOMIAL + [] + [stress_yy] + order = CONSTANT + family = MONOMIAL + [] + [stress_xy] + order = CONSTANT + family = MONOMIAL + [] + [vonmises] + order = CONSTANT + family = MONOMIAL + [] +[] + +[Kernels] + [TensorMechanics] + use_displaced_mesh = false + displacements = 'disp_x disp_y' + add_variables = true + incremental = true + generate_output = 'stress_xx stress_yy' + [] +[] + +[AuxKernels] + [stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + execute_on = timestep_end + [] + [stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + execute_on = timestep_end + [] + [stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + execute_on = timestep_end + [] + [vonmises] + type = RankTwoScalarAux + rank_two_tensor = stress + variable = vonmises + scalar_type = vonmisesStress + execute_on = timestep_end + [] +[] + +[Functions] + [rampConstant] + type = PiecewiseLinear + x = '0. 1e-12 1.0' + y = '1.0 1.0 1.0' + scale_factor = -1e1 + [] +[] + +[BCs] + [fix] + type = DirichletBC + boundary = 'top bottom fix' + variable = disp_x + value = 0.0 + [] + [fixy] + type = DirichletBC + boundary = 'fix' + variable = disp_y + value = 0.0 + [] + [Pressure] + [Side1] + boundary = 'top' + function = rampConstant + [] + [Side2] + boundary = 'bottom' + function = rampConstant + [] + [] +[] + +[Materials] + [elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 200000 + poissons_ratio = 0.3 + [] + [strain] + type = ComputeCrackTipEnrichmentIncrementalStrain + displacements = 'disp_x disp_y' + crack_front_definition = crackFrontDefinition + enrichment_displacements = 'enrich1_x enrich2_x enrich3_x enrich4_x enrich1_y enrich2_y enrich3_y enrich4_y' + [] + [stress] + type = ComputeFiniteStrainElasticStress + [] +[] + +[Executioner] + type = Transient + + solve_type = 'Newton' + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + automatic_scaling = true + [Quadrature] + type = GAUSS + order = SECOND + [] + + l_max_its = 50 + nl_max_its = 100 + + nl_abs_tol = 1e-5 + nl_rel_tol = 1e-5 + + start_time = 0.0 + dt = 0.1 + end_time = 0.3 + max_xfem_update = 1 +[] + +[Preconditioning] + [smp] + type = SMP + full = true + [] +[] diff --git a/modules/xfem/test/tests/crack_tip_test/gold/output_JIntegral_out.csv b/modules/xfem/test/tests/crack_tip_test/gold/output_JIntegral_out.csv new file mode 100644 index 000000000000..a8369d028f5b --- /dev/null +++ b/modules/xfem/test/tests/crack_tip_test/gold/output_JIntegral_out.csv @@ -0,0 +1,4 @@ +time,II_KII_1_1,II_KI_1_1,J_1_1 +0.1,-2.6724828472357e-15,15.309312344397,0.0011232903683779 +0.2,-1.7276278305883e-14,32.315814385054,0.0061163481112729 +0.3,5.1824461981461e-16,24.811530197213,0.0051774888774892 diff --git a/modules/xfem/test/tests/crack_tip_test/gold/output_cutter_mesh_xfemcutter_XFEMCutMeshOutput.e-s0003 b/modules/xfem/test/tests/crack_tip_test/gold/output_cutter_mesh_xfemcutter_XFEMCutMeshOutput.e-s0003 new file mode 100644 index 000000000000..c9c3aad33bdb Binary files /dev/null and b/modules/xfem/test/tests/crack_tip_test/gold/output_cutter_mesh_xfemcutter_XFEMCutMeshOutput.e-s0003 differ diff --git a/modules/xfem/test/tests/crack_tip_test/line0p45.e b/modules/xfem/test/tests/crack_tip_test/line0p45.e new file mode 100644 index 000000000000..6c1a23a7daa1 Binary files /dev/null and b/modules/xfem/test/tests/crack_tip_test/line0p45.e differ diff --git a/modules/xfem/test/tests/crack_tip_test/output_JIntegral.i b/modules/xfem/test/tests/crack_tip_test/output_JIntegral.i new file mode 100644 index 000000000000..aee4e4ad3d17 --- /dev/null +++ b/modules/xfem/test/tests/crack_tip_test/output_JIntegral.i @@ -0,0 +1,9 @@ +[Outputs] + #exodus = true + csv = true + execute_on = timestep_end + [./console] + type = Console + output_linear = true + [../] +[] diff --git a/modules/xfem/test/tests/crack_tip_test/output_cutter_mesh.i b/modules/xfem/test/tests/crack_tip_test/output_cutter_mesh.i new file mode 100644 index 000000000000..2d9e121461c5 --- /dev/null +++ b/modules/xfem/test/tests/crack_tip_test/output_cutter_mesh.i @@ -0,0 +1,11 @@ +[Outputs] + [xfemcutter] + type = XFEMCutMeshOutput + xfem_cutter_uo = cut_mesh2 + [] + execute_on = timestep_end + [./console] + type = Console + output_linear = true + [../] +[] diff --git a/modules/xfem/test/tests/crack_tip_test/square3_3.e b/modules/xfem/test/tests/crack_tip_test/square3_3.e new file mode 100644 index 000000000000..b99ac36ee570 Binary files /dev/null and b/modules/xfem/test/tests/crack_tip_test/square3_3.e differ diff --git a/modules/xfem/test/tests/crack_tip_test/tests b/modules/xfem/test/tests/crack_tip_test/tests new file mode 100644 index 000000000000..229b488ac110 --- /dev/null +++ b/modules/xfem/test/tests/crack_tip_test/tests @@ -0,0 +1,22 @@ +[Tests] + issues = '#31729' + design = 'ComputeCrackTipEnrichmentIncrementalStrain.md' + requirement = 'The XFEM module shall represent a propagating crack in a 2D mechanics model with XFEM, ' + 'where the crack tip is arbitrarily inside elements' + [tip_location] + type = Exodiff + input = 'elastic_test.i output_cutter_mesh.i' + exodiff = 'output_cutter_mesh_xfemcutter_XFEMCutMeshOutput.e-s0003' + abs_zero = 1e-8 + map = false + capabilities = 'unique_id' + requirement = 'test crack tip location' + [] + [JIntegral] + type = CSVDiff + input = 'elastic_test.i output_JIntegral.i' + csvdiff = 'output_JIntegral_out.csv' + capabilities = 'unique_id' + requirement = 'test J integral' + [] +[]