NeoPZ
tpzbiharmonicestimator.cpp
Go to the documentation of this file.
1 
7 #include "pzbiharmonic.h"
8 #include "pzbndcond.h"
9 
10 #include <cmath>
11 using namespace std;
12 
15 {
16  this->SetExactSolutions(NULL, NULL);
17 }
18 
20 {
21 }
22 
24  void (*fp)(TPZVec<REAL> &loc,TPZVec<STATE> &val,TPZFMatrix<STATE> &deriv),
25  void (*fd)(TPZVec<REAL> &locdual,TPZVec<STATE> &valdual,TPZFMatrix<STATE> &derivdual)){
26  this->fPrimalExactSol = fp;
27  this->fDualExactSol = fd;
28 }
29 
30 //static STATE Pi = 4.*atan(1.);
32 std::ofstream CE ( "debugContributeErros.txt" );
33 // std::ofstream CoutPontual ( "solPontual.txt" );
35  REAL weight,
36  TPZVec<REAL> &nk){
37  TPZVec<STATE> u_exactp(1);
38  TPZFMatrix<STATE> du_exactp(9,1);
39  TPZVec<STATE> u_exactd(1);
40  TPZFMatrix<STATE> du_exactd(9,1);
41  this->fPrimalExactSol(data.x,u_exactp,du_exactp);
42  this->fDualExactSol(data.x,u_exactd,du_exactd);
43 
44  int numbersol = data.sol.size();
45  if (numbersol != 1) {
46  DebugStop();
47  }
48 
49  TPZManVector<REAL,3> x = data.x;
50  TPZVec<STATE> sol = data.sol[0];
51  TPZFMatrix<STATE> dsol = data.dsol[0];
52  // REAL Laplac = dsol(2,0);
53  STATE LaplacLaplac = dsol(8,0);
54  // REAL DuDx = dsol(0,0);
55  // REAL DuDy = dsol(1,0);
56 
57  if (this->fForcingFunction)
58  {
59  TPZManVector<STATE, 1> result(1);
60  this->fForcingFunction->Execute(x, result);
61  this->fXf = result[0];
62  }
63  CE << "x: "<< x[0]<<","<<x[1]<<std::endl;
64  CE << "sol: "<< sol[0]<<","<<sol[1]<<","<<sol[2]<<","<<sol[3]<<std::endl;
65  CE << "weight: "<< weight <<std::endl;
66  CE << "LaplacLaplac: "<< LaplacLaplac <<std::endl;
67 
68  // std::cout.precision(16);
69  // if(fabs(x[0]-0.5)<0.001 && fabs(x[1]-0.5)<0.001){
70  // CoutPontual << "x: "<< x[0]<<","<<x[1]<<std::endl;
71  // CoutPontual << "sol: "<< sol[0]<<std::endl;
72  // }
73  STATE f = this->fXf;
74  CE << "f: "<<f <<std::endl;
75  CE << " "<<std::endl;
76  STATE Z = sol[2];
77  STATE Z1 = sol[3];
78  nk[0] += weight *( f - LaplacLaplac )*(Z1-Z);
79  nk[1] += weight *( f - LaplacLaplac )*(u_exactd[0]-Z);//segundo a teoria deve dar o erro exato do funcional(a menos do erro na integracao numerica)!
80 }
81 
83  REAL weight,
84  TPZVec<STATE> &nkL,
85  TPZVec<STATE> &nkR ){
86  TPZManVector<REAL,3> normal = data.normal;
87  TPZManVector<REAL,3> x = data.x;
88  int LeftPOrder = dataleft.p;
89  int RightPOrder = dataright.p;
90  int numbersol = dataleft.sol.size();
91  if (numbersol != 1) {
92  DebugStop();
93  }
94 
95  TPZVec<STATE> solL=dataleft.sol[0];
96  TPZVec<STATE> solR=dataright.sol[0];
97  TPZFMatrix<STATE> dsolL=dataleft.dsol[0];
98  TPZFMatrix<STATE> dsolR=dataright.dsol[0];
99  STATE faceSize=data.HSize;
100 
101  STATE alpha=gSigmaA*(pow((STATE)LeftPOrder,(STATE)gL_alpha)+pow((STATE)RightPOrder,(STATE)gL_alpha)) /(2.*pow(faceSize,(STATE)gM_alpha));
102  STATE betta=gSigmaB*(pow((STATE)LeftPOrder,(STATE)gL_betta)+pow((STATE)RightPOrder,(STATE)gL_betta)) /(2.*pow(faceSize,(STATE)gM_betta));
103 
104  // solucoes exatas
105  TPZVec<STATE> u_exactp(1);
106  TPZFMatrix<STATE> du_exactp(9,1);
107  TPZVec<STATE> u_exactd(1);
108  TPZFMatrix<STATE> du_exactd(9,1);
109  this->fPrimalExactSol(data.x,u_exactp,du_exactp);
110  this->fDualExactSol(data.x,u_exactd,du_exactd);
111 
112 
113  // u -> left e right
114  STATE uL = solL[0];
115  STATE uR = solR[0];
116  STATE graduL[2];
117  graduL[0] = dsolL(0,0);
118  graduL[1] = dsolL(1,0);
119  STATE graduR[2];
120  graduR[0] = dsolR(0,0);
121  graduR[1] = dsolR(1,0);
122  STATE LaplacianoUL =dsolL(2,0);
123  STATE LaplacianoUR =dsolR(2,0);
124  STATE GRADLaplacianoUL[2];
125  GRADLaplacianoUL[0] =dsolL(3,0);
126  GRADLaplacianoUL[1] =dsolL(4,0);
127  STATE GRADLaplacianoUR[2];
128  GRADLaplacianoUR[0] =dsolR(3,0);
129  GRADLaplacianoUR[1] =dsolR(4,0);
130 
131  // dualhp -> left e right
132  STATE dualhpL = solL[2];
133  STATE dualhpR = solR[2];
134  STATE graddualhpL[2];
135  graddualhpL[0] = dsolL(0,2);
136  graddualhpL[1] = dsolL(1,2);
137  STATE graddualhpR[2];
138  graddualhpR[0] = dsolR(0,2);
139  graddualhpR[1] = dsolR(1,2);
140  STATE LaplacianodualhpL = dsolL(2,2);
141  STATE LaplacianodualhpR = dsolR(2,2);
142  STATE GRADLaplacianodualhpL[2];
143  GRADLaplacianodualhpL[0] = dsolL(3,2);
144  GRADLaplacianodualhpL[1] = dsolL(4,2);
145  STATE GRADLaplacianodualhpR[2];
146  GRADLaplacianodualhpR[0] = dsolR(3,2);
147  GRADLaplacianodualhpR[1] = dsolR(4,2);
148 
149  // dualhpPlus -> left e right
150  STATE dualhpPlusL = solL[3];
151  STATE dualhpPlusR = solR[3];
152  STATE graddualhpPlusL[2];
153  graddualhpPlusL[0] = dsolL(0,3);
154  graddualhpPlusL[1] = dsolL(1,3);
155  STATE graddualhpPlusR[2];
156  graddualhpPlusR[0] = dsolR(0,3);
157  graddualhpPlusR[1] = dsolR(1,3);
158  STATE LaplacianodualhpPlusL = dsolL(2,3);
159  STATE LaplacianodualhpPlusR = dsolR(2,3);
160  STATE GRADLaplacianodualhpPlusL[2];
161  GRADLaplacianodualhpPlusL[0] = dsolL(3,3);
162  GRADLaplacianodualhpPlusL[1] = dsolL(4,3);
163  STATE GRADLaplacianodualhpPlusR[2];
164  GRADLaplacianodualhpPlusR[0] = dsolR(3,3);
165  GRADLaplacianodualhpPlusR[1] = dsolR(4,3);
166 
167  // sao seis termos com left e right
168 
169  // primeiro termo (s/ penalizacao)
170  nkL[0]+=-1.0*weight*
171  (LaplacianoUL-LaplacianoUR)*
172  0.5*(normal[0]*(graddualhpPlusL[0]-graddualhpL[0])+normal[1]*(graddualhpPlusL[1]-graddualhpL[1]));
173  nkR[0]+=-1.0*weight*
174  (LaplacianoUL-LaplacianoUR)*
175  (/*-*/0.5)*(normal[0]*(graddualhpPlusR[0]-graddualhpR[0])+normal[1]*(graddualhpPlusR[1]-graddualhpR[1]));
176  // segundo termo (s/ penalizacao)
177  nkL[0]+=weight*((GRADLaplacianoUL[0]-GRADLaplacianoUR[0])*normal[0]+(GRADLaplacianoUL[1]-GRADLaplacianoUR[1])*normal[1])*0.5*(dualhpPlusL-dualhpL);
178  nkR[0]+=/*-*/weight*((GRADLaplacianoUL[0]-GRADLaplacianoUR[0])*normal[0]+(GRADLaplacianoUL[1]-GRADLaplacianoUR[1])*normal[1])*0.5*(dualhpPlusR-dualhpR);
179  // terceiro termo (s/ penalizacao)
180  nkL[0]+=-weight*(uL-uR)*0.5*(normal[0]*(GRADLaplacianodualhpPlusL[0]-GRADLaplacianodualhpL[0])+normal[1]*(GRADLaplacianodualhpPlusL[1]-GRADLaplacianodualhpL[1]));
181  nkR[0]+=-weight*(uL-uR)*0.5*(/*-*/normal[0]*(GRADLaplacianodualhpPlusR[0]-GRADLaplacianodualhpR[0])/*-*/+normal[1]*(GRADLaplacianodualhpPlusR[1]-GRADLaplacianodualhpR[1]));
182  // quarto termo ( s/ penalizacao)
183  nkL[0]+=weight*((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1])*0.5*(LaplacianodualhpPlusL-LaplacianodualhpL);
184  nkR[0]+=/*-*/weight*((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1])*0.5*(LaplacianodualhpPlusR-LaplacianodualhpR);
185 
186  // quinto termo
187  nkL[0]+=-weight*alpha*(uL-uR)*(dualhpPlusL-dualhpL);
188  nkR[0]+=-weight*alpha*(uL-uR)*(-1.)*(dualhpPlusR-dualhpR);
189  //sexto termo
190  nkL[0]+=-weight*betta*((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1]) * (normal[0]*(graddualhpPlusL[0]-graddualhpL[0])+normal[1]*(graddualhpPlusL[1]-graddualhpL[1]));
191  nkR[0]+=-weight*betta*((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1])*
192  (-1.)*(normal[0]*(graddualhpPlusR[0]-graddualhpR[0])+normal[1]*(graddualhpPlusR[1]-graddualhpR[1]));
193 
194  // com sol's exatas
195 
196  // primeiro termo
197  nkL[1]+=-weight*(LaplacianoUL-LaplacianoUR)*0.5*(normal[0]*(du_exactd[0]-graddualhpL[0])+normal[1]*(du_exactd[1]-graddualhpL[1]));
198  nkR[1]+=-weight*(LaplacianoUL-LaplacianoUR)*0.5*(normal[0]*(du_exactd[0]-graddualhpR[0])+normal[1]*(du_exactd[1]-graddualhpR[1]));
199  // segundo termo
200  nkL[1]+=weight*((GRADLaplacianoUL[0]-GRADLaplacianoUR[0])*normal[0]+(GRADLaplacianoUL[1]-GRADLaplacianoUR[1])*normal[1])*0.5*(u_exactd[0]-dualhpL);
201  nkR[1]+=weight*((GRADLaplacianoUL[0]-GRADLaplacianoUR[0])*normal[0]+(GRADLaplacianoUL[1]-GRADLaplacianoUR[1])*normal[1])*0.5*(u_exactd[0]-dualhpR);
202  // terceiro termo
203  nkL[1]+=-weight*(uL-uR)*0.5*(normal[0]*(du_exactd[3]-GRADLaplacianodualhpL[0])+normal[1]*(du_exactd[4]-GRADLaplacianodualhpL[1]));
204  nkR[1]+=-weight*(uL-uR)*0.5*(normal[0]*(du_exactd[3]-GRADLaplacianodualhpR[0])+normal[1]*(du_exactd[4]-GRADLaplacianodualhpR[1]));
205  // quarto termo
206  nkL[1]+=weight*((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1])*0.5*(du_exactd[2]-LaplacianodualhpL);
207  nkR[1]+=weight*((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1])*0.5*(du_exactd[2]-LaplacianodualhpR);
208  // quinto termo
209  nkL[1]+=-weight*alpha*(uL-uR)*(u_exactd[0]-dualhpL);
210  nkR[1]+=-weight*alpha*(uL-uR)*(-1.)*(u_exactd[0]-dualhpR);
211  //sexto termo
212  nkL[1]+=-weight* /*betaL*/ ((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1]) * (normal[0]*(du_exactd[0]-graddualhpL[0])+normal[1]*(du_exactd[1]-graddualhpL[1]));
213  nkR[1]+=-weight* /*betaR*/ ((graduL[0]-graduR[0])*normal[0]+(graduL[1]-graduR[1])*normal[1]) * (-1.)*(normal[0]*(du_exactd[0]-graddualhpR[0])+normal[1]*(du_exactd[1]-graddualhpR[1]));//acertar betas!!!!!!!!
214 
215 }
216 
218  REAL weight,
219  TPZVec<STATE> &nk,
220  TPZBndCond &bc){
221  TPZFMatrix<REAL> dphi = data.dphix;
222  TPZFMatrix<REAL> dphiL = dataleft.dphix;
223  TPZFMatrix<REAL> phi = data.phi;
224  TPZFMatrix<REAL> phiL = dataleft.phi;
225  TPZManVector<REAL,3> normal = data.normal;
226  TPZManVector<REAL,3> x = data.x;
227  int LeftPOrder=dataleft.p;
228  int numbersol = dataleft.sol.size();
229  if (numbersol != 1) {
230  DebugStop();
231  }
232 
233  TPZVec<STATE> solL=dataleft.sol[0];
234  TPZFMatrix<STATE> dsolL=dataleft.dsol[0];
235  STATE faceSize=data.HSize;
236 
237  TPZVec<STATE> u_exactp(1);
238  TPZFMatrix<STATE> du_exactp(9,1);
239  TPZVec<STATE> u_exactd(1);
240  TPZFMatrix<STATE> du_exactd(9,1);
241  this->fDualExactSol(data.x,u_exactd,du_exactd);
242  this->fPrimalExactSol(data.x,u_exactp,du_exactp);
243 
244  STATE alpha = gSigmaA*pow((STATE)LeftPOrder, (STATE)gL_alpha) / pow(faceSize, (STATE)gM_alpha);
245  STATE betta = gSigmaB*pow((STATE)LeftPOrder, (STATE)gL_betta) / pow(faceSize, (STATE)gM_betta);
246 
247  STATE u = dataleft.sol[0][0];
248  STATE gradu[2];
249  gradu[0] = dataleft.dsol[0](0,0);
250  gradu[1] = dataleft.dsol[0](1,0);
251 
252  STATE dualhp = dataleft.sol[0][2];//pq a dual hp está na terceira coluna
253  STATE GRADdualhp[2];//
254  STATE GRADLaplacdualhp[2];
255  GRADdualhp[0] = dataleft.dsol[0](0,2);
256  GRADdualhp[1] = dataleft.dsol[0](1,2);
257  STATE Laplacdualhp = dataleft.dsol[0](2,2);
258  GRADLaplacdualhp[0] = dataleft.dsol[0](3,2);
259  GRADLaplacdualhp[1] = dataleft.dsol[0](4,2);
260 
261  STATE dualhpPlus = dataleft.sol[0][3];//
262  STATE GRADdualhpPlus[2];//
263  STATE GRADLaplacdualhpPlus[2];//
264  GRADdualhpPlus[0]= dataleft.dsol[0](0,3);//
265  GRADdualhpPlus[1]= dataleft.dsol[0](1,3);//
266  STATE LaplacdualhpPlus = dataleft.dsol[0](2,3);
267  GRADLaplacdualhpPlus[0] = dataleft.dsol[0](3,3);
268  GRADLaplacdualhpPlus[1] = dataleft.dsol[0](4,3);
269 
270 
271  STATE gradZ[2];
272  gradZ[0] = GRADdualhp[0]; gradZ[1] = GRADdualhp[1];
273  STATE gradZPlus[2];
274  gradZPlus[0] = GRADdualhpPlus[0]; gradZPlus[1] = GRADdualhpPlus[1];
275 
276  STATE ValBC_D;
277  STATE ValBC_N;
278  ValBC_D = bc.Val2()(0,0);// Dirichlet
279  ValBC_N = bc.Val2()(1,0);// Neumann
280 
281  STATE ResD = ValBC_D - u;
282  STATE ResN = ValBC_N - (gradu[0]*normal[0]+gradu[1]*normal[1]);
283 
284  nk[0] += weight * ResD * ( (GRADLaplacdualhpPlus[0]*normal[0]+GRADLaplacdualhpPlus[1]*normal[1]) - (GRADLaplacdualhp[0]*normal[0]+GRADLaplacdualhp[1]*normal[1]) );
285  nk[0] -= weight*ResN*(LaplacdualhpPlus-Laplacdualhp);
286  nk[0] += weight*alpha*ResD*(dualhpPlus-dualhp);
287  nk[0] += weight*betta*ResN*( (GRADdualhpPlus[0]*normal[0]+GRADdualhpPlus[1]*normal[1]) - (GRADdualhp[0]*normal[0]+GRADdualhp[1]*normal[1]) );
288 
289 
290  nk[1] += weight * ResD * ( (du_exactd[3]*normal[0]+du_exactd[4]*normal[1]) - (GRADLaplacdualhp[0]*normal[0]+GRADLaplacdualhp[1]*normal[1]) );
291  nk[1] -= weight*ResN*(du_exactd[2]-Laplacdualhp);
292  nk[1] += weight*alpha*ResD*(u_exactd[0]-dualhp);
293  nk[1] += weight*betta*ResN*( (du_exactd[0]*normal[0]+du_exactd[1]*normal[1]) - (GRADdualhp[0]*normal[0]+GRADdualhp[1]*normal[1]) );
294 }
295 
297 STATE L2ErrorPrimal = 0.;
299 STATE L2ErrorDual = 0.;
300 
302  REAL weight,
303  TPZVec<REAL> &nk){
304  int phis = data.phi.Rows();
305  int dphir = data.dphix.Rows();
306  OrderSolution(data);
307  TPZFNMatrix<100,REAL> axes(3,3,0.),jacinv(2,2,0.);
308  TPZFNMatrix<200,STATE> ek(data.phi.Rows(),data.phi.Rows(),0.),ef(data.phi.Rows(),1,0.);
309  this->Contribute(data,weight,ek,ef);
310  data.phi.Resize(phis,1);
311  data.dphix.Resize(dphir,phis);
312  // std:: cout<< "x=("<< data.x[0]<<","<< data.x[1]<<")"<< " "<< "xsim=(" << 1.-data.x[0]<<","<< 1.-data.x[1]<< ")";
313  // ek.Print("Teste", std::cout, EFormatted);
314 
315  nk.Resize(10);
316 
317  // L2ErrorPrimal+= ((dphi(1,3)-dphi(1,1))*(dphi(1,3)-dphi(1,1))+(dphi(0,3)-dphi(0,1))*(dphi(0,3)-dphi(0,1))+(phi(3,0)-phi(1,0))*(phi(3,0)-phi(1,0)))*weight;
318  // L2ErrorDual += ((dphi(1,2)-dphi(1,0))*(dphi(1,2)-dphi(1,0))+(dphi(0,2)-dphi(0,0))*(dphi(0,2)-dphi(0,0))+(phi(2,0)-phi(0,0))*(phi(2,0)-phi(0,0)))*weight;
319 
320  nk[0] += ek(0,0);//norma da energia do erro aprox. dual
321  nk[3] += ek(1,1);// norma da energia do erro aprox. primal
322  nk[6] += ek(0,1);// estimador de erro (p e p+2)
323  nk[1] += ek(2,2);// norma da energia do erro dual exato
324  nk[4] += ek(3,3);// norma da energia do erro primal exato
325  nk[7] += ek(2,3);// estimador de erro (p e exato)
326  nk[2] += ek(4,4);// norma da energia do erro dual exato da sol. enriquecida
327  nk[5] += ek(5,5);// norma da energia do erro primal exato da sol. enriquecida
328  nk[8] += ek(4,5);// estimador de erro (p+2 e exato)
329  nk[9] += 0.; // indice de efetividade
330  // std:: cout<< " "<< nk[0]<<" "<< nk[1]<<" "<< nk[2]<< " "<< nk[3]<< " " <<nk[4]<< " "<< nk[5]<< std::endl;
331 }
332 
333 
334 
336  REAL weight,
337  TPZVec<STATE> &nkL,
338  TPZVec<STATE> &nkR){
339  // return;
340  OrderSolutionLeft(data,dataleft);
341  OrderSolutionRight(data,dataright);
342  int sz = dataleft.phi.Rows()+dataright.phi.Rows();
343  int ss = dataleft.phi.Rows(); //inclui isso pq a matriz eh szXsz mas os blocos são sz/2 X sz/2
344  TPZFNMatrix<100,REAL> axesleft(3,3,0.),axesright;
345  TPZFNMatrix<100,STATE> ef(sz,1,0.),ek(sz,sz,0.);
346  this->ContributeInterface(data,dataleft,dataright,weight,ek,ef);
347 
348  // std:: cout<< "x=("<< data.x[0]<<","<< data.x[1]<<")"; //" "<< "xsim=(" << 1.-data.x[0]<<","<< 1.-data.x[1]<< ")";
349  // ek.Print("Teste", std::cout, EFormatted);
350  nkL.Resize(10);
351  nkR.Resize(10);
352 
353  nkL[0] += 0.5*(ek(0,0)+ek(ss,0)+ek(0,ss)+ek(ss,ss)); // norma da energia do erro aprox. dual-Left
354  nkL[3] += 0.5*(ek(1,1)+ek(ss+1,1)+ek(1,1+ss)+ek(1+ss,1+ss)); // norma da energia do erro aprox. primal-Left
355  nkL[6] += 0.5*(ek(0,1)+ek(ss,ss+1)+ek(1,ss)+ek(ss+1,0)); // estimador (p e p+2)-Left // O erro está no segundo termo
356  nkL[1] += 0.5*(ek(2,2)+ek(2+ss,2+ss)+ek(2+ss,2)+ek(2,2+ss)); // norma da energia do erro exato da solucao dual-Left
357  nkL[4] += 0.5*(ek(3,3)+ek(3+ss,3+ss)+ek(3+ss,3)+ek(3,3+ss)); // norma da energia do erro exato da solucao primal-Left
358  nkL[7] += 0.5*(ek(2,3)+ek(ss+2,ss+3)+ek(3,2+ss)+ek(3+ss,2)); // estimador(p e exata)-Left
359  nkL[2] += 0.5*(ek(4,4)+ek(4+ss,4+ss)+ek(4+ss,4)+ek(4,4+ss));// norma da energia do erro dual-left exato da sol. enriquecida
360  nkL[5] += 0.5*(ek(5,5)+ek(ss+5,ss+5)+ek(5+ss,5)+ek(5,5+ss));// norma da energia do erro primal-left exato da sol. enriquecida
361  nkL[8] += 0.5*(ek(4,5)+ek(ss+4,ss+5)+ek(5,4+ss)+ek(5+ss,4)); // estimador(p+2 e exata)-Left
362  nkL[9] += 0.; // indíce de efetividade
363 
364  nkR = nkL;
365 
366  // nkR[0] += ek(ss,ss) +0.5*(ek(0,ss)+ek(ss,0)); // norma da energia do erro aprox. dual-Right
367  // nkR[3] += ek(1+ss,1+ss)+0.5*(ek(1,ss+1)+ek(ss+1,1));// norma da energia do erro aprox. dual-Right
368  // nkR[6] += ek(ss,ss+1)+(ek(ss+1,0) +ek(1,ss))/2.; // estimador(p e p+2)-Right // O erro está no segundo termo
369  // nkR[1] += ek(ss+2,ss+2)+0.5*(ek(2,2+ss)+ek(2+ss,2));// norma da energia do erro exato da sol. dual Right
370  // nkR[4] += ek(ss+3,ss+3)+0.5*(ek(3,3+ss)+ek(3+ss,3));// norma da energia do erro exato da sol. primal Right
371  // nkR[7] += ek(ss+2,ss+3)+(ek(3,2+ss)+ek(3+ss,2))/2.;// estimador (p e exato) - Right
372  // nkR[2] += ek(ss+4,ss+4)+0.5*(ek(4,4+ss)+ek(4+ss,4));// norma da energia do erro dual-right exato da sol. enriquecida
373  // nkR[5] += ek(ss+5,ss+5)+0.5*(ek(5,5+ss)+ek(5+ss,5));// norma da energia do erro primal-right exato da sol. enriquecida
374  // nkR[8] += ek(ss+4,ss+5)+(ek(5,4+ss)+ek(5+ss,4))/2.;// estimador (p+2 e exato) - Right
375  // nkR[9] += 0.; // indice de efetividade
376 }
377 
378 
379 
380 
382  REAL weight,
383  TPZVec<STATE> &nk,
384  TPZBndCond &bc){
385  // return;
386  OrderSolutionLeft(data,dataleft);
387 // OrderSolutionRight(data,dataright);
388  TPZFNMatrix<100,STATE> ek(dataleft.phi.Rows(),dataleft.phi.Rows(),0.),ef(dataleft.phi.Rows(),1,0.);
389 
390  this->ContributeBCInterface(data,dataleft,weight,ek,ef,bc);
391  // std:: cout<< "x=("<< data.x[0]<<","<< data.x[1]<<")";
392  // ek.Print("Teste", std::cout, EFormatted);
393  nk.Resize(10);
394 
395  nk[0] += ek(0,0);//norma da energia do erro aprox. dual
396  nk[1] += ek(1,1);// norma da energia do erro aprox.primal
397  nk[2] += ek(0,1);// estimador de erro (p e p+2)
398  nk[3] += ek(2,2);// norma da energia do erro dual exato
399  nk[4] += ek(3,3);// norma da energia do erro primal exato
400  nk[5] += ek(2,3);// estimador de erro (p e exato)
401  nk[6] += ek(4,4);// norma da energia do erro dual exato enriquecido
402  nk[7] += ek(5,5);// norma da energia do erro primal exato enriquecido
403  nk[8] += ek(4,5);// estimador de erro (p+2 e exato)
404  nk[9] += 0.; // indice de efetividade
405  // std:: cout<< " "<< nk[0]<<" "<< nk[1]<<" "<< nk[2]<< " "<< nk[3]<< " " <<nk[4]<< " "<< nk[5]<< std::endl;
406 }
407 
408 
409 
411  TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
412  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,
413  TPZVec<REAL> &values) {
414  TPZVec<STATE> sol(1), dsol(9,0.);
415 
416  Solution(u,dudx,axes,1,sol);
417  Solution(u,dudx,axes,2,dsol);
418  TPZVec<STATE> tmp_1(1,0.);
419  Psi(x, tmp_1);
420  values[0] = sol[0]*tmp_1[0];// parte da implentação do funcional, funcional com núcleo
421  values[1] = (sol[0] - u_exact[0])*(sol[0] - u_exact[0]); // error em norma L2
422  values[2] = 0.; // erro em semi norma H1
423  for(int id=0; id<2; id++) {
424  values[2] += (dsol[id] - du_exact(id,0))*(dsol[id] - du_exact(id,0));// erro em semi norma H1
425  }//(-grad(u).n)v
426  values[3] = values[1]+values[2]; // erro em norma H1
427 }
428 
429 
431  pisci[0]/*=8*exp(10*(2 - 2*x[0] + pow(x[0],2) - 2*x[1] + pow(x[1],2)))*
432  (1 - 46*x[1] + 1129*pow(x[1],2) - 120000*pow(x[0],7)*pow(-1 + x[1],2)*pow(x[1],2) + 20000*pow(x[0],8)*pow(-1 + x[1],2)*pow(x[1],2) -
433  2526*pow(x[1],3) + 2043*pow(x[1],4) - 800*pow(x[1],5) + 200*pow(x[1],6) +
434  x[0]*(-46 + 2116*x[1] - 31836*pow(x[1],2) + 76000*pow(x[1],3) - 73880*pow(x[1],4) + 36800*pow(x[1],5) - 9200*pow(x[1],6)) -
435  800*pow(x[0],5)*(1 - 46*x[1] + 996*pow(x[1],2) - 2260*pow(x[1],3) + 1910*pow(x[1],4) - 800*pow(x[1],5) + 200*pow(x[1],6)) +
436  200*pow(x[0],6)*(1 - 46*x[1] + 1986*pow(x[1],2) - 4240*pow(x[1],3) + 2900*pow(x[1],4) - 800*pow(x[1],5) + 200*pow(x[1],6)) +
437  pow(x[0],2)*(1129 - 31836*x[1] + 322916*pow(x[1],2) - 838160*pow(x[1],3) + 1045960*pow(x[1],4) - 796800*pow(x[1],5) + 397200*pow(x[1],6) -
438  120000*pow(x[1],7) + 20000*pow(x[1],8)) - 2*pow(x[0],3)*
439  (1263 - 38000*x[1] + 419080*pow(x[1],2) - 1066400*pow(x[1],3) + 1264600*pow(x[1],4) - 904000*pow(x[1],5) + 424000*pow(x[1],6) -
440  120000*pow(x[1],7) + 20000*pow(x[1],8)) + pow(x[0],4)*
441  (2043 - 73880*x[1] + 1045960*pow(x[1],2) - 2529200*pow(x[1],3) + 2604400*pow(x[1],4) - 1528000*pow(x[1],5) + 580000*pow(x[1],6) -
442  120000*pow(x[1],7) + 20000*pow(x[1],8)));*/
443 
444  // Experimento 1
445  =pow(2.,8)*( 24.*(1.-x[0])*(1.-x[0])*x[0]*x[0]
446  + 24.*(1.-x[1])*(1.-x[1])*x[1]*x[1]
447  + 2*(2*(1-x[0])*(1-x[0])-8*(1-x[0])*x[0]+2*x[0]*x[0])*(1-x[1])*(1-x[1])
448  - 8*(2*(1-x[0])*(1-x[0])-8*(1-x[0])*x[0]+2*x[0]*x[0])*(1-x[1])* x[1]
449  + 2* (2* (1 - x[0])* (1 - x[0]) - 8* (1 - x[0])* x[0] + 2* x[0]* x[0])* x[1]* x[1]
450  + (2*(1-x[0])*(1-x[0])-8*(1-x[0])*x[0]+2*x[0]*x[0])*(2* (1 - x[1])* (1 - x[1])
451  - 8* (1 - x[1])* x[1] + 2* x[1]* x[1]));//validacao
452 }// função núcleo do funcional
453 
455 {
456  int numbersol = data.dsol.size();
457  if (numbersol != 1) {
458  DebugStop();
459  }
460 
461  data.phi.Resize(6,1);
462  data.dphix.Resize(data.dsol[0].Rows(),6);
463  data.phi(0,0) = data.sol[0][3]- data.sol[0][2];//erro aprox. dual
464  data.phi(1,0) = data.sol[0][1]- data.sol[0][0];//erro aprox. primal
465  data.phi(2,0) = 0.;//erro exato dual p
466  data.phi(3,0) = 0.;//erro exato primal p
467  data.phi(4,0) = 0.;//erro exato dual p enriquecido
468  data.phi(5,0) = 0.;//erro exato primal p enriquecido
469 
470  int nd = data.dphix.Rows();
471  int id;
472  for(id=0; id<nd; id++)
473  {
474  data.dphix(id,0) = data.dsol[0](id,3)-data.dsol[0](id,2);// derivada do erro aprox dual
475  data.dphix(id,1) = data.dsol[0](id,1)-data.dsol[0](id,0);// derivada do erro aprox. primal
476  data.dphix(id,2) = 0. ;// derivada do erro exato dual
477  data.dphix(id,3) = 0. ;// derivada do erro exato primal
478  data.dphix(id,4) = 0. ;// derivada do erro exato dual enriquecido
479  data.dphix(id,5) = 0. ;// derivada do erro exato primal enriquecido
480  }
481  TPZVec<STATE> u_exactp(1);
482  TPZFMatrix<STATE> du_exactp(nd,1);
483  TPZVec<STATE> u_exactd(1);
484  TPZFMatrix<STATE> du_exactd(nd,1);
485  // erro exato primal
486  if (this->fPrimalExactSol )
487  { this->fPrimalExactSol(data.x,u_exactp,du_exactp);
488  data.phi(3,0)=u_exactp[0]-data.sol[0][0];
489  for(id=0; id<nd; id++)
490  { data.dphix(id,3) = du_exactp(id,0)- data.dsol[0](id,0);// derivada do erro primal exato
491  data.dphix(id,5) = du_exactp(id,0)- data.dsol[0](id,1);// derivada do erro primal exato enriquecido
492  }
493  }
494 
495  //erro exato dual
496  if (this->fDualExactSol )
497  { this->fDualExactSol(data.x,u_exactd,du_exactd);
498  data.phi(2,0)=u_exactd[0]-data.sol[0][2];
499  // std::cout<< "dual exato= "<< u_exactd[0]<< std::endl;
500  // std::cout<< "dual aprx p2= "<< data.sol[2]<< std::endl;
501  // std::cout<< "dual aprx q3= "<< data.sol[3]<< std::endl;
502  for(id=0; id<nd; id++)
503  { data.dphix(id,2) = du_exactd(id,0)- data.dsol[0](id,2);// derivada do erro dual exato
504  data.dphix(id,4) = du_exactd(id,0)- data.dsol[0](id,3);// derivada do erro dual exato enriquecido
505  }
506  }
507 }
508 
510 {
511  int numbersol = dataleft.dsol.size();
512  if (numbersol != 1) {
513  DebugStop();
514  }
515 
516  dataleft.phi.Resize(6,1);
517  dataleft.dphix.Resize(dataleft.dsol[0].Rows(),6);
518  dataleft.phi(0,0) = dataleft.sol[0][3]-dataleft.sol[0][2];//erro aprox. dual
519  dataleft.phi(1,0) = dataleft.sol[0][1]-dataleft.sol[0][0];//erro aprox. primal
520  dataleft.phi(2,0) = 0.;//erro exato dual
521  dataleft.phi(3,0) = 0.;//erro exato primal
522  dataleft.phi(4,0) = 0.;//erro exato dual enriquecido
523  dataleft.phi(5,0) = 0.;//erro exato primal enriquecido
524 
525  int nd = dataleft.dphix.Rows();
526  int id;
527  for(id=0; id<nd; id++)
528  {
529  dataleft.dphix(id,0) = dataleft.dsol[0](id,3)-dataleft.dsol[0](id,2);// derivada do erro aprox dual
530  dataleft.dphix(id,1) = dataleft.dsol[0](id,1)-dataleft.dsol[0](id,0);// derivada do erro aprox. primal
531  dataleft.dphix(id,2) = 0. ;// derivada do erro exato dual
532  dataleft.dphix(id,3) = 0. ;// derivada do erro exato primal
533  dataleft.dphix(id,4) = 0. ;// derivada do erro exato dual enriquecido
534  dataleft.dphix(id,5) = 0. ;// derivada do erro exato primal enriquecido
535  }
536 
537  TPZVec<STATE> u_exactp(1);
538  TPZFMatrix<STATE> du_exactp(nd,1);
539  TPZVec<STATE> u_exactd(1);
540  TPZFMatrix<STATE> du_exactd(nd,1);
541 
542  // erro exato primal
543  if (this->fPrimalExactSol )
544  {
545  this->fPrimalExactSol(data.x,u_exactp,du_exactp);
546  dataleft.phi(3,0)=u_exactp[0]-dataleft.sol[0][0];
547  dataleft.phi(5,0)=u_exactp[0]-dataleft.sol[0][1];
548  for(id=0; id<nd; id++)
549  { dataleft.dphix(id,3) = du_exactp(id,0)- dataleft.dsol[0](id,0);// derivada do erro primal exato
550  dataleft.dphix(id,5) = du_exactp(id,0)- dataleft.dsol[0](id,1);// derivada do erro primal exato enriquecido
551  }
552  }
553 
554  //erro exato dual
555  if (this->fDualExactSol )
556  {
557  this->fDualExactSol(data.x,u_exactd,du_exactd);
558  dataleft.phi(2,0)=u_exactd[0]-dataleft.sol[0][2];
559  dataleft.phi(4,0)=u_exactd[0]-dataleft.sol[0][3];
560  for(id=0; id<nd; id++)
561  { dataleft.dphix(id,2) = du_exactd(id,0)- dataleft.dsol[0](id,2);// derivada do erro dual exato
562  dataleft.dphix(id,4) = du_exactd(id,0)- dataleft.dsol[0](id,3);// derivada do erro dual exato enriquecido
563  }
564  }
565 }
566 
568 {
569  int numbersol = dataright.dsol.size();
570  if (numbersol != 1) {
571  DebugStop();
572  }
573 
574  dataright.phi.Resize(6,1);
575  dataright.dphix.Resize(dataright.dsol[0].Rows(),6);
576  dataright.phi(0,0) = dataright.sol[0][3]-dataright.sol[0][2];//erro aprox. dual
577  dataright.phi(1,0) = dataright.sol[0][1]-dataright.sol[0][0];//erro aprox. primal
578  dataright.phi(2,0) = 0.;//erro exato dual
579  dataright.phi(3,0) = 0.;//erro exato primal
580  dataright.phi(4,0) = 0.;//erro exato dual enriquecido
581  dataright.phi(5,0) = 0.;//erro exato primal enriquecido
582 
583  int nd = dataright.dphix.Rows();
584  int id;
585  for(id=0; id<nd; id++)
586  {
587  dataright.dphix(id,0) = dataright.dsol[0](id,3)-dataright.dsol[0](id,2);// derivada do erro aprox dual
588  dataright.dphix(id,1) = dataright.dsol[0](id,1)-dataright.dsol[0](id,0);// derivada do erro aprox. primal
589  dataright.dphix(id,2) = 0. ;// derivada do erro exato dual
590  dataright.dphix(id,3) = 0. ;// derivada do erro exato primal
591  dataright.dphix(id,4) = 0. ;// derivada do erro exato dual enriquecido
592  dataright.dphix(id,5) = 0. ;// derivada do erro exato primal enriquecido
593 
594  }
595 
596  TPZVec<STATE> u_exactp(1);
597  TPZFMatrix<STATE> du_exactp(8,1);
598  TPZVec<STATE> u_exactd(1);
599  TPZFMatrix<STATE> du_exactd(8,1);
600 
601  // erro exato primal
602  if (this->fPrimalExactSol )
603  { this->fPrimalExactSol(data.x,u_exactp,du_exactp);
604  dataright.phi(3,0)=u_exactp[0]-dataright.sol[0][0];
605  dataright.phi(5,0)=u_exactp[0]-dataright.sol[0][1];
606  for(id=0; id<nd; id++)
607  { dataright.dphix(id,3) = du_exactp(id,0)- dataright.dsol[0](id,0);// derivada do erro primal exato
608  dataright.dphix(id,5) = du_exactp(id,0)- dataright.dsol[0](id,1);// derivada do erro primal exato enriquecido
609  }
610  }
611 
612  //erro exato dual
613  if (this->fDualExactSol )
614  { this->fDualExactSol(data.x,u_exactd,du_exactd);
615  dataright.phi(2,0)=u_exactd[0]-dataright.sol[0][2];
616  dataright.phi(4,0)=u_exactd[0]-dataright.sol[0][3];
617  for(id=0; id<nd; id++)
618  { dataright.dphix(id,2) = du_exactd(id,0)- dataright.dsol[0](id,2);// derivada do erro dual exato
619  dataright.dphix(id,4) = du_exactd(id,0)- dataright.dsol[0](id,3);// derivada do erro dual exato enriquecido
620  }
621  }
622 
623 }
624 
626  return Hash("TPZBiharmonicEstimator") ^ TPZBiharmonic::ClassId() << 1;
627 }
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
virtual void ContributeInterfaceBCErrorsSimple(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZVec< STATE > &nk, TPZBndCond &bc)
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...
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
TPZBiharmonicEstimator(int nummat, STATE f)
Constructor.
Contains the TPZBiharmonic class which implements a discontinuous Galerkin formulation for the bi-har...
int ClassId() const override
Unique identifier for serialization purposes.
clarg::argBool bc("-bc", "binary checkpoints", false)
void OrderSolution(TPZMaterialData &data)
Computes the primal and dual exact error.
void Psi(TPZVec< REAL > &x, TPZVec< STATE > &pisci)
Kernel of the functional.
virtual int ClassId() const override
Unique identifier for serialization purposes.
std::ofstream CE("debugContributeErros.txt")
Output file to contribute erros.
static REAL gSigmaA
Definition: pzbiharmonic.h:26
static REAL gM_alpha
Definition: pzbiharmonic.h:26
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
static REAL gM_betta
Definition: pzbiharmonic.h:26
Implements discontinuous Galerkin formulation for the bi-harmonic equation.
Definition: pzbiharmonic.h:19
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
static REAL gL_alpha
Definition: pzbiharmonic.h:26
static REAL gSigmaB
Definition: pzbiharmonic.h:26
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Implements integral over element&#39;s volume.
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
int p
maximum polinomial order of the shape functions
virtual void ContributeErrorsSimple(TPZMaterialData &data, REAL weight, TPZVec< REAL > &nk)
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
void SetExactSolutions(void(*fp)(TPZVec< REAL > &loc, TPZVec< STATE > &val, TPZFMatrix< STATE > &deriv), void(*fd)(TPZVec< REAL > &locdual, TPZVec< STATE > &valdual, TPZFMatrix< STATE > &derivdual))
Set the pointer of the solution function.
void(* fPrimalExactSol)(TPZVec< REAL > &loc, TPZVec< STATE > &val, TPZFMatrix< STATE > &deriv)
Attributes required for goal oriented error estimation validation.
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
STATE L2ErrorDual
Initializing variable for computing L2 error to dual form.
virtual void ContributeInterfaceErrorsDual(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZVec< STATE > &nkL, TPZVec< STATE > &nkR)
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
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
f
Definition: test.py:287
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
virtual void ContributeInterfaceBCErrorsDual(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZVec< STATE > &nk, TPZBndCond &bc)
Implements integration of the boundary interface part of an error estimator.
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
void OrderSolutionLeft(TPZMaterialData &data, TPZMaterialData &dataleft)
Computes the primal and dual exact error over dataleft considering the solution in data...
REAL HSize
measure of the size of the element
virtual void ContributeInterfaceErrorsSimple(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZVec< STATE > &nkL, TPZVec< STATE > &nkR)
virtual void ContributeErrorsDual(TPZMaterialData &data, REAL weight, TPZVec< REAL > &nk)
Implements integration of the internal part of an error estimator.
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...
STATE L2ErrorPrimal
Initializing variable for computing L2 error.
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
void(* fDualExactSol)(TPZVec< REAL > &loc, TPZVec< STATE > &val, TPZFMatrix< STATE > &deriv)
void Errors(TPZVec< REAL > &x, TPZVec< STATE > &u, TPZFMatrix< STATE > &dudx, TPZFMatrix< REAL > &axes, TPZVec< STATE > &, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &values) override
Compute the error due to the difference between the interpolated flux and the flux computed based on...
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
static REAL gL_betta
Definition: pzbiharmonic.h:26
void OrderSolutionRight(TPZMaterialData &data, TPZMaterialData &dataright)
Computes the primal and dual exact error over dataright considering the solution in data...
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
Estimates error to biharmonic problem. Also computes the contributions on elements and interfaces...
def values
Definition: rdt.py:119
TPZSolVec sol
vector of the solutions at the integration point
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
Contains TPZBiharmonicEstimator class estimates error to biharmonic problem.