NeoPZ
TPZMatLaplacian.cpp
Go to the documentation of this file.
1 
6 #include "TPZMatLaplacian.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 #include "pzaxestools.h"
16 #include "pzextractval.h"
17 
18 #include <cmath>
19 
20 #ifdef LOG4CXX
21 static LoggerPtr logger(Logger::getLogger("pz.material.TPZMatLaplacian"));
22 #endif
23 
24 
25 using namespace std;
26 
27 TPZMatLaplacian::TPZMatLaplacian(int nummat, int dim) :
28 TPZRegisterClassId(&TPZMatLaplacian::ClassId), TPZDiscontinuousGalerkin(nummat), fXf(0.), fDim(dim), fTensorK(dim,dim,0.), fInvK(dim,dim,0.)
29 {
30  fK = 1.;
31  for (int i=0; i<dim; i++) {
32  fTensorK(i,i) = 1.;
33  fInvK(i,i) = 1.;
34  }
35  fPenaltyConstant = 1000.;
36  this->SetNonSymmetric();
37  this->SetNoPenalty();
38 }
39 
40 TPZMatLaplacian::TPZMatLaplacian()
41 : TPZRegisterClassId(&TPZMatLaplacian::ClassId), TPZDiscontinuousGalerkin(), fXf(0.), fDim(1), fTensorK(1,1,1.), fInvK(1,1,1.){
42  fK = 1.;
43  fPenaltyConstant = 1000.;
44  this->SetNonSymmetric();
45  this->SetNoPenalty();
46 }
47 
48 TPZMatLaplacian::TPZMatLaplacian(const TPZMatLaplacian &copy)
49 : TPZRegisterClassId(&TPZMatLaplacian::ClassId), TPZDiscontinuousGalerkin(copy)
50 {
51  this->operator =(copy);
52 }
53 
54 TPZMatLaplacian & TPZMatLaplacian::operator=(const TPZMatLaplacian &copy){
55  TPZDiscontinuousGalerkin::operator = (copy);
56  fXf = copy.fXf;
57  fDim = copy.fDim;
58  fK = copy.fK;
59  fTensorK = copy.fTensorK;
60  fInvK = copy.fInvK;
61  fSymmetry = copy.fSymmetry;
62  fPenaltyConstant = copy.fPenaltyConstant;
63  this->fPenaltyType = copy.fPenaltyType;
64  this->fPermeabilityFunction = copy.fPermeabilityFunction;
65  return *this;
66 }
67 
68 void TPZMatLaplacian::SetParameters(STATE diff, STATE f) {
69  fK = diff;
70  fTensorK.Zero();
71  fInvK.Zero();
72  for (int i=0; i<fDim; i++) {
73  fTensorK(i,i) = diff;
74  fInvK(i,i) = STATE(1.)/diff;
75  }
76  fXf = f;
77 }
78 
79 //void TPZMatLaplacian::GetParameters(STATE &diff, STATE &f) {
80 // diff = fK;
81 // f = fXf;
82 //}
83 
84 TPZMatLaplacian::~TPZMatLaplacian() {
85 }
86 
87 int TPZMatLaplacian::NStateVariables() const {
88  return 1;
89 }
90 
91 void TPZMatLaplacian::Print(std::ostream &out) {
92  out << "name of material : " << Name() << "\n";
93  out << "Laplace operator multiplier fK "<< fK << endl;
94  out << "Forcing vector fXf " << fXf << endl;
95  out << "Penalty constant " << fPenaltyConstant << endl;
96  out << "Base Class properties :";
97  TPZMaterial::Print(out);
98  out << "\n";
99 }
100 
101 void TPZMatLaplacian::Contribute(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
102 {
103 
104  if(data.numberdualfunctions)
105  {
106  ContributeHDiv(data , weight , ek, ef);
107 
108  return;
109  }
110 
111  TPZFMatrix<REAL> &phi = data.phi;
112  TPZFMatrix<REAL> &dphi = data.dphix;
113  TPZVec<REAL> &x = data.x;
114  // TPZFMatrix<REAL> &axes = data.axes;
115  // TPZFMatrix<REAL> &jacinv = data.jacinv;
116  int phr = phi.Rows();
117 
118  STATE fXfLoc = fXf;
119 
120  if(fForcingFunction) { // phi(in, 0) = phi_in
121  TPZManVector<STATE,1> res(1);
122  //TPZFMatrix<STATE> dres(Dimension(),1);
123  //fForcingFunction->Execute(x,res,dres); // dphi(i,j) = dphi_j/dxi
124  fForcingFunction->Execute(x,res);
125  fXfLoc = res[0];
126  }
127 
128  STATE KPerm = fK;
129  if (fPermeabilityFunction) {
130  TPZFNMatrix<9,STATE> perm, invperm;
131  TPZManVector<STATE,3> func;
132  TPZFNMatrix<18,STATE> dfunc(6,3,0.);
133  fPermeabilityFunction->Execute(x, func, dfunc);
134  KPerm = dfunc(0,0);
135  }
136 
137  //Equacao de Poisson
138  for( int in = 0; in < phr; in++ ) {
139  int kd;
140  ef(in, 0) += (STATE)weight * fXfLoc * (STATE)phi(in,0);
141  for( int jn = 0; jn < phr; jn++ ) {
142  //ek(in,jn) += (STATE)weight*((STATE)(phi(in,0)*phi(jn,0)));
143  for(kd=0; kd<fDim; kd++) {
144  ek(in,jn) += (STATE)weight*(KPerm*(STATE)(dphi(kd,in)*dphi(kd,jn)));
145  }
146  }
147  }
148 
149  if (this->IsSymetric()){
150  if ( !ek.VerifySymmetry(1.e-10) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
151  }
152 
153 }
154 
155 void TPZMatLaplacian::Contribute(TPZMaterialData &data,REAL weight, TPZFMatrix<STATE> &ef)
156 {
157 
158  if(data.numberdualfunctions)
159  {
160  DebugStop();
161 // ContributeHDiv(data , weight , ef);
162 
163  return;
164  }
165 
166  TPZFMatrix<REAL> &phi = data.phi;
167  TPZFMatrix<REAL> &dphi = data.dphix;
168  TPZVec<REAL> &x = data.x;
169  // TPZFMatrix<REAL> &axes = data.axes;
170  // TPZFMatrix<REAL> &jacinv = data.jacinv;
171  int phr = phi.Rows();
172 
173  STATE fXfLoc = fXf;
174 
175  if(fForcingFunction) { // phi(in, 0) = phi_in
176  TPZManVector<STATE,1> res(1);
177  //TPZFMatrix<STATE> dres(Dimension(),1);
178  //fForcingFunction->Execute(x,res,dres); // dphi(i,j) = dphi_j/dxi
179  fForcingFunction->Execute(x,res);
180  fXfLoc = res[0];
181  }
182 
183  //Equacao de Poisson
184  for( int in = 0; in < phr; in++ ) {
185  int kd;
186  ef(in, 0) += (STATE)weight * fXfLoc * (STATE)phi(in,0);
187  for(kd=0; kd<fDim; kd++) {
188  ef(in,0) -= (STATE)weight*(fK*(STATE)(dphi(kd,in)*data.dsol[0](kd,0)));
189  }
190  }
191 }
193 void TPZMatLaplacian::ContributeHDiv(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
194 {
201  //TPZVec<REAL> &x = data.x;
202  STATE fXfLoc = fXf;
203  if(fForcingFunction) { // phi(in, 0) = phi_in
204  TPZManVector<STATE> res(1);
205  fForcingFunction->Execute(data.x,res); // dphi(i,j) = dphi_j/dxi
206  fXfLoc = res[0];
207  }
208  int numvec = data.fVecShapeIndex.NElements();
209  int numdual = data.numberdualfunctions;
210  int numprimalshape = data.phi.Rows()-numdual;
211 
212  int i,j;
213  REAL kreal = 0.;
214 #ifdef STATE_COMPLEX
215  kreal = fK.real();
216 #else
217  kreal = fK;
218 #endif
219  REAL ratiok = 1./kreal;
220  for(i=0; i<numvec; i++)
221  {
222  int ivecind = data.fVecShapeIndex[i].first;
223  int ishapeind = data.fVecShapeIndex[i].second;
224  for (j=0; j<numvec; j++) {
225  int jvecind = data.fVecShapeIndex[j].first;
226  int jshapeind = data.fVecShapeIndex[j].second;
227  REAL prod = data.fNormalVec(0,ivecind)*data.fNormalVec(0,jvecind)+
228  data.fNormalVec(1,ivecind)*data.fNormalVec(1,jvecind)+
229  data.fNormalVec(2,ivecind)*data.fNormalVec(2,jvecind);//faz o produto escalar entre u e v--> Matriz A
230  ek(i,j) += weight*ratiok*data.phi(ishapeind,0)*data.phi(jshapeind,0)*prod;
231 
232 
233 
234  }
235  TPZFNMatrix<3> ivec(3,1);
236  ivec(0,0) = data.fNormalVec(0,ivecind);
237  ivec(1,0) = data.fNormalVec(1,ivecind);
238  ivec(2,0) = data.fNormalVec(2,ivecind);
239  TPZFNMatrix<3> axesvec(3,1);
240  data.axes.Multiply(ivec,axesvec);
241  int iloc;
242  REAL divwq = 0.;
243  for(iloc=0; iloc<fDim; iloc++)
244  {
245  divwq += axesvec(iloc,0)*data.dphix(iloc,ishapeind);
246  }
247  for (j=0; j<numdual; j++) {
248  REAL fact = (-1.)*weight*data.phi(numprimalshape+j,0)*divwq;//calcula o termo da matriz B^T e B
249  ek(i,numvec+j) += fact;
250  ek(numvec+j,i) += fact;//-div
251  }
252  }
253  for(i=0; i<numdual; i++)
254  {
255  ef(numvec+i,0) += (STATE)(weight*data.phi(numprimalshape+i,0))*fXfLoc;//calcula o termo da matriz f
256  }
257 
258 }
259 
260 void TPZMatLaplacian::ContributeBCHDiv(TPZMaterialData &data,REAL weight,
261  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
262  int numvec = data.phi.Rows();
263 
264  TPZFMatrix<REAL> &phi = data.phi;
265  TPZManVector<STATE,1> v2(1);
266  v2[0] = bc.Val2()(0,0);
267  if (bc.HasForcingFunction()) {
268  bc.ForcingFunction()->Execute(data.x, v2);
269  }
270 
271  switch (bc.Type()) {
272  case 1 : // Neumann condition
273  int i,j;
274  for(i=0; i<numvec; i++)
275  {
276  //int ishapeind = data.fVecShapeIndex[i].second;
277  ef(i,0)+= (STATE)(gBigNumber * phi(i,0) * weight) * v2[0];
278  for (j=0; j<numvec; j++) {
279  //int jshapeind = data.fVecShapeIndex[j].second;
280  ek(i,j) += gBigNumber * phi(i,0) * phi(j,0) * weight;
281  }
282  }
283  break;
284  case 0 :
285  {// Dirichlet condition
286  int in;
287  for(in = 0 ; in < numvec; in++) {
288  ef(in,0) += (STATE)((-1.)* phi(in,0) * weight)*v2[0];
289  }
290  }
291  break;
292  case 2 : // mixed condition
293  {
294  int in,jn;
295  for(in = 0 ; in < numvec; in++) {
296  //int ishapeind = data.fVecShapeIndex[in].second;
297  ef(in,0) += v2[0] * (STATE)(phi(in,0) * weight);
298  for (jn = 0; jn < numvec; jn++) {
299  // int jshapeind = data.fVecShapeIndex[jn].second;
300  ek(in,jn) += (STATE)(weight*phi(in,0)*phi(jn,0)) *bc.Val1()(0,0);
301  }
302  }
303  }
304  break;
305  case 3: // outflow condition
306  break;
307  }
308 
309  if (this->IsSymetric()) {//only 1.e-3 because of bignumbers.
310  if ( !ek.VerifySymmetry( 1.e-3 ) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
311  }
312 
313 
314 }
315 
316 void TPZMatLaplacian::ContributeBC(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
317 
318  ContributeBC(datavec[0],weight,ek,ef,bc);
319 }
320 
321 void TPZMatLaplacian::ContributeBC(TPZMaterialData &data,REAL weight,
322  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
323 
324  if(data.fShapeType == TPZMaterialData::EVecandShape || data.fShapeType == TPZMaterialData::EVecShape)
325  {
326 
327  ContributeBCHDiv(data , weight , ek, ef, bc);
328 
329 
330  return;
331  }
332 
333  TPZFMatrix<REAL> &phi = data.phi;
334  // TPZFMatrix<REAL> &axes = data.axes;
335  int phr = phi.Rows();
336  short in,jn;
337  STATE v2[1];
338  v2[0] = bc.Val2()(0,0);
339 
340  if(bc.HasForcingFunction()) { // phi(in, 0) = phi_in // JORGE 2013 01 26
341  TPZManVector<STATE> res(1);
342  TPZFNMatrix<3,STATE> dres(3,1);
343  bc.ForcingFunction()->Execute(data.x,res,dres); // dphi(i,j) = dphi_j/dxi
344  v2[0] = res[0];
345  }
346 
347  switch (bc.Type()) {
348  case 0 : // Dirichlet condition
349  for(in = 0 ; in < phr; in++) {
350  ef(in,0) += (STATE)(gBigNumber* phi(in,0) * weight) * v2[0];
351  for (jn = 0 ; jn < phr; jn++) {
352  ek(in,jn) += gBigNumber * phi(in,0) * phi(jn,0) * weight;
353  }
354  }
355  break;
356  case 1 : // Neumann condition
357  for(in = 0 ; in < phi.Rows(); in++) {
358  ef(in,0) += v2[0] * (STATE)(phi(in,0) * weight);
359  }
360  break;
361  case 2 : // mixed condition
362  for(in = 0 ; in < phi.Rows(); in++) {
363  ef(in, 0) += v2[0] * (STATE)(phi(in, 0) * weight);
364  for (jn = 0 ; jn < phi.Rows(); jn++) {
365  ek(in,jn) += bc.Val1()(0,0) * (STATE)(phi(in,0) * phi(jn,0) * weight); // peso de contorno => integral de contorno
366  }
367  }
368  break;
369  }
370 
371  if (this->IsSymetric()) {//only 1.e-3 because of bignumbers.
372  if ( !ek.VerifySymmetry( 1.e-3 ) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
373  }
374 }
375 
377 int TPZMatLaplacian::VariableIndex(const std::string &name){
378  if(!strcmp("Solution",name.c_str())) return 1;
379  if(!strcmp("Derivative",name.c_str())) return 2;
380  if(!strcmp("KDuDx",name.c_str())) return 3;
381  if(!strcmp("KDuDy",name.c_str())) return 4;
382  if(!strcmp("KDuDz",name.c_str())) return 5;
383  if(!strcmp("NormKDu",name.c_str())) return 6;
384  if(!strcmp("MinusKGradU",name.c_str())) return 7;
385  if(!strcmp("POrder",name.c_str())) return 8;
386  if(!strcmp("Laplac",name.c_str())) return 9;
387  if(!strcmp("Stress",name.c_str())) return 10;
388  if(!strcmp("Flux",name.c_str())) return 10;
389  if(!strcmp("Pressure",name.c_str())) return 11;
390 
391  if(!strcmp("ExactSolution",name.c_str())) return 12;
392  if(!strcmp("ExactFlux",name.c_str())) return 13;
393  if(!strcmp("Divergence",name.c_str())) return 14;
394  if(!strcmp("ExactDiv",name.c_str())) return 15;
395 
396  if(!strcmp("PressureOmega1",name.c_str())) return 16;
397  if(!strcmp("PressureOmega2",name.c_str())) return 17;
398  if(!strcmp("FluxOmega1",name.c_str())) return 18;
399 
400  if(!strcmp("GradFluxX",name.c_str())) return 19;
401  if(!strcmp("GradFluxY",name.c_str())) return 20;
402  if(!strcmp("FluxL2",name.c_str())) return 21;//Only To calculate l2 error
403  if(!strcmp("Permeability",name.c_str())) return 22; // output the permeability
404  return TPZMaterial::VariableIndex(name);
405 }
406 
407 int TPZMatLaplacian::NSolutionVariables(int var){
408  if(var == 1) return 1;
409  if(var == 2) return 3;//arrumar o fluxo de hdiv para ser fdim tbem enquanto isso faco isso
410  if ((var == 3) || (var == 4) || (var == 5) || (var == 6)) return 1;
411  if (var == 7) return fDim;
412  if (var == 8) return 1;
413  if (var == 9) return 1;
414  if (var==10) return fDim;
415  if (var==11) return 1;
416 
417  if (var==12) return 1;
418  if (var==13) return fDim;
419  if (var==14) return 1;
420  if (var==15) return 1;
421  //teste de acoplamento
422  if (var==16) return 1;
423  if (var==17) return 1;
424  if (var==18) return 3;
425  if (var==19) return 3;
426  if (var==20) return 3;
427  if (var==21) return fDim;
428  if (var==22) return 1; // number of permeabilities
429 
430 
431  return TPZMaterial::NSolutionVariables(var);
432 }
433 
434 void TPZMatLaplacian::Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout){
435 
436  TPZVec<STATE> pressure(1);
437  TPZVec<REAL> pto(3);
438  TPZFMatrix<STATE> flux(3,1);
439 
440  int numbersol = data.sol.size();
441  if (numbersol != 1) {
442  DebugStop();
443  }
444  STATE perm = fK;
445  if(fPermeabilityFunction)
446  {
447  TPZManVector<STATE,3> f;
448  TPZFNMatrix<18,STATE> df(6,3);
449  fPermeabilityFunction->Execute(data.x, f, df);
450  perm = df(0,0);
451  }
452 
453 #ifndef STATE_COMPLEX
454 
455  switch (var) {
456  /* case 7:
457  {
458  // { //MinusKGradU
459  int id;
460  TPZManVector<STATE> dsolp(3,0);
461  dsolp[0] = data.dsol[0](0,0)*data.axes(0,0)+data.dsol[0](1,0)*data.axes(1,0);
462  dsolp[1] = data.dsol[0](0,0)*data.axes(0,1)+data.dsol[0](1,0)*data.axes(1,1);
463  dsolp[2] = data.dsol[0](0,0)*data.axes(0,2)+data.dsol[0](1,0)*data.axes(1,2);
464  for(id=0 ; id<fDim; id++)
465  {
466  Solout[id] = -1. * this->fK * dsolp[id];
467  }
468  }
469  break;*/
470  case 8:
471  Solout[0] = data.p;
472  break;
473  case 10:
474  if (data.numberdualfunctions) {
475 
476  Solout[0]=data.sol[0][0];
477  Solout[1]=data.sol[0][1];
478 
479  }
480  else {
481  TPZFNMatrix<3,STATE> dsolxy(3,0);
482  TPZAxesTools<STATE>::Axes2XYZ(data.dsol[0], dsolxy, data.axes);
483  for (int i=0; i<fDim; i++) {
484  Solout[i] = -perm*dsolxy(i,0);
485  }
486  }
487 
488  break;
489 
490  case 21:
491  Solout[0]=data.sol[0][0];
492  Solout[1]=data.sol[0][1];
493  break;
494 
495  case 11:
496  if (data.numberdualfunctions) {
497  Solout[0]=data.sol[0][2];
498  }
499  else{
500  Solout[0]=data.sol[0][0];
501  }
502  break;
503 
504  case 12:
505  fForcingFunctionExact->Execute(data.x,pressure,flux);
506 
507  Solout[0]=pressure[0];
508  break;
509  case 13:
510  fForcingFunctionExact->Execute(data.x,pressure,flux);
511 
512  Solout[0]=flux(0,0);
513  Solout[1]=flux(1,0);
514  break;
515 
516  case 14:
517  {
518  if (data.numberdualfunctions){
519  Solout[0]=data.sol[0][data.sol[0].NElements()-1];
520  }else{
521  Solout[0]=data.dsol[0](0,0)+data.dsol[0](1,1);
522  }
523 
524  }
525  break;
526 
527  case 15:
528  {
529  fForcingFunctionExact->Execute(data.x,pressure,flux);
530  Solout[0]=flux(2,0);
531  }
532  break;
533 
534 
535  case 16:
536  if (data.numberdualfunctions) {
537  Solout[0]=data.sol[0][2];
538  }
539  else {
540  std::cout<<"Pressao somente em Omega1"<<std::endl;
541  Solout[0]=0;//NULL;
542  }
543 
544  break;
545 
546  case 17:
547  if (!data.numberdualfunctions) {
548  Solout[0]=data.sol[0][0];
549  }
550  else {
551  std::cout<<"Pressao somente em omega2"<<std::endl;
552  Solout[0]=0;//NULL;
553  }
554 
555  break;
556  case 18:
557  if( data.numberdualfunctions){
558  Solout[0]=data.sol[0][0];//fluxo de omega1
559  Solout[1]=data.sol[0][1];
560  // Solout[2]=data.sol[2];
561  return;
562  }
563 
564  case 19:
565  if(data.numberdualfunctions){
566  Solout[0]=data.dsol[0](0,0);//fluxo de omega1
567  Solout[1]=data.dsol[0](1,0);
568  Solout[2]=data.dsol[0](2,0);
569  return;
570  }
571  case 20:
572  if( data.numberdualfunctions){
573  Solout[0]=data.dsol[0](0,1);//fluxo de omega1
574  Solout[1]=data.dsol[0](1,1);
575  Solout[2]=data.dsol[0](2,1);
576  return;
577  }
578  else {
579  std::cout<<"Pressao somente em omega2"<<std::endl;
580  Solout[0]=0;//NULL;
581  }
582  break;
583  case 22:
584  // output the permeability
585  Solout[0] = perm;
586  break;
587  default:
588  if (data.sol[0].size() == 4) {
589  data.sol[0][0] = data.sol[0][2];
590  }
591 
592  this->Solution(data.sol[0], data.dsol[0], data.axes, var, Solout);
593  break;
594  }
595 #endif
596 }
597 
598 void TPZMatLaplacian::Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,TPZVec<STATE> &Solout){
599 
600 #ifndef STATE_COMPLEX
601  Solout.Resize( this->NSolutionVariables( var ) );
602 
603  if(var == 1){
604  Solout[0] = Sol[0];//function
605  return;
606  }
607  if(var == 2) {
608  int id;
609  TPZFNMatrix<50,STATE> dsoldx;
610  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
611  for(id=0 ; id<3; id++) {
612  Solout[id] = dsoldx(id,0);//derivate
613  }
614  return;
615  }//var == 2
616  if (var == 3){ //KDuDx
617  TPZFNMatrix<9,STATE> dsoldx;
618  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
619  Solout[0] = dsoldx(0,0) * this->fK;
620  return;
621  }//var ==3
622  if (var == 4){ //KDuDy
623  TPZFNMatrix<9,STATE> dsoldx;
624  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
625  Solout[0] = dsoldx(1,0) * this->fK;
626  return;
627  }//var == 4
628  if (var == 5){ //KDuDz
629  TPZFNMatrix<9,STATE> dsoldx;
630  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
631  Solout[0] = dsoldx(2,0) * this->fK;
632  return;
633  }//var == 5
634  if (var == 6){ //NormKDu
635  int id;
636  REAL val = 0.;
637  for(id=0 ; id<fDim; id++){
638  val += (DSol(id,0) * this->fK) * (DSol(id,0) * this->fK);
639  }
640  Solout[0] = sqrt(val);
641  return;
642  }//var == 6
643  if (var == 7){ //MinusKGradU
644  int id;
645  //REAL val = 0.;
646  TPZFNMatrix<9,STATE> dsoldx;
647  TPZAxesTools<STATE>::Axes2XYZ(DSol, dsoldx, axes);
648  for(id=0 ; id<fDim; id++) {
649  Solout[id] = -1. * this->fK * dsoldx(id,0);
650  }
651  return;
652  }//var == 7
653  if(var == 9){//Laplac
654  Solout.Resize(1);
655  Solout[0] = DSol(2,0);
656  return;
657  }//Laplac
658 
659 #endif
660  TPZMaterial::Solution(Sol, DSol, axes, var, Solout);
661 
662 }//method
663 
664 void TPZMatLaplacian::Flux(TPZVec<REAL> &/*x*/, TPZVec<STATE> &/*Sol*/, TPZFMatrix<STATE> &/*DSol*/, TPZFMatrix<REAL> &/*axes*/, TPZVec<STATE> &/*flux*/) {
665  //Flux(TPZVec<REAL> &x, TPZVec<REAL> &Sol, TPZFMatrix<REAL> &DSol, TPZFMatrix<REAL> &axes, TPZVec<REAL> &flux)
666 }
667 void TPZMatLaplacian::ErrorsHdiv(TPZMaterialData &data,TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values){
668 
669  TPZVec<STATE> sol(1),dsol(fDim),div(1);
670  if(data.numberdualfunctions) Solution(data,11,sol);//pressao
671  Solution(data,21,dsol);//fluxo
672  Solution(data,14,div);//divergente
673 
674  TPZFNMatrix<18,STATE> perm(2*fDim,fDim);
675  if (fPermeabilityFunction) {
676  TPZManVector<STATE,3> dumf(3,0.);
677  fPermeabilityFunction->Execute(data.x, dumf, perm);
678  }
679  else
680  {
681  for (int i=0; i<fDim; i++) {
682  perm(i,i) = fK;
683  perm(i+fDim,i) = STATE(1.)/fK;
684  }
685  }
686 
687 
688 #ifdef LOG4CXX
689  {
690  std::stringstream sout;
691  sout<< "\n";
692  sout << " Pto " << data.x << std::endl;
693  sout<< " pressao exata " <<u_exact <<std::endl;
694  sout<< " pressao aprox " <<sol <<std::endl;
695  sout<< " ---- "<<std::endl;
696  sout<< " fluxo exato " <<du_exact<<std::endl;
697  sout<< " fluxo aprox " <<dsol<<std::endl;
698  sout<< " ---- "<<std::endl;
699  if(du_exact.Rows()>fDim) sout<< " div exato " <<du_exact(2,0)<<std::endl;
700  sout<< " div aprox " <<div<<std::endl;
701  LOGPZ_DEBUG(logger,sout.str())
702  }
703 #endif
704 
705 
706  //values[0] : pressure error using L2 norm
707  if(data.numberdualfunctions){
708  REAL diffP = abs(u_exact[0]-sol[0]);
709  values[0] = diffP*diffP;
710  }
711  //values[1] : flux error using L2 norm
712  for(int id=0; id<fDim; id++) {
713  REAL diffFlux = abs(dsol[id] - du_exact(id,0));
714  values[1] += abs(fK)*diffFlux*diffFlux;
715  }
716 
717  if(du_exact.Rows()>fDim){
718  //values[2] : divergence using L2 norm
719  REAL diffDiv = abs(div[0] - du_exact(2,0));
720  values[2]=diffDiv*diffDiv;
721  //values[3] : Hdiv norm => values[1]+values[2];
722  values[3]= values[1]+values[2];
723  }
724 
725 #ifdef LOG4CXX
726  {
727  std::stringstream sout;
728  sout << " Erro pressao " << values[0]<< std::endl;
729  sout<< "Erro fluxo " <<values[1]<<std::endl;
730  sout<< " Erro div " <<values[2] <<std::endl;
731  LOGPZ_DEBUG(logger,sout.str())
732  }
733 #endif
734 
735 
736 }
737 
738 void TPZMatLaplacian::Errors(TPZVec<REAL> &x,TPZVec<STATE> &u,
739  TPZFMatrix<STATE> &dudxaxes, TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
740  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) {
741 
742  values.Resize(3);
743  values.Fill(0.0);
744 
745  TPZFNMatrix<9,STATE> perm(2*fDim,fDim);
746  TPZManVector<STATE,3> val(fDim);
747  if (fPermeabilityFunction) {
748  fPermeabilityFunction->Execute(x, val, perm);
749  }
750  else
751  {
752  for (int i=0; i<fDim; i++) {
753  for (int j=0; j<fDim; j++)
754  {
755  perm(i,j) = this->fTensorK(i,j);
756  perm(fDim+i,j) = this->fInvK(i,j);
757  }
758  }
759  }
760  TPZFNMatrix<3,STATE> dudx(3,1,0.);
761  TPZAxesTools<STATE>::Axes2XYZ(dudxaxes, dudx, axes);
762 
764  values[1] = TPZExtractVal::val((u[0] - u_exact[0])*(u[0] - u_exact[0]));
765 
767  values[2] = 0.;
768  for(int i = 0; i < du_exact.Rows(); i++){
769  values[2] += TPZExtractVal::val( (dudx(i,0) - du_exact(i,0))*(dudx(i,0) - du_exact(i,0)));
770  }
771  // Energy Norm
772  values[0] = 0.;
773  for (int i=0; i<fDim; i++) {
774  for (int j=0; j<fDim; j++) {
775  values[0] += TPZExtractVal::val((dudx(i,0) - du_exact(i,0))*perm(i,j)*(dudx(j,0) - du_exact(j,0)));
776  }
777  }
779 
780 }
781 
782 void TPZMatLaplacian::BCInterfaceJump(TPZVec<REAL> &x, TPZSolVec &leftu,TPZBndCond &bc,TPZSolVec & jump){
783  int numbersol = leftu.size();
784  for (int is=0; is<numbersol ; is++) {
785  jump[is].Resize(1);
786  if(bc.Type() == 0){ //DIRICHLET
787  STATE f = bc.Val2()(0,0);
788  jump[is][0] = leftu[is][0] - f;
789  }
790  else{
791  jump[is].Fill(0.);
792  }
793  }
794 }//method
795 
796 
797 void TPZMatLaplacian::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright,
798  REAL weight,
799  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef){
800 
801  TPZFMatrix<REAL> &dphiLdAxes = dataleft.dphix;
802  TPZFMatrix<REAL> &dphiRdAxes = dataright.dphix;
803  TPZFMatrix<REAL> &phiL = dataleft.phi;
804  TPZFMatrix<REAL> &phiR = dataright.phi;
805  TPZManVector<REAL,3> &normal = data.normal;
806 
807  TPZFNMatrix<660> dphiL, dphiR;
808  TPZAxesTools<REAL>::Axes2XYZ(dphiLdAxes, dphiL, dataleft.axes);
809  TPZAxesTools<REAL>::Axes2XYZ(dphiRdAxes, dphiR, dataright.axes);
810 
811  int &LeftPOrder=dataleft.p;
812  int &RightPOrder=dataright.p;
813 
814  REAL &faceSize=data.HSize;
815 
816 
817  int nrowl = phiL.Rows();
818  int nrowr = phiR.Rows();
819  int il,jl,ir,jr,id;
820 
821  //diffusion term
822  STATE leftK, rightK;
823  leftK = this->fK;
824  rightK = this->fK;
825 
826  // 1) phi_I_left, phi_J_left
827  for(il=0; il<nrowl; il++) {
828  REAL dphiLinormal = 0.;
829  for(id=0; id<fDim; id++) {
830  dphiLinormal += dphiL(id,il)*normal[id];
831  }
832  for(jl=0; jl<nrowl; jl++) {
833  REAL dphiLjnormal = 0.;
834  for(id=0; id<fDim; id++) {
835  dphiLjnormal += dphiL(id,jl)*normal[id];
836  }
837  ek(il,jl) += (STATE)(weight * ( this->fSymmetry * (0.5)*dphiLinormal*phiL(jl,0)-(0.5)*dphiLjnormal*phiL(il,0))) * leftK;
838  }
839  }
840 
841  // 2) phi_I_right, phi_J_right
842  for(ir=0; ir<nrowr; ir++) {
843  REAL dphiRinormal = 0.;
844  for(id=0; id<fDim; id++) {
845  dphiRinormal += dphiR(id,ir)*normal[id];
846  }
847  for(jr=0; jr<nrowr; jr++) {
848  REAL dphiRjnormal = 0.;
849  for(id=0; id<fDim; id++) {
850  dphiRjnormal += dphiR(id,jr)*normal[id];
851  }
852  ek(ir+nrowl,jr+nrowl) += (STATE)(weight * (this->fSymmetry * ((-0.5) * dphiRinormal * phiR(jr) ) + (0.5) * dphiRjnormal * phiR(ir))) * rightK;
853  }
854  }
855 
856  // 3) phi_I_left, phi_J_right
857  for(il=0; il<nrowl; il++) {
858  REAL dphiLinormal = 0.;
859  for(id=0; id<fDim; id++) {
860  dphiLinormal += dphiL(id,il)*normal[id];
861  }
862  for(jr=0; jr<nrowr; jr++) {
863  REAL dphiRjnormal = 0.;
864  for(id=0; id<fDim; id++) {
865  dphiRjnormal += dphiR(id,jr)*normal[id];
866  }
867  ek(il,jr+nrowl) += (STATE)weight * ((STATE)fSymmetry * ((STATE)((-0.5) * dphiLinormal * phiR(jr)) * leftK ) - (STATE)((0.5) * dphiRjnormal * phiL(il))* rightK );
868  }
869  }
870 
871  // 4) phi_I_right, phi_J_left
872  for(ir=0; ir<nrowr; ir++) {
873  REAL dphiRinormal = 0.;
874  for(id=0; id<fDim; id++) {
875  dphiRinormal += dphiR(id,ir)*normal[id];
876  }
877  for(jl=0; jl<nrowl; jl++) {
878  REAL dphiLjnormal = 0.;
879  for(id=0; id<fDim; id++) {
880  dphiLjnormal += dphiL(id,jl)*normal[id];
881  }
882  ek(ir+nrowl,jl) += (STATE)weight * (
883  (STATE)(fSymmetry * (0.5) * dphiRinormal * phiL(jl)) * rightK + (STATE)((0.5) * dphiLjnormal * phiR(ir)) * leftK
884  );
885  }
886  }
887 
888  if (this->IsSymetric()){
889  if ( !ek.VerifySymmetry(1.e-10) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
890  }
891 
892  if (this->fPenaltyConstant == 0.) return;
893 
894  leftK = this->fK;
895  rightK = this->fK;
896 
897 
898 
899  //penalty = <A p^2>/h
900  REAL penalty = fPenaltyConstant * (0.5 * (abs(leftK)*LeftPOrder*LeftPOrder + abs(rightK)*RightPOrder*RightPOrder)) / faceSize;
901 
902  if (this->fPenaltyType == ESolutionPenalty || this->fPenaltyType == EBoth){
903 
904  // 1) left i / left j
905  for(il=0; il<nrowl; il++) {
906  for(jl=0; jl<nrowl; jl++) {
907  ek(il,jl) += weight * penalty * phiL(il,0) * phiL(jl,0);
908  }
909  }
910 
911  // 2) right i / right j
912  for(ir=0; ir<nrowr; ir++) {
913  for(jr=0; jr<nrowr; jr++) {
914  ek(ir+nrowl,jr+nrowl) += weight * penalty * phiR(ir,0) * phiR(jr,0);
915  }
916  }
917 
918  // 3) left i / right j
919  for(il=0; il<nrowl; il++) {
920  for(jr=0; jr<nrowr; jr++) {
921  ek(il,jr+nrowl) += -1.0 * weight * penalty * phiR(jr,0) * phiL(il,0);
922  }
923  }
924 
925  // 4) right i / left j
926  for(ir=0; ir<nrowr; ir++) {
927  for(jl=0; jl<nrowl; jl++) {
928  ek(ir+nrowl,jl) += -1.0 * weight * penalty * phiL(jl,0) * phiR(ir,0);
929  }
930  }
931 
932  }
933 
934  if (this->fPenaltyType == EFluxPenalty || this->fPenaltyType == EBoth){
935 
936  REAL NormalFlux_i = 0.;
937  REAL NormalFlux_j = 0.;
938 
939  // 1) left i / left j
940  for(il=0; il<nrowl; il++) {
941  NormalFlux_i = 0.;
942  for(id=0; id<fDim; id++) {
943  NormalFlux_i += dphiL(id,il)*normal[id];
944  }
945  for(jl=0; jl<nrowl; jl++) {
946  NormalFlux_j = 0.;
947  for(id=0; id<fDim; id++) {
948  NormalFlux_j += dphiL(id,jl)*normal[id];
949  }
950  ek(il,jl) += (STATE)(weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * leftK;
951  }
952  }
953 
954  // 2) right i / right j
955  for(ir=0; ir<nrowr; ir++) {
956  NormalFlux_i = 0.;
957  for(id=0; id<fDim; id++) {
958  NormalFlux_i += dphiR(id,ir)*normal[id];
959  }
960  for(jr=0; jr<nrowr; jr++) {
961  NormalFlux_j = 0.;
962  for(id=0; id<fDim; id++) {
963  NormalFlux_j += dphiR(id,jr)*normal[id];
964  }
965  ek(ir+nrowl,jr+nrowl) += (STATE)(weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * rightK;
966  }
967  }
968 
969  // 3) left i / right j
970  for(il=0; il<nrowl; il++) {
971  NormalFlux_i = 0.;
972  for(id=0; id<fDim; id++) {
973  NormalFlux_i += dphiL(id,il)*normal[id];
974  }
975  for(jr=0; jr<nrowr; jr++) {
976  NormalFlux_j = 0.;
977  for(id=0; id<fDim; id++) {
978  NormalFlux_j += dphiR(id,jr)*normal[id];
979  }
980  ek(il,jr+nrowl) += (STATE)((-1.) * weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * rightK;
981  }
982  }
983 
984  // 4) right i / left j
985  for(ir=0; ir<nrowr; ir++) {
986  NormalFlux_i = 0.;
987  for(id=0; id<fDim; id++) {
988  NormalFlux_i += dphiR(id,ir)*normal[id];
989  }
990  for(jl=0; jl<nrowl; jl++) {
991  NormalFlux_j = 0.;
992  for(id=0; id<fDim; id++) {
993  NormalFlux_j += dphiL(id,jl)*normal[id];
994  }
995  ek(ir+nrowl,jl) += (STATE)((-1.) * weight * ((1.)/penalty) * NormalFlux_i * NormalFlux_j) * leftK;
996  }
997  }
998 
999  }
1000 
1001 }
1002 
1004 void TPZMatLaplacian::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
1005  REAL weight,
1006  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc){
1007 
1008  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
1009  TPZFMatrix<REAL> &phiL = dataleft.phi;
1010  TPZManVector<REAL,3> &normal = data.normal;
1011  int POrder= dataleft.p;
1012  REAL faceSize=data.HSize;
1013 
1014  STATE v2[1];
1015  v2[0] = bc.Val2()(0,0);
1016 
1017  if(bc.HasForcingFunction()) { // phi(in, 0) = phi_in
1018  TPZManVector<STATE> res(1);
1019  bc.ForcingFunction()->Execute(data.x,res); // dphi(i,j) = dphi_j/dxi
1020  v2[0] = res[0];
1021  }
1022 
1023  // cout << "Material Id " << bc.Id() << " normal " << normal << "\n";
1024  int il,jl,nrowl,id;
1025  nrowl = phiL.Rows();
1026  switch(bc.Type()) {
1027  case 0: // DIRICHLET
1028 
1029  //Diffusion
1030  for(il=0; il<nrowl; il++) {
1031  REAL dphiLinormal = 0.;
1032  for(id=0; id<fDim; id++) {
1033  dphiLinormal += dphiL(id,il)*normal[id];
1034  }
1035  ef(il,0) += (STATE)(weight*dphiLinormal*fSymmetry)*fK*v2[0];
1036  for(jl=0; jl<nrowl; jl++) {
1037  REAL dphiLjnormal = 0.;
1038  for(id=0; id<fDim; id++) {
1039  dphiLjnormal += dphiL(id,jl)*normal[id];
1040  }
1041  ek(il,jl) += (STATE)(weight*(fSymmetry * dphiLinormal * phiL(jl,0) - dphiLjnormal * phiL(il,0)))*fK;
1042  }
1043  }
1044 
1045  case 1: // Neumann
1046  for(il=0; il<nrowl; il++) {
1047  ef(il,0) += (STATE)(weight*phiL(il,0))*v2[0];
1048  }
1049  break;
1050 
1051  default:
1052  PZError << __PRETTY_FUNCTION__ << " - Wrong boundary condition type\n";
1053  break;
1054  }
1055  if (this->IsSymetric()){
1056  if ( !ek.VerifySymmetry(1.e-10) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
1057  }
1058 
1059  if (this->fPenaltyConstant == 0.) return;
1060 
1061  if (this->fPenaltyType == ESolutionPenalty || this->fPenaltyType == EBoth){
1062  nrowl = phiL.Rows();
1063  const REAL penalty = fPenaltyConstant * abs(fK) * POrder * POrder / faceSize; //Ap^2/h
1064  REAL outflow = 0.;
1065 
1066  switch(bc.Type()) {
1067  case 0: // DIRICHLET
1068  for(il=0; il<nrowl; il++) {
1069  ef(il,0) += (STATE)(weight*penalty*phiL(il,0))*v2[0];
1070  for(jl=0; jl<nrowl; jl++) {
1071  ek(il,jl) += weight * penalty * phiL(il,0) * phiL(jl,0);
1072  }
1073  }
1074 
1075  break;
1076  case 1: // Neumann
1077  if(outflow > 0.)
1078  {
1079  for(il=0; il<nrowl; il++)
1080  {
1081  for(jl=0; jl<nrowl; jl++)
1082  {
1083  ek(il,jl) += weight * outflow * phiL(il,0) * phiL(jl,0);
1084  }
1085  }
1086  }
1087  //nothing to be done
1088  break;
1089  default:
1090  PZError << "TPZMatLaplacian::Wrong boundary condition type\n";
1091  break;
1092  }
1093 
1094  }
1095 
1096 }
1097 
1098 
1099 void TPZMatLaplacian::FillDataRequirements(TPZMaterialData &data)
1100 {
1101  data.SetAllRequirements(false);
1102  data.fNeedsNeighborSol = false;
1103  data.fNeedsNeighborCenter = false;
1104  data.fNeedsNormal = true;
1105  data.fNeedsSol = true;
1106 }
1107 
1108 void TPZMatLaplacian::FillDataRequirementsInterface(TPZMaterialData &data)
1109 {
1110  data.SetAllRequirements(false);
1111  data.fNeedsNeighborSol = false;
1112  data.fNeedsNeighborCenter = false;
1113  data.fNeedsNormal = true;
1114  data.fNeedsHSize = true;
1115 }
1116 
1117 REAL TPZMatLaplacian::ComputeSquareResidual(TPZVec<REAL>& X, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol){
1118  // residual = -fK Laplac(u) - (fXf)
1119  STATE fXfLoc = fXf;
1120  if(fForcingFunction) {
1121  TPZManVector<STATE> res(1);
1122  fForcingFunction->Execute(X,res);
1123  fXfLoc = res[0];
1124  }
1125 
1126  STATE laplacU = (STATE)0;
1127  if(this->Dimension() == 1){
1128  laplacU = dsol(1,0);
1129  }
1130  if(this->Dimension() == 2){
1131  laplacU = dsol(2,0);
1132  }
1133 
1134  REAL result = abs(-this->fK * laplacU - (fXfLoc));
1135  return (result*result);
1136 }
1137 
1138 void TPZMatLaplacian::Write(TPZStream &buf, int withclassid) const {
1139  TPZDiscontinuousGalerkin::Write(buf, withclassid);
1140  buf.Write(&fXf, 1);
1141  buf.Write(&fDim, 1);
1142  buf.Write(&fK, 1);
1143  buf.Write(&fSymmetry, 1);
1144  buf.Write(&fPenaltyConstant,1);
1145 }
1146 
1147 void TPZMatLaplacian::Read(TPZStream &buf, void *context) {
1148  TPZDiscontinuousGalerkin::Read(buf, context);
1149  buf.Read(&fXf, 1);
1150  buf.Read(&fDim, 1);
1151  buf.Read(&fK, 1);
1152  buf.Read(&fSymmetry, 1);
1153  buf.Read(&fPenaltyConstant,1);
1154 }
1155 
1156 int TPZMatLaplacian::ClassId() const{
1157  return Hash("TPZMatLaplacian") ^ TPZDiscontinuousGalerkin::ClassId() << 1;
1158 }
1159 
1160 #ifndef BORLAND
1161 template class TPZRestoreClass<TPZMatLaplacian>;
1162 #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
TPZMatLaplacian::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: TPZMatLaplacian.cpp:377
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
TPZMatLaplacian::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
Definition: TPZMatLaplacian.cpp:598
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
TPZMatLaplacian::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: TPZMatLaplacian.cpp:664
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZDiscontinuousGalerkin::Read
void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzdiscgal.cpp:118
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
TPZMatLaplacian::fXf
STATE fXf
Forcing function value.
Definition: TPZMatLaplacian.h:25
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZMatLaplacian::ContributeBCInterface
virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Definition: TPZMatLaplacian.cpp:1004
TPZMaterial::operator=
TPZMaterial & operator=(const TPZMaterial &copy)
operator =
Definition: TPZMaterial.cpp:62
TPZMatLaplacian::ErrorsHdiv
void ErrorsHdiv(TPZMaterialData &data, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &values) override
Definition: TPZMatLaplacian.cpp:667
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.
TPZMaterialData::SetAllRequirements
void SetAllRequirements(bool set)
Set all flags at once.
Definition: pzmaterialdata.cpp:90
TPZMatLaplacian::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
It computes a contribution to the stiffness matrix and load vector at one integration point...
Definition: TPZMatLaplacian.cpp:101
TPZMatLaplacian::fPenaltyConstant
REAL fPenaltyConstant
Constant multiplier of penalty term, when required is set.
Definition: TPZMatLaplacian.h:55
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZMaterialData::fShapeType
MShapeFunctionType fShapeType
Definition: pzmaterialdata.h:43
abs
TinyFad< 8, T > abs(const TinyFad< 8, T > &in)
Definition: tinyfadeight.h:846
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZMatLaplacian::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: TPZMatLaplacian.cpp:738
TPZMatLaplacian::fDim
int fDim
Problem dimension.
Definition: TPZMatLaplacian.h:28
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
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
TPZMatLaplacian::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: TPZMatLaplacian.cpp:1156
TPZMaterialData::p
int p
maximum polinomial order of the shape functions
Definition: pzmaterialdata.h:75
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
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
TPZMatLaplacian::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: TPZMatLaplacian.cpp:1138
test.f
f
Definition: test.py:287
TPZMatLaplacian::Read
virtual void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: TPZMatLaplacian.cpp:1147
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZDiscontinuousGalerkin::Write
void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzdiscgal.cpp:114
TPZMatLaplacian::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: TPZMatLaplacian.cpp:1117
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
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
TPZMatLaplacian::SetNonSymmetric
void SetNonSymmetric()
Set material elliptic term as the Baumann&#39;s formulation, i.e. the non-symmetrical formulation...
Definition: TPZMatLaplacian.h:94
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZMatLaplacian::fPenaltyType
EPenaltyType fPenaltyType
Penalty term definition.
Definition: TPZMatLaplacian.h:47
TPZMatLaplacian::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: TPZMatLaplacian.cpp:87
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
TPZMatLaplacian::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data) override
Fill material data parameter with necessary requirements for the.
Definition: TPZMatLaplacian.cpp:1099
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZMatLaplacian::Name
virtual std::string Name() override
Returns the name of the material.
Definition: TPZMatLaplacian.h:189
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
TPZExtractVal::val
static REAL val(const int number)
Definition: pzextractval.h:21
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.
TPZMatLaplacian::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: TPZMatLaplacian.cpp:91
TPZMatLaplacian::BCInterfaceJump
virtual void BCInterfaceJump(TPZVec< REAL > &x, TPZSolVec &leftu, TPZBndCond &bc, TPZSolVec &jump) override
Computes interface jump from element to Dirichlet boundary condition.
Definition: TPZMatLaplacian.cpp:782
TPZMatLaplacian::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, 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...
Definition: TPZMatLaplacian.cpp:321
TPZMatLaplacian::operator=
TPZMatLaplacian & operator=(const TPZMatLaplacian &copy)
Definition: TPZMatLaplacian.cpp:54
TPZMatLaplacian::~TPZMatLaplacian
virtual ~TPZMatLaplacian()
Definition: TPZMatLaplacian.cpp:84
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
TPZMatLaplacian::fInvK
TPZFNMatrix< 9, STATE > fInvK
Definition: TPZMatLaplacian.h:34
TPZMatLaplacian::fK
STATE fK
Coeficient which multiplies the Laplacian operator.
Definition: TPZMatLaplacian.h:31
TPZMatLaplacian::SetNoPenalty
void SetNoPenalty()
Defines no penalty terms in ContributeInterface.
Definition: TPZMatLaplacian.h:58
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
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
TPZMatLaplacian::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: TPZMatLaplacian.cpp:797
pzaxestools.h
Contains declaration of the TPZAxesTools class which implements verifications over axes...
TPZMatLaplacian::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMatLaplacian.cpp:407
TPZMatLaplacian::fPermeabilityFunction
TPZAutoPointer< TPZFunction< STATE > > fPermeabilityFunction
Pointer to forcing function, it is the Permeability and its inverse.
Definition: TPZMatLaplacian.h:50
TPZMatLaplacian::Dimension
virtual int Dimension() const override
Returns the integrable dimension of the material.
Definition: TPZMatLaplacian.h:134
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZMatLaplacian.h
Contains the TPZMatLaplacian class.
TPZMatLaplacian::fSymmetry
REAL fSymmetry
Symmetry coefficient of elliptic term.
Definition: TPZMatLaplacian.h:41
rdt.values
def values
Definition: rdt.py:119
TPZMatLaplacian::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: TPZMatLaplacian.cpp:193
TPZMaterialData::fVecShapeIndex
TPZManVector< std::pair< int, int64_t > > fVecShapeIndex
correspondence between normal vector index and index of the shape functions
Definition: pzmaterialdata.h:99
TPZMatLaplacian::SetParameters
void SetParameters(STATE diff, STATE f)
Set a uniform diffusion constant and external flux.
Definition: TPZMatLaplacian.cpp:68
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZMatLaplacian::ContributeBCHDiv
virtual void ContributeBCHDiv(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: TPZMatLaplacian.cpp:260
TPZMatLaplacian::fTensorK
TPZFNMatrix< 9, STATE > fTensorK
Tensor de permeabilidade.
Definition: TPZMatLaplacian.h:34
TPZMatLaplacian::FillDataRequirementsInterface
virtual void FillDataRequirementsInterface(TPZMaterialData &data) override
Fill material data parameter with necessary requirements for the ContributeInterface method...
Definition: TPZMatLaplacian.cpp:1108
TPZRestoreClass
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
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
Generated on Sun Aug 16 2020 11:17:04 for NeoPZ by   doxygen 1.8.13