diff --git a/src/common/aka_csr.hh b/src/common/aka_csr.hh index 7227cceee..f9f0ca805 100644 --- a/src/common/aka_csr.hh +++ b/src/common/aka_csr.hh @@ -1,258 +1,289 @@ /** * @file aka_csr.hh * * @author Nicolas Richart * * @date creation: Wed Apr 20 2011 * @date last modification: Sun Dec 03 2017 * * @brief A compresed sparse row structure based on akantu Arrays * * @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_array.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_CSR_HH__ #define __AKANTU_AKA_CSR_HH__ namespace akantu { /** * This class can be used to store the structure of a sparse matrix or for * vectors with variable number of component per element * * @param nb_rows number of rows of a matrix or size of a vector. */ template class CSR { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: explicit CSR(UInt nb_rows = 0) : nb_rows(nb_rows), rows_offsets(nb_rows + 1, 1, "rows_offsets"), rows(0, 1, "rows") { rows_offsets.clear(); }; virtual ~CSR() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// does nothing inline void beginInsertions(){}; /// insert a new entry val in row row inline UInt insertInRow(UInt row, const T & val) { UInt pos = rows_offsets(row)++; rows(pos) = val; return pos; } /// access an element of the matrix inline const T & operator()(UInt row, UInt col) const { AKANTU_DEBUG_ASSERT(rows_offsets(row + 1) - rows_offsets(row) > col, "This element is not present in this CSR"); return rows(rows_offsets(row) + col); } /// access an element of the matrix inline T & operator()(UInt row, UInt col) { AKANTU_DEBUG_ASSERT(rows_offsets(row + 1) - rows_offsets(row) > col, "This element is not present in this CSR"); return rows(rows_offsets(row) + col); } inline void endInsertions() { for (UInt i = nb_rows; i > 0; --i) rows_offsets(i) = rows_offsets(i - 1); rows_offsets(0) = 0; } inline void countToCSR() { for (UInt i = 1; i < nb_rows; ++i) rows_offsets(i) += rows_offsets(i - 1); for (UInt i = nb_rows; i >= 1; --i) rows_offsets(i) = rows_offsets(i - 1); rows_offsets(0) = 0; } inline void clearRows() { rows_offsets.clear(); rows.resize(0); }; inline void resizeRows(UInt nb_rows) { this->nb_rows = nb_rows; rows_offsets.resize(nb_rows + 1); rows_offsets.clear(); } inline void resizeCols() { rows.resize(rows_offsets(nb_rows)); } inline void copy(Array & offsets, Array & values) { offsets.copy(rows_offsets); values.copy(rows); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// returns the number of rows inline UInt getNbRows() const { return rows_offsets.size() - 1; }; /// returns the number of non-empty columns in a given row inline UInt getNbCols(UInt row) const { return rows_offsets(row + 1) - rows_offsets(row); }; /// returns the offset (start of columns) for a given row inline UInt & rowOffset(UInt row) { return rows_offsets(row); }; /// iterator on a row template class iterator_internal : public std::iterator { public: using _parent = std::iterator; using pointer = typename _parent::pointer; using reference = typename _parent::reference; explicit iterator_internal(pointer x = nullptr) : pos(x){}; iterator_internal(const iterator_internal & it) : pos(it.pos){}; iterator_internal & operator++() { ++pos; return *this; }; iterator_internal operator++(int) { iterator tmp(*this); operator++(); return tmp; }; iterator_internal & operator--() { --pos; return *this; }; iterator_internal operator--(int) { iterator_internal tmp(*this); operator--(); return tmp; }; bool operator==(const iterator_internal & rhs) { return pos == rhs.pos; }; bool operator!=(const iterator_internal & rhs) { return pos != rhs.pos; }; reference operator*() { return *pos; }; pointer operator->() const { return pos; }; private: pointer pos; }; +#ifndef SWIG + template class CSRRow { + public: + CSRRow(R * begin, R * end) : begin_(begin), end_(end) {} + + inline auto begin() const { return begin_; } + inline auto end() const { return end_; } + + private: + iterator_internal begin_, end_; + }; +#endif + using iterator = iterator_internal; using const_iterator = iterator_internal; inline iterator begin(UInt row) { return iterator(rows.storage() + rows_offsets(row)); }; inline iterator end(UInt row) { return iterator(rows.storage() + rows_offsets(row + 1)); }; inline const_iterator begin(UInt row) const { return const_iterator(rows.storage() + rows_offsets(row)); }; inline const_iterator end(UInt row) const { return const_iterator(rows.storage() + rows_offsets(row + 1)); }; +#ifndef SWIG +private: + template decltype(auto) make_row(R * begin, R * end) { + return CSRRow(begin, end); + } + +public: + inline decltype(auto) getRow(UInt row) { + return make_row(rows.storage() + rows_offsets(row), + rows.storage() + rows_offsets(row + 1)); + } + + inline decltype(auto) getRow(UInt row) const { + return make_row(rows.storage() + rows_offsets(row), + rows.storage() + rows_offsets(row + 1)); + } +#endif + inline iterator rbegin(UInt row) { return iterator(rows.storage() + rows_offsets(row + 1) - 1); }; inline iterator rend(UInt row) { return iterator(rows.storage() + rows_offsets(row) - 1); }; inline const Array & getRowsOffset() const { return rows_offsets; }; inline const Array & getRows() const { return rows; }; inline Array & getRows() { return rows; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: UInt nb_rows; /// array of size nb_rows containing the offset where the values are stored in Array rows_offsets; /// compressed row values, values of row[i] are stored between rows_offsets[i] /// and rows_offsets[i+1] Array rows; }; /* -------------------------------------------------------------------------- */ /* Data CSR */ /* -------------------------------------------------------------------------- */ /** * Inherits from CSR and can contain information such as matrix values * where the mother class would be a CSR structure for row and cols * * @return nb_rows */ template class DataCSR : public CSR { public: DataCSR(UInt nb_rows = 0) : CSR(nb_rows), data(0, 1){}; inline void resizeCols() { CSR::resizeCols(); data.resize(rows_offsets(nb_rows)); } inline const Array & getData() const { return data; }; private: Array data; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "aka_csr_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const CSR & _this) // { // _this.printself(stream); // return stream; // } } // akantu #endif /* __AKANTU_AKA_CSR_HH__ */ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index 6358292e8..19e4192a7 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,1845 +1,1844 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Leonardo Snozzi * @author Marco Vocialta * * @date creation: Fri Aug 20 2010 * @date last modification: Wed Feb 21 2018 * * @brief All mesh utils necessary for various tasks * * @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 "mesh_utils.hh" #include "element_synchronizer.hh" #include "fe_engine.hh" #include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); /// count number of occurrence of each node UInt nb_nodes = mesh.getNbNodes(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); // AKANTU_DEBUG_ASSERT( // mesh.firstType(spatial_dimension) != mesh.lastType(spatial_dimension), // "Some elements must be found in right dimension to compute facets!"); for (auto && ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type, ghost_type); auto conn_it = mesh.getConnectivity(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) for (UInt n = 0; n < conn_it->size(); ++n) ++node_to_elem.rowOffset((*conn_it)(n)); } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list Element e; node_to_elem.beginInsertions(); for (auto && ghost_type : ghost_types) { e.ghost_type = ghost_type; for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); auto conn_it = mesh.getConnectivity(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; for (UInt n = 0; n < conn_it->size(); ++n) node_to_elem.insertInRow((*conn_it)(n), e); } } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh, CSR & node_to_elem, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, ghost_type).size(); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) ++node_to_elem.rowOffset(conn_val[el_offset + n]); } /// convert the occurrence array in a csr one node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(spatial_dimension - 1, ghost_type)) { mesh.getConnectivity(type, ghost_type).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt to_dimension) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt from_dimension, UInt to_dimension) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT( mesh_facets.isMeshFacets(), "The mesh_facets should be initialized with initMeshFacets"); /// generate facets buildFacetsDimension(mesh, mesh_facets, false, from_dimension); /// sort facets and generate sub-facets for (UInt i = from_dimension - 1; i > to_dimension; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension) { AKANTU_DEBUG_IN(); // save the current parent of mesh_facets and set it temporarly to mesh since // mesh is the one containing the elements for which mesh_facets has the // sub-elements // example: if the function is called with mesh = mesh_facets const Mesh * mesh_facets_parent = nullptr; try { mesh_facets_parent = &mesh_facets.getMeshParent(); } catch (...) { } mesh_facets.defineMeshParent(mesh); MeshAccessor mesh_accessor(mesh_facets); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const auto mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; std::vector connected_elements; // init the SubelementToElement data to improve performance for (auto && ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(dimension, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type); auto facet_types = mesh.getAllFacetTypes(type); for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); mesh_accessor.getElementToSubelement(facet_type, ghost_type); mesh_accessor.getConnectivity(facet_type, ghost_type); } } } const ElementSynchronizer * synchronizer = nullptr; if (mesh.isDistributed()) { synchronizer = &(mesh.getElementSynchronizer()); } Element current_element; for (auto && ghost_type : ghost_types) { GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; for (auto && type : mesh.elementTypes(dimension, ghost_type)) { auto facet_types = mesh.getAllFacetTypes(type); current_element.type = type; for (auto && ft : arange(facet_types.size())) { auto facet_type = facet_types(ft); auto nb_element = mesh.getNbElement(type, ghost_type); auto element_to_subelement = &mesh_facets.getElementToSubelement(facet_type, ghost_type); auto connectivity_facets = &mesh_facets.getConnectivity(facet_type, ghost_type); auto nb_facet_per_element = mesh.getNbFacetsPerElement(type, ft); const auto & element_connectivity = mesh.getConnectivity(type, ghost_type); Matrix facet_local_connectivity( mesh.getFacetLocalConnectivity(type, ft)); auto nb_nodes_per_facet = connectivity_facets->getNbComponent(); Vector facet(nb_nodes_per_facet); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; for (UInt f = 0; f < nb_facet_per_element; ++f) { for (UInt n = 0; n < nb_nodes_per_facet; ++n) facet(n) = element_connectivity(el, facet_local_connectivity(f, n)); UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0)); counter.resize(first_node_nb_elements); counter.clear(); // loop over the other nodes to search intersecting elements, // which are the elements that share another node with the // starting element after first_node UInt local_el = 0; auto first_node_elements = node_to_elem.begin(facet(0)); auto first_node_elements_end = node_to_elem.end(facet(0)); for (; first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (UInt n = 1; n < nb_nodes_per_facet; ++n) { auto node_elements_begin = node_to_elem.begin(facet(n)); auto node_elements_end = node_to_elem.end(facet(n)); counter(local_el) += std::count(node_elements_begin, node_elements_end, *first_node_elements); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; connected_elements.clear(); for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) { Element real_el = node_to_elem(facet(0), el_f); if (not(counter(el_f) == nb_nodes_per_facet - 1)) continue; ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } if (minimum_el != current_element) continue; bool full_ghost_facet = false; UInt n = 0; while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n))) ++n; if (n == nb_nodes_per_facet) full_ghost_facet = true; if (full_ghost_facet) continue; if (boundary_only and nb_element_connected_to_facet != 1) continue; std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(current_element); if (nb_element_connected_to_facet == 1) { /// boundary facet elements.push_back(ElementNull); } else if (nb_element_connected_to_facet == 2) { /// internal facet elements.push_back(connected_elements[1]); /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2] = {_not_ghost, _not_ghost}; for (UInt el = 0; el < connected_elements.size(); ++el) gt[el] = connected_elements[el].ghost_type; if ((gt[0] == _not_ghost) xor (gt[1] == _not_ghost)) { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { prank[el] = synchronizer->getRank(connected_elements[el]); } // ugly trick from Marco detected :P bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) facet_ghost_type = _not_ghost; else facet_ghost_type = _ghost; connectivity_facets = &mesh_facets.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_facets.getElementToSubelement( facet_type, facet_ghost_type); } } else { /// facet of facet for (UInt i = 1; i < nb_element_connected_to_facet; ++i) { elements.push_back(connected_elements[i]); } } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->size() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.size(); ++elem) { Element loc_el = elements[elem]; if (loc_el.type == _not_defined) continue; Array & subelement_to_element = mesh_facets.getSubelementToElement(loc_el.type, loc_el.ghost_type); UInt nb_facet_per_loc_element = subelement_to_element.getNbComponent(); for (UInt f_in = 0; f_in < nb_facet_per_loc_element; ++f_in) { auto & el = subelement_to_element(loc_el.element, f_in); if (el.type != _not_defined) continue; el.type = facet_type; el.element = current_facet; el.ghost_type = facet_ghost_type; break; } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = &mesh_accessor.getConnectivity(facet_type, facet_ghost_type); element_to_subelement = &mesh_accessor.getElementToSubelement( facet_type, facet_ghost_type); } } } } } } // restore the parent of mesh_facet if (mesh_facets_parent) mesh_facets.defineMeshParent(*mesh_facets_parent); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, Array & local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes_numbers) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; for (UInt i = 0; i < old_nodes_numbers.size(); ++i) { renumbering_map[old_nodes_numbers(i)] = i; } /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element) * nb_nodes_per_element, renumbering_map); old_nodes_numbers.resize(renumbering_map.size()); for (auto & renumber_pair : renumbering_map) { old_nodes_numbers(renumber_pair.second) = renumber_pair.first; } renumbering_map.clear(); MeshAccessor mesh_accessor(mesh); /// copy the renumbered connectivity to the right place auto & local_conn = mesh_accessor.getConnectivity(type); local_conn.resize(nb_local_element); memcpy(local_conn.storage(), local_connectivities.storage(), nb_local_element * nb_nodes_per_element * sizeof(UInt)); auto & ghost_conn = mesh_accessor.getConnectivity(type, _ghost); ghost_conn.resize(nb_ghost_element); std::memcpy(ghost_conn.storage(), local_connectivities.storage() + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); auto & ghost_counter = mesh_accessor.getGhostsCounters(type, _ghost); ghost_counter.resize(nb_ghost_element, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity( Array & list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes.storage(); UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; auto it = renumbering_map.find(node); if (it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Array & nodes_removed = remove_nodes.getList(); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { auto ghost_type = (GhostType)gt; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); Array & connectivity = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity.size()); renumberNodesInConnectivity( connectivity, nb_element * nb_nodes_per_element, renumbering_map); } } Array & new_numbering = remove_nodes.getNewNumbering(); std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); auto it = renumbering_map.begin(); auto end = renumbering_map.end(); for (; it != end; ++it) { new_numbering(it->first) = it->second; } for (UInt i = 0; i < new_numbering.size(); ++i) { if (new_numbering(i) == UInt(-1)) nodes_removed.push_back(i); } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_COHESIVE_ELEMENT) /* -------------------------------------------------------------------------- */ UInt MeshUtils::insertCohesiveElements( Mesh & mesh, Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion, Array & doubled_nodes, Array & new_elements, bool only_double_facets) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt elements_to_insert = updateFacetToDouble(mesh_facets, facet_insertion); if (elements_to_insert > 0) { if (spatial_dimension == 1) { doublePointFacet(mesh, mesh_facets, doubled_nodes); } else { doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes, true); findSubfacetToDouble(mesh_facets); if (spatial_dimension == 2) { doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes); } else if (spatial_dimension == 3) { doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false); findSubfacetToDouble(mesh_facets); doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes); } } if (!only_double_facets) updateCohesiveData(mesh, mesh_facets, new_elements); } return elements_to_insert; } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doubleNodes(Mesh & mesh, const std::vector & old_nodes, Array & doubled_nodes) { AKANTU_DEBUG_IN(); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt old_nb_nodes = position.size(); UInt new_nb_nodes = old_nb_nodes + old_nodes.size(); UInt old_nb_doubled_nodes = doubled_nodes.size(); UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.size(); position.resize(new_nb_nodes); doubled_nodes.resize(new_nb_doubled_nodes); Array::iterator> position_begin = position.begin(spatial_dimension); for (UInt n = 0; n < old_nodes.size(); ++n) { UInt new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes[n]; doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node; /// update position std::copy(position_begin + old_nodes[n], position_begin + old_nodes[n] + 1, position_begin + new_node); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets, UInt facet_dimension, Array & doubled_nodes, bool facet_mode) { AKANTU_DEBUG_IN(); for (auto gt_facet : ghost_types) { for (auto && type_facet : mesh_facets.elementTypes(facet_dimension, gt_facet)) { auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); if (nb_facet_to_double == 0) continue; // this while fail if multiple facet types // \TODO handle multiple sub-facet types auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(type_facet); auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto nb_nodes_per_facet = conn_facet.getNbComponent(); auto old_nb_facet = conn_facet.size(); auto new_nb_facet = old_nb_facet + nb_facet_to_double; #ifndef AKANTU_NDEBUG // memory initialization are slow but help debug conn_facet.resize(new_nb_facet, UInt(-1)); #else conn_facet.resize(new_nb_facet); #endif auto conn_facet_begin = conn_facet.begin(nb_nodes_per_facet); auto & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); #ifndef AKANTU_NDEBUG subfacet_to_facet.resize(new_nb_facet, ElementNull); #else subfacet_to_facet.resize(new_nb_facet); #endif auto subfacet_to_facet_begin = subfacet_to_facet.begin(nb_subfacet_per_facet); Element new_facet{type_facet, old_nb_facet, gt_facet}; auto conn_facet_new_it = conn_facet_begin + new_facet.element; auto subfacet_to_facet_new_it = subfacet_to_facet_begin + new_facet.element; for (UInt facet = 0; facet < nb_facet_to_double; ++facet, ++new_facet.element, ++conn_facet_new_it, ++subfacet_to_facet_new_it) { UInt old_facet = facets_to_double(facet); /// adding a new facet by copying original one /// copy connectivity in new facet *conn_facet_new_it = conn_facet_begin[old_facet]; /// update subfacet_to_facet *subfacet_to_facet_new_it = subfacet_to_facet_begin[old_facet]; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; /// update facet_to_subfacet array mesh_facets.getElementToSubelement(subfacet).push_back(new_facet); } } /// update facet_to_subfacet and _segment_3 facets if any if (not facet_mode) { updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true); updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true); updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } else updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt MeshUtils::updateFacetToDouble( Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); UInt nb_facets_to_double = 0.; for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { const auto & f_insertion = facet_insertion(type_facet, gt_facet); auto & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); ElementType el_type = _not_defined; GhostType el_gt = _casper; UInt nb_facet_per_element = 0; Element old_facet_el{type_facet, 0, gt_facet}; Array * facet_to_element = nullptr; for (UInt f = 0; f < f_insertion.size(); ++f) { if (f_insertion(f) == false) continue; ++nb_facets_to_double; if (element_to_facet(f)[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || element_to_facet(f)[1].kind() == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } f_to_double.push_back(f); UInt new_facet = mesh_facets.getNbElement(type_facet, gt_facet) + f_to_double.size() - 1; old_facet_el.element = f; /// update facet_to_element vector Element & elem_to_update = element_to_facet(f)[1]; UInt el = elem_to_update.element; if (elem_to_update.ghost_type != el_gt || elem_to_update.type != el_type) { el_type = elem_to_update.type; el_gt = elem_to_update.ghost_type; facet_to_element = &mesh_facets.getSubelementToElement(el_type, el_gt); nb_facet_per_element = facet_to_element->getNbComponent(); } auto begin = facet_to_element->storage() + el * nb_facet_per_element; auto f_update = std::find(begin, begin + nb_facet_per_element, old_facet_el); AKANTU_DEBUG_ASSERT(f_update != begin + nb_facet_per_element, "Facet not found"); f_update->element = new_facet; /// update elements connected to facet std::vector first_facet_list = element_to_facet(f); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1]; element_to_facet(f)[1] = ElementNull; element_to_facet(new_facet)[1] = ElementNull; } } } AKANTU_DEBUG_OUT(); return nb_facets_to_double; } /* -------------------------------------------------------------------------- */ void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { - for (auto type : - mesh_facets.elementTypes(_all_dimensions, gt)) { + for (auto type : mesh_facets.elementTypes(_all_dimensions, gt)) { mesh_facets.getDataPointer("facet_to_double", type, gt, 1, false); mesh_facets.getDataPointer>( "facets_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "elements_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer>( "subfacets_to_subsubfacet_double", type, gt, 1, false); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::findSubfacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); if (spatial_dimension == 1) return; for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); GhostType gt_subfacet = _casper; ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet); GhostType gt_subsubfacet = _casper; Array * conn_subfacet = nullptr; Array * sf_to_double = nullptr; Array> * sf_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; Array> * el_to_subfacet_double = nullptr; UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet); UInt nb_subsubfacet; UInt nb_nodes_per_sf_el; if (subsubfacet_mode) { nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet); nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet); } else nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array * subsubfacet_to_subfacet = nullptr; UInt old_nb_facet = subfacet_to_facet.size() - nb_facet_to_double; Element current_facet{type_facet, 0, gt_facet}; std::vector element_list; std::vector facet_list; std::vector * subfacet_list; if (subsubfacet_mode) subfacet_list = new std::vector; /// map to filter subfacets Array> * facet_to_subfacet = nullptr; /// this is used to make sure that both new and old facets are /// checked UInt facets[2]; /// loop on every facet for (UInt f_index = 0; f_index < 2; ++f_index) { for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { facets[bool(f_index)] = facets_to_double(facet); facets[!bool(f_index)] = old_nb_facet + facet; UInt old_facet = facets[0]; UInt new_facet = facets[1]; Element & starting_element = element_to_facet(new_facet)[0]; current_facet.element = old_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; if (subsubfacet_mode) { subsubfacet_to_subfacet = &mesh_facets.getSubelementToElement( type_subfacet, gt_subfacet); } else { conn_subfacet = &mesh_facets.getConnectivity(type_subfacet, gt_subfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subfacet, gt_subfacet); } } if (subsubfacet_mode) { /// loop on every subsubfacet for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) { Element & subsubfacet = (*subsubfacet_to_subfacet)(subfacet.element, ssf); if (subsubfacet == ElementNull) continue; if (gt_subsubfacet != subsubfacet.ghost_type) { gt_subsubfacet = subsubfacet.ghost_type; conn_subfacet = &mesh_facets.getConnectivity(type_subsubfacet, gt_subsubfacet); sf_to_double = &mesh_facets.getData( "facet_to_double", type_subsubfacet, gt_subsubfacet); sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subsubfacet, gt_subsubfacet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subsubfacet, gt_subsubfacet); el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_subsubfacet, gt_subsubfacet); facet_to_subfacet = &mesh_facets.getElementToSubelement( type_subsubfacet, gt_subsubfacet); } UInt global_ssf = subsubfacet.element; Vector subsubfacet_connectivity( conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet( mesh_facets, starting_element, current_facet, subsubfacet_connectivity, element_list, facet_list, subfacet_list) == false && removeElementsInVector(*subfacet_list, (*facet_to_subfacet)(global_ssf)) == false) { sf_to_double->push_back(global_ssf); sf_to_subfacet_double->push_back(*subfacet_list); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } else { const UInt global_sf = subfacet.element; Vector subfacet_connectivity( conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subfacet is to be doubled if (findElementsAroundSubfacet( mesh_facets, starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false && removeElementsInVector( facet_list, (*facet_to_subfacet)(global_sf)) == false) { sf_to_double->push_back(global_sf); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } } } if (subsubfacet_mode) delete subfacet_list; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets, Array & new_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; MeshAccessor mesh_facets_accessor(mesh_facets); for (auto gt_facet : ghost_types) { for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); if (nb_facet_to_double == 0) continue; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); auto & facet_to_coh_element = mesh_facets_accessor.getSubelementToElement(type_cohesive, gt_facet); auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto & conn_cohesive = mesh.getConnectivity(type_cohesive, gt_facet); UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt old_nb_cohesive_elements = conn_cohesive.size(); UInt new_nb_cohesive_elements = conn_cohesive.size() + nb_facet_to_double; UInt old_nb_facet = element_to_facet.size() - nb_facet_to_double; facet_to_coh_element.resize(new_nb_cohesive_elements); conn_cohesive.resize(new_nb_cohesive_elements); UInt new_elements_old_size = new_elements.size(); new_elements.resize(new_elements_old_size + nb_facet_to_double); Element c_element{type_cohesive, 0, gt_facet}; Element f_element{type_facet, 0, gt_facet}; UInt facets[2]; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension ? 1 : 0] = f_to_double(facet); facets[third_dimension ? 0 : 1] = old_nb_facet + facet; UInt cohesive_element = old_nb_cohesive_elements + facet; /// store doubled facets f_element.element = facets[0]; facet_to_coh_element(cohesive_element, 0) = f_element; f_element.element = facets[1]; facet_to_coh_element(cohesive_element, 1) = f_element; /// modify cohesive elements' connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) { conn_cohesive(cohesive_element, n) = conn_facet(facets[0], n); conn_cohesive(cohesive_element, n + nb_nodes_per_facet) = conn_facet(facets[1], n); } /// update element_to_facet vectors c_element.element = cohesive_element; element_to_facet(facets[0])[1] = c_element; element_to_facet(facets[1])[1] = c_element; /// add cohesive element to the element event list new_elements(new_elements_old_size + facet) = c_element; } } } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); if (spatial_dimension != 1) return; auto & position = mesh.getNodes(); for (auto gt_facet : ghost_types) { - for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { - auto & conn_facet = - mesh_facets.getConnectivity(type_facet, gt_facet); + for (auto type_facet : + mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { + auto & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); const auto & facets_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); auto nb_facet_to_double = facets_to_double.size(); auto new_nb_facet = element_to_facet.size(); auto old_nb_facet = element_to_facet.size() - nb_facet_to_double; auto old_nb_nodes = position.size(); auto new_nb_nodes = old_nb_nodes + nb_facet_to_double; position.resize(new_nb_nodes); conn_facet.resize(new_nb_facet); auto old_nb_doubled_nodes = doubled_nodes.size(); doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double); - for(auto && data_facet : enumerate(facets_to_double)) { + for (auto && data_facet : enumerate(facets_to_double)) { const auto & old_facet = std::get<1>(data_facet); auto facet = std::get<0>(data_facet); auto new_facet = old_nb_facet + facet; auto el = element_to_facet(new_facet)[0]; auto old_node = conn_facet(old_facet); auto new_node = old_nb_nodes + facet; /// update position position(new_node) = position(old_node); conn_facet(new_facet) = new_node; Vector conn_segment = mesh.getConnectivity(el); /// update facet connectivity auto it = std::find(conn_segment.begin(), conn_segment.end(), old_node); *it = new_node; doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node; doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, GhostType gt_facet, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (type_facet != _segment_3) return; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.size(); UInt old_nb_facet = mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); /// this ones matter only for segments in 3D Array> * el_to_subfacet_double = nullptr; Array> * f_to_subfacet_double = nullptr; if (third_dim_segments) { el_to_subfacet_double = &mesh_facets.getData>( "elements_to_subfacet_double", type_facet, gt_facet); f_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_facet, gt_facet); } std::vector middle_nodes; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt node = conn_facet(old_facet, 2); if (!mesh.isPureGhostNode(node)) middle_nodes.push_back(node); } UInt n = doubled_nodes.size(); doubleNodes(mesh, middle_nodes, doubled_nodes); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt old_node = conn_facet(old_facet, 2); if (mesh.isPureGhostNode(old_node)) continue; UInt new_node = doubled_nodes(n, 1); UInt new_facet = old_nb_facet + facet; conn_facet(new_facet, 2) = new_node; if (third_dim_segments) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_facet(new_facet)); updateElementalConnectivity(mesh, old_node, new_node, (*el_to_subfacet_double)(facet), &(*f_to_subfacet_double)(facet)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_facet(new_facet)); } ++n; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); /// update subfacet_to_facet vector ElementType type_facet = _not_defined; GhostType gt_facet = _casper; Array * subfacet_to_facet = nullptr; UInt nb_subfacet_per_facet = 0; UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) - nb_subfacet_to_double; Array> * facet_list = nullptr; if (facet_mode) facet_list = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); else facet_list = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); Element old_subfacet_el{type_subfacet, 0, gt_subfacet}; Element new_subfacet_el{type_subfacet, 0, gt_subfacet}; for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { old_subfacet_el.element = sf_to_double(sf); new_subfacet_el.element = old_nb_subfacet + sf; for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) { Element & facet = (*facet_list)(sf)[f]; if (facet.type != type_facet || facet.ghost_type != gt_facet) { type_facet = facet.type; gt_facet = facet.ghost_type; subfacet_to_facet = &mesh_facets.getSubelementToElement(type_facet, gt_facet); nb_subfacet_per_facet = subfacet_to_facet->getNbComponent(); } Element * sf_update = std::find( subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet, subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, old_subfacet_el); AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet != subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, "Subfacet not found"); *sf_update = new_subfacet_el; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); Array> & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); Array> * facet_to_subfacet_double = nullptr; if (facet_mode) { facet_to_subfacet_double = &mesh_facets.getData>( "facets_to_subfacet_double", type_subfacet, gt_subfacet); } else { facet_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); } UInt old_nb_subfacet = facet_to_subfacet.size(); facet_to_subfacet.resize(old_nb_subfacet + nb_subfacet_to_double); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) facet_to_subfacet(old_nb_subfacet + sf) = (*facet_to_subfacet_double)(sf); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (spatial_dimension == 1) return; for (auto gt_subfacet : ghost_types) { for (auto type_subfacet : mesh_facets.elementTypes(0, gt_subfacet)) { auto & sf_to_double = mesh_facets.getData( "facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.size(); if (nb_subfacet_to_double == 0) continue; AKANTU_DEBUG_ASSERT( type_subfacet == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); auto & conn_subfacet = mesh_facets.getConnectivity(type_subfacet, gt_subfacet); UInt old_nb_subfacet = conn_subfacet.size(); UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double; conn_subfacet.resize(new_nb_subfacet); std::vector nodes_to_double; UInt old_nb_doubled_nodes = doubled_nodes.size(); /// double nodes for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_subfacet = sf_to_double(sf); nodes_to_double.push_back(conn_subfacet(old_subfacet)); } doubleNodes(mesh, nodes_to_double, doubled_nodes); /// add new nodes in connectivity for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt new_subfacet = old_nb_subfacet + sf; UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); conn_subfacet(new_subfacet) = new_node; } /// update facet and element connectivity Array> & f_to_subfacet_double = mesh_facets.getData>("facets_to_subfacet_double", type_subfacet, gt_subfacet); Array> & el_to_subfacet_double = mesh_facets.getData>( "elements_to_subfacet_double", type_subfacet, gt_subfacet); Array> * sf_to_subfacet_double = nullptr; if (spatial_dimension == 3) sf_to_subfacet_double = &mesh_facets.getData>( "subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0); UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); updateElementalConnectivity(mesh, old_node, new_node, el_to_subfacet_double(sf), &f_to_subfacet_double(sf)); updateElementalConnectivity(mesh_facets, old_node, new_node, f_to_subfacet_double(sf)); if (spatial_dimension == 3) updateElementalConnectivity(mesh_facets, old_node, new_node, (*sf_to_subfacet_double)(sf)); } if (spatial_dimension == 2) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true); } else if (spatial_dimension == 3) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::flipFacets( Mesh & mesh_facets, const ElementTypeMapArray & remote_global_connectivities, GhostType gt_facet) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); /// get global connectivity for local mesh ElementTypeMapArray local_global_connectivities( "local_global_connectivity", mesh_facets.getID(), mesh_facets.getMemoryID()); local_global_connectivities.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _ghost_type = gt_facet, _with_nb_nodes_per_element = true, _with_nb_element = true); mesh_facets.getGlobalConnectivity(local_global_connectivities); /// loop on every facet for (auto type_facet : mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) { auto & connectivity = mesh_facets.getConnectivity(type_facet, gt_facet); auto & local_global_connectivity = local_global_connectivities(type_facet, gt_facet); const auto & remote_global_connectivity = remote_global_connectivities(type_facet, gt_facet); auto & element_per_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); auto & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); auto nb_nodes_per_facet = connectivity.getNbComponent(); auto nb_nodes_per_P1_facet = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet)); for (auto && data : zip(make_view(connectivity, nb_nodes_per_facet), make_view(local_global_connectivity, nb_nodes_per_facet), make_view(remote_global_connectivity, nb_nodes_per_facet), make_view(subfacet_to_facet, subfacet_to_facet.getNbComponent()), make_view(element_per_facet))) { auto & conn = std::get<0>(data); auto & local_gconn = std::get<1>(data); const auto & remote_gconn = std::get<2>(data); /// skip facet if connectivities are the same if (local_gconn == remote_gconn) continue; /// re-arrange connectivity auto conn_tmp = conn; auto begin = local_gconn.begin(); auto end = local_gconn.end(); std::transform(remote_gconn.begin(), remote_gconn.end(), conn.begin(), [&](auto && gnode) { auto it = std::find(begin, end, gnode); AKANTU_DEBUG_ASSERT(it != end, "Node not found"); return conn_tmp(it - begin); }); /// if 3D, check if facets are just rotated if (spatial_dimension == 3) { auto begin = remote_gconn.storage(); /// find first node auto it = std::find(begin, begin + remote_gconn.size(), local_gconn(0)); UInt n, start = it - begin; /// count how many nodes in the received connectivity follow /// the same order of those in the local connectivity for (n = 1; n < nb_nodes_per_P1_facet && local_gconn(n) == remote_gconn((start + n) % nb_nodes_per_P1_facet); ++n) ; /// skip the facet inversion if facet is just rotated if (n == nb_nodes_per_P1_facet) continue; } /// update data to invert facet auto & element_per_facet = std::get<4>(data); std::swap(element_per_facet[0], element_per_facet[1]); auto & subfacets_of_facet = std::get<3>(data); std::swap(subfacets_of_facet(0), subfacets_of_facet(1)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if (mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call " "the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID(), mesh.getMemoryID()); barycenters.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = _all_dimensions); // mesh.initElementTypeMapArray(barycenters, spatial_dimension, // _all_dimensions); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(_all_dimensions, ghost_type)) { element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & barycenters_arr = barycenters(type, ghost_type); barycenters_arr.resize(nb_element); auto bary = barycenters_arr.begin(spatial_dimension); auto bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { element.element = el; mesh.getBarycenter(element, *bary); } } } MeshAccessor mesh_accessor(mesh); for (Int sp(spatial_dimension); sp >= 1; --sp) { if (mesh.getNbElement(sp) == 0) continue; for (auto ghost_type : ghost_types) { for (auto & type : mesh.elementTypes(sp, ghost_type)) { mesh_accessor.getSubelementToElement(type, ghost_type) - .resize(mesh.getNbElement(type, ghost_type), ElementNull); + .resize(mesh.getNbElement(type, ghost_type)); + mesh_accessor.getSubelementToElement(type, ghost_type).set(ElementNull); } for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { mesh_accessor.getElementToSubelement(type, ghost_type) .resize(mesh.getNbElement(type, ghost_type)); mesh.getElementToSubelement(type, ghost_type).clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (auto ghost_type : ghost_types) { facet_element.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(sp - 1, ghost_type)) { facet_element.type = type; auto & element_to_subelement = mesh.getElementToSubelement(type, ghost_type); const auto & connectivity = mesh.getConnectivity(type, ghost_type); - auto fit = connectivity.begin(mesh.getNbNodesPerElement(type)); - auto fend = connectivity.end(mesh.getNbNodesPerElement(type)); + for (auto && data : enumerate( + make_view(connectivity, mesh.getNbNodesPerElement(type)))) { + const auto & facet = std::get<1>(data); + facet_element.element = std::get<0>(data); - UInt fid = 0; - for (; fit != fend; ++fit, ++fid) { - const Vector & facet = *fit; - facet_element.element = fid; std::map element_seen_counter; - UInt nb_nodes_per_facet = + auto nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(type)); - for (UInt n(0); n < nb_nodes_per_facet; ++n) { - - auto eit = nodes_to_elements.begin(facet(n)); - auto eend = nodes_to_elements.end(facet(n)); - for (; eit != eend; ++eit) { - auto & elem = *eit; + // count the number of node in common between the facet and the other + // element connected to the nodes of the facet + for (auto node : arange(nb_nodes_per_facet)) { + for (auto & elem : nodes_to_elements.getRow(facet(node))) { auto cit = element_seen_counter.find(elem); if (cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } + // check which are the connected elements std::vector connected_elements; for (auto && cit : element_seen_counter) { if (cit.second == nb_nodes_per_facet) connected_elements.push_back(cit.first); } - for (auto & connected_element : connected_elements) - element_to_subelement(fid).push_back(connected_element); + // add the connected elements as sub-elements + for (auto & connected_element : connected_elements) { + element_to_subelement(facet_element.element) + .push_back(connected_element); + } - for (UInt ce = 0; ce < connected_elements.size(); ++ce) { - Element & elem = connected_elements[ce]; - Array & subelement_to_element = - mesh_accessor.getSubelementToElement(elem.type, - elem.ghost_type); + // add the element as sub-element to the connected elements + for (auto & connected_element : connected_elements) { + Vector subelements_to_element = + mesh.getSubelementToElement(connected_element); - UInt f(0); - for (; f < mesh.getNbFacetsPerElement(elem.type) && - subelement_to_element(elem.element, f) != ElementNull; - ++f) - ; + // find the position where to insert the element + auto it = std::find(subelements_to_element.begin(), + subelements_to_element.end(), ElementNull); AKANTU_DEBUG_ASSERT( - f < mesh.getNbFacetsPerElement(elem.type), - "The element " << elem << " seems to have too many facets!! (" - << f << " < " - << mesh.getNbFacetsPerElement(elem.type) << ")"); - - subelement_to_element(elem.element, f) = facet_element; + it != subelements_to_element.end(), + "The element " + << connected_element << " seems to have too many facets!! (" + << (it - subelements_to_element.begin()) << " < " + << mesh.getNbFacetsPerElement(connected_element.type) + << ")"); + + *it = facet_element; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MeshUtils::findElementsAroundSubfacet( const Mesh & mesh_facets, const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & element_list, std::vector & facet_list, std::vector * subfacet_list) { AKANTU_DEBUG_IN(); bool facet_matched = false; element_list.clear(); facet_list.clear(); if (third_dim_points) { subfacet_list->clear(); } element_list.push_back(starting_element); std::queue elements_to_check; elements_to_check.push(starting_element); /// keep going as long as there are elements to check while (not elements_to_check.empty()) { /// check current element Element & current_element = elements_to_check.front(); const Vector facets_to_element = mesh_facets.getSubelementToElement(current_element); // for every facet of the element for (auto & current_facet : facets_to_element) { if (current_facet == ElementNull) continue; if (current_facet == end_facet) facet_matched = true; // facet already listed if (std::find(facet_list.begin(), facet_list.end(), current_facet) != facet_list.end()) continue; // subfacet_connectivity is not in the connectivity of current_facet; if ((std::find(facet_list.begin(), facet_list.end(), current_facet) != facet_list.end()) or not hasElement(mesh_facets.getConnectivity(current_facet), subfacet_connectivity)) continue; facet_list.push_back(current_facet); if (third_dim_points) { const Vector subfacets_of_facet = mesh_facets.getSubelementToElement(current_facet); /// check subfacets for (const auto & current_subfacet : subfacets_of_facet) { if (current_subfacet == ElementNull) continue; if ((std::find(subfacet_list->begin(), subfacet_list->end(), current_subfacet) == subfacet_list->end()) and hasElement(mesh_facets.getConnectivity(current_subfacet), subfacet_connectivity)) subfacet_list->push_back(current_subfacet); } } /// consider opposing element const auto & elements_to_facet = mesh_facets.getElementToSubelement(current_facet); UInt opposing = 0; if (elements_to_facet[0] == current_element) opposing = 1; auto & opposing_element = elements_to_facet[opposing]; /// skip null elements since they are on a boundary if (opposing_element == ElementNull) continue; /// skip this element if already added if (std::find(element_list.begin(), element_list.end(), opposing_element) != element_list.end()) continue; /// only regular elements have to be checked if (opposing_element.kind() == _ek_regular) elements_to_check.push(opposing_element); element_list.push_back(opposing_element); AKANTU_DEBUG_ASSERT( hasElement(mesh_facets.getMeshParent().getConnectivity(opposing_element), subfacet_connectivity), "Subfacet doesn't belong to this element"); } /// erased checked element from the list elements_to_check.pop(); } AKANTU_DEBUG_OUT(); return facet_matched; } /* -------------------------------------------------------------------------- */ void MeshUtils::updateElementalConnectivity( Mesh & mesh, UInt old_node, UInt new_node, const std::vector & element_list, const std::vector * #if defined(AKANTU_COHESIVE_ELEMENT) facet_list #endif ) { AKANTU_DEBUG_IN(); for (auto & element : element_list) { if (element.type == _not_defined) continue; Vector connectivity = mesh.getConnectivity(element); #if defined(AKANTU_COHESIVE_ELEMENT) if (element.kind() == _ek_cohesive) { AKANTU_DEBUG_ASSERT( facet_list != nullptr, "Provide a facet list in order to update cohesive elements"); const Vector facets = mesh.getMeshFacets().getSubelementToElement(element); auto facet_nb_nodes = connectivity.size() / 2; /// loop over cohesive element's facets for (const auto & facet : enumerate(facets)) { /// skip facets if not present in the list if (std::find(facet_list->begin(), facet_list->end(), std::get<1>(facet)) == facet_list->end()) { continue; } auto n = std::get<0>(facet); auto begin = connectivity.begin() + n * facet_nb_nodes; auto end = begin + facet_nb_nodes; auto it = std::find(begin, end, old_node); AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element"); *it = new_node; } } else #endif { auto it = std::find(connectivity.begin(), connectivity.end(), old_node); AKANTU_DEBUG_ASSERT(it != connectivity.end(), "Node not found in current element"); /// update connectivity *it = new_node; } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 10ab5ef1b..608b782ec 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,899 +1,881 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Thu Feb 08 2018 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2015-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 "dof_manager_default.hh" #include "communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver_default.hh" #include "sparse_matrix_aij.hh" #include "terms_to_assemble.hh" #include "time_step_solver_default.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = i; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { if (c_jcn >= c_irn) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addElementalMatrixToUnsymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), global_equation_number(0, 1, "global_equation_number"), synchronizer(nullptr) {} /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), jacobian_release(0), global_equation_number(0, 1, "global_equation_number"), first_global_dof_id(0), synchronizer(nullptr) { if (this->mesh->isDistributed()) this->synchronizer = std::make_unique( *this, this->id + ":dof_synchronizer", this->memory_id); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() = default; /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = array_to_assemble.size() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.size() == nb_degree_of_freedoms, "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); typename Array::const_scalar_iterator arr_it = array_to_assemble.begin_reinterpret(nb_degree_of_freedoms); Array::const_scalar_iterator equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++arr_it, ++equ_it) { global_array(*equ_it) += scale_factor * (*arr_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id) : DOFData(dof_id), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManagerDefault::getNewDOFData(const ID & dof_id) { this->dofs[dof_id] = std::make_unique(dof_id); return *this->dofs[dof_id]; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs) { // auto prank = this->communicator.whoAmI(); // auto psize = this->communicator.getNbProc(); // access the relevant data to update auto & dof_data = this->getDOFDataTyped(dof_id); const auto & support_type = dof_data.support_type; const auto & group = dof_data.group_support; if (support_type == _dst_nodal and group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_dofs, nb_pure_local_dofs, [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs, [](UInt node) -> UInt { return node; }); } // update the synchronizer if needed if (this->synchronizer) this->synchronizer->registerDOFs(dof_id); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, support_type); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, group_support); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; std::unique_ptr sm = std::make_unique(*this, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrixAIJ & sm_to_copy = this->getMatrix(matrix_to_copy_id); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) { SparseMatrix & matrix = DOFManager::getMatrix(id); return dynamic_cast(matrix); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerDefault::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; std::unique_ptr nls; switch (type) { #if defined(AKANTU_IMPLICIT) case _nls_newton_raphson: case _nls_newton_raphson_modified: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } case _nls_linear: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } #endif case _nls_lumped: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } default: AKANTU_EXCEPTION("The asked type of non linear solver is not supported by " "this dof manager"); } return this->registerNonLinearSolver(non_linear_solver_id, nls); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique( *this, type, non_linear_solver, time_step_solver_id, this->memory_id); return this->registerTimeStepSolver(time_step_solver_id, tss); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearResidual() { this->residual.resize(this->local_system_size); this->residual.clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearMatrix(const ID & mtx) { this->getMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateGlobalBlockedDofs() { auto it = this->dofs.begin(); auto end = this->dofs.end(); this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.resize(this->local_system_size); this->global_blocked_dofs.clear(); for (; it != end; ++it) { if (!this->hasBlockedDOFs(it->first)) continue; DOFData & dof_data = *it->second; this->assembleToGlobalArray(it->first, *dof_data.blocked_dofs, this->global_blocked_dofs, true); } } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = equation_number.size(); local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent()); auto loc_it = local_array.begin_reinterpret(nb_degree_of_freedoms); auto equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++loc_it, ++equ_it) { (*loc_it) = global_array(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->global_solution, solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToResidual( const ID & dof_id, const Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); this->assembleToGlobalArray(dof_id, array_to_assemble, this->residual, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToLumpedMatrix( const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); Array & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { SparseMatrixAIJ & A = this->getMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, 1.); Array tmp_residual(this->residual.size(), 1, 0.); A.matVecMul(data_cache, tmp_residual, scale_factor, 1.); this->residual += tmp_residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleLumpedMatMulVectToResidual( const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { const Array & A = this->getLumpedMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, scale_factor); auto A_it = A.begin(); auto A_end = A.end(); auto x_it = data_cache.begin(); auto r_it = this->residual.begin(); for (; A_it != A_end; ++A_it, ++x_it, ++r_it) { *r_it += *A_it * *x_it; } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); DOFData & dof_data = this->getDOFData(dof_id); AKANTU_DEBUG_ASSERT(dof_data.support_type == _dst_nodal, "This function applies only on Nodal dofs"); this->addToProfile(matrix_id, dof_id, type, ghost_type); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt nb_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); nb_element = group_elements.size(); filter_it = group_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } } AKANTU_DEBUG_ASSERT(elementary_mat.size() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = dof_data.dof->getNbComponent(); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); auto conn_begin = connectivity.begin(nb_nodes_per_element); auto conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (filter_it) ++filter_it; else ++conn_it; if (A.getMatrixType() == _symmetric) if (elemental_matrix_type == _symmetric) this->addSymmetricElementalMatrixToSymmetric(A, *el_mat_it, element_eq_nb, A.size()); else this->addUnsymmetricElementalMatrixToSymmetric(A, *el_mat_it, element_eq_nb, A.size()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb, A.size()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { const Array & equation_number_m = this->getLocalEquationNumbers(dof_id_m); const Array & equation_number_n = this->getLocalEquationNumbers(dof_id_n); SparseMatrixAIJ & A = this->getMatrix(matrix_id); for (const auto & term : terms) { UInt gi = this->localToGlobalEquationNumber(equation_number_m(term.i())); UInt gj = this->localToGlobalEquationNumber(equation_number_n(term.j())); A.add(gi, gj, term); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const DOFData & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) return; auto mat_dof = std::make_pair(matrix_id, dof_id); auto type_pair = std::make_pair(type, ghost_type); auto prof_it = this->matrix_profiled_dofs.find(mat_dof); if (prof_it != this->matrix_profiled_dofs.end() && std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) != prof_it->second.end()) return; UInt nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt size = A.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto & connectivity = this->mesh->getConnectivity(type, ghost_type); auto cbegin = connectivity.begin(nb_nodes_per_element); auto cit = cbegin; UInt nb_elements = connectivity.size(); UInt * ge_it = nullptr; if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); ge_it = group_elements.storage(); nb_elements = group_elements.size(); } UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; Vector element_eq_nb(size_mat); for (UInt e = 0; e < nb_elements; ++e) { if (ge_it) cit = cbegin + *ge_it; this->extractElementEquationNumber( equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (ge_it) ++ge_it; else ++cit; for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = element_eq_nb(i); if (c_irn < size) { for (UInt j = 0; j < size_mat; ++j) { UInt c_jcn = element_eq_nb(j); if (c_jcn < size) { A.add(c_irn, c_jcn); } } } } } this->matrix_profiled_dofs[mat_dof].push_back(type_pair); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::applyBoundary(const ID & matrix_id) { SparseMatrixAIJ & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getValueRelease()) { auto are_equal = std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (not are_equal) J.applyBoundary(); previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getValueRelease(); } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getGlobalResidual() { if (this->synchronizer) { if (not this->global_residual) { this->global_residual = std::make_unique>(0, 1, "global_residual"); } if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->global_residual->resize(this->system_size); this->synchronizer->gather(this->residual, *this->global_residual); } else { this->synchronizer->gather(this->residual); } return *this->global_residual; } else { return this->residual; } } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getResidual() const { return this->residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::setGlobalSolution(const Array & solution) { if (this->synchronizer) { if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->synchronizer->scatter(this->global_solution, solution); } else { this->synchronizer->scatter(this->global_solution); } } else { AKANTU_DEBUG_ASSERT(solution.size() == this->global_solution.size(), "Sequential call to this function needs the solution " "to be the same size as the global_solution"); this->global_solution.copy(solution); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { DOFManager::onNodesAdded(nodes_list, event); if (this->synchronizer) this->synchronizer->onNodesAdded(nodes_list); } /* -------------------------------------------------------------------------- */ std::pair DOFManagerDefault::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { UInt nb_new_local_dofs, nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = DOFManager::updateNodalDOFs(dof_id, nodes_list); auto & dof_data = this->getDOFDataTyped(dof_id); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor(DOFManagerDefault::DOFDataDefault & dof_data, DOFManagerDefault & dof_manager) : dof_data(dof_data), dof_manager(dof_manager) { for (auto && pair : zip(dof_data.local_equation_number, dof_data.associated_nodes)) { UInt node, dof; std::tie(dof, node) = pair; dofs_per_node[node].push_back(dof); } } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_ask_nodes) { return nodes.size() * dof_data.dof->getNbComponent() * sizeof(UInt); } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_ask_nodes) { for (auto & node : nodes) { auto & dofs = dofs_per_node.at(node); for (auto & dof : dofs) { buffer << dof_manager.global_equation_number(dof); } } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { if (tag == _gst_ask_nodes) { for (auto & node : nodes) { auto & dofs = dofs_per_node[node]; for (auto dof : dofs) { UInt global_dof; buffer >> global_dof; dof_manager.global_equation_number(dof) = global_dof; dof_manager.global_to_local_mapping[global_dof] = dof; } } } } protected: std::unordered_map> dofs_per_node; DOFManagerDefault::DOFDataDefault & dof_data; DOFManagerDefault & dof_manager; }; /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateDOFsData( DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, const std::function & getNode) { auto prank = this->communicator.whoAmI(); auto psize = this->communicator.getNbProc(); // access the relevant data to update const auto & support_type = dof_data.support_type; // resize all relevant arrays this->residual.resize(this->local_system_size, 0.); this->dofs_type.resize(this->local_system_size); this->global_solution.resize(this->local_system_size, 0.); this->global_blocked_dofs.resize(this->local_system_size, true); this->previous_global_blocked_dofs.resize(this->local_system_size, true); this->global_equation_number.resize(this->local_system_size, -1); for (auto & lumped_matrix : lumped_matrices) lumped_matrix.second->resize(this->local_system_size); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); // determine the first local/global dof id to use Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = nb_new_pure_local; this->communicator.allGather(nb_dofs_per_proc); this->first_global_dof_id += std::accumulate( nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0); UInt first_dof_id = this->local_system_size - nb_new_local_dofs; if (support_type == _dst_nodal) { dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_new_local_dofs); } // update per dof info for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number.push_back(local_eq_num); bool is_local_dof = true; // determine the dof type switch (support_type) { case _dst_nodal: { UInt node = getNode(d / dof_data.dof->getNbComponent()); this->dofs_type(local_eq_num) = this->mesh->getNodeType(node); dof_data.associated_nodes.push_back(node); + is_local_dof = this->mesh->isLocalOrMasterNode(node); break; } case _dst_generic: { this->dofs_type(local_eq_num) = _nt_normal; break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } // update global id for local dofs if (is_local_dof) { this->global_equation_number(local_eq_num) = this->first_global_dof_id; this->global_to_local_mapping[this->first_global_dof_id] = local_eq_num; ++this->first_global_dof_id; } else { this->global_equation_number(local_eq_num) = UInt(-1); } } // synchronize the global numbering for slaves if (support_type == _dst_nodal && this->synchronizer) { - //auto nb_dofs_per_node = dof_data.dof->getNbComponent(); - GlobalDOFInfoDataAccessor data_accessor(dof_data, *this); auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, _gst_ask_nodes); - - // std::vector counters(nb_new_local_dofs / nb_dofs_per_node); - - // for (UInt d = 0; d < nb_new_local_dofs; ++d) { - // UInt local_eq_num = first_dof_id + d; - // if (not this->isSlaveDOF(local_eq_num)) - // continue; - - // UInt node = d / nb_dofs_per_node; - // UInt dof_count = counters[node]++; - // node = getNode(node); - - // UInt global_dof_id = data_accessor.getNthDOFForNode(dof_count, node); - - // this->global_equation_number(local_eq_num) = global_dof_id; - // this->global_to_local_mapping[global_dof_id] = local_eq_num; - // } } this->first_global_dof_id += std::accumulate( nb_dofs_per_proc.begin() + prank + 1, nb_dofs_per_proc.end(), 0); } /* -------------------------------------------------------------------------- */ // register in factory static bool default_dof_manager_is_registered[[gnu::unused]] = DefaultDOFManagerFactory::getInstance().registerAllocator( "default", [](const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(id, mem_id); }); static bool dof_manager_is_registered[[gnu::unused]] = DOFManagerFactory::getInstance().registerAllocator( "default", [](Mesh & mesh, const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(mesh, id, mem_id); }); } // namespace akantu diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc index 029055f20..035ccfb53 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc @@ -1,333 +1,334 @@ /** * @file test_solid_mechanics_model_dynamics.cc * * @author Guillaume Anciaux * * @date creation: Wed Nov 29 2017 * @date last modification: Wed Feb 21 2018 * * @brief test of the class SolidMechanicsModel on the 3d cube * * @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 "boundary_condition_functor.hh" #include "test_solid_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { const Real A = 1e-1; const Real E = 1.; const Real poisson = 3. / 10; const Real lambda = E * poisson / (1 + poisson) / (1 - 2 * poisson); const Real mu = E / 2 / (1. + poisson); const Real rho = 1.; const Real cp = std::sqrt((lambda + 2 * mu) / rho); const Real cs = std::sqrt(mu / rho); const Real c = std::sqrt(E / rho); const Vector k = {.5, 0., 0.}; const Vector psi1 = {0., 0., 1.}; const Vector psi2 = {0., 1., 0.}; const Real knorm = k.norm(); /* -------------------------------------------------------------------------- */ template struct Verification {}; /* -------------------------------------------------------------------------- */ template <> struct Verification<1> { void displacement(Vector & disp, const Vector & coord, Real current_time) { const auto & x = coord(_x); const Real omega = c * k[0]; disp(_x) = A * cos(k[0] * x - omega * current_time); } void velocity(Vector & vel, const Vector & coord, Real current_time) { const auto & x = coord(_x); const Real omega = c * k[0]; vel(_x) = omega * A * sin(k[0] * x - omega * current_time); } }; /* -------------------------------------------------------------------------- */ template <> struct Verification<2> { void displacement(Vector & disp, const Vector & X, Real current_time) { Vector kshear = {k[1], k[0]}; Vector kpush = {k[0], k[1]}; const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; disp = A * (kpush * cos(phase_p) + kshear * cos(phase_s)); } void velocity(Vector & vel, const Vector & X, Real current_time) { Vector kshear = {k[1], k[0]}; Vector kpush = {k[0], k[1]}; const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; vel = A * (kpush * omega_p * cos(phase_p) + kshear * omega_s * cos(phase_s)); } }; /* -------------------------------------------------------------------------- */ template <> struct Verification<3> { void displacement(Vector & disp, const Vector & coord, Real current_time) { const auto & X = coord; Vector kpush = k; Vector kshear1(3); Vector kshear2(3); kshear1.crossProduct(k, psi1); kshear2.crossProduct(k, psi2); const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; disp = A * (kpush * cos(phase_p) + kshear1 * cos(phase_s) + kshear2 * cos(phase_s)); } void velocity(Vector & vel, const Vector & coord, Real current_time) { const auto & X = coord; Vector kpush = k; Vector kshear1(3); Vector kshear2(3); kshear1.crossProduct(k, psi1); kshear2.crossProduct(k, psi2); const Real omega_p = knorm * cp; const Real omega_s = knorm * cs; Real phase_p = X.dot(kpush) - omega_p * current_time; Real phase_s = X.dot(kpush) - omega_s * current_time; vel = A * (kpush * omega_p * cos(phase_p) + kshear1 * omega_s * cos(phase_s) + kshear2 * omega_s * cos(phase_s)); } }; /* -------------------------------------------------------------------------- */ template class SolutionFunctor : public BC::Dirichlet::DirichletFunctor { public: SolutionFunctor(Real current_time, SolidMechanicsModel & model) : current_time(current_time), model(model) {} public: static constexpr UInt dim = ElementClass<_type>::getSpatialDimension(); inline void operator()(UInt node, Vector & flags, Vector & primal, const Vector & coord) const { flags(0) = true; auto & vel = model.getVelocity(); auto it = vel.begin(model.getSpatialDimension()); Vector v = it[node]; Verification verif; verif.displacement(primal, coord, current_time); verif.velocity(v, coord, current_time); } private: Real current_time; SolidMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ // This fixture uses somewhat finer meshes representing bars. template class TestSMMFixtureBar : public TestSMMFixture { using parent = TestSMMFixture; public: void SetUp() override { this->mesh_file = "../patch_tests/data/bar" + aka::to_string(this->type) + ".msh"; parent::SetUp(); getStaticParser().parse("test_solid_mechanics_model_" "dynamics_material.dat"); auto analysis_method = analysis_method_::value; this->model->initFull(_analysis_method = analysis_method); if (this->dump_paraview) { std::stringstream base_name; base_name << "bar" << analysis_method << this->type; this->model->setBaseName(base_name.str()); this->model->addDumpFieldVector("displacement"); + this->model->addDumpFieldVector("blocked_dofs"); this->model->addDumpField("mass"); this->model->addDumpField("velocity"); this->model->addDumpField("acceleration"); this->model->addDumpFieldVector("external_force"); this->model->addDumpFieldVector("internal_force"); this->model->addDumpField("stress"); this->model->addDumpField("strain"); } time_step = this->model->getStableTimeStep() / 10.; this->model->setTimeStep(time_step); // std::cout << "timestep: " << time_step << std::endl; const auto & position = this->mesh->getNodes(); auto & displacement = this->model->getDisplacement(); auto & velocity = this->model->getVelocity(); constexpr auto dim = parent::spatial_dimension; Verification verif; for (auto && tuple : zip(make_view(position, dim), make_view(displacement, dim), make_view(velocity, dim))) { verif.displacement(std::get<1>(tuple), std::get<0>(tuple), 0.); verif.velocity(std::get<2>(tuple), std::get<0>(tuple), 0.); } if (dump_paraview) this->model->dump(); /// boundary conditions this->model->applyBC(SolutionFunctor(0., *this->model), "BC"); wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1 simulation_time = 5 / wave_velocity; max_steps = simulation_time / time_step; // 100 auto ndump = 50; dump_freq = max_steps / ndump; } void solveStep() { constexpr auto dim = parent::spatial_dimension; Real current_time = 0.; const auto & position = this->mesh->getNodes(); const auto & displacement = this->model->getDisplacement(); UInt nb_nodes = this->mesh->getNbNodes(); UInt nb_global_nodes = this->mesh->getNbGlobalNodes(); max_error = -1.; Array displacement_solution(nb_nodes, dim); Verification verif; for (UInt s = 0; s < this->max_steps; ++s, current_time += this->time_step) { if (s % this->dump_freq == 0 && this->dump_paraview) this->model->dump(); /// boundary conditions this->model->applyBC( SolutionFunctor(current_time, *this->model), "BC"); // compute the disp solution for (auto && tuple : zip(make_view(position, dim), make_view(displacement_solution, dim))) { verif.displacement(std::get<1>(tuple), std::get<0>(tuple), current_time); } // compute the error solution Real disp_error = 0.; auto n = 0; for (auto && tuple : zip(make_view(displacement, dim), make_view(displacement_solution, dim))) { if (this->mesh->isLocalOrMasterNode(n)) { auto diff = std::get<1>(tuple) - std::get<0>(tuple); disp_error += diff.dot(diff); } ++n; } this->mesh->getCommunicator().allReduce(disp_error, SynchronizerOperation::_sum); disp_error = sqrt(disp_error) / nb_global_nodes; max_error = std::max(disp_error, max_error); this->model->solveStep(); } } protected: Real time_step; Real wave_velocity; Real simulation_time; UInt max_steps; UInt dump_freq; bool dump_paraview{false}; Real max_error{-1}; }; template using analysis_method_t = std::integral_constant; using TestTypes = gtest_list_t; template using TestSMMFixtureBarExplicit = TestSMMFixtureBar>; TYPED_TEST_CASE(TestSMMFixtureBarExplicit, TestTypes); /* -------------------------------------------------------------------------- */ TYPED_TEST(TestSMMFixtureBarExplicit, Dynamics) { this->solveStep(); EXPECT_NEAR(this->max_error, 0., 2e-3); // std::cout << "max error: " << max_error << std::endl; } /* -------------------------------------------------------------------------- */ template using TestSMMFixtureBarImplicit = TestSMMFixtureBar>; TYPED_TEST_CASE(TestSMMFixtureBarImplicit, TestTypes); TYPED_TEST(TestSMMFixtureBarImplicit, Dynamics) { this->solveStep(); EXPECT_NEAR(this->max_error, 0., 2e-3); } }