NeoPZ
TPZConsLawTest.cpp
Go to the documentation of this file.
1 
5 #include "TPZConsLawTest.h"
6 #include "pzelmat.h"
7 #include "pzbndcond.h"
8 #include "pzmatrix.h"
9 #include "pzfmatrix.h"
10 #include "pzerror.h"
11 #include <math.h>
12 #include <cmath>
13 
14 using namespace std;
15 
16 
17 TPZConsLawTest::TPZConsLawTest(int nummat, TPZVec<STATE> B,int artdiff,STATE delta_t,int dim,STATE delta,int test) : TPZConservationLaw(nummat,delta_t,dim), fXf(1,1,0.), fB(dim) {
18 
19  if(artdiff<0 || artdiff>3){
20  PZError << "TPZConsLawTest::TPZConsLawTest artificial diffusion parameter, default 1\n";
21  artdiff = 1;
22  }
23  if(B.NElements() != dim){
24  PZError << "TPZConsLawTest::TPZConsLawTest error in dimension of B, default dimension " << dim << endl;
25  B.Resize(dim);
26  int i;
27  for(i=0;i<dim;i++) B[i] = 0.;
28  B[dim-1] = 1.0;
29  }
30  int i;
31  for(i=0;i<dim;i++) fB[i] = B[i];
32  fArtificialDiffusion = artdiff;//SUPG, SL or Bornhauss
33  fDelta = delta;
34  fTest = test;
35 }
36 
38 }
39 
41  return 1;
42 }
43 
44 void TPZConsLawTest::Print(std::ostream &out) {
45  out << "name of material : " << Name() << "\n";
46  out << "properties : \n";
47  TPZMaterial::Print(out);
48 }
49 
51  REAL weight,
53  TPZFMatrix<STATE> &ef) {
54 
55  TPZFMatrix<REAL> &dphi = data.dphix;
56  // TPZFMatrix<REAL> &dphiL = data.dphixl;
57  // TPZFMatrix<REAL> &dphiR = data.dphixr;
58  TPZFMatrix<REAL> &phi = data.phi;
59  // TPZFMatrix<REAL> &phiL = data.phil;
60  // TPZFMatrix<REAL> &phiR = data.phir;
61  // TPZManVector<REAL,3> &normal = data.normal;
62  TPZManVector<REAL,3> &x = data.x;
63  // int &POrder=data.p;
64  // int &LeftPOrder=data.leftp;
65  // int &RightPOrder=data.rightp;
66  int numbersol = data.sol.size();
67  if (numbersol != 1) {
68  DebugStop();
69  }
70 
71  TPZVec<STATE> &sol=data.sol[0];
72  // TPZVec<REAL> &solL=data.soll;
73  // TPZVec<REAL> &solR=data.solr;
74  // TPZFMatrix<REAL> &dsol=data.dsol;
75  // TPZFMatrix<REAL> &dsolL=data.dsoll;
76  // TPZFMatrix<REAL> &dsolR=data.dsolr;
77  // REAL &faceSize=data.HSize;
78  // TPZFMatrix<REAL> &daxesdksi=data.daxesdksi;
79  // TPZFMatrix<REAL> &axes=data.axes;
80 
81  int phr = phi.Rows();// phi(in, 0) = phi_in , dphi(i,jn) = dphi_jn/dxi
82 
83  if(fForcingFunction) {
85  fForcingFunction->Execute(x,res);//fXf(0,0) = res[0];//if(!sol[0])
86  sol[0] = res[0];//solu�o inicial, na itera�o 0 sol = 0 e res != 0
87  } //nas itera�es > 0 sol != 0 e res = 0
88  int dim = dphi.Rows();
89  if(Dimension() != dim)
90  PZError << "TPZConsLawTest::Contribute dimension error, dimension = " << dim;
91 
92  STATE sum1 = 0.,sum2=0.;
93  int i;
94  STATE time = TimeStep();
95  for( int in = 0; in < phr; in++ ) {
96  //primeira parcela
97  sum1 = phi(in, 0) * sol[0];
98  //segunda parcela
99  sum2 = 0.;
100  for(i=0;i<dim;i++) sum2 += B(i,x) * dphi(i,in);//F = (B0*u.B1*u,B2*u)
101  sum2 *= time * sol[0];
102  //terceira parcela: contribui� do elemento interface: feita na classe TPZInterfaceElement
103  //contribui� final
104  ef(in, 0) += weight *(sum1+sum2);//* fXf(0,0);
105 
106  for( int jn = 0; jn < phr; jn++ ) {
107  //primeira parcela
108  ek(in,jn) += weight * phi(in,0) * phi(jn,0);
109  for(int ki=0; ki<dim; ki++) {
110  for(int kj=0; kj<dim; kj++) {
111  //segunda parcela
112  ek(in,jn) += weight * time * Delta() * T(ki,x) * B(kj,x) * ( dphi(kj,in) * dphi(ki,jn) );
113  }
114  }
115  }
116  }
117 }
118 
120 
121  return fDelta;
122 }
123 
125 
126  //REAL cfl = CFL(1);
127  //REAL delta = ( (10./3.)*cfl*cfl - (2./3.)*cfl + 1./10. );
128  return 0.0;//delta;
129 }
130 
132 
133  return (1.0/(2.0*degree+1.0));
134 }
135 
137 
138  if(fTest==0) return fB[i];
139 
140  if(fTest==1){
141  if(i==0) return -x[1];
142  if(i==1) return x[0];
143  if(i==2) return 0.0;
144  }
145  if(fTest==2){
146  if(i==0) return -x[1];
147  if(i==1) return x[0];
148  if(i==2) return 0.;
149  }
150  return 0.;
151 }
152 
154 
155  //SUPG
156  if(fArtificialDiffusion == 1){
157  int i;
158  STATE norm = 0;
159  int dim = Dimension();
160  for(i=0;i<dim;i++) norm += B(i,x)*B(i,x);
161  norm = sqrt(norm);
162  if(jn==0) return B(0,x)/norm;
163  if(jn==1) return B(1,x)/norm;
164  if(jn==2) return B(2,x)/norm;
165  }
166  //LS
167  if(fArtificialDiffusion == 2){
168  cout << "TPZConsLawTest::T artificial diffusion LS not implemented\n";
169  return 0;
170  }
171  //Bornhaus
172  if(fArtificialDiffusion == 3){
173  cout << "TPZConsLawTest::T artificial diffusion Bornhaus not implemented\n";
174  return 0;
175  }
176  cout << "TPZConsLawTest::T case not implemented, case = " << jn << endl;
177  return 0.0;
178 }
179 
181  REAL weight,
182  TPZFMatrix<STATE> &ek,
183  TPZFMatrix<STATE> &ef){
184 
185  // TPZFMatrix<REAL> &dphi = data.dphix;
186  // TPZFMatrix<REAL> &dphiL = data.dphixl;
187  // TPZFMatrix<REAL> &dphiR = data.dphixr;
188  // TPZFMatrix<REAL> &phi = data.phi;
189  TPZFMatrix<REAL> &phiL = dataleft.phi;
190  TPZFMatrix<REAL> &phiR = dataright.phi;
191  TPZManVector<REAL,3> &normal = data.normal;
192  TPZManVector<REAL,3> &x = data.x;
193  // int &POrder=data.p;
194  // int &LeftPOrder=data.leftp;
195  // int &RightPOrder=data.rightp;
196  // TPZVec<REAL> &sol=data.sol;
197  int numbersol = dataleft.sol.size();
198  if (numbersol != 1) {
199  DebugStop();
200  }
201 
202  TPZVec<STATE> &solL=dataleft.sol[0];
203  TPZVec<STATE> &solR=dataright.sol[0];
204  // TPZFMatrix<REAL> &dsol=data.dsol;
205  // TPZFMatrix<REAL> &dsolL=data.dsoll;
206  // TPZFMatrix<REAL> &dsolR=data.dsolr;
207  // REAL &faceSize=data.HSize;
208  // TPZFMatrix<REAL> &daxesdksi=data.daxesdksi;
209  // TPZFMatrix<REAL> &axes=data.axes;
210 
211 
212  int phrl = phiL.Rows();
213  int phrr = phiR.Rows();
214 
215  if(fForcingFunction) { // phi(in, 0) = phi_in
216  TPZManVector<STATE> res(1);// dphi(i,j) = dphi_j/dxi
217  fForcingFunction->Execute(x,res);
218  if(phrl) solL[0] = res[0];//solu�o inicial interior (t=0), na itera�o 0 sol = 0 e res != 0
219  if(phrr) solR[0] = res[0];//nas itera�es > 0 sol != 0 e res = 0
220  }
221 
222  STATE Bn=0.0;
223  //a normal est�imersa no mesmo espao da malha: R, R, R
224  int dim = Dimension(),i;
225  for(i=0;i<dim;i++) Bn += B(i,x)*normal[i];
226  //contribui� do elemento interface
227  REAL time = TimeStep();
228  int efc = 0;
229  if(Bn > 0){
230  //a normal �interface sempre aponta do seu elemento esquerdo para o seu direito
231  for( int in = 0; in < phrl; in++ ) {
232  ef(efc , 0) += -time * weight * solL[0] * phiL(in, 0) * Bn;//0 2 4 ..
233  efc++;
234  }
235  for( int in = 0; in < phrr; in++ ) {
236  ef(efc, 0) += -time * weight * solL[0] * phiR(in, 0) * -Bn;//1 3 5 ..
237  efc++;
238  }
239  } else {
240  for( int in = 0; in < phrl; in++ ) {
241  ef(efc , 0) += -time * weight * solR[0] * phiL(in, 0) * Bn;//0 2 4 ..
242  efc++;
243  }
244  for( int in = 0; in < phrr; in++ ) {
245  ef(efc, 0) += -time * weight * solR[0] * phiR(in, 0) * -Bn;//1 3 5 ..
246  efc++;
247  }
248  }
249 }
250 
252  REAL weight,
253  TPZFMatrix<STATE> &ek,
254  TPZFMatrix<STATE> &ef,
255  TPZBndCond &bc) {
256 
257  // TPZFMatrix<REAL> &dphi = data.dphix;
258  // TPZFMatrix<REAL> &dphiL = data.dphixl;
259  // TPZFMatrix<REAL> &dphiR = data.dphixr;
260  TPZFMatrix<REAL> &phi = data.phi;
261  // TPZFMatrix<REAL> &phiL = data.phil;
262  // TPZFMatrix<REAL> &phiR = data.phir;
263  // TPZManVector<REAL,3> &normal = data.normal;
264  // TPZManVector<REAL,3> &x = data.x;
265  // int &POrder=data.p;
266  // int &LeftPOrder=data.leftp;
267  // int &RightPOrder=data.rightp;
268  // TPZVec<REAL> &sol=data.sol;
269  // TPZVec<REAL> &solL=data.soll;
270  // TPZVec<REAL> &solR=data.solr;
271  // TPZFMatrix<REAL> &dsol=data.dsol;
272  // TPZFMatrix<REAL> &dsolL=data.dsoll;
273  // TPZFMatrix<REAL> &dsolR=data.dsolr;
274  // REAL &faceSize=data.HSize;
275  // TPZFMatrix<REAL> &daxesdksi=data.daxesdksi;
276  // TPZFMatrix<REAL> &axes=data.axes;
277 
278  int phr = phi.Rows();
279  short in,jn;
280  STATE v2[1];
281  v2[0] = bc.Val2()(0,0);
282 
283  switch (bc.Type()) {
284  case 0 : // Dirichlet condition
285  for(in = 0 ; in < phr; in++) {
286  ef(in,0) += gBigNumber * v2[0] * phi(in,0) * weight;
287  for (jn = 0 ; jn < phr; jn++) {
288  ek(in,jn) += gBigNumber * phi(in,0) * phi(jn,0) * weight;
289  }
290  }
291  break;
292  case 1 : // Neumann condition
293  for(in = 0 ; in < phi.Rows(); in++) {
294  ef(in,0) += v2[0] * phi(in,0) * weight;
295  }
296  break;
297  case 2 : // condi�o mista
298  for(in = 0 ; in < phi.Rows(); in++) {
299  ef(in, 0) += v2[0] * phi(in, 0) * weight;
300  for (jn = 0 ; jn < phi.Rows(); jn++) {
301  ek(in,jn) += bc.Val1()(0,0) * phi(in,0) *
302  phi(jn,0) * weight; // peso de contorno => integral de contorno
303  }
304  }
305  }
306 }
307 
309 int TPZConsLawTest::VariableIndex(const std::string &name){
310  if(!strcmp("Solution",name.c_str())) return 1;
311  if(!strcmp("Derivate",name.c_str())) return 2;
312  cout << "TPZConsLawTest::VariableIndex Error\n";
313  return -1;
314 }
315 
317 
318  if(var == 1) return 1;
319  if(var == 2) return Dimension();
320  cout << "TPZConsLawTest::NSolutionVariables Error\n";
321  return 0;
322 }
323 
325 
326  if(var == 0 || var == 1) Solout[0] = Sol[0];//function
327  if(var == 2) {
328  Solout[0] = DSol(0,0);//derivate
329  Solout[1] = DSol(1,0);//derivate
330  Solout[2] = DSol(2,0);//derivate
331  }
332 }
333 
334 void TPZConsLawTest::Flux(TPZVec<REAL> &/*x*/, TPZVec<STATE> &/*Sol*/, TPZFMatrix<STATE> &/*DSol*/, TPZFMatrix<REAL> &/*axes*/, TPZVec<STATE> &/*flux*/) {
335  //Flux(TPZVec<REAL> &x, TPZVec<REAL> &Sol, TPZFMatrix<REAL> &DSol, TPZFMatrix<REAL> &axes, TPZVec<REAL> &flux)
336 }
337 
339  TPZFMatrix<STATE> &dudx, TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
340  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) {
341 
342  TPZVec<STATE> sol(1),dsol(3);
343  Solution(u,dudx,axes,1,sol);
344  Solution(u,dudx,axes,2,dsol);
345  REAL dx = dsol[0]*axes(0,0)+dsol[1]*axes(1,0)+dsol[2]*axes(2,0);
346  REAL dy = dsol[0]*axes(0,1)+dsol[1]*axes(1,1)+dsol[2]*axes(2,1);
347  REAL dz = dsol[0]*axes(0,2)+dsol[1]*axes(1,2)+dsol[2]*axes(2,2);
348  //values[1] : eror em norma L2
349  values[1] = pow(sol[0] - u_exact[0],(STATE)2.0);
350  //values[2] : erro em semi norma H1
351  values[2] = (dx - du_exact(0,0))*(dx - du_exact(0,0));
352  values[2] += ((dy - du_exact(1,0))*(dy - du_exact(1,0)));
353  values[2] += ((dz - du_exact(2,0))*dz - du_exact(2,0));
354  //values[0] : erro em norma H1 <=> norma Energia
355  values[0] = values[1]+values[2];
356 }
357 
359 
360  if(!bcleft){
361  PZError << "TPZConsLawTest::ComputeSolLeft null bundary condition return\n";
362  return;
363  }
364  //int nstate = NStateVariables();
365  TPZFMatrix<REAL> jacinv(0,0),axes(0,0);
366  switch (bcleft->Type()){
367  case 0://Dirichlet
368  case 1://Neumann
369  case 2://Mista
370  PZError << "TPZConsLawTest::ComputeSolLeft boundary condition error\n";
371  break;
372  case 3://Dirichlet: nada a fazer a CC �a correta
373  break;
374  case 4://recuperar valor da solu� MEF direita: saida livre
375  soll[0] = solr[0];
376  break;
377  default:
378  PZError << "TPZConsLawTest::ContributeInterface Boundary Condition Type Not Exists\n";
379  }
380 }
381 
382 
384 
385  if(!bcright){
386  PZError << "TPZConsLawTest::ComputeSolLeft null bundary condition return\n";
387  return;
388  }
389  //int nstate = NStateVariables();
390  TPZFMatrix<REAL> jacinv(0,0),axes(0,0);
391  switch (bcright->Type()){
392  case 0://Dirichlet
393  case 1://Neumann
394  case 2://Mista
395  PZError << "TPZConsLawTest::ComputeSolLeft boundary condition error\n";
396  break;
397  case 3://Dirichlet: nada a fazer a CC �a correta
398  break;
399  case 4://recuperar valor da solu� MEF esquerda: saida livre
400  solr[0] = soll[0];
401  break;
402  default:
403  PZError << "TPZConsLawTest::ContributeInterface Boundary Condition Type Not Exists\n";
404  }
405 }
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
int fTest
Integer for integration degree of the initial solution.
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
clarg::argBool bc("-bc", "binary checkpoints", false)
int Type()
Definition: pzbndcond.h:249
Definition: test.py:1
int Dimension() const override
Returns the dimension of the problem.
Definition: pzconslaw.h:276
Defines PZError.
void ComputeSolRight(TPZVec< STATE > &solr, TPZVec< STATE > &soll, TPZVec< REAL > &normal, TPZBndCond *bcright)
void degree(int root, int adj_num, int adj_row[], int adj[], int mask[], int deg[], int *iccsze, int ls[], int node_num)
Definition: rcm.cpp:875
virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Contributes to the residual vector and tangent matrix the face-based quantities.
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Contributes to the residual vector and tangent matrix the volume-based quantities.
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)
virtual int NStateVariables() const override
Number of state variables according to the dimension.
virtual int VariableIndex(const std::string &name) override
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
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...
virtual void Print(std::ostream &out) override
Prints the state of internal variables.
REAL TimeStep()
Returns the value of the time step.
Definition: pzconslaw.h:307
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
void ComputeSolLeft(TPZVec< STATE > &solr, TPZVec< STATE > &soll, TPZVec< REAL > &normal, TPZBndCond *bcleft)
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
expr_ dx(i) *cos(expr_.val())
virtual void Flux(TPZVec< REAL > &x, TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, TPZVec< STATE > &flux) override
Compute the value of the flux function to be used by ZZ error estimator.
STATE B(int i, TPZVec< REAL > &x)
Contains TPZMatrixclass which implements full matrix (using column major representation).
#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
TPZConsLawTest(int nummat, TPZVec< STATE > B, int artdiff, STATE delta_t, int dim, STATE delta, int test=0)
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
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
string res
Definition: test.py:151
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
static REAL gBigNumber
Big number to penalization method, used for Dirichlet conditions.
Definition: TPZMaterial.h:83
Contains TPZMatrix<TVar>class, root matrix class.
virtual ~TPZConsLawTest()
STATE T(int jn, TPZVec< REAL > &x)
Implements the interface for conservation laws, keeping track of the timestep as well.
Definition: pzconslaw.h:71
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...
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
REAL CFL()
Returns the CFL number.
Definition: pzconslaw.h:281
Contains the TPZConsLawTest class for test. Material as conservation law.
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZVec< STATE > fB
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
def values
Definition: rdt.py:119
TPZSolVec sol
vector of the solutions at the integration point
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
virtual std::string Name() override
Returns the material name.
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Contributes to the residual vector the boundary conditions.
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15