Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/examples/cluster_statistics.py
17e997c1486afeatures/netcdf-mpi

Last Change
Raw File

cluster_statistics.py
View Options

#!/usr/bin/python
# -*- coding: utf-8 -*-
##
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
#
# @section LICENSE
#
# Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
# Solides)
#
# Tamaas 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.
#
# Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
#
#
################################################################
import numpy as np
################################################################
def plotSurface(surf):
fig = plt.figure()
axe = fig.add_subplot(111)
img = axe.imshow(surf)
cbar = fig.colorbar(img)
################################################################
# surface generation
################################################################
from tamaas import *
n = 128
SG = SurfaceGeneratorFilterFFT()
SG.getGridSize().assign(n)
SG.getHurst().assign(0.8)
SG.getRMS().assign(1.);
SG.getQ0().assign(4);
SG.getQ1().assign(4);
SG.getQ2().assign(64);
SG.getRandomSeed().assign(20);
SG.Init()
s = SG.buildSurface()
rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(s)
s *= 1./rms_slopes_spectral
s -= s.max()
################################################################
# surface plot
################################################################
import matplotlib.pyplot as plt
#fig1 = plotSurface(s)
################################################################
# contact contact
################################################################
bem = BemPolonski(s)
bem.setEffectiveModulus(1.)
load = 0.1
bem.computeEquilibrium(1e-13,load)
tractions = bem.getTractions()
displacements = bem.getDisplacements()
#fig2 = plotSurface(tractions)
#fig3 = plotSurface(displacements)
################################################################
# cluster statistics
################################################################
area = ContactArea(tractions.shape[0],1.);
area.getSurface()[tractions > 0.] = 1
# Cluster detection
print "==============================="
print "Detect contact clusters: "
area.detectContactClusters()
print "Found " + str(area.getNumClusters())
clusters = ContactClusterCollection(area);
cluster_list = area.getContactClusters();
cluster_areas = [float(c.getA())/float(n)**2 for c in cluster_list]
cluster_perimeters = [c.getP() for c in cluster_list]
print "cluster_areas",cluster_areas
print "cluster_perimeters",cluster_perimeters
print "cluster_total_area",clusters.getTotalArea()
print "cluster_total_perimeter",clusters.getTotalPerimeter()
print "nb_clusters_with_hole",clusters.getNbClustersWithHole()
print "nb_clusters",clusters.getNbClusters()
################################################################
value,bins = np.histogram(cluster_areas,bins=50)
fig = plt.figure()
axe = fig.add_subplot(111)
axe.loglog(bins[:-1],value,'-o')
################################################################
plt.show()
c4science · Help