← Back to branch summary

~ubuntu-branches/debian/jessie/yade/jessie

« back to all changes in this revision

Viewing changes to pkg/dem/PeriIsoCompressor.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2014-08-04 19:34:58 UTC
  • mfrom: (1.1.11)
  • Revision ID: package-import@ubuntu.com-20140804193458-cw8qhnujxe9wzi15
Tags: 1.11.0-1
* [a0600ae] Imported Upstream version 1.11.0
* [a3055e0] Do not use parallel build on kfreebsd-amd64 and s390x.
* [f86b405] Remove applied patches.

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/dem/PeriIsoCompressor.cpp
8
8
#include<yade/pkg/dem/DemXDofGeom.hpp>
9
9
#include<yade/lib/pyutil/gil.hpp>
10
10
 
11
 
using namespace std;
12
 
 
13
11
YADE_PLUGIN((PeriIsoCompressor)(PeriTriaxController)(Peri3dController))
14
12
 
15
13
CREATE_LOGGER(PeriIsoCompressor);
50
48
                Real sigAvg=(sigma[0]+sigma[1]+sigma[2])/3., avgArea=(cellArea[0]+cellArea[1]+cellArea[2])/3., avgSize=(cellSize[0]+cellSize[1]+cellSize[2])/3.;
51
49
                Real avgGrow=1e-4*(sigmaGoal-sigAvg)*avgArea/(avgStiffness>0?avgStiffness:1);
52
50
                Real maxToAvg=maxSize/avgSize;
53
 
                if(abs(maxToAvg*avgGrow)>maxDisplPerStep) avgGrow=Mathr::Sign(avgGrow)*maxDisplPerStep/maxToAvg;
 
51
                if(std::abs(maxToAvg*avgGrow)>maxDisplPerStep) avgGrow=Mathr::Sign(avgGrow)*maxDisplPerStep/maxToAvg;
54
52
                Real okGrow=-(minSize-2.1*maxSpan)/maxToAvg;
55
53
                if(avgGrow<okGrow) throw runtime_error("Unable to shring cell due to maximum body size (although required by stress condition). Increase particle rigidity, increase total sample dimensions, or decrease goal stress.");
56
54
                // avgGrow=max(avgGrow,-(minSize-2.1*maxSpan)/maxToAvg);
57
55
                if(avgStiffness>0) { sigma+=(avgGrow*avgStiffness)*Vector3r::Ones(); sigAvg+=avgGrow*avgStiffness; }
58
 
                if(abs((sigAvg-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
 
56
                if(std::abs((sigAvg-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
59
57
                cellGrow=(avgGrow/avgSize)*cellSize;
60
58
        }
61
59
        else{ // handle each dimension separately
65
63
                        // FIXME: that is why the fixup 1e-4 is needed here
66
64
                        // FIXME: or perhaps maxDisplaPerStep=1e-2*charLen is too big??
67
65
                        cellGrow[axis]=1e-4*(sigmaGoal-sigma[axis])*cellArea[axis]/(avgStiffness>0?avgStiffness:1);  // FIXME: wrong dimensions? See PeriTriaxController
68
 
                        if(abs(cellGrow[axis])>maxDisplPerStep) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxDisplPerStep;
 
66
                        if(std::abs(cellGrow[axis])>maxDisplPerStep) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxDisplPerStep;
69
67
                        cellGrow[axis]=max(cellGrow[axis],-(cellSize[axis]-2.1*maxSpan));
70
68
                        // crude way of predicting sigma, for steps when it is not computed from intrs
71
69
                        if(avgStiffness>0) sigma[axis]+=cellGrow[axis]*avgStiffness; // FIXME: dimensions
72
 
                        if(abs((sigma[axis]-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
 
70
                        if(std::abs((sigma[axis]-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
73
71
                }
74
72
        }
75
73
        TRVAR4(cellGrow,sigma,sigmaGoal,avgStiffness);
115
113
                Vector3r branch=Body::byId(I->getId2(),scene)->state->pos + scene->cell->hSize*I->cellDist.cast<Real>() -Body::byId(I->getId1(),scene)->state->pos;
116
114
                stressTensor+=f*branch.transpose();
117
115
                if( !dynCell ){
118
 
                        for ( int i=0; i<3; i++ ) sumStiff[i]+=abs ( gsc->normal[i] ) *nsi->kn+ ( 1-abs ( gsc->normal[i] ) ) *nsi->ks;
 
116
                        for ( int i=0; i<3; i++ ) sumStiff[i]+=std::abs ( gsc->normal[i] ) *nsi->kn+ ( 1-std::abs ( gsc->normal[i] ) ) *nsi->ks;
119
117
                        n++;}
120
118
        }
121
119
        // Divide by volume as in stressTensor=sum(fi*lj)/Volume (Love equation)
184
182
                // steady evolution with fluctuations; see TriaxialStressController
185
183
                if (!dynCell) strain_rate=(1-growDamping)*strain_rate+.8*prevGrow[axis];
186
184
                // limit maximum strain rate
187
 
                if (abs(strain_rate)>maxStrainRate[axis]) strain_rate = Mathr::Sign(strain_rate)*maxStrainRate[axis];
 
185
                if (std::abs(strain_rate)>maxStrainRate[axis]) strain_rate = Mathr::Sign(strain_rate)*maxStrainRate[axis];
188
186
                // do not shrink below minimum cell size (periodic collider condition), although it is suboptimal WRT resulting stress
189
187
                strain_rate=max(strain_rate,-(cellSize[axis]-2.1*maxBodySpan[axis])/scene->dt);
190
188
 
195
193
                // signal if condition not satisfied
196
194
                if(stressMask&(1<<axis)){
197
195
                        Real curr=stress[axis];
198
 
                        if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>relStressTol) || abs(curr-goal[axis])>absStressTol) allOk=false;
 
196
                        if((goal[axis]!=0 && std::abs((curr-goal[axis])/goal[axis])>relStressTol) || std::abs(curr-goal[axis])>absStressTol) allOk=false;
199
197
                }else{
200
198
                        Real curr=strain[axis];
201
199
                        // since strain is prescribed exactly, tolerances need just to accomodate rounding issues
202
 
                        if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>1e-6) || abs(curr-goal[axis])>1e-6){
 
200
                        if((goal[axis]!=0 && std::abs((curr-goal[axis])/goal[axis])>1e-6) || std::abs(curr-goal[axis])>1e-6){
203
201
                                allOk=false;
204
 
                                if(doUpdate) LOG_DEBUG("Strain not OK; "<<abs(curr-goal[axis])<<">1e-6");}
 
202
                                if(doUpdate) LOG_DEBUG("Strain not OK; "<<std::abs(curr-goal[axis])<<">1e-6");}
205
203
                }
206
204
        }
207
205
        // update stress and strain
286
284
                                // convert relative progress values of ##Path to absolute values
287
285
                                PATH_OP_OP(i,j,0) *= 1./PATH_OP_OP(i,pathSizes[i]-1,0);
288
286
                                // convert relative stress/strain values of ##Path to absolute stress strain values
289
 
                                if (abs(PATH_OP_OP(i,pathSizes[i]-1,1)) >= 1e-9) { // the final value is not 0 (otherwise always absolute values are considered)
 
287
                                if (std::abs(PATH_OP_OP(i,pathSizes[i]-1,1)) >= 1e-9) { // the final value is not 0 (otherwise always absolute values are considered)
290
288
                                        PATH_OP_OP(i,j,1) *= goal(i)/PATH_OP_OP(i,pathSizes[i]-1,1);
291
289
                                }
292
290
                        }
387
385
                }
388
386
        }
389
387
 
390
 
        // correction coefficient ix strainRate.maxabs() > maxStrainRate
 
388
        // correction coefficient ix strainRate.maxstd::abs() > maxStrainRate
391
389
        Real srCorr = (strainRate.cwiseAbs().maxCoeff() > maxStrainRate)? (maxStrainRate/strainRate.cwiseAbs().maxCoeff()):1.;
392
390
        strainRate *= srCorr;
393
391