diff --git a/src/mesh_utils/mesh_utils_pbc.cc b/src/mesh_utils/mesh_utils_pbc.cc index ab934b111..ea00ebf40 100644 --- a/src/mesh_utils/mesh_utils_pbc.cc +++ b/src/mesh_utils/mesh_utils_pbc.cc @@ -1,291 +1,298 @@ /** * @file mesh_utils_pbc.cc * * @author Guillaume Anciaux * @author David Simon Kammer * * @date creation: Wed Feb 09 2011 * @date last modification: Tue Sep 15 2015 * * @brief periodic boundary condition connectivity tweak * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /// class that sorts a set of nodes of same coordinates in 'dir' direction class CoordinatesComparison { public: CoordinatesComparison(const UInt dimension, const UInt dir_1, const UInt dir_2, - Real normalization, Real tolerance) - : dim(dimension), dir_1(dir_1), dir_2(dir_2), normalization(normalization), - tolerance(tolerance) {} + Real normalization, Real tolerance, const Array & coords) + : dim(dimension), dir_1(dir_1), dir_2(dir_2), normalization(normalization), + tolerance(tolerance), coords_it(coords.begin(dim)) {} // answers the question whether n2 is larger or equal to n1 - bool operator()(const Vector & n1, const Vector & n2) { - Real diff = n1(dir_1) - n2(dir_1);; + bool operator()(const UInt n1, const UInt n2) { + Vector coords_n1 = coords_it[n1]; + Vector coords_n2 = coords_it[n2]; + return this->operator()(coords_n1, coords_n2); + } + + bool operator()(const Vector & coords_n1, const Vector & coords_n2) { + Real diff = coords_n1(dir_1) - coords_n2(dir_1);; if (dim == 2 || std::abs(diff) / normalization > tolerance) return diff > 0. ? false : true; else if (dim > 2) { - diff = n1(dir_2) - n2(dir_2);; - return diff > 0 ? false : true; + diff = coords_n1(dir_2) - coords_n2(dir_2);; + return diff > 0 ? false : true; } return true; } private: UInt dim; UInt dir_1; UInt dir_2; Real normalization; Real tolerance; + const Array::const_vector_iterator coords_it; }; /* -------------------------------------------------------------------------- */ void MeshUtils::computePBCMap(const Mesh & mesh, const UInt dir, std::map & pbc_pair) { Array selected_left; Array selected_right; const UInt dim = mesh.getSpatialDimension(); Array::const_vector_iterator it = mesh.getNodes().begin(dim); Array::const_vector_iterator end = mesh.getNodes().end(dim); if (dim <= dir) return; const Vector & lower_bounds = mesh.getLowerBounds(); const Vector & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_INFO("min " << lower_bounds(dir)); AKANTU_DEBUG_INFO("max " << upper_bounds(dir)); for (UInt node = 0; it != end; ++it, ++node) { const Vector & coords = *it; AKANTU_DEBUG_TRACE("treating " << coords(dir)); if (Math::are_float_equal(coords(dir), lower_bounds(dir))) { AKANTU_DEBUG_TRACE("pushing node " << node << " on the left side"); selected_left.push_back(node); } else if (Math::are_float_equal(coords(dir), upper_bounds(dir))) { selected_right.push_back(node); AKANTU_DEBUG_TRACE("pushing node " << node << " on the right side"); } } AKANTU_DEBUG_INFO( "found " << selected_left.getSize() << " and " << selected_right.getSize() << " nodes at each boundary for direction " << dir); // match pairs MeshUtils::matchPBCPairs(mesh, dir, selected_left, selected_right, pbc_pair); } /* -------------------------------------------------------------------------- */ void MeshUtils::computePBCMap(const Mesh & mesh, const SurfacePair & surface_pair, std::map & pbc_pair) { Array selected_first; Array selected_second; // find nodes on surfaces const ElementGroup & first_surf = mesh.getElementGroup(surface_pair.first); const ElementGroup & second_surf = mesh.getElementGroup(surface_pair.second); // if this surface pair is not on this proc if (first_surf.getNbNodes() == 0 || second_surf.getNbNodes() == 0) { AKANTU_DEBUG_WARNING("computePBCMap has at least one surface without any " "nodes. I will ignore it."); return; } // copy nodes from element group selected_first.copy(first_surf.getNodeGroup().getNodes()); selected_second.copy(second_surf.getNodeGroup().getNodes()); // coordinates const Array & coords = mesh.getNodes(); const UInt dim = mesh.getSpatialDimension(); // variables to find min and max of surfaces Real first_max[3], first_min[3]; Real second_max[3], second_min[3]; for (UInt i = 0; i < dim; ++i) { first_min[i] = std::numeric_limits::max(); second_min[i] = std::numeric_limits::max(); first_max[i] = -std::numeric_limits::max(); second_max[i] = -std::numeric_limits::max(); } // find min and max of surface nodes for (Array::scalar_iterator it = selected_first.begin(); it != selected_first.end(); ++it) { for (UInt i = 0; i < dim; ++i) { if (first_min[i] > coords(*it, i)) first_min[i] = coords(*it, i); if (first_max[i] < coords(*it, i)) first_max[i] = coords(*it, i); } } for (Array::scalar_iterator it = selected_second.begin(); it != selected_second.end(); ++it) { for (UInt i = 0; i < dim; ++i) { if (second_min[i] > coords(*it, i)) second_min[i] = coords(*it, i); if (second_max[i] < coords(*it, i)) second_max[i] = coords(*it, i); } } // find direction of pbc Int first_dir = -1; #ifndef AKANTU_NDEBUG Int second_dir = -2; #endif for (UInt i = 0; i < dim; ++i) { if (Math::are_float_equal(first_min[i], first_max[i])) { first_dir = i; } #ifndef AKANTU_NDEBUG if (Math::are_float_equal(second_min[i], second_max[i])) { second_dir = i; } #endif } AKANTU_DEBUG_ASSERT(first_dir == second_dir, "Surface pair has not same direction. Surface " << surface_pair.first << " dir=" << first_dir << " ; Surface " << surface_pair.second << " dir=" << second_dir); UInt dir = first_dir; // match pairs if (first_min[dir] < second_min[dir]) MeshUtils::matchPBCPairs(mesh, dir, selected_first, selected_second, pbc_pair); else MeshUtils::matchPBCPairs(mesh, dir, selected_second, selected_first, pbc_pair); } /* -------------------------------------------------------------------------- */ void MeshUtils::matchPBCPairs(const Mesh & mesh, const UInt dir, Array & selected_left, Array & selected_right, std::map & pbc_pair) { // tolerance is that large because most meshers generate points coordinates // with single precision only (it is the case of GMSH for instance) Real tolerance = 1e-7; const UInt dim = mesh.getSpatialDimension(); Real normalization = mesh.getUpperBounds()(dir) - mesh.getLowerBounds()(dir); AKANTU_DEBUG_ASSERT(std::abs(normalization) > Math::getTolerance(), "In matchPBCPairs: The normalization is zero. " << "Did you compute the bounding box of the mesh?"); UInt odir_1 = UInt(-1), odir_2 = UInt(-1); if (dim == 3) { if (dir == _x) { odir_1 = _y; odir_2 = _z; } else if (dir == _y) { odir_1 = _x; odir_2 = _z; } else if (dir == _z) { odir_1 = _x; odir_2 = _y; } } else if (dim == 2) { if (dir == _x) { odir_1 = _y; } else if (dir == _y) { odir_1 = _x; } } CoordinatesComparison compare_nodes(dim, odir_1, odir_2, normalization, - tolerance); + tolerance, mesh.getNodes()); std::sort(selected_left.begin(), selected_left.end(), compare_nodes); std::sort(selected_right.begin(), selected_right.end(), compare_nodes); Array::scalar_iterator it_left = selected_left.begin(); Array::scalar_iterator end_left = selected_left.end(); Array::scalar_iterator it_right = selected_right.begin(); Array::scalar_iterator end_right = selected_right.end(); Array::const_vector_iterator nit = mesh.getNodes().begin(dim); while ((it_left != end_left) && (it_right != end_right)) { UInt i1 = *it_left; UInt i2 = *it_right; Vector coords1 = nit[i1]; Vector coords2 = nit[i2]; AKANTU_DEBUG_TRACE("do I pair? " << i1 << "(" << coords1 << ") with" << i2 << "(" << coords2 << ") in direction " << dir); Real dx = 0.0; Real dy = 0.0; if (dim >= 2) dx = coords1(odir_1) - coords2(odir_1); if (dim == 3) dy = coords1(odir_2) - coords2(odir_2); if (std::abs(dx * dx + dy * dy) / normalization < tolerance) { // then i match these pairs if (pbc_pair.count(i2)) { i2 = pbc_pair[i2]; } pbc_pair[i1] = i2; AKANTU_DEBUG_TRACE("pairing " << i1 << "(" << coords1 << ") with" << i2 << "(" << coords2 << ") in direction " << dir); ++it_left; ++it_right; } else if (compare_nodes(coords1, coords2)) { ++it_left; } else { ++it_right; } } AKANTU_DEBUG_INFO("found " << pbc_pair.size() << " pairs for direction " << dir); } __END_AKANTU__