Homec4science
Diffusion akantu (features/11-computation-of-contact-energy-with-getenergy-contact) Browse

rAKA/test/test_synchronizer/test_grid_synchronizer.cc
a51a521437e3features/11-computation…

Last Change
Raw File

test_grid_synchronizer.cc
View Options

/**
* @file test_grid_synchronizer.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Nov 25 2011
* @date last modification: Tue Jun 24 2014
*
* @brief test the GridSynchronizer object
*
* @section LICENSE
*
* Copyright (©) 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 "aka_common.hh"
#include "aka_grid_dynamic.hh"
#include "mesh.hh"
#include "grid_synchronizer.hh"
#include "mesh_partition.hh"
#include "synchronizer_registry.hh"
#include "test_data_accessor.hh"
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.hh"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
const UInt spatial_dimension = 3;
typedef std::map<std::pair<Element, Element>, Real> pair_list;
#include "test_grid_tools.hh"
static void updatePairList(const ElementTypeMapArray<Real> & barycenter,
const SpatialGrid<Element> & grid,
Real radius,
pair_list & neighbors,
neighbors_map_t<spatial_dimension>::type & neighbors_map) {
AKANTU_DEBUG_IN();
GhostType ghost_type = _not_ghost;
Element e;
e.ghost_type = ghost_type;
// generate the pair of neighbor depending of the cell_list
ElementTypeMapArray<Real>::type_iterator it = barycenter.firstType(_all_dimensions, ghost_type);
ElementTypeMapArray<Real>::type_iterator last_type = barycenter.lastType(0, ghost_type);
for(; it != last_type; ++it) {
// loop over quad points
e.type = *it;
e.element = 0;
const Array<Real> & barycenter_vect = barycenter(*it, ghost_type);
UInt sp = barycenter_vect.getNbComponent();
Array<Real>::const_iterator< Vector<Real> > bary =
barycenter_vect.begin(sp);
Array<Real>::const_iterator< Vector<Real> > bary_end =
barycenter_vect.end(sp);
for(;bary != bary_end; ++bary, e.element++) {
#if !defined(AKANTU_NDEBUG)
Point<spatial_dimension> pt1(*bary);
#endif
SpatialGrid<Element>::CellID cell_id = grid.getCellID(*bary);
SpatialGrid<Element>::neighbor_cells_iterator first_neigh_cell =
grid.beginNeighborCells(cell_id);
SpatialGrid<Element>::neighbor_cells_iterator last_neigh_cell =
grid.endNeighborCells(cell_id);
// loop over neighbors cells of the one containing the current element
for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) {
SpatialGrid<Element>::Cell::const_iterator first_neigh_el =
grid.beginCell(*first_neigh_cell);
SpatialGrid<Element>::Cell::const_iterator last_neigh_el =
grid.endCell(*first_neigh_cell);
// loop over the quadrature point in the current cell of the cell list
for (;first_neigh_el != last_neigh_el; ++first_neigh_el){
const Element & elem = *first_neigh_el;
Array<Real>::const_iterator< Vector<Real> > neigh_it =
barycenter(elem.type, elem.ghost_type).begin(sp);
const Vector<Real> & neigh_bary = neigh_it[elem.element];
Real distance = bary->distance(neigh_bary);
if(distance <= radius) {
#if !defined(AKANTU_NDEBUG)
Point<spatial_dimension> pt2(neigh_bary);
neighbors_map[pt1].push_back(pt2);
#endif
std::pair<Element, Element> pair = std::make_pair(e, elem);
pair_list::iterator p = neighbors.find(pair);
if(p != neighbors.end()) {
AKANTU_DEBUG_ERROR("Pair already registered [" << e << " " << elem << "] -> " << p->second << " " << distance);
} else {
neighbors[pair] = distance;
}
}
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Main */
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[]) {
akantu::initialize(argc, argv);
Real radius = 0.2;
Mesh mesh(spatial_dimension);
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
Int psize = comm.getNbProc();
Int prank = comm.whoAmI();
DistributedSynchronizer * dist = NULL;
if(prank == 0) {
mesh.read("bar3d.msh");
MeshPartition * partition = new MeshPartitionScotch(mesh, spatial_dimension);
partition->partitionate(psize);
dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition);
delete partition;
} else {
dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL);
}
mesh.computeBoundingBox();
const Vector<Real> & lower_bounds = mesh.getLowerBounds();
const Vector<Real> & upper_bounds = mesh.getUpperBounds();
Vector<Real> center = 0.5 * (upper_bounds + lower_bounds);
Vector<Real> spacing(spatial_dimension);
for (UInt i = 0; i < spatial_dimension; ++i) {
spacing[i] = radius * 1.2;
}
SpatialGrid<Element> grid(spatial_dimension, spacing, center);
GhostType ghost_type = _not_ghost;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
ElementTypeMapArray<Real> barycenters("", "");
mesh.initElementTypeMapArray(barycenters, spatial_dimension, spatial_dimension);
Element e;
e.ghost_type = ghost_type;
for(; it != last_type; ++it) {
UInt nb_element = mesh.getNbElement(*it, ghost_type);
e.type = *it;
Array<Real> & barycenter = barycenters(*it, ghost_type);
barycenter.resize(nb_element);
Array<Real>::iterator< Vector<Real> > bary_it = barycenter.begin(spatial_dimension);
for (UInt elem = 0; elem < nb_element; ++elem) {
mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type);
e.element = elem;
grid.insert(e, *bary_it);
++bary_it;
}
}
std::stringstream sstr; sstr << "mesh_" << prank << ".msh";
mesh.write(sstr.str());
Mesh grid_mesh(spatial_dimension, "grid_mesh", 0);
std::stringstream sstr_grid; sstr_grid << "grid_mesh_" << prank << ".msh";
grid.saveAsMesh(grid_mesh);
grid_mesh.write(sstr_grid.str());
std::cout << "Pouet 1" << std::endl;
GridSynchronizer * grid_communicator = GridSynchronizer::createGridSynchronizer(mesh, grid);
std::cout << "Pouet 2" << std::endl;
ghost_type = _ghost;
it = mesh.firstType(spatial_dimension, ghost_type);
last_type = mesh.lastType(spatial_dimension, ghost_type);
e.ghost_type = ghost_type;
for(; it != last_type; ++it) {
UInt nb_element = mesh.getNbElement(*it, ghost_type);
e.type = *it;
Array<Real> & barycenter = barycenters(*it, ghost_type);
barycenter.resize(nb_element);
Array<Real>::iterator< Vector<Real> > bary_it = barycenter.begin(spatial_dimension);
for (UInt elem = 0; elem < nb_element; ++elem) {
mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type);
e.element = elem;
grid.insert(e, *bary_it);
++bary_it;
}
}
Mesh grid_mesh_ghost(spatial_dimension, "grid_mesh_ghost", 0);
std::stringstream sstr_gridg; sstr_gridg << "grid_mesh_ghost_" << prank << ".msh";
grid.saveAsMesh(grid_mesh_ghost);
grid_mesh_ghost.write(sstr_gridg.str());
std::cout << "Pouet 3" << std::endl;
neighbors_map_t<spatial_dimension>::type neighbors_map;
pair_list neighbors;
updatePairList(barycenters, grid, radius, neighbors, neighbors_map);
pair_list::iterator nit = neighbors.begin();
pair_list::iterator nend = neighbors.end();
std::stringstream sstrp; sstrp << "pairs_" << prank;
std::ofstream fout(sstrp.str().c_str());
for(;nit != nend; ++nit) {
fout << "[" << nit->first.first << "," << nit->first.second << "] -> "
<< nit->second << std::endl;
}
std::string file = "neighbors_ref";
std::stringstream sstrf; sstrf << file << "_" << prank;
file = sstrf.str();
std::ofstream nout;
nout.open(file.c_str());
neighbors_map_t<spatial_dimension>::type::iterator it_n = neighbors_map.begin();
neighbors_map_t<spatial_dimension>::type::iterator end_n = neighbors_map.end();
for(;it_n != end_n; ++it_n) {
std::sort(it_n->second.begin(), it_n->second.end());
std::vector< Point<spatial_dimension> >::iterator it_v = it_n->second.begin();
std::vector< Point<spatial_dimension> >::iterator end_v = it_n->second.end();
nout << "####" << std::endl;
nout << it_n->second.size() << std::endl;
nout << it_n->first << std::endl;
nout << "#" << std::endl;
for(;it_v != end_v; ++it_v) {
nout << *it_v << std::endl;
}
}
fout.close();
AKANTU_DEBUG_INFO("Creating TestAccessor");
TestAccessor test_accessor(mesh, barycenters);
SynchronizerRegistry synch_registry(test_accessor);
synch_registry.registerSynchronizer(*dist, _gst_smm_mass);
synch_registry.registerSynchronizer(*grid_communicator, _gst_test);
AKANTU_DEBUG_INFO("Synchronizing tag on Dist");
synch_registry.synchronize(_gst_smm_mass);
AKANTU_DEBUG_INFO("Synchronizing tag on Grid");
synch_registry.synchronize(_gst_test);
delete grid_communicator;
delete dist;
akantu::finalize();
return EXIT_SUCCESS;
}
c4science · Help