NeoPZ
pznonlinbiharmonic.cpp
Go to the documentation of this file.
1 
6 #include "pznonlinbiharmonic.h"
7 #include "pzelmat.h"
8 #include "pzbndcond.h"
9 #include "pzmatrix.h"
10 #include "pzfmatrix.h"
11 #include "pzerror.h"
12 #include <math.h>
13 #include <cmath>
14 
15 // NIPG SIPG SSIPG1 SSIPG2
16 STATE TPZNonLinBiharmonic::gLambda1 = 1.0; // -1 1 -1 1
17 STATE TPZNonLinBiharmonic::gLambda2 = 1.0; // -1 1 1 -1
18 STATE TPZNonLinBiharmonic::gSigmaA = 10.0;// 10 10 10 10
19 STATE TPZNonLinBiharmonic::gSigmaB = 10.0;// 10 10 10 10
20 STATE TPZNonLinBiharmonic::gL_alpha = 6.0; // [0, 6
21 STATE TPZNonLinBiharmonic::gM_alpha = 3.0; // IGUAL
22 STATE TPZNonLinBiharmonic::gL_betta = 4.0; // [-2, 4
23 STATE TPZNonLinBiharmonic::gM_betta = 1.0; // IGUAL
24 STATE TPZNonLinBiharmonic::g_teta = 0.5; // Parametro da parte advectiva.
25 STATE TPZNonLinBiharmonic::Re = 50.0; //
26 int TPZNonLinBiharmonic::NorP = 1; // Constante. Se for 1, entao Metodo de Newton
27 // Se for 0, entao Metodo de Picard
28 
29 using namespace std;
30 
31 TPZNonLinBiharmonic::TPZNonLinBiharmonic(int nummat, STATE f) : TPZRegisterClassId(&TPZNonLinBiharmonic::ClassId),
32 TPZDiscontinuousGalerkin(nummat),
33 fXf(f){}
34 
35 TPZNonLinBiharmonic::~TPZNonLinBiharmonic() {
36 }
37 
38 void TPZNonLinBiharmonic::Print(std::ostream &out) {
39  out << "name of material : " << Name() << "\n";
40  out << "properties : \n";
41  TPZMaterial::Print(out);
42 
43 }
44 
45 // Fazer outro metodo igual a esse sem os ... ek ... e com os parametros adequados
46 
47 void TPZNonLinBiharmonic::Contribute(TPZMaterialData &data,
48  REAL weight,
49  TPZFMatrix<STATE> &ek,
50  TPZFMatrix<STATE> &ef) {
51  TPZFMatrix<REAL> &dphi = data.dphix;
52  TPZFMatrix<REAL> &phi = data.phi;
53  TPZManVector<REAL,3> &x = data.x;
54  int numbersol = data.dsol.size();
55  if (numbersol != 1) {
56  DebugStop();
57  }
58 
59  TPZFMatrix<STATE> &dsol=data.dsol[0];
60 
61  int phr = phi.Rows();
62  STATE Re_1 = 1./Re;
63 
64  if(fForcingFunction) { // phi(in, 0) = phi_in
65  TPZManVector<STATE> res(1);
66  fForcingFunction->Execute(x,res); // dphi(i,j) = dphi_j/dxi
67  fXf = res[0];
68  }
69  //Equa�o de non linear biharmonic
70  for( int in = 0; in < phr; in++ ) {
71  ef(in, 0) += weight * fXf*phi(in,0)
72  - weight*(Re_1*dsol(2,0)*dphi(2,in)
73  + dsol(2,0)*(dsol(1,0)*dphi(0,in) - dsol(0,0)*dphi(1,in)) ) ;
74 
75  for( int jn = 0; jn < phr; jn++ ) {
76  ek(in,jn) += weight * ( Re_1*dphi(2,in) * dphi(2,jn)
77 
78  + dphi(2,jn)*(dsol(1,0)*dphi(0,in) - dsol(0,0)*dphi(1,in)) );
79  if(NorP==1)
80  ek(in,jn) += weight * ( dsol(2,0)*(dphi(1,jn)*dphi(0,in) - dphi(0,jn)*dphi(1,in) ) );
81  }
82  }
83 }
84 
85 
86 void TPZNonLinBiharmonic::ContributeBC(TPZMaterialData &data,
87  REAL weight,
88  TPZFMatrix<STATE> &ek,
89  TPZFMatrix<STATE> &ef,
90  TPZBndCond &bc) {
91 
92  //NOT TO BE DONE HERE
93  PZError << "TPZBiHarminic::ContributeBC - It should never be called.";
94 }
95 
97 int TPZNonLinBiharmonic::VariableIndex(const std::string &name){
98  if(!strcmp("Displacement6",name.c_str())) return 0;
99  if(!strcmp("Solution",name.c_str())) return 1;
100  if(!strcmp("Derivate",name.c_str())) return 2;
101  if(!strcmp("POrder",name.c_str())) return 10;
102  cout << "TPZNonLinBiharmonic::VariableIndex Error\n";
103  return -1;
104 }
105 
106 int TPZNonLinBiharmonic::NSolutionVariables(int var){
107  if(var == 0) return 6;
108  if(var == 1) return 1;
109  if(var == 2) return 2;
110  if(var == 10) return 1;
111  return TPZMaterial::NSolutionVariables(var);
112 }
113 
114 void TPZNonLinBiharmonic::Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &/*axes*/,
115  int var,TPZVec<STATE> &Solout){
116  if(var == 0 || var == 1) Solout[0] = Sol[0];//function
117  if(var == 2) {
118  Solout.Resize(DSol.Rows());
119  int id;
120  for(id=0 ; id < DSol.Rows(); id++) {
121  Solout[id] = DSol(id,0);//derivate
122  }
123  }
124 }
125 
126 void TPZNonLinBiharmonic::Flux(TPZVec<REAL> &/*x*/, TPZVec<STATE> &/*Sol*/,
127  TPZFMatrix<STATE> &/*DSol*/, TPZFMatrix<REAL> &/*axes*/,
128  TPZVec<STATE> &/*flux*/) {
129  //Flux(TPZVec<REAL> &x, TPZVec<STATE> &Sol, TPZFMatrix<STATE> &DSol, TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux)
130 }
131 
132 void TPZNonLinBiharmonic::Errors(TPZVec<REAL> &/*x*/,TPZVec<STATE> &u, TPZFMatrix<STATE> &dudx,
133  TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
134  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,
135  TPZVec<REAL> &values) {
136 
137  TPZVec<STATE> sol(1), dsol(8,0.);
138  Solution(u,dudx,axes,1,sol);
139  Solution(u,dudx,axes,2,dsol);
140  //values[1] : error em norma L2
141  values[1] = (sol[0] - u_exact[0])*(sol[0] - u_exact[0]);
142 
143  //values[2] : erro em semi norma H1
144  values[2] = 0.;
145  for(int id=0; id<2; id++) {
146  values[2] += (dsol[id] - du_exact(id,0))*(dsol[id] - du_exact(id,0));
147  }
148  //values[3] : erro em semi norma Laplace
149  values[3] = (dsol[2] - du_exact(2,0))*(dsol[2] - du_exact(2,0));
150 
151  //values[0] : erro em norma H1
152  values[0] = values[1]+values[2];
153  // dxx
154  values[5] = (dsol[5] - du_exact(5,0))*(dsol[5] - du_exact(5,0));
155  // dyy
156  values[6] = (dsol[6] - du_exact(6,0))*(dsol[6] - du_exact(6,0));
157  // dxy
158  values[7] = (dsol[7] - du_exact(7,0))*(dsol[7] - du_exact(7,0));
159  //values[4] : erro em norma H2
160  values[4] = values[5]+values[6]+values[7]+values[0];
161 
162 }
163 
164 void TPZNonLinBiharmonic::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright,
165  REAL weight,
166  TPZFMatrix<STATE> &ek,
167  TPZFMatrix<STATE> &ef){
168  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
169  TPZFMatrix<REAL> &dphiR = dataright.dphix;
170  TPZFMatrix<REAL> &phiL = dataleft.phi;
171  TPZFMatrix<REAL> &phiR = dataright.phi;
172  TPZManVector<REAL,3> &normal = data.normal;
173 
174  int LeftPOrder=dataleft.p;
175  int RightPOrder=dataright.p;
176  int numbersol = dataleft.sol.size();
177  if (numbersol != 1) {
178  DebugStop();
179  }
180 
181  TPZVec<STATE> &solL=dataleft.sol[0];
182  TPZVec<STATE> &solR=dataright.sol[0];
183  TPZFMatrix<STATE> &dsolL=dataleft.dsol[0];
184  TPZFMatrix<STATE> &dsolR=dataright.dsol[0];
185  STATE faceSize=data.HSize;
186 
187  int nrowl = phiL.Rows();
188  int nrowr = phiR.Rows();
189  int il,jl,ir,jr,id;
190  STATE Re_1 = 1./Re;
191 
192  STATE alpha=gSigmaA*(pow(((STATE)LeftPOrder),gL_alpha)+pow(((STATE)RightPOrder), gL_alpha) ) /
193  (2. * pow(faceSize, gM_alpha) );
194 
195 
196  STATE betta=gSigmaB*(pow(((STATE)LeftPOrder),gL_betta)+pow(((STATE)RightPOrder), gL_betta) ) /
197  (2. * pow(faceSize, gM_betta) );
198 
199  //
200  // advectivo
201  // 'Unica Integral
202  STATE norFnorF = normal[0]*normal[0];
203  STATE norSnorS = normal[1]*normal[1];
204 
205  STATE uLnorF= dsolL(1,0)*normal[0];
206  STATE vLnorS= -dsolL(0,0)*normal[1];
207 
208  STATE uRnorF= dsolR(1,0)*normal[0]; //
209  STATE vRnorS= -dsolR(0,0)*normal[1];
210 
211  STATE absUnorL = fabs(uLnorF + vLnorS);
212  STATE absUnorR = fabs(uRnorF + vRnorS);
213 
214 
215  for(il=0; il<nrowl; il++) {
216  ef(il,0) += weight*0.5*( dsolL(2,0)*phiL(il,0)*(uLnorF + vLnorS))
217  +weight*g_teta*absUnorL*(dsolL(2,0)*phiL(il,0)*norFnorF +
218  dsolL(2,0)*phiL(il,0)*norSnorS);
219 
220  for(jl=0; jl<nrowl; jl++) {
221  STATE deriv_uLnorF = dphiL(1,jl)*normal[0];
222  STATE deriv_vLnorS = -dphiL(0,jl)*normal[1];
223 
224  STATE deriv_absUnorL = 0.;
225  if(uLnorF+vLnorS > 0.) deriv_absUnorL += (deriv_uLnorF+deriv_vLnorS);
226  else deriv_absUnorL -= (deriv_uLnorF+deriv_vLnorS);
227 
228  ek(il,jl) -= weight*0.5*( dphiL(2,jl)*phiL(il,0)*(uLnorF + vLnorS))
229  +weight*g_teta*absUnorL*(dphiL(2,jl)*phiL(il,0)*norFnorF +
230  dphiL(2,jl)*phiL(il,0)*norSnorS) ;
231 
232  if(NorP==1)
233  ek(il,jl) -= + weight*0.5*(dsolL(2,0)*phiL(il,0)*(deriv_uLnorF+deriv_vLnorS)) +
234  weight*g_teta*deriv_absUnorL*(dsolL(2,0)*phiL(il,0)*norFnorF +
235  dsolL(2,0)*phiL(il,0)*norSnorS );
236  }
237  }
238 
239  for(ir=0; ir<nrowr; ir++) {
240  //
241  ef(ir+nrowl,0) += weight*0.5*(
242  - dsolR(2,0)*phiR(ir,0)*(uRnorF + vRnorS))
243  +weight*g_teta*absUnorR*(dsolR(2,0)*phiR(ir,0)*norFnorF +
244  dsolR(2,0)*phiR(ir,0)*norSnorS);
245 
246  for(jr=0; jr<nrowr; jr++) {
247  STATE deriv_uRnorF = dphiR(1,jr)*normal[0];
248  STATE deriv_vRnorS = -dphiR(0,jr)*normal[1];
249 
250  STATE deriv_absUnorR = 0.;
251  if(uRnorF + vRnorS > 0.) deriv_absUnorR += (deriv_uRnorF+deriv_vRnorS);
252  else deriv_absUnorR -= (deriv_uRnorF+deriv_vRnorS);
253 
254  ek(ir+nrowl,jr+nrowl) -= weight*0.5*(
255  - dphiR(2,jr)*phiR(ir,0)*(uRnorF + vRnorS))
256  +weight*g_teta*absUnorR*(dphiR(2,jr)*phiR(ir,0)*norFnorF +
257  dphiR(2,jr)*phiR(ir,0)*norSnorS );
258  if(NorP==1)
259  ek(ir+nrowl,jr+nrowl) -= + weight*0.5*(
260  - dsolR(2,0)*phiR(ir,0)*(deriv_uRnorF + deriv_vRnorS))
261  +weight*g_teta*deriv_absUnorR*(dsolR(2,0)*phiR(ir,0)*norFnorF +
262  dsolR(2,0)*phiR(ir,0)*norSnorS);
263  }
264  }
265 
266 
267  for(il=0; il<nrowl; il++) {
268  //
269  ef(il,0) += weight*0.5*(
270  dsolR(2,0)*phiL(il,0)*(uRnorF + vRnorS))
271  -weight*g_teta*absUnorR*(dsolR(2,0)*phiL(il,0)*norFnorF +
272  dsolR(2,0)*phiL(il,0)*norSnorS);
273 
274 
275  for(jr=0; jr<nrowr; jr++) {
276  STATE deriv_uRnorF = dphiR(1,jr)*normal[0];
277  STATE deriv_vRnorS = -dphiR(0,jr)*normal[1];
278 
279  STATE deriv_absUnorR = 0.;
280  if(uRnorF + vRnorS > 0.) deriv_absUnorR += (deriv_uRnorF+deriv_vRnorS);
281  else deriv_absUnorR -= (deriv_uRnorF+deriv_vRnorS);
282 
283  ek(il,jr+nrowl) -= weight*0.5*(
284  dphiR(2,jr)*phiL(il,0)*(uRnorF + vRnorS))
285  -weight*g_teta*absUnorR*(dphiR(2,jr)*phiL(il,0)*norFnorF +
286  dphiR(2,jr)*phiL(il,0)*norSnorS) ;
287 
288  if(NorP==1)
289  ek(il,jr+nrowl) -= +weight*0.5*(
290  dsolR(2,0)*phiL(il,0)*(deriv_uRnorF + deriv_vRnorS))
291  -weight*g_teta*deriv_absUnorR*(dsolR(2,0)*phiL(il,0)*norFnorF +
292  dsolR(2,0)*phiL(il,0)*norSnorS) ;
293  }
294  }
295 
296  for(ir=0; ir<nrowr; ir++) {
297  ef(ir+nrowl,0) += weight*0.5*(
298  - dsolL(2,0)*phiR(ir,0)*(uLnorF + vLnorS))
299  -weight*g_teta*absUnorL*(dsolL(2,0)*phiR(ir,0)*norFnorF +
300  dsolL(2,0)*phiR(ir,0)*norSnorS);
301 
302  for(jl=0; jl<nrowl; jl++) {
303  STATE deriv_uLnorF = dphiL(1,jl)*normal[0];
304  STATE deriv_vLnorS = -dphiL(0,jl)*normal[1];
305 
306  STATE deriv_absUnorL = 0.;
307  if(uLnorF+vLnorS > 0.) deriv_absUnorL += (deriv_uLnorF+deriv_vLnorS);
308  else deriv_absUnorL -= (deriv_uLnorF+deriv_vLnorS);
309 
310  ek(ir+nrowl,jl) -= weight*0.5*(
311  - dphiL(2,jl)*phiR(ir,0)*(uLnorF + vLnorS))
312  -weight*g_teta*absUnorL*(dphiL(2,jl)*phiR(ir,0)*norFnorF +
313  dphiL(2,jl)*phiR(ir,0)*norSnorS );
314  if(NorP==1)
315  ek(ir+nrowl,jl) -= +weight*0.5*(
316  - dsolL(2,0)*phiR(ir,0)*(deriv_uLnorF + deriv_vLnorS))
317 
318  -weight*g_teta*deriv_absUnorL*(dsolL(2,0)*phiR(ir,0)*norFnorF +
319  dsolL(2,0)*phiR(ir,0)*norSnorS);
320 
321  }
322  }
323 
324  /*
325  * biharmonic
326  */
327  /* Primeira Integral */
328  for(il=0; il<nrowl; il++) {
329  STATE dphiLinormal = 0.;
330  for(id=0; id<2; id++) {
331  dphiLinormal += dphiL(3+id,il)*normal[id];
332  }
333  for(jl=0; jl<nrowl; jl++) {
334  STATE dphiLjnormal = 0.;
335  for(id=0; id<2; id++) {
336  dphiLjnormal += dphiL(3+id,jl)*normal[id];
337  }
338 
339  ek(il,jl) += Re_1*weight*0.5*(+gLambda1*dphiLinormal*phiL(jl,0) + dphiLjnormal*phiL(il,0));
340  }
341  STATE dsolLnormal = 0.;
342  for(id=0; id<2; id++) {
343  dsolLnormal += dsolL(3+id,0)*normal[id];
344  }
345  ef(il,0) -= Re_1*weight*0.5*(+gLambda1*dphiLinormal*solL[0] + dsolLnormal*phiL(il,0));
346  }
347 
348  for(ir=0; ir<nrowr; ir++) {
349  STATE dphiRinormal = 0.;
350  for(id=0; id<2; id++) {
351  dphiRinormal += dphiR(3+id,ir)*normal[id];
352  }
353  for(jr=0; jr<nrowr; jr++) {
354  STATE dphiRjnormal = 0.;
355  for(id=0; id<2; id++) {
356  dphiRjnormal += dphiR(3+id,jr)*normal[id];
357  }
358  ek(ir+nrowl,jr+nrowl) += Re_1*weight*0.5*(
359  -gLambda1*dphiRinormal*phiR(jr,0) - dphiRjnormal*phiR(ir,0));
360  }
361  STATE dsolRnormal = 0.;
362  for(id=0; id<2; id++) {
363  dsolRnormal += dsolR(3+id,0)*normal[id];
364  }
365  ef(ir+nrowl,0) -= Re_1*weight*0.5*(
366  -gLambda1*dphiRinormal*solR[0] - dsolRnormal*phiR(ir,0));
367 
368  }
369 
370  for(il=0; il<nrowl; il++) {
371  STATE dphiLinormal = 0.;
372  for(id=0; id<2; id++) {
373  dphiLinormal += dphiL(3+id,il)*normal[id];
374  }
375  for(jr=0; jr<nrowr; jr++) {
376  STATE dphiRjnormal = 0.;
377  for(id=0; id<2; id++) {
378  dphiRjnormal += dphiR(3+id,jr)*normal[id];
379  }
380  ek(il,jr+nrowl) += Re_1*weight*0.5*(
381  - gLambda1*dphiLinormal*phiR(jr,0) + dphiRjnormal*phiL(il,0));
382  }
383  STATE dsolRnormal = 0.;
384  for(id=0; id<2; id++) {
385  dsolRnormal += dsolR(3+id,0)*normal[id];
386  }
387  ef(il,0) -= Re_1*weight*0.5*(
388  - gLambda1*dphiLinormal*solR[0] + dsolRnormal*phiL(il,0));
389 
390  }
391  for(ir=0; ir<nrowr; ir++) {
392  STATE dphiRinormal = 0.;
393  for(id=0; id<2; id++) {
394  dphiRinormal += dphiR(3+id,ir)*normal[id];
395  }
396  for(jl=0; jl<nrowl; jl++) {
397  STATE dphiLjnormal = 0.;
398  for(id=0; id<2; id++) {
399  dphiLjnormal += dphiL(3+id,jl)*normal[id];
400  }
401  ek(ir+nrowl,jl) += Re_1*weight*0.5*(
402  + gLambda1*dphiRinormal*phiL(jl,0) - dphiLjnormal*phiR(ir,0));
403  }
404  STATE dsolLnormal = 0.;
405  for(id=0; id<2; id++) {
406  dsolLnormal += dsolL(3+id,0)*normal[id];
407  }
408  ef(ir+nrowl,0) -= Re_1*weight*0.5*(
409  + gLambda1*dphiRinormal*solL[0] - dsolLnormal*phiR(ir,0));
410  }
411 
412  /* Segunda Integral */
413  for(il=0; il<nrowl; il++) {
414  STATE dphiLinormal = 0.;
415  for(id=0; id<2; id++) {
416  dphiLinormal += dphiL(id,il)*normal[id];
417  }
418  for(jl=0; jl<nrowl; jl++) {
419  STATE dphiLjnormal = 0.;
420  for(id=0; id<2; id++) {
421  dphiLjnormal += dphiL(id,jl)*normal[id];
422  }
423 
424  ek(il,jl) += Re_1*weight*0.5*(
425  - dphiLinormal*dphiL(2,jl) - gLambda2*dphiLjnormal*dphiL(2,il));
426  }
427  STATE dsolLnormal = 0.;
428  for(id=0; id<2; id++) {
429  dsolLnormal += dsolL(id,0)*normal[id];
430  }
431  ef(il,0) -= Re_1*weight*0.5*(
432  - dphiLinormal*dsolL(2,0) - gLambda2*dsolLnormal*dphiL(2,il) );
433 
434  }
435  for(ir=0; ir<nrowr; ir++) {
436  STATE dphiRinormal = 0.;
437  for(id=0; id<2; id++) {
438  dphiRinormal += dphiR(id,ir)*normal[id];
439  }
440  for(jr=0; jr<nrowr; jr++) {
441  STATE dphiRjnormal = 0.;
442  for(id=0; id<2; id++) {
443  dphiRjnormal += dphiR(id,jr)*normal[id];
444  }
445  ek(ir+nrowl,jr+nrowl) += Re_1*weight*0.5*(
446  + dphiRinormal*dphiR(2,jr) + gLambda2*dphiRjnormal*dphiR(2,ir));
447  }
448  STATE dsolRnormal = 0.;
449  for(id=0; id<2; id++) {
450  dsolRnormal += dsolR(id,0)*normal[id];
451  }
452  ef(ir+nrowl,0) -= Re_1*weight*0.5*(
453  + dphiRinormal*dsolR(2,0) + gLambda2*dsolRnormal*dphiR(2,ir));
454  }
455  for(il=0; il<nrowl; il++) {
456  STATE dphiLinormal = 0.;
457  for(id=0; id<2; id++) {
458  dphiLinormal += dphiL(id,il)*normal[id];
459  }
460  for(jr=0; jr<nrowr; jr++) {
461  STATE dphiRjnormal = 0.;
462  for(id=0; id<2; id++) {
463  dphiRjnormal += dphiR(id,jr)*normal[id];
464  }
465  ek(il,jr+nrowl) += Re_1*weight*0.5*(
466  - dphiLinormal*dphiR(2,jr) + gLambda2*dphiRjnormal*dphiL(2,il));
467  }
468  STATE dsolRnormal = 0.;
469  for(id=0; id<2; id++) {
470  dsolRnormal += dsolR(id,0)*normal[id];
471  }
472  ef(il,0) -= Re_1*weight*0.5*(
473  - dphiLinormal*dsolR(2,0) + gLambda2*dsolRnormal*dphiL(2,il));
474  }
475 
476  for(ir=0; ir<nrowr; ir++) {
477  STATE dphiRinormal = 0.;
478  for(id=0; id<2; id++) {
479  dphiRinormal += dphiR(id,ir)*normal[id];
480  }
481  for(jl=0; jl<nrowl; jl++) {
482  STATE dphiLjnormal = 0.;
483  for(id=0; id<2; id++) {
484  dphiLjnormal += dphiL(id,jl)*normal[id];
485  }
486  ek(ir+nrowl,jl) += Re_1*weight*0.5*(
487  + dphiRinormal*dphiL(2,jl) - gLambda2*dphiLjnormal*dphiR(2,ir));
488  }
489  STATE dsolLnormal = 0.;
490  for(id=0; id<2; id++) {
491  dsolLnormal += dsolL(id,0)*normal[id];
492  }
493  ef(ir+nrowl,0) -= Re_1*weight*0.5*(
494  + dphiRinormal*dsolL(2,0) - gLambda2*dsolLnormal*dphiR(2,ir));
495 
496  }
497 
498  /* Terceira Integral */
499  for(il=0; il<nrowl; il++) {
500  STATE dphiLinormal = 0.;
501  for(id=0; id<2; id++) {
502  dphiLinormal += dphiL(id,il)*normal[id];
503  }
504  for(jl=0; jl<nrowl; jl++) {
505  STATE dphiLjnormal = 0.;
506  for(id=0; id<2; id++) {
507  dphiLjnormal += dphiL(id,jl)*normal[id];
508  }
509 
510  ek(il,jl) += Re_1*weight*(
511  alpha * phiL(jl,0)*phiL(il,0) +
512  betta * dphiLinormal*dphiLjnormal);
513  }
514  STATE dsolLnormal = 0.;
515  for(id=0; id<2; id++) {
516  dsolLnormal += dsolL(id,0)*normal[id];
517  }
518  ef(il,0) -= Re_1*weight*(
519  alpha * solL[0]*phiL(il,0) +
520  betta * dphiLinormal*dsolLnormal);
521  }
522 
523  for(ir=0; ir<nrowr; ir++) {
524  STATE dphiRinormal = 0.;
525  for(id=0; id<2; id++) {
526  dphiRinormal += dphiR(id,ir)*normal[id];
527  }
528  for(jr=0; jr<nrowr; jr++) {
529  STATE dphiRjnormal = 0.;
530  for(id=0; id<2; id++) {
531  dphiRjnormal += dphiR(id,jr)*normal[id];
532  }
533  ek(ir+nrowl,jr+nrowl) += Re_1*weight*(
534  alpha * phiR(jr,0)*phiR(ir,0) +
535  betta * dphiRinormal*dphiRjnormal );
536  }
537  STATE dsolRnormal = 0.;
538  for(id=0; id<2; id++) {
539  dsolRnormal += dsolR(id,0)*normal[id];
540  }
541  ef(ir+nrowl,0) -= Re_1*weight*(
542  alpha * solR[0]*phiR(ir,0) +
543  betta * dphiRinormal*dsolRnormal );
544  }
545 
546  for(il=0; il<nrowl; il++) {
547  STATE dphiLinormal = 0.;
548  for(id=0; id<2; id++) {
549  dphiLinormal += dphiL(id,il)*normal[id];
550  }
551  for(jr=0; jr<nrowr; jr++) {
552  STATE dphiRjnormal = 0.;
553  for(id=0; id<2; id++) {
554  dphiRjnormal += dphiR(id,jr)*normal[id];
555  }
556  ek(il,jr+nrowl) += Re_1*weight*(
557  - alpha * phiR(jr,0)*phiL(il,0)
558  - betta * dphiLinormal*dphiRjnormal);
559  }
560  STATE dsolRnormal = 0.;
561  for(id=0; id<2; id++) {
562  dsolRnormal += dsolR(id,0)*normal[id];
563  }
564  ef(il,0) -= Re_1*weight*(
565  - alpha * solR[0]*phiL(il,0)
566  - betta * dphiLinormal*dsolRnormal);
567  }
568 
569  for(ir=0; ir<nrowr; ir++) {
570  STATE dphiRinormal = 0.;
571  for(id=0; id<2; id++) {
572  dphiRinormal += dphiR(id,ir)*normal[id];
573  }
574  for(jl=0; jl<nrowl; jl++) {
575  STATE dphiLjnormal = 0.;
576  for(id=0; id<2; id++) {
577  dphiLjnormal += dphiL(id,jl)*normal[id];
578  }
579  ek(ir+nrowl,jl) += Re_1*weight*(
580  - alpha * phiL(jl,0)*phiR(ir,0)
581  - betta * dphiRinormal*dphiLjnormal);
582  }
583  STATE dsolLnormal = 0.;
584  for(id=0; id<2; id++) {
585  dsolLnormal += dsolL(id,0)*normal[id];
586  }
587  ef(ir+nrowl,0) -= Re_1*weight*(
588  - alpha * solL[0]*phiR(ir,0)
589  - betta * dphiRinormal*dsolLnormal);
590  }
591 }
592 
593 void TPZNonLinBiharmonic::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
594  REAL weight,
595  TPZFMatrix<STATE> &ek,
596  TPZFMatrix<STATE> &ef,
597  TPZBndCond &bc) {
598  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
599  TPZFMatrix<REAL> &phiL = dataleft.phi;
600  TPZManVector<REAL,3> &normal = data.normal;
601  int POrder=data.p;
602  int numbersol = dataleft.sol.size();
603  if (numbersol != 1) {
604  DebugStop();
605  }
606 
607  TPZVec<STATE> &solL=dataleft.sol[0];
608  TPZFMatrix<STATE> &dsolL=dataleft.dsol[0];
609  STATE faceSize=data.HSize;
610 
611  STATE alpha = gSigmaA*pow(((STATE)POrder), gL_alpha) / pow(faceSize, gM_alpha);
612  STATE betta = gSigmaB*pow(((STATE)POrder), gL_betta) / pow(faceSize, gM_betta);
613 
614  int il,jl,nrowl,id;
615  nrowl = phiL.Rows();
616  STATE Re_1 = 1./Re;
617 
618  /* Primeira Integral */
619  for(il=0; il<nrowl; il++) {
620 
621  STATE dphiLinormal = 0.;
622  for(id=0; id<2; id++) {
623  dphiLinormal += dphiL(3+id,il)*normal[id];
624  }
625 
626  // Termos de Dirichlet - em Val2()(0,0)
627  ef(il,0) += + Re_1*gLambda1*weight*dphiLinormal*bc.Val2()(0,0)
628  + Re_1*alpha * weight*bc.Val2()(0,0)*phiL(il) ;
629 
630  for(jl=0; jl<nrowl; jl++) {
631  STATE dphiLjnormal = 0.;
632  for(id=0; id<2; id++) {
633  dphiLjnormal += dphiL(3+id,jl)*normal[id];
634  }
635  ek(il,jl) += Re_1*weight*(+ gLambda1*dphiLinormal*phiL(jl,0)+ dphiLjnormal*phiL(il,0)); //2 1
636  }
637  STATE dsolLnormal = 0.;
638  for(id=0; id<2; id++) {
639  dsolLnormal += dsolL(3+id,0)*normal[id];
640  }
641  ef(il,0) -= Re_1*weight*(+ gLambda1*dphiLinormal*solL[0]+ dsolLnormal*phiL(il,0));
642 
643  }
644 
645  /* Segunda Integral */
646  for(il=0; il<nrowl; il++) {
647  STATE dphiLinormal = 0.;
648  for(id=0; id<2; id++) {
649  dphiLinormal += dphiL(id,il)*normal[id];
650  }
651 
652  // Termos de Neuwmann - em Val2()(1,0)
653  ef(il,0) += - Re_1*gLambda2*weight*bc.Val2()(1,0)*dphiL(2,il)
654  + Re_1*betta * weight*bc.Val2()(1,0)*dphiLinormal ;
655 
656  for(jl=0; jl<nrowl; jl++) {
657  STATE dphiLjnormal = 0.;
658  for(id=0; id<2; id++) {
659  dphiLjnormal += dphiL(id,jl)*normal[id];
660  }
661  ek(il,jl) += Re_1*weight*(- dphiLinormal*dphiL(2,jl) - gLambda2*dphiLjnormal*dphiL(2,il) );
662  }
663 
664  STATE dsolLnormal = 0.;
665  for(id=0; id<2; id++) {
666  dsolLnormal += dsolL(id,0)*normal[id];
667  }
668  ef(il,0) -= Re_1*weight*(- dphiLinormal*dsolL(2,0) - gLambda2*dsolLnormal*dphiL(2,il) );
669 
670  }
671 
672  for(il=0; il<nrowl; il++) {
673  STATE dphiLinormal = 0.;
674  for(id=0; id<2; id++) {
675  dphiLinormal += dphiL(id,il)*normal[id];
676  }
677  for(jl=0; jl<nrowl; jl++) {
678  STATE dphiLjnormal = 0.;
679  for(id=0; id<2; id++) {
680  dphiLjnormal += dphiL(id,jl)*normal[id];
681  }
682  ek(il,jl) += Re_1*weight*(
683  alpha * phiL(jl,0)*phiL(il,0) +
684  betta * dphiLinormal*dphiLjnormal);
685  }
686  STATE dsolLnormal = 0.;
687  for(id=0; id<2; id++) {
688  dsolLnormal += dsolL(id,0)*normal[id];
689  }
690  ef(il,0) -= Re_1*weight*(
691  alpha * solL[0]*phiL(il,0) +
692  betta * dphiLinormal*dsolLnormal);
693 
694  }
695 
696  /*
697  * Termo advectivo
698  */
699 
700  STATE ULnormal = dsolL(1,0)*normal[0] - dsolL(0,0)*normal[1];
701  for(il=0; il<nrowl; il++){
702  for(jl=0; jl<nrowl; jl++){
703  ek(il,jl) -= weight*dphiL(2,jl)*phiL(il,0)*ULnormal;
704 
705  if(NorP==1)
706  ek(il,jl) -= +weight*dsolL(2,0)*phiL(il,0)*( dphiL(1,jl)*normal[0] - dphiL(0,jl)*normal[1]);
707  }
708 
709  ef(il,0) += weight*dsolL(2,0)*phiL(il,0)*ULnormal;
710 
711  }
712 
713 }
714 
715 int TPZNonLinBiharmonic::ClassId() const{
716  return Hash("TPZNonLinBiharmonic") ^ TPZDiscontinuousGalerkin::ClassId() << 1;
717 }
TPZDiscontinuousGalerkin
Defines the interface which material objects need to implement for discontinuous Galerkin formulation...
Definition: pzdiscgal.h:20
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZNonLinBiharmonic::ContributeBCInterface
virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to stiffness matrix and load vector at one BC integration point...
Definition: pznonlinbiharmonic.cpp:593
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
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
TPZNonLinBiharmonic::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Implements boundary conditions for continuous Galerkin.
Definition: pznonlinbiharmonic.cpp:86
pznonlinbiharmonic.h
Contains the TPZNonLinBiharmonic class which implements a discontinuous Galerkin formulation for the ...
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZNonLinBiharmonic::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: pznonlinbiharmonic.cpp:126
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
pzerror.h
Defines PZError.
TPZNonLinBiharmonic
This class implements discontinuous Galerkin formulation for the non-linear bi-harmonic equation...
Definition: pznonlinbiharmonic.h:19
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZNonLinBiharmonic::Name
virtual std::string Name() override
Returns the name of the material.
Definition: pznonlinbiharmonic.h:51
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
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...
TPZNonLinBiharmonic::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
Definition: pznonlinbiharmonic.cpp:114
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
TPZNonLinBiharmonic::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: pznonlinbiharmonic.cpp:106
test.f
f
Definition: test.py:287
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
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
test.res
string res
Definition: test.py:151
TPZMaterialData::HSize
REAL HSize
measure of the size of the element
Definition: pzmaterialdata.h:83
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZNonLinBiharmonic::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Implements integral over element&#39;s volume.
Definition: pznonlinbiharmonic.cpp:47
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZNonLinBiharmonic::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: pznonlinbiharmonic.cpp:38
TPZNonLinBiharmonic::TPZNonLinBiharmonic
TPZNonLinBiharmonic(int nummat, STATE f)
Inicialisation of biharmonic material.
Definition: pznonlinbiharmonic.cpp:31
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZNonLinBiharmonic::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: pznonlinbiharmonic.cpp:164
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZNonLinBiharmonic::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: pznonlinbiharmonic.cpp:715
TPZNonLinBiharmonic::Errors
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: pznonlinbiharmonic.cpp:132
rdt.values
def values
Definition: rdt.py:119
TPZNonLinBiharmonic::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: pznonlinbiharmonic.cpp:97
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13