~ubuntu-branches/ubuntu/trusty/yade/trusty

Vector3r bbMin=tuple2vec(python::extract<python::tuple>(Aabb[0])()), bbMax=tuple2vec(python::extract<python::tuple>(Aabb[1])()); Vector3r box=bbMax-bbMin;

shared_ptr<Scene> rb=Omega::instance().getScene();

Real E=0;

FOREACH(const shared_ptr<Interaction>&i, *rb->interactions){

if(!i->phys) continue;

shared_ptr<CpmPhys> bc=dynamic_pointer_cast<CpmPhys>(i->phys); if(!bc) continue;

const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb), b2=Body::byId(i->getId2(),rb);

bool isIn1=isInBB(b1->state->pos,bbMin,bbMax), isIn2=isInBB(b2->state->pos,bbMin,bbMax);

if(!isIn1 && !isIn2) continue;

Real weight=1.;

if((!isIn1 && isIn2) || (isIn1 && !isIn2)){

//shared_ptr<Body> bIn=isIn1?b1:b2, bOut=isIn2?b2:b1;

Vector3r vIn=(isIn1?b1:b2)->state->pos, vOut=(isIn2?b1:b2)->state->pos;

#define _WEIGHT_COMPONENT(axis) if(vOut[axis]<bbMin[axis]) weight=min(weight,abs((vOut[axis]-bbMin[axis])/(vOut[axis]-vIn[axis]))); else if(vOut[axis]>bbMax[axis]) weight=min(weight,abs((vOut[axis]-bbMax[axis])/(vOut[axis]-vIn[axis])));

_WEIGHT_COMPONENT(0); _WEIGHT_COMPONENT(1); _WEIGHT_COMPONENT(2);

assert(weight>=0 && weight<=1);

}

E+=.5*bc->E*bc->crossSection*pow(bc->epsN,2)+.5*bc->G*bc->crossSection*pow(bc->epsT.norm(),2);

}

return E/(box[0]*box[1]*box[2]);

}

#endif

/*! Set velocity of all dynamic particles so that if their motion were unconstrained,

* axis given by axisPoint and axisDirection would be reached in timeToAxis

* (or, point at distance subtractDist from axis would be reached).

* The code is analogous to AxialGravityEngine and is intended to give initial motion

* to particles subject to axial compaction to speed up the process. */

void velocityTowardsAxis(const Vector3r& axisPoint, const Vector3r& axisDirection, Real timeToAxis, Real subtractDist=0., Real perturbation=0.1){

FOREACH(const shared_ptr<Body>&b, *(Omega::instance().getScene()->bodies)){

if(b->state->blockedDOFs==State::DOF_ALL) continue;

const Vector3r& x0=b->state->pos;

const Vector3r& x1=axisPoint;

const Vector3r x2=axisPoint+axisDirection;

Vector3r closestAxisPoint=(x2-x1) * /* t */ (-(x1-x0).dot(x2-x1))/((x2-x1).squaredNorm());

Vector3r toAxis=closestAxisPoint-x0;

if(subtractDist>0) toAxis*=(toAxis.norm()-subtractDist)/toAxis.norm();

b->state->vel=toAxis/timeToAxis;

Vector3r ppDiff=perturbation*(1./sqrt(3.))*Vector3r(Mathr::UnitRandom(),Mathr::UnitRandom(),Mathr::UnitRandom())*b->state->vel.norm();

b->state->vel+=ppDiff;

}

BOOST_PYTHON_FUNCTION_OVERLOADS(velocityTowardsAxis_overloads,velocityTowardsAxis,3,5);

void particleConfinement(){

CpmStateUpdater csu; csu.update();

}

// makes linker error out with monolithic build..

#if 0

python::dict testNumpy(){

Scene* scene=Omega::instance().getScene().get();

int dim1[]={scene->bodies->size()};

int dim2[]={scene->bodies->size(),3};

numpy_boost<Real,1> mass(dim1);

numpy_boost<Real,2> vel(dim2);

FOREACH(const shared_ptr<Body>& b, *scene->bodies){

if(!b) continue;

mass[b->getId()]=b->state->mass;

VECTOR3R_TO_NUMPY(vel[b->getId()],b->state->vel);

}

py::dict ret;

ret["mass"]=mass;

ret["vel"]=vel;

return ret;

}

#endif

#if 0

100

/* Compute stress tensor on each particle */

101

void particleMacroStress(void){

102

Scene* scene=Omega::instance().getScene().get();

103

// find interactions of each particle

104

std::vector<std::list<Body::id_t> > bIntr(scene->bodies->size());

105

FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){

106

if(!i->isReal) continue;

107

// only contacts between 2 spheres

108

Sphere* s1=dynamic_cast<Sphere*>(Body::byId(i->getId1(),scene)->shape.get())

109

Sphere* s2=dynamic_cast<Sphere*>(Body::byId(i->getId2(),scene)->shape.get())

110

if(!s1 || !s2) continue;

111

bIntr[i.getId1()].push_back(i.getId2());

112

bIntr[i.getId2()].push_back(i.getId1());

113

}

114

for(Body::id_t id1=0; id1<(Body::id_t)bIntr.size(); id1++){

115

if(bIntr[id1].size()==0) continue;

116

Matrix3r ss(Matrix3r::Zero()); // stress tensor on current particle

117

FOREACH(Body::id_t id2, bIntr[id1]){

118

const shared_ptr<Interaction> i=scene->interactions->find(id1,id2);

119

assert(i);

120

Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(i->geom);

121

CpmPhys* phys=YADE_CAST<CpmPhys*>(i->phys);

122

Real d=(geom->se31->pos-geom->se32->pos).norm(); // current contact length

123

const Vector3r& n=geom->normal;

124

const Real& A=phys->crossSection;

125

const Vector3r& sigmaT=phys->sigmaT;

126

const Real& sigmaN=phys->sigmaN;

127

for(int i=0; i<3; i++) for(int j=0;j<3; j++){

128

ss[i][j]=d*A*(sigmaN*n[i]*n[j]+.5*(sigmaT[i]*n[j]+sigmaT[j]*n[i]));

129

}

130

}

131

// divide ss by V of the particle

132

// FIXME: for now, just use 2*(4/3)*π*r³ (.5 porosity)

133

ss*=1/(2*(4/3)*Mathr::PI*);

134

}

135

}

136

#endif

137

#include<yade/lib/smoothing/WeightedAverage2d.hpp>

138

/* Fastly locate interactions within given distance from a point in 2d (projected to plane) */

139

struct HelixInteractionLocator2d{

140

struct FlatInteraction{ Real r,h,theta; shared_ptr<Interaction> i; FlatInteraction(Real _r, Real _h, Real _theta, const shared_ptr<Interaction>& _i): r(_r), h(_h), theta(_theta), i(_i){}; };

141

shared_ptr<GridContainer<FlatInteraction> > grid;

142

Real thetaSpan;

143

int axis;

144

HelixInteractionLocator2d(Real dH_dTheta, int _axis, Real periodStart, Real theta0, Real thetaMin, Real thetaMax): axis(_axis){

145

Scene* scene=Omega::instance().getScene().get();

146

Real inf=std::numeric_limits<Real>::infinity();

147

Vector2r lo=Vector2r(inf,inf), hi(-inf,-inf);

148

Real minD0(inf),maxD0(-inf);

149

Real minTheta(inf), maxTheta(-inf);

150

std::list<FlatInteraction> intrs;

151

// first pass: find extrema for positions and interaction lengths, build interaction list

152

FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){

153

Dem3DofGeom* ge=dynamic_cast<Dem3DofGeom*>(i->geom.get());

154

CpmPhys* ph=dynamic_cast<CpmPhys*>(i->phys.get());

155

if(!ge || !ph) continue;

156

Real r,h,theta;

157

boost::tie(r,h,theta)=Shop::spiralProject(ge->contactPoint,dH_dTheta,axis,periodStart,theta0);

158

if(!isnan(thetaMin) && theta<thetaMin) continue;

159

if(!isnan(thetaMax) && theta>thetaMax) continue;

160

lo=lo.cwiseMin(Vector2r(r,h)); hi=hi.cwiseMax(Vector2r(r,h));

161

minD0=min(minD0,ge->refLength); maxD0=max(maxD0,ge->refLength);

162

minTheta=min(minTheta,theta); maxTheta=max(maxTheta,theta);

163

intrs.push_back(FlatInteraction(r,h,theta,i));

164

}

165

// create grid, put interactions inside

166

Vector2i nCells=Vector2i(max(10,(int)((hi[0]-lo[0])/(2*minD0))),max(10,(int)((hi[1]-lo[1])/(2*minD0))));

167

Vector2r hair=1e-2*Vector2r((hi[0]-lo[0])/nCells[0],(hi[1]-lo[1])/nCells[1]); // avoid rounding issue on boundary, enlarge by 1/100 cell size on each side

168

grid=shared_ptr<GridContainer<FlatInteraction> >(new GridContainer<FlatInteraction>(lo-hair,hi+hair,nCells));

169

FOREACH(const FlatInteraction& fi, intrs){

170

grid->add(fi,Vector2r(fi.r,fi.h));

171

}

172

thetaSpan=maxTheta-minTheta;

173

}

174

py::list intrsAroundPt(const Vector2r& pt, Real radius){

175

py::list ret;

176

FOREACH(const Vector2i& v, grid->circleFilter(pt,radius)){

177

FOREACH(const FlatInteraction& fi, grid->grid[v[0]][v[1]]){

178

if((pow(fi.r-pt[0],2)+pow(fi.h-pt[1],2))>radius*radius) continue; // too far

179

ret.append(fi.i);

180

}

181

}

182

return ret;

183

}

184

// return macroscopic stress around interactions that project around given point and their average omega

185

// stresses are rotated around axis back by theta angle

186

python::tuple macroAroundPt(const Vector2r& pt, Real radius){

187

Matrix3r ss(Matrix3r::Zero());

188

Real omegaCumm=0, kappaCumm=0; int nIntr=0;

189

FOREACH(const Vector2i& v, grid->circleFilter(pt,radius)){

190

FOREACH(const FlatInteraction& fi, grid->grid[v[0]][v[1]]){

191

if((pow(fi.r-pt[0],2)+pow(fi.h-pt[1],2))>radius*radius) continue; // too far

192

Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(fi.i->geom.get());

193

CpmPhys* phys=YADE_CAST<CpmPhys*>(fi.i->phys.get());

194

// transformation matrix, to rotate to the plane

195

Vector3r ax(Vector3r::Zero()); ax[axis]=1.;

196

Quaternionr q(AngleAxisr(fi.theta,ax)); q=q.conjugate();

197

Matrix3r TT=q.toRotationMatrix();

198

199

Real d=(geom->se31.position-geom->se32.position).norm(); // current contact length

200

const Vector3r& n=TT*geom->normal;

201

const Real& A=phys->crossSection;

202

const Vector3r& sigmaT=TT*phys->sigmaT;

203

const Real& sigmaN=phys->sigmaN;

204

for(int i=0; i<3; i++) for(int j=0;j<3; j++){

205

ss(i,j)+=d*A*(sigmaN*n[i]*n[j]+.5*(sigmaT[i]*n[j]+sigmaT[j]*n[i]));

206

}

207

omegaCumm+=phys->omega; kappaCumm+=phys->kappaD;

208

nIntr++;

209

}

210

}

211

// divide by approx spatial volume over which we averaged:

212

// spiral cylinder with two half-spherical caps at ends

213

ss*=1/((4/3.)*Mathr::PI*pow(radius,3)+Mathr::PI*pow(radius,2)*(thetaSpan*pt[0]-2*radius));

214

return python::make_tuple(nIntr,ss,omegaCumm/nIntr,kappaCumm/nIntr);

215

}

216

Vector2r getLo(){ return grid->getLo(); }

217

Vector2r getHi(){ return grid->getHi(); }

218

219

};

220

221

#ifdef YADE_VTK

222

223

/* Fastly locate interactions within given distance from a given point. See python docs for details. */

224

class InteractionLocator{

225

// object doing the work, see http://www.vtk.org/doc/release/5.2/html/a01048.html

226

vtkPointLocator *locator;

227

// used by locator

228

vtkUnstructuredGrid* grid;

229

vtkPoints* points;

230

// maps vtk's id to Interaction objects

231

vector<shared_ptr<Interaction> > intrs;

232

// axis-aligned bounds of all interactions

233

Vector3r mn,mx;

234

// count of interactions we hold

235

int cnt;

236

public:

237

InteractionLocator(){

238

// traverse all real interactions in the current simulation

239

// add them to points, create grid with those points

240

// store shared_ptr's to interactions in intrs separately

241

Scene* scene=Omega::instance().getScene().get();

242

locator=vtkPointLocator::New();

243

points=vtkPoints::New();

244

int count=0;

245

FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){

246

if(!i->isReal()) continue;

247

Dem3DofGeom* ge=dynamic_cast<Dem3DofGeom*>(i->geom.get());

248

if(!ge) continue;

249

const Vector3r& cp(ge->contactPoint);

250

int id=points->InsertNextPoint(cp[0],cp[1],cp[2]);

251

if(intrs.size()<=(size_t)id){intrs.resize(id+1000);}

252

intrs[id]=i;

253

count++;

254

}

255

if(count==0) throw std::runtime_error("Zero interactions when constructing InteractionLocator!");

256

double bb[6];

257

points->ComputeBounds(); points->GetBounds(bb);

258

mn=Vector3r(bb[0],bb[2],bb[4]); mx=Vector3r(bb[1],bb[3],bb[5]);

259

cnt=points->GetNumberOfPoints();

260

261

grid=vtkUnstructuredGrid::New();

262

grid->SetPoints(points);

263

locator->SetDataSet(grid);

264

locator->BuildLocator();

265

}

266

// cleanup

267

~InteractionLocator(){ locator->Delete(); points->Delete(); grid->Delete(); }

268

269

py::list intrsAroundPt(const Vector3r& pt, Real radius){

270

vtkIdList *ids=vtkIdList::New();

271

locator->FindPointsWithinRadius(radius,(const double*)(&pt),ids);

272

int numIds=ids->GetNumberOfIds();

273

py::list ret;

274

for(int i=0; i<numIds; i++){

275

// LOG_TRACE("Add "<<i<<"th id "<<ids->GetId(i));

276

ret.append(intrs[ids->GetId(i)]);

277

}

278

return ret;

279

}

280

python::tuple macroAroundPt(const Vector3r& pt, Real radius, Real forceVolume=-1){

281

Matrix3r ss(Matrix3r::Zero());

282

vtkIdList *ids=vtkIdList::New();

283

locator->FindPointsWithinRadius(radius,(const double*)(&pt),ids);

284

int numIds=ids->GetNumberOfIds();

285

Real omegaCumm=0, kappaCumm=0;

286

for(int k=0; k<numIds; k++){

287

const shared_ptr<Interaction>& I(intrs[ids->GetId(k)]);

288

Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->geom.get());

289

CpmPhys* phys=YADE_CAST<CpmPhys*>(I->phys.get());

290

Real d=(geom->se31.position-geom->se32.position).norm(); // current contact length

291

const Vector3r& n=geom->normal;

292

const Real& A=phys->crossSection;

293

const Vector3r& sigmaT=phys->sigmaT;

294

const Real& sigmaN=phys->sigmaN;

295

for(int i=0; i<3; i++) for(int j=0;j<3; j++){

296

ss(i,j)+=d*A*(sigmaN*n[i]*n[j]+.5*(sigmaT[i]*n[j]+sigmaT[j]*n[i]));

297

}

298

omegaCumm+=phys->omega; kappaCumm+=phys->kappaD;

299

}

300

Real volume=(forceVolume>0?forceVolume:(4/3.)*Mathr::PI*pow(radius,3));

301

ss*=1/volume;

302

return py::make_tuple(ss,omegaCumm/numIds,kappaCumm/numIds);

303

}

304

py::tuple getBounds(){ return py::make_tuple(mn,mx);}

305

int getCnt(){ return cnt; }

306

};

307

#endif

308

309

BOOST_PYTHON_MODULE(_eudoxos){

310

YADE_SET_DOCSTRING_OPTS;

311

py::def("velocityTowardsAxis",velocityTowardsAxis,velocityTowardsAxis_overloads(py::args("axisPoint","axisDirection","timeToAxis","subtractDist","perturbation")));

312

// def("spiralSphereStresses2d",spiralSphereStresses2d,(python::arg("dH_dTheta"),python::arg("axis")=2));

313

py::def("particleConfinement",particleConfinement);

314

#if 0

315

import_array();

316

py::def("testNumpy",testNumpy);

317

#endif

318

#ifdef YADE_VTK

319

py::class_<InteractionLocator>("InteractionLocator","Locate all (real) interactions in space by their :yref:`contact point<Dem3DofGeom::contactPoint>`. When constructed, all real interactions are spatially indexed (uses `vtkPointLocator <http://www.vtk.org/doc/release/5.4/html/a01247.html>`_ internally). Use instance methods to use those data. \n\n.. note::\n\tData might become inconsistent with real simulation state if simulation is being run between creation of this object and spatial queries.")

320

.def("intrsAroundPt",&InteractionLocator::intrsAroundPt,((python::arg("point"),python::arg("maxDist"))),"Return list of real interactions that are not further than *maxDist* from *point*.")

321

.def("macroAroundPt",&InteractionLocator::macroAroundPt,((python::arg("point"),python::arg("maxDist"),python::arg("forceVolume")=-1)),"Return tuple of averaged stress tensor (as Matrix3), average omega and average kappa values. *forceVolume* can be used (if positive) rather than the sphere (with *maxDist* radius) volume for the computation. (This is useful if *point* and *maxDist* encompass empty space that you want to avoid.)")

322

.add_property("bounds",&InteractionLocator::getBounds,"Return coordinates of lower and uppoer corner of axis-aligned abounding box of all interactions")

323

.add_property("count",&InteractionLocator::getCnt,"Number of interactions held")

324

;

325

#endif

326

py::class_<HelixInteractionLocator2d>("HelixInteractionLocator2d",

327

"Locate all real interactions in 2d plane (reduced by spiral projection from 3d, using ``Shop::spiralProject``, which is the same as :yref:`yade.utils.spiralProject`) using their :yref:`contact points<Dem3DofGeom::contactPoint>`. \n\n.. note::\n\tDo not run simulation while using this object.",

328

python::init<Real,int,Real,Real,Real,Real>((python::arg("dH_dTheta"),python::arg("axis")=0,python::arg("periodStart")=NaN,python::arg("theta0")=0,python::arg("thetaMin")=NaN,python::arg("thetaMax")=NaN),":param float dH_dTheta: Spiral inclination, i.e. height increase per 1 radian turn;\n:param int axis: axis of rotation (0=x,1=y,2=z)\n:param float theta: spiral angle at zero height (theta intercept)\n:param float thetaMin: only interactions with $\\theta$≥\\ *thetaMin* will be considered (NaN to deactivate)\n:param float thetaMax: only interactions with $\\theta$≤\\ *thetaMax* will be considered (NaN to deactivate)\n\nSee :yref:`yade.utils.spiralProject`.")

329

)

330

.def("intrsAroundPt",&HelixInteractionLocator2d::intrsAroundPt,(python::arg("pt2d"),python::arg("radius")),"Return list of interaction objects that are not further from *pt2d* than *radius* in the projection plane")

331

.def("macroAroundPt",&HelixInteractionLocator2d::macroAroundPt,(python::arg("pt2d"),python::arg("radius")),"Compute macroscopic stress around given point; the interaction ($n$ and $\\sigma^T$ are rotated to the projection plane by $\\theta$ (as given by :yref:`yade.utils.spiralProject`) first, but no skew is applied). The formula used is\n\n.. math::\n\n \\sigma_{ij}=\\frac{1}{V}\\sum_{IJ}d^{IJ}A^{IJ}\\left[\\sigma^{N,IJ}n_i^{IJ}n_j^{IJ}+\\frac{1}{2}\\left(\\sigma_i^{T,IJ}n_j^{IJ}+\\sigma_j^{T,IJ}n_i^{IJ}\\right)\\right]\n\nwhere the sum is taken over volume $V$ containing interactions $IJ$ between spheres $I$ and $J$;\n\n* $i$, $j$ indices denote Cartesian components of vectors and tensors,\n* $d^{IJ}$ is current distance between spheres $I$ and $J$,\n* $A^{IJ}$ is area of contact $IJ$,\n* $n$ is ($\\theta$-rotated) interaction normal (unit vector pointing from center of $I$ to the center of $J$)\n* $\\sigma^{N,IJ}$ is normal stress (as scalar) in contact $IJ$,\n* $\\sigma^{T,IJ}$ is shear stress in contact $IJ$ in global coordinates and $\\theta$-rotated. \n\nAdditionally, computes average of :yref:`CpmPhys.omega` ($\\bar\\omega$) and :yref:`CpmPhys.kappaD` ($\\bar\\kappa_D$). *N* is the number of interactions in the volume given.\n\n:return: tuple of (*N*, $\\tens{\\sigma}$, $\\bar\\omega$, $\\bar\\kappa_D$).\n")

332

.add_property("lo",&HelixInteractionLocator2d::getLo,"Return lower corner of the rectangle containing all interactions.")

333

.add_property("hi",&HelixInteractionLocator2d::getHi,"Return upper corner of the rectangle containing all interactions.");

334

335

}

Older »