NeoPZ
pznonlinearpoisson3d.cpp
Go to the documentation of this file.
1 
6 #include "pznonlinearpoisson3d.h"
7 #include "pzbndcond.h"
8 
9 using namespace std;
10 
11 TPZNonLinearPoisson3d::TPZNonLinearPoisson3d(int nummat, int dim):
12 TPZRegisterClassId(&TPZNonLinearPoisson3d::ClassId),
13 TPZMatPoisson3dReferred(nummat, dim){
14  this->fIsReferred = true;
15  this->SetNoStabilizationTerm();
16 }
17 
18 TPZNonLinearPoisson3d::TPZNonLinearPoisson3d(const TPZNonLinearPoisson3d &cp):
19 TPZRegisterClassId(&TPZNonLinearPoisson3d::ClassId),
20 TPZMatPoisson3dReferred(cp){
21  this->fIsReferred = cp.fIsReferred;
22  this->fStabilizationType = cp.fStabilizationType;
23 }
24 
25 TPZNonLinearPoisson3d::~TPZNonLinearPoisson3d(){
26 
27 }
28 
29 void TPZNonLinearPoisson3d::SetSUPGStab(STATE sd){
30  this->fStabilizationType = ESUPG;
31  this->fSD = sd;
32 }
33 
34 void TPZNonLinearPoisson3d::SetGradientStab(STATE sd){
35  this->fStabilizationType = EGradient;
36  this->fSD = sd;
37 }
38 
39 void TPZNonLinearPoisson3d::SetNoStabilizationTerm(){
40  this->fStabilizationType = ENoStabilization;
41  this->fSD = 0.;
42 }
43 
44 void TPZNonLinearPoisson3d::Contribute(TPZMaterialData &data,
45  REAL weight,
46  TPZFMatrix<STATE> &ek,
47  TPZFMatrix<STATE> &ef){
48 
49  TPZFMatrix<REAL> &dphi = data.dphix;
50  TPZFMatrix<REAL> &phi = data.phi;
51  TPZManVector<REAL,3> &x = data.x;
52  int numbersol = data.sol.size();
53  if (numbersol != 1) {
54  DebugStop();
55  }
56 
57  TPZVec<STATE> &sol=data.sol[0];
58  TPZFMatrix<STATE> &dsol=data.dsol[0];
59  TPZFMatrix<REAL> &axes=data.axes;
60  TPZFMatrix<REAL> &jacinv=data.jacinv;
61 
62  if (this->IsReferred()){
63  this->SetConvectionTerm(dsol, axes);
64  }
65 
66  int phr = phi.Rows();
67 
68  if(fForcingFunction) {
69  TPZManVector<STATE> res(1);
70  fForcingFunction->Execute(x,res);
71  fXf = res[0];
72  }
73  STATE delx = 0.;
74  STATE ConvDirAx[3] = {0.};
75  if(fC != 0.0) {
76  int di,dj;
77  delx = 0.;
78  for(di=0; di<fDim; di++) {
79  for(dj=0; dj<fDim; dj++) {
80  delx = (delx<fabs(jacinv(di,dj))) ? fabs(jacinv(di,dj)) : delx;
81  }
82  }
83  delx = 2./delx;
84 
85 
86  switch(fDim) {
87  case 1:
88  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
89  break;
90  case 2:
91  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
92  ConvDirAx[1] = axes(1,0)*fConvDir[0]+axes(1,1)*fConvDir[1]+axes(1,2)*fConvDir[2];
93  break;
94  case 3:
95  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
96  ConvDirAx[1] = axes(1,0)*fConvDir[0]+axes(1,1)*fConvDir[1]+axes(1,2)*fConvDir[2];
97  ConvDirAx[2] = axes(2,0)*fConvDir[0]+axes(2,1)*fConvDir[1]+axes(2,2)*fConvDir[2];
98  break;
99  default:
100  PZError << "TPZNonLinearPoisson3d::Contribute dimension error " << fDim << endl;
101  }
102  }
103 
104  for( int in = 0; in < phr; in++ ) {
105  int kd;
106  STATE dphiic = 0;
107  for(kd = 0; kd<fDim; kd++) dphiic += ConvDirAx[kd]*dphi(kd,in);
108  ef(in, 0) += - weight * ( fXf*phi(in,0) + 0.5*fSD*delx*fC*dphiic*fXf );
109  for(kd = 0; kd < fDim; kd++){
110  ef(in, 0) += -1. * weight * ( +fK * ( dphi(kd,in) * dsol(kd,0) )
111  -fC * ( ConvDirAx[kd]* dphi(kd,in) * sol[0] ) );
112  }//kd
113 
114  for( int jn = 0; jn < phr; jn++ ) {
115  for(kd=0; kd<fDim; kd++) {
116  ek(in,jn) += weight * (
117  +fK * ( dphi(kd,in) * dphi(kd,jn) )
118  -fC * ( ConvDirAx[kd]* dphi(kd,in) * phi(jn) ) );
119  }
120  }
121  }//in
122 
123  if (fStabilizationType == ESUPG){
124 
125  for( int in = 0; in < phr; in++ ) {
126  int kd;
127  STATE dphiic = 0;
128  for(kd = 0; kd<fDim; kd++) dphiic += ConvDirAx[kd]*dphi(kd,in);
129  ef(in, 0) += - weight * ( + 0.5*fSD*delx*fC*dphiic*fXf );
130  for(kd = 0; kd < fDim; kd++){
131  ef(in, 0) += -1. * weight * ( +0.5 * fSD * delx * fC * dphiic * dsol(kd,0) * ConvDirAx[kd] );
132  }//kd
133 
134  for( int jn = 0; jn < phr; jn++ ) {
135  for(kd=0; kd<fDim; kd++) {
136  ek(in,jn) += weight * (
137  +0.5 * fSD * delx * fC * dphiic * dphi(kd,jn)* ConvDirAx[kd]
138  );
139  }
140  }
141  }//in
142 
143  }//SUPG
144 
145  if (fStabilizationType == EGradient){
146 
147  //computing norm of solution gradient
148  STATE dsolNorm = 0.;
149  for(int d = 0; d < fDim; d++) dsolNorm += dsol(d,0)*dsol(d,0);
150  dsolNorm = sqrt(dsolNorm);
151  if (dsolNorm < 1e-16) dsolNorm = 1.;
152 
153  //loop over i shape functions
154  int kd;
155  for( int in = 0; in < phr; in++ ){
156 
157  //computing gradV.gradU/Norm(gradU)
158  STATE dphiic = 0.;
159  for(kd = 0; kd<fDim; kd++) dphiic += dsol(kd,0) * dphi(kd,in) / dsolNorm;
160 
161  ef(in, 0) += - weight * ( 0.5*fSD*delx* /*dsolNorm**/ dphiic* fXf );
162 
163  double aux = 0.;
164  for(kd = 0; kd < fDim; kd++){
165  aux += dphiic * dsol(kd,0)*(fC*ConvDirAx[kd]);
166  }//kd
167  ef(in,0) += -1.* ( +0.5 * fSD * delx * aux * weight );
168 
169  for( int jn = 0; jn < phr; jn++ ) {
170  STATE DdphiicDalpha = 0.;
171  for(kd = 0; kd<fDim; kd++) DdphiicDalpha += dphi(kd,jn)*dphi(kd,in)/dsolNorm;
172  double aux = 0.;
173  for(kd=0; kd<fDim; kd++) {
174  aux += (fC*ConvDirAx[kd]) * ( dphiic*dphi(kd,jn) + dsol(kd,0)*DdphiicDalpha );
175  }//kd
176  ek(in,jn) += +0.5 * fSD * delx * aux * weight;
177  }//jn
178 
179  }//in
180 
181  }//EGradiente
182 
183  if (this->fC == 0.){
184  if ( !ek.VerifySymmetry() ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
185  }
186 
187 }
188 
189 void TPZNonLinearPoisson3d::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ef) {
190  TPZFMatrix<REAL> &dphi = data.dphix;
191  TPZFMatrix<REAL> &phi = data.phi;
192  TPZManVector<REAL,3> &x = data.x;
193  int numbersol = data.sol.size();
194  if (numbersol != 1) {
195  DebugStop();
196  }
197 
198  TPZVec<STATE> &sol=data.sol[0];
199  TPZFMatrix<STATE> &dsol=data.dsol[0];
200  TPZFMatrix<REAL> &axes=data.axes;
201  TPZFMatrix<REAL> &jacinv=data.jacinv;
202 
203  if (this->IsReferred()){
204  this->SetConvectionTerm(dsol, axes);
205  }
206 
207  int phr = phi.Rows();
208 
209  if(fForcingFunction) {
210  TPZManVector<STATE> res(1);
211  fForcingFunction->Execute(x,res);
212  fXf = res[0];
213  }
214  STATE delx = 0.;
215  STATE ConvDirAx[3] = {0.};
216  if(fC != 0.0) {
217  int di,dj;
218  delx = 0.;
219  for(di=0; di<fDim; di++) {
220  for(dj=0; dj<fDim; dj++) {
221  delx = (delx<fabs(jacinv(di,dj))) ? fabs(jacinv(di,dj)) : delx;
222  }
223  }
224  delx = 2./delx;
225 
226 
227  switch(fDim) {
228  case 1:
229  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
230  break;
231  case 2:
232  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
233  ConvDirAx[1] = axes(1,0)*fConvDir[0]+axes(1,1)*fConvDir[1]+axes(1,2)*fConvDir[2];
234  break;
235  case 3:
236  ConvDirAx[0] = axes(0,0)*fConvDir[0]+axes(0,1)*fConvDir[1]+axes(0,2)*fConvDir[2];
237  ConvDirAx[1] = axes(1,0)*fConvDir[0]+axes(1,1)*fConvDir[1]+axes(1,2)*fConvDir[2];
238  ConvDirAx[2] = axes(2,0)*fConvDir[0]+axes(2,1)*fConvDir[1]+axes(2,2)*fConvDir[2];
239  break;
240  default:
241  PZError << "TPZNonLinearPoisson3d::Contribute dimension error " << fDim << endl;
242  }
243  }
244 
245  for( int in = 0; in < phr; in++ ) {
246  int kd;
247  STATE dphiic = 0;
248  for(kd = 0; kd<fDim; kd++) dphiic += ConvDirAx[kd]*dphi(kd,in);
249  ef(in, 0) += - weight * ( fXf*phi(in,0) + 0.5*fSD*delx*fC*dphiic*fXf );
250  for(kd = 0; kd < fDim; kd++){
251  ef(in, 0) += -1. * weight * ( +fK * ( dphi(kd,in) * dsol(kd,0) )
252  -fC * ( ConvDirAx[kd]* dphi(kd,in) * sol[0] ) );
253  }//kd
254 
255  }//in
256 
257  if (fStabilizationType == ESUPG){
258 
259  for( int in = 0; in < phr; in++ ) {
260  int kd;
261  STATE dphiic = 0;
262  for(kd = 0; kd<fDim; kd++) dphiic += ConvDirAx[kd]*dphi(kd,in);
263  ef(in, 0) += - weight * ( + 0.5*fSD*delx*fC*dphiic*fXf );
264  for(kd = 0; kd < fDim; kd++){
265  ef(in, 0) += -1. * weight * ( +0.5 * fSD * delx * fC * dphiic * dsol(kd,0) * ConvDirAx[kd] );
266  }//kd
267 
268  }//in
269 
270  }//SUPG
271 
272  if (fStabilizationType == EGradient){
273 
274  //computing norm of solution gradient
275  STATE dsolNorm = 0.;
276  for(int d = 0; d < fDim; d++) dsolNorm += dsol(d,0)*dsol(d,0);
277  dsolNorm = sqrt(dsolNorm);
278  if (dsolNorm < 1e-16) dsolNorm = 1.;
279 
280  //loop over i shape functions
281  int kd;
282  for( int in = 0; in < phr; in++ ){
283 
284  //computing gradV.gradU/Norm(gradU)
285  STATE dphiic = 0.;
286  for(kd = 0; kd<fDim; kd++) dphiic += dsol(kd,0) * dphi(kd,in) / dsolNorm;
287 
288  ef(in, 0) += - weight * ( 0.5*fSD*delx* /*dsolNorm**/ dphiic* fXf );
289 
290  double aux = 0.;
291  for(kd = 0; kd < fDim; kd++){
292  aux += dphiic * dsol(kd,0)*(fC*ConvDirAx[kd]);
293  }//kd
294  ef(in,0) += -1.* ( +0.5 * fSD * delx * aux * weight );
295 
296  }//in
297 
298  }//EGradiente
299 
300 }//void
301 
302 void TPZNonLinearPoisson3d::ContributeBC(TPZMaterialData &data,
303  REAL weight,
304  TPZFMatrix<STATE> &ek,
305  TPZFMatrix<STATE> &ef,
306  TPZBndCond &bc){
307  TPZFMatrix<REAL> &phi = data.phi;
308  int numbersol = data.sol.size();
309  if (numbersol != 1) {
310  DebugStop();
311  }
312 
313  TPZVec<STATE> &sol=data.sol[0];
314  TPZFMatrix<REAL> &axes=data.axes;
315 
316  int phr = phi.Rows();
317  short in,jn;
318  STATE v2[1];
319  v2[0] = bc.Val2()(0,0);
320 
321  switch (bc.Type()) {
322  case 0 : { // Dirichlet condition
323  for(in = 0 ; in < phr; in++) {
324  ef(in,0) += weight * ( gBigNumber * v2[0] * phi(in,0) - gBigNumber * phi(in,0) * sol[0] );
325  for (jn = 0 ; jn < phr; jn++) {
326  ek(in,jn) += gBigNumber * phi(in,0) * phi(jn,0) * weight;
327  }
328  }
329  }
330  break;
331 
332  case 1 : { // Neumann condition
333  for(in = 0 ; in < phi.Rows(); in++) {
334  ef(in,0) += v2[0] * phi(in,0) * weight;
335  }
336  }
337  break;
338 
339  case 2 :{ // condicao mista
340  cout << __PRETTY_FUNCTION__ << " at line " << __LINE__ << " not implemented\n";
341  }
342  break;
343 
344  case 3 : { // outflow condition
345 
346  if (this->IsReferred()){
347  PZError << "Error at " << __PRETTY_FUNCTION__
348  << " - the outflow boundary condition can not be implemented for referred elements derived from TPZInterpolatedElement\n";
349  }
350 
351  int id, il, jl;
352  REAL normal[3];
353  if (fDim == 1) PZError << __PRETTY_FUNCTION__ << " - ERROR! The normal vector is not available for 1D TPZInterpolatedElement\n";
354  if (fDim == 2){
355  normal[0] = axes(0,1);
356  normal[1] = axes(1,1);
357  }
358  if (fDim == 3){
359  normal[0] = axes(0,2);
360  normal[1] = axes(1,2);
361  normal[2] = axes(2,2);
362  }
363  STATE ConvNormal = 0.;
364  for(id=0; id<fDim; id++) ConvNormal += fC*fConvDir[id]*normal[id];
365  if(ConvNormal > 0.) {
366  for(il=0; il<phr; il++) {
367  for(jl=0; jl<phr; jl++) {
368  ek(il,jl) += weight * ConvNormal * phi(il)*phi(jl);
369  }
370  ef(il,0) += -1. * weight * ConvNormal * phi(il) * sol[0];
371  }
372  }
373  else{
374  if (ConvNormal < 0.) std::cout << "Boundary condition error: inflow detected in outflow boundary condition: ConvNormal = " << ConvNormal << "\n";
375  }
376  }
377  break;
378 
379  default :{
380  PZError << __PRETTY_FUNCTION__ << " at line " << __LINE__ << " - Error! Wrong boundary condition type\n";
381  }
382  break;
383  }
384 
385  if (this->IsSymetric()) {//only 1.e-3 because of bignumbers.
386  if ( !ek.VerifySymmetry( 1.e-3 ) ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
387  }
388 }
389 
390 void TPZNonLinearPoisson3d::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright,
391  REAL weight,
392  TPZFMatrix<STATE> &ek,
393  TPZFMatrix<STATE> &ef){
394  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
395  TPZFMatrix<REAL> &dphiR = dataright.dphix;
396  TPZFMatrix<REAL> &phiL = dataleft.phi;
397  TPZFMatrix<REAL> &phiR = dataright.phi;
398  TPZManVector<REAL,3> &normal = data.normal;
399  int numbersol = dataleft.sol.size();
400  if (numbersol != 1) {
401  DebugStop();
402  }
403 
404  TPZVec<STATE> &solL=dataleft.sol[0];
405  TPZVec<STATE> &solR=dataright.sol[0];
406  TPZFMatrix<STATE> &dsolL=dataleft.dsol[0];
407  TPZFMatrix<STATE> &dsolR=dataright.dsol[0];
408 
409  if (this->IsReferred()){
410  this->SetConvectionTermInterface(dsolL, dsolR);
411  }
412 
413  const int nrowl = phiL.Rows();
414  const int nrowr = phiR.Rows();
415 
416  //Convection term
417  STATE ConvNormal = 0.;
418  for(int id=0; id<fDim; id++) ConvNormal += fC * fConvDir[id]*normal[id];
419  if(ConvNormal > 0.) {
420  for(int il=0; il<nrowl; il++) {
421  ef(il, 0) += -1. * weight * ConvNormal * phiL(il) * solL[0];
422  for(int jl=0; jl<nrowl; jl++) {
423  ek(il,jl) += weight * ConvNormal * phiL(il)*phiL(jl);
424  }
425  }
426  for(int ir=0; ir<nrowr; ir++) {
427  ef(ir+nrowl,0) += -1. * (-1. * weight * ConvNormal * phiR(ir) * solL[0] );
428  for(int jl=0; jl<nrowl; jl++) {
429  ek(ir+nrowl,jl) -= weight * ConvNormal * phiR(ir) * phiL(jl);
430  }
431  }
432  } else {
433  for(int ir=0; ir<nrowr; ir++) {
434  ef(ir+nrowl,0) += -1. * (-1. * weight * ConvNormal * phiR(ir) * solR[0] );
435  for(int jr=0; jr<nrowr; jr++) {
436  ek(ir+nrowl,jr+nrowl) -= weight * ConvNormal * phiR(ir) * phiR(jr);
437  }
438  }
439  for(int il=0; il<nrowl; il++) {
440  ef(il,0) += -1. * weight * ConvNormal * phiL(il) * solR[0];
441  for(int jr=0; jr<nrowr; jr++) {
442  ek(il,jr+nrowl) += weight * ConvNormal * phiL(il) * phiR(jr);
443  }
444  }
445  }
446 
447 
448  //diffusion term
449  STATE leftK, rightK;
450  leftK = this->fK;
451  rightK = this->fK;
452 
453  //Compute GradSol . normal
454  STATE DSolLNormal = 0.;
455  STATE DSolRNormal = 0.;
456  for(int id=0; id<fDim; id++) {
457  DSolLNormal += dsolL(id,0)*normal[id];
458  DSolRNormal += dsolR(id,0)*normal[id];
459  }//for
460 
461  // 1) phi_I_left, phi_J_left
462  for(int il=0; il<nrowl; il++) {
463  STATE dphiLinormal = 0.;
464  for(int id=0; id<fDim; id++) {
465  dphiLinormal += dphiL(id,il)*normal[id];
466  }
467 
468  //ef = F - K u
469  ef(il,0) += -1. * (weight * leftK * (this->fSymmetry * 0.5 * dphiLinormal*solL[0]-0.5*DSolLNormal*phiL(il,0)));
470 
471  for(int jl=0; jl<nrowl; jl++) {
472  STATE dphiLjnormal = 0.;
473  for(int id=0; id<fDim; id++) {
474  dphiLjnormal += dphiL(id,jl)*normal[id];
475  }
476  ek(il,jl) += weight * leftK * (
477  this->fSymmetry * 0.5*dphiLinormal*phiL(jl,0)-0.5*dphiLjnormal*phiL(il,0)
478  );
479  }
480  }
481 
482  // 2) phi_I_right, phi_J_right
483  for(int ir=0; ir<nrowr; ir++) {
484  STATE dphiRinormal = 0.;
485  for(int id=0; id<fDim; id++) {
486  dphiRinormal += dphiR(id,ir)*normal[id];
487  }
488 
489  //ef = F - K u
490  ef(ir+nrowl,0) += -1. * weight * rightK * ( this->fSymmetry * (-0.5 * dphiRinormal * solR[0] ) + 0.5 * DSolRNormal * phiR(ir) );
491 
492  for(int jr=0; jr<nrowr; jr++) {
493  STATE dphiRjnormal = 0.;
494  for(int id=0; id<fDim; id++) {
495  dphiRjnormal += dphiR(id,jr)*normal[id];
496  }
497  ek(ir+nrowl,jr+nrowl) += weight * rightK * (
498  this->fSymmetry * (-0.5 * dphiRinormal * phiR(jr) ) + 0.5 * dphiRjnormal * phiR(ir)
499  );
500  }
501  }
502 
503  // 3) phi_I_left, phi_J_right
504  for(int il=0; il<nrowl; il++) {
505  STATE dphiLinormal = 0.;
506  for(int id=0; id<fDim; id++) {
507  dphiLinormal += dphiL(id,il)*normal[id];
508  }
509 
510  //ef = F - K u
511  ef(il,0) += -1. * weight * ( this->fSymmetry * (-0.5 * dphiLinormal * leftK * solR[0] ) - 0.5 * DSolRNormal * rightK * phiL(il) );
512 
513  for(int jr=0; jr<nrowr; jr++) {
514  STATE dphiRjnormal = 0.;
515  for(int id=0; id<fDim; id++) {
516  dphiRjnormal += dphiR(id,jr)*normal[id];
517  }
518  ek(il,jr+nrowl) += weight * (
519  this->fSymmetry * (-0.5 * dphiLinormal * leftK * phiR(jr) ) - 0.5 * dphiRjnormal * rightK * phiL(il)
520  );
521  }
522  }
523 
524  // 4) phi_I_right, phi_J_left
525  for(int ir=0; ir<nrowr; ir++) {
526  STATE dphiRinormal = 0.;
527  for(int id=0; id<fDim; id++) {
528  dphiRinormal += dphiR(id,ir)*normal[id];
529  }
530 
531  //ef = F - K u
532  ef(ir+nrowl,0) += -1. * weight * (this->fSymmetry * 0.5 * dphiRinormal * rightK * solL[0] + 0.5 * DSolLNormal * leftK * phiR(ir));
533 
534  for(int jl=0; jl<nrowl; jl++) {
535  STATE dphiLjnormal = 0.;
536  for(int id=0; id<fDim; id++) {
537  dphiLjnormal += dphiL(id,jl)*normal[id];
538  }
539  ek(ir+nrowl,jl) += weight * (
540  this->fSymmetry * 0.5 * dphiRinormal * rightK * phiL(jl) + 0.5 * dphiLjnormal * leftK * phiR(ir)
541  );
542  }
543  }
544 
545  if (this->IsSymetric()){
546  if ( !ek.VerifySymmetry() ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
547  }
548 }
549 
550 void TPZNonLinearPoisson3d::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
551  REAL weight,
552  TPZFMatrix<STATE> &ek,
553  TPZFMatrix<STATE> &ef,
554  TPZBndCond &bc) {
555  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
556  TPZFMatrix<REAL> &phiL = dataleft.phi;
557  TPZManVector<REAL,3> &normal = data.normal;
558  int numbersol = dataleft.sol.size();
559  if (numbersol != 1) {
560  DebugStop();
561  }
562 
563  TPZVec<STATE> &solL=dataleft.sol[0];
564  TPZFMatrix<STATE> &dsolL=dataleft.dsol[0];
565 
566  if (this->IsReferred()){
567  this->SetConvectionTermInterface(dsolL, dsolL);
568  }
569 
570  int il,jl,nrowl,id;
571  nrowl = phiL.Rows();
572  STATE ConvNormal = 0.;
573  for(id=0; id<fDim; id++) ConvNormal += fC * fConvDir[id]*normal[id];
574 
575  //Compute GradSol . normal
576  STATE DSolLNormal = 0.;
577  for(id=0; id<fDim; id++) {
578  DSolLNormal += dsolL(id,0)*normal[id];
579  }//for
580 
581  switch(bc.Type()) {
582  case 0: // DIRICHLET
583 
584  //Diffusion
585  for(il=0; il<nrowl; il++) {
586  STATE dphiLinormal = 0.;
587  for(id=0; id<fDim; id++) {
588  dphiLinormal += dphiL(id,il)*normal[id];
589  }
590  ef(il,0) += weight*fK*dphiLinormal*bc.Val2()(0,0) * this->fSymmetry;
591 
592  //ef = F - K u
593  ef(il,0) += -1. * weight*fK*(this->fSymmetry * dphiLinormal * solL[0] - DSolLNormal * phiL(il,0));
594 
595  for(jl=0; jl<nrowl; jl++) {
596  STATE dphiLjnormal = 0.;
597  for(id=0; id<fDim; id++) {
598  dphiLjnormal += dphiL(id,jl)*normal[id];
599  }
600  ek(il,jl) += weight*fK*(this->fSymmetry * dphiLinormal * phiL(jl,0) - dphiLjnormal * phiL(il,0));
601  }
602  }
603 
604  //Convection
605  if(ConvNormal > 0.) {
606  for(il=0; il<nrowl; il++) {
607  ef(il,0) += -1. * weight * ConvNormal * phiL(il) * solL[0];
608  for(jl=0; jl<nrowl; jl++) {
609  ek(il,jl) += weight * ConvNormal * phiL(il)*phiL(jl);
610  }
611  }
612  } else {
613  for(il=0; il<nrowl; il++) {
614  ef(il,0) -= weight * ConvNormal * bc.Val2()(0,0) * phiL(il);
615  }
616  }
617 
618  break;
619  case 1: // Neumann
620  for(il=0; il<nrowl; il++) {
621  ef(il,0) += weight*phiL(il,0)*bc.Val2()(0,0);
622  }
623  break;
624 
625  case 3: // outflow condition
626  if(ConvNormal > 0.) {
627  for(il=0; il<nrowl; il++) {
628  for(jl=0; jl<nrowl; jl++) {
629  ek(il,jl) += weight * ConvNormal * phiL(il)*phiL(jl);
630  }
631  ef(il,0) += -1. * weight * ConvNormal * phiL(il) * solL[0];
632  }
633  }
634  else {
635  if (ConvNormal < 0.){
636  std::cout << "Boundary condition error: inflow detected in outflow boundary condition: ConvNormal = " << ConvNormal << "\n";
637  }
638  }
639  break;
640 
641  default:
642  PZError << __PRETTY_FUNCTION__ << " - Wrong boundary condition type\n";
643  break;
644  }
645  if (this->IsSymetric()){
646  if ( !ek.VerifySymmetry() ) cout << __PRETTY_FUNCTION__ << "\nMATRIZ NAO SIMETRICA" << endl;
647  }
648 }
649 
650 int TPZNonLinearPoisson3d::ClassId() const{
651  return Hash("TPZNonLinearPoisson3d") ^ TPZMatPoisson3dReferred::ClassId() << 1;
652 }
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
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
pznonlinearpoisson3d.h
Contains the TPZNonLinearPoisson3d class.
TPZMaterialData::normal
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
Definition: pzmaterialdata.h:69
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZNonLinearPoisson3d::SetGradientStab
void SetGradientStab(STATE sd=1.0)
Define gradient stabilization term.
Definition: pznonlinearpoisson3d.cpp:34
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZMatPoisson3d::fSD
STATE fSD
Multiplication value for the streamline diffusion term.
Definition: pzpoisson3d.h:53
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZNonLinearPoisson3d::ContributeBCInterface
virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Definition: pznonlinearpoisson3d.cpp:550
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
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
TPZMatPoisson3dReferred::SetConvectionTermInterface
void SetConvectionTermInterface(TPZFMatrix< STATE > &dsolL, TPZFMatrix< STATE > &dsolR)
Sets convection term for ContributeInterface methods.
Definition: pzpoisson3dreferred.cpp:50
TPZMatPoisson3d::fC
REAL fC
Variable which multiplies the convection term of the equation.
Definition: pzpoisson3d.h:40
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
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
TPZMatPoisson3d::fK
STATE fK
Coeficient which multiplies the Laplacian operator.
Definition: pzpoisson3d.h:37
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZNonLinearPoisson3d::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: pznonlinearpoisson3d.cpp:390
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
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
TPZNonLinearPoisson3d::TPZNonLinearPoisson3d
TPZNonLinearPoisson3d(int nummat, int dim)
Definition: pznonlinearpoisson3d.cpp:11
test.res
string res
Definition: test.py:151
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
TPZNonLinearPoisson3d::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: pznonlinearpoisson3d.cpp:302
TPZMatrix::VerifySymmetry
virtual int VerifySymmetry(REAL tol=1.e-13) const
Checks if current matrix value is symmetric.
Definition: pzmatrix.cpp:1501
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZNonLinearPoisson3d::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: pznonlinearpoisson3d.cpp:44
TPZMatPoisson3d::fConvDir
REAL fConvDir[3]
Direction of the convection operator.
Definition: pzpoisson3d.h:43
TPZMatPoisson3d::fDim
int fDim
Problem dimension.
Definition: pzpoisson3d.h:34
TPZMatPoisson3dReferred::SetConvectionTerm
void SetConvectionTerm(TPZFMatrix< STATE > &dsol, TPZFMatrix< REAL > &axes)
Sets convection term.
Definition: pzpoisson3dreferred.cpp:17
TPZNonLinearPoisson3d::SetSUPGStab
void SetSUPGStab(STATE sd=1.0)
Define SUPG stabilization term.
Definition: pznonlinearpoisson3d.cpp:29
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
TPZMatPoisson3dReferred
This class implements a version of TPZMatPoisson3d where the convection term is given at each integr...
Definition: pzpoisson3dreferred.h:22
TPZNonLinearPoisson3d::SetNoStabilizationTerm
void SetNoStabilizationTerm()
Define no stabilization term.
Definition: pznonlinearpoisson3d.cpp:39
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZNonLinearPoisson3d::fStabilizationType
int fStabilizationType
Stabilization term definition.
Definition: pznonlinearpoisson3d.h:24
TPZNonLinearPoisson3d::ClassId
int ClassId() const override
Unique identifier for serialization purposes.
Definition: pznonlinearpoisson3d.cpp:650
TPZMatPoisson3dReferred::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzpoisson3dreferred.cpp:109
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13