currentFE.hpp 15.5 KB
Newer Older
1
/* -*- mode: c++ -*-
2
 This file is part of the LifeV library
passerini's avatar
passerini committed
3
 Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politecnico di Milano
4
5
6
7
8
9
10
11
12
13
14
15
16
17

 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
prudhomm's avatar
prudhomm committed
18
*/
19
20
21
#ifndef _CURRENTFE_H
#define _CURRENTFE_H

prudhomm's avatar
prudhomm committed
22
23
24
#include <life/lifecore/life.hpp>
#include <life/lifefem/geoMap.hpp>
#include <life/lifefem/refFE.hpp>
simone's avatar
simone committed
25
//#include <life/lifefem/geoMap.hpp>
26
27
28
29
30
/*!
  \file currentFE.h
  \brief Structure for the current finite element
*/

prudhomm's avatar
prudhomm committed
31
32
namespace LifeV
{
33
34
/*!
  \class CurrentFE
prudhomm's avatar
prudhomm committed
35
  \brief The class for a finite element
36
  \author J.-F. Gerbeau
prudhomm's avatar
prudhomm committed
37
  \date 04/2002
prudhomm's avatar
prudhomm committed
38

39
40
  modified: the update methods to have
            them dependent on the dimension (lifev becomes multiscale...)
41
     (nbCoor = Nb Dimension = 1, 2 or 3).
prudhomm's avatar
prudhomm committed
42

43
44
45
     I removed the Macro "#if defined(THREEDIM)" to
     use a normal switch.
     Question: IS IT TOO SLOW???
prudhomm's avatar
prudhomm committed
46

47
48
49
     I also factorized some code, and
     postponed the implementation of template methods
     after the class declaration.
50
51
  \author Vincent Martin
  \date 07/2004
52
53
*/

54
55
class CurrentFE
{
56
private:
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
    //! compute only the jacobian matrix
    void _comp_jacobian();
    //! call _comp_jacobian() and compute its determinant
    void _comp_jacobian_and_det();
    //! call _comp_jacobian() and compute its inverse and its determinant
    void _comp_inv_jacobian_and_det();
    void _comp_quad_point_coor();

    //! update the definition of the geo points  (VM 07/2004)
    //! contains a switch over nbCoor (= dimension of pb): slow???
    template <class GEOELE>
    void _update_point( const GEOELE& geoele );

    //! compute phiDer
    void _comp_phiDer();
    //! compute the second derivative phiDer2
    void _comp_phiDer2();
    //! compute phiDer and phiDer2
    void _comp_phiDerDer2();


    UInt _currentId;
simone's avatar
simone committed
79
80
81
82
    UInt _currentLocalId;



83
#ifdef TEST_PRE
84
85
86
87
88

    bool _hasJac;
    bool _hasFirstDeriv;
    bool _hasSecondDeriv;
    bool _hasQuadPtCoor;
89
90
#endif
public:
91
    CurrentFE( const RefFE& _refFE, const GeoMap& _geoMap, const QuadRule& _qr );
92
93
94
95
private:
    CurrentFE( );
    CurrentFE( const CurrentFE& );
public:
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
    const int nbGeoNode;
    const int nbNode;
    const int nbCoor;
    const int nbQuadPt;
    const int nbDiag;
    const int nbUpper;
    const int nbPattern;
    KNM<Real> point;
    const RefFE& refFE;
    const GeoMap& geoMap;
    const QuadRule& qr;
    KNM<Real> phi;
    KNMK<Real> dPhiRef;
    KNMKL<Real> dPhiRef2;
    //
    KNMK<Real> phiDer;
    KNMK<Real> jacobian;
    KNMK<Real> tInvJac;
    KNM<Real> phiGeo;
    KNMK<Real> dPhiGeo;
    KN<Real> weightDet;
    KN<Real> detJac;
    KNM<Real> quadPt;  //!< Coordinates of the quadrature points on the current element
    //second derivatives not yet done (four dimensions KNMKL ?)
    KNMKL<Real> phiDer2;
    //--------------------------------------------------------------------------
    /*!
      return the diameter of the element in the 1-norm
    */
    Real diameter() const;
grandper's avatar
grandper committed
126
127
128
129
    /*!
      return the diameter of the element in the 2-norm
    */
    Real diameter2() const;
130
131
132
133
134
135
136
    /*!
      return the id of the current element (updated with the update* functions)
     */
    inline UInt currentId() const
    {
        return _currentId;
    }
simone's avatar
simone committed
137
138
139
140
141

    inline UInt currentLocalId() const
    {
        return _currentLocalId;
    }
142
#ifdef TEST_PRE
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
    /*!
      return true if the determinant has been updated
      (can ONLY be used if TEST_PRE is defined at compilation time)
     */
    inline bool hasJac() const
    {
        return _hasJac;
    }
    /*!
      return true if the first derivatives have been updated
      (can ONLY be used if TEST_PRE is defined at compilation time)
     */
    inline bool hasFirstDeriv() const
    {
        return _hasFirstDeriv;
    }
    /*!
      return true if the second derivatives have been updated
      (can ONLY be used if TEST_PRE is defined at compilation time)
     */
    inline bool hasSecondDeriv() const
    {
        return _hasSecondDeriv;
    }
    /*!
      return true if the coordinate of the quadrature points have been updated
      (can ONLY be used if TEST_PRE is defined at compilation time)
     */
    inline bool hasQuadPtCoor() const
    {
        return _hasQuadPtCoor;
    }
175
#endif
176
177
178
179
180
181
182
183
184
    /*!
      compute the coordinate (x,y,z)= F(xi,eta,zeta), (F: geo mappping)
      where (xi,eta,zeta) are the coor in the ref element
      and   (x,y,z) the coor in the current element
      (if the code is compiled in 2D mode then z=0 and zeta is disgarded)
    */
    void coorMap( Real& x, Real& y, Real& z,
                  const Real & xi, const Real & eta, const Real &
                  zeta ) const;
quinodoz's avatar
   
quinodoz committed
185
186
187
188
    /*!
      compute the coordinate (xi,eta,zeta)=inv(F)(x,y,z)
    */
    void coorBackMap(const Real& x, const Real& y, const Real& z,
grandper's avatar
grandper committed
189
190
		     Real & xi, Real & eta, Real& zeta) const;

quinodoz's avatar
   
quinodoz committed
191
192
193
194
195
    /*!
      compute the jacobian at a given point : d x_compx / d zeta_compzeta
    */
    Real pointJacobian(const Real& hat_x, const Real& hat_y, const Real& hat_z,
		  int compx, int compzeta) const;
grandper's avatar
grandper committed
196

quinodoz's avatar
   
quinodoz committed
197
198
199
200
201
202
203
    /*!
      compute the inverse jacobian
     */

    Real pointInverseJacobian(const Real& hat_x, const Real& hat_y, const Real& hat_z,
		  int compx, int compzeta) const;

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
    /*!
      return the barycenter of the element
     */
    void barycenter( Real& x, Real& y, Real& z );
    /*!  return (x,y,z) = the global coordinates of the quadrature point ig
      in the current element. \warning this function is almost obsolete since if
      you call the function updateFirstDerivQuadPt rather than updateFirstDeriv
      (for example), the coordinates of the quadrature points have already
      been computed and may be obtained via quadPt(ig,icoor). This is usually
      much less expensive since it avoids many calls to coorQuadPt
    */
    inline void coorQuadPt( Real& x, Real& y, Real& z, const int ig ) const
    {
        coorMap( x, y, z, qr.quadPointCoor( ig, 0 ), qr.quadPointCoor( ig, 1 ),
                 qr.quadPointCoor( ig, 2 ) );
    }
    //!  patternFirst(i): row index in the element matrix of the i-th term of the pattern
    inline int patternFirst( int i ) const
    {
        return refFE.patternFirst( i );
    }
    //! patternSecond(i): column index in the element matrix of the i-th term of the pattern
    inline int patternSecond( int i ) const
    {
        return refFE.patternSecond( i );
    }
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
    /*!
      Return the measure of the current element
    */
    Real measure() const;

    //---------------------------------------
    //! DECLARATION of the update methods
    //---------------------------------------
    /*!
      minimal update: we just identify the id of the current element
    */
    template <class GEOELE>
    void update( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian on
      the current element
    */
    template <class GEOELE>
    void updateJac( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian and quadPt
      on the current element
    */
    template <class GEOELE>
    void updateJacQuadPt( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian,
      tInvJac, phiDer on the current element
    */
    template <class GEOELE>
    void updateFirstDeriv( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian,
      tInvJac, phiDer and quadPt on the current element
    */
    template <class GEOELE>
    void updateFirstDerivQuadPt( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian,
      tInvJac, phiDer2 on the current element
    */
    template <class GEOELE>
    void updateSecondDeriv( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian,
      tInvJac, phiDer2 on the current element
    */
    template <class GEOELE>
    void updateSecondDerivQuadPt( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian,
      tInvJac, phiDer, phiDer2 on the current element
    */
    template <class GEOELE>
    void updateFirstSecondDeriv( const GEOELE& geoele );
    /*!
      compute the arrays detJac, weightDet, jacobian,
      tInvJac, phiDer, phiDer2 on the current element
    */
    template <class GEOELE>
    void updateFirstSecondDerivQuadPt( const GEOELE& geoele );
292
293
294
295
296

};


//----------------------------------------------------------------------
prudhomm's avatar
prudhomm committed
297
//! IMPLEMENTATION
298
299
//----------------------------------------------------------------------
//! update the definition of the geo points (VM 07/04)
300
301
template <class GEOELE>
void CurrentFE::_update_point( const GEOELE& geoele )
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
    //! old way: it was practical, but is it still working???
    // for(int icoor=0;icoor<nbCoor;icoor++)
    // point(i,icoor) =  geoele.coor(i+1,icoor+1);

    switch ( nbCoor )
    {
    case 1:    //! 1D
        for ( int i = 0;i < nbGeoNode;i++ )
        {
            point( i, 0 ) = geoele.point( i + 1 ).x();
        }
        break;
    case 2:    //! 2D
        for ( int i = 0;i < nbGeoNode;i++ )
        {
            point( i, 0 ) = geoele.point( i + 1 ).x();
            point( i, 1 ) = geoele.point( i + 1 ).y();
        }
        break;
    case 3:    //! 3D
        for ( int i = 0;i < nbGeoNode;i++ )
        {
            point( i, 0 ) = geoele.point( i + 1 ).x();
            point( i, 1 ) = geoele.point( i + 1 ).y();
            point( i, 2 ) = geoele.point( i + 1 ).z();
        }
        break;
    default:
winkelma's avatar
winkelma committed
331
332
333
        std::ostringstream err_msg;
        err_msg << "Dimension " << nbCoor << " not available.";

334
        ERROR_MSG( err_msg.str().c_str() );
335
336
    }

337
338
339
340
341
342
343
344
}

//---------------------------------------
//! IMPLEMENTATION of the CurrentFE::update methods
//---------------------------------------

/*!
    minimal update: we just identify the id of the current element
345
  */
346
347
template <class GEOELE>
void CurrentFE::update( const GEOELE& geoele )
348
{
349
#ifdef TEST_PRE
350
351
352
353
    _hasJac = false;
    _hasFirstDeriv = false;
    _hasSecondDeriv = false;
    _hasQuadPtCoor = false;
354
#endif
355

simone's avatar
simone committed
356
357
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
358
}
359

360
361
362
363
/*!
    compute the arrays detJac, weightDet, jacobian on
    the current element
*/
364
365
template <class GEOELE>
void CurrentFE::updateJac( const GEOELE& geoele )
366
367
{
#ifdef TEST_PRE
368
369
370
371
    _hasJac = true;
    _hasFirstDeriv = false;
    _hasSecondDeriv = false;
    _hasQuadPtCoor = false;
372
#endif
373

simone's avatar
simone committed
374
375
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
376
377
378
379
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the jacobian and its determinant...
    _comp_jacobian_and_det();
380
}
381

382
383
384
385
/*!
    compute the arrays detJac, weightDet, jacobian and quadPt
    on the current element
*/
386
387
template <class GEOELE>
void CurrentFE::updateJacQuadPt( const GEOELE& geoele )
388
389
{
#ifdef TEST_PRE
390
391
392
393
    _hasJac = true;
    _hasFirstDeriv = false;
    _hasSecondDeriv = false;
    _hasQuadPtCoor = true;
394
#endif
395

simone's avatar
simone committed
396
397
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
398
399
400
401
402
403
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the jacobian and its determinant...
    _comp_jacobian_and_det();
    //! and the coordinates of the quadrature points
    _comp_quad_point_coor();
404
405
406
407
408
409
}

/*!
    compute the arrays detJac, weightDet, jacobian,
    tInvJac, phiDer on the current element
*/
410
411
template <class GEOELE>
void CurrentFE::updateFirstDeriv( const GEOELE& geoele )
412
413
{
#ifdef TEST_PRE
414
415
416
417
    _hasJac = true;
    _hasFirstDeriv = true;
    _hasSecondDeriv = false;
    _hasQuadPtCoor = false;
418
#endif
419

simone's avatar
simone committed
420
421
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
422
423
424
425
426
427
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the inverse jacobian...
    _comp_inv_jacobian_and_det();
    //! product InvJac by dPhiRef to compute phiDer
    _comp_phiDer();
428
429
430
431
432
433
}

/*!
    compute the arrays detJac, weightDet, jacobian,
    tInvJac, phiDer and quadPt on the current element
*/
434
435
template <class GEOELE>
void CurrentFE::updateFirstDerivQuadPt( const GEOELE& geoele )
436
437
{
#ifdef TEST_PRE
438
439
440
441
    _hasJac = true;
    _hasFirstDeriv = true;
    _hasSecondDeriv = false;
    _hasQuadPtCoor = true;
442
#endif
443

simone's avatar
simone committed
444
445
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
446
447
448
449
450
451
452
453
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the inverse jacobian...
    _comp_inv_jacobian_and_det();
    //! product InvJac by dPhiRef to compute phiDer
    _comp_phiDer();
    //! and the coordinates of the quadrature points
    _comp_quad_point_coor();
454
455
456
457
458
459
460
}

// A. Veneziani, October 30, 2002
/*!
  compute the arrays detJac, weightDet, jacobian,
  tInvJac, phiDer2 on the current element
*/
461
462
template <class GEOELE>
void CurrentFE::updateSecondDeriv( const GEOELE& geoele )
463
464
{
#ifdef TEST_PRE
465
466
467
468
    _hasJac = true;
    _hasFirstDeriv = false;
    _hasSecondDeriv = true;
    _hasQuadPtCoor = false;
469
#endif
470

simone's avatar
simone committed
471
472
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
473
474
475
476
477
478
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the inverse jacobian...
    _comp_inv_jacobian_and_det();
    //! compute the second derivative
    _comp_phiDer2();
479
480
481
}

/*!
482
483
484
    compute the arrays detJac, weightDet, jacobian,
    tInvJac, phiDer2 on the current element
  */
485
486
template <class GEOELE>
void CurrentFE::updateSecondDerivQuadPt( const GEOELE& geoele )
487
{
488
#ifdef TEST_PRE
489
490
491
492
    _hasJac = true;
    _hasFirstDeriv = false;
    _hasSecondDeriv = true;
    _hasQuadPtCoor = true;
493
#endif
494

simone's avatar
simone committed
495
496
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
497
498
499
500
501
502
503
504
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the inverse jacobian...
    _comp_inv_jacobian_and_det();
    //! compute the second derivative
    _comp_phiDer2();
    //! and the coordinates of the quadrature points
    _comp_quad_point_coor();
505
506
}
/*!
507
508
509
    compute the arrays detJac, weightDet, jacobian,
    tInvJac, phiDer, phiDer2 on the current element
  */
510
511
template <class GEOELE>
void CurrentFE::updateFirstSecondDeriv( const GEOELE& geoele )
512
{
513
#ifdef TEST_PRE
514
515
516
517
    _hasJac = true;
    _hasFirstDeriv = true;
    _hasSecondDeriv = true;
    _hasQuadPtCoor = false;
518
#endif
519

simone's avatar
simone committed
520
521
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
522
523
524
525
526
527
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the inverse jacobian...
    _comp_inv_jacobian_and_det();
    //! compute phiDer and phiDer2
    _comp_phiDerDer2();
528
529
}
/*!
530
531
532
    compute the arrays detJac, weightDet, jacobian,
    tInvJac, phiDer, phiDer2 on the current element
  */
533
534
template <class GEOELE>
void CurrentFE::updateFirstSecondDerivQuadPt( const GEOELE& geoele )
535
{
536
#ifdef TEST_PRE
537
538
539
540
    _hasJac = true;
    _hasFirstDeriv = true;
    _hasSecondDeriv = true;
    _hasQuadPtCoor = true;
541
#endif
542

simone's avatar
simone committed
543
544
    _currentId      = geoele.id();
    _currentLocalId = geoele.localId();
545
546
547
548
549
550
551
552
    //! update the definition of the geo points
    _update_point( geoele );
    //! compute the inverse jacobian...
    _comp_inv_jacobian_and_det();
    //! compute phiDer and phiDer2
    _comp_phiDerDer2();
    //! and the coordinates of the quadrature points
    _comp_quad_point_coor();
553
}
prudhomm's avatar
prudhomm committed
554
}
555
556

#endif