diff --git a/src/fem/mesh_inline_impl.cc b/src/fem/mesh_inline_impl.cc index 9e47d1061..eba644049 100644 --- a/src/fem/mesh_inline_impl.cc +++ b/src/fem/mesh_inline_impl.cc @@ -1,531 +1,535 @@ /** * @file mesh_inline_impl.cc * * @author Guillaume Anciaux * @author Marco Vocialta * @author Nicolas Richart * * @date Thu Jul 15 00:41:12 2010 * * @brief Implementation of the inline functions of the mesh 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 . * */ __END_AKANTU__ #if defined(AKANTU_COHESIVE_ELEMENT) # include "cohesive_element.hh" #endif #include "static_communicator.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh) : new_numbering(mesh.getNbNodes(), 1, "new_numbering") { } /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh) : new_numbering("new_numbering", mesh.getID()) { } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedElementsEvent & event) { connectivities.onElementsRemoved(event.getNewNumbering()); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedNodesEvent & event) { if(created_nodes) removeNodesFromArray(*nodes , event.getNewNumbering()); if(nodes_global_ids) removeNodesFromArray(*nodes_global_ids, event.getNewNumbering()); if(nodes_type) removeNodesFromArray(*nodes_type , event.getNewNumbering()); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template inline void Mesh::removeNodesFromArray(Array & vect, const Array & new_numbering) { Array tmp(vect.getSize(), vect.getNbComponent()); UInt nb_component = vect.getNbComponent(); UInt new_nb_nodes = 0; for (UInt i = 0; i < new_numbering.getSize(); ++i) { UInt new_i = new_numbering(i); if(new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i * nb_component, nb_component * sizeof(T)); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::elementToLinearized(const Element & elem) const { AKANTU_DEBUG_ASSERT(elem.type < _max_element_type && elem.element < types_offsets.values[elem.type+1], "The element " << elem << "does not exists in the mesh " << id); return types_offsets.values[elem.type] + elem.element; } /* -------------------------------------------------------------------------- */ inline Element Mesh::linearizedToElement (UInt linearized_element) const { UInt t; for (t = _not_defined; t != _max_element_type && linearized_element >= types_offsets(t); ++t); AKANTU_DEBUG_ASSERT(linearized_element < types_offsets(t), "The linearized element " << linearized_element << "does not exists in the mesh " << id); --t; ElementType type = ElementType(t); return Element(type, linearized_element - types_offsets.values[t], _not_ghost, getKind(type)); } /* -------------------------------------------------------------------------- */ inline void Mesh::updateTypesOffsets(const GhostType & ghost_type) { + Array * types_offsets_ptr = &this->types_offsets; + if(ghost_type == _ghost) types_offsets_ptr = &this->ghost_types_offsets; + Array & types_offsets = *types_offsets_ptr; + types_offsets.clear(); type_iterator it = firstType(_all_dimensions, ghost_type, _ek_not_defined); type_iterator last = lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != last; ++it) types_offsets(*it) = connectivities(*it, ghost_type).getSize(); for (UInt t = _not_defined + 1; t < _max_element_type; ++t) types_offsets(t) += types_offsets(t - 1); for (UInt t = _max_element_type; t > _not_defined; --t) types_offsets(t) = types_offsets(t - 1); types_offsets(0) = 0; } /* -------------------------------------------------------------------------- */ inline const Mesh::ConnectivityTypeList & Mesh::getConnectivityTypeList(const GhostType & ghost_type) const { if (ghost_type == _not_ghost) return type_set; else return ghost_type_set; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbBoundaries() const { return boundaries.getNbBoundaries(); } /* -------------------------------------------------------------------------- */ inline Array * Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if(nodes_global_ids == NULL) { std::stringstream sstr; sstr << id << ":nodes_global_ids"; nodes_global_ids = &(alloc(sstr.str(), nodes->getSize(), 1)); } AKANTU_DEBUG_OUT(); return nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Array * Mesh::getNodesTypePointer() { AKANTU_DEBUG_IN(); if(nodes_type == NULL) { std::stringstream sstr; sstr << id << ":nodes_type"; nodes_type = &(alloc(sstr.str(), nodes->getSize(), 1, -1)); } AKANTU_DEBUG_OUT(); return nodes_type; } /* -------------------------------------------------------------------------- */ inline Array * Mesh::getConnectivityPointer(const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Array * tmp; if(!connectivities.exists(type, ghost_type)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); tmp = &(connectivities.alloc(0, nb_nodes_per_element, type, ghost_type)); AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); if (ghost_type == _not_ghost) type_set.insert(type); else ghost_type_set.insert(type); updateTypesOffsets(ghost_type); } else { tmp = &connectivities(type, ghost_type); } AKANTU_DEBUG_OUT(); return tmp; } /* -------------------------------------------------------------------------- */ inline Array< std::vector > * Mesh::getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type) { return getDataPointer< std::vector >(type, "element_to_subelement", ghost_type); } /* -------------------------------------------------------------------------- */ inline Array * Mesh::getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type) { return getDataPointer(type, "subelement_to_element", ghost_type, getNbFacetsPerElement(type)); } /* -------------------------------------------------------------------------- */ inline const Array< std::vector > & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) const { return getData< std::vector >(type, "element_to_subelement", ghost_type); } /* -------------------------------------------------------------------------- */ inline Array< std::vector > & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) { return getData< std::vector >(type, "element_to_subelement", ghost_type); } /* -------------------------------------------------------------------------- */ inline const Array & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) const { return getData(type, "subelement_to_element", ghost_type); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) { return getData(type, "subelement_to_element", ghost_type); } /* -------------------------------------------------------------------------- */ template inline Array * Mesh::getDataPointer(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type, UInt nb_component) { Array & tmp = mesh_data.getElementalDataArrayAlloc(data_name, el_type, ghost_type, nb_component); tmp.resize(getNbElement(el_type, ghost_type)); return &tmp; } /* -------------------------------------------------------------------------- */ template inline const Array & Mesh::getData(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type) const { return mesh_data.getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getData(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type) { return mesh_data.getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline const ByElementTypeArray & Mesh::getData(const std::string & data_name) const { return mesh_data.getElementalData(data_name); } /* -------------------------------------------------------------------------- */ template inline ByElementTypeArray & Mesh::getData(const std::string & data_name) { return mesh_data.getElementalData(data_name); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const ElementType & type, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); try { const Array & conn = connectivities(type, ghost_type); AKANTU_DEBUG_OUT(); return conn.getSize(); } catch (...) { AKANTU_DEBUG_OUT(); return 0; } } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const UInt spatial_dimension, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt nb_element = 0; type_iterator it = firstType(spatial_dimension, ghost_type, _ek_not_defined); type_iterator last = lastType(spatial_dimension, ghost_type, _ek_not_defined); for (; it != last; ++it) nb_element += getNbElement(*it, ghost_type); AKANTU_DEBUG_OUT(); return nb_element; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type) const { AKANTU_DEBUG_IN(); UInt * conn_val = getConnectivity(type, ghost_type).values; UInt nb_nodes_per_element = getNbNodesPerElement(type); Real local_coord[spatial_dimension * nb_nodes_per_element]; UInt offset = element * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { memcpy(local_coord + n * spatial_dimension, nodes->storage() + conn_val[offset + n] * spatial_dimension, spatial_dimension*sizeof(Real)); } Math::barycenter(local_coord, nb_nodes_per_element, spatial_dimension, barycenter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(const Element & element, Vector & barycenter) const { getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElement(const ElementType & type) { UInt nb_nodes_per_element = 0; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass::getNbNodesPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT); #undef GET_NB_NODES_PER_ELEMENT return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getP1ElementType(const ElementType & type) { ElementType p1_type = _not_defined; #define GET_P1_TYPE(type) \ p1_type = ElementClass::getP1ElementType() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_P1_TYPE); #undef GET_P1_TYPE return p1_type; } /* -------------------------------------------------------------------------- */ inline ElementKind Mesh::getKind(const ElementType & type) { ElementKind kind = _ek_not_defined; #define GET_KIND(type) \ kind = ElementClass::getKind() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_KIND); #undef GET_KIND return kind; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getSpatialDimension(const ElementType & type) { UInt spatial_dimension = 0; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass::getSpatialDimension() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION); #undef GET_SPATIAL_DIMENSION return spatial_dimension; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getFacetType(const ElementType & type) { ElementType surface_type = _not_defined; #define GET_FACET_TYPE(type) \ surface_type = ElementClass::getFacetType() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_TYPE); #undef GET_FACET_TYPE return surface_type; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline Matrix Mesh::getFacetLocalConnectivity(const ElementType & type) { AKANTU_DEBUG_IN(); Matrix mat; #define GET_FACET_CON(type) \ mat = ElementClass::getFacetLocalConnectivityPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON); #undef GET_FACET_CON AKANTU_DEBUG_OUT(); return mat; } /* -------------------------------------------------------------------------- */ inline Matrix Mesh::getFacetConnectivity(UInt element, const ElementType & type, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); Matrix local_facets = getFacetLocalConnectivity(type); Matrix facets(local_facets.rows(), local_facets.cols()); const Array & conn = connectivities(type, ghost_type); for (UInt f = 0; f < facets.rows(); ++f) { for (UInt n = 0; n < facets.cols(); ++n) { facets(f, n) = conn(element, local_facets(f, n)); } } AKANTU_DEBUG_OUT(); return facets; } /* -------------------------------------------------------------------------- */ template inline void Mesh::extractNodalValuesFromElement(const Array & nodal_values, T * local_coord, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const { for (UInt n = 0; n < n_nodes; ++n) { memcpy(local_coord + n * nb_degree_of_freedom, nodal_values.storage() + connectivity[n] * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(T)); } } /* -------------------------------------------------------------------------- */ #define DECLARE_GET_BOUND(Var, var) \ inline void Mesh::get##Var##Bounds(Real * var) const { \ for (UInt i = 0; i < spatial_dimension; ++i) { \ var[i] = var##_bounds[i]; \ } \ } \ DECLARE_GET_BOUND(Lower, lower) DECLARE_GET_BOUND(Upper, upper) DECLARE_GET_BOUND(LocalLower, local_lower) DECLARE_GET_BOUND(LocalUpper, local_upper) #undef DECLARE_GET_BOUND /* -------------------------------------------------------------------------- */ inline void Mesh::addConnectivityType(const ElementType & type, const GhostType & ghost_type){ getConnectivityPointer(type, ghost_type); } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPureGhostNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -3 : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -2 || (*nodes_type)(n) == -1 : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -1 : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -2 : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return nodes_type ? (*nodes_type)(n) >= 0 : false; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->getSize(); } /* -------------------------------------------------------------------------- */ inline Int Mesh::getNodeType(UInt local_id) const { return nodes_type ? (*nodes_type)(local_id) : -1; } /* -------------------------------------------------------------------------- */ inline const SubBoundary & Mesh::getSubBoundary(const std::string & name) const{ Boundary::const_iterator it = getBoundary().find(name); if(it == getBoundary().end()) { AKANTU_EXCEPTION("No sub-boundary named " << name << "!"); } return *it; } /* -------------------------------------------------------------------------- */ inline SubBoundary & Mesh::getSubBoundary(const std::string & name) { Boundary::iterator it = getBoundary().find(name); if(it == getBoundary().end()) { AKANTU_EXCEPTION("No sub-boundary named " << name << "!"); } return *it; } diff --git a/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh b/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh index 41d9ee9bb..3f238501a 100644 --- a/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh +++ b/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh @@ -1,254 +1,256 @@ /** * @file dumper_iohelper_tmpl_homogenizing_field.hh * * @author Nicolas Richart * * @date Fri Oct 26 21:52:40 2012 * * @brief description of field homogenizing field * * @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 . * */ /* -------------------------------------------------------------------------- */ -__END_AKANTU__ - -#include "solid_mechanics_model.hh" +namespace akantu { + class SolidMechanicsModel; + class FEM; +} __BEGIN_AKANTU__ - /* -------------------------------------------------------------------------- */ /* Homogenizing containers */ /* -------------------------------------------------------------------------- */ template class sub_type> class DumperIOHelper::PaddingHomogenizingFunctor { public: PaddingHomogenizingFunctor(Container & cont) : cont(cont){ } UInt getNbComponent() { typename Container::iterator it = cont.begin(); typename Container::iterator end = cont.end(); UInt nb_comp = 0; for (; it != end; ++it) nb_comp = std::max(nb_comp, (*it).size()); return nb_comp; } sub_type operator()(const sub_type & vect, __attribute__((unused)) const ElementType & type) { return vect; } protected: Container & cont; }; /* -------------------------------------------------------------------------- */ template class sub_type> class DumperIOHelper::AvgHomogenizingFunctor { public: AvgHomogenizingFunctor(Container & cont) : cont(cont){ } UInt getNbComponent() { typename Container::iterator it = cont.begin(); typename Container::iterator end = cont.end(); UInt nb_comp = 0; UInt i = 0; for (; it != end; ++it) { ElementType type = it.getType(); UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type); nb_comp = std::max(nb_comp, (*it).size() / nb_quad); ++i; } return nb_comp; } sub_type operator()(const sub_type & vect, const ElementType & type) { UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type); if(nb_quad == 1) { return vect; } else { UInt s = vect.size() / nb_quad; sub_type v(s); for (UInt i = 0; i < s; ++i) { for (UInt q = 0; q < nb_quad; ++q) { v[i] += vect[i + q*s] / Real(nb_quad); } } return v; } } protected: Container & cont; }; /* specialization */ template class DumperIOHelper::AvgHomogenizingFunctor { public: AvgHomogenizingFunctor(Container & cont) : cont(cont){ } UInt getNbComponent() { typename Container::iterator it = cont.begin(); typename Container::iterator end = cont.end(); UInt nb_comp = 0; UInt i = 0; for (; it != end; ++it) { ElementType type = it.getType(); UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type); nb_comp = std::max(nb_comp, (*it).size() / nb_quad); ++i; } return nb_comp; } Matrix operator()(const Matrix & vect, const ElementType & type) { UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type); if(nb_quad == 1) { return vect; } else { UInt m = vect.rows(); UInt n = vect.cols() / nb_quad; Matrix v(m, n); for (UInt q = 0; q < nb_quad; ++q) for (UInt i = 0; i < m; ++i) for (UInt j = 0; j < n; ++j) v(i, j) += vect(i, j + q*n) / Real(nb_quad); return v; } } protected: Container & cont; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class, template class, bool> class, bool> class Container, template class, bool> class int_sub_iterator, template class> class Funct, template class ret_type, bool filtered> class DumperIOHelper::HomogenizedField : public Field { protected: typedef typename Container::iterator sub_iterator; typedef Funct, ret_type> Functor; public: class iterator { public: iterator(const sub_iterator & it, Functor & func) : it(it), func(func) {} bool operator!=(const iterator & it) { return it.it != this->it; } iterator operator++() { ++this->it; return *this; } ret_type operator*() { return func(*it, it.getType()); } protected: sub_iterator it; Functor & func; }; HomogenizedField(const SolidMechanicsModel & model, const std::string & field_id, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : cont(model, field_id, spatial_dimension, ghost_type, element_kind, filter), funct(cont) { nb_component = funct.getNbComponent(); } HomogenizedField(const SolidMechanicsModel & model, const std::string & field_id, UInt n, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : cont(model, field_id, n, spatial_dimension, ghost_type, element_kind, filter), funct(cont) { nb_component = funct.getNbComponent(); } HomogenizedField(const FEM & fem, const ByElementTypeArray & field, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : cont(fem, field, spatial_dimension, ghost_type, element_kind, filter), funct(cont) { nb_component = funct.getNbComponent(); } HomogenizedField(const FEM & fem, const ByElementTypeArray & field, UInt n, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : cont(fem, field, n, spatial_dimension, ghost_type, element_kind, filter), funct(cont) { nb_component = funct.getNbComponent(); } iterator begin() { return iterator(cont.begin(), funct); } iterator end () { return iterator(cont.end(), funct); } virtual void registerToDumper(const std::string & id, iohelper::Dumper & dupmer) { dupmer.addElemDataField(id, *this); } UInt getDim() { if(padding_n && padding_m) return padding_m*padding_n; else return nb_component; } UInt size() { return cont.size(); } iohelper::DataType getDataType() { return iohelper::getDataType(); } UInt isHomogeneous() { return true; } void setPadding(UInt n, UInt m = 1) { Field::setPadding(n, m); cont.setPadding(n, m); } protected: Container cont; Functor funct; UInt nb_component; }; + +__END_AKANTU__ diff --git a/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh b/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh index 0f2be8511..c0d043c6f 100644 --- a/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh +++ b/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh @@ -1,384 +1,390 @@ /** * @file dumper_iohelper_tmpl_material_internal_field.hh * * @author Nicolas Richart * * @date Fri Oct 26 21:52:40 2012 * * @brief description of material internal field * * @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__ + /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class R> class DumperIOHelper::MaterialPaddingHelper : public PaddingHelper { public: MaterialPaddingHelper(const SolidMechanicsModel & model) : model(model) { } inline R pad(const R & in, UInt padding_m, UInt padding_n, UInt nb_data, __attribute__((unused)) UInt material_id) { return PaddingHelper::pad(in, padding_m, padding_n, nb_data); } protected: const SolidMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ template class R> class DumperIOHelper::StressPaddingHelper : public MaterialPaddingHelper { public: StressPaddingHelper(const SolidMechanicsModel & model) : MaterialPaddingHelper(model) { AKANTU_DEBUG_ERROR("This class exists only to pad stress (Matrix) to 3D"); } inline R pad(const R & in, UInt padding_m, UInt padding_n, UInt nb_data, UInt material_id) {} protected: const Material & material; }; template<> class DumperIOHelper::StressPaddingHelper : public MaterialPaddingHelper { public: StressPaddingHelper(const SolidMechanicsModel & model) : MaterialPaddingHelper(model) { } inline Matrix pad(const Matrix & in, UInt padding_m, UInt padding_n, UInt nb_data, UInt material_id) { if(padding_m <= in.rows() && padding_n * nb_data <= in.cols()) return in; else { AKANTU_DEBUG_ASSERT(padding_m == 3 && padding_n == 3, "This function can only pad to 3D"); Matrix stress = MaterialPaddingHelper::pad(in, 3, 3, nb_data, material_id); if(in.rows() == 2 && in.cols() == 2 * nb_data) { const Material & material = model.getMaterial(material_id); bool plane_strain = !material.getParam("Plane_Stress"); if(plane_strain) { Real nu = material.getParam("nu"); for (UInt d = 0; d < nb_data; ++d) { stress(2, 2 + 3*d) = nu * (stress(0, 0 + 3*d) + stress(1, 1 + 3*d)); } } } return stress; } } }; /* -------------------------------------------------------------------------- */ template class R> class DumperIOHelper::StrainPaddingHelper : public MaterialPaddingHelper { public: StrainPaddingHelper(const SolidMechanicsModel & model) : MaterialPaddingHelper(model) { AKANTU_DEBUG_ERROR("This class exists only to pad strain (Matrix) to 3D"); } inline R pad(const R & in, UInt padding_m, UInt padding_n, UInt nb_data, UInt material_id) {} }; template<> class DumperIOHelper::StrainPaddingHelper : public MaterialPaddingHelper { public: StrainPaddingHelper(const SolidMechanicsModel & model) : MaterialPaddingHelper(model) { } inline Matrix pad(const Matrix & in, UInt padding_m, UInt padding_n, UInt nb_data, UInt material_id) { if(padding_m <= in.rows() && padding_n * nb_data <= in.cols()) return in; else { AKANTU_DEBUG_ASSERT(padding_m == 3 && padding_n == 3, "This function can only pad to 3D"); Matrix strain = MaterialPaddingHelper::pad(in, 3, 3, nb_data, material_id); if(in.rows() == 2 && in.cols() == 2 * nb_data) { const Material & material = model.getMaterial(material_id); bool plane_stress = material.getParam("Plane_Stress"); if(plane_stress) { Real nu = material.getParam("nu"); for (UInt d = 0; d < nb_data; ++d) { strain(2, 2 + 3*d) = nu / (nu - 1) * (strain(0, 0 + 3*d) + strain(1, 1 + 3*d)); } } } return strain; } } }; /* -------------------------------------------------------------------------- */ /* Element material field iterator/container */ /* -------------------------------------------------------------------------- */ template class ret_type, template class> class padding_helper_type, template class, bool> class int_iterator, bool filtered> class DumperIOHelper::generic_internal_material_field_iterator : public generic_quadrature_point_iterator, filtered> { public: typedef generic_quadrature_point_iterator, filtered> parent; typedef typename parent::it_type it_type; typedef typename parent::data_type data_type; typedef typename parent::return_type return_type; typedef typename parent::field_type field_type; typedef typename parent::internal_iterator internal_iterator; typedef typename Array::template const_iterator< ret_type > internal_material_iterator; public: generic_internal_material_field_iterator(const field_type & element_material, UInt n, const typename field_type::type_iterator & t_it, const typename field_type::type_iterator & t_it_end, const internal_iterator & it, ElementType element_type, const GhostType ghost_type = _not_ghost, const ByElementTypeArray * filter = NULL, UInt * fit = NULL) : parent(element_material, 2, t_it, t_it_end, it, element_type, ghost_type, filter, fit), out_n(n), model(NULL), padding_helper(NULL) { } ~generic_internal_material_field_iterator() { delete padding_helper; } generic_internal_material_field_iterator(const generic_internal_material_field_iterator & it) : parent(it) { out_n = it.out_n; field_id = it.field_id; if(it.model) { model = it.model; padding_helper = new padding_helper_type(*model); } } return_type operator*() { const ret_type & material_id = *this->vit; UInt nb_data = this->fem->getNbQuadraturePoints(*this->tit); try { const Array & vect = model->getMaterial(material_id[1]).getArray(field_id, *this->tit, this->ghost_type); if(vect.getSize() == 0 || vect.getSize() < material_id[0]) // vector exists but has a wrong size return return_type(); UInt ln = out_n; if(out_n == 0) ln = vect.getNbComponent(); internal_material_iterator it = iterator_helper::begin(vect, ln, vect.getNbComponent() / ln * nb_data, vect.getSize() / nb_data); it += material_id[0]; return padding_helper->pad(*it, this->padding_m, this->padding_n, nb_data, material_id[1]); } catch (...) { return return_type(); } } void setModel(const SolidMechanicsModel & model) { this->model = &model; this->padding_helper = new padding_helper_type(model); } void setFieldID(const std::string & field_id) { this->field_id = field_id; } protected: virtual UInt getNbDataPerElem(__attribute__((unused)) const ElementType & type) { return 1; } protected: UInt out_n; const SolidMechanicsModel * model; std::string field_id; padding_helper_type * padding_helper; }; /* -------------------------------------------------------------------------- */ template class ret_type, bool filtered> class DumperIOHelper::internal_material_field_iterator : public generic_internal_material_field_iterator { public: typedef generic_internal_material_field_iterator parent; typedef typename parent::it_type it_type; typedef typename parent::data_type data_type; typedef typename parent::return_type return_type; typedef typename parent::field_type field_type; typedef typename parent::internal_iterator internal_iterator; typedef typename Array::template const_iterator< ret_type > internal_material_iterator; public: internal_material_field_iterator(const field_type & element_material, UInt n, const typename field_type::type_iterator & t_it, const typename field_type::type_iterator & t_it_end, const internal_iterator & it, ElementType element_type, const GhostType ghost_type = _not_ghost, const ByElementTypeArray * filter = NULL, UInt * fit = NULL) : parent(element_material, n, t_it, t_it_end, it, element_type, ghost_type, filter, fit) { } }; /* -------------------------------------------------------------------------- */ template class ret_type, bool filtered> class DumperIOHelper::material_stress_field_iterator : public generic_internal_material_field_iterator { public: typedef generic_internal_material_field_iterator parent; typedef typename parent::it_type it_type; typedef typename parent::data_type data_type; typedef typename parent::return_type return_type; typedef typename parent::field_type field_type; typedef typename parent::internal_iterator internal_iterator; typedef typename Array::template const_iterator< ret_type > internal_material_iterator; public: material_stress_field_iterator(const field_type & element_material, UInt n, const typename field_type::type_iterator & t_it, const typename field_type::type_iterator & t_it_end, const internal_iterator & it, ElementType element_type, const GhostType ghost_type = _not_ghost, const ByElementTypeArray * filter = NULL, UInt * fit = NULL) : parent(element_material, n, t_it, t_it_end, it, element_type, ghost_type, filter, fit) { } }; /* -------------------------------------------------------------------------- */ template class ret_type, bool filtered> class DumperIOHelper::material_strain_field_iterator : public generic_internal_material_field_iterator { public: typedef generic_internal_material_field_iterator parent; typedef typename parent::it_type it_type; typedef typename parent::data_type data_type; typedef typename parent::return_type return_type; typedef typename parent::field_type field_type; typedef typename parent::internal_iterator internal_iterator; typedef typename Array::template const_iterator< ret_type > internal_material_iterator; public: material_strain_field_iterator(const field_type & element_material, UInt n, const typename field_type::type_iterator & t_it, const typename field_type::type_iterator & t_it_end, const internal_iterator & it, ElementType element_type, const GhostType ghost_type = _not_ghost, const ByElementTypeArray * filter = NULL, UInt * fit = NULL) : parent(element_material, n, t_it, t_it_end, it, element_type, ghost_type, filter, fit) { } }; /* -------------------------------------------------------------------------- */ template class ret_type, template class, bool> class iterator_type, bool filtered> class DumperIOHelper::InternalMaterialField : public GenericQuadraturePointsField, ret_type, filtered> { typedef GenericQuadraturePointsField, ret_type, filtered> parent; public: typedef typename parent::iterator iterator; public: /* ------------------------------------------------------------------------ */ InternalMaterialField(const SolidMechanicsModel & model, const std::string & field_id, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : parent(model.getFEM(), model.getElementIndexByMaterial(), 0, spatial_dimension, ghost_type, element_kind, filter), model(model), field_id(field_id) { init(); } InternalMaterialField(const SolidMechanicsModel & model, const std::string & field_id, UInt n, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : parent(model.getFEM(), model.getElementIndexByMaterial(), n, spatial_dimension, ghost_type, element_kind, filter), model(model), field_id(field_id) { init(); } virtual iterator begin() { iterator it = parent::begin(); it.setItSize(2,1); it.setModel(model); it.setFieldID(field_id); return it; } virtual iterator end () { iterator it = parent::end(); it.setItSize(2,1); it.setModel(model); it.setFieldID(field_id); return it; } protected: void init() { typename ByElementTypeArray::type_iterator tit = this->field.firstType(this->spatial_dimension, this->ghost_type, this->element_kind); typename ByElementTypeArray::type_iterator end = this->field.lastType (this->spatial_dimension, this->ghost_type, this->element_kind); UInt nb_materials = model.getNbMaterials(); bool found = false; for(;tit != end; ++tit) { // UInt nb_quad = this->fem.getNbQuadraturePoints(*tit); for (UInt ma = 0; ma < nb_materials; ++ma) { const Material & mat = model.getMaterial(ma); try { const Array & vect __attribute__((unused)) = mat.getArray(field_id, *tit, this->ghost_type); found = true; } catch (...) { }; } } if(!found) AKANTU_DEBUG_ERROR("The field " << field_id << " does not exists in any materials"); } virtual UInt getNbDataPerElem(__attribute__((unused)) const ElementType & type) { return 1; } protected: const SolidMechanicsModel & model; std::string field_id; }; + + +__END_AKANTU__ + diff --git a/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh b/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh index 2fb4da5c7..8de6516ed 100644 --- a/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh +++ b/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh @@ -1,180 +1,180 @@ /** * @file dumper_iohelper_tmpl_quadrature_points_field.hh * * @author Nicolas Richart * * @date Fri Oct 26 21:52:40 2012 * * @brief Description of quadrature points fields * * @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 . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "fem.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class ret_type, class daughter, bool filtered> class DumperIOHelper::generic_quadrature_point_iterator : public element_iterator { public: typedef element_iterator parent; typedef typename parent::it_type it_type; typedef typename parent::data_type data_type; typedef typename parent::return_type return_type; typedef typename parent::field_type field_type; typedef typename parent::internal_iterator internal_iterator; public: generic_quadrature_point_iterator(const field_type & field, UInt n, const typename field_type::type_iterator & t_it, const typename field_type::type_iterator & t_it_end, const internal_iterator & it, ElementType element_type, const GhostType ghost_type = _not_ghost, const ByElementTypeArray * filter = NULL, UInt * fit = NULL) : parent(field, n, t_it, t_it_end, it, element_type, ghost_type, filter, fit) { } void setFEM(const FEM & fem) { this->fem = &fem; } protected: virtual UInt getNbDataPerElem(const ElementType & type) { return fem->getNbQuadraturePoints(type); } protected: const FEM * fem; }; /* -------------------------------------------------------------------------- */ template class ret_type, bool filtered> class DumperIOHelper::quadrature_point_iterator : public generic_quadrature_point_iterator, filtered > { public: typedef generic_quadrature_point_iterator, filtered> parent; typedef typename parent::it_type it_type; typedef typename parent::data_type data_type; typedef typename parent::return_type return_type; typedef typename parent::field_type field_type; typedef typename parent::internal_iterator internal_iterator; public: quadrature_point_iterator(const field_type & field, UInt n, const typename field_type::type_iterator & t_it, const typename field_type::type_iterator & t_it_end, const internal_iterator & it, ElementType element_type, const GhostType ghost_type = _not_ghost, const ByElementTypeArray * filter = NULL, UInt * fit = NULL) : parent(field, n, t_it, t_it_end, it, element_type, ghost_type, filter, fit) { } return_type operator*() { UInt nb_data = this->fem->getNbQuadraturePoints(*this->tit); return PaddingHelper::pad(*this->vit, this->padding_m, this->padding_n, nb_data); } }; /* -------------------------------------------------------------------------- */ /* Fields type description */ /* -------------------------------------------------------------------------- */ template class ret_type, bool filtered> class DumperIOHelper::GenericQuadraturePointsField : public GenericElementalField { public: typedef iterator_type iterator; typedef GenericElementalField parent; GenericQuadraturePointsField(const FEM & fem, const ByElementTypeArray & field, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : parent(field, spatial_dimension, - ghost_type, element_kind, filtered), fem(fem) { } + ghost_type, element_kind, filter), fem(fem) { } GenericQuadraturePointsField(const FEM & fem, const ByElementTypeArray & field, UInt n, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : parent(field, n, spatial_dimension, ghost_type, element_kind, filter), fem(fem) { } /* ------------------------------------------------------------------------ */ virtual iterator begin() { iterator it = parent::begin(); it.setFEM(fem); return it; } virtual iterator end () { iterator it = parent::end(); it.setFEM(fem); return it; } const FEM & getFEM() { return fem; } protected: virtual UInt getNbDataPerElem(const ElementType & type) { return fem.getNbQuadraturePoints(type); } protected: const FEM & fem; }; /* -------------------------------------------------------------------------- */ template class ret_type, template class, bool> class iterator_type, bool filtered> class DumperIOHelper::QuadraturePointsField : public GenericQuadraturePointsField, ret_type, filtered> { public: typedef iterator_type iterator; typedef GenericQuadraturePointsField parent; QuadraturePointsField(const FEM & fem, const ByElementTypeArray & field, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : parent(fem, field, spatial_dimension, ghost_type, element_kind, filter) { } QuadraturePointsField(const FEM & fem, const ByElementTypeArray & field, UInt n, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined, const ByElementTypeArray * filter = NULL) : parent(fem, field, n, spatial_dimension, ghost_type, element_kind, filter) { } }; diff --git a/src/io/dumper/dumper_paraview.hh b/src/io/dumper/dumper_paraview.hh index d99d66578..dd5797f2f 100644 --- a/src/io/dumper/dumper_paraview.hh +++ b/src/io/dumper/dumper_paraview.hh @@ -1,75 +1,75 @@ /** * @file dumper_paraview.hh * * @author Nicolas Richart * * @date Fri Oct 26 21:52:40 2012 * * @brief Dumper Paraview using IOHelper * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DUMPER_PARAVIEW_HH__ #define __AKANTU_DUMPER_PARAVIEW_HH__ #include "dumper_iohelper.hh" __BEGIN_AKANTU__ class DumperParaview : public DumperIOHelper { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - DumperParaview(const std::string & filename, - const std::string & directory = "./paraview", - bool parallel = true); + DumperParaview(const std::string & filename, + const std::string & directory = "./paraview", + bool parallel = true); virtual ~DumperParaview(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: // void dump(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: void setBaseName(const std::string & basename) { DumperIOHelper::setBaseName(basename); static_cast(dumper)->setVTUSubDirectory(filename + "-VTU"); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; __END_AKANTU__ #endif /* __AKANTU_DUMPER_PARAVIEW_HH__ */ diff --git a/src/io/mesh_io/mesh_io_diana.cc b/src/io/mesh_io/mesh_io_diana.cc index f07573338..8e1092c8b 100644 --- a/src/io/mesh_io/mesh_io_diana.cc +++ b/src/io/mesh_io/mesh_io_diana.cc @@ -1,500 +1,500 @@ /** * @file mesh_io_diana.cc * * @author David Simon Kammer * @author Alodie Schneuwly * * @date Sat Mar 26 20:43:38 2011 * * @brief handles diana meshes * * @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 "mesh_io_diana.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIODiana::MeshIODiana() { canReadSurface = true; canReadExtendedData = true; _diana_to_akantu_element_types["T9TM"] = _triangle_3; _diana_to_akantu_element_types["Q12TM"] = _quadrangle_4; _diana_to_akantu_element_types["TP18L"] = _pentahedron_6; _diana_to_akantu_element_types["TE12L"] = _tetrahedron_4; _diana_to_akantu_element_types["HX24L"] = _hexahedron_8; _diana_to_akantu_mat_prop["YOUNG"] = "E"; _diana_to_akantu_mat_prop["DENSIT"] = "rho"; _diana_to_akantu_mat_prop["POISON"] = "nu"; } /* -------------------------------------------------------------------------- */ MeshIODiana::~MeshIODiana() { std::map *>::iterator ng_it; std::map *>::iterator eg_it; for (ng_it = node_groups.begin(); ng_it != node_groups.end(); ++ng_it) { delete ng_it->second; } for (eg_it = element_groups.begin(); eg_it != element_groups.end(); ++eg_it) { delete eg_it->second; } } /* -------------------------------------------------------------------------- */ inline void my_getline(std::ifstream & infile, std::string & line) { std::getline(infile, line); //read the line size_t pos = line.find("\r"); /// remove the extra \r if needed line = line.substr(0, pos); } /* -------------------------------------------------------------------------- */ void MeshIODiana::read(const std::string & filename, Mesh & mesh) { AKANTU_DEBUG_IN(); std::ifstream infile; infile.open(filename.c_str()); std::string line; UInt first_node_number = std::numeric_limits::max(); std::vector global_to_local_index; if(!infile.good()) { AKANTU_DEBUG_ERROR("Cannot open file " << filename); } while(infile.good()) { my_getline(infile, line); /// read all nodes if(line == "'COORDINATES'") { line = readCoordinates(infile, mesh, first_node_number); } /// read all elements if (line == "'ELEMENTS'") { line = readElements(infile, mesh, global_to_local_index, first_node_number); } /// read the material properties and write a .dat file if (line == "'MATERIALS'") { line = readMaterial(infile, filename); } /// read the material properties and write a .dat file if (line == "'GROUPS'") { line = readGroups(infile, global_to_local_index, first_node_number); } } infile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshIODiana::write(__attribute__((unused)) const std::string & filename, __attribute__((unused)) const Mesh & mesh) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readCoordinates(std::ifstream & infile, Mesh & mesh, UInt & first_node_number) { AKANTU_DEBUG_IN(); Array & nodes = const_cast &>(mesh.getNodes()); std::string line; UInt index; Real coord[3]; do { my_getline(infile, line); if("'ELEMENTS'" == line) break; //end = true; //else { /// for each node, read the coordinates std::stringstream sstr_node(line); sstr_node >> index >> coord[0] >> coord[1] >> coord[2]; //if (!sstr_node.fail()) //break; first_node_number = first_node_number < index ? first_node_number : index; nodes.push_back(coord); // } } while(true);//!end); AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ UInt MeshIODiana::readInterval(std::stringstream & line, std::set & interval) { UInt first; line >> first; if(line.fail()) { return 0; } interval.insert(first); UInt second; int dash; dash = line.get(); if(dash == '-') { line >> second; interval.insert(second); return 2; } if(line.fail()) line.clear(std::ios::eofbit); // in case of get at end of the line else line.unget(); return 1; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readGroups(std::ifstream & infile, std::vector & global_to_local_index, UInt first_node_number) { AKANTU_DEBUG_IN(); std::string line; my_getline(infile, line); bool reading_nodes_group = false; while(line != "'SUPPORTS'") { if(line == "NODES") { reading_nodes_group = true; my_getline(infile, line); } if(line == "ELEMEN") { reading_nodes_group = false; my_getline(infile, line); } std::stringstream *str = new std::stringstream(line); UInt id; std::string name; char c; *str >> id >> name >> c; Array * list_ids = new Array(0, 1, name); UInt s = 1; bool end = false; while(!end) { while(!str->eof() && s != 0) { std::set interval; s = readInterval(*str, interval); std::set::iterator it = interval.begin(); if(s == 1) list_ids->push_back(*it); if(s == 2) { UInt first = *it; ++it; UInt second = *it; for(UInt i = first; i <= second; ++i) { list_ids->push_back(i); } } } if(str->fail()) end = true; else { my_getline(infile, line); delete str; str = new std::stringstream(line); } } delete str; if(reading_nodes_group) { for (UInt i = 0; i < list_ids->getSize(); ++i) { (*list_ids)(i) -= first_node_number; } node_groups[name] = list_ids; } else { std::vector * elem = new std::vector; elem->reserve(list_ids->getSize()); for (UInt i = 0; i < list_ids->getSize(); ++i) { Element & e = global_to_local_index[(*list_ids)(i)-1]; if(e.type != _not_defined) elem->push_back(e); } element_groups[name] = elem; delete list_ids; } my_getline(infile, line); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readElements(std::ifstream & infile, Mesh & mesh, std::vector & global_to_local_index, UInt first_node_number) { AKANTU_DEBUG_IN(); std::string line; my_getline(infile, line); if("CONNECTIVITY" == line) { line = readConnectivity(infile, mesh, global_to_local_index, first_node_number); } /// read the line corresponding to the materials if ("MATERIALS" == line) { line = readMaterialElement(infile, mesh, global_to_local_index); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readConnectivity(std::ifstream & infile, Mesh & mesh, std::vector & global_to_local_index, UInt first_node_number) { AKANTU_DEBUG_IN(); Int index; std::string lline; std::string diana_type; ElementType akantu_type, akantu_type_old = _not_defined; Array *connectivity = NULL; UInt node_per_element = 0; Element elem; bool end = false; do { my_getline(infile, lline); std::stringstream sstr_elem(lline); if(lline == "MATERIALS") end = true; else { /// traiter les coordonnees sstr_elem >> index; sstr_elem >> diana_type; akantu_type = _diana_to_akantu_element_types[diana_type]; if(akantu_type != _not_defined) { if(akantu_type != akantu_type_old) { connectivity = mesh.getConnectivityPointer(akantu_type); node_per_element = connectivity->getNbComponent(); akantu_type_old = akantu_type; } UInt local_connect[node_per_element]; for(UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; node_index -= first_node_number; local_connect[j] = node_index; } connectivity->push_back(local_connect); elem.type = akantu_type; elem.element = connectivity->getSize() - 1; } else { elem.type = _not_defined; elem.element = UInt(-1); } global_to_local_index.push_back(elem); } } while(!end); AKANTU_DEBUG_OUT(); return lline; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readMaterialElement(std::ifstream & infile, Mesh & mesh, std::vector & global_to_local_index) { AKANTU_DEBUG_IN(); std::string line; std::stringstream sstr_tag_name; sstr_tag_name << "tag_" << 0; Mesh::type_iterator it = mesh.firstType(); Mesh::type_iterator end = mesh.lastType(); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); - mesh.getDataPointer(*it, "material", _not_ghost)->resize(nb_element); + mesh.getDataPointer(*it, "material", _not_ghost, 1)->resize(nb_element); } my_getline(infile, line); while(line != "'MATERIALS'") { line = line.substr(line.find('/') + 1, std::string::npos); // erase the first slash / of the line char tutu[250]; strcpy(tutu, line.c_str()); Array temp_id(0, 2); UInt mat; while(true){ std::stringstream sstr_intervals_elements(line); UInt id[2]; char temp; while(sstr_intervals_elements.good()) { sstr_intervals_elements >> id[0] >> temp >> id[1]; // >> "/" >> mat; if(!sstr_intervals_elements.fail()) temp_id.push_back(id); } if (sstr_intervals_elements.fail()) { sstr_intervals_elements.clear(); sstr_intervals_elements.ignore(); sstr_intervals_elements >> mat; break; } my_getline(infile, line); } // loop over elements // UInt * temp_id_val = temp_id.values; for (UInt i = 0; i < temp_id.getSize(); ++i) - for (UInt j=temp_id(i,0); j<=temp_id(i,1); ++j) { + for (UInt j = temp_id(i, 0); j <= temp_id(i, 1); ++j) { Element & element = global_to_local_index[j - 1]; if(element.type == _not_defined) continue; UInt elem = element.element; ElementType type = element.type; Array & data = *(mesh.getDataPointer(type, "material", _not_ghost)); data(elem) = mat; } my_getline(infile, line); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readMaterial(std::ifstream & infile, const std::string & filename) { AKANTU_DEBUG_IN(); std::stringstream mat_file_name; mat_file_name << "material_" << filename; std::ofstream material_file; material_file.open(mat_file_name.str().c_str());//mat_file_name.str()); UInt mat_index; std::string line; bool first_mat = true; bool end = false; typedef std::map MatProp; MatProp mat_prop; do{ my_getline(infile, line); std::stringstream sstr_material(line); if("'GROUPS'" == line) { if(!mat_prop.empty()) { material_file << "material elastic [" << std::endl; for(MatProp::iterator it = mat_prop.begin(); it != mat_prop.end(); ++it) material_file << "\t" << it->first << " = " << it->second << std::endl; material_file << "]" << std::endl; mat_prop.clear(); } end = true; } else { /// traiter les caractéristiques des matériaux sstr_material >> mat_index; if(!sstr_material.fail()) { if(!first_mat) { if(!mat_prop.empty()) { material_file << "material elastic [" << std::endl; for(MatProp::iterator it = mat_prop.begin(); it != mat_prop.end(); ++it) material_file << "\t" << it->first << " = " << it->second << std::endl; material_file << "]" << std::endl; mat_prop.clear(); } } first_mat = false; } else { sstr_material.clear(); } std::string prop_name; sstr_material >> prop_name; std::map::iterator it; it = _diana_to_akantu_mat_prop.find(prop_name); if(it != _diana_to_akantu_mat_prop.end()) { Real value; sstr_material >> value; mat_prop[it->second] = value; } else { AKANTU_DEBUG_INFO("In material reader, property " << prop_name << "not recognized"); } } } while (!end); AKANTU_DEBUG_OUT(); return line; } __END_AKANTU__ diff --git a/src/mesh_utils/mesh_partition.cc b/src/mesh_utils/mesh_partition.cc index 741c9f2b7..199d06e61 100644 --- a/src/mesh_utils/mesh_partition.cc +++ b/src/mesh_utils/mesh_partition.cc @@ -1,377 +1,380 @@ /** * @file mesh_partition.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date Tue Aug 17 16:19:43 2010 * * @brief implementation of common part of all partitioner * * @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 "mesh_partition.hh" #include "mesh_utils.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MeshPartition::MeshPartition(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id) : Memory(memory_id), id("MeshPartitioner"), mesh(mesh), spatial_dimension(spatial_dimension), partitions ("partition" , id, memory_id), ghost_partitions ("ghost_partition" , id, memory_id), ghost_partitions_offset("ghost_partition_offset", id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MeshPartition::~MeshPartition() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * conversion in c++ of the GENDUALMETIS (mesh.c) function wrote by George in * Metis (University of Minnesota) */ void MeshPartition::buildDualGraph(Array & dxadj, Array & dadjncy, Array & edge_loads, const EdgeLoadFunctor & edge_load_func, const Array & pairs) { AKANTU_DEBUG_IN(); // tweak mesh; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_element_p1[nb_types]; UInt magic_number[nb_types]; // UInt * conn_val[nb_types]; UInt nb_element[nb_types]; Array * conn[nb_types]; Array * conn_tmp[nb_types]; Array lin_to_element; Element elem; elem.ghost_type = _not_ghost; + const_cast(mesh).updateTypesOffsets(_not_ghost); + const_cast(mesh).updateTypesOffsets(_ghost); + Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != mesh.getSpatialDimension()) continue; elem.type = type; ElementType type_p1 = Mesh::getP1ElementType(type); nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); nb_nodes_per_element_p1[nb_good_types] = Mesh::getNbNodesPerElement(type_p1); nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).getSize(); magic_number[nb_good_types] = Mesh::getNbNodesPerElement(Mesh::getFacetType(type_p1)); conn[nb_good_types] = &const_cast &>(mesh.getConnectivity(type, _not_ghost)); for (UInt i = 0; i < nb_element[nb_good_types]; ++i) { elem.element = i; lin_to_element.push_back(elem); } if(pairs.getSize() != 0) { conn_tmp[nb_good_types] = new Array(mesh.getConnectivity(type, _not_ghost)); for (UInt i = 0; i < pairs.getSize(); ++i) { for (UInt el = 0; el < nb_element[nb_good_types]; ++el) { for (UInt n = 0; n < nb_nodes_per_element[nb_good_types]; ++n) { if(pairs(i, 1) == (*conn[nb_good_types])(el, n)) (*conn[nb_good_types])(el, n) = pairs(i, 0); } } } } nb_good_types++; } CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem); UInt nb_total_element = 0; UInt nb_total_node_element = 0; for (UInt t = 0; t < nb_good_types; ++t) { nb_total_element += nb_element[t]; nb_total_node_element += nb_element[t]*nb_nodes_per_element_p1[t]; } dxadj.resize(nb_total_element + 1); /// initialize the dxadj array for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) dxadj(linerized_el) = nb_nodes_per_element_p1[t]; /// convert the dxadj_val array in a csr one for (UInt i = 1; i < nb_total_element; ++i) dxadj(i) += dxadj(i-1); for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i-1); dxadj(0) = 0; dadjncy.resize(2*dxadj(nb_total_element)); edge_loads.resize(2*dxadj(nb_total_element)); /// weight map to determine adjacency unordered_map::type weight_map; for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) { /// fill the weight map for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n) { UInt node = (*conn[t])(el, n); CSR::iterator k; for (k = node_to_elem.rbegin(node); k != node_to_elem.rend(node); --k) { UInt current_el = *k; if(current_el <= linerized_el) break; unordered_map::type::iterator it_w; it_w = weight_map.find(current_el); if(it_w == weight_map.end()) { weight_map[current_el] = 1; } else { it_w->second++; } } } /// each element with a weight of the size of a facet are adjacent unordered_map::type::iterator it_w; for(it_w = weight_map.begin(); it_w != weight_map.end(); ++it_w) { if(it_w->second == magic_number[t]) { UInt adjacent_el = it_w->first; #if defined(AKANTU_COHESIVE_ELEMENT) /// Patch in order to prevent neighboring cohesive elements /// from detecting each other ElementKind linearized_el_kind = mesh.linearizedToElement(linerized_el).kind; ElementKind adjacent_el_kind = mesh.linearizedToElement(adjacent_el).kind; if (linearized_el_kind == adjacent_el_kind && linearized_el_kind == _ek_cohesive) continue; #endif UInt index_adj = dxadj(adjacent_el )++; UInt index_lin = dxadj(linerized_el)++; dadjncy(index_lin) = adjacent_el; dadjncy(index_adj) = linerized_el; } } weight_map.clear(); } } if(pairs.getSize() != 0) { for (UInt i = 0; i < nb_good_types; ++i) { conn[i]->copy(*conn_tmp[i]); delete conn_tmp[i]; } } Int k_start = 0; for (UInt t = 0, linerized_el = 0, j = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) { for (Int k = k_start; k < dxadj(linerized_el); ++k, ++j) dadjncy(j) = dadjncy(k); dxadj(linerized_el) = j; k_start += nb_nodes_per_element_p1[t]; } for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i - 1); dxadj(0) = 0; UInt adj = 0; for (UInt i = 0; i < nb_total_element; ++i) { UInt nb_adj = dxadj(i + 1) - dxadj(i); for (UInt j = 0; j < nb_adj; ++j, ++adj) { Int el_adj_id = dadjncy(dxadj(i) + j); Element el = lin_to_element(i); Element el_adj = lin_to_element(el_adj_id); Int load = edge_load_func(el, el_adj); edge_loads(adj) = load; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::fillPartitionInformation(const Mesh & mesh, const Int * linearized_partitions) { AKANTU_DEBUG_IN(); CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined); UInt linearized_el = 0; for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); partitions .alloc(nb_element, 1, type, _not_ghost); CSR & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); ghost_partitions .alloc(0, 1, type, _ghost); const Array & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el, ++linearized_el) { UInt part = linearized_partitions[linearized_el]; partitions(type, _not_ghost)(el) = part; std::list list_adj_part; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity.values[el * nb_nodes_per_element + n]; CSR::iterator ne; for (ne = node_to_elem.begin(node); ne != node_to_elem.end(node); ++ne) { UInt adj_el = *ne; UInt adj_part = linearized_partitions[adj_el]; if(part != adj_part) { list_adj_part.push_back(adj_part); } } } list_adj_part.sort(); list_adj_part.unique(); for(std::list::iterator adj_it = list_adj_part.begin(); adj_it != list_adj_part.end(); ++adj_it) { ghost_part_csr.getRows().push_back(*adj_it); ghost_part_csr.rowOffset(el)++; ghost_partitions(type, _ghost).push_back(*adj_it); ghost_partitions_offset(type, _ghost)(el)++; } } ghost_part_csr.countToCSR(); /// convert the ghost_partitions_offset array in an offset array Array & ghost_partitions_offset_ptr = ghost_partitions_offset(type, _ghost); for (UInt i = 1; i < nb_element; ++i) ghost_partitions_offset_ptr(i) += ghost_partitions_offset_ptr(i-1); for (UInt i = nb_element; i > 0; --i) ghost_partitions_offset_ptr(i) = ghost_partitions_offset_ptr(i-1); ghost_partitions_offset_ptr(0) = 0; } // Facets Mesh::type_iterator fit = mesh.firstType(spatial_dimension - 1, _not_ghost, _ek_not_defined); Mesh::type_iterator fend = mesh.lastType(spatial_dimension - 1, _not_ghost, _ek_not_defined); for(; fit != fend; ++fit) { ElementType type = *fit; UInt nb_element = mesh.getNbElement(type); partitions .alloc(nb_element, 1, type, _not_ghost); AKANTU_DEBUG_WARNING("Allocating partitions for " << type); CSR & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); ghost_partitions .alloc(0, 1, type, _ghost); AKANTU_DEBUG_WARNING("Allocating ghost_partitions for " << type); const Array< std::vector > & elem_to_subelem = mesh.getElementToSubelement(type, _not_ghost); for(UInt i(0); i < mesh.getNbElement(type, _not_ghost); ++i) { // Facet loop const std::vector & adjacent_elems = elem_to_subelem(i); Element min_elem; UInt min_part(std::numeric_limits::max()); std::set adjacent_parts; for(UInt j(0); j < adjacent_elems.size(); ++j) { UInt adjacent_elem_id = adjacent_elems[j].element; UInt adjacent_elem_part = partitions(adjacent_elems[j].type, adjacent_elems[j].ghost_type)(adjacent_elem_id); if(adjacent_elem_part < min_part) { min_part = adjacent_elem_part; min_elem = adjacent_elems[j]; } adjacent_parts.insert(adjacent_elem_part); } partitions(type, _not_ghost)(i) = min_part; CSR::iterator git = ghost_partitions_csr(min_elem.type, _not_ghost).begin(min_elem.element); CSR::iterator gend = ghost_partitions_csr(min_elem.type, _not_ghost).end(min_elem.element); for(; git != gend; ++git) { adjacent_parts.insert(*git); } adjacent_parts.erase(min_part); std::set::const_iterator pit = adjacent_parts.begin(); std::set::const_iterator pend = adjacent_parts.end(); for(; pit != pend; ++pit) { ghost_part_csr.getRows().push_back(*pit); ghost_part_csr.rowOffset(i)++; ghost_partitions(type, _ghost).push_back(*pit); } ghost_partitions_offset(type, _ghost)(i+1) = ghost_partitions_offset(type, _ghost)(i+1) + adjacent_elems.size(); } ghost_part_csr.countToCSR(); } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index eac2ce660..e75ffeb41 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,513 +1,516 @@ /** * @file material.hh * * @author Marco Vocialta * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Mother class for all materials * * @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 "aka_common.hh" #include "aka_memory.hh" #include "parser.hh" #include "data_accessor.hh" #include "material_parameters.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } __BEGIN_AKANTU__ /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : protected Memory, public DataAccessor, public Parsable, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register a material parameter template void registerParam(std::string name, T & variable, T default_value, ParamAccessType type, std::string description = ""); template void registerParam(std::string name, T & variable, ParamAccessType type, std::string description = ""); template void registerInternal(__attribute__((unused)) ByElementTypeArray & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// read parameter from file virtual bool parseParam(const std::string & key, const std::string & value, const ID & id); /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleResidual(GhostType ghost_type); virtual void UpdateStressesAtT(GhostType ghost_type){}; /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {}; virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(__attribute__((unused)) Real h, __attribute__((unused)) const Element & element = ElementNull) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ virtual void interpolateStress(ByElementTypeReal & result); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ virtual void initElementalFieldInterpolation(ByElementTypeReal & interpolation_points_coordinates); protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } template void computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type); /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(const Matrix & B, Matrix & Bvoigt, UInt nb_nodes_per_element) const; inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElement(); /// compute the coordinates of the quadrature points void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const; /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Array & field, Array & result, const GhostType ghost_type); /// template function to initialize the elemental field interpolation template void initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const GhostType ghost_type); /// build the coordinate matrix for the interpolation on elemental field template inline void buildElementalFieldInterpolationCoodinates(const Matrix & coordinates, Matrix & coordMatrix); /// get the size of the coordiante matrix used in the interpolation template inline UInt getSizeElementalFieldInterpolationCoodinates(GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for a quadrature point inline void computePotentialEnergyOnQuad(Matrix & grad_u, Matrix & sigma, Real & epot); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(ElementType type, UInt index, Vector & epot_on_quad_points); protected: /// allocate an internal vector template void initInternalArray(ByElementTypeArray & vect, UInt nb_component, bool temporary = false, ElementKind element_kind = _ek_regular); public: /// resize an internal vector template void resizeInternalArray(ByElementTypeArray & vect, ElementKind element_kind = _ek_regular) const; /* ------------------------------------------------------------------------ */ template inline void gradUToF(const Matrix & grad_u, Matrix & F); inline void rightCauchy(const Matrix & F, Matrix & C); inline void leftCauchy (const Matrix & F, Matrix & B); inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); template inline void packElementDataHelper(const ByElementTypeArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ByElementTypeArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; public: /* ------------------------------------------------------------------------ */ virtual inline void onElementsAdded(const Array & element_list, const NewElementsEvent & event); virtual inline void onElementsRemoved(const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event); protected: template void removeQuadraturePointsFromArrays(ByElementTypeArray & data, const ByElementTypeUInt & new_numbering); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); /// return the potential energy for the subset of elements contained by the material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained by the material virtual Real getEnergy(std::string energy_id); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); + AKANTU_GET_MACRO(Strain, strain, const ByElementTypeReal &); + AKANTU_GET_MACRO(Stress, stress, const ByElementTypeReal &); + AKANTU_GET_MACRO(ElementFilter, element_filter, const ByElementTypeUInt &); bool isNonLocal() const { return is_non_local; } const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); template inline T getParam(const ID & param) const; template inline void setParam(const ID & param, T value); bool isFiniteDeformation() const { return finite_deformation; } protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material ID id; /// Finite deformation bool finite_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// is the vector for potential energy initialized // bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// tell if using in non local mode or not bool is_non_local; /// spatial dimension UInt spatial_dimension; /// elemental field interpolation coordinates ByElementTypeReal interpolation_inverse_coordinates; /// elemental field interpolation points ByElementTypeReal interpolation_points_matrices; /// list of the paramters MaterialParameters params; private: /// boolean to know if the material has been initialized bool is_init; std::map internal_vectors_real; std::map internal_vectors_uint; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "material_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ this->stress(el_type, \ ghost_type).resize(this->strain(el_type, \ ghost_type).getSize()); \ \ Array::iterator< Matrix > stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ tangent_mat.resize(this->strain(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Array::iterator< Matrix > tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } /* -------------------------------------------------------------------------- */ #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> /* -------------------------------------------------------------------------- */ /* Material list */ /* -------------------------------------------------------------------------- */ // elastic materials #include "material_elastic.hh" #include "material_elastic_orthotropic.hh" #define AKANTU_CORE_MATERIAL_LIST \ ((2, (elastic , MaterialElastic ))) \ ((2, (elastic_orthotropic, MaterialElasticOrthotropic))) #if defined(AKANTU_EXTRA_MATERIALS) # include "material_extra_includes.hh" #else # define AKANTU_EXTRA_MATERIAL_LIST #endif #if defined(AKANTU_COHESIVE_ELEMENT) # include "material_cohesive_includes.hh" #else # define AKANTU_COHESIVE_MATERIAL_LIST #endif #if defined(AKANTU_DAMAGE_NON_LOCAL) # include "material_non_local_includes.hh" #else # define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #endif #define AKANTU_MATERIAL_LIST \ AKANTU_CORE_MATERIAL_LIST \ AKANTU_EXTRA_MATERIAL_LIST \ AKANTU_COHESIVE_MATERIAL_LIST \ AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 9ac64ce32..7db341daa 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1746 +1,1752 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the SolidMechanicsModel 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif -/* -------------------------------------------------------------------------- */ -__BEGIN_AKANTU__ - #ifdef AKANTU_USE_IOHELPER # include "dumper_iohelper_tmpl_homogenizing_field.hh" # include "dumper_iohelper_tmpl_material_internal_field.hh" #endif +/* -------------------------------------------------------------------------- */ +__BEGIN_AKANTU__ + /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), Dumpable(), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), element_index_by_material("element index by material", id), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->displacement_t = NULL; materials.clear(); mesh.registerEventHandler(*this); this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; if(solver) delete solver; if(mass_matrix) delete mass_matrix; if(velocity_damping_matrix) delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; if(jacobian_matrix) delete jacobian_matrix; if(dof_synchronizer) delete dof_synchronizer; if(synch_parallel) delete synch_parallel; AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(std::string material_file, AnalysisMethod analysis_method) { method = analysis_method; // initialize the model initModel(); // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; } // initialize the materials if(material_file != "") { readMaterials(material_file); initMaterials(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; synch_parallel = &createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); // synch_registry->registerSynchronizer(synch_parallel, _gst_smm_for_strain); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEMBoundary() { FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnControlPoints(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit() { AKANTU_DEBUG_IN(); // method = _explicit_dynamic; if (integrator) delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); std::stringstream sstr; sstr << id << ":increment_acceleration"; increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysFiniteDeformation() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":displacement_t"; displacement_t = &(alloc (sstr_disp_t.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); element_index_by_material.alloc(nb_element, 2, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); initBC(*this, *displacement, *force); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); changeLocalEquationNumberforPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->values; Real * position_val = mesh.getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); // f = f_ext - f_int // f = f_ext if (need_initialize) initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); updateStresses(); std::vector::iterator mat_it; #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant non local stresses"); synch_registry->waitEndSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant stresses"); // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateStresses() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::UpdateStressesAtT() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.UpdateStressesAtT(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (isExplicit()) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials updateStresses(); /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * accel_val = acceleration->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } boundary_val++; res_val++; mass_val++; accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv if(velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_lumped_mass) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *boundary, f_m2a); } else if (method == _explicit_consistent_mass) { solveDynamic(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & boundary, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * boundary_val = boundary.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*boundary_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; boundary_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { memcpy(increment->values, displacement->values, displacement->getSize()*displacement->getNbComponent()*sizeof(Real)); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_nodes = displacement->getSize(); for (UInt n = 0; n < nb_nodes; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *boundary, *increment_acceleration); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else - UInt nb_global_node = mesh.getNbGlobalNodes(); + UInt nb_global_nodes = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":jacobian_matrix"; - jacobian_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric, + jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); if (!isExplicit()) { - std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); - } #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS - solver->initialize(options); + if(solver) + solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; + if (!increment) setIncrementFlagOn(); + initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*boundary); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveDynamic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); - if(!increment) setIncrementFlagOn(); - updateResidualInternal(); solveDynamic(*increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic() { AKANTU_DEBUG_IN(); - UInt nb_nodes = displacement->getSize(); - UInt nb_degree_of_freedom = displacement->getNbComponent(); + AKANTU_DEBUG_INFO("Solving Ku = f"); + AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, + "You should first initialize the implicit solver and assemble the stiffness matrix"); - Array * normal_boundary = new Array (nb_nodes, nb_degree_of_freedom, "normal_boundary"); - normal_boundary->clear(); - UInt n_angles(0); - if(nb_degree_of_freedom==2) { - n_angles=1; - } - else if(nb_degree_of_freedom==3) { - n_angles=3; - } else { - AKANTU_DEBUG_ERROR("Invalid number of degrees if freedom in solveStatic."); - } + UInt nb_nodes = displacement->getSize(); + UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; + + jacobian_matrix->copyContent(*stiffness_matrix); + jacobian_matrix->applyBoundary(*boundary); - Array * Euler_angles = new Array (nb_nodes, n_angles, "Euler_angles"); - Euler_angles->clear(); - solveStatic(*normal_boundary, *Euler_angles); + increment->clear(); + + solver->setRHS(*residual); + solver->solve(*increment); + + Real * increment_val = increment->values; + Real * displacement_val = displacement->values; + bool * boundary_val = boundary->values; + + for (UInt j = 0; j < nb_degree_of_freedom; + ++j, ++displacement_val, ++increment_val, ++boundary_val) { + if (!(*boundary_val)) { + *displacement_val += *increment_val; + } + } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic(Array & boundary_normal, Array & EulerAngles) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); // if(method != _static) jacobian_matrix->copyContent(*stiffness_matrix); Array * residual_rotated = new Array (nb_nodes, nb_degree_of_freedom, "residual_rotated"); //stiffness_matrix->saveMatrix("stiffness_original.out"); jacobian_matrix->applyBoundaryNormal(boundary_normal, EulerAngles, *residual, (*stiffness_matrix).getA(), *residual_rotated); //jacobian_matrix->saveMatrix("stiffness_rotated_dir.out"); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual_rotated); delete residual_rotated; if (!increment) setIncrementFlagOn(); solver->solve(*increment); Matrix T(nb_degree_of_freedom, nb_degree_of_freedom); Matrix small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Matrix T_small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool constrain_ij = false; for (UInt j = 0; j < nb_degree_of_freedom; j++) { if (boundary_normal(n, j)) { constrain_ij = true; break; } } if (constrain_ij) { if (nb_degree_of_freedom == 2) { Real Theta = EulerAngles(n, 0); T(0, 0) = cos(Theta); T(0, 1) = -sin(Theta); T(1, 1) = cos(Theta); T(1, 0) = sin(Theta); } else if (nb_degree_of_freedom == 3) { Real Theta_x = EulerAngles(n, 0); Real Theta_y = EulerAngles(n, 1); Real Theta_z = EulerAngles(n, 2); T(0, 0) = cos(Theta_y) * cos(Theta_z); T(0, 1) = -cos(Theta_y) * sin(Theta_z); T(0, 2) = sin(Theta_y); T(1, 0) = cos(Theta_x) * sin(Theta_z) + cos(Theta_z) * sin(Theta_x) * sin(Theta_y); T(1, 1) = cos(Theta_x) * cos(Theta_z) - sin(Theta_x) * sin(Theta_y) * sin(Theta_z); T(1, 2) = -cos(Theta_y) * sin(Theta_x); T(2, 0) = sin(Theta_x) * sin(Theta_z) - cos(Theta_x) * cos(Theta_z) * sin(Theta_y); T(2, 1) = cos(Theta_z) * sin(Theta_x) + cos(Theta_x) * sin(Theta_y) * sin(Theta_z); T(2, 2) = cos(Theta_x) * cos(Theta_y); } small_rhs.clear(); T_small_rhs.clear(); for (UInt j = 0; j < nb_degree_of_freedom; j++) if(!(boundary_normal(n,j)) ) small_rhs(j,j)=increment_val[j]; T_small_rhs.mul(T,small_rhs); for (UInt j = 0; j < nb_degree_of_freedom; j++){ if(!(boundary_normal(n,j))){ for (UInt k = 0; k < nb_degree_of_freedom; k++) displacement_val[k]+=T_small_rhs(k,j); } } displacement_val += nb_degree_of_freedom; boundary_val += nb_degree_of_freedom; increment_val += nb_degree_of_freedom; } else { for (UInt j = 0; j < nb_degree_of_freedom; j++, ++displacement_val, ++increment_val, ++boundary_val) { if (!(*boundary_val)) { *displacement_val += *increment_val; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { if(!velocity_damping_matrix) velocity_damping_matrix = new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) { Real error; bool tmp = testConvergenceIncrement(tolerance, error); AKANTU_DEBUG_INFO("Norm of increment : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) { AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * boundary_val = boundary->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } boundary_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); error = norm[0] / norm[1]; return (error < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) { Real error; bool tmp = testConvergenceResidual(tolerance, error); AKANTU_DEBUG_INFO("Norm of residual : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *boundary, *increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; - increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0)); + increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0.)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::iterator< Vector > eibm = element_index_by_material(*it, ghost_type).begin(2); Array X(0, nb_nodes_per_element*spatial_dimension); FEM::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::iterator< Matrix > X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) { elem.element = el; Real el_size = getFEM().getElementInradius(*X_el, *it); Real el_dt = mat_val[(*eibm)(1)]->getStableTimeStep(el_size, elem); min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEM().interpolateOnQuadraturePoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::iterator< Vector > vit = vel_on_quad.begin(spatial_dimension); Array::iterator< Vector > vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[element_index_by_material(type)(index, 1)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEM().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = boundary->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work"){ return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(energy_id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index){ AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } std::vector::iterator mat_it; Vector mat = element_index_by_material(type, _not_ghost).begin(2)[index]; Real energy = materials[mat(1)]->getEnergy(energy_id, type, mat(0)); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if(displacement) displacement->resize(nb_nodes); if(mass ) mass ->resize(nb_nodes); if(velocity ) velocity ->resize(nb_nodes); if(acceleration) acceleration->resize(nb_nodes); if(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(boundary ) boundary ->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); if(current_position) current_position->resize(nb_nodes); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } if (method != _explicit_lumped_mass) { delete stiffness_matrix; delete jacobian_matrix; delete solver; SolverOptions solver_options; initImplicit((method == _implicit_dynamic), solver_options); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); /// \todo have rules to choose the correct material UInt mat_id = 0; UInt * mat_id_vect = NULL; try { const NewMaterialElementsEvent & event_mat = dynamic_cast(event); mat_id_vect = event_mat.getMaterialList().storage(); } catch(...) { } for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; // element_material(elem.type, elem.ghost_type).push_back(UInt(0)); if(mat_id_vect) mat_id = mat_id_vect[el]; Material & mat = *materials[mat_id]; UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type); UInt id[2]; id[0] = mat_index; id[1] = 0; element_index_by_material(elem.type, elem.ghost_type).push_back(id); } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } if(method != _explicit_lumped_mass) AKANTU_DEBUG_TO_IMPLEMENT(); assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event) { // MeshUtils::purifyMesh(mesh); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering); if(mass ) mesh.removeNodesFromArray(*mass , new_numbering); if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if(force ) mesh.removeNodesFromArray(*force , new_numbering); if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); if(boundary ) mesh.removeNodesFromArray(*boundary , new_numbering); if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromArray(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ internalAddDumpFieldToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField<>(mesh, spatial_dimension, _not_ghost, _ek_regular)); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, spatial_dimension, _not_ghost, _ek_regular)); } else { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::HomogenizedField(*this, field_id, spatial_dimension, _not_ghost, _ek_regular)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpBoundaryField(const std::string & field_id, const std::string & boundary_name) { SubBoundary & boundary_ref = mesh.getSubBoundary(boundary_name); this->addDumpBoundaryFieldToDumper(boundary_ref.getDefaultDumperName(), field_id, boundary_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpBoundaryFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & boundary_name) { #ifdef AKANTU_USE_IOHELPER SubBoundary & boundary_ref = mesh.getSubBoundary(boundary_name); const Array * nodal_filter = &(boundary_ref.getNodes()); const ByElementTypeArray * elemental_filter = &(boundary_ref.getElements()); #define ADD_FIELD(field, type) \ boundary_ref.registerField(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { boundary_ref.registerField(dumper_name, field_id, new DumperIOHelper::ElementPartitionField(mesh, spatial_dimension, _not_ghost, _ek_regular, elemental_filter)); } else if(field_id == "element_index_by_material") { boundary_ref.registerField(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, spatial_dimension, _not_ghost, _ek_regular, elemental_filter)); } else { boundary_ref.registerField(dumper_name, field_id, new DumperIOHelper::HomogenizedField(*this, field_id, spatial_dimension, _not_ghost, _ek_regular, elemental_filter)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpBoundaryField(const std::string & field_id, const std::string & boundary_name) { SubBoundary & boundary_ref = mesh.getSubBoundary(boundary_name); this->removeDumpBoundaryFieldFromDumper(boundary_ref.getDefaultDumperName(), field_id, boundary_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpBoundaryFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & boundary_name) { SubBoundary & boundary_ref = mesh.getSubBoundary(boundary_name); boundary_ref.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpBoundaryFieldVector(const std::string & field_id, const std::string & boundary_name) { SubBoundary & boundary_ref = mesh.getSubBoundary(boundary_name); this->addDumpBoundaryFieldVectorToDumper(boundary_ref.getDefaultDumperName(), field_id, boundary_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpBoundaryFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & boundary_name) { #ifdef AKANTU_USE_IOHELPER SubBoundary & boundary_ref = mesh.getSubBoundary(boundary_name); const Array * nodal_filter = &(boundary_ref.getNodes()); const ByElementTypeArray * elemental_filter = &(boundary_ref.getElements()); #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter); \ f->setPadding(3); \ boundary_ref.registerField(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular, elemental_filter); field->setPadding(3); boundary_ref.registerField(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER if(field_id == "stress") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else if(field_id == "strain") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEM().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; element_index_by_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/solver/solver_mumps.cc b/src/solver/solver_mumps.cc index 638882948..bbcf73d66 100644 --- a/src/solver/solver_mumps.cc +++ b/src/solver/solver_mumps.cc @@ -1,314 +1,335 @@ /** * @file solver_mumps.cc * * @author Nicolas Richart * * @date Mon Dec 13 10:48:06 2010 * * @brief implem of SolverMumps 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 . * * @section DESCRIPTION * * @subsection Ctrl_param Control parameters * * ICNTL(1), * ICNTL(2), * ICNTL(3) : output streams for error, diagnostics, and global messages * * ICNTL(4) : verbose level : 0 no message - 4 all messages * * ICNTL(5) : type of matrix, 0 assembled, 1 elementary * * ICNTL(6) : control the permutation and scaling(default 7) see mumps doc for * more information * * ICNTL(7) : determine the pivot order (default 7) see mumps doc for more * information * * ICNTL(8) : describe the scaling method used * * ICNTL(9) : 1 solve A x = b, 0 solve At x = b * * ICNTL(10) : number of iterative refinement when NRHS = 1 * * ICNTL(11) : > 0 return statistics * * ICNTL(12) : only used for SYM = 2, ordering strategy * * ICNTL(13) : * * ICNTL(14) : percentage of increase of the estimated working space * * ICNTL(15-17) : not used * * ICNTL(18) : only used if ICNTL(5) = 0, 0 matrix centralized, 1 structure on * host and mumps give the mapping, 2 structure on host and distributed matrix * for facto, 3 distributed matrix * * ICNTL(19) : > 0, Shur complement returned * * ICNTL(20) : 0 rhs dense, 1 rhs sparse * * ICNTL(21) : 0 solution in rhs, 1 solution distributed in ISOL_loc and SOL_loc * allocated by user * * ICNTL(22) : 0 in-core, 1 out-of-core * * ICNTL(23) : maximum memory allocatable by mumps pre proc * * ICNTL(24) : controls the detection of "null pivot rows" * * ICNTL(25) : * * ICNTL(26) : * * ICNTL(27) : * * ICNTL(28) : 0 automatic choice, 1 sequential analysis, 2 parallel analysis * * ICNTL(29) : 0 automatic choice, 1 PT-Scotch, 2 ParMetis */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_MPI #include "static_communicator_mpi.hh" #endif #include "solver_mumps.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolverMumps::SolverMumps(SparseMatrix & matrix, const ID & id, const MemoryID & memory_id) : Solver(matrix, id, memory_id), is_mumps_data_initialized(false), rhs_is_local(true) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_MPI parallel_method = SolverMumpsOptions::_fully_distributed; #else //AKANTU_USE_MPI parallel_method = SolverMumpsOptions::_master_slave_distributed; #endif //AKANTU_USE_MPI CommunicatorEventHandler & comm_event_handler = *this; communicator.registerEventHandler(comm_event_handler); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolverMumps::~SolverMumps() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::destroyMumpsData() { AKANTU_DEBUG_IN(); if(is_mumps_data_initialized) { mumps_data.job = _smj_destroy; // destroy dmumps_c(&mumps_data); is_mumps_data_initialized = false; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::onCommunicatorFinalize(const StaticCommunicator & comm) { AKANTU_DEBUG_IN(); try{ #if defined(AKANTU_USE_MPI) const StaticCommunicatorMPI & comm_mpi = dynamic_cast(comm.getRealStaticCommunicator()); if(mumps_data.comm_fortran == MPI_Comm_c2f(comm_mpi.getMPICommunicator())) #endif destroyMumpsData(); } catch(...) {} AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +std::ostream & operator <<(std::ostream & stream, const DMUMPS_STRUC_C & _this) { + stream << "DMUMPS Data [" << std::endl; + stream << " + job : " << _this.job << std::endl; + stream << " + par : " << _this.par << std::endl; + stream << " + sym : " << _this.sym << std::endl; + stream << " + comm_fortran : " << _this.comm_fortran << std::endl; + stream << " + nz : " << _this.nz << std::endl; + stream << " + irn : " << _this.irn << std::endl; + stream << " + jcn : " << _this.jcn << std::endl; + stream << " + nz_loc : " << _this.nz_loc << std::endl; + stream << " + irn_loc : " << _this.irn_loc << std::endl; + stream << " + jcn_loc : " << _this.jcn_loc << std::endl; + stream << "]"; + return stream; +} + + /* -------------------------------------------------------------------------- */ void SolverMumps::initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method) { switch(parallel_method) { case SolverMumpsOptions::_fully_distributed: icntl(18) = 3; //fully distributed icntl(28) = 0; //automatic choice mumps_data.nz_loc = matrix->getNbNonZero(); mumps_data.irn_loc = matrix->getIRN().values; mumps_data.jcn_loc = matrix->getJCN().values; break; case SolverMumpsOptions::_master_slave_distributed: if(prank == 0) { mumps_data.nz = matrix->getNbNonZero(); mumps_data.irn = matrix->getIRN().values; mumps_data.jcn = matrix->getJCN().values; } else { mumps_data.nz = 0; mumps_data.irn = NULL; mumps_data.jcn = NULL; icntl(18) = 0; //centralized icntl(28) = 0; //sequential analysis } break; + default: + AKANTU_DEBUG_ERROR("This case should not happen!!"); } } /* -------------------------------------------------------------------------- */ void SolverMumps::initialize(SolverOptions & options) { AKANTU_DEBUG_IN(); mumps_data.par = 1; if(SolverMumpsOptions * opt = dynamic_cast(&options)) { if(opt->parallel_method == SolverMumpsOptions::_master_slave_distributed) { mumps_data.par = 0; } } mumps_data.sym = 2 * (matrix->getSparseMatrixType() == _symmetric); prank = communicator.whoAmI(); #ifdef AKANTU_USE_MPI mumps_data.comm_fortran = MPI_Comm_c2f(dynamic_cast(communicator.getRealStaticCommunicator()).getMPICommunicator()); #endif if(AKANTU_DEBUG_TEST(dblTrace)) { icntl(1) = 2; icntl(2) = 2; icntl(3) = 2; icntl(4) = 4; } mumps_data.job = _smj_initialize; //initialize dmumps_c(&mumps_data); + + is_mumps_data_initialized = true; /* ------------------------------------------------------------------------ */ UInt size = matrix->getSize(); if(prank == 0) { std::stringstream sstr_rhs; sstr_rhs << id << ":rhs"; - rhs = &(alloc(sstr_rhs.str(), size, 1, REAL_INIT_VALUE)); + rhs = &(alloc(sstr_rhs.str(), size, 1, 0.)); } else { rhs = NULL; } /// No outputs icntl(1) = 0; icntl(2) = 0; icntl(3) = 0; icntl(4) = 0; - mumps_data.nz_alloc = 0; if (AKANTU_DEBUG_TEST(dblDump)) icntl(4) = 4; mumps_data.n = size; if(AKANTU_DEBUG_TEST(dblDump)) { strcpy(mumps_data.write_problem, "mumps_matrix.mtx"); } /* ------------------------------------------------------------------------ */ // Default Scaling icntl(8) = 77; icntl(5) = 0; // Assembled matrix SolverMumpsOptions * opt = dynamic_cast(&options); if(opt) parallel_method = opt->parallel_method; initMumpsData(parallel_method); mumps_data.job = _smj_analyze; //analyze dmumps_c(&mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::setRHS(Array & rhs) { if(prank == 0) { matrix->getDOFSynchronizer().gather(rhs, 0, this->rhs); } else { matrix->getDOFSynchronizer().gather(rhs, 0); } } /* -------------------------------------------------------------------------- */ void SolverMumps::solve() { AKANTU_DEBUG_IN(); if(parallel_method == SolverMumpsOptions::_fully_distributed) mumps_data.a_loc = matrix->getA().values; else if(prank == 0) { mumps_data.a = matrix->getA().values; } if(prank == 0) { mumps_data.rhs = rhs->values; } /// Default centralized dense second member icntl(20) = 0; icntl(21) = 0; mumps_data.job = _smj_factorize_solve; //solve dmumps_c(&mumps_data); AKANTU_DEBUG_ASSERT(info(1) != -10, "Singular matrix"); AKANTU_DEBUG_ASSERT(info(1) == 0, "Error in mumps during solve process, check mumps user guide INFO(1) =" << info(1)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::solve(Array & solution) { AKANTU_DEBUG_IN(); solve(); if(prank == 0) { matrix->getDOFSynchronizer().scatter(solution, 0, this->rhs); } else { matrix->getDOFSynchronizer().scatter(solution, 0); } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/solver/solver_mumps.hh b/src/solver/solver_mumps.hh index cb26dd4f9..505fc832d 100644 --- a/src/solver/solver_mumps.hh +++ b/src/solver/solver_mumps.hh @@ -1,165 +1,165 @@ /** * @file solver_mumps.hh * * @author Nicolas Richart * * @date Mon Dec 13 10:48:06 2010 * * @brief Solver class implementation for the mumps solver * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_MUMPS_HH__ #define __AKANTU_SOLVER_MUMPS_HH__ #include #include "solver.hh" #include "static_communicator.hh" __BEGIN_AKANTU__ class SolverMumpsOptions : public SolverOptions { public: enum ParallelMethod { _fully_distributed, _master_slave_distributed }; SolverMumpsOptions(ParallelMethod parallel_method = _fully_distributed) : SolverOptions(), parallel_method(parallel_method) { } virtual void niceFunctionWhichDoesNothing() { SolverOptions::niceFunctionWhichDoesNothing(); AKANTU_DEBUG_ERROR("Nothing!!! (TWICE)"); }; private: friend class SolverMumps; ParallelMethod parallel_method; }; class SolverMumps : public Solver, public CommunicatorEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SolverMumps(SparseMatrix & sparse_matrix, const ID & id = "solver_mumps", const MemoryID & memory_id = 0); virtual ~SolverMumps(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile and do the analysis part void initialize(SolverOptions & options = _solver_no_options); void initializeSlave(SolverOptions & options = _solver_no_options); /// factorize and solve the system void solve(Array & solution); void solve(); void solveSlave(); virtual void setRHS(Array & rhs); /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; virtual void onCommunicatorFinalize(const StaticCommunicator & communicator); private: void destroyMumpsData(); inline Int & icntl(UInt i) { return mumps_data.icntl[i - 1]; } inline Int & info(UInt i) { return mumps_data.info[i - 1]; } void initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// mumps data DMUMPS_STRUC_C mumps_data; /// specify if the mumps_data are initialized or not bool is_mumps_data_initialized; UInt prank; /* ------------------------------------------------------------------------ */ /* Local types */ /* ------------------------------------------------------------------------ */ private: SolverMumpsOptions::ParallelMethod parallel_method; bool rhs_is_local; enum SolverMumpsJob { _smj_initialize = -1, _smj_analyze = 1, _smj_factorize = 2, _smj_solve = 3, _smj_analyze_factorize = 4, _smj_factorize_solve = 5, _smj_complete = 6, // analyze, factorize, solve _smj_destroy = -2 }; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "solver_mumps_inline_impl.cc" - +std::ostream & operator <<(std::ostream & stream, const DMUMPS_STRUC_C & _this); /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const SolverMumps & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_SOLVER_MUMPS_HH__ */ diff --git a/src/synchronizer/distributed_synchronizer.cc b/src/synchronizer/distributed_synchronizer.cc index 919936693..c08f93fa1 100644 --- a/src/synchronizer/distributed_synchronizer.cc +++ b/src/synchronizer/distributed_synchronizer.cc @@ -1,1162 +1,1168 @@ /** * @file distributed_synchronizer.cc * * @author Nicolas Richart * * @date Thu Jun 16 16:36:52 2011 * * @brief implementation of a communicator using a static_communicator for real * send/receive * * @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 "aka_common.hh" #include "distributed_synchronizer.hh" #include "static_communicator.hh" #include "mesh_utils.hh" #include "mesh_data.hh" /* -------------------------------------------------------------------------- */ #include #include #include #if defined(AKANTU_DEBUG_TOOLS) # include "aka_debug_tools.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ DistributedSynchronizer::DistributedSynchronizer(Mesh & mesh, SynchronizerID id, MemoryID memory_id) : Synchronizer(id, memory_id), mesh(mesh), static_communicator(&StaticCommunicator::getStaticCommunicator()) { AKANTU_DEBUG_IN(); nb_proc = static_communicator->getNbProc(); rank = static_communicator->whoAmI(); send_element = new Array[nb_proc]; recv_element = new Array[nb_proc]; mesh.registerEventHandler(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DistributedSynchronizer::~DistributedSynchronizer() { AKANTU_DEBUG_IN(); for (UInt p = 0; p < nb_proc; ++p) { send_element[p].clear(); recv_element[p].clear(); } delete [] send_element; delete [] recv_element; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DistributedSynchronizer * DistributedSynchronizer:: createDistributedSynchronizerMesh(Mesh & mesh, const MeshPartition * partition, UInt root, SynchronizerID id, MemoryID memory_id) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); DistributedSynchronizer & communicator = *(new DistributedSynchronizer(mesh, id, memory_id)); if(nb_proc == 1) return &communicator; UInt * local_connectivity = NULL; UInt * local_partitions = NULL; Array * old_nodes = mesh.getNodesGlobalIdsPointer(); old_nodes->resize(0); Array * nodes = mesh.getNodesPointer(); UInt spatial_dimension = nodes->getNbComponent(); /* ------------------------------------------------------------------------ */ /* Local (rank == root) */ /* ------------------------------------------------------------------------ */ if(my_rank == root) { + mesh.nb_global_nodes = mesh.nodes->getSize(); + AKANTU_DEBUG_ASSERT(partition->getNbPartition() == nb_proc, "The number of partition does not match the number of processors: " << partition->getNbPartition() << " != " << nb_proc); /** * connectivity and communications scheme construction */ Mesh::type_iterator it = mesh.firstType(_all_dimensions, _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, _not_ghost, _ek_not_defined); UInt count = 0; /* --- MAIN LOOP ON TYPES --- */ for(; it != end; ++it) { ElementType type = *it; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(*it); UInt nb_local_element[nb_proc]; UInt nb_ghost_element[nb_proc]; UInt nb_element_to_send[nb_proc]; memset(nb_local_element, 0, nb_proc*sizeof(UInt)); memset(nb_ghost_element, 0, nb_proc*sizeof(UInt)); memset(nb_element_to_send, 0, nb_proc*sizeof(UInt)); const Array & partition_num = partition->getPartition(type, _not_ghost); const CSR & ghost_partition = partition->getGhostPartitionCSR()(type, _not_ghost); // const Array & ghost_partition = partition->getGhostPartition(type, _ghost); // const Array & ghost_partition_offset = partition->getGhostPartitionOffset(type, _ghost); /* -------------------------------------------------------------------- */ /// constructing the reordering structures for (UInt el = 0; el < nb_element; ++el) { nb_local_element[partition_num(el)]++; for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { nb_ghost_element[*part]++; } nb_element_to_send[partition_num(el)] += ghost_partition.getNbCols(el) + 1; } /// allocating buffers UInt * buffers[nb_proc]; UInt * buffers_tmp[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { UInt size = nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]); buffers[p] = new UInt[size]; buffers_tmp[p] = buffers[p]; } /// copying the local connectivity UInt * conn_val = mesh.getConnectivity(type, _not_ghost).values; for (UInt el = 0; el < nb_element; ++el) { memcpy(buffers_tmp[partition_num(el)], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[partition_num(el)] += nb_nodes_per_element; } /// copying the connectivity of ghost element for (UInt el = 0; el < nb_element; ++el) { for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { UInt proc = *part; memcpy(buffers_tmp[proc], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[proc] += nb_nodes_per_element; } } /// tag info std::vector tag_names; mesh.getMeshData().getTagNames(tag_names, type); // Make sure the tags are sorted (or at least not in random order), // because they come from a map !! std::sort(tag_names.begin(), tag_names.end()); UInt nb_tags = tag_names.size(); /* -------->>>>-SIZE + CONNECTIVITY------------------------------------ */ /// send all connectivity and ghost information to all processors std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { UInt size[5]; size[0] = (UInt) type; size[1] = nb_local_element[p]; size[2] = nb_ghost_element[p]; size[3] = nb_element_to_send[p]; size[4] = nb_tags; AKANTU_DEBUG_INFO("Sending connectivities informations to proc " << p << " TAG("<< Tag::genTag(my_rank, count, TAG_SIZES) <<")"); comm.send(size, 5, p, Tag::genTag(my_rank, count, TAG_SIZES)); AKANTU_DEBUG_INFO("Sending connectivities to proc " << p << " TAG("<< Tag::genTag(my_rank, count, TAG_CONNECTIVITY) <<")"); requests.push_back(comm.asyncSend(buffers[p], nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]), p, Tag::genTag(my_rank, count, TAG_CONNECTIVITY))); } else { local_connectivity = buffers[p]; } } /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element[root], nb_ghost_element[root], type, *old_nodes); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete [] buffers[p]; } /* -------------------------------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { buffers[p] = new UInt[nb_ghost_element[p] + nb_element_to_send[p]]; buffers_tmp[p] = buffers[p]; } /// splitting the partition information to send them to processors UInt count_by_proc[nb_proc]; memset(count_by_proc, 0, nb_proc*sizeof(UInt)); for (UInt el = 0; el < nb_element; ++el) { *(buffers_tmp[partition_num(el)]++) = ghost_partition.getNbCols(el); UInt i(0); for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { *(buffers_tmp[partition_num(el)]++) = *part; } } for (UInt el = 0; el < nb_element; ++el) { UInt i(0); for (CSR::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { *(buffers_tmp[*part]++) = partition_num(el); } } /* -------->>>>-PARTITIONS--------------------------------------------- */ /// last data to compute the communication scheme for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { AKANTU_DEBUG_INFO("Sending partition informations to proc " << p << " TAG("<< Tag::genTag(my_rank, count, TAG_PARTITIONS) <<")"); requests.push_back(comm.asyncSend(buffers[p], nb_element_to_send[p] + nb_ghost_element[p], p, Tag::genTag(my_rank, count, TAG_PARTITIONS))); } else { local_partitions = buffers[p]; } } if(Mesh::getSpatialDimension(type) == mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); communicator.fillCommunicationScheme(local_partitions, nb_local_element[root], nb_ghost_element[root], type); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete [] buffers[p]; } /* -------------------------------------------------------------------- */ /// send data assossiated to the mesh /* -------->>>>-TAGS------------------------------------------------- */ if(nb_tags != 0) { UInt mesh_data_sizes_buffer_length; MeshData & mesh_data = mesh.getMeshData(); { // Sending information about the tags in mesh_data: name, data type and // number of components of the underlying array associated to the current type DynamicCommunicationBuffer mesh_data_sizes_buffer; std::vector::const_iterator names_it = tag_names.begin(); std::vector::const_iterator names_end = tag_names.end(); for(;names_it != names_end; ++names_it) { mesh_data_sizes_buffer << *names_it; mesh_data_sizes_buffer << mesh_data.getTypeCode(*names_it); mesh_data_sizes_buffer << mesh_data.getNbComponent(*names_it, type); } mesh_data_sizes_buffer_length = mesh_data_sizes_buffer.getSize(); AKANTU_DEBUG_INFO("Broadcasting the size of the information about the mesh data tags: (" << mesh_data_sizes_buffer_length << ")." ); comm.broadcast(&mesh_data_sizes_buffer_length, 1, root); AKANTU_DEBUG_INFO("Broadcasting the information about the mesh data tags, addr " << (void*)mesh_data_sizes_buffer.storage()); if(mesh_data_sizes_buffer_length !=0) comm.broadcast(mesh_data_sizes_buffer.storage(), mesh_data_sizes_buffer.getSize(), root); } if(mesh_data_sizes_buffer_length !=0) { //Sending the actual data to each processor DynamicCommunicationBuffer buffers[nb_proc]; std::vector::const_iterator names_it = tag_names.begin(); std::vector::const_iterator names_end = tag_names.end(); // Loop over each tag for the current type for(;names_it != names_end; ++names_it) { // Type code of the current tag (i.e. the tag named *names_it) communicator.fillTagBuffer(mesh_data, buffers, *names_it, type, partition_num, ghost_partition); - - // Reinitializing the mesh data on the master - communicator.populateMeshData(mesh_data, - buffers[root], - *names_it, - type, - mesh_data.getTypeCode(*names_it), - mesh_data.getNbComponent(*names_it, type), - nb_local_element[root], - nb_ghost_element[root]); } std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { AKANTU_DEBUG_INFO("Sending " << buffers[p].getSize() << " bytes of mesh data to proc " << p << " TAG("<< Tag::genTag(my_rank, count, TAG_MESH_DATA) <<")"); requests.push_back(comm.asyncSend(buffers[p].storage(), buffers[p].getSize(), p, Tag::genTag(my_rank, count, TAG_MESH_DATA))); } } + names_it = tag_names.begin(); + // Loop over each tag for the current type + for(;names_it != names_end; ++names_it) { + // Reinitializing the mesh data on the master + communicator.populateMeshData(mesh_data, + buffers[root], + *names_it, + type, + mesh_data.getTypeCode(*names_it), + mesh_data.getNbComponent(*names_it, type), + nb_local_element[root], + nb_ghost_element[root]); + } + // Nécsesaire ? oui comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); } } ++count; } /* -------->>>>-SIZE----------------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { UInt size[5]; size[0] = (UInt) _not_defined; size[1] = 0; size[2] = 0; size[3] = 0; size[4] = 0; AKANTU_DEBUG_INFO("Sending empty connectivities informations to proc " << p << " TAG("<< Tag::genTag(my_rank, count, TAG_SIZES) <<")"); comm.send(size, 5, p, Tag::genTag(my_rank, count, TAG_SIZES)); } } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ /** * Nodes coordinate construction and synchronization */ std::multimap< UInt, std::pair > nodes_to_proc; /// get the list of nodes to send and send them Real * local_nodes = NULL; UInt nb_nodes_per_proc[nb_proc]; UInt * nodes_per_proc[nb_proc]; /* --------<<<<-NB_NODES + NODES----------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { UInt nb_nodes = 0; // UInt * buffer; if(p != root) { AKANTU_DEBUG_INFO("Receiving number of nodes from proc " << p << " TAG("<< Tag::genTag(p, 0, TAG_NB_NODES) <<")"); comm.receive(&nb_nodes, 1, p, Tag::genTag(p, 0, TAG_NB_NODES)); nodes_per_proc[p] = new UInt[nb_nodes]; nb_nodes_per_proc[p] = nb_nodes; AKANTU_DEBUG_INFO("Receiving list of nodes from proc " << p << " TAG("<< Tag::genTag(p, 0, TAG_NODES) <<")"); comm.receive(nodes_per_proc[p], nb_nodes, p, Tag::genTag(p, 0, TAG_NODES)); } else { nb_nodes = old_nodes->getSize(); nb_nodes_per_proc[p] = nb_nodes; nodes_per_proc[p] = old_nodes->values; } /// get the coordinates for the selected nodes Real * nodes_to_send = new Real[nb_nodes * spatial_dimension]; Real * nodes_to_send_tmp = nodes_to_send; for (UInt n = 0; n < nb_nodes; ++n) { memcpy(nodes_to_send_tmp, nodes->values + spatial_dimension * nodes_per_proc[p][n], spatial_dimension * sizeof(Real)); // nodes_to_proc.insert(std::make_pair(buffer[n], std::make_pair(p, n))); nodes_to_send_tmp += spatial_dimension; } /* -------->>>>-COORDINATES-------------------------------------------- */ if(p != root) { /// send them for distant processors AKANTU_DEBUG_INFO("Sending coordinates to proc " << p << " TAG("<< Tag::genTag(my_rank, 0, TAG_COORDINATES) <<")"); comm.send(nodes_to_send, nb_nodes * spatial_dimension, p, Tag::genTag(my_rank, 0, TAG_COORDINATES)); delete [] nodes_to_send; } else { /// save them for local processor local_nodes = nodes_to_send; } } /// construct the local nodes coordinates UInt nb_nodes = old_nodes->getSize(); nodes->resize(nb_nodes); memcpy(nodes->values, local_nodes, nb_nodes * spatial_dimension * sizeof(Real)); delete [] local_nodes; Array * nodes_type_per_proc[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { nodes_type_per_proc[p] = new Array(nb_nodes_per_proc[p]); } communicator.fillNodesType(mesh); /* --------<<<<-NODES_TYPE-1--------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { AKANTU_DEBUG_INFO("Receiving first nodes types from proc " << p << " TAG("<< Tag::genTag(my_rank, count, TAG_NODES_TYPE) <<")"); comm.receive(nodes_type_per_proc[p]->values, nb_nodes_per_proc[p], p, Tag::genTag(p, 0, TAG_NODES_TYPE)); } else { nodes_type_per_proc[p]->copy(mesh.getNodesType()); } for (UInt n = 0; n < nb_nodes_per_proc[p]; ++n) { if((*nodes_type_per_proc[p])(n) == -2) nodes_to_proc.insert(std::make_pair(nodes_per_proc[p][n], std::make_pair(p, n))); } } std::multimap< UInt, std::pair >::iterator it_node; std::pair< std::multimap< UInt, std::pair >::iterator, std::multimap< UInt, std::pair >::iterator > it_range; for (UInt i = 0; i < mesh.nb_global_nodes; ++i) { it_range = nodes_to_proc.equal_range(i); if(it_range.first == nodes_to_proc.end() || it_range.first->first != i) continue; UInt node_type = (it_range.first)->second.first; for (it_node = it_range.first; it_node != it_range.second; ++it_node) { UInt proc = it_node->second.first; UInt node = it_node->second.second; if(proc != node_type) nodes_type_per_proc[proc]->values[node] = node_type; } } /* -------->>>>-NODES_TYPE-2--------------------------------------------- */ std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { AKANTU_DEBUG_INFO("Sending nodes types to proc " << p << " TAG("<< Tag::genTag(my_rank, 0, TAG_NODES_TYPE) <<")"); requests.push_back(comm.asyncSend(nodes_type_per_proc[p]->values, nb_nodes_per_proc[p], p, Tag::genTag(my_rank, 0, TAG_NODES_TYPE))); } else { mesh.getNodesTypePointer()->copy(*nodes_type_per_proc[p]); } } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { if(p != root) delete [] nodes_per_proc[p]; delete nodes_type_per_proc[p]; } /* ---------------------------------------------------------------------- */ /* Distant (rank != root) */ /* ---------------------------------------------------------------------- */ } else { /** * connectivity and communications scheme construction on distant processors */ ElementType type = _not_defined; UInt count = 0; do { /* --------<<<<-SIZE--------------------------------------------------- */ UInt size[5] = { 0 }; comm.receive(size, 5, root, Tag::genTag(root, count, TAG_SIZES)); type = (ElementType) size[0]; UInt nb_local_element = size[1]; UInt nb_ghost_element = size[2]; UInt nb_element_to_send = size[3]; UInt nb_tags = size[4]; if(type != _not_defined) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /* --------<<<<-CONNECTIVITY----------------------------------------- */ local_connectivity = new UInt[(nb_local_element + nb_ghost_element) * nb_nodes_per_element]; AKANTU_DEBUG_INFO("Receiving connectivities from proc " << root); comm.receive(local_connectivity, nb_nodes_per_element * (nb_local_element + nb_ghost_element), root, Tag::genTag(root, count, TAG_CONNECTIVITY)); AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element, nb_ghost_element, type, *old_nodes); delete [] local_connectivity; /* --------<<<<-PARTITIONS--------------------------------------------- */ local_partitions = new UInt[nb_element_to_send + nb_ghost_element * 2]; AKANTU_DEBUG_INFO("Receiving partition informations from proc " << root); comm.receive(local_partitions, nb_element_to_send + nb_ghost_element * 2, root, Tag::genTag(root, count, TAG_PARTITIONS)); if(Mesh::getSpatialDimension(type) == mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); communicator.fillCommunicationScheme(local_partitions, nb_local_element, nb_ghost_element, type); } delete [] local_partitions; /* --------<<<<-TAGS------------------------------------------------- */ if(nb_tags != 0) { UInt mesh_data_sizes_buffer_length; CommunicationBuffer mesh_data_sizes_buffer; MeshData & mesh_data = mesh.getMeshData(); AKANTU_DEBUG_INFO("Receiving the size of the information about the mesh data tags."); comm.broadcast(&mesh_data_sizes_buffer_length, 1, root); if(mesh_data_sizes_buffer_length != 0) { mesh_data_sizes_buffer.resize(mesh_data_sizes_buffer_length); AKANTU_DEBUG_INFO("Receiving the information about the mesh data tags, addr " << (void*)mesh_data_sizes_buffer.storage()); comm.broadcast(mesh_data_sizes_buffer.storage(), mesh_data_sizes_buffer_length, root); AKANTU_DEBUG_INFO("Size of the information about the mesh data: " << mesh_data_sizes_buffer_length); std::vector tag_names; std::vector tag_type_codes; std::vector tag_nb_component; tag_names.resize(nb_tags); tag_type_codes.resize(nb_tags); tag_nb_component.resize(nb_tags); CommunicationBuffer mesh_data_buffer; UInt type_code_int; for(UInt i(0); i < nb_tags; ++i) { mesh_data_sizes_buffer >> tag_names[i]; mesh_data_sizes_buffer >> type_code_int; tag_type_codes[i] = static_cast(type_code_int); mesh_data_sizes_buffer >> tag_nb_component[i]; } std::vector::const_iterator names_it = tag_names.begin(); std::vector::const_iterator names_end = tag_names.end(); CommunicationStatus mesh_data_comm_status; AKANTU_DEBUG_INFO("Checking size of data to receive for mesh data TAG(" << Tag::genTag(root, count, TAG_MESH_DATA) << ")"); comm.probe(root, Tag::genTag(root, count, TAG_MESH_DATA), mesh_data_comm_status); UInt mesh_data_buffer_size(mesh_data_comm_status.getSize()); AKANTU_DEBUG_INFO("Receiving " << mesh_data_buffer_size << " bytes of mesh data TAG(" << Tag::genTag(root, count, TAG_MESH_DATA) << ")"); mesh_data_buffer.resize(mesh_data_buffer_size); comm.receive(mesh_data_buffer.storage(), mesh_data_buffer_size, root, Tag::genTag(root, count, TAG_MESH_DATA)); // Loop over each tag for the current type UInt k(0); for(; names_it != names_end; ++names_it, ++k) { communicator.populateMeshData(mesh_data, mesh_data_buffer, *names_it, type, tag_type_codes[k], tag_nb_component[k], nb_local_element, nb_ghost_element); } } } } ++count; } while(type != _not_defined); /** * Nodes coordinate construction and synchronization on distant processors */ /* -------->>>>-NB_NODES + NODES----------------------------------------- */ AKANTU_DEBUG_INFO("Sending list of nodes to proc " << root); UInt nb_nodes = old_nodes->getSize(); comm.send(&nb_nodes, 1, root, Tag::genTag(my_rank, 0, TAG_NB_NODES)); comm.send(old_nodes->values, nb_nodes, root, Tag::genTag(my_rank, 0, TAG_NODES)); /* --------<<<<-COORDINATES---------------------------------------------- */ nodes->resize(nb_nodes); AKANTU_DEBUG_INFO("Receiving coordinates from proc " << root); comm.receive(nodes->values, nb_nodes * spatial_dimension, root, Tag::genTag(root, 0, TAG_COORDINATES)); communicator.fillNodesType(mesh); /* --------<<<<-NODES_TYPE-2--------------------------------------------- */ Int * nodes_types = mesh.getNodesTypePointer()->values; AKANTU_DEBUG_INFO("Sending first nodes types to proc " << root); comm.send(nodes_types, nb_nodes, root, Tag::genTag(my_rank, 0, TAG_NODES_TYPE)); /* --------<<<<-NODES_TYPE-2--------------------------------------------- */ AKANTU_DEBUG_INFO("Receiving nodes types from proc " << root); comm.receive(nodes_types, nb_nodes, root, Tag::genTag(root, 0, TAG_NODES_TYPE)); } comm.broadcast(&(mesh.nb_global_nodes), 1, root); MeshUtils::fillElementToSubElementsData(mesh); AKANTU_DEBUG_OUT(); return &communicator; } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::fillTagBuffer(const MeshData & mesh_data, DynamicCommunicationBuffer * buffers, const std::string & tag_name, const ElementType & el_type, const Array & partition_num, const CSR & ghost_partition) { #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem) : { \ fillTagBufferTemplated(mesh_data, buffers, tag_name, el_type, partition_num, ghost_partition); \ break; \ } \ MeshDataTypeCode data_type_code = mesh_data.getTypeCode(tag_name); switch(data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default : AKANTU_DEBUG_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::fillNodesType(Mesh & mesh) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); Int * nodes_type = mesh.getNodesTypePointer()->values; UInt * nodes_set = new UInt[nb_nodes]; std::fill_n(nodes_set, nb_nodes, 0); const UInt NORMAL_SET = 1; const UInt GHOST_SET = 2; bool * already_seen = new bool[nb_nodes]; for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; UInt set = NORMAL_SET; if (gt == _ghost) set = GHOST_SET; std::fill_n(already_seen, nb_nodes, false); Mesh::type_iterator it = mesh.firstType(_all_dimensions, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, gt, _ek_not_defined); for(; it != end; ++it) { ElementType type = *it; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type, gt); Array::iterator< Vector > conn_it = mesh.getConnectivity(type, gt).begin(nb_nodes_per_element); for (UInt e = 0; e < nb_element; ++e, ++conn_it) { Vector & conn = *conn_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { AKANTU_DEBUG_ASSERT(conn(n) < nb_nodes, "Node " << conn(n) << " bigger than number of nodes " << nb_nodes); if(!already_seen[conn(n)]) { nodes_set[conn(n)] += set; already_seen[conn(n)] = true; } } } } } delete [] already_seen; for (UInt i = 0; i < nb_nodes; ++i) { if(nodes_set[i] == NORMAL_SET) nodes_type[i] = -1; else if(nodes_set[i] == GHOST_SET) nodes_type[i] = -3; else if(nodes_set[i] == (GHOST_SET + NORMAL_SET)) nodes_type[i] = -2; } delete [] nodes_set; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::fillCommunicationScheme(const UInt * partition, UInt nb_local_element, UInt nb_ghost_element, ElementType type) { AKANTU_DEBUG_IN(); Element element; element.type = type; element.kind = Mesh::getKind(type); const UInt * part = partition; part = partition; for (UInt lel = 0; lel < nb_local_element; ++lel) { UInt nb_send = *part; part++; element.element = lel; element.ghost_type = _not_ghost; for (UInt p = 0; p < nb_send; ++p) { UInt proc = *part; part++; AKANTU_DEBUG(dblAccessory, "Must send : " << element << " to proc " << proc); (send_element[proc]).push_back(element); } } for (UInt gel = 0; gel < nb_ghost_element; ++gel) { UInt proc = *part; part++; element.element = gel; element.ghost_type = _ghost; AKANTU_DEBUG(dblAccessory, "Must recv : " << element << " from proc " << proc); recv_element[proc].push_back(element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::asynchronousSynchronize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); if (communications.find(tag) == communications.end()) { communications[tag].resize(nb_proc); computeBufferSize(data_accessor, tag); } Communication & communication = communications[tag]; AKANTU_DEBUG_ASSERT(communication.send_requests.size() == 0, "There must be some pending sending communications. Tag is " << tag); std::map::iterator t_it = tag_counter.find(tag); UInt counter = 0; if(t_it == tag_counter.end()) { tag_counter[tag] = 0; } else { counter = ++(t_it->second); } for (UInt p = 0; p < nb_proc; ++p) { UInt ssize = communication.size_to_send[p]; if(p == rank || ssize == 0) continue; CommunicationBuffer & buffer = communication.send_buffer[p]; buffer.resize(ssize); #ifndef AKANTU_NDEBUG UInt nb_elements = send_element[p].getSize(); AKANTU_DEBUG_INFO("Packing data for proc " << p << " (" << ssize << "/" << nb_elements <<" data to send/elements)"); #endif data_accessor.packElementData(buffer, send_element[p], tag); AKANTU_DEBUG_ASSERT(buffer.getPackedSize() == ssize, "a problem have been introduced with " << "false sent sizes declaration " << buffer.getPackedSize() << " != " << ssize); AKANTU_DEBUG_INFO("Posting send to proc " << p << " (tag: " << tag << " - " << ssize << " data to send)" << " [" << Tag::genTag(rank, counter, tag) << "]"); communication.send_requests.push_back(static_communicator->asyncSend(buffer.storage(), ssize, p, Tag::genTag(rank, counter, tag))); } AKANTU_DEBUG_ASSERT(communication.recv_requests.size() == 0, "There must be some pending receive communications"); for (UInt p = 0; p < nb_proc; ++p) { UInt rsize = communication.size_to_receive[p]; if(p == rank || rsize == 0) continue; CommunicationBuffer & buffer = communication.recv_buffer[p]; buffer.resize(rsize); AKANTU_DEBUG_INFO("Posting receive from proc " << p << " (tag: " << tag << " - " << rsize << " data to receive) " << " [" << Tag::genTag(p, counter, tag) << "]"); communication.recv_requests.push_back(static_communicator->asyncReceive(buffer.storage(), rsize, p, Tag::genTag(p, counter, tag))); } #if defined(AKANTU_DEBUG_TOOLS) && defined(AKANTU_CORE_CXX11) static std::set tags; if(tags.find(tag) == tags.end()) { debug::element_manager.print(debug::_dm_synch, [&send_element, rank, nb_proc, tag, id](const Element & el)->std::string { std::stringstream out; UInt elp = 0; for (UInt p = 0; p < nb_proc; ++p) { UInt pos = send_element[p].find(el); if(pos != UInt(-1)) { if(elp > 0) out << std::endl; out << id << " send (" << pos << "/" << send_element[p].getSize() << ") to proc " << p << " tag:" << tag; ++elp; } } return out.str(); }); debug::element_manager.print(debug::_dm_synch, [&recv_element, rank, nb_proc, tag, id](const Element & el)->std::string { std::stringstream out; UInt elp = 0; for (UInt p = 0; p < nb_proc; ++p) { if(p == rank) continue; UInt pos = recv_element[p].find(el); if(pos != UInt(-1)) { if(elp > 0) out << std::endl; out << id << " recv (" << pos << "/" << recv_element[p].getSize() << ") from proc " << p << " tag:" << tag; ++elp; } } return out.str(); }); tags.insert(tag); } #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::waitEndSynchronize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(communications.find(tag) != communications.end(), "No communication with the tag \"" << tag <<"\" started"); Communication & communication = communications[tag]; std::vector req_not_finished; std::vector * req_not_finished_tmp = &req_not_finished; std::vector * recv_requests_tmp = &(communication.recv_requests); // static_communicator->waitAll(recv_requests); while(!recv_requests_tmp->empty()) { for (std::vector::iterator req_it = recv_requests_tmp->begin(); req_it != recv_requests_tmp->end() ; ++req_it) { CommunicationRequest * req = *req_it; if(static_communicator->testRequest(req)) { UInt proc = req->getSource(); AKANTU_DEBUG_INFO("Unpacking data coming from proc " << proc); CommunicationBuffer & buffer = communication.recv_buffer[proc]; data_accessor.unpackElementData(buffer, recv_element[proc], tag); buffer.resize(0); AKANTU_DEBUG_ASSERT(buffer.getLeftToUnpack() == 0, "all data have not been unpacked: " << buffer.getLeftToUnpack() << " bytes left"); static_communicator->freeCommunicationRequest(req); } else { req_not_finished_tmp->push_back(req); } } std::vector * swap = req_not_finished_tmp; req_not_finished_tmp = recv_requests_tmp; recv_requests_tmp = swap; req_not_finished_tmp->clear(); } AKANTU_DEBUG_INFO("Waiting that every send requests are received"); static_communicator->waitAll(communication.send_requests); for (std::vector::iterator req_it = communication.send_requests.begin(); req_it != communication.send_requests.end() ; ++req_it) { CommunicationRequest & req = *(*req_it); if(static_communicator->testRequest(&req)) { UInt proc = req.getDestination(); CommunicationBuffer & buffer = communication.send_buffer[proc]; buffer.resize(0); static_communicator->freeCommunicationRequest(&req); } } communication.send_requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::computeBufferSize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); for (UInt p = 0; p < nb_proc; ++p) { UInt ssend = 0; UInt sreceive = 0; if(p != rank) { if(send_element[p].getSize() != 0) { ssend = data_accessor.getNbDataForElements(send_element[p], tag); AKANTU_DEBUG_INFO("I have " << ssend << "(" << ssend / 1024. << "kB - "<< send_element[p].getSize() <<" element(s)) data to send to " << p << " for tag " << tag); } if(recv_element[p].getSize() != 0) { sreceive = data_accessor.getNbDataForElements(recv_element[p], tag); AKANTU_DEBUG_INFO("I have " << sreceive << "(" << sreceive / 1024. << "kB - "<< recv_element[p].getSize() <<" element(s)) data to receive for tag " << tag); } } communications[tag].size_to_send [p] = ssend; communications[tag].size_to_receive[p] = sreceive; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); Int prank = StaticCommunicator::getStaticCommunicator().whoAmI(); Int psize = StaticCommunicator::getStaticCommunicator().getNbProc(); stream << "[" << prank << "/" << psize << "]" << space << "DistributedSynchronizer [" << std::endl; for (UInt p = 0; p < nb_proc; ++p) { if (p == UInt(prank)) continue; stream << "[" << prank << "/" << psize << "]" << space << " + Communication to proc " << p << " [" << std::endl; if(AKANTU_DEBUG_TEST(dblDump)) { stream << "[" << prank << "/" << psize << "]" << space << " - Element to send to proc " << p << " [" << std::endl; Array::iterator it_el = send_element[p].begin(); Array::iterator end_el = send_element[p].end(); for(;it_el != end_el; ++it_el) stream << "[" << prank << "/" << psize << "]" << space << " " << *it_el << std::endl; stream << "[" << prank << "/" << psize << "]" << space << " ]" << std::endl; stream << "[" << prank << "/" << psize << "]" << space << " - Element to recv from proc " << p << " [" << std::endl; it_el = recv_element[p].begin(); end_el = recv_element[p].end(); for(;it_el != end_el; ++it_el) stream << "[" << prank << "/" << psize << "]" << space << " " << *it_el << std::endl;stream << "[" << prank << "/" << psize << "]" << space << " ]" << std::endl; } std::map< SynchronizationTag, Communication>::const_iterator it = communications.begin(); std::map< SynchronizationTag, Communication>::const_iterator end = communications.end(); for (; it != end; ++it) { const SynchronizationTag & tag = it->first; const Communication & communication = it->second; UInt ssend = communication.size_to_send[p]; UInt sreceive = communication.size_to_receive[p]; stream << "[" << prank << "/" << psize << "]" << space << " - Tag " << tag << " -> " << ssend << "byte(s) -- <- " << sreceive << "byte(s)" << std::endl; } } } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::onElementsRemoved(const Array & element_to_remove, const ByElementTypeUInt & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt psize = comm.getNbProc(); UInt prank = comm.whoAmI(); std::vector isend_requests; Array * list_of_el = new Array[nb_proc]; // Handling ghost elements for (UInt p = 0; p < psize; ++p) { if (p == prank) continue; Array & recv = recv_element[p]; if(recv.getSize() == 0) continue; Array::iterator recv_begin = recv.begin(); Array::iterator recv_end = recv.end(); Array::const_iterator er_it = element_to_remove.begin(); Array::const_iterator er_end = element_to_remove.end(); Array & list = list_of_el[p]; for (UInt i = 0; recv_begin != recv_end; ++i, ++recv_begin) { const Element & el = *recv_begin; Array::const_iterator pos = std::find(er_it, er_end, el); if(pos == er_end) { list.push_back(i); } } if(list.getSize() == recv.getSize()) list.push_back(UInt(0)); else list.push_back(UInt(-1)); AKANTU_DEBUG_INFO("Sending a message of size " << list.getSize() << " to proc " << p << " TAG(" << Tag::genTag(prank, 0, 0) << ")"); isend_requests.push_back(comm.asyncSend(list.storage(), list.getSize(), p, Tag::genTag(prank, 0, 0))); list.erase(list.getSize() - 1); if(list.getSize() == recv.getSize()) continue; Array new_recv; for (UInt nr = 0; nr < list.getSize(); ++nr) { Element & el = recv(list(nr)); el.element = new_numbering(el.type, el.ghost_type)(el.element); new_recv.push_back(el); } AKANTU_DEBUG_INFO("I had " << recv.getSize() << " elements to recv from proc " << p << " and " << list.getSize() << " elements to keep. I have " << new_recv.getSize() << " elements left."); recv.copy(new_recv); } for (UInt p = 0; p < psize; ++p) { if (p == prank) continue; Array & send = send_element[p]; if(send.getSize() == 0) continue; CommunicationStatus status; AKANTU_DEBUG_INFO("Getting number of elements of proc " << p << " not needed anymore TAG("<< Tag::genTag(p, 0, 0) <<")"); comm.probe(p, Tag::genTag(p, 0, 0), status); Array list(status.getSize()); AKANTU_DEBUG_INFO("Receiving list of elements (" << status.getSize() - 1 << " elements) no longer needed by proc " << p << " TAG("<< Tag::genTag(p, 0, 0) <<")"); comm.receive(list.storage(), list.getSize(), p, Tag::genTag(p, 0, 0)); if(list.getSize() == 1 && list(0) == 0) continue; list.erase(list.getSize() - 1); Array new_send; for (UInt ns = 0; ns < list.getSize(); ++ns) { new_send.push_back(send(list(ns))); } AKANTU_DEBUG_INFO("I had " << send.getSize() << " elements to send to proc " << p << " and " << list.getSize() << " elements to keep. I have " << new_send.getSize() << " elements left."); send.copy(new_send); } comm.waitAll(isend_requests); comm.freeCommunicationRequest(isend_requests); delete [] list_of_el; AKANTU_DEBUG_OUT(); } // void DistributedSynchronizer::checkCommunicationScheme() { // for (UInt p = 0; p < psize; ++p) { // if (p == prank) continue; // for(UInt e(0), e < recv_element.getSize()) // } // } /* -------------------------------------------------------------------------- */ void DistributedSynchronizer::buildPrankToElement(ByElementTypeUInt & prank_to_element) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt psize = comm.getNbProc(); UInt prank = comm.whoAmI(); UInt spatial_dimension = mesh.getSpatialDimension(); mesh.initByElementTypeArray(prank_to_element, 1, spatial_dimension, false, _ek_not_defined, true); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined); /// assign prank to all not ghost elements for (; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); Array & prank_to_el = prank_to_element(*it); for (UInt el = 0; el < nb_element; ++el) { prank_to_el(el) = prank; } } /// assign prank to all ghost elements for (UInt p = 0; p < psize; ++p) { UInt nb_ghost_element = recv_element[p].getSize(); for (UInt el = 0; el < nb_ghost_element; ++el) { UInt element = recv_element[p](el).element; ElementType type = recv_element[p](el).type; GhostType ghost_type = recv_element[p](el).ghost_type; Array & prank_to_el = prank_to_element(type, ghost_type); prank_to_el(element) = p; } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/test/test_synchronizer/test_grid_synchronizer.cc b/test/test_synchronizer/test_grid_synchronizer.cc index 8d31bc95d..9904e71f7 100644 --- a/test/test_synchronizer/test_grid_synchronizer.cc +++ b/test/test_synchronizer/test_grid_synchronizer.cc @@ -1,296 +1,296 @@ /** * @file test_grid_synchronizer.cc * * @author Nicolas Richart * * @date Fri Nov 25 17:00:17 2011 * * @brief test the GridSynchronizer object * * @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 "aka_common.hh" #include "aka_grid_dynamic.hh" #include "mesh.hh" #include "grid_synchronizer.hh" #include "mesh_partition.hh" #include "synchronizer_registry.hh" #include "test_data_accessor.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" #endif //AKANTU_USE_IOHELPER using namespace akantu; const UInt spatial_dimension = 3; typedef std::map, Real> pair_list; #include "test_grid_tools.hh" static void updatePairList(const ByElementTypeReal & barycenter, const SpatialGrid & grid, Real radius, pair_list & neighbors, neighbors_map_t::type & neighbors_map) { AKANTU_DEBUG_IN(); GhostType ghost_type = _not_ghost; Element e; e.ghost_type = ghost_type; // generate the pair of neighbor depending of the cell_list ByElementTypeReal::type_iterator it = barycenter.firstType(_all_dimensions, ghost_type); ByElementTypeReal::type_iterator last_type = barycenter.lastType(0, ghost_type); for(; it != last_type; ++it) { // loop over quad points e.type = *it; e.element = 0; const Array & barycenter_vect = barycenter(*it, ghost_type); UInt sp = barycenter_vect.getNbComponent(); Array::const_iterator< Vector > bary = barycenter_vect.begin(sp); Array::const_iterator< Vector > bary_end = barycenter_vect.end(sp); for(;bary != bary_end; ++bary, e.element++) { Point pt1(*bary); SpatialGrid::CellID cell_id = grid.getCellID(*bary); SpatialGrid::neighbor_cells_iterator first_neigh_cell = grid.beginNeighborCells(cell_id); SpatialGrid::neighbor_cells_iterator last_neigh_cell = grid.endNeighborCells(cell_id); // loop over neighbors cells of the one containing the current element for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { SpatialGrid::Cell::const_iterator first_neigh_el = grid.beginCell(*first_neigh_cell); SpatialGrid::Cell::const_iterator last_neigh_el = grid.endCell(*first_neigh_cell); // loop over the quadrature point in the current cell of the cell list for (;first_neigh_el != last_neigh_el; ++first_neigh_el){ const Element & elem = *first_neigh_el; Array::const_iterator< Vector > neigh_it = barycenter(elem.type, elem.ghost_type).begin(sp); const Vector & neigh_bary = neigh_it[elem.element]; Real distance = bary->distance(neigh_bary); if(distance <= radius) { Point pt2(neigh_bary); neighbors_map[pt1].push_back(pt2); std::pair pair = std::make_pair(e, elem); pair_list::iterator p = neighbors.find(pair); if(p != neighbors.end()) { AKANTU_DEBUG_ERROR("Pair already registered [" << e << " " << elem << "] -> " << p->second << " " << distance); } else { neighbors[pair] = distance; } } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Main */ /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(argc, argv); Real radius = 0.2; Mesh mesh(spatial_dimension); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); DistributedSynchronizer * dist = NULL; if(prank == 0) { mesh.read("bar3d.msh"); MeshPartition * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLowerBounds(lower_bounds); mesh.getUpperBounds(upper_bounds); Vector spacing(spatial_dimension); Vector center(spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { spacing[i] = radius * 1.2; center[i] = (upper_bounds[i] + lower_bounds[i]) / 2.; } SpatialGrid grid(spatial_dimension, spacing, center); GhostType ghost_type = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); - ByElementTypeReal barycenters("", "", 0); + ByElementTypeReal barycenters("", ""); mesh.initByElementTypeArray(barycenters, spatial_dimension, spatial_dimension); Element e; e.ghost_type = ghost_type; for(; it != last_type; ++it) { UInt nb_element = mesh.getNbElement(*it, ghost_type); e.type = *it; Array & barycenter = barycenters(*it, ghost_type); barycenter.resize(nb_element); Array::iterator< Vector > bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem) { mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type); e.element = elem; grid.insert(e, *bary_it); ++bary_it; } } std::stringstream sstr; sstr << "mesh_" << prank << ".msh"; mesh.write(sstr.str()); Mesh grid_mesh(spatial_dimension, "grid_mesh", 0); std::stringstream sstr_grid; sstr_grid << "grid_mesh_" << prank << ".msh"; grid.saveAsMesh(grid_mesh); grid_mesh.write(sstr_grid.str()); std::cout << "Pouet 1" << std::endl; GridSynchronizer * grid_communicator = GridSynchronizer::createGridSynchronizer(mesh, grid); std::cout << "Pouet 2" << std::endl; ghost_type = _ghost; it = mesh.firstType(spatial_dimension, ghost_type); last_type = mesh.lastType(spatial_dimension, ghost_type); e.ghost_type = ghost_type; for(; it != last_type; ++it) { UInt nb_element = mesh.getNbElement(*it, ghost_type); e.type = *it; Array & barycenter = barycenters(*it, ghost_type); barycenter.resize(nb_element); Array::iterator< Vector > bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem) { mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type); e.element = elem; grid.insert(e, *bary_it); ++bary_it; } } Mesh grid_mesh_ghost(spatial_dimension, "grid_mesh_ghost", 0); std::stringstream sstr_gridg; sstr_gridg << "grid_mesh_ghost_" << prank << ".msh"; grid.saveAsMesh(grid_mesh_ghost); grid_mesh_ghost.write(sstr_gridg.str()); std::cout << "Pouet 3" << std::endl; neighbors_map_t::type neighbors_map; pair_list neighbors; updatePairList(barycenters, grid, radius, neighbors, neighbors_map); pair_list::iterator nit = neighbors.begin(); pair_list::iterator nend = neighbors.end(); std::stringstream sstrp; sstrp << "pairs_" << prank; std::ofstream fout(sstrp.str().c_str()); for(;nit != nend; ++nit) { fout << "[" << nit->first.first << "," << nit->first.second << "] -> " << nit->second << std::endl; } std::string file = "neighbors_ref"; std::stringstream sstrf; sstrf << file << "_" << prank; file = sstrf.str(); std::ofstream nout; nout.open(file); neighbors_map_t::type::iterator it_n = neighbors_map.begin(); neighbors_map_t::type::iterator end_n = neighbors_map.end(); for(;it_n != end_n; ++it_n) { std::sort(it_n->second.begin(), it_n->second.end()); std::vector< Point >::iterator it_v = it_n->second.begin(); std::vector< Point >::iterator end_v = it_n->second.end(); nout << "####" << std::endl; nout << it_n->second.size() << std::endl; nout << it_n->first << std::endl; nout << "#" << std::endl; for(;it_v != end_v; ++it_v) { nout << *it_v << std::endl; } } fout.close(); AKANTU_DEBUG_INFO("Creating TestAccessor"); TestAccessor test_accessor(mesh, barycenters); SynchronizerRegistry synch_registry(test_accessor); synch_registry.registerSynchronizer(*dist, _gst_smm_mass); synch_registry.registerSynchronizer(*grid_communicator, _gst_test); AKANTU_DEBUG_INFO("Synchronizing tag on Dist"); synch_registry.synchronize(_gst_smm_mass); AKANTU_DEBUG_INFO("Synchronizing tag on Grid"); synch_registry.synchronize(_gst_test); delete grid_communicator; delete dist; akantu::finalize(); return EXIT_SUCCESS; }