NavierStokesHandler_miguel.hpp 10.2 KB
Newer Older
prudhomm's avatar
prudhomm committed
1
2
/*
  This file is part of the LifeV library
3
  Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politechnico di Milano
prudhomm's avatar
prudhomm committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
  
  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
/*!
  \file NavierStokesHandler.h
  \author M.A. Fernandez
  \date 06/2003 
  \version 1.0

  \brief This file contains an abstract class for NavierStokes solvers.  

*/

#ifndef _NAVIERSTOKESHANDLER_H_
#define _NAVIERSTOKESHANDLER_H_

#include <cmath>
#include <sstream> 
#include <cstdlib>


#include "dataNavierStokes.hpp"
#include "dataAztec.hpp"
#include "refFE.hpp"
#include "dof.hpp"
#include "lifeV.hpp"
#include "medit_wrtrs.hpp"
#include "bcCond.hpp"


/*! 
  \class NavierStokesHandler

  Abstract class which defines the general structure of a NavierStokes solver.
  For each new NavierStokes solver  we have to implement the corresponding timeAdvance and an iterate methods 
  
*/

template <typename Mesh>
class NavierStokesHandler:
public DataNavierStokes<Mesh> { 
 
 public:

  typedef Real (*Function)(const Real&, const Real&, const Real&, const Real&, const ID&);

  //! Constructor   
  /*!
    \param data_file GetPot data file
    \param refFE_u reference FE for the velocity
    \param refFE_p reference FE for the pressure
    \param Qr_u volumic quadrature rule for the velocity
    \param bdQr_u surface quadrature rule for the velocity 
    \param Qr_p volumic quadrature rule for the pressure
    \param bdQr_p surface quadrature rule for the pressure
    \param BCh_u boundary conditions for the velocity
  */
  NavierStokesHandler(const GetPot& data_file,  const RefFE& refFE_u, 
		const RefFE& refFE_p, const QuadRule& Qr_u, const QuadRule& bdQr_u, 
		const QuadRule& Qr_p, const QuadRule& bdQr_p, BC_Handler& BCh_u);
    
  //! Sets initial condition for the velocity
  void initialize(const Function& u0); 

  //! Update the right  hand side  for time advancing   
  /*! 
    \param source volumic source  
    \param time present time
  */
  virtual void timeAdvance(const Function source, const Real& time) = 0; 

  //! Update convective term, bc treatment and solve the linearized ns system
  virtual void iterate() = 0;

  //! Returns the velocity vector
  PhysVectUnknown<Vector>& u();
  
  //! Returns the pressure
  ScalUnknown<Vector>& p();

  //! Returns the velocity Dof 
  const Dof& uDof() const;

  //! Returns the pressure Dof 
  const Dof& pDof() const;

  //! Postprocessing
  void postProcess();

  //! Do nothing destructor
  virtual ~NavierStokesHandler() {}
  
 protected:
  
  //! Reference FE for the velocity
  const RefFE& _refFE_u;
  
  //! Reference FE for the pressure
  const RefFE& _refFE_p;
 
  //! The Dof object associated with the velocity
  Dof _dof_u;

  //! The Dof object associated with the pressure
  Dof _dof_p;

  //! The number of total velocity dofs  
  UInt _dim_u;

  //! The number of total pressure dofs  
  UInt _dim_p;

  //! Quadrature rule for velocity volumic elementary computations
  const QuadRule& _Qr_u;
  
  //! Quadrature rule for velocity surface elementary computations
  const QuadRule& _bdQr_u;

  //! Quadrature rule for pressure volumic elementary computations
  const QuadRule& _Qr_p;
  
  //! Quadrature rule for pressure surface elementary computations
  const QuadRule& _bdQr_p;

  //! Current FE for the velocity u
  CurrentFE _fe_u;
  CurrentBdFE _feBd_u;

  //! Current FE for the pressure p
  CurrentFE _fe_p;

  //! The velocity
  PhysVectUnknown<Vector> _u;

  //! The pressure
  ScalUnknown<Vector> _p;
   
  //! The BC handler
  BC_Handler& _BCh_u;

  //! The actual time
  Real _time;

  //! Aux. var. for PostProc
  UInt _count;
};



//
// IMPLEMENTATION
//


// Constructor
template <typename Mesh> 
NavierStokesHandler<Mesh>::
NavierStokesHandler(const GetPot& data_file,  const RefFE& refFE_u, 
		    const RefFE& refFE_p, const QuadRule& Qr_u, const QuadRule& bdQr_u, 
		    const QuadRule& Qr_p, const QuadRule& bdQr_p, BC_Handler& BCh_u):
     DataNavierStokes<Mesh>(data_file),    
     _refFE_u(refFE_u),
     _refFE_p(refFE_p),
     _dof_u(_mesh,_refFE_u),
     _dof_p(_mesh,_refFE_p),
     _dim_u(_dof_u.numTotalDof()),
     _dim_p(_dof_p.numTotalDof()), 
     _Qr_u(Qr_u),
     _bdQr_u(bdQr_u),
     _Qr_p(Qr_p),
     _bdQr_p(bdQr_p),
     _fe_u(_refFE_u,_mesh.getGeoMap(),_Qr_u),
     _feBd_u(_refFE_u.boundaryFE(),_mesh.getGeoMap().boundaryMap(),_bdQr_u),
     _fe_p(_refFE_p,_mesh.getGeoMap(),_Qr_p),
     _u(_dim_u),
     _p(_dim_p),
     _BCh_u(BCh_u),
     _time(0),
     _count(0) {}


// Returns the velocity vector
template<typename Mesh> PhysVectUnknown<Vector>& 
NavierStokesHandler<Mesh>::u() {
  return _u;
}

// Returns the pressure
template<typename Mesh> ScalUnknown<Vector>&
NavierStokesHandler<Mesh>::p() {
  return _p;
}

// Returns the velocity Dof 
template<typename Mesh> const Dof& 
NavierStokesHandler<Mesh>::uDof() const {
  return _dof_u;
}

// Returns the pressure Dof 
template<typename Mesh>  const Dof& 
NavierStokesHandler<Mesh>::pDof() const {
  return _dof_p;
}

// Postprocessing 
template<typename Mesh>  void 
NavierStokesHandler<Mesh>::postProcess() {
  ostringstream index;
  string name;
 
  ++_count;
 
  if (fmod(float(_count),float(_verbose)) == 0.0) {
    cout << "  o-  Post-processing \n";
    index << (_count/_verbose);
 
    switch( index.str().size() ) {
    case 1:
      name = "00"+index.str();
      break;
    case 2:
      name = "0"+index.str();
      break;
    case 3:
      name = index.str();
      break;
    }
 
 
    wr_medit_ascii_scalar("press."+name+".bb",_p.giveVec(),_p.size());
    wr_medit_ascii_scalar("vel_x."+name+".bb",_u.giveVec(),_mesh.numVertices());
    wr_medit_ascii_scalar("vel_y."+name+".bb",_u.giveVec() + _dim_u,_mesh.numVertices());
    wr_medit_ascii_scalar("vel_z."+name+".bb",_u.giveVec() + 2*_dim_u,_mesh.numVertices());
    // wr_medit_ascii_vector("veloc."+name+".bb",_u.giveVec(),_mesh.numVertices(),_dim_u);
    system(("ln -s "+_mesh_file+" press."+name+".mesh").data());
    system(("ln -s "+_mesh_file+" vel_x."+name+".mesh").data());
    system(("ln -s "+_mesh_file+" vel_y."+name+".mesh").data());
    system(("ln -s "+_mesh_file+" vel_z."+name+".mesh").data());
    // system(("ln -s "+_mesh_file+" veloc."+name+".mesh").data());
  }
}


// Sets the initial condition
template<typename Mesh> void 
NavierStokesHandler<Mesh>::initialize(const Function& u0) {
  
  // Initialize pressure
  _p.vec()=0.0;

  // Initialize velocity

  typedef  typename Mesh::VolumeShape GeoShape; // Element shape
   
  UInt nDofpV  = _refFE_u.nbDofPerVertex; // number of Dof per vertex
  UInt nDofpE  = _refFE_u.nbDofPerEdge;   // number of Dof per edge
  UInt nDofpF  = _refFE_u.nbDofPerFace;   // number of Dof per face
  UInt nDofpEl = _refFE_u.nbDofPerVolume; // number of Dof per Volume
  
  UInt nElemV = GeoShape::numVertices; // Number of element's vertices 
  UInt nElemE = GeoShape::numEdges;    // Number of element's edges
  UInt nElemF = GeoShape::numFaces;    // Number of element's faces
  
  UInt nDofElemV = nElemV*nDofpV; // number of vertex's Dof on a Element
  UInt nDofElemE = nElemE*nDofpE; // number of edge's Dof on a Element
  UInt nDofElemF = nElemF*nDofpF; // number of face's Dof on a Element
    
  ID nbComp = _u.nbcomp(); // Number of components of the mesh velocity

  Real x, y, z;
 
  ID lDof;

  // Loop on elements of the mesh
  for (ID iElem=1; iElem <= _mesh.numVolumes(); ++iElem) {
       
    _fe_u.updateJac( _mesh.volume(iElem) );

    // Vertex based Dof 
    if ( nDofpV ) { 
      
      // loop on element vertices 
      for (ID iVe=1; iVe<=nElemV; ++iVe){
	
	// Loop number of Dof per vertex
	for (ID l=1; l<=nDofpV; ++l) {
	  lDof = (iVe-1)*nDofpV + l; // Local dof in this element
	  
	  // Nodal coordinates
	  _fe_u.coorMap(x, y, z, _fe_u.refFE.xi(lDof-1), _fe_u.refFE.eta(lDof-1), _fe_u.refFE.zeta(lDof-1));  
	  
	  // Loop on data vector components
	  for (UInt icmp=0; icmp < nbComp; ++icmp)
       	    _u.vec()( icmp*_dim_u + _dof_u.localToGlobal(iElem,lDof) - 1 ) = u0(0.0,x,y,z,icmp+1);  
	}
      }
    }
    
    // Edge based Dof 
    if (nDofpE) { 
	
      // loop on element edges 
      for (ID iEd=1; iEd <=nElemE; ++iEd) {
	
	// Loop number of Dof per edge
	for (ID l=1; l<=nDofpE; ++l) {
	  lDof = nDofElemV + (iEd-1)*nDofpE + l; // Local dof in the adjacent Element
	
	  // Nodal coordinates
	  _fe_u.coorMap(x, y, z, _fe_u.refFE.xi(lDof-1), _fe_u.refFE.eta(lDof-1), _fe_u.refFE.zeta(lDof-1));
	
	  // Loop on data vector components
	  for (UInt icmp=0; icmp < nbComp; ++icmp) 
	    _u.vec()( icmp*_dim_u + _dof_u.localToGlobal(iElem,lDof) - 1 ) = u0(0.0,x,y,z,icmp+1);
	}
      }
    }  
    
    // Face based Dof 
    if (nDofpF) { 
      
      // loop on element faces
      for (ID iFa=1; iFa <=nElemF; ++iFa) {
	
	// Loop on number of Dof per face
	for (ID l=1; l<=nDofpF; ++l) {
	  
	  lDof = nDofElemE + nDofElemV + (iFa-1)*nDofpF + l; // Local dof in the adjacent Element
     
	  // Nodal coordinates
	  _fe_u.coorMap(x, y, z, _fe_u.refFE.xi(lDof-1), _fe_u.refFE.eta(lDof-1), _fe_u.refFE.zeta(lDof-1));
		  
	  // Loop on data vector components
	  for (UInt icmp=0; icmp < nbComp; ++icmp) 
	    _u.vec()( icmp*_dim_u + _dof_u.localToGlobal(iElem,lDof) - 1) = u0(0.0,x,y,z,icmp+1);   
	}
      }
    }

    // Element based Dof 
    // Loop on number of Dof per Element
    for (ID l=1; l<=nDofpEl; ++l) {
      lDof = nDofElemF + nDofElemE + nDofElemV + l; // Local dof in the Element
        
      // Nodal coordinates
      _fe_u.coorMap(x, y, z, _fe_u.refFE.xi(lDof-1), _fe_u.refFE.eta(lDof-1), _fe_u.refFE.zeta(lDof-1));
	      
      // Loop on data vector components
      for (UInt icmp=0; icmp < nbComp; ++icmp)
	_u.vec()( icmp*_dim_u + _dof_u.localToGlobal(iElem,lDof) - 1) =   u0(0.0,x,y,z,icmp+1);      
    }
  }
}


#endif