This file is indexed.

/usr/include/rheolef/field_nonlinear_expr.h is in librheolef-dev 6.5-1build1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
#ifndef _RHEOLEF_FIELD_NONLINEAR_EXPR_H
#define _RHEOLEF_FIELD_NONLINEAR_EXPR_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef 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 General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
///
/// =========================================================================
//
// field_nonlinear_expr: separate from field_expr
//
// it is a prototype: only compose(f,uh) is still handled
//
// OUTLINE:
// 1) main wrapper class
// 2) unary function call: (f expr)
// 3) binary function call: (f expr1 expr2)
// 4) jump of a field
#include "rheolef/field_nonlinear_expr_terminal.h"
#include "rheolef/field_evaluate.h"
#include "rheolef/tensor4.h"
#include "rheolef/field_functor.h"
#include "rheolef/operators.h"

#include <boost/functional.hpp>

namespace rheolef {

// ---------------------------------------------------------------------------
// 1) main wrapper class
// ---------------------------------------------------------------------------
template<class RawExpr>
class field_nonlinear_expr {
public:
// typedefs:

  typedef geo_element::size_type         size_type;
  typedef typename RawExpr::memory_type  memory_type;
  typedef typename RawExpr::result_type  result_type;
  typedef typename RawExpr::value_type   value_type;
  typedef typename RawExpr::float_type   float_type;
  typedef typename RawExpr::scalar_type  scalar_type;

// alocators:

  field_nonlinear_expr (const RawExpr& raw_expr) 
    : _raw_expr(raw_expr) {}

// accessors:

  static const space_constant::valued_type valued_hint = RawExpr::valued_hint;

  space_constant::valued_type valued_tag() const { return _raw_expr.valued_tag(); }

  void initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x) const {
    _raw_expr.initialize (omega, hat_x);
  }
  bool initialize (const space_basic<float_type,memory_type>& Xh) const {
    return _raw_expr.initialize (Xh);
  }
  template<class Result>
  void evaluate (const geo_element& K, std::vector<Result>& value) const {
    _raw_expr.evaluate (K, value);
  }
  void evaluate_on_side (const geo_element& K, const side_information_type& sid, std::vector<result_type>& value) const {
    _raw_expr.evaluate_on_side (K, sid, value);
  }
  template<class Result>
  void valued_check() const {
    _raw_expr.valued_check<Result>();
  }

protected:
// data:
  RawExpr   _raw_expr;
};
// ---------------------------------------------------------------------------
// 2) unary function call: (f expr)
// ---------------------------------------------------------------------------
template<class UnaryFunction, class Expr>
class field_nonlinear_expr_uf {
public:
// typedefs:

  typedef geo_element::size_type                      size_type;
  typedef typename Expr::memory_type                  memory_type;
  typedef typename details::generic_unary_traits<UnaryFunction>::template result_hint<typename Expr::result_type>::type
						      result_type;
  typedef result_type                                 value_type;
  typedef typename scalar_traits<value_type>::type    scalar_type;
  typedef typename  float_traits<value_type>::type    float_type;
  typedef field_nonlinear_expr_uf<UnaryFunction,Expr> self_type;

// alocators:

  field_nonlinear_expr_uf (const UnaryFunction& f, const field_nonlinear_expr<Expr>& expr) 
    : _f(f), _expr(expr) {}

// accessors:

  static const space_constant::valued_type valued_hint = space_constant::valued_tag_traits<result_type>::value;

  space_constant::valued_type valued_tag() const {
    return details::generic_unary_traits<UnaryFunction>::valued_tag (_expr.valued_tag());
  }

// initializators:

  void initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x) const {
    _expr.initialize (omega, hat_x);
  }
  bool initialize (const space_basic<float_type,memory_type>& Xh) const {
    return _expr.initialize (Xh);
  }
// evaluator:

  template<class Result, class Arg, class Status>
  struct evaluate_call_check {
    void operator() (const self_type& obj, const geo_element& K, std::vector<Result>& value) const {
      fatal_macro ("invalid type resolution: Result="<<typename_macro(Result)
          << ", Arg="<<typename_macro(Arg)
          << ", UnaryFunction="<<typename_macro(UnaryFunction)
      );
    }
  };
  template<class Result, class Arg>
  struct evaluate_call_check<Result,Arg,mpl::true_> {
    void operator() (const self_type& obj, const geo_element& K, std::vector<Result>& value) const {
      std::vector<Arg> tmp_value;
      obj._expr.evaluate (K, tmp_value);
      value.resize(tmp_value.size());
      typename std::vector<Arg>::const_iterator tmp = tmp_value.begin();
      for (typename std::vector<Result>::iterator
	iter = value.begin(),
	last = value.end(); iter != last; ++iter, ++tmp) {
          *iter = obj._f(*tmp);
      }
    }
  };
  template<class Result, class Arg>
  void evaluate_call (const geo_element& K, std::vector<Result>& value) const {
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<Arg,Result>::result_type result_type;
    typedef typename mpl::and_<
		typename details::is_equal<Result,result_type>::type
	       ,typename mpl::not_<typename details::is_error<Arg>::type>::type
	      >::type
      status_t;
    evaluate_call_check<Result,Arg,status_t> eval;
    eval (*this, K, value);
  }
  // when arg is known at run-time:
  template<class This, class Result, class Arg, space_constant::valued_type ArgTag>
  struct evaluate_switch {
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      typedef typename scalar_traits<Arg>::type T;
      space_constant::valued_type arg_valued_tag = obj._expr.valued_tag();
      switch (arg_valued_tag) {
        case space_constant::scalar:
	  obj.template evaluate_call<Result,T> (K, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, point_basic<T> > (K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, tensor_basic<T> > (K, value); break;
        default: { error_macro ("unexpected valued tag="<<arg_valued_tag); }
      }
    }
  };
  // when arg is known at compile-time:
  template<class This, class Result, class Arg>
  struct evaluate_switch <This, Result, Arg, space_constant::scalar> {
    typedef typename scalar_traits<Arg>::type T;
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result,T> (K, value);
    }
  };
  template<class This, class Result, class Arg>
  struct evaluate_switch <This, Result, Arg, space_constant::vector> {
    typedef typename scalar_traits<Arg>::type T;
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result, point_basic<T> > (K, value);
    }
  };
  template<class This, class Result, class Arg>
  struct evaluate_switch <This, Result, Arg, space_constant::tensor> {
    typedef typename scalar_traits<Arg>::type T;
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result, tensor_basic<T> > (K, value);
    }
  };
  template<class Result>
  void evaluate (const geo_element& K, std::vector<Result>& value) const
  {
   typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
                     A1;
    typedef field_nonlinear_expr_uf<UnaryFunction, Expr> This;
    static const space_constant::valued_type argument_tag = space_constant::valued_tag_traits<A1>::value;
    evaluate_switch <This, Result, A1, argument_tag> helper;
    helper.evaluate (*this, K, value);
  }
  template<class Result>
  void evaluate_on_side (const geo_element& K, const side_information_type& sid, std::vector<Result>& value) const { 
    // TODO: group with evaluate
    fatal_macro ("uf::evaluate_on_side: not yet");
  }

  template<class Result>
  void valued_check() const {
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
                     A1;
    if (! is_undeterminated<A1>::value) _expr.valued_check<A1>();
  }

protected:
// data:
  UnaryFunction                 _f;
  field_nonlinear_expr<Expr>    _expr;
};
// ---------------------------------------------------------------------------
// 3) binary function call: (f expr1 expr2)
// ---------------------------------------------------------------------------
template<class BinaryFunction, class Expr1, class Expr2>
class field_nonlinear_expr_bf {
public:
// typedefs:

  typedef geo_element::size_type                   size_type;
  typedef typename details::generic_binary_traits<BinaryFunction>::template result_hint<typename Expr1::result_type,typename Expr2::result_type>::type result_type;
  typedef result_type                              value_type;
  typedef typename scalar_traits<value_type>::type scalar_type;
  typedef typename  float_traits<value_type>::type float_type;
  typedef typename Expr1::memory_type              memory_type;

// alocators:

  field_nonlinear_expr_bf (const BinaryFunction&             f, 
		    const Expr1& expr1,
                    const Expr2& expr2)
    : _f(f), _expr1(expr1), _expr2(expr2)
  {
  }

// accessors:
 
  static const space_constant::valued_type valued_hint = space_constant::valued_tag_traits<result_type>::value;

  space_constant::valued_type valued_tag() const {
    return details::generic_binary_traits<BinaryFunction>::valued_tag(_expr1.valued_tag(), _expr2.valued_tag());
  }

// initializers:

  void initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x) const
  {
    _expr1.initialize (omega, hat_x); 
    _expr2.initialize (omega, hat_x);
  }
  bool initialize (const space_basic<float_type,memory_type>& Xh) const {
    bool is_homogeneous1 = _expr1.initialize (Xh);
    bool is_homogeneous2 = _expr2.initialize (Xh);
    return is_homogeneous1 && is_homogeneous2;
  }

// evaluators:

  template<class Result, class Arg1, class Arg2>
  void evaluate_internal2 (const geo_element& K, std::vector<Result>& value) const {
#ifdef TO_CLEAN
    _check<Result,Arg1,Arg2> ();
#endif // TO_CLEAN
    std::vector<Arg1> tmp1_value;
    std::vector<Arg2> tmp2_value;
    _expr1.evaluate (K, tmp1_value);
    _expr2.evaluate (K, tmp2_value);
    value.resize(tmp1_value.size());
    typename std::vector<Arg1>::const_iterator tmp1 = tmp1_value.begin();
    typename std::vector<Arg2>::const_iterator tmp2 = tmp2_value.begin();
    for (typename std::vector<Result>::iterator
        iter = value.begin(),
	last = value.end(); iter != last; ++iter, ++tmp1, ++tmp2) {
      *iter = _f (*tmp1, *tmp2);
    }
  }
  template<class This, class Result, class ReturnType, class Arg1, class Arg2>
  struct evaluate_internal {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      fatal_macro ("unexpected return type "	
	<< pretty_typename_macro(ReturnType) << ": "
	<< pretty_typename_macro(Result) << " was expected for function "
	<< pretty_typename_macro(BinaryFunction) << "("
	<< pretty_typename_macro(Arg1) << ","
	<< pretty_typename_macro(Arg2) << ")");
    }
  };
  template<class This, class Result, class Arg1, class Arg2>
  struct evaluate_internal<This,Result,Result,Arg1,Arg2> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_internal2 <Result,Arg1,Arg2> (K, value);
    }
  };
  template<class Result, class Arg1, class Arg2>
  void evaluate_call (const geo_element& K, std::vector<Result>& value) const {
    typedef typename details::generic_binary_traits<BinaryFunction>::template result_hint<Arg1,Arg2>::type ReturnType;
    typedef field_nonlinear_expr_bf<BinaryFunction, Expr1, Expr2> This;
    evaluate_internal<This,Result,ReturnType,Arg1,Arg2> eval_int;
    eval_int (*this, K, value);
  }
  // when both args are defined at compile time:
  template<class This, class Result,
	class Arg1,        space_constant::valued_type Arg1Tag,
  	class Arg2,        space_constant::valued_type Arg2Tag>
  struct evaluate_switch {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result, Arg1, Arg2> (K, value);
    }
  };
  // when both args are undefined at compile time:
  template<class This, class Result,
	class Arg1,
  	class Arg2>
  struct evaluate_switch<This, Result,
	Arg1,   space_constant::last_valued,
        Arg2,   space_constant::last_valued> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      typedef typename scalar_traits<Arg1>::type T1;
      typedef typename scalar_traits<Arg2>::type T2;
      space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
      space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
      switch (arg1_valued_tag) {
        case space_constant::scalar: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      obj.template evaluate_call<Result, T1, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, T1, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, T1, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, T1, tensor3_basic<T2> > (K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        case space_constant::vector: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      obj.template evaluate_call<Result, point_basic<T1>, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, point_basic<T1>, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, point_basic<T1>, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, point_basic<T1>, tensor3_basic<T2> > (K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      obj.template evaluate_call<Result, tensor_basic<T1>, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, tensor_basic<T1>, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, tensor_basic<T1>, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, tensor_basic<T1>, tensor3_basic<T2> > (K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        case space_constant::tensor3: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, tensor3_basic<T2> > (K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
      }
    }
  };
  // when only first arg is defined at compile time:
  template<class This, class Result,
	class Arg1,        space_constant::valued_type Arg1Tag,
  	class Arg2>
  struct evaluate_switch<This, Result,
	Arg1,   Arg1Tag,
        Arg2,   space_constant::last_valued> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      typedef typename scalar_traits<Arg2>::type T2;
      space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
      switch (arg2_valued_tag) {
        case space_constant::scalar:
	  obj.template evaluate_call<Result, Arg1, T2>                (K, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, Arg1, point_basic<T2> >  (K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, Arg1, tensor_basic<T2> > (K, value); break;
        case space_constant::tensor3:
	  obj.template evaluate_call<Result, Arg1, tensor3_basic<T2> > (K, value); break;
        default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
      }
    }
  };
  // when only second arg is defined at compile time:
  template<class This, class Result,
	class Arg1,     
  	class Arg2,        space_constant::valued_type Arg2Tag>
  struct evaluate_switch<This, Result,
	Arg1,              space_constant::last_valued,
        Arg2,              Arg2Tag> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      typedef typename scalar_traits<Arg1>::type T1;
      space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
      switch (arg1_valued_tag) {
        case space_constant::scalar:
	  obj.template evaluate_call<Result, T1, Arg2>               (K, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, point_basic<T1>, Arg2>  (K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, tensor_basic<T1>, Arg2> (K, value); break;
        case space_constant::tensor3:
	  obj.template evaluate_call<Result, tensor3_basic<T1>, Arg2> (K, value); break;
        default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
      }
    }
  };
  template<class Result>
  void evaluate (const geo_element& K, std::vector<Result>& value) const
  {
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::first_argument_type   A1;
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::second_argument_type A2;
    static const space_constant::valued_type  first_argument_tag = space_constant::valued_tag_traits<A1>::value;
    static const space_constant::valued_type second_argument_tag = space_constant::valued_tag_traits<A2>::value;
    typedef field_nonlinear_expr_bf<BinaryFunction, Expr1, Expr2> This;
    evaluate_switch <This, Result, A1, first_argument_tag, A2, second_argument_tag>   eval;
    eval (*this, K, value);
  }
  template<class Result>
  void valued_check() const {
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::first_argument_type   A1;
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::second_argument_type A2;
    if (! is_undeterminated<A1>::value) _expr1.valued_check<A1>();
    if (! is_undeterminated<A2>::value) _expr2.valued_check<A2>();
  }

protected:
// data:
  BinaryFunction  _f;
  Expr1           _expr1;
  Expr2           _expr2;
};
// ----------------------------------------------------------------------------
// 4) jump of a field
// ----------------------------------------------------------------------------
template<class Expr>
class field_nonlinear_expr_dg {
public:
// typedefs:

  typedef geo_element::size_type          size_type;
  typedef typename Expr::memory_type      memory_type;
  typedef typename Expr::result_type      result_type;
  typedef typename Expr::value_type       value_type;
  typedef typename Expr::scalar_type      scalar_type;
  typedef typename Expr::float_type       float_type;

// alocators:

  field_nonlinear_expr_dg (const Expr& expr, const float_type& c0, const float_type& c1);

// accessors:

  static const space_constant::valued_type valued_hint = Expr::valued_hint;

  space_constant::valued_type valued_tag() const { return _expr0.valued_tag(); }

  void initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& quad) const;
  bool initialize (const space_basic<float_type,memory_type>& Xh) const;

  template<class ValueType>
  void evaluate (const geo_element& K, std::vector<ValueType>& value) const;

  template<class ValueType>
  void valued_check() const;

// data:
protected:
  Expr                                      _expr0, _expr1;
  float_type                                _c0,    _c1;
  mutable geo_basic<float_type,memory_type> _omega;
};

template<class Expr>
inline
field_nonlinear_expr_dg<Expr>::field_nonlinear_expr_dg (
  const Expr&       expr, 
  const float_type& c0,
  const float_type& c1)
 : _expr0(expr), _expr1(expr),
   _c0(c0),      _c1(c1),
   _omega()
{
}
template<class Expr>
inline
void
field_nonlinear_expr_dg<Expr>::initialize (
    const geo_basic<float_type,memory_type>& side_dom,
    const quadrature<float_type>&            quad) const
{
  _expr0.initialize (side_dom, quad);
  _expr1.initialize (side_dom, quad);
  _omega = side_dom.get_background_geo();
}
template<class Expr>
inline
bool
field_nonlinear_expr_dg<Expr>::initialize (
    const space_basic<float_type,memory_type>& Xh) const
{
  _expr0.initialize (Xh);
  _expr1.initialize (Xh);
  _omega = Xh.get_geo().get_background_geo();
  return true;
}
template<class Expr>
template<class ValueType>
inline
void 
field_nonlinear_expr_dg<Expr>::valued_check() const 
{
  space_constant::valued_type valued_tag = space_constant::valued_tag_traits<ValueType>::value;
  check_macro (_expr0.valued_tag() == valued_tag,
 	              "unexpected "<<  space_constant::valued_name(_expr0.valued_tag())
        << "-valued field while a " << space_constant::valued_name(valued_tag)
        << "-valued one is expected in expression");
}
template<class Expr>
template<class Result>
void
field_nonlinear_expr_dg<Expr>::evaluate (
    const geo_element&   K,
    std::vector<Result>& value) const
{
  size_type L_map_d = K.dimension() + 1;
  size_type L_dis_ie0, L_dis_ie1;
  side_information_type sid0, sid1;

  L_dis_ie0 = K.master(0);
  L_dis_ie1 = K.master(1);
  check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
      "unexpected isolated mesh side K="<<K);
  if (L_dis_ie1 == std::numeric_limits<size_type>::max()) {
    // K is a boundary side
    const geo_element& L0 = _omega.dis_get_geo_element (L_map_d, L_dis_ie0);
    L0.get_side_informations (K, sid0);
    _expr0.evaluate_on_side (L0, sid0, value);
    // average (i.e. _c0==0.5): fix it on the boundary where c0=1 : average(v)=v on the boundary
    Float c0 = (_c0 != 0.5) ? _c0 : 1;
    for (size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
      value[loc_idof] = c0*value[loc_idof];
    }
    return;
  }
  // K is an internal side
  std::vector<Result> value1;
  const geo_element& L0 = _omega.dis_get_geo_element (L_map_d, L_dis_ie0);
  const geo_element& L1 = _omega.dis_get_geo_element (L_map_d, L_dis_ie1);
  L0.get_side_informations (K, sid0);
  L1.get_side_informations (K, sid1);
  _expr0.evaluate_on_side (L0, sid0, value);
  _expr1.evaluate_on_side (L1, sid1, value1);
  for (size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
    value[loc_idof] = _c0*value[loc_idof] + _c1*value1[loc_idof];
  }
}

} // namespace rheolef
#endif // _RHEOLEF_FIELD_nONLINEAR_EXPR_H