diff --git a/examples/phase_field/phase_field_notch.cc b/examples/phase_field/phase_field_notch.cc index 14deb3da7..daa52d11c 100644 --- a/examples/phase_field/phase_field_notch.cc +++ b/examples/phase_field/phase_field_notch.cc @@ -1,120 +1,182 @@ /** * @file phase_field_notch.cc * * @author Mohit Pundir * * @date creation: Tue Oct 02 2018 * @date last modification: Wed Apr 07 2021 * * @brief Example of phase field model * * * @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 "phase_field_model.hh" #include "solid_mechanics_model.hh" +#include "group_manager.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; using clk = std::chrono::high_resolution_clock; using second = std::chrono::duration; using millisecond = std::chrono::duration; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ +class PhaseFieldElementFilter : public GroupManager::ClusteringFilter { +public: + PhaseFieldElementFilter(const PhaseFieldModel & model, + const Real max_damage = 1.) + : model(model), is_unbroken(max_damage) {} + + + bool operator()(const Element & el) const override { + + const Array & mat_indexes = + model.getPhaseFieldByElement(el.type, el.ghost_type); + const Array & mat_loc_num = + model.getPhaseFieldLocalNumbering(el.type, el.ghost_type); + + const auto & mat = model.getPhaseField(mat_indexes(el.element)); + + UInt el_index = mat_loc_num(el.element); + UInt nb_quad_per_element = + model.getFEEngine("PhaseFieldFEEngine") + .getNbIntegrationPoints(el.type, el.ghost_type); + + const Array & damage_array = mat.getDamage(el.type, el.ghost_type); + + AKANTU_DEBUG_ASSERT(nb_quad_per_element * el_index < damage_array.size(), + "This quadrature point is out of range"); + + const Real * element_damage = + damage_array.storage() + nb_quad_per_element * el_index; + + UInt unbroken_quads = std::count_if( + element_damage, element_damage + nb_quad_per_element, is_unbroken); + + return (unbroken_quads > 0); + } + +private: + struct IsUnbrokenFunctor { + IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {} + bool operator()(const Real & x) const { return x > max_damage; } + const Real max_damage; + }; + + const PhaseFieldModel & model; + const IsUnbrokenFunctor is_unbroken; + +}; + int main(int argc, char * argv[]) { initialize("material_notch.dat", argc, argv); // create mesh Mesh mesh(spatial_dimension); mesh.read("square_notch.msh"); CouplerSolidPhaseField coupler(mesh); auto & model = coupler.getSolidMechanicsModel(); auto & phase = coupler.getPhaseFieldModel(); model.initFull(_analysis_method = _static); auto && mat_selector = std::make_shared>("physical_names", model); model.setMaterialSelector(mat_selector); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); model.applyBC(BC::Dirichlet::FixedValue(0., _y), "bottom"); model.applyBC(BC::Dirichlet::FixedValue(0., _x), "left"); model.setBaseName("phase_notch"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpFieldVector("displacement"); model.addDumpField("damage"); model.dump(); UInt nbSteps = 1500; Real increment = 1e-5; auto start_time = clk::now(); for (UInt s = 1; s < nbSteps; ++s) { if (s >= 500) { increment = 1.e-6; } if (s % 10 == 0) { constexpr char wheel[] = "/-\\|"; auto elapsed = clk::now() - start_time; auto time_per_step = elapsed / s; std::cout << "\r[" << wheel[(s / 10) % 4] << "] " << std::setw(5) << s << "/" << nbSteps << " (" << std::setprecision(2) << std::fixed << std::setw(8) << millisecond(time_per_step).count() << "ms/step - elapsed: " << std::setw(8) << second(elapsed).count() << "s - ETA: " << std::setw(8) << second((nbSteps - s) * time_per_step).count() << "s)" << std::string(' ', 20) << std::flush; } model.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "top"); - + coupler.solve(); + auto energy = phase.getEnergy(); + if (s % 100 == 0) { model.dump(); } } + Real damage_limit = 0.08; + auto global_nb_clusters = + mesh.createClusters(spatial_dimension, "crack", PhaseFieldElementFilter(phase, damage_limit)); + + + auto nb_fragment = mesh.getNbElementGroups(spatial_dimension); + + model.dumpGroup("crack_0"); + + std::cout << global_nb_clusters << std::endl; + std::cout << nb_fragment << std::endl; + finalize(); return EXIT_SUCCESS; } diff --git a/src/model/phase_field/phasefield.cc b/src/model/phase_field/phasefield.cc index 207e02915..0264c9e32 100644 --- a/src/model/phase_field/phasefield.cc +++ b/src/model/phase_field/phasefield.cc @@ -1,383 +1,383 @@ /** * @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, const ID & id) : Parsable(ParserType::_phasefield, id), id(id), fem(model.getFEEngine()), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id), damage_on_qpoints("damage", *this), phi("phi", *this), strain("strain", *this), gradd("grad_d", *this), driving_force("driving_force", *this), dissipated_energy("dissipated_energy", *this), damage_energy("damage_energy", *this), damage_energy_density("damage_energy_density", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseField::PhaseField(PhaseFieldModel & model, UInt 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()), element_filter("element_filter", id), damage_on_qpoints("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), gradd("grad_d", *this, dim, fe_engine, this->element_filter), driving_force("driving_force", *this, dim, fe_engine, this->element_filter), dissipated_energy("dissipated_energy", *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() = 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_on_qpoints.initialize(1); phi.initialize(1); driving_force.initialize(1); gradd.initialize(spatial_dimension); strain.initialize(spatial_dimension * spatial_dimension); dissipated_energy.initialize(1); 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_uint.begin(); it != internal_vectors_uint.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(); UInt 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::computeDissipatedEnergyByElements() { AKANTU_DEBUG_IN(); const Array & damage = model.getDamage(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { Array & elem_filter = element_filter(type, _not_ghost); if (elem_filter.empty()) { continue; } Array & damage_interpolated = damage_on_qpoints(type, _not_ghost); // compute the damage on quadrature points fem.interpolateOnIntegrationPoints(damage, damage_interpolated, 1, type, _not_ghost); Array & gradd_vect = gradd(type, _not_ghost); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(damage, gradd_vect, - spatial_dimension, type, _not_ghost, + 1, type, _not_ghost, elem_filter); computeDissipatedEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::computeDissipatedEnergy(ElementType /*unused*/) { AKANTU_DEBUG_IN(); AKANTU_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real PhaseField::getEnergy(){ AKANTU_DEBUG_IN(); Real edis = 0.; computeDissipatedEnergyByElements(); /// integrate the dissipated energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { edis += fem.integrate(dissipated_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real PhaseField::getEnergy(ElementType type, UInt index) { Real edis = 0.; Vector edis_on_quad_points(fem.getNbIntegrationPoints(type)); computeDissipatedEnergyByElement(type, index, edis_on_quad_points); edis = fem.integrate(edis_on_quad_points, type, element_filter(type)(index)); return edis; } /* -------------------------------------------------------------------------- */ 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