NeoPZ
TPZEuler.cpp
Go to the documentation of this file.
1 
6 #include "TPZEuler.h"
7 #include "pzerror.h"
8 #include "pzbndcond.h"
9 
10 using namespace std;
11 TEulerDiffusivity TPZEuler::gEul;
12 
13 void TPZEuler::SetData(istream &data) {
14  TPZMaterial::SetData(data);
15  data >> fDeltaT;
16 }
17 TPZMaterial * TPZEuler::NewMaterial() {
18  TPZEuler *result = new TPZEuler(*this);
19 
20  return result;
21 }
22 void TPZEuler::Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,
23  TPZVec<STATE> &Solout){
24  if(var == 1) {
25  Solout.Resize(1);
26  Solout[0] = gEul.Pressure(Sol);
27  } else if(var ==2) {
28  Solout.Resize(1);
29  Solout[0] = Sol[0];
30  } else if(var == 3) {
31  Solout.Resize(2);
32  Solout[0] = Sol[1];
33  Solout[1] = Sol[2];
34  } else TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
35 }
36 int TPZEuler::NSolutionVariables(int index) {
37  if(index == 1) return 1;
38  if(index == 2) return 1;
39  if(index == 3) return 2;
40  return TPZMaterial::NSolutionVariables(index);
41 }
42 int TPZEuler::VariableIndex(const std::string &name) {
43  if(!strcmp(name.c_str(),"pressure")) return 1;
44  if(!strcmp(name.c_str(),"density")) return 2;
45  if(!strcmp(name.c_str(),"velocity")) return 3;
46  return TPZMaterial::VariableIndex(name);
47 }
48 void TPZEuler::Print(std::ostream & out) {
49  TPZMaterial::Print(out);
50 }
51 void TPZEuler::ContributeBC(TPZMaterialData &data,REAL weight,
52  TPZFMatrix<STATE> &ek,
53  TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
54  // TPZFMatrix<REAL> &dphi = data.dphix;
55  // TPZFMatrix<REAL> &dphiL = data.dphixl;
56  // TPZFMatrix<REAL> &dphiR = data.dphixr;
57  TPZFMatrix<REAL> &phi = data.phi;
58  // TPZFMatrix<REAL> &phiL = data.phil;
59  // TPZFMatrix<REAL> &phiR = data.phir;
60  // TPZManVector<REAL,3> &normal = data.normal;
61  // TPZManVector<REAL,3> &x = data.x;
62  // int &POrder=data.p;
63  // int &LeftPOrder=data.leftp;
64  // int &RightPOrder=data.rightp;
65  int numbersol = data.sol.size();
66  if (numbersol != 1) {
67  DebugStop();
68  }
69 
70  TPZVec<STATE> &sol=data.sol[0];
71  // TPZVec<REAL> &solL=data.soll;
72  // TPZVec<REAL> &solR=data.solr;
73  // TPZFMatrix<REAL> &dsol=data.dsol;
74  // TPZFMatrix<REAL> &dsolL=data.dsoll;
75  // TPZFMatrix<REAL> &dsolR=data.dsolr;
76  // REAL &faceSize=data.HSize;
77  // TPZFMatrix<REAL> &daxesdksi=data.daxesdksi;
78  TPZFMatrix<REAL> &axes=data.axes;
79 
80  if(fState == 0) return;
81  if(bc.Material() != this){
82  PZError << "TPZMat1dLin.apply_bc warning : this material didn't create the boundary condition!\n";
83  }
84 
85  if(bc.Type() < 0 && bc.Type() > 3){
86  PZError << "TPZEuler.aplybc, unknown boundary condition type :" <<
87  bc.Type() << " boundary condition ignored\n";
88  return;
89  }
90 
91 
92  int numdof = NStateVariables();
93  int numnod = ek.Rows()/numdof;
94  int r = numdof;
95 
96  TPZVec<STATE> flux(8);
97  gEul.Flux(sol,flux);
98  REAL normal[2] = {axes(0,1),-axes(0,0)};
99  /*
100  int i;
101  cout << " flux = " << x[0] << ' ' << x[1] << ' ' << "normal" << normal[0] << ' ' << normal[1] << ' ';
102  for(i=0; i<4; i++) cout << flux[i+4] << ' ';
103  cout << endl;
104  */
105  int idf,jdf,in,jn;
106  switch(bc.Type()){
107 
108  case 0:
109  for(in=0 ; in<numnod ; ++in){
110  for(idf = 0;idf<r;idf++) {
111  (ef)(in*r+idf,0) += gBigNumber*phi(in,0)*bc.Val2()(idf,0)*weight;
112  }
113  for(jn=0 ; jn<numnod ; ++jn) {
114  for(idf = 0;idf<r;idf++) {
115  ek(in*r+idf,jn*r+idf) += gBigNumber*phi(in,0)*phi(jn,0)*weight;
116  }
117  }
118  }
119  break;
120 
121  case 1:
122  for(in=0 ; in<numnod ; ++in){
123  for(idf = 0;idf<r;idf++) {
124  (ef)(in*r+idf,0) += phi(in,0)*bc.Val2()(idf,0)*weight;
125  }
126  }
127  break;
128 
129  case 2:
130  for(in=0 ; in<numnod ; ++in){
131  for(idf = 0;idf<r;idf++) {
132  (ef)(in*r+idf,0) += phi(in,0)*bc.Val2()(idf,0)*weight;
133  }
134  for(jn=0 ; jn<numnod ; ++jn) {
135  for(idf = 0;idf<r;idf++) {
136  for(jdf = 0;jdf<r;jdf++) {
137  ek(in*r+idf,jn*r+jdf) += bc.Val1()(idf,jdf)*phi(in,0)*phi(jn,0)*weight;
138  }
139  }
140  }
141  case 3: {
142  TPZFMatrix<STATE> A(4,4),B(4,4);
143  gEul.JacobFlux(sol,A,B);
144  for(in=0; in<numnod; in++) {
145  for(idf=0; idf<4; idf++) {
146  /*
147  ef(4*in+idf) += weight*fDeltaT*(
148  -phi(in,0)*flux[idf]*normal[0]
149  -phi(in,0)*flux[idf+4]*normal[1]
150  );
151  */
152  for(jn=0; jn<numnod; jn++) {
153  for(jdf=0; jdf<4; jdf++) {
154  ek(4*in+idf,4*jn+jdf) += weight*fDeltaT*
155  (phi(in,0)*A(idf,jdf)*phi(jn,0)*normal[0]
156  + phi(in,0)*B(idf,jdf)*phi(jn,0)*normal[1]);
157 
158  }
159  }
160  }
161  }
162  }
163  }//fim switch
164  }
165 }
166 void TPZEuler::Contribute(TPZMaterialData &data, REAL weight,TPZFMatrix<STATE> &ek,
167  TPZFMatrix<STATE> &ef) {
168  TPZFMatrix<REAL> &dphi = data.dphix;
169  // TPZFMatrix<REAL> &dphiL = data.dphixl;
170  // TPZFMatrix<REAL> &dphiR = data.dphixr;
171  TPZFMatrix<REAL> &phi = data.phi;
172  // TPZFMatrix<REAL> &phiL = data.phil;
173  // TPZFMatrix<REAL> &phiR = data.phir;
174  // TPZManVector<REAL,3> &normal = data.normal;
175  TPZManVector<REAL,3> &x = data.x;
176  // int &POrder=data.p;
177  // int &LeftPOrder=data.leftp;
178  // int &RightPOrder=data.rightp;
179  int numbersol = data.sol.size();
180  if (numbersol != 1) {
181  DebugStop();
182  }
183 
184  TPZVec<STATE> &sol=data.sol[0];
185  // TPZVec<REAL> &solL=data.soll;
186  // TPZVec<REAL> &solR=data.solr;
187  // TPZFMatrix<REAL> &dsol=data.dsol;
188  // TPZFMatrix<REAL> &dsolL=data.dsoll;
189  // TPZFMatrix<REAL> &dsolR=data.dsolr;
190  // REAL &faceSize=data.HSize;
191  TPZFMatrix<REAL> &daxesdksi=data.jacinv;
192  TPZFMatrix<REAL> &axes=data.axes;
193 
194  int nshape = phi.Rows();
195  REAL dphix[2];
196  int in,jn,idf,jdf;
197  for(in=0; in<nshape; in++) {
198  dphix[0] = dphi(0,in)*axes(0,0)+dphi(1,in)*axes(1,0);
199  dphix[1] = dphi(1,in)*axes(0,1)+dphi(1,in)*axes(1,1);
200  dphi(0,in) = dphix[0];
201  dphi(1,in) = dphix[1];
202  }
203 
204  if(fState == 0 && ! fForcingFunction) {
205  cout << "TPZEuler failed to suply initial conditions\n";
206  return;
207  }
208  if(fState == 0) {
209  TPZVec<STATE> force(4);
210  fForcingFunction->Execute(x,force);
211  for(in=0; in<nshape; in++) {
212  for(idf=0; idf<4; idf++) {
213  ef(4*in+idf) += weight*phi(in,0)*gBigNumber*force[idf];
214  for(jn=0; jn<nshape; jn++) {
215  ek(4*in+idf,4*jn+idf) += weight*phi(in,0)*phi(jn,0)*gBigNumber;
216  }
217  }
218  }
219  return;
220  }
221  TPZVec<STATE> flux(8);
222  gEul.Flux(sol,flux);
223  /*
224  int i;
225  cout << " flux = " << x[0] << ' ' << x[1] << ' ';
226  for(i=0; i<4; i++) cout << flux[i+4] << ' ';
227  cout << endl << "dsol";
228  for(i=0; i<4; i++) cout << dsol(1,i) << ' ';
229  cout << endl;
230  */
231  TPZFMatrix<STATE> KXX(4,4),KXY(4,4),KYX(4,4),KYY(4,4),A(4,4),B(4,4);
232  TPZFMatrix<REAL> jacinv(2,2);
233  REAL jacdet = daxesdksi(0,0)*daxesdksi(1,1)-daxesdksi(0,1)*daxesdksi(1,0);
234  jacinv(0,0) = daxesdksi(1,1)/jacdet;
235  jacinv(1,1) = daxesdksi(0,0)/jacdet;
236  jacinv(0,1) = -daxesdksi(0,1)/jacdet;
237  jacinv(1,0) = -daxesdksi(1,0)/jacdet;
238  gEul.MatrixDiff(sol,axes,jacinv,KXX,KXY,KYX,KYY);
239  gEul.JacobFlux(sol,A,B);
240  for(in=0; in<nshape; in++) {
241  for(idf=0; idf<4; idf++) {
242  ef(4*in+idf) += weight*phi(in,0)*sol[idf];
243 
244  //+fDeltaT*dphi(0,in)*flux[idf]
245  //+fDeltaT*dphi(1,in)*flux[idf+4]);
246  /*
247  for(jdf=0; jdf<4; jdf++) {
248  ef(in*4+idf,0) += weight*fDeltaT*(
249  -dphi(0,in)*KXX(idf,jdf)*dsol(0,jdf)
250  -dphi(0,in)*KXY(idf,jdf)*dsol(1,jdf)
251  -dphi(1,in)*KYX(idf,jdf)*dsol(0,jdf)
252  -dphi(1,in)*KYY(idf,jdf)*dsol(1,jdf)
253  );
254  }
255  */
256  for(jn=0; jn<nshape; jn++) {
257  for(jdf=0; jdf<4;jdf++) {
258 
259  ek(4*in+idf,4*jn+jdf) += fDeltaT*weight*(
260  dphi(0,in)*KXX(idf,jdf)*dphi(0,jn)
261  +dphi(0,in)*KXY(idf,jdf)*dphi(1,jn)
262  +dphi(1,in)*KYX(idf,jdf)*dphi(0,jn)
263  +dphi(1,in)*KYY(idf,jdf)*dphi(1,jn)
264 
265  -dphi(0,in)*A(idf,jdf)*phi(jn)
266  -dphi(1,in)*B(idf,jdf)*phi(jn)
267  );
268  }
269  ek(4*in+idf,4*jn+idf) += weight*phi(in,0)*phi(jn,0);
270  }
271  }
272  }
273 }
274 int TPZEuler::NStateVariables() const {
275  return 4;
276 }
277 int TPZEuler::Dimension() const {
278  return 2;
279 }
280 TPZEuler::TPZEuler(TPZEuler & copy)
281 : TPZRegisterClassId(&TPZEuler::ClassId), TPZMaterial(copy){
282  fDeltaT = copy.fDeltaT;
283  fState = copy.fState;
284 }
285 TPZEuler::TPZEuler(int id, STATE deltat)
286 : TPZRegisterClassId(&TPZEuler::ClassId), TPZMaterial(id){
287  fDeltaT = deltat;
288  fState = 0;
289 }
290 
291 int TPZEuler::ClassId() const{
292  return Hash("TPZEuler") ^ TPZMaterial::ClassId() << 1;
293 }
TPZMaterial::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout)
Returns the solution associated with the var index based on the finite element approximation.
Definition: TPZMaterial.cpp:200
TPZEuler::gEul
static TEulerDiffusivity gEul
Definition: TPZEuler.h:106
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZEuler::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: TPZEuler.cpp:36
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZMaterial::VariableIndex
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Definition: TPZMaterial.cpp:141
pzerror.h
Defines PZError.
TPZEuler::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Computes contribution to the stiffness matrix and right hand side at an integration point...
Definition: TPZEuler.cpp:166
TPZEuler.h
Contains the TPZEuler class which implements a a linear scalar convection equation.
TPZEuler::TPZEuler
TPZEuler(TPZEuler &copy)
Copy constructor.
Definition: TPZEuler.cpp:280
TPZMaterialData::jacinv
TPZFNMatrix< 9, REAL > jacinv
value of the inverse of the jacobian at the integration point
Definition: pzmaterialdata.h:67
TPZEuler::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
Definition: TPZEuler.cpp:22
TPZEuler::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: TPZEuler.cpp:291
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::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
TPZMaterial
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
TPZEuler::VariableIndex
virtual int VariableIndex(const std::string &name) override
Returns the variable index associated with the name.
Definition: TPZEuler.cpp:42
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
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
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
TPZEuler
This class implements a linear scalar convection equation with modified SUPG difusion.
Definition: TPZEuler.h:22
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TEulerDiffusivity
Implements a numerical diffusivity coeficient for the SUPG method. Analysis: Solving process Analysis...
Definition: eulerdif.h:40
TPZMaterialData::axes
TPZFNMatrix< 9, REAL > axes
axes indicating the directions of the derivatives of the shapefunctions
Definition: pzmaterialdata.h:63
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZEuler::SetData
virtual void SetData(std::istream &data) override
Reads data of the material from a istream (file data)
Definition: TPZEuler.cpp:13
TPZEuler::Print
virtual void Print(std::ostream &out=std::cout) override
Print out the data associated with the material.
Definition: TPZEuler.cpp:48
TPZMaterial::ClassId
int ClassId() const override
Unique identifier for serialization purposes.
Definition: TPZMaterial.cpp:382
idf
static int idf[6]
Definition: TPZCompMeshTools.cpp:135
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZEuler::NewMaterial
virtual TPZMaterial * NewMaterial() override
To create another material of the same type.
Definition: TPZEuler.cpp:17
TPZEuler::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to the stiffness matrix and load vector at one BC integration point...
Definition: TPZEuler.cpp:51
TPZEuler::fState
int fState
Definition: TPZEuler.h:108
TPZEuler::Dimension
virtual int Dimension() const override
Returns the integrable dimension of the material.
Definition: TPZEuler.cpp:277
TPZEuler::fDeltaT
STATE fDeltaT
Definition: TPZEuler.h:107
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
TPZBndCond::Material
TPZMaterial * Material() const
Definition: pzbndcond.h:263
TPZMaterial::SetData
virtual void SetData(std::istream &data)
Reads data of the material from a istream (file data)
Definition: TPZMaterial.cpp:286
TPZEuler::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: TPZEuler.cpp:274
Generated on Sun Aug 16 2020 11:17:04 for NeoPZ by   doxygen 1.8.13