yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00695
[svn] r1534 - in trunk: gui/py pkg/dem/PreProcessor
Author: eudoxos
Date: 2008-09-30 20:38:52 +0200 (Tue, 30 Sep 2008)
New Revision: 1534
Modified:
trunk/gui/py/_utils.cpp
trunk/gui/py/eudoxos.py
trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
trunk/pkg/dem/PreProcessor/TriaxialTest.hpp
Log:
1. TriaxialTest now takes readiusMean argument (by default negative, i.e. disabled); if >0, box size is scaled so that both requested porosity and mean size can be preserved.
2. Attempt to fix computation of stress from stored elastic energy.
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2008-09-30 10:49:22 UTC (rev 1533)
+++ trunk/gui/py/_utils.cpp 2008-09-30 18:38:52 UTC (rev 1534)
@@ -17,20 +17,20 @@
bool isInBB(Vector3r p, Vector3r bbMin, Vector3r bbMax){return p[0]>bbMin[0] && p[0]<bbMax[0] && p[1]>bbMin[1] && p[1]<bbMax[1] && p[2]>bbMin[2] && p[2]<bbMax[2];}
/* \todo implement groupMask */
-python::tuple aabbExtrema(Real cutoff=0.0){
+python::tuple aabbExtrema(Real cutoff=0.0, bool centers=false){
if(cutoff<0. || cutoff>1.) throw invalid_argument("Cutoff must be >=0 and <=1.");
Real inf=std::numeric_limits<Real>::infinity();
Vector3r minimum(inf,inf,inf),maximum(-inf,-inf,-inf);
FOREACH(const shared_ptr<Body>& b, *Omega::instance().getRootBody()->bodies){
shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->geometricalModel); if(!s) continue;
Vector3r rrr(s->radius,s->radius,s->radius);
- minimum=componentMinVector(minimum,b->physicalParameters->se3.position-rrr);
- maximum=componentMaxVector(maximum,b->physicalParameters->se3.position+rrr);
+ minimum=componentMinVector(minimum,b->physicalParameters->se3.position-(centers?0:rrr));
+ maximum=componentMaxVector(maximum,b->physicalParameters->se3.position+(centers?0:rrr));
}
Vector3r dim=maximum-minimum;
return python::make_tuple(vec2tuple(minimum+.5*cutoff*dim),vec2tuple(maximum-.5*cutoff*dim));
}
-BOOST_PYTHON_FUNCTION_OVERLOADS(aabbExtrema_overloads,aabbExtrema,0,1);
+BOOST_PYTHON_FUNCTION_OVERLOADS(aabbExtrema_overloads,aabbExtrema,0,2);
python::tuple negPosExtremeIds(int axis, Real distFactor=1.1){
python::tuple extrema=aabbExtrema();
@@ -163,7 +163,7 @@
BOOST_PYTHON_MODULE(_utils){
def("PWaveTimeStep",PWaveTimeStep);
- def("aabbExtrema",aabbExtrema,aabbExtrema_overloads(args("cutoff")));
+ def("aabbExtrema",aabbExtrema,aabbExtrema_overloads(args("cutoff","centers")));
def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor")));
def("coordsAndDisplacements",coordsAndDisplacements,coordsAndDisplacements_overloads(args("AABB")));
def("setRefSe3",setRefSe3);
Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py 2008-09-30 10:49:22 UTC (rev 1533)
+++ trunk/gui/py/eudoxos.py 2008-09-30 18:38:52 UTC (rev 1534)
@@ -9,7 +9,7 @@
# E=(1/2)σεAl # global stored energy
# σ=EE/(.5εAl)=EE/(.5εV)
from yade import utils
- dim=utils.aabbDim(cutoff)
+ dim=utils.aabbDim(cutoff,centers=True)
return utils.elasticEnergy(utils.aabbExtrema(cutoff))/(.5*strain*dim[0]*dim[1]*dim[2])
def estimatePoissonYoung(principalAxis,stress=0,plot=False,cutoff=0.):
Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TriaxialTest.cpp 2008-09-30 10:49:22 UTC (rev 1533)
+++ trunk/pkg/dem/PreProcessor/TriaxialTest.cpp 2008-09-30 18:38:52 UTC (rev 1534)
@@ -81,15 +81,12 @@
#include <boost/random/normal_distribution.hpp>
+CREATE_LOGGER(TriaxialTest);
-
using namespace boost;
using namespace std;
-typedef pair<Vector3r, Real> BasicSphere;
-//! generate a list of non-overlapping spheres
-string GenerateCloud(vector<BasicSphere>& sphere_list, Vector3r lowerCorner, Vector3r upperCorner, long number, Real rad_std_dev, Real porosity);
TriaxialTest::TriaxialTest () : FileGenerator()
@@ -161,6 +158,7 @@
biaxial2dTest=false;
radiusStdDev=0.3;
+ radiusMean=-1; // no radius specified
// wall_top_id =0;
// wall_bottom_id =0;
@@ -192,6 +190,7 @@
REGISTER_ATTRIBUTE(maxMultiplier);
REGISTER_ATTRIBUTE(finalMaxMultiplier);
REGISTER_ATTRIBUTE(radiusStdDev);
+ REGISTER_ATTRIBUTE(radiusMean);
REGISTER_ATTRIBUTE(sphereYoungModulus);
REGISTER_ATTRIBUTE(spherePoissonRatio);
@@ -277,6 +276,25 @@
rootBody->bodies = shared_ptr<BodyContainer>(new BodyRedirectionVector);
shared_ptr<Body> body;
+
+
+
+ /* if _mean_radius is not given (i.e. <=0), then calculate it from box size;
+ * OTOH, if it is specified, scale the box preserving its ratio and lowerCorner so that the radius can be as requested
+ */
+ Vector3r dimensions=upperCorner-lowerCorner; Real volume=dimensions.X()*dimensions.Y()*dimensions.Z();
+ Real porosity=.75;
+ Real really_radiusMean;
+
+ if(radiusMean<=0) really_radiusMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*numberOfGrains),1/3.);
+ else {
+ Real boxScaleFactor=radiusMean*pow((4/3.)*Mathr::PI*numberOfGrains/(volume*(1-porosity)),1/3.);
+ LOG_INFO("Mean radius value of "<<radiusMean<<" requested, scaling box by "<<boxScaleFactor);
+ dimensions*=boxScaleFactor;
+ upperCorner=lowerCorner+dimensions;
+ really_radiusMean=radiusMean;
+ }
+
@@ -386,27 +404,18 @@
vector<BasicSphere> sphere_list;
if(importFilename!="") sphere_list=Shop::loadSpheresFromFile(importFilename,lowerCorner,upperCorner);
- else message+=GenerateCloud(sphere_list, lowerCorner, upperCorner, numberOfGrains, radiusStdDev, 0.75);
-
+ else message+=GenerateCloud(sphere_list, lowerCorner, upperCorner, numberOfGrains, radiusStdDev, really_radiusMean, porosity);
vector<BasicSphere>::iterator it = sphere_list.begin();
vector<BasicSphere>::iterator it_end = sphere_list.end();
-
- for (;it!=it_end; ++it)
- {
- cerr << "sphere (" << it->first << " " << it->second << ")"<<endl;
- createSphere(body,it->first,it->second,false,true);
+ FOREACH(const BasicSphere& it, sphere_list){
+ LOG_DEBUG("sphere (" << it.first << " " << it.second << ")");
+ createSphere(body,it.first,it.second,false,true);
rootBody->bodies->insert(body);
}
+
return true;
-// return "Generated a sample inside box of dimensions: ("
-// + lexical_cast<string>(lowerCorner[0]) + ","
-// + lexical_cast<string>(lowerCorner[1]) + ","
-// + lexical_cast<string>(lowerCorner[2]) + ") and ("
-// + lexical_cast<string>(upperCorner[0]) + ","
-// + lexical_cast<string>(upperCorner[1]) + ","
-// + lexical_cast<string>(upperCorner[2]) + ").";
}
@@ -609,10 +618,12 @@
//cerr << "fin de section triaxialcompressionEngine = shared_ptr<TriaxialCompressionEngine> (new TriaxialCompressionEngine);" << std::endl;
// recording global stress
- triaxialStateRecorder = shared_ptr<TriaxialStateRecorder>(new TriaxialStateRecorder);
- triaxialStateRecorder-> outputFile = WallStressRecordFile + Key;
- triaxialStateRecorder-> interval = recordIntervalIter;
- //triaxialStateRecorderer-> thickness = thickness;
+ if(recordIntervalIter>0){
+ triaxialStateRecorder = shared_ptr<TriaxialStateRecorder>(new TriaxialStateRecorder);
+ triaxialStateRecorder-> outputFile = WallStressRecordFile + Key;
+ triaxialStateRecorder-> interval = recordIntervalIter;
+ //triaxialStateRecorderer-> thickness = thickness;
+ }
shared_ptr<MakeItFlat> makeItFlat(new MakeItFlat);
makeItFlat->direction=2;
@@ -699,7 +710,7 @@
}
-string GenerateCloud(vector<BasicSphere>& sphere_list, Vector3r lowerCorner, Vector3r upperCorner, long number, Real rad_std_dev, Real porosity)
+string TriaxialTest::GenerateCloud(vector<BasicSphere>& sphere_list, Vector3r lowerCorner, Vector3r upperCorner, long number, Real rad_std_dev, Real mean_radius, Real porosity)
{
typedef boost::minstd_rand StdGenerator;
static StdGenerator generator;
@@ -711,11 +722,8 @@
sphere_list.clear();
long tries = 1000; //nb of tries for positionning the next sphere
Vector3r dimensions = upperCorner - lowerCorner;
-
- Real mean_radius = std::pow(dimensions.X()*dimensions.Y()*dimensions.Z()*(1-porosity)/(3.1416*1.3333*number),0.333333);
- //cerr << mean_radius;
- std::cerr << "generating aggregates ... ";
+ LOG_INFO("Generating aggregates ...");
long t, i;
for (i=0; i<number; ++i) {
Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.hpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TriaxialTest.hpp 2008-09-30 10:49:22 UTC (rev 1533)
+++ trunk/pkg/dem/PreProcessor/TriaxialTest.hpp 2008-09-30 18:38:52 UTC (rev 1534)
@@ -73,9 +73,9 @@
StabilityCriterion,
maxMultiplier, ///max multiplier of diameters during internal compaction
finalMaxMultiplier,
- wallOversizeFactor,
- radiusStdDev; // make walls bigger (/smaller) than necessary by this factor
-
+ wallOversizeFactor, // make walls bigger (/smaller) than necessary by this factor
+ radiusStdDev,
+ radiusMean;
bool wall_top
,wall_bottom
,wall_1
@@ -127,6 +127,11 @@
void createSphere(shared_ptr<Body>& body, Vector3r position, Real radius,bool big,bool dynamic);
void createActors(shared_ptr<MetaBody>& rootBody);
void positionRootBody(shared_ptr<MetaBody>& rootBody);
+
+ typedef pair<Vector3r, Real> BasicSphere;
+ //! generate a list of non-overlapping spheres
+ string GenerateCloud(vector<BasicSphere>& sphere_list, Vector3r lowerCorner, Vector3r upperCorner, long number, Real rad_std_dev, Real mean_radius, Real porosity);
+
public :
TriaxialTest ();
@@ -137,6 +142,7 @@
void registerAttributes();
REGISTER_CLASS_NAME(TriaxialTest);
REGISTER_BASE_CLASS_NAME(FileGenerator);
+ DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(TriaxialTest,false);
