diff --git a/src/common/aka_math.cc b/src/common/aka_math.cc index 8eaab1f64..bb9b74e7f 100644 --- a/src/common/aka_math.cc +++ b/src/common/aka_math.cc @@ -1,250 +1,252 @@ /** * @file aka_math.cc * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author David Simon Kammer * @author Nicolas Richart * @author Leonardo Snozzi * @author Peter Spijker * @author Marco Vocialta * * @date creation: Wed Aug 04 2010 * @date last modification: Sun Aug 13 2017 * * @brief Implementation of the math toolbox * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_array.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void Math::matrix_vector(UInt m, UInt n, const Array & A, const Array & x, Array & y, Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.size() == x.size(), "The vector A(" << A.getID() << ") and the vector x(" << x.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * n, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT( x.getNbComponent() == n, "The vector x(" << x.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT( y.getNbComponent() == n, "The vector y(" << y.getID() << ") do not the good number of component."); UInt nb_element = A.size(); UInt offset_A = A.getNbComponent(); UInt offset_x = x.getNbComponent(); y.resize(nb_element); Real * A_val = A.storage(); Real * x_val = x.storage(); Real * y_val = y.storage(); for (UInt el = 0; el < nb_element; ++el) { matrix_vector(m, n, A_val, x_val, y_val, alpha); A_val += offset_A; x_val += offset_x; y_val += offset_x; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Math::matrix_matrix(UInt m, UInt n, UInt k, const Array & A, const Array & B, Array & C, Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.size() == B.size(), "The vector A(" << A.getID() << ") and the vector B(" << B.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * k, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT( B.getNbComponent() == k * n, "The vector B(" << B.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT( C.getNbComponent() == m * n, "The vector C(" << C.getID() << ") do not the good number of component."); UInt nb_element = A.size(); UInt offset_A = A.getNbComponent(); UInt offset_B = B.getNbComponent(); UInt offset_C = C.getNbComponent(); C.resize(nb_element); Real * A_val = A.storage(); Real * B_val = B.storage(); Real * C_val = C.storage(); for (UInt el = 0; el < nb_element; ++el) { matrix_matrix(m, n, k, A_val, B_val, C_val, alpha); A_val += offset_A; B_val += offset_B; C_val += offset_C; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Array & A, const Array & B, Array & C, Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.size() == B.size(), "The vector A(" << A.getID() << ") and the vector B(" << B.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * k, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT( B.getNbComponent() == k * n, "The vector B(" << B.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT( C.getNbComponent() == m * n, "The vector C(" << C.getID() << ") do not the good number of component."); UInt nb_element = A.size(); UInt offset_A = A.getNbComponent(); UInt offset_B = B.getNbComponent(); UInt offset_C = C.getNbComponent(); C.resize(nb_element); Real * A_val = A.storage(); Real * B_val = B.storage(); Real * C_val = C.storage(); for (UInt el = 0; el < nb_element; ++el) { matrix_matrixt(m, n, k, A_val, B_val, C_val, alpha); A_val += offset_A; B_val += offset_B; C_val += offset_C; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Math::compute_tangents(const Array & normals, Array & tangents) { AKANTU_DEBUG_IN(); UInt spatial_dimension = normals.getNbComponent(); UInt tangent_components = spatial_dimension * (spatial_dimension - 1); + if (tangent_components == 0) return; + AKANTU_DEBUG_ASSERT( tangent_components == tangents.getNbComponent(), "Cannot compute the tangents, the storage array for tangents" << " does not have the good amount of components."); UInt nb_normals = normals.size(); tangents.resize(nb_normals, 0.); Real * normal_it = normals.storage(); Real * tangent_it = tangents.storage(); /// compute first tangent for (UInt q = 0; q < nb_normals; ++q) { /// if normal is orthogonal to xy plane, arbitrarly define tangent if (Math::are_float_equal(Math::norm2(normal_it), 0)) tangent_it[0] = 1; else Math::normal2(normal_it, tangent_it); normal_it += spatial_dimension; tangent_it += tangent_components; } /// compute second tangent (3D case) if (spatial_dimension == 3) { normal_it = normals.storage(); tangent_it = tangents.storage(); for (UInt q = 0; q < nb_normals; ++q) { Math::normal3(normal_it, tangent_it, tangent_it + spatial_dimension); normal_it += spatial_dimension; tangent_it += tangent_components; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Math::reduce(Array & array) { UInt nb_values = array.size(); if (nb_values == 0) return 0.; UInt nb_values_to_sum = nb_values >> 1; std::sort(array.begin(), array.end()); // as long as the half is not empty while (nb_values_to_sum) { UInt remaining = (nb_values - 2 * nb_values_to_sum); if (remaining) array(nb_values - 2) += array(nb_values - 1); // sum to consecutive values and store the sum in the first half for (UInt i = 0; i < nb_values_to_sum; ++i) { array(i) = array(2 * i) + array(2 * i + 1); } nb_values = nb_values_to_sum; nb_values_to_sum >>= 1; } return array(0); } } // akantu diff --git a/src/mesh/mesh_inline_impl.cc b/src/mesh/mesh_inline_impl.cc index 684c76019..af4c9e1ef 100644 --- a/src/mesh/mesh_inline_impl.cc +++ b/src/mesh/mesh_inline_impl.cc @@ -1,767 +1,767 @@ /** * @file mesh_inline_impl.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Thu Jul 15 2010 * @date last modification: Mon Dec 18 2017 * * @brief Implementation of the inline functions of the mesh class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "element_class.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_INLINE_IMPL_CC__ #define __AKANTU_MESH_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline ElementKind Element::kind() const { return Mesh::getKind(type); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template Mesh::ElementTypesIteratorHelper Mesh::elementTypes(pack &&... _pack) const { return connectivities.elementTypes(_pack...); } /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh) : new_numbering(mesh.getNbNodes(), 1, "new_numbering") {} /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh, const ID & new_numbering_id) : new_numbering(new_numbering_id, mesh.getID(), mesh.getMemoryID()) {} /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewElementsEvent & event) { this->nodes_to_elements.resize(nodes->size()); for (const auto & elem : event.getList()) { const Array & conn = connectivities(elem.type, elem.ghost_type); UInt nb_nodes_per_elem = this->getNbNodesPerElement(elem.type); for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = conn(elem.element, n); if (not nodes_to_elements[node]) nodes_to_elements[node] = std::make_unique>(); nodes_to_elements[node]->insert(elem); } } EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(NewNodesEvent & event) { this->computeBoundingBox(); this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedElementsEvent & event) { this->connectivities.onElementsRemoved(event.getNewNumbering()); this->fillNodesToElements(); this->computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedNodesEvent & event) { const auto & new_numbering = event.getNewNumbering(); this->removeNodesFromArray(*nodes, new_numbering); if (nodes_global_ids and not is_mesh_facets) this->removeNodesFromArray(*nodes_global_ids, new_numbering); if (not is_mesh_facets) this->removeNodesFromArray(*nodes_flags, new_numbering); if (not nodes_to_elements.empty()) { std::vector>> tmp( nodes_to_elements.size()); auto it = nodes_to_elements.begin(); UInt new_nb_nodes = 0; for (auto new_i : new_numbering) { if (new_i != UInt(-1)) { tmp[new_i] = std::move(*it); ++new_nb_nodes; } ++it; } tmp.resize(new_nb_nodes); std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin()); } computeBoundingBox(); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template inline void Mesh::removeNodesFromArray(Array & vect, const Array & new_numbering) { Array tmp(vect.size(), vect.getNbComponent()); UInt nb_component = vect.getNbComponent(); UInt new_nb_nodes = 0; for (UInt i = 0; i < new_numbering.size(); ++i) { UInt new_i = new_numbering(i); if (new_i != UInt(-1)) { T * to_copy = vect.storage() + i * nb_component; std::uninitialized_copy(to_copy, to_copy + nb_component, tmp.storage() + new_i * nb_component); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if (not nodes_global_ids) { nodes_global_ids = std::make_shared>( nodes->size(), 1, getID() + ":nodes_global_ids"); for (auto && global_ids : enumerate(*nodes_global_ids)) { std::get<1>(global_ids) = std::get<0>(global_ids); } } AKANTU_DEBUG_OUT(); return *nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getConnectivityPointer(const ElementType & type, const GhostType & ghost_type) { if (connectivities.exists(type, ghost_type)) return connectivities(type, ghost_type); if (ghost_type != _not_ghost) { ghosts_counters.alloc(0, 1, type, ghost_type, 1); } AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); return connectivities.alloc(0, nb_nodes_per_element, type, ghost_type); } /* -------------------------------------------------------------------------- */ inline Array> & Mesh::getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type) { return getDataPointer>("element_to_subelement", type, ghost_type, 1, true); } /* -------------------------------------------------------------------------- */ inline Array & Mesh::getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type) { auto & array = getDataPointer( "subelement_to_element", type, ghost_type, getNbFacetsPerElement(type), true, is_mesh_facets, ElementNull); return array; } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement() const { return getData>("element_to_subelement"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) const { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) { return getData>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getElementToSubelement(const Element & element) const { return getData>("element_to_subelement")(element); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getElementToSubelement(const Element & element) { return getData>("element_to_subelement")(element); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement() const { return getData("subelement_to_element"); } /* -------------------------------------------------------------------------- */ inline const auto & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) const { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline auto & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) { return getData("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElement(const Element & element) const { const auto & sub_to_element = this->getSubelementToElement(element.type, element.ghost_type); auto it = sub_to_element.begin(sub_to_element.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getSubelementToElement(const Element & element) { auto & sub_to_element = this->getSubelementToElement(element.type, element.ghost_type); auto it = sub_to_element.begin(sub_to_element.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent) { Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) tmp.resize(mesh_parent->getNbElement(el_type, ghost_type)); else tmp.resize(this->getNbElement(el_type, ghost_type)); } else { tmp.resize(0); } return tmp; } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_) { Array & tmp = this->getElementalDataArrayAlloc( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) tmp.resize(mesh_parent->getNbElement(el_type, ghost_type), defaul_); else tmp.resize(this->getNbElement(el_type, ghost_type), defaul_); } else { tmp.resize(0); } return tmp; } /* -------------------------------------------------------------------------- */ template inline const Array & Mesh::getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type) const { return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline Array & Mesh::getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type) { return this->getElementalDataArray(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template inline const ElementTypeMapArray & Mesh::getData(const ID & data_name) const { return this->getElementalData(data_name); } /* -------------------------------------------------------------------------- */ template inline ElementTypeMapArray & Mesh::getData(const ID & data_name) { return this->getElementalData(data_name); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const ElementType & type, const GhostType & ghost_type) const { try { const Array & conn = connectivities(type, ghost_type); return conn.size(); } catch (...) { return 0; } } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & kind) const { AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == UInt(-1), "spatial_dimension is " << spatial_dimension << " and is greater than 3 !"); UInt nb_element = 0; for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) nb_element += getNbElement(type, ghost_type); return nb_element; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(const Element & element, Vector & barycenter) const { Vector conn = getConnectivity(element); Matrix local_coord(spatial_dimension, conn.size()); auto node_begin = make_view(*nodes, spatial_dimension).begin(); for (auto && node : enumerate(conn)) { local_coord(std::get<0>(node)) = Vector(node_begin[std::get<1>(node)]); } Math::barycenter(local_coord.storage(), conn.size(), spatial_dimension, barycenter.storage()); } /* -------------------------------------------------------------------------- */ 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 UInt Mesh::getNbFacetTypes(const ElementType & type, __attribute__((unused)) UInt t) { UInt nb = 0; #define GET_NB_FACET_TYPE(type) nb = ElementClass::getNbFacetTypes() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET_TYPE); #undef GET_NB_FACET_TYPE return nb; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getFacetType(const ElementType & type, UInt t) { #define GET_FACET_TYPE(type) return ElementClass::getFacetType(t); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return _not_defined; } /* -------------------------------------------------------------------------- */ inline constexpr auto Mesh::getAllFacetTypes(const ElementType & type) { #define GET_FACET_TYPE(type) return ElementClass::getFacetTypes(); AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE); #undef GET_FACET_TYPE return ElementClass<_not_defined>::getFacetTypes(); } /* -------------------------------------------------------------------------- */ 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 UInt Mesh::getNbFacetsPerElement(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetLocalConnectivity(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); #define GET_FACET_CON(type) \ AKANTU_DEBUG_OUT(); \ return ElementClass::getFacetLocalConnectivityPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON); #undef GET_FACET_CON AKANTU_DEBUG_OUT(); return ElementClass<_not_defined>::getFacetLocalConnectivityPerElement(0); // This avoid a compilation warning but will certainly // also cause a segfault if reached } /* -------------------------------------------------------------------------- */ inline auto Mesh::getFacetConnectivity(const Element & element, UInt t) const { AKANTU_DEBUG_IN(); Matrix local_facets(getFacetLocalConnectivity(element.type, t)); Matrix facets(local_facets.rows(), local_facets.cols()); const Array & conn = connectivities(element.type, element.ghost_type); for (UInt f = 0; f < facets.rows(); ++f) { for (UInt n = 0; n < facets.cols(); ++n) { facets(f, n) = conn(element.element, local_facets(f, n)); } } AKANTU_DEBUG_OUT(); return facets; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivity(const Element & element) const { const auto & conn = connectivities(element.type, element.ghost_type); auto it = conn.begin(conn.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ inline VectorProxy Mesh::getConnectivity(const Element & element) { auto & conn = connectivities(element.type, element.ghost_type); auto it = conn.begin(conn.getNbComponent()); return it[element.element]; } /* -------------------------------------------------------------------------- */ 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)); } } /* -------------------------------------------------------------------------- */ inline void Mesh::addConnectivityType(const ElementType & type, const GhostType & ghost_type) { getConnectivityPointer(type, ghost_type); } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPureGhostNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_pure_ghost; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_local_master_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_normal; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_master; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_slave; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPeriodicSlave(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) == NodeFlag::_periodic_slave; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPeriodicMaster(UInt n) const { return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) == NodeFlag::_periodic_master; } /* -------------------------------------------------------------------------- */ inline NodeFlag Mesh::getNodeFlag(UInt local_id) const { return (*nodes_flags)(local_id); } /* -------------------------------------------------------------------------- */ inline Int Mesh::getNodePrank(UInt local_id) const { auto it = nodes_prank.find(local_id); return it == nodes_prank.end() ? -1 : it->second; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeLocalId(UInt global_id) const { - AKANTU_DEBUG_ASSERT(nodes_global_ids != nullptr, - "The global ids are note set."); + if (nodes_global_ids == nullptr) + return global_id; return nodes_global_ids->find(global_id); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->size(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElementList(const Array & elements) { UInt nb_nodes_per_element = 0; UInt nb_nodes = 0; ElementType current_element_type = _not_defined; for (const auto & el : elements) { if (el.type != current_element_type) { current_element_type = el.type; nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type); } nb_nodes += nb_nodes_per_element; } return nb_nodes; } /* -------------------------------------------------------------------------- */ inline Mesh & Mesh::getMeshFacets() { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshFacets() const { if (!this->mesh_facets) AKANTU_SILENT_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshParent() const { if (!this->mesh_parent) AKANTU_SILENT_EXCEPTION( "No parent mesh is defined! This is only valid in a mesh_facets"); return *this->mesh_parent; } /* -------------------------------------------------------------------------- */ void Mesh::addPeriodicSlave(UInt slave, UInt master) { if (master == slave) return; // if pair already registered auto master_slaves = periodic_master_slave.equal_range(master); auto slave_it = std::find_if(master_slaves.first, master_slaves.second, [&](auto & pair) { return pair.second == slave; }); if (slave_it == master_slaves.second) { // no duplicates periodic_master_slave.insert(std::make_pair(master, slave)); AKANTU_DEBUG_INFO("adding periodic slave, slave gid:" << getNodeGlobalId(slave) << " [lid: " << slave << "]" << ", master gid:" << getNodeGlobalId(master) << " [lid: " << master << "]"); // std::cout << "adding periodic slave, slave gid:" << getNodeGlobalId(slave) // << " [lid: " << slave << "]" // << ", master gid:" << getNodeGlobalId(master) // << " [lid: " << master << "]" << std::endl; } periodic_slave_master[slave] = master; auto set_flag = [&](auto node, auto flag) { (*nodes_flags)[node] &= ~NodeFlag::_periodic_mask; // clean periodic flags (*nodes_flags)[node] |= flag; }; set_flag(slave, NodeFlag::_periodic_slave); set_flag(master, NodeFlag::_periodic_master); } /* -------------------------------------------------------------------------- */ UInt Mesh::getPeriodicMaster(UInt slave) const { return periodic_slave_master.at(slave); } /* -------------------------------------------------------------------------- */ class Mesh::PeriodicSlaves { using internal_iterator = std::unordered_multimap::const_iterator; std::pair pair; public: PeriodicSlaves(const Mesh & mesh, UInt master) : pair(mesh.getPeriodicMasterSlaves().equal_range(master)) {} PeriodicSlaves(const PeriodicSlaves & other) = default; PeriodicSlaves(PeriodicSlaves && other) = default; PeriodicSlaves & operator=(const PeriodicSlaves & other) = default; class const_iterator { internal_iterator it; public: const_iterator(internal_iterator it) : it(std::move(it)) {} const_iterator operator++() { ++it; return *this; } bool operator!=(const const_iterator & other) { return other.it != it; } auto operator*() { return it->second; } }; auto begin() { return const_iterator(pair.first); } auto end() { return const_iterator(pair.second); } }; /* -------------------------------------------------------------------------- */ inline decltype(auto) Mesh::getPeriodicSlaves(UInt master) const { return PeriodicSlaves(*this, master); } /* -------------------------------------------------------------------------- */ inline Vector Mesh::getConnectivityWithPeriodicity(const Element & element) const { Vector conn = getConnectivity(element); if (not isPeriodic()) { return conn; } for (auto && node : conn) { if (isPeriodicSlave(node)) { node = getPeriodicMaster(node); } } return conn; } } // namespace akantu #endif /* __AKANTU_MESH_INLINE_IMPL_CC__ */ diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic.cc b/test/test_model/test_model_solver/test_model_solver_dynamic.cc index a0e7f651b..0b0e826cb 100644 --- a/test/test_model/test_model_solver/test_model_solver_dynamic.cc +++ b/test/test_model/test_model_solver/test_model_solver_dynamic.cc @@ -1,226 +1,226 @@ /** * @file test_model_solver_dynamic.cc * * @author Nicolas Richart * * @date creation: Wed Apr 13 2016 * @date last modification: Tue Feb 20 2018 * * @brief Test default dof manager * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "data_accessor.hh" #include "dof_manager.hh" #include "dof_manager_default.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "model_solver.hh" #include "non_linear_solver.hh" #include "sparse_matrix.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #include "dumper_element_partition.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "test_model_solver_my_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef EXPLICIT #define EXPLICIT true #endif using namespace akantu; static void genMesh(Mesh & mesh, UInt nb_nodes); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); UInt prank = Communicator::getStaticCommunicator().whoAmI(); UInt global_nb_nodes = 201; UInt max_steps = 400; Real time_step = 0.001; Mesh mesh(1); Real F = -9.81; bool _explicit = EXPLICIT; const Real pulse_width = 0.2; const Real A = 0.01; if (prank == 0) genMesh(mesh, global_nb_nodes); mesh.distribute(); - mesh.makePeriodic(_x); + //mesh.makePeriodic(_x); MyModel model(F, mesh, _explicit); model.forces.clear(); model.blocked.clear(); for (auto && n : arange(mesh.getNbNodes())) { Real x = mesh.getNodes()(n) - 0.2; // Sinus * Gaussian Real L = pulse_width; Real k = 2 * M_PI / L; //model.displacement(n) = A * sin(k * x) * exp(-(k * x) * (k * x) / (L * L)); model.displacement(n) = A * cos(k * x); model.velocity(n) = k * A * sin(k * x); } if (!_explicit) { model.getNewSolver("dynamic", _tsst_dynamic, _nls_newton_raphson); model.setIntegrationScheme("dynamic", "disp", _ist_trapezoidal_rule_2, IntegrationScheme::_displacement); } else { model.getNewSolver("dynamic", _tsst_dynamic_lumped, _nls_lumped); model.setIntegrationScheme("dynamic", "disp", _ist_central_difference, IntegrationScheme::_acceleration); } model.setTimeStep(time_step); // #if EXPLICIT == true // std::ofstream output("output_dynamic_explicit.csv"); // #else // std::ofstream output("output_dynamic_implicit.csv"); // #endif if (prank == 0) { std::cout << std::scientific; std::cout << std::setw(14) << "time" << "," << std::setw(14) << "wext" << "," << std::setw(14) << "epot" << "," << std::setw(14) << "ekin" << "," << std::setw(14) << "total" << "," << std::setw(14) << "max_disp" << "," << std::setw(14) << "min_disp" << std::endl; } Real wext = 0.; model.getDOFManager().clearResidual(); model.assembleResidual(); Real epot = 0; // model.getPotentialEnergy(); Real ekin = 0; // model.getKineticEnergy(); Real einit = ekin + epot; Real etot = ekin + epot - wext - einit; Real max_disp = 0., min_disp = 0.; for (auto && disp : model.displacement) { max_disp = std::max(max_disp, disp); min_disp = std::min(min_disp, disp); } if (prank == 0) { std::cout << std::setw(14) << 0. << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << "," << std::setw(14) << max_disp << "," << std::setw(14) << min_disp << std::endl; } #if EXPLICIT == false NonLinearSolver & solver = model.getDOFManager().getNonLinearSolver("dynamic"); solver.set("max_iterations", 20); #endif auto * dumper = new DumperParaview("dynamic", "./paraview"); mesh.registerExternalDumper(*dumper, "dynamic", true); mesh.addDumpMesh(mesh); mesh.addDumpFieldExternalToDumper("dynamic", "displacement", model.displacement); mesh.addDumpFieldExternalToDumper("dynamic", "velocity", model.velocity); mesh.addDumpFieldExternalToDumper("dynamic", "forces", model.forces); mesh.addDumpFieldExternalToDumper("dynamic", "acceleration", model.acceleration); mesh.dump(); for (UInt i = 1; i < max_steps + 1; ++i) { model.solveStep("dynamic"); mesh.dump(); #if EXPLICIT == false // int nb_iter = solver.get("nb_iterations"); // Real error = solver.get("error"); // bool converged = solver.get("converged"); // if (prank == 0) // std::cerr << error << " " << nb_iter << " -> " << converged << std::endl; #endif epot = model.getPotentialEnergy(); ekin = model.getKineticEnergy(); wext += model.getExternalWorkIncrement(); etot = ekin + epot - wext - einit; Real max_disp = 0., min_disp = 0.; for (auto && disp : model.displacement) { max_disp = std::max(max_disp, disp); min_disp = std::min(min_disp, disp); } if (prank == 0) { std::cout << std::setw(14) << time_step * i << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << "," << std::setw(14) << max_disp << "," << std::setw(14) << min_disp << std::endl; } // if (std::abs(etot) > 1e-1) { // AKANTU_ERROR("The total energy of the system is not conserved!"); //} } // output.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); } conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } mesh_accessor.makeReady(); }