NeoPZ
pzpoisson3d.cpp
Go to the documentation of this file.
1 
6 #include "pzpoisson3d.h"
7 #include "pzelmat.h"
8 #include "pzbndcond.h"
9 #include "pzmatrix.h"
10 #include "pzfmatrix.h"
11 #include "pzerror.h"
12 #include "pzmaterialdata.h"
13 #include <math.h>
14 #include "pzlog.h"
15 
16 #include <cmath>
17 
18 #ifdef LOG4CXX
19 static LoggerPtr logger(Logger::getLogger("pz.material.poisson3d"));
20 #endif
21 
22 
23 using namespace std;
24 STATE TPZMatPoisson3d::gAlfa = 0.5;
25 
26 TPZMatPoisson3d::TPZMatPoisson3d(int nummat, int dim) : TPZRegisterClassId(&TPZMatPoisson3d::ClassId),
27 TPZDiscontinuousGalerkin(nummat), fXf(0.), fDim(dim), fSD(0.) {
28  if(dim < 1)
29  {
30  DebugStop();
31  }
32  fK = 1.;
33  fC = 0.;
34  fConvDir[0] = 0.;
35  fConvDir[1] = 0.;
36  fConvDir[2] = 0.;
37  fPenaltyConstant = 1000.;
38  this->SetNonSymmetric();
39  this->SetNoPenalty();
40  fShapeHdiv=false;
41  fNeumann=false;
42 }
43 
44 TPZMatPoisson3d::TPZMatPoisson3d():TPZRegisterClassId(&TPZMatPoisson3d::ClassId),
45 TPZDiscontinuousGalerkin(), fXf(0.), fDim(1), fSD(0.){
46  fK = 1.;
47  fC = 0.;
48  fConvDir[0] = 0.;
49  fConvDir[1] = 0.;
50  fConvDir[2] = 0.;
51  fPenaltyConstant = 1000.;
52  this->SetNonSymmetric();
53  this->SetNoPenalty();
54  fShapeHdiv=false;
55  fNeumann=false;
56 }
57 
58 TPZMatPoisson3d::TPZMatPoisson3d(const TPZMatPoisson3d &copy):TPZRegisterClassId(&TPZMatPoisson3d::ClassId),
59 TPZDiscontinuousGalerkin(copy){
60  this->operator =(copy);
61 }
62 
63 TPZMatPoisson3d & TPZMatPoisson3d::operator=(const TPZMatPoisson3d &copy){
64  TPZDiscontinuousGalerkin::operator = (copy);
65  fXf = copy.fXf;
66  fDim = copy.fDim;
67  fK = copy.fK;
68  fC = copy.fC;
69  fShapeHdiv= copy.fShapeHdiv;
70  for (int i = 0; i < 3; i++) fConvDir[i] = copy.fConvDir[i];
71  fSymmetry = copy.fSymmetry;
72  fSD = copy.fSD;
73  fPenaltyConstant = copy.fPenaltyConstant;
74  this->fPenaltyType = copy.fPenaltyType;
75  return *this;
76 }
77 
78 void TPZMatPoisson3d::SetParameters(STATE diff, REAL conv, TPZVec<REAL> &convdir) {
79  fK = diff;
80  fC = conv;
81  int d;
82  for(d=0; d<fDim; d++) fConvDir[d] = convdir[d];
83 }
84 
85 void TPZMatPoisson3d::GetParameters(STATE &diff, REAL &conv, TPZVec<REAL> &convdir) {
86  diff = fK;
87  conv = fC;
88  int d;
89  for(d=0; d<fDim; d++) convdir[d] = fConvDir[d];
90 }
91 
92 TPZMatPoisson3d::~TPZMatPoisson3d() {
93 }
94 
95 //int TPZMatPoisson3d::NStateVariables() const {
96 // return 1;
97 //}
98 
99 void TPZMatPoisson3d::Print(std::ostream &out) {
100  out << "name of material : " << Name() << "\n";
101  out << "Laplace operator multiplier fK "<< fK << endl;
102  out << "Convection coeficient fC " << fC << endl;
103  out << "Convection direction " << fConvDir[0] << ' ' << fConvDir[1] << ' ' << fConvDir[2] << endl;
104  out << "Forcing vector fXf " << fXf << endl;
105  out << "Penalty constant " << fPenaltyConstant << endl;
106  out << "Base Class properties :";
107  TPZMaterial::Print(out);
108  out << "\n";
109 }
110 
111 void TPZMatPoisson3d::Contribute(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef) {
112 
113  if(data.numberdualfunctions)
114  {
115  ContributeHDiv(data , weight , ek, ef);
116 
117  return;
118  }
119 
120  if(fNeumann){
121 
122  LocalNeumanContribute(data , weight , ek, ef);
123 
124  return;
125  }
126 
127  TPZFMatrix<REAL> &phi = data.phi;
128  TPZFMatrix<REAL> &dphi = data.dphix;
129  TPZVec<REAL> &x = data.x;
130  TPZFMatrix<REAL> &axes = data.axes;
131  TPZFMatrix<REAL> &jacinv = data.jacinv;
132  int phr = phi.Rows();
133 
134  STATE fXfLoc = fXf;
135 
136  if(fForcingFunction) { // phi(in, 0) = phi_in
137  TPZManVector<STATE,1> res(1);
138  TPZFMatrix<STATE> dres(Dimension(),1);
139  fForcingFunction->Execute(x,res,dres); // dphi(i,j) = dphi_j/dxi
140  fXfLoc = res[0];
141  }
142  REAL delx = 0.;
143  STATE ConvDirAx[3] = {0.};
144  if(fC != 0.0) {
145  int di,dj;
146  delx = 0.;
147  for(di=0; di<fDim; di++) {
148  for(dj=0; dj<fDim; dj++) {
149  delx = (delx<fabs(jacinv(di,dj))) ? fabs(jacinv(di,dj)) : delx;
150  }
151  }
152  delx = 2./delx;
153 
154 
155  switch(fDim) {
156  case 1:
157  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
158  break;
159  case 2:
160  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
161  ConvDirAx[1] = axes(1,0)*fConvDir[0]+axes(1,1)*fConvDir[1]+axes(1,2)*fConvDir[2];
162  break;
163  case 3:
164  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
165  ConvDirAx[1] = axes(1,0)*fConvDir[0]+axes(1,1)*fConvDir[1]+axes(1,2)*fConvDir[2];
166  ConvDirAx[2] = axes(2,0)*fConvDir[0]+axes(2,1)*fConvDir[1]+axes(2,2)*fConvDir[2];
167  break;
168  default:
169  PZError << "TPZMatPoisson3d::Contribute dimension error " << fDim << endl;
170  }
171  }
172 
173 // std::cout << " fxfloc " << fXfLoc << " weight " << weight << " prod " << fXfLoc*weight << std::endl;
174  //Equacao de Poisson
175  for( int in = 0; in < phr; in++ ) {
176  int kd;
177  STATE dphiic = 0;
178  for(kd = 0; kd<fDim; kd++) dphiic += ConvDirAx[kd]*(STATE)dphi(kd,in);
179  ef(in, 0) += - (STATE)weight * fXfLoc * ( (STATE)phi(in,0) + (STATE)(0.5*delx*fC)*fSD*dphiic );
180  for( int jn = 0; jn < phr; jn++ ) {
181  for(kd=0; kd<fDim; kd++) {
182  ek(in,jn) += (STATE)weight * (
183  +fK * (STATE)( dphi(kd,in) * dphi(kd,jn) )
184  - (STATE)(fC* dphi(kd,in) * phi(jn)) * ConvDirAx[kd]
185  + (STATE)(0.5 * delx * fC * dphi(kd,jn)) * fSD * dphiic * ConvDirAx[kd]
186  );
187  }
188  }
189  }
190 
191  if (this->IsSymetric()){
192  if ( !ek.VerifySymmetry() ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
193  }
194 
195 
196 }
197 
199 void TPZMatPoisson3d::ContributeHDiv(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
200 {
207  //TPZVec<REAL> &x = data.x;
208  STATE fXfLoc = fXf;
209  if(fForcingFunction) { // phi(in, 0) = phi_in
210  TPZManVector<STATE> res(1);
211  fForcingFunction->Execute(data.x,res); // dphi(i,j) = dphi_j/dxi
212  fXfLoc = res[0];
213  }
214  int numvec = data.fVecShapeIndex.NElements();
215  int numdual = data.numberdualfunctions;
216  int numprimalshape = data.phi.Rows()-numdual;
217 
218  int i,j;
219  REAL kreal = 0.;
220 #ifdef STATE_COMPLEX
221  kreal = fK.real();
222 #else
223  kreal = fK;
224 #endif
225  REAL ratiok = 1./kreal;
226  for(i=0; i<numvec; i++)
227  {
228  int ivecind = data.fVecShapeIndex[i].first;
229  int ishapeind = data.fVecShapeIndex[i].second;
230  for (j=0; j<numvec; j++) {
231  int jvecind = data.fVecShapeIndex[j].first;
232  int jshapeind = data.fVecShapeIndex[j].second;
233  REAL prod = data.fNormalVec(0,ivecind)*data.fNormalVec(0,jvecind)+
234  data.fNormalVec(1,ivecind)*data.fNormalVec(1,jvecind)+
235  data.fNormalVec(2,ivecind)*data.fNormalVec(2,jvecind);//faz o produto escalar entre u e v--> Matriz A
236  ek(i,j) += weight*ratiok*data.phi(ishapeind,0)*data.phi(jshapeind,0)*prod;
237 
238 
239 
240  }
241  TPZFNMatrix<3> ivec(3,1);
242  ivec(0,0) = data.fNormalVec(0,ivecind);
243  ivec(1,0) = data.fNormalVec(1,ivecind);
244  ivec(2,0) = data.fNormalVec(2,ivecind);
245  TPZFNMatrix<3> axesvec(3,1);
246  data.axes.Multiply(ivec,axesvec);
247  int iloc;
248  REAL divwq = 0.;
249  for(iloc=0; iloc<fDim; iloc++)
250  {
251  divwq += axesvec(iloc,0)*data.dphix(iloc,ishapeind);
252  }
253  for (j=0; j<numdual; j++) {
254  REAL fact = (-1.)*weight*data.phi(numprimalshape+j,0)*divwq;//calcula o termo da matriz B^T e B
255  ek(i,numvec+j) += fact;
256  ek(numvec+j,i) += fact;//-div
257  }
258  }
259  for(i=0; i<numdual; i++)
260  {
261  ef(numvec+i,0) += (STATE)((-1.)*weight*data.phi(numprimalshape+i,0))*fXfLoc;//calcula o termo da matriz f
262 
263 #ifdef LOG4CXX
264  if (logger->isDebugEnabled())
265  {
266  std::stringstream sout;
267  sout<< "Verificando termo fonte\n";
268  sout << " pto " <<data.x << std::endl;
269  sout<< " fpto " <<fXfLoc <<std::endl;
270  LOGPZ_DEBUG(logger,sout.str())
271  }
272 #endif
273  }
274 
275 }
276 
277 //
278 void TPZMatPoisson3d::LocalNeumanContribute(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
279 {
285  int newRows=ek.Rows()+1;
286  int newCols=ek.Cols()+1;
287  ek.Resize( newRows, newCols );
288  ef.Resize(newRows,1);
289 
290  STATE fXfLoc = fXf;
291  if(fForcingFunction) { // phi(in, 0) = phi_in
292  TPZManVector<STATE> res(1);
293  fForcingFunction->Execute(data.x,res); // dphi(i,j) = dphi_j/dxi
294  fXfLoc = res[0];
295  }
296 
297 
298  TPZFMatrix<REAL> &phi = data.phi;
299  TPZFMatrix<REAL> &dphi = data.dphix;
300  TPZVec<REAL> &x = data.x;
301  TPZFMatrix<REAL> &axes = data.axes;
302  TPZFMatrix<REAL> &jacinv = data.jacinv;
303  TPZVec<STATE> &sol=data.sol[0];
304  TPZGradSolVec &dsol=data.dsol;
305 
306 
307  int phr = phi.Rows();
308  int nlinhaek=ek.Rows();
309  int ncolek=ek.Cols();
310 
311 
312  int i,j;
313  REAL kreal = 0.;
314 #ifdef STATE_COMPLEX
315  kreal = fK.real();
316 #else
317  kreal = fK;
318 #endif
319  REAL ratiok = 1./kreal;
320  //--
321  for( int in = 0; in < phr; in++ ){
322  for( int jn = 0; jn < phr; jn++ ) {
323  for(int kd=0; kd<fDim; kd++) {
324  ek(in,jn) += (STATE)weight * (fK * (STATE)( dphi(kd,in) * dphi(kd,jn) ));//gradphi_i.gradphi_j
325  }
326 
327  }
328  ek(in,ncolek-1) += (STATE)weight * phi(in);
329  ek(nlinhaek-1,in) += (STATE)weight * phi(in);
330 
331  }
332 
333  for( int in = 0; in < phr; in++ ) {
334  for( int jn = 0; jn < phr; jn++ ) {
335  for(int kd=0; kd<fDim; kd++) {
336  ef(in,0) += (STATE)weight * ( dphi(kd,in) * dsol[0](kd,jn) );
337  }
338 
339  }
340  }
341  ef(phr,0)+=(STATE)weight * (sol[0]);
342 
343 }
344 
345 
347 
348 void TPZMatPoisson3d::ContributeBCHDiv(TPZMaterialData &data,REAL weight,
349  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
350  int numvec = data.phi.Rows();
351 
352  TPZFMatrix<REAL> &phi = data.phi;
353  TPZManVector<STATE,1> v2(1);
354  v2[0] = bc.Val2()(0,0);
355  if (bc.HasForcingFunction()) {
356  bc.ForcingFunction()->Execute(data.x, v2);
357  }
358 
359  switch (bc.Type()) {
360  case 1 : // Neumann condition
361  int i,j;
362  for(i=0; i<numvec; i++)
363  {
364  //int ishapeind = data.fVecShapeIndex[i].second;
365  ef(i,0)+= (STATE)(gBigNumber * phi(i,0) * weight) * v2[0];
366  for (j=0; j<numvec; j++) {
367  //int jshapeind = data.fVecShapeIndex[j].second;
368  ek(i,j) += gBigNumber * phi(i,0) * phi(j,0) * weight;
369  }
370  }
371  break;
372  case 0 :
373  {// Dirichlet condition
374  int in;
375  for(in = 0 ; in < numvec; in++) {
376  ef(in,0) += (STATE)((-1.)* phi(in,0) * weight)*v2[0];
377  }
378  }
379  break;
380  case 2 : // mixed condition
381  {
382  int in,jn;
383  for(in = 0 ; in < numvec; in++) {
384  //int ishapeind = data.fVecShapeIndex[in].second;
385  ef(in,0) += v2[0] * (STATE)(phi(in,0) * weight);
386  for (jn = 0; jn < numvec; jn++) {
387  // int jshapeind = data.fVecShapeIndex[jn].second;
388  ek(in,jn) += (STATE)(weight*phi(in,0)*phi(jn,0)) *bc.Val1()(0,0);
389  }
390  }
391  }
392  break;
393  case 3: // outflow condition
394  break;
395  }
396 
397  if (this->IsSymetric()) {//only 1.e-3 because of bignumbers.
398  if ( !ek.VerifySymmetry( 1.e-3 ) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
399  }
400 
401 
402 }
403 void TPZMatPoisson3d::ContributeBC(TPZMaterialData &data,REAL weight,
404  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
405 
406  if(fShapeHdiv==true)//data.fShapeType == TPZMaterialData::EVecandShape || data.fShapeType == TPZMaterialData::EVecShape)
407  {
408 
409  ContributeBCHDiv(data , weight , ek, ef, bc);
410 
411 
412  return;
413  }
414 
415  TPZFMatrix<REAL> &phi = data.phi;
416  TPZFMatrix<REAL> &axes = data.axes;
417  int phr = phi.Rows();
418  short in,jn;
419  STATE v2[1];
420  v2[0] = bc.Val2()(0,0);
421 
422  if(bc.HasForcingFunction()) { // phi(in, 0) = phi_in // JORGE 2013 01 26
423  TPZManVector<STATE,1> res(1);
424  bc.ForcingFunction()->Execute(data.x,res); // dphi(i,j) = dphi_j/dxi
425 // if(fabs(res[0]) > 1.e-6)
426 // {
427 // bc.ForcingFunction()->Execute(data.x,res); // dphi(i,j) = dphi_j/dxi
428 // DebugStop();
429 // }
430  v2[0] = res[0];
431  }
432 
433  switch (bc.Type()) {
434  case 0 : // Dirichlet condition
435  for(in = 0 ; in < phr; in++) {
436  ef(in,0) += (STATE)(gBigNumber* phi(in,0) * weight) * v2[0];
437  for (jn = 0 ; jn < phr; jn++) {
438  ek(in,jn) += gBigNumber * phi(in,0) * phi(jn,0) * weight;
439  }
440  }
441  break;
442  case 1 : // Neumann condition
443  for(in = 0 ; in < phi.Rows(); in++) {
444  ef(in,0) += v2[0] * (STATE)(phi(in,0) * weight);
445  }
446  break;
447  case 2 : // mixed condition
448  for(in = 0 ; in < phi.Rows(); in++) {
449  ef(in, 0) += v2[0] * (STATE)(phi(in, 0) * weight);
450  for (jn = 0 ; jn < phi.Rows(); jn++) {
451  ek(in,jn) += bc.Val1()(0,0) * (STATE)(phi(in,0) * phi(jn,0) * weight); // peso de contorno => integral de contorno
452  }
453  }
454  break;
455  case 3: // outflow condition
456  int id, il, jl;
457  REAL normal[3];
458  if (fDim == 1) PZError << __PRETTY_FUNCTION__ << " - ERROR! The normal vector is not available for 1D TPZInterpolatedElement\n";
459  if (fDim == 2){
460  normal[0] = axes(0,1);
461  //normal[1] = axes(1,1);
462  }
463  if (fDim == 3){
464  normal[0] = axes(0,2);
465  normal[1] = axes(1,2);
466  normal[2] = axes(2,2);
467  }
468  REAL ConvNormal = 0.;
469  for(id=0; id<fDim; id++) ConvNormal += fC*fConvDir[id]*normal[id];
470  if(ConvNormal > 0.) {
471  for(il=0; il<phr; il++) {
472  for(jl=0; jl<phr; jl++) {
473  ek(il,jl) += weight * ConvNormal * phi(il)*phi(jl);
474  }
475  }
476  }
477  else{
478  if (ConvNormal < 0.) std::cout << "Boundary condition error: inflow detected in outflow boundary condition: ConvNormal = " << ConvNormal << "\n";
479  }
480  break;
481  }
482 
483  if (this->IsSymetric()) {//only 1.e-3 because of bignumbers.
484  if ( !ek.VerifySymmetry( 1.e-3 ) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
485  }
486 }
487 
489 int TPZMatPoisson3d::VariableIndex(const std::string &name){
490  if(!strcmp("Solution",name.c_str())) return 1;
491  if(!strcmp("Derivative",name.c_str())) return 2;
492  if(!strcmp("KDuDx",name.c_str())) return 3;
493  if(!strcmp("KDuDy",name.c_str())) return 4;
494  if(!strcmp("KDuDz",name.c_str())) return 5;
495  if(!strcmp("NormKDu",name.c_str())) return 6;
496  if(!strcmp("MinusKGradU",name.c_str())) return 7;
497  if(!strcmp("POrder",name.c_str())) return 8;
498  if(!strcmp("Laplac",name.c_str())) return 9;
499  if(!strcmp("Stress",name.c_str())) return 10;
500  if(!strcmp("Flux",name.c_str())) return 10;
501  if(!strcmp("Pressure",name.c_str())) return 11;
502 
503  if (!strcmp("ExactPressure", name.c_str())) return 12;
504  if(!strcmp("ExactSolution",name.c_str())) return 12;
505  if(!strcmp("ExactFlux",name.c_str())) return 13;
506  if(!strcmp("Divergence",name.c_str())) return 14;
507  if(!strcmp("ExactDiv",name.c_str())) return 15;
508 
509  if(!strcmp("PressureOmega1",name.c_str())) return 16;
510  if(!strcmp("PressureOmega2",name.c_str())) return 17;
511  if(!strcmp("FluxOmega1",name.c_str())) return 18;
512 
513  if(!strcmp("GradFluxX",name.c_str())) return 19;
514  if(!strcmp("GradFluxY",name.c_str())) return 20;
515  if(!strcmp("FluxL2",name.c_str())) return 21;//Only To calculate l2 error
516  if(!strcmp("OrdemP",name.c_str())) return 99;
517  return TPZMaterial::VariableIndex(name);
518 }
519 
520 int TPZMatPoisson3d::NSolutionVariables(int var){
521  if(var == 1) return 1;
522  if(var == 2) return fDim;//arrumar o fluxo de hdiv para ser fdim tbem enquanto isso faco isso
523  if ((var == 3) || (var == 4) || (var == 5) || (var == 6)) return 1;
524  if (var == 7) return fDim;
525  if (var == 8) return 1;
526  if (var == 9) return 1;
527  if (var==10) return fDim;
528  if (var==11) return 1;
529 
530  if (var==12) return 1;
531  if (var==13) return fDim;
532  if (var==14) return 1;
533  if (var==15) return 1;
534  //teste de acoplamento
535  if (var==16) return 1;
536  if (var==17) return 1;
537  if (var==18) return 3;
538  if (var==19) return 3;
539  if (var==20) return 3;
540  if (var==21) return fDim;
541 
542 
543  return TPZMaterial::NSolutionVariables(var);
544 }
545 
546 void TPZMatPoisson3d::Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout){
547 
548  TPZVec<STATE> pressure(1);
549  TPZVec<REAL> pto(3);
550  TPZFMatrix<STATE> flux(3,1);
551 
552  int numbersol = data.sol.size();
553  if (numbersol != 1) {
554  DebugStop();
555  }
556 
557  // Solout.Resize(this->NSolutionVariables(var));
558 
559 #ifndef STATE_COMPLEX
560 
561  switch (var) {
562  /* case 7:
563  {
564  // { //MinusKGradU
565  int id;
566  TPZManVector<STATE> dsolp(3,0);
567  dsolp[0] = data.dsol[0](0,0)*data.axes(0,0)+data.dsol[0](1,0)*data.axes(1,0);
568  dsolp[1] = data.dsol[0](0,0)*data.axes(0,1)+data.dsol[0](1,0)*data.axes(1,1);
569  dsolp[2] = data.dsol[0](0,0)*data.axes(0,2)+data.dsol[0](1,0)*data.axes(1,2);
570  for(id=0 ; id<fDim; id++)
571  {
572  Solout[id] = -1. * this->fK * dsolp[id];
573  }
574  }
575  break;*/
576  case 8:
577  Solout[0] = data.p;
578  break;
579  case 10:
580  if (data.numberdualfunctions) {
581 
582  Solout[0]=data.sol[0][0];
583  Solout[1]=data.sol[0][1];
584  Solout[2]=data.sol[0][2];
585 
586  }
587  else {
588  this->Solution(data.sol[0], data.dsol[0], data.axes, 2, Solout);
589  }
590 
591  break;
592 
593  case 21:
594  for(int k=0;k<fDim;k++){
595  Solout[k]=data.sol[0][k];
596  }
597  break;
598 
599  case 11:
600  if (data.numberdualfunctions) {
601  Solout[0]=data.sol[0][2];
602  }
603  else{
604  Solout[0]=data.sol[0][0];
605  }
606  break;
607 
608  case 12:
609  fForcingFunctionExact->Execute(data.x,pressure,flux);
610 
611  Solout[0]=pressure[0];
612  break;
613  case 13:
614  fForcingFunctionExact->Execute(data.x,pressure,flux);
615 
616  Solout[0]=flux(0,0);
617  Solout[1]=flux(1,0);
618  break;
619 
620  case 14:
621  {
622  if (data.numberdualfunctions){
623  Solout[0]=data.sol[0][data.sol[0].NElements()-1];
624  }else{
625  //Solout[0]=data.dsol[0](0,0)+data.dsol[0](1,1)+data.dsol[0](2,2);
626  STATE val = 0.;
627  for(int i=0; i<fDim; i++){
628  val += data.dsol[0](i,i);
629  }
630  Solout[0] = val;
631  }
632  }
633  break;
634 
635  case 15:
636  {
637  fForcingFunctionExact->Execute(data.x,pressure,flux);
638  Solout[0]=flux(fDim,0);
639  }
640  break;
641 
642 
643  case 16:
644  if (data.numberdualfunctions) {
645  Solout[0]=data.sol[0][2];
646  }
647  else {
648  std::cout<<"Pressao somente em Omega1"<<std::endl;
649  Solout[0]=0;//NULL;
650  }
651 
652  break;
653 
654  case 17:
655  if (!data.numberdualfunctions) {
656  Solout[0]=data.sol[0][0];
657  }
658  else {
659  std::cout<<"Pressao somente em omega2"<<std::endl;
660  Solout[0]=0;//NULL;
661  }
662 
663  break;
664  case 18:
665  if( data.numberdualfunctions){
666  Solout[0]=data.sol[0][0];//fluxo de omega1
667  Solout[1]=data.sol[0][1];
668  // Solout[2]=data.sol[2];
669  return;
670  }
671 
672  case 19:
673  if(data.numberdualfunctions){
674  Solout[0]=data.dsol[0](0,0);//fluxo de omega1
675  Solout[1]=data.dsol[0](1,0);
676  Solout[2]=data.dsol[0](2,0);
677  return;
678  }
679  case 20:
680  if( data.numberdualfunctions){
681  Solout[0]=data.dsol[0](0,1);//fluxo de omega1
682  Solout[1]=data.dsol[0](1,1);
683  Solout[2]=data.dsol[0](2,1);
684  return;
685  }
686  else {
687  std::cout<<"Pressao somente em omega2"<<std::endl;
688  Solout[0]=0;//NULL;
689  }
690  break;
691  default:
692 
693  if (data.sol[0].size() == 4) {
694 
695  data.sol[0][0] = data.sol[0][2];
696  }
697 
698  this->Solution(data.sol[0], data.dsol[0], data.axes, var, Solout);
699  break;
700  }
701 #endif
702 }
703 
704 #include "pzaxestools.h"
705 void TPZMatPoisson3d::Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,TPZVec<STATE> &Solout){
706 
707 #ifndef STATE_COMPLEX
708  Solout.Resize( this->NSolutionVariables( var ) );
709 
710  if(var == 1){
711  Solout[0] = Sol[0];//function
712  return;
713  }
714  if(var == 2) {
715  int id;
716  for(id=0 ; id<fDim; id++) {
717  TPZFNMatrix<9,STATE> dsoldx;
718  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
719  Solout[id] = dsoldx(id,0);//derivate
720  }
721  return;
722  }//var == 2
723  if (var == 3){ //KDuDx
724  TPZFNMatrix<9,STATE> dsoldx;
725  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
726  Solout[0] = dsoldx(0,0) * this->fK;
727  return;
728  }//var ==3
729  if (var == 4){ //KDuDy
730  TPZFNMatrix<9,STATE> dsoldx;
731  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
732  Solout[0] = dsoldx(1,0) * this->fK;
733  return;
734  }//var == 4
735  if (var == 5){ //KDuDz
736  TPZFNMatrix<9,STATE> dsoldx;
737  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
738  Solout[0] = dsoldx(2,0) * this->fK;
739  return;
740  }//var == 5
741  if (var == 6){ //NormKDu
742  int id;
743  REAL val = 0.;
744  for(id=0 ; id<fDim; id++){
745  val += (DSol(id,0) * this->fK) * (DSol(id,0) * this->fK);
746  }
747  Solout[0] = sqrt(val);
748  return;
749  }//var == 6
750  if (var == 7){ //MinusKGradU
751  int id;
752  //REAL val = 0.;
753  TPZFNMatrix<9,STATE> dsoldx;
754  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
755  for(id=0 ; id<fDim; id++) {
756  Solout[id] = -1. * this->fK * dsoldx(id,0);
757  }
758  return;
759  }//var == 7
760  if(var == 9){//Laplac
761  Solout.Resize(1);
762  Solout[0] = DSol(2,0);
763  return;
764  }//Laplac
765 
766 #endif
767  TPZMaterial::Solution(Sol, DSol, axes, var, Solout);
768 
769 }//method
770 
771 void TPZMatPoisson3d::Flux(TPZVec<REAL> &/*x*/, TPZVec<STATE> &/*Sol*/, TPZFMatrix<STATE> &/*DSol*/, TPZFMatrix<REAL> &/*axes*/, TPZVec<STATE> &/*flux*/) {
772  //Flux(TPZVec<REAL> &x, TPZVec<REAL> &Sol, TPZFMatrix<REAL> &DSol, TPZFMatrix<REAL> &axes, TPZVec<REAL> &flux)
773 }
774 void TPZMatPoisson3d::ErrorsHdiv(TPZMaterialData &data,TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values){
775 
776  values.Fill(0.0);
777  TPZVec<STATE> sol(1,0.),dsol(fDim,0.),div(1,0.);
778  if(data.numberdualfunctions) Solution(data,11,sol);//pressao
779  Solution(data,21,dsol);//fluxo
780  Solution(data,14,div);//divergente
781 
782 #ifdef LOG4CXX
783  if(logger->isDebugEnabled()){
784  std::stringstream sout;
785  sout<< "\n";
786  sout << " Pto " << data.x << std::endl;
787  sout<< " pressao exata " <<u_exact <<std::endl;
788  sout<< " pressao aprox " <<sol <<std::endl;
789  sout<< " ---- "<<std::endl;
790  sout<< " fluxo exato " <<du_exact(0,0)<<", " << du_exact(1,0)<<std::endl;
791  sout<< " fluxo aprox " <<dsol<<std::endl;
792  sout<< " ---- "<<std::endl;
793  if(du_exact.Rows()>fDim) sout<< " div exato " <<du_exact(2,0)<<std::endl;
794  sout<< " div aprox " <<div<<std::endl;
795  LOGPZ_DEBUG(logger,sout.str())
796  }
797 #endif
798 
799 
800  //values[0] : pressure error using L2 norm
801  if(data.numberdualfunctions){
802  REAL diffP = abs(u_exact[0]-sol[0]);
803  values[0] = diffP*diffP;
804  }
805  //values[1] : flux error using L2 norm
806  for(int id=0; id<fDim; id++) {
807  REAL diffFlux = abs(dsol[id] - du_exact(id,0));
808  values[1] += diffFlux*diffFlux;
809  }
810  if(du_exact.Rows()>fDim){
811  //values[2] : divergence using L2 norm
812  REAL diffDiv = abs(div[0] - du_exact(fDim,0));
813  values[2]=diffDiv*diffDiv;
814  //values[3] : Hdiv norm => values[1]+values[2];
815  values[3]= values[1]+values[2];
816  }
817 }
818 
819 void TPZMatPoisson3d::Errors(TPZVec<REAL> &x,TPZVec<STATE> &u,
820  TPZFMatrix<STATE> &dudx, TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
821  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) {
822 
823  values.Resize(NEvalErrors());
824  values.Fill(0.0);
825  TPZManVector<STATE> dudxEF(1,0.), dudyEF(1,0.),dudzEF(1,0.);
826  this->Solution(u,dudx,axes,VariableIndex("KDuDx"), dudxEF);
827  STATE fraq = dudxEF[0]/fK;
828  fraq = fraq - du_exact(0,0);
829  REAL diff = fabs(fraq);
830  values[3] = diff*diff;
831  if(fDim > 1) {
832  this->Solution(u,dudx,axes, this->VariableIndex("KDuDy"), dudyEF);
833  fraq = dudyEF[0]/fK;
834  fraq = fraq - du_exact(1,0);
835  diff = fabs(fraq);
836  values[4] = diff*diff;
837  if(fDim > 2) {
838  this->Solution(u,dudx,axes, this->VariableIndex("KDuDz"), dudzEF);
839  fraq = dudzEF[0]/fK;
840  fraq = fraq - du_exact(2,0);
841  diff = fabs(fraq);
842  values[5] = diff*diff;
843  }
844  }
845 
846  TPZManVector<STATE,3> sol(1),dsol(3,0.);
847  Solution(u,dudx,axes,1,sol);
848  Solution(u,dudx,axes,2,dsol);
849  int id;
850  //values[1] : eror em norma L2
851  diff = fabs(sol[0] - u_exact[0]);
852  values[1] = diff*diff;
853  //values[2] : erro em semi norma H1
854  values[2] = 0.;
855  for(id=0; id<fDim; id++) {
856  diff = fabs(dsol[id] - du_exact(id,0));
857  values[2] += abs(fK)*diff*diff;
858  }
859  //values[0] : erro em norma H1 <=> norma Energia
860  values[0] = values[1]+values[2];
861 }
862 
863 void TPZMatPoisson3d::BCInterfaceJump(TPZVec<REAL> &x, TPZSolVec &leftu,TPZBndCond &bc,TPZSolVec & jump){
864  int numbersol = leftu.size();
865  for (int is=0; is<numbersol ; is++) {
866  jump[is].Resize(1);
867  if(bc.Type() == 0){ //DIRICHLET
868  STATE f = bc.Val2()(0,0);
869  jump[is][0] = leftu[is][0] - f;
870  }
871  else{
872  jump[is].Fill(0.);
873  }
874  }
875 }//method
876 
877 #ifdef _AUTODIFF
878 void TPZMatPoisson3d::ContributeEnergy(TPZVec<REAL> &x,
879  TPZVec<FADFADREAL> &sol,
880  TPZVec<FADFADREAL> &dsol,
881  FADFADREAL &U,
882  REAL weight)
883 {
884  int dim = dsol.NElements()/sol.NElements();
885 
886  //Equa�o de Poisson
887  if(sol.NElements() != 1) PZError << "";
888  REAL vartocast = 0.;
889 #ifdef STATE_COMPLEX
890  vartocast = fXf.real();
891 #else
892  vartocast = fXf;
893 #endif
894  U+= sol[0] * FADREAL(weight * vartocast);
895 
896 #ifdef STATE_COMPLEX
897  vartocast = fK.real();
898 #else
899  vartocast = fK;
900 #endif
901  switch(dim)
902  {
903  case 1:
904  U+=vartocast*(dsol[0] * dsol[0])*FADREAL(weight/2.); // U=((du/dx)^2)/2
905 
906  break;
907  case 2:
908  U+=vartocast*(dsol[0] * dsol[0] +
909  dsol[1] * dsol[1])*(weight/2.); // U=((du/dx)^2+(du/dy)^2)/2
910  /*Buff = dsol[0] * dsol[0];
911  Buff += dsol[1] * dsol[1];
912  U += Buff * FADREAL(weight/2.); // U=((du/dx)^2+(du/dy)^2)/2*/
913  break;
914  case 3:
915  U+=vartocast*(dsol[0] * dsol[0] + dsol[1] * dsol[1] +
916  dsol[2] * dsol[2])*(weight/2.); // U=((du/dx)^2+(du/dy)^2+(du/dz)^2)/2*/
917  /*Buff = dsol[0] * dsol[0];
918  Buff += dsol[1] * dsol[1];
919  Buff += dsol[2] * dsol[2];
920  U += Buff * FADREAL(weight/2.); // U=((du/dx)^2+(du/dy)^2+(du/dz)^2)/2*/
921  break;
922  }
923 }
924 
925 void TPZMatPoisson3d::ContributeBCEnergy(TPZVec<REAL> & x,TPZVec<FADFADREAL> & sol, FADFADREAL &U, REAL weight, TPZBndCond &bc) {
926  REAL vartocast = 0.;
927 #ifdef STATE_COMPLEX
928  vartocast = bc.Val2()(0,0).real();
929 #else
930  vartocast = bc.Val2()(0,0);
931 #endif
932  FADFADREAL solMinBC = sol[0] - FADREAL(vartocast);
933 
934 
935  switch (bc.Type()) {
936  case 0 : // Dirichlet condition
937  // U += 1/2* Big * weight * Integral((u - u0)^2 dOmega)
938  U += (solMinBC * solMinBC) * FADREAL(weight * gBigNumber / 2.);
939  break;
940  case 1 : // Neumann condition
941  // U -= weight * Integral([g].u dOmega)
942  U -= sol[0] * FADREAL(vartocast*weight);
943  break;
944  case 2 : // condi�o mista
945 #ifdef STATE_COMPLEX
946  vartocast = bc.Val1()(0,0).real();
947 #else
948  vartocast = bc.Val1()(0,0);
949 #endif
950  // U += 1/2 * weight * Integral(<(u-u0), [g].(u-u0)> dOmega)
951  U += ( solMinBC * /*scalar*/ FADREAL(vartocast) * /*matrix oprt*/ solMinBC ) * FADREAL(weight / 2.);
952  break;
953 
954  }
955 }
956 
957 #endif
958 
959 
960 void TPZMatPoisson3d::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright,
961  REAL weight,
962  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef){
963 
964  TPZFMatrix<REAL> &dphiLdAxes = dataleft.dphix;
965  TPZFMatrix<REAL> &dphiRdAxes = dataright.dphix;
966  TPZFMatrix<REAL> &phiL = dataleft.phi;
967  TPZFMatrix<REAL> &phiR = dataright.phi;
968  TPZManVector<REAL,3> &normal = data.normal;
969 
970  TPZFNMatrix<660> dphiL, dphiR;
971  TPZAxesTools<REAL>::Axes2XYZ(dphiLdAxes, dphiL, dataleft.axes);
972  TPZAxesTools<REAL>::Axes2XYZ(dphiRdAxes, dphiR, dataright.axes);
973 
974  int &LeftPOrder=dataleft.p;
975  int &RightPOrder=dataright.p;
976 
977  REAL &faceSize=data.HSize;
978 
979 
980  int nrowl = phiL.Rows();
981  int nrowr = phiR.Rows();
982  int il,jl,ir,jr,id;
983 
984  //Convection term
985  REAL ConvNormal = 0.;
986  for(id=0; id<fDim; id++) ConvNormal += fC * fConvDir[id] * normal[id];
987  if(ConvNormal > 0.) {
988  for(il=0; il<nrowl; il++) {
989  for(jl=0; jl<nrowl; jl++) {
990  ek(il,jl) += weight * ConvNormal * phiL(il)*phiL(jl);
991  }
992  }
993  for(ir=0; ir<nrowr; ir++) {
994  for(jl=0; jl<nrowl; jl++) {
995  ek(ir+nrowl,jl) -= weight * ConvNormal * phiR(ir) * phiL(jl);
996  }
997  }
998  } else {
999  for(ir=0; ir<nrowr; ir++) {
1000  for(jr=0; jr<nrowr; jr++) {
1001  ek(ir+nrowl,jr+nrowl) -= weight * ConvNormal * phiR(ir) * phiR(jr);
1002  }
1003  }
1004  for(il=0; il<nrowl; il++) {
1005  for(jr=0; jr<nrowr; jr++) {
1006  ek(il,jr+nrowl) += weight * ConvNormal * phiL(il) * phiR(jr);
1007  }
1008  }
1009  }
1010 
1011  if(IsZero(fK)) return;
1012  //diffusion term
1013  STATE leftK, rightK;
1014  leftK = this->fK;
1015  rightK = this->fK;
1016 
1017  // 1) phi_I_left, phi_J_left
1018  for(il=0; il<nrowl; il++) {
1019  REAL dphiLinormal = 0.;
1020  for(id=0; id<fDim; id++) {
1021  dphiLinormal += dphiL(id,il)*normal[id];
1022  }
1023  for(jl=0; jl<nrowl; jl++) {
1024  REAL dphiLjnormal = 0.;
1025  for(id=0; id<fDim; id++) {
1026  dphiLjnormal += dphiL(id,jl)*normal[id];
1027  }
1028  ek(il,jl) += (STATE)(weight * ( this->fSymmetry * (0.5)*dphiLinormal*phiL(jl,0)-(0.5)*dphiLjnormal*phiL(il,0))) * leftK;
1029  }
1030  }
1031 
1032  // 2) phi_I_right, phi_J_right
1033  for(ir=0; ir<nrowr; ir++) {
1034  REAL dphiRinormal = 0.;
1035  for(id=0; id<fDim; id++) {
1036  dphiRinormal += dphiR(id,ir)*normal[id];
1037  }
1038  for(jr=0; jr<nrowr; jr++) {
1039  REAL dphiRjnormal = 0.;
1040  for(id=0; id<fDim; id++) {
1041  dphiRjnormal += dphiR(id,jr)*normal[id];
1042  }
1043  ek(ir+nrowl,jr+nrowl) += (STATE)(weight * (this->fSymmetry * ((-0.5) * dphiRinormal * phiR(jr) ) + (0.5) * dphiRjnormal * phiR(ir))) * rightK;
1044  }
1045  }
1046 
1047  // 3) phi_I_left, phi_J_right
1048  for(il=0; il<nrowl; il++) {
1049  REAL dphiLinormal = 0.;
1050  for(id=0; id<fDim; id++) {
1051  dphiLinormal += dphiL(id,il)*normal[id];
1052  }
1053  for(jr=0; jr<nrowr; jr++) {
1054  REAL dphiRjnormal = 0.;
1055  for(id=0; id<fDim; id++) {
1056  dphiRjnormal += dphiR(id,jr)*normal[id];
1057  }
1058  ek(il,jr+nrowl) += (STATE)weight * ((STATE)fSymmetry * ((STATE)((-0.5) * dphiLinormal * phiR(jr)) * leftK ) - (STATE)((0.5) * dphiRjnormal * phiL(il))* rightK );
1059  }
1060  }
1061 
1062  // 4) phi_I_right, phi_J_left
1063  for(ir=0; ir<nrowr; ir++) {
1064  REAL dphiRinormal = 0.;
1065  for(id=0; id<fDim; id++) {
1066  dphiRinormal += dphiR(id,ir)*normal[id];
1067  }
1068  for(jl=0; jl<nrowl; jl++) {
1069  REAL dphiLjnormal = 0.;
1070  for(id=0; id<fDim; id++) {
1071  dphiLjnormal += dphiL(id,jl)*normal[id];
1072  }
1073  ek(ir+nrowl,jl) += (STATE)weight * (
1074  (STATE)(fSymmetry * (0.5) * dphiRinormal * phiL(jl)) * rightK + (STATE)((0.5) * dphiLjnormal * phiR(ir)) * leftK
1075  );
1076  }
1077  }
1078 
1079  if (this->IsSymetric()){
1080  if ( !ek.VerifySymmetry() ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
1081  }
1082 
1083  if (this->fPenaltyConstant == 0.) return;
1084 
1085  leftK = this->fK;
1086  rightK = this->fK;
1087 
1088 
1089 
1090  //penalty = <A p^2>/h
1091  REAL penalty = fPenaltyConstant * (0.5 * (abs(leftK)*LeftPOrder*LeftPOrder + abs(rightK)*RightPOrder*RightPOrder)) / faceSize;
1092 
1093  if (this->fPenaltyType == ESolutionPenalty || this->fPenaltyType == EBoth){
1094 
1095  // 1) left i / left j
1096  for(il=0; il<nrowl; il++) {
1097  for(jl=0; jl<nrowl; jl++) {
1098  ek(il,jl) += weight * penalty * phiL(il,0) * phiL(jl,0);
1099  }
1100  }
1101 
1102  // 2) right i / right j
1103  for(ir=0; ir<nrowr; ir++) {
1104  for(jr=0; jr<nrowr; jr++) {
1105  ek(ir+nrowl,jr+nrowl) += weight * penalty * phiR(ir,0) * phiR(jr,0);
1106  }
1107  }
1108 
1109  // 3) left i / right j
1110  for(il=0; il<nrowl; il++) {
1111  for(jr=0; jr<nrowr; jr++) {
1112  ek(il,jr+nrowl) += -1.0 * weight * penalty * phiR(jr,0) * phiL(il,0);
1113  }
1114  }
1115 
1116  // 4) right i / left j
1117  for(ir=0; ir<nrowr; ir++) {
1118  for(jl=0; jl<nrowl; jl++) {
1119  ek(ir+nrowl,jl) += -1.0 * weight * penalty * phiL(jl,0) * phiR(ir,0);
1120  }
1121  }
1122 
1123  }
1124 
1125  if (this->fPenaltyType == EFluxPenalty || this->fPenaltyType == EBoth){
1126 
1127  REAL NormalFlux_i = 0.;
1128  REAL NormalFlux_j = 0.;
1129 
1130  // 1) left i / left j
1131  for(il=0; il<nrowl; il++) {
1132  NormalFlux_i = 0.;
1133  for(id=0; id<fDim; id++) {
1134  NormalFlux_i += dphiL(id,il)*normal[id];
1135  }
1136  for(jl=0; jl<nrowl; jl++) {
1137  NormalFlux_j = 0.;
1138  for(id=0; id<fDim; id++) {
1139  NormalFlux_j += dphiL(id,jl)*normal[id];
1140  }
1141  ek(il,jl) += (STATE)(weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * leftK;
1142  }
1143  }
1144 
1145  // 2) right i / right j
1146  for(ir=0; ir<nrowr; ir++) {
1147  NormalFlux_i = 0.;
1148  for(id=0; id<fDim; id++) {
1149  NormalFlux_i += dphiR(id,ir)*normal[id];
1150  }
1151  for(jr=0; jr<nrowr; jr++) {
1152  NormalFlux_j = 0.;
1153  for(id=0; id<fDim; id++) {
1154  NormalFlux_j += dphiR(id,jr)*normal[id];
1155  }
1156  ek(ir+nrowl,jr+nrowl) += (STATE)(weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * rightK;
1157  }
1158  }
1159 
1160  // 3) left i / right j
1161  for(il=0; il<nrowl; il++) {
1162  NormalFlux_i = 0.;
1163  for(id=0; id<fDim; id++) {
1164  NormalFlux_i += dphiL(id,il)*normal[id];
1165  }
1166  for(jr=0; jr<nrowr; jr++) {
1167  NormalFlux_j = 0.;
1168  for(id=0; id<fDim; id++) {
1169  NormalFlux_j += dphiR(id,jr)*normal[id];
1170  }
1171  ek(il,jr+nrowl) += (STATE)((-1.) * weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * rightK;
1172  }
1173  }
1174 
1175  // 4) right i / left j
1176  for(ir=0; ir<nrowr; ir++) {
1177  NormalFlux_i = 0.;
1178  for(id=0; id<fDim; id++) {
1179  NormalFlux_i += dphiR(id,ir)*normal[id];
1180  }
1181  for(jl=0; jl<nrowl; jl++) {
1182  NormalFlux_j = 0.;
1183  for(id=0; id<fDim; id++) {
1184  NormalFlux_j += dphiL(id,jl)*normal[id];
1185  }
1186  ek(ir+nrowl,jl) += (STATE)((-1.) * weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * leftK;
1187  }
1188  }
1189 
1190  }
1191 
1192 }
1193 
1195 void TPZMatPoisson3d::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
1196  REAL weight,
1197  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc){
1198 
1199  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
1200  TPZFMatrix<REAL> &phiL = dataleft.phi;
1201  TPZManVector<REAL,3> &normal = data.normal;
1202  int POrder= dataleft.p;
1203  REAL faceSize=data.HSize;
1204 
1205  // cout << "Material Id " << bc.Id() << " normal " << normal << "\n";
1206  int il,jl,nrowl,id;
1207  nrowl = phiL.Rows();
1208  REAL ConvNormal = 0.;
1209  for(id=0; id<fDim; id++) ConvNormal += fC*fConvDir[id]*normal[id];
1210  switch(bc.Type()) {
1211  case 0: // DIRICHLET
1212 
1213  //Diffusion
1214  for(il=0; il<nrowl; il++) {
1215  REAL dphiLinormal = 0.;
1216  for(id=0; id<fDim; id++) {
1217  dphiLinormal += dphiL(id,il)*normal[id];
1218  }
1219  ef(il,0) += (STATE)(weight*dphiLinormal*fSymmetry)*fK*bc.Val2()(0,0);
1220  for(jl=0; jl<nrowl; jl++) {
1221  REAL dphiLjnormal = 0.;
1222  for(id=0; id<fDim; id++) {
1223  dphiLjnormal += dphiL(id,jl)*normal[id];
1224  }
1225  ek(il,jl) += (STATE)(weight*(fSymmetry * dphiLinormal * phiL(jl,0) - dphiLjnormal * phiL(il,0)))*fK;
1226  }
1227  }
1228 
1229  //Convection
1230  if(ConvNormal > 0.) {
1231  for(il=0; il<nrowl; il++) {
1232  for(jl=0; jl<nrowl; jl++) {
1233  ek(il,jl) += weight * ConvNormal * phiL(il)*phiL(jl);
1234  }
1235  }
1236  } else {
1237  for(il=0; il<nrowl; il++) {
1238  ef(il,0) -= (STATE)(weight * ConvNormal * phiL(il)) * bc.Val2()(0,0);
1239  }
1240  }
1241 
1242  break;
1243 
1244  case 1: // Neumann
1245  for(il=0; il<nrowl; il++) {
1246  ef(il,0) += (STATE)(weight*phiL(il,0))*bc.Val2()(0,0);
1247  }
1248  break;
1249 
1250  case 3: // outflow condition
1251  if(ConvNormal > 0.) {
1252  for(il=0; il<nrowl; il++) {
1253  for(jl=0; jl<nrowl; jl++) {
1254  ek(il,jl) += weight * ConvNormal * phiL(il)*phiL(jl);
1255  }
1256  }
1257  }
1258  else {
1259  if (ConvNormal < 0.) std::cout << "Boundary condition error: inflow detected in outflow boundary condition: ConvNormal = " << ConvNormal << "\n";
1260  }
1261  break;
1262 
1263  default:
1264  PZError << __PRETTY_FUNCTION__ << " - Wrong boundary condition type\n";
1265  break;
1266  }
1267  if (this->IsSymetric()){
1268  if ( !ek.VerifySymmetry() ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
1269  }
1270 
1271  if (this->fPenaltyConstant == 0.) return;
1272 
1273  if (this->fPenaltyType == ESolutionPenalty || this->fPenaltyType == EBoth){
1274  nrowl = phiL.Rows();
1275  const REAL penalty = fPenaltyConstant * abs(fK) * POrder * POrder / faceSize; //Ap^2/h
1276  REAL outflow = 0.;
1277  for(il=0; il<fDim; il++) outflow += fC * fConvDir[il] * normal[il];
1278 
1279 
1280  switch(bc.Type()) {
1281  case 0: // DIRICHLET
1282  for(il=0; il<nrowl; il++) {
1283  ef(il,0) += (STATE)(weight * penalty * phiL(il,0)) * bc.Val2()(0,0);
1284  for(jl=0; jl<nrowl; jl++) {
1285  ek(il,jl) += weight * penalty * phiL(il,0) * phiL(jl,0);
1286  }
1287  }
1288 
1289  break;
1290  case 1: // Neumann
1291  if(outflow > 0.)
1292  {
1293  for(il=0; il<nrowl; il++)
1294  {
1295  for(jl=0; jl<nrowl; jl++)
1296  {
1297  ek(il,jl) += weight * outflow * phiL(il,0) * phiL(jl,0);
1298  }
1299  }
1300  }
1301  //nothing to be done
1302  break;
1303  default:
1304  PZError << "TPZMatPoisson3d::Wrong boundary condition type\n";
1305  break;
1306  }
1307 
1308  }
1309 
1310 }
1311 
1312 void TPZMatPoisson3d::InterfaceErrors(TPZVec<REAL> &/*x*/,
1313  TPZVec<STATE> &leftu, TPZFMatrix<STATE> &leftdudx, /* TPZFMatrix<REAL> &leftaxes,*/
1314  TPZVec<STATE> &rightu, TPZFMatrix<STATE> &rightdudx, /* TPZFMatrix<REAL> &rightaxes,*/
1315  TPZVec<STATE> &/*flux*/,
1316  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values,
1317  TPZVec<STATE> normal, STATE elsize) {
1318  // #warning Metodo nao funcional
1319  TPZManVector<STATE,3> Lsol(1), Ldsol(3,0.), Rsol(1), Rdsol(3,0.);
1320 
1321  TPZFMatrix<REAL> fake_axes(fDim,fDim,0.);
1322 
1323  Solution(leftu,leftdudx,fake_axes,1,Lsol);
1324  Solution(leftu,leftdudx,fake_axes,2,Ldsol);
1325 
1326  Solution(rightu,rightdudx,fake_axes,1,Rsol);
1327  Solution(rightu,rightdudx,fake_axes,2,Rdsol);
1328 
1329 #ifdef PZDEBUG
1330  if ( (leftdudx.Rows() != rightdudx.Rows()) || (leftdudx.Rows() != du_exact.Rows()) ){
1331  PZError << "TPZMatPoisson3d::InterfaceErrors - Left and right matrices should have"
1332  << endl
1333  << "same sizes in internal boundaries."
1334  << endl;
1335  exit (-1);
1336  }
1337 #endif
1338 
1339  STATE Ldsolnormal = 0., Rdsolnormal = 0., ExactDNormal = 0.;
1340  for(int id = 0; id < fDim; id++) {
1341  Ldsolnormal += Ldsol[id] * normal[id];
1342  Rdsolnormal += Rdsol[id] * normal[id];
1343  ExactDNormal += du_exact(id, 0) * normal[id];
1344  }
1345 
1346  values.Resize(3);
1347  STATE aux;
1348 
1349  //values[1] : eror em norma L2
1350 
1351  //Jump aprox. solution - jump of exact solution i.e. zero
1352  aux = (Lsol[0] - Rsol[0]);
1353 
1354  //*= h ^ -gAlfa
1355  aux *= pow(elsize, (STATE(-1.)) * gAlfa);
1356  REAL auxnorm = abs(aux);
1357  values[1] = auxnorm * auxnorm;
1358 
1359  //values[2] : erro em semi norma H1
1360  values[2] = 0.;
1361 
1362  for(int id=0; id<fDim; id++) {
1363  //Normal gradient average <grad V> = 0.5 * (grad_left.n + grad_right.n)
1364  aux = STATE(0.5) * (Ldsolnormal + Rdsolnormal);
1365  //<grad V> - <grad exact> = <grad V> - grad exact
1366  aux = aux - ExactDNormal;
1367  //*= h ^ gAlfa
1368  aux *= pow(elsize, gAlfa);
1369  auxnorm = abs(aux);
1370  values[2] += auxnorm * auxnorm;
1371  }
1372  //values[0] : erro em norma H1 <=> norma Energia
1373  values[0] = values[1]+values[2];
1374 }
1375 
1376 REAL TPZMatPoisson3d::ComputeSquareResidual(TPZVec<REAL>& X, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol){
1377  // residual = -fK Laplac(u) + fC * div(fConvDir*u) - (-fXf)
1378  STATE fXfLoc = fXf;
1379  if(fForcingFunction) {
1380  TPZManVector<STATE> res(1);
1381  fForcingFunction->Execute(X,res);
1382  fXfLoc = res[0];
1383  }
1384 
1385  STATE laplacU = (STATE)0;
1386  STATE divBetaU = (STATE)0;
1387  if(this->Dimension() == 1){
1388  laplacU = dsol(1,0);
1389  divBetaU = (STATE)(this->fC * this->fConvDir[0]) * dsol(0,0);
1390  }
1391  if(this->Dimension() == 2){
1392  laplacU = dsol(2,0);
1393  divBetaU = (STATE)fC * ( (STATE)fConvDir[0] * dsol(0,0) + (STATE)fConvDir[1] * dsol(1,0) );
1394  }
1395 
1396  REAL result = abs(-this->fK * laplacU + divBetaU - (-fXfLoc));
1397  return (result*result);
1398 }
1399 
1400 void TPZMatPoisson3d::Write(TPZStream &buf, int withclassid) const{
1401  TPZDiscontinuousGalerkin::Write(buf, withclassid);
1402  buf.Write(&fXf, 1);
1403  buf.Write(&fDim, 1);
1404  buf.Write(&fK, 1);
1405  buf.Write(&fC, 1);
1406  buf.Write(fConvDir, 3);
1407  buf.Write(&fSymmetry, 1);
1408  buf.Write(&fSD, 1);
1409  buf.Write(&fPenaltyConstant,1);
1410  buf.Write(&gAlfa, 1);
1411 }
1412 
1413 void TPZMatPoisson3d::Read(TPZStream &buf, void *context){
1414  TPZDiscontinuousGalerkin::Read(buf, context);
1415  buf.Read(&fXf, 1);
1416  buf.Read(&fDim, 1);
1417  buf.Read(&fK, 1);
1418  buf.Read(&fC, 1);
1419  buf.Read(fConvDir, 3);
1420  buf.Read(&fSymmetry, 1);
1421  buf.Read(&fSD, 1);
1422  buf.Read(&fPenaltyConstant,1);
1423  buf.Read(&gAlfa, 1);
1424 }
1425 
1426 int TPZMatPoisson3d::ClassId() const{
1427  return Hash("TPZMatPoisson3d") ^ TPZDiscontinuousGalerkin::ClassId() << 1;
1428 }
1429 
1430 #ifndef BORLAND
1431 template class TPZRestoreClass<TPZMatPoisson3d>;
1432 #endif
TPZDiscontinuousGalerkin
Defines the interface which material objects need to implement for discontinuous Galerkin formulation...
Definition: pzdiscgal.h:20
TPZMaterial::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout)
Returns the solution associated with the var index based on the finite element approximation.
Definition: TPZMaterial.cpp:200
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZMatPoisson3d::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzpoisson3d.cpp:1400
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
TPZMatPoisson3d::ComputeSquareResidual
virtual REAL ComputeSquareResidual(TPZVec< REAL > &X, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol) override
Compute square of residual of the differential equation at one integration point. ...
Definition: pzpoisson3d.cpp:1376
TPZMatPoisson3d::Read
virtual void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzpoisson3d.cpp:1413
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
TPZMatPoisson3d::ContributeBCHDiv
virtual void ContributeBCHDiv(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzpoisson3d.cpp:348
IsZero
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
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.
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
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)
TPZMatPoisson3d::fSD
STATE fSD
Multiplication value for the streamline diffusion term.
Definition: pzpoisson3d.h:53
TPZDiscontinuousGalerkin::Read
void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzdiscgal.cpp:118
TPZMatPoisson3d::fPenaltyConstant
REAL fPenaltyConstant
Constant multiplyer of penalty term, when required is set.
Definition: pzpoisson3d.h:67
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
TPZMatPoisson3d::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: pzpoisson3d.cpp:960
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZMaterial::operator=
TPZMaterial & operator=(const TPZMaterial &copy)
operator =
Definition: TPZMaterial.cpp:62
TPZMatPoisson3d::GetParameters
void GetParameters(STATE &diff, REAL &conv, TPZVec< REAL > &convdir)
Definition: pzpoisson3d.cpp:85
TPZMaterial::VariableIndex
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Definition: TPZMaterial.cpp:141
pzerror.h
Defines PZError.
TPZMatPoisson3d::NEvalErrors
virtual int NEvalErrors() override
Returns the number of norm errors. Default is 3: energy, L2 and H1.
Definition: pzpoisson3d.h:290
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZMaterialData::jacinv
TPZFNMatrix< 9, REAL > jacinv
value of the inverse of the jacobian at the integration point
Definition: pzmaterialdata.h:67
abs
TinyFad< 8, T > abs(const TinyFad< 8, T > &in)
Definition: tinyfadeight.h:846
TPZMatPoisson3d::~TPZMatPoisson3d
virtual ~TPZMatPoisson3d()
Definition: pzpoisson3d.cpp:92
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZMatPoisson3d::Flux
virtual void Flux(TPZVec< REAL > &x, TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, TPZVec< STATE > &flux) override
Computes the value of the flux function to be used by ZZ error estimator.
Definition: pzpoisson3d.cpp:771
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
pzpoisson3d.h
Contains the TPZMatPoisson3d class.
TPZBndCond::Val2
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZMaterialData::phi
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
Definition: pzmaterialdata.h:55
TPZMaterialData::p
int p
maximum polinomial order of the shape functions
Definition: pzmaterialdata.h:75
TPZMatPoisson3d::fC
REAL fC
Variable which multiplies the convection term of the equation.
Definition: pzpoisson3d.h:40
TPZMatPoisson3d::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzpoisson3d.cpp:1426
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZMaterial::Print
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
Definition: TPZMaterial.cpp:103
pzelmat.h
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
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
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZMatPoisson3d::fXf
STATE fXf
Forcing function value.
Definition: pzpoisson3d.h:31
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
TPZMatPoisson3d::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ef) override
It computes a contribution to the residual vector at one integration point.
Definition: pzpoisson3d.h:199
test.f
f
Definition: test.py:287
TPZMatPoisson3d::ContributeBCInterface
virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Definition: pzpoisson3d.cpp:1195
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZMatPoisson3d::fK
STATE fK
Coeficient which multiplies the Laplacian operator.
Definition: pzpoisson3d.h:37
TPZMatPoisson3d::SetParameters
virtual void SetParameters(STATE diff, REAL conv, TPZVec< REAL > &convdir)
Definition: pzpoisson3d.cpp:78
TPZMatPoisson3d::ContributeHDiv
virtual void ContributeHDiv(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Compute the contribution at an integration point to the stiffness matrix of the HDiv formulation...
Definition: pzpoisson3d.cpp:199
TPZDiscontinuousGalerkin::Write
void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzdiscgal.cpp:114
TPZMatPoisson3d::InterfaceErrors
void InterfaceErrors(TPZVec< REAL > &, TPZVec< STATE > &leftu, TPZFMatrix< STATE > &leftdudx, TPZVec< STATE > &rightu, TPZFMatrix< STATE > &rightdudx, TPZVec< STATE > &, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &values, TPZVec< STATE > normal, STATE elsize)
Definition: pzpoisson3d.cpp:1312
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
TPZMaterialData::fNormalVec
TPZFNMatrix< 180 > fNormalVec
list of normal vectors
Definition: pzmaterialdata.h:101
TPZMatPoisson3d::ErrorsHdiv
void ErrorsHdiv(TPZMaterialData &data, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &values) override
Definition: pzpoisson3d.cpp:774
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
TPZMatPoisson3d::LocalNeumanContribute
virtual void LocalNeumanContribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzpoisson3d.cpp:278
TPZMaterialData::numberdualfunctions
int numberdualfunctions
number of dual function (e.g. pressure in HDiv approximations)
Definition: pzmaterialdata.h:93
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZMaterialData::axes
TPZFNMatrix< 9, REAL > axes
axes indicating the directions of the derivatives of the shapefunctions
Definition: pzmaterialdata.h:63
test.res
string res
Definition: test.py:151
TPZMatPoisson3d::SetNoPenalty
void SetNoPenalty()
Defines no penalty terms in ContributeInterface.
Definition: pzpoisson3d.h:70
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZMatPoisson3d::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: pzpoisson3d.cpp:489
pzmaterialdata.h
Contains the TPZMaterialData class which implements an interface between TPZCompEl::CalcStiff and TPZ...
TPZMaterialData::HSize
REAL HSize
measure of the size of the element
Definition: pzmaterialdata.h:83
TPZMaterial::HasForcingFunction
virtual int HasForcingFunction()
Directive that gives true if the material has a forcing function.
Definition: TPZMaterial.h:472
TPZMatPoisson3d::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: pzpoisson3d.h:207
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
TPZMatPoisson3d::gAlfa
static STATE gAlfa
Using in InterfaceErrors.
Definition: pzpoisson3d.h:82
TPZMatPoisson3d::fPenaltyType
EPenaltyType fPenaltyType
Penalty term definition.
Definition: pzpoisson3d.h:59
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
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZMatPoisson3d::BCInterfaceJump
virtual void BCInterfaceJump(TPZVec< REAL > &x, TPZSolVec &leftu, TPZBndCond &bc, TPZSolVec &jump) override
Computes interface jump from element to Dirichlet boundary condition.
Definition: pzpoisson3d.cpp:863
TPZMatPoisson3d::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: pzpoisson3d.cpp:520
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZMatPoisson3d
DESCRIBE PLEASE.
Definition: pzpoisson3d.h:26
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZMatrix::VerifySymmetry
virtual int VerifySymmetry(REAL tol=1.e-13) const
Checks if current matrix value is symmetric.
Definition: pzmatrix.cpp:1501
TPZMatPoisson3d::Errors
virtual void Errors(TPZVec< REAL > &x, TPZVec< STATE > &u, TPZFMatrix< STATE > &dudx, TPZFMatrix< REAL > &axes, TPZVec< STATE > &flux, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &values) override
Computes the error due to the difference between the interpolated flux and the flux computed based o...
Definition: pzpoisson3d.cpp:819
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::SetNonSymmetric
void SetNonSymmetric()
Set material elliptic term as the Baumann&#39;s formulation, i.e. the non-symmetrical formulation...
Definition: pzpoisson3d.h:117
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
TPZMatPoisson3d::fConvDir
REAL fConvDir[3]
Direction of the convection operator.
Definition: pzpoisson3d.h:43
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
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
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
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
rdt.values
def values
Definition: rdt.py:119
TPZMaterialData::fVecShapeIndex
TPZManVector< std::pair< int, int64_t > > fVecShapeIndex
correspondence between normal vector index and index of the shape functions
Definition: pzmaterialdata.h:99
TPZMatPoisson3d::fSymmetry
REAL fSymmetry
Symmetry coefficient of elliptic term.
Definition: pzpoisson3d.h:50
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZRestoreClass
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
TPZMatPoisson3d::Name
virtual std::string Name() override
Returns the name of the material.
Definition: pzpoisson3d.h:192
Generated on Sun Aug 16 2020 11:17:01 for NeoPZ by   doxygen 1.8.13