diff --git a/examples/python/phase-field/material_static.dat b/examples/python/phase-field/material_static.dat index 3ae8da692..744e85b9d 100644 --- a/examples/python/phase-field/material_static.dat +++ b/examples/python/phase-field/material_static.dat @@ -1,16 +1,16 @@ material phasefield [ name = plate rho = 1. E = 210.0 nu = 0.3 - eta = 0.0 + eta = 1e-8 Plane_Stress = false ] phasefield exponential [ name = plate l0 = 0.0075 gc = 2.7e-3 E = 210.0 nu = 0.3 ] diff --git a/examples/python/phase-field/phasefield-static.py b/examples/python/phase-field/phasefield-static.py index c0e6efb2f..765a58422 100644 --- a/examples/python/phase-field/phasefield-static.py +++ b/examples/python/phase-field/phasefield-static.py @@ -1,132 +1,132 @@ #!/usr/bin/env python # coding: utf-8 """ phasefield-static.py: Static phase field example""" __author__ = "Mohit Pundir" __credits__ = [ "Mohit Pundir ", ] __copyright__ = ( "Copyright (©) 2018-2021 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 aka.parseInput("material_static.dat") dim = 2 mesh = aka.Mesh(dim) mesh.read("plate_static.msh") model = aka.CouplerSolidPhaseField(mesh) solid = model.getSolidMechanicsModel() phase = model.getPhaseFieldModel() solid.initFull(_analysis_method=aka._static) solver = solid.getNonLinearSolver("static") solver.set("max_iterations", 100) solver.set("threshold", 1e-8) solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) solid.getNewSolver( "linear_static", aka.TimeStepSolverType.static, aka.NonLinearSolverType.linear ) solid.setIntegrationScheme( "linear_static", "displacement", aka.IntegrationSchemeType.pseudo_time ) phase.initFull(_analysis_method=aka._static) phase.getNewSolver( "nonlinear_static", aka.TimeStepSolverType.static, aka.NonLinearSolverType.newton_raphson, ) phase.setIntegrationScheme( "nonlinear_static", "damage", aka.IntegrationSchemeType.pseudo_time ) solver = phase.getNonLinearSolver("nonlinear_static") solver.set("max_iterations", 100) solver.set("threshold", 1e-4) solver.set("convergence_type", aka.SolveConvergenceCriteria.solution) solid.applyBC(aka.FixedValue(0, aka._y), "bottom") solid.applyBC(aka.FixedValue(0, aka._x), "left") # Initialization for bulk vizualisation solid.setBaseName("phasefield-static") solid.addDumpFieldVector("displacement") solid.addDumpFieldVector("external_force") solid.addDumpField("strain") solid.addDumpField("stress") solid.addDumpField("damage") solid.addDumpField("blocked_dofs") nb_dofs = solid.getMesh().getNbNodes() * dim increment = solid.getIncrement() displacement = solid.getDisplacement() displacement = displacement.reshape(nb_dofs) blocked_dofs = solid.getBlockedDOFs() blocked_dofs = blocked_dofs.reshape(nb_dofs) damage = phase.getDamage() -tolerance = 1e-6 +tolerance = 1e-5 -steps = 1500 -increment = 1e-5 +steps = 1000 +increment = 5e-6 for n in range(steps): print("Computing iteration " + str(n + 1) + "/" + str(steps)) solid.applyBC(aka.IncrementValue(increment, aka._y), "top") mask = blocked_dofs == False # NOQA: E712 iiter = 0 error_disp = 1 error_dam = 1 displacement_prev = displacement[mask].copy() damage_prev = damage.copy() damage_prev = damage_prev # solve using staggered scheme while error_disp > tolerance or error_dam > tolerance: model.solve("linear_static", "") displacement_new = displacement[mask] damage_new = damage delta_disp = displacement_new - displacement_prev delta_dam = damage_new - damage_prev error_disp = np.linalg.norm(delta_disp) error_dam = np.linalg.norm(delta_dam) iiter += 1 displacement_prev = displacement_new.copy() damage_prev = damage_new.copy() print(error_dam, error_disp) if iiter > 500: raise Exception("Convergence not reached") if n % 50 == 0: solid.dump() solid.dump() diff --git a/packages/phase_field.cmake b/packages/phase_field.cmake index c23b2e38d..048e416e3 100644 --- a/packages/phase_field.cmake +++ b/packages/phase_field.cmake @@ -1,53 +1,54 @@ #=============================================================================== # @file phase_field.cmake # # @author Mohit Pundir # # @date creation: Fri Sep 03 2010 # @date last modification: Wed Jun 30 2021 # # @brief package description for phase field model # # # @section LICENSE # # Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . # #=============================================================================== package_declare(phase_field DEPENDS model_couplers DESCRIPTION "Use Phase Field package of Akantu") package_declare_sources(phase_field model/phase_field/phasefield.cc model/phase_field/phasefield.hh model/phase_field/phasefield_inline_impl.cc model/phase_field/phasefield_selector.hh model/phase_field/phasefield_selector_tmpl.hh model/phase_field/phasefields/phasefield_exponential.hh model/phase_field/phasefields/phasefield_exponential.cc + model/phase_field/phasefields/phasefield_exponential_inline_impl.hh model/phase_field/phase_field_model.cc model/phase_field/phase_field_model.hh model/phase_field/phase_field_model_inline_impl.cc model/model_couplers/coupler_solid_phasefield.hh model/model_couplers/coupler_solid_phasefield.cc ) diff --git a/packages/solid_mechanics.cmake b/packages/solid_mechanics.cmake index 113376256..c7b8776ac 100644 --- a/packages/solid_mechanics.cmake +++ b/packages/solid_mechanics.cmake @@ -1,118 +1,118 @@ #=============================================================================== # @file solid_mechanics.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Dec 04 2017 # @date last modification: Fri Mar 26 2021 # # @brief package description for core # # # @section LICENSE # # Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . # #=============================================================================== package_declare(solid_mechanics DEFAULT ON DESCRIPTION "Solid mechanics model" DEPENDS core lapack ) package_declare_sources(solid_mechanics model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.hh model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.hh model/solid_mechanics/solid_mechanics_model_io.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.hh model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.hh model/solid_mechanics/materials/material_damage/material_anisotropic_damage.cc model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.hh model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_phasefield.cc model/solid_mechanics/materials/material_damage/material_phasefield.hh - model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.cc + model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.hh model/solid_mechanics/materials/material_damage/material_von_mises_mazars.cc model/solid_mechanics/materials/material_damage/material_von_mises_mazars.hh model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh model/solid_mechanics/materials/material_non_local.hh model/solid_mechanics/materials/material_non_local_tmpl.hh model/solid_mechanics/materials/material_non_local_includes.hh ) package_declare_material_infos(solid_mechanics LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_extra_files_to_package(solid_mechanics SOURCES model/solid_mechanics/material_list.hh.in ) diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc index 458a4feb7..540882a72 100644 --- a/python/py_dof_manager.cc +++ b/python/py_dof_manager.cc @@ -1,216 +1,216 @@ /* * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "py_dof_manager.hh" #include "py_aka_array.hh" #include "py_akantu_pybind11_compatibility.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { class PySolverCallback : public SolverCallback { public: using SolverCallback::SolverCallback; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(MatrixType, SolverCallback, getMatrixType, matrix_id); } /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleMatrix, matrix_id); } /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleLumpedMatrix, matrix_id); } /// callback to assemble the residual (rhs) void assembleResidual() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleResidual, ); } /// callback for the predictor (in case of dynamic simulation) void predictor() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, predictor, ); } /// callback for the corrector (in case of dynamic simulation) void corrector() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, corrector, ); } void beforeSolveStep() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, beforeSolveStep, ); } void afterSolveStep(bool converged) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, afterSolveStep, converged); } }; class PyInterceptSolverCallback : public InterceptSolverCallback { public: using InterceptSolverCallback::InterceptSolverCallback; MatrixType getMatrixType(const ID & matrix_id) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(MatrixType, InterceptSolverCallback, getMatrixType, matrix_id); } void assembleMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleMatrix, matrix_id); } /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleLumpedMatrix, matrix_id); } void assembleResidual() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleResidual, ); } void predictor() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, predictor, ); } void corrector() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, corrector, ); } void beforeSolveStep() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, beforeSolveStep, ); } void afterSolveStep(bool converged) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, afterSolveStep, converged); } }; } // namespace /* -------------------------------------------------------------------------- */ void register_dof_manager(py::module & mod) { py::class_>(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) .def( "hasMatrix", [](DOFManager & self, const ID & name) -> bool { return self.hasMatrix(name); }, py::arg("name")) .def("assembleToResidual", &DOFManager::assembleToResidual, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("scale_factor") = 1.) .def("assembleToLumpedMatrix", &DOFManager::assembleToLumpedMatrix, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("lumped_mtx"), py::arg("scale_factor") = 1.) .def("assemblePreassembledMatrix", &DOFManager::assemblePreassembledMatrix, py::arg("matrix_id"), py::arg("terms")) .def("zeroResidual", &DOFManager::zeroResidual); - py::class_(mod, "NonLinearSolver") + py::class_(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); py::class_(mod, "TimeStepSolver") .def("getIntegrationScheme", &TimeStepSolver::getIntegrationScheme); py::class_(mod, "SolverCallback") .def(py::init_alias()) .def("getMatrixType", &SolverCallback::getMatrixType) .def("assembleMatrix", &SolverCallback::assembleMatrix) .def("assembleLumpedMatrix", &SolverCallback::assembleLumpedMatrix) .def("assembleResidual", [](SolverCallback & self) { self.assembleResidual(); }) .def("predictor", &SolverCallback::predictor) .def("corrector", &SolverCallback::corrector) .def("beforeSolveStep", &SolverCallback::beforeSolveStep) .def("afterSolveStep", &SolverCallback::afterSolveStep) .def_property_readonly("dof_manager", &SolverCallback::getSCDOFManager, py::return_value_policy::reference); py::class_(mod, "InterceptSolverCallback") .def(py::init_alias()); } } // namespace akantu diff --git a/python/py_mesh.cc b/python/py_mesh.cc index 6928a6cdc..6fbfccafc 100644 --- a/python/py_mesh.cc +++ b/python/py_mesh.cc @@ -1,259 +1,259 @@ /** * @file py_mesh.cc * * @author Guillaume Anciaux * @author Philip Mueller * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Mon Mar 15 2021 * * @brief pybind11 interface to Mesh * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { /* ------------------------------------------------------------------------ */ template auto register_element_type_map_array(py::module & mod, const std::string & name) { auto element_type_map_class = py::class_, std::shared_ptr>>( mod, ("ElementTypeMapArray" + name).c_str()) .def(py::init(), py::arg("id") = "by_element_type_array", py::arg("parent_id") = "no_parent") .def( "__call__", [](ElementTypeMapArray & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference, py::keep_alive<0, 1>()) .def( "elementTypes", [](ElementTypeMapArray & self, UInt _dim, GhostType _ghost_type, ElementKind _kind) -> std::vector { auto types = self.elementTypes(_dim, _ghost_type, _kind); std::vector _types; for (auto && t : types) { _types.push_back(t); } return _types; }, py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular) .def( "initialize", [](ElementTypeMapArray & self, const Mesh & mesh, GhostType ghost_type, Int nb_component, Int spatial_dimension, ElementKind element_kind, bool with_nb_element, bool with_nb_nodes_per_element, T default_value, bool do_not_default) { self.initialize( mesh, _ghost_type = ghost_type, _nb_component = nb_component, _spatial_dimension = (spatial_dimension == -2 ? mesh.getSpatialDimension() : spatial_dimension), _element_kind = element_kind, _with_nb_element = with_nb_element, _with_nb_nodes_per_element = with_nb_nodes_per_element, _default_value = default_value, _do_not_default = do_not_default); }, py::arg("mesh"), py::arg("ghost_type") = _casper, py::arg("nb_component") = 1, py::arg("spatial_dimension") = -2, py::arg("element_kind") = _ek_not_defined, py::arg("with_nb_element") = false, py::arg("with_nb_nodes_per_element") = false, py::arg("default_value") = T{}, py::arg("do_not_default") = false); return element_type_map_class; } } // namespace /* -------------------------------------------------------------------------- */ void register_mesh(py::module & mod) { py::class_(mod, "PeriodicSlaves") .def( "__iter__", [](Mesh::PeriodicSlaves & _this) { - py::make_iterator(_this.begin(), _this.end()); + return py::make_iterator(_this.begin(), _this.end()); }, py::keep_alive<0, 1>()); py::class_(mod, "MeshData") .def( "getElementalDataUInt", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference) .def( "getElementalDataReal", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference); py::class_(mod, "Mesh", py::multiple_inheritance()) .def(py::init(), py::arg("spatial_dimension"), py::arg("id") = "mesh") .def("read", &Mesh::read, py::arg("filename"), py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file") .def( "getNodes", [](Mesh & self) -> decltype(auto) { return self.getNodes(); }, py::return_value_policy::reference) .def("getNbNodes", &Mesh::getNbNodes) .def( "getConnectivity", [](Mesh & self, ElementType type) -> decltype(auto) { return self.getConnectivity(type); }, py::return_value_policy::reference) .def( "getConnectivities", [](Mesh & self) -> decltype(auto) { return self.getConnectivities(); }, py::return_value_policy::reference) .def( "addConnectivityType", [](Mesh & self, ElementType type, GhostType ghost_type) -> void { self.addConnectivityType(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("distribute", [](Mesh & self) { self.distribute(); }) .def("isDistributed", [](const Mesh & self) { return self.isDistributed(); }) .def("fillNodesToElements", &Mesh::fillNodesToElements, py::arg("dimension") = _all_dimensions) .def("getAssociatedElements", [](Mesh & self, const Idx & node, py::list list) { Array elements; self.getAssociatedElements(node, elements); for (auto && element : elements) { list.append(element); } }) .def("makePeriodic", [](Mesh & self, const SpatialDirection & direction) { self.makePeriodic(direction); }) .def( "getNbElement", [](Mesh & self, const Int spatial_dimension, GhostType ghost_type, ElementKind kind) { return self.getNbElement(spatial_dimension, ghost_type, kind); }, py::arg("spatial_dimension") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined) .def( "getNbElement", [](Mesh & self, ElementType type, GhostType ghost_type) { return self.getNbElement(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def_static( "getSpatialDimension", [](ElementType & type) { return Mesh::getSpatialDimension(type); }) .def( "getDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> decltype(auto) { return _this.getData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "hasDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> bool { return _this.hasData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); }) .def("getPeriodicMaster", &Mesh::getPeriodicMaster) .def("getPeriodicSlaves", &Mesh::getPeriodicSlaves) .def("isPeriodicSlave", &Mesh::isPeriodicSlave) .def("isPeriodicMaster", &Mesh::isPeriodicMaster) .def( "getMeshFacets", [](const Mesh & self) -> const Mesh & { return self.getMeshFacets(); }, py::return_value_policy::reference) .def("initMeshFacets", &Mesh::initMeshFacets, py::arg("id") = "mesh_facets", py::return_value_policy::reference); /* ------------------------------------------------------------------------ */ py::class_(mod, "MeshUtils") .def_static("buildFacets", &MeshUtils::buildFacets); py::class_(mod, "MeshAccessor") .def(py::init(), py::arg("mesh")) .def( "resizeConnectivity", [](MeshAccessor & self, Int new_size, ElementType type, GhostType gt) -> void { self.resizeConnectivity(new_size, type, gt); }, py::arg("new_size"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "resizeNodes", [](MeshAccessor & self, Int new_size) -> void { self.resizeNodes(new_size); }, py::arg("new_size")) .def("makeReady", &MeshAccessor::makeReady); register_element_type_map_array(mod, "Real"); register_element_type_map_array(mod, "UInt"); register_element_type_map_array(mod, "Int"); register_element_type_map_array(mod, "bool"); // register_element_type_map_array(mod, "String"); } } // namespace akantu diff --git a/src/model/phase_field/phasefield.cc b/src/model/phase_field/phasefield.cc index 4bc87e499..61302178f 100644 --- a/src/model/phase_field/phasefield.cc +++ b/src/model/phase_field/phasefield.cc @@ -1,282 +1,285 @@ /** * @file phasefield.cc * * @author Mohit Pundir * * @date creation: Fri Jun 19 2020 * @date last modification: Fri May 14 2021 * * @brief Implementation of the common part of the phasefield class * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "phasefield.hh" #include "phase_field_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PhaseField::PhaseField(PhaseFieldModel & model, Int dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Parsable(ParserType::_phasefield, id), id(id), fem(fe_engine), - model(model), spatial_dimension(this->model.getSpatialDimension()), + model(model), g_c("g_c", *this), + spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id), damage("damage", *this, dim, fe_engine, this->element_filter), phi("phi", *this, dim, fe_engine, this->element_filter), strain("strain", *this, dim, fe_engine, this->element_filter), driving_force("driving_force", *this, dim, fe_engine, this->element_filter), damage_energy("damage_energy", *this, dim, fe_engine, this->element_filter), damage_energy_density("damage_energy_density", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseField::PhaseField(PhaseFieldModel & model, const ID & id) : PhaseField(model, model.getSpatialDimension(), model.getMesh(), model.getFEEngine(), id) {} /* -------------------------------------------------------------------------- */ PhaseField::~PhaseField() = default; /* -------------------------------------------------------------------------- */ void PhaseField::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("l0", l0, Real(0.), _pat_parsable | _pat_readable, "length scale parameter"); registerParam("gc", g_c, _pat_parsable | _pat_readable, "critical local fracture energy density"); registerParam("E", E, _pat_parsable | _pat_readable, "Young's modulus"); registerParam("nu", nu, _pat_parsable | _pat_readable, "Poisson ratio"); damage.initialize(1); phi.initialize(1); driving_force.initialize(1); + g_c.initialize(1); + strain.initialize(spatial_dimension * spatial_dimension); damage_energy_density.initialize(1); damage_energy.initialize(spatial_dimension * spatial_dimension); } /* -------------------------------------------------------------------------- */ void PhaseField::initPhaseField() { AKANTU_DEBUG_IN(); this->phi.initializeHistory(); this->resizeInternals(); updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { it->second->resize(); } for (auto it = internal_vectors_int.begin(); it != internal_vectors_int.end(); ++it) { it->second->resize(); } for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) { it->second->resize(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::updateInternalParameters() { this->lambda = this->nu * this->E / ((1 + this->nu) * (1 - 2 * this->nu)); this->mu = this->E / (2 * (1 + this->nu)); } /* -------------------------------------------------------------------------- */ void PhaseField::computeAllDrivingForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Int spatial_dimension = model.getSpatialDimension(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { auto & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { continue; } computeDrivingForce(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) { auto & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { continue; } auto nb_element = elem_filter.size(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array nt_driving_force(nb_quadrature_points, nb_nodes_per_element); fem.computeNtb(driving_force(type, ghost_type), nt_driving_force, type, ghost_type, elem_filter); Array int_nt_driving_force(nb_element, nb_nodes_per_element); fem.integrate(nt_driving_force, int_nt_driving_force, nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalArrayLocalArray( int_nt_driving_force, internal_force, type, ghost_type, 1, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { auto & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { AKANTU_DEBUG_OUT(); return; } auto nb_element = elem_filter.size(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nt_b_n = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "N^t*b*N"); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); // damage_energy_density_on_qpoints = gc/l0 + phi = scalar auto & damage_energy_density_vect = damage_energy_density(type, ghost_type); // damage_energy_on_qpoints = gc*l0 = scalar auto & damage_energy_vect = damage_energy(type, ghost_type); fem.computeBtDB(damage_energy_vect, *bt_d_b, 2, type, ghost_type, elem_filter); fem.computeNtbN(damage_energy_density_vect, *nt_b_n, type, ghost_type, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{N}^t * \mathbf{w} * /// \mathbf{N}@f$ auto K_n = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_n"); fem.integrate(*nt_b_n, *K_n, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_n, type, _not_ghost, _symmetric, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{B}^t * \mathbf{W} * /// \mathbf{B}@f$ auto K_b = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_b"); fem.integrate(*bt_d_b, *K_b, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_b, type, _not_ghost, _symmetric, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::beforeSolveStep() { this->savePreviousState(); this->computeAllDrivingForces(_not_ghost); } /* -------------------------------------------------------------------------- */ void PhaseField::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void PhaseField::savePreviousState() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) { if (pair.second->hasHistory()) { pair.second->saveCurrentValues(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "PhaseField Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/model/phase_field/phasefield.hh b/src/model/phase_field/phasefield.hh index 1571ff642..e4922c661 100644 --- a/src/model/phase_field/phasefield.hh +++ b/src/model/phase_field/phasefield.hh @@ -1,299 +1,309 @@ /** * @file phasefield.hh * * @author Mohit Pundir * * @date creation: Fri Jun 19 2020 * @date last modification: Wed Jun 23 2021 * * @brief Mother class for all phasefield laws * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_factory.hh" #include "data_accessor.hh" #include "parsable.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASEFIELD_HH__ #define __AKANTU_PHASEFIELD_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class PhaseFieldModel; class PhaseField; } // namespace akantu namespace akantu { template using InternalPhaseField = InternalFieldTmpl; +template <> +inline void +ParameterTyped>::setAuto( + const ParserParameter & in_param) { + Parameter::setAuto(in_param); + RandomParameter random_param = in_param; + param.setRandomDistribution(random_param); +} + using PhaseFieldFactory = Factory; class PhaseField : public DataAccessor, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PhaseField(const PhaseField & phase) = delete; PhaseField & operator=(const PhaseField & phase) = delete; /// Initialize phasefield with defaults PhaseField(PhaseFieldModel & model, const ID & id = ""); /// Initialize phasefield with custom mesh & fe_engine PhaseField(PhaseFieldModel & model, Int dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor ~PhaseField() override; protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(InternalPhaseField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } template void unregisterInternal(InternalPhaseField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } /// initialize the phasefield computed parameter virtual void initPhaseField(); /// virtual void beforeSolveStep(); /// virtual void afterSolveStep(); /// assemble the residual for this phasefield virtual void assembleInternalForces(GhostType ghost_type); /// assemble the stiffness matrix for this phasefield virtual void assembleStiffnessMatrix(GhostType ghost_type); /// compute the driving force for this phasefield virtual void computeAllDrivingForces(GhostType ghost_type = _not_ghost); /// save the phi in the phi internal field if needed virtual void savePreviousState(); /// add an element to the local mesh filter inline UInt addElement(ElementType type, UInt element, GhostType ghost_type); inline UInt addElement(const Element & element); /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; protected: /// resize the internals arrrays virtual void resizeInternals(); /// function called to updatet the internal parameters when the /// modifiable parameters are modified virtual void updateInternalParameters(); // constitutive law for driving force virtual void computeDrivingForce(ElementType /* el_type */, GhostType /* ghost_type */ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline Int getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(Model, model, const PhaseFieldModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO(Strain, strain, const ElementTypeMapArray &); AKANTU_GET_MACRO_NOT_CONST(Strain, strain, ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); AKANTU_GET_MACRO_NOT_CONST(Damage, damage, ElementTypeMapArray &); AKANTU_GET_MACRO(Damage, damage, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, Idx); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); template const InternalPhaseField & getInternal(const ID & id) const; template InternalPhaseField & getInternal(const ID & id); template inline bool isInternal(const ID & id, ElementKind element_kind) const; template inline void setParam(const ID & param, T value); inline const Parameter & getParam(const ID & param) const; template void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_int; std::map *> internal_vectors_bool; protected: ID id; /// Link to the fem object in the model FEEngine & fem; /// phasefield name std::string name; /// The model to whch the phasefield belong PhaseFieldModel & model; /// length scale parameter Real l0; /// critical energy release rate - Real g_c; + // Real g_c; + RandomInternalField g_c; /// Young's modulus Real E; /// Poisson ratio Real nu; /// Lame's first parameter Real lambda; /// Lame's second paramter Real mu; /// spatial dimension Int spatial_dimension; /// list of element handled by the phasefield ElementTypeMapArray element_filter; /// damage arrays ordered by element types InternalPhaseField damage; /// phi arrays ordered by element types InternalPhaseField phi; /// strain arrays ordered by element types InternalPhaseField strain; /// driving force ordered by element types InternalPhaseField driving_force; /// damage energy ordered by element types InternalPhaseField damage_energy; /// damage energy density ordered by element types InternalPhaseField damage_energy_density; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const PhaseField & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "phasefield_inline_impl.cc" #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ #define PHASEFIELD_DEFAULT_ALLOCATOR(id, phase_name) \ [](const ID &, PhaseFieldModel & model, \ const ID & id) -> std::unique_ptr { \ return std::make_unique(model, id); \ } #define INSTANTIATE_PHASEFIELD(id, phase_name) \ static bool phasefield_is_allocated_##id [[gnu::unused]] = \ PhaseFieldFactory::getInstance().registerAllocator( \ #id, PHASEFIELD_DEFAULT_ALLOCATOR(id, phase_name)) #endif diff --git a/src/model/phase_field/phasefields/phasefield_exponential.cc b/src/model/phase_field/phasefields/phasefield_exponential.cc index e1f5b5be1..5d428cc66 100644 --- a/src/model/phase_field/phasefields/phasefield_exponential.cc +++ b/src/model/phase_field/phasefields/phasefield_exponential.cc @@ -1,71 +1,82 @@ /** * @file phasefield_exponential.cc * * @author Mohit Pundir * * @date creation: Fri Jun 19 2020 * @date last modification: Wed Jun 23 2021 * * @brief Specialization of the phasefield law class for exponential type * law * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "phasefield_exponential.hh" +#include "aka_common.hh" +#include namespace akantu { /* -------------------------------------------------------------------------- */ PhaseFieldExponential::PhaseFieldExponential(PhaseFieldModel & model, const ID & id) : PhaseField(model, id) {} /* -------------------------------------------------------------------------- */ void PhaseFieldExponential::updateInternalParameters() { PhaseField::updateInternalParameters(); - Matrix d = - Matrix::Identity(spatial_dimension, spatial_dimension) * this->g_c * - this->l0; - damage_energy.set(d); + for (const auto & type : + element_filter.elementTypes(spatial_dimension, _not_ghost)) { + for (auto && tuple : zip(make_view(this->damage_energy(type, _not_ghost), + spatial_dimension, spatial_dimension), + this->g_c(type, _not_ghost))) { + Matrix d = + Matrix::Identity(spatial_dimension, spatial_dimension) * + std::get<1>(tuple) * this->l0; + std::get<0>(tuple) = d; + } + } } /* -------------------------------------------------------------------------- */ void PhaseFieldExponential::computeDrivingForce(ElementType el_type, GhostType ghost_type) { for (auto && tuple : zip(this->phi(el_type, ghost_type), this->phi.previous(el_type, ghost_type), this->driving_force(el_type, ghost_type), this->damage_energy_density(el_type, ghost_type), make_view(this->strain(el_type, ghost_type), - spatial_dimension, spatial_dimension))) { + spatial_dimension, spatial_dimension), + this->g_c(el_type, ghost_type))) { computePhiOnQuad(std::get<4>(tuple), std::get<0>(tuple), std::get<1>(tuple)); - computeDamageEnergyDensityOnQuad(std::get<0>(tuple), std::get<3>(tuple)); + computeDamageEnergyDensityOnQuad(std::get<0>(tuple), std::get<3>(tuple), + std::get<5>(tuple)); computeDrivingForceOnQuad(std::get<0>(tuple), std::get<2>(tuple)); } } INSTANTIATE_PHASEFIELD(exponential, PhaseFieldExponential); } // namespace akantu diff --git a/src/model/phase_field/phasefields/phasefield_exponential.hh b/src/model/phase_field/phasefields/phasefield_exponential.hh index 969373b49..99cef22c4 100644 --- a/src/model/phase_field/phasefields/phasefield_exponential.hh +++ b/src/model/phase_field/phasefields/phasefield_exponential.hh @@ -1,133 +1,76 @@ /** * @file phasefield_exponential.hh * * @author Mohit Pundir * * @date creation: Fri Jun 19 2020 * @date last modification: Wed Jun 23 2021 * * @brief Phasefield law for approximating discrete crack as an exponential * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "phasefield.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ #define __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ namespace akantu { class PhaseFieldExponential : public PhaseField { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PhaseFieldExponential(PhaseFieldModel & model, const ID & id = ""); ~PhaseFieldExponential() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: void computePhiOnQuad(const Matrix & /*strain_quad*/, Real & /*phi_quad*/, Real & /*phi_hist_quad*/); void computeDrivingForce(ElementType /*el_type*/, GhostType /*ghost_type*/) override; inline void computeDrivingForceOnQuad(const Real & /*phi_quad*/, Real & /*driving_force_quad*/); inline void computeDamageEnergyDensityOnQuad(const Real & /*phi_quad*/, - Real & /*dam_energy_quad*/); + Real & /*dam_energy_quad*/, + const Real & /*g_c_quad*/); public: void updateInternalParameters() override; }; -/* -------------------------------------------------------------------------- */ -inline void -PhaseFieldExponential::computeDrivingForceOnQuad(const Real & phi_quad, - Real & driving_force_quad) { - driving_force_quad = 2.0 * phi_quad; -} +} // namespace akantu /* -------------------------------------------------------------------------- */ -inline void PhaseFieldExponential::computeDamageEnergyDensityOnQuad( - const Real & phi_quad, Real & dam_energy_quad) { - dam_energy_quad = 2.0 * phi_quad + this->g_c / this->l0; -} - /* -------------------------------------------------------------------------- */ -inline void -PhaseFieldExponential::computePhiOnQuad(const Matrix & strain_quad, - Real & phi_quad, Real & phi_hist_quad) { - - Matrix strain_plus(spatial_dimension, spatial_dimension); - // Matrix strain_minus(spatial_dimension, spatial_dimension); - Matrix strain_dir(spatial_dimension, spatial_dimension); - Matrix strain_diag_plus(spatial_dimension, spatial_dimension); - // Matrix strain_diag_minus(spatial_dimension, spatial_dimension); - - Vector strain_values(spatial_dimension); - - Real trace_plus; - - strain_plus.zero(); - // strain_minus.zero(); - strain_dir.zero(); - strain_values.zero(); - strain_diag_plus.zero(); - // strain_diag_minus.zero(); - - strain_quad.eig(strain_values, strain_dir); - - for (Int i = 0; i < spatial_dimension; i++) { - strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); - // strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); - } - - strain_plus = strain_dir * strain_diag_plus * strain_dir.transpose(); - - // Nicolas: @mohit is the second transpose really here ? - // strain_minus = strain_dir * strain_diag_minus * strain_dir.transpose(); - - trace_plus = std::max(Real(0.), strain_quad.trace()); - // trace_minus = std::min(Real(0.), strain_quad.trace()); - - auto I = Matrix::Identity(spatial_dimension, spatial_dimension); - - Matrix sigma_plus = I * lambda * trace_plus + 2 * mu * strain_plus; - // Matrix sigma_minus = I * lambda * trace_minus + 2 * mu * - // strain_minus; - - phi_quad = 0.5 * sigma_plus.doubleDot(strain_quad); - if (phi_quad < phi_hist_quad) { - phi_quad = phi_hist_quad; - } -} - -} // namespace akantu +#include "phasefield_exponential_inline_impl.hh" #endif diff --git a/src/model/phase_field/phasefields/phasefield_exponential_inline_impl.hh b/src/model/phase_field/phasefields/phasefield_exponential_inline_impl.hh new file mode 100644 index 000000000..2187170fc --- /dev/null +++ b/src/model/phase_field/phasefields/phasefield_exponential_inline_impl.hh @@ -0,0 +1,71 @@ +#include "phasefield_exponential.hh" + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +inline void +PhaseFieldExponential::computeDrivingForceOnQuad(const Real & phi_quad, + Real & driving_force_quad) { + driving_force_quad = 2.0 * phi_quad; +} + +/* -------------------------------------------------------------------------- */ +inline void PhaseFieldExponential::computeDamageEnergyDensityOnQuad( + const Real & phi_quad, Real & dam_energy_quad, const Real & g_c_quad) { + dam_energy_quad = 2.0 * phi_quad + g_c_quad / this->l0; +} + +/* -------------------------------------------------------------------------- */ +inline void +PhaseFieldExponential::computePhiOnQuad(const Matrix & strain_quad, + Real & phi_quad, Real & phi_hist_quad) { + Matrix strain_plus(spatial_dimension, spatial_dimension); + Matrix strain_minus(spatial_dimension, spatial_dimension); + Matrix strain_dir(spatial_dimension, spatial_dimension); + Matrix strain_diag_plus(spatial_dimension, spatial_dimension); + Matrix strain_diag_minus(spatial_dimension, spatial_dimension); + + Vector strain_values(spatial_dimension); + + Real trace_plus, trace_minus; + + strain_plus.zero(); + strain_minus.zero(); + strain_dir.zero(); + strain_values.zero(); + strain_diag_plus.zero(); + strain_diag_minus.zero(); + + strain_quad.eig(strain_values, strain_dir); + + for (Int i = 0; i < spatial_dimension; i++) { + strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); + strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); + } + + Matrix mat_tmp(spatial_dimension, spatial_dimension); + Matrix sigma_plus(spatial_dimension, spatial_dimension); + Matrix sigma_minus(spatial_dimension, spatial_dimension); + + strain_plus = strain_dir * strain_diag_plus * strain_dir.transpose(); + strain_minus = strain_dir * strain_diag_minus * strain_dir.transpose(); + + trace_plus = std::max(Real(0.), strain_quad.trace()); + trace_minus = std::min(Real(0.), strain_quad.trace()); + + for (Int i = 0; i < spatial_dimension; i++) { + for (Int j = 0; j < spatial_dimension; j++) { + sigma_plus(i, j) = static_cast(i == j) * lambda * trace_plus + + 2 * mu * strain_plus(i, j); + sigma_minus(i, j) = static_cast(i == j) * lambda * trace_minus + + 2 * mu * strain_minus(i, j); + } + } + + phi_quad = sigma_plus.doubleDot(strain_plus) / 2.; + if (phi_quad < phi_hist_quad) { + phi_quad = phi_hist_quad; + } +} + +} // namespace akantu diff --git a/src/model/solid_mechanics/materials/internal_field.hh b/src/model/solid_mechanics/materials/internal_field.hh index 574618d3c..410461ab6 100644 --- a/src/model/solid_mechanics/materials/internal_field.hh +++ b/src/model/solid_mechanics/materials/internal_field.hh @@ -1,248 +1,250 @@ /** * @file internal_field.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Mar 26 2021 * * @brief Material internal properties * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_common.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_INTERNAL_FIELD_HH_ #define AKANTU_INTERNAL_FIELD_HH_ namespace akantu { class Material; class FEEngine; /** * class for the internal fields of materials * to store values for each quadrature */ -template +template class InternalFieldTmpl : public ElementTypeMapArray { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: + using Material = Material_; + InternalFieldTmpl(const ID & id, Material & material); ~InternalFieldTmpl() override; /// This constructor is only here to let cohesive elements compile InternalFieldTmpl(const ID & id, Material & material, FEEngine & fem, const ElementTypeMapArray & element_filter); /// More general constructor InternalFieldTmpl(const ID & id, Material & material, Int dim, FEEngine & fem, const ElementTypeMapArray & element_filter); InternalFieldTmpl(const ID & id, const InternalFieldTmpl & other); auto operator=(const InternalFieldTmpl &) -> InternalFieldTmpl = delete; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to reset the FEEngine for the internal field virtual void setFEEngine(FEEngine & fe_engine); /// function to reset the element kind for the internal virtual void setElementKind(ElementKind element_kind); /// initialize the field to a given number of component virtual void initialize(UInt nb_component); /// activate the history of this field virtual void initializeHistory(); /// resize the arrays and set the new element to 0 virtual void resize(); /// set the field to a given value v virtual void setDefaultValue(const T & v); /// reset all the fields to the default value virtual void reset(); /// save the current values in the history virtual void saveCurrentValues(); /// restore the previous values from the history virtual void restorePreviousValues(); /// remove the quadrature points corresponding to suppressed elements virtual void removeIntegrationPoints(const ElementTypeMapArray & new_numbering); /// print the content void printself(std::ostream & stream, int /*indent*/ = 0) const override; /// get the default value inline operator T() const; virtual auto getFEEngine() -> FEEngine & { return *fem; } virtual auto getFEEngine() const -> const FEEngine & { return *fem; } protected: /// initialize the arrays in the ElementTypeMapArray void internalInitialize(Int nb_component); /// set the values for new internals virtual void setArrayValues(T * begin, T * end); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get filter types for range loop decltype(auto) elementTypes(GhostType ghost_type = _not_ghost) const { return ElementTypeMapArray::elementTypes( _spatial_dimension = this->spatial_dimension, _element_kind = this->element_kind, _ghost_type = ghost_type); } /// get filter types for range loop decltype(auto) filterTypes(GhostType ghost_type = _not_ghost) const { return this->element_filter.elementTypes( _spatial_dimension = this->spatial_dimension, _element_kind = this->element_kind, _ghost_type = ghost_type); } /// get the array for a given type of the element_filter decltype(auto) getFilter(ElementType type, GhostType ghost_type = _not_ghost) const { return (this->element_filter(type, ghost_type)); } /// get the Array corresponding to the type en ghost_type specified virtual auto operator()(ElementType type, GhostType ghost_type = _not_ghost) -> Array & { return ElementTypeMapArray::operator()(type, ghost_type); } virtual auto operator()(ElementType type, GhostType ghost_type = _not_ghost) const -> const Array & { return ElementTypeMapArray::operator()(type, ghost_type); } virtual auto previous(ElementType type, GhostType ghost_type = _not_ghost) -> Array & { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return this->previous_values->operator()(type, ghost_type); } virtual auto previous(ElementType type, GhostType ghost_type = _not_ghost) const -> const Array & { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return this->previous_values->operator()(type, ghost_type); } virtual auto previous() -> InternalFieldTmpl & { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return *(this->previous_values); } virtual auto previous() const -> const InternalFieldTmpl & { AKANTU_DEBUG_ASSERT(previous_values != nullptr, "The history of the internal " << this->getID() << " has not been activated"); return *(this->previous_values); } /// check if the history is used or not auto hasHistory() const -> bool { return (previous_values != nullptr); } /// get the kind treated by the internal AKANTU_GET_MACRO_AUTO(ElementKind, element_kind); /// return the number of components AKANTU_GET_MACRO_AUTO(NbComponent, nb_component); /// return the spatial dimension corresponding to the internal element type /// loop filter AKANTU_GET_MACRO_AUTO(SpatialDimension, spatial_dimension); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the material for which this is an internal parameter Material & material; /// the fem containing the mesh and the element informations FEEngine * fem{nullptr}; /// Element filter if needed const ElementTypeMapArray & element_filter; /// default value T default_value{}; /// spatial dimension of the element to consider Int spatial_dimension{0}; /// ElementKind of the element to consider ElementKind element_kind{_ek_regular}; /// Number of component of the internal field Int nb_component{0}; /// Is the field initialized bool is_init{false}; /// previous values std::unique_ptr> previous_values; }; /// standard output stream operator template inline auto operator<<(std::ostream & stream, const InternalFieldTmpl & _this) -> std::ostream & { _this.printself(stream); return stream; } template using InternalField = InternalFieldTmpl; } // namespace akantu #endif /* AKANTU_INTERNAL_FIELD_HH_ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc b/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc index 0d3d4ca9b..69efd7bd9 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc @@ -1,84 +1,89 @@ /** * @file material_phasefield.cc * * @author Mohit Pundir * * @date creation: Mon Dec 13 2010 * @date last modification: Fri Apr 02 2021 * * @brief Specialization of the material class for the phasefield material * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "material_phasefield.hh" +#include "aka_common.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialPhaseField::MaterialPhaseField(SolidMechanicsModel & model, const ID & id) - : Parent(model, id) { - - AKANTU_DEBUG_IN(); - + : Parent(model, id), effective_damage("effective_damage", *this) { this->registerParam("eta", eta, Real(0.), _pat_parsable, "eta"); this->damage.initialize(0); - - AKANTU_DEBUG_OUT(); + this->effective_damage.initialize(1); } -/* -------------------------------------------------------------------------- */ template void MaterialPhaseField::computeStress(ElementType el_type, GhostType ghost_type) { - for (auto && args : Parent::getArguments(el_type, ghost_type)) { + computeEffectiveDamage(el_type, ghost_type); + + for (auto && args : getArguments(el_type, ghost_type)) { computeStressOnQuad(args); } } /* -------------------------------------------------------------------------- */ template void MaterialPhaseField::computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type) { - AKANTU_DEBUG_IN(); + computeEffectiveDamage(el_type, ghost_type); for (auto && args : - Parent::getArgumentsTangent(tangent_matrix, el_type, ghost_type)) { + getArgumentsTangent(tangent_matrix, el_type, ghost_type)) { computeTangentModuliOnQuad(args); } +} - AKANTU_DEBUG_OUT(); +/* -------------------------------------------------------------------------- */ +template +void MaterialPhaseField::computeEffectiveDamage(ElementType el_type, + GhostType ghost_type) { + for (auto && args : getArguments(el_type, ghost_type)) { + computeEffectiveDamageOnQuad(args); + } } /* -------------------------------------------------------------------------- */ template class MaterialPhaseField<1>; template class MaterialPhaseField<2>; template class MaterialPhaseField<3>; static bool material_is_allocated_phasefield = instantiateMaterial("phasefield"); } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh b/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh index 0f3639a94..fa08b2679 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh @@ -1,85 +1,112 @@ /** * @file material_phasefield.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Apr 02 2021 * * @brief Phasefield damage law * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_common.hh" #include "material.hh" #include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_PHASEFIELD_HH__ #define __AKANTU_MATERIAL_PHASEFIELD_HH__ namespace akantu { template class MaterialPhaseField : public MaterialDamage { using Parent = MaterialDamage; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialPhaseField(SolidMechanicsModel & model, const ID & id = ""); ~MaterialPhaseField() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// compute the tangent stiffness matrix for an element type void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override; + /* ------------------------------------------------------------------------ */ + decltype(auto) getArguments(ElementType el_type, + GhostType ghost_type = _not_ghost) { + return zip_append(Parent::getArguments(el_type, ghost_type), + "effective_damage"_n = make_view( + this->effective_damage(el_type, ghost_type))); + } + + decltype(auto) getArgumentsTangent(Array & tangent_matrix, + ElementType el_type, + GhostType ghost_type) { + return zip_append( + Parent::getArgumentsTangent(tangent_matrix, el_type, ghost_type), + "effective_damage"_n = + make_view(this->effective_damage(el_type, ghost_type))); + } + protected: /// constitutive law for a given quadrature point template inline void computeStressOnQuad(Args && args); /// compute the tangent stiffness matrix for a given quadrature point template inline void computeTangentModuliOnQuad(Args && args); - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ -public: + /// Compute the effective damage + void computeEffectiveDamage(ElementType el_type, + GhostType ghost_type = _not_ghost); + + template inline void computeEffectiveDamageOnQuad(Args && args); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: Real eta; + + // effective damage to conserve stiffness in compression + InternalField effective_damage; }; } // namespace akantu -#include "material_phasefield_inline_impl.cc" +/* -------------------------------------------------------------------------- */ +/* inline functions */ +/* -------------------------------------------------------------------------- */ + +#include "material_phasefield_inline_impl.hh" + #endif /* __AKANTU_MATERIAL_PHASEFIELD_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.cc b/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh similarity index 52% rename from src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.cc rename to src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh index ebd3eb00d..87995613f 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh @@ -1,53 +1,97 @@ /** * @file material_phasefield_inline_impl.cc * * @author Mohit Pundir * * @date creation: Mon Dec 13 2010 * @date last modification: Fri Apr 02 2021 * * @brief Implementation of the inline functions of the material phasefield * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "material_phasefield.hh" + +#ifndef AKANTU_MATERIAL_PHASEFIELD_INLINE_IMPL_HH_ +#define AKANTU_MATERIAL_PHASEFIELD_INLINE_IMPL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { template template inline void MaterialPhaseField::computeStressOnQuad(Args && args) { MaterialElastic::computeStressOnQuad(args); - auto && dam = args["damage"_n]; + auto && dam = args["effective_damage"_n]; args["sigma"_n] *= (1 - dam) * (1 - dam) + eta; } /* -------------------------------------------------------------------------- */ template template void MaterialPhaseField::computeTangentModuliOnQuad(Args && args) { MaterialElastic::computeTangentModuliOnQuad(args); - auto dam = args["damage"_n]; + auto dam = args["effective_damage"_n]; args["tangent_moduli"_n] *= (1 - dam) * (1 - dam) + eta; } +/* -------------------------------------------------------------------------- */ +template +template +inline void +MaterialPhaseField::computeEffectiveDamageOnQuad(Args && args) { + using Mat = Matrix; + + auto strain = this->template gradUToEpsilon(args["grad_u"_n]); + + Mat strain_dir; + Vector strain_values; + strain.eig(strain_values, strain_dir); + + Mat strain_diag_plus; + Mat strain_diag_minus; + + for (UInt i = 0; i < dim; i++) { + strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); + strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); + } + + Mat strain_plus = strain_dir * strain_diag_plus * strain_dir.transpose(); + Mat strain_minus = strain_dir * strain_diag_minus * strain_dir.transpose(); + + auto trace_plus = std::max(Real(0.), strain.trace()); + auto trace_minus = std::min(Real(0.), strain.trace()); + + Mat sigma_plus = + Mat::Identity() * trace_plus * this->lambda + 2. * this->mu * strain_plus; + Mat sigma_minus = Mat::Identity() * trace_minus * this->lambda + + 2. * this->mu * strain_minus; + + auto strain_energy_plus = sigma_plus.doubleDot(strain_plus) / 2.; + auto strain_energy_minus = sigma_minus.doubleDot(strain_minus) / 2.; + + args["effective_damage"_n] = + args["damage"_n] * (strain_energy_minus < strain_energy_plus); +} + } // namespace akantu +#endif diff --git a/src/model/solid_mechanics/materials/random_internal_field.hh b/src/model/solid_mechanics/materials/random_internal_field.hh index 58fa09ff6..78f9995cd 100644 --- a/src/model/solid_mechanics/materials/random_internal_field.hh +++ b/src/model/solid_mechanics/materials/random_internal_field.hh @@ -1,106 +1,108 @@ /** * @file random_internal_field.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Mar 26 2021 * * @brief Random internal material parameter * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_common.hh" #include "aka_random_generator.hh" #include "internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_RANDOM_INTERNAL_FIELD_HH_ #define AKANTU_RANDOM_INTERNAL_FIELD_HH_ namespace akantu { /** * class for the internal fields of materials with a random * distribution */ template class BaseField = InternalField, template class Generator = RandomGenerator> class RandomInternalField : public BaseField { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - RandomInternalField(const ID & id, Material & material); + using ParentMaterial = typename BaseField::Material; + + RandomInternalField(const ID & id, ParentMaterial & material); ~RandomInternalField() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ RandomInternalField operator=(const RandomInternalField &) = delete; public: AKANTU_GET_MACRO(RandomParameter, random_parameter, const RandomParameter &); /// initialize the field to a given number of component void initialize(UInt nb_component) override; /// set the field to a given value void setDefaultValue(const T & value) override; /// set the specified random distribution to a given parameter void setRandomDistribution(const RandomParameter & param); /// print the content void printself(std::ostream & stream, int indent = 0) const override; protected: void setArrayValues(T * begin, T * end) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: inline operator Real() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// random parameter containing the distribution and base value RandomParameter random_parameter; }; /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const RandomInternalField & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AKANTU_RANDOM_INTERNAL_FIELD_HH_ */ diff --git a/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh b/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh index dc94496a6..e79be2cb5 100644 --- a/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh +++ b/src/model/solid_mechanics/materials/random_internal_field_tmpl.hh @@ -1,126 +1,126 @@ /** * @file random_internal_field_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Fri Mar 26 2021 * * @brief Random internal material parameter implementation * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_common.hh" #include "aka_random_generator.hh" #include "internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_RANDOM_INTERNAL_FIELD_TMPL_HH_ #define AKANTU_RANDOM_INTERNAL_FIELD_TMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> RandomInternalField::RandomInternalField( - const ID & id, Material & material) + const ID & id, ParentMaterial & material) : BaseField(id, material), random_parameter(T()) {} /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> RandomInternalField::~RandomInternalField() = default; /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> void RandomInternalField::initialize( UInt nb_component) { this->internalInitialize(nb_component); } /* ------------------------------------------------------------------------ */ template class BaseField, template class Generator> void RandomInternalField::setDefaultValue( const T & value) { random_parameter.setBaseValue(value); this->reset(); } /* ------------------------------------------------------------------------ */ template class BaseField, template class Generator> void RandomInternalField::setRandomDistribution( const RandomParameter & param) { random_parameter = param; this->reset(); } /* ------------------------------------------------------------------------ */ template class BaseField, template class Generator> void RandomInternalField::printself( std::ostream & stream, int indent [[gnu::unused]]) const { stream << "RandomInternalField [ "; random_parameter.printself(stream); stream << " ]"; #if !defined(AKANTU_NDEBUG) if (AKANTU_DEBUG_TEST(dblDump)) { stream << std::endl; - InternalField::printself(stream, indent); + BaseField::printself(stream, indent); } #endif } /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> void RandomInternalField::setArrayValues(T * begin, T * end) { random_parameter.template setValues(begin, end); } /* -------------------------------------------------------------------------- */ template class BaseField, template class Generator> inline RandomInternalField::operator Real() const { return random_parameter.getBaseValue(); } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped>::setAuto( const ParserParameter & in_param) { Parameter::setAuto(in_param); RandomParameter r = in_param; param.setRandomDistribution(r); } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* AKANTU_RANDOM_INTERNAL_FIELD_TMPL_HH_ */ diff --git a/src/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc index fbe6a4a8b..9e9036286 100644 --- a/src/solver/sparse_matrix_aij.cc +++ b/src/solver/sparse_matrix_aij.cc @@ -1,302 +1,307 @@ /** * @file sparse_matrix_aij.cc * * @author Nicolas Richart * * @date creation: Fri Aug 21 2015 * @date last modification: Fri Jul 24 2020 * * @brief Implementation of the AIJ sparse matrix * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "sparse_matrix_aij.hh" #include "aka_iterators.hh" #include "aka_math.hh" #include "dof_manager_default.hh" #include "dof_synchronizer.hh" #include "solver_vector_default.hh" #include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(DOFManagerDefault & dof_manager, const MatrixType & matrix_type, const ID & id) : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager), irn(0, 1, id + ":irn"), jcn(0, 1, id + ":jcn"), a(0, 1, id + ":a") {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id) : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager), irn(matrix.irn, id + ":irn"), jcn(matrix.jcn, id + ":jcn"), a(matrix.a, id + ":a") {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::~SparseMatrixAIJ() = default; /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::applyBoundary(Real block_val) { AKANTU_DEBUG_IN(); const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); auto begin = blocked_dofs.begin(); auto end = blocked_dofs.end(); auto is_blocked = [&](auto && i) -> bool { auto il = this->dof_manager.globalToLocalEquationNumber(i); return std::binary_search(begin, end, il); }; for (auto && ij_a : zip(irn, jcn, a)) { auto ni = std::get<0>(ij_a) - 1; auto nj = std::get<1>(ij_a) - 1; if (is_blocked(ni) or is_blocked(nj)) { std::get<2>(ij_a) = std::get<0>(ij_a) != std::get<1>(ij_a) ? 0. : this->dof_manager.isLocalOrMasterDOF( this->dof_manager.globalToLocalEquationNumber(ni)) ? block_val : 0.; } } this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); auto m = this->size_; auto & comm = dof_manager.getCommunicator(); // write header if (comm.whoAmI() == 0) { outfile << "%%MatrixMarket matrix coordinate pattern"; if (this->matrix_type == _symmetric) { outfile << " symmetric"; } else { outfile << " general"; } outfile << std::endl; outfile << m << " " << m << " " << this->nb_non_zero << std::endl; } for (auto p : arange(comm.getNbProc())) { // write content if (comm.whoAmI() == p) { for (auto && data : zip(this->irn, this->jcn)) { outfile << std::get<0>(data) << " " << std::get<1>(data) << " 1" << std::endl; } } comm.barrier(); } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); auto & comm = dof_manager.getCommunicator(); // open and set the properties of the stream std::ofstream outfile; if (0 == comm.whoAmI()) { outfile.open(filename.c_str()); } else { outfile.open(filename.c_str(), std::ios_base::app); } outfile.precision(std::numeric_limits::digits10); // write header decltype(nb_non_zero) nnz = this->nb_non_zero; comm.allReduce(nnz); if (comm.whoAmI() == 0) { outfile << "%%MatrixMarket matrix coordinate real"; if (this->matrix_type == _symmetric) { outfile << " symmetric"; } else { outfile << " general"; } outfile << std::endl; outfile << this->size_ << " " << this->size_ << " " << nnz << std::endl; } for (auto p : arange(comm.getNbProc())) { // write content if (comm.whoAmI() == p) { for (auto && data : zip(this->irn, this->jcn, this->a)) { outfile << std::get<0>(data) << " " << std::get<1>(data) << " " << std::get<2>(data) << std::endl; } } comm.barrier(); } // time to end outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const Array & x, Array & y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); + Array tmp(y); + tmp.zero(); y *= beta; auto x_it = make_view(x).begin(); auto y_it = make_view(y).begin(); for (auto && data : zip(this->irn, this->jcn, this->a)) { auto i = this->dof_manager.globalToLocalEquationNumber(std::get<0>(data) - 1); auto j = this->dof_manager.globalToLocalEquationNumber(std::get<1>(data) - 1); const auto & A = std::get<2>(data); y_it[i] += alpha * A * x_it[j]; if ((this->matrix_type == _symmetric) && (i != j)) { y_it[j] += alpha * A * x_it[i]; } } if (this->dof_manager.hasSynchronizer()) { - this->dof_manager.getSynchronizer().reduceSynchronizeArray(y); + this->dof_manager.getSynchronizer().reduceSynchronizeArray( + tmp); } + y += tmp; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const SolverVector & _x, SolverVector & _y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); auto && x = aka::as_type(_x).getVector(); auto && y = aka::as_type(_y).getVector(); this->matVecMul(x, y, alpha, beta); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyContent(const SparseMatrix & matrix) { AKANTU_DEBUG_IN(); const auto & mat = aka::as_type(matrix); AKANTU_DEBUG_ASSERT(nb_non_zero == mat.getNbNonZero(), "The to matrix don't have the same profiles"); memcpy(a.data(), mat.getA().data(), nb_non_zero * sizeof(Real)); this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyProfile(const SparseMatrix & other) { const auto & A = aka::as_type(other); SparseMatrix::clearProfile(); this->irn.copy(A.irn); this->jcn.copy(A.jcn); this->irn_jcn_k.clear(); Idx i; Idx j; Idx k; for (auto && data : enumerate(irn, jcn)) { std::tie(k, i, j) = data; this->irn_jcn_k[this->key(i - 1, j - 1)] = k; } this->nb_non_zero = this->irn.size(); this->a.resize(this->nb_non_zero); this->a.set(0.); this->size_ = A.size_; this->profile_release = A.profile_release; this->value_release++; } /* -------------------------------------------------------------------------- */ template void SparseMatrixAIJ::addMeToTemplated(MatrixType & B, Real alpha) const { Idx i; Idx j; Real A_ij; for (auto && tuple : zip(irn, jcn, a)) { std::tie(i, j, A_ij) = tuple; B.add(i - 1, j - 1, alpha * A_ij); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::addMeTo(SparseMatrix & B, Real alpha) const { if (aka::is_of_type(B)) { this->addMeToTemplated(aka::as_type(B), alpha); } else { // this->addMeToTemplated(*this, alpha); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::mul(Real alpha) { this->a *= alpha; this->value_release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::set(Real val) { a.set(val); this->value_release++; } /* -------------------------------------------------------------------------- */ bool SparseMatrixAIJ::isFinite() const { return this->a.isFinite(); } } // namespace akantu diff --git a/test/test_model/test_phase_field_model/CMakeLists.txt b/test/test_model/test_phase_field_model/CMakeLists.txt index 5d962ac21..8de9302cb 100644 --- a/test/test_model/test_phase_field_model/CMakeLists.txt +++ b/test/test_model/test_phase_field_model/CMakeLists.txt @@ -1,55 +1,61 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Fri Sep 03 2010 # @date last modification: Tue Jun 08 2021 # # @brief test for the phase field model # # # @section LICENSE # # Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . # #=============================================================================== register_test(test_phasefield_selector SOURCES test_phasefield_selector.cc FILES_TO_COPY phasefield_selector.dat phasefield_selector.msh PACKAGE phase_field ) register_test(test_phase_solid_coupling SOURCES test_phase_solid_coupling.cc FILES_TO_COPY material_coupling.dat test_one_element.msh PACKAGE phase_field ) +register_test(test_phase_field_anisotropic + SOURCES test_phase_field_anisotropic.cc + FILES_TO_COPY material_coupling.dat test_one_element.msh + PACKAGE phase_field + ) + register_test(test_phase_solid_explicit SOURCES test_phase_solid_explicit.cc FILES_TO_COPY material_coupling.dat test_one_element.msh PACKAGE phase_field ) register_test(test_multi_material SOURCES test_multi_material.cc FILES_TO_COPY material_multiple.dat test_two_element.msh PACKAGE phase_field ) diff --git a/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc b/test/test_model/test_phase_field_model/test_phase_field_anisotropic.cc similarity index 61% copy from test/test_model/test_phase_field_model/test_phase_solid_explicit.cc copy to test/test_model/test_phase_field_model/test_phase_field_anisotropic.cc index a3bb93478..f83dec4ec 100644 --- a/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc +++ b/test/test_model/test_phase_field_model/test_phase_field_anisotropic.cc @@ -1,154 +1,148 @@ -/** - * @file test_phase_solid_explicit.cc - * - * @author Mohit Pundir - * - * @date creation: Sun Feb 28 2021 - * @date last modification: Fri Jun 25 2021 - * - * @brief test of the class PhaseFieldModel on the 2d square - * - * - * @section LICENSE - * - * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * 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_common.hh" #include "coupler_solid_phasefield.hh" +#include "material.hh" +#include "material_phasefield.hh" #include "non_linear_solver.hh" +#include "phase_field_model.hh" +#include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ +#include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; -const Int spatial_dimension = 2; +const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { - std::ofstream os("data-explicit.csv"); - os << "#strain stress damage analytical_sigma analytical_damage error_stress " - "error_damage" - << std::endl; + std::ofstream os("data.csv"); + os << "#strain stress damage analytical_sigma analytical_damage" << std::endl; initialize("material_coupling.dat", argc, argv); Mesh mesh(spatial_dimension); mesh.read("test_one_element.msh"); CouplerSolidPhaseField coupler(mesh); auto & model = coupler.getSolidMechanicsModel(); auto & phase = coupler.getPhaseFieldModel(); - model.initFull(_analysis_method = _explicit_lumped_mass); - - Real time_factor = 0.8; - Real stable_time_step = model.getStableTimeStep() * time_factor; - model.setTimeStep(stable_time_step); + model.initFull(_analysis_method = _static); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); model.setBaseName("phase_solid"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpFieldVector("displacement"); model.addDumpField("damage"); model.dump(); - Int nbSteps = 1000; + UInt nbSteps = 1500; Real increment = 1e-4; auto & stress = model.getMaterial(0).getArray("stress", _quadrangle_4); auto & damage = model.getMaterial(0).getArray("damage", _quadrangle_4); Real analytical_damage{0.}; Real analytical_sigma{0.}; auto & phasefield = phase.getPhaseField(0); const Real E = phasefield.getParam("E"); const Real nu = phasefield.getParam("nu"); Real c22 = E * (1 - nu) / ((1 + nu) * (1 - 2 * nu)); const Real gc = phasefield.getParam("gc"); const Real l0 = phasefield.getParam("l0"); Real error_stress{0.}; + Real error_damage{0.}; - for (Int s = 0; s < nbSteps; ++s) { - Real axial_strain = increment * s; + Real max_strain{0.}; + + for (UInt s = 0; s < nbSteps; ++s) { + Real axial_strain{0.}; + if (s < 500) { + axial_strain = increment * s; + } else if (s < 1000) { + axial_strain = (1500 - 2 * double(s)) * increment; + } else { + axial_strain = (3 * double(s) - 3500) * increment; + } applyDisplacement(model, axial_strain); - coupler.solve("explicit_lumped", "static"); + if (axial_strain > max_strain) { + max_strain = axial_strain; + } + + coupler.solve("static", "static"); - analytical_damage = axial_strain * axial_strain * c22 / - (gc / l0 + axial_strain * axial_strain * c22); - analytical_sigma = - c22 * axial_strain * (1 - analytical_damage) * (1 - analytical_damage); + analytical_damage = max_strain * max_strain * c22 / + (gc / l0 + max_strain * max_strain * c22); + if (axial_strain < 0.) { + analytical_sigma = c22 * axial_strain; + } else { + analytical_sigma = c22 * axial_strain * (1 - analytical_damage) * + (1 - analytical_damage); + } error_stress = std::abs(analytical_sigma - stress(0, 3)) / analytical_sigma; error_damage = std::abs(analytical_damage - damage(0)) / analytical_damage; - if (error_damage > 1e-8 and error_stress > 1e-8) { + if ((error_damage > 1e-8 or error_stress > 1e-8) and + std::abs(axial_strain) < 1e-13) { + std::cerr << std::left << std::setw(15) + << "Error damage: " << error_damage << std::endl; + std::cerr << std::left << std::setw(15) + << "Error stress: " << error_stress << std::endl; return EXIT_FAILURE; } os << axial_strain << " " << stress(0, 3) << " " << damage(0) << " " << analytical_sigma << " " << analytical_damage << " " << error_stress << " " << error_damage << std::endl; model.dump(); } os.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel & model, Real & increment) { auto & displacement = model.getDisplacement(); auto & positions = model.getMesh().getNodes(); auto & blocked_dofs = model.getBlockedDOFs(); - for (Int n = 0; n < model.getMesh().getNbNodes(); ++n) { + for (UInt n = 0; n < model.getMesh().getNbNodes(); ++n) { if (positions(n, 1) == -0.5) { displacement(n, 0) = 0; displacement(n, 1) = 0; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } else { displacement(n, 0) = 0; displacement(n, 1) = increment; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } } } + +/* -------------------------------------------------------------------------- */ diff --git a/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc b/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc index 7cbd2f82c..d422887cb 100644 --- a/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc +++ b/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc @@ -1,273 +1,277 @@ /** * @file test_phase_solid_coupling.cc * * @author Mohit Pundir * * @date creation: Sun Jan 06 2019 * @date last modification: Wed Mar 03 2021 * * @brief test of the class PhaseFieldModel on the 2d square * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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_common.hh" #include "material.hh" #include "material_phasefield.hh" #include "non_linear_solver.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; const Int spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); void computeStrainOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, GhostType); void computeDamageOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, GhostType); void gradUToEpsilon(const Matrix &, Matrix &); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { std::ofstream os("data.csv"); os << "#strain stress damage analytical_sigma analytical_damage" << std::endl; initialize("material_coupling.dat", argc, argv); Mesh mesh(spatial_dimension); mesh.read("test_one_element.msh"); SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _static); PhaseFieldModel phase(mesh); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); model.setBaseName("phase_solid"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpFieldVector("displacement"); model.addDumpField("damage"); model.dump(); Int nbSteps = 1000; Real increment = 1e-4; auto & stress = model.getMaterial(0).getArray("stress", _quadrangle_4); auto & damage = model.getMaterial(0).getArray("damage", _quadrangle_4); Real analytical_damage{0.}; Real analytical_sigma{0.}; auto & phasefield = phase.getPhaseField(0); const Real E = phasefield.getParam("E"); const Real nu = phasefield.getParam("nu"); Real c22 = E * (1 - nu) / ((1 + nu) * (1 - 2 * nu)); const Real gc = phasefield.getParam("gc"); const Real l0 = phasefield.getParam("l0"); Real error_stress{0.}; Real error_damage{0.}; for (Int s = 0; s < nbSteps; ++s) { Real axial_strain = increment * s; applyDisplacement(model, axial_strain); model.solveStep(); computeStrainOnQuadPoints(model, phase, _not_ghost); phase.solveStep(); computeDamageOnQuadPoints(model, phase, _not_ghost); model.assembleInternalForces(); analytical_damage = axial_strain * axial_strain * c22 / (gc / l0 + axial_strain * axial_strain * c22); analytical_sigma = c22 * axial_strain * (1 - analytical_damage) * (1 - analytical_damage); error_stress = std::abs(analytical_sigma - stress(0, 3)) / analytical_sigma; error_damage = std::abs(analytical_damage - damage(0)) / analytical_damage; - if (error_damage > 1e-8 and error_stress > 1e-8) { + if (error_damage > 1e-8 or error_stress > 1e-8) { + std::cerr << std::left << std::setw(15) + << "Error damage: " << error_damage << std::endl; + std::cerr << std::left << std::setw(15) + << "Error stress: " << error_stress << std::endl; return EXIT_FAILURE; } os << axial_strain << " " << stress(0, 3) << " " << damage(0) << " " << analytical_sigma << " " << analytical_damage << " " << error_stress << " " << error_damage << std::endl; model.dump(); } os.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel & model, Real & increment) { auto & displacement = model.getDisplacement(); auto & positions = model.getMesh().getNodes(); auto & blocked_dofs = model.getBlockedDOFs(); for (Int n = 0; n < model.getMesh().getNbNodes(); ++n) { if (positions(n, 1) == -0.5) { displacement(n, 0) = 0; displacement(n, 1) = 0; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } else { displacement(n, 0) = 0; displacement(n, 1) = increment; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } } } /* -------------------------------------------------------------------------- */ void computeStrainOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, GhostType ghost_type) { auto & mesh = solid.getMesh(); auto nb_materials = solid.getNbMaterials(); auto nb_phasefields = phase.getNbPhaseFields(); AKANTU_DEBUG_ASSERT( nb_phasefields == nb_materials, "The number of phasefields and materials should be equal"); for (auto index : arange(nb_materials)) { auto & material = solid.getMaterial(index); for (auto index2 : arange(nb_phasefields)) { auto & phasefield = phase.getPhaseField(index2); if (phasefield.getName() == material.getName()) { auto & strain_on_qpoints = phasefield.getStrain(); auto & gradu_on_qpoints = material.getGradU(); for (const auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type); auto & gradu_on_qpoints_vect = gradu_on_qpoints(type, ghost_type); for (auto && values : zip(make_view(strain_on_qpoints_vect, spatial_dimension, spatial_dimension), make_view(gradu_on_qpoints_vect, spatial_dimension, spatial_dimension))) { auto & strain = std::get<0>(values); auto & grad_u = std::get<1>(values); Material::gradUToEpsilon(grad_u, strain); } } break; } } } } /* -------------------------------------------------------------------------- */ void computeDamageOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, GhostType ghost_type) { auto & fem = phase.getFEEngine(); auto & mesh = phase.getMesh(); auto nb_materials = solid.getNbMaterials(); auto nb_phasefields = phase.getNbPhaseFields(); AKANTU_DEBUG_ASSERT( nb_phasefields == nb_materials, "The number of phasefields and materials should be equal"); for (auto index : arange(nb_materials)) { auto & material = solid.getMaterial(index); for (auto index2 : arange(nb_phasefields)) { auto & phasefield = phase.getPhaseField(index2); if (phasefield.getName() == material.getName()) { switch (spatial_dimension) { case 1: { auto & mat = dynamic_cast &>(material); auto & solid_damage = mat.getDamage(); for (const auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); fem.interpolateOnIntegrationPoints( phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } case 2: { auto & mat = dynamic_cast &>(material); auto & solid_damage = mat.getDamage(); for (const auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); fem.interpolateOnIntegrationPoints( phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } default: break; } } } } } /* -------------------------------------------------------------------------- */ void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (Int i = 0; i < spatial_dimension; ++i) { for (Int j = 0; j < spatial_dimension; ++j) epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } } diff --git a/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc b/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc index a3bb93478..076e67d32 100644 --- a/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc +++ b/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc @@ -1,154 +1,158 @@ /** * @file test_phase_solid_explicit.cc * * @author Mohit Pundir * * @date creation: Sun Feb 28 2021 * @date last modification: Fri Jun 25 2021 * * @brief test of the class PhaseFieldModel on the 2d square * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 "coupler_solid_phasefield.hh" #include "non_linear_solver.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; const Int spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { std::ofstream os("data-explicit.csv"); os << "#strain stress damage analytical_sigma analytical_damage error_stress " "error_damage" << std::endl; initialize("material_coupling.dat", argc, argv); Mesh mesh(spatial_dimension); mesh.read("test_one_element.msh"); CouplerSolidPhaseField coupler(mesh); auto & model = coupler.getSolidMechanicsModel(); auto & phase = coupler.getPhaseFieldModel(); model.initFull(_analysis_method = _explicit_lumped_mass); Real time_factor = 0.8; Real stable_time_step = model.getStableTimeStep() * time_factor; model.setTimeStep(stable_time_step); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); model.setBaseName("phase_solid"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpFieldVector("displacement"); model.addDumpField("damage"); model.dump(); Int nbSteps = 1000; Real increment = 1e-4; auto & stress = model.getMaterial(0).getArray("stress", _quadrangle_4); auto & damage = model.getMaterial(0).getArray("damage", _quadrangle_4); Real analytical_damage{0.}; Real analytical_sigma{0.}; auto & phasefield = phase.getPhaseField(0); const Real E = phasefield.getParam("E"); const Real nu = phasefield.getParam("nu"); Real c22 = E * (1 - nu) / ((1 + nu) * (1 - 2 * nu)); const Real gc = phasefield.getParam("gc"); const Real l0 = phasefield.getParam("l0"); Real error_stress{0.}; Real error_damage{0.}; for (Int s = 0; s < nbSteps; ++s) { Real axial_strain = increment * s; applyDisplacement(model, axial_strain); coupler.solve("explicit_lumped", "static"); analytical_damage = axial_strain * axial_strain * c22 / (gc / l0 + axial_strain * axial_strain * c22); analytical_sigma = c22 * axial_strain * (1 - analytical_damage) * (1 - analytical_damage); error_stress = std::abs(analytical_sigma - stress(0, 3)) / analytical_sigma; error_damage = std::abs(analytical_damage - damage(0)) / analytical_damage; if (error_damage > 1e-8 and error_stress > 1e-8) { + std::cerr << std::left << std::setw(15) + << "Error damage: " << error_damage << std::endl; + std::cerr << std::left << std::setw(15) + << "Error stress: " << error_stress << std::endl; return EXIT_FAILURE; } os << axial_strain << " " << stress(0, 3) << " " << damage(0) << " " << analytical_sigma << " " << analytical_damage << " " << error_stress << " " << error_damage << std::endl; model.dump(); } os.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel & model, Real & increment) { auto & displacement = model.getDisplacement(); auto & positions = model.getMesh().getNodes(); auto & blocked_dofs = model.getBlockedDOFs(); for (Int n = 0; n < model.getMesh().getNbNodes(); ++n) { if (positions(n, 1) == -0.5) { displacement(n, 0) = 0; displacement(n, 1) = 0; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } else { displacement(n, 0) = 0; displacement(n, 1) = increment; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } } }