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);
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
