Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78365021
contact_cluster.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Aug 20, 05:53
Size
9 KB
Mime Type
text/x-c
Expires
Thu, Aug 22, 05:53 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20033894
Attached To
rTAMAAS tamaas
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"
/* -------------------------------------------------------------------------- */
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;
}
/* -------------------------------------------------------------------------- */
Event Timeline
Log In to Comment