Page MenuHomec4science
Files F63573127

mesh_utils_pbc.cc
No OneTemporary
Actions

File Metadata

Created
Tue, May 21, 02:23

mesh_utils_pbc.cc
View Options

/**
* @file mesh_utils_pbc.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
*
* @date creation: Wed Feb 09 2011
* @date last modification: Fri Aug 01 2014
*
* @brief periodic boundary condition connectivity tweak
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <map>
/* -------------------------------------------------------------------------- */
#include "mesh_utils.hh"
#include "element_group.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 dirx, const UInt diry,
Real normalization,
Real tolerance,
Real * coords):
dim(dimension),dir_x(dirx),dir_y(diry),normalization(normalization),
tolerance(tolerance),coordinates(coords){}
// answers the question whether n2 is larger or equal to n1
bool operator() (UInt n1, UInt n2){
Real p1_x = coordinates[dim*n1+dir_x];
Real p2_x = coordinates[dim*n2+dir_x];
Real diff_x = p1_x - p2_x;
if (dim == 2 || std::abs(diff_x)/normalization > tolerance)
return diff_x > 0.0 ? false : true;
else if (dim > 2){
Real p1_y = coordinates[dim*n1+dir_y];
Real p2_y = coordinates[dim*n2+dir_y];
Real diff_y = p1_y - p2_y;
return diff_y > 0 ? false : true;
}
return true;
}
private:
UInt dim;
UInt dir_x;
UInt dir_y;
Real normalization;
Real tolerance;
Real * coordinates;
};
/* -------------------------------------------------------------------------- */
void MeshUtils::computePBCMap(const Mesh & mymesh,
const UInt dir,
std::map<UInt,UInt> & pbc_pair){
Array<UInt> selected_left;
Array<UInt> selected_right;
Real * coords = mymesh.nodes->storage();
const UInt nb_nodes = mymesh.nodes->getSize();
const UInt dim = mymesh.getSpatialDimension();
if (dim <= dir) return;
AKANTU_DEBUG_INFO("min " << mymesh.lower_bounds[dir]);
AKANTU_DEBUG_INFO("max " << mymesh.upper_bounds[dir]);
for (UInt i = 0; i < nb_nodes; ++i) {
AKANTU_DEBUG_TRACE("treating " << coords[dim*i+dir]);
if(Math::are_float_equal(coords[dim*i+dir], mymesh.lower_bounds[dir])){
AKANTU_DEBUG_TRACE("pushing node " << i << " on the left side");
selected_left.push_back(i);
}
else if(Math::are_float_equal(coords[dim*i+dir], mymesh.upper_bounds[dir])){
selected_right.push_back(i);
AKANTU_DEBUG_TRACE("pushing node " << i << " 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(mymesh, dir, selected_left, selected_right, pbc_pair);
}
/* -------------------------------------------------------------------------- */
void MeshUtils::computePBCMap(const Mesh & mymesh,
const SurfacePair & surface_pair,
std::map<UInt,UInt> & pbc_pair) {
Array<UInt> selected_first;
Array<UInt> selected_second;
// find nodes on surfaces
const ElementGroup & first_surf = mymesh.getElementGroup(surface_pair.first);
const ElementGroup & second_surf = mymesh.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<Real> & coords = mymesh.getNodes();
const UInt dim = mymesh.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<Real>::max();
second_min[i] = std::numeric_limits<Real>::max();
first_max[i] = -std::numeric_limits<Real>::max();
second_max[i] = -std::numeric_limits<Real>::max();
}
// find min and max of surface nodes
for (Array<UInt>::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<UInt>::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(mymesh, dir, selected_first, selected_second, pbc_pair);
else
MeshUtils::matchPBCPairs(mymesh, dir, selected_second, selected_first, pbc_pair);
}
/* -------------------------------------------------------------------------- */
void MeshUtils::matchPBCPairs(const Mesh & mymesh,
const UInt dir,
Array<UInt> & selected_left,
Array<UInt> & selected_right,
std::map<UInt,UInt> & 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;
Real * coords = mymesh.nodes->storage();
const UInt dim = mymesh.getSpatialDimension();
Real normalization = mymesh.upper_bounds[dir]-mymesh.lower_bounds[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 dir_x = UInt(-1) ,dir_y = UInt(-1);
if (dim == 3){
if (dir == 0){
dir_x = 1;dir_y = 2;
}
else if (dir == 1){
dir_x = 0;dir_y = 2;
}
else if (dir == 2){
dir_x = 0;dir_y = 1;
}
}
else if (dim == 2){
if (dir == 0){
dir_x = 1;
}
else if (dir == 1){
dir_x = 0;
}
}
CoordinatesComparison compare_nodes(dim,dir_x,dir_y,normalization,tolerance,coords);
std::sort(selected_left.begin(),selected_left.end(),compare_nodes);
std::sort(selected_right.begin(),selected_right.end(),compare_nodes);
Array<UInt>::scalar_iterator it_left = selected_left.begin();
Array<UInt>::scalar_iterator end_left = selected_left.end();
Array<UInt>::scalar_iterator it_right = selected_right.begin();
Array<UInt>::scalar_iterator end_right = selected_right.end();
while ((it_left != end_left) && (it_right != end_right) ){
UInt i1 = *it_left;
UInt i2 = *it_right;
AKANTU_DEBUG_TRACE("do I pair? " << i1 << "("
<< coords[dim*i1] << "," << coords[dim*i1+1] << ","
<< coords[dim*i1+2]
<< ") with"
<< i2 << "("
<< coords[dim*i2] << "," << coords[dim*i2+1] << ","
<< coords[dim*i2+2]
<< ") in direction " << dir);
Real dx = 0.0;
Real dy = 0.0;
if (dim >= 2) dx = coords[dim*i1 + dir_x] - coords[dim*i2 + dir_x];
if (dim == 3) dy = coords[dim*i1 + dir_y] - coords[dim*i2 + dir_y];
if (fabs(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 << "("
<< coords[dim*i1] << "," << coords[dim*i1+1] << ","
<< coords[dim*i1+2]
<< ") with"
<< i2 << "("
<< coords[dim*i2] << "," << coords[dim*i2+1] << ","
<< coords[dim*i2+2]
<< ") in direction " << dir);
++it_left;
++it_right;
} else if (compare_nodes(i1, i2)) {
++it_left;
} else {
++it_right;
}
}
AKANTU_DEBUG_INFO("found " << pbc_pair.size() << " pairs for direction " << dir);
}
__END_AKANTU__

Event Timeline

c4science · Help