NeoPZ
pzbiharmonic.cpp
Go to the documentation of this file.
1 
6 #include "pzbiharmonic.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 
14 #include <cmath>
15 
16 // NIPG SIPG SSIPG1 SSIPG2
17 REAL TPZBiharmonic::gLambda1 = 1.0; // -1 1 -1 1
18 REAL TPZBiharmonic::gLambda2 = 1.0; // -1 1 1 -1
19 REAL TPZBiharmonic::gSigmaA = 10.0;// 10 10 10 10
20 REAL TPZBiharmonic::gSigmaB = 10.0;// 10 10 10 10
21 REAL TPZBiharmonic::gL_alpha = 6.0; // [0, 6]
22 REAL TPZBiharmonic::gM_alpha = 3.0; // IGUAL
23 REAL TPZBiharmonic::gL_betta = 4.0; // [-2, 4]
24 REAL TPZBiharmonic::gM_betta = 1.0; // IGUAL
25 using namespace std;
26 
27 
28 TPZBiharmonic::TPZBiharmonic(int nummat, REAL f)
30 fXf(f){}
31 
33 }
34 
35 void TPZBiharmonic::Print(std::ostream &out) {
36  out << "name of material : " << Name() << "\n";
37  out << "properties : \n";
38  TPZMaterial::Print(out);
39 
40 }
41 
43  REAL weight,
45  TPZFMatrix<STATE> &ef) {
46  TPZFMatrix<REAL> &dphi = data.dphix;
47  TPZFMatrix<REAL> &phi = data.phi;
48  TPZManVector<REAL,3> &x = data.x;
49 
50  int phr = phi.Rows();
51 
52  if(fForcingFunction) { // phi(in, 0) = phi_in
54  TPZFMatrix<STATE> grad;
55  fForcingFunction->Execute(x,res,grad); // dphi(i,j) = dphi_j/dxi
56  fXf = res[0];
57  }
58  //Equa�o de Poisson
59  for( int in = 0; in < phr; in++ ) {
60  ef(in, 0) += weight * fXf*phi(in,0);
61  for( int jn = 0; jn < phr; jn++ ) {
62  ek(in,jn) += weight * ( dphi(2,in) * dphi(2,jn) );
63  }
64  }
65 }
66 
67 
69  REAL weight,
72  TPZBndCond &bc) {
73 
74  //NOT TO BE DONE HERE
75  PZError << "TPZBiHarminic::ContributeBC - It should never be called.";
76 }
77 
79 int TPZBiharmonic::VariableIndex(const std::string &name){
80  if(!strcmp("Displacement6",name.c_str())) return 0;
81  if(!strcmp("Solution",name.c_str())) return 1;
82  if(!strcmp("Derivate",name.c_str())) return 2;
83  if(!strcmp("POrder",name.c_str())) return 10;
84  cout << "TPZBiharmonic::VariableIndex Error\n";
85  return -1;
86 }
87 
89  if(var == 0) return 6;
90  if(var == 1) return 1;
91  if(var == 2) return 2;
92  if(var == 10) return 1;
94  // cout << "TPZBiharmonic::NSolutionVariables Error\n";
95  // return 0;
96 }
97 
99  int var,TPZVec<STATE> &Solout){
100  if(var == 0 || var == 1) Solout[0] = Sol[0];//function
101  if(var == 2) {
102  Solout.Resize(DSol.Rows());
103  int id;
104  for(id=0 ; id < DSol.Rows(); id++) {
105  Solout[id] = DSol(id,0);//derivate
106  }
107  }
108 }
109 
111  TPZFMatrix<STATE> &/*DSol*/, TPZFMatrix<REAL> &/*axes*/,
112  TPZVec<STATE> &/*flux*/) {
113  //Flux(TPZVec<REAL> &x, TPZVec<REAL> &Sol, TPZFMatrix<REAL> &DSol, TPZFMatrix<REAL> &axes, TPZVec<REAL> &flux)
114 }
115 
117  TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
118  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,
119  TPZVec<REAL> &values) {
120  TPZVec<STATE> sol(1), dsol(8,0.);
121  Solution(u,dudx,axes,1,sol);
122  Solution(u,dudx,axes,2,dsol);
123  //values[1] : error em norma L2
124  values[1] = (sol[0] - u_exact[0])*(sol[0] - u_exact[0]);
125 
126  //values[2] : erro em semi norma H1
127  values[2] = 0.;
128  for(int id=0; id<2; id++) {
129  values[2] += (dsol[id] - du_exact(id,0))*(dsol[id] - du_exact(id,0));
130  }
131  //values[3] : erro em semi norma Laplace
132  values[3] = (dsol[2] - du_exact(2,0))*(dsol[2] - du_exact(2,0));
133 
134  //values[0] : erro em norma H1
135  values[0] = values[1]+values[2];
136  // dxx
137  values[5] = (dsol[5] - du_exact(5,0))*(dsol[5] - du_exact(5,0));
138  // dyy
139  values[6] = (dsol[6] - du_exact(6,0))*(dsol[6] - du_exact(6,0));
140  // dxy
141  values[7] = (dsol[7] - du_exact(7,0))*(dsol[7] - du_exact(7,0));
142  //values[4] : erro em norma H2
143  values[4] = values[5]+values[6]+values[7]+values[0];
144 }
145 
147  REAL weight,
148  TPZFMatrix<STATE> &ek,
149  TPZFMatrix<STATE> &ef){
150  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
151  TPZFMatrix<REAL> &dphiR = dataright.dphix;
152  TPZFMatrix<REAL> &phiL = dataleft.phi;
153  TPZFMatrix<REAL> &phiR = dataright.phi;
154  TPZManVector<REAL,3> &normal = data.normal;
155  int LeftPOrder=dataleft.p;
156  int RightPOrder=dataright.p;
157  REAL faceSize=data.HSize;
158 
159  int nrowl = phiL.Rows();
160  int nrowr = phiR.Rows();
161  int il,jl,ir,jr,id;
162 
163  REAL alpha=gSigmaA*(pow(((REAL)LeftPOrder),gL_alpha)+pow(((REAL)RightPOrder), gL_alpha) ) /
164  (2. * pow(faceSize, gM_alpha) );
165 
166  REAL betta=gSigmaB*(pow(((REAL)LeftPOrder),gL_betta)+pow(((REAL)RightPOrder), gL_betta) ) /
167  (2. * pow(faceSize, gM_betta) );
168 
169  /* Primeira Integral */
170  for(il=0; il<nrowl; il++) {
171  REAL dphiLinormal = 0.;
172  for(id=0; id<2; id++) {
173  dphiLinormal += dphiL(3+id,il)*normal[id];
174  }
175  for(jl=0; jl<nrowl; jl++) {
176  REAL dphiLjnormal = 0.;
177  for(id=0; id<2; id++) {
178  dphiLjnormal += dphiL(3+id,jl)*normal[id];
179  }
180 
181  ek(il,jl) += weight*0.5*(+gLambda1*dphiLinormal*phiL(jl,0) + dphiLjnormal*phiL(il,0));
182  }
183  }
184  for(ir=0; ir<nrowr; ir++) {
185  REAL dphiRinormal = 0.;
186  for(id=0; id<2; id++) {
187  dphiRinormal += dphiR(3+id,ir)*normal[id];
188  }
189  for(jr=0; jr<nrowr; jr++) {
190  REAL dphiRjnormal = 0.;
191  for(id=0; id<2; id++) {
192  dphiRjnormal += dphiR(3+id,jr)*normal[id];
193  }
194  ek(ir+nrowl,jr+nrowl) += weight*0.5*(
195  -gLambda1*dphiRinormal*phiR(jr,0) - dphiRjnormal*phiR(ir,0)
196  );
197  }
198  }
199  for(il=0; il<nrowl; il++) {
200  REAL dphiLinormal = 0.;
201  for(id=0; id<2; id++) {
202  dphiLinormal += dphiL(3+id,il)*normal[id];
203  }
204  for(jr=0; jr<nrowr; jr++) {
205  REAL dphiRjnormal = 0.;
206  for(id=0; id<2; id++) {
207  dphiRjnormal += dphiR(3+id,jr)*normal[id];
208  }
209  ek(il,jr+nrowl) += weight*0.5*(
210  - gLambda1*dphiLinormal*phiR(jr,0) + dphiRjnormal*phiL(il,0)
211  );
212  }
213  }
214  for(ir=0; ir<nrowr; ir++) {
215  REAL dphiRinormal = 0.;
216  for(id=0; id<2; id++) {
217  dphiRinormal += dphiR(3+id,ir)*normal[id];
218  }
219  for(jl=0; jl<nrowl; jl++) {
220  REAL dphiLjnormal = 0.;
221  for(id=0; id<2; id++) {
222  dphiLjnormal += dphiL(3+id,jl)*normal[id];
223  }
224  ek(ir+nrowl,jl) += weight*0.5*(
225  + gLambda1*dphiRinormal*phiL(jl,0) - dphiLjnormal*phiR(ir,0)
226  );
227  }
228  }
229 
230  /* Segunda Integral */
231  for(il=0; il<nrowl; il++) {
232  REAL dphiLinormal = 0.;
233  for(id=0; id<2; id++) {
234  dphiLinormal += dphiL(id,il)*normal[id];
235  }
236  for(jl=0; jl<nrowl; jl++) {
237  REAL dphiLjnormal = 0.;
238  for(id=0; id<2; id++) {
239  dphiLjnormal += dphiL(id,jl)*normal[id];
240  }
241 
242  ek(il,jl) += weight*0.5*(
243  - dphiLinormal*dphiL(2,jl) - gLambda2*dphiLjnormal*dphiL(2,il)
244  );
245  }
246  }
247  for(ir=0; ir<nrowr; ir++) {
248  REAL dphiRinormal = 0.;
249  for(id=0; id<2; id++) {
250  dphiRinormal += dphiR(id,ir)*normal[id];
251  }
252  for(jr=0; jr<nrowr; jr++) {
253  REAL dphiRjnormal = 0.;
254  for(id=0; id<2; id++) {
255  dphiRjnormal += dphiR(id,jr)*normal[id];
256  }
257  ek(ir+nrowl,jr+nrowl) += weight*0.5*(
258  + dphiRinormal*dphiR(2,jr) + gLambda2*dphiRjnormal*dphiR(2,ir)
259  );
260  }
261  }
262  for(il=0; il<nrowl; il++) {
263  REAL dphiLinormal = 0.;
264  for(id=0; id<2; id++) {
265  dphiLinormal += dphiL(id,il)*normal[id];
266  }
267  for(jr=0; jr<nrowr; jr++) {
268  REAL dphiRjnormal = 0.;
269  for(id=0; id<2; id++) {
270  dphiRjnormal += dphiR(id,jr)*normal[id];
271  }
272  ek(il,jr+nrowl) += weight*0.5*(
273  - dphiLinormal*dphiR(2,jr) + gLambda2*dphiRjnormal*dphiL(2,il)
274  );
275  }
276  }
277  for(ir=0; ir<nrowr; ir++) {
278  REAL dphiRinormal = 0.;
279  for(id=0; id<2; id++) {
280  dphiRinormal += dphiR(id,ir)*normal[id];
281  }
282  for(jl=0; jl<nrowl; jl++) {
283  REAL dphiLjnormal = 0.;
284  for(id=0; id<2; id++) {
285  dphiLjnormal += dphiL(id,jl)*normal[id];
286  }
287  ek(ir+nrowl,jl) += weight*0.5*(
288  + dphiRinormal*dphiL(2,jl) - gLambda2*dphiLjnormal*dphiR(2,ir)
289  );
290  }
291  }
292 
293  /* Terceira Integral */
294  for(il=0; il<nrowl; il++) {
295  REAL dphiLinormal = 0.;
296  for(id=0; id<2; id++) {
297  dphiLinormal += dphiL(id,il)*normal[id];
298  }
299  for(jl=0; jl<nrowl; jl++) {
300  REAL dphiLjnormal = 0.;
301  for(id=0; id<2; id++) {
302  dphiLjnormal += dphiL(id,jl)*normal[id];
303  }
304 
305  ek(il,jl) += weight*(
306  alpha * phiL(jl,0)*phiL(il,0) +
307  betta * dphiLinormal*dphiLjnormal
308  );
309  }
310  }
311  for(ir=0; ir<nrowr; ir++) {
312  REAL dphiRinormal = 0.;
313  for(id=0; id<2; id++) {
314  dphiRinormal += dphiR(id,ir)*normal[id];
315  }
316  for(jr=0; jr<nrowr; jr++) {
317  REAL dphiRjnormal = 0.;
318  for(id=0; id<2; id++) {
319  dphiRjnormal += dphiR(id,jr)*normal[id];
320  }
321  ek(ir+nrowl,jr+nrowl) += weight*(
322  alpha * phiR(jr,0)*phiR(ir,0) +
323  betta * dphiRinormal*dphiRjnormal
324  );
325  }
326  }
327  for(il=0; il<nrowl; il++) {
328  REAL dphiLinormal = 0.;
329  for(id=0; id<2; id++) {
330  dphiLinormal += dphiL(id,il)*normal[id];
331  }
332  for(jr=0; jr<nrowr; jr++) {
333  REAL dphiRjnormal = 0.;
334  for(id=0; id<2; id++) {
335  dphiRjnormal += dphiR(id,jr)*normal[id];
336  }
337  ek(il,jr+nrowl) += weight*(
338  - alpha * phiR(jr,0)*phiL(il,0)
339  - betta * dphiLinormal*dphiRjnormal
340  );
341  }
342  }
343  for(ir=0; ir<nrowr; ir++) {
344  REAL dphiRinormal = 0.;
345  for(id=0; id<2; id++) {
346  dphiRinormal += dphiR(id,ir)*normal[id];
347  }
348  for(jl=0; jl<nrowl; jl++) {
349  REAL dphiLjnormal = 0.;
350  for(id=0; id<2; id++) {
351  dphiLjnormal += dphiL(id,jl)*normal[id];
352  }
353  ek(ir+nrowl,jl) += weight*(
354  - alpha * phiL(jl,0)*phiR(ir,0)
355  - betta * dphiRinormal*dphiLjnormal
356  );
357  }
358  }
359 }
360 
362  REAL weight,
363  TPZFMatrix<STATE> &ek,
364  TPZFMatrix<STATE> &ef,
365  TPZBndCond &bc) {
366 
367  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
368  TPZFMatrix<REAL> &phiL = dataleft.phi;
369  TPZManVector<REAL,3> &normal = data.normal;
370  int POrder=data.p; // I need some explains, why you use reference & - Jorge
371  REAL faceSize=data.HSize;
372 
373  REAL alpha = gSigmaA*pow(((REAL)POrder), gL_alpha) / pow(faceSize, gM_alpha);
374  REAL betta = gSigmaB*pow(((REAL)POrder), gL_betta) / pow(faceSize, gM_betta);
375 
376  int il,jl,nrowl,id;
377  nrowl = phiL.Rows();
378 
379  /* Primeira Integral */
380  for(il=0; il<nrowl; il++) {
381 
382  REAL dphiLinormal = 0.;
383  for(id=0; id<2; id++) {
384  dphiLinormal += dphiL(3+id,il)*normal[id];
385  }
386 
387  // Termos de Dirichlet - em Val2()(0,0)
388  ef(il,0) += + gLambda1*weight*dphiLinormal*bc.Val2()(0,0)
389  + alpha * weight*bc.Val2()(0,0)*phiL(il) ;
390 
391  for(jl=0; jl<nrowl; jl++) {
392  REAL dphiLjnormal = 0.;
393  for(id=0; id<2; id++) {
394  dphiLjnormal += dphiL(3+id,jl)*normal[id];
395  }
396  ek(il,jl) += weight*(+ gLambda1*dphiLinormal*phiL(jl,0)+ dphiLjnormal*phiL(il,0)); //2 1
397  }
398  }
399 
400  /* Segunda Integral */
401  for(il=0; il<nrowl; il++) {
402  REAL dphiLinormal = 0.;
403  for(id=0; id<2; id++) {
404  dphiLinormal += dphiL(id,il)*normal[id];
405  }
406 
407  // Termos de Neuwmann - em Val2()(1,0)
408  ef(il,0) += - gLambda2*weight*bc.Val2()(1,0)*dphiL(2,il)
409  + betta * weight*bc.Val2()(1,0)*dphiLinormal ;
410 
411  for(jl=0; jl<nrowl; jl++) {
412  REAL dphiLjnormal = 0.;
413  for(id=0; id<2; id++) {
414  dphiLjnormal += dphiL(id,jl)*normal[id];
415  }
416  ek(il,jl) += weight*(- dphiLinormal*dphiL(2,jl) - gLambda2*dphiLjnormal*dphiL(2,il) );
417  }
418  }
419 
420 
421  for(il=0; il<nrowl; il++) {
422  REAL dphiLinormal = 0.;
423  for(id=0; id<2; id++) {
424  dphiLinormal += dphiL(id,il)*normal[id];
425  }
426  for(jl=0; jl<nrowl; jl++) {
427  REAL dphiLjnormal = 0.;
428  for(id=0; id<2; id++) {
429  dphiLjnormal += dphiL(id,jl)*normal[id];
430  }
431  ek(il,jl) += weight*(
432  alpha * phiL(jl,0)*phiL(il,0) +
433  betta * dphiLinormal*dphiLjnormal
434  );
435  }
436  }
437 }
438 
440  return Hash("TPZBiharmonic") ^ TPZDiscontinuousGalerkin::ClassId() << 1;
441 }
Defines the interface which material objects need to implement for discontinuous Galerkin formulation...
Definition: pzdiscgal.h:20
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
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...
int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzdiscgal.cpp:110
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
static REAL gLambda1
Definition: pzbiharmonic.h:26
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Contains the TPZBiharmonic class which implements a discontinuous Galerkin formulation for the bi-har...
clarg::argBool bc("-bc", "binary checkpoints", false)
virtual int ClassId() const override
Unique identifier for serialization purposes.
static REAL gLambda2
Definition: pzbiharmonic.h:26
static REAL gSigmaA
Definition: pzbiharmonic.h:26
static REAL gM_alpha
Definition: pzbiharmonic.h:26
Defines PZError.
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
static REAL gM_betta
Definition: pzbiharmonic.h:26
Implements discontinuous Galerkin formulation for the bi-harmonic equation.
Definition: pzbiharmonic.h:19
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
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
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
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Implements boundary conditions for continuous Galerkin.
Contains TPZMatrixclass which implements full matrix (using column major representation).
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
string res
Definition: test.py:151
REAL HSize
measure of the size of the element
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...
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
Contains TPZMatrix<TVar>class, root matrix class.
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
Compute the error due to the difference between the interpolated flux and the flux computed based on ...
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.
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
virtual ~TPZBiharmonic()
Destructor.
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
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
def values
Definition: rdt.py:119
TPZBiharmonic(int nummat, REAL f)
Inicialisation of biharmonic material.
virtual int VariableIndex(const std::string &name) override
virtual std::string Name() override
Returns the name of the material.
Definition: pzbiharmonic.h:49
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15