Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/contact_cluster.cpp
e96d4bebeabcfeatures/netcdf-mpi

Last Change
Raw File

contact_cluster.cpp
View Options

//==================================================================//
// Percolation detection and characterization of contact clusters //
// To be integrated in FFT-BEM code //
// //
// V.A. Yastrebov, 2014 //
// Centre des Materiaux, CNRS, MINES ParisTech //
// //
// G. Anciaux, 2014 //
// EPFL, ENAC, IIC, LSMS //
//==================================================================//
#include "contact_cluster.hh"
#include <iostream>
#include <cmath>
#include "contact_area.hh"
#include "cluster_grow.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
void ContactCluster::computeBoundingBox(ContactArea & contact_area) {
int N = contact_area.size();
int i0 = -1;
int j0 = -1;
int i1 = -1;
int j1 = -1;
for (int i=0; i<N; i++) {
for (int j=0; j<N; j++) {
if (contact_area(i,j) == this->id && i0 < 0) i0 = i;
if (contact_area(j,i) == this->id && j0 < 0) j0 = i;
if (contact_area(N-i-1,N-j-1) == this->id && i1 < 0) i1 = N-i-1;
if (contact_area(N-j-1,N-i-1) == this->id && j1 < 0) j1 = N-i-1;
// if all extremum are found then break
if (i0 >= 0 && j0 >= 0 && i1 >= 0 && j1 >= 0) break;
}
// if all extremum are found then break
if (i0 >= 0 && j0 >= 0 && i1 >= 0 && j1 >= 0) break;
}
if (i0 < 0 || j0 < 0) SURFACE_FATAL("bounds not found " << id << " " << contact_area(0,0));
if (i1 < 0 || j1 < 0) SURFACE_FATAL("bounds not found");
lb_corner[0] = i0;
lb_corner[1] = j0;
rt_corner[0] = i1;
rt_corner[1] = j1;
bbox_size[0] = i1-i0+1;
bbox_size[1] = j1-j0+1;
}
/* -------------------------------------------------------------------------- */
void ContactCluster::checkLines(UInt dir,ContactArea & contact_area){
UInt N = contact_area.size();
UInt i0 = lb_corner[dir];
start[dir] = -1;
stop[dir] = -1;
for (UInt i = 0; i < N; i++) {
bool fail = false;
for (UInt j = 0; j < N; j++){
int val = -1;
if (dir == 0) val = contact_area(i,j);
if (dir == 1) val = contact_area(j,i);
if (val == id) {
fail = true;
break;
}
}
if (i0 == 0 && !fail && start[dir] < 0) {
start[dir] = i;
if (stop[dir] >= 0) SURFACE_FATAL("Error occured at "<<__LINE__);
}
if (i0 == 0 && fail && start[dir] >= 0 && stop[dir] < 0)
stop[dir] = i;
}
}
/* -------------------------------------------------------------------------- */
void ContactCluster::minimizeBoundingBox(ContactArea & contact_area) {
this->split_cluster = false;
UInt N = contact_area.size();
if (bbox_size[0] != N && bbox_size[1] != N) return;
// Minimize bounding box
this->split_cluster = true;
// Check if there are vertical lines that do not contact "id"
this->checkLines(0,contact_area);
// Check if there are horizontal lines that do not contact "id"
this->checkLines(1,contact_area);
if (start[0] >= 0) bbox_size[0] -= (stop[0]-start[0])-1;
if (start[1] >= 0) bbox_size[1] -= (stop[1]-start[1])-1;
}
/* -------------------------------------------------------------------------- */
void ContactCluster::extractCluster(ContactArea & contact_area) {
UInt Nx = this->size(0);
UInt Ny = this->size(1);
UInt i0 = lb_corner[0];
UInt j0 = lb_corner[1];
UInt N = contact_area.size();
mass_center[0] = 0;
mass_center[1] = 0;
for (UInt i = 0; i < Nx; i++) {
for (UInt j = 0; j < Ny; j++) {
if (i+i0 >= N || j+j0 >= N)
SURFACE_FATAL(" i0=" << i0
<< ",j0=" << j0
<< ", i=" << i
<< ", j=" << j);
int val = contact_area(i+i0,j+j0);
if (val == id) {
this->at(i,j) = 1;
mass_center[0] += i+i0;
mass_center[1] += j+j0;
}
}
}
}
/* -------------------------------------------------------------------------- */
void ContactCluster::extractClusterPeriodic(ContactArea & contact_area) {
// if the cluster is split by periodicity
mass_center[0] = 0;
mass_center[1] = 0;
UInt Nx = this->size(0);
UInt Ny = this->size(1);
UInt N = contact_area.size();
UInt i0 = lb_corner[0];
UInt j0 = lb_corner[1];
UInt Nxmax = start[0]<0?Nx:N;
UInt Nymax = start[1]<0?Ny:N;
for (UInt i = 0; i < Nxmax; i++) {
for (UInt j = 0; j < Nymax; j++) {
if (i+i0>=0 && i+i0<N &&
j+j0>=0 && j+j0<N) {
int val = contact_area(i+i0,j+j0);
if (val == id) {
if (start[0] >= 0 && start[1] < 0) {
mass_center[1] += j+j0;
if ((int)i < start[0]) {
this->at(i,j) = 1;
mass_center[0]+=i;
}
else {
this->at(i-(stop[0]-start[0])+1,j) = 1;
mass_center[0] -= (N-i);
}
}
else if (start[1] >= 0 && start[0] < 0) {
mass_center[0]+=i+i0;
if ((int)j < start[1]) {
this->at(i,j) = 1;
mass_center[1]+=j;
}
else {
this->at(i,j-(stop[1]-start[1])+1) = 1;
mass_center[1] -= (N-j);
}
}
else if (start[1] >= 0 && start[0] >= 0) {
if ((int)j < start[1] && (int)i < start[0]) {
this->at(i,j) = 1;
mass_center[0] += i;
mass_center[1] += j;
}
else if ((int)j < start[1] && (int)i >= start[0]) {
this->at(i-(stop[0]-start[0])+1,j) = 1;
mass_center[1]+=j;
mass_center[0]-=(N-i);
}
else if ((int)j >= start[1] && (int)i >= start[0]) {
this->at(i-(stop[0]-start[0])+1,j-(stop[1]-start[1])+1) = 1;
mass_center[0] -= (N-i);
mass_center[1] -= (N-j);
}
else if ((int)j >= start[1] && (int)i < start[0]) {
this->at(i,j-(stop[1]-start[1])+1) = 1;
mass_center[0] += i;
mass_center[1] -= (N-j);
}
}
else if (start[0] < 0 && start[1] < 0) {
this->at(i,j) = 1;
mass_center[0] += i+i0;
mass_center[1]+=j+j0;
}
}
}
else {
std::cout<<"Error i0="<<i0<<",j0="<<j0
<<",i="<<i<<",j="<<j
<<",Nxmax="<<Nxmax << ",Nymax" << Nymax
<<std::endl;
abort();
}
}
}
}
/* -------------------------------------------------------------------------- */
ContactCluster::ContactCluster(int color, ContactArea & contact_area)
: Map2d(0,1.) {
concave = false;
id = 0;
lb_corner[0] = 0.;
lb_corner[1] = 0.;
rt_corner[0] = 0.;
rt_corner[1] = 0.;
start[0] = 0.;
stop[1] = 0.;
bbox_size[0] = 0.;
bbox_size[1] = 0.;
A = 0.;
P = 0.;
IP = 0.;
mass_center[0] = 0;
mass_center[1] = 0;
Ps = 0.;
Cd = 0.;
//RI = 0.;
comp = 0.;
split_cluster = 0.;
this->id = color;
this->computeBoundingBox(contact_area);
if (bbox_size[0] < 0 || bbox_size[1] < 0) {
SURFACE_FATAL("Something wrong with the detection of the bounding box of a contact cluster");
}
this->minimizeBoundingBox(contact_area);
this->setGridSize(bbox_size);
if (!split_cluster) this->extractCluster(contact_area);
else this->extractClusterPeriodic(contact_area);
}
/* -------------------------------------------------------------------------- */
void ContactCluster::computeStatistics() {
bool start = false;
// UInt start_i, start_j;
A=0;
P=0;
IP=0;
// UInt Nmax = bbox_size[0] > bbox_size[1] ? bbox_size[0]:bbox_size[1];
UInt intersection[2] = {0,0};
for (UInt i = 0; i < bbox_size[0]; i++) {
intersection[0] = 0;
intersection[1] = 0;
for (UInt j = 0; j < bbox_size[1]; j++) {
// Area
A += this->at(i,j);
// Perimeter at edges
if ( (i==0 || i== bbox_size[0] - 1) && this->at(i,j) == 1) {
if (!start) {
start=true;
// start_i=i;
//start_j=j;
}
P++;
}
if ( (j==0 || j== bbox_size[0] - 1) &&
this->at(i,j) == 1) {
if (!start) {
start=true;
//start_i=i;
//start_j=j;
}
P++;
}
// Perimeter along lines
if (i>0 && this->at(i,j) != this->at(i-1,j)) {
if (!start) {
start=true;
//start_i=i;
//start_j=j;
}
P++;
}
// Perimeter along columns
if (j>0 && this->at(i,j) != this->at(i,j-1)) {
P++;
}
// Interpixel perimeter along lines
if (i>0 && this->at(i,j) == this->at(i-1,j) && this->at(i,j) > 0)
IP++;
// Interpixel perimeter along columns
if (j>0 && this->at(i,j) == this->at(i,j-1) && this->at(i,j) > 0)
IP++;
}
}
for (UInt k = 0; k < 2; ++k) {
mass_center[k] /= Real(A);
//if (mass_center[k] < 0) {
// mass_center[k] += N;
//}
}
this->comp = P/sqrt(Real(A));
Cd = (A-P/4.)/(A-sqrt(1.*A));
if (intersection[0] > 2 || intersection[1] > 2)
concave=true;
this->findHoles();
}
/* -------------------------------------------------------------------------- */
void ContactCluster::findHoles() {
UInt Nx = this->size(0);
UInt Ny = this->size(1);
ClusterGrow clg(*this,id,-1,0,id);
for (UInt i0 = 0; i0 < Nx; ++i0) {
for (UInt j0 = 0; j0 < Ny; ++j0) {
// check if the point touch the border and is not in contact
if ( (i0 == 0 || j0 == 0 ||
i0 == Nx-1 || j0 == Ny-1) &&
this->at(i0,j0) == 0) {
clg.grow<true>(i0,j0);
}
}
}
// find the holes
std::map<UInt,Real> holes_areas;
for (UInt i=0; i<Nx; i++)
for (UInt j=0; j<Ny; j++)
if (this->at(i,j) > id){
if (holes_areas.count(this->at(i,j)) == 0)
holes_areas[this->at(i,j)] = 0;
holes_areas[this->at(i,j)] += 1;
}
// this->RI = (Nx*Ny-Ai-A)/Real(A);
}
/* -------------------------------------------------------------------------- */
void ContactCluster::printself(std::ostream & stream, int indent) const {
stream << mass_center[0] << " " << mass_center[1] << " "
<< bbox_size[0] << " " << bbox_size[1] << " "
<< sqrt(A) << " "
<< P << " "
<< comp << " "
<< IP << " "
<< Cd << " "
<< this->getNbHoles() << " "
// << RI << " "
<< id;
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help