diff --git a/test/test_model/test_phase_field_model/CMakeLists.txt b/test/test_model/test_phase_field_model/CMakeLists.txt index 06f932900..b2fa84d97 100644 --- a/test/test_model/test_phase_field_model/CMakeLists.txt +++ b/test/test_model/test_phase_field_model/CMakeLists.txt @@ -1,39 +1,46 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Thu Dec 20 2018 # @date last modification: Thu Dec 20 2018 # # @brief test for the phase field model # # @section LICENSE # # Copyright (©) 2010-2018 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 core ) register_test(test_phase_solid_coupling SOURCES test_phase_solid_coupling.cc FILES_TO_COPY material_coupling.dat test_one_element.msh PACKAGE core ) +register_test(test_phase_solid_explicit + SOURCES test_phase_solid_explicit.cc + FILES_TO_COPY material_coupling.dat test_one_element.msh + PACKAGE core + ) + + register_test(test_multi_material SOURCES test_multi_material.cc FILES_TO_COPY material_multiple.dat test_two_element.msh PACKAGE core ) diff --git a/test/test_model/test_phase_field_model/test_multi_material.cc b/test/test_model/test_phase_field_model/test_multi_material.cc index adc892777..612c13114 100644 --- a/test/test_model/test_phase_field_model/test_multi_material.cc +++ b/test/test_model/test_phase_field_model/test_multi_material.cc @@ -1,258 +1,130 @@ /** * @file tets_phase_field_2d.cc * * @author Mohit Pundir * * @date creation: Mon Oct 1 2018 * * @brief test of the class PhaseFieldModel on the 2d square * * @section LICENSE * * Copyright (©) 2015 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 "non_linear_solver.hh" +#include "coupler_solid_phasefield.hh" #include "solid_mechanics_model.hh" #include "phase_field_model.hh" #include "material.hh" #include "material_phasefield.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); -void computeStrainOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const GhostType &); -void computeDamageOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const GhostType &); -void gradUToEpsilon(const Matrix &, Matrix &); /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { initialize("material_multiple.dat", argc, argv); Mesh mesh(spatial_dimension); mesh.read("test_two_element.msh"); - SolidMechanicsModel model(mesh); + CouplerSolidPhaseField coupler(mesh); + auto & model = coupler.getSolidMechanicsModel(); + auto & phase = coupler.getPhaseFieldModel(); + auto && mat_selector = std::make_shared>( "physical_names", model); model.setMaterialSelector(mat_selector); model.initFull(_analysis_method = _explicit_lumped_mass); Real time_step = model.getStableTimeStep(); time_step *= 0.8; model.setTimeStep(time_step); - - PhaseFieldModel phase(mesh); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); model.setBaseName("multi_material"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpField("damage"); model.addDumpFieldVector("displacement"); model.addDumpField("blocked_dofs"); model.dump(); - UInt nbSteps = 10000; - Real increment = 1e-5; + UInt nbSteps = 1000; + Real increment = 1e-4; for (UInt 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(); + coupler.solve(); + model.dump(); } 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 (UInt n = 0; n < model.getMesh().getNbNodes(); ++n) { if (positions(n, 1) == -1) { displacement(n, 1) = 0; blocked_dofs(n, 1) = true; displacement(n, 0) = 0; blocked_dofs(n ,0) = true; } else if (positions(n, 1) == 1) { displacement(n, 0) = 0; displacement(n, 1) = increment; blocked_dofs(n, 0) = true; blocked_dofs(n ,1) = true; } else { displacement(n, 0) = 0; blocked_dofs(n, 0) = true; } } } - -/* -------------------------------------------------------------------------- */ -void computeStrainOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, - const 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 (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); - gradUToEpsilon(grad_u, strain); - } - } - - break; - } - - } - } -} - - -/* -------------------------------------------------------------------------- */ -void computeDamageOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, - const 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 = static_cast &>(material); - auto & solid_damage = mat.getDamage(); - auto & phase_damage = phasefield.getDamage(); - - for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { - auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); - auto & phase_damage_on_qpoints_vect = phase_damage(type, ghost_type); - - fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, - 1, type, ghost_type); - - } - - break; - } - case 2: { - auto & mat = static_cast &>(material); - auto & solid_damage = mat.getDamage(); - auto & phase_damage = phasefield.getDamage(); - - for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { - auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); - auto & phase_damage_on_qpoints_vect = phase_damage(type, ghost_type); - - fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, - 1, type, ghost_type); - - } - break; - } - default: - auto & mat = static_cast &>(material); - break; - } - - } - } - } - -} - - -/* -------------------------------------------------------------------------- */ -void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { - for (UInt i=0; i < spatial_dimension; ++i) { - for (UInt 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_coupling.cc b/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc index 453c96809..03cc09ad6 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,268 +1,280 @@ /** * @file test_phase_field_coupling.cc * * @author Mohit Pundir * * @date creation: Thu Feb 25 2021 * * @brief test of the class PhaseFieldModel on the 2d square * * @section LICENSE * * Copyright (©) 2015 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 "non_linear_solver.hh" #include "solid_mechanics_model.hh" #include "phase_field_model.hh" #include "material.hh" #include "material_phasefield.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); void computeStrainOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const GhostType &); void computeDamageOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const 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(); UInt 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 & mat = model.getMaterial(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 (UInt 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 > 0.01) { + return EXIT_FAILURE; + } + os << axial_strain << " " << stress(0, 3) << " " << damage(0) << " " << analytical_sigma << " " << analytical_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 (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; } } } /* -------------------------------------------------------------------------- */ void computeStrainOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, const 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 (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); gradUToEpsilon(grad_u, strain); } } break; } } } } /* -------------------------------------------------------------------------- */ void computeDamageOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, const 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 = static_cast &>(material); auto & solid_damage = mat.getDamage(); auto & phase_damage = phasefield.getDamage(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); auto & phase_damage_on_qpoints_vect = phase_damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } case 2: { auto & mat = static_cast &>(material); auto & solid_damage = mat.getDamage(); auto & phase_damage = phasefield.getDamage(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); auto & phase_damage_on_qpoints_vect = phase_damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } default: auto & mat = static_cast &>(material); break; } } } } } /* -------------------------------------------------------------------------- */ void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i=0; i < spatial_dimension; ++i) { for (UInt 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 new file mode 100644 index 000000000..f9fc9eef9 --- /dev/null +++ b/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc @@ -0,0 +1,161 @@ +/** + * @file test_phase_field_explicit.cc + * + * @author Mohit Pundir + * + * @date creation: Thu Feb 28 2021 + * + * @brief test of the class PhaseFieldModel on the 2d square + * + * @section LICENSE + * + * Copyright (©) 2015 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 "non_linear_solver.hh" +#include "coupler_solid_phasefield.hh" +#include "solid_mechanics_model.hh" +#include "phase_field_model.hh" +#include "material.hh" +#include "material_phasefield.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ + +using namespace akantu; +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; + + 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(); + + UInt 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 & mat = model.getMaterial(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 (UInt s = 0; s < nbSteps; ++s) { + Real axial_strain = increment * s; + applyDisplacement(model, axial_strain); + + coupler.solve(); + + 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 > 0.01) { + 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 (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; + } + } +}