Law2_ScGeom_CapillaryPhys_Capillarity.cpp 24.1 KB
Newer Older
1
/*************************************************************************
2
3
*  Copyright (C) 2006 by luc Scholtes                                    *
*  luc.scholtes@hmg.inpg.fr                                              *
4
5
6
7
8
*                                                                        *
*  This program is free software; it is licensed under the terms of the  *
*  GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/

9
#include "Law2_ScGeom_CapillaryPhys_Capillarity.hpp"
10
11
#include <pkg/common/ElastMat.hpp>
#include <pkg/dem/ScGeom.hpp>
12

13
14
15
16
17
#include <pkg/dem/CapillaryPhys.hpp>
#include <pkg/dem/HertzMindlin.hpp>
#include <core/Omega.hpp>
#include <core/Scene.hpp>
#include <lib/base/Math.hpp>
18

19
YADE_PLUGIN((Law2_ScGeom_CapillaryPhys_Capillarity));
20

21
void Law2_ScGeom_CapillaryPhys_Capillarity::postLoad(Law2_ScGeom_CapillaryPhys_Capillarity&){
22

23
  capillary = shared_ptr<capillarylaw>(new capillarylaw);
24
25
26
27
28
29
30
31
32
33
34
  capillary->fill("M(r=1)");
  capillary->fill("M(r=1.1)");
  capillary->fill("M(r=1.25)");
  capillary->fill("M(r=1.5)");
  capillary->fill("M(r=1.75)");
  capillary->fill("M(r=2)");
  capillary->fill("M(r=3)");
  capillary->fill("M(r=4)");
  capillary->fill("M(r=5)");
  capillary->fill("M(r=10)");
}
35
36


Bruno Chareyre's avatar
   
Bruno Chareyre committed
37
MeniscusParameters::MeniscusParameters()
38
39
40
41
42
{
        V = 0;
        F = 0;
        delta1 = 0;
        delta2 = 0;
43
}
44

Bruno Chareyre's avatar
   
Bruno Chareyre committed
45
MeniscusParameters::MeniscusParameters(const MeniscusParameters &source)
46
47
48
49
50
51
52
{
        V = source.V;
        F = source.F;
        delta1 = source.delta1;
        delta2 = source.delta2;
}

Bruno Chareyre's avatar
   
Bruno Chareyre committed
53
MeniscusParameters::~MeniscusParameters()
54
{}
Bruno Chareyre's avatar
   
Bruno Chareyre committed
55

56
void Law2_ScGeom_CapillaryPhys_Capillarity::action()
57
{
Anton Gladky's avatar
Anton Gladky committed
58
	if (!scene) cerr << "scene not defined!";
59
	if (!capillary) postLoad(*this);//when the script does not define arguments, postLoad is never called
60
	shared_ptr<BodyContainer>& bodies = scene->bodies;
61
62
63
64
65
66
67
68
	
	//check for contact model once (assuming that contact model does not change)
	if (!hertzInitialized){//NOTE: We are assuming that only one type is used in one simulation here
		FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
			if (I->isReal()) {
				if (CapillaryPhys::getClassIndexStatic()==I->phys->getClassIndex()) hertzOn=false;
				else if (MindlinCapillaryPhys::getClassIndexStatic()==I->phys->getClassIndex()) hertzOn=true;
				else LOG_ERROR("The capillary law is not implemented for interactions using"<<I->phys->getClassName());
69
70
				bodiesMenisciiList.initialized = false;//must be re-initialized after creation of first real contact in the model
				hertzInitialized = true;
71
72
73
74
75
76
				break;
			}
		}
	}
	
	if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene,hertzOn);
77

78
	FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){ // could be done in parallel as soon as OpenMPVector class (lib/base/openmp-accu.hpp) is extended See http://www.mail-archive.com/yade-dev@lists.launchpad.net/msg10842.html and msg11238.html
79
		/// interaction is real
80
		if (interaction->isReal()) {
81
82
83
84
85
86
87
			CapillaryPhys* cundallContactPhysics=NULL;
			MindlinCapillaryPhys* mindlinContactPhysics=NULL;

			/// contact physics depends on the contact law, that is used (either linear model or hertz model)
			if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys*>(interaction->phys.get());//use CapillaryPhys for linear model
			else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>(interaction->phys.get());//use MindlinCapillaryPhys for hertz model
			
88
89
90
91
92
			unsigned int id1 = interaction->getId1();
			unsigned int id2 = interaction->getId2();
			Body* b1 = (*bodies)[id1].get();
			Body* b2 = (*bodies)[id2].get();

93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
			/// interaction geometry search (this test is to compute capillarity only between spheres (probably a better way to do that)
			int geometryIndex1 = (*bodies)[id1]->shape->getClassIndex(); // !!!
			int geometryIndex2 = (*bodies)[id2]->shape->getClassIndex();
			if (!(geometryIndex1 == geometryIndex2)) continue;

			/// definition of interacting objects (not necessarily in contact)
			ScGeom* currentContactGeometry = static_cast<ScGeom*>(interaction->geom.get());

			/// Capillary components definition:
			Real liquidTension = 0.073; 	// superficial water tension N/m (0.073 is water tension at 20 Celsius degrees)

			/// Interacting Grains:
			// If you want to define a ratio between YADE sphere size and real sphere size
			Real alpha=1;
			Real R1 = 0;
			R1=alpha*std::min(currentContactGeometry->radius2,currentContactGeometry->radius1) ;
			Real R2 = 0;
			R2=alpha*std::max(currentContactGeometry->radius2,currentContactGeometry->radius1) ;

			/// intergranular distance
113
			Real D = - alpha * currentContactGeometry->penetrationDepth;
114

115
116
			if ( D<0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere
				D=max(0.,D); // defines fCap when spheres interpenetrate. D<0 leads to wrong interpolation as D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
117
				if (!hertzOn) {
118
					if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert(interaction);
119
120
					cundallContactPhysics->meniscus=true;
				} else {
121
					if (fusionDetection && !mindlinContactPhysics->meniscus) bodiesMenisciiList.insert(interaction);
122
					mindlinContactPhysics->meniscus=true;
123
				}
124
125
126
127
128
			}
			Real Dinterpol = D/R2;

			/// Suction (Capillary pressure):
			Real Pinterpol = 0;
129
130
131
132
			if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : capillaryPressure*(R2/liquidTension);
			else Pinterpol = mindlinContactPhysics->isBroken ? 0 : capillaryPressure*(R2/liquidTension);
			if (!hertzOn) cundallContactPhysics->capillaryPressure = capillaryPressure;
			else mindlinContactPhysics->capillaryPressure = capillaryPressure;
133
134
135
136
137
138

			/// Capillary solution finder:
			if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) {
				int* currentIndexes =  hertzOn? mindlinContactPhysics->currentIndexes : cundallContactPhysics->currentIndexes;
				//If P=0, we use null solution
				MeniscusParameters
139
				solution(Pinterpol? capillary->interpolate(R1,R2,Dinterpol, Pinterpol, currentIndexes) : MeniscusParameters());
140
141
				/// capillary adhesion force
				Real Finterpol = solution.F;
142
				Vector3r fCap = - Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
143
144
				if (!hertzOn) cundallContactPhysics->fCap = fCap;
				else mindlinContactPhysics->fCap = fCap;
145
146
147
				/// meniscus volume
				Real Vinterpol = solution.V * pow(R2/alpha,3);
				if (!hertzOn) {
148
					cundallContactPhysics->vMeniscus = Vinterpol;
149
150
151
					if (Vinterpol != 0) cundallContactPhysics->meniscus = true;
					else cundallContactPhysics->meniscus = false;
				} else {
152
					mindlinContactPhysics->vMeniscus = Vinterpol;
153
154
155
156
					if (Vinterpol != 0) mindlinContactPhysics->meniscus = true;
					else mindlinContactPhysics->meniscus = false;
				}
				if (!Vinterpol) {
157
					if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove(interaction);
158
					if (D>0) scene->interactions->requestErase(interaction);
159
					else if ((Pinterpol > 0) && (showError)) {
160
						LOG_ERROR("No meniscus found at a contact. capillaryPressure may be too large wrt. the loaded data files."); // V=0 at a contact reveals a problem if and only if uc* > 0
161
162
						showError = false;//show error message once / avoid console spam
					}
163
164
165
166
167
168
169
170
				}
				/// wetting angles
				if (!hertzOn) {
					cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
					cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
				} else {
					mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
					mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
171
				}
172
			}
173
		///interaction is not real	//If the interaction is not real, it should not be in the list
174
		} else if (fusionDetection) bodiesMenisciiList.remove(interaction);
175
176
	}
	if (fusionDetection) checkFusion();
Bruno Chareyre's avatar
   
Bruno Chareyre committed
177

178
179
180
181
182
183
184
185
        #ifdef YADE_OPENMP
        const long size=scene->interactions->size();
        #pragma omp parallel for schedule(guided) num_threads(ompThreads>0 ? min(ompThreads,omp_get_max_threads()) : omp_get_max_threads())
        for(long i=0; i<size; i++){
            const shared_ptr<Interaction>& interaction=(*scene->interactions)[i];
        #else
        FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){
        #endif
186
		if (interaction->isReal()) {
187
188
			CapillaryPhys* cundallContactPhysics=NULL;
			MindlinCapillaryPhys* mindlinContactPhysics=NULL;
189
190
			if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys*>(interaction->phys.get());//use CapillaryPhys for linear model
			else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>(interaction->phys.get());//use MindlinCapillaryPhys for hertz model
191

192
193
			if ((hertzOn && mindlinContactPhysics->meniscus) || (!hertzOn && cundallContactPhysics->meniscus)) {
				if (fusionDetection) {//version with effect of fusion
Bruno Chareyre's avatar
Bruno Chareyre committed
194
195
					//BINARY VERSION : if fusionNumber!=0 then no capillary force
					short int& fusionNumber = hertzOn?mindlinContactPhysics->fusionNumber:cundallContactPhysics->fusionNumber;
196
					if (binaryFusion) {
Anton Gladky's avatar
Anton Gladky committed
197
						if (fusionNumber!=0) {	//cerr << "fusion" << endl;
198
							hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap = Vector3r::Zero();
Bruno Chareyre's avatar
Bruno Chareyre committed
199
							continue;
200
201
						}
					}
Bruno Chareyre's avatar
Bruno Chareyre committed
202
					//LINEAR VERSION : capillary force is divided by (fusionNumber + 1) - NOTE : any decreasing function of fusionNumber can be considered in fact
203
					else if (fusionNumber !=0) hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap /= (fusionNumber+1.);
204
				}
205
206
				scene->forces.addForce(interaction->getId1(),-(hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap));
				scene->forces.addForce(interaction->getId2(),  hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap );
207
			}
208
		}
209
	}
Bruno Chareyre's avatar
   
Bruno Chareyre committed
210
}
211
212
213

capillarylaw::capillarylaw()
{}
214
215
216
217
218
219
220

void capillarylaw::fill(const char* filename)
{
        data_complete.push_back(Tableau(filename));

}

221
void Law2_ScGeom_CapillaryPhys_Capillarity::checkFusion()
Bruno Chareyre's avatar
   
Bruno Chareyre committed
222
223
{
	//Reset fusion numbers
224
	FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){ // same remark for parallel loops, the most problematic part ?
225
226
227
228
		if ( interaction->isReal()) {
			if (!hertzOn) static_cast<CapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
			else static_cast<MindlinCapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
		}
229
	}
230

Bruno Chareyre's avatar
   
Bruno Chareyre committed
231
	list< shared_ptr<Interaction> >::iterator firstMeniscus, lastMeniscus, currentMeniscus;
232
	Real angle1 = -1.0; Real angle2 = -1.0;
Anton Gladky's avatar
Anton Gladky committed
233
	
234
	for ( int i=0; i< bodiesMenisciiList.size(); ++i ) { // i is the index (or id) of the body being tested
235
236
237
238
		CapillaryPhys* cundallInteractionPhysics1=NULL;
		MindlinCapillaryPhys* mindlinInteractionPhysics1=NULL;
		CapillaryPhys* cundallInteractionPhysics2=NULL;
		MindlinCapillaryPhys* mindlinInteractionPhysics2=NULL;
239
		if ( !bodiesMenisciiList[i].empty() ) {
Bruno Chareyre's avatar
   
Bruno Chareyre committed
240
			lastMeniscus = bodiesMenisciiList[i].end();
241
			for ( firstMeniscus=bodiesMenisciiList[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus ) { //FOR EACH MENISCUS ON THIS BODY...
Bruno Chareyre's avatar
   
Bruno Chareyre committed
242
				currentMeniscus = firstMeniscus; ++currentMeniscus;
243
				if (!hertzOn) {
244
					cundallInteractionPhysics1 = YADE_CAST<CapillaryPhys*>((*firstMeniscus)->phys.get());
245
246
					if (i == (*firstMeniscus)->getId1()) angle1=cundallInteractionPhysics1->Delta1;//get angle of meniscus1 on body i
					else angle1=cundallInteractionPhysics1->Delta2;
247
				}
248
				else {
249
					mindlinInteractionPhysics1 = YADE_CAST<MindlinCapillaryPhys*>((*firstMeniscus)->phys.get());
250
251
					if (i == (*firstMeniscus)->getId1()) angle1=mindlinInteractionPhysics1->Delta1;//get angle of meniscus1 on body i
					else angle1=mindlinInteractionPhysics1->Delta2;
252
				}
Bruno Chareyre's avatar
   
Bruno Chareyre committed
253
				for ( ;currentMeniscus!= lastMeniscus; ++currentMeniscus) {//... CHECK FUSION WITH ALL OTHER MENISCII ON THE BODY
254
					if (!hertzOn) {
255
						cundallInteractionPhysics2 = YADE_CAST<CapillaryPhys*>((*currentMeniscus)->phys.get());
256
257
						if (i == (*currentMeniscus)->getId1()) angle2=cundallInteractionPhysics2->Delta1;//get angle of meniscus2 on body i
						else angle2=cundallInteractionPhysics2->Delta2;
258
					}
259
					else {
260
						mindlinInteractionPhysics2 = YADE_CAST<MindlinCapillaryPhys*>((*currentMeniscus)->phys.get());
261
262
						if (i == (*currentMeniscus)->getId1()) angle2=mindlinInteractionPhysics2->Delta1;//get angle of meniscus2 on body i
						else angle2=mindlinInteractionPhysics2->Delta2;
263
					}
Anton Gladky's avatar
Anton Gladky committed
264
					if (angle1==0 || angle2==0) cerr << "THIS SHOULD NOT HAPPEN!!"<< endl;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
265

Anton Gladky's avatar
Anton Gladky committed
266
					//cerr << "angle1 = " << angle1 << " | angle2 = " << angle2 << endl;
267

268
269
					Vector3r normalFirstMeniscus = YADE_CAST<ScGeom*>((*firstMeniscus)->geom.get())->normal;
					Vector3r normalCurrentMeniscus = YADE_CAST<ScGeom*>((*currentMeniscus)->geom.get())->normal;
270

Bruno Chareyre's avatar
   
Bruno Chareyre committed
271
					Real normalDot = 0;
272
					if ((*firstMeniscus)->getId1() ==  (*currentMeniscus)->getId1() ||  (*firstMeniscus)->getId2()  == (*currentMeniscus)->getId2()) normalDot = normalFirstMeniscus.dot(normalCurrentMeniscus);
273
					else normalDot = - (normalFirstMeniscus.dot(normalCurrentMeniscus));
274

Bruno Chareyre's avatar
   
Bruno Chareyre committed
275
276
277
					Real normalAngle = 0;
					if (normalDot >= 0 ) normalAngle = Mathr::FastInvCos0(normalDot);
					else normalAngle = ((Mathr::PI) - Mathr::FastInvCos0(-(normalDot)));
278

279
					if ((angle1+angle2)*Mathr::DEG_TO_RAD > normalAngle) {
280
281
						if (!hertzOn) {++(cundallInteractionPhysics1->fusionNumber); ++(cundallInteractionPhysics2->fusionNumber);}//count +1 if 2 meniscii are overlaping
						else {++(mindlinInteractionPhysics1->fusionNumber); ++(mindlinInteractionPhysics2->fusionNumber);}
282
					}
Bruno Chareyre's avatar
   
Bruno Chareyre committed
283
284
285
286
287
				}
			}
		}
	}
}
288

289
MeniscusParameters capillarylaw::interpolate(Real R1, Real R2, Real D, Real P, int* index)
Anton Gladky's avatar
Anton Gladky committed
290
{	//cerr << "interpolate" << endl;
291
292
293
294
295
296
297
        if (R1 > R2) {
                Real R3 = R1;
                R1 = R2;
                R2 = R3;
        }

        Real R = R2/R1;
Anton Gladky's avatar
Anton Gladky committed
298
        //cerr << "R = " << R << endl;
299

Bruno Chareyre's avatar
   
Bruno Chareyre committed
300
301
302
        MeniscusParameters result_inf;
        MeniscusParameters result_sup;
        MeniscusParameters result;
303
        int i = 0;
304

305
        for ( ; i < (NB_R_VALUES); i++)
306
307
        {
                Real data_R = data_complete[i].R;
Anton Gladky's avatar
Anton Gladky committed
308
                //cerr << "i = " << i << endl;
309

310
                if (data_R > R)	// Attention a l'ordre ds lequel vont etre ranges les tableau R (croissant)
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
                {
                        Tableau& tab_inf=data_complete[i-1];
                        Tableau& tab_sup=data_complete[i];

                        Real r=(R-tab_inf.R)/(tab_sup.R-tab_inf.R);

                        result_inf = tab_inf.Interpolate2(D,P,index[0], index[1]);
                        result_sup = tab_sup.Interpolate2(D,P,index[2], index[3]);

                        result.V = result_inf.V*(1-r) + r*result_sup.V;
                        result.F = result_inf.F*(1-r) + r*result_sup.F;
                        result.delta1 = result_inf.delta1*(1-r) + r*result_sup.delta1;
                        result.delta2 = result_inf.delta2*(1-r) + r*result_sup.delta2;

                        i=NB_R_VALUES;
Anton Gladky's avatar
Anton Gladky committed
326
                        //cerr << "i = " << i << endl;
327
328
329
330
331
332

                }
                else if (data_complete[i].R == R)
                {
                        result = data_complete[i].Interpolate2(D,P, index[0], index[1]);
                        i=NB_R_VALUES;
Anton Gladky's avatar
Anton Gladky committed
333
                        //cerr << "i = " << i << endl;
334
335
336
337
338
339
340
341
342
343
344
                }
        }
        return result;
}

Tableau::Tableau()
{}

Tableau::Tableau(const char* filename)

{
Anton Gladky's avatar
Anton Gladky committed
345
        ifstream file (filename);
346
        file >> R;
Anton Gladky's avatar
Anton Gladky committed
347
        //cerr << "r = " << R << endl;
348
349
350
351
        int n_D;	//number of D values
        file >> n_D;

        if (!file.is_open())
352
353
354
355
	{
		static bool first=true;
		if(first)
		{
Anton Gladky's avatar
Anton Gladky committed
356
	                cout << "WARNING: cannot open files used for capillary law, all forces will be null. Instructions on how to download and install them is found here : https://yade-dem.org/wiki/CapillaryTriaxialTest." << endl;
357
358
359
360
			first=false;
		}
		return;
	}
361
362
363
        for (int i=0; i<n_D; i++)
                full_data.push_back(TableauD(file));
        file.close();
364

365
366
367
368
369
}

Tableau::~Tableau()
{}

Bruno Chareyre's avatar
   
Bruno Chareyre committed
370
MeniscusParameters Tableau::Interpolate2(Real D, Real P, int& index1, int& index2)
371

Anton Gladky's avatar
Anton Gladky committed
372
{	//cerr << "interpolate2" << endl;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
373
374
375
        MeniscusParameters result;
        MeniscusParameters result_inf;
        MeniscusParameters result_sup;
376

377
        for ( unsigned int i=0; i < full_data.size(); ++i)
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
        {
                if (full_data[i].D > D )	// ok si D rang�s ds l'ordre croissant

                {
                        Real rD = (D-full_data[i-1].D)/(full_data[i].D-full_data[i-1].D);

                        result_inf = full_data[i-1].Interpolate3(P, index1);
                        result_sup = full_data[i].Interpolate3(P, index2);

                        result.V = result_inf.V*(1-rD) + rD*result_sup.V;
                        result.F = result_inf.F*(1-rD) + rD*result_sup.F;
                        result.delta1 = result_inf.delta1*(1-rD) + rD*result_sup.delta1;
                        result.delta2 = result_inf.delta2*(1-rD) + rD*result_sup.delta2;

                        i = full_data.size();
                }
                else if (full_data[i].D == D)
                {
                        result=full_data[i].Interpolate3(P, index1);

                        i=full_data.size();
                }

        }
        return result;
}

TableauD::TableauD()
406
{}
407

Anton Gladky's avatar
Anton Gladky committed
408
TableauD::TableauD(ifstream& file)
409
410
411
412
413
{
        int i=0;
        Real x;
        int n_lines;	//pb: n_lines is real!!!
        file >> n_lines;
Anton Gladky's avatar
Anton Gladky committed
414
        //cout << n_lines << endl;
415
416
417
418

        file.ignore(200, '\n'); // saute les caract�res (200 au maximum) jusque au caract�re \n (fin de ligne)*_

        if (n_lines!=0)
419
                for (; i<n_lines; ++i) {
420
421
422
423
424
425
426
427
428
429
                        data.push_back(vector<Real> ());
                        for (int j=0; j < 6; ++j)	// [D,P,V,F,delta1,delta2]
                        {
                                file >> x;
                                data[i].push_back(x);
                        }
                }
        D = data[i-1][0];
}

Bruno Chareyre's avatar
   
Bruno Chareyre committed
430
MeniscusParameters TableauD::Interpolate3(Real P, int& index)
Anton Gladky's avatar
Anton Gladky committed
431
{	//cerr << "interpolate3" << endl;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
432
        MeniscusParameters result;
433
        int dataSize = data.size();
434

435
436
        if (index < dataSize && index>0)
        {
Bruno Chareyre's avatar
   
Bruno Chareyre committed
437
        	if (data[index][1] >= P && data[index-1][1] < P)
438
        	{
439
        		//compteur1+=1;
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
        		Real Pinf=data[index-1][1];
                        Real Finf=data[index-1][3];
                        Real Vinf=data[index-1][2];
                        Real Delta1inf=data[index-1][4];
                        Real Delta2inf=data[index-1][5];

                        Real Psup=data[index][1];
                        Real Fsup=data[index][3];
                        Real Vsup=data[index][2];
                        Real Delta1sup=data[index][4];
                        Real Delta2sup=data[index][5];

                        result.V = Vinf+((Vsup-Vinf)/(Psup-Pinf))*(P-Pinf);
                        result.F = Finf+((Fsup-Finf)/(Psup-Pinf))*(P-Pinf);
                        result.delta1 = Delta1inf+((Delta1sup-Delta1inf)/(Psup-Pinf))*(P-Pinf);
                        result.delta2 = Delta2inf+((Delta2sup-Delta2inf)/(Psup-Pinf))*(P-Pinf);
                        return result;
457

458
459
        	}
        }
460
	//compteur2+=1;
461
462
        for (int k=1; k < dataSize; ++k) 	// Length(data) ??

Anton Gladky's avatar
Anton Gladky committed
463
        {	//cerr << "k = " << k << endl;
464
465
                if ( data[k][1] > P) 	// OK si P rangés ds l'ordre croissant

Anton Gladky's avatar
Anton Gladky committed
466
                {	//cerr << "if" << endl;
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
                        Real Pinf=data[k-1][1];
                        Real Finf=data[k-1][3];
                        Real Vinf=data[k-1][2];
                        Real Delta1inf=data[k-1][4];
                        Real Delta2inf=data[k-1][5];

                        Real Psup=data[k][1];
                        Real Fsup=data[k][3];
                        Real Vsup=data[k][2];
                        Real Delta1sup=data[k][4];
                        Real Delta2sup=data[k][5];

                        result.V = Vinf+((Vsup-Vinf)/(Psup-Pinf))*(P-Pinf);
                        result.F = Finf+((Fsup-Finf)/(Psup-Pinf))*(P-Pinf);
                        result.delta1 = Delta1inf+((Delta1sup-Delta1inf)/(Psup-Pinf))*(P-Pinf);
                        result.delta2 = Delta2inf+((Delta2sup-Delta2inf)/(Psup-Pinf))*(P-Pinf);
                        index = k;

                        k=dataSize;
                }
                else if (data[k][1] == P)

Anton Gladky's avatar
Anton Gladky committed
489
                {	//cerr << "elseif" << endl;
490
491
492
493
494
495
496
497
498
499
500
501
502
503
                        result.V = data[k][2];
                        result.F = data[k][3];
                        result.delta1 = data[k][4];
                        result.delta2 = data[k][5];
                        index = k;

                        k=dataSize;
                }

        }
        return result;
}

TableauD::~TableauD()
504
{}
505
506
507

std::ostream& operator<<(std::ostream& os, Tableau& T)
{
Anton Gladky's avatar
Anton Gladky committed
508
        os << "Tableau : R=" << T.R << endl;
509
        for (unsigned int i=0; i<T.full_data.size(); i++) {
Anton Gladky's avatar
Anton Gladky committed
510
                os << "TableauD : D=" << T.full_data[i].D << endl;
511
512
                for (unsigned int j=0; j<T.full_data[i].data.size();j++) {
                        for (unsigned int k=0; k<T.full_data[i].data[j].size(); k++)
513
                                os << T.full_data[i].data[j][k] << " ";
Anton Gladky's avatar
Anton Gladky committed
514
                        os << endl;
515
516
                }
        }
Anton Gladky's avatar
Anton Gladky committed
517
        os << endl;
518
519
        return os;
}
Bruno Chareyre's avatar
   
Bruno Chareyre committed
520

521
522
//what is this function for? it is never called...
//TODO: remove?
523
BodiesMenisciiList::BodiesMenisciiList(Scene * scene, bool hertzOn)
Bruno Chareyre's avatar
   
Bruno Chareyre committed
524
525
{
	initialized=false;
526
	prepare(scene, hertzOn);
Bruno Chareyre's avatar
   
Bruno Chareyre committed
527
528
}

529

530
bool BodiesMenisciiList::prepare(Scene * scene, bool hertzOn)
Bruno Chareyre's avatar
   
Bruno Chareyre committed
531
{
Anton Gladky's avatar
Anton Gladky committed
532
	//cerr << "preparing bodiesInteractionsList" << endl;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
533
	interactionsOnBody.clear();
534
	shared_ptr<BodyContainer>& bodies = scene->bodies;
535

536
	Body::id_t MaxId = -1;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
537
538
539
540
541
542
543
544
545
546
547
	BodyContainer::iterator bi    = bodies->begin();
	BodyContainer::iterator biEnd = bodies->end();
	for(  ; bi!=biEnd ; ++bi )
	{
		MaxId=max(MaxId, (*bi)->getId());
	}
	interactionsOnBody.resize(MaxId+1);
	for ( unsigned int i=0; i<interactionsOnBody.size(); ++i )
	{
		interactionsOnBody[i].clear();
	}
548
	
549
550
	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ // parallel version using Engine::ompThreads variable not accessible, this function is not one of the Engine..
                if (I->isReal()) {
551
552
			if (!hertzOn) {if (static_cast<CapillaryPhys*>(I->phys.get())->meniscus) insert(I);}
			else {if (static_cast<MindlinCapillaryPhys*>(I->phys.get())->meniscus) insert(I);}
Bruno Chareyre's avatar
   
Bruno Chareyre committed
553
554
                }
        }
555

Bruno Chareyre's avatar
   
Bruno Chareyre committed
556
557
558
559
560
	return initialized=true;
}

bool BodiesMenisciiList::insert(const shared_ptr< Interaction >& interaction)
{
561
  checkLengthBuffer(interaction);
Bruno Chareyre's avatar
   
Bruno Chareyre committed
562
563
	interactionsOnBody[interaction->getId1()].push_back(interaction);
	interactionsOnBody[interaction->getId2()].push_back(interaction);
564
	return true;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
565
566
567
568
569
}


bool BodiesMenisciiList::remove(const shared_ptr< Interaction >& interaction)
{
570
  checkLengthBuffer(interaction);
Bruno Chareyre's avatar
   
Bruno Chareyre committed
571
	interactionsOnBody[interaction->getId1()].remove(interaction);
572
	interactionsOnBody[interaction->getId2()].remove(interaction);
Bruno Chareyre's avatar
   
Bruno Chareyre committed
573
574
575
	return true;
}

576
577
578
579
580
581
582
void BodiesMenisciiList::checkLengthBuffer(const shared_ptr<Interaction>& interaction) {
	Body::id_t maxBodyId = std::max(interaction->getId1(), interaction->getId2());
	if (unsigned(maxBodyId) >= interactionsOnBody.size()) {
		interactionsOnBody.resize(maxBodyId+1);
	}
}

Bruno Chareyre's avatar
   
Bruno Chareyre committed
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
list< shared_ptr<Interaction> >&  BodiesMenisciiList::operator[] (int index)
{
	return interactionsOnBody[index];
}

int BodiesMenisciiList::size()
{
	return interactionsOnBody.size();
}

void BodiesMenisciiList::display()
{
	list< shared_ptr<Interaction> >::iterator firstMeniscus;
	list< shared_ptr<Interaction> >::iterator lastMeniscus;
	for ( unsigned int i=0; i<interactionsOnBody.size(); ++i )
	{
		if ( !interactionsOnBody[i].empty() )
		{
			lastMeniscus = interactionsOnBody[i].end();
Anton Gladky's avatar
Anton Gladky committed
602
			//cerr << "size = "<<interactionsOnBody[i].size() << " empty="<<interactionsOnBody[i].empty() <<endl;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
603
604
			for ( firstMeniscus=interactionsOnBody[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus )
			{
605
				if ( *firstMeniscus ){
Bruno Chareyre's avatar
   
Bruno Chareyre committed
606
					if ( firstMeniscus->get() )
Anton Gladky's avatar
Anton Gladky committed
607
608
						cerr << "(" << ( *firstMeniscus )->getId1() << ", " << ( *firstMeniscus )->getId2() <<") ";
					else cerr << "(void)";
609
				}
Bruno Chareyre's avatar
   
Bruno Chareyre committed
610
			}
Anton Gladky's avatar
Anton Gladky committed
611
			cerr << endl;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
612
		}
Anton Gladky's avatar
Anton Gladky committed
613
		else cerr << "empty" << endl;
Bruno Chareyre's avatar
   
Bruno Chareyre committed
614
615
616
617
618
619
620
	}
}

BodiesMenisciiList::BodiesMenisciiList()
{
	initialized=false;
}
621