yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10301
[Branch ~yade-pkg/yade/git-trunk] Rev 3771: Merge branch 'master' of github.com:yade/trunk
Merge authors:
Anton Gladky (gladky-anton)
jduriez (jduriez)
Raphael Maurin <raph_maurin@xxxxxxxxxxx>
------------------------------------------------------------
revno: 3771 [merge]
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Wed 2013-12-04 23:20:23 +0100
message:
Merge branch 'master' of github.com:yade/trunk
added:
examples/stl-gts/
examples/stl-gts/README
examples/stl-gts/cone.geo
examples/stl-gts/convert2stl.sh
examples/stl-gts/gts-stl.py
scripts/checks-and-tests/checks/checkViscElEng.py
modified:
examples/clumps/addToClump-example.py
pkg/dem/ViscoelasticPM.cpp
pkg/dem/ViscoelasticPM.hpp
py/_utils.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'examples/clumps/addToClump-example.py'
--- examples/clumps/addToClump-example.py 2013-08-21 06:35:10 +0000
+++ examples/clumps/addToClump-example.py 2013-12-03 07:53:08 +0000
@@ -3,6 +3,8 @@
'''This example shows usage of addToClump() and appendClumped().'''
+from yade import pack,export,qt
+
#define material for all bodies:
id_Mat=O.materials.append(FrictMat(young=1e6,poisson=0.3,density=1000,frictionAngle=1))
Mat=O.materials[id_Mat]
@@ -81,3 +83,5 @@
O.dt=1e-6
print '\nPress Play button ... '
+renderer = qt.Renderer()
+qt.View()
=== added directory 'examples/stl-gts'
=== added file 'examples/stl-gts/README'
--- examples/stl-gts/README 1970-01-01 00:00:00 +0000
+++ examples/stl-gts/README 2013-12-03 07:47:09 +0000
@@ -0,0 +1,7 @@
+This example shows, how to perform the following conversion:
+GEO->STL->GTS->Yade(clumps)
+
+Then GTS-file (mesh) is importing into Yade-Script and new spheres
+are adding as a clump within the volume of this mesh.
+
+The script requires gmsh and libgts-bin installed.
=== added file 'examples/stl-gts/cone.geo'
--- examples/stl-gts/cone.geo 1970-01-01 00:00:00 +0000
+++ examples/stl-gts/cone.geo 2013-12-03 07:47:09 +0000
@@ -0,0 +1,38 @@
+acc = 20.0/1000.0;
+D = 280.0/2.0/1000.0;
+d = 180.0/2.0/1000.0;
+h = 180.0/1000.0;
+Point(1) = {0.0, 0.0, 0.0, acc};
+Point(2) = {0.0, D, 0.0, acc};
+Point(3) = {0.0, -D, 0.0, acc};
+Point(4) = {D, 0.0, 0.0, acc};
+Point(5) = {-D, 0.0, 0.0, acc};
+Point(6) = {0.0, 0.0, h, acc};
+Circle(1) = {2, 1, 4};
+Circle(2) = {4, 1, 3};
+Circle(3) = {3, 1, 5};
+Circle(4) = {5, 1, 2};
+Line(5) = {2, 1};
+Line(6) = {2, 6};
+Line(7) = {4, 1};
+Line(8) = {4, 6};
+Line(9) = {3, 1};
+Line(10) = {3, 6};
+Line(11) = {5, 1};
+Line(12) = {5, 6};
+Line Loop(13) = {1, 7, -5};
+Plane Surface(14) = {13};
+Line Loop(15) = {2, 9, -7};
+Plane Surface(16) = {15};
+Line Loop(17) = {3, 11, -9};
+Plane Surface(18) = {17};
+Line Loop(19) = {4, 5, -11};
+Plane Surface(20) = {19};
+Line Loop(21) = {3, 12, -10};
+Ruled Surface(22) = {21};
+Line Loop(23) = {4, 6, -12};
+Ruled Surface(24) = {23};
+Line Loop(25) = {1, 8, -6};
+Ruled Surface(26) = {25};
+Line Loop(27) = {2, 10, -8};
+Ruled Surface(28) = {27};
=== added file 'examples/stl-gts/convert2stl.sh'
--- examples/stl-gts/convert2stl.sh 1970-01-01 00:00:00 +0000
+++ examples/stl-gts/convert2stl.sh 2013-12-03 07:47:09 +0000
@@ -0,0 +1,3 @@
+#!/bin/bash
+gmsh -2 cone.geo -o cone.stl
+stl2gts -r < cone.stl > cone.gts
=== added file 'examples/stl-gts/gts-stl.py'
--- examples/stl-gts/gts-stl.py 1970-01-01 00:00:00 +0000
+++ examples/stl-gts/gts-stl.py 2013-12-03 07:47:09 +0000
@@ -0,0 +1,42 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+# © 2009 Václav Šmilauer <eudoxos@xxxxxxxx>
+# © 2013 Anton Gladky <gladk@xxxxxxxxxx>
+
+from yade import pack
+import gts, os.path, locale
+
+surf=gts.read(open('cone.gts'))
+
+if surf.is_closed():
+ pred=pack.inGtsSurface(surf)
+ aabb=pred.aabb()
+ dim0=aabb[1][0]-aabb[0][0]; radius=dim0/70. # get some characteristic dimension, use it for radius
+ O.bodies.appendClumped(pack.regularHexa(pred,radius=radius,gap=radius/4.))
+ surf.translate(0,-(aabb[1][1]-aabb[0][1])/2.0,-(aabb[1][2]-aabb[0][2])) # move surface down so that facets are underneath the falling spheres
+O.bodies.append(pack.gtsSurface2Facets(surf,wire=True))
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],label='collider'),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
+ [Ip2_FrictMat_FrictMat_FrictPhys()],
+ [Law2_L3Geom_FrictPhys_ElPerfPl()],
+ ),
+ NewtonIntegrator(damping=.1,gravity=[0,0,-500.0]),
+ PyRunner(iterPeriod=1000,command='timing.stats(); O.pause();'),
+ PyRunner(iterPeriod=10,command='addPlotData()')
+]
+O.dt=.7*PWaveTimeStep()
+O.saveTmp()
+O.timingEnabled=True
+O.trackEnergy=True
+from yade import plot
+plot.plots={'i':('total',O.energy.keys,)}
+def addPlotData(): plot.addData(i=O.iter,total=O.energy.total(),**O.energy)
+plot.plot(subPlots=False)
+
+from yade import timing
+from yade import qt
+qt.View()
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp 2013-08-28 13:31:12 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2013-12-04 19:25:29 +0000
@@ -14,8 +14,22 @@
/* ViscElPhys */
ViscElPhys::~ViscElPhys(){}
+/* Contact parameter calculation function */
+Real Ip2_ViscElMat_ViscElMat_ViscElPhys::contactParameterCalculation(const Real& l1, const Real& l2, const bool& massMultiply){
+ if (massMultiply) {
+ // If one of paramaters > 0. we DO NOT return 0
+ Real a = (l1?1/l1:0) + (l2?1/l2:0);
+ if (a) return 1/a;
+ else return 0;
+ } else {
+ // If one of paramaters > 0, we return 0
+ return (l1>0 or l2>0)?l1*l2/(l1+l2):0;
+ }
+}
+
/* Ip2_ViscElMat_ViscElMat_ViscElPhys */
void Ip2_ViscElMat_ViscElMat_ViscElPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
+
// no updates of an existing contact
if(interaction->phys) return;
ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
@@ -36,27 +50,18 @@
const Real ks2 = mat2->ks*mass2; const Real cs2 = mat2->cs*mass2;
ViscElPhys* phys = new ViscElPhys();
-
- if ((kn1>0) or (kn2>0)) {
- phys->kn = 1/( ((kn1>0)?1/kn1:0) + ((kn2>0)?1/kn2:0) );
- } else {
- phys->kn = 0;
- }
- if ((ks1>0) or (ks2>0)) {
- phys->ks = 1/( ((ks1>0)?1/ks1:0) + ((ks2>0)?1/ks2:0) );
- } else {
- phys->ks = 0;
- }
-
- if ((mR1>0) or (mR2>0)) {
+
+ phys->kn = contactParameterCalculation(kn1,kn2, mat1->massMultiply&&mat2->massMultiply);
+ phys->ks = contactParameterCalculation(ks1,ks2, mat1->massMultiply&&mat2->massMultiply);
+ phys->cn = contactParameterCalculation(cn1,cn2, mat1->massMultiply&&mat2->massMultiply);
+ phys->cs = contactParameterCalculation(cs1,cs2, mat1->massMultiply&&mat2->massMultiply);
+
+ if ((mR1>0) or (mR2>0)) {
phys->mR = 2.0/( ((mR1>0)?1/mR1:0) + ((mR2>0)?1/mR2:0) );
} else {
phys->mR = 0;
}
-
- phys->cn = (cn1?1/cn1:0) + (cn2?1/cn2:0); phys->cn = phys->cn?1/phys->cn:0;
- phys->cs = (cs1?1/cs1:0) + (cs2?1/cs2:0); phys->cs = phys->cs?1/phys->cs:0;
-
+
phys->tangensOfFrictionAngle = std::tan(std::min(mat1->frictionAngle, mat2->frictionAngle));
phys->shearForce = Vector3r(0,0,0);
@@ -209,3 +214,4 @@
addTorque(id2, c2x.cross(f)-momentResistance,scene);
}
}
+
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2013-08-28 12:00:42 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2013-12-04 19:25:29 +0000
@@ -65,6 +65,8 @@
virtual void go(const shared_ptr<Material>& b1,
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
+ private :
+ Real contactParameterCalculation(const Real& l1,const Real& l2, const bool& massMultiply);
YADE_CLASS_BASE_DOC(Ip2_ViscElMat_ViscElMat_ViscElPhys,IPhysFunctor,"Convert 2 instances of :yref:`ViscElMat` to :yref:`ViscElPhys` using the rule of consecutive connection.");
FUNCTOR2D(ViscElMat,ViscElMat);
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2013-11-13 07:37:51 +0000
+++ py/_utils.cpp 2013-11-28 11:05:20 +0000
@@ -537,7 +537,7 @@
py::def("normalShearStressTensors",Shop__normalShearStressTensors,(py::args("compressionPositive")=false,py::args("splitNormalTensor")=false,py::args("thresholdForce")=NaN),"Compute overall stress tensor of the periodic cell decomposed in 2 parts, one contributed by normal forces, the other by shear forces. The formulation can be found in [Thornton2000]_, eq. (3):\n\n.. math:: \\tens{\\sigma}_{ij}=\\frac{2}{V}\\sum R N \\vec{n}_i \\vec{n}_j+\\frac{2}{V}\\sum R T \\vec{n}_i\\vec{t}_j\n\nwhere $V$ is the cell volume, $R$ is \"contact radius\" (in our implementation, current distance between particle centroids), $\\vec{n}$ is the normal vector, $\\vec{t}$ is a vector perpendicular to $\\vec{n}$, $N$ and $T$ are norms of normal and shear forces.\n\n:param bool splitNormalTensor: if true the function returns normal stress tensor split into two parts according to the two subnetworks of strong an weak forces.\n\n:param Real thresholdForce: threshold value according to which the normal stress tensor can be split (e.g. a zero value would make distinction between tensile and compressive forces).");
py::def("fabricTensor",Shop__fabricTensor,(py::args("splitTensor")=false,py::args("revertSign")=false,py::args("thresholdForce")=NaN),"Compute the fabric tensor of the periodic cell. The original paper can be found in [Satake1982]_.\n\n:param bool splitTensor: split the fabric tensor into two parts related to the strong and weak contact forces respectively.\n\n:param bool revertSign: it must be set to true if the contact law's convention takes compressive forces as positive.\n\n:param Real thresholdForce: if the fabric tensor is split into two parts, a threshold value can be specified otherwise the mean contact force is considered by default. It is worth to note that this value has a sign and the user needs to set it according to the convention adopted for the contact law. To note that this value could be set to zero if one wanted to make distinction between compressive and tensile forces.");
py::def("bodyStressTensors",Shop__getStressLWForEachBody,"Compute and return a table with per-particle stress tensors. Each tensor represents the average stress in one particle, obtained from the contour integral of applied load as detailed below. This definition is considering each sphere as a continuum. It can be considered exact in the context of spheres at static equilibrium, interacting at contact points with negligible volume changes of the solid phase (this last assumption is not restricting possible deformations and volume changes at the packing scale).\n\nProof:\n\nFirst, we remark the identity: $\\sigma_{ij}=\\delta_{ik}\\sigma_{kj}=x_{i,k}\\sigma_{kj}=(x_{i}\\sigma_{kj})_{,k}-x_{i}\\sigma_{kj,k}$.\n\nAt equilibrium, the divergence of stress is null: $\\sigma_{kj,k}=\\vec{0}$. Consequently, after divergence theorem: $\\frac{1}{V}\\int_V \\sigma_{ij}dV = \\frac{1}{V}\\int_V (x_{i}\\sigma_{kj})_{,k}dV = \\frac{1}{V}\\int_{\\partial V}x_i\\sigma_{kj}n_kdS = \\frac{1}{V}\\sum_bx_i^bf_j^b$.\n\nThe last equality is implicitely based on the representation of external loads as Dirac distributions whose zeros are the so-called *contact points*: 0-sized surfaces on which the *contact forces* are applied, located at $x_i$ in the deformed configuration.\n\nA weighted average of per-body stresses will give the average stress inside the solid phase. There is a simple relation between the stress inside the solid phase and the stress in an equivalent continuum in the absence of fluid pressure. For porosity $n$, the relation reads: $\\sigma_{ij}^{equ.}=(1-n)\\sigma_{ij}^{solid}$.");
- py::def("getStress",Shop::getStress,(py::args("volume")=0),"Compute and return Love-Weber stress tensor:\n\n $\\sigma_{ij}=\\frac{1}{V}\\sum_b l_i^b f_j^b$, where the sum is over all interactions, with $l$ the branch vector (joining centers of the bodies) and $f$ is the contact force. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases.");
+ py::def("getStress",Shop::getStress,(py::args("volume")=0),"Compute and return Love-Weber stress tensor:\n\n $\\sigma_{ij}=\\frac{1}{V}\\sum_b f_i^b l_j^b$, where the sum is over all interactions, with $f$ the contact force and $l$ the branch vector (joining centers of the bodies). Stress is negativ for repulsive contact forces, i.e. compression. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases.");
py::def("getCapillaryStress",Shop::getCapillaryStress,(py::args("volume")=0),"Compute and return Love-Weber capillary stress tensor:\n\n $\\sigma^{cap}_{ij}=\\frac{1}{V}\\sum_b l_i^b f^{cap,b}_j$, where the sum is over all interactions, with $l$ the branch vector (joining centers of the bodies) and $f^{cap}$ is the capillary force. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases. Only the CapillaryPhys interaction type is supported presently.");
py::def("getBodyIdsContacts",Shop__getBodyIdsContacts,(py::args("bodyID")=0),"Get a list of body-ids, which contacts the given body.");
py::def("maxOverlapRatio",maxOverlapRatio,"Return maximum overlap ration in interactions (with :yref:`ScGeom`) of two :yref:`spheres<Sphere>`. The ratio is computed as $\\frac{u_N}{2(r_1 r_2)/r_1+r_2}$, where $u_N$ is the current overlap distance and $r_1$, $r_2$ are radii of the two spheres in contact.");
=== added file 'scripts/checks-and-tests/checks/checkViscElEng.py'
--- scripts/checks-and-tests/checks/checkViscElEng.py 1970-01-01 00:00:00 +0000
+++ scripts/checks-and-tests/checks/checkViscElEng.py 2013-12-04 15:43:55 +0000
@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+# 2 spheres. One is fixed, another one failing from up.
+# Calculate en (coefficient of restitution) and compare with precalculated value
+# Coefficient of restitution is the ratio of speeds after and before an impact
+# This check-simulation checks the correctness of ViscoElasticEngine
+
+from yade import utils, plot
+o = Omega()
+fr = 0.5;rho=2000
+tc = 0.001; en = 0.7; et = 0.7; o.dt = 0.0002*tc
+
+
+r1 = 0.003
+r2 = 0.002
+
+
+param = getViscoelasticFromSpheresInteraction(tc,en,et)
+
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr,density=rho,**param))
+
+
+id1 = O.bodies.append(sphere(center=[0,0,0],radius=r1,material=mat1,fixed=True))
+id2 = O.bodies.append(sphere(center=[0,0,(r1+r2*2.0)],radius=r2,material=mat1,fixed=False))
+
+
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=(r1+r2)*5.0),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom()],
+ [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+ [Law2_ScGeom_ViscElPhys_Basic()],
+ ),
+ NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
+ PyRunner(command='addPlotData()',iterPeriod=100),
+]
+
+v0 = 0
+en = 0
+tolerance = 0.0001
+
+def addPlotData():
+ global v0, en
+ v = O.bodies[id2].state.vel[2]
+
+ if v0<=0 and v>0:
+ en=-v/v0
+ print ("Precalculated en value %f" % 0.736356797441)
+ print ("Obtained en value %f" % en)
+ O.pause()
+ v0=v
+
+O.run(1000000)
+O.wait()
+
+if ((abs(0.736356797441-en)/en)>tolerance):
+ resultStatus += 1
