NeoPZ
pzmultiphase.cpp
Go to the documentation of this file.
1 //
2 // pzmultiphase.h
3 // PZ
4 //
5 // Created by Omar Duran on 19/08/2013.
6 // Copyright (c) 2012 LabMec-Unicamp. All rights reserved.
7 //
8 
9 #include "pzmultiphase.h"
10 #include "pzlog.h"
11 #include "pzbndcond.h"
12 #include "pzfmatrix.h"
13 #include "pzaxestools.h"
14 #include <math.h>
15 
16 #include <iostream>
17 
18 #ifdef LOG4CXX
19 static LoggerPtr logger(Logger::getLogger("pz.multiphase"));
20 #endif
21 
22 #ifdef LOG4CXX
23 static LoggerPtr logdata(Logger::getLogger("pz.material.multiphase.data"));
24 #endif
25 
26 TPZMultiphase::TPZMultiphase():
27 TPZRegisterClassId(&TPZMultiphase::ClassId),
28 TPZDiscontinuousGalerkin()
29 {
30  fDim = 2;
31  fTheta = 1.0;
32  fGamma = 1.0;
33  fDeltaT = 1.0;
34  fmatId = 0;
35  fLref = 1.0;
36  fPref= 1.0;
37  fKref = 1.0;
38  fRhoref = 1.0;
39  fEtaref = 1.0;
40  fxi = 1.0;
41 }
42 
43 TPZMultiphase::TPZMultiphase(int matid, int dim):
44 TPZRegisterClassId(&TPZMultiphase::ClassId),
45 TPZDiscontinuousGalerkin(matid)
46 {
47  // Two-dimensional analysis
48 
49  if(dim<0 || dim >2){
50  DebugStop();
51  }
52 
53  fDim = dim;
54  fTheta = 1.0;
55  fGamma = 1.0;
56  fDeltaT = 1.0;
57  fmatId = matid;
58  fLref = 1.0;
59  fPref= 1.0;
60  fKref = 1.0;
61  fRhoref = 1.0;
62  fEtaref = 1.0;
63  fxi = 1.0;
64 }
65 
66 
67 TPZMultiphase::~TPZMultiphase()
68 {
69 }
70 
71 int TPZMultiphase::Dimension() const {return fDim;};
72 
73 int TPZMultiphase::MatId() {return fmatId;}
74 
75 int TPZMultiphase::NStateVariables() const {return 8;}
76 
77 void TPZMultiphase::Print(std::ostream &out) {
78  out << "name of material : " << Name() << "\n";
79  out << "Coeficient which multiplies the gradient operator "<< "my var" << std::endl;
80  out << "Base Class properties :";
81  TPZMaterial::Print(out);
82  out << "\n";
83 }
84 
85 
86 
87 // Data set
88 
93 STATE TPZMultiphase::LameLambda()
94 {
95  REAL lamelambda = 1.0;
96  return lamelambda;
97 }
98 
103 STATE TPZMultiphase::LameLambdaU()
104 {
105  REAL lamelambdaU = 1.0;
106  return lamelambdaU;
107 }
108 
113 STATE TPZMultiphase::LameMu()
114 {
115  REAL lameMu = 0.5;
116  return lameMu;
117 }
118 
123 STATE TPZMultiphase::BiotAlpha()
124 {
125  REAL alpha = 0.0;
126  return alpha;
127 }
128 
133 STATE TPZMultiphase::Se()
134 {
135  REAL Se = 0.0;
136  return Se;
137 }
138 
139 
141 void TPZMultiphase::CapillaryPressure(REAL Sw, REAL &pc, REAL &DpcDSo){
142 // Bentsen and Anli Capillary Pressure model
143  REAL PsiToPa = 6894.76;
144  REAL Pct = 1.0, Pcs = 0.5, Swi = 0.0;
145 
146 // pc = 0.0;
147 // DpcDSo = 0.0;
148 
149  pc = (Pct - Pcs * log ((0.000000001 + Sw - Swi)/(1.0 - Swi)))*PsiToPa/this->fPref;
150  DpcDSo = -(Pcs/(0.000000001 + Sw - Swi))*PsiToPa/this->fPref;
151 }
152 
154 void TPZMultiphase::Kro(REAL Sw, REAL &Kro, REAL &dKroDSw){
155 
156 
157 // // Zero model
158 // Kro = 0.0;
159 // dKroDSw = 0.0;
160 
161  // linear model
162  Kro = 1.0-Sw;
163  dKroDSw = -1.0;
164 
165  // // Non-linear model
166  // Kro = (1-Sw)*(1-Sw);
167  // dKroDSw = 2.0*(1-Sw)*(-1.0);
168 
169  // // Non-linear model
170  // Kro = (1-Sw)*(1-Sw)*(1-Sw);
171  // dKroDSw = 3.0*(1-Sw)*(1-Sw)*(-1.0);
172 }
173 
174 
176 void TPZMultiphase::Krw(REAL Sw, REAL &Krw, REAL &dKrwDSw){
177 
178 // // Zero model
179 // Krw = 1;
180 // dKrwDSw = 0.0;
181 
182  // linear model
183  Krw = Sw;
184  dKrwDSw = 1.0;
185 
186  // // Non-linear model
187  // Krw = (Sw)*(Sw);
188  // dKrwDSw = 2.0*(Sw)*(1.0);
189 
190  // // Non-linear model
191  // Krw = (Sw)*(Sw)*(Sw);
192  // dKrwDSw = 3.0*(Sw)*(Sw)*(1.0);
193 
194 }
195 
197 void TPZMultiphase::SWaterstar(REAL &Swstar, REAL &Po, REAL &Sw)
198 {
199 
200  REAL waterdensity,oildensity;
201  REAL dwaterdensitydp,doildensitydp;
202  REAL waterviscosity, oilviscosity;
203  REAL dwaterviscositydp, doilviscositydp;
204 // REAL SUser=0.98;
205 
206  RhoWater(Po,waterdensity,dwaterdensitydp);
207  RhoOil(Po,oildensity,doildensitydp);
208  WaterViscosity(Po,waterviscosity,dwaterviscositydp);
209  OilViscosity(Po,oilviscosity,doilviscositydp);
210 
211  // Zero model
212  Swstar = Sw;
213 
214  // linear model
215  Swstar = (waterviscosity*oildensity)/(oilviscosity*waterdensity+waterviscosity*oildensity);
216 
217 // if(Swstar <= SUser)
218 // {
219 // Swstar = SUser;
220 // }
221 
222 // // Non-linear model
223 // Swstar = (sqrt(oilviscosity*waterviscosity*waterdensity*oildensity)-(waterviscosity*oildensity))/
224 // (oilviscosity*waterdensity-waterviscosity*oildensity);
225 
226 
227  // // Non-linear model cubic it is a large simbolic expression, so we got a problem here!!!!!
228  // Krw = (Sw)*(Sw)*(Sw);
229  // dKrwDSw = 3.0*(Sw)*(Sw)*(1.0);
230 
231 }
232 
234 void TPZMultiphase::Porosity(REAL po, REAL &poros, REAL &dPorosDpo){
235  const REAL comp = (0.0e-10)*(fPref);
236  const REAL pref = (1.0e6)/(fPref);
237  const REAL porosRef = 0.1;
238  poros = porosRef*exp(comp*(po-pref));
239  dPorosDpo = comp*porosRef*exp(comp*(po-pref));
240 }
241 
242 
244 void TPZMultiphase::RhoOil(REAL po, REAL &RhoOil, REAL &dRhoOilDpo){
245  const REAL Oilcomp = (0.0e-8)*(fPref);
246  const REAL pref = (1.0e6)/(fPref);
247  RhoOil = RhoOilSC()*exp(Oilcomp*(po-pref))/(fRhoref);
248  dRhoOilDpo = Oilcomp*RhoOilSC()*exp(Oilcomp*(po-pref))/(fRhoref);
249 }
250 
252 void TPZMultiphase::RhoWater(REAL po, REAL &RhoWater, REAL &dRhoWaterDpo){
253  const REAL Watercomp = (0.0e-9)*(fPref);
254  const REAL pref = (1.0e6)/(fPref);
255  RhoWater = RhoWaterSC()*exp(Watercomp*(po-pref))/(fRhoref);
256  dRhoWaterDpo = Watercomp*RhoWaterSC()*exp(Watercomp*(po-pref))/(fRhoref);
257 }
258 
259 
261 void TPZMultiphase::OilViscosity(REAL po, REAL &OilViscosity, REAL &dOilViscosityDpo){
262  const REAL OilViscRef = (1.0e-3)/(fEtaref);
263  OilViscosity = OilViscRef;
264  dOilViscosityDpo = 0;
265 }
266 
268 void TPZMultiphase::WaterViscosity(REAL po, REAL &WaterViscosity, REAL &dWaterViscosityDpo){
269  const REAL WaterViscRef = (1.0e-3)/(fEtaref);
270  WaterViscosity = WaterViscRef;
271  dWaterViscosityDpo = 0;
272 }
273 
275 void TPZMultiphase::OilLabmda(REAL &OilLabmda, REAL Po, REAL Sw, REAL &dOilLabmdaDPo, REAL &dOilLabmdaDSw){
276 
277  REAL krOil,Oilviscosity,OilDensity;
278  REAL dKroDSw,DOilviscosityDp,DOilDensityDp;
279 
280  Kro(Sw, krOil, dKroDSw);
281  OilViscosity(Po, Oilviscosity,DOilviscosityDp);
282  RhoOil(Po, OilDensity,DOilDensityDp);
283 
284  OilLabmda = ((krOil)/(Oilviscosity))*(OilDensity);
285  dOilLabmdaDPo = (DOilDensityDp/Oilviscosity)*(krOil);
286  dOilLabmdaDSw = (OilDensity/Oilviscosity)*(dKroDSw);
287 }
288 
289 
291 void TPZMultiphase::WaterLabmda(REAL &WaterLabmda, REAL Pw, REAL Sw, REAL &dWaterLabmdaDPw, REAL &dWaterLabmdaDSw){
292 
293  REAL krWater,Waterviscosity,WaterDensity;
294  REAL dKrwDSw,DWaterviscosityDp,DWaterDensityDp;
295 
296  Krw(Sw, krWater, dKrwDSw);
297  WaterViscosity(Pw, Waterviscosity,DWaterviscosityDp);
298  RhoWater(Pw, WaterDensity,DWaterDensityDp);
299 
300  WaterLabmda = ((krWater)/(Waterviscosity))*(WaterDensity);
301  dWaterLabmdaDPw = (DWaterDensityDp/Waterviscosity)*(krWater);
302  dWaterLabmdaDSw = (WaterDensity/Waterviscosity)*(dKrwDSw);
303 }
304 
305 
306 
308 void TPZMultiphase::Labmda(REAL &Labmda, REAL Pw, REAL Sw, REAL &dLabmdaDPw, REAL &dLabmdaDSw){
309 
310  REAL OilLabmdav;
311  REAL dOilLabmdaDpo,dOilLabmdaDSw;
312 
313  REAL WaterLabmdav;
314  REAL dWaterLabmdaDpo,dWaterLabmdaDSw;
315 
316  OilLabmda(OilLabmdav, Pw, Sw, dOilLabmdaDpo, dOilLabmdaDSw);
317  WaterLabmda(WaterLabmdav, Pw, Sw, dWaterLabmdaDpo, dWaterLabmdaDSw);
318 
319  Labmda = OilLabmdav + WaterLabmdav;
320  dLabmdaDPw = dOilLabmdaDpo+dWaterLabmdaDpo;
321  dLabmdaDSw = dOilLabmdaDSw+dWaterLabmdaDSw;
322 }
323 
325 void TPZMultiphase::fOil(REAL &fOil, REAL Pw, REAL Sw, REAL &dfOilDPw, REAL &dfOilDSw){
326 
327  REAL OilLabmdab;
328  REAL dOilLabmdaDpo,dOilLabmdaDSw;
329 
330  // REAL WaterLabmdav;
331  // REAL dWaterLabmdaDpo,dWaterLabmdaDSw;
332 
333  REAL Lambdab;
334  REAL dLabmdaDPw, dLabmdaDSw;
335 
336  OilLabmda(OilLabmdab, Pw, Sw, dOilLabmdaDpo, dOilLabmdaDSw);
337  // WaterLabmda(WaterLabmdav, Pw, Sw, dWaterLabmdaDpo, dWaterLabmdaDSw);
338  Labmda(Lambdab,Pw,Sw,dLabmdaDPw,dLabmdaDSw);
339 
340  fOil = ((OilLabmdab)/(Lambdab));
341  dfOilDPw = ((dOilLabmdaDpo)/(Lambdab))-((OilLabmdab)/(Lambdab*Lambdab))*dLabmdaDPw;
342  dfOilDSw = ((dOilLabmdaDSw)/(Lambdab))-((OilLabmdab)/(Lambdab*Lambdab))*dLabmdaDSw;
343 
344 }
345 
346 
348 void TPZMultiphase::fWater(REAL &fWater, REAL Pw, REAL Sw, REAL &dfWaterDPw, REAL &dfWaterDSw){
349 
350  // REAL OilLabmdab;
351  // REAL dOilLabmdaDpo,dOilLabmdaDSw;
352 
353  REAL WaterLabmdab;
354  REAL dWaterLabmdaDpo,dWaterLabmdaDSw;
355 
356  REAL Lambdab;
357  REAL dLabmdaDPw, dLabmdaDSw;
358 
359  // OilLabmda(OilLabmdab, Pw, Sw, dOilLabmdaDpo, dOilLabmdaDSw);
360  WaterLabmda(WaterLabmdab, Pw, Sw, dWaterLabmdaDpo, dWaterLabmdaDSw);
361  Labmda(Lambdab,Pw,Sw,dLabmdaDPw,dLabmdaDSw);
362 
363  fWater = ((WaterLabmdab)/(Lambdab));
364  dfWaterDPw = ((dWaterLabmdaDpo)/(Lambdab))-((WaterLabmdab)/(Lambdab*Lambdab))*dLabmdaDPw;
365  dfWaterDSw = ((dWaterLabmdaDSw)/(Lambdab))-((WaterLabmdab)/(Lambdab*Lambdab))*dLabmdaDSw;
366 
367 }
368 
373 void TPZMultiphase::fWaterstar(REAL &fWstar, REAL Pw, REAL Sw, REAL &dfWstarDPw, REAL &dfWstarDSw)
374 {
375  REAL fwater, foil;
376  REAL dfwaterdP, dfoildP;
377  REAL dfwaterdSw, dfoildSw;
378 
379  REAL Swcutoff;
380  SWaterstar(Swcutoff,Pw,Sw);
381 
382  if(Sw <= Swcutoff)
383  {
384  fOil(foil,Pw,Sw,dfoildP,dfoildSw);
385  fWater(fwater,Pw,Sw,dfwaterdP,dfwaterdSw);
386  fWstar = foil*fwater;
387  dfWstarDPw = (dfoildP)*fwater + (dfwaterdP)*foil;
388  dfWstarDSw = (dfoildSw)*fwater + (dfwaterdSw)*foil;
389  }
390  else
391  {
392  fOil(foil,Pw,Swcutoff,dfoildP,dfoildSw);
393  fWater(fwater,Pw,Swcutoff,dfwaterdP,dfwaterdSw);
394  fWstar = foil*fwater;
395  dfWstarDPw = 0.0*((dfoildP)*fwater + (dfwaterdP)*foil);
396  dfWstarDSw = 0.0*((dfoildSw)*fwater + (dfwaterdSw)*foil);
397  }
398 
399 
400 }
401 
406 void TPZMultiphase::fOilstar(REAL &fOstar, REAL Pw, REAL Sw, REAL &dfOstarDPw, REAL &dfOstarDSw)
407 {
408  REAL fwater, foil;
409  REAL dfwaterdP, dfoildP;
410  REAL dfwaterdSw, dfoildSw;
411  REAL Swcutoff;
412  SWaterstar(Swcutoff,Pw,Sw);
413 
414 
415  if(Sw >= Swcutoff)
416  {
417  fOil(foil,Pw,Sw,dfoildP,dfoildSw);
418  fWater(fwater,Pw,Sw,dfwaterdP,dfwaterdSw);
419  fOstar = foil*fwater;
420  dfOstarDPw = (dfoildP)*fwater + (dfwaterdP)*foil;
421  dfOstarDSw = (dfoildSw)*fwater + (dfwaterdSw)*foil;
422  }
423  else
424  {
425  fOil(foil,Pw,Swcutoff,dfoildP,dfoildSw);
426  fWater(fwater,Pw,Swcutoff,dfwaterdP,dfwaterdSw);
427  fOstar = foil*fwater;
428  dfOstarDPw = 0.0*((dfoildP)*fwater + (dfwaterdP)*foil);
429  dfOstarDSw = 0.0*((dfoildSw)*fwater + (dfwaterdSw)*foil);
430  }
431 }
432 
437 void TPZMultiphase::fstar(REAL &fStar, REAL Pw, REAL Sw, REAL Gdotn, REAL &dfstarDPw, REAL &dfstarDSw)
438 {
439  REAL fwaterStar, foilStar;
440  REAL dfwaterStardP, dfoilStardP;
441  REAL dfwaterStardSw, dfoilStardSw;
442 
443  if(Gdotn > 0.0)
444  {
445  fWaterstar(fwaterStar,Pw,Sw,dfwaterStardP,dfwaterStardSw);
446  fStar = fwaterStar;
447  dfstarDPw = dfwaterStardP;
448  dfstarDSw = dfwaterStardSw;
449  }
450  else
451  {
452  fOilstar(foilStar,Pw,Sw,dfoilStardP,dfoilStardSw);
453  fStar = foilStar;
454  dfstarDPw = dfoilStardP;
455  dfstarDSw = dfoilStardSw;
456 // if(fabs(Gdotn) <= 1.0e-14)
457 // {
458 // fStar = 0.0;
459 // dfstarDPw = 0.0;
460 // dfstarDSw = 0.0;
461 // }
462 
463  }
464 
465 }
466 
467 
468 // Fad Methods ///////////////////////////////////////////////////////////////////////////////////////////////////////
469 
470 #ifdef _AUTODIFF
471 
473 void TPZMultiphase::CapillaryPressure(BFadREAL So, BFadREAL &pc){
474  pc = 0.0;
475 }
476 
478 void TPZMultiphase::Kro(BFadREAL Sw, BFadREAL &Kro){
479  Kro = 1.0-Sw.val();
480 }
481 
482 
484 void TPZMultiphase::Krw(BFadREAL Sw, BFadREAL &Krw){
485  Krw = Sw;
486 }
487 
488 
490 void TPZMultiphase::Porosity(BFadREAL po, BFadREAL &poros){
491  const REAL comp = 1.0e-10;
492  const REAL pref = 1.0e6;
493  const REAL porosRef = 0.3;
494  poros = porosRef*exp(comp*((po.val())-pref));
495 }
496 
497 
499 void TPZMultiphase::RhoOil(BFadREAL po, BFadREAL &RhoOil){
500  const REAL Oilcomp = 0.0e-8;
501  const REAL pref = 1.0e6;
502  RhoOil = RhoOilSC()*exp(Oilcomp*((po.val())-pref));
503 }
504 
506 void TPZMultiphase::RhoWater(BFadREAL po, BFadREAL &RhoWater){
507  const REAL Watercomp = 0.0e-9;
508  const REAL pref = 1.0e6;
509  RhoWater = RhoWaterSC()*exp(Watercomp*((po.val())-pref));
510 }
511 
512 
514 void TPZMultiphase::OilViscosity(BFadREAL po, BFadREAL &OilViscosity){
515  const REAL OilViscRef = 1.0e-3;
516  OilViscosity = OilViscRef;
517 }
518 
520 void TPZMultiphase::WaterViscosity(BFadREAL po, BFadREAL &WaterViscosity){
521  const REAL WaterViscRef = 1.0e-3;
522  WaterViscosity = WaterViscRef;
523 }
524 
526 void TPZMultiphase::OilLabmda(BFadREAL OilLabmda, BFadREAL Po, BFadREAL &Sw){
527 
528  BFadREAL krOil,Oilviscosity,OilDensity;
529 
530  Kro(Sw, krOil);
531  OilViscosity(Po,Oilviscosity);
532  RhoOil(Po, OilDensity);
533 
534  OilLabmda = ((krOil)/(Oilviscosity))*(OilDensity);
535 
536 }
537 
538 
540 void TPZMultiphase::WaterLabmda(BFadREAL WaterLabmda, BFadREAL Pw, BFadREAL &Sw){
541 
542  BFadREAL krWater,Waterviscosity,WaterDensity;
543 
544  Krw(Sw, krWater);
545  WaterViscosity(Pw,Waterviscosity);
546  RhoWater(Pw,WaterDensity);
547 
548  WaterLabmda = ((krWater)/(Waterviscosity))*(WaterDensity);
549 
550 }
551 
552 
554 void TPZMultiphase::Labmda(BFadREAL Labmda, BFadREAL Pw, BFadREAL &Sw){
555 
556  BFadREAL OilLabmdaf, WaterLabmdaf;
557 
558  OilLabmda(OilLabmdaf, Pw, Sw);
559  WaterLabmda(WaterLabmdaf, Pw, Sw);
560 
561  Labmda = OilLabmdaf + WaterLabmdaf;
562 
563 }
564 
565 
567 void TPZMultiphase::fOil(BFadREAL fOil, BFadREAL Pw, BFadREAL &Sw)
568 {
569 }
570 
571 
573 void TPZMultiphase::fWater(BFadREAL fWater, BFadREAL Pw, BFadREAL &Sw)
574 {
575 
576 }
577 
578 #endif
579 
580 // Fad Methods ///////////////////////////////////////////////////////////////////////////////////////////////////////
581 
582 
583 // Constant data
584 
586 REAL TPZMultiphase::RhoOilSC(){
587  return 1000.0;
588 }
589 
591 REAL TPZMultiphase::RhoWaterSC(){
592  return 1000.0;
593 }
594 
596 TPZFMatrix<REAL> TPZMultiphase::Gravity(){
597  TPZFMatrix<REAL> gravity(3,1,0.0);
598  gravity(0,0) = 0.0;
599  gravity(1,0) = 0.0;
600  gravity(2,0) = 0.0;
601  gravity *= (1.0/(fPref/(fLref*fRhoref)));
602  return gravity;
603 }
604 
606 void TPZMultiphase::K(TPZFMatrix<REAL> &Kab){
607  Kab.Resize(3,3);
608  Kab.Zero();
609  Kab(0,0) = 1.0e-13;
610  Kab(1,1) = 1.0e-13;
611  Kab(2,2) = 0.0e-13;
612  Kab *= 1.0/(fKref);
613 
614 }
615 
617 TPZFMatrix<REAL> TPZMultiphase::Kinv(TPZFMatrix<REAL> &Kab){
618  TPZFMatrix<REAL> Kinverse(3,3,0.0);
619 
620  REAL Constant = (-1.0*Kab(0,1)*Kab(1,0)+Kab(0,0)*Kab(1,1));
621 
622  Kinverse(0,0) = 1.0*Kab(1,1)/(Constant);
623  Kinverse(0,1) = - 1.0*Kab(0,1)/(Constant);
624  Kinverse(1,0) = - 1.0*Kab(1,0)/(Constant);
625  Kinverse(1,1) = 1.0*Kab(0,0)/(Constant);
626  return Kinverse;
627 }
628 
629 
631 void TPZMultiphase::LoadKMap(std::string MaptoRead)
632 {
633 
634  std::string FileName;
635  std::string stringTemp;
636  FileName = MaptoRead;
637 
638  // Definitions
639  int64_t numelements=0;
640  int64_t numKData=0;
641 
642  // Scanning for total Number geometric elements
643  int64_t NumEntitiestoRead;
644  TPZStack <std::string> SentinelString;
645  {
646 
647  // reading a general mesh information by filter
648  std::ifstream read (FileName.c_str());
649  std::string FlagString;
650  // int flag = -1;
651  while(read)
652  {
653  char buf[1024];
654  read.getline(buf, 1024);
655  std::string str(buf);
656  if(str == "KABSOLUTE" || str == "KABSOLUTE\r")
657  {
658  SentinelString.push_back(str);
659  }
660 
661  }
662 
663  FlagString = "EndReading";
664  SentinelString.push_back(FlagString);
665  }
666 
667  NumEntitiestoRead = SentinelString.size();
668 
669  TPZStack <int> GeneralData(NumEntitiestoRead,0);
670  TPZStack <int> DataToProcess(NumEntitiestoRead,-1);
671 
672  {
673  // reading a general map information by filter
674  std::ifstream read (FileName.c_str());
675  std::string FlagString;
676  int64_t cont = 0;
677 
678  while(read)
679  {
680  char buf[1024];
681  read.getline(buf, 1024);
682  std::string str(buf);
683 
684  // Reading General Data
685  if(str == SentinelString[cont])
686  {
687  FlagString = str;
688  }
689  if(SentinelString[cont] == "" || SentinelString[cont] == "\r")
690  {
691  cont++;
692  }
693 
694  if(SentinelString[cont] == "EndReading")
695  {
696  break;
697  }
698 
699  if( (str != "" || str != "\r") && FlagString == SentinelString[cont])
700  {
701  // Data scaning
702  while (read) {
703  char buftemp[1024];
704  read.getline(buftemp, 1024);
705  std::string strtemp(buftemp);
706  GeneralData[cont]++;
707  if(strtemp == "" || strtemp == "\r")
708  {
709  FlagString = "";
710  GeneralData[cont]--;
711  cont++;
712  break;
713  }
714  }
715 
716  }
717 
718 
719  }
720  }
721 
722 
723  for (int64_t i = 0 ; i < NumEntitiestoRead; i++ )
724  {
725 
726  if(SentinelString[i] == "KABSOLUTE" || SentinelString[i] == "KABSOLUTE\r")
727  {
728  numKData=GeneralData[i];
729  DataToProcess[i]=0;
730  }
731  }
732 
733  numelements=numKData;
734  TPZFMatrix<REAL> Kabsolute(3,3);
735  Kabsolute.Resize(3,3);
736  Kabsolute.Zero();
737  // TPZStack<TPZFMatrix<REAL> > KabsoluteMap(numelements,0);
738  TPZStack<TPZFMatrix<REAL> > KabsoluteMap(100,0);
739 
740  int64_t elementId = 0;
741  int64_t ContOfKs = 0;
742 
743  REAL kxx , kxy, kxz;
744  REAL kyx , kyy, kyz;
745  REAL kzx , kzy, kzz;
746 
747 
748 
749  {
750 
751  // reading a general mesh information by filter
752  std::ifstream read (FileName.c_str());
753  std::string FlagString;
754  int64_t cont = 0;
755  // int dim = 0;
756  // int flag = 0;
757  while(read)
758  {
759  char buf[1024];
760  read.getline(buf, 1024);
761  std::string str(buf);
762  std::string strtemp="InitialState";
763 
764  // Reading General Data
765  if(str == SentinelString[cont])
766  {
767  FlagString = str;
768  }
769 
770  if(SentinelString[cont] == "" || SentinelString[cont] == "\r")
771  {
772  cont++;
773  }
774 
775  if(SentinelString[cont] == "EndReading")
776  {
777  break;
778  }
779 
780  if( (str != "" || str != "\r" )&& FlagString == SentinelString[cont])
781  {
782  // Data scaning
783  while (read)
784  {
785 
786  switch (DataToProcess[cont])
787  {
788  case 0:
789  {
790  //"KABSOLUTE"
791  if (GeneralData[cont] != 0)
792  {
793  read >> elementId;
794  read >> kxx;
795  read >> kxy;
796  read >> kxz;
797  read >> kyx;
798  read >> kyy;
799  read >> kyz;
800  read >> kzx;
801  read >> kzy;
802  read >> kzz;
803 
804  Kabsolute(0,0)=(1.0/fKref)*kxx;
805  Kabsolute(0,1)=(1.0/fKref)*kxy;
806  Kabsolute(0,2)=(1.0/fKref)*kxz;
807  Kabsolute(1,0)=(1.0/fKref)*kyx;
808  Kabsolute(1,1)=(1.0/fKref)*kyy;
809  Kabsolute(1,2)=(1.0/fKref)*kyz;
810  Kabsolute(2,0)=(1.0/fKref)*kzx;
811  Kabsolute(2,1)=(1.0/fKref)*kzy;
812  Kabsolute(2,2)=(1.0/fKref)*kzz;
813 
814  KabsoluteMap[elementId]=Kabsolute;
815 
816  ContOfKs++;
817  }
818  if(ContOfKs == numKData)
819  {
820  strtemp = "";
821  }
822  }
823  break;
824  default:
825  {
826  strtemp = "";
827  }
828  break;
829  }
830 
831  if(strtemp == "" || strtemp == "\r")
832  {
833  FlagString = "";
834  cont++;
835  break;
836  }
837  }
838 
839  }
840 
841 
842  }
843  }
844 
845 
846  SetKMap(KabsoluteMap);
847 
848 }
849 
850 
851 
852 // Contribute methods
853 
854 void TPZMultiphase::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef) {
855  DebugStop();
856 }
857 
858 
859 void TPZMultiphase::Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
860 {
861 
862 #ifdef PZDEBUG
863  int nref = datavec.size();
864  if (nref != 5 )
865  {
866  std::cout << " Error. datavec size is different from 5 \n";
867  DebugStop();
868  }
869 #endif
870 
871 
872  // Getting weight functions
873  TPZFMatrix<REAL> &phiU = datavec[0].phi;
874  TPZFMatrix<REAL> &phiQ = datavec[1].phi;
875  TPZFMatrix<REAL> &phiP = datavec[2].phi;
876  TPZFMatrix<REAL> &phiS = datavec[3].phi;
877 // TPZFMatrix<REAL> &phiQG = datavec[4].phi;
878 
879  TPZFMatrix<REAL> &dphiU = datavec[0].dphix;
880  // TPZFMatrix<REAL> &dphiQ = datavec[0].dphix;
881  TPZFMatrix<REAL> &dphiP = datavec[2].dphix;
882  TPZFMatrix<REAL> &dphiS = datavec[3].dphix;
883 
884  // number of test functions for each state variable
885  int phrU, phrQ, phrP, phrS, phrQG;
886  phrU = phiU.Rows();
887  phrQ = datavec[1].fVecShapeIndex.NElements();
888  phrP = phiP.Rows();
889  phrS = phiS.Rows();
890  phrQG = datavec[4].fVecShapeIndex.NElements();
891 
892  // blocks
893  int UStateVar = 2;
894  int FirstU = 0;
895  int FirstQ = phrU * UStateVar + FirstU;
896  int FirstP = phrQ + FirstQ;
897  int FirstS = phrP + FirstP;
898 // int FirstQG = FirstS;
899 
900  // Getting and computing another required data
901  REAL TimeStep = this->fDeltaT;
902  REAL Theta = this->fTheta;
903  REAL Gamma = this->fGamma;
904  int GeoID = datavec[1].gelElId;
905  TPZFMatrix<REAL> Kabsolute;
906  TPZFMatrix<REAL> Kinverse;
907  TPZFMatrix<REAL> Gfield;
908  if (fYorN)
909  {
910  Kabsolute=this->fKabsoluteMap[GeoID];
911  }
912  else
913  {
914  this->K(Kabsolute);
915  }
916 
917  Kinverse=this->Kinv(Kabsolute);
918  Gfield = this->Gravity();
919 
920  TPZManVector<STATE,3> sol_u = datavec[0].sol[0];
921  TPZManVector<STATE,3> sol_q = datavec[1].sol[0];
922  TPZManVector<STATE,3> sol_p = datavec[2].sol[0];
923  TPZManVector<STATE,3> sol_s = datavec[3].sol[0];
924  TPZManVector<STATE,3> sol_qg = datavec[4].sol[0];
925 
926  TPZFMatrix<STATE> dsol_u = datavec[0].dsol[0];
927  TPZFMatrix<STATE> dsol_q =datavec[1].dsol[0];
928  TPZFMatrix<STATE> dsol_p =datavec[2].dsol[0];
929  TPZFMatrix<STATE> dsol_s =datavec[3].dsol[0];
930 
931  REAL LambdaL, LambdaLU, MuL;
932  REAL Balpha, Sestr;
933 
934  REAL rockporosity, oildensity, waterdensity;
935  REAL drockporositydp, doildensitydp, dwaterdensitydp;
936 
937  REAL oilviscosity, waterviscosity;
938  REAL doilviscositydp, dwaterviscositydp;
939 
940  REAL bulklambda, oillambda, waterlambda;
941  REAL dbulklambdadp, doillambdadp, dwaterlambdadp;
942  REAL dbulklambdads, doillambdads, dwaterlambdads;
943 
944  REAL bulkfoil, bulkfwater;
945  REAL dbulkfoildp, dbulkfwaterdp;
946  REAL dbulkfoilds, dbulkfwaterds;
947 
948  // Functions computed at point x_{k} for each integration point
949  LambdaL = this->LameLambda();
950  LambdaLU = this->LameLambdaU();
951  MuL = this->LameMu();
952  Balpha = this->BiotAlpha();
953  Sestr = this->Se();
954 
955  REAL Pressure = sol_p[0];
956 
957  int VecPos= 0;
958  this->Porosity(sol_p[VecPos], rockporosity, drockporositydp);
959  this->RhoOil(sol_p[VecPos], oildensity, doildensitydp);
960  this->RhoWater(sol_p[VecPos], waterdensity, dwaterdensitydp);
961  this->OilViscosity(sol_p[VecPos], oilviscosity, doilviscositydp);
962  this->WaterViscosity(sol_p[VecPos], waterviscosity, dwaterviscositydp);
963  this->OilLabmda(oillambda, sol_p[VecPos], sol_s[VecPos], doillambdadp, doillambdads);
964  this->WaterLabmda(waterlambda, sol_p[VecPos], sol_s[VecPos], dwaterlambdadp, dwaterlambdads);
965  this->Labmda(bulklambda, sol_p[VecPos], sol_s[VecPos], dbulklambdadp, dbulklambdads);
966  this->fOil(bulkfoil, sol_p[VecPos], sol_s[VecPos], dbulkfoildp, dbulkfoilds);
967  this->fWater(bulkfwater, sol_p[VecPos], sol_s[VecPos], dbulkfwaterdp, dbulkfwaterds);
968 
969  // ////////////////////////// Jacobian Matrix ///////////////////////////////////
970  // Contribution of domain integrals for Jacobian matrix
971  // n+1 time step
972  if(gState == ECurrentState)
973  {
974 
975  // Elasticity Block (Equation for elasticity )
976  // Elastic equation
977  // Linear strain operator
978  // Ke Matrix
979  TPZFMatrix<REAL> du(2,2);
980  for(int iu = 0; iu < phrU; iu++ )
981  {
982  // Derivative for Vx
983  du(0,0) = dphiU(0,iu)*datavec[0].axes(0,0)+dphiU(1,iu)*datavec[0].axes(1,0);
984  // Derivative for Vy
985  du(1,0) = dphiU(0,iu)*datavec[0].axes(0,1)+dphiU(1,iu)*datavec[0].axes(1,1);
986 
987  for(int ju = 0; ju < phrU; ju++)
988  {
989  // Derivative for Ux
990  du(0,1) = dphiU(0,ju)*datavec[0].axes(0,0)+dphiU(1,ju)*datavec[0].axes(1,0);
991  // Derivative for Uy
992  du(1,1) = dphiU(0,ju)*datavec[0].axes(0,1)+dphiU(1,ju)*datavec[0].axes(1,1);
993 
994  if (this->fPlaneStress == 1)
995  {
996  /* Plain stress state */
997  ek(2*iu + FirstU, 2*ju + FirstU) += weight*((4*(MuL)*(LambdaL+MuL)/(LambdaL+2*MuL))*du(0,0)*du(0,1) + (2*MuL)*du(1,0)*du(1,1));
998 
999  ek(2*iu + FirstU, 2*ju+1 + FirstU) += weight*((2*(MuL)*(LambdaL)/(LambdaL+2*MuL))*du(0,0)*du(1,1) + (2*MuL)*du(1,0)*du(0,1));
1000 
1001  ek(2*iu+1 + FirstU, 2*ju + FirstU) += weight*((2*(MuL)*(LambdaL)/(LambdaL+2*MuL))*du(1,0)*du(0,1) + (2*MuL)*du(0,0)*du(1,1));
1002 
1003  ek(2*iu+1 + FirstU, 2*ju+1 + FirstU) += weight*((4*(MuL)*(LambdaL+MuL)/(LambdaL+2*MuL))*du(1,0)*du(1,1) + (2*MuL)*du(0,0)*du(0,1));
1004  }
1005  else
1006  {
1007  /* Plain Strain State */
1008  ek(2*iu + FirstU,2*ju + FirstU) += weight* ((LambdaL + 2*MuL)*du(0,0)*du(0,1) + (MuL)*du(1,0)*du(1,1));
1009 
1010  ek(2*iu + FirstU,2*ju+1 + FirstU) += weight* (LambdaL*du(0,0)*du(1,1) + (MuL)*du(1,0)*du(0,1));
1011 
1012  ek(2*iu+1 + FirstU,2*ju + FirstU) += weight* (LambdaL*du(1,0)*du(0,1) + (MuL)*du(0,0)*du(1,1));
1013 
1014  ek(2*iu+1 + FirstU,2*ju+1 + FirstU) += weight* ((LambdaL + 2*MuL)*du(1,0)*du(1,1) + (MuL)*du(0,0)*du(0,1));
1015 
1016  }
1017  }
1018  }
1019 
1020  // Matrix Qc
1021  // Coupling matrix for Elastic equation
1022  for(int in = 0; in < phrU; in++ )
1023  {
1024  du(0,0) = dphiU(0,in)*datavec[0].axes(0,0)+dphiU(1,in)*datavec[0].axes(1,0);
1025  du(1,0) = dphiU(0,in)*datavec[0].axes(0,1)+dphiU(1,in)*datavec[0].axes(1,1);
1026 
1027  for(int jp = 0; jp < phrP; jp++)
1028  {
1029  ek(2*in + FirstU,jp + FirstP) += (-1.0)*Balpha*weight*(phiP(jp,0)*du(0,0));
1030  ek(2*in+1 + FirstU,jp + FirstP) += (-1.0)*Balpha*weight*(phiP(jp,0)*du(1,0));
1031  }
1032  }
1033 
1034 
1035  REAL SaturationAtnplusOne = sol_s[0];
1036  // First Block (Equation One) constitutive law
1037  // Integrate[(Kinv/bulklambda)*dot(v,v), Omega_{e} ] (Equation One)
1038  REAL OneOverLambda = 1.0/bulklambda;
1039  for(int iq=0; iq<phrQ; iq++)
1040  {
1041 
1042  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1043  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1044  for (int jq=0; jq<phrQ; jq++)
1045  {
1046 
1047  int jvectorindex = datavec[1].fVecShapeIndex[jq].first;
1048  int jshapeindex = datavec[1].fVecShapeIndex[jq].second;
1049 
1050  if (fnewWS)
1051  {
1052  REAL vec1 = (Kinverse(0,0)*datavec[1].fNormalVec(0,ivectorindex)+Kinverse(0,1)*datavec[1].fNormalVec(1,ivectorindex));
1053  REAL vec2 = (Kinverse(1,0)*datavec[1].fNormalVec(0,ivectorindex)+Kinverse(1,1)*datavec[1].fNormalVec(1,ivectorindex));
1054 
1055  REAL dotprod =
1056  (phiQ(ishapeindex,0)*vec1) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(0,jvectorindex)) +
1057  (phiQ(ishapeindex,0)*vec2) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(1,jvectorindex)) ; // dot(K q,v)
1058 
1059  ek(iq + FirstQ,jq + FirstQ) += weight * OneOverLambda * dotprod;
1060  }
1061  else
1062  {
1063  REAL dotprod =
1064  (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex)) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(0,jvectorindex)) +
1065  (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex)) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(1,jvectorindex)) ; // dot(q,v)
1066 
1067  ek(iq + FirstQ,jq + FirstQ) += weight * OneOverLambda * dotprod;
1068  }
1069 
1070  }
1071 
1072  }
1073 
1074 
1075  // First Block (Equation One) constitutive law
1076  // Integrate[(d(1/bulklambdal)/dS)*dot(q,v), Omega_{e} ] (Equation One)
1077  /*REAL OneOverLambda = 1/bulklambda;*/
1078  for(int iq=0; iq<phrQ; iq++)
1079  {
1080 
1081  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1082  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1083  for (int jsat=0; jsat<phrS; jsat++)
1084  {
1085  if (fnewWS)
1086  {
1087  REAL dotprod =
1088  (Kinverse(0,0)*sol_q[0]+Kinverse(0,1)*sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex)) +
1089  (Kinverse(1,0)*sol_q[0]+Kinverse(1,1)*sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1090 
1091  ek(iq + FirstQ,jsat + FirstS) -= weight * dbulklambdads * OneOverLambda * OneOverLambda * phiS(jsat,0) * dotprod;
1092  }
1093  else
1094  {
1095  REAL dotprod =
1096  (sol_q[0]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex)) +
1097  (sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex)) ;
1098 
1099  ek(iq + FirstQ,jsat + FirstS) -= weight * dbulklambdads * OneOverLambda * OneOverLambda * phiS(jsat,0) * dotprod;
1100  }
1101 
1102  }
1103 
1104  }
1105 
1106 
1107  // First Block (Equation One) constitutive law
1108  // Integrate[(d(1/bulklambdal)/dP)*dot(q,v), Omega_{e} ] (Equation One)
1109  for(int iq=0; iq<phrQ; iq++)
1110  {
1111 
1112  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1113  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1114  for (int jp=0; jp<phrP; jp++)
1115  {
1116 
1117  if (fnewWS)
1118  {
1119  REAL dotprod =
1120  (Kinverse(0,0)*sol_q[0]+Kinverse(0,1)*sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex)) +
1121  (Kinverse(1,0)*sol_q[0]+Kinverse(1,1)*sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex)) ;
1122 
1123  ek(iq + FirstQ,jp + FirstP) -= weight * dbulklambdadp * OneOverLambda * OneOverLambda * phiP(jp,0) * dotprod;
1124  }
1125  else
1126  {
1127  REAL dotprod =
1128  (sol_q[0]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex)) +
1129  (sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex)) ;
1130 
1131  ek(iq + FirstQ,jp + FirstP) -= weight * dbulklambdadp * OneOverLambda * OneOverLambda * phiP(jp,0) * dotprod;
1132  }
1133  }
1134  }
1135 
1136 
1137  // This block was verified
1138  // First Block (Equation One) constitutive law
1139  // Integrate [ K dot(v,grad(P)) , Omega_{e}] (Equation One)
1140  for(int iq=0; iq<phrQ; iq++)
1141  {
1142 
1143  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1144  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1145 
1146  for (int jp=0; jp<phrP; jp++)
1147  {
1148 
1149 
1150  // Compute grad(W)
1151  TPZManVector<STATE> dsolp(2,0);
1152  dsolp[0] = dphiP(0,jp)*datavec[1].axes(0,0)+dphiP(1,jp)*datavec[1].axes(1,0);
1153  dsolp[1] = dphiP(0,jp)*datavec[1].axes(0,1)+dphiP(1,jp)*datavec[1].axes(1,1);
1154 
1155  if (fnewWS)
1156  {
1157 
1158  REAL e1e1 = (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex))*(dsolp[0]);
1159  REAL e2e2 = (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex))*(dsolp[1]);
1160 
1161  ek(iq + FirstQ,jp + FirstP) += weight * ( e1e1 + e2e2 );
1162 
1163  }
1164  else
1165  {
1166  REAL e1e1 = (Kabsolute(0,0)*(dsolp[0])+
1167  Kabsolute(0,1)*(dsolp[1]))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1168 
1169  REAL e2e2 = (Kabsolute(1,0)*(dsolp[0])+
1170  Kabsolute(1,1)*(dsolp[1]))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1171 
1172  ek(iq + FirstQ,jp + FirstP) += weight * ( e1e1 + e2e2 );
1173  }
1174  }
1175 
1176  }
1177 
1178  // This block was verified
1179  // First Block (Equation One) constitutive law
1180  // Integrate [ K (dot(f1*rho1*grad(g*z),v)+dot(f1*rho1*grad(g*z),v)) , Omega_{e}] (Equation One)
1181  // dS/dPalpha;
1182  for(int iq=0; iq<phrQ; iq++)
1183  {
1184 
1185  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1186  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1187 
1188  if (fnewWS)
1189  {
1190 
1191  REAL e1e1 = (Gfield(0,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1192  REAL e2e2 = (Gfield(1,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1193 
1194  for (int jp=0; jp<phrP; jp++)
1195  {
1196  ek(iq + FirstQ,jp + FirstP) -= weight * ( ( bulkfwater * dwaterdensitydp + bulkfoil * doildensitydp ) * phiP(jp,0) * (e1e1 + e2e2));
1197  }
1198 
1199  }
1200  else
1201  {
1202  REAL e1e1 = (Kabsolute(0,0)*(Gfield(0,0))+
1203  Kabsolute(0,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1204 
1205  REAL e2e2 = (Kabsolute(1,0)*(Gfield(0,0))+
1206  Kabsolute(1,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1207 
1208  for (int jp=0; jp<phrP; jp++)
1209  {
1210  ek(iq + FirstQ,jp + FirstP) -= weight * ( ( bulkfwater * dwaterdensitydp + bulkfoil * doildensitydp ) * phiP(jp,0) * (e1e1 + e2e2));
1211  }
1212  }
1213 
1214  }
1215 
1216  // This block was verified
1217  // First Block (Equation One) constitutive law
1218  // Integrate [ K (dot(f1*rho1*grad(g*z),v)+dot(f1*rho1*grad(g*z),v)) , Omega_{e}] (Equation One)
1219  // dP/dPalpha;
1220  for(int iq=0; iq<phrQ; iq++)
1221  {
1222 
1223  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1224  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1225 
1226  if (fnewWS)
1227  {
1228 
1229  REAL e1e1 = (Gfield(0,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1230  REAL e2e2 = (Gfield(1,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1231 
1232  for (int jp=0; jp<phrP; jp++)
1233  {
1234  ek(iq + FirstQ,jp + FirstP) -= weight * ( ( dbulkfwaterdp * waterdensity + dbulkfoildp * oildensity ) * phiP(jp,0) * (e1e1 + e2e2));
1235  }
1236 
1237  }
1238  else
1239  {
1240  REAL e1e1 = (Kabsolute(0,0)*(Gfield(0,0))+
1241  Kabsolute(0,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1242 
1243  REAL e2e2 = (Kabsolute(1,0)*(Gfield(0,0))+
1244  Kabsolute(1,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1245 
1246  for (int jp=0; jp<phrP; jp++)
1247  {
1248  ek(iq + FirstQ,jp + FirstP) -= weight * ( ( dbulkfwaterdp * waterdensity + dbulkfoildp * oildensity ) * phiP(jp,0) * (e1e1 + e2e2));
1249  }
1250  }
1251 
1252  }
1253 
1254  // This block was verified
1255  // First Block (Equation One) constitutive law
1256  // Integrate [ K (dot(f1*rho1*grad(g*z),v)+dot(f1*rho1*grad(g*z),v)) , Omega_{e}] (Equation One)
1257  // dS/dSalpha;
1258  for(int iq=0; iq < phrQ; iq++)
1259  {
1260 
1261  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1262  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1263 
1264  if (fnewWS)
1265  {
1266 
1267  REAL e1e1 = (Gfield(0,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1268  REAL e2e2 = (Gfield(1,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1269 
1270  for (int jsat=0; jsat < phrS; jsat++)
1271  {
1272  ek(iq+ FirstQ,jsat + FirstS) -= weight * ( ( dbulkfwaterds * waterdensity + dbulkfoilds * oildensity ) * phiS(jsat,0) * (e1e1 + e2e2));
1273  }
1274  }
1275  else
1276  {
1277  REAL e1e1 = (Kabsolute(0,0)*(Gfield(0,0))+
1278  Kabsolute(0,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1279 
1280  REAL e2e2 = (Kabsolute(1,0)*(Gfield(0,0))+
1281  Kabsolute(1,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1282 
1283  for (int jsat=0; jsat<phrS; jsat++)
1284  {
1285  ek(iq + FirstQ,jsat + FirstS) -= weight * ( ( dbulkfwaterds * waterdensity + dbulkfoilds * oildensity ) * phiS(jsat,0) * (e1e1 + e2e2));
1286  }
1287  }
1288 
1289  }
1290 
1291  // Poroelastic Contribution
1292  REAL divphiU, divU, dsolU[2][2];
1293 
1294  dsolU[0][0] = dsol_u(0,0)*datavec[0].axes(0,0)+dsol_u(1,0)*datavec[0].axes(1,0); // dUx/dx
1295  dsolU[1][0] = dsol_u(0,0)*datavec[0].axes(0,1)+dsol_u(1,0)*datavec[0].axes(1,1); // dUx/dy
1296 
1297  dsolU[0][1] = dsol_u(0,1)*datavec[0].axes(0,0)+dsol_u(1,1)*datavec[0].axes(1,0); // dUy/dx
1298  dsolU[1][1] = dsol_u(0,1)*datavec[0].axes(0,1)+dsol_u(1,1)*datavec[0].axes(1,1); // dUy/dy
1299  divU = dsolU[0][0]+dsolU[1][1]+0.0;
1300 
1301  REAL PhiStar = Balpha * divU + Sestr * Pressure;
1302 
1303  // Second Block (Equation Two) Bulk flux equation
1304  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1305  // d(porosity)/dPalpha
1306  for(int ip=0; ip < phrP; ip++)
1307  {
1308  for (int ju = 0; ju < phrU; ju++)
1309  {
1310  TPZFMatrix<REAL> du(2,2);
1311 
1312  du(0,0) = dphiU(0,ju)*datavec[0].axes(0,0); // dPhiU/dx dalphax
1313  du(1,0) = dphiU(0,ju)*datavec[0].axes(0,1); // dPhiU/dy dalphax
1314 
1315  du(0,1) = dphiU(1,ju)*datavec[0].axes(1,0); // dPhiU/dx dalphay
1316  du(1,1) = dphiU(1,ju)*datavec[0].axes(1,1); // dPhiU/dy dalphay
1317 
1318  divphiU = du(0,0)+du(1,0)+0.0;
1319 
1320  REAL dPhiStardUx = Balpha * (du(0,0)+du(1,0));
1321  REAL dPhiStardUy = Balpha * (du(0,1)+du(1,1));
1322 
1323  REAL Integratingx = phiP(ip,0) * dPhiStardUx * (waterdensity * (SaturationAtnplusOne) + oildensity * (1 - SaturationAtnplusOne));
1324  REAL Integratingy = phiP(ip,0) * dPhiStardUy * (waterdensity * (SaturationAtnplusOne) + oildensity * (1 - SaturationAtnplusOne));
1325 
1326  ek(ip + FirstP,2*ju + FirstU) += (-1.0) * weight * Integratingx;
1327  ek(ip + FirstP,2*ju + 1 + FirstU) += (-1.0) * weight * Integratingy;
1328  }
1329  }
1330 
1331  // Second Block (Equation Two) Bulk flux equation
1332  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1333  // d(porosity)/dPalpha
1334  for(int ip=0; ip < phrP; ip++)
1335  {
1336  for (int jp=0; jp < phrP; jp++)
1337  {
1338  REAL dPhiStardP = Sestr * phiP(jp,0);
1339  REAL Integrating = phiP(ip,0) * dPhiStardP * (waterdensity * (SaturationAtnplusOne) + oildensity * (1 - SaturationAtnplusOne));
1340  ek(ip + FirstP,jp + FirstP) += (-1.0) * weight * Integrating;
1341  }
1342  }
1343 
1344  // Second Block (Equation Two) Bulk flux equation
1345  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1346  // d(porosity)/dPalpha
1347  for(int ip=0; ip < phrP; ip++)
1348  {
1349  for (int jp=0; jp < phrP; jp++)
1350  {
1351  REAL Integrating = phiP(ip,0) * PhiStar * (dwaterdensitydp * (SaturationAtnplusOne) + doildensitydp * (1 - SaturationAtnplusOne)) * phiP(jp,0);
1352  ek(ip + FirstP,jp + FirstP) -= weight * Integrating;
1353  }
1354  }
1355 
1356  // Second Block (Equation Two) Bulk flux equation
1357  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1358  // d(porosity)/dPalpha
1359  for(int ip=0; ip < phrP; ip++)
1360  {
1361  for (int jp=0; jp < phrP; jp++)
1362  {
1363  REAL Integrating = phiP(ip,0) * drockporositydp * phiP(jp,0) * (waterdensity * (SaturationAtnplusOne) + oildensity * (1 - SaturationAtnplusOne));
1364  ek(ip + FirstP,jp + FirstP) += (-1.0) * weight * Integrating;
1365  }
1366  }
1367 
1368 
1369  // Second Block (Equation Two) Bulk flux equation
1370  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1371  // d(porosity)/dPalpha
1372  for(int ip=0; ip < phrP; ip++)
1373  {
1374  for (int jp=0; jp < phrP; jp++)
1375  {
1376  REAL Integrating = phiP(ip,0) * rockporosity * (dwaterdensitydp * (SaturationAtnplusOne) + doildensitydp * (1 - SaturationAtnplusOne)) * phiP(jp,0);
1377  ek(ip + FirstP,jp + FirstP) -= weight * Integrating;
1378  }
1379  }
1380 
1381 
1382  // Second Block (Equation Two) Bulk flux equation
1383  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1384  // d(S1)/dSalpha
1385  for(int ip=0; ip<phrP; ip++)
1386  {
1387  for (int jsat=0; jsat<phrS; jsat++)
1388  {
1389  REAL Integrating = phiP(ip,0) * PhiStar * (waterdensity - oildensity) * phiS(jsat,0);
1390  ek(ip + FirstP,jsat + FirstS) -= weight * Integrating;
1391  }
1392 
1393  }
1394 
1395  // Second Block (Equation Two) Bulk flux equation
1396  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1397  // d(S1)/dSalpha
1398  for(int ip=0; ip<phrP; ip++)
1399  {
1400  for (int jsat=0; jsat<phrS; jsat++)
1401  {
1402  REAL Integrating = phiP(ip,0) * rockporosity * (waterdensity - oildensity) * phiS(jsat,0);
1403  ek(ip + FirstP,jsat + FirstS) -= weight * Integrating;
1404  }
1405 
1406  }
1407 
1408  // Second Block (Equation Two) Bulk flux equation
1409  // Integrate[dot(grad(w),v), Omega_{e}] (Equation Two)
1410  for(int ip=0; ip<phrP; ip++)
1411  {
1412  // Compute grad(W)
1413  TPZManVector<STATE> dsolp(2,0);
1414  dsolp[0] = dphiP(0,ip)*datavec[2].axes(0,0)+dphiP(1,ip)*datavec[2].axes(1,0);
1415  dsolp[1] = dphiP(0,ip)*datavec[2].axes(0,1)+dphiP(1,ip)*datavec[2].axes(1,1);
1416 
1417  for (int jq=0; jq<phrQ; jq++)
1418  {
1419 
1420  int jvectorindex = datavec[1].fVecShapeIndex[jq].first;
1421  int jshapeindex = datavec[1].fVecShapeIndex[jq].second;
1422 
1423  REAL dotprod =
1424  (dsolp[0]) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(0,jvectorindex)) +
1425  (dsolp[1]) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(1,jvectorindex)) ;
1426 
1427  ek(ip + FirstP,jq + FirstQ) += (Gamma) * (TimeStep) * weight * dotprod;
1428  }
1429 
1430  }
1431 
1432  // Third Vector Block (Equation three)
1433  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1434  for(int isat=0; isat<phrS; isat++)
1435  {
1436 
1437  for (int jsat=0; jsat<phrS; jsat++)
1438  {
1439  ek(isat + FirstS,jsat+ FirstS) += (PhiStar * waterdensity) * weight * phiS(isat,0) * phiS(jsat,0);
1440  }
1441 
1442  }
1443 
1444  // Third Vector Block (Equation three)
1445  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1446  for(int isat=0; isat<phrS; isat++)
1447  {
1448 
1449  for (int jsat=0; jsat<phrS; jsat++)
1450  {
1451  ek(isat + FirstS,jsat+ FirstS) += (rockporosity * waterdensity) * weight * phiS(isat,0) * phiS(jsat,0);
1452  }
1453 
1454  }
1455 
1456  // Poroelastic contribution
1457  // Third Vector Block (Equation three)
1458  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1459  for(int isat=0; isat<phrS; isat++)
1460  {
1461  for (int ju=0; ju<phrU; ju++)
1462  {
1463  TPZFMatrix<REAL> du(2,2);
1464 
1465  du(0,0) = dphiU(0,ju)*datavec[0].axes(0,0); // dPhiU/dx dalphax
1466  du(1,0) = dphiU(0,ju)*datavec[0].axes(0,1); // dPhiU/dy dalphax
1467 
1468  du(0,1) = dphiU(1,ju)*datavec[0].axes(1,0); // dPhiU/dx dalphay
1469  du(1,1) = dphiU(1,ju)*datavec[0].axes(1,1); // dPhiU/dy dalphay
1470 
1471  divphiU = du(0,0)+du(1,0)+0.0;
1472 
1473  REAL dPhiStardUx = Balpha * (du(0,0)+du(1,0));
1474  REAL dPhiStardUy = Balpha * (du(0,1)+du(1,1));
1475 
1476  ek(isat + FirstS,2*ju + FirstU) += weight * (dPhiStardUx * waterdensity) * SaturationAtnplusOne * phiS(isat,0);
1477  ek(isat + FirstS,2*ju + 1 + FirstU) += weight * (dPhiStardUy * waterdensity) * SaturationAtnplusOne * phiS(isat,0);
1478  }
1479  }
1480 
1481  // Poroelastic contribution
1482  // Third Vector Block (Equation three)
1483  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1484  for(int isat=0; isat<phrS; isat++)
1485  {
1486  for (int jp=0; jp<phrP; jp++)
1487  {
1488  REAL dPhiStardP = Sestr * phiP(jp,0);
1489  ek(isat + FirstS,jp + FirstP) += weight * (dPhiStardP * waterdensity + PhiStar * dwaterdensitydp * phiP(jp,0)) * SaturationAtnplusOne * phiS(isat,0);
1490  }
1491  }
1492 
1493  // Third Vector Block (Equation three)
1494  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1495  for(int isat=0; isat<phrS; isat++)
1496  {
1497 
1498  for (int jp=0; jp<phrP; jp++)
1499  {
1500  ek(isat + FirstS,jp + FirstP) += weight * (drockporositydp * waterdensity + rockporosity * dwaterdensitydp) * SaturationAtnplusOne * phiS(isat,0) * phiP(jp,0);
1501  }
1502 
1503  }
1504 
1505  // Third Vector Block (Equation three)
1506  // Integrate[dot(d(f1(S1)/ds)q,grad(L)), Omega_{e}] (Equation three)
1507  for(int isat=0; isat<phrS; isat++)
1508  {
1509 
1510  // Compute grad(L)
1511  TPZManVector<STATE> Gradphis(2,0);
1512  Gradphis[0] = dphiS(0,isat)*datavec[3].axes(0,0)+dphiS(1,isat)*datavec[3].axes(1,0);
1513  Gradphis[1] = dphiS(0,isat)*datavec[3].axes(0,1)+dphiS(1,isat)*datavec[3].axes(1,1);
1514 
1515 
1516  REAL dotprod =
1517  (Gradphis[0]) * (sol_q[0]) +
1518  (Gradphis[1]) * (sol_q[1]);
1519 
1520  for (int jsat=0; jsat<phrS; jsat++)
1521  {
1522  ek(isat + FirstS,jsat + FirstS) -= (Theta) * (TimeStep) * weight * dbulkfwaterds * phiS(jsat,0) * dotprod;
1523  }
1524  }
1525 
1526  // Third Vector Block (Equation three)
1527  // Integrate[dot(f1(S1)v,grad(L)), Omega_{e}] (Equation three)
1528  for(int isat=0; isat<phrS; isat++)
1529  {
1530  // Compute grad(L)
1531  TPZManVector<STATE> Gradphis(2,0);
1532  Gradphis[0] = dphiS(0,isat)*datavec[3].axes(0,0)+dphiS(1,isat)*datavec[3].axes(1,0);
1533  Gradphis[1] = dphiS(0,isat)*datavec[3].axes(0,1)+dphiS(1,isat)*datavec[3].axes(1,1);
1534 
1535  for (int jq=0; jq<phrQ; jq++)
1536  {
1537 
1538  int jvectorindex = datavec[1].fVecShapeIndex[jq].first;
1539  int jshapeindex = datavec[1].fVecShapeIndex[jq].second;
1540  REAL dotprod =
1541  (Gradphis[0]) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(0,jvectorindex)) +
1542  (Gradphis[1]) * (phiQ(jshapeindex,0)*datavec[1].fNormalVec(1,jvectorindex));
1543 
1544  ek(isat + FirstS,jq + FirstQ) -= (Theta) * (TimeStep) * weight * bulkfwater * dotprod;
1545 
1546  }
1547 
1548  }
1549 
1550  // ////////////////////////// Jacobian Matrix ///////////////////////////////////
1551  // End of contribution of domain integrals for Jacobian matrix
1552 
1553  }
1554 
1555 
1556  // n time step
1557  // This values are constant in Newton iteration
1558  if(gState == ELastState)
1559  {
1560 
1561  }
1562 
1563 
1564 
1565  // Compute Residual contribution ef for domain integrals
1566  this->Contribute(datavec, weight, ef);
1567 
1568 
1569 }
1570 
1571 // Residual vector contribution
1572 void TPZMultiphase::Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ef)
1573 {
1574  // Getting weight functions
1575  TPZFMatrix<REAL> &phiU = datavec[0].phi;
1576  TPZFMatrix<REAL> &phiQ = datavec[1].phi;
1577  TPZFMatrix<REAL> &phiP = datavec[2].phi;
1578  TPZFMatrix<REAL> &phiS = datavec[3].phi;
1579 // TPZFMatrix<REAL> &phiQG = datavec[4].phi;
1580 
1581  TPZFMatrix<REAL> &dphiU = datavec[0].dphix;
1582  // TPZFMatrix<REAL> &dphiQ = datavec[0].dphix;
1583  TPZFMatrix<REAL> &dphiP = datavec[2].dphix;
1584  TPZFMatrix<REAL> &dphiS = datavec[3].dphix;// This a null value since S is constant by element.
1585 
1586 
1587  // number of test functions for each state variable
1588  int phrU, phrQ, phrP, phrS, phrQG;
1589  phrU = phiU.Rows();
1590  phrQ = datavec[1].fVecShapeIndex.NElements(); // phiQ.Rows();
1591  phrP = phiP.Rows();
1592  phrS = phiS.Rows();
1593  phrQG = datavec[4].fVecShapeIndex.NElements();
1594 
1595 
1596  // blocks
1597  int UStateVar = 2;
1598  int FirstU = 0;
1599  int FirstQ = phrU * UStateVar + FirstU;
1600  int FirstP = phrQ + FirstQ;
1601  int FirstS = phrP + FirstP;
1602 
1603  // Getting and computing another required data
1604  REAL TimeStep = this->fDeltaT;
1605  REAL Theta = this->fTheta;
1606  REAL Gamma = this->fGamma;
1607  int GeoID = datavec[1].gelElId;
1608  TPZFMatrix<REAL> Kabsolute;
1609  TPZFMatrix<REAL> Kinverse;
1610  TPZFMatrix<REAL> Gfield;
1611  if (fYorN)
1612  {
1613  Kabsolute=this->fKabsoluteMap[GeoID];
1614  }
1615  else
1616  {
1617  this->K(Kabsolute);
1618  }
1619 
1620  Kinverse=this->Kinv(Kabsolute);
1621  Gfield = this->Gravity();
1622  TPZManVector<STATE,3> sol_u =datavec[0].sol[0];
1623  TPZManVector<STATE,3> sol_q =datavec[1].sol[0];
1624  TPZManVector<STATE,3> sol_p =datavec[2].sol[0];
1625  TPZManVector<STATE,3> sol_s =datavec[3].sol[0];
1626  TPZManVector<STATE,3> sol_qg =datavec[4].sol[0];
1627 
1628  TPZFMatrix<STATE> dsol_u = datavec[0].dsol[0];
1629  TPZFMatrix<STATE> dsol_q =datavec[1].dsol[0];
1630  TPZFMatrix<STATE> dsol_p =datavec[2].dsol[0];
1631  TPZFMatrix<STATE> dsol_s =datavec[3].dsol[0];
1632 
1633  TPZFMatrix<> axesQ, axesP;
1634  axesQ=datavec[1].axes;
1635 
1636  REAL Pressure = sol_p[0];
1637 
1638  REAL LambdaL, LambdaLU, MuL;
1639  REAL Balpha, Sestr;
1640 
1641  REAL rockporosity, oildensity, waterdensity;
1642  REAL drockporositydp, doildensitydp, dwaterdensitydp;
1643 
1644  REAL oilviscosity, waterviscosity;
1645  REAL doilviscositydp, dwaterviscositydp;
1646 
1647  REAL bulklambda, oillambda, waterlambda;
1648  REAL dbulklambdadp, doillambdadp, dwaterlambdadp;
1649  REAL dbulklambdads, doillambdads, dwaterlambdads;
1650 
1651  REAL bulkfoil, bulkfwater;
1652  REAL dbulkfoildp, dbulkfwaterdp;
1653  REAL dbulkfoilds, dbulkfwaterds;
1654 
1655  // Functions computed at point x_{k} for each integration point
1656  LambdaL = this->LameLambda();
1657  LambdaLU = this->LameLambdaU();
1658  MuL = this->LameMu();
1659  Balpha = this->BiotAlpha();
1660  Sestr = this->Se();
1661 
1662  int VecPos= 0;
1663  // REAL PressureRef = 1.0e6;
1664  this->Porosity(sol_p[VecPos], rockporosity, drockporositydp);
1665  this->RhoOil(sol_p[VecPos], oildensity, doildensitydp);
1666  this->RhoWater(sol_p[VecPos], waterdensity, dwaterdensitydp);
1667  this->OilViscosity(sol_p[VecPos], oilviscosity, doilviscositydp);
1668  this->WaterViscosity(sol_p[VecPos], waterviscosity, dwaterviscositydp);
1669  this->OilLabmda(oillambda, sol_p[VecPos], sol_s[VecPos], doillambdadp, doillambdads);
1670  this->WaterLabmda(waterlambda, sol_p[VecPos], sol_s[VecPos], dwaterlambdadp, dwaterlambdads);
1671  this->Labmda(bulklambda, sol_p[VecPos], sol_s[VecPos], dbulklambdadp, dbulklambdads);
1672  this->fOil(bulkfoil, sol_p[VecPos], sol_s[VecPos], dbulkfoildp, dbulkfoilds);
1673  this->fWater(bulkfwater, sol_p[VecPos], sol_s[VecPos], dbulkfwaterdp, dbulkfwaterds);
1674 
1675  // ////////////////////////// Residual Vector ///////////////////////////////////
1676  // Contribution of domain integrals for Residual Vector
1677 
1678  // n+1 time step
1679  if(gState == ECurrentState)
1680  {
1681  // Elastic equation
1682  // Linear strain operator
1683  // Ke Matrix
1684  TPZFMatrix<REAL> du(2,2);
1685  TPZFMatrix<REAL> dux(2,2);
1686  TPZFMatrix<REAL> duy(2,2);
1687  // Required check out of this implementation
1688  // Derivative for Ux
1689  dux(0,1) = dsol_u(0,0)*datavec[0].axes(0,0)+dsol_u(1,0)*datavec[0].axes(1,0); // dUx/dx
1690  dux(1,1) = dsol_u(0,0)*datavec[0].axes(0,1)+dsol_u(1,0)*datavec[0].axes(1,1); // dUx/dy
1691 
1692  // Derivative for Uy
1693  duy(0,1) = dsol_u(0,1)*datavec[0].axes(0,0)+dsol_u(1,1)*datavec[0].axes(1,0); // dUy/dx
1694  duy(1,1) = dsol_u(0,1)*datavec[0].axes(0,1)+dsol_u(1,1)*datavec[0].axes(1,1); // dUy/dy
1695 
1696  for(int iu = 0; iu < phrU; iu++ )
1697  {
1698  // Derivative for Vx
1699  du(0,0) = dphiU(0,iu)*datavec[0].axes(0,0)+dphiU(1,iu)*datavec[0].axes(1,0);
1700  // Derivative for Vy
1701  du(1,0) = dphiU(0,iu)*datavec[0].axes(0,1)+dphiU(1,iu)*datavec[0].axes(1,1);
1702 
1703  // Fu Vector Force right hand term check the gravity term
1704 // ef(2*iu + FirstU) += weight*Gfield(0,0)*phiU(iu, 0);
1705 // ef(2*iu+1 + FirstU) += weight*Gfield(1,0)*phiU(iu, 0);
1706 
1707  if (fPlaneStress == 1)
1708  {
1709  /* Plain stress state */
1710  ef(2*iu + FirstU) += weight*((4*(MuL)*(LambdaL+MuL)/(LambdaL+2*MuL))*du(0,0)*dux(0,1) + (2*MuL)*du(1,0)*dux(1,1));
1711 
1712  ef(2*iu + FirstU) += weight*((2*(MuL)*(LambdaL)/(LambdaL+2*MuL))*du(0,0)*duy(1,1) + (2*MuL)*du(1,0)*duy(0,1));
1713 
1714  ef(2*iu+1 + FirstU) += weight*((2*(MuL)*(LambdaL)/(LambdaL+2*MuL))*du(1,0)*dux(0,1) + (2*MuL)*du(0,0)*dux(1,1));
1715 
1716  ef(2*iu+1 + FirstU) += weight*((4*(MuL)*(LambdaL+MuL)/(LambdaL+2*MuL))*du(1,0)*duy(1,1) + (2*MuL)*du(0,0)*duy(0,1));
1717  }
1718  else
1719  {
1720 
1721 // /* Plain Strain State */
1722 // ek(2*iu + FirstU,2*ju + FirstU) += weight* ((LambdaL + 2*MuL)*du(0,0)*du(0,1) + (MuL)*du(1,0)*du(1,1));
1723 //
1724 // ek(2*iu + FirstU,2*ju+1 + FirstU) += weight* (LambdaL*du(0,0)*du(1,1) + (MuL)*du(1,0)*du(0,1));
1725 //
1726 // ek(2*iu+1 + FirstU,2*ju + FirstU) += weight* (LambdaL*du(1,0)*du(0,1) + (MuL)*du(0,0)*du(1,1));
1727 //
1728 // ek(2*iu+1 + FirstU,2*ju+1 + FirstU) += weight* ((LambdaL + 2*MuL)*du(1,0)*du(1,1) + (MuL)*du(0,0)*du(0,1));
1729 
1730  /* Plain Strain State */
1731  ef(2*iu + FirstU) += weight* ((LambdaL + 2*MuL)*du(0,0)*dux(0,1) + (MuL)*du(1,0)*(dux(1,1)));
1732 
1733  ef(2*iu + FirstU) += weight* (LambdaL*du(0,0)*duy(1,1) + (MuL)*du(1,0)*(duy(0,1)));
1734 
1735  ef(2*iu+1 + FirstU) += weight* (LambdaL*du(1,0)*dux(0,1) + (MuL)*du(0,0)*(dux(1,1)));
1736 
1737  ef(2*iu+1 + FirstU) += weight* ((LambdaL + 2*MuL)*du(1,0)*duy(1,1) + (MuL)*du(0,0)*(duy(0,1)));
1738  }
1739  }
1740 
1741  // Matrix Qc
1742  // Coupling matrix for Elastic equation
1743  for(int in = 0; in < phrU; in++ )
1744  {
1745  du(0,0) = dphiU(0,in)*datavec[0].axes(0,0)+dphiU(1,in)*datavec[0].axes(1,0);
1746  du(1,0) = dphiU(0,in)*datavec[0].axes(0,1)+dphiU(1,in)*datavec[0].axes(1,1);
1747 
1748  ef(2*in + FirstU) += (-1.0)*Balpha*weight*(Pressure*du(0,0));
1749  ef(2*in+1 + FirstU) += (-1.0)*Balpha*weight*(Pressure*du(1,0));
1750  }
1751 
1752  // This block was verified
1753  REAL SaturationAtnplusOne = sol_s[0];
1754  // First Block (Equation One) constitutive law
1755  // Integrate[(viscosity/density)*dot(q,v), Omega_{e}] (Equation One)
1756  // REAL ViscOverdensity = waterviscosity/waterdensity;
1757  REAL OneOverLambda = 1.0/bulklambda;
1758  for(int iq=0; iq<phrQ; iq++)
1759  {
1760 
1761  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1762  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1763 
1764  if (fnewWS)
1765  {
1766  REAL dotprod =
1767  (Kinverse(0,0)*sol_q[0]+Kinverse(0,1)*sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex)) +
1768  (Kinverse(1,0)*sol_q[0]+Kinverse(1,1)*sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex)) ;
1769 
1770  ef(iq + FirstQ) += OneOverLambda * weight * dotprod;
1771  }
1772  else
1773  {
1774  REAL dotprod =
1775  (sol_q[0]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex)) +
1776  (sol_q[1]) * (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex)) ;
1777 
1778  ef(iq + FirstQ) += OneOverLambda * weight * dotprod;
1779  }
1780 
1781  }
1782 
1783 
1784  // This block was verified
1785  // First Block (Equation One) constitutive law
1786  // Integrate [ K dot(v,grad(P)) , Omega_{e}] (Equation One)
1787  // Compute grad(P)
1788  TPZManVector<STATE> dsolp(2,0);
1789  dsolp[0] = dsol_p(0,0)*datavec[2].axes(0,0)+dsol_p(1,0)*datavec[2].axes(1,0);
1790  dsolp[1] = dsol_p(0,0)*datavec[2].axes(0,1)+dsol_p(1,0)*datavec[2].axes(1,1);
1791 
1792  for(int iq=0; iq<phrQ; iq++)
1793  {
1794 
1795  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1796  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1797 
1798  if (fnewWS)
1799  {
1800  REAL e1e1 = (phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex))*(dsolp[0]);
1801  REAL e2e2 = (phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex))*(dsolp[1]);
1802 
1803  ef(iq + FirstQ) += weight * (e1e1 + e2e2);
1804  }
1805  else
1806  {
1807  REAL e1e1 = (Kabsolute(0,0)*(dsolp[0])+
1808  Kabsolute(0,1)*(dsolp[1]))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1809 
1810  REAL e2e2 = (Kabsolute(1,0)*(dsolp[0])+
1811  Kabsolute(1,1)*(dsolp[1]))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1812 
1813  ef(iq + FirstQ) += weight * (e1e1 + e2e2);
1814  }
1815  }
1816 
1817 
1818  // This block was verified
1819  // First Block (Equation One) constitutive law
1820  // Integrate [ K (dot(f1*rho1*grad(g*z),v)+dot(f1*rho1*grad(g*z),v)) , Omega_{e}] (Equation One)
1821 
1822  for(int iq=0; iq<phrQ; iq++)
1823  {
1824 
1825  int ivectorindex = datavec[1].fVecShapeIndex[iq].first;
1826  int ishapeindex = datavec[1].fVecShapeIndex[iq].second;
1827 
1828  if (fnewWS)
1829  {
1830 
1831  REAL e1e1 = (Gfield(0,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1832  REAL e2e2 = (Gfield(1,0))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1833 
1834  ef(iq + FirstQ) -= weight * ( ( bulkfwater * waterdensity + bulkfoil * oildensity )* (e1e1 + e2e2));
1835  }
1836  else
1837  {
1838  REAL e1e1 = (Kabsolute(0,0)*(Gfield(0,0))+
1839  Kabsolute(0,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(0,ivectorindex));
1840 
1841  REAL e2e2 = (Kabsolute(1,0)*(Gfield(0,0))+
1842  Kabsolute(1,1)*(Gfield(1,0)))*(phiQ(ishapeindex,0)*datavec[1].fNormalVec(1,ivectorindex));
1843 
1844  ef(iq + FirstQ) -= weight * ((bulkfwater * waterdensity + bulkfoil * oildensity)* (e1e1 + e2e2));
1845  }
1846 
1847  }
1848 
1849  // Poroelastic Contribution
1850  REAL divU, dsolU[2][2];
1851  dsolU[0][0] = dsol_u(0,0)*datavec[0].axes(0,0)+dsol_u(1,0)*datavec[0].axes(1,0); // dUx/dx
1852  dsolU[1][0] = dsol_u(0,0)*datavec[0].axes(0,1)+dsol_u(1,0)*datavec[0].axes(1,1); // dUx/dy
1853 
1854  dsolU[0][1] = dsol_u(0,1)*datavec[0].axes(0,0)+dsol_u(1,1)*datavec[0].axes(1,0); // dUy/dx
1855  dsolU[1][1] = dsol_u(0,1)*datavec[0].axes(0,1)+dsol_u(1,1)*datavec[0].axes(1,1); // dUy/dy
1856 
1857  divU = dsolU[0][0]+dsolU[1][1]+0.0;
1858  REAL PhiStar = Balpha * divU + Sestr * Pressure;
1859 
1860  // Second Block (Equation Two) Bulk flux equation
1861  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1862  for(int ip=0; ip<phrP; ip++)
1863  {
1864  // Here rockporosity is the poroelastic contribution
1865  REAL Integrating = phiP(ip,0) * PhiStar * (waterdensity * SaturationAtnplusOne + oildensity * (1 - SaturationAtnplusOne));
1866  ef(ip + FirstP) += (-1.0) * weight * Integrating;
1867  }
1868 
1869  // Second Block (Equation Two) Bulk flux equation
1870  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1871  for(int ip=0; ip<phrP; ip++)
1872  {
1873  // Here rockporosity is the poroelastic contribution
1874  REAL Integrating = phiP(ip,0) * rockporosity * (waterdensity * SaturationAtnplusOne + oildensity * (1 - SaturationAtnplusOne));
1875  ef(ip + FirstP) += (-1.0) * weight * Integrating;
1876  }
1877 
1878  // Second Block (Equation Two) Bulk flux equation
1879  // Integrate[dot(grad(W),q), Omega_{e}] (Equation Two)
1880  // This block was verified
1881  for(int ip=0; ip<phrP; ip++)
1882  {
1883  // Compute grad(W)
1884  TPZManVector<STATE> dsolp(2,0);
1885  dsolp[0] = dphiP(0,ip)*datavec[2].axes(0,0)+dphiP(1,ip)*datavec[2].axes(1,0);
1886  dsolp[1] = dphiP(0,ip)*datavec[2].axes(0,1)+dphiP(1,ip)*datavec[2].axes(1,1);
1887 
1888  REAL dotprod =
1889  (dsolp[0]) * (sol_q[0]) +
1890  (dsolp[1]) * (sol_q[1]);
1891 
1892  ef(ip + FirstP) += (Gamma) * (TimeStep) * weight * dotprod;
1893  }
1894 
1895  // Poroelastic Contribution
1896  // Third Vector Block (Equation three)
1897  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1898  for(int isat=0; isat<phrS; isat++)
1899  {
1900  // Here rockporosity is the poroelastic contribution
1901  ef(isat + FirstS) += weight * (PhiStar * waterdensity) * phiS(isat,0) * SaturationAtnplusOne;
1902  }
1903 
1904 
1905  // This block was verified
1906  // Third Vector Block (Equation three)
1907  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1908  for(int isat=0; isat<phrS; isat++)
1909  {
1910  // Here rockporosity is the poroelastic contribution
1911  ef(isat + FirstS) += weight * (rockporosity * waterdensity) * phiS(isat,0) * SaturationAtnplusOne;
1912  }
1913 
1914 
1915  // Third Vector Block (Equation three)
1916  // Integrate[dot(f1(S1)q,grad(L)), Omega_{e}] (Equation three)
1917  // std::cout << "phrS: " << phrS << " \n" << std::endl;
1918  for(int isat=0; isat<phrS; isat++)
1919  {
1920  // Compute grad(L)
1921  TPZManVector<STATE> Gradphis(2,0);
1922  Gradphis[0] = dphiS(0,isat)*datavec[3].axes(0,0)+dphiS(1,isat)*datavec[3].axes(1,0);
1923  Gradphis[1] = dphiS(0,isat)*datavec[3].axes(0,1)+dphiS(1,isat)*datavec[3].axes(1,1);
1924 
1925  REAL dotprod =
1926  (Gradphis[0]) * (sol_q[0]) +
1927  (Gradphis[1]) * (sol_q[1]);
1928 
1929  ef(isat + FirstS) -= (Theta) * (TimeStep) * weight * bulkfwater * dotprod;
1930  }
1931 
1932  }
1933 
1934 
1935  // n time step
1936  if(gState == ELastState)
1937  {
1938  REAL SaturationAtnTimeStep = sol_s[0]; // Gettin Saturation at n time step
1939 
1940  // Poroelastic contribution
1941  REAL divU, dsolU[2][2];
1942  dsolU[0][0] = dsol_u(0,0)*datavec[0].axes(0,0)+dsol_u(1,0)*datavec[0].axes(1,0); // dUx/dx
1943  dsolU[1][0] = dsol_u(0,0)*datavec[0].axes(0,1)+dsol_u(1,0)*datavec[0].axes(1,1); // dUx/dy
1944 
1945  dsolU[0][1] = dsol_u(0,1)*datavec[0].axes(0,0)+dsol_u(1,1)*datavec[0].axes(1,0); // dUy/dx
1946  dsolU[1][1] = dsol_u(0,1)*datavec[0].axes(0,1)+dsol_u(1,1)*datavec[0].axes(1,1); // dUy/dy
1947 
1948  divU = dsolU[0][0]+dsolU[1][1]+0.0;
1949  REAL PhiStar = Balpha * divU + Sestr * Pressure;
1950 
1951  // Second Block (Equation Two) Bulk flux equation
1952  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1953  for(int ip = 0; ip < phrP; ip++)
1954  {
1955  // Here rockporosity is the poroelastic contribution
1956  REAL Integrating = phiP(ip,0) * PhiStar * (waterdensity * SaturationAtnTimeStep + oildensity * (1 - SaturationAtnTimeStep));
1957  ef(ip + FirstP) += weight * Integrating;
1958  }
1959 
1960  // Second Block (Equation Two) Bulk flux equation
1961  // Integrate[W*(d(\phi*(rho1 * S1 + rho2 * S2)/dt)), Omega_{e}] (Equation Two)
1962  for(int ip = 0; ip < phrP; ip++)
1963  {
1964  // Here rockporosity is the poroelastic contribution
1965  REAL Integrating = phiP(ip,0) * rockporosity * (waterdensity * SaturationAtnTimeStep + oildensity * (1 - SaturationAtnTimeStep));
1966  ef(ip + FirstP) += weight * Integrating;
1967  }
1968 
1969  // Second Block (Equation Two) Bulk flux equation
1970  // Integrate[dot(grad(w),q), Omega_{e}] (Equation Two)
1971  // This block was verified
1972  for(int ip=0; ip<phrP; ip++)
1973  {
1974  // Compute grad(W)
1975  TPZManVector<STATE> dsolp(2,0);
1976  dsolp[0] = dphiP(0,ip)*datavec[2].axes(0,0)+dphiP(1,ip)*datavec[2].axes(1,0);
1977  dsolp[1] = dphiP(0,ip)*datavec[2].axes(0,1)+dphiP(1,ip)*datavec[2].axes(1,1);
1978 
1979  REAL dotprod =
1980  (dsolp[0]) * (sol_q[0]) +
1981  (dsolp[1]) * (sol_q[1]);
1982 
1983  ef(ip + FirstP) += (1-Gamma) * (TimeStep) * weight * dotprod;
1984  }
1985 
1986 
1987  // Poroelastic Contribution
1988  // Third Vector Block (Equation three)
1989  // Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1990  for(int isat=0; isat<phrS; isat++)
1991  {
1992  // Here rockporosity is the poroelastic contribution
1993  ef(isat + FirstS) += (-1.0) * weight * (PhiStar * waterdensity) * phiS(isat,0) * SaturationAtnTimeStep;
1994  }
1995 
1996  // Third Vector Block (Equation three)
1997  // (-1.0) * Integrate[porosity*density*L*S, Omega_{e}] (Equation three)
1998  for(int isat=0; isat<phrS; isat++)
1999  {
2000  // Here rockporosity is the poroelastic contribution
2001  ef(isat + FirstS) += (-1.0) * weight * (rockporosity * waterdensity) * phiS(isat,0) * SaturationAtnTimeStep;
2002  }
2003 
2004  // Third Vector Block (Equation three)
2005  // Integrate[dot(f1(S1)q,grad(L)), Omega_{e}] (Equation three)
2006  for(int isat=0; isat<phrS; isat++)
2007  {
2008  // Compute grad(L)
2009  TPZManVector<STATE> Gradphis(2,0);
2010  Gradphis[0] = dphiS(0,isat)*datavec[3].axes(0,0)+dphiS(1,isat)*datavec[3].axes(1,0);
2011  Gradphis[1] = dphiS(0,isat)*datavec[3].axes(0,1)+dphiS(1,isat)*datavec[3].axes(1,1);
2012 
2013  REAL dotprod =
2014  (Gradphis[0]) * (sol_q[0]) +
2015  (Gradphis[1]) * (sol_q[1]) ;
2016 
2017  ef(isat + FirstS) -= (1-Theta) * (TimeStep) * weight * bulkfwater * dotprod;
2018  }
2019  }
2020 
2021  // ////////////////////////// Residual Vector ///////////////////////////////////
2022  // End of contribution of domain integrals for Residual Vector
2023 
2024 
2025 }
2026 
2027 
2028 void TPZMultiphase::ContributeInterface(TPZVec<TPZMaterialData> &datavec,TPZVec<TPZMaterialData> &dataleftvec,TPZVec<TPZMaterialData> &datarightvec,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
2029 {
2030  DebugStop();
2031 }
2032 
2033 void TPZMultiphase::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
2034 {
2035  DebugStop();
2036 }
2037 
2038 void TPZMultiphase::ContributeInterface(TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, TPZVec<TPZMaterialData> &dataright, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
2039 {
2040 
2041  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
2042  TPZFMatrix<REAL> &phiUR = dataright[0].phi;
2043 
2044  TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
2045  TPZFMatrix<REAL> &phiQR = dataright[1].phi;
2046 
2047  TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
2048  TPZFMatrix<REAL> &phiPR = dataright[2].phi;
2049 
2050  TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
2051  TPZFMatrix<REAL> &phiSR = dataright[3].phi;
2052 
2053  TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
2054 // TPZFMatrix<REAL> &phiQGR = dataright[4].phi;
2055 
2056 
2057  TPZManVector<REAL,3> &normal = data.normal;
2058 
2059  REAL n1 = normal[0];
2060  REAL n2 = normal[1];
2061  // REAL n3 = normal[2];
2062 
2063  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
2064  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
2065  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
2066  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
2067  TPZManVector<STATE,3> sol_qgL =dataleft[4].sol[0];
2068 
2069  TPZManVector<STATE,3> sol_uR =dataright[0].sol[0];
2070  TPZManVector<STATE,3> sol_qR =dataright[1].sol[0];
2071  TPZManVector<STATE,3> sol_pR =dataright[2].sol[0];
2072  TPZManVector<STATE,3> sol_sR =dataright[3].sol[0];
2073  TPZManVector<STATE,3> sol_qgR =dataright[4].sol[0];
2074 
2075  // Getting Q solution for left and right side
2076  REAL qxL = sol_qL[0];
2077  REAL qyL = sol_qL[1];
2078  REAL qxR = sol_qR[0];
2079  REAL qyR = sol_qR[1];
2080  REAL dotqnL = (qxL*n1) + (qyL*n2);
2081  REAL dotqnR = (qxR*n1) + (qyR*n2);
2082 
2083  // Getting QG solution for left and right side
2084 // REAL qgxL = sol_qgL[0];
2085 // REAL qgyL = sol_qgL[1];
2086 // REAL qgxR = sol_qgR[0];
2087 // REAL qgyR = sol_qgR[1];
2088 // REAL dotqgnL = (qgxL*n1) + (qgyL*n2);
2089 // REAL dotqgnR = (qgxR*n1) + (qgyR*n2);
2090 
2091  // Getting S solution for left and right side
2092  REAL SaturationL = sol_sL[0];
2093  REAL SaturationR = sol_sR[0];
2094 
2095  // Getting another required data
2096  REAL TimeStep = this->fDeltaT;
2097  REAL Theta = this->fTheta;
2098  REAL Gamma = this->fGamma;
2099  this->fnewWS=true;
2100 
2101  // Getting Harmonic mean of permeabilities
2102  TPZFMatrix<REAL> KabsoluteLeft;
2103  TPZFMatrix<REAL> KabsoluteRight;
2104  TPZFMatrix<REAL> Gfield;
2105  TPZFMatrix<REAL> Kmean(3,3,0);
2106 
2107  int GeoIDLeft = dataleft[1].gelElId;
2108  int GeoIDRight = dataright[1].gelElId;
2109 
2110  REAL rockporosityl, oildensityl, waterdensityl;
2111  REAL drockporositydpl, doildensitydpl, dwaterdensitydpl;
2112  REAL rockporosityr, oildensityr, waterdensityr;
2113  REAL drockporositydpr, doildensitydpr, dwaterdensitydpr;
2114 
2115  REAL oilviscosityl, waterviscosityl;
2116  REAL doilviscositydpl, dwaterviscositydpl;
2117  REAL oilviscosityr, waterviscosityr;
2118  REAL doilviscositydpr, dwaterviscositydpr;
2119 
2120  REAL bulklambdal, oillambdal, waterlambdal;
2121  REAL dbulklambdadpl, doillambdadpl, dwaterlambdadpl;
2122  REAL dbulklambdadsl, doillambdadsl, dwaterlambdadsl;
2123  REAL bulklambdar, oillambdar, waterlambdar;
2124  REAL dbulklambdadpr, doillambdadpr, dwaterlambdadpr;
2125  REAL dbulklambdadsr, doillambdadsr, dwaterlambdadsr;
2126 
2127  REAL bulkfoill, bulkfwaterl;
2128  REAL dbulkfoildpl, dbulkfwaterdpl;
2129  REAL dbulkfoildsl, dbulkfwaterdsl;
2130 
2131  REAL bulkfoilr, bulkfwaterr;
2132  REAL dbulkfoildpr, dbulkfwaterdpr;
2133  REAL dbulkfoildsr, dbulkfwaterdsr;
2134 
2135  REAL bulkfStarl;
2136  REAL dbulkfStardpl;
2137  REAL dbulkfStardsl;
2138 
2139  REAL bulkfStarr;
2140  REAL dbulkfStardpr;
2141  REAL dbulkfStardsr;
2142 
2143  REAL Pcl;
2144  REAL dPcdSl;
2145 
2146  REAL Pcr;
2147  REAL dPcdSr;
2148 
2149  // Functions computed at point x_{k} for each integration point
2150  int VecPos= 0;
2151  // REAL PressureRef = 1.0e6;
2152  this->Porosity(sol_pL[VecPos], rockporosityl, drockporositydpl);
2153  this->RhoOil(sol_pL[VecPos], oildensityl, doildensitydpl);
2154  this->RhoWater(sol_pL[VecPos], waterdensityl, dwaterdensitydpl);
2155  this->OilViscosity(sol_pL[VecPos], oilviscosityl, doilviscositydpl);
2156  this->WaterViscosity(sol_pL[VecPos], waterviscosityl, dwaterviscositydpl);
2157  this->OilLabmda(oillambdal, sol_pL[VecPos], sol_sL[VecPos], doillambdadpl, doillambdadsl);
2158  this->WaterLabmda(waterlambdal, sol_pL[VecPos], sol_sL[VecPos], dwaterlambdadpl, dwaterlambdadsl);
2159  this->Labmda(bulklambdal, sol_pL[VecPos], sol_sL[VecPos], dbulklambdadpl, dbulklambdadsl);
2160  this->fOil(bulkfoill, sol_pL[VecPos], sol_sL[VecPos], dbulkfoildpl, dbulkfoildsl);
2161  this->fWater(bulkfwaterl, sol_pL[VecPos], sol_sL[VecPos], dbulkfwaterdpl, dbulkfwaterdsl);
2162  this->CapillaryPressure(sol_sL[VecPos],Pcl,dPcdSl);
2163 
2164  this->Porosity(sol_pR[VecPos], rockporosityr, drockporositydpr);
2165  this->RhoOil(sol_pR[VecPos], oildensityr, doildensitydpr);
2166  this->RhoWater(sol_pR[VecPos], waterdensityr, dwaterdensitydpr);
2167  this->OilViscosity(sol_pR[VecPos], oilviscosityr, doilviscositydpr);
2168  this->WaterViscosity(sol_pR[VecPos], waterviscosityr, dwaterviscositydpr);
2169  this->OilLabmda(oillambdar, sol_pR[VecPos], sol_sR[VecPos], doillambdadpr, doillambdadsr);
2170  this->WaterLabmda(waterlambdar, sol_pR[VecPos], sol_sR[VecPos], dwaterlambdadpr, dwaterlambdadsr);
2171  this->Labmda(bulklambdar, sol_pR[VecPos], sol_sR[VecPos], dbulklambdadpr, dbulklambdadsr);
2172  this->fOil(bulkfoilr, sol_pR[VecPos], sol_sR[VecPos], dbulkfoildpr, dbulkfoildsr);
2173  this->fWater(bulkfwaterr, sol_pR[VecPos], sol_sR[VecPos], dbulkfwaterdpr, dbulkfwaterdsr);
2174  this->CapillaryPressure(sol_sR[VecPos],Pcr,dPcdSr);
2175 
2176  if (fYorN)
2177  {
2178 
2179  KabsoluteLeft=fKabsoluteMap[GeoIDLeft];
2180  KabsoluteRight=fKabsoluteMap[GeoIDRight];
2181 
2182  }
2183  else
2184  {
2185  this->K(KabsoluteLeft);
2186  this->K(KabsoluteRight);
2187  }
2188 
2189  Gfield = this->Gravity();
2190  REAL Gravitydotnl = Gfield(0,0)*(n1) + Gfield(1,0)*(n2);
2191  REAL Gravitydotnr = Gfield(0,0)*(n1) + Gfield(1,0)*(n2);
2192  this->fstar(bulkfStarl,sol_pL[VecPos], sol_sL[VecPos],1.0*Gravitydotnl,dbulkfStardpl,dbulkfStardsl);
2193  this->fstar(bulkfStarr,sol_pR[VecPos], sol_sR[VecPos],-1.0*Gravitydotnr,dbulkfStardpr,dbulkfStardsr);
2194 
2195  // Min of fstar at Gamma
2196  REAL bulkfstar = 0.0;
2197  REAL dbulkfstardp = 0.0;
2198  REAL dbulkfstards = 0.0;
2199  if(fabs(bulkfStarl) >= fabs(bulkfStarr))
2200  {
2201  bulkfstar = bulkfStarr;
2202  dbulkfstardp = dbulkfStardpr;
2203  dbulkfstards = dbulkfStardsr;
2204  }
2205  else
2206  {
2207  bulkfstar = bulkfStarl;
2208  dbulkfstardp = dbulkfStardpl;
2209  dbulkfstards = dbulkfStardsl;
2210  }
2211 
2212 
2213  for (int in=0; in < 3; in++) {
2214  for (int jn=0; jn < 3; jn++) {
2215  if (KabsoluteLeft(in,jn)+KabsoluteRight(in,jn)==0) {
2216  Kmean(in,jn)= 0.0;
2217  }
2218  else
2219  {
2220  Kmean(in,jn)= (2.0*KabsoluteLeft(in,jn)*KabsoluteRight(in,jn))/(KabsoluteLeft(in,jn)+KabsoluteRight(in,jn));
2221  }
2222 
2223  }
2224  }
2225 
2226  TPZFMatrix<STATE> KGL(3,3,0.0),KGR(3,3,0.0), KG(3,3,0.0);
2227  KGL(0,0) = KabsoluteLeft(0,0)*Gfield(0,0)+KabsoluteLeft(0,1)*Gfield(1,0);
2228  KGL(1,0) = KabsoluteLeft(1,0)*Gfield(0,0)+KabsoluteLeft(1,1)*Gfield(1,0);
2229  KGR(0,0) = KabsoluteRight(0,0)*Gfield(0,0)+KabsoluteRight(0,1)*Gfield(1,0);
2230  KGR(1,0) = KabsoluteRight(1,0)*Gfield(0,0)+KabsoluteRight(1,1)*Gfield(1,0);
2231 
2232  KG(0,0) = Kmean(0,0)*Gfield(0,0)+Kmean(0,1)*Gfield(1,0);
2233  KG(1,0) = Kmean(1,0)*Gfield(0,0)+Kmean(1,1)*Gfield(1,0);
2234 
2235 // REAL GravityFluxL = (KabsoluteLeft(0,0)*Gfield(0,0) + KabsoluteLeft(0,1)*Gfield(1,0))*(n1)+
2236 // (KabsoluteLeft(1,0)*Gfield(0,0) + KabsoluteLeft(1,1)*Gfield(1,0))*(n2);
2237 
2238 // REAL GravityFluxR = (KabsoluteRight(0,0)*Gfield(0,0) + KabsoluteRight(0,1)*Gfield(1,0))*(n1)+
2239 // (KabsoluteRight(1,0)*Gfield(0,0) + KabsoluteRight(1,1)*Gfield(1,0))*(n2);
2240 
2241 
2242  int URowsleft = phiUL.Rows();
2243  int URowsRight = phiUR.Rows();
2244 
2245  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
2246  int QRowsRight = dataright[1].fVecShapeIndex.NElements();
2247 
2248  int PRowsleft = phiPL.Rows();
2249  int PRowsRight = phiPR.Rows();
2250 
2251  int SRowsleft = phiSL.Rows();
2252  int SRowsRight = phiSR.Rows();
2253 
2254  int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
2255 // int QGRowsRight = dataright[4].fVecShapeIndex.NElements();
2256 
2257  int UStateVar = 2;
2258 
2259  int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
2260  int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
2261 
2262  int FirstUL = 0;
2263  int FirstQL = URowsleft * UStateVar + FirstUL;
2264  int FirstPL = QRowsleft + FirstQL;
2265  int FirstSL = PRowsleft + FirstPL;
2266  int FirstQGL = SRowsleft + FirstSL;
2267 
2268  int FirstUR = 0;
2269  int FirstQR = URowsRight * UStateVar + FirstUR;
2270  int FirstPR = QRowsRight + FirstQR;
2271  int FirstSR = PRowsRight + FirstPR;
2272 // int FirstQGR = SRowsRight + FirstSR;
2273 
2274 
2275  if(gState == ECurrentState)
2276  {
2277 
2278  // ////////////////////////// Jacobian Matrix ///////////////////////////////////
2279  // Contribution of contour integrals for Jacobian matrix
2280 
2281  // First Block (Equation One) constitutive law
2282  // Integrate[L dot(K v, n), Gamme_{e}] (Equation One) Left-Left part
2283 
2284  for (int iq=0; iq < QRowsleft; iq++)
2285  {
2286 
2287  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
2288  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
2289 
2290  for (int jp=0; jp < PRowsleft; jp++)
2291  {
2292 
2293  if (fnewWS)
2294  {
2295  REAL e1e1 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))*(n1);
2296  REAL e2e2 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex))*(n2);
2297  ek(iq + FirstQL, jp + FirstPL) += (-1.0) * weight * (e1e1 + e2e2 ) * phiPL(jp,0);
2298  }
2299  else
2300  {
2301  REAL e1e1 = (Kmean(0,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
2302  Kmean(0,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n1);
2303 
2304  REAL e2e2 = (Kmean(1,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
2305  Kmean(1,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n2);
2306  ek(iq + FirstQL, jp + FirstPL) += (-1.0) * weight * (e1e1 + e2e2 ) * phiPL(jp,0);
2307  }
2308 
2309  }
2310 
2311  }
2312 
2313  // First Block (Equation One) constitutive law
2314  // Integrate[L dot(K v, n), Gamme_{e}] (Equation One) Right-Right Part
2315  for (int iq=0; iq < QRowsRight; iq++)
2316  {
2317  int iRvectorindex = dataright[1].fVecShapeIndex[iq].first;
2318  int iRshapeindex = dataright[1].fVecShapeIndex[iq].second;
2319 
2320  for (int jp=0; jp < PRowsRight; jp++)
2321  {
2322 
2323  if (fnewWS)
2324  {
2325 
2326  REAL e1e1 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))*(n1);
2327  REAL e2e2 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex))*(n2);
2328 
2329  ek(iq + FirstQR + iRightInterfaceBlock,jp + FirstPR + jRightInterfaceBlock) += (1.0) * weight * (e1e1 + e2e2 ) * phiPR(jp,0) ;
2330  }
2331  else
2332  {
2333  REAL e1e1 = (Kmean(0,0)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))+
2334  Kmean(0,1)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex)))*(n1);
2335 
2336  REAL e2e2 = (Kmean(1,0)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))+
2337  Kmean(1,1)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex)))*(n2);
2338 
2339  ek(iq + FirstQR + iRightInterfaceBlock,jp + FirstPR + jRightInterfaceBlock) += (1.0) * weight * (e1e1 + e2e2 ) * phiPR(jp,0) ;
2340  }
2341 
2342 
2343 
2344  }
2345 
2346  }
2347 
2348  REAL dSwPcdSL = SaturationL * dPcdSl + Pcl;
2349  REAL dSwPcdSR = SaturationR * dPcdSr + Pcr;
2350 
2351  //This block was verified
2352  // First Block (Equation One) constitutive law Capillary Pressure
2353  // Integrate[Sw Pc dot( v, n), Gamme_{e}] (Equation One) Left-Left part
2354  for (int iq=0; iq < QRowsleft; iq++)
2355  {
2356  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
2357  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
2358  for (int jsat=0; jsat < SRowsleft; jsat++)
2359  {
2360  REAL e1e1 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))*(n1);
2361  REAL e2e2 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex))*(n2);
2362  ek(iq + FirstQL, jsat + FirstSL) += (-1.0) * (-1.0) * weight * dSwPcdSL * (e1e1 + e2e2 ) * phiSL(jsat,0);
2363  }
2364  }
2365 
2366  //This block was verified
2367  // First Block (Equation One) constitutive law Capillary Pressure
2368  // Integrate[Sw Pc dot(v, n), Gamme_{e}] (Equation One) Right-Right Part
2369  for (int iq=0; iq < QRowsRight; iq++)
2370  {
2371  int iRvectorindex = dataright[1].fVecShapeIndex[iq].first;
2372  int iRshapeindex = dataright[1].fVecShapeIndex[iq].second;
2373 
2374  for (int jsat=0; jsat < SRowsleft; jsat++)
2375  {
2376  REAL e1e1 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))*(n1);
2377  REAL e2e2 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex))*(n2);
2378  ek(iq + FirstQR + iRightInterfaceBlock,jsat + FirstSR + jRightInterfaceBlock) += (-1.0) * (1.0) * weight * dSwPcdSR * (e1e1 + e2e2 ) * phiSR(jsat,0) ;
2379  }
2380 
2381  }
2382 
2383  // Second Block (Equation Two) Bulk flux equation
2384  // Integrate[L dot(v, n), Gamma_{e}] (Equation Two) Left-Left Part
2385  for (int ip=0; ip < PRowsleft; ip++)
2386  {
2387 
2388  for (int jq=0; jq<QRowsleft; jq++)
2389  {
2390 
2391  int jvectorindex = dataleft[1].fVecShapeIndex[jq].first;
2392  int jshapeindex = dataleft[1].fVecShapeIndex[jq].second;
2393 
2394  REAL dotprod =
2395  (n1) * (phiQL(jshapeindex,0)*dataleft[1].fNormalVec(0,jvectorindex)) +
2396  (n2) * (phiQL(jshapeindex,0)*dataleft[1].fNormalVec(1,jvectorindex)) ;
2397 
2398  ek(ip + FirstPL,jq + FirstQL) += (-1.0) * (Gamma) * (TimeStep) * weight * dotprod * phiPL(ip,0);
2399 
2400  }
2401 
2402  }
2403 
2404 
2405  // Second Block (Equation Two) Bulk flux equation
2406  // Integrate[L dot(v, n), Gamma_{e}] (Equation Two) Right-Right Part
2407  for (int ip=0; ip < PRowsRight; ip++)
2408  {
2409 
2410  for (int jq=0; jq<QRowsRight; jq++)
2411  {
2412 
2413  int jvectorindex = dataright[1].fVecShapeIndex[jq].first;
2414  int jshapeindex = dataright[1].fVecShapeIndex[jq].second;
2415 
2416  REAL dotprod =
2417  (n1) * (phiQR(jshapeindex,0)*dataright[1].fNormalVec(0,jvectorindex)) +
2418  (n2) * (phiQR(jshapeindex,0)*dataright[1].fNormalVec(1,jvectorindex)) ;
2419 
2420  ek(ip + FirstPR + iRightInterfaceBlock,jq + FirstQR + jRightInterfaceBlock) += (1.0) * (Gamma) * (TimeStep) * weight * dotprod * phiPR(ip,0);
2421 
2422  }
2423 
2424  }
2425 
2426 
2427 
2428  // Upwind scheme
2429  // Third Vector Block (Equation three) Saturation equation
2430 
2431  REAL UpwindSaturation = 0.0;
2432 
2433  if (dotqnL > 0.0)
2434  {
2435  UpwindSaturation = bulkfwaterl;
2436 
2437  // Theta * TimeStep * Integrate[L L^{upwind} dot(v, n), Gamme_{e}] (Equation three) Left-Left Part
2438  for(int isat=0; isat<SRowsleft; isat++)
2439  {
2440  for(int jsat=0; jsat<SRowsleft; jsat++)
2441  {
2442  ek(isat + FirstSL ,jsat + FirstSL)
2443  += weight * (Theta) * (TimeStep) * phiSL(isat,0) * dbulkfwaterdsl * phiSL(jsat,0) * dotqnL;
2444  }
2445  }
2446 
2447  // Theta * TimeStep * Integrate[L L^{upwind} dot(v, n), Gamme_{e}] (Equation three) Right-Left Part
2448  for(int isat=0; isat<SRowsRight; isat++)
2449  {
2450 
2451  for(int jsat=0; jsat<SRowsleft; jsat++)
2452  {
2453  ek(isat + FirstSR + iRightInterfaceBlock,jsat + FirstSL)
2454  -= weight * (Theta) * (TimeStep) * phiSR(isat,0) * dbulkfwaterdsl * phiSL(jsat,0) * dotqnL;
2455  }
2456  }
2457 
2458  }
2459 
2460  else
2461 
2462  {
2463  UpwindSaturation = bulkfwaterr;
2464 
2465  // Theta * TimeStep * Integrate[L L^{upwind} dot(v, n), Gamme_{e}] (Equation three) Left-Right Part
2466  for(int isat=0; isat<SRowsleft; isat++)
2467  {
2468  for(int jsat=0; jsat<SRowsRight; jsat++)
2469  {
2470  ek(isat + FirstSL,jsat + FirstSR + jRightInterfaceBlock)
2471  += weight * (Theta) * (TimeStep) * phiSL(isat,0) * dbulkfwaterdsr * phiSR(jsat,0) * dotqnR;
2472  }
2473  }
2474  // Theta * TimeStep * Integrate[L L^{upwind} dot(v, n), Gamme_{e}] (Equation three) Right-Right Part
2475  for(int isat=0; isat<SRowsRight; isat++)
2476  {
2477 
2478  for(int jsat=0; jsat<SRowsRight; jsat++)
2479  {
2480  ek(isat + FirstSR + iRightInterfaceBlock,jsat + FirstSR + jRightInterfaceBlock)
2481  -= weight * (Theta) * (TimeStep) * phiSR(isat,0) * dbulkfwaterdsr * phiSR(jsat,0) * dotqnR;
2482  }
2483  }
2484 
2485  }
2486 
2487 
2488  // Third Vector Block (Equation three) Saturation equation
2489  // Theta * TimeStep * Integrate[L S dot(v, n), Gamme_{e}] (Equation three) Left-Left Part
2490  for (int isat=0; isat < SRowsleft; isat++) {
2491 
2492  for (int jq=0; jq < QRowsleft; jq++)
2493  {
2494  int jLvectorindex = dataleft[1].fVecShapeIndex[jq].first;
2495  int jLshapeindex = dataleft[1].fVecShapeIndex[jq].second;
2496 
2497  REAL dotprodL =
2498  (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(0,jLvectorindex)) * (n1) +
2499  (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(1,jLvectorindex)) * (n2) ;//+
2500  // (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(2,jLvectorindex)) * (n3) ; // dot(q,n) left
2501 
2502  ek(isat + FirstSL ,jq + FirstQL) += weight * (Theta) * (TimeStep) * phiSL(isat,0) * UpwindSaturation * dotprodL;
2503 
2504  }
2505  }
2506 
2507  // Theta * TimeStep * Integrate[L S dot(v, n), Gamme_{e}] (Equation three) Right-Left Part
2508  for (int isat=0; isat < SRowsRight; isat++)
2509  {
2510  for (int jq=0; jq < QRowsleft; jq++)
2511  {
2512  int jLvectorindex = dataleft[1].fVecShapeIndex[jq].first;
2513  int jLshapeindex = dataleft[1].fVecShapeIndex[jq].second;
2514 
2515  REAL dotprodL =
2516  (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(0,jLvectorindex)) * (n1) +
2517  (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(1,jLvectorindex)) * (n2) ;//+
2518  // (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(2,jLvectorindex)) * (n3) ; // dot(q,n) left
2519 
2520  ek(isat+ FirstSR + iRightInterfaceBlock,jq + FirstQL) -= weight * (Theta) * (TimeStep) * phiSR(isat,0) * UpwindSaturation * dotprodL;
2521 
2522  }
2523  }
2524 
2525  REAL QGgstarL = bulklambdal * (waterdensityl - oildensityl) * (KGL(0,0)*n1 + KGL(1,0)*n2);
2526  REAL QGgstarR = bulklambdar * (waterdensityr - oildensityr) * (KGR(0,0)*n1 + KGR(1,0)*n2);
2527 
2528  // It is needed to implement the complete derivative of P
2529  REAL dQGgstarLdP = ((bulklambdal * (dwaterdensitydpl - doildensitydpl)) +
2530  (dbulklambdadpl * (waterdensityl - oildensityl))) * (KGL(0,0)*n1 + KGL(1,0)*n2);
2531 
2532  REAL dQGgstarRdP = ((bulklambdar * (dwaterdensitydpr - doildensitydpr)) +
2533  (dbulklambdadpr * (waterdensityr - oildensityr))) * (KGR(0,0)*n1 + KGR(1,0)*n2);
2534 
2535  REAL dQGgstarLdS = (dbulklambdadsl) * (waterdensityl - oildensityl) * (KGL(0,0)*n1 + KGL(1,0)*n2);
2536  REAL dQGgstarRdS = (dbulklambdadsr) * (waterdensityr - oildensityr) * (KGR(0,0)*n1 + KGR(1,0)*n2);
2537 
2538  REAL QGstar = 0.0;
2539  REAL dQGstardP = 0.0;
2540  REAL dQGstardS = 0.0;
2541 
2542 // if(Gravitydotnl > 0.0 )
2543 // {
2544  if(fabs(1.0*bulkfStarl*QGgstarL) < fabs(1.0*bulkfStarr*QGgstarR))
2545  {
2546  QGstar = 1.0 * bulkfStarl * QGgstarL;
2547  dQGstardS = 1.0*(dbulkfStardsl * QGgstarL + bulkfStarl * dQGgstarLdS);
2548  dQGstardP = 1.0*(dbulkfStardpl * QGgstarL + bulkfStarl * dQGgstarLdP);
2549 
2550  // Gravitational segregation scheme
2551  // Four Block (Equation Four) gravitational flux constitutive law
2552  // Integrate[L dot(K v, n), Gamma_{e}] (Equation One) Left-Left part
2553  for (int iqg=0; iqg < QGRowsleft; iqg++)
2554  {
2555  int iLvectorindex = dataleft[4].fVecShapeIndex[iqg].first;
2556  int iLshapeindex = dataleft[4].fVecShapeIndex[iqg].second;
2557  REAL e1e1i = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(0,iLvectorindex))*(n1);
2558  REAL e2e2i = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(1,iLvectorindex))*(n2);
2559 
2560  for(int jsat=0; jsat<SRowsleft; jsat++)
2561  {
2562  // This degree of freedom is equal for Left and Right part
2563  ek(iqg + FirstQGL, jsat + FirstSL)
2564  -= weight * (dQGstardS * phiSL(jsat,0)) * (e1e1i + e2e2i);
2565  }
2566 
2567  for(int jp=0; jp<PRowsleft; jp++)
2568  {
2569  // This degree of freedom is equal for Left and Right part
2570  ek(iqg + FirstQGL, jp + FirstPL)
2571  -= weight * (dQGstardP * phiPL(jp,0)) * (e1e1i + e2e2i);
2572  }
2573 
2574  }
2575 
2576  }
2577  else
2578  {
2579  QGstar = 1.0 *bulkfStarr*QGgstarR;
2580  dQGstardS = 1.0*(dbulkfStardsr * QGgstarR + bulkfStarr * dQGgstarRdS);
2581  dQGstardP = 1.0*(dbulkfStardpr * QGgstarR + bulkfStarr * dQGgstarRdP);
2582 
2583  // Gravitational segregation scheme
2584  // Four Block (Equation Four) gravitational flux constitutive law
2585  // Integrate[L dot(K v, n), Gamma_{e}] (Equation One) Left-Left part
2586  for (int iqg=0; iqg < QGRowsleft; iqg++)
2587  {
2588  int iLvectorindex = dataleft[4].fVecShapeIndex[iqg].first;
2589  int iLshapeindex = dataleft[4].fVecShapeIndex[iqg].second;
2590  REAL e1e1i = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(0,iLvectorindex))*(n1);
2591  REAL e2e2i = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(1,iLvectorindex))*(n2);
2592 
2593  for(int jsat=0; jsat<SRowsRight; jsat++)
2594  {
2595  // This degree of freedom is equal for Left and Right part
2596  ek(iqg + FirstQGL, jRightInterfaceBlock + jsat + FirstSR)
2597  -= weight * (dQGstardS * phiSR(jsat,0)) * (e1e1i + e2e2i);
2598  }
2599 
2600  for(int jp=0; jp<PRowsRight; jp++)
2601  {
2602  // This degree of freedom is equal for Left and Right part
2603  ek(iqg + FirstQGL, jRightInterfaceBlock + jp + FirstPR )
2604  -= weight * (dQGstardP * phiPR(jp,0)) * (e1e1i + e2e2i);
2605  }
2606 
2607  }
2608 
2609  }
2610 
2611 
2612  // Gravitational segregation scheme
2613  // Four Block (Equation Four) gravitational flux constitutive law
2614  // Integrate[L dot(K v, n), Gamma_{e}] (Equation One) Left-Left part
2615  for (int iqg=0; iqg < QGRowsleft; iqg++)
2616  {
2617  int iLvectorindex = dataleft[4].fVecShapeIndex[iqg].first;
2618  int iLshapeindex = dataleft[4].fVecShapeIndex[iqg].second;
2619  REAL e1e1i = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(0,iLvectorindex))*(n1);
2620  REAL e2e2i = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(1,iLvectorindex))*(n2);
2621 
2622  for (int jqg=0; jqg < QGRowsleft; jqg++)
2623  {
2624  int jLvectorindex = dataleft[4].fVecShapeIndex[jqg].first;
2625  int jLshapeindex = dataleft[4].fVecShapeIndex[jqg].second;
2626  REAL e1e1j = (phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(0,jLvectorindex))*(n1);
2627  REAL e2e2j = (phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(1,jLvectorindex))*(n2);
2628 
2629  // This degree of freedom is equal for Left and Right part
2630  ek(iqg + FirstQGL, jqg + FirstQGL)
2631  += weight * (e1e1j + e2e2j) * (e1e1i + e2e2i);
2632  }
2633 
2634  }
2635 
2636 
2637  // Gravitational segregation scheme
2638  // (Theta) * deltat * Integrate[L*dot(qg,n), Gamma_{e}] (Equation three) Left-Left Part
2639  for(int isat=0; isat < SRowsleft; isat++)
2640  {
2641  for (int jqg=0; jqg < QGRowsleft; jqg++)
2642  {
2643  int jLvectorindex = dataleft[4].fVecShapeIndex[jqg].first;
2644  int jLshapeindex = dataleft[4].fVecShapeIndex[jqg].second;
2645  REAL e1e1j = (phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(0,jLvectorindex))*(n1);
2646  REAL e2e2j = (phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(1,jLvectorindex))*(n2);
2647 
2648  REAL ResidualPart = (Theta) * (TimeStep) * ( e1e1j + e2e2j );
2649  ek(isat + FirstSL,jqg + FirstQGL)
2650  += weight * phiSL(isat,0) * ResidualPart;
2651  }
2652  }
2653 
2654  // (Theta) * deltat * Integrate[L* dot(qg,n), Gamma_{e}] (Equation three) Right-Left Part
2655  for(int isat=0; isat < SRowsRight; isat++)
2656  {
2657  for (int jqg=0; jqg < QGRowsleft; jqg++)
2658  {
2659  int jLvectorindex = dataleft[4].fVecShapeIndex[jqg].first;
2660  int jLshapeindex = dataleft[4].fVecShapeIndex[jqg].second;
2661  REAL e1e1j = (phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(0,jLvectorindex))*(n1);
2662  REAL e2e2j = (phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(1,jLvectorindex))*(n2);
2663 
2664  REAL ResidualPart = (Theta) * (TimeStep) * ( e1e1j + e2e2j );
2665  ek(isat + FirstSR + iRightInterfaceBlock,jqg + FirstQGL)
2666  -= weight * phiSR(isat,0) * ResidualPart;
2667  }
2668  }
2669 
2670 
2671  }
2672 
2673 
2674  // ////////////////////////// Jacobian Matrix ///////////////////////////////////
2675  // End of contribution of countour integrals for Jacobian matrix
2676 
2677 
2678 
2679  // n time step
2680  // This values are constant in Newton iteration
2681  if(gState == ELastState)
2682  {
2683 
2684  }
2685 
2686  // if (fnewWS) {
2687  // Compute Residual contribution ef for contour integrals
2688  this->ContributeInterface(data, dataleft, dataright, weight, ef);
2689  // }
2690 
2691 
2692 }
2693 
2694 
2695 
2696 void TPZMultiphase::ContributeInterface(TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, TPZVec<TPZMaterialData> &dataright, REAL weight, TPZFMatrix<STATE> &ef)
2697 {
2698 
2699  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
2700  TPZFMatrix<REAL> &phiUR = dataright[0].phi;
2701 
2702  TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
2703  TPZFMatrix<REAL> &phiQR = dataright[1].phi;
2704 
2705  TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
2706  TPZFMatrix<REAL> &phiPR = dataright[2].phi;
2707 
2708  TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
2709  TPZFMatrix<REAL> &phiSR = dataright[3].phi;
2710 
2711  TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
2712 // TPZFMatrix<REAL> &phiQGR = dataright[4].phi;
2713 
2714  TPZManVector<REAL,3> &normal = data.normal;
2715 
2716 
2717  REAL n1 = normal[0];
2718  REAL n2 = normal[1];
2719 
2720  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
2721  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
2722  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
2723  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
2724  TPZManVector<STATE,3> sol_qgL =dataleft[4].sol[0];
2725 
2726  TPZManVector<STATE,3> sol_uR =dataright[0].sol[0];
2727  TPZManVector<STATE,3> sol_qR =dataright[1].sol[0];
2728  TPZManVector<STATE,3> sol_pR =dataright[2].sol[0];
2729  TPZManVector<STATE,3> sol_sR =dataright[3].sol[0];
2730  TPZManVector<STATE,3> sol_qgR =dataright[4].sol[0];
2731 
2732 
2733  // Getting Q solution for left and right side
2734  REAL qxL = sol_qL[0];
2735  REAL qyL = sol_qL[1];
2736  REAL qxR = sol_qR[0];
2737  REAL qyR = sol_qR[1];
2738  // REAL qz = sol_qR[2];
2739 
2740  // Getting Qg solution for left and right side
2741  REAL qgxL = sol_qgL[0];
2742  REAL qgyL = sol_qgL[1];
2743  REAL qgxR = sol_qgR[0];
2744  REAL qgyR = sol_qgR[1];
2745 
2746  // Getting P solution for left and right side
2747  REAL PseudoPressureL = sol_pL[0];
2748  REAL PseudoPressureR = sol_pR[0];
2749 
2750  // Getting S solution for left and right side
2751  REAL SaturationL = sol_sL[0];
2752  REAL SaturationR = sol_sR[0];
2753 
2754  REAL dotqnL = (qxL*n1) + (qyL*n2);
2755  REAL dotqnR = (qxR*n1) + (qyR*n2);
2756 
2757  REAL dotqgnL = (qgxL*n1) + (qgyL*n2);
2758  REAL dotqgnR = (qgxR*n1) + (qgyR*n2);
2759 
2760  // Getting another required data
2761  REAL TimeStep = this->fDeltaT;
2762  REAL Theta = this->fTheta;
2763  REAL Gamma = this->fGamma;
2764  this->fnewWS=true;
2765 
2766  // Getting Harmonic mean of permeabilities
2767  TPZFMatrix<REAL> KabsoluteLeft;
2768  TPZFMatrix<REAL> KabsoluteRight;
2769  TPZFMatrix<REAL> Gfield;
2770  TPZFMatrix<REAL> Kmean(3,3,0);
2771 
2772  int GeoIDLeft = dataleft[1].gelElId;
2773  int GeoIDRight = dataright[1].gelElId;
2774 
2775  REAL rockporosityl, oildensityl, waterdensityl;
2776  REAL drockporositydpl, doildensitydpl, dwaterdensitydpl;
2777  REAL rockporosityr, oildensityr, waterdensityr;
2778  REAL drockporositydpr, doildensitydpr, dwaterdensitydpr;
2779 
2780  REAL oilviscosityl, waterviscosityl;
2781  REAL doilviscositydpl, dwaterviscositydpl;
2782  REAL oilviscosityr, waterviscosityr;
2783  REAL doilviscositydpr, dwaterviscositydpr;
2784 
2785  REAL bulklambdal, oillambdal, waterlambdal;
2786  REAL dbulklambdadpl, doillambdadpl, dwaterlambdadpl;
2787  REAL dbulklambdadsl, doillambdadsl, dwaterlambdadsl;
2788  REAL bulklambdar, oillambdar, waterlambdar;
2789  REAL dbulklambdadpr, doillambdadpr, dwaterlambdadpr;
2790  REAL dbulklambdadsr, doillambdadsr, dwaterlambdadsr;
2791 
2792 
2793  REAL bulkfoill, bulkfwaterl;
2794  REAL dbulkfoildpl, dbulkfwaterdpl;
2795  REAL dbulkfoildsl, dbulkfwaterdsl;
2796 
2797  REAL bulkfoilr, bulkfwaterr;
2798  REAL dbulkfoildpr, dbulkfwaterdpr;
2799  REAL dbulkfoildsr, dbulkfwaterdsr;
2800 
2801  REAL bulkfStarl;
2802  REAL dbulkfStardpl;
2803  REAL dbulkfStardsl;
2804 
2805  REAL bulkfStarr;
2806  REAL dbulkfStardpr;
2807  REAL dbulkfStardsr;
2808 
2809  REAL Pcl;
2810  REAL dPcdSl;
2811 
2812  REAL Pcr;
2813  REAL dPcdSr;
2814 
2815  // Functions computed at point x_{k} for each integration point
2816  int VecPos= 0;
2817  // REAL PressureRef = 1.0e6;
2818  this->Porosity(sol_pL[VecPos], rockporosityl, drockporositydpl);
2819  this->RhoOil(sol_pL[VecPos], oildensityl, doildensitydpl);
2820  this->RhoWater(sol_pL[VecPos], waterdensityl, dwaterdensitydpl);
2821  this->OilViscosity(sol_pL[VecPos], oilviscosityl, doilviscositydpl);
2822  this->WaterViscosity(sol_pL[VecPos], waterviscosityl, dwaterviscositydpl);
2823  this->OilLabmda(oillambdal, sol_pL[VecPos], sol_sL[VecPos], doillambdadpl, doillambdadsl);
2824  this->WaterLabmda(waterlambdal, sol_pL[VecPos], sol_sL[VecPos], dwaterlambdadpl, dwaterlambdadsl);
2825  this->Labmda(bulklambdal, sol_pL[VecPos], sol_sL[VecPos], dbulklambdadpl, dbulklambdadsl);
2826  this->fOil(bulkfoill, sol_pL[VecPos], sol_sL[VecPos], dbulkfoildpl, dbulkfoildsl);
2827  this->fWater(bulkfwaterl, sol_pL[VecPos], sol_sL[VecPos], dbulkfwaterdpl, dbulkfwaterdsl);
2828  this->CapillaryPressure(sol_sL[VecPos],Pcl,dPcdSl);
2829 
2830  this->Porosity(sol_pR[VecPos], rockporosityr, drockporositydpr);
2831  this->RhoOil(sol_pR[VecPos], oildensityr, doildensitydpr);
2832  this->RhoWater(sol_pR[VecPos], waterdensityr, dwaterdensitydpr);
2833  this->OilViscosity(sol_pR[VecPos], oilviscosityr, doilviscositydpr);
2834  this->WaterViscosity(sol_pR[VecPos], waterviscosityr, dwaterviscositydpr);
2835  this->OilLabmda(oillambdar, sol_pR[VecPos], sol_sR[VecPos], doillambdadpr, doillambdadsr);
2836  this->WaterLabmda(waterlambdar, sol_pR[VecPos], sol_sR[VecPos], dwaterlambdadpr, dwaterlambdadsr);
2837  this->Labmda(bulklambdar, sol_pR[VecPos], sol_sR[VecPos], dbulklambdadpr, dbulklambdadsr);
2838  this->fOil(bulkfoilr, sol_pR[VecPos], sol_sR[VecPos], dbulkfoildpr, dbulkfoildsr);
2839  this->fWater(bulkfwaterr, sol_pR[VecPos], sol_sR[VecPos], dbulkfwaterdpr, dbulkfwaterdsr);
2840  this->CapillaryPressure(sol_sR[VecPos],Pcr,dPcdSr);
2841 
2842  if (fYorN)
2843  {
2844 
2845  KabsoluteLeft=fKabsoluteMap[GeoIDLeft];
2846 
2847  KabsoluteRight=fKabsoluteMap[GeoIDRight];
2848  }
2849  else
2850  {
2851  this->K(KabsoluteLeft);
2852  this->K(KabsoluteRight);
2853  }
2854 
2855  Gfield = this->Gravity();
2856  REAL Gravitydotnl = Gfield(0,0)*(n1) + Gfield(1,0)*(n2);
2857  REAL Gravitydotnr = Gfield(0,0)*(n1) + Gfield(1,0)*(n2);
2858  this->fstar(bulkfStarl,sol_pL[VecPos], sol_sL[VecPos],1.0*Gravitydotnl,dbulkfStardpl,dbulkfStardsl);
2859  this->fstar(bulkfStarr,sol_pR[VecPos], sol_sR[VecPos],-1.0*Gravitydotnr,dbulkfStardpr,dbulkfStardsr);
2860 
2861  // Min of fstar at Gamma
2862  REAL bulkfstar = 0.0;
2863  REAL dbulkfstardp = 0.0;
2864  REAL dbulkfstards = 0.0;
2865  if(fabs(bulkfStarl) >= fabs(bulkfStarr))
2866  {
2867  bulkfstar = bulkfStarr;
2868  dbulkfstardp = dbulkfStardpr;
2869  dbulkfstards = dbulkfStardsr;
2870  }
2871  else
2872  {
2873  bulkfstar = bulkfStarl;
2874  dbulkfstardp = dbulkfStardpl;
2875  dbulkfstards = dbulkfStardsl;
2876  }
2877 
2878 
2879  for (int in=0; in < 3; in++) {
2880  for (int jn=0; jn < 3; jn++) {
2881  if (KabsoluteLeft(in,jn)+KabsoluteRight(in,jn)<=0) {
2882  Kmean(in,jn)= 0.0;
2883  }
2884  else
2885  {
2886  Kmean(in,jn)= (2.0*KabsoluteLeft(in,jn)*KabsoluteRight(in,jn))/(KabsoluteLeft(in,jn)+KabsoluteRight(in,jn));
2887  }
2888 
2889  }
2890  }
2891 
2892  TPZFMatrix<STATE> KGL(3,3,0.0),KGR(3,3,0.0), KG(3,3,0.0);
2893  KGL(0,0) = KabsoluteLeft(0,0)*Gfield(0,0)+KabsoluteLeft(0,1)*Gfield(1,0);
2894  KGL(1,0) = KabsoluteLeft(1,0)*Gfield(0,0)+KabsoluteLeft(1,1)*Gfield(1,0);
2895  KGR(0,0) = KabsoluteRight(0,0)*Gfield(0,0)+KabsoluteRight(0,1)*Gfield(1,0);
2896  KGR(1,0) = KabsoluteRight(1,0)*Gfield(0,0)+KabsoluteRight(1,1)*Gfield(1,0);
2897 
2898  KG(0,0) = Kmean(0,0)*Gfield(0,0)+Kmean(0,1)*Gfield(1,0);
2899  KG(1,0) = Kmean(1,0)*Gfield(0,0)+Kmean(1,1)*Gfield(1,0);
2900 
2901 // REAL GravityFluxL = (KabsoluteLeft(0,0)*Gfield(0,0) + KabsoluteLeft(0,1)*Gfield(1,0))*(n1)+
2902 // (KabsoluteLeft(1,0)*Gfield(0,0) + KabsoluteLeft(1,1)*Gfield(1,0))*(n2);
2903 //
2904 // REAL GravityFluxR = (KabsoluteRight(0,0)*Gfield(0,0) + KabsoluteRight(0,1)*Gfield(1,0))*(n1)+
2905 // (KabsoluteRight(1,0)*Gfield(0,0) + KabsoluteRight(1,1)*Gfield(1,0))*(n2);
2906 
2907  int URowsleft = phiUL.Rows();
2908  int URowsRight = phiUR.Rows();
2909  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
2910  int QRowsRight = dataright[1].fVecShapeIndex.NElements();
2911  int PRowsleft = phiPL.Rows();
2912  int PRowsRight = phiPR.Rows();
2913  int SRowsleft = phiSL.Rows();
2914  int SRowsRight = phiSR.Rows();
2915 
2916  int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
2917 // int QGRowsRight = dataright[4].fVecShapeIndex.NElements();
2918 
2919  int UStateVar = 2;
2920 
2921  int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
2922 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
2923 
2924  int FirstUL = 0;
2925  int FirstQL = URowsleft * UStateVar + FirstUL;
2926  int FirstPL = QRowsleft + FirstQL;
2927  int FirstSL = PRowsleft + FirstPL;
2928  int FirstQGL = SRowsleft + FirstSL;
2929 
2930  int FirstUR = 0;
2931  int FirstQR = URowsRight * UStateVar + FirstUR;
2932  int FirstPR = QRowsRight + FirstQR;
2933  int FirstSR = PRowsRight + FirstPR;
2934 // int FirstQGR = SRowsRight + FirstSR;
2935 
2936 
2937  // ////////////////////////// Residual Vector ///////////////////////////////////
2938  // Contribution of contour integrals for Residual Vector
2939 
2940  // n+1 time step
2941  if(gState == ECurrentState)
2942  {
2943 
2944 
2945  //This block was verified
2946  // First Block (Equation One) constitutive law
2947  // Integrate[P dot(K v, n), Gamme_{e}] (Equation One) Left-Left part
2948  for (int iq=0; iq < QRowsleft; iq++)
2949  {
2950 
2951  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
2952  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
2953 
2954  if (fnewWS)
2955  {
2956 
2957  REAL e1e1 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))*(n1);
2958  REAL e2e2 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex))*(n2);
2959 
2960  ef(iq + FirstQL) += (-1.0) * weight * (e1e1 + e2e2 ) * PseudoPressureL;
2961  }
2962  else
2963  {
2964  REAL e1e1 = (Kmean(0,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
2965  Kmean(0,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n1);
2966 
2967  REAL e2e2 = (Kmean(1,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
2968  Kmean(1,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n2);
2969 
2970  ef(iq + FirstQL) += (-1.0) * weight * (e1e1 + e2e2 ) * PseudoPressureL;
2971  }
2972 
2973 
2974  }
2975 
2976 
2977 
2978  //This block was verified
2979  // First Block (Equation One) constitutive law
2980  // Integrate[P dot(K v, n), Gamme_{e}] (Equation One) Right-Right Part
2981  for (int iq=0; iq < QRowsRight; iq++)
2982  {
2983 
2984  int iRvectorindex = dataright[1].fVecShapeIndex[iq].first;
2985  int iRshapeindex = dataright[1].fVecShapeIndex[iq].second;
2986 
2987  if (fnewWS)
2988  {
2989 
2990  REAL e1e1 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))*(n1);
2991  REAL e2e2 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex))*(n2);
2992 
2993 
2994  ef(iq + iRightInterfaceBlock + FirstQR) += (1.0) * weight * (e1e1 + e2e2 ) * PseudoPressureR;
2995  }
2996  else
2997  {
2998  REAL e1e1 = (Kmean(0,0)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))+
2999  Kmean(0,1)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex)))*(n1);
3000 
3001  REAL e2e2 = (Kmean(1,0)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))+
3002  Kmean(1,1)*(phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex)))*(n2);
3003 
3004  ef(iq + iRightInterfaceBlock + FirstQR) += (1.0) * weight * (e1e1 + e2e2 ) * PseudoPressureR;
3005  }
3006 
3007  }
3008 
3009 
3010  REAL SwPcL = SaturationL * Pcl;
3011  REAL SwPcR = SaturationR * Pcr;
3012 
3013  //This block was verified
3014  // First Block (Equation One) constitutive law Capillary Pressure
3015  // Integrate[Sw Pc dot( v, n), Gamme_{e}] (Equation One) Left-Left part
3016  for (int iq=0; iq < QRowsleft; iq++)
3017  {
3018  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
3019  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
3020  REAL e1e1 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))*(n1);
3021  REAL e2e2 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex))*(n2);
3022  ef(iq + FirstQL) += (-1.0) * (-1.0) * weight * (e1e1 + e2e2 ) * SwPcL;
3023 
3024  }
3025 
3026  //This block was verified
3027  // First Block (Equation One) constitutive law Capillary Pressure
3028  // Integrate[Sw Pc dot(v, n), Gamme_{e}] (Equation One) Right-Right Part
3029  for (int iq=0; iq < QRowsRight; iq++)
3030  {
3031  int iRvectorindex = dataright[1].fVecShapeIndex[iq].first;
3032  int iRshapeindex = dataright[1].fVecShapeIndex[iq].second;
3033  REAL e1e1 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(0,iRvectorindex))*(n1);
3034  REAL e2e2 = (phiQR(iRshapeindex,0)*dataright[1].fNormalVec(1,iRvectorindex))*(n2);
3035 
3036  ef(iq + iRightInterfaceBlock + FirstQR) += (-1.0) * (1.0) * weight * (e1e1 + e2e2 ) * SwPcR;
3037  }
3038 
3039  // Second Block (Equation Two) Bulk flux equation
3040  // Integrate[L dot(q, n), Gamme_{e}] (Equation Two) Left-Left Part
3041  for (int ip=0; ip < PRowsleft; ip++)
3042  {
3043  ef(ip + FirstPL) += (-1.0) * (Gamma) * (TimeStep) * weight * dotqnL * phiPL(ip,0);
3044  }
3045 
3046  // Second Block (Equation Two) Bulk flux equation
3047  // Integrate[L dot(q, n), Gamme_{e}] (Equation Two) Right-Right Part
3048  for (int ip=0; ip < PRowsRight; ip++)
3049  {
3050  ef(ip + FirstPR + iRightInterfaceBlock) += (1.0) * (Gamma) * (TimeStep) * weight * dotqnR * phiPR(ip,0);
3051  }
3052 
3053 
3054  REAL UpwindSaturation = 0.0;
3055 
3056  if (dotqnL > 0.0)
3057  {
3058  UpwindSaturation = bulkfwaterl;
3059 
3060  }
3061  else
3062  {
3063  UpwindSaturation = bulkfwaterr;
3064 
3065  }
3066 
3067  // Third Vector Block (Equation three)
3068  // (Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Left-Left Part
3069  for(int isat=0; isat < SRowsleft; isat++)
3070  {
3071  REAL ResidualPart = (Theta) * (TimeStep) * ( UpwindSaturation * dotqnL );
3072  ef(isat + FirstSL) += weight * phiSL(isat,0) * ResidualPart;
3073  }
3074 
3075  // (Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Right-Left Part
3076  for(int isat=0; isat < SRowsRight; isat++)
3077  {
3078  REAL ResidualPart = (Theta) * (TimeStep) * ( UpwindSaturation * dotqnR );
3079  ef(isat + iRightInterfaceBlock + FirstSR) -= weight * phiSR(isat,0) * ResidualPart;
3080  }
3081 
3082  REAL QGgstarL = bulklambdal * (waterdensityl - oildensityl) * (KGL(0,0)*n1 + KGL(1,0)*n2);
3083  REAL QGgstarR = bulklambdar * (waterdensityr - oildensityr) * (KGR(0,0)*n1 + KGR(1,0)*n2);
3084  REAL QGstar = 0.0;
3085 
3086 // if(Gravitydotnl > 0.0 )
3087 // {
3088  if(fabs(1.0*bulkfStarl*QGgstarL) < fabs(1.0*bulkfStarr*QGgstarR))
3089  {
3090  QGstar = 1.0*bulkfStarl*QGgstarL;
3091  }
3092  else
3093  {
3094  QGstar = 1.0*bulkfStarr*QGgstarR;
3095  }
3096 
3097 // #ifdef LOG4CXX
3098 // if(logdata->isDebugEnabled())
3099 // {
3100 // std::stringstream sout;
3101 // sout << "SLeft =" << sol_sL[VecPos] << " bulkfStarl =" << bulkfStarl << " Gravitydotnl =" << Gravitydotnl << std::endl;
3102 // sout << "SRight =" << sol_sR[VecPos] << " bulkfStarr =" << bulkfStarr << " Gravitydotnr =" << -1.0*Gravitydotnr << std::endl;
3103 // sout << "QGstarL =" << 1.0*bulkfStarl*QGgstarL << " Gravitydotnl =" << Gravitydotnl << std::endl;
3104 // sout << "QGstarR =" << 1.0*bulkfStarr*QGgstarR << " Gravitydotnl =" << -1.0*Gravitydotnr << std::endl;
3105 // sout << "QGstar =" << QGstar << std::endl;
3106 // LOGPZ_DEBUG(logdata,sout.str());
3107 // }
3108 // #endif
3109 
3110  // Gravitational segregation scheme
3111  // Four Block (Equation Four) gravitational flux constitutive law
3112  // Integrate[L dot(K v, n), Gamma{e}] (Equation One) Left-Left part
3113  for (int iqg=0; iqg < QGRowsleft; iqg++)
3114  {
3115  int iLvectorindex = dataleft[4].fVecShapeIndex[iqg].first;
3116  int iLshapeindex = dataleft[4].fVecShapeIndex[iqg].second;
3117  REAL e1e1 = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(0,iLvectorindex))*(n1);
3118  REAL e2e2 = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(1,iLvectorindex))*(n2);
3119  // This degree of freedom is the same for Left and Right part
3120  ef(iqg + FirstQGL) += weight * (dotqgnL-QGstar) * (e1e1 + e2e2);
3121 
3122  }
3123 
3124  // Gravitational segregation scheme
3125  // (Theta) * deltat * Integrate[L*dot(qg,n), Gamma_{e}] (Equation three) Left-Left Part
3126  for(int isat=0; isat < SRowsleft; isat++)
3127  {
3128  REAL ResidualPart = (Theta) * (TimeStep) * ( dotqgnL );
3129  ef(isat + FirstSL) += weight * phiSL(isat,0) * ResidualPart;
3130  }
3131 
3132  // (Theta) * deltat * Integrate[L* dot(qg,n), Gamma_{e}] (Equation three) Right-Left Part
3133  for(int isat=0; isat < SRowsRight; isat++)
3134  {
3135  REAL ResidualPart = (Theta) * (TimeStep) * ( dotqgnL );
3136  ef(isat + iRightInterfaceBlock + FirstSR) -= weight * phiSR(isat,0) * ResidualPart;
3137  }
3138 
3139 // // Four Block (Equation Four) gravitational flux constitutive law
3140 // // Integrate[L dot(K v, n), Gamme_{e}] (Equation One) Left-Left part
3141 // for (int iqg=0; iqg < QGRowsleft; iqg++)
3142 // {
3143 //
3144 // int iLvectorindex = dataleft[3].fVecShapeIndex[iqg].first;
3145 // int iLshapeindex = dataleft[3].fVecShapeIndex[iqg].second;
3146 //
3147 // if (fnewWS)
3148 // {
3149 // REAL e1e1 = (phiQGL(iLshapeindex,0)*dataleft[3].fNormalVec(0,iLvectorindex))*(n1);
3150 // REAL e2e2 = (phiQGL(iLshapeindex,0)*dataleft[3].fNormalVec(1,iLvectorindex))*(n2);
3151 // ef(iqg+QRowsleft+PRowsleft+SRowsleft)
3152 // += (-1.0) * weight * (e1e1 + e2e2 ) * (waterdensityl - oildensityl);
3153 // }
3154 // else
3155 // {
3156 // REAL e1e1 = (Kmean(0,0)*(phiQGL(iLshapeindex,0)*dataleft[3].fNormalVec(0,iLvectorindex))+
3157 // Kmean(0,1)*(phiQGL(iLshapeindex,0)*dataleft[3].fNormalVec(1,iLvectorindex)))*(n1);
3158 //
3159 // REAL e2e2 = (Kmean(1,0)*(phiQGL(iLshapeindex,0)*dataleft[3].fNormalVec(0,iLvectorindex))+
3160 // Kmean(1,1)*(phiQGL(iLshapeindex,0)*dataleft[3].fNormalVec(1,iLvectorindex)))*(n2);
3161 // ef(iqg+QRowsleft+PRowsleft+SRowsleft)
3162 // += (-1.0) * weight * (e1e1 + e2e2 ) * (waterdensityl - oildensityl);
3163 // }
3164 //
3165 // }
3166 //
3167 // // Four Block (Equation Four) gravitational flux constitutive law
3168 // // Integrate[L dot(K v, n), Gamme_{e}] (Equation One) Right-Right Part
3169 // for (int iqg=0; iqg < QGRowsRight; iqg++)
3170 // {
3171 // int iRvectorindex = dataright[3].fVecShapeIndex[iqg].first;
3172 // int iRshapeindex = dataright[3].fVecShapeIndex[iqg].second;
3173 //
3174 // if (fnewWS)
3175 // {
3176 //
3177 // REAL e1e1 = (phiQGR(iRshapeindex,0)*dataright[3].fNormalVec(0,iRvectorindex))*(n1);
3178 // REAL e2e2 = (phiQGR(iRshapeindex,0)*dataright[3].fNormalVec(1,iRvectorindex))*(n2);
3179 //
3180 // ef(iqg + iRightInterfaceBlock + QGRowsRight + PRowsRight + SRowsRight)
3181 // += (1.0) * weight * (e1e1 + e2e2 ) * (waterdensityr - oildensityr);
3182 // }
3183 // else
3184 // {
3185 // REAL e1e1 = (Kmean(0,0)*(phiQGR(iRshapeindex,0)*dataright[3].fNormalVec(0,iRvectorindex))+
3186 // Kmean(0,1)*(phiQGR(iRshapeindex,0)*dataright[3].fNormalVec(1,iRvectorindex)))*(n1);
3187 //
3188 // REAL e2e2 = (Kmean(1,0)*(phiQGR(iRshapeindex,0)*dataright[3].fNormalVec(0,iRvectorindex))+
3189 // Kmean(1,1)*(phiQGR(iRshapeindex,0)*dataright[3].fNormalVec(1,iRvectorindex)))*(n2);
3190 //
3191 // ef(iqg + iRightInterfaceBlock + QGRowsRight + PRowsRight + SRowsRight)
3192 // += (1.0) * weight * (e1e1 + e2e2 ) * (waterdensityr - oildensityr);
3193 // }
3194 //
3195 // }
3196 //
3197 
3198 
3199  }
3200 
3201 
3202 
3203 
3204  // n time step
3205  if(gState == ELastState)
3206  {
3207 
3208  // Second Block (Equation Two) Bulk flux equation
3209  // Integrate[L dot(q, n), Gamme_{e}] (Equation Two) Left-Left Part
3210  for (int ip=0; ip < PRowsleft; ip++)
3211  {
3212  ef(ip + FirstPL) -= (1-Gamma) * (TimeStep) * weight * dotqnL * phiPL(ip,0);
3213  }
3214 
3215  // Second Block (Equation Two) Bulk flux equation
3216  // Integrate[L dot(q, n), Gamme_{e}] (Equation Two) Right-Right Part
3217  for (int ip=0; ip < PRowsRight; ip++)
3218  {
3219  ef(ip + iRightInterfaceBlock + FirstPR) += (1-Gamma) * (TimeStep) * weight * dotqnR * phiPR(ip,0);
3220  }
3221 
3222  REAL dotqnL = (qxL*n1) + (qyL*n2);
3223  // REAL dotqnR = (qxR*n1) + (qyR*n2);
3224  REAL UpwindSaturation = 0.0;
3225 
3226  if (dotqnL > 0.0)
3227  {
3228  UpwindSaturation = bulkfwaterl;
3229 
3230  }
3231  else
3232  {
3233  UpwindSaturation = bulkfwaterr;
3234 
3235  }
3236 
3237 
3238  // Third Vector Block (Equation three)
3239  // (1-Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Left-Left Part
3240  for(int isat=0; isat<SRowsleft; isat++)
3241  {
3242  REAL ResidualPart = (1-Theta) * (TimeStep) * ( UpwindSaturation * dotqnL );
3243  ef(isat + FirstSL) += weight * phiSL(isat,0) * ResidualPart;
3244  }
3245 
3246  // (1-Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Right-Left Part
3247  for(int isat=0; isat<SRowsRight; isat++)
3248  {
3249  REAL ResidualPart = (1-Theta) * (TimeStep) * ( UpwindSaturation * dotqnL );
3250  ef(isat + iRightInterfaceBlock + FirstSR) -= weight * phiSR(isat,0) * ResidualPart;
3251  }
3252 
3253  // Gravitational segregation scheme
3254  // (Theta) * deltat * Integrate[L*dot(qg,n), Gamma_{e}] (Equation three) Left-Left Part
3255  for(int isat=0; isat<SRowsleft; isat++)
3256  {
3257  REAL ResidualPart = (1-Theta) * (TimeStep) * ( dotqgnL );
3258  ef(isat + FirstSL) += weight * phiSL(isat,0) * ResidualPart;
3259  }
3260 
3261  // (Theta) * deltat * Integrate[L* dot(qg,n), Gamma_{e}] (Equation three) Right-Left Part
3262  for(int isat=0; isat<SRowsRight; isat++)
3263  {
3264  REAL ResidualPart = (1-Theta) * (TimeStep) * ( dotqgnR);
3265  ef(isat + iRightInterfaceBlock + FirstSR) -= weight * phiSR(isat,0) * ResidualPart;
3266  }
3267 
3268  }
3269 
3270  // ////////////////////////// Residual Vector ///////////////////////////////////
3271  // End of contribution of contour integrals for Residual Vector
3272 
3273 }
3274 
3275 
3276 
3277 void TPZMultiphase::ContributeBC(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ef, TPZBndCond &bc)
3278 {
3279  DebugStop();
3280 }
3281 
3282 void TPZMultiphase::ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, TPZBndCond &bc)
3283 {
3284  DebugStop();
3285 }
3286 
3287 
3288 void TPZMultiphase::ContributeBC(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc){
3289 
3290 
3291 #ifdef PZDEBUG
3292  int nref = datavec.size();
3293  if (nref != 5 ) {
3294  std::cout << " Error.!! datavec size not equal to 4 \n";
3295  DebugStop();
3296  }
3297 #endif
3298 
3299 
3300 
3301 }
3302 
3303 void TPZMultiphase::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix<STATE> &ef,TPZBndCond &bc)
3304 {
3305  DebugStop();
3306 }
3307 
3308 void TPZMultiphase::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
3309 {
3310  DebugStop();
3311 }
3312 
3313 void TPZMultiphase::ContributeBCInterface(TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
3314 {
3315 
3316  int nref = dataleft.size();
3317  if (nref != 5) {
3318  std::cout << " Error:: datavec size must to be equal to 4 \n" << std::endl;
3319  DebugStop();
3320  }
3321 
3322  if (bc.Val2().Rows() != 6){
3323  std::cout << " Error:: This material need boundary conditions for qx, qy, p (pore pressure) and s (Saturation) .\n" << std::endl;
3324  std::cout << " give me one matrix with this form Val2(3,1).\n" << std::endl;
3325  DebugStop();
3326  }
3327 
3328  if (bc.Val1().Rows() != 6){
3329  std::cout << " Error:: This material need boundary conditions for ux, uy, qx, qy, p (pore pressure) and s (Saturation) .\n" << std::endl;
3330  DebugStop();
3331  }
3332 
3333  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
3334  TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
3335  TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
3336  TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
3337  TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
3338 
3339  int URowsleft = phiUL.Rows();
3340  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
3341  // int QRowsleft1 = phiQL.Rows();
3342  int PRowsleft = phiPL.Rows();
3343  int SRowsleft = phiSL.Rows();
3344  int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
3345 
3346  int UStateVar = 2;
3347 
3348 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
3349 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
3350 
3351  int FirstUL = 0;
3352  int FirstQL = URowsleft * UStateVar + FirstUL;
3353  int FirstPL = QRowsleft + FirstQL;
3354  int FirstSL = PRowsleft + FirstPL;
3355  int FirstQGL = SRowsleft + FirstSL;
3356 
3357  TPZManVector<REAL,3> &normal = data.normal;
3358  REAL n1 = normal[0];
3359  REAL n2 = normal[1];
3360 
3361  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
3362  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
3363  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
3364  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
3365 
3366  // Getting Q solution for left and right side
3367 // REAL uxL = sol_uL[0];
3368 // REAL uyL = sol_uL[1];
3369  REAL qxL = sol_qL[0];
3370  REAL qyL = sol_qL[1];
3371  REAL dotqnL = (qxL*n1) + (qyL*n2);
3372 
3373  // Getting P solution for left side
3374  REAL PseudoPressureL = sol_pL[0];
3375 
3376  // Getting S solution for left side
3377  // REAL SaturationL = sol_sL[0];
3378 
3379  // Getting another required data
3380 
3381  REAL TimeStep = this->fDeltaT;
3382  REAL Theta = this->fTheta;
3383  REAL Gamma = this->fGamma;
3384  this->fnewWS=true;
3385 
3386  // Getting Harmonic mean of permeabilities
3387  int GeoIDLeft = dataleft[1].gelElId;
3388  TPZFMatrix<REAL> Kabsolute;
3389  TPZFMatrix<REAL> Kinverse;
3390  TPZFMatrix<REAL> Gfield;
3391 
3392  if (fYorN)
3393  {
3394  Kabsolute=this->fKabsoluteMap[GeoIDLeft];
3395  }
3396  else
3397  {
3398  this->K(Kabsolute);
3399  }
3400 
3401  Kinverse = this->Kinv(Kabsolute);
3402  Gfield = this->Gravity();
3403 
3404  REAL GravityFluxL = (Kabsolute(0,0)*Gfield(0,0) + Kabsolute(0,1)*Gfield(1,0))*(n1)+
3405  (Kabsolute(1,0)*Gfield(0,0) + Kabsolute(1,1)*Gfield(1,0))*(n2);
3406 
3407  REAL rockporosityl, oildensityl, waterdensityl;
3408  REAL drockporositydpl, doildensitydpl, dwaterdensitydpl;
3409 
3410  REAL oilviscosityl, waterviscosityl;
3411  REAL doilviscositydpl, dwaterviscositydpl;
3412 
3413  REAL bulklambdal, oillambdal, waterlambdal;
3414  REAL dbulklambdadpl, doillambdadpl, dwaterlambdadpl;
3415  REAL dbulklambdadsl, doillambdadsl, dwaterlambdadsl;
3416 
3417  REAL bulkfoill, bulkfwaterl;
3418  REAL dbulkfoildpl, dbulkfwaterdpl;
3419  REAL dbulkfoildsl, dbulkfwaterdsl;
3420 
3421 
3422  // Functions computed at point x_{k} for each integration point
3423  int VecPos= 0;
3424 
3425  this->Porosity(sol_pL[VecPos], rockporosityl, drockporositydpl);
3426  this->RhoOil(sol_pL[VecPos], oildensityl, doildensitydpl);
3427  this->RhoWater(sol_pL[VecPos], waterdensityl, dwaterdensitydpl);
3428  this->OilViscosity(sol_pL[VecPos], oilviscosityl, doilviscositydpl);
3429  this->WaterViscosity(sol_pL[VecPos], waterviscosityl, dwaterviscositydpl);
3430  this->OilLabmda(oillambdal, sol_pL[VecPos], sol_sL[VecPos], doillambdadpl, doillambdadsl);
3431  this->WaterLabmda(waterlambdal, sol_pL[VecPos], sol_sL[VecPos], dwaterlambdadpl, dwaterlambdadsl);
3432  this->Labmda(bulklambdal, sol_pL[VecPos], sol_sL[VecPos], dbulklambdadpl, dbulklambdadsl);
3433  this->fOil(bulkfoill, sol_pL[VecPos], sol_sL[VecPos], dbulkfoildpl, dbulkfoildsl);
3434  this->fWater(bulkfwaterl, sol_pL[VecPos], sol_sL[VecPos], dbulkfwaterdpl, dbulkfwaterdsl);
3435 
3437  // Regular Controur integrals
3438 
3439  // n+1 time step
3440  if(gState == ECurrentState)
3441  {
3442 
3443  //First Block (Equation One) constitutive law
3444  //Integrate[L dot(K v, n), Gamma_{e}] (Equation One) Left-Left part
3445  for (int iq=0; iq < QRowsleft; iq++)
3446  {
3447 
3448  for (int jp=0; jp < PRowsleft; jp++)
3449  {
3450  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
3451  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
3452 
3453  if (fnewWS)
3454  {
3455  REAL e1e1 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))*(n1);
3456  REAL e2e2 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex))*(n2);
3457 
3458  ek(iq + FirstQL, jp + FirstPL) += (-1.0) * weight * (e1e1 + e2e2 ) * phiPL(jp,0);
3459 
3460  }
3461  else
3462  {
3463  REAL e1e1 = (Kabsolute(0,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
3464  Kabsolute(0,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n1);
3465 
3466  REAL e2e2 = (Kabsolute(1,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
3467  Kabsolute(1,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n2);
3468 
3469  ek(iq + FirstQL, jp + FirstPL) += (-1.0) * weight * (e1e1 + e2e2 ) * phiPL(jp,0);
3470 
3471  }
3472  }
3473  }
3474 
3475  // Second Block (Equation Two) Bulk flux equation
3476  // Integrate[L dot(v, n), Gamme_{e}] (Equation Two) Left-Left Part
3477  for (int ip=0; ip < PRowsleft; ip++)
3478  {
3479 
3480  for (int jq=0; jq< QRowsleft; jq++)
3481  {
3482 
3483  int jvectorindex = dataleft[1].fVecShapeIndex[jq].first;
3484  int jshapeindex = dataleft[1].fVecShapeIndex[jq].second;
3485 
3486  REAL dotprod =
3487  (n1) * (phiQL(jshapeindex,0)*dataleft[1].fNormalVec(0,jvectorindex)) +
3488  (n2) * (phiQL(jshapeindex,0)*dataleft[1].fNormalVec(1,jvectorindex)) ;
3489 
3490  ek(ip + FirstPL,jq + FirstQL) += (-1.0) * (Gamma) * (TimeStep) * weight * dotprod * phiPL(ip,0);
3491 
3492  }
3493 
3494  }
3495 
3496 
3497  // This block was verified
3498  // First Block (Equation One) constitutive law
3499  // Integrate[P dot(K v, n), Gamme_{e}] (Equation One) Left-Left part
3500  for (int iq=0; iq < QRowsleft; iq++)
3501  {
3502 
3503  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
3504  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
3505 
3506  if (fnewWS)
3507  {
3508 
3509  REAL e1e1 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))*(n1);
3510  REAL e2e2 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex))*(n2);
3511 
3512  ef(iq + FirstQL) += (-1.0) * weight * (e1e1 + e2e2) * PseudoPressureL;
3513 
3514  }
3515  else
3516  {
3517  REAL e1e1 = (Kabsolute(0,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
3518  Kabsolute(0,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n1);
3519 
3520  REAL e2e2 = (Kabsolute(1,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
3521  Kabsolute(1,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n2);
3522 
3523 
3524  ef(iq + FirstQL) += (-1.0) * weight * (e1e1 + e2e2) * PseudoPressureL;
3525 
3526  }
3527 
3528  }
3529 
3530 
3531 
3532  for (int ip=0; ip < PRowsleft; ip++)
3533  {
3534  ef(ip+FirstPL) += (-1.0) * (Gamma) * (TimeStep) * weight * dotqnL * phiPL(ip,0);
3535  }
3536 
3537  // Null gravitational flux on boundaries
3538  for(int iqg=0; iqg < QGRowsleft; iqg++)
3539  {
3540  int iLvectorindex = dataleft[4].fVecShapeIndex[iqg].first;
3541  int iLshapeindex = dataleft[4].fVecShapeIndex[iqg].second;
3542 
3543  REAL vni = (phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(0,iLvectorindex)*n1)+(phiQGL(iLshapeindex,0)*dataleft[4].fNormalVec(1,iLvectorindex)*n2);
3544 
3545  for (int jqg=0; jqg < QGRowsleft; jqg++)
3546  {
3547  int jLvectorindex = dataleft[4].fVecShapeIndex[jqg].first;
3548  int jLshapeindex = dataleft[4].fVecShapeIndex[jqg].second;
3549 
3550  REAL vnj = (phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(0,jLvectorindex)*n1)+(phiQGL(jLshapeindex,0)*dataleft[4].fNormalVec(1,jLvectorindex)*n2);
3551 
3552  ek(iqg + FirstQGL,jqg + FirstQGL) += weight * (gBigNumber * ( vnj ) * vni );
3553  }
3554  }
3555 
3556 
3557  }
3558 
3559 
3560  // n time step
3561  if(gState == ELastState)
3562  {
3563 
3564  // Second Block (Equation Two) Bulk flux equation
3565  // Integrate[L dot(q, n), Gamme_{e}] (Equation Two) Left-Left Part
3566  for (int ip=0; ip < PRowsleft; ip++)
3567  {
3568  ef(ip + FirstPL) -= (1-Gamma) * (TimeStep) * weight * dotqnL * phiPL(ip,0);
3569  }
3570 
3571  REAL UpwindFSaturation = 0.0;
3572  REAL signofG = 1.0;
3573 
3574  if (dotqnL > 0.0)
3575  {
3576  UpwindFSaturation = bulkfwaterl;
3577 
3578  }
3579  else
3580  {
3581  UpwindFSaturation = 0.0;
3582 
3583  }
3584 
3585  // Third Vector Block (Equation three)
3586  // (Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Left-Left Part
3587  for(int isat=0; isat<SRowsleft; isat++)
3588  {
3589  REAL ResidualPart = (1-Theta) * (TimeStep) * ( UpwindFSaturation * bulkfoill * bulklambdal * (waterdensityl - oildensityl) * GravityFluxL );
3590  ef(isat + FirstSL) += signofG * (-1.0) * weight * phiSL(isat,0) * ResidualPart;
3591  }
3592 
3593  }
3594 
3595 
3596  // Regular Controur integrals
3598 
3599  STATE v2[6];
3600  v2[0] = bc.Val2()(0,0); // ux
3601  v2[1] = bc.Val2()(1,0); // uy
3602  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
3603  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
3604  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
3605  v2[5] = bc.Val2()(5,0); // Saturation
3606 // REAL qN = (v2[2]*n1 + v2[3]*n2); // Normal Flux
3607 
3608 
3609  switch (bc.Type()) {
3610  case 1 : // inflow BC bc: Ux, Uy, Qn, Sin
3611  {
3612  if(gState == ECurrentState)
3613  {
3614  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3615  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3616  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3617  ApplySin(data,dataleft,weight,ek,ef,bc);
3618  }
3619 
3620  if(gState == ELastState)
3621  {
3622  ApplySin(data,dataleft,weight,ef,bc);
3623  }
3624  }
3625  break;
3626 
3627  case 2 : // inflow BC bc: Tx, Ty, Qn, Sin
3628  {
3629  if(gState == ECurrentState)
3630  {
3631  ApplySigmaN(data,dataleft,weight,ek,ef,bc);
3632  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3633  ApplySin(data,dataleft,weight,ek,ef,bc);
3634  }
3635 
3636  if(gState == ELastState)
3637  {
3638  ApplySin(data,dataleft,weight,ef,bc);
3639  }
3640 
3641  }
3642  break;
3643 
3644  case 3 : // inflow BC bc: Ux, Qn, Sin
3645  {
3646  if(gState == ECurrentState)
3647  {
3648  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3649  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3650  ApplySin(data,dataleft,weight,ek,ef,bc);
3651  }
3652 
3653  if(gState == ELastState)
3654  {
3655  ApplySin(data,dataleft,weight,ef,bc);
3656  }
3657  }
3658  break;
3659 
3660  case 4 : // inflow BC bc: Uy, Qn, Sin
3661  {
3662  if(gState == ECurrentState)
3663  {
3664  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3665  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3666  ApplySin(data,dataleft,weight,ek,ef,bc);
3667  }
3668 
3669  if(gState == ELastState)
3670  {
3671  ApplySin(data,dataleft,weight,ef,bc);
3672  }
3673  }
3674  break;
3675 
3676  case 5 : // outflow BC bc: Ux, Uy, Qn, Sout
3677  {
3678  if(gState == ECurrentState)
3679  {
3680  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3681  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3682  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3683  ApplySout(data,dataleft,weight,ek,ef,bc);
3684  }
3685 
3686  if(gState == ELastState)
3687  {
3688  ApplySout(data,dataleft,weight,ef,bc);
3689  }
3690  }
3691  break;
3692 
3693  case 6 : // outflow BC bc: Tx, Ty, Qn, Sout
3694  {
3695  if(gState == ECurrentState)
3696  {
3697  ApplySigmaN(data,dataleft,weight,ek,ef,bc);
3698  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3699  ApplySout(data,dataleft,weight,ek,ef,bc);
3700  }
3701 
3702  if(gState == ELastState)
3703  {
3704  ApplySout(data,dataleft,weight,ef,bc);
3705  }
3706 
3707  }
3708  break;
3709 
3710  case 7 : // outflow BC bc: Ux, Qn, Sout
3711  {
3712  if(gState == ECurrentState)
3713  {
3714  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3715  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3716  ApplySout(data,dataleft,weight,ek,ef,bc);
3717  }
3718 
3719  if(gState == ELastState)
3720  {
3721  ApplySout(data,dataleft,weight,ef,bc);
3722  }
3723  }
3724  break;
3725 
3726  case 8 : // outflow BC bc: Uy, Qn, Sout
3727  {
3728  if(gState == ECurrentState)
3729  {
3730  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3731  ApplyQnD(data,dataleft,weight,ek,ef,bc);
3732  ApplySout(data,dataleft,weight,ek,ef,bc);
3733  }
3734 
3735  if(gState == ELastState)
3736  {
3737  ApplySout(data,dataleft,weight,ef,bc);
3738  }
3739  }
3740  break;
3741 
3742  case 9 : // inflow BC bc: Ux, Uy, PN, Sin
3743  {
3744  if(gState == ECurrentState)
3745  {
3746  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3747  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3748  ApplyPN(data,dataleft,weight,ek,ef,bc);
3749  ApplySin(data,dataleft,weight,ek,ef,bc);
3750  }
3751 
3752  if(gState == ELastState)
3753  {
3754  ApplySin(data,dataleft,weight,ef,bc);
3755  }
3756  }
3757  break;
3758 
3759  case 10 : // inflow BC bc: Tx, Ty, PN, Sin
3760  {
3761  if(gState == ECurrentState)
3762  {
3763  ApplySigmaN(data,dataleft,weight,ek,ef,bc);
3764  ApplyPN(data,dataleft,weight,ek,ef,bc);
3765  ApplySin(data,dataleft,weight,ek,ef,bc);
3766  }
3767 
3768  if(gState == ELastState)
3769  {
3770  ApplySin(data,dataleft,weight,ef,bc);
3771  }
3772 
3773  }
3774  break;
3775 
3776  case 11 : // inflow BC bc: Ux, PN, Sin
3777  {
3778  if(gState == ECurrentState)
3779  {
3780  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3781  ApplyPN(data,dataleft,weight,ek,ef,bc);
3782  ApplySin(data,dataleft,weight,ek,ef,bc);
3783  }
3784 
3785  if(gState == ELastState)
3786  {
3787  ApplySin(data,dataleft,weight,ef,bc);
3788  }
3789  }
3790  break;
3791 
3792  case 12 : // inflow BC bc: Uy, PN, Sin
3793  {
3794  if(gState == ECurrentState)
3795  {
3796  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3797  ApplyPN(data,dataleft,weight,ek,ef,bc);
3798  ApplySin(data,dataleft,weight,ek,ef,bc);
3799  }
3800 
3801  if(gState == ELastState)
3802  {
3803  ApplySin(data,dataleft,weight,ef,bc);
3804  }
3805  }
3806  break;
3807 
3808  case 13 : // outflow BC bc: Ux, Uy, PN, Sout
3809  {
3810  if(gState == ECurrentState)
3811  {
3812  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3813  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3814  ApplyPN(data,dataleft,weight,ek,ef,bc);
3815  ApplySout(data,dataleft,weight,ek,ef,bc);
3816  }
3817 
3818  if(gState == ELastState)
3819  {
3820  ApplySout(data,dataleft,weight,ef,bc);
3821  }
3822  }
3823  break;
3824 
3825  case 14 : // outflow BC bc: Tx, Ty, PN, Sout
3826  {
3827  if(gState == ECurrentState)
3828  {
3829  ApplySigmaN(data,dataleft,weight,ek,ef,bc);
3830  ApplyPN(data,dataleft,weight,ek,ef,bc);
3831  ApplySout(data,dataleft,weight,ek,ef,bc);
3832  }
3833 
3834  if(gState == ELastState)
3835  {
3836  ApplySout(data,dataleft,weight,ef,bc);
3837  }
3838 
3839  }
3840  break;
3841 
3842  case 15 : // outflow BC bc: Ux, PN, Sout
3843  {
3844  if(gState == ECurrentState)
3845  {
3846  ApplyUxD(data,dataleft,weight,ek,ef,bc);
3847  ApplyPN(data,dataleft,weight,ek,ef,bc);
3848  ApplySout(data,dataleft,weight,ek,ef,bc);
3849  }
3850 
3851  if(gState == ELastState)
3852  {
3853  ApplySout(data,dataleft,weight,ef,bc);
3854  }
3855  }
3856  break;
3857 
3858  case 16 : // outflow BC bc: Uy, PN, Sout
3859  {
3860  if(gState == ECurrentState)
3861  {
3862  ApplyUyD(data,dataleft,weight,ek,ef,bc);
3863  ApplyPN(data,dataleft,weight,ek,ef,bc);
3864  ApplySout(data,dataleft,weight,ek,ef,bc);
3865  }
3866 
3867  if(gState == ELastState)
3868  {
3869  ApplySout(data,dataleft,weight,ef,bc);
3870  }
3871  }
3872  break;
3873 
3874  default: std::cout << "This BC doesn't exist." << std::endl;
3875  {
3876 
3877  DebugStop();
3878  }
3879  break;
3880  }
3881 
3882 }
3883 
3884  void TPZMultiphase::ApplyUxD (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
3885  {
3886 
3887  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
3888 // TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
3889 // TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
3890 // TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
3891 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
3892 
3893  int URowsleft = phiUL.Rows();
3894 // int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
3895  // int QRowsleft1 = phiQL.Rows();
3896 // int PRowsleft = phiPL.Rows();
3897 // int SRowsleft = phiSL.Rows();
3898 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
3899 
3900 // int UStateVar = 2;
3901 
3902 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
3903 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
3904 
3905 // int FirstUL = 0;
3906 // int FirstQL = URowsleft * UStateVar + FirstUL;
3907 // int FirstPL = QRowsleft + FirstQL;
3908 // int FirstSL = PRowsleft + FirstPL;
3909 // int FirstQGL = SRowsleft + FirstSL;
3910 
3911 
3912  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
3913 
3914  // Getting Q solution for left and right side
3915  REAL uxL = sol_uL[0];
3916 // REAL uyL = sol_uL[1];
3917 
3918  STATE v2[6];
3919  v2[0] = bc.Val2()(0,0); // ux
3920  v2[1] = bc.Val2()(1,0); // uy
3921  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
3922  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
3923  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
3924  v2[5] = bc.Val2()(5,0); // Saturation
3925 
3926  // Dirichlet condition for each state variable
3927  // Elasticity Equation
3928  for(int iu = 0 ; iu < URowsleft; iu++)
3929  {
3930  // Contribution for load Vector
3931  ef(2*iu) += this->fxi * gBigNumber * weight * (uxL - v2[0]) * phiUL(iu,0); // X displacement Value
3932 // ef(2*iu+1) += gBigNumber * weight * (uyL - v2[1]) * phiUL(iu,0); // y displacement Value
3933 
3934  for(int ju = 0 ; ju < URowsleft; ju++)
3935  {
3936  // Contribution for Stiffness Matrix
3937  ek(2*iu,2*ju) += this->fxi * gBigNumber * weight * phiUL(iu,0) * phiUL(ju,0); // X displacement
3938 // ek(2*iu+1,2*ju+1) += gBigNumber * weight * phiUL(iu,0) * phiUL(ju,0); // Y displacement
3939  }
3940  }
3941 
3942 
3943 
3944  }
3945 
3946  void TPZMultiphase::ApplyUyD (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
3947  {
3948 
3949  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
3950 // TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
3951 // TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
3952 // TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
3953 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
3954 
3955  int URowsleft = phiUL.Rows();
3956 // int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
3957  // int QRowsleft1 = phiQL.Rows();
3958 // int PRowsleft = phiPL.Rows();
3959 // int SRowsleft = phiSL.Rows();
3960 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
3961 
3962 // int UStateVar = 2;
3963 
3964 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
3965 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
3966 
3967 // int FirstUL = 0;
3968 // int FirstQL = URowsleft * UStateVar + FirstUL;
3969 // int FirstPL = QRowsleft + FirstQL;
3970 // int FirstSL = PRowsleft + FirstPL;
3971 // int FirstQGL = SRowsleft + FirstSL;
3972 
3973 
3974  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
3975 
3976  // Getting Q solution for left and right side
3977 // REAL uxL = sol_uL[0];
3978  REAL uyL = sol_uL[1];
3979 
3980  STATE v2[6];
3981  v2[0] = bc.Val2()(0,0); // ux
3982  v2[1] = bc.Val2()(1,0); // uy
3983  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
3984  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
3985  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
3986  v2[5] = bc.Val2()(5,0); // Saturation
3987 
3988  // Dirichlet condition for each state variable
3989  // Elasticity Equation
3990  for(int iu = 0 ; iu < URowsleft; iu++)
3991  {
3992  // Contribution for load Vector
3993 // ef(2*iu) += gBigNumber * weight * (uxL - v2[0]) * phiUL(iu,0); // X displacement Value
3994  ef(2*iu+1) += this->fxi * gBigNumber * weight * (uyL - v2[1]) * phiUL(iu,0); // y displacement Value
3995 
3996  for(int ju = 0 ; ju < URowsleft; ju++)
3997  {
3998  // Contribution for Stiffness Matrix
3999 // ek(2*iu,2*ju) += gBigNumber * weight * phiUL(iu,0) * phiUL(ju,0); // X displacement
4000  ek(2*iu+1,2*ju+1) += this->fxi * gBigNumber * weight * phiUL(iu,0) * phiUL(ju,0); // Y displacement
4001  }
4002  }
4003 
4004 
4005  }
4006 
4007  void TPZMultiphase::ApplySigmaN (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
4008  {
4009 
4010  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
4011 // TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
4012 // TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
4013 // TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
4014 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
4015 
4016  int URowsleft = phiUL.Rows();
4017 // int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
4018 // // int QRowsleft1 = phiQL.Rows();
4019 // int PRowsleft = phiPL.Rows();
4020 // int SRowsleft = phiSL.Rows();
4021 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
4022 //
4023 // int UStateVar = 2;
4024 
4025 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4026 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4027 
4028 // int FirstUL = 0;
4029 // int FirstQL = URowsleft * UStateVar + FirstUL;
4030 // int FirstPL = QRowsleft + FirstQL;
4031 // int FirstSL = PRowsleft + FirstPL;
4032 // int FirstQGL = SRowsleft + FirstSL;
4033 
4034 
4035  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
4036 
4037  // Getting Q solution for left and right side
4038 // REAL uxL = sol_uL[0];
4039 // REAL uyL = sol_uL[1];
4040 
4041  STATE v2[6];
4042  v2[0] = bc.Val2()(0,0); // ux
4043  v2[1] = bc.Val2()(1,0); // uy
4044  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
4045  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
4046  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
4047  v2[5] = bc.Val2()(5,0); // Saturation
4048 
4049  // Dirichlet condition for each state variable
4050  // Elasticity Equation
4051  for(int iu = 0 ; iu < URowsleft; iu++)
4052  {
4053  // Contribution for load Vector
4054  ef(2*iu) += weight * ( v2[0]) * phiUL(iu,0); // X Traction Value
4055  ef(2*iu+1) += weight * ( v2[1]) * phiUL(iu,0); // y Traction Value
4056 
4057  }
4058 
4059  }
4060 
4061  void TPZMultiphase::ApplyQnD (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
4062  {
4063 
4064  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
4065  TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
4066 // TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
4067 // TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
4068 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
4069 
4070  int URowsleft = phiUL.Rows();
4071  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
4072 // int PRowsleft = phiPL.Rows();
4073 // int SRowsleft = phiSL.Rows();
4074 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
4075 
4076  int UStateVar = 2;
4077 
4078 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4079 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4080 
4081  int FirstUL = 0;
4082  int FirstQL = URowsleft * UStateVar + FirstUL;
4083 // int FirstPL = QRowsleft + FirstQL;
4084 // int FirstSL = PRowsleft + FirstPL;
4085 // int FirstQGL = SRowsleft + FirstSL;
4086 
4087  TPZManVector<REAL,3> &normal = data.normal;
4088  REAL n1 = normal[0];
4089  REAL n2 = normal[1];
4090 
4091  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
4092 
4093  // Getting Q solution for left and right side
4094  REAL qxL = sol_qL[0];
4095  REAL qyL = sol_qL[1];
4096  REAL dotqnL = (qxL*n1) + (qyL*n2);
4097 
4098 
4099  STATE v2[6];
4100  v2[0] = bc.Val2()(0,0); // ux
4101  v2[1] = bc.Val2()(1,0); // uy
4102  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
4103  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
4104  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
4105  v2[5] = bc.Val2()(5,0); // Saturation
4106  REAL qN = (v2[2]*n1 + v2[3]*n2); // Normal Flux
4107 
4108 
4109  // Phil's Hint
4110  for(int iq=0; iq < QRowsleft; iq++)
4111  {
4112  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
4113  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
4114 
4115  REAL vni = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex)*n1)+(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)*n2);
4116  // The coefficient 0.0001 is required for balance the residual contribution
4117  ef(iq + FirstQL) += weight * ( this->fxi * (gBigNumber * ( dotqnL - qN ) * vni ) );
4118  // ef(iq) += weight * ( (gBigNumber * ( dotqnL - qN ) * ( dotqnL - qN ) * vni ) );
4119  // ef(iq) += weight * ( (gBigNumber * ( dotqnL - qN ) * ( dotqnL - qN ) * ( dotqnL - qN ) * ( dotqnL - qN ) * vni ) );
4120 
4121  for (int jq=0; jq < QRowsleft; jq++)
4122  {
4123  int jLvectorindex = dataleft[1].fVecShapeIndex[jq].first;
4124  int jLshapeindex = dataleft[1].fVecShapeIndex[jq].second;
4125 
4126  REAL vnj = (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(0,jLvectorindex)*n1)+(phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(1,jLvectorindex)*n2);
4127  ek(iq + FirstQL,jq + FirstQL) += weight * ( this->fxi * (gBigNumber * ( vnj ) * vni ) );
4128  // ek(iq,jq) += weight * ( 2.0 * (gBigNumber * ( dotqnL - qN ) * ( vnj ) * vni ) );
4129  // ek(iq,jq) += weight * ( 4.0 * (gBigNumber * ( dotqnL - qN ) * ( dotqnL - qN ) * ( dotqnL - qN ) * ( vnj ) * vni ) );
4130  }
4131  }
4132  }
4133 
4134 
4135  void TPZMultiphase::ApplyPN (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
4136  {
4137 
4138  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
4139  TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
4140 // TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
4141 // TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
4142 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
4143 
4144  int URowsleft = phiUL.Rows();
4145  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
4146 // int PRowsleft = phiPL.Rows();
4147 // int SRowsleft = phiSL.Rows();
4148 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
4149 
4150  int UStateVar = 2;
4151 
4152 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4153 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4154 
4155  int FirstUL = 0;
4156  int FirstQL = URowsleft * UStateVar + FirstUL;
4157 // int FirstPL = QRowsleft + FirstQL;
4158 // int FirstSL = PRowsleft + FirstPL;
4159 // int FirstQGL = SRowsleft + FirstSL;
4160 
4161  TPZManVector<REAL,3> &normal = data.normal;
4162  REAL n1 = normal[0];
4163  REAL n2 = normal[1];
4164 
4165  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
4166  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
4167  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
4168  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
4169 
4170  // Getting Q solution for left and right side
4171 // REAL uxL = sol_uL[0];
4172 // REAL uyL = sol_uL[1];
4173 // REAL qxL = sol_qL[0];
4174 // REAL qyL = sol_qL[1];
4175 // REAL dotqnL = (qxL*n1) + (qyL*n2);
4176 
4177  // Getting P solution for left side
4178 // REAL PseudoPressureL = sol_pL[0];
4179 
4180  // Getting another required data
4181 
4182 // REAL TimeStep = this->fDeltaT;
4183 // REAL Theta = this->fTheta;
4184 // REAL Gamma = this->fGamma;
4185  this->fnewWS=true;
4186 
4187  // Getting Harmonic mean of permeabilities
4188  int GeoIDLeft = dataleft[1].gelElId;
4189  TPZFMatrix<REAL> Kabsolute;
4190  TPZFMatrix<REAL> Kinverse;
4191  TPZFMatrix<REAL> Gfield;
4192 
4193  if (fYorN)
4194  {
4195  Kabsolute=this->fKabsoluteMap[GeoIDLeft];
4196  }
4197  else
4198  {
4199  this->K(Kabsolute);
4200  }
4201 
4202  Kinverse = this->Kinv(Kabsolute);
4203  Gfield = this->Gravity();
4204 
4205  STATE v2[6];
4206  v2[0] = bc.Val2()(0,0); // ux
4207  v2[1] = bc.Val2()(1,0); // uy
4208  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
4209  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
4210  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
4211  v2[5] = bc.Val2()(5,0); // Saturation
4212 // REAL qN = (v2[2]*n1 + v2[3]*n2); // Normal Flux
4213 
4214  if(fBCForcingFunction)
4215  {
4216  TPZManVector<STATE> PValue(1);
4217  fBCForcingFunction->Execute(dataleft[2].x,PValue);
4218  v2[4] = PValue[0];
4219  }
4220 
4221  // This block was verified
4222  // First Block (Equation One) constitutive law
4223  // Integrate[P dot(K v, n), Gamme_{e}] (Equation One) Left-Left part
4224  for (int iq=0; iq < QRowsleft; iq++)
4225  {
4226 
4227  int iLvectorindex = dataleft[1].fVecShapeIndex[iq].first;
4228  int iLshapeindex = dataleft[1].fVecShapeIndex[iq].second;
4229 
4230  if (fnewWS)
4231  {
4232 
4233  REAL e1e1 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))*(n1);
4234  REAL e2e2 = (phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex))*(n2);
4235 
4236  ef(iq + FirstQL) += (1.0) * weight * (e1e1 + e2e2 ) * (v2[4]);
4237 
4238  }
4239  else
4240  {
4241  REAL e1e1 = (Kabsolute(0,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
4242  Kabsolute(0,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n1);
4243 
4244  REAL e2e2 = (Kabsolute(1,0)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(0,iLvectorindex))+
4245  Kabsolute(1,1)*(phiQL(iLshapeindex,0)*dataleft[1].fNormalVec(1,iLvectorindex)))*(n2);
4246 
4247  ef(iq + FirstQL) += (1.0) * weight * (e1e1 + e2e2 ) * (v2[4]);
4248 
4249  }
4250  }
4251 
4252  }
4253 
4254  void TPZMultiphase::ApplySin (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
4255  {
4256 
4257  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
4258 // TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
4259  TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
4260  TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
4261 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
4262 
4263  int URowsleft = phiUL.Rows();
4264  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
4265  int PRowsleft = phiPL.Rows();
4266  int SRowsleft = phiSL.Rows();
4267 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
4268 
4269  int UStateVar = 2;
4270 
4271 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4272 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4273 
4274  int FirstUL = 0;
4275  int FirstQL = URowsleft * UStateVar + FirstUL;
4276  int FirstPL = QRowsleft + FirstQL;
4277  int FirstSL = PRowsleft + FirstPL;
4278 // int FirstQGL = SRowsleft + FirstSL;
4279 
4280  TPZManVector<REAL,3> &normal = data.normal;
4281  REAL n1 = normal[0];
4282  REAL n2 = normal[1];
4283 
4284  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
4285  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
4286  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
4287  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
4288 
4289  // Getting Q solution for left and right side
4290 // REAL uxL = sol_uL[0];
4291 // REAL uyL = sol_uL[1];
4292 // REAL qxL = sol_qL[0];
4293 // REAL qyL = sol_qL[1];
4294 // REAL dotqnL = (qxL*n1) + (qyL*n2);
4295 
4296  // Getting P solution for left side
4297 // REAL PseudoPressureL = sol_pL[0];
4298 
4299  // Getting another required data
4300 
4301  REAL TimeStep = this->fDeltaT;
4302  REAL Theta = this->fTheta;
4303 // REAL Gamma = this->fGamma;
4304  this->fnewWS=true;
4305 
4306  // Getting Harmonic mean of permeabilities
4307  int GeoIDLeft = dataleft[1].gelElId;
4308  TPZFMatrix<REAL> Kabsolute;
4309  TPZFMatrix<REAL> Kinverse;
4310  TPZFMatrix<REAL> Gfield;
4311 
4312  if (fYorN)
4313  {
4314  Kabsolute=this->fKabsoluteMap[GeoIDLeft];
4315  }
4316  else
4317  {
4318  this->K(Kabsolute);
4319  }
4320 
4321  Kinverse = this->Kinv(Kabsolute);
4322  Gfield = this->Gravity();
4323 
4324 // REAL GravityFluxL = (Kabsolute(0,0)*Gfield(0,0) + Kabsolute(0,1)*Gfield(1,0))*(n1)+
4325 // (Kabsolute(1,0)*Gfield(0,0) + Kabsolute(1,1)*Gfield(1,0))*(n2);
4326 
4327  REAL rockporosityl, oildensityl, waterdensityl;
4328  REAL drockporositydpl, doildensitydpl, dwaterdensitydpl;
4329 
4330  REAL oilviscosityl, waterviscosityl;
4331  REAL doilviscositydpl, dwaterviscositydpl;
4332 
4333  REAL bulklambdal, oillambdal, waterlambdal;
4334  REAL dbulklambdadpl, doillambdadpl, dwaterlambdadpl;
4335  REAL dbulklambdadsl, doillambdadsl, dwaterlambdadsl;
4336 
4337  REAL bulkfoill, bulkfwaterl;
4338  REAL dbulkfoildpl, dbulkfwaterdpl;
4339  REAL dbulkfoildsl, dbulkfwaterdsl;
4340 
4341 
4342  // Functions computed at point x_{k} for each integration point
4343  int VecPos= 0;
4344 
4345  this->Porosity(sol_pL[VecPos], rockporosityl, drockporositydpl);
4346  this->RhoOil(sol_pL[VecPos], oildensityl, doildensitydpl);
4347  this->RhoWater(sol_pL[VecPos], waterdensityl, dwaterdensitydpl);
4348  this->OilViscosity(sol_pL[VecPos], oilviscosityl, doilviscositydpl);
4349  this->WaterViscosity(sol_pL[VecPos], waterviscosityl, dwaterviscositydpl);
4350  this->OilLabmda(oillambdal, sol_pL[VecPos], sol_sL[VecPos], doillambdadpl, doillambdadsl);
4351  this->WaterLabmda(waterlambdal, sol_pL[VecPos], sol_sL[VecPos], dwaterlambdadpl, dwaterlambdadsl);
4352  this->Labmda(bulklambdal, sol_pL[VecPos], sol_sL[VecPos], dbulklambdadpl, dbulklambdadsl);
4353  this->fOil(bulkfoill, sol_pL[VecPos], sol_sL[VecPos], dbulkfoildpl, dbulkfoildsl);
4354  this->fWater(bulkfwaterl, sol_pL[VecPos], sol_sL[VecPos], dbulkfwaterdpl, dbulkfwaterdsl);
4355 
4356  STATE v2[6];
4357  v2[0] = bc.Val2()(0,0); // ux
4358  v2[1] = bc.Val2()(1,0); // uy
4359  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
4360  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
4361  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
4362  v2[5] = bc.Val2()(5,0); // Saturation
4363  REAL qN = (v2[2]*n1 + v2[3]*n2); // Normal Flux
4364 
4365  // Upwind scheme
4366  // Third Vector Block (Equation three) Saturation equation
4367  REAL UpwindSaturation = 0.0;
4368 
4369 // if (dotqnL > 0.0)
4370 // {
4371 // this->fWater(bulkfwaterl, sol_pL[VecPos], v2[5], dbulkfwaterdpl, dbulkfwaterdsl);
4372 // UpwindSaturation = bulkfwaterl;
4373 //
4374 // // Theta * TimeStep * Integrate[L L^{upwind} dot(v, n), Gamme_{e}] (Equation three) Bc-Left Part
4375 // for(int isat=0; isat<SRowsleft; isat++)
4376 // {
4377 // for(int jsat=0; jsat<SRowsleft; jsat++)
4378 // {
4379 // ek(isat+QRowsleft+PRowsleft,jsat+QRowsleft+PRowsleft) -= (-1.0) * weight *
4380 // (Theta) * (TimeStep) * phiSL(isat,0)* dbulkfwaterdsl * phiSL(jsat,0) * dotqnL;
4381 // }
4382 // }
4383 //
4384 //
4385 // }
4386 //
4387 // else
4388 // {
4389  this->fWater(bulkfwaterl, sol_pL[VecPos], v2[5], dbulkfwaterdpl, dbulkfwaterdsl);
4390  UpwindSaturation = bulkfwaterl;
4391 
4392 // }
4393 
4394 
4395 // // Theta * TimeStep * Integrate[L S dot(v, n), Gamme_{e}] (Equation three) Right-Left Part
4396 // for (int isat=0; isat < SRowsleft; isat++) {
4397 //
4398 // for (int jq=0; jq<QRowsleft; jq++)
4399 // {
4400 // int jLvectorindex = dataleft[1].fVecShapeIndex[jq].first;
4401 // int jLshapeindex = dataleft[1].fVecShapeIndex[jq].second;
4402 //
4403 // REAL dotprodL =
4404 // (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(0,jLvectorindex)) * (n1) +
4405 // (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(1,jLvectorindex)) * (n2) ;//+
4406 //
4407 // ek(isat+QRowsleft+PRowsleft,jq) -= (-1.0) * weight * (Theta) * (TimeStep) * phiSL(isat,0) * (UpwindSaturation) * dotprodL;
4408 //
4409 // }
4410 //
4411 // }
4412 
4413  // (Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Right-Left Part
4414  for(int isat=0; isat < SRowsleft; isat++)
4415  {
4416  REAL ResidualPart = (Theta) * (TimeStep) * ( (UpwindSaturation) * qN );
4417  ef(isat + FirstSL) -= (-1.0) * weight * phiSL(isat,0) * ResidualPart;
4418  }
4419 
4420  }
4421 
4422  void TPZMultiphase::ApplySout (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
4423  {
4424 
4425  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
4426  TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
4427  TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
4428  TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
4429 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
4430 
4431  int URowsleft = phiUL.Rows();
4432  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
4433  int PRowsleft = phiPL.Rows();
4434  int SRowsleft = phiSL.Rows();
4435 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
4436 
4437  int UStateVar = 2;
4438 
4439 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4440 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4441 
4442  int FirstUL = 0;
4443  int FirstQL = URowsleft * UStateVar + FirstUL;
4444  int FirstPL = QRowsleft + FirstQL;
4445  int FirstSL = PRowsleft + FirstPL;
4446 // int FirstQGL = SRowsleft + FirstSL;
4447 
4448  TPZManVector<REAL,3> &normal = data.normal;
4449  REAL n1 = normal[0];
4450  REAL n2 = normal[1];
4451 
4452  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
4453  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
4454  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
4455  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
4456 
4457  // Getting Q solution for left and right side
4458 // REAL uxL = sol_uL[0];
4459 // REAL uyL = sol_uL[1];
4460  REAL qxL = sol_qL[0];
4461  REAL qyL = sol_qL[1];
4462  REAL dotqnL = (qxL*n1) + (qyL*n2);
4463 
4464  // Getting P solution for left side
4465 // REAL PseudoPressureL = sol_pL[0];
4466 
4467  // Getting another required data
4468 
4469  REAL TimeStep = this->fDeltaT;
4470  REAL Theta = this->fTheta;
4471 // REAL Gamma = this->fGamma;
4472  this->fnewWS=true;
4473 
4474  // Getting Harmonic mean of permeabilities
4475  int GeoIDLeft = dataleft[1].gelElId;
4476  TPZFMatrix<REAL> Kabsolute;
4477  TPZFMatrix<REAL> Kinverse;
4478  TPZFMatrix<REAL> Gfield;
4479 
4480  if (fYorN)
4481  {
4482  Kabsolute=this->fKabsoluteMap[GeoIDLeft];
4483  }
4484  else
4485  {
4486  this->K(Kabsolute);
4487  }
4488 
4489  Kinverse = this->Kinv(Kabsolute);
4490  Gfield = this->Gravity();
4491 
4492 // REAL GravityFluxL = (Kabsolute(0,0)*Gfield(0,0) + Kabsolute(0,1)*Gfield(1,0))*(n1)+
4493 // (Kabsolute(1,0)*Gfield(0,0) + Kabsolute(1,1)*Gfield(1,0))*(n2);
4494 
4495  REAL rockporosityl, oildensityl, waterdensityl;
4496  REAL drockporositydpl, doildensitydpl, dwaterdensitydpl;
4497 
4498  REAL oilviscosityl, waterviscosityl;
4499  REAL doilviscositydpl, dwaterviscositydpl;
4500 
4501  REAL bulklambdal, oillambdal, waterlambdal;
4502  REAL dbulklambdadpl, doillambdadpl, dwaterlambdadpl;
4503  REAL dbulklambdadsl, doillambdadsl, dwaterlambdadsl;
4504 
4505  REAL bulkfoill, bulkfwaterl;
4506  REAL dbulkfoildpl, dbulkfwaterdpl;
4507  REAL dbulkfoildsl, dbulkfwaterdsl;
4508 
4509 
4510  // Functions computed at point x_{k} for each integration point
4511  int VecPos= 0;
4512 
4513  this->Porosity(sol_pL[VecPos], rockporosityl, drockporositydpl);
4514  this->RhoOil(sol_pL[VecPos], oildensityl, doildensitydpl);
4515  this->RhoWater(sol_pL[VecPos], waterdensityl, dwaterdensitydpl);
4516  this->OilViscosity(sol_pL[VecPos], oilviscosityl, doilviscositydpl);
4517  this->WaterViscosity(sol_pL[VecPos], waterviscosityl, dwaterviscositydpl);
4518  this->OilLabmda(oillambdal, sol_pL[VecPos], sol_sL[VecPos], doillambdadpl, doillambdadsl);
4519  this->WaterLabmda(waterlambdal, sol_pL[VecPos], sol_sL[VecPos], dwaterlambdadpl, dwaterlambdadsl);
4520  this->Labmda(bulklambdal, sol_pL[VecPos], sol_sL[VecPos], dbulklambdadpl, dbulklambdadsl);
4521  this->fOil(bulkfoill, sol_pL[VecPos], sol_sL[VecPos], dbulkfoildpl, dbulkfoildsl);
4522  this->fWater(bulkfwaterl, sol_pL[VecPos], sol_sL[VecPos], dbulkfwaterdpl, dbulkfwaterdsl);
4523 
4524  STATE v2[6];
4525  v2[0] = bc.Val2()(0,0); // ux
4526  v2[1] = bc.Val2()(1,0); // uy
4527  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
4528  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
4529  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
4530  v2[5] = bc.Val2()(5,0); // Saturation
4531 // REAL qN = (v2[2]*n1 + v2[3]*n2); // Normal Flux
4532 
4533 
4534  // Upwind scheme
4535  // Third Vector Block (Equation three) Saturation equation
4536  REAL UpwindSaturation = 0.0;
4537 
4538  if (dotqnL > 0.0)
4539  {
4540 
4541  UpwindSaturation = bulkfwaterl;
4542  // Theta * TimeStep * Integrate[L L^{upwind} dot(v, n), Gamme_{e}] (Equation three) Bc-Left Part
4543  for(int isat=0; isat<SRowsleft; isat++)
4544  {
4545  for(int jsat=0; jsat<SRowsleft; jsat++)
4546  {
4547  ek(isat + FirstSL,jsat + FirstSL) -= (-1.0) * weight
4548  * (Theta) * (TimeStep) * phiSL(isat,0) * dbulkfwaterdsl * phiSL(jsat,0) * dotqnL;
4549  }
4550  }
4551  }
4552 
4553  else
4554 
4555  {
4556  UpwindSaturation = bulkfwaterl;
4557  if (dotqnL < 0.0 && fabs(dotqnL) > 1.0e-10) { std::cout << "Boundary condition error: inflow detected in outflow boundary condition: dotqnL = " << dotqnL << "\n";}
4558  }
4559 
4560  // Theta * TimeStep * Integrate[L S dot(v, n), Gamme_{e}] (Equation three) Right-Left Part
4561  for (int isat=0; isat < SRowsleft; isat++) {
4562 
4563  for (int jq=0; jq<QRowsleft; jq++)
4564  {
4565  int jLvectorindex = dataleft[1].fVecShapeIndex[jq].first;
4566  int jLshapeindex = dataleft[1].fVecShapeIndex[jq].second;
4567 
4568 
4569  REAL dotprodL =
4570  (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(0,jLvectorindex)) * (n1) +
4571  (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(1,jLvectorindex)) * (n2) ;//+
4572  // (phiQL(jLshapeindex,0)*dataleft[1].fNormalVec(2,jLvectorindex)) * (n3) ; // dot(q,n) left
4573 
4574  ek(isat + FirstSL,jq + FirstQL) -= (-1.0) * weight * (Theta) * (TimeStep) * phiSL(isat,0) * (UpwindSaturation) * dotprodL;
4575 
4576  }
4577  }
4578 
4579 
4580  // (Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Right-Left Part
4581  for(int isat=0; isat<SRowsleft; isat++)
4582  {
4583  REAL ResidualPart = (Theta) * (TimeStep) * ( (UpwindSaturation) * dotqnL );
4584  ef(isat + FirstSL) -= (-1.0) * weight * phiSL(isat,0) * ResidualPart;
4585  }
4586 
4587  }
4588 
4589  void TPZMultiphase::ApplySin (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
4590  {
4591 
4592  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
4593 // TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
4594  TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
4595  TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
4596 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
4597 
4598  int URowsleft = phiUL.Rows();
4599  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
4600  int PRowsleft = phiPL.Rows();
4601  int SRowsleft = phiSL.Rows();
4602 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
4603 
4604  int UStateVar = 2;
4605 
4606 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4607 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4608 
4609  int FirstUL = 0;
4610  int FirstQL = URowsleft * UStateVar + FirstUL;
4611  int FirstPL = QRowsleft + FirstQL;
4612  int FirstSL = PRowsleft + FirstPL;
4613 // int FirstQGL = SRowsleft + FirstSL;
4614 
4615  TPZManVector<REAL,3> &normal = data.normal;
4616  REAL n1 = normal[0];
4617  REAL n2 = normal[1];
4618 
4619  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
4620  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
4621  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
4622  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
4623 
4624  // Getting Q solution for left and right side
4625 // REAL uxL = sol_uL[0];
4626 // REAL uyL = sol_uL[1];
4627 // REAL qxL = sol_qL[0];
4628 // REAL qyL = sol_qL[1];
4629 // REAL dotqnL = (qxL*n1) + (qyL*n2);
4630 
4631  // Getting P solution for left side
4632 // REAL PseudoPressureL = sol_pL[0];
4633 
4634  // Getting another required data
4635 
4636  REAL TimeStep = this->fDeltaT;
4637  REAL Theta = this->fTheta;
4638 // REAL Gamma = this->fGamma;
4639  this->fnewWS=true;
4640 
4641  // Getting Harmonic mean of permeabilities
4642  int GeoIDLeft = dataleft[1].gelElId;
4643  TPZFMatrix<REAL> Kabsolute;
4644  TPZFMatrix<REAL> Kinverse;
4645  TPZFMatrix<REAL> Gfield;
4646 
4647  if (fYorN)
4648  {
4649  Kabsolute=this->fKabsoluteMap[GeoIDLeft];
4650  }
4651  else
4652  {
4653  this->K(Kabsolute);
4654  }
4655 
4656  Kinverse = this->Kinv(Kabsolute);
4657  Gfield = this->Gravity();
4658 
4659 // REAL GravityFluxL = (Kabsolute(0,0)*Gfield(0,0) + Kabsolute(0,1)*Gfield(1,0))*(n1)+
4660 // (Kabsolute(1,0)*Gfield(0,0) + Kabsolute(1,1)*Gfield(1,0))*(n2);
4661 
4662  REAL rockporosityl, oildensityl, waterdensityl;
4663  REAL drockporositydpl, doildensitydpl, dwaterdensitydpl;
4664 
4665  REAL oilviscosityl, waterviscosityl;
4666  REAL doilviscositydpl, dwaterviscositydpl;
4667 
4668  REAL bulklambdal, oillambdal, waterlambdal;
4669  REAL dbulklambdadpl, doillambdadpl, dwaterlambdadpl;
4670  REAL dbulklambdadsl, doillambdadsl, dwaterlambdadsl;
4671 
4672  REAL bulkfoill, bulkfwaterl;
4673  REAL dbulkfoildpl, dbulkfwaterdpl;
4674  REAL dbulkfoildsl, dbulkfwaterdsl;
4675 
4676 
4677  // Functions computed at point x_{k} for each integration point
4678  int VecPos= 0;
4679 
4680  this->Porosity(sol_pL[VecPos], rockporosityl, drockporositydpl);
4681  this->RhoOil(sol_pL[VecPos], oildensityl, doildensitydpl);
4682  this->RhoWater(sol_pL[VecPos], waterdensityl, dwaterdensitydpl);
4683  this->OilViscosity(sol_pL[VecPos], oilviscosityl, doilviscositydpl);
4684  this->WaterViscosity(sol_pL[VecPos], waterviscosityl, dwaterviscositydpl);
4685  this->OilLabmda(oillambdal, sol_pL[VecPos], sol_sL[VecPos], doillambdadpl, doillambdadsl);
4686  this->WaterLabmda(waterlambdal, sol_pL[VecPos], sol_sL[VecPos], dwaterlambdadpl, dwaterlambdadsl);
4687  this->Labmda(bulklambdal, sol_pL[VecPos], sol_sL[VecPos], dbulklambdadpl, dbulklambdadsl);
4688  this->fOil(bulkfoill, sol_pL[VecPos], sol_sL[VecPos], dbulkfoildpl, dbulkfoildsl);
4689  this->fWater(bulkfwaterl, sol_pL[VecPos], sol_sL[VecPos], dbulkfwaterdpl, dbulkfwaterdsl);
4690 
4691  STATE v2[6];
4692  v2[0] = bc.Val2()(0,0); // ux
4693  v2[1] = bc.Val2()(1,0); // uy
4694  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
4695  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
4696  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
4697  v2[5] = bc.Val2()(5,0); // Saturation
4698  REAL qN = (v2[2]*n1 + v2[3]*n2); // Normal Flux
4699 
4700  REAL UpwindSaturation = 0.0;
4701 
4702 //
4703 // if (dotqnL > 0.0)
4704 // {
4705 // this->fWater(bulkfwaterl, sol_pL[VecPos], v2[3], dbulkfwaterdpl, dbulkfwaterdsl);
4706 // UpwindSaturation = bulkfwaterl;
4707 // }
4708 //
4709 // else
4710 // {
4711  this->fWater(bulkfwaterl, sol_pL[VecPos], v2[5], dbulkfwaterdpl, dbulkfwaterdsl);
4712  UpwindSaturation = bulkfwaterl;
4713 
4714 //
4715 // }
4716 //
4717 
4718  // (Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Right-Left Part
4719  for(int isat=0; isat<SRowsleft; isat++)
4720  {
4721 
4722  REAL ResidualPart = (1-Theta) * (TimeStep) * ( (UpwindSaturation) * qN );
4723  ef(isat + FirstSL) -= (-1.0) * weight * phiSL(isat,0) * ResidualPart;
4724  }
4725 
4726  }
4727 
4728  void TPZMultiphase::ApplySout (TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
4729  {
4730 
4731  TPZFMatrix<REAL> &phiUL = dataleft[0].phi;
4732 // TPZFMatrix<REAL> &phiQL = dataleft[1].phi;
4733  TPZFMatrix<REAL> &phiPL = dataleft[2].phi;
4734  TPZFMatrix<REAL> &phiSL = dataleft[3].phi;
4735 // TPZFMatrix<REAL> &phiQGL = dataleft[4].phi;
4736 
4737  int URowsleft = phiUL.Rows();
4738  int QRowsleft = dataleft[1].fVecShapeIndex.NElements();
4739  int PRowsleft = phiPL.Rows();
4740  int SRowsleft = phiSL.Rows();
4741 // int QGRowsleft = dataleft[4].fVecShapeIndex.NElements();
4742 
4743  int UStateVar = 2;
4744 
4745 // int iRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4746 // int jRightInterfaceBlock = URowsleft * UStateVar + QRowsleft + PRowsleft + SRowsleft + QGRowsleft;
4747 
4748  int FirstUL = 0;
4749  int FirstQL = URowsleft * UStateVar + FirstUL;
4750  int FirstPL = QRowsleft + FirstQL;
4751  int FirstSL = PRowsleft + FirstPL;
4752 // int FirstQGL = SRowsleft + FirstSL;
4753 
4754  TPZManVector<REAL,3> &normal = data.normal;
4755  REAL n1 = normal[0];
4756  REAL n2 = normal[1];
4757 
4758  TPZManVector<STATE,3> sol_uL =dataleft[0].sol[0];
4759  TPZManVector<STATE,3> sol_qL =dataleft[1].sol[0];
4760  TPZManVector<STATE,3> sol_pL =dataleft[2].sol[0];
4761  TPZManVector<STATE,3> sol_sL =dataleft[3].sol[0];
4762 
4763  // Getting Q solution for left and right side
4764 // REAL uxL = sol_uL[0];
4765 // REAL uyL = sol_uL[1];
4766  REAL qxL = sol_qL[0];
4767  REAL qyL = sol_qL[1];
4768  REAL dotqnL = (qxL*n1) + (qyL*n2);
4769 
4770  // Getting P solution for left side
4771 // REAL PseudoPressureL = sol_pL[0];
4772 
4773  // Getting another required data
4774 
4775  REAL TimeStep = this->fDeltaT;
4776  REAL Theta = this->fTheta;
4777 // REAL Gamma = this->fGamma;
4778  this->fnewWS=true;
4779 
4780  // Getting Harmonic mean of permeabilities
4781  int GeoIDLeft = dataleft[1].gelElId;
4782  TPZFMatrix<REAL> Kabsolute;
4783  TPZFMatrix<REAL> Kinverse;
4784  TPZFMatrix<REAL> Gfield;
4785 
4786  if (fYorN)
4787  {
4788  Kabsolute=this->fKabsoluteMap[GeoIDLeft];
4789  }
4790  else
4791  {
4792  this->K(Kabsolute);
4793  }
4794 
4795  Kinverse = this->Kinv(Kabsolute);
4796  Gfield = this->Gravity();
4797 
4798 // REAL GravityFluxL = (Kabsolute(0,0)*Gfield(0,0) + Kabsolute(0,1)*Gfield(1,0))*(n1)+
4799 // (Kabsolute(1,0)*Gfield(0,0) + Kabsolute(1,1)*Gfield(1,0))*(n2);
4800 
4801  REAL rockporosityl, oildensityl, waterdensityl;
4802  REAL drockporositydpl, doildensitydpl, dwaterdensitydpl;
4803 
4804  REAL oilviscosityl, waterviscosityl;
4805  REAL doilviscositydpl, dwaterviscositydpl;
4806 
4807  REAL bulklambdal, oillambdal, waterlambdal;
4808  REAL dbulklambdadpl, doillambdadpl, dwaterlambdadpl;
4809  REAL dbulklambdadsl, doillambdadsl, dwaterlambdadsl;
4810 
4811  REAL bulkfoill, bulkfwaterl;
4812  REAL dbulkfoildpl, dbulkfwaterdpl;
4813  REAL dbulkfoildsl, dbulkfwaterdsl;
4814 
4815 
4816  // Functions computed at point x_{k} for each integration point
4817  int VecPos= 0;
4818 
4819  this->Porosity(sol_pL[VecPos], rockporosityl, drockporositydpl);
4820  this->RhoOil(sol_pL[VecPos], oildensityl, doildensitydpl);
4821  this->RhoWater(sol_pL[VecPos], waterdensityl, dwaterdensitydpl);
4822  this->OilViscosity(sol_pL[VecPos], oilviscosityl, doilviscositydpl);
4823  this->WaterViscosity(sol_pL[VecPos], waterviscosityl, dwaterviscositydpl);
4824  this->OilLabmda(oillambdal, sol_pL[VecPos], sol_sL[VecPos], doillambdadpl, doillambdadsl);
4825  this->WaterLabmda(waterlambdal, sol_pL[VecPos], sol_sL[VecPos], dwaterlambdadpl, dwaterlambdadsl);
4826  this->Labmda(bulklambdal, sol_pL[VecPos], sol_sL[VecPos], dbulklambdadpl, dbulklambdadsl);
4827  this->fOil(bulkfoill, sol_pL[VecPos], sol_sL[VecPos], dbulkfoildpl, dbulkfoildsl);
4828  this->fWater(bulkfwaterl, sol_pL[VecPos], sol_sL[VecPos], dbulkfwaterdpl, dbulkfwaterdsl);
4829 
4830  STATE v2[6];
4831  v2[0] = bc.Val2()(0,0); // ux
4832  v2[1] = bc.Val2()(1,0); // uy
4833  v2[2] = bc.Val2()(2,0)/this->fQref; // qx
4834  v2[3] = bc.Val2()(3,0)/this->fQref; // qy
4835  v2[4] = bc.Val2()(4,0)/this->fPref; // Pressure
4836  v2[5] = bc.Val2()(5,0); // Saturation
4837 // REAL qN = (v2[2]*n1 + v2[3]*n2); // Normal Flux
4838 
4839  REAL UpwindSaturation = 0.0;
4840 
4841  if (dotqnL > 0.0)
4842  {
4843  UpwindSaturation = bulkfwaterl;
4844  }
4845  else
4846  {
4847  UpwindSaturation = bulkfwaterl;
4848  if (dotqnL < 0.0 && fabs(dotqnL) > 1.0e-12) {std::cout << "Boundary condition error: inflow detected in outflow boundary condition: dotqnL = " << dotqnL << "\n";}
4849  }
4850 
4851  // (Theta) * deltat * Integrate[L*S dot(q,n), Gamma_{e}] (Equation three) Right-Left Part
4852  for(int isat=0; isat<SRowsleft; isat++)
4853  {
4854  REAL ResidualPart = (1-Theta) * (TimeStep) * ( (UpwindSaturation) * dotqnL );
4855  ef(isat + FirstSL) -= (-1.0) * weight * phiSL(isat,0) * ResidualPart;
4856  }
4857 
4858  }
4859 
4861 int TPZMultiphase::VariableIndex(const std::string &name){
4862  if(!strcmp("BulkVelocity",name.c_str())) return 1;
4863  if(!strcmp("WaterVelocity",name.c_str())) return 2;
4864  if(!strcmp("OilVelocity",name.c_str())) return 3;
4865  if(!strcmp("WeightedPressure",name.c_str())) return 4;
4866  if(!strcmp("WaterSaturation",name.c_str())) return 5;
4867  if(!strcmp("OilSaturation",name.c_str())) return 6;
4868  if(!strcmp("WaterDensity",name.c_str())) return 7;
4869  if(!strcmp("OilDensity",name.c_str())) return 8;
4870  if(!strcmp("RockPorosity",name.c_str())) return 9;
4871  if(!strcmp("GravityVelocity",name.c_str())) return 10;
4872  if(!strcmp("Kabsolute",name.c_str())) return 11;
4873  if(!strcmp("SwExact",name.c_str())) return 12;
4874  if(!strcmp("Displacement",name.c_str())) return 13;
4875  if(!strcmp("SigmaX",name.c_str())) return 14;
4876  if(!strcmp("SigmaY",name.c_str())) return 15;
4877 
4878  return TPZMaterial::VariableIndex(name);
4879 }
4880 
4881 int TPZMultiphase::NSolutionVariables(int var){
4882  if(var == 1) return 2;
4883  if(var == 2) return 2;
4884  if(var == 3) return 2;
4885  if(var == 4) return 1;
4886  if(var == 5) return 1;
4887  if(var == 6) return 1;
4888  if(var == 7) return 1;
4889  if(var == 8) return 1;
4890  if(var == 9) return 1;
4891  if(var == 10) return 2;
4892  if(var == 11) return 3;
4893  if(var == 12) return 1;
4894  if(var == 13) return 3;
4895  if(var == 14) return 1;
4896  if(var == 15) return 1;
4897 
4898  return TPZMaterial::NSolutionVariables(var);
4899 }
4900 
4901 void TPZMultiphase::Solution(TPZVec<TPZMaterialData> &datavec, int var, TPZVec<STATE> &Solout){
4902 
4903  Solout.Resize( this->NSolutionVariables(var));
4904  TPZVec<STATE> SolU, SolQ, SolP, SolS, SolQG, SolSExact(1);
4905  TPZFMatrix<STATE> SolSExactD;
4906  SolU = datavec[0].sol[0];
4907  SolQ = datavec[1].sol[0];
4908  SolP = datavec[2].sol[0];
4909  SolS = datavec[3].sol[0];
4910  SolQG = datavec[4].sol[0];
4911 
4912  TPZFNMatrix <6,STATE> dSolU;
4913  dSolU = datavec[0].dsol[0];
4914  TPZFNMatrix <9> axesU;
4915  axesU = datavec[0].axes;
4916 
4917  REAL epsx;
4918  REAL epsy;
4919  REAL epsxy;
4920  REAL SigX;
4921  REAL SigY;
4922  REAL SigZ;
4923  REAL Tau, DSolxy[2][2];
4924  REAL divu;
4925 
4926  if(var == 1){ //function (state variable Q)
4927  Solout[0] = SolQ[0];
4928  Solout[1] = SolQ[1];
4929  return;
4930  }
4931 
4932  if(var == 4){
4933  Solout[0] = SolP[0];//function (state variable p)
4934  return;
4935  }
4936 
4937  if(var == 5){
4938  Solout[0] = SolS[0];//function (state variable S)
4939  return;
4940  }
4941 
4942  if(var == 6){
4943  Solout[0] = 1.0-SolS[0];//function (state variable S)
4944  return;
4945  }
4946 
4947  if(var == 10){ //function (state variable QG)
4948  Solout[0] = SolQG[0];
4949  Solout[1] = SolQG[1];
4950  return;
4951  }
4952 
4953  if(var == 11){ //function (state variable QG)
4954 
4955  TPZFMatrix<REAL> Kabsolute;
4956  TPZFMatrix<REAL> Kinverse;
4957  TPZFMatrix<REAL> Gfield;
4958  if (fYorN)
4959  {
4960  int iel = datavec[0].gelElId;
4961  Kabsolute=this->fKabsoluteMap[iel];
4962  }
4963  else
4964  {
4965  this->K(Kabsolute);
4966  }
4967 
4968  Kinverse=this->Kinv(Kabsolute);
4969  Gfield = this->Gravity();
4970 
4971  Solout[0] = Kabsolute(0,0);
4972  Solout[1] = Kabsolute(1,1);
4973  Solout[2] = Kabsolute(2,2);
4974 
4975  return;
4976  }
4977 
4978 
4979  if(var == 12){ // ExactSolution
4980  fTimedependentFunctionExact->Execute(datavec[2].x, fTime, SolSExact,SolSExactD);
4981  Solout[0] = SolSExact[0];
4982  return;
4983  }
4984 
4985  if(var == 13){ //function (state variable U)
4986  Solout[0] = SolU[0];
4987  Solout[1] = SolU[1];
4988  Solout[2] = 0.0;
4989  return;
4990  }
4991 
4992 
4993 
4994  DSolxy[0][0] = dSolU(0,0)*axesU(0,0)+dSolU(1,0)*axesU(1,0); // dUx/dx
4995  DSolxy[1][0] = dSolU(0,0)*axesU(0,1)+dSolU(1,0)*axesU(1,1); // dUx/dy
4996 
4997  DSolxy[0][1] = dSolU(0,1)*axesU(0,0)+dSolU(1,1)*axesU(1,0); // dUy/dx
4998  DSolxy[1][1] = dSolU(0,1)*axesU(0,1)+dSolU(1,1)*axesU(1,1); // dUy/dy
4999 
5000  divu = DSolxy[0][0]+DSolxy[1][1]+0.0;
5001 
5002  REAL lambda,lambdau, mu, Balpha, Sestr;
5003  // Functions computed at point x_{k} for each integration point
5004  lambda = this->LameLambda();
5005  lambdau = this->LameLambdaU();
5006  mu = this->LameMu();
5007  Balpha = this->BiotAlpha();
5008  Sestr = this->Se();
5009 
5010  epsx = DSolxy[0][0];// du/dx
5011  epsy = DSolxy[1][1];// dv/dy
5012  epsxy = 0.5*(DSolxy[1][0]+DSolxy[0][1]);
5013  REAL C11 = 4*(mu)*(lambda+mu)/(lambda+2*mu);
5014  REAL C22 = 2*(mu)*(lambda)/(lambda+2*mu);
5015 
5016  if (this->fPlaneStress==1)
5017  {
5018  SigX = C11*epsx+C22*epsy;
5019  SigY = C11*epsy+C22*epsx;
5020  SigZ = 0.0;
5021  Tau = 2.0*mu*epsxy;
5022  }
5023  else
5024  {
5025  SigX = ((lambda + 2*mu)*(epsx) + (lambda)*epsy - Balpha*SolP[0]);
5026  SigY = ((lambda + 2*mu)*(epsy) + (lambda)*epsx - Balpha*SolP[0]);
5027  SigZ = lambda*divu - Balpha*SolP[0];
5028  Tau = 2.0*mu*epsxy;
5029  }
5030 
5031  if(var == 14){ // (SigmaX)
5032  Solout[0] = SigX;
5033  return;
5034  }
5035 
5036  if(var == 15){ // (SigmaY)
5037  Solout[0] = SigY;
5038  return;
5039  }
5040 
5041 }
5042 
5043 
5044 void TPZMultiphase::FillDataRequirements(TPZVec<TPZMaterialData > &datavec)
5045 {
5046  int nref = datavec.size();
5047  for(int i = 0; i<nref; i++ )
5048  {
5049  datavec[i].SetAllRequirements(false);
5050  datavec[i].fNeedsSol = true;
5051  datavec[i].fNeedsNeighborSol = true;
5052  datavec[i].fNeedsNeighborCenter = false;
5053  datavec[i].fNeedsNormal = true;
5054  }
5055 
5056 }
5057 
5058 void TPZMultiphase::FillBoundaryConditionDataRequirement(int type,TPZVec<TPZMaterialData > &datavec){
5059  int nref = datavec.size();
5060  for(int i = 0; i<nref; i++)
5061  {
5062  datavec[i].fNeedsSol = true;
5063  datavec[i].fNeedsNormal = true;
5064  datavec[i].fNeedsNeighborSol = true;
5065  }
5066 }
5067 
5068 int TPZMultiphase::ClassId() const{
5069  return Hash("TPZMultiphase") ^ TPZDiscontinuousGalerkin::ClassId() << 1;
5070 }
TPZDiscontinuousGalerkin
Defines the interface which material objects need to implement for discontinuous Galerkin formulation...
Definition: pzdiscgal.h:20
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZMultiphase::OilLabmda
void OilLabmda(REAL &OilLabmda, REAL Po, REAL Sw, REAL &dOilLabmdaDPo, REAL &dOilLabmdaDSw)
Oil mobility. .
Definition: pzmultiphase.cpp:275
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZMultiphase::Labmda
void Labmda(REAL &Labmda, REAL Pw, REAL Sw, REAL &dLabmdaDPw, REAL &dLabmdaDSw)
Bulk mobility. .
Definition: pzmultiphase.cpp:308
TPZMultiphase::fGamma
REAL fGamma
Parameter representing temporal scheme for conservation equation.
Definition: pzmultiphase.h:180
TPZMultiphase::Contribute
virtual void Contribute(TPZVec< TPZMaterialData > &datavec, REAL weight, TPZFMatrix< STATE > &ef) override
It computes a contribution to the stiffness matrix and load vector at one integration point to multip...
Definition: pzmultiphase.cpp:1572
TPZDiscontinuousGalerkin::ClassId
int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzdiscgal.cpp:110
TPZMaterialData::normal
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
Definition: pzmaterialdata.h:69
TPZMultiphase::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: pzmultiphase.cpp:75
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZMultiphase::MatId
virtual int MatId()
Definition: pzmultiphase.cpp:73
TPZMultiphase::Name
virtual std::string Name() override
Returns the name of the material.
Definition: pzmultiphase.h:73
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZMultiphase::OilViscosity
void OilViscosity(REAL po, REAL &OilViscosity, REAL &dOilViscosityDpo)
Oil viscosity. .
Definition: pzmultiphase.cpp:261
TPZMaterial::fBCForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fBCForcingFunction
Pointer to bc forcing function, it is a variable boundary condition at differential equation...
Definition: TPZMaterial.h:59
TPZMultiphase::RhoWaterSC
REAL RhoWaterSC()
Water density on standard conditions - kg/m3.
Definition: pzmultiphase.cpp:591
TPZMultiphase::FillDataRequirements
virtual void FillDataRequirements(TPZVec< TPZMaterialData > &datavec) override
Fill material data parameter with necessary requirements for the Contribute method. Here, in base class, all requirements are considered as necessary. Each derived class may optimize performance by selecting only the necessary data.
Definition: pzmultiphase.cpp:5044
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZMultiphase::CapillaryPressure
void CapillaryPressure(REAL So, REAL &pc, REAL &DpcDSo)
Capilar pressure. .
Definition: pzmultiphase.cpp:141
TPZMultiphase::WaterLabmda
void WaterLabmda(REAL &WaterLabmda, REAL Pw, REAL Sw, REAL &dWaterLabmdaDPw, REAL &dWaterLabmdaDSw)
Water mobility. .
Definition: pzmultiphase.cpp:291
TPZMaterial::VariableIndex
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Definition: TPZMaterial.cpp:141
TPZMultiphase::Kro
void Kro(REAL Sw, REAL &Kro, REAL &dKroDSw)
Oil relative permeability. .
Definition: pzmultiphase.cpp:154
TPZMultiphase::Dimension
virtual int Dimension() const override
Returns the integrable dimension of the material.
Definition: pzmultiphase.cpp:71
TPZMultiphase::RhoOil
void RhoOil(REAL po, REAL &RhoOil, REAL &dRhoOilDpo)
Definition: pzmultiphase.cpp:244
TPZMultiphase::fOil
void fOil(REAL &fOil, REAL Pw, REAL Sw, REAL &dfOilDPw, REAL &dfOilDSw)
Fractional oil flux. .
Definition: pzmultiphase.cpp:325
TPZMultiphase::BiotAlpha
STATE BiotAlpha()
Biot parameter Parameter. .
Definition: pzmultiphase.cpp:123
TPZMultiphase::fRhoref
REAL fRhoref
density reference - kg/m3
Definition: pzmultiphase.h:155
TPZMultiphase::fKabsoluteMap
TPZStack< TPZFMatrix< REAL > > fKabsoluteMap
K map.
Definition: pzmultiphase.h:162
TPZMultiphase::WaterViscosity
void WaterViscosity(REAL po, REAL &WaterViscosity, REAL &dWaterViscosityDpo)
Water viscosity. .
Definition: pzmultiphase.cpp:268
TPZMultiphase::ApplyUyD
virtual void ApplyUyD(TPZMaterialData &data, TPZVec< TPZMaterialData > &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzmultiphase.cpp:3946
TPZMultiphase::fstar
void fstar(REAL &fStar, REAL Pw, REAL Sw, REAL Gdotn, REAL &dfstarDPw, REAL &dfstarDSw)
Fractional product upwind function, Gdotn > 0 means water decrease (dSw/dt < 0), Gdotn < 0 means wate...
Definition: pzmultiphase.cpp:437
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
TPZMultiphase::fQref
REAL fQref
Pressure reference - Pa.
Definition: pzmultiphase.h:149
TPZMultiphase::ApplySin
virtual void ApplySin(TPZMaterialData &data, TPZVec< TPZMaterialData > &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzmultiphase.cpp:4254
TPZMultiphase::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzmultiphase.cpp:5068
TPZBndCond::Val2
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZMultiphase::fDim
int fDim
Problem dimension.
Definition: pzmultiphase.h:35
TPZMultiphase::fLref
REAL fLref
Characteristic length - m.
Definition: pzmultiphase.h:143
TPZMultiphase::fmatId
int fmatId
Material id.
Definition: pzmultiphase.h:38
TPZMultiphase::ApplySout
virtual void ApplySout(TPZMaterialData &data, TPZVec< TPZMaterialData > &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzmultiphase.cpp:4422
TPZMaterial::Print
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
Definition: TPZMaterial.cpp:103
TPZMultiphase::ApplySigmaN
virtual void ApplySigmaN(TPZMaterialData &data, TPZVec< TPZMaterialData > &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzmultiphase.cpp:4007
TPZMultiphase::LameMu
STATE LameMu()
Lame Second Parameter. .
Definition: pzmultiphase.cpp:113
TPZMultiphase::fYorN
bool fYorN
Use or not K map.
Definition: pzmultiphase.h:168
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZMultiphase::Solution
virtual void Solution(TPZVec< TPZMaterialData > &datavec, int var, TPZVec< STATE > &Solout) override
Returns the solution associated with the var index based on the finite element approximation.
Definition: pzmultiphase.cpp:4901
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
TPZMultiphase::fEtaref
REAL fEtaref
viscosity reference - Pa s
Definition: pzmultiphase.h:158
TPZMultiphase::ApplyPN
virtual void ApplyPN(TPZMaterialData &data, TPZVec< TPZMaterialData > &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzmultiphase.cpp:4135
TPZMultiphase::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: pzmultiphase.cpp:4861
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
stats.read
def read(filename)
Definition: stats.py:13
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZMultiphase::fPref
REAL fPref
Pressure reference - Pa.
Definition: pzmultiphase.h:152
TPZMultiphase::fWater
void fWater(REAL &fWater, REAL Pw, REAL Sw, REAL &dfWaterDPw, REAL &dfWaterDSw)
Fractional water flux. .
Definition: pzmultiphase.cpp:348
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZMultiphase::ContributeBCInterface
virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to residual vector at one BC integration point.
Definition: pzmultiphase.cpp:3303
TPZMultiphase
Material to solve a 2d multiphase transport problem by multiphysics simulation.
Definition: pzmultiphase.h:30
TPZMultiphase::fDeltaT
REAL fDeltaT
Simulation time step.
Definition: pzmultiphase.h:171
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZMultiphase::SetKMap
void SetKMap(TPZStack< TPZFMatrix< REAL > > KabsoluteMap)
Definition: pzmultiphase.h:237
TPZMultiphase::Se
STATE Se()
//Se o 1/M coeficiente poroelastico de armazenamento a volume constante.
Definition: pzmultiphase.cpp:133
TPZMultiphase::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: pzmultiphase.cpp:77
TPZMaterial::gBigNumber
static REAL gBigNumber
Big number to penalization method, used for Dirichlet conditions.
Definition: TPZMaterial.h:83
TPZMultiphase::ApplyUxD
virtual void ApplyUxD(TPZMaterialData &data, TPZVec< TPZMaterialData > &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzmultiphase.cpp:3884
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZMaterial::fTimedependentFunctionExact
TPZAutoPointer< TPZFunction< STATE > > fTimedependentFunctionExact
Pointer to time dependent exact solution function, needed to calculate exact error.
Definition: TPZMaterial.h:56
TPZMultiphase::LameLambdaU
STATE LameLambdaU()
Undrained Lame First Parameter. .
Definition: pzmultiphase.cpp:103
TPZMultiphase::fOilstar
void fOilstar(REAL &fOstar, REAL Pw, REAL Sw, REAL &dfOstarDPw, REAL &dfOstarDSw)
Fractional product function, oil decrease direction (dSw/dt > 0). .
Definition: pzmultiphase.cpp:406
TPZMultiphase::Porosity
void Porosity(REAL po, REAL &poros, REAL &dPorosDp)
Rock porosity. .
Definition: pzmultiphase.cpp:234
TPZMultiphase::fWaterstar
void fWaterstar(REAL &fWstar, REAL Pw, REAL Sw, REAL &dfWstarDPw, REAL &dfWstarDSw)
Fractional product function, water decrease direction (dSw/dt < 0). .
Definition: pzmultiphase.cpp:373
TPZMultiphase::fKref
REAL fKref
Permeability reference - m2.
Definition: pzmultiphase.h:146
TPZMultiphase::~TPZMultiphase
virtual ~TPZMultiphase()
Definition: pzmultiphase.cpp:67
TPZMultiphase::ContributeInterface
virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
It computes a contribution to stiffness matrix and load vector at one integration point...
Definition: pzmultiphase.cpp:2033
log
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ log
Definition: tfadfunc.h:130
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZMultiphase::RhoOilSC
REAL RhoOilSC()
Oil density on standard conditions - kg/m3.
Definition: pzmultiphase.cpp:586
TPZStack
This class implements a stack object. Utility.
Definition: pzcheckmesh.h:14
TPZMultiphase::Krw
void Krw(REAL Sw, REAL &Krw, REAL &dKrwSo)
Water relative permeability. .
Definition: pzmultiphase.cpp:176
TPZMultiphase::FillBoundaryConditionDataRequirement
virtual void FillBoundaryConditionDataRequirement(int type, TPZVec< TPZMaterialData > &datavec) override
This method defines which parameters need to be initialized in order to compute the contribution of t...
Definition: pzmultiphase.cpp:5058
pzaxestools.h
Contains declaration of the TPZAxesTools class which implements verifications over axes...
TPZMultiphase::ContributeBC
virtual void ContributeBC(TPZVec< TPZMaterialData > &datavec, REAL weight, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to the stiffness matrix and load vector at one BC integration point to mul...
Definition: pzmultiphase.cpp:3277
exp
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ exp
Definition: tfadfunc.h:125
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZMultiphase::ApplyQnD
virtual void ApplyQnD(TPZMaterialData &data, TPZVec< TPZMaterialData > &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzmultiphase.cpp:4061
TPZMultiphase::LoadKMap
void LoadKMap(std::string MaptoRead)
Absolute permeability.
Definition: pzmultiphase.cpp:631
TPZMultiphase::K
void K(TPZFMatrix< REAL > &Kab)
Absolute permeability.
Definition: pzmultiphase.cpp:606
TPZMultiphase::fTheta
REAL fTheta
Parameter representing temporal scheme for transport equation.
Definition: pzmultiphase.h:177
TPZMultiphase::fTime
REAL fTime
Simulation current time.
Definition: pzmultiphase.h:174
TPZStack::push_back
void push_back(const T object)
Definition: pzstack.h:48
TPZMultiphase::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: pzmultiphase.cpp:4881
TPZMultiphase::LameLambda
STATE LameLambda()
Lame First Parameter. .
Definition: pzmultiphase.cpp:93
TPZMultiphase::SWaterstar
void SWaterstar(REAL &Swstar, REAL &Po, REAL &Sw)
Water saturation maximum value of the fractional flow product function. .
Definition: pzmultiphase.cpp:197
TPZMultiphase::Gravity
TPZFMatrix< REAL > Gravity()
Gravity.
Definition: pzmultiphase.cpp:596
TPZMultiphase::Kinv
TPZFMatrix< REAL > Kinv(TPZFMatrix< REAL > &Kab)
Absolute permeability inverse.
Definition: pzmultiphase.cpp:617
TPZMultiphase::fxi
REAL fxi
Big number balance.
Definition: pzmultiphase.h:44
TPZMultiphase::fPlaneStress
int fPlaneStress
plane stress condition
Definition: pzmultiphase.h:165
TPZMultiphase::RhoWater
void RhoWater(REAL pw, REAL &RhoWater, REAL &dRhoWaterDpo)
Definition: pzmultiphase.cpp:252
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13