← Back to branch summary

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

« back to all changes in this revision

Viewing changes to pkg/pfv/FlowEngine.ipp.in

  • 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/pfv/FlowEngine.ipp.in
275
275
        buffer.clear();
276
276
        buffer.resize(scene->bodies->size());
277
277
        shared_ptr<Sphere> sph ( new Sphere );
278
 
        const int Sph_Index = sph->getClassIndexStatic();
 
278
        const int Sph_Index = sph->getClassIndexStatic();
279
279
        FOREACH ( const shared_ptr<Body>& b, *scene->bodies ) {
280
 
                if (!b || ignoredBody==b->getId()) continue;
281
 
                posData& dat = buffer[b->getId()];
 
280
                if (!b || ignoredBody==b->getId() || b->isClumpMember()) continue;
 
281
                posData& dat = buffer[b->getId()];
282
282
                dat.id=b->getId();
283
283
                dat.pos=b->state->pos;
284
284
                dat.isSphere= (b->shape->getClassIndex() ==  Sph_Index);
 
285
                dat.isClump = b->isClump();
285
286
                if (dat.isSphere) dat.radius = YADE_CAST<Sphere*>(b->shape.get())->radius;
 
287
                if (dat.isClump) {
 
288
                        const shared_ptr<Clump>& clump = YADE_PTR_CAST<Clump>(b->shape);
 
289
                        const shared_ptr<Body>& member = Body::byId(clump->members.begin()->first,scene);
 
290
                        dat.radius = pow( (3*b->state->mass)/(4*Mathr::PI*member->material->density) , 1.0/3.0); //use equivalent radius of clump (just valid for nearly spherical clumps)
 
291
                }
286
292
                dat.exists=true;
287
293
        }
288
294
}
293
299
        solver->xMin = Mathr::MAX_REAL, solver->xMax = -Mathr::MAX_REAL, solver->yMin = Mathr::MAX_REAL, solver->yMax = -Mathr::MAX_REAL, solver->zMin = Mathr::MAX_REAL, solver->zMax = -Mathr::MAX_REAL;
294
300
        FOREACH ( const posData& b, buffer ) {
295
301
                if ( !b.exists ) continue;
296
 
                if ( b.isSphere ) {
 
302
                if ( b.isSphere || b.isClump ) {
297
303
                        const Real& rad = b.radius;
298
304
                        const Real& x = b.pos[0];
299
305
                        const Real& y = b.pos[1];
358
364
///Using one-by-one insertion
359
365
        vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
360
366
        FOREACH ( const posData& b, buffer ) {
361
 
                if ( !b.exists ) continue;
362
 
                if ( b.isSphere ) {
363
 
                        if (b.id==ignoredBody) continue;
364
 
                        flow.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
365
 
        }
 
367
                if ( !b.exists || b.id==ignoredBody ) continue;
 
368
                if ( b.isSphere || b.isClump ) flow.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );
 
369
        }
366
370
        flow.tesselation().redirected=true;//By inserting one-by-one, we already redirected
367
371
        flow.shearLubricationForces.resize ( flow.tesselation().maxId+1 );
368
372
        flow.shearLubricationTorques.resize ( flow.tesselation().maxId+1 );
391
395
                        case ( 3 ) : cell->info().volume() = volumeCellTripleFictious ( cell ); break;
392
396
                        default: break; 
393
397
                }
394
 
                if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow.volumeSolidPore(cell) ); }
 
398
                if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( std::abs(cell->info().volume()) - flow.volumeSolidPore(cell) ); }
395
399
        }
396
400
        if (debug) cout << "Volumes initialised." << endl;
397
401
}
486
490
                }
487
491
        }
488
492
        Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
489
 
        return abs ( Volume );
 
493
        return std::abs ( Volume );
490
494
}
491
495
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
492
496
template<class Cellhandle>
521
525
        Real Vol1=0.5* ( ( A-BS ).cross ( B-BS ) ) [coord[1]]* ( 0.333333333* ( 2*B[coord[1]]+A[coord[1]] )-Wall_coordinate[1] );
522
526
        //second pyramid with triangular base (A,AS,BS)
523
527
        Real Vol2=0.5* ( ( AS-BS ).cross ( A-BS ) ) [coord[1]]* ( 0.333333333* ( B[coord[1]]+2*A[coord[1]] )-Wall_coordinate[1] );
524
 
        return abs ( Vol1+Vol2 );
 
528
        return std::abs ( Vol1+Vol2 );
525
529
}
526
530
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
527
531
template<class Cellhandle>
549
553
                }
550
554
        }
551
555
        Real Volume = ( A[coord[0]]-Wall_coordinate[0] ) * ( A[coord[1]]-Wall_coordinate[1] ) * ( A[coord[2]]-Wall_coordinate[2] );
552
 
        return abs ( Volume );
 
556
        return std::abs ( Volume );
553
557
}
554
558
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
555
559
template<class Cellhandle>
620
624
                                        bool v1fictious = Tes.vertex ( id1 )->info().isFictious;
621
625
                                        int bnd = v1fictious? id1 : id2;
622
626
                                        int coord = flow.boundary(bnd).coordinate;
623
 
                                        O1O2 = v1fictious ? abs((sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]) : abs((sph1->state->pos + makeVector3r(Tes.vertex(id1)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]);
 
627
                                        O1O2 = v1fictious ? std::abs((sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]) : std::abs((sph1->state->pos + makeVector3r(Tes.vertex(id1)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]);
624
628
                                        if(v1fictious)
625
629
                                                normal = makeVector3r(flow.boundary(id1).normal);
626
630
                                        else
725
729
//                 if ( !cell->info().isReal() ) continue;
726
730
                if ( cell->info().isGhost ) continue;
727
731
                for ( int i=0;i<3;i++ )
728
 
                        meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
729
 
                volume+=abs ( cell->info().volume() );
 
732
                        meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * std::abs ( cell->info().volume() ) );
 
733
                volume+=std::abs ( cell->info().volume() );
730
734
        }
731
735
        return ( meanVel/volume );
732
736
}