diff --git a/python/py_fragment_manager.cc b/python/py_fragment_manager.cc index 21ee72939..cd2a48e9b 100644 --- a/python/py_fragment_manager.cc +++ b/python/py_fragment_manager.cc @@ -1,47 +1,50 @@ /* -------------------------------------------------------------------------- */ +#include "py_aka_array.hh" #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](FragmentManager & self) -> decltype(auto) { \ return self.func_name(); \ }) void register_fragment_manager(py::module & mod) { - py::class_(mod, "FragmentManager") + py::class_(mod, "FragmentManager") .def(py::init(), py::arg("model"), py::arg("dump_data") = true, py::arg("ID") = "fragment_manager") - .def_function(buildFragments) + .def("buildFragments", &FragmentManager::buildFragments, + py::arg("damage_limit") = 1.) .def_function(computeCenterOfMass) .def_function(computeVelocity) .def_function(computeInertiaMoments) - .def_function(computeAllData) + .def("computeAllData", &FragmentManager::computeAllData, + py::arg("damage_limit") = 1.) .def_function(computeNbElementsPerFragment) .def_function(getNbFragment) .def_function(getMass) .def_function(getVelocity) .def_function(getMomentsOfInertia) .def_function(getPrincipalDirections) .def_function(getNbElementsPerFragment); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc index a308bd7c0..ebc197151 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc @@ -1,544 +1,545 @@ /** * @file fragment_manager.cc * * @author Aurelia Isabel Cuba Ramos * @author Marco Vocialta * * @date creation: Thu Jan 23 2014 * @date last modification: Tue Feb 20 2018 * * @brief Group manager to handle fragments * * * Copyright (©) 2014-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 . * */ /* -------------------------------------------------------------------------- */ #include "fragment_manager.hh" #include "aka_iterators.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "material_cohesive.hh" #include "mesh_iterators.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ FragmentManager::FragmentManager(SolidMechanicsModelCohesive & model, bool dump_data, const ID & id) : GroupManager(model.getMesh(), id), model(model), mass_center(0, model.getSpatialDimension(), "mass_center"), mass(0, model.getSpatialDimension(), "mass"), velocity(0, model.getSpatialDimension(), "velocity"), inertia_moments(0, model.getSpatialDimension(), "inertia_moments"), principal_directions( 0, model.getSpatialDimension() * model.getSpatialDimension(), "principal_directions"), quad_coordinates("quad_coordinates", id), mass_density("mass_density", id), nb_elements_per_fragment(0, 1, "nb_elements_per_fragment"), dump_data(dump_data) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); /// compute quadrature points' coordinates quad_coordinates.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(quad_coordinates, spatial_dimension, // spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getMesh().getNodes(), quad_coordinates); /// store mass density per quadrature point mass_density.initialize(mesh, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(mass_density, 1, spatial_dimension, // _not_ghost); storeMassDensityPerIntegrationPoint(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ class CohesiveElementFilter : public GroupManager::ClusteringFilter { public: CohesiveElementFilter(const SolidMechanicsModelCohesive & model, const Real max_damage = 1.) : model(model), is_unbroken(max_damage) {} bool operator()(const Element & el) const override { if (Mesh::getKind(el.type) == _ek_regular) { return true; } const Array & mat_indexes = model.getMaterialByElement(el.type, el.ghost_type); const Array & mat_loc_num = model.getMaterialLocalNumbering(el.type, el.ghost_type); const auto & mat = static_cast( model.getMaterial(mat_indexes(el.element))); UInt el_index = mat_loc_num(el.element); UInt nb_quad_per_element = model.getFEEngine("CohesiveFEEngine") .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; } + bool operator()(const Real & x) const { return x < max_damage; } const Real max_damage; }; const SolidMechanicsModelCohesive & model; const IsUnbrokenFunctor is_unbroken; }; /* -------------------------------------------------------------------------- */ void FragmentManager::buildFragments(Real damage_limit) { AKANTU_DEBUG_IN(); if (mesh.isDistributed()) { auto & cohesive_synchronizer = model.getCohesiveSynchronizer(); cohesive_synchronizer.synchronize(model, SynchronizationTag::_smmc_damage); } auto & mesh_facets = mesh.getMeshFacets(); UInt spatial_dimension = model.getSpatialDimension(); std::string fragment_prefix("fragment"); /// generate fragments global_nb_fragment = createClusters(spatial_dimension, mesh_facets, fragment_prefix, CohesiveElementFilter(model, damage_limit)); nb_fragment = getNbElementGroups(spatial_dimension); fragment_index.resize(nb_fragment); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { auto name = std::get<0>(data).getName(); /// get fragment index std::string fragment_index_string = name.substr(fragment_prefix.size() + 1); std::get<1>(data) = std::stoul(fragment_index_string); } /// compute fragments' mass computeMass(); if (dump_data) { createDumpDataArray(fragment_index, "fragments", true); createDumpDataArray(mass, "fragments mass"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeMass() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// create a unit field per quadrature point, since to compute mass /// it's neccessary to integrate only density ElementTypeMapArray unit_field("unit_field", id); unit_field.initialize(model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); integrateFieldOnFragments(unit_field, mass); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeCenterOfMass() { AKANTU_DEBUG_IN(); /// integrate position multiplied by density integrateFieldOnFragments(quad_coordinates, mass_center); /// divide it by the fragments' mass Real * mass_storage = mass.storage(); Real * mass_center_storage = mass_center.storage(); UInt total_components = mass_center.size() * mass_center.getNbComponent(); for (UInt i = 0; i < total_components; ++i) { mass_center_storage[i] /= mass_storage[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeVelocity() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// compute velocity per quadrature point ElementTypeMapArray velocity_field("velocity_field", id); velocity_field.initialize( model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getVelocity(), velocity_field); /// integrate on fragments integrateFieldOnFragments(velocity_field, velocity); /// divide it by the fragments' mass Real * mass_storage = mass.storage(); Real * velocity_storage = velocity.storage(); UInt total_components = velocity.size() * velocity.getNbComponent(); for (UInt i = 0; i < total_components; ++i) { velocity_storage[i] /= mass_storage[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Given the distance @f$ \mathbf{r} @f$ between a quadrature point * and its center of mass, the moment of inertia is computed as \f[ * I_\mathrm{CM} = \mathrm{tr}(\mathbf{r}\mathbf{r}^\mathrm{T}) * \mathbf{I} - \mathbf{r}\mathbf{r}^\mathrm{T} \f] for more * information check Wikipedia * (http://en.wikipedia.org/wiki/Moment_of_inertia#Identities_for_a_skew-symmetric_matrix) * */ void FragmentManager::computeInertiaMoments() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); computeCenterOfMass(); /// compute local coordinates products with respect to the center of match ElementTypeMapArray moments_coords("moments_coords", id); moments_coords.initialize(model.getFEEngine(), _nb_component = spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); /// loop over fragments for (auto && data : zip(iterateElementGroups(), make_view(mass_center, spatial_dimension))) { const auto & el_list = std::get<0>(data).getElements(); auto & mass_center = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); auto & moments_coords_array = moments_coords(type); const auto & quad_coordinates_array = quad_coordinates(type); const auto & el_list_array = el_list(type); auto moments_begin = moments_coords_array.begin(spatial_dimension, spatial_dimension); auto quad_coordinates_begin = quad_coordinates_array.begin(spatial_dimension); Vector relative_coords(spatial_dimension); for (UInt el = 0; el < el_list_array.size(); ++el) { UInt global_el = el_list_array(el); /// loop over quadrature points for (UInt q = 0; q < nb_quad_per_element; ++q) { UInt global_q = global_el * nb_quad_per_element + q; Matrix moments_matrix = moments_begin[global_q]; const Vector & quad_coord_vector = quad_coordinates_begin[global_q]; /// to understand this read the documentation written just /// before this function relative_coords = quad_coord_vector; relative_coords -= mass_center; moments_matrix.outerProduct(relative_coords, relative_coords); Real trace = moments_matrix.trace(); moments_matrix *= -1.; moments_matrix += Matrix::eye(spatial_dimension, trace); } } } } /// integrate moments Array integrated_moments(global_nb_fragment, spatial_dimension * spatial_dimension); integrateFieldOnFragments(moments_coords, integrated_moments); /// compute and store principal moments inertia_moments.resize(global_nb_fragment); principal_directions.resize(global_nb_fragment); auto integrated_moments_it = integrated_moments.begin(spatial_dimension, spatial_dimension); auto inertia_moments_it = inertia_moments.begin(spatial_dimension); auto principal_directions_it = principal_directions.begin(spatial_dimension, spatial_dimension); for (UInt frag = 0; frag < global_nb_fragment; ++frag, ++integrated_moments_it, ++inertia_moments_it, ++principal_directions_it) { integrated_moments_it->eig(*inertia_moments_it, *principal_directions_it); } if (dump_data) { createDumpDataArray(inertia_moments, "moments of inertia"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void FragmentManager::computeAllData() { +void FragmentManager::computeAllData(Real damage_limit) { AKANTU_DEBUG_IN(); - buildFragments(); + buildFragments(damage_limit); computeVelocity(); computeInertiaMoments(); + computeNbElementsPerFragment(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::storeMassDensityPerIntegrationPoint() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : mesh.elementTypes(_spatial_dimension = spatial_dimension, _element_kind = _ek_regular)) { Array & mass_density_array = mass_density(type); UInt nb_element = mesh.getNbElement(type); UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); mass_density_array.resize(nb_element * nb_quad_per_element); const Array & mat_indexes = model.getMaterialByElement(type); Real * mass_density_it = mass_density_array.storage(); /// store mass_density for each element and quadrature point for (UInt el = 0; el < nb_element; ++el) { Material & mat = model.getMaterial(mat_indexes(el)); for (UInt q = 0; q < nb_quad_per_element; ++q, ++mass_density_it) { *mass_density_it = mat.getRho(); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::integrateFieldOnFragments( ElementTypeMapArray & field, Array & output) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); UInt nb_component = output.getNbComponent(); /// integration part output.resize(global_nb_fragment); output.zero(); auto output_begin = output.begin(nb_component); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & el_list = std::get<0>(data).getElements(); auto fragment_index = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); const Array & density_array = mass_density(type); Array & field_array = field(type); const Array & elements = el_list(type); /// generate array to be integrated by filtering fragment's elements Array integration_array(elements.size() * nb_quad_per_element, nb_component); auto field_array_begin = field_array.begin_reinterpret( nb_quad_per_element, nb_component, field_array.size() / nb_quad_per_element); auto density_array_begin = density_array.begin_reinterpret( nb_quad_per_element, density_array.size() / nb_quad_per_element); for (auto && data : enumerate(make_view( integration_array, nb_quad_per_element, nb_component))) { UInt global_el = elements(std::get<0>(data)); auto & int_array = std::get<1>(data); int_array = field_array_begin[global_el]; /// multiply field by density const Vector & density_vector = density_array_begin[global_el]; for (UInt i = 0; i < nb_quad_per_element; ++i) { for (UInt j = 0; j < nb_component; ++j) { int_array(i, j) *= density_vector(i); } } } /// integrate the field over the fragment Array integrated_array(elements.size(), nb_component); model.getFEEngine().integrate(integration_array, integrated_array, nb_component, type, _not_ghost, elements); /// sum over all elements and store the result Vector output_tmp(output_begin[fragment_index]); output_tmp += std::accumulate(integrated_array.begin(nb_component), integrated_array.end(nb_component), Vector(nb_component)); } } /// sum output over all processors const Communicator & comm = mesh.getCommunicator(); comm.allReduce(output, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeNbElementsPerFragment() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); nb_elements_per_fragment.resize(global_nb_fragment); nb_elements_per_fragment.zero(); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & el_list = std::get<0>(data).getElements(); auto fragment_index = std::get<1>(data); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_element = el_list(type).size(); nb_elements_per_fragment(fragment_index) += nb_element; } } /// sum values over all processors const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_elements_per_fragment, SynchronizerOperation::_sum); if (dump_data) { createDumpDataArray(nb_elements_per_fragment, "elements per fragment"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FragmentManager::createDumpDataArray(Array & data, std::string name, bool fragment_index_output) { AKANTU_DEBUG_IN(); if (data.empty()) { return; } auto & mesh_not_const = const_cast(mesh); auto && spatial_dimension = mesh.getSpatialDimension(); auto && nb_component = data.getNbComponent(); auto && data_begin = data.begin(nb_component); /// loop over fragments for (auto && data : zip(iterateElementGroups(), fragment_index)) { const auto & fragment = std::get<0>(data); auto fragment_idx = std::get<1>(data); /// loop over cluster types for (auto && type : fragment.elementTypes(spatial_dimension)) { /// init mesh data auto & mesh_data = mesh_not_const.getDataPointer( name, type, _not_ghost, nb_component); auto mesh_data_begin = mesh_data.begin(nb_component); /// fill mesh data for (const auto & elem : fragment.getElements(type)) { Vector md_tmp = mesh_data_begin[elem]; if (fragment_index_output) { md_tmp(0) = fragment_idx; } else { md_tmp = data_begin[fragment_idx]; } } } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh index 0deb462da..3ef723a11 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh @@ -1,167 +1,167 @@ /** * @file fragment_manager.hh * * @author Marco Vocialta * * @date creation: Thu Jan 23 2014 * @date last modification: Thu Jul 06 2017 * * @brief Group manager to handle fragments * * * Copyright (©) 2014-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 . * */ /* -------------------------------------------------------------------------- */ #include "group_manager.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_FRAGMENT_MANAGER_HH_ #define AKANTU_FRAGMENT_MANAGER_HH_ namespace akantu { class SolidMechanicsModelCohesive; } namespace akantu { /* -------------------------------------------------------------------------- */ class FragmentManager : public GroupManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FragmentManager(SolidMechanicsModelCohesive & model, bool dump_data = true, const ID & id = "fragment_manager"); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ private: /// store mass density per integration point void storeMassDensityPerIntegrationPoint(); /// integrate an elemental field multiplied by density on global /// fragments void integrateFieldOnFragments(ElementTypeMapArray & field, Array & output); /// compute fragments' mass void computeMass(); /// create dump data for a single array template void createDumpDataArray(Array & data, std::string name, bool fragment_index_output = false); public: /// build fragment list (cohesive elements are considered broken if /// damage >= damage_limit) void buildFragments(Real damage_limit = 1.); /// compute fragments' center of mass void computeCenterOfMass(); /// compute fragments' velocity void computeVelocity(); /// computes principal moments of inertia with respect to the center /// of mass of each fragment void computeInertiaMoments(); /// compute all fragments' data - void computeAllData(); + void computeAllData(Real damage_limit = 1.); /// compute number of elements per fragment void computeNbElementsPerFragment(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get number of fragments AKANTU_GET_MACRO(NbFragment, global_nb_fragment, UInt); /// get fragments' mass AKANTU_GET_MACRO(Mass, mass, const Array &); /// get fragments' center of mass AKANTU_GET_MACRO(CenterOfMass, mass_center, const Array &); /// get fragments' velocity AKANTU_GET_MACRO(Velocity, velocity, const Array &); /// get fragments' principal moments of inertia AKANTU_GET_MACRO(MomentsOfInertia, inertia_moments, const Array &); /// get fragments' principal directions AKANTU_GET_MACRO(PrincipalDirections, principal_directions, const Array &); /// get number of elements per fragment AKANTU_GET_MACRO(NbElementsPerFragment, nb_elements_per_fragment, const Array &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// local_fragment index Array fragment_index; /// global number of fragments (parallel simulations) UInt global_nb_fragment; /// number of fragments UInt nb_fragment; /// cohesive solid mechanics model associated SolidMechanicsModelCohesive & model; /// fragments' center of mass Array mass_center; /// fragments' mass Array mass; /// fragments' velocity Array velocity; /// fragments' principal moments of inertia with respect to the /// center of mass Array inertia_moments; /// fragments' principal directions Array principal_directions; /// quadrature points' coordinates ElementTypeMapArray quad_coordinates; /// mass density per quadrature point ElementTypeMapArray mass_density; /// fragment filter Array nb_elements_per_fragment; /// dump data bool dump_data; }; } // namespace akantu #endif /* AKANTU_FRAGMENT_MANAGER_HH_ */