diff --git a/examples/python/cohesive/custom_material.py b/examples/python/cohesive/custom_material.py index 04edfa027..04b943051 100755 --- a/examples/python/cohesive/custom_material.py +++ b/examples/python/cohesive/custom_material.py @@ -1,193 +1,194 @@ #!/usr/bin/env python3 # pylint: disable=missing-module-docstring # pylint: disable=missing-function-docstring __copyright__ = ( "Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)" "Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)" ) __license__ = "LGPLv3" import numpy as np import akantu as aka spatial_dimension = 2 class LinearCohesive(aka.MaterialCohesive): """Material Cohesive Linear.""" def __init__(self, model, _id): """Construct the material and register the parameters to the parser.""" super().__init__(model, _id) super().registerParamReal( "G_c", aka._pat_readable | aka._pat_parsable, "Fracture energy" ) super().registerParamReal( "beta", aka._pat_readable | aka._pat_parsable, "" ) self.registerInternalReal("delta_max", 1) self.beta = 0. self.sigma_c = 0. self.delta_c = 0. self.G_c = 0. def initMaterial(self): """Initialize the parameters for the material.""" super().initMaterial() self.sigma_c = self.getReal("sigma_c") self.G_c = self.getReal("G_c") self.beta = self.getReal("beta") self.delta_c = 2 * self.G_c / self.sigma_c def checkInsertion(self, _check_only): """Check if need to insert a cohesive element.""" model = self.getModel() facets = self.getFacetFilter() inserter = model.getElementInserter() for type_facet in facets.elementTypes(dim=(spatial_dimension - 1)): facet_filter = facets(type_facet) nb_facet = facet_filter.shape[0] if nb_facet == 0: continue fe_engine = model.getFEEngine("FacetsFEEngine") facets_check = inserter.getCheckFacets(type_facet) insertion = inserter.getInsertionFacets(type_facet) nb_quad_facet = fe_engine.getNbIntegrationPoints(type_facet) normals = fe_engine.getNormalsOnIntegrationPoints(type_facet) facets_stresses = model.getStressOnFacets(type_facet).reshape( normals.shape[0] // nb_quad_facet, nb_quad_facet, 2, spatial_dimension, spatial_dimension, ) tangents = model.getTangents(type_facet) for facet, facet_stresses in zip(facet_filter, facets_stresses): if not facets_check[facet]: continue ref_stress = 0 for q, quad_stresses in enumerate(facet_stresses): current_quad = facet * nb_quad_facet + q normal = normals[current_quad, :].ravel() tangent = tangents[current_quad, :].ravel() stress_1 = quad_stresses.T[0] stress_2 = quad_stresses.T[1] avg_stress = stress_1 + stress_2 / 2.0 traction = avg_stress.dot(normal) n = traction.dot(normal) n = max(0, n) t = traction.dot(tangent) ref_stress = max( ref_stress, np.sqrt(n * n + t * t / (self.beta**2)) ) if ref_stress > self.sigma_c: insertion[facet] = True # constitutive law - def computeTraction(self, normals, el_type, ghost_type): + def computeTraction(self, el_type, ghost_type): """Compute the traction for a given opening.""" + normals = self.getNormals(el_type, ghost_type) openings = self.getOpening(el_type, ghost_type) tractions = self.getTraction(el_type, ghost_type) delta_max = self.getInternalReal("delta_max")(el_type) for el in range(normals.shape[0]): normal = normals[el].ravel() opening = openings[el].ravel() delta_n = opening.dot(normal) * normal delta_s = opening - delta_n delta = self.beta * np.linalg.norm(delta_s) ** 2 + \ np.linalg.norm(delta_n) ** 2 delta_max[el] = max(delta, delta_max[el]) tractions[el, :] = ( (delta * delta_s + delta_n) * self.sigma_c / delta * (1 - delta / self.delta_c) ) def allocator(_dim, _unused, model, _id): # NOQA """Register the material to the material factory.""" return LinearCohesive(model, _id) mat_factory = aka.MaterialFactory.getInstance() mat_factory.registerAllocator("local_cohesive", allocator) # ------------------------------------------------------------------------- # Initialization # ------------------------------------------------------------------------- aka.parseInput("local_material.dat") mesh = aka.Mesh(spatial_dimension) mesh.read("plate.msh") model = aka.SolidMechanicsModelCohesive(mesh) model.initFull(_analysis_method=aka._static, _is_extrinsic=True) model.initNewSolver(aka._explicit_lumped_mass) model.setBaseName("plate") model.addDumpFieldVector("displacement") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("blocked_dofs") model.setBaseNameToDumper("cohesive elements", "cohesive") model.addDumpFieldVectorToDumper("cohesive elements", "displacement") model.addDumpFieldToDumper("cohesive elements", "damage") model.addDumpFieldVectorToDumper("cohesive elements", "tractions") model.addDumpFieldVectorToDumper("cohesive elements", "opening") # ------------------------------------------------------------------------- # Boundary conditions # ------------------------------------------------------------------------- model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked") model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked") trac = np.zeros(spatial_dimension) traction = 0.095 trac[int(aka._y)] = traction model.getExternalForce()[:] = 0 model.applyBC(aka.FromTraction(trac), "Traction") print("Solve for traction ", traction) solver = model.getNonLinearSolver("static") solver.set("max_iterations", 100) solver.set("threshold", 1e-10) solver.set("convergence_type", aka.SolveConvergenceCriteria.residual) model.solveStep("static") model.dump() model.dump("cohesive elements") model.setTimeStep(model.getStableTimeStep() * 0.1) maxsteps = 100 for i in range(0, maxsteps): print("{0}/{1}".format(i, maxsteps)) model.checkCohesiveStress() model.solveStep("explicit_lumped") if i % 10 == 0: model.dump() model.dump("cohesive elements") diff --git a/src/model/phase_field/phase_field_model_inline_impl.hh b/src/model/phase_field/phase_field_model_inline_impl.hh index afc009270..0c7f2aaa6 100644 --- a/src/model/phase_field/phase_field_model_inline_impl.hh +++ b/src/model/phase_field/phase_field_model_inline_impl.hh @@ -1,90 +1,90 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . */ /* -------------------------------------------------------------------------- */ #include "aka_named_argument.hh" #include "phasefield_selector.hh" #include "phasefield_selector_tmpl.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__ #define __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline decltype(auto) PhaseFieldModel::getPhaseFields() { return make_dereference_adaptor(phasefields); } /* -------------------------------------------------------------------------- */ inline decltype(auto) PhaseFieldModel::getPhaseFields() const { return make_dereference_adaptor(phasefields); } /* -------------------------------------------------------------------------- */ inline PhaseField & PhaseFieldModel::getPhaseField(Idx mat_index) { - AKANTU_DEBUG_ASSERT(mat_index < phasefields.size(), + AKANTU_DEBUG_ASSERT(mat_index < Idx(phasefields.size()), "The model " << id << " has no phasefield no " << mat_index); return *phasefields[mat_index]; } /* -------------------------------------------------------------------------- */ inline const PhaseField & PhaseFieldModel::getPhaseField(Idx mat_index) const { - AKANTU_DEBUG_ASSERT(mat_index < phasefields.size(), + AKANTU_DEBUG_ASSERT(mat_index < Idx(phasefields.size()), "The model " << id << " has no phasefield no " << mat_index); return *phasefields[mat_index]; } /* -------------------------------------------------------------------------- */ inline PhaseField & PhaseFieldModel::getPhaseField(const std::string & name) { auto it = phasefields_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(), "The model " << id << " has no phasefield named " << name); return *phasefields[it->second]; } /* -------------------------------------------------------------------------- */ inline Int PhaseFieldModel::getPhaseFieldIndex(const std::string & name) const { auto it = phasefields_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(), "The model " << id << " has no phasefield named " << name); return it->second; } /* -------------------------------------------------------------------------- */ inline const PhaseField & PhaseFieldModel::getPhaseField(const std::string & name) const { auto it = phasefields_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(), "The model " << id << " has no phasefield named " << name); return *phasefields[it->second]; } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh b/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh index a5efd9ab0..e56aa8ced 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh @@ -1,367 +1,366 @@ /** * Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "material_anisotropic_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_ #define AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_ namespace akantu { struct EmptyIteratorContainer { using size_type = std::size_t; struct iterator { auto & operator++() { return *this; } Real operator*() { return 0; } bool operator!=(const iterator & /*unused*/) const { return true; } bool operator==(const iterator & /*unused*/) const { return false; } }; auto begin() const // NOLINT(readability-convert-member-functions-to-static) { return iterator(); } auto end() const // NOLINT(readability-convert-member-functions-to-static) { return iterator(); } }; } // namespace akantu namespace std { template <> struct iterator_traits<::akantu::EmptyIteratorContainer::iterator> { using iterator_category = forward_iterator_tag; using value_type = akantu::Real; using difference_type = std::ptrdiff_t; using pointer = akantu::Real *; using reference = akantu::Real &; }; } // namespace std namespace akantu { namespace { template void tensorPlus_(const Eigen::MatrixBase & A, Op && oper) { Vector A_eigs; A.eig(A_eigs); for (auto & ap : A_eigs) { oper(ap); } } template auto tensorPlus2(const Eigen::MatrixBase & A) { Real square = 0; tensorPlus_(A, [&](Real eig) { eig = std::max(eig, 0.); square += eig * eig; }); return square; } template auto tensorPlusTrace(const Eigen::MatrixBase & A) { Real trace_plus = 0.; Real trace_minus = 0.; tensorPlus_(A, [&](Real eig) { trace_plus += std::max(eig, 0.); trace_minus += std::min(eig, 0.); }); - return std::make_pair(trace_plus, trace_minus); + return std::pair(trace_plus, trace_minus); } template auto tensorPlusOp(const Eigen::MatrixBase & A, Eigen::MatrixBase & A_directions, Op && oper) { Vector A_eigs; Matrix A_diag; A.eig(A_eigs, A_directions); for (auto && [i, eig] : enumerate(A_eigs)) { A_diag(i, i) = oper(std::max(eig, 0.), i); } Matrix res = A_directions * A_diag * A_directions.transpose(); return res; } template auto tensorPlus(const Eigen::MatrixBase & A, Eigen::MatrixBase & A_directions) { return tensorPlusOp(A, A_directions, [](Real x, Real) { return x; }); } template auto tensorPlusOp(const Eigen::MatrixBase & A, Op && oper) { Matrix A_directions; A_directions.zero(); return tensorPlusOp(A, A_directions, std::forward(oper)); } template auto tensorPlus(const Eigen::MatrixBase & A) { return tensorPlusOp(A, [](Real x, Real) { return x; }); } template auto tensorSqrt(const Eigen::MatrixBase & A) { return tensorPlusOp(A, [](Real x, Int) { return std::sqrt(x); }); } } // namespace /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> MaterialAnisotropicDamage:: MaterialAnisotropicDamage(SolidMechanicsModel & model, const ID & id) : Parent(model, id), damage("damage_tensor", *this), elastic_stress("elastic_stress", *this), equivalent_strain("equivalent_strain", *this), trace_damage("trace_damage", *this), equivalent_strain_function(*this), damage_threshold_function(*this) { this->registerParam("Dc", Dc, _pat_parsable, "Critical damage"); this->damage.initialize(dim * dim); this->elastic_stress.initialize(dim * dim); this->equivalent_strain.initialize(1); this->trace_damage.initialize(1); this->trace_damage.initializeHistory(); } /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> template void MaterialAnisotropicDamage:: damageStress(Eigen::MatrixBase & sigma, const Eigen::MatrixBase & sigma_el, const Eigen::MatrixBase & D, Real TrD) { // σ_(n + 1) = (1 − D_(n + 1))^(1/2) σ~_(n + 1) (1 − D_(n + 1))^(1 / 2) // - ((1 − D_(n + 1)) : σ~_(n + 1))/ (3 - Tr(D_(n+1))) (1 − D_(n + 1)) // + 1/3 (1 - Tr(D_(n+1)) _+ + _-) I Matrix one_D = Matrix::Identity() - D; auto sqrt_one_D = tensorSqrt(one_D); Real Tr_sigma_plus; Real Tr_sigma_minus; std::tie(Tr_sigma_plus, Tr_sigma_minus) = tensorPlusTrace(sigma_el); auto I = Matrix::Identity(); sigma = sqrt_one_D * sigma_el * sqrt_one_D - (one_D.doubleDot(sigma_el) / (dim - TrD) * one_D) + 1. / dim * ((1 - TrD) * Tr_sigma_plus - Tr_sigma_minus) * I; } /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> template void MaterialAnisotropicDamage::computeStressOnQuad(Args && args) { auto & sigma = args["sigma"_n]; auto & grad_u = args["grad_u"_n]; auto & sigma_el = args["sigma_el"_n]; auto & epsilon_hat = args["epsilon_hat"_n]; auto & D = args["damage"_n]; auto & TrD_n_1 = args["TrD_n_1"_n]; auto & TrD = args["TrD"_n]; auto & equivalent_strain_data = args["equivalent_strain_data"_n]; auto & damage_threshold_data = args["damage_threshold_data"_n]; Matrix Dtmp; Real TrD_n_1_tmp; Matrix epsilon; // yes you read properly this is a label for a goto auto computeDamage = [&]() { MaterialElastic::computeStressOnQuad(args); epsilon = this->template gradUToEpsilon(grad_u); // evaluate the damage criteria epsilon_hat = equivalent_strain_function(epsilon, equivalent_strain_data); // evolve the damage if needed auto K_TrD = damage_threshold_function.K(TrD, damage_threshold_data); auto f = epsilon_hat - K_TrD; // if test function > 0 evolve the damage if (f > 0) { TrD_n_1_tmp = damage_threshold_function.K_inv(epsilon_hat, damage_threshold_data); auto epsilon_plus = tensorPlus(epsilon); auto delta_lambda = (TrD_n_1_tmp - TrD) / (epsilon_hat * epsilon_hat); Dtmp = D + delta_lambda * epsilon_plus; return true; } return false; }; // compute a temporary version of the new damage auto is_damage_updated = computeDamage(); if (is_damage_updated) { /// Check and correct for broken case if (Dtmp.trace() > Dc) { if (epsilon.trace() > 0) { // tensile loading auto kpa = this->kpa; auto lambda = this->lambda; // change kappa to Kappa_broken = (1-Dc) Kappa kpa = (1 - Dc) * kpa; this->E = 9 * kpa * (kpa - lambda) / (3 * kpa - lambda); this->nu = lambda / (3 * kpa - lambda); this->updateInternalParameters(); computeDamage(); } else if (std::abs(epsilon.trace()) < 1e-10) { // deviatoric case Matrix n; std::vector ns; tensorPlusOp(Dtmp, n, [&](Real x, Int i) { if (x > this->Dc) { ns.push_back(i); return this->Dc; } return x; }); } } TrD_n_1 = TrD_n_1_tmp; D = Dtmp; } else { TrD_n_1 = TrD; } // apply the damage to the stress damageStress(sigma, sigma_el, D, TrD_n_1); } /* -------------------------------------------------------------------------- */ template class EquivalentStrain, template class DamageThreshold, template class Parent> void MaterialAnisotropicDamage::computeStress(ElementType type, GhostType ghost_type) { for (auto && args : getArguments(type, ghost_type)) { computeStressOnQuad(args); } } /* -------------------------------------------------------------------------- */ /* EquivalentStrain functions */ /* -------------------------------------------------------------------------- */ template class EquivalentStrainMazars : public EmptyIteratorContainer { public: EquivalentStrainMazars(Material & /*mat*/) {} template Real operator()(const Eigen::MatrixBase & epsilon, Other &&... /*other*/) { - Real epsilon_hat = 0.; - std::tie(epsilon_hat, std::ignore) = tensorPlusTrace(epsilon); + auto && [epsilon_hat, _] = tensorPlusTrace(epsilon); return std::sqrt(epsilon_hat); } }; template class EquivalentStrainMazarsDruckerPrager : public EquivalentStrainMazars { public: EquivalentStrainMazarsDruckerPrager(Material & mat) : EquivalentStrainMazars(mat) { mat.registerParam("k", k, _pat_parsable, "k"); } template Real operator()(const Eigen::MatrixBase & epsilon, Real) { Real epsilon_hat = EquivalentStrainMazars::operator()(epsilon); epsilon_hat += k * epsilon.trace(); return epsilon_hat; } protected: Real k; }; /* -------------------------------------------------------------------------- */ /* DamageThreshold functions */ /* -------------------------------------------------------------------------- */ template class DamageThresholdLinear : public EmptyIteratorContainer { public: DamageThresholdLinear(Material & mat) : mat(mat) { mat.registerParam("A", A, _pat_parsable, "A"); mat.registerParam("K0", K0, _pat_parsable, "K0"); } template Real K(Real x, Other &&... /*other*/) { return 1. / A * x + K0; } template Real K_inv(Real x, Other &&... /*other*/) { return A * (x - K0); } private: Material & mat; Real A; Real K0; }; template class DamageThresholdTan : public EmptyIteratorContainer { public: DamageThresholdTan(Material & mat) : mat(mat) { mat.registerParam("a", a, _pat_parsable, "a"); mat.registerParam("A", A, _pat_parsable, "A"); mat.registerParam("K0", K0, _pat_parsable, "K0"); } template Real K(Real x, Other &&... /*other*/) { return a * std::tan(std::atan2(x, a) - std::atan2(K0, a)); } template Real K_inv(Real x, Other &&... /*other*/) { return a * A * (std::atan2(x, a) - std::atan2(K0, a)); } private: Material & mat; Real a{2.93e-4}; Real A{5e3}; Real K0{5e-5}; }; } // namespace akantu #endif /* AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_ */