Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Jan Stránský posted a new comment: Just for future reference, this is exactly why I insisted on providing the code. I think you can leave the distribution of the points, but change the condition a bit, something like: %%% not tested, just an idea ia = 0.5*(i+ic); % average i ja = 0.5*(j+jc); % average j d = pi/18; % or something like that A(hang,5)=length(find(PHI(:)>ia-d&PHI(:)ja-d&Theta(:)https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Needs information => Solved jsonscript confirmed that the question is solved: Hello Jan, >Is this what you are observing? Yes, your explanation clarifies my query. In my code, I just tried to extend the 2D Rose diagram idea to spherical distribution plot. I haven't thought about the area distribution factor. It makes sense now. I have been looking at the wrong distribution plot. So, a sphere divided with equal areas, should give a more uniform distribution. I would try potting that. Thank you so much. This answers my question. -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Open => Needs information Jan Stránský requested more information: Thanks for the code > for i=pi/18:2*(pi/36):(2*pi) > for j=pi/18:(pi/18):(pi) > A(hang,5)=length(find(PHI(:)>ic&PHI(:)jc&Theta(:) %Distribution to different bins I don't think uniform separate distribution of theta and phi are correct here. Such division has higher "area" around equator and lower "area" near poles (just have a look at the "globus earth model" with meridians and parallels). This implies, that there is higher probability of normals to fall into equator bin than to the polar bin. I.e., you should have bins with more normals near equator and less normals near poles. Is this what you are observing? cheers Jan -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed:
https://answers.launchpad.net/yade/+question/691560
Status: Answered => Open
jsonscript is still having a problem:
Hello Jan,
Thank you so much for your advice. I am trying to understand the method
you mentioned for evaluation of isotropy.
Meanwhile I wish to send you my matlab code that I used to plot the
contact-normal distribution.
Matlab code:
#---#
clear all
clc
[id1,id2,nx,ny,nz,fn] = textread('normals.txt','%f %f %f %f %f %f');
c=zeros(2*length(nx),3);
for i=1:length(nx)
c(i,1)=nx(i);
c(i,2)=ny(i);
c(i,3)=nz(i);
end
%Also considering the opposite direction of every contact normal
for i=(length(nx)+1):(2*length(nx))
c(i,1)=-nx(i-length(nx));
c(i,2)=-ny(i-length(nx));
c(i,3)=-nz(i-length(nx));
end
[R,Theta,PHI]=cart2sph_luxp2(c); %Cartesian corrdiantes to polar
coordinates
hang=1;
ic=0;
jc=0;
A=zeros(650,5);
for i=pi/18:2*(pi/36):(2*pi)
for j=pi/18:(pi/18):(pi)
A(hang,5)=length(find(PHI(:)>ic&PHI(:)jc&Theta(:)3.14
jc=0;
else
jc=j;
end
end
ic=i;
end
M=max(A(:,5));
polar(nan,M*0.905);
hold all;
M
sum=sum(A(:,5))
%20:705
for k=1:648
theta1=A(k,3);
theta2=A(k,4);
phi1=A(k,1);
phi2=A(k,2);
r=A(k,5);
x1=r*sin(theta1)*cos(phi1);
y1=r*sin(theta1)*sin(phi1);
z1=r*cos(theta1);
x2=r*sin(theta2)*cos(phi1);
y2=r*sin(theta2)*sin(phi1);
z2=r*cos(theta2);
x3=r*sin(theta2)*cos(phi2);
y3=r*sin(theta2)*sin(phi2);
z3=r*cos(theta2);
x4=r*sin(theta1)*cos(phi2);
y4=r*sin(theta1)*sin(phi2);
z4=r*cos(theta1);
X1=[x1,x2,0]; %Faces
Y1=[y1,y2,0];
Z1=[z1,z2,0];
X2=[x2,x3,0];
Y2=[y2,y3,0];
Z2=[z2,z3,0];
X3=[x3,x4,0];
Y3=[y3,y4,0];
Z3=[z3,z4,0];
X4=[x4,x1,0];
Y4=[y4,y1,0];
Z4=[z4,z1,0];
X5=[x1,x2,x3,x4];
Y5=[y1,y2,y3,y4];
Z5=[z1,z2,z3,z4];
a1=r/M;
if a1>1
a1=1;
end
C=[a1, a1, a1]; %Colour according to Number
ColorSpec = a1*[1,0.5,0.5];
fill3(X1,Y1,Z1,ColorSpec); %Plane-1,2,3,4 and Top Face
fill3(X2,Y2,Z2,ColorSpec);
fill3(X3,Y3,Z3,ColorSpec);
fill3(X4,Y4,Z4,ColorSpec);
fill3(X5,Y5,Z5,ColorSpec);
end
[Xs,Ys,Zs] = sphere(18);
hs = surf(M*Xs,M*Ys,M*Zs);
set(hs,'facecolor','none','edgecolor','b','edgealpha',0.2)
set(gca,{'xtick' 'ytick' 'ztick' 'vis' 'clim'},{[] [] [] 'on' [0 M]});
axis equal vis3d;
view(3);
drawnow;
function [R,Theta,PHI]=cart2sph_luxp2(ShapeData)
% This function is used to transport Cartisian coordinates To Spherical
% coordinates
% According the notes in my PHD report.
% input a n-by-3 matrix including the points associating x,y,z
% Edit By LUXP
% 2010-11-30
x=ShapeData(:,1);
y=ShapeData(:,2);
z=ShapeData(:,3);
R=sqrt(x.^2+y.^2+z.^2);
Theta=acos(z./R);
PHI=acos(x./sqrt(x.^2+y.^2));
Temp=PHI(y<0);
Temp=2*pi-Temp;
PHI(y<0)=Temp;
end
#--#
The plot of distribution actually confused me. The plot is not a
homogeneous vector distribution in the contact normals. Vertical
direction seems to have less contact normals distributed.
Do you think there is a mistake in my way of visualisation (matlab
code)? I checked it several times and am convinced its correct.
Your suggestions will be valuable. Thank you.
--
You received this question notification because your team yade-users is
an answer contact for Yade.
___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed:
https://answers.launchpad.net/yade/+question/691560
Status: Open => Answered
Jan Stránský proposed the following answer:
Hello,
thanks for the code. I have tried it and the results are IMHO nicely isotropic.
Evaluation:
I chose 200 random unit vectors and for these vectors I counted normals that
are inclined from the vector maximally by given angle.
>From these counts I computed relative standard deviation, whose value is:
1.1 % for angle 30 deg
3.2 % for angle 15 deg
11 % for angle 5 deg
Plotting the values gives relatively nice semi-sphere.
Of course it does not give the same number for every direction, but I did not
find any systematic problem ("there is little contact in the vertical
direction, and the middle is recessed")
###
# Extract normals from file
with open("normals.txt") as f:
lines = f.readlines()
lines = [l.split()[2:] for l in lines]
lines = [[float(v) for v in l] for l in lines]
normals = [Vector3(*l) for l in lines]
print("Number of normals:",len(normals))
# Function to test isotropy.
# For random direction counts normals "close" to this direction
def countSimilarNormals(limAngle=pi/6):
# random unit vector with positive z coordinate
from random import gauss
v = [gauss(0, 1) for i in range(3)]
mag = sqrt(sum(pow(x,2) for x in v))
v = Vector3(*[x/mag for x in v])
if v[2] < 0: v *= -1
# counting "similar" normals
ret = 0
for n in normals:
# compute angle
dot = v.dot(n)
dot = abs(dot)
angle = acos(dot)
# if it is within the limit, count it
if angle <= limAngle:
ret += 1
return v,ret
# count similar normals for a few random directions
data = [countSimilarNormals(pi/6/2/3) for i in range(200)]
directions,counts = zip(*data)
# some statistics
cmax = max(counts)
import statistics
cmean = statistics.mean(counts)
cstddev = statistics.stdev(counts)
print("Isotropy fulfilled with",cstddev/cmean*100,"% relative standard
deviation")
# plot
dvals = [d*c/cmax for (d,c) in zip(directions,counts)]
x,y,z = [[d[i] for d in dvals] for i in (0,1,2)]
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes(projection="3d")
ax.auto_scale_xyz((-1,1),(-1,1),(0,1))
ax.scatter3D(x,y,z)
plt.show()
# sorry, I don't know how to set aspect ratio, so the result is displayed as
semi-ellipsoid
###
> I plotted the spherical histogram using Matlab.
> Am I doing something wrong here?
please provide the code. As discussed, this stage may contain a mistake.
> I couldn't attach my histogram plot here
yes, thanks for respecting the rule, but you can add the data (i.e. a
few lines of e.g. angle-number pairs) or the code producing the plot or
the data.
cheers
Jan
--
You received this question notification because your team yade-users is
an answer contact for Yade.
___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed:
https://answers.launchpad.net/yade/+question/691560
Status: Answered => Open
jsonscript is still having a problem:
Hello Jan,
I realise I wasn't clear enough in my previous comments. I apologise for
that.
Here is the MWE script (Runtime 2-3 mins) :
#---#
from yade import pack
from yade import plot
O.periodic=True
O.cell.setBox(.3,.3,.3)
sp=pack.SpherePack()
radius=5e-3
num=sp.makeCloud(minCorner=(0,0,0),maxCorner=(.3,.3,.3),rMean=radius,rRelFuzz=0.5,num=1,periodic=True)
O.bodies.append([sphere(s[0],s[1]) for s in sp])
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
PeriTriaxController(dynCell=True,maxUnbalanced=2e-3,relStressTol=2e-4,goal=[-1e5,-1e5,-1e5],stressMask=7,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
NewtonIntegrator(damping=.2)
]
def triaxDone():
print ('Here we are:
stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff)
for i in O.interactions:
n=i.geom.normal
f=open("normals.txt","a")
print ("%s %s %0.9f %0.9f %0.9f" % (i.id1, i.id2, n[0], n[1], n[2]), file=f)
f.close()
O.pause()
O.dt=PWaveTimeStep()
O.run()
O.wait()
#--#
I tried to consolidate a specimen isotropically with 100kPa and later checked
for the isotropy in contact-normal distribution.
>From the contact-normals recorded in the text file "normals.txt", I plotted
>the spherical histogram using Matlab. It seems to me, the distribution is not
>even near to isotropic.
I couldn't attach my histogram plot here, respecting "No-external link"
rule of this forum. But, I see the distribution is not isotropic.
Am I doing something wrong here? if not, is there some other preparation
method to generate isotropic sample w.r.t. both stress and contact
normals?
--
You received this question notification because your team yade-users is
an answer contact for Yade.
___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Jan Stránský proposed the following answer: > I just want to know how to prepare uniform samples to finish triaxial test.. Out of mny options, see [2] figure 2.3 created by [3]. But as I said, it is possible (and very likely IMO) that your packing already is isotropic and there is a mistake in your evaluation. cheers Jan [2] https://github.com/stranskyjan/phd-thesis/blob/master/pdfs/PhD_thesis_Stransky_2018_v2.pdf [3] https://github.com/stranskyjan/phd-thesis/blob/master/codes/scripts/randomPeriPack/packingIsotropy.sh -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Open => Answered Jan Stránský proposed the following answer: Yes, PeriTriaxController cannot guarantee a contact normal isotropy. Which strongly depends on "contact normal isotropy" definition. The contacts are random and discrete with finite number of contacts, so the isotropy in the strict sense is simply impossible. On the other hand, the contact normal distribution should be "very close" to isotropy (depending on definition and evaluation). Mechanical properties should be "very close" to isotropy, too. If it is not like that in your case and you want a help, please provide a MWE. cheers Jan -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Needs information => Open jsonscript gave more information on the question: Hello, Peritriaxcontroller we use is to produce isotropic stress state.. but they cannot guarantee us that we would get a contact normal isotropy. Thanks! -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Open => Needs information Jan Stránský requested more information: Hello, please read [1] and provide a MWE illustrating your result and how you measured "contact normal". Is it a problem? Of course it is possible to give you some code to prepare uniform samples, but it can be very different from your needs and therefore useless. It is also possible that your sample is already isotropic and you just have some mistake in the isotropy-evaluation. So a MWE provided by you is the most reasonable way to start a reasonable discussion.. cheers Jan [1] https://www.yade-dem.org/wiki/Howtoask -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Needs information => Open jsonscript gave more information on the question: hello. Sorry, I just want to know how to prepare uniform samples to finish triaxial test.. Because my sample after consolidation is always anisotropic.. Thanks! -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Open => Needs information Jan Stránský requested more information: > some codes are listed here: instant replay: please read [1] and provide a MWE* illustrating result and how you measured "contact normal". To not-just-guess (out of mny possible reasons), we need your code. cheers Jan *M = minimal *W = working (complete is a necessary condition) [1] https://www.yade-dem.org/wiki/Howtoask -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Open => Needs information Jan Stránský requested more information: Hello, please provide a MWE [1] illustrating result, how you measured "contact normal". To not-just-guess (out of mny possible reasons), we need your code. cheers Jan [1] https://www.yade-dem.org/wiki/Howtoask -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Needs information => Open jsonscript gave more information on the question: Hello， Triaxial compression test with sand.. I meaurse the contact normal by i.geom.normal in O.interactions..then draw the rose diagram.. The distribution diagram is not a sphere, there is little contact in the vertical direction, and the middle is recessed... Here is part of my code: PeriTriaxController(dynCell=True,maxUnbalanced=2e-3,relStressTol=2e-5,goal=[-1e5,-1e5,-1e5],stressMask=7,globUpdate=5,maxStrainRate=[0.1,0.1,0.1],doneHook='triaxDone()',label='triax'), Thanks! -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
Re: [Yade-users] [Question #691560]: contact normal is not spherical
Question #691560 on Yade changed: https://answers.launchpad.net/yade/+question/691560 Status: Open => Needs information Jan Stránský requested more information: Hello, please provide a MWE [1] illustrating result, how you measured "contact normal" and what exactly is "spherical". cheers Jan [1] https://www.yade-dem.org/wiki/Howtoask -- You received this question notification because your team yade-users is an answer contact for Yade. ___ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
