diff --git a/src/io/dumper/dumpable_inline_impl.hh b/src/io/dumper/dumpable_inline_impl.hh index 2597ebb59..5b1c4e1da 100644 --- a/src/io/dumper/dumpable_inline_impl.hh +++ b/src/io/dumper/dumpable_inline_impl.hh @@ -1,369 +1,369 @@ /** * @file dumpable_inline_impl.hh * * @author Nicolas Richart * * @date Sat Jul 13 11:35:22 2013 * * @brief Implementation of the Dumpable class * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline Dumpable::Dumpable() : default_dumper("") { } /* -------------------------------------------------------------------------- */ inline Dumpable::~Dumpable() { - DumperMap::iterator it = this->dumpers.begin(); - DumperMap::iterator end = this->dumpers.end(); + // DumperMap::iterator it = this->dumpers.begin(); + // DumperMap::iterator end = this->dumpers.end(); // for (; it != end; ++it) { // DumperSet::iterator fit = this->external_dumpers.find(it->first); // DumperSet::iterator fend = this->external_dumpers.end(); // if (fit == fend) // delete it->second; // } } /* -------------------------------------------------------------------------- */ template inline void Dumpable::registerDumper(const std::string & dumper_name, const std::string & file_name, const bool is_default) { AKANTU_DEBUG_ASSERT(this->dumpers.find(dumper_name) == this->dumpers.end(), "Dumper " + dumper_name + "is already registered."); std::string name = file_name; if (name == "") name = dumper_name; this->dumpers[dumper_name] = new T(name); if (is_default) this->default_dumper = dumper_name; } /* -------------------------------------------------------------------------- */ inline void Dumpable::registerExternalDumper(DumperIOHelper & dumper, const std::string & dumper_name, const bool is_default) { this->dumpers[dumper_name] = &dumper; if (is_default) this->default_dumper = dumper_name; } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpMesh(const Mesh & mesh, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { this->addDumpMeshToDumper(this->default_dumper, mesh, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpMeshToDumper(const std::string & dumper_name, const Mesh & mesh, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerMesh(mesh, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFilteredMesh(const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { this->addDumpFilteredMeshToDumper(this->default_dumper, mesh, elements_filter, nodes_filter, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFilteredMeshToDumper(const std::string & dumper_name, const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerFilteredMesh(mesh, elements_filter, nodes_filter, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpField(const std::string & field_id) { this->addDumpFieldToDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_TO_IMPLEMENT(); }; /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFieldExternal(const std::string & field_id, DumperIOHelper::Field * field) { this->addDumpFieldExternalToDumper(this->default_dumper, field_id, field); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFieldExternalToDumper(const std::string & dumper_name, const std::string & field_id, DumperIOHelper::Field * field) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field); } /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternal(const std::string & field_id, const Array & field) { this->addDumpFieldExternalToDumper(this->default_dumper, field_id, field); }; /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternalToDumper(const std::string & dumper_name, const std::string & field_id, const Array & field) { DumperIOHelper::Field * field_cont = new DumperIOHelper::NodalField(field); DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field_cont); } /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternal(const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { this->addDumpFieldExternalToDumper(this->default_dumper, field_id, field, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternalToDumper(const std::string & dumper_name, const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { DumperIOHelper::Field * field_cont = new DumperIOHelper::ElementalField(field, spatial_dimension, ghost_type, element_kind); DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field_cont); } /* -------------------------------------------------------------------------- */ inline void Dumpable::removeDumpField(const std::string & field_id) { this->removeDumpFieldFromDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ inline void Dumpable::removeDumpFieldFromDumper(const std::string & dumper_name, const std::string & field_id) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.unRegisterField(field_id); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFieldVector(const std::string & field_id) { this->addDumpFieldVectorToDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFieldTensor(const std::string & field_id) { this->addDumpFieldTensorToDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ inline void Dumpable::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ inline void Dumpable::setDirectory(const std::string & directory) { this->setDirectoryToDumper(this->default_dumper, directory); }; /* -------------------------------------------------------------------------- */ inline void Dumpable::setDirectoryToDumper(const std::string & dumper_name, const std::string & directory) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.setDirectory(directory); } /* -------------------------------------------------------------------------- */ inline void Dumpable::setBaseName(const std::string & basename) { this->setBaseNameToDumper(this->default_dumper, basename); } /* -------------------------------------------------------------------------- */ inline void Dumpable::setBaseNameToDumper(const std::string & dumper_name, const std::string & basename) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.setBaseName(basename); } /* -------------------------------------------------------------------------- */ inline void Dumpable::setTimeStepToDumper(Real time_step) { this->setTimeStepToDumper(this->default_dumper, time_step); } /* -------------------------------------------------------------------------- */ inline void Dumpable::setTimeStepToDumper(const std::string & dumper_name, Real time_step) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.setTimeStep(time_step); }; /* -------------------------------------------------------------------------- */ inline void Dumpable::dump(const std::string & dumper_name) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.dump(); } /* -------------------------------------------------------------------------- */ inline void Dumpable::dump() { this->dump(this->default_dumper); } /* -------------------------------------------------------------------------- */ inline void Dumpable::dump(const std::string & dumper_name, UInt step) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.dump(step); } /* -------------------------------------------------------------------------- */ inline void Dumpable::dump(UInt step) { this->dump(this->default_dumper, step); } /* -------------------------------------------------------------------------- */ inline void Dumpable::dump(const std::string & dumper_name, Real time, UInt step) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.dump(time, step); } /* -------------------------------------------------------------------------- */ inline void Dumpable::dump(Real time, UInt step) { this->dump(this->default_dumper, time, step); } /* -------------------------------------------------------------------------- */ inline void Dumpable::internalAddDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id, DumperIOHelper::Field * field) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field); } /* -------------------------------------------------------------------------- */ inline DumperIOHelper & Dumpable::getDumper() { return this->getDumper(this->default_dumper); } /* -------------------------------------------------------------------------- */ inline DumperIOHelper & Dumpable::getDumper(const std::string & dumper_name) { DumperMap::iterator it = this->dumpers.find(dumper_name); DumperMap::iterator end = this->dumpers.end(); if (it == end) AKANTU_EXCEPTION("Dumper " << dumper_name << "has not been registered, yet."); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline T & Dumpable::getDumper(const std::string & dumper_name) { DumperIOHelper & dumper = this->getDumper(dumper_name); try { T & templated_dumper = dynamic_cast(dumper); return templated_dumper; } catch (...) { AKANTU_EXCEPTION("Dumper " << dumper_name << " is not of type: " << debug::demangle(typeid(T).name())); } } /* -------------------------------------------------------------------------- */ inline std::string Dumpable::getDefaultDumperName() const { return this->default_dumper; } __END_AKANTU__ diff --git a/src/model/solid_mechanics/fragment_manager.cc b/src/model/solid_mechanics/fragment_manager.cc index e666e9dac..027ca9ea1 100644 --- a/src/model/solid_mechanics/fragment_manager.cc +++ b/src/model/solid_mechanics/fragment_manager.cc @@ -1,646 +1,654 @@ /** * @file fragment_manager.cc * @author Marco Vocialta * @date Mon Jan 20 14:33:49 2014 * * @brief Group manager to handle fragments * * @section LICENSE * * Copyright (©) 2010-2011 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 "material_cohesive.hh" #include +#include +#include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ FragmentManager::FragmentManager(SolidMechanicsModelCohesive & model, bool dump_data, const ID & id, const MemoryID & memory_id) : GroupManager(model.getMesh(), id, memory_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"), quad_coordinates("quad_coordinates", id), mass_density("mass_density", id), fragment_filter(0, 1, "fragment_filter"), dump_data(dump_data) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); /// compute quadrature points' coordinates mesh.initElementTypeMapArray(quad_coordinates, spatial_dimension, spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnQuadraturePoints(model.getMesh().getNodes(), quad_coordinates); /// store mass density per quadrature point mesh.initElementTypeMapArray(mass_density, 1, spatial_dimension, _not_ghost); storeMassDensityPerQuadraturePoint(); 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 { if (el.kind == _ek_regular) return true; const Array & el_id_by_mat = model.getElementIndexByMaterial(el.type, el.ghost_type); const MaterialCohesive & mat = static_cast (model.getMaterial(el_id_by_mat(el.element, 0))); UInt el_index = el_id_by_mat(el.element, 1); UInt nb_quad_per_element = model.getFEEngine("CohesiveFEEngine").getNbQuadraturePoints(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.getSize(), "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); if (unbroken_quads > 0) return true; return false; } private: struct IsUnbrokenFunctor { IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {} bool operator() (const Real & x) {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 defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) DistributedSynchronizer * cohesive_synchronizer = const_cast(model.getCohesiveSynchronizer()); if (cohesive_synchronizer) { cohesive_synchronizer->computeBufferSize(model, _gst_smmc_damage); cohesive_synchronizer->asynchronousSynchronize(model, _gst_smmc_damage); cohesive_synchronizer->waitEndSynchronize(model, _gst_smmc_damage); } #endif DistributedSynchronizer & synchronizer = const_cast(model.getSynchronizer()); Mesh & mesh_facets = const_cast(model.getMeshFacets()); UInt spatial_dimension = model.getSpatialDimension(); std::string fragment_prefix("fragment"); /// generate fragments global_nb_fragment = createClusters(spatial_dimension, fragment_prefix, CohesiveElementFilter(model, damage_limit), &synchronizer, &mesh_facets); nb_fragment = getNbElementGroups(spatial_dimension); fragment_index.resize(nb_fragment); UInt * fragment_index_it = fragment_index.storage(); /// loop over fragments for(const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++fragment_index_it) { /// get fragment index std::string fragment_index_string = it->first.substr(fragment_prefix.size() + 1); std::stringstream sstr(fragment_index_string.c_str()); sstr >> *fragment_index_it; AKANTU_DEBUG_ASSERT(!sstr.fail(), "fragment_index is not an integer"); } /// 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); mesh.initElementTypeMapArray(unit_field, spatial_dimension, spatial_dimension, _not_ghost); ElementTypeMapArray::type_iterator it = unit_field.firstType(spatial_dimension, _not_ghost, _ek_regular); ElementTypeMapArray::type_iterator end = unit_field.lastType(spatial_dimension, _not_ghost, _ek_regular); for (; it != end; ++it) { ElementType type = *it; Array & field_array = unit_field(type); UInt nb_element = mesh.getNbElement(type); UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type); field_array.resize(nb_element * nb_quad_per_element); field_array.set(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.getSize() * 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); mesh.initElementTypeMapArray(velocity_field, spatial_dimension, spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnQuadraturePoints(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.getSize() * 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); mesh.initElementTypeMapArray(moments_coords, spatial_dimension * spatial_dimension, spatial_dimension, _not_ghost); /// resize the by element type ElementTypeMapArray::type_iterator it = moments_coords.firstType(spatial_dimension, _not_ghost, _ek_regular); ElementTypeMapArray::type_iterator end = moments_coords.lastType(spatial_dimension, _not_ghost, _ek_regular); for (; it != end; ++it) { ElementType type = *it; Array & field_array = moments_coords(type); UInt nb_element = mesh.getNbElement(type); UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type); field_array.resize(nb_element * nb_quad_per_element); } /// compute coordinates Array::const_iterator< Vector > mass_center_it = mass_center.begin(spatial_dimension); /// loop over fragments for(const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++mass_center_it) { const ElementTypeMapArray & el_list = it->second->getElements(); ElementTypeMapArray::type_iterator type_it = el_list.firstType(spatial_dimension, _not_ghost, _ek_regular); ElementTypeMapArray::type_iterator type_end = el_list.lastType(spatial_dimension, _not_ghost, _ek_regular); /// loop over elements of the fragment for (; type_it != type_end; ++type_it) { ElementType type = *type_it; UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type); Array & moments_coords_array = moments_coords(type); const Array & quad_coordinates_array = quad_coordinates(type); const Array & el_list_array = el_list(type); Array::iterator< Matrix > moments_begin = moments_coords_array.begin(spatial_dimension, spatial_dimension); Array::const_iterator< Vector > quad_coordinates_begin = quad_coordinates_array.begin(spatial_dimension); Vector relative_coords(spatial_dimension); for (UInt el = 0; el < el_list_array.getSize(); ++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_it; 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); Array::iterator< Matrix > integrated_moments_it = integrated_moments.begin(spatial_dimension, spatial_dimension); Array::iterator< Vector > inertia_moments_it = inertia_moments.begin(spatial_dimension); for (UInt frag = 0; frag < global_nb_fragment; ++frag, - ++integrated_moments_it, ++inertia_moments_it) + ++integrated_moments_it, ++inertia_moments_it) { integrated_moments_it->eig(*inertia_moments_it); + // sort moments + std::sort(inertia_moments_it->storage(), + inertia_moments_it->storage() + spatial_dimension, + std::greater()); + } + if (dump_data) createDumpDataArray(inertia_moments, "moments of inertia"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeAllData() { AKANTU_DEBUG_IN(); buildFragments(); computeVelocity(); computeInertiaMoments(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::storeMassDensityPerQuadraturePoint() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator end = mesh.lastType(spatial_dimension); for (; it != end; ++it) { ElementType type = *it; Array & mass_density_array = mass_density(type); UInt nb_element = mesh.getNbElement(type); UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type); mass_density_array.resize(nb_element * nb_quad_per_element); Array & el_index_by_mat = model.getElementIndexByMaterial(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(el_index_by_mat(el, 0)); 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.clear(); UInt * fragment_index_it = fragment_index.storage(); Array::iterator< Vector > output_begin = output.begin(nb_component); /// loop over fragments for(const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++fragment_index_it) { const ElementTypeMapArray & el_list = it->second->getElements(); ElementTypeMapArray::type_iterator type_it = el_list.firstType(spatial_dimension, _not_ghost, _ek_regular); ElementTypeMapArray::type_iterator type_end = el_list.lastType(spatial_dimension, _not_ghost, _ek_regular); /// loop over elements of the fragment for (; type_it != type_end; ++type_it) { ElementType type = *type_it; UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type); const Array & density_array = mass_density(type); Array & field_array = field(type); const Array & elements = el_list(type); UInt nb_element = elements.getSize(); /// generate array to be integrated by filtering fragment's elements Array integration_array(elements.getSize() * nb_quad_per_element, nb_component); Array::iterator< Matrix > int_array_it = integration_array.begin_reinterpret(nb_quad_per_element, nb_component, nb_element); Array::iterator< Matrix > int_array_end = integration_array.end_reinterpret(nb_quad_per_element, nb_component, nb_element); Array::iterator< Matrix > field_array_begin = field_array.begin_reinterpret(nb_quad_per_element, nb_component, field_array.getSize() / nb_quad_per_element); Array::const_iterator< Vector > density_array_begin = density_array.begin_reinterpret(nb_quad_per_element, density_array.getSize() / nb_quad_per_element); for (UInt el = 0; int_array_it != int_array_end; ++int_array_it, ++el) { UInt global_el = elements(el); *int_array_it = 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_it)(i, j) *= density_vector(i); } } } /// integrate the field over the fragment Array integrated_array(elements.getSize(), nb_component); model.getFEEngine().integrate(integration_array, integrated_array, nb_component, type, _not_ghost, elements); /// sum over all elements and store the result output_begin[*fragment_index_it] += std::accumulate(integrated_array.begin(nb_component), integrated_array.end(nb_component), Vector(nb_component)); } } /// sum output over all processors StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); comm.allReduce(output.storage(), global_nb_fragment * nb_component, _so_sum); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt FragmentManager::getNbBigFragments(UInt minimum_nb_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); Array nb_element_per_fragment(global_nb_fragment); nb_element_per_fragment.clear(); UInt * fragment_index_it = fragment_index.storage(); /// loop over fragments for(const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++fragment_index_it) { const ElementTypeMapArray & el_list = it->second->getElements(); ElementTypeMapArray::type_iterator type_it = el_list.firstType(spatial_dimension, _not_ghost, _ek_regular); ElementTypeMapArray::type_iterator type_end = el_list.lastType(spatial_dimension, _not_ghost, _ek_regular); /// loop over elements of the fragment for (; type_it != type_end; ++type_it) { ElementType type = *type_it; UInt nb_element = el_list(type).getSize(); nb_element_per_fragment(*fragment_index_it) += nb_element; } } /// sum values over all processors StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); comm.allReduce(nb_element_per_fragment.storage(), global_nb_fragment, _so_sum); /// count big enough fragments fragment_filter.resize(global_nb_fragment); fragment_filter.clear(); UInt nb_big_fragment = 0; for (UInt frag = 0; frag < global_nb_fragment; ++frag) { if (nb_element_per_fragment(frag) >= minimum_nb_elements) { fragment_filter(frag) = true; ++nb_big_fragment; } } if (dump_data) createDumpDataArray(nb_element_per_fragment, "elements per fragment"); AKANTU_DEBUG_OUT(); return nb_big_fragment; } /* -------------------------------------------------------------------------- */ template void FragmentManager::createDumpDataArray(Array & data, std::string name, bool fragment_index_output) { AKANTU_DEBUG_IN(); if (data.getSize() == 0) return; typedef typename Array::template iterator< Vector > data_iterator; Mesh & mesh_not_const = const_cast(mesh); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_component = data.getNbComponent(); UInt * fragment_index_it = fragment_index.storage(); data_iterator data_begin = data.begin(nb_component); /// loop over fragments for(const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++fragment_index_it) { const ElementGroup & fragment = *(it->second); /// loop over cluster types ElementGroup::type_iterator type_it = fragment.firstType(spatial_dimension); ElementGroup::type_iterator type_end = fragment.lastType(spatial_dimension); for(; type_it != type_end; ++type_it) { ElementType type = *type_it; /// init mesh data Array * mesh_data = mesh_not_const.getDataPointer(name, type, _not_ghost, nb_component); data_iterator mesh_data_begin = mesh_data->begin(nb_component); ElementGroup::const_element_iterator el_it = fragment.element_begin(type); ElementGroup::const_element_iterator el_end = fragment.element_end(type); /// fill mesh data if (fragment_index_output) { for (; el_it != el_end; ++el_it) mesh_data_begin[*el_it](0) = *fragment_index_it; } else { for (; el_it != el_end; ++el_it) mesh_data_begin[*el_it] = data_begin[*fragment_index_it]; } } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc index e24f88e9e..0931e6a79 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc @@ -1,416 +1,426 @@ /** * @file test_cohesive_parallel_buildfragments.cc * @author Marco Vocialta * @date Tue Dec 17 12:03:03 2013 * * @brief Test to build fragments in parallel * * @section LICENSE * * Copyright (©) 2010-2011 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 #include #include +#include +#include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" #include "material_cohesive.hh" #include "fragment_manager.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; void verticalInsertionLimit(SolidMechanicsModelCohesive &); void displaceElements(SolidMechanicsModelCohesive &, const Real, const Real); bool isInertiaEqual(const Vector &, const Vector &); void rotateArray(Array & array, Real angle); const UInt spatial_dimension = 3; const UInt total_nb_fragment = 4; const Real rotation_angle = M_PI / 4.; const Real global_tolerance = 1.e-9; int main(int argc, char *argv[]) { initialize("material.dat", argc, argv); Math::setTolerance(global_tolerance); Mesh mesh(spatial_dimension); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); akantu::MeshPartition * partition = NULL; if(prank == 0) { // Read the mesh mesh.read("mesh.msh"); /// partition the mesh MeshUtils::purifyMesh(mesh); partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); } SolidMechanicsModelCohesive model(mesh); model.initParallel(partition, NULL, true); delete partition; /// model initialization model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); mesh.computeBoundingBox(); Real L = mesh.getUpperBounds()(0) - mesh.getLowerBounds()(0); Real h = mesh.getUpperBounds()(1) - mesh.getLowerBounds()(1); Real rho = model.getMaterial("bulk").getParam("rho"); Real theoretical_mass = L * h * h * rho; Real frag_theo_mass = theoretical_mass / total_nb_fragment; UInt nb_element = mesh.getNbElement(spatial_dimension, _not_ghost, _ek_regular); comm.allReduce(&nb_element, 1, _so_sum); UInt nb_element_per_fragment = nb_element / total_nb_fragment; FragmentManager fragment_manager(model); fragment_manager.computeAllData(); fragment_manager.getNbBigFragments(nb_element_per_fragment + 1); model.setBaseName("extrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("stress"); model.addDumpField("partitions"); model.addDumpField("fragments"); model.addDumpField("fragments mass"); model.addDumpField("moments of inertia"); model.addDumpField("elements per fragment"); model.dump(); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_test"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump("cohesive elements"); /// set check facets verticalInsertionLimit(model); model.assembleMassLumped(); model.synchronizeBoundaries(); /// impose initial displacement Array & displacement = model.getDisplacement(); Array & velocity = model.getVelocity(); const Array & position = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { displacement(n, 0) = position(n, 0) * 0.1; velocity(n, 0) = position(n, 0); } rotateArray(mesh.getNodes(), rotation_angle); // rotateArray(displacement, rotation_angle); // rotateArray(velocity, rotation_angle); model.updateResidual(); model.checkCohesiveStress(); model.dump(); model.dump("cohesive elements"); const Array & fragment_mass = fragment_manager.getMass(); const Array & fragment_center = fragment_manager.getCenterOfMass(); Real el_size = L / total_nb_fragment; Real lim = -L/2 + el_size * 0.99; /// define theoretical inertia moments Vector small_frag_inertia(spatial_dimension); small_frag_inertia(0) = frag_theo_mass * (h*h + h*h) / 12.; small_frag_inertia(1) = frag_theo_mass * (el_size*el_size + h*h) / 12.; small_frag_inertia(2) = frag_theo_mass * (el_size*el_size + h*h) / 12.; + std::sort(small_frag_inertia.storage(), + small_frag_inertia.storage() + spatial_dimension, + std::greater()); + const Array & inertia_moments = fragment_manager.getMomentsOfInertia(); Array::const_iterator< Vector > inertia_moments_begin = inertia_moments.begin(spatial_dimension); /// displace one fragment each time for (UInt frag = 1; frag <= total_nb_fragment; ++frag) { if (prank == 0) std::cout << "Generating fragment: " << frag << std::endl; fragment_manager.computeAllData(); /// check number of big fragments UInt nb_big_fragment = fragment_manager.getNbBigFragments(nb_element_per_fragment + 1); model.dump(); model.dump("cohesive elements"); if (frag < total_nb_fragment) { if (nb_big_fragment != 1) { AKANTU_DEBUG_ERROR("The number of big fragments is wrong: " << nb_big_fragment); return EXIT_FAILURE; } } else { if (nb_big_fragment != 0) { AKANTU_DEBUG_ERROR("The number of big fragments is wrong: " << nb_big_fragment); return EXIT_FAILURE; } } /// check number of fragments UInt nb_fragment_num = fragment_manager.getNbFragment(); if (nb_fragment_num != frag) { AKANTU_DEBUG_ERROR("The number of fragments is wrong! Numerical: " << nb_fragment_num << " Theoretical: " << frag); return EXIT_FAILURE; } /// check mass computation if (frag < total_nb_fragment) { Real total_mass = 0.; UInt small_fragments = 0; for (UInt f = 0; f < nb_fragment_num; ++f) { const Vector & current_inertia = inertia_moments_begin[f]; if (Math::are_float_equal(fragment_mass(f, 0), frag_theo_mass)) { /// check center of mass if (fragment_center(f, 0) > (L * frag / total_nb_fragment - L / 2)) { AKANTU_DEBUG_ERROR("Fragment center is wrong!"); return EXIT_FAILURE; } /// check moment of inertia if (!isInertiaEqual(current_inertia, small_frag_inertia)) { AKANTU_DEBUG_ERROR("Inertia moments are wrong"); return EXIT_FAILURE; } ++small_fragments; total_mass += frag_theo_mass; } else { /// check the moment of inertia for the biggest fragment Real big_frag_mass = frag_theo_mass * (total_nb_fragment - frag + 1); Real big_frag_size = el_size * (total_nb_fragment - frag + 1); Vector big_frag_inertia(spatial_dimension); big_frag_inertia(0) = big_frag_mass * (h*h + h*h) / 12.; big_frag_inertia(1) = big_frag_mass * (big_frag_size*big_frag_size + h*h) / 12.; big_frag_inertia(2) = big_frag_mass * (big_frag_size*big_frag_size + h*h) / 12.; + std::sort(big_frag_inertia.storage(), + big_frag_inertia.storage() + spatial_dimension, + std::greater()); + if (!isInertiaEqual(current_inertia, big_frag_inertia)) { AKANTU_DEBUG_ERROR("Inertia moments are wrong"); return EXIT_FAILURE; } } } if (small_fragments != nb_fragment_num - 1) { AKANTU_DEBUG_ERROR("The number of small fragments is wrong!"); return EXIT_FAILURE; } if (!Math::are_float_equal(total_mass, small_fragments * frag_theo_mass)) { AKANTU_DEBUG_ERROR("The mass of small fragments is wrong!"); return EXIT_FAILURE; } } /// displace fragments rotateArray(mesh.getNodes(), -rotation_angle); // rotateArray(displacement, -rotation_angle); displaceElements(model, lim, el_size * 2); rotateArray(mesh.getNodes(), rotation_angle); // rotateArray(displacement, rotation_angle); model.updateResidual(); lim += el_size; } model.dump(); model.dump("cohesive elements"); /// check centers const Array & fragment_velocity = fragment_manager.getVelocity(); Real initial_position = -L / 2. + el_size / 2.; for (UInt frag = 0; frag < total_nb_fragment; ++frag) { Real theoretical_center = initial_position + el_size * frag; UInt f_index = 0; while (f_index < total_nb_fragment && !Math::are_float_equal(fragment_center(f_index, 0), theoretical_center)) ++f_index; if (f_index == total_nb_fragment) { AKANTU_DEBUG_ERROR("The fragments' center is wrong!"); return EXIT_FAILURE; } f_index = 0; while (f_index < total_nb_fragment && !Math::are_float_equal(fragment_velocity(f_index, 0), theoretical_center)) ++f_index; if (f_index == total_nb_fragment) { AKANTU_DEBUG_ERROR("The fragments' velocity is wrong!"); return EXIT_FAILURE; } } finalize(); if (prank == 0) std::cout << "OK: test_cohesive_buildfragments was passed!" << std::endl; return EXIT_SUCCESS; } void verticalInsertionLimit(SolidMechanicsModelCohesive & model) { UInt spatial_dimension = model.getSpatialDimension(); const Mesh & mesh_facets = model.getMeshFacets(); const Array & position = mesh_facets.getNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type); for(; it != end; ++it) { ElementType type = *it; Array & check_facets = model.getElementInserter().getCheckFacets(type, ghost_type); const Array & connectivity = mesh_facets.getConnectivity(type, ghost_type); UInt nb_nodes_per_facet = connectivity.getNbComponent(); for (UInt f = 0; f < check_facets.getSize(); ++f) { if (!check_facets(f)) continue; UInt nb_aligned_nodes = 1; Real first_node_pos = position(connectivity(f, 0), 0); for (; nb_aligned_nodes < nb_nodes_per_facet; ++nb_aligned_nodes) { Real other_node_pos = position(connectivity(f, nb_aligned_nodes), 0); if (! Math::are_float_equal(first_node_pos, other_node_pos)) break; } if (nb_aligned_nodes != nb_nodes_per_facet) { check_facets(f) = false; } } } } } void displaceElements(SolidMechanicsModelCohesive & model, const Real lim, const Real amount) { UInt spatial_dimension = model.getSpatialDimension(); Array & displacement = model.getDisplacement(); Mesh & mesh = model.getMesh(); UInt nb_nodes = mesh.getNbNodes(); Array displaced(nb_nodes); displaced.clear(); Vector barycenter(spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { ElementType type = *it; const Array & connectivity = mesh.getConnectivity(type, ghost_type); UInt nb_element = connectivity.getSize(); UInt nb_nodes_per_element = connectivity.getNbComponent(); Array::const_vector_iterator conn_el = connectivity.begin(nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter(el, type, barycenter.storage(), ghost_type); if (barycenter(0) < lim) { const Vector & conn = conn_el[el]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = conn(n); if (!displaced(node)) { displacement(node, 0) -= amount; displaced(node) = true; } } } } } } } bool isInertiaEqual(const Vector & a, const Vector & b) { UInt nb_terms = a.size(); UInt equal_terms = 0; while (equal_terms < nb_terms && std::abs(a(equal_terms) - b(equal_terms)) / a(equal_terms) < Math::getTolerance()) ++equal_terms; return equal_terms == nb_terms; } void rotateArray(Array & array, Real angle) { UInt spatial_dimension = array.getNbComponent(); Real rotation_values[] = {std::cos(angle), std::sin(angle), 0, -std::sin(angle), std::cos(angle), 0, 0, 0, 1}; Matrix rotation(rotation_values, spatial_dimension, spatial_dimension); RVector displaced_node(spatial_dimension); Array::vector_iterator node_it = array.begin(spatial_dimension); Array::vector_iterator node_end = array.end(spatial_dimension); for (; node_it != node_end; ++node_it) { displaced_node.mul(rotation, *node_it); *node_it = displaced_node; } }