← Back to branch summary

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

« back to all changes in this revision

Viewing changes to pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2014-10-20 21:31:37 UTC
  • mfrom: (1.1.13)
  • Revision ID: package-import@ubuntu.com-20141020213137-wt0doalgtw493ohd
Tags: 1.12.0-1
http://bugs.debian.org/763259
* [3bba065] Imported Upstream version 1.12.0. (Closes: #763259)
* [1ac4c00] Remove patch applied by upstream.
* [0b5f802] Set Standards-Version: 3.9.6. No changes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp
12
12
//FIXME : needs "requestErase" somewhere
13
13
 
14
14
#include "Law2_ScGeom_CapillaryPhys_Capillarity.hpp"
15
 
#include <yade/pkg/common/ElastMat.hpp>
16
 
#include <yade/pkg/dem/ScGeom.hpp>
 
15
#include <pkg/common/ElastMat.hpp>
 
16
#include <pkg/dem/ScGeom.hpp>
17
17
 
18
 
#include <yade/pkg/dem/CapillaryPhys.hpp>
19
 
#include <yade/pkg/dem/HertzMindlin.hpp>
20
 
#include <yade/core/Omega.hpp>
21
 
#include <yade/core/Scene.hpp>
22
 
#include <yade/lib/base/Math.hpp>
 
18
#include <pkg/dem/CapillaryPhys.hpp>
 
19
#include <pkg/dem/HertzMindlin.hpp>
 
20
#include <core/Omega.hpp>
 
21
#include <core/Scene.hpp>
 
22
#include <lib/base/Math.hpp>
23
23
 
24
24
YADE_PLUGIN((Law2_ScGeom_CapillaryPhys_Capillarity));
25
25
 
65
65
        shared_ptr<BodyContainer>& bodies = scene->bodies;
66
66
        if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene);
67
67
 
68
 
        InteractionContainer::iterator ii    = scene->interactions->begin();
69
 
        InteractionContainer::iterator iiEnd = scene->interactions->end();
70
68
        bool hertzInitialized = false;
71
 
        for (; ii!=iiEnd ; ++ii) {
72
 
 
 
69
        #ifdef YADE_OPENMP
 
70
        const long size=scene->interactions->size();
 
71
        #pragma omp parallel for schedule(guided) num_threads(ompThreads>0 ? min(ompThreads,omp_get_max_threads()) : omp_get_max_threads())
 
72
        for(long i=0; i<size; i++){
 
73
                const shared_ptr<Interaction>& interaction=(*scene->interactions)[i];
 
74
        #else
 
75
        FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){
 
76
        #endif
73
77
                /// interaction is real
74
 
                if ((*ii)->isReal()) {
75
 
                        const shared_ptr<Interaction>& interaction = *ii;
 
78
                if (interaction->isReal()) {
76
79
                        if (!hertzInitialized) {//NOTE: We are assuming that only one type is used in one simulation here
77
80
                                if (CapillaryPhys::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=false;
78
81
                                else if (MindlinCapillaryPhys::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=true;
116
119
                        if ((currentContactGeometry->penetrationDepth>=0)|| D<=0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere
117
120
                                D=0; // defines fCap when spheres interpenetrate. D<0 leads to wrong interpolation has D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
118
121
                                if (!hertzOn) {
119
 
                                        if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
 
122
                                        if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert(interaction);
120
123
                                        cundallContactPhysics->meniscus=true;
121
124
                                } else {
122
 
                                        if (fusionDetection && !mindlinContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
 
125
                                        if (fusionDetection && !mindlinContactPhysics->meniscus) bodiesMenisciiList.insert(interaction);
123
126
                                        mindlinContactPhysics->meniscus=true;
124
127
                                }
125
128
                        }
155
158
                                        else mindlinContactPhysics->meniscus = false;
156
159
                                }
157
160
                                if (!Vinterpol) {
158
 
                                        if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove((*ii));
 
161
                                        if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove(interaction);
159
162
                                        if (D>0) scene->interactions->requestErase(interaction);
160
163
                                }
161
164
                                /// wetting angles
168
171
                                }
169
172
                        }
170
173
                ///interaction is not real      //If the interaction is not real, it should not be in the list
171
 
                } else if (fusionDetection) bodiesMenisciiList.remove((*ii));
 
174
                } else if (fusionDetection) bodiesMenisciiList.remove(interaction);
172
175
        }
173
176
        if (fusionDetection) checkFusion();
174
177
 
175
 
        for (ii= scene->interactions->begin(); ii!=iiEnd ; ++ii) {
176
 
                if ((*ii)->isReal()) {
 
178
        #ifdef YADE_OPENMP
 
179
        #pragma omp parallel for schedule(guided) num_threads(ompThreads>0 ? min(ompThreads,omp_get_max_threads()) : omp_get_max_threads())
 
180
        for(long i=0; i<size; i++){
 
181
                const shared_ptr<Interaction>& interaction=(*scene->interactions)[i];
 
182
        #else
 
183
        FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){
 
184
        #endif
 
185
                if (interaction->isReal()) {
177
186
                        CapillaryPhys* cundallContactPhysics=NULL;
178
187
                        MindlinCapillaryPhys* mindlinContactPhysics=NULL;
179
 
                        if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys*>((*ii)->phys.get());//use CapillaryPhys for linear model
180
 
                        else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>((*ii)->phys.get());//use MindlinCapillaryPhys for hertz model
 
188
                        if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys*>(interaction->phys.get());//use CapillaryPhys for linear model
 
189
                        else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>(interaction->phys.get());//use MindlinCapillaryPhys for hertz model
181
190
 
182
191
                        if ((hertzOn && mindlinContactPhysics->meniscus) || (!hertzOn && cundallContactPhysics->meniscus)) {
183
192
                                if (fusionDetection) {//version with effect of fusion
192
201
                                        //LINEAR VERSION : capillary force is divided by (fusionNumber + 1) - NOTE : any decreasing function of fusionNumber can be considered in fact
193
202
                                        else if (fusionNumber !=0) hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap /= (fusionNumber+1.);
194
203
                                }
195
 
                                scene->forces.addForce((*ii)->getId1(), hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap);
196
 
                                scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap));
 
204
                                scene->forces.addForce(interaction->getId1(), hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap);
 
205
                                scene->forces.addForce(interaction->getId2(),-(hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap));
197
206
                        }
198
207
                }
199
208
        }
211
220
void Law2_ScGeom_CapillaryPhys_Capillarity::checkFusion()
212
221
{
213
222
        //Reset fusion numbers
214
 
        InteractionContainer::iterator ii    = scene->interactions->begin();
215
 
        InteractionContainer::iterator iiEnd = scene->interactions->end();
216
 
        for( ; ii!=iiEnd ; ++ii ) {
217
 
                if ((*ii)->isReal()) {
218
 
                        if (!hertzOn) static_cast<CapillaryPhys*>((*ii)->phys.get())->fusionNumber=0;
219
 
                        else static_cast<MindlinCapillaryPhys*>((*ii)->phys.get())->fusionNumber=0;
220
 
                }
 
223
        #ifdef YADE_OPENMP
 
224
        const long size=scene->interactions->size();
 
225
        #pragma omp parallel for schedule(guided) num_threads(ompThreads>0 ? min(ompThreads,omp_get_max_threads()) : omp_get_max_threads())
 
226
        for(long i=0; i<size; i++){
 
227
                const shared_ptr<Interaction>& interaction=(*scene->interactions)[i];
 
228
        #else
 
229
        FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){
 
230
        #endif
 
231
                if ( !interaction->isReal()) continue;
 
232
                if (!hertzOn) static_cast<CapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
 
233
                else static_cast<MindlinCapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
221
234
        }
222
235
 
223
236
        list< shared_ptr<Interaction> >::iterator firstMeniscus, lastMeniscus, currentMeniscus;
535
548
                interactionsOnBody[i].clear();
536
549
        }
537
550
 
538
 
        InteractionContainer::iterator ii    = scene->interactions->begin();
539
 
        InteractionContainer::iterator iiEnd = scene->interactions->end();
540
 
        for(  ; ii!=iiEnd ; ++ii ) {
541
 
                if ((*ii)->isReal()) {
542
 
                        if (static_cast<CapillaryPhys*>((*ii)->phys.get())->meniscus) insert(*ii);
 
551
        FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ // parallel version using Engine::ompThreads variable not accessible, this function is not one of the Engine..
 
552
                if (I->isReal()) {
 
553
                    if (static_cast<CapillaryPhys*>(I->phys.get())->meniscus) insert(I);
543
554
                }
544
555
        }
545
556