NeoPZ
mixedpoisson.cpp
Go to the documentation of this file.
1 
8 #include "mixedpoisson.h"
9 #include "pzlog.h"
10 #include "pzbndcond.h"
11 #include "pzfmatrix.h"
12 #include "pzaxestools.h"
13 
14 
15 #include <iostream>
16 
17 #ifdef LOG4CXX
18 static LoggerPtr logdata(Logger::getLogger("pz.mixedpoisson.data"));
19 static LoggerPtr logerror(Logger::getLogger("pz.mixedpoisson.error"));
20 #endif
21 
22 TPZMixedPoisson::TPZMixedPoisson(): TPZRegisterClassId(&TPZMixedPoisson::ClassId), TPZMatPoisson3d() {
23  fvisc = 1.;
24  ff = 0.;
25  fIsStabilized = false;
26  fdelta1 = 0.;
27  fdelta2 = 0.;
28  fUseHdois = false;
29  fh2 = 1.;
30 
31  fInvK.Redim(3, 3);
32  fTensorK.Redim(3, 3);
33  fTensorK.Identity();
34  fInvK.Identity();
35  fPermeabilityFunction = NULL;
36 }
37 
38 TPZMixedPoisson::TPZMixedPoisson(int matid, int dim): TPZRegisterClassId(&TPZMixedPoisson::ClassId), TPZMatPoisson3d(matid,dim) {
39  if (dim < 1) {
40  DebugStop();
41  }
42  fvisc = 1.;
43  ff = 0.;
44  fIsStabilized = false;
45  fdelta1 = 0.;
46  fdelta2 = 0.;
47  fUseHdois = false;
48  fh2 = 1.;
49 
50  fInvK.Redim(3, 3);
51  fTensorK.Resize(3, 3);
52  fInvK.Identity();
53  fTensorK.Identity();
54  fPermeabilityFunction = NULL;
55 }
56 
57 TPZMixedPoisson::~TPZMixedPoisson() {
58 }
59 
60 TPZMixedPoisson::TPZMixedPoisson(const TPZMixedPoisson &cp) :TPZRegisterClassId(&TPZMixedPoisson::ClassId), TPZMatPoisson3d(cp) {
61  fvisc = cp.fvisc;
62  ff = cp.ff;
63  fIsStabilized = cp.fIsStabilized;
64  fdelta1 = cp.fdelta1;
65  fdelta2 = cp.fdelta2;
66  fUseHdois = cp.fUseHdois;
67  fh2 = cp.fh2;
68 
69  fInvK = cp.fInvK;
70  fTensorK = cp.fTensorK;
71  fPermeabilityFunction = cp.fPermeabilityFunction;
72 }
73 
74 TPZMixedPoisson & TPZMixedPoisson::operator=(const TPZMixedPoisson &copy){
75  TPZMatPoisson3d::operator=(copy);
76  fvisc = copy.fvisc;
77  ff = copy.ff;
78  fIsStabilized = copy.fIsStabilized;
79  fdelta1 = copy.fdelta1;
80  fdelta2 = copy.fdelta2;
81  fUseHdois = copy.fUseHdois;
82  fh2 = copy.fh2;
83 
84  fInvK = copy.fInvK;
85  fTensorK = copy.fTensorK;
86  fPermeabilityFunction = copy.fPermeabilityFunction;
87  return *this;
88 }
89 
90 int TPZMixedPoisson::NStateVariables() const {
91  return 1;
92 }
93 
94 void TPZMixedPoisson::Print(std::ostream &out) {
95  out << "name of material : " << Name() << "\n";
96  out << "Base Class properties :";
97  TPZMatPoisson3d::Print(out);
98  out << "\n";
99 }
100 
102 void TPZMixedPoisson::GetPermeability(TPZVec<REAL> &x, TPZFMatrix<REAL> &PermTensor, TPZFMatrix<REAL> &InvPermTensor)
103 {
104 
105 
106  PermTensor = this->fTensorK;
107  InvPermTensor = this->fInvK;
108  //int rtens = 2*fDim;
109  if(fPermeabilityFunction){
110  PermTensor.Zero();
111  InvPermTensor.Zero();
112  TPZFNMatrix<18,STATE> resultMat(2*3,3,0.);
113  TPZManVector<STATE> res;
114  fPermeabilityFunction->Execute(x,res,resultMat);
115 
116  for(int id=0; id<3; id++){
117  for(int jd=0; jd<3; jd++){
118 
119  PermTensor(id,jd) = resultMat(id,jd);
120  InvPermTensor(id,jd) = resultMat(id+fDim,jd);
121  }
122  }
123  }
124 }
125 
126 
127 void TPZMixedPoisson::Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef) {
128 
129 //#ifdef PZDEBUG
130 // int nref = datavec.size();
131 //
132 // if (nref != 2) {
133 // std::cout << " Error. The size of the datavec is different from 2." << std::endl;
134 // DebugStop();
135 // }
136 //#endif
137 
138 
139 
140  STATE force = ff;
141  if(fForcingFunction) {
142  TPZManVector<STATE> res(1);
143  fForcingFunction->Execute(datavec[1].x,res);
144  force = res[0];
145  }
146 
147  TPZFNMatrix<9,REAL> PermTensor;
148  TPZFNMatrix<9,REAL> InvPermTensor;
149 
150  GetPermeability(datavec[1].x, PermTensor, InvPermTensor);
151 
152  // Setting the phis
153  TPZFMatrix<REAL> &phiQ = datavec[0].phi;
154  TPZFMatrix<REAL> &phip = datavec[1].phi;
155  TPZFMatrix<REAL> &dphiQ = datavec[0].dphix;
156  TPZFMatrix<REAL> &dphiP = datavec[1].dphix;
157  TPZFNMatrix<9,REAL> dphiPXY(3,dphiP.Cols());
158  TPZAxesTools<REAL>::Axes2XYZ(dphiP, dphiPXY, datavec[1].axes);
159 
160 
161  REAL &faceSize = datavec[0].HSize;
162  if(fUseHdois==true){
163  fh2 = faceSize*faceSize;
164  }else fh2 = 1.;
165 
166  int phrq, phrp;
167  phrp = phip.Rows();
168  phrq = datavec[0].fVecShapeIndex.NElements();
169 
170  int nactive = 0;
171  for (int i=0; i<datavec.size(); i++) {
172  if (datavec[i].fActiveApproxSpace) {
173  nactive++;
174  }
175  }
176 #ifdef PZDEBUG
177  if(nactive == 4)
178  {
179  int phrgb = datavec[2].phi.Rows();
180  int phrub = datavec[3].phi.Rows();
181  if(phrp+phrq+phrgb+phrub != ek.Rows())
182  {
183  DebugStop();
184  }
185  }else
186  {
187  if(phrp+phrq != ek.Rows())
188  {
189  DebugStop();
190  }
191  }
192 #endif
193  //Calculate the matrix contribution for flux. Matrix A
194  for(int iq=0; iq<phrq; iq++)
195  {
196  //ef(iq, 0) += 0.;
197  int ivecind = datavec[0].fVecShapeIndex[iq].first;
198  int ishapeind = datavec[0].fVecShapeIndex[iq].second;
199  TPZFNMatrix<3,REAL> ivec(3,1,0.);
200  for(int id=0; id<3; id++){
201  ivec(id,0) = datavec[0].fNormalVec(id,ivecind);
202  }
203 
204  //Inserindo termo de estabilizacao no termo de fonte
205  REAL divqi = 0.;
206  if(fIsStabilized)
207  {
208  //calculando div(qi)
209  TPZFNMatrix<3,REAL> axesvec(3,1,0.);
210  datavec[0].axes.Multiply(ivec,axesvec);
211  for(int iloc=0; iloc<fDim; iloc++)
212  {
213  divqi += axesvec(iloc,0)*dphiQ(iloc,ishapeind);
214  }
215  ef(iq, 0) += weight*(fdelta2*divqi*force);
216  }
217 
218  TPZFNMatrix<3,REAL> ivecZ(3,1,0.);
219  TPZFNMatrix<3,REAL> jvecZ(3,1,0.);
220  for (int jq=0; jq<phrq; jq++)
221  {
222  TPZFNMatrix<3,REAL> jvec(3,1,0.);
223  int jvecind = datavec[0].fVecShapeIndex[jq].first;
224  int jshapeind = datavec[0].fVecShapeIndex[jq].second;
225 
226  for(int id=0; id<3; id++){
227  jvec(id,0) = datavec[0].fNormalVec(id,jvecind);
228  }
229 
230  //dot product between Kinv[u]v
231  jvecZ.Zero();
232  for(int id=0; id<fDim; id++){
233  for(int jd=0; jd<fDim; jd++){
234  jvecZ(id,0) += InvPermTensor(id,jd)*jvec(jd,0);
235  }
236  }
237  //jvecZ.Print("mat1 = ");
238  REAL prod1 = ivec(0,0)*jvecZ(0,0) + ivec(1,0)*jvecZ(1,0) + ivec(2,0)*jvecZ(2,0);
239  ek(iq,jq) += fvisc*weight*phiQ(ishapeind,0)*phiQ(jshapeind,0)*prod1;
240 
241 
242  //Inserindo termos de estabilizacao na matriz do fluxo
243  if(fIsStabilized==true)
244  {
245  //termos de delta1
246  //dot product between uKinv[v]
247  ivecZ.Zero();
248  for(int id=0; id<fDim; id++){
249  for(int jd=0; jd<fDim; jd++){
250  ivecZ(id,0) += InvPermTensor(id,jd)*ivec(jd,0);
251  }
252  }
253  //ivecZ.Print("mat2 = ");
254  REAL prod2 = ivecZ(0,0)*jvec(0,0) + ivecZ(1,0)*jvec(1,0) + ivecZ(2,0)*jvec(2,0);
255  ek(iq,jq) += (-1.)*weight*fh2*fdelta1*fvisc*phiQ(ishapeind,0)*phiQ(jshapeind,0)*prod2;
256 
257 
258  //termos de delta2: //dot product between divQu.divQv
259  REAL divqj = 0.;
260  TPZFNMatrix<3,REAL> axesvec(3,1,0.);
261  datavec[0].axes.Multiply(jvec,axesvec);
262  //calculando div(qj)
263  for(int jloc=0; jloc<fDim; jloc++)
264  {
265  divqj += axesvec(jloc,0)*dphiQ(jloc,jshapeind);
266  }
267  ek(iq,jq) += weight*fdelta2*divqi*divqj;
268  }
269 
270  }
271  }
272 
273 
274  // Coupling terms between flux and pressure. Matrix B
275  for(int iq=0; iq<phrq; iq++)
276  {
277  int ivecind = datavec[0].fVecShapeIndex[iq].first;
278  int ishapeind = datavec[0].fVecShapeIndex[iq].second;
279 
280  TPZFNMatrix<3,REAL> ivec(3,1,0.);
281  for(int id=0; id<3; id++){
282  ivec(id,0) = datavec[0].fNormalVec(id,ivecind);
283  //ivec(1,0) = datavec[0].fNormalVec(1,ivecind);
284  //ivec(2,0) = datavec[0].fNormalVec(2,ivecind);
285  }
286  TPZFNMatrix<3,REAL> axesvec(3,1,0.);
287  datavec[0].axes.Multiply(ivec,axesvec);
288 
289  REAL divwq = 0.;
290  for(int iloc=0; iloc<fDim; iloc++)
291  {
292  divwq += axesvec(iloc,0)*dphiQ(iloc,ishapeind);
293  }
294  for (int jp=0; jp<phrp; jp++) {
295 
296  REAL fact = (-1.)*weight*phip(jp,0)*divwq;
297  // Matrix B
298  ek(iq, phrq+jp) += fact;
299 
300  // Matrix B^T
301  ek(phrq+jp,iq) += fact;
302 
303 
304  //Inserindo termo de estabilizacao: delta1
305  if(fIsStabilized==true)
306  {
307  //produto gardPu.Qv
308  REAL dotVGradP = 0.;
309 
310  for(int k =0; k<fDim; k++)
311  {
312  dotVGradP += ivec(k,0)*phiQ(ishapeind,0)*dphiP(k,jp);
313  }
314 
315  REAL integration = (-1.)*weight*fh2*fdelta1*dotVGradP;
316 
317  // Estabilizacao delta1 na Matrix B
318  ek(iq, phrq+jp) += integration;
319 
320  // Estabilizacao delta1 na Matrix BˆT
321  ek(phrq+jp,iq) += integration;
322  }
323  }
324  }
325 
326  //termo fonte referente a equacao da pressao
327  for(int ip=0; ip<phrp; ip++){
328  ef(phrq+ip,0) += (-1.)*weight*force*phip(ip,0);
329  }
330 
331  //Contribution for estabilization delta1 for gradPu*gradPv. Matrix D
332  if(fIsStabilized==true)
333  {
334  //produto KgradPu x KgradPv
335  TPZFNMatrix<3,REAL> dphiPuZ(dphiP.Rows(),dphiP.Cols(),0.);
336  PermTensor.Multiply(dphiPXY, dphiPuZ);
337 
338  REAL umSobreVisc = 1./fvisc;
339  for(int ip=0; ip<phrp; ip++)
340  {
341  for(int jp=0; jp<phrp; jp++)
342  {
343  for(int k =0; k<3; k++)
344  {
345  ek(phrq+ip, phrq+jp) += (-1.)*weight*fh2*fdelta1*umSobreVisc*dphiPuZ(k,ip)*dphiPXY(k,jp);
346  }
347 
348  }
349  }
350  }
351  if(nactive == 4)
352  {
353  for(int ip=0; ip<phrp; ip++)
354  {
355  ek(phrq+ip,phrq+phrp) += phip(ip,0)*weight;
356  ek(phrq+phrp,phrq+ip) += phip(ip,0)*weight;
357  }
358  ek(phrp+phrq+1,phrq+phrp) += -weight;
359  ek(phrq+phrp,phrp+phrq+1) += -weight;
360  }
361  //
362 // #ifdef LOG4CXX
363 // if(logdata->isDebugEnabled())
364 // {
365 // std::stringstream sout;
366 // sout<<"\n\n Matriz ek e vetor fk \n ";
367 // ek.Print("ekmph = ",sout,EMathematicaInput);
368 // ef.Print("efmph = ",sout,EMathematicaInput);
369 // LOGPZ_DEBUG(logdata,sout.str());
370 // }
371 // #endif
372 
373 }
374 
375 
376 //void TPZMixedPoisson::Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef) {
377 //
378 //#ifdef PZDEBUG
379 // int nref = datavec.size();
380 // if (nref != 2 ) {
381 // std::cout << " Erro. The size of the datavec is different from 2 \n";
382 // DebugStop();
383 // }
384 //#endif
385 //
386 // STATE force = ff;
387 // if(fForcingFunction) {
388 // TPZManVector<STATE> res(1);
389 // fForcingFunction->Execute(datavec[1].x,res);
390 // force = res[0];
391 // }
392 //
393 // // Setting the phis
394 // TPZFMatrix<REAL> &phiQ = datavec[0].phi;
395 // TPZFMatrix<REAL> &phip = datavec[1].phi;
396 // TPZFMatrix<REAL> &dphiQ = datavec[0].dphix;
397 // TPZFMatrix<REAL> &dphiP = datavec[1].dphix;
398 //
399 // REAL &faceSize = datavec[0].HSize;
400 // if(fUseHdois==true){
401 // fh2 = faceSize*faceSize;
402 // }else fh2 = 1.;
403 //
404 // int phrq, phrp;
405 // phrp = phip.Rows();
406 // phrq = datavec[0].fVecShapeIndex.NElements();
407 //
408 // //Calculate the matrix contribution for flux. Matrix A
409 // REAL InvK = fvisc/fk;
410 // for(int iq=0; iq<phrq; iq++)
411 // {
412 // //ef(iq, 0) += 0.;
413 // int ivecind = datavec[0].fVecShapeIndex[iq].first;
414 // int ishapeind = datavec[0].fVecShapeIndex[iq].second;
415 // TPZFNMatrix<3> ivec(3,1);
416 // ivec(0,0) = datavec[0].fNormalVec(0,ivecind);
417 // ivec(1,0) = datavec[0].fNormalVec(1,ivecind);
418 // ivec(2,0) = datavec[0].fNormalVec(2,ivecind);
419 //
420 // //Inserindo termo de estabilizacao no termo de fonte
421 // REAL divqi = 0.;
422 // if(fIsStabilized==true)
423 // {
424 // //calculando div(qi)
425 // TPZFNMatrix<3> axesvec(3,1);
426 // datavec[0].axes.Multiply(ivec,axesvec);
427 // for(int iloc=0; iloc<fDim; iloc++)
428 // {
429 // divqi += axesvec(iloc,0)*dphiQ(iloc,ishapeind);
430 // }
431 // ef(iq, 0) += weight*(fdelta2*divqi*force);
432 // }
433 //
434 // for (int jq=0; jq<phrq; jq++)
435 // {
436 // TPZFNMatrix<3> jvec(3,1);
437 // int jvecind = datavec[0].fVecShapeIndex[jq].first;
438 // int jshapeind = datavec[0].fVecShapeIndex[jq].second;
439 // jvec(0,0) = datavec[0].fNormalVec(0,jvecind);
440 // jvec(1,0) = datavec[0].fNormalVec(1,jvecind);
441 // jvec(2,0) = datavec[0].fNormalVec(2,jvecind);
442 //
443 // //dot product between u and v
444 // REAL prod = ivec(0,0)*jvec(0,0) + ivec(1,0)*jvec(1,0) + ivec(2,0)*jvec(2,0);
445 // ek(iq,jq) += InvK*weight*phiQ(ishapeind,0)*phiQ(jshapeind,0)*prod;
446 //
447 //
448 // //Inserindo termos de estabilizacao na matriz do fluxo
449 // if(fIsStabilized==true)
450 // {
451 // //termos de delta1
452 // ek(iq,jq) += (-1.)*weight*fh2*fdelta1*InvK*phiQ(ishapeind,0)*phiQ(jshapeind,0)*prod;
453 //
454 //
455 // //termos de delta2
456 // REAL divqj = 0.;
457 // TPZFNMatrix<3> axesvec(3,1);
458 // datavec[0].axes.Multiply(jvec,axesvec);
459 // //calculando div(qj)
460 // for(int jloc=0; jloc<fDim; jloc++)
461 // {
462 // divqj += axesvec(jloc,0)*dphiQ(jloc,jshapeind);
463 // }
464 // ek(iq,jq) += weight*fdelta2*divqi*divqj;
465 // }
466 //
467 // }
468 // }
469 //
470 //
471 // // Coupling terms between flux and pressure. Matrix B
472 // for(int iq=0; iq<phrq; iq++)
473 // {
474 // int ivecind = datavec[0].fVecShapeIndex[iq].first;
475 // int ishapeind = datavec[0].fVecShapeIndex[iq].second;
476 //
477 // TPZFNMatrix<3> ivec(3,1);
478 // ivec(0,0) = datavec[0].fNormalVec(0,ivecind);
479 // ivec(1,0) = datavec[0].fNormalVec(1,ivecind);
480 // ivec(2,0) = datavec[0].fNormalVec(2,ivecind);
481 // TPZFNMatrix<3> axesvec(3,1);
482 // datavec[0].axes.Multiply(ivec,axesvec);
483 //
484 // REAL divwq = 0.;
485 // for(int iloc=0; iloc<fDim; iloc++)
486 // {
487 // divwq += axesvec(iloc,0)*dphiQ(iloc,ishapeind);
488 // }
489 // for (int jp=0; jp<phrp; jp++) {
490 //
491 // REAL fact = (-1.)*weight*phip(jp,0)*divwq;
492 // // Matrix B
493 // ek(iq, phrq+jp) += fact;
494 //
495 // // Matrix B^T
496 // ek(phrq+jp,iq) += fact;
497 //
498 //
499 // //Inserindo termo de estabilizacao: delta1
500 // if(fIsStabilized==true)
501 // {
502 // REAL dotVGradP = 0.;
503 // for(int k =0; k<fDim; k++)
504 // {
505 // dotVGradP += ivec(k,0)*phiQ(ishapeind,0)*dphiP(k,jp);
506 // }
507 //
508 // REAL integration = (-1.)*weight*fh2*fdelta1*dotVGradP;
509 //
510 // // Estabilizacao delta1 na Matrix B
511 // ek(iq, phrq+jp) += integration;
512 //
513 // // Estabilizacao delta1 na Matrix BˆT
514 // ek(phrq+jp,iq) += integration;
515 // }
516 // }
517 // }
518 //
519 // //termo fonte referente a equacao da pressao
520 // for(int ip=0; ip<phrp; ip++){
521 // ef(phrq+ip,0) += (-1.)*weight*force*phip(ip,0);
522 // }
523 //
524 // //Contribution for estabilization delta1 for gradPu*gradPv. Matrix D
525 // if(fIsStabilized==true)
526 // {
527 // REAL mK = fk/fvisc;
528 // for(int ip=0; ip<phrp; ip++)
529 // {
530 // for(int jp=0; jp<phrp; jp++)
531 // {
532 // for(int k =0; k<fDim; k++)
533 // {
534 // ek(phrq+ip, phrq+jp) += (-1.)*weight*fh2*fdelta1*mK*dphiP(k,ip)*dphiP(k,jp);
535 // }
536 //
537 // }
538 // }
539 // }
551 //
552 //}
553 
554 void TPZMixedPoisson::ContributeBC(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc){
555 
556 #ifdef PZDEBUG
557  int nref = datavec.size();
558 // if (nref != 2 ) {
559 // std::cout << " Erro.!! datavec tem que ser de tamanho 2 \n";
560 // DebugStop();
561 // }
562  if (bc.Type() > 2 ) {
563  std::cout << " Erro.!! Neste material utiliza-se apenas condicoes de Neumann e Dirichlet\n";
564  DebugStop();
565  }
566 #endif
567 
568 
569  int dim = Dimension();
570 
571  TPZFMatrix<REAL> &phiQ = datavec[0].phi;
572  int phrq = phiQ.Rows();
573 
574  REAL v2 = bc.Val2()(0,0);
575  REAL v1 = bc.Val1()(0,0);
576  if(bc.HasForcingFunction())
577  {
578  TPZManVector<STATE> res(3);
579  TPZFNMatrix<9,STATE> gradu(3,1,0.);
580  bc.ForcingFunction()->Execute(datavec[0].x,res,gradu);
581  if(bc.Type() == 0)
582  {
583  v2 = res[0];
584  }
585  else if(bc.Type() == 1 || bc.Type() == 2)
586  {
587  TPZFNMatrix<9,REAL> PermTensor, InvPermTensor;
588  GetPermeability(datavec[0].x, PermTensor, InvPermTensor);
589  REAL normflux = 0.;
590  for(int i=0; i<3; i++)
591  {
592  for(int j=0; j<3; j++)
593  {
594  normflux -= datavec[0].normal[i]*PermTensor(i,j)*gradu(j,0);
595  }
596  }
597  v2 = -normflux;
598  if(bc.Type() ==2)
599  {
600  v2 = -res[0]+v2/v1;
601  }
602  }
603  else
604  {
605  DebugStop();
606  }
607  }else
608  {
609  v2 = bc.Val2()(0,0);
610  }
611 
612  switch (bc.Type()) {
613  case 0 : // Dirichlet condition
614  //primeira equacao
615  for(int iq=0; iq<phrq; iq++)
616  {
617  //the contribution of the Dirichlet boundary condition appears in the flow equation
618  ef(iq,0) += (-1.)*v2*phiQ(iq,0)*weight;
619  }
620  break;
621 
622  case 1 : // Neumann condition
623  //primeira equacao
624  for(int iq=0; iq<phrq; iq++)
625  {
626  ef(iq,0)+= gBigNumber*v2*phiQ(iq,0)*weight;
627  for (int jq=0; jq<phrq; jq++) {
628 
629  ek(iq,jq)+= gBigNumber*phiQ(iq,0)*phiQ(jq,0)*weight;
630  }
631  }
632  break;
633 
634  case 2 : // mixed condition
635  for(int iq = 0; iq < phrq; iq++) {
636 
637  ef(iq,0) += v2*phiQ(iq,0)*weight;
638  for (int jq = 0; jq < phrq; jq++) {
639  ek(iq,jq) += weight/v1*phiQ(iq,0)*phiQ(jq,0);
640  }
641  }
642  break;
643  }
644 
645 }
646 
648 int TPZMixedPoisson::VariableIndex(const std::string &name){
649  if(!strcmp("Flux",name.c_str())) return 31;
650  if(!strcmp("Pressure",name.c_str())) return 32;
651  if(!strcmp("GradFluxX",name.c_str())) return 33;
652  if(!strcmp("GradFluxY",name.c_str())) return 34;
653  if(!strcmp("DivFlux",name.c_str())) return 35;
654 
655  if(!strcmp("ExactPressure",name.c_str())) return 36;
656  if(!strcmp("ExactFlux",name.c_str())) return 37;
657 
658  if(!strcmp("POrder",name.c_str())) return 38;
659  if(!strcmp("GradPressure",name.c_str())) return 39;
660  if(!strcmp("Divergence",name.c_str())) return 40;
661  if(!strcmp("ExactDiv",name.c_str())) return 41;
662  if (!strcmp("Derivative",name.c_str())) {
663  return 42;
664  }
665  if (!strcmp("Permeability",name.c_str())) {
666  return 43;
667  }
668  if(!strcmp("g_average",name.c_str())) return 44;
669  if(!strcmp("u_average",name.c_str())) return 45;
670  return TPZMatPoisson3d::VariableIndex(name);
671 
672 }
673 
674 int TPZMixedPoisson::NSolutionVariables(int var){
675  if(var == 31) return 3;
676  if(var == 32 || var==8) return 1;
677  if(var == 33) return 3;
678  if(var == 34) return 3;
679  if(var == 35) return 1;
680  if(var == 36) return 1;
681  if(var == 37) return fDim;
682  if(var == 38) return 1;
683  if(var == 39) return fDim;
684  if(var == 40 || var == 41) return 1;
685  if(var == 42) return 3;
686  if(var == 43) return 1;
687  if(var == 44) return 1;
688  if(var == 45) return 1;
689  return TPZMatPoisson3d::NSolutionVariables(var);
690 }
691 
692 void TPZMixedPoisson::Solution(TPZVec<TPZMaterialData> &datavec, int var, TPZVec<STATE> &Solout){
693 
694  Solout.Resize( this->NSolutionVariables(var));
695 
696  TPZVec<STATE> SolP, SolQ;
697 
698  TPZFNMatrix<9,REAL> PermTensor = fTensorK;
699  TPZFNMatrix<9,REAL> InvPermTensor = fInvK;
700  //int rtens = 2*fDim;
701  if(fPermeabilityFunction){
702  PermTensor.Redim(3,3);
703  InvPermTensor.Redim(3,3);
704  TPZFNMatrix<18,STATE> resultMat(2*3,3);
705  TPZManVector<STATE> res;
706  fPermeabilityFunction->Execute(datavec[1].x,res,resultMat);
707 
708  for(int id=0; id<3; id++){
709  for(int jd=0; jd<3; jd++){
710 
711  PermTensor(id,jd) = resultMat(id,jd);
712  InvPermTensor(id,jd) = resultMat(id+3,jd);
713  }
714  }
715  }
716 
717 
718 
719  // SolQ = datavec[0].sol[0];
720  SolP = datavec[1].sol[0];
721 
722  if(var == 31){ //function (state variable Q)
723  for (int i=0; i<3; i++)
724  {
725  Solout[i] = datavec[0].sol[0][i];
726 
727  }
728  return;
729  }
730 
731  if(var == 32){
732  Solout[0] = SolP[0];//function (state variable p)
733  return;
734  }
735 
736  if(var==33){
737  Solout[0]=datavec[0].dsol[0](0,0);
738  Solout[1]=datavec[0].dsol[0](1,0);
739  Solout[2]=datavec[0].dsol[0](2,0);
740  return;
741  }
742 
743  if(var==34){
744  Solout[0]=datavec[0].dsol[0](0,1);
745  Solout[1]=datavec[0].dsol[0](1,1);
746  Solout[2]=datavec[0].dsol[0](2,1);
747  return;
748  }
749 
750  if(var==35){
751  Solout[0]=datavec[0].dsol[0](0,0)+datavec[0].dsol[0](1,1);
752  return;
753  }
754 
755  TPZVec<REAL> ptx(3);
756  TPZVec<STATE> solExata(1);
757  TPZFNMatrix<3,STATE> flux(fDim+1,1);//pq colocar fdim +1?
758  TPZFNMatrix<3,STATE> gradu(fDim+1,1);
759 
760  //Exact solution
761  if(var == 36){
762  if(fForcingFunctionExact)
763  {
764  fForcingFunctionExact->Execute(datavec[0].x, solExata,flux);
765  }
766  Solout[0] = solExata[0];
767  return;
768  }//var6
769 
770  if(var == 37){
771 
772 
773  if(fForcingFunctionExact)
774  {
775  fForcingFunctionExact->Execute(datavec[0].x, solExata,gradu);
776 
777  }
778 
779  TPZFNMatrix<3, REAL> fluxtmp(fDim + 1, 1);
780  TPZFNMatrix<3, REAL> gradutmp(fDim + 1, 1);
781  for (int idi = 0; idi < gradu.Rows(); idi++)
782  for (int jdi = 0; jdi < gradu.Cols(); jdi++)
783  gradutmp(idi, jdi) = gradu(idi, jdi);
784 
785  PermTensor.Multiply(gradutmp, fluxtmp);
786 
787  for (int i=0; i<fDim; i++)
788  {
789  Solout[i] = fluxtmp(i,0);
790  }
791 
792  return;
793  }//var7
794 
795  if(var==38){
796  Solout[0] = datavec[1].p;
797  return;
798  }
799 
800  if(var==39){
801  TPZFNMatrix<3,REAL> dsoldx;
802  TPZFMatrix<REAL> dsoldaxes(fDim,1);
803 // dsoldaxes(0,0) = datavec[1].dsol[0][0];
804 // dsoldaxes(1,0) = datavec[1].dsol[0][1];
805 
806  for (int i=0; i<fDim; i++)
807  {
808  dsoldaxes(i,0) = datavec[1].dsol[0][i];;
809  }
810 
811  TPZAxesTools<REAL>::Axes2XYZ(dsoldaxes, dsoldx, datavec[1].axes);
812 // Solout[0] = dsoldx(0,0);
813 // Solout[1] = dsoldx(1,0);
814 
815  for (int i=0; i<fDim; i++)
816  {
817  Solout[i] = dsoldx(i,0);
818  }
819 
820 
821  return;
822  }
823 
824  if(var==40){
825  Solout[0]=datavec[0].dsol[0](0,0)+datavec[0].dsol[0](1,1);
826  return;
827  }
828 
829  if(var==41){
830  fForcingFunctionExact->Execute(datavec[0].x,solExata,flux);
831  Solout[0]=flux(2,0);
832  return;
833  }
834  if(var == 42)
835  {
836  for(int i=0; i<fDim; i++)
837  {
838  Solout[i] = 0.;
839  }
840  for (int i=0; i<fDim; i++) {
841  for (int j=0; j<fDim; j++) {
842  Solout[i] -= InvPermTensor(i,j)*datavec[0].sol[0][i];
843  }
844  }
845  return;
846  }
847  if(var ==43)
848  {
849  Solout[0] = PermTensor(0,0);
850  return;
851  }
852 
853  if(datavec.size() == 4)
854  {
855  if(var ==44)
856  {
857  Solout[0] = datavec[2].sol[0][0];
858  return;
859  }
860  if(var ==45)
861  {
862  Solout[0] = datavec[3].sol[0][0];
863  return;
864  }
865 
866  }
867 
868  TPZMatPoisson3d::Solution(datavec,var,Solout);
869 }
870 
871 
872 void TPZMixedPoisson::FillDataRequirements(TPZVec<TPZMaterialData > &datavec)
873 {
874  int nref = datavec.size();
875  for(int i = 0; i<nref; i++ )
876  {
877  datavec[i].SetAllRequirements(false);
878  datavec[i].fNeedsNeighborSol = false;
879  datavec[i].fNeedsNeighborCenter = false;
880  datavec[i].fNeedsNormal = true;
881  datavec[i].fNeedsHSize = true;
882  }
883 }
884 
885 
886 void TPZMixedPoisson::Errors(TPZVec<TPZMaterialData> &data, TPZVec<STATE> &u_exact, TPZFMatrix<STATE> &du_exact, TPZVec<REAL> &errors)
887 {
888 
889 
890 // TPZVec<REAL> &x,TPZVec<STATE> &u,
891 // TPZFMatrix<STATE> &dudx, TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
892 // TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) {
893 
894  errors.Resize(NEvalErrors());
895  errors.Fill(0.0);
896  TPZManVector<STATE,3> deriv(3,0.), pressure(1,0.);
897  this->Solution(data,VariableIndex("Derivative"), deriv);
898  this->Solution(data,VariableIndex("Pressure"), pressure);
899 
900  TPZFNMatrix<9,STATE> perm(2*3,3);
901  TPZManVector<STATE,3> val(3);
902  if (fPermeabilityFunction) {
903  fPermeabilityFunction->Execute(data[0].x, val, perm);
904  }
905  else
906  {
907  for (int i=0; i<3; i++) {
908  for (int j=0; j<3; j++)
909  {
910  perm(i,j) = this->fTensorK(i,j);
911  perm(3+i,j) = this->fInvK(i,j);
912  }
913  }
914  }
915 
916 #ifdef LOG4CXX
917  if(logerror->isDebugEnabled())
918  {
919  std::stringstream sout;
920  sout.precision(14);
921  sout << "x " << data[0].x << std::endl;
922  sout << "u_exact " << u_exact << std::endl;
923  du_exact.Print("du_exact",sout);
924  sout << "deriv " << deriv << std::endl;
925  sout << "pressure " << pressure << std::endl;
926  LOGPZ_DEBUG(logerror, sout.str())
927  }
928 #endif
929  int id;
930  //values[1] : eror em norma L2
931  STATE diff = pressure[0] - u_exact[0];
932  errors[1] = diff*diff;
933  //values[2] : erro em semi norma H1
934  errors[2] = 0.;
935  for(id=0; id<3; id++) {
936  diff = deriv[id] - du_exact(id,0);
937  errors[2] += fK*diff*diff;
938  }
939  //values[0] : erro em norma H1 <=> norma Energia
940 
941  errors[0] = 0.;
942  for (int i=0; i<3; i++) {
943  for (int j=0; j<3; j++) {
944  errors[0] += (deriv[i] - du_exact(i,0))*perm(i,j)*(deriv[i] - du_exact(j,0));
945  }
946  }
947 }
948 
949 int TPZMixedPoisson::ClassId() const{
950  return Hash("TPZMixedPoisson") ^ TPZMatPoisson3d::ClassId() << 1;
951 }
952 
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZMixedPoisson::NEvalErrors
virtual int NEvalErrors() override
Returns the number of norm errors. Default is 3: energy, L2 and H1.
Definition: mixedpoisson.h:199
TPZMixedPoisson::fInvK
TPZFNMatrix< 9, REAL > fInvK
inverse of the permeability tensor.
Definition: mixedpoisson.h:44
TPZMixedPoisson::fvisc
REAL fvisc
fluid viscosity
Definition: mixedpoisson.h:47
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.
TPZManVector
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZAxesTools::Axes2XYZ
static void Axes2XYZ(const TPZFMatrix< TVar > &dudaxes, TPZFMatrix< TVar > &dudx, const TPZFMatrix< REAL > &axesv, bool colMajor=true)
Makes the basis transformation from axes basis to euclidian basis.
Definition: pzaxestools.h:61
TPZMatPoisson3d::Dimension
int Dimension() const override
Returns the integrable dimension of the material.
Definition: pzpoisson3d.h:158
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZMatrix::Identity
virtual void Identity()
Converts the matrix in an identity matrix.
Definition: pzmatrix.cpp:199
TPZMixedPoisson::GetPermeability
void GetPermeability(TPZVec< REAL > &x, TPZFMatrix< REAL > &K, TPZFMatrix< REAL > &invK)
return the permeability and compute it if there is permeability function
Definition: mixedpoisson.cpp:102
TPZMixedPoisson::ContributeBC
virtual void ContributeBC(TPZVec< TPZMaterialData > &datavec, REAL weight, TPZFMatrix< STATE > &ek, 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: mixedpoisson.cpp:554
TPZMixedPoisson::fh2
REAL fh2
Coeficient that multiplies the Stabilization term fdelta1.
Definition: mixedpoisson.h:57
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZMixedPoisson::Name
virtual std::string Name() override
Returns the name of the material.
Definition: mixedpoisson.h:83
TPZBndCond::Val2
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZMixedPoisson::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: mixedpoisson.cpp:648
TPZMixedPoisson::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: mixedpoisson.cpp:872
TPZMatPoisson3d::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzpoisson3d.cpp:1426
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZMaterial::ForcingFunction
TPZAutoPointer< TPZFunction< STATE > > & ForcingFunction()
Returns a procedure as source function for the material.
Definition: TPZMaterial.h:397
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
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
TPZMixedPoisson::Solution
virtual void Solution(TPZVec< TPZMaterialData > &datavec, int var, TPZVec< STATE > &Solout) override
It return a solution to multiphysics simulation.
Definition: mixedpoisson.cpp:692
TPZMatPoisson3d::fK
STATE fK
Coeficient which multiplies the Laplacian operator.
Definition: pzpoisson3d.h:37
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
TPZMaterial::fForcingFunctionExact
TPZAutoPointer< TPZFunction< STATE > > fForcingFunctionExact
Pointer to exact solution function, needed to calculate exact error.
Definition: TPZMaterial.h:50
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
test.res
string res
Definition: test.py:151
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZMixedPoisson::ff
REAL ff
Forcing function value.
Definition: mixedpoisson.h:38
TPZMatPoisson3d::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: pzpoisson3d.cpp:489
TPZMaterial::HasForcingFunction
virtual int HasForcingFunction()
Directive that gives true if the material has a forcing function.
Definition: TPZMaterial.h:472
TPZMatPoisson3d::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: pzpoisson3d.h:264
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZMaterial::gBigNumber
static REAL gBigNumber
Big number to penalization method, used for Dirichlet conditions.
Definition: TPZMaterial.h:83
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZMixedPoisson
Material to solve a mixed poisson problem 2d by multiphysics simulation.
Definition: mixedpoisson.h:34
TPZMixedPoisson::Contribute
virtual void Contribute(TPZVec< TPZMaterialData > &datavec, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
It computes a contribution to the stiffness matrix and load vector at one integration point to multip...
Definition: mixedpoisson.cpp:127
TPZMatPoisson3d::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: pzpoisson3d.cpp:520
TPZMatPoisson3d
DESCRIBE PLEASE.
Definition: pzpoisson3d.h:26
TPZMixedPoisson::fTensorK
TPZFNMatrix< 9, REAL > fTensorK
permeability tensor. Coeficient which multiplies the gradient operator
Definition: mixedpoisson.h:41
TPZMixedPoisson::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: mixedpoisson.cpp:90
TPZMixedPoisson::fPermeabilityFunction
TPZAutoPointer< TPZFunction< STATE > > fPermeabilityFunction
post-processing procedure for error estimation as Ainsworth
Definition: mixedpoisson.h:63
TPZMatrix::Multiply
virtual void Multiply(const TPZFMatrix< TVar > &A, TPZFMatrix< TVar > &res, int opt=0) const
It mutiplies itself by TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:1916
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZMatPoisson3d::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: pzpoisson3d.cpp:99
TPZVec::Fill
void Fill(const T &copy, const int64_t from=0, const int64_t numelem=-1)
Will fill the elements of the vector with a copy object.
Definition: pzvec.h:460
pzaxestools.h
Contains declaration of the TPZAxesTools class which implements verifications over axes...
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZMixedPoisson::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: mixedpoisson.cpp:674
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZMatPoisson3d::fDim
int fDim
Problem dimension.
Definition: pzpoisson3d.h:34
TPZMatPoisson3d::operator=
TPZMatPoisson3d & operator=(const TPZMatPoisson3d &copy)
Definition: pzpoisson3d.cpp:63
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZMixedPoisson::~TPZMixedPoisson
virtual ~TPZMixedPoisson()
Definition: mixedpoisson.cpp:57
TPZMixedPoisson::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: mixedpoisson.cpp:94
TPZMixedPoisson::fIsStabilized
bool fIsStabilized
Choose Stabilized method.
Definition: mixedpoisson.h:50
TPZMixedPoisson::fdelta1
REAL fdelta1
Coeficient of Stabilization.
Definition: mixedpoisson.h:53
TPZMixedPoisson::Errors
virtual void Errors(TPZVec< TPZMaterialData > &data, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &errors) override
Definition: mixedpoisson.cpp:886
TPZMixedPoisson::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: mixedpoisson.cpp:949
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
TPZMixedPoisson::operator=
TPZMixedPoisson & operator=(const TPZMixedPoisson &copy)
Definition: mixedpoisson.cpp:74
Generated on Sun Aug 16 2020 11:17:01 for NeoPZ by   doxygen 1.8.13