447-dealII Report

**447-dealII**
Exec Time (%)	Codelet Name	Error (%)
20.4	mapping_q1__ZNK9MappingQ1ILi3EE12comput...	14.5
1.3	fe_poly__ZNK7FE_PolyI24TensorProductPol...	100.0
1.4	mapping_q__ZNK8MappingQILi3EE20apply_la...	16.2
2.6	sparse_matrix_double_include_lac_sparse...	2.5
11.4	fe__ZNK17FiniteElementBaseILi3EE11compu...	100.0
2.3	fe_poly__ZNK7FE_PolyI24TensorProductPol...	95.9
20.4	mapping_q1__ZNK9MappingQ1ILi3EE12comput...	14.5
1.4	mapping_q__ZNK8MappingQILi3EE20apply_la...	16.2
7.6	full_matrix_double_include_lac_full_mat...	2.3
1.0	step_14__usr_lib_gcc_x86_64_linux_gnu_4...	100.0

Exec Time : 3.343264e+10 cycles

Matching : 30.6%

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivofe_values_include_base_tensor_baseZN6TensorILi1ELi3EEpLERKS0448

Real time: 6.815835e+08 -- Predicted time: 1.007903e+10

Callcount: 28375652

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5578538	### 1	28391622.6666667	3.550000e+02	2.400000e+01	93.2

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q1_include_base_tensor_baseZN6TensorILi1ELi3EEaSERKS0398

Real time: 4.319898e+08 -- Predicted time: 6.303630e+09

Callcount: 16561775

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
3255978	### 1	14062807.1386129	4.070000e+02	2.400000e+01	94.1
5996782	### 2	1882970.80145248	2.570000e+02	3.600000e+01	85.9
7017149	### 3	628388.995525038	1.530000e+02	4.000000e+01	73.8

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivosparse_matrix_double_include_lac_sparse_matrix_templates__ZNK12SparseMatrixIdE5vmultI6VectorIdES3_EEvRT_RKT0__348

Real time: 8.875573e+08 -- Predicted time: 9.106381e+08

Callcount: 668

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
34	### 1	196.44175047068	8.501700e+04	8.710800e+04	2.4
86	### 2	45.0061481034199	3.706780e+05	3.643400e+05	1.7
127	### 3	53.2374927957941	5.833710e+05	5.691120e+05	2.4
205	### 4	54.02055341834	7.125660e+05	7.008080e+05	1.6
270	### 5	53.0298434026489	9.108890e+05	8.942680e+05	1.8
362	### 6	60.077322866498	1.386346e+06	1.364564e+06	1.5
423	### 7	64.0997806150684	2.035361e+06	2.023840e+06	0.5
518	### 8	72.0241569557941	2.993134e+06	2.980508e+06	0.4
595	### 9	80.0419319365667	4.123683e+06	3.870272e+06	6.1

                //----------------------------  polynomial.cc  -----------------------
//      $Id: polynomial.cc,v 1.3 2006/01/23 23:49:36 wolf Exp $   
//    Version: $Name:  $
//
//    Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//----------------------------   polynomial.cc  ----------------------


#include <base/polynomial.h>
#include <base/exceptions.h>
#include <base/thread_management.h>

#include <cmath>
#include <stdlib.h>


// have a lock that guarantees that at most one thread is changing and
// accessing the @p{coefficients} arrays of classes implementing
// polynomials with tables. make this lock local to this file.
//
// having only one lock for all of these classes is probably not going
// to be a problem since we only need it on very rare occasions. if
// someone finds this is a bottleneck, feel free to replace it by a
// more fine-grained solution
namespace 
{
  Threads::ThreadMutex coefficients_lock;
}



namespace Polynomials
{

// -------------------- class Polynomial ---------------- //


  template <typename number>
  Polynomial<number>::Polynomial (const std::vector<number> &a)
                  :
                  coefficients(a)
  {}



  template <typename number>
  number
  Polynomial<number>::value (const number x) const
  {
    Assert (coefficients.size() > 0, ExcVoidPolynomial());
    const unsigned int m=coefficients.size();

                                     // Horner scheme
    number value = coefficients.back();
    for (int k=m-2; k>=0; --k)
      value = value*x + coefficients[k];

    return value;
  }



  template <typename number>
  void
  Polynomial<number>::value (const number         x,
                             std::vector<number> &values) const
  {
    Assert (coefficients.size() > 0, ExcVoidPolynomial());
    Assert (values.size() > 0, ExcEmptyArray());
    const unsigned int values_size=values.size();
  
  
                                     // if we only need the value, then
                                     // call the other function since
                                     // that is significantly faster
                                     // (there is no need to allocate
                                     // and free memory, which is really
                                     // expensive compared to all the
                                     // other operations!)
    if (values_size == 1)
      {
        values[0] = value(x);
        return;
      };

                                     // if there are derivatives needed,
                                     // then do it properly by the
                                     // full Horner scheme
    const unsigned int m=coefficients.size();
    std::vector<number> a(coefficients);
    unsigned int j_faculty=1;

                                     // loop over all requested
                                     // derivatives. note that
                                     // derivatives @p{j>m} are
                                     // necessarily zero, as they
                                     // differentiate the polynomial
                                     // more often than the highest
                                     // power is
    const unsigned int min_valuessize_m=std::min(values_size, m);
    for (unsigned int j=0; j<min_valuessize_m; ++j)
      {
        for (int k=m-2; k>=static_cast<int>(j); --k)
          a[k]+=x*a[k+1];
        values[j]=j_faculty*a[j];

        j_faculty*=j+1;
      }

                                     // fill higher derivatives by zero
    for (unsigned int j=min_valuessize_m; j<values_size; ++j)
      values[j] = 0;
  }


  template <typename number>
  void
  Polynomial<number>::scale (std::vector<number> &coefficients,
                             const number         factor)
  {
    double f = 1.;
    for (typename std::vector<number>::iterator c = coefficients.begin();
         c != coefficients.end(); ++c)
      {
        *c *= f;
        f *= factor;
      }  
  }



  template <typename number>
  void
  Polynomial<number>::scale (const number factor)
  {
    scale (coefficients, factor);
  }



  template <typename number>
  void
  Polynomial<number>::multiply (std::vector<number> &coefficients,
                                const number         factor)
  {
    for (typename std::vector<number>::iterator c = coefficients.begin();
         c != coefficients.end(); ++c)
      *c *= factor;
  }



  template <typename number>
  Polynomial<number>&
  Polynomial<number>::operator *= (const double s)
  {
    for (typename std::vector<number>::iterator c = coefficients.begin();
         c != coefficients.end(); ++c)
      *c *= s;
    return *this;
  }

  
  template <typename number>
  Polynomial<number>&
  Polynomial<number>::operator += (const Polynomial<number>& p)
  {
				     // if necessary expand the number
				     // of coefficients we store
    if (p.coefficients.size() > coefficients.size())
      coefficients.resize (p.coefficients.size(), 0.);

    for (unsigned int i=0; i<p.coefficients.size(); ++i)
      coefficients[i] += p.coefficients[i];

    return *this;
  }

  
  template <typename number>
  Polynomial<number>&
  Polynomial<number>::operator -= (const Polynomial<number>& p)
  {
				     // if necessary expand the number
				     // of coefficients we store
    if (p.coefficients.size() > coefficients.size())
      coefficients.resize (p.coefficients.size(), 0.);

    for (unsigned int i=0; i<p.coefficients.size(); ++i)
      coefficients[i] -= p.coefficients[i];

    return *this;
  }

  
  template <typename number>
  template <typename number2>
  void
  Polynomial<number>::shift(std::vector<number>& coefficients,
                            const number2 offset)
  {  
                                     // Copy coefficients to a vector of
                                     // accuracy given by the argument
    std::vector<number2> new_coefficients(coefficients.begin(),
                                          coefficients.end());
  
                                     // Traverse all coefficients from
                                     // c_1. c_0 will be modified by
                                     // higher degrees, only.
    for (unsigned int d=1; d<new_coefficients.size(); ++d)
      {
        const unsigned int n = d;
                                         // Binomial coefficients are
                                         // needed for the
                                         // computation. The rightmost
                                         // value is unity.
        unsigned int binomial_coefficient = 1;

                                         // Powers of the offset will be
                                         // needed and computed
                                         // successively.
        number2 offset_power = offset;
      
                                         // Compute (x+offset)^d
                                         // and modify all values c_k
                                         // with k<d.
                                         // The coefficient in front of
                                         // x^d is not modified in this step.
        for (unsigned int k=0;k<d;++k)
          {
                                             // Recursion from Bronstein
                                             // Make sure no remainders
                                             // occur in integer
                                             // division.
            binomial_coefficient = (binomial_coefficient*(n-k))/(k+1);

            new_coefficients[d-k-1] += new_coefficients[d]
                                       * binomial_coefficient
                                       * offset_power;
            offset_power *= offset;
          }
                                         // The binomial coefficient
                                         // should have gone through a
                                         // whole row of Pascal's
                                         // triangle.
        Assert (binomial_coefficient == 1, ExcInternalError());
      }

                                     // copy new elements to old vector
    coefficients.assign(new_coefficients.begin(), new_coefficients.end());
  }


  template <typename number>
  template <typename number2>
  void
  Polynomial<number>::shift (const number2 offset)
  {
    shift (coefficients, offset);
  }


  
  template <typename number>
  Polynomial<number>
  Polynomial<number>::derivative () const
  {
    if (degree() == 0)
      return Monomial<number>(0, 0.);

    std::vector<number> newcoefficients (coefficients.size()-1);
    for (unsigned int i=1 ; i<coefficients.size() ; ++i)
      newcoefficients[i-1] = i * coefficients[i];

    return Polynomial<number> (newcoefficients);
  }
  

  template <typename number>
  void
  Polynomial<number>::print (std::ostream& out) const
  {
    for (int i=degree();i>=0;--i)
      {
        out << static_cast<double>(coefficients[i])
            << " x^" << i << std::endl;
      }
  }


// ------------------ class Monomial -------------------------- //

  template <typename number>
  std::vector<number>
  Monomial<number>::make_vector(unsigned int n,
				double coefficient)
  {
    std::vector<number> result(n+1, 0.);
    result[n] = coefficient;
    return result;
  }
  
  
  template <typename number>
  Monomial<number>::Monomial (unsigned int n,
			      double coefficient)
		  : Polynomial<number>(make_vector(n, coefficient))
  {}
  
  
  template <typename number>
  std::vector<Polynomial<number> >
  Monomial<number>::generate_complete_basis (const unsigned int degree)
  {
    std::vector<Polynomial<number> > v;
    v.reserve(degree+1);
    for (unsigned int i=0; i<=degree; ++i)
      v.push_back (Monomial<number>(i));
    return v;
  }
// ------------------ class LagrangeEquidistant --------------- //

  LagrangeEquidistant::LagrangeEquidistant (const unsigned int n,
                                            const unsigned int support_point)
		  :
                  Polynomial<double>(compute_coefficients(n,support_point))
  {}



  std::vector<double> 
  LagrangeEquidistant::compute_coefficients (const unsigned int n,
                                             const unsigned int support_point)
  {
    std::vector<double> a (n+1);
    Assert(support_point<n+1, ExcIndexRange(support_point, 0, n+1));

    unsigned int n_functions=n+1;
    Assert(support_point<n_functions,
           ExcIndexRange(support_point, 0, n_functions));
    double const *x=0;
  
    switch (n)
      {
        case 1:
        {
          static const double x1[4]=
            {
                  1.0, -1.0,
                  0.0, 1.0
            };
          x=&x1[0];
          break;	
        }
        case 2:
        {
          static const double x2[9]=
            {
                  1.0, -3.0, 2.0,
                  0.0, 4.0, -4.0,
                  0.0, -1.0, 2.0
            };
          x=&x2[0];
          break;
        }
        case 3:
        {
          static const double x3[16]=
            {
                  1.0, -11.0/2.0, 9.0, -9.0/2.0,
                  0.0, 9.0, -45.0/2.0, 27.0/2.0,
                  0.0, -9.0/2.0, 18.0, -27.0/2.0,
                  0.0, 1.0, -9.0/2.0, 9.0/2.0
            };
          x=&x3[0];
          break;
        }
        case 4:
        {
          static const double x4[25]=
            {
                  1.0, -25.0/3.0, 70.0/3.0, -80.0/3.0, 32.0/3.0,
                  0.0, 16.0, -208.0/3.0, 96.0, -128.0/3.0,
                  0.0, -12.0, 76.0, -128.0, 64.0,
                  0.0, 16.0/3.0, -112.0/3.0, 224.0/3.0, -128.0/3.0,
                  0.0, -1.0, 22.0/3.0, -16.0, 32.0/3.0
            };
          x=&x4[0];
          break;	
        }
        case 5:
        {
          static const double x5[36]=
            {
                  1.0, -137.0/12.0, 375.0/8.0, -2125.0/24.0, 625.0/8.0, -625.0/24.0,
                  0.0, 25.0, -1925.0/12.0, 8875.0/24.0, -4375.0/12.0, 3125.0/24.0,
                  0.0, -25.0, 2675.0/12.0, -7375.0/12.0, 8125.0/12.0, -3125.0/12.0,
                  0.0, 50.0/3.0, -325.0/2.0, 6125.0/12.0, -625.0, 3125.0/12.0,
                  0.0, -25.0/4.0, 1525.0/24.0, -5125.0/24.0, 6875.0/24.0, -3125.0/24.0,
                  0.0, 1.0, -125.0/12.0, 875.0/24.0, -625.0/12.0, 625.0/24.0
            };
          x=&x5[0];
          break;
        }
        case 6:
        {
          static const double x6[49]=
            {
                  1.0, -147.0/10.0, 406.0/5.0, -441.0/2.0, 315.0, -1134.0/5.0,
                  324.0/5.0, 0.0, 36.0, -1566.0/5.0, 1044.0, -1674.0, 1296.0,
                  -1944.0/5.0, 0.0, -45.0, 1053.0/2.0, -4149.0/2.0, 3699.0, -3078.0,
                  972.0, 0.0, 40.0, -508.0, 2232.0, -4356.0, 3888.0, -1296.0, 0.0,
                  -45.0/2.0, 297.0, -2763.0/2.0, 2889.0, -2754.0, 972.0, 0.0,
                  36.0/5.0, -486.0/5.0, 468.0, -1026.0, 5184.0/5.0, -1944.0/5.0, 0.0,
                  -1.0, 137.0/10.0, -135.0/2.0, 153.0, -162.0, 324.0/5.0
            };
          x=&x6[0];
          break;
        }
        case 7:
        {
          static const double x7[64]=
            {
                  1.0, -363.0/20.0, 22981.0/180.0, -331681.0/720.0, 16807.0/18.0,
                  -386561.0/360.0, 117649.0/180.0, -117649.0/720.0, 0.0, 49.0,
                  -10927.0/20.0, 109417.0/45.0, -88837.0/16.0, 991613.0/144.0,
                  -352947.0/80.0, 823543.0/720.0, 0.0, -147.0/2.0, 43071.0/40.0,
                  -1347647.0/240.0, 170471.0/12.0, -151263.0/8.0, 1529437.0/120.0,
                  -823543.0/240.0, 0.0, 245.0/3.0, -46501.0/36.0, 133427.0/18.0,
                  -2926819.0/144.0, 4151329.0/144.0, -2941225.0/144.0,
                  823543.0/144.0, 0.0, -245.0/4.0, 2009.0/2.0, -872935.0/144.0,
                  52822.0/3.0, -1899191.0/72.0, 117649.0/6.0, -823543.0/144.0, 0.0,
                  147.0/5.0, -9849.0/20.0, 45962.0/15.0, -444185.0/48.0,
                  1159683.0/80.0, -2705927.0/240.0, 823543.0/240.0, 0.0, -49.0/6.0,
                  49931.0/360.0, -634207.0/720.0, 98441.0/36.0, -319333.0/72.0,
                  1294139.0/360.0, -823543.0/720.0, 0.0, 1.0, -343.0/20.0,
                  9947.0/90.0, -16807.0/48.0, 84035.0/144.0, -117649.0/240.0,
                  117649.0/720.0
            };
          x=&x7[0];
          break;
        }
        case 8:
        {
          static const double x8[81]=
            {
                  1.0, -761.0/35.0, 59062.0/315.0, -4272.0/5.0, 34208.0/15.0,
                  -18432.0/5.0, 53248.0/15.0, -65536.0/35.0, 131072.0/315.0, 0.0,
                  64.0, -30784.0/35.0, 44672.0/9.0, -673792.0/45.0, 235520.0/9.0,
                  -1196032.0/45.0, 131072.0/9.0, -1048576.0/315.0, 0.0, -112.0,
                  9936.0/5.0, -587296.0/45.0, 1956992.0/45.0, -733184.0/9.0,
                  3915776.0/45.0, -2228224.0/45.0, 524288.0/45.0, 0.0, 448.0/3.0,
                  -128192.0/45.0, 102016.0/5.0, -1097728.0/15.0, 145408.0,
                  -2441216.0/15.0, 1441792.0/15.0, -1048576.0/45.0, 0.0, -140.0,
                  2764.0, -186496.0/9.0, 703552.0/9.0, -1466368.0/9.0, 1712128.0/9.0,
                  -1048576.0/9.0, 262144.0/9.0, 0.0, 448.0/5.0, -9024.0/5.0,
                  626048.0/45.0, -2443264.0/45.0, 5285888.0/45.0, -6406144.0/45.0,
                  4063232.0/45.0, -1048576.0/45.0, 0.0, -112.0/3.0, 34288.0/45.0,
                  -5984.0, 358784.0/15.0, -53248.0, 999424.0/15.0, -131072.0/3.0,
                  524288.0/45.0, 0.0, 64.0/7.0, -6592.0/35.0, 67456.0/45.0,
                  -274432.0/45.0, 124928.0/9.0, -802816.0/45.0, 3801088.0/315.0,
                  -1048576.0/315.0, 0.0, -1.0, 726.0/35.0, -7504.0/45.0, 30944.0/45.0,
                  -14336.0/9.0, 94208.0/45.0, -65536.0/45.0, 131072.0/315.0
            };
          x=&x8[0];
          break;
        }
        case 9:
        {
          static const double x9[100]=
            {
                  1.0, -7129.0/280.0, 58635.0/224.0, -40707.0/28.0, 623295.0/128.0,
                  -6589431.0/640.0, 885735.0/64.0, -5137263.0/448.0, 4782969.0/896.0,
                  -4782969.0/4480.0, 0.0, 81.0, -373329.0/280.0, 10307331.0/1120.0,
                  -5589243.0/160.0, 51221727.0/640.0, -4546773.0/40.0,
                  31355019.0/320.0, -52612659.0/1120.0, 43046721.0/4480.0, 0.0,
                  -162.0, 475389.0/140.0, -15190173.0/560.0, 18152829.0/160.0,
                  -44529507.0/160.0, 33244587.0/80.0, -3720087.0/10.0,
                  205667667.0/1120.0, -43046721.0/1120.0, 0.0, 252.0, -56601.0/10.0,
                  1959363.0/40.0, -8776431.0/40.0, 91020753.0/160.0,
                  -71035947.0/80.0, 16474671.0/20.0, -33480783.0/80.0,
                  14348907.0/160.0, 0.0, -567.0/2.0, 526419.0/80.0, -4752351.0/80.0,
                  89119521.0/320.0, -241241409.0/320.0, 195629337.0/160.0,
                  -187598673.0/160.0, 196101729.0/320.0, -43046721.0/320.0, 0.0,
                  1134.0/5.0, -21465.0/4.0, 795339.0/16.0, -3844017.0/16.0,
                  215023653.0/320.0, -18009945.0/16.0, 35606547.0/32.0,
                  -4782969.0/8.0, 43046721.0/320.0, 0.0, -126.0, 60381.0/20.0,
                  -2276289.0/80.0, 22480173.0/160.0, -64448703.0/160.0,
                  55447011.0/80.0, -28166373.0/40.0, 62178597.0/160.0,
                  -14348907.0/160.0, 0.0, 324.0/7.0, -78327.0/70.0, 2989629.0/280.0,
                  -2142531.0/40.0, 25043337.0/160.0, -22025277.0/80.0,
                  80247591.0/280.0, -90876411.0/560.0, 43046721.0/1120.0, 0.0,
                  -81.0/8.0, 275967.0/1120.0, -1328967.0/560.0, 7712091.0/640.0,
                  -22878207.0/640.0, 20490003.0/320.0, -21789081.0/320.0,
                  176969853.0/4480.0, -43046721.0/4480.0, 0.0, 1.0, -6849.0/280.0,
                  265779.0/1120.0, -194643.0/160.0, 2337903.0/640.0, -531441.0/80.0,
                  2302911.0/320.0, -4782969.0/1120.0, 4782969.0/4480.0
            };
          x=&x9[0];
          break;
        }
        case 10:
        {
          static const double x10[121]=
            {
                  1.0, -7381.0/252.0, 177133.0/504.0, -10511875.0/4536.0,
                  42711625.0/4536.0, -5369375.0/216.0, 4695625.0/108.0,
                  -9453125.0/189.0, 6875000.0/189.0, -8593750.0/567.0,
                  1562500.0/567.0, 0.0, 100.0, -121525.0/63.0, 1997825.0/126.0,
                  -82992625.0/1134.0, 3775625.0/18.0, -20965625.0/54.0,
                  4187500.0/9.0, -65937500.0/189.0, 3125000.0/21.0,
                  -15625000.0/567.0, 0.0, -225.0, 153025.0/28.0, -2898075.0/56.0,
                  33095875.0/126.0, -57981875.0/72.0, 56396875.0/36.0,
                  -17546875.0/9.0, 94843750.0/63.0, -41406250.0/63.0, 7812500.0/63.0,
                  0.0, 400.0, -654100.0/63.0, 20028950.0/189.0, -108434750.0/189.0,
                  16686250.0/9.0, -33868750.0/9.0, 43625000.0/9.0, -242500000.0/63.0,
                  325000000.0/189.0, -62500000.0/189.0, 0.0, -525.0, 168775.0/12.0,
                  -1792225.0/12.0, 91073375.0/108.0, -102070625.0/36.0,
                  107321875.0/18.0, -71281250.0/9.0, 19375000.0/3.0, -26562500.0/9.0,
                  15625000.0/27.0, 0.0, 504.0, -13754.0, 149625.0, -7818625.0/9.0,
                  27074375.0/9.0, -58608125.0/9.0, 80000000.0/9.0, -66875000.0/9.0,
                  31250000.0/9.0, -6250000.0/9.0, 0.0, -350.0, 174025.0/18.0,
                  -11544725.0/108.0, 34178875.0/54.0, -80666875.0/36.0,
                  89384375.0/18.0, -62468750.0/9.0, 5937500.0, -76562500.0/27.0,
                  15625000.0/27.0, 0.0, 1200.0/7.0, -100300.0/21.0, 1121950.0/21.0,
                  -60659750.0/189.0, 10401250.0/9.0, -7831250.0/3.0,
                  234625000.0/63.0, -205000000.0/63.0, 100000000.0/63.0,
                  -62500000.0/189.0, 0.0, -225.0/4.0, 88325.0/56.0, -996675.0/56.0,
                  54486625.0/504.0, -28405625.0/72.0, 32584375.0/36.0,
                  -11828125.0/9.0, 73750000.0/63.0, -36718750.0/63.0, 7812500.0/63.0,
                  0.0, 100.0/9.0, -6575.0/21.0, 4033825.0/1134.0, -24717625.0/1134.0,
                  4341875.0/54.0, -10090625.0/54.0, 7437500.0/27.0,
                  -47187500.0/189.0, 71875000.0/567.0, -15625000.0/567.0, 0.0, -1.0,
                  7129.0/252.0, -162875.0/504.0, 1130750.0/567.0, -59375.0/8.0,
                  1883125.0/108.0, -78125.0/3.0, 4531250.0/189.0, -781250.0/63.0,
                  1562500.0/567.0
            };
          x=&x10[0];
          break;
        }
        default:
              Assert(false, ExcNotImplemented());
      }

    Assert(x!=0, ExcInternalError());
    for (unsigned int i=0; i<n_functions; ++i)
      a[i]=x[support_point*n_functions+i];
  
    return a;
  }


  std::vector<Polynomial<double> >
  LagrangeEquidistant::
  generate_complete_basis (const unsigned int degree)
  {
    if (degree==0)
                                       // create constant polynomial
      return std::vector<Polynomial<double> >
        (1, Polynomial<double> (std::vector<double> (1,1.)));
    else
      {
                                         // create array of Lagrange
                                         // polynomials
        std::vector<Polynomial<double> > v;
        for (unsigned int i=0; i<=degree; ++i)
          v.push_back(LagrangeEquidistant(degree,i));
        return v;
      };
  }



// ------------------ class Legendre --------------- //


//TODO:[?] This class leaks memory, but only at the very end of a program.
// Since it expands the Legendre<number>::coefficients array, the elements
// of this static variable are not destroyed at the end of the program
// run. While this is not a problem (since the returned memory could
// not be used anyway then), it is a little confusing when looking at
// a memory checker such as "purify". Maybe, this should be handled somehow
// to avoid this confusion in future.

// Reserve space for polynomials up to degree 19. Should be sufficient
// for the start.
  std::vector<const std::vector<double> *>
  Legendre::recursive_coefficients(20,
				   static_cast<const std::vector<double>*>(0));
  std::vector<const std::vector<double> *>
  Legendre::shifted_coefficients(20,
				 static_cast<const std::vector<double>*>(0));


  Legendre::Legendre (const unsigned int k)
                  :
                  Polynomial<double> (get_coefficients(k))
  {}


  
  void
  Legendre::compute_coefficients (const unsigned int k_)
  {
                                     // this function is unused in SPEC and is
                                     // removed because the original required
                                     // work arounds for different compilers
    abort ();
  }



  const std::vector<double> &
  Legendre::get_coefficients (const unsigned int k)
  {
                                     // first make sure the coefficients
                                     // get computed if so necessary
    compute_coefficients (k);

                                     // then get a pointer to the array
                                     // of coefficients. do that in a MT
                                     // safe way
    Threads::ThreadMutex::ScopedLock lock (coefficients_lock);
    return *shifted_coefficients[k];
  }



  std::vector<Polynomial<double> >
  Legendre::generate_complete_basis (const unsigned int degree)
  {
    std::vector<Polynomial<double> > v;
    v.reserve(degree+1);
    for (unsigned int i=0; i<=degree; ++i)
      v.push_back (Legendre(i));
    return v;
  }



// ------------------ class Hierarchical --------------- //



// Reserve space for polynomials up to degree 19. Should be sufficient
// for the start.
  std::vector<const std::vector<double> *>
  Hierarchical::recursive_coefficients(
     20, static_cast<const std::vector<double>*>(0));



  Hierarchical::Hierarchical (const unsigned int k)
                      :
                      Polynomial<double> (get_coefficients(k))
  {}



  void
  Hierarchical::compute_coefficients (const unsigned int k_)
  {
    unsigned int k = k_;

                                     // first make sure that no other
                                     // thread intercepts the operation
                                     // of this function
                                     // for this, acquire the lock
                                     // until we quit this function
    Threads::ThreadMutex::ScopedLock lock(coefficients_lock);

                                     // The first 2 coefficients 
                                     // are hard-coded
    if (k==0)
      k=1;
                                     // check: does the information
                                     // already exist?
    if (  (recursive_coefficients.size() < k+1) ||
	  ((recursive_coefficients.size() >= k+1) && 
           (recursive_coefficients[k] == 0)) )
				           // no, then generate the
				           // respective coefficients
      {
	recursive_coefficients.resize (k+1, 0);
      
	if (k<=1)
	  {
                                             // create coefficients
                                             // vectors for k=0 and k=1
                                             //
                                             // allocate the respective
                                             // amount of memory and
                                             // later assign it to the
                                             // coefficients array to
                                             // make it const
	    std::vector<double> *c0 = new std::vector<double>(2);
	    (*c0)[0] =  1.;
	    (*c0)[1] = -1.;

	    std::vector<double> *c1 = new std::vector<double>(2);
	    (*c1)[0] = 0.;
	    (*c1)[1] = 1.;

                                             // now make these arrays
                                             // const
	    recursive_coefficients[0] = c0;
	    recursive_coefficients[1] = c1;
	  }
	else if (k==2)
	  {
            coefficients_lock.release ();
	    compute_coefficients(1);
	    coefficients_lock.acquire ();

	    std::vector<double> *c2 = new std::vector<double>(3);

	    const double a = 1.; //1./8.;

	    (*c2)[0] =   0.*a;
	    (*c2)[1] =  -4.*a;
	    (*c2)[2] =   4.*a;
	    
	    recursive_coefficients[2] = c2;
	  }
	else
	  {
                                             // for larger numbers,
                                             // compute the coefficients
                                             // recursively. to do so,
                                             // we have to release the
                                             // lock temporarily to
                                             // allow the called
                                             // function to acquire it
                                             // itself
	    coefficients_lock.release ();
	    compute_coefficients(k-1);
	    coefficients_lock.acquire ();

	    std::vector<double> *ck = new std::vector<double>(k+1);
	   
	    const double a = 1.; //1./(2.*k);

	    (*ck)[0] = - a*(*recursive_coefficients[k-1])[0];
	  
	    for (unsigned int i=1; i<=k-1; ++i)
		(*ck)[i] = a*( 2.*(*recursive_coefficients[k-1])[i-1]
			       - (*recursive_coefficients[k-1])[i] );
	  
	    (*ck)[k] = a*2.*(*recursive_coefficients[k-1])[k-1];
	                                  // for even degrees, we need
	                                  // to add a multiple of
	                                  // basis fcn phi_2
	    if ( (k%2) == 0 )
	      {
		double b = 1.; //8.;
		//for (unsigned int i=1; i<=k; i++)
		//  b /= 2.*i;

		(*ck)[1] += b*(*recursive_coefficients[2])[1];
		(*ck)[2] += b*(*recursive_coefficients[2])[2];
	      }	  
                                             // finally assign the newly
                                             // created vector to the
                                             // const pointer in the
                                             // coefficients array
	    recursive_coefficients[k] = ck;
	  };
      };
  }



  const std::vector<double> &
  Hierarchical::get_coefficients (const unsigned int k)
  {
				   // first make sure the coefficients
				   // get computed if so necessary
    compute_coefficients (k);

				   // then get a pointer to the array
				   // of coefficients. do that in a MT
                                     // safe way
    coefficients_lock.acquire ();
    const std::vector<double> *p = recursive_coefficients[k];
    coefficients_lock.release ();

				   // return the object pointed
				   // to. since this object does not
				   // change any more once computed,
				   // this is MT safe
    return *p;
  }



  std::vector<Polynomial<double> >
  Hierarchical::generate_complete_basis (const unsigned int degree)
  {
    if (degree==0)
                                       // create constant
                                       // polynomial. note that we
                                       // can't use the other branch
                                       // of the if-statement, since
                                       // calling the constructor of
                                       // this class with argument
                                       // zero does _not_ create the
                                       // constant polynomial, but
                                       // rather 1-x
      return std::vector<Polynomial<double> >
        (1, Polynomial<double> (std::vector<double> (1,1.)));
    else
      {
        std::vector<Polynomial<double> > v;
        v.reserve(degree+1);
        for (unsigned int i=0; i<=degree; ++i)
          v.push_back (Hierarchical(i));
        return v;
      }
  }
}



// ------------------ explicit instantiations --------------- //

namespace Polynomials
{
  template class Polynomial<float>;
  template class Polynomial<double>;
  template class Polynomial<long double>;

  template void Polynomial<float>::shift(const float offset);
  template void Polynomial<float>::shift(const double offset);
  template void Polynomial<double>::shift(const double offset);
  template void Polynomial<long double>::shift(const long double offset);
  template void Polynomial<float>::shift(const long double offset);
  template void Polynomial<double>::shift(const long double offset);

  template class Monomial<float>;
  template class Monomial<double>;
  template class Monomial<long double>;
}

invivopolynomial__ZNK11Polynomials10PolynomialIdE5valueEd_61

Real time: 1.230907e+09 -- Predicted time: 1.202738e+10

Callcount: 43046364

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
8462740	### 1	32072995.9571816	2.820000e+02	2.400000e+01	91.4
9168246	### 2	11006641.7509872	2.710000e+02	2.800000e+01	89.6

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivofe_values_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470

Real time: 6.813907e+08 -- Predicted time: 1.391051e+10

Callcount: 28375652

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5578538	### 1	28384365.4820374	4.900000e+02	2.400000e+01	95.1
13323871	### 2	5353.44396795557	4.060000e+02	2.800000e+01	93.1

                //----------------------------  mapping_q1.cc  ---------------------------
//    $Id: mapping_q1.cc,v 1.1 2004/09/14 00:51:24 wolf Exp $
//    Version: $Name:  $
//
//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//----------------------------  mapping_q1.cc  ---------------------------

#include <base/tensor.h>
#include <base/quadrature.h>
#include <base/memory_consumption.h>
#include <lac/full_matrix.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <fe/fe_values.h>
#include <fe/mapping_q1.h>

#include <cmath>
#include <algorithm>
#include <memory>




template <int dim>
const unsigned int MappingQ1<dim>::n_shape_functions;



template<int dim>
MappingQ1<dim>::InternalData::InternalData (const unsigned int n_shape_functions)
		:
		is_mapping_q1_data(true),
		n_shape_functions (n_shape_functions)
{}



template <int dim>
unsigned int
MappingQ1<dim>::InternalData::memory_consumption () const
{
  return (Mapping<dim>::InternalDataBase::memory_consumption() +
	  MemoryConsumption::memory_consumption (shape_values) +
	  MemoryConsumption::memory_consumption (shape_derivatives) +
	  MemoryConsumption::memory_consumption (covariant) +
	  MemoryConsumption::memory_consumption (contravariant) +
	  MemoryConsumption::memory_consumption (unit_tangentials) +
	  MemoryConsumption::memory_consumption (aux) +
	  MemoryConsumption::memory_consumption (mapping_support_points) +
	  MemoryConsumption::memory_consumption (cell_of_current_support_points) +
	  MemoryConsumption::memory_consumption (is_mapping_q1_data) +
	  MemoryConsumption::memory_consumption (n_shape_functions));
}



template <int dim>
MappingQ1<dim>::MappingQ1 ()
{}



template<int dim>
void
MappingQ1<dim>::compute_shapes (const std::vector<Point<dim> > &unit_points,
				InternalData &data) const
{
  if (data.is_mapping_q1_data)
    MappingQ1<dim>::compute_shapes_virtual(unit_points, data);
  else
    compute_shapes_virtual(unit_points, data);
}


#if (deal_II_dimension == 1)

template<>
const unsigned int MappingQ1<deal_II_dimension>::vertex_mapping[2] =
{ 0, 1 };


template<>
void
MappingQ1<1>::compute_shapes_virtual (const std::vector<Point<1> > &unit_points,
				      InternalData& data) const
{
  const unsigned int n_points=unit_points.size();
  for (unsigned int k = 0 ; k < n_points ; ++k)
    {
      double x = unit_points[k](0);
      
      if (data.shape_values.size()!=0)
	{
	  Assert(data.shape_values.size()==n_shape_functions*n_points,
		 ExcInternalError());
	  data.shape(k,0) = 1.-x;
	  data.shape(k,1) = x;
	}
      if (data.shape_derivatives.size()!=0)
	{
	  Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
		 ExcInternalError());
	  data.derivative(k,0)[0] = -1.;
	  data.derivative(k,1)[0] = 1.;
	}
    }
}

#endif

#if (deal_II_dimension == 2)

template<> const unsigned int
MappingQ1<2>::vertex_mapping[4] =
{ 0, 1, 3, 2 };


template<>
void
MappingQ1<2>::compute_shapes_virtual (const std::vector<Point<2> > &unit_points,
				      InternalData &data) const
{
  const unsigned int n_points=unit_points.size();
  for (unsigned int k = 0 ; k < n_points ; ++k)
    {
      double x = unit_points[k](0);
      double y = unit_points[k](1);
      
      if (data.shape_values.size()!=0)
	{
	  Assert(data.shape_values.size()==n_shape_functions*n_points,
		 ExcInternalError());
	  data.shape(k,0) = (1.-x)*(1.-y);
	  data.shape(k,1) = x*(1.-y);
	  data.shape(k,2) = (1.-x)*y;
	  data.shape(k,3) = x*y;
	}
      if (data.shape_derivatives.size()!=0)
	{
	  Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
		 ExcInternalError());
	  data.derivative(k,0)[0] = (y-1.);
	  data.derivative(k,1)[0] = (1.-y);
	  data.derivative(k,2)[0] = -y;
	  data.derivative(k,3)[0] = y;
	  data.derivative(k,0)[1] = (x-1.);
	  data.derivative(k,1)[1] = -x;
	  data.derivative(k,2)[1] = (1.-x);
	  data.derivative(k,3)[1] = x;
	}
    }
}

#endif

#if (deal_II_dimension == 3)

template<>
const unsigned int MappingQ1<deal_II_dimension>::vertex_mapping[8] =
{ 0, 1, 4, 5, 3, 2, 7, 6 };


template<>
void
MappingQ1<3>::compute_shapes_virtual (const std::vector<Point<3> > &unit_points,
				      InternalData &data) const
{
  const unsigned int n_points=unit_points.size();
  for (unsigned int k = 0 ; k < n_points ; ++k)
    {
      double x = unit_points[k](0);
      double y = unit_points[k](1);
      double z = unit_points[k](2);
      
      if (data.shape_values.size()!=0)
	{
	  Assert(data.shape_values.size()==n_shape_functions*n_points,
		 ExcInternalError());
	  data.shape(k,0) = (1.-x)*(1.-y)*(1.-z);
	  data.shape(k,1) = x*(1.-y)*(1.-z);
	  data.shape(k,2) = (1.-x)*y*(1.-z);
	  data.shape(k,3) = x*y*(1.-z);
	  data.shape(k,4) = (1.-x)*(1.-y)*z;
	  data.shape(k,5) = x*(1.-y)*z;
	  data.shape(k,6) = (1.-x)*y*z;
	  data.shape(k,7) = x*y*z;
	}
      if (data.shape_derivatives.size()!=0)
	{
	  Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
		 ExcInternalError());
	  data.derivative(k,0)[0] = (y-1.)*(1.-z);
	  data.derivative(k,1)[0] = (1.-y)*(1.-z);
	  data.derivative(k,2)[0] = -y*(1.-z);
	  data.derivative(k,3)[0] = y*(1.-z);
	  data.derivative(k,4)[0] = (y-1.)*z;
	  data.derivative(k,5)[0] = (1.-y)*z;
	  data.derivative(k,6)[0] = -y*z;
	  data.derivative(k,7)[0] = y*z;
	  data.derivative(k,0)[1] = (x-1.)*(1.-z);
	  data.derivative(k,1)[1] = -x*(1.-z);
	  data.derivative(k,2)[1] = (1.-x)*(1.-z);
	  data.derivative(k,3)[1] = x*(1.-z);
	  data.derivative(k,4)[1] = (x-1.)*z;
	  data.derivative(k,5)[1] = -x*z;
	  data.derivative(k,6)[1] = (1.-x)*z;
	  data.derivative(k,7)[1] = x*z;
	  data.derivative(k,0)[2] = (x-1)*(1.-y);
	  data.derivative(k,1)[2] = x*(y-1.);
	  data.derivative(k,2)[2] = (x-1.)*y;
	  data.derivative(k,3)[2] = -x*y;
	  data.derivative(k,4)[2] = (1.-x)*(1.-y);
	  data.derivative(k,5)[2] = x*(1.-y);
	  data.derivative(k,6)[2] = (1.-x)*y;
	  data.derivative(k,7)[2] = x*y;
	}
    }
}

#endif

template <int dim>
UpdateFlags
MappingQ1<dim>::update_once (const UpdateFlags in) const
{
  UpdateFlags out = UpdateFlags(in & (update_transformation_values
				      | update_transformation_gradients));

				   // Shape function values
  if (in & update_q_points)
    out |= update_transformation_values;

				   // Shape function gradients
  if (in & (update_covariant_transformation
	    | update_contravariant_transformation
	    | update_JxW_values
	    | update_boundary_forms
	    | update_normal_vectors))
    out |= update_transformation_gradients;

  return out;
}



template <int dim>
UpdateFlags
MappingQ1<dim>::update_each (const UpdateFlags in) const
{
				   // Select flags of concern for the
				   // transformation.
  UpdateFlags out = UpdateFlags(in & (update_q_points
				      | update_covariant_transformation
				      | update_contravariant_transformation
				      | update_JxW_values
				      | update_boundary_forms
				      | update_normal_vectors));

				   // The following is a little incorrect:
				   // If not applied on a face,
				   // update_boundary_forms does not
				   // make sense. On the other hand,
				   // it is necessary on a
				   // face. Currently,
				   // update_boundary_forms is simply
				   // ignored for the interior of a
				   // cell.
  if (out & (update_JxW_values
	     |update_normal_vectors))
    out |= update_boundary_forms;
  
  if (out & (update_covariant_transformation
	     | update_JxW_values
	     | update_boundary_forms
	     | update_normal_vectors))
    out |= update_contravariant_transformation;

  return out;
}


template <int dim>
void
MappingQ1<dim>::compute_data (const UpdateFlags      update_flags,
			      const Quadrature<dim> &q,
			      const unsigned int     n_original_q_points,
			      InternalData          &data) const
{
  const unsigned int n_q_points = q.n_quadrature_points;

  data.update_once = update_once(update_flags);
  data.update_each = update_each(update_flags);
  data.update_flags = data.update_once | data.update_each;

  const UpdateFlags flags(data.update_flags);
  
  if (flags & update_transformation_values)
    data.shape_values.resize(data.n_shape_functions * n_q_points);

  if (flags & update_transformation_gradients)
    data.shape_derivatives.resize(data.n_shape_functions * n_q_points);

  if (flags & update_covariant_transformation)
    data.covariant.resize(n_original_q_points);

  if (flags & update_contravariant_transformation)
    data.contravariant.resize(n_original_q_points);

  compute_shapes (q.get_points(), data);
}



template <int dim>
typename Mapping<dim>::InternalDataBase *
MappingQ1<dim>::get_data (const UpdateFlags update_flags,
			  const Quadrature<dim>& q) const
{
  InternalData* data = new InternalData(n_shape_functions);
  compute_data (update_flags, q, q.n_quadrature_points, *data);
  return data;
}


#if deal_II_dimension == 1

template <>
void
MappingQ1<1>::compute_face_data (const UpdateFlags update_flags,
                                 const Quadrature<1>& q,
                                 const unsigned int n_original_q_points,
                                 InternalData& data) const
{
  compute_data (update_flags, q, n_original_q_points, data);
}

#else

template <int dim>
void
MappingQ1<dim>::compute_face_data (const UpdateFlags update_flags,
				   const Quadrature<dim>& q,
				   const unsigned int n_original_q_points,
				   InternalData& data) const
{
  compute_data (update_flags, q, n_original_q_points, data);

  if (data.update_flags & update_boundary_forms)
    {
      data.aux.reinit(dim-1, n_original_q_points);      
      
				       // Compute tangentials to the
				       // unit cell.
      const unsigned int nfaces = GeometryInfo<dim>::faces_per_cell;
      data.unit_tangentials.reinit(nfaces*(dim-1),
				   n_original_q_points);
      for (unsigned int i=0; i<nfaces; ++i)
	{
					   // Base index of the
					   // non-zero entry of the
					   // tangential vector.  Add
					   // dim so we can subtract 1
					   // without getting negative
					   // values.
	  unsigned int nindex = this->normal_directions[i]/2 + dim;

					   // First tangential has a
					   // non-zero in component
					   // (i+1)%dim, if normal is
					   // non-zero in i.
	  Tensor<1,dim> tangential;
	  tangential[(nindex+1)%dim] = (this->normal_directions[i]%2) ? -1 : 1;
	  std::fill (data.unit_tangentials[i].begin(),
		     data.unit_tangentials[i].end(),
		     tangential);
	  
	  if (dim>2)
	    {
					       // Second tangential
					       // has a non-zero in
					       // component (i-1)%dim,
					       // if normal is
					       // non-zero in
					       // i. Creates a
					       // right-handed system.
	      Tensor<1,dim> tangential;
	      tangential[(nindex-1)%dim] = 1.;
	      std::fill (data.unit_tangentials[i+nfaces].begin(),
			 data.unit_tangentials[i+nfaces].end(),
			 tangential);
	    }
	}
    }
}

#endif

  

template <int dim>
typename Mapping<dim>::InternalDataBase *
MappingQ1<dim>::get_face_data (const UpdateFlags        update_flags,
			       const Quadrature<dim-1> &quadrature) const
{
  InternalData* data = new InternalData(n_shape_functions);
  compute_face_data (update_flags,
		     QProjector<dim>::project_to_all_faces(quadrature),
		     quadrature.n_quadrature_points,
		     *data);

  return data;
}



template <int dim>
typename Mapping<dim>::InternalDataBase *
MappingQ1<dim>::get_subface_data (const UpdateFlags update_flags,
				  const Quadrature<dim-1>& quadrature) const
{
  InternalData* data = new InternalData(n_shape_functions);
  compute_face_data (update_flags,
		     QProjector<dim>::project_to_all_subfaces(quadrature),
		     quadrature.n_quadrature_points,
		     *data);

  return data;
}




template <int dim>
void
MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cell,
			      const unsigned int       n_q_points,
			      const DataSetDescriptor  data_set,
			      InternalData            &data,
			      std::vector<Point<dim> > &quadrature_points) const
{
  const UpdateFlags update_flags(data.current_update_flags());

  if (update_flags & update_q_points)
    {
      Assert (quadrature_points.size() == n_q_points,
	      ExcDimensionMismatch(quadrature_points.size(), n_q_points));
      std::fill(quadrature_points.begin(),
		quadrature_points.end(),
		Point<dim>());
    }

  if (update_flags & update_covariant_transformation)
    {
      Assert (data.covariant.size() == n_q_points,
	      ExcDimensionMismatch(data.covariant.size(), n_q_points));
    }
  
  if (update_flags & update_contravariant_transformation)
    {
      Assert (data.contravariant.size() == n_q_points,
	      ExcDimensionMismatch(data.contravariant.size(), n_q_points));
      std::fill(data.contravariant.begin(),
		data.contravariant.end(),
		Tensor<2,dim>());
    }
  
  if (update_flags & update_jacobian_grads)
    {
      Assert(false, ExcNotImplemented());
//      Assert (covariant_grads.size () == n_q_points,
//	      ExcDimensionMismatch(covariant_grads.size(), n_q_points));
    }
  
				   // if necessary, recompute the
				   // support points of the
				   // transformation of this cell
				   // (note that we need to first
				   // check the triangulation pointer,
				   // since otherwise the second test
				   // might trigger an exception if
				   // the triangulations are not the
				   // same)
  if ((data.mapping_support_points.size() == 0)
      ||
      (&cell->get_triangulation() !=
       &data.cell_of_current_support_points->get_triangulation())
      ||
      (cell != data.cell_of_current_support_points))
    {
      compute_mapping_support_points(cell, data.mapping_support_points);
      data.cell_of_current_support_points = cell;
    }

                                   // first compute quadrature points
  if (update_flags & update_q_points)
    for (unsigned int point=0; point<n_q_points; ++point)
      for (unsigned int k=0; k<data.n_shape_functions; ++k)
        quadrature_points[point]
          += data.shape(point+data_set,k)
             * data.mapping_support_points[k];
  
                                   // then Jacobians
  if (update_flags & update_contravariant_transformation)
    for (unsigned int point=0; point<n_q_points; ++point)
      for (unsigned int k=0; k<data.n_shape_functions; ++k)
        for (unsigned int i=0; i<dim; ++i)
          for (unsigned int j=0; j<dim; ++j)
            data.contravariant[point][i][j]
              += (data.derivative(point+data_set, k)[j]
                  *
                  data.mapping_support_points[k][i]);

				       // invert contravariant for
				       // covariant transformation
				       // matrices
  if (update_flags & update_covariant_transformation)
    for (unsigned int point=0; point<n_q_points; ++point)
      data.covariant[point] = invert(data.contravariant[point]);
}


template <int dim>
void
MappingQ1<dim>::compute_mapping_support_points(
  const typename Triangulation<dim>::cell_iterator &cell,
  std::vector<Point<dim> > &a) const
{
  a.resize(GeometryInfo<dim>::vertices_per_cell);

  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
    a[i] = cell->vertex(vertex_mapping[i]);
}



template <int dim>
void
MappingQ1<dim>::fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
				const Quadrature<dim>                &q,
				typename Mapping<dim>::InternalDataBase      &mapping_data,
				std::vector<Point<dim> >                  &quadrature_points,
				std::vector<double>                       &JxW_values) const
{
  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
  Assert(data_ptr!=0, ExcInternalError());
  InternalData &data=*data_ptr;

  const unsigned int n_q_points=q.n_quadrature_points;
  
  compute_fill (cell, n_q_points, DataSetDescriptor::cell (),
                data, quadrature_points);

  
  const UpdateFlags update_flags(data.current_update_flags());
  const std::vector<double> &weights=q.get_weights();

				   // Multiply quadrature weights
				   // by Jacobian determinants
  if (update_flags & update_JxW_values)
    {      
      Assert (JxW_values.size() == n_q_points,
	      ExcDimensionMismatch(JxW_values.size(), n_q_points));
      for (unsigned int point=0; point<n_q_points; ++point)
	JxW_values[point]
	  = determinant(data.contravariant[point])*weights[point];
    }
}



template <int dim>
void
MappingQ1<dim>::compute_fill_face (const typename DoFHandler<dim>::cell_iterator &cell,
				   const unsigned int      face_no,
				   const bool              is_subface,
				   const unsigned int      n_q_points,
				   const DataSetDescriptor data_set,
				   const std::vector<double>   &weights,
				   InternalData           &data,
				   std::vector<Point<dim> >    &quadrature_points,
				   std::vector<double>         &JxW_values,
				   std::vector<Tensor<1,dim> > &boundary_forms,
				   std::vector<Point<dim> >    &normal_vectors) const
{
  compute_fill (cell, n_q_points, data_set, data, quadrature_points);

  const UpdateFlags update_flags(data.current_update_flags());
  
  if (update_flags & update_boundary_forms)
    {
      Assert (boundary_forms.size()==n_q_points,
	      ExcDimensionMismatch(boundary_forms.size(), n_q_points));
      if (update_flags & update_normal_vectors)
	Assert (normal_vectors.size()==n_q_points,
		ExcDimensionMismatch(normal_vectors.size(), n_q_points));
      if (update_flags & update_JxW_values)
	Assert (JxW_values.size() == n_q_points,
		ExcDimensionMismatch(JxW_values.size(), n_q_points));
      
      
      Assert (data.aux[0].size() <= data.unit_tangentials[face_no].size(),
	      ExcInternalError());
      transform_contravariant(data.aux[0].begin(),
			      data.aux[0].end(),
			      data.unit_tangentials[face_no].begin(),
			      data);

      typename std::vector<Tensor<1,dim> >::iterator
	result = boundary_forms.begin();
      const typename std::vector<Tensor<1,dim> >::iterator
	end = boundary_forms.end();
      const Tensor<1,dim> *
	tang1 = data.aux[0].begin();
      
      switch (dim)
	{
	  case 2:
          {
            for (; result != end; ++result, ++tang1)
              cross_product (*result, *tang1);
            break;
          };

	  case 3:
          {
            Assert (face_no+GeometryInfo<dim>::faces_per_cell <
                    data.unit_tangentials.n_rows(),
                    ExcInternalError());
            Assert (data.aux[1].size() <=
                    data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell].size(),
                    ExcInternalError());
            transform_contravariant(data.aux[1].begin(),
                                    data.aux[1].end(),
                                    data.unit_tangentials[
                                      face_no+GeometryInfo<dim>::faces_per_cell].begin(),
                                    data);
            const Tensor<1,dim> *tang2 = data.aux[1].begin();
            for (;result != end; ++result, ++tang1, ++tang2)
              cross_product (*result, *tang1, *tang2);
            break;
          };

	  default:
                Assert(false, ExcNotImplemented());
	}
      
      if (update_flags & (update_normal_vectors
			  | update_JxW_values))
	for (unsigned int i=0;i<boundary_forms.size();++i)
	  {
	    const double f = std::sqrt(contract(boundary_forms[i],
                                                boundary_forms[i]));
	    if (update_flags & update_JxW_values)
	      {
		JxW_values[i] = f * weights[i];
		if (is_subface)
		  JxW_values[i] /= GeometryInfo<dim>::subfaces_per_face;
	      }
	    if (update_flags & update_normal_vectors)
              normal_vectors[i] = boundary_forms[i] / f;
	  }
    }
}


template <int dim>
void
MappingQ1<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
				     const unsigned int       face_no,
				     const Quadrature<dim-1> &q,
				     typename Mapping<dim>::InternalDataBase &mapping_data,
				     std::vector<Point<dim> >     &quadrature_points,
				     std::vector<double>          &JxW_values,
				     std::vector<Tensor<1,dim> >  &boundary_forms,
				     std::vector<Point<dim> >     &normal_vectors) const
{
  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
  Assert(data_ptr!=0, ExcInternalError());
  InternalData &data=*data_ptr;

  const unsigned int n_q_points = q.n_quadrature_points;
  
  compute_fill_face (cell, face_no, false,
		     n_q_points,
		     DataSetDescriptor::face (face_no,
                                              cell->face_orientation(face_no),
                                              n_q_points),
		     q.get_weights(),
		     data,
		     quadrature_points,
		     JxW_values,
		     boundary_forms,
		     normal_vectors);
}


template <int dim>
void
MappingQ1<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
					const unsigned int       face_no,
					const unsigned int       sub_no,
					const Quadrature<dim-1> &q,
					typename Mapping<dim>::InternalDataBase &mapping_data,
					std::vector<Point<dim> >     &quadrature_points,
					std::vector<double>          &JxW_values,
					std::vector<Tensor<1,dim> >  &boundary_forms,
					std::vector<Point<dim> >     &normal_vectors) const
{
  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
  Assert(data_ptr!=0, ExcInternalError());
  InternalData &data=*data_ptr;

  const unsigned int n_q_points = q.n_quadrature_points;
  
  compute_fill_face (cell, face_no, true,
		     n_q_points,
		     DataSetDescriptor::sub_face (face_no, sub_no,
                                                  cell->face_orientation(face_no),
                                                  n_q_points),
		     q.get_weights(),
		     data,
		     quadrature_points,
		     JxW_values,
		     boundary_forms,
		     normal_vectors);
}


#if (deal_II_dimension == 1)

template <>
void
MappingQ1<1>::compute_fill_face (const DoFHandler<1>::cell_iterator &,
				 const unsigned int,
				 const bool,
				 const unsigned int,
				 const DataSetDescriptor,
				 const std::vector<double> &,
				 InternalData &,
				 std::vector<Point<1> > &,
				 std::vector<double> &,
				 std::vector<Tensor<1,1> > &,
				 std::vector<Point<1> > &) const
{
  Assert(false, ExcNotImplemented());
}


template <>
void
MappingQ1<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &,
				   const unsigned,
				   const Quadrature<0>&,
				   Mapping<1>::InternalDataBase&,
				   std::vector<Point<1> >&,
				   std::vector<double>&,
				   std::vector<Tensor<1,1> >&,
				   std::vector<Point<1> >&) const
{
  Assert(false, ExcNotImplemented());
}


template <>
void
MappingQ1<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &,
				      const unsigned,
				      const unsigned,
				      const Quadrature<0>&,
				      Mapping<1>::InternalDataBase&,
				      std::vector<Point<1> >&,
				      std::vector<double>&,
				      std::vector<Tensor<1,1> >&,
				      std::vector<Point<1> >&) const
{
  Assert(false, ExcNotImplemented());
}
#endif


template <int dim>
void
MappingQ1<dim>::transform_covariant (Tensor<1,dim>       *begin,
				     Tensor<1,dim>       *end,
				     const Tensor<1,dim> *src,
				     const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const InternalData *data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
  Assert(data_ptr!=0, ExcInternalError());
  const InternalData &data=*data_ptr;

  Assert (data.update_flags & update_covariant_transformation,
	  ExcAccessToUninitializedField());

//TODO: [GK] Can we do a similar assertion?  
//  Assert (dst.size() == data.contravariant.size(),
//	  ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size()));

  typename std::vector<Tensor<2,dim> >::const_iterator
    tensor = data.covariant.begin();
  
  while (begin!=end)
    {
      contract (*(begin++), *(src++), *(tensor++));
    }
}



template <int dim>
void
MappingQ1<dim>::transform_covariant (Tensor<2,dim>       *begin,
				     Tensor<2,dim>       *end,
				     const Tensor<2,dim> *src,
				     const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const InternalData *data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
  Assert(data_ptr!=0, ExcInternalError());
  const InternalData &data=*data_ptr;

  Assert (data.update_flags & update_covariant_transformation,
	  ExcAccessToUninitializedField());

//TODO: [GK] Can we do a similar assertion?  
//  Assert (dst.size() == data.contravariant.size(),
//	  ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size()));

  typename std::vector<Tensor<2,dim> >::const_iterator
    tensor = data.covariant.begin();
  
  while (begin!=end)
    {
      contract (*(begin++), *(src++), *(tensor++));
    }
}



template <int dim>
void
MappingQ1<dim>::
transform_contravariant (Tensor<1,dim>       *begin,
                         Tensor<1,dim>       *end,
                         const Tensor<1,dim> *src,
                         const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
  Assert(data_ptr!=0, ExcInternalError());
  const InternalData &data=*data_ptr;

  Assert (data.update_flags & update_contravariant_transformation,
	  ExcAccessToUninitializedField());

  typename std::vector<Tensor<2,dim> >::const_iterator
    tensor = data.contravariant.begin();
  
  while (begin!=end)
    contract (*(begin++), *(tensor++), *(src++));
}



template <int dim>
void
MappingQ1<dim>::
transform_contravariant (Tensor<2,dim>       *begin,
                         Tensor<2,dim>       *end,
                         const Tensor<2,dim> *src,
                         const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
  Assert(data_ptr!=0, ExcInternalError());
  const InternalData &data=*data_ptr;

  Assert (data.update_flags & update_contravariant_transformation,
	  ExcAccessToUninitializedField());

  typename std::vector<Tensor<2,dim> >::const_iterator
    tensor = data.contravariant.begin();
  
  while (begin!=end)
    contract (*(begin++), *(tensor++), *(src++));
}



template <int dim>
Point<dim>
MappingQ1<dim>::
transform_unit_to_real_cell (const typename Triangulation<dim>::cell_iterator &cell,
                             const Point<dim>                                 &p) const
{
				   // Use the get_data function to
				   // create an InternalData with data
				   // vectors of the right size and
				   // transformation shape values
				   // already computed at point p.
  const Quadrature<dim> point_quadrature(p);

  std::auto_ptr<InternalData>
    mdata (dynamic_cast<InternalData *> (
             get_data(update_transformation_values, point_quadrature)));

				   // compute the mapping support
				   // points
  compute_mapping_support_points(cell, mdata->mapping_support_points);

  return transform_unit_to_real_cell_internal(*mdata);
}



template <int dim>
Point<dim>
MappingQ1<dim>::
transform_unit_to_real_cell_internal (const InternalData &data) const
{
  const unsigned int n_mapping_points=data.mapping_support_points.size();
  Assert(data.shape_values.size()==n_mapping_points, ExcInternalError());
  
				   // use now the InternalData to
				   // compute the point in real space.
  Point<dim> p_real;
  for (unsigned int i=0; i<data.mapping_support_points.size(); ++i)
    p_real += data.mapping_support_points[i] * data.shape(0,i);

  return p_real;
}



template <int dim>
Point<dim>
MappingQ1<dim>::
transform_real_to_unit_cell (const typename Triangulation<dim>::cell_iterator &cell,
                             const Point<dim>                                 &p) const
{
				   // Let the start value of the
				   // newton iteration be the center
				   // of the unit cell
  Point<dim> p_unit;
  for (unsigned int i=0; i<dim; ++i)
    p_unit(i)=0.5;

				   // Use the get_data function to
				   // create an InternalData with data
				   // vectors of the right size and
				   // transformation shape values and
				   // derivatives already computed at
				   // point p_unit.
  const Quadrature<dim> point_quadrature(p_unit);
  std::auto_ptr<InternalData>
    mdata(dynamic_cast<InternalData *> (
            MappingQ1<dim>::get_data(update_transformation_values
                                     | update_transformation_gradients,
                                     point_quadrature)));

  MappingQ1<dim>::compute_mapping_support_points (cell,
                                                  mdata->mapping_support_points);
  Assert(mdata->mapping_support_points.size() ==
         GeometryInfo<dim>::vertices_per_cell,
	 ExcInternalError());

				   // perform the newton iteration.
  transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
  
  return p_unit;
}



template <int dim>
void
MappingQ1<dim>::
transform_real_to_unit_cell_internal
(const typename Triangulation<dim>::cell_iterator &cell,
 const Point<dim>                                 &p,
 InternalData                                     &mdata,
 Point<dim>                                       &p_unit) const
{
  const unsigned int n_shapes=mdata.shape_values.size();
  Assert(n_shapes!=0, ExcInternalError());
  Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
  
  std::vector<Point<dim> > &points=mdata.mapping_support_points;
  Assert(points.size()==n_shapes, ExcInternalError());
  
				   // Newton iteration to solve
				   // f(x)=p(x)-p=0
				   // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
  
				   // The start value is set to be the
				   // center of the unit cell.
  
				   // The shape values and derivatives
				   // of the mapping at this point are
				   // previously computed.

  
				   // f(x)
  Point<dim> p_real(transform_unit_to_real_cell_internal(mdata));
  Point<dim> f = p_real-p;

  const double eps=1e-15*cell->diameter();
  unsigned int loop=0;
  while (f.square()>eps*eps && loop++<10)
    {      
				       // f'(x)
      Tensor<2,dim> df;
      for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
	{
	  const Tensor<1,dim> &grad_transform=mdata.derivative(0,k);
	  const Point<dim> &point=points[k];
	  
	  for (unsigned int i=0; i<dim; ++i)
	    for (unsigned int j=0; j<dim; ++j)
	      df[i][j]+=point[i]*grad_transform[j];
	}
      
				       // Solve  [f'(x)]d=f(x)
      Tensor<1,dim> d;
      Tensor<2,dim> df_1;

      df_1 = invert(df);
      contract (d, df_1, static_cast<const Tensor<1,dim>&>(f));

				       // update of p_unit
      p_unit -= d;
				       // shape values and derivatives
				       // at new p_unit point
      compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
      
				       // f(x)
      p_real = transform_unit_to_real_cell_internal(mdata);
      f = p_real-p;
    }
}

//----------------------------------------------------------------------//

template class MappingQ1<deal_II_dimension>;

invivomapping_q1__ZNK9MappingQ1ILi3EE12compute_fillERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjN10QProjectorILi3EE17DataSetDescriptorERNS0_12InternalDataERSt6vectorI5PointILi3EESaISE_EE_511

Real time: 6.839042e+09 -- Predicted time: 8.004070e+09

Callcount: 200996

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
44406	### 1	51860.3272562573	2.132800e+04	1.752400e+04	17.8
57682	### 2	85055.8160029189	8.199000e+03	5.932000e+03	27.6
85151	### 3	40484.6544719299	2.358000e+03	6.640000e+02	71.8
38981	### 4	15112.1653793958	2.985340e+05	2.627720e+05	11.9
110385	### 5	15905.4443647815	1.001960e+05	8.802000e+04	12.1

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q_include_fe_mapping_q__ZNK8MappingQILi3EE19transform_covariantEP6TensorILi1ELi3EES3_PKS2_RKN7MappingILi3EE16InternalDataBaseE_1159

Real time: 1.061856e+09 -- Predicted time: 2.817437e+09

Callcount: 3736156

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
2983302	### 1	1520812.69879824	6.950000e+02	1.640000e+02	76.4
447132	### 2	1815389.03416444	9.240000e+02	4.320000e+02	53.2
2953975	### 3	384503.332439062	2.160000e+02	3.200000e+01	85.1

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivofe_include_base_tensor_base__ZN6TensorILi1ELi3EEC2Eb_330

Real time: 7.223488e+08 -- Predicted time: 4.667876e+09

Callcount: 30017177

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5901255	### 1	30115327.1666667	1.550000e+02	2.400000e+01	84.5

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivocompressed_sparsity_pattern__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_iterator_base_funcsZSt9advanceISt23_Rb_tree_const_iteratorIjElEvRT_T0_St26bidirectional_iterator_tag_141

Real time: 4.813426e+09 -- Predicted time: 2.089787e+10

Callcount: 14910185

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
11332995	### 1	791270.895231783	2.015000e+03	7.360000e+02	63.4
10836508	### 2	4659182.78182895	1.422000e+03	2.640000e+02	81.4
9065929	### 3	2217825.33303805	8.450000e+02	7.600000e+01	91.0
3100098	### 4	1658068.29023802	1.482000e+03	3.880000e+02	73.8
2618622	### 5	2906766.95223727	9.140000e+02	1.720000e+02	81.1
3313228	### 6	2074487.91319344	1.906000e+03	5.040000e+02	73.5
6866998	### 7	629424.047609547	1.999000e+03	9.680000e+02	51.5
7654632	### 8	154030.852206614	3.102000e+03	1.660000e+03	46.4

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q_include_base_tensor__Z8contractILi3EEvR6TensorILi1EXT_EERKS1_RKS0_ILi2EXT_EE_502

Real time: 1.722230e+09 -- Predicted time: 1.163764e+10

Callcount: 63821468

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
12547041	### 1	32462663.1522979	1.740000e+02	2.800000e+01	83.9
45723359	### 2	31521785.1556524	1.900000e+02	2.400000e+01	87.3

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivostep_14__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_construct__ZNSt12_Destroy_auxILb0EE9__destroyIPSt3setIjSt4lessIjESaIjEEEEvT_S8__102

Real time: 3.677175e+08 -- Predicted time: 0.000000e+00

Callcount: 18

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
13	### 1	59.7087661500223	0.000000e+00	2.047056e+06	100.0
16	### 2	2.98472930493437	0.000000e+00	8.268019e+07	100.0

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivocompressed_sparsity_pattern__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_tree__ZNSt8_Rb_treeIjjSt9_IdentityIjESt4lessIjESaIjEE24_M_get_insert_unique_posERKj_1321

Real time: 1.393543e+09 -- Predicted time: 5.988079e+09

Callcount: 13395685

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
2633538	### 1	3764550.77860929	3.480000e+02	1.120000e+02	67.8
10894808	### 2	4799196.97803069	4.820000e+02	4.000000e+01	91.7
8203242	### 3	3829849.27776982	4.740000e+02	1.440000e+02	69.6
9048125	### 4	1645070.02782194	3.340000e+02	2.040000e+02	38.9

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivostep_14_include_base_tensor_baseZNK6TensorILi1ELi3EEmlERKS0493

Real time: 3.439522e+09 -- Predicted time: 3.286173e+10

Callcount: 143234298

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
28159280	### 1	143146545.582015	2.290000e+02	2.400000e+01	89.5
67347440	### 2	196540.389767921	4.130000e+02	3.600000e+01	91.2

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q_include_base_tensor_baseZN6TensorILi1ELi3EEC2ERKS0350

Real time: 6.715702e+09 -- Predicted time: 1.143546e+11

Callcount: 279315570

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
54912304	### 1	278441208.180053	4.100000e+02	2.400000e+01	94.1
271599465	### 2	1087963.41709714	1.780000e+02	2.800000e+01	84.2

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q_include_base_tensor__Z8contractILi3EEdRK6TensorILi1EXT_EES3__460

Real time: 5.469959e+08 -- Predicted time: 4.856730e+09

Callcount: 22781250

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
4478701	### 1	22801549.8333333	2.130000e+02	2.400000e+01	88.7

                //----------------------------  fe.cc  ---------------------------
//    $Id: fe.cc,v 1.2 2006/01/23 23:49:36 wolf Exp $
//    Version: $Name:  $
//
//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//----------------------------  fe.cc  ---------------------------


#include <fe/fe.h>
#include <base/memory_consumption.h>
#include <fe/mapping.h>
#include <fe/fe_values.h>
#include <base/quadrature.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <grid/tria_boundary.h>

#include <algorithm>
#include <functional>
#include <numeric>



/*------------------------------- FiniteElementBase ----------------------*/


template <int dim>
const double FiniteElementBase<dim>::fd_step_length = 1.0e-6;


template <int dim>
void
FiniteElementBase<dim>::
InternalDataBase::initialize_2nd (const FiniteElement<dim> *element,
				  const Mapping<dim>       &mapping,
				  const Quadrature<dim>    &quadrature)
{
				   // if we shall compute second
				   // derivatives, then we do so by
				   // finite differencing the
				   // gradients. that we do by
				   // evaluating the gradients of
				   // shape values at points shifted
				   // star-like a little in each
				   // coordinate direction around each
				   // quadrature point.
				   //
				   // therefore generate 2*dim (the
				   // number of evaluation points)
				   // FEValues objects with slightly
				   // shifted positions
  std::vector<Point<dim> > diff_points (quadrature.n_quadrature_points);
  
  differences.resize(2*dim);
  for (unsigned int d=0; d<dim; ++d)
    {
      Point<dim> shift;
      shift (d) = fd_step_length;

				       // generate points and FEValues
				       // objects shifted in
				       // plus-direction. note that
				       // they only need to compute
				       // gradients, not more
      for (unsigned int i=0; i<diff_points.size(); ++i)
	diff_points[i] = quadrature.point(i) + shift;
      const Quadrature<dim> plus_quad (diff_points);
      differences[d] = new FEValues<dim> (mapping, *element,
					  plus_quad, update_gradients);

				       // now same in minus-direction
      for (unsigned int i=0; i<diff_points.size(); ++i)
	diff_points[i] = quadrature.point(i) - shift;
      const Quadrature<dim> minus_quad (diff_points);
      differences[d+dim] = new FEValues<dim> (mapping, *element,
					      minus_quad, update_gradients); 
    }
}




template <int dim>
FiniteElementBase<dim>::InternalDataBase::~InternalDataBase ()
{
  for (unsigned int i=0; i<differences.size (); ++i)
    if (differences[i] != 0)
      {
					 // delete pointer and set it
					 // to zero to avoid
					 // inadvertant use
	delete differences[i];
	differences[i] = 0;
      };
}




template <int dim>
FiniteElementBase<dim>::
FiniteElementBase (const FiniteElementData<dim> &fe_data,
                   const std::vector<bool> &restriction_is_additive_flags,
                   const std::vector<std::vector<bool> > &nonzero_components)
		:
		FiniteElementData<dim> (fe_data),
                system_to_component_table (this->dofs_per_cell),
                face_system_to_component_table(this->dofs_per_face),
                system_to_base_table(this->dofs_per_cell),
                face_system_to_base_table(this->dofs_per_face),		
                component_to_base_table (this->components,
                                         std::make_pair(0U, 0U)),
                restriction_is_additive_flags(restriction_is_additive_flags),
		nonzero_components (nonzero_components),
		n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)),
		cached_primitivity (std::find_if (n_nonzero_components_table.begin(),
                                                  n_nonzero_components_table.end(),
                                                  std::bind2nd(std::not_equal_to<unsigned int>(),
                                                               1U))
                                    ==
                                    n_nonzero_components_table.end())
{
  Assert (restriction_is_additive_flags.size() == this->dofs_per_cell,
	  ExcDimensionMismatch(restriction_is_additive_flags.size(),
			       this->dofs_per_cell));
  Assert (nonzero_components.size() == this->dofs_per_cell,
	  ExcInternalError());
  for (unsigned int i=0; i<nonzero_components.size(); ++i)
    {
      Assert (nonzero_components[i].size() == this->n_components(),
	      ExcInternalError());
      Assert (std::count (nonzero_components[i].begin(),
			  nonzero_components[i].end(),
			  true)
	      >= 1,
	      ExcInternalError());
      Assert (n_nonzero_components_table[i] >= 1,
	      ExcInternalError());
      Assert (n_nonzero_components_table[i] <= this->n_components(),
	      ExcInternalError());      
    };
  
                                   // initialize some tables in the
				   // default way, i.e. if there is
				   // only one (vector-)component; if
				   // there are several, then the
				   // constructor of the derived class
				   // needs to overwrite these arrays
  for (unsigned int j=0 ; j<this->dofs_per_cell ; ++j)
    {
      system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
      system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);      
    }
  for (unsigned int j=0 ; j<this->dofs_per_face ; ++j)
    {
      face_system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
      face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);      
    }
}


template <int dim>
double
FiniteElementBase<dim>::shape_value (const unsigned int,
				     const Point<dim> &) const
{
  AssertThrow(false, ExcUnitShapeValuesDoNotExist());
  return 0.;
}



template <int dim>
double
FiniteElementBase<dim>::shape_value_component (const unsigned int,
					       const Point<dim> &,
					       const unsigned int) const
{
  AssertThrow(false, ExcUnitShapeValuesDoNotExist());
  return 0.;
}



template <int dim>
Tensor<1,dim>
FiniteElementBase<dim>::shape_grad (const unsigned int,
				    const Point<dim> &) const
{
  AssertThrow(false, ExcUnitShapeValuesDoNotExist());
  return Tensor<1,dim> ();
}



template <int dim>
Tensor<1,dim>
FiniteElementBase<dim>::shape_grad_component (const unsigned int,
					      const Point<dim> &,
					      const unsigned int) const
{
  AssertThrow(false, ExcUnitShapeValuesDoNotExist());
  return Tensor<1,dim> ();
}



template <int dim>
Tensor<2,dim>
FiniteElementBase<dim>::shape_grad_grad (const unsigned int,
					 const Point<dim> &) const
{
  AssertThrow(false, ExcUnitShapeValuesDoNotExist());
  return Tensor<2,dim> ();
}



template <int dim>
Tensor<2,dim>
FiniteElementBase<dim>::shape_grad_grad_component (const unsigned int,
						   const Point<dim> &,
						   const unsigned int) const
{
  AssertThrow(false, ExcUnitShapeValuesDoNotExist());
  return Tensor<2,dim> ();
}


template <int dim>
const FullMatrix<double> &
FiniteElementBase<dim>::get_restriction_matrix (const unsigned int child) const
{
  Assert (child<GeometryInfo<dim>::children_per_cell,
	  ExcIndexRange(child, 0, GeometryInfo<dim>::children_per_cell));
  Assert (restriction[child].n() != 0, ExcProjectionVoid());
  return restriction[child];
}



template <int dim>
const FullMatrix<double> &
FiniteElementBase<dim>::get_prolongation_matrix (const unsigned int child) const
{
  Assert (child<GeometryInfo<dim>::children_per_cell,
	  ExcIndexRange(child, 0, GeometryInfo<dim>::children_per_cell));
  Assert (prolongation[child].n() != 0, ExcEmbeddingVoid());
  return prolongation[child];
}


template <int dim>
bool
FiniteElementBase<dim>::prolongation_is_implemented () const
{
  for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
    {
      Assert ((prolongation[c].m() == this->dofs_per_cell) ||
              (prolongation[c].m() == 0),
              ExcInternalError());
      Assert ((prolongation[c].n() == this->dofs_per_cell) ||
              (prolongation[c].n() == 0),
              ExcInternalError());
      if ((prolongation[c].m() == 0) ||
          (prolongation[c].n() == 0))
        return false;
    }
  return true;
}



template <int dim>
bool
FiniteElementBase<dim>::constraints_are_implemented () const
{
  return (this->dofs_per_face  == 0) || (interface_constraints.m() != 0);
}



template <int dim>
const FullMatrix<double> &
FiniteElementBase<dim>::constraints () const
{
  Assert ((this->dofs_per_face  == 0) || (interface_constraints.m() != 0),
          ExcConstraintsVoid());
  
  if (dim==1)
    Assert ((interface_constraints.m()==0) && (interface_constraints.n()==0),
	    ExcWrongInterfaceMatrixSize(interface_constraints.m(),
					interface_constraints.n()));
  
  return interface_constraints;
}



template <int dim>
TableIndices<2>
FiniteElementBase<dim>::interface_constraints_size () const 
{
  switch (dim)
    {
      case 1:
            return TableIndices<2> (0U, 0U);
      case 2:
            return TableIndices<2> (this->dofs_per_vertex +
                                    2*this->dofs_per_line,
                                    this->dofs_per_face);
      case 3:
            return TableIndices<2> (5*this->dofs_per_vertex +
                                    12*this->dofs_per_line  +
                                    4*this->dofs_per_quad,
                                    this->dofs_per_face);
      default:
            Assert (false, ExcNotImplemented());
    };
  return TableIndices<2> (deal_II_numbers::invalid_unsigned_int,
                          deal_II_numbers::invalid_unsigned_int);
}



template <int dim>
void
FiniteElementBase<dim>::
get_interpolation_matrix (const FiniteElementBase<dim> &,
			  FullMatrix<double>           &) const
{
				   // by default, no interpolation
				   // implemented. so throw exception,
				   // as documentation says
  AssertThrow (false,
               typename FiniteElementBase<dim>::
               ExcInterpolationNotImplemented());
}

                                   


template <int dim>
bool
FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &f) const
{
  return ((static_cast<const FiniteElementData<dim>&>(*this) ==
	   static_cast<const FiniteElementData<dim>&>(f)) &&
	  (interface_constraints == f.interface_constraints));
}



template <int dim>
const std::vector<Point<dim> > &
FiniteElementBase<dim>::get_unit_support_points () const
{
				   // a finite element may define
				   // support points, but only if
				   // there are as many as there are
				   // degrees of freedom
  Assert ((unit_support_points.size() == 0) ||
	  (unit_support_points.size() == this->dofs_per_cell),
	  ExcInternalError());
  return unit_support_points;
}



template <int dim>
bool
FiniteElementBase<dim>::has_support_points () const
{
  return (unit_support_points.size() != 0);
}



template <int dim>
Point<dim>
FiniteElementBase<dim>::unit_support_point (const unsigned index) const
{
  Assert (index < this->dofs_per_cell,
          ExcIndexRange (index, 0, this->dofs_per_cell));
  Assert (unit_support_points.size() == this->dofs_per_cell,
          ExcFEHasNoSupportPoints ());
  return unit_support_points[index];
}



template <int dim>
const std::vector<Point<dim-1> > &
FiniteElementBase<dim>::get_unit_face_support_points () const
{
				   // a finite element may define
				   // support points, but only if
				   // there are as many as there are
				   // degrees of freedom on a face
  Assert ((unit_face_support_points.size() == 0) ||
	  (unit_face_support_points.size() == this->dofs_per_face),
	  ExcInternalError());
  return unit_face_support_points;
}



template <int dim>
bool
FiniteElementBase<dim>::has_face_support_points () const
{
  return (unit_face_support_points.size() != 0);
}



template <int dim>
Point<dim-1>
FiniteElementBase<dim>::unit_face_support_point (const unsigned index) const
{
  Assert (index < this->dofs_per_face,
          ExcIndexRange (index, 0, this->dofs_per_face));
  Assert (unit_face_support_points.size() == this->dofs_per_face,
          ExcFEHasNoSupportPoints ());
  return unit_face_support_points[index];
}



template <int dim>
unsigned int
FiniteElementBase<dim>::memory_consumption () const
{
  return (sizeof(FiniteElementData<dim>) +
	  MemoryConsumption::
	  memory_consumption<FullMatrix<double>, sizeof(restriction)/sizeof(restriction[0])>
	  (restriction)+
	  MemoryConsumption::memory_consumption
	  <FullMatrix<double>, sizeof(prolongation)/sizeof(prolongation[0])>
	  (prolongation) +
	  MemoryConsumption::memory_consumption (interface_constraints) +
	  MemoryConsumption::memory_consumption (system_to_component_table) +
	  MemoryConsumption::memory_consumption (face_system_to_component_table) +
	  MemoryConsumption::memory_consumption (system_to_base_table) +
	  MemoryConsumption::memory_consumption (face_system_to_base_table) +	  
	  MemoryConsumption::memory_consumption (component_to_base_table) +
	  MemoryConsumption::memory_consumption (restriction_is_additive_flags) +
	  MemoryConsumption::memory_consumption (nonzero_components) +
	  MemoryConsumption::memory_consumption (n_nonzero_components_table));
}



template <int dim>
void
FiniteElementBase<dim>::
compute_2nd (const Mapping<dim>                   &mapping,
	     const typename DoFHandler<dim>::cell_iterator &cell,
	     const unsigned int,
	     typename Mapping<dim>::InternalDataBase &mapping_internal,
	     InternalDataBase                     &fe_internal,
	     FEValuesData<dim>                    &data) const
{
//TODO[GK]: This function presently has a flaw: it ignores the offset parameter. This means that if it is called for faces or subfaces, it computes the data for _all_ faces or subfaces, even though we only need to have this for one of them. given that computing second derivatives is expensive, this is significant. What should be done is to honor the offset parameter and of course only compute what we actually need
  Assert ((fe_internal.update_each | fe_internal.update_once)
	  & update_second_derivatives,
	  ExcInternalError());

  const unsigned int total_nonzero_components
    = std::accumulate (n_nonzero_components_table.begin(),
                       n_nonzero_components_table.end(),
                       0U);
  
  Assert (data.shape_2nd_derivatives.n_rows() == total_nonzero_components,
	  ExcInternalError());
				   // Number of quadrature points
  const unsigned int n_q_points = data.shape_2nd_derivatives.n_cols();

				   // first reinit the fe_values
				   // objects used for the finite
				   // differencing stuff
  for (unsigned int d=0; d<dim; ++d)
    {
      fe_internal.differences[d]->reinit(cell);
      fe_internal.differences[d+dim]->reinit(cell);
    }

				   // collection of difference
				   // quotients of gradients in each
				   // direction (first index) and at
				   // all q-points (second index)
  Table<2,Tensor<1,dim> >    diff_quot (dim, n_q_points);
  std::vector<Tensor<1,dim> > diff_quot2 (n_q_points);

				   // for all nonzero components of
				   // all shape functions at all
				   // quadrature points and difference
				   // quotients in all directions:
  unsigned int total_index = 0;
  for (unsigned int shape_index=0; shape_index<this->dofs_per_cell; ++shape_index)
    for (unsigned int n=0; n<n_nonzero_components(shape_index); ++n, ++total_index)
      {
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int q=0; q<n_q_points; ++q)
            {
                                               // get gradient at points
                                               // shifted slightly to
                                               // the right and to the
                                               // left in the present
                                               // coordinate direction
                                               //
                                               // note that things
                                               // might be more
                                               // difficult if the
                                               // shape function has
                                               // more than one
                                               // non-zero component,
                                               // so find out about
                                               // the actual component
                                               // if necessary
              Tensor<1,dim> right, left;
              if (is_primitive(shape_index))
                {
                  right = fe_internal.differences[d1]->shape_grad(shape_index, q);
                  left  = fe_internal.differences[d1+dim]->shape_grad(shape_index, q);
                }
              else
                {
                                                   // get the
                                                   // component index
                                                   // of the n-th
                                                   // nonzero
                                                   // compoment
                  unsigned int component=0;
                  for (unsigned int nonzero_comp=0; component<this->n_components();
                       ++component)
                    if (nonzero_components[shape_index][component] == true)
                      {
                        ++nonzero_comp;
                                                         // check
                                                         // whether we
                                                         // have found
                                                         // the
                                                         // component
                                                         // we are
                                                         // looking
                                                         // for. note
                                                         // that
                                                         // nonzero_comp
                                                         // is 1-based
                                                         // by the way
                                                         // we compute
                                                         // it
                        if (nonzero_comp == n+1)
                          break;
                      }
                  Assert (component < this->n_components(),
                          ExcInternalError());

                  right = fe_internal.differences[d1]
                          ->shape_grad_component(shape_index, q, component);
                  left  = fe_internal.differences[d1+dim]
                          ->shape_grad_component(shape_index, q, component);
                };
              
                                               // compute the second
                                               // derivative from a
                                               // symmetric difference
                                               // approximation
              for (unsigned int d=0; d<dim; ++d)
                diff_quot[d][q][d1] = 1./(2*fd_step_length) * (right[d]-left[d]);
            }
        
                                         // up to now we still have
                                         // difference quotients on the
                                         // unit cell, so transform it
                                         // to something on the real
                                         // cell
        for (unsigned int d=0; d<dim; ++d)
          {
            Assert (diff_quot2.size() <=
                    diff_quot[d].size(),
                    ExcInternalError());
            mapping.transform_covariant (&*diff_quot2.begin(), &*diff_quot2.begin()+diff_quot2.size(),
                                         diff_quot[d].begin(),
                                         mapping_internal);
            
            for (unsigned int q=0; q<n_q_points; ++q)
              for (unsigned int d1=0; d1<dim; ++d1)
                data.shape_2nd_derivatives[total_index][q][d][d1] = diff_quot2[q][d1];
          }
      }
}



template <int dim>
std::vector<unsigned int>
FiniteElementBase<dim>::
compute_n_nonzero_components (const std::vector<std::vector<bool> > &nonzero_components)
{
  std::vector<unsigned int> retval (nonzero_components.size());
  for (unsigned int i=0; i<nonzero_components.size(); ++i)
    retval[i] = std::count (nonzero_components[i].begin(),
			    nonzero_components[i].end(),
			    true);
  return retval;
}



/*------------------------------- FiniteElement ----------------------*/

template <int dim>
FiniteElement<dim>::FiniteElement (const FiniteElementData<dim> &fe_data,
				   const std::vector<bool> &restriction_is_additive_flags,
				   const std::vector<std::vector<bool> > &nonzero_components)
                :
		FiniteElementBase<dim> (fe_data,
					restriction_is_additive_flags,
					nonzero_components)
{}



template <int dim>
FiniteElement<dim>::FiniteElement (const FiniteElement<dim> &)
                :
                FiniteElementBase<dim> (FiniteElementData<dim>(),
                                        std::vector<bool> (),
                                        std::vector<std::vector<bool> >())
{
  Assert (false,
          ExcMessage ("Finite element objects don't support copying "
                      "semantics through the copy constructor. If "
                      "you want to copy a finite element, use the "
                      "clone() function."));
}



template <int dim>
FiniteElement<dim>::~FiniteElement ()
{}



template <int dim>
typename Mapping<dim>::InternalDataBase *
FiniteElement<dim>::get_face_data (const UpdateFlags       flags,
				   const Mapping<dim>      &mapping,
				   const Quadrature<dim-1> &quadrature) const
{
  return get_data (flags, mapping,
		   QProjector<dim>::project_to_all_faces(quadrature));
}



template <int dim>
typename Mapping<dim>::InternalDataBase *
FiniteElement<dim>::get_subface_data (const UpdateFlags        flags,
				      const Mapping<dim>      &mapping,
				      const Quadrature<dim-1> &quadrature) const
{
  return get_data (flags, mapping,
		   QProjector<dim>::project_to_all_subfaces(quadrature));
}



template <int dim>
unsigned int
FiniteElement<dim>::memory_consumption () const
{
  return FiniteElementBase<dim>::memory_consumption ();
}




template <int dim>
const FiniteElement<dim>&
FiniteElement<dim>::base_element(const unsigned index) const
{
  Assert (index==0, ExcIndexRange(index,0,1));
  return *this;
}

/*------------------------------- Explicit Instantiations -------------*/

template class FiniteElementBase<deal_II_dimension>;
template class FiniteElement<deal_II_dimension>;

invivofe__ZNK17FiniteElementBaseILi3EE11compute_2ndERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjRNS2_16InternalDataBaseERNS0_16InternalDataBaseER12FEValuesDataILi3EE_489

Real time: 3.816822e+09 -- Predicted time: 0.000000e+00

Callcount: 6664

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
520	### 1	5174.71340529709	0.000000e+00	2.115880e+05	100.0
4115	### 2	1498.47457923449	0.000000e+00	1.807420e+06	100.0

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivofe_include_base_tensor_baseZN6TensorILi1ELi3EEaSERKS0398

Real time: 7.204908e+08 -- Predicted time: 5.554267e+09

Callcount: 30004893

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5898840	### 1	30023063.8333333	1.850000e+02	2.400000e+01	87.0

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivofe_values_include_base_tensor_baseZN6TensorILi1ELi3EEC2ERKS0350

Real time: 6.820756e+08 -- Predicted time: 1.105129e+10

Callcount: 28405766

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5579034	### 1	28409487	3.890000e+02	2.400000e+01	93.8

                //----------------------  tensor_product_polynomials.cc  ------------
//    $Id: tensor_product_polynomials.cc,v 1.1 2004/09/14 00:51:25 wolf Exp $
//    Version: $Name:  $
//
//    Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//----------------------  tensor_product_polynomials.cc  ------------


#include <base/tensor_product_polynomials.h>
#include <base/exceptions.h>
#include <base/table.h>



/* ------------------- TensorProductPolynomials -------------- */

template <>
void
TensorProductPolynomials<1>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[1]) const
{
  Assert (i<polynomials.size(), ExcInternalError());
  indices[0] = index_map[i];
}



template <>
void
TensorProductPolynomials<2>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[2]) const
{
  const unsigned int n_pols = polynomials.size();
  Assert (i<n_pols*n_pols, ExcInternalError());
  const unsigned int n=index_map[i];

  indices[0] = n % n_pols;
  indices[1] = n / n_pols;
}



template <>
void
TensorProductPolynomials<3>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[3]) const
{
  const unsigned int n_pols = polynomials.size();
  Assert (i<n_pols*n_pols*n_pols, ExcInternalError());
  const unsigned int n=index_map[i];

  indices[0] = n % n_pols;
  indices[1] = (n/n_pols) % n_pols;
  indices[2] = n / (n_pols*n_pols);
}


template <int dim>
void
TensorProductPolynomials<dim>::output_indices(std::ostream &out) const
{
  unsigned int ix[dim];
  for (unsigned int i=0; i<n_tensor_pols; ++i)
    {
      compute_index(i,ix);
      out << i << "\t";
      for (unsigned int d=0; d<dim; ++d)
	out << ix[d] << " ";
      out << std::endl;
    }
}



template <int dim>
void
TensorProductPolynomials<dim>::set_numbering(
  const std::vector<unsigned int> &renumber)
{
  Assert(renumber.size()==index_map.size(),
	 ExcDimensionMismatch(renumber.size(), index_map.size()));

  index_map=renumber;
  for (unsigned int i=0; i<index_map.size(); ++i)
    index_map_inverse[index_map[i]]=i;
}



template <int dim>
double
TensorProductPolynomials<dim>::compute_value (const unsigned int i,
                                              const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);
  
  double value=1.;
  for (unsigned int d=0; d<dim; ++d)
    value *= polynomials[indices[d]].value(p(d));
  
  return value;
}

  
template <int dim>
Tensor<1,dim>
TensorProductPolynomials<dim>::compute_grad (const unsigned int i,
                                             const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

                                   // compute values and
                                   // uni-directional derivatives at
                                   // the given point in each
                                   // co-ordinate direction
  std::vector<std::vector<double> > v(dim, std::vector<double> (2));
  for (unsigned int d=0; d<dim; ++d)
    polynomials[indices[d]].value(p(d), v[d]);
  
  Tensor<1,dim> grad;
  for (unsigned int d=0; d<dim; ++d)
    {
      grad[d] = 1.;
      for (unsigned int x=0; x<dim; ++x)
        grad[d] *= v[x][d==x];
    }
  
  return grad;
}


template <int dim>
Tensor<2,dim>
TensorProductPolynomials<dim>::compute_grad_grad (const unsigned int i,
                                                  const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

  std::vector<std::vector<double> > v(dim, std::vector<double> (3));
  for (unsigned int d=0; d<dim; ++d)
    polynomials[indices[d]].value(p(d), v[d]);
  
  Tensor<2,dim> grad_grad;
  for (unsigned int d1=0; d1<dim; ++d1)
    for (unsigned int d2=0; d2<dim; ++d2)
      {
        grad_grad[d1][d2] = 1.;
        for (unsigned int x=0; x<dim; ++x)
          {
            unsigned int derivative=0;
            if (d1==x || d2==x)
              {
                if (d1==d2)
                  derivative=2;
                else
                  derivative=1;
              } 
            grad_grad[d1][d2] *= v[x][derivative];
          }
      }

  return grad_grad;
}




template <int dim>
void
TensorProductPolynomials<dim>::
compute (const Point<dim>            &p,
         std::vector<double>         &values,
         std::vector<Tensor<1,dim> > &grads,
         std::vector<Tensor<2,dim> > &grad_grads) const
{
  Assert (values.size()==n_tensor_pols || values.size()==0,
          ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
  Assert (grads.size()==n_tensor_pols|| grads.size()==0,
          ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
  Assert (grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
          ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));

  const bool update_values     = (values.size() == n_tensor_pols),
             update_grads      = (grads.size()==n_tensor_pols),
             update_grad_grads = (grad_grads.size()==n_tensor_pols);

                                   // check how many
                                   // values/derivatives we have to
                                   // compute
  unsigned int n_values_and_derivatives = 0;
  if (update_values)
    n_values_and_derivatives = 1;
  if (update_grads)
    n_values_and_derivatives = 2;
  if (update_grad_grads)
    n_values_and_derivatives = 3;


                                   // compute the values (and
                                   // derivatives, if necessary) of
                                   // all polynomials at this
                                   // evaluation point
  Table<2,std::vector<double> > v(dim, polynomials.size());
  for (unsigned int d=0; d<dim; ++d)
    for (unsigned int i=0; i<polynomials.size(); ++i)
      {
        v(d,i).resize (n_values_and_derivatives, 0.);
        polynomials[i].value(p(d), v(d,i));
      };
  
  for (unsigned int i=0; i<n_tensor_pols; ++i)
    {
                                       // first get the
                                       // one-dimensional indices of
                                       // this particular tensor
                                       // product polynomial
      unsigned int indices[dim];
      compute_index (i, indices);
      
      if (update_values)
        {
          values[i] = 1;
          for (unsigned int x=0; x<dim; ++x)
            values[i] *= v(x,indices[x])[0];
        }
  
      if (update_grads)
        for (unsigned int d=0; d<dim; ++d)
          {
            grads[i][d] = 1.;            
            for (unsigned int x=0; x<dim; ++x)
              grads[i][d] *= v(x,indices[x])[d==x];
          }

      if (update_grad_grads)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            {
              grad_grads[i][d1][d2] = 1.;
              for (unsigned int x=0; x<dim; ++x)
                {
                  unsigned int derivative=0;
                  if (d1==x || d2==x)
                    {
                      if (d1==d2)
                        derivative=2;
                      else
                        derivative=1;
                    } 
                  grad_grads[i][d1][d2]
                    *= v(x,indices[x])[derivative];
                }
            }
    }
}



template<int dim>
unsigned int
TensorProductPolynomials<dim>::n() const
{
  return n_tensor_pols;
}




/* ------------------- AnisotropicPolynomials -------------- */


template <int dim>
AnisotropicPolynomials<dim>::
AnisotropicPolynomials(const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
		:
		polynomials (pols),
		n_tensor_pols(get_n_tensor_pols(pols))
{
  Assert (pols.size() == dim, ExcInvalidDim(pols.size()));
  for (unsigned int d=0; d<dim; ++d)
    Assert (pols[d].size() > 0,
            ExcMessage ("The number of polynomials must be larger than zero "
                        "for all coordinate directions."));
}




template <>
void
AnisotropicPolynomials<1>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[1]) const
{
  Assert (i<polynomials[0].size(), ExcInternalError());
  indices[0] = i;
}



template <>
void
AnisotropicPolynomials<2>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[2]) const
{
  Assert (i < polynomials[0].size()*polynomials[1].size(),
          ExcInternalError());

  indices[0] = i % polynomials[0].size();
  indices[1] = i / polynomials[0].size();
}



template <>
void
AnisotropicPolynomials<3>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[3]) const
{
  Assert (i < polynomials[0].size()*polynomials[1].size()*polynomials[2].size(),
          ExcInternalError());

  indices[0] = i % polynomials[0].size();
  indices[1] = (i/polynomials[0].size()) % polynomials[1].size();
  indices[2] = i / (polynomials[0].size()*polynomials[1].size());
}



template <int dim>
double
AnisotropicPolynomials<dim>::compute_value (const unsigned int i,
                                              const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);
  
  double value=1.;
  for (unsigned int d=0; d<dim; ++d)
    value *= polynomials[d][indices[d]].value(p(d));
  
  return value;
}

  
template <int dim>
Tensor<1,dim>
AnisotropicPolynomials<dim>::compute_grad (const unsigned int i,
                                             const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

                                   // compute values and
                                   // uni-directional derivatives at
                                   // the given point in each
                                   // co-ordinate direction
  std::vector<std::vector<double> > v(dim, std::vector<double> (2));
  for (unsigned int d=0; d<dim; ++d)
    polynomials[d][indices[d]].value(p(d), v[d]);
  
  Tensor<1,dim> grad;
  for (unsigned int d=0; d<dim; ++d)
    {
      grad[d] = 1.;
      for (unsigned int x=0; x<dim; ++x)
        grad[d] *= v[x][d==x];
    }
  
  return grad;
}


template <int dim>
Tensor<2,dim>
AnisotropicPolynomials<dim>::compute_grad_grad (const unsigned int i,
                                                  const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

  std::vector<std::vector<double> > v(dim, std::vector<double> (3));
  for (unsigned int d=0; d<dim; ++d)
    polynomials[d][indices[d]].value(p(d), v[d]);
  
  Tensor<2,dim> grad_grad;
  for (unsigned int d1=0; d1<dim; ++d1)
    for (unsigned int d2=0; d2<dim; ++d2)
      {
        grad_grad[d1][d2] = 1.;
        for (unsigned int x=0; x<dim; ++x)
          {
            unsigned int derivative=0;
            if (d1==x || d2==x)
              {
                if (d1==d2)
                  derivative=2;
                else
                  derivative=1;
              } 
            grad_grad[d1][d2] *= v[x][derivative];
          }
      }

  return grad_grad;
}




template <int dim>
void
AnisotropicPolynomials<dim>::
compute (const Point<dim>            &p,
         std::vector<double>         &values,
         std::vector<Tensor<1,dim> > &grads,
         std::vector<Tensor<2,dim> > &grad_grads) const
{
  Assert (values.size()==n_tensor_pols || values.size()==0,
          ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
  Assert (grads.size()==n_tensor_pols|| grads.size()==0,
          ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
  Assert (grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
          ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));

  const bool update_values     = (values.size() == n_tensor_pols),
             update_grads      = (grads.size()==n_tensor_pols),
             update_grad_grads = (grad_grads.size()==n_tensor_pols);

                                   // check how many
                                   // values/derivatives we have to
                                   // compute
  unsigned int n_values_and_derivatives = 0;
  if (update_values)
    n_values_and_derivatives = 1;
  if (update_grads)
    n_values_and_derivatives = 2;
  if (update_grad_grads)
    n_values_and_derivatives = 3;


                                   // compute the values (and
                                   // derivatives, if necessary) of
                                   // all polynomials at this
                                   // evaluation point
  std::vector<std::vector<std::vector<double> > > v(dim);
  for (unsigned int d=0; d<dim; ++d)
    {
      v[d].resize (polynomials[d].size());
      for (unsigned int i=0; i<polynomials[d].size(); ++i)
        {
          v[d][i].resize (n_values_and_derivatives, 0.);
          polynomials[d][i].value(p(d), v[d][i]);
        };
    }
  
  for (unsigned int i=0; i<n_tensor_pols; ++i)
    {
                                       // first get the
                                       // one-dimensional indices of
                                       // this particular tensor
                                       // product polynomial
      unsigned int indices[dim];
      compute_index (i, indices);
      
      if (update_values)
        {
          values[i] = 1;
          for (unsigned int x=0; x<dim; ++x)
            values[i] *= v[x][indices[x]][0];
        }
  
      if (update_grads)
        for (unsigned int d=0; d<dim; ++d)
          {
            grads[i][d] = 1.;            
            for (unsigned int x=0; x<dim; ++x)
              grads[i][d] *= v[x][indices[x]][d==x ? 1 : 0];
          }

      if (update_grad_grads)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            {
              grad_grads[i][d1][d2] = 1.;
              for (unsigned int x=0; x<dim; ++x)
                {
                  unsigned int derivative=0;
                  if (d1==x || d2==x)
                    {
                      if (d1==d2)
                        derivative=2;
                      else
                        derivative=1;
                    } 
                  grad_grads[i][d1][d2]
                    *= v[x][indices[x]][derivative];
                }
            }
    }
}



template<int dim>
unsigned int
AnisotropicPolynomials<dim>::n() const
{
  return n_tensor_pols;
}


template <int dim>
unsigned int
AnisotropicPolynomials<dim>::
get_n_tensor_pols (const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
{
  unsigned int y = 1;
  for (unsigned int d=0; d<dim; ++d)
    y *= pols[d].size();
  return y;
}



/* ------------------- explicit instantiations -------------- */
template class TensorProductPolynomials<1>;
template class TensorProductPolynomials<2>;
template class TensorProductPolynomials<3>;

template class AnisotropicPolynomials<1>;
template class AnisotropicPolynomials<2>;
template class AnisotropicPolynomials<3>;

invivotensor_product_polynomials__ZNK24TensorProductPolynomialsILi3EE13compute_valueEjRK5PointILi3EE_108

Real time: 1.509820e+09 -- Predicted time: 9.303152e+09

Callcount: 14347366

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
2820634	### 1	6601505.80206089	6.600000e+02	1.040000e+02	84.2
5403367	### 2	164143.481944912	7.090000e+02	1.200000e+02	83.0
11668817	### 3	6639089.34226384	6.300000e+02	1.000000e+02	84.1
2452473	### 4	263184.571769708	7.310000e+02	9.200000e+01	87.4
9649768	### 5	366080.895411981	6.460000e+02	9.600000e+01	85.1
4921876	### 6	323853.870032694	6.740000e+02	1.080000e+02	83.9

                //----------------------------------------------------------------
//    $Id: fe_poly.cc,v 1.1 2004/09/14 00:51:24 wolf Exp $
//    Version: $Name:  $
//
//    Copyright (C) 2004 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//----------------------------------------------------------------

#include <base/tensor_product_polynomials.h>
#include <base/polynomials_p.h>
#include <fe/fe_poly.h>
#include <fe/fe_values.h>


template <class POLY, int dim>
FE_Poly<POLY,dim>::FE_Poly (unsigned int degree,
			    const POLY& poly_space,
			    const FiniteElementData<dim> &fe_data,
			    const std::vector<bool> &restriction_is_additive_flags,
			    const std::vector<std::vector<bool> > &nonzero_components):
		FiniteElement<dim> (fe_data,
				    restriction_is_additive_flags,
				    nonzero_components),
                degree(degree),
                poly_space(poly_space)
{}


template <class POLY, int dim>
unsigned int
FE_Poly<POLY,dim>::get_degree () const
{
  return degree;
}


template <class POLY, int dim>
double
FE_Poly<POLY,dim>::shape_value (const unsigned int i,
				const Point<dim> &p) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  return poly_space.compute_value(i, p);
}


template <class POLY, int dim>
double
FE_Poly<POLY,dim>::shape_value_component (const unsigned int i,
					   const Point<dim> &p,
					   const unsigned int component) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  Assert (component == 0, ExcIndexRange (component, 0, 1));
  return poly_space.compute_value(i, p);
}



template <class POLY, int dim>
Tensor<1,dim>
FE_Poly<POLY,dim>::shape_grad (const unsigned int i,
				const Point<dim> &p) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  return poly_space.compute_grad(i, p);
}



template <class POLY, int dim>
Tensor<1,dim>
FE_Poly<POLY,dim>::shape_grad_component (const unsigned int i,
					  const Point<dim> &p,
					  const unsigned int component) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  Assert (component == 0, ExcIndexRange (component, 0, 1));
  return poly_space.compute_grad(i, p);
}



template <class POLY, int dim>
Tensor<2,dim>
FE_Poly<POLY,dim>::shape_grad_grad (const unsigned int i,
				     const Point<dim> &p) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  return poly_space.compute_grad_grad(i, p);
}



template <class POLY, int dim>
Tensor<2,dim>
FE_Poly<POLY,dim>::shape_grad_grad_component (const unsigned int i,
					       const Point<dim> &p,
					       const unsigned int component) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  Assert (component == 0, ExcIndexRange (component, 0, 1));
  return poly_space.compute_grad_grad(i, p);
}



//----------------------------------------------------------------------
// Auxiliary functions
//----------------------------------------------------------------------




template <class POLY, int dim>
UpdateFlags
FE_Poly<POLY,dim>::update_once (const UpdateFlags flags) const
{
				   // for this kind of elements, only
				   // the values can be precomputed
				   // once and for all. set this flag
				   // if the values are requested at
				   // all
  return (update_default | (flags & update_values));
}



template <class POLY, int dim>
UpdateFlags
FE_Poly<POLY,dim>::update_each (const UpdateFlags flags) const
{
  UpdateFlags out = update_default;

  if (flags & update_gradients)
    out |= update_gradients | update_covariant_transformation;
  if (flags & update_second_derivatives)
    out |= update_second_derivatives | update_covariant_transformation;

  return out;
}



//----------------------------------------------------------------------
// Data field initialization
//----------------------------------------------------------------------

template <class POLY, int dim>
typename Mapping<dim>::InternalDataBase *
FE_Poly<POLY,dim>::get_data (const UpdateFlags      update_flags,
			     const Mapping<dim>    &mapping,
			     const Quadrature<dim> &quadrature) const
{
				   // generate a new data object and
				   // initialize some fields
  InternalData* data = new InternalData;

				   // check what needs to be
				   // initialized only once and what
				   // on every cell/face/subface we
				   // visit
  data->update_once = update_once(update_flags);
  data->update_each = update_each(update_flags);
  data->update_flags = data->update_once | data->update_each;

  const UpdateFlags flags(data->update_flags);
  const unsigned int n_q_points = quadrature.n_quadrature_points;

				   // some scratch arrays
  std::vector<double> values(0);
  std::vector<Tensor<1,dim> > grads(0);
  std::vector<Tensor<2,dim> > grad_grads(0);

				   // initialize fields only if really
				   // necessary. otherwise, don't
				   // allocate memory
  if (flags & update_values)
    {
      values.resize (this->dofs_per_cell);
      data->shape_values.reinit (this->dofs_per_cell,
				 n_q_points);
    }

  if (flags & update_gradients)
    {
      grads.resize (this->dofs_per_cell);
      data->shape_gradients.reinit (this->dofs_per_cell,
				    n_q_points);
    }

				   // if second derivatives through
				   // finite differencing is required,
				   // then initialize some objects for
				   // that
  if (flags & update_second_derivatives)
    data->initialize_2nd (this, mapping, quadrature);

				   // next already fill those fields
				   // of which we have information by
				   // now. note that the shape
				   // gradients are only those on the
				   // unit cell, and need to be
				   // transformed when visiting an
				   // actual cell
  if (flags & (update_values | update_gradients))
    for (unsigned int i=0; i<n_q_points; ++i)
      {
	poly_space.compute(quadrature.point(i),
			   values, grads, grad_grads);
	
	if (flags & update_values)
	  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
	    data->shape_values[k][i] = values[k];
	
	if (flags & update_gradients)
	  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
	    data->shape_gradients[k][i] = grads[k];
      }
  return data;
}




//----------------------------------------------------------------------
// Fill data of FEValues
//----------------------------------------------------------------------

template <class POLY, int dim>
void
FE_Poly<POLY,dim>::fill_fe_values (const Mapping<dim>                   &mapping,
			   const typename DoFHandler<dim>::cell_iterator &cell,
			   const Quadrature<dim>                &quadrature,
			   typename Mapping<dim>::InternalDataBase &mapping_data,
			   typename Mapping<dim>::InternalDataBase &fedata,
			   FEValuesData<dim>                    &data) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
  
  const UpdateFlags flags(fe_data.current_update_flags());

  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
    {
      if (flags & update_values)
	for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
	  data.shape_values(k,i) = fe_data.shape_values[k][i];
      
      if (flags & update_gradients)
	{
	  Assert (data.shape_gradients[k].size() <=
		  fe_data.shape_gradients[k].size(),
		  ExcInternalError());	  
	  mapping.transform_covariant(data.shape_gradients[k].begin(),
				      data.shape_gradients[k].end(),
				      fe_data.shape_gradients[k].begin(),
				      mapping_data);
	};
    }

  const typename QProjector<dim>::DataSetDescriptor dsd;
  if (flags & update_second_derivatives)
    this->compute_2nd (mapping, cell, dsd.cell(),
		       mapping_data, fe_data, data);
}



template <class POLY, int dim>
void
FE_Poly<POLY,dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
				const typename DoFHandler<dim>::cell_iterator &cell,
				const unsigned int                    face,
				const Quadrature<dim-1>              &quadrature,
				typename Mapping<dim>::InternalDataBase       &mapping_data,
				typename Mapping<dim>::InternalDataBase       &fedata,
				FEValuesData<dim>                    &data) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);

				   // offset determines which data set
				   // to take (all data sets for all
				   // faces are stored contiguously)
  
  const typename QProjector<dim>::DataSetDescriptor dsd;
  const typename QProjector<dim>::DataSetDescriptor offset
    = dsd.face (face, cell->face_orientation(face),
		quadrature.n_quadrature_points);
  
  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);

  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
    {
      if (flags & update_values)
        for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
	  data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
      
      if (flags & update_gradients)
	{
	  Assert (data.shape_gradients[k].size() + offset <=
		  fe_data.shape_gradients[k].size(),
		  ExcInternalError());
	  mapping.transform_covariant(data.shape_gradients[k].begin(),
				      data.shape_gradients[k].end(),
				      fe_data.shape_gradients[k].begin()+offset,
				      mapping_data);
	};
    }

  if (flags & update_second_derivatives)
    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
}



template <class POLY, int dim>
void
FE_Poly<POLY,dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
				   const typename DoFHandler<dim>::cell_iterator &cell,
				   const unsigned int                    face,
				   const unsigned int                    subface,
				   const Quadrature<dim-1>              &quadrature,
				   typename Mapping<dim>::InternalDataBase       &mapping_data,
				   typename Mapping<dim>::InternalDataBase       &fedata,
				   FEValuesData<dim>                    &data) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);

				   // offset determines which data set
				   // to take (all data sets for all
				 // sub-faces are stored contiguously)

  const typename QProjector<dim>::DataSetDescriptor dsd;
  const typename QProjector<dim>::DataSetDescriptor offset
    = dsd.sub_face (face, subface, cell->face_orientation(face),
		    quadrature.n_quadrature_points);

  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);

  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
    {
      if (flags & update_values)
        for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
	  data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
      
      if (flags & update_gradients)
	{
	  Assert (data.shape_gradients[k].size() + offset <=
		  fe_data.shape_gradients[k].size(),
		  ExcInternalError());
	  mapping.transform_covariant(data.shape_gradients[k].begin(),
				      data.shape_gradients[k].end(),
				      fe_data.shape_gradients[k].begin()+offset,
				      mapping_data);
	};
    }
  
  if (flags & update_second_derivatives)
    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
}



template <class POLY, int dim>
unsigned int
FE_Poly<POLY,dim>::n_base_elements () const
{
  return 1;
}



template <class POLY, int dim>
const FiniteElement<dim> &
FE_Poly<POLY,dim>::base_element (const unsigned int index) const
{
  Assert (index==0, ExcIndexRange(index, 0, 1));
  return *this;
}



template <class POLY, int dim>
unsigned int
FE_Poly<POLY,dim>::element_multiplicity (const unsigned int index) const
{
  Assert (index==0, ExcIndexRange(index, 0, 1));
  return 1;
}



template class FE_Poly<TensorProductPolynomials<deal_II_dimension>, deal_II_dimension>;
template class FE_Poly<PolynomialSpace<deal_II_dimension>, deal_II_dimension>;
template class FE_Poly<PolynomialsP<deal_II_dimension>, deal_II_dimension>;

invivofe_poly__ZNK7FE_PolyI24TensorProductPolynomialsILi3EELi3EE14fill_fe_valuesERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEERK10QuadratureILi3EERNS4_16InternalDataBaseESI_R12FEValuesDataILi3EE_252

Real time: 7.792083e+08 -- Predicted time: 3.152913e+07

Callcount: 167658

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
18864	### 1	50191.4205494091	0.000000e+00	1.008000e+03	100.0
35671	### 2	47026.5991374805	0.000000e+00	1.481200e+04	100.0
136355	### 3	82971.4035043435	3.800000e+02	4.800000e+01	87.3
83336	### 4	6941.47286045352	0.000000e+00	4.116000e+03	100.0

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q_include_base_tensor_baseZN6TensorILi1ELi3EEpLERKS0448

Real time: 3.261883e+09 -- Predicted time: 1.046865e+10

Callcount: 135858852

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
26709298	### 1	135956499.333333	7.700000e+01	2.400000e+01	68.8

                //----------------------------  mapping_q.cc  ---------------------------
//    $Id: mapping_q.cc,v 1.3 2006/01/23 23:49:36 wolf Exp $
//    Version: $Name:  $
//
//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//----------------------------  mapping_q.cc  ---------------------------

#include <fe/mapping_q.h>
#include <fe/fe_q.h>
#include <base/polynomial.h>
#include <base/quadrature.h>
#include <base/quadrature_lib.h>
#include <base/memory_consumption.h>
#include <base/tensor_product_polynomials.h>
#include <lac/full_matrix.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary.h>
#include <dofs/dof_accessor.h>
#include <fe/fe_tools.h>

#include <numeric>
#include <memory>


template <int dim>
const bool MappingQ<dim>::use_mapping_q_on_all_cells;


template<int dim>
MappingQ<dim>::InternalData::InternalData (const unsigned int n_shape_functions)
		:
		MappingQ1<dim>::InternalData(n_shape_functions),
				use_mapping_q1_on_current_cell(false),
				mapping_q1_data(1 << dim)
{
  this->is_mapping_q1_data=false;
}



template<int dim>
unsigned int
MappingQ<dim>::InternalData::memory_consumption () const 
{
  return (MappingQ1<dim>::InternalData::memory_consumption () +
	  MemoryConsumption::memory_consumption (unit_normals) +
	  MemoryConsumption::memory_consumption (use_mapping_q1_on_current_cell) +
	  MemoryConsumption::memory_consumption (mapping_q1_data));
}



#if deal_II_dimension == 1

// in 1d, it is irrelevant which polynomial degree to use, since all
// cells are scaled linearly
template<>
MappingQ<1>::MappingQ (const unsigned int):
		degree(1),
		n_inner(0),
		n_outer(0),
		tensor_pols(0),
		n_shape_functions(2),
		renumber(0)
{}


template<>
MappingQ<1>::~MappingQ ()
{}

#endif



template<typename number>
number power(const number x, const unsigned int y)
{
				   // since the power to which x is
				   // raised is usually the space
				   // dimension, and since this is
				   // rarely larger than three, the
				   // following code is optimal and
				   // cannot be further optimized by
				   // grouping of operands to reduce
				   // the number of multiplications
				   // from O(x) to O(log x)
  number value=1;
  for (unsigned int i=0; i<y; ++i)
    value *= x;
  return value;
}


template<int dim>
MappingQ<dim>::MappingQ (const unsigned int p)
                :
		degree(p),
		n_inner(power(degree-1, dim)),
		n_outer((dim==2) ? 4+4*(degree-1)
			:8+12*(degree-1)+6*(degree-1)*(degree-1)),
		tensor_pols(0),
		n_shape_functions(power(degree+1,dim)),
		renumber(0)
{
				   // Construct the tensor product
				   // polynomials used as shape
				   // functions for the Qp mapping of
				   // cells at the boundary.
  std::vector<Polynomials::LagrangeEquidistant> v;
  for (unsigned int i=0; i<=degree; ++i)
    v.push_back(Polynomials::LagrangeEquidistant(degree,i));

  tensor_pols = new TensorProductPolynomials<dim> (v);
  Assert (n_shape_functions==tensor_pols->n(),
	  ExcInternalError());
  Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
  
				   // build the renumbering of the
				   // shape functions of the Qp
				   // mapping.
  renumber.resize(n_shape_functions,0);
  FETools::lexicographic_to_hierarchic_numbering (FE_Q<dim>(degree),
						  renumber);

				   // build laplace_on_quad_vector
  if (degree>1)
    {
      if (dim >= 2)
	set_laplace_on_quad_vector(laplace_on_quad_vector);
      if (dim >= 3)
	set_laplace_on_hex_vector(laplace_on_hex_vector);
    }
}


template<int dim>
MappingQ<dim>::~MappingQ ()
{
  delete tensor_pols;
}



#if deal_II_dimension == 1

template<>
void
MappingQ<1>::compute_shapes_virtual (const std::vector<Point<1> > &unit_points,
				     MappingQ1<1>::InternalData   &data) const
{
  MappingQ1<1>::compute_shapes_virtual(unit_points, data);
}

#endif



template<int dim>
void
MappingQ<dim>::compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
				       typename MappingQ1<dim>::InternalData &data) const
{
  const unsigned int n_points=unit_points.size();
  std::vector<double> values;
  std::vector<Tensor<1,dim> > grads;
  if (data.shape_values.size()!=0)
    {
      Assert(data.shape_values.size()==n_shape_functions*n_points,
	     ExcInternalError());
      values.resize(n_shape_functions);
    }
  if (data.shape_derivatives.size()!=0)
    {
      Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
	     ExcInternalError());
      grads.resize(n_shape_functions);
    }
  
				   // dummy variable of size 0
  std::vector<Tensor<2,dim> > grad2;

  
  if (data.shape_values.size()!=0 || data.shape_derivatives.size()!=0)
    for (unsigned int point=0; point<n_points; ++point)
      {
	tensor_pols->compute(unit_points[point], values, grads, grad2);
	
	if (data.shape_values.size()!=0)
	  for (unsigned int i=0; i<n_shape_functions; ++i)
	    data.shape(point,renumber[i]) = values[i];
	
	if (data.shape_derivatives.size()!=0)
	  for (unsigned int i=0; i<n_shape_functions; ++i)
	    data.derivative(point,renumber[i]) = grads[i];
      }
}



template <int dim>
typename Mapping<dim>::InternalDataBase *
MappingQ<dim>::get_data (const UpdateFlags update_flags,
			 const Quadrature<dim> &quadrature) const
{
  InternalData *data = new InternalData(n_shape_functions);
  this->compute_data (update_flags, quadrature,
                      quadrature.n_quadrature_points, *data);
  if (!use_mapping_q_on_all_cells)
    this->compute_data (update_flags, quadrature,
                        quadrature.n_quadrature_points, data->mapping_q1_data);
  return data;
}



template <int dim>
typename Mapping<dim>::InternalDataBase *
MappingQ<dim>::get_face_data (const UpdateFlags update_flags,
			      const Quadrature<dim-1>& quadrature) const
{
  InternalData *data = new InternalData(n_shape_functions);
  const Quadrature<dim> q (QProjector<dim>::project_to_all_faces(quadrature));
  this->compute_face_data (update_flags, q,
                           quadrature.n_quadrature_points, *data);
  if (!use_mapping_q_on_all_cells)
    this->compute_face_data (update_flags, q,
                             quadrature.n_quadrature_points,
                             data->mapping_q1_data);
  return data;
}



template <int dim>
typename Mapping<dim>::InternalDataBase *
MappingQ<dim>::get_subface_data (const UpdateFlags update_flags,
				 const Quadrature<dim-1>& quadrature) const
{
  InternalData *data = new InternalData(n_shape_functions);
  const Quadrature<dim> q (QProjector<dim>::project_to_all_subfaces(quadrature));
  this->compute_face_data (update_flags, q,
                           quadrature.n_quadrature_points, *data);
  if (!use_mapping_q_on_all_cells)
    this->compute_face_data (update_flags, q,
                             quadrature.n_quadrature_points,
                             data->mapping_q1_data);
  return data;
}



template <int dim>
void
MappingQ<dim>::fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
			       const Quadrature<dim>                &q,
			       typename Mapping<dim>::InternalDataBase       &mapping_data,
			       std::vector<Point<dim> >             &quadrature_points,
			       std::vector<double>                  &JxW_values) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &data = dynamic_cast<InternalData&> (mapping_data);

				   // check whether this cell needs
				   // the full mapping or can be
				   // treated by a reduced Q1 mapping,
				   // e.g. if the cell is in the
				   // interior of the domain
  data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
					  || cell->has_boundary_lines());

				   // depending on this result, use
				   // this or the other data object
				   // for the mapping
  typename MappingQ1<dim>::InternalData *p_data=0;
  if (data.use_mapping_q1_on_current_cell)
    p_data=&data.mapping_q1_data;
  else
    p_data=&data;
  
  MappingQ1<dim>::fill_fe_values(cell, q, *p_data,
				 quadrature_points, JxW_values);
}



template <int dim>
void
MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
				    const unsigned int       face_no,
				    const Quadrature<dim-1> &q,
				    typename Mapping<dim>::InternalDataBase &mapping_data,
				    std::vector<Point<dim> >     &quadrature_points,
				    std::vector<double>          &JxW_values,
				    std::vector<Tensor<1,dim> >  &exterior_forms,
				    std::vector<Point<dim> >     &normal_vectors) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &data = dynamic_cast<InternalData&> (mapping_data);
  
				   // check whether this cell needs
				   // the full mapping or can be
				   // treated by a reduced Q1 mapping,
				   // e.g. if the cell is entirely in
				   // the interior of the domain. note
				   // that it is not sufficient to ask
				   // whether the present _face_ is in
				   // the interior, as the mapping on
				   // the face depends on the mapping
				   // of the cell, which in turn
				   // depends on the fact whether
				   // _any_ of the faces of this cell
				   // is at the boundary, not only the
				   // present face
  data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
					|| cell->has_boundary_lines());

				   // depending on this result, use
				   // this or the other data object
				   // for the mapping
  typename MappingQ1<dim>::InternalData *p_data=0;
  if (data.use_mapping_q1_on_current_cell)
    p_data=&data.mapping_q1_data;
  else
    p_data=&data;

  const unsigned int n_q_points=q.n_quadrature_points;
  this->compute_fill_face (cell, face_no, false,
                           n_q_points,
                           QProjector<dim>::DataSetDescriptor::
                           face (face_no, cell->face_orientation(face_no),
                                 n_q_points),
                           q.get_weights(),
                           *p_data,
                           quadrature_points, JxW_values,
                           exterior_forms, normal_vectors);
}


template <int dim>
void
MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
				       const unsigned int       face_no,
				       const unsigned int       sub_no,
				       const Quadrature<dim-1> &q,
				       typename Mapping<dim>::InternalDataBase &mapping_data,
				       std::vector<Point<dim> >     &quadrature_points,
				       std::vector<double>          &JxW_values,
				       std::vector<Tensor<1,dim> >  &exterior_forms,
				       std::vector<Point<dim> >     &normal_vectors) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &data = dynamic_cast<InternalData&> (mapping_data);

				   // check whether this cell needs
				   // the full mapping or can be
				   // treated by a reduced Q1 mapping,
				   // e.g. if the cell is entirely in
				   // the interior of the domain. note
				   // that it is not sufficient to ask
				   // whether the present _face_ is in
				   // the interior, as the mapping on
				   // the face depends on the mapping
				   // of the cell, which in turn
				   // depends on the fact whether
				   // _any_ of the faces of this cell
				   // is at the boundary, not only the
				   // present face
  data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
					|| cell->has_boundary_lines());

				   // depending on this result, use
				   // this or the other data object
				   // for the mapping
  typename MappingQ1<dim>::InternalData *p_data=0;
  if (data.use_mapping_q1_on_current_cell)
    p_data=&data.mapping_q1_data;
  else
    p_data=&data;

  const unsigned int n_q_points=q.n_quadrature_points;
  this->compute_fill_face (cell, face_no, true,
                           n_q_points,
                           QProjector<dim>::DataSetDescriptor::
                           sub_face (face_no, sub_no,
                                     cell->face_orientation(face_no),
                                     n_q_points),
                           q.get_weights(),
                           *p_data,
                           quadrature_points, JxW_values,
                           exterior_forms, normal_vectors);
}


#if deal_II_dimension==1

template <>
void
MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const
{
  Assert(false, ExcInternalError());
}

#else

template <int dim>
void
MappingQ<dim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
{
  Assert(degree>1, ExcInternalError());
  const unsigned int n_inner_2d=(degree-1)*(degree-1);
  const unsigned int n_outer_2d=4+4*(degree-1);

				   // first check whether we have precomputed
				   // the values for some polynomial degree;
				   // the sizes of arrays is
				   // n_inner_2d*n_outer_2d
  double const *loqv_ptr=0;
  switch (degree)
    {
                                       // for degree==1, we shouldn't have to
                                       // compute any support points, since
                                       // all of them are on the vertices
      
      case 2:
      {
                                         // (checked these values against the
                                         // output of compute_laplace_vector
                                         // again, and found they're indeed
                                         // right -- just in case someone
                                         // wonders where they come from --
                                         // WB)
	static const double loqv2[1*8]
	  ={1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.};
	loqv_ptr=&loqv2[0];
        Assert (sizeof(loqv2)/sizeof(loqv2[0]) ==
                n_inner_2d * n_outer_2d,
                ExcInternalError());
	
	break;
      }

      case 3:
      {
                                         // (same as above)
	static const double loqv3[4*12]
	  ={80/1053., 1/81., 11/1053., 1/81., 25/117., 44/351.,
	    7/117., 16/351., 7/117., 16/351., 25/117., 44/351.,
	    1/81., 80/1053., 1/81., 11/1053., 44/351., 25/117.,
	    25/117., 44/351., 16/351., 7/117., 7/117., 16/351.,
	    1/81., 11/1053., 1/81., 80/1053., 7/117., 16/351.,
	    16/351., 7/117., 25/117., 44/351., 44/351., 25/117.,
	    11/1053., 1/81., 80/1053., 1/81., 16/351., 7/117.,
	    44/351., 25/117., 44/351., 25/117., 16/351., 7/117.};
        Assert (sizeof(loqv3)/sizeof(loqv3[0]) ==
                n_inner_2d * n_outer_2d,
                ExcInternalError());
        
	loqv_ptr=&loqv3[0];
	
	break;
      }


      case 4:
      {
	static const double loqv4[9*16]
	  ={0.074059218503115712, -0.0010757446289059922,
            0.0019142922390714627, -0.0010757446289060069,
            0.22312738654318912, 0.13468513060151868,
            0.03812914216116723, 0.029131600026332517,
            0.02200737428129396, 0.016008355644312244,
            0.029131600026332527, 0.022007374281293915,
            0.016008355644312224, 0.22312738654318917,
            0.13468513060151874, 0.038129142161167237,
            
            0.0066480315133420603, 0.0066480315133427186,
            0.0028734528616576258, 0.0028734528616574584,
            0.15277168188202372, 0.23481527607092728,
            0.15277168188202397, 0.079035726825841868,
            0.059692382812499903, 0.036198648174158153,
            0.024962693117977778, 0.040819489554071289,
            0.024962693117977889, 0.079035726825843783,
            0.059692382812500312, 0.036198648174158236,
            
            -0.0010757446289069229, 0.074059218503115892,
            -0.0010757446289058842, 0.0019142922390713395,
            0.038129142161167293, 0.13468513060151846,
            0.22312738654318981, 0.22312738654318776,
            0.13468513060151835, 0.038129142161167202,
            0.016008355644312171, 0.022007374281293943,
            0.029131600026332617, 0.029131600026335094,
            0.02200737428129395, 0.016008355644312293,
            
            0.0066480315133420733, 0.0028734528616574727,
            0.0028734528616576362, 0.0066480315133427255,
            0.079035726825843755, 0.059692382812500257,
            0.036198648174158243, 0.024962693117977792,
            0.04081948955407131, 0.024962693117977903,
            0.079035726825841868, 0.059692382812499799,
            0.036198648174158098, 0.15277168188202384,
            0.23481527607092734, 0.15277168188202397,
            
            0.011067708333333018, 0.011067708333333358,
            0.011067708333333736, 0.01106770833333337,
            0.067708333333334217, 0.1035156250000009,
            0.067708333333334356, 0.067708333333333759,
            0.10351562499999903, 0.067708333333333995,
            0.067708333333333703, 0.10351562499999885,
            0.067708333333333898, 0.067708333333334245,
            0.10351562500000108, 0.067708333333334397,
            
            0.0028734528616571847, 0.0066480315133423621,
            0.0066480315133430378, 0.0028734528616573343,
            0.036198648174158195, 0.059692382812500278,
            0.079035726825844074, 0.15277168188202336,
            0.23481527607092598, 0.15277168188202367,
            0.036198648174158042, 0.059692382812499861,
            0.079035726825842201, 0.024962693117977788,
            0.040819489554074009, 0.024962693117977879,
            
            -0.0010757446289069133, 0.0019142922390713401,
            -0.0010757446289058651, 0.07405921850311592,
            0.029131600026335087, 0.022007374281293911,
            0.016008355644312279, 0.016008355644312175,
            0.022007374281293957, 0.029131600026332641,
            0.22312738654318776, 0.13468513060151827,
            0.038129142161167182, 0.038129142161167286,
            0.13468513060151852, 0.22312738654318992,
            
            0.0028734528616571756, 0.0028734528616573209,
            0.0066480315133430369, 0.0066480315133423742,
            0.024962693117977757, 0.040819489554073919,
            0.024962693117977847, 0.036198648174158049,
            0.059692382812499924, 0.079035726825842215,
            0.15277168188202328, 0.23481527607092575,
            0.15277168188202356, 0.036198648174158167,
            0.059692382812500319, 0.079035726825844088,
            
            0.0019142922390712369, -0.0010757446289068027,
            0.07405921850311617, -0.0010757446289067784,
            0.016008355644312283, 0.022007374281293967,
            0.02913160002633523, 0.038129142161167258,
            0.13468513060151821, 0.22312738654318864,
            0.038129142161167265, 0.13468513060151813,
            0.22312738654318859, 0.016008355644312276,
            0.022007374281294009, 0.029131600026335244
          };
        
        Assert (sizeof(loqv4)/sizeof(loqv4[0]) ==
                n_inner_2d * n_outer_2d,
                ExcInternalError());
        
	loqv_ptr=&loqv4[0];
	
	break;
      }
      
				       // no other cases implemented,
				       // so simply fall through
      default:
            break;
    }
  
  if (loqv_ptr!=0)
    {
				       // precomputed. copy values to
				       // the loqvs array
      loqvs.reinit(n_inner_2d, n_outer_2d);
      for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
	for (unsigned int k=0; k<n_outer_2d; ++k)
	  loqvs[unit_point][k]=loqv_ptr[unit_point*n_outer_2d+k];
    }
  else
    {
				       // not precomputed, then do so now
      if (dim==2)
	compute_laplace_vector(loqvs);
      else
					 // computing the Laplace vector for
					 // faces is not supported in 3d at
					 // present. presumably, doing so
					 // would not be so hard: we would
					 // only have to call the function in
					 // 2d, i.e. the quad(=face) values in
					 // 3d are equal to the quad(=cell)
					 // values in 2d. however, that would
					 // require us to link in the 2d
					 // library, which is kind of awkward
					 // (note that compute_laplace_vector
					 // really makes use of a lot of 2d
					 // stuff, such as FEValues etc). an
					 // alternative would be to precompute
					 // the values of this array for a
					 // couple of higher mapping orders,
					 // pin down their values and insert
					 // them into the array above.
	Assert (false, ExcNotImplemented());
  }

				   // the sum of weights of the points
				   // at the outer rim should be
				   // one. check this
  for (unsigned int unit_point=0; unit_point<loqvs.n_rows(); ++unit_point)
    Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
				     loqvs[unit_point].end(),0.)-1)<1e-13,
	   ExcInternalError());
}

#endif


#if deal_II_dimension==3

template <>
void
MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
{
  Assert(degree>1, ExcInternalError());

				   // first check whether we have
				   // precomputed the values for some
				   // polynomial degree
  double const *lohv_ptr=0;
  if (degree==2)
    {
      static const double loqv2[26]
  	={1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128.,
	  7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192.,
	  7/192., 7/192., 7/192., 7/192.,
	  1/12., 1/12., 1/12., 1/12., 1/12., 1/12.};
      
      lohv_ptr=&loqv2[0];
    }
  
  if (lohv_ptr!=0)
    {
				       // precomputed. copy values to
				       // the lohvs array
      lohvs.reinit(n_inner, n_outer);
      for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
	for (unsigned int k=0; k<n_outer; ++k)
	  lohvs[unit_point][k]=lohv_ptr[unit_point*n_outer+k];
    }
  else
				     // not precomputed, then do so now
    compute_laplace_vector(lohvs);
    
				   // the sum of weights of the points
				   // at the outer rim should be
				   // one. check this
  for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
    Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
				     lohvs[unit_point].end(),0.) - 1)<1e-13,
	   ExcInternalError());
}

#endif


template <int dim>
void
MappingQ<dim>::set_laplace_on_hex_vector(Table<2,double> &) const
{
  Assert(false, ExcInternalError());
}




#if deal_II_dimension==1

template <>
void
MappingQ<1>::compute_laplace_vector(Table<2,double> &) const
{
  Assert(false, ExcInternalError());
}

#else


template <int dim>
void
MappingQ<dim>::compute_laplace_vector(Table<2,double> &lvs) const
{
  Assert(lvs.n_rows()==0, ExcInternalError());
  Assert(dim==2 || dim==3, ExcNotImplemented());

                                   // for degree==1, we shouldn't have to
                                   // compute any support points, since all of
                                   // them are on the vertices
  Assert(degree>1, ExcInternalError());

				   // compute the shape
				   // gradients at the quadrature
				   // points on the unit cell
  const QGauss<dim> quadrature(degree+1);
  const unsigned int n_q_points=quadrature.n_quadrature_points;
  
  InternalData quadrature_data(n_shape_functions);
  quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
  this->compute_shapes(quadrature.get_points(), quadrature_data);
  
				   // Compute the stiffness matrix of
				   // the inner dofs
  FullMatrix<double> S(n_inner);
  for (unsigned int point=0; point<n_q_points; ++point)
    for (unsigned int i=0; i<n_inner; ++i)
      for (unsigned int j=0; j<n_inner; ++j)
	S(i,j) += contract(quadrature_data.derivative(point, n_outer+i),
			   quadrature_data.derivative(point, n_outer+j))
		  * quadrature.weight(point);
  
				   // Compute the components of T to be the
				   // product of gradients of inner and
				   // outer shape functions.
  FullMatrix<double> T(n_inner, n_outer);
  for (unsigned int point=0; point<n_q_points; ++point)
    for (unsigned int i=0; i<n_inner; ++i)
      for (unsigned int k=0; k<n_outer; ++k)
	T(i,k) += contract(quadrature_data.derivative(point, n_outer+i),
			   quadrature_data.derivative(point, k))
		  *quadrature.weight(point);
  
  FullMatrix<double> S_1(n_inner);
  S_1.invert(S);
  
  FullMatrix<double> S_1_T(n_inner, n_outer);
  
				   // S:=S_1*T
  S_1.mmult(S_1_T,T);
  
				   // Resize and initialize the
				   // lvs
  lvs.reinit (n_inner, n_outer);
  for (unsigned int i=0; i<n_inner; ++i)
    for (unsigned int k=0; k<n_outer; ++k)
      lvs(i,k) = -S_1_T(i,k);
}

#endif



template <int dim>
void
MappingQ<dim>::apply_laplace_vector(const Table<2,double> &lvs,
				    std::vector<Point<dim> > &a) const
{
                                   // check whether the data we need
                                   // is really available. if you fail
                                   // here and if
                                   // lvs==laplace_on_quad_vector in
                                   // the calling function, then we
                                   // didn't compute the quad laplace
                                   // vector. this is mentioned in the
                                   // constructor of this class,
                                   // although I don't understand the
                                   // reason for not aborting there
                                   // any more [WB]
  Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree));
  
  const unsigned int n_inner_apply=lvs.n_rows();
  Assert(n_inner_apply==n_inner || n_inner_apply==(degree-1)*(degree-1),
	 ExcInternalError());
  const unsigned int n_outer_apply=lvs.n_cols();
  Assert(a.size()==n_outer_apply,
	 ExcDimensionMismatch(a.size(), n_outer_apply));

				   // compute each inner point as
				   // linear combination of the outer
				   // points. the weights are given by
				   // the lvs entries, the outer
				   // points are the first (existing)
				   // elements of a
  for (unsigned int unit_point=0; unit_point<n_inner_apply; ++unit_point)
    {
      Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
      Point<dim> p;
      for (unsigned int k=0; k<n_outer_apply; ++k)
	p+=lvs[unit_point][k]*a[k];

      a.push_back(p);
    }
}


template <int dim>
void
MappingQ<dim>::compute_mapping_support_points(
  const typename Triangulation<dim>::cell_iterator &cell,
  std::vector<Point<dim> > &a) const
{
				   // if this is a cell for which we
				   // want to compute the full
				   // mapping, then get them from the
				   // following function
  if (use_mapping_q_on_all_cells || cell->has_boundary_lines())
    compute_support_points_laplace(cell, a);
  else
				     // otherwise: use a Q1 mapping
				     // for which the mapping shape
				     // function support points are
				     // simply the vertices of the
				     // cell
    {
      a.resize(GeometryInfo<dim>::vertices_per_cell);
      
      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
	a[i] = cell->vertex(this->vertex_mapping[i]);
    }
}

  
template <int dim>
void
MappingQ<dim>::compute_support_points_laplace(const typename Triangulation<dim>::cell_iterator &cell,
					      std::vector<Point<dim> > &a) const
{
				   // in any case, we need the
				   // vertices first
  a.resize(GeometryInfo<dim>::vertices_per_cell);
  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
    a[i] = cell->vertex(i);
  
  if (degree>1)
    switch (dim)
      {
	case 2:
					   // in 2d, add the
					   // points on the four
					   // bounding lines to
					   // the exterior (outer)
					   // points
	  add_line_support_points (cell, a);
	  apply_laplace_vector (laplace_on_quad_vector,a);
	  break;

	case 3:
					   // in 3d also add the
					   // points located on
					   // the boundary faces
	  add_line_support_points (cell, a);
	  add_quad_support_points (cell, a);
	  apply_laplace_vector (laplace_on_hex_vector, a);
	  break;
	      
	default:
	  Assert(false, ExcNotImplemented());
	  break;
      };
}





#if deal_II_dimension==1

template <>
void
MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
				      std::vector<Point<1> > &) const
{
				   // there are no points on bounding
				   // lines which are to be added
  const unsigned int dim=1;
  Assert (dim > 1, ExcImpossibleInDim(dim));
}

#endif


template <int dim>
void
MappingQ<dim>::add_line_support_points (const typename Triangulation<dim>::cell_iterator &cell,
					std::vector<Point<dim> > &a) const
{
  static const StraightBoundary<dim> straight_boundary;
				   // if we only need the midpoint,
				   // then ask for it.
  if (degree==2)
    {
      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
	{
	  const typename Triangulation<dim>::line_iterator line = cell->line(line_no);
	  const Boundary<dim> * const boundary
	    = (line->at_boundary() ?
	       &line->get_triangulation().get_boundary(line->boundary_indicator()) :
	       &straight_boundary);
	  
	  a.push_back(boundary->get_new_point_on_line(line));
	};
    }
  else
				     // otherwise call the more
				     // complicated functions and ask
				     // for inner points from the
				     // boundary description
    {
      std::vector<Point<dim> > line_points (degree-1);
      
				       // loop over each of the lines,
				       // and if it is at the
				       // boundary, then first get the
				       // boundary description and
				       // second compute the points on
				       // it
      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
	{
	  const typename Triangulation<dim>::line_iterator line = cell->line(line_no);
	  
	  const Boundary<dim> * const boundary
	    = (line->at_boundary() ?
	       &line->get_triangulation().get_boundary(line->boundary_indicator()) :
	       &straight_boundary);
	  
	  boundary->get_intermediate_points_on_line (line, line_points);
	  a.insert (a.end(), line_points.begin(), line_points.end());
	}
    }
}




#if deal_II_dimension==3


template<>
void
MappingQ<3>::
add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
                        std::vector<Point<3> >                &a) const
{
  const unsigned int faces_per_cell    = GeometryInfo<3>::faces_per_cell,
		     vertices_per_face = GeometryInfo<3>::vertices_per_face,
		     lines_per_face    = GeometryInfo<3>::lines_per_face,
		     vertices_per_cell = GeometryInfo<3>::vertices_per_cell;

                                   // mapping of faces to vertex
                                   // indices for properly oriented
                                   // faces
  static const unsigned int face_vertex_to_cell_vertex1
    [faces_per_cell][vertices_per_face]={{0,1,2,3},
					 {4,5,6,7},
					 {0,1,5,4},
					 {1,5,6,2},
					 {3,2,6,7},
					 {0,4,7,3}};
                                   // same if the face is in opposite
                                   // orientation
  static const unsigned int face_vertex_to_cell_vertex2
    [faces_per_cell][vertices_per_face]={{0,3,2,1},
					 {4,7,6,5},
					 {0,4,5,1},
					 {1,2,6,5},
					 {3,7,6,2},
					 {0,3,7,4}};

                                   // indices of the lines that bound
                                   // a properly oriented face
  static const unsigned int face_line_to_cell_line1
    [faces_per_cell][lines_per_face]={{0,1,2,3},
				      {4,5,6,7},
				      {0,9,4,8},
				      {9,5,10,1},
				      {2,10,6,11},
				      {8,7,11,3}};
                                   // same for faces in wrong
                                   // orientation
  static const unsigned int face_line_to_cell_line2
    [faces_per_cell][lines_per_face]={{3,2,1,0},
				      {7,6,5,4},
				      {8,4,9,0},
				      {1,10,5,9},
				      {11,6,10,2},
				      {3,11,7,8}};

  static const StraightBoundary<3> straight_boundary;
				   // used if face quad at boundary or
				   // entirely in the interior of the
				   // domain
  std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
				   // used if only one line of face
				   // quad is at boundary
  std::vector<Point<3> > b(4*degree);
  
  
				   // loop over all faces and collect
				   // points on them
  for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
    {
      const Triangulation<3>::face_iterator face = cell->face(face_no);

                                       // select the correct mappings
                                       // for the present face
      const bool face_orientation = cell->face_orientation(face_no);

      const unsigned int *face_vertex_to_cell_vertex
        = (face_orientation ?
           face_vertex_to_cell_vertex1[face_no] :
           face_vertex_to_cell_vertex2[face_no]);
      const unsigned int *face_line_to_cell_line
        = (face_orientation ?
           face_line_to_cell_line1[face_no] :
           face_line_to_cell_line2[face_no]);

                                       // some sanity checks up front
      for (unsigned int i=0; i<vertices_per_face; ++i)
        Assert(face->vertex_index(i)==
               cell->vertex_index(face_vertex_to_cell_vertex[i]),
               ExcInternalError());
      
      for (unsigned int i=0; i<lines_per_face; ++i)
        Assert(face->line(i)==
               cell->line(face_line_to_cell_line[i]),
               ExcInternalError());
      
				       // if face at boundary, then
				       // ask boundary object to
				       // return intermediate points
				       // on it
      if (face->at_boundary())
	{
	  face->get_triangulation().get_boundary(face->boundary_indicator())
	    .get_intermediate_points_on_quad (face, quad_points);
	  a.insert (a.end(), quad_points.begin(), quad_points.end());
	}
      else
	{
					   // face is not at boundary,
					   // but maybe some of its
					   // lines are. count them
	  unsigned int lines_at_boundary=0;
	  for (unsigned int i=0; i<lines_per_face; ++i)
	    if (face->line(i)->at_boundary())
	      ++lines_at_boundary;
	  
	  Assert(lines_at_boundary<lines_per_face, ExcInternalError());

					   // if at least one of the
					   // lines bounding this quad
					   // is at the boundary, then
					   // collect points
					   // separately
	  if (lines_at_boundary>0)
	    {
					       // call of function
					       // apply_laplace_vector
					       // increases size of b
					       // about 1. There
					       // resize b for the
					       // case the mentioned
					       // function was already
					       // called.
	      b.resize(4*degree);
	      
					       // b is of size
					       // 4*degree, make sure
					       // that this is the
					       // right size
	      Assert(b.size()==vertices_per_face+lines_per_face*(degree-1),
		     ExcDimensionMismatch(b.size(),
                                          vertices_per_face+lines_per_face*(degree-1)));
	      
					       // sort the points into b
              for (unsigned int i=0; i<vertices_per_face; ++i)
                b[i]=a[face_vertex_to_cell_vertex[i]];
		      
              for (unsigned int i=0; i<lines_per_face; ++i)
                for (unsigned int j=0; j<degree-1; ++j)
                  b[vertices_per_face+i*(degree-1)+j]=
                    a[vertices_per_cell+face_line_to_cell_line[i]*(degree-1)+j];

					       // Now b includes the
					       // right order of
					       // support points on
					       // the quad to apply
					       // the laplace vector
	      apply_laplace_vector(laplace_on_quad_vector, b);
	      Assert(b.size()==4*degree+(degree-1)*(degree-1),
		     ExcDimensionMismatch(b.size(), 4*degree+(degree-1)*(degree-1)));
	      
	      for (unsigned int i=0; i<(degree-1)*(degree-1); ++i)
		a.push_back(b[4*degree+i]);
	    }
	  else
	    {
					       // face is entirely in
					       // the interior. get
					       // intermediate points
					       // from a straight
					       // boundary object
	      straight_boundary.get_intermediate_points_on_quad (face, quad_points);
	      a.insert (a.end(), quad_points.begin(), quad_points.end());
	    }
	}
    }
}

#endif


template<int dim>
void
MappingQ<dim>::
add_quad_support_points(const typename Triangulation<dim>::cell_iterator &,
                        std::vector<Point<dim> > &) const
{
  Assert (dim > 2, ExcImpossibleInDim(dim));
}



template <int dim>
void
MappingQ<dim>::
transform_covariant (Tensor<1,dim>       *begin,
                     Tensor<1,dim>       *end,
                     const Tensor<1,dim> *src,
                     const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const typename MappingQ1<dim>::InternalData *q1_data =
    dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
  Assert(q1_data!=0, ExcInternalError());
  
  typename std::vector<Tensor<2,dim> >::const_iterator tensor;

  if (q1_data->is_mapping_q1_data)
    tensor = q1_data->covariant.begin();
  else
    {
      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
      Assert(data!=0, ExcInternalError());

      if (data->use_mapping_q1_on_current_cell)
	tensor = data->mapping_q1_data.covariant.begin();
      else
	tensor = data->covariant.begin();    
    }

  while (begin!=end)
    contract (*(begin++), *(src++), *(tensor++));
}



template <int dim>
void
MappingQ<dim>::transform_covariant (Tensor<2,dim>       *begin,
				    Tensor<2,dim>       *end,
				    const Tensor<2,dim> *src,
				    const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const typename MappingQ1<dim>::InternalData *q1_data =
    dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
  Assert(q1_data!=0, ExcInternalError());
  
  typename std::vector<Tensor<2,dim> >::const_iterator tensor;

  if (q1_data->is_mapping_q1_data)
    tensor = q1_data->covariant.begin();
  else
    {
      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
      Assert(data!=0, ExcInternalError());

      if (data->use_mapping_q1_on_current_cell)
	tensor = data->mapping_q1_data.covariant.begin();
      else
	tensor = data->covariant.begin();
    }

  while (begin!=end)
    contract (*(begin++), *(src++), *(tensor++));
}



template <int dim>
void
MappingQ<dim>::
transform_contravariant (Tensor<1,dim>       *begin,
                         Tensor<1,dim>       *end,
                         const Tensor<1,dim> *src,
                         const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const typename MappingQ1<dim>::InternalData *q1_data =
    dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
  Assert(q1_data!=0, ExcInternalError());
  
  typename std::vector<Tensor<2,dim> >::const_iterator tensor;

  if (q1_data->is_mapping_q1_data)
    tensor = q1_data->contravariant.begin();
  else
    {
      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
      Assert(data!=0, ExcInternalError());

      if (data->use_mapping_q1_on_current_cell)
	tensor = data->mapping_q1_data.contravariant.begin();
      else
	tensor = data->contravariant.begin();    
    }
  
  while (begin!=end)
    contract (*(begin++), *(tensor++), *(src++));
}



template <int dim>
void
MappingQ<dim>::transform_contravariant (Tensor<2,dim>       *begin,
					Tensor<2,dim>       *end,
					const Tensor<2,dim> *src,
					const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
  const typename MappingQ1<dim>::InternalData *q1_data =
    dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
  Assert(q1_data!=0, ExcInternalError());
  
  typename std::vector<Tensor<2,dim> >::const_iterator tensor;

  if (q1_data->is_mapping_q1_data)
    tensor = q1_data->contravariant.begin();
  else
    {
      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
      Assert(data!=0, ExcInternalError());

      if (data->use_mapping_q1_on_current_cell)
	tensor = data->mapping_q1_data.contravariant.begin();
      else
	tensor = data->contravariant.begin();    
    }

  while (begin!=end)
    contract (*(begin++), *(tensor++), *(src++));
}



template <int dim>
Point<dim>
MappingQ<dim>::
transform_unit_to_real_cell (const typename Triangulation<dim>::cell_iterator &cell,
                             const Point<dim>                                 &p) const
{
				   // Use the get_data function to
				   // create an InternalData with data
				   // vectors of the right size and
				   // transformation shape values
				   // already computed at point p.
  const Quadrature<dim> point_quadrature(p);
  std::auto_ptr<InternalData>
    mdata (dynamic_cast<InternalData *> (
             get_data(update_transformation_values, point_quadrature)));
  
  mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
					    || cell->has_boundary_lines());

  typename MappingQ1<dim>::InternalData
    *p_data = (mdata->use_mapping_q1_on_current_cell ?
               &mdata->mapping_q1_data :
               &*mdata);

  compute_mapping_support_points(cell, p_data->mapping_support_points);
  
  return this->transform_unit_to_real_cell_internal(*p_data);
}



template <int dim>
Point<dim>
MappingQ<dim>::
transform_real_to_unit_cell (const typename Triangulation<dim>::cell_iterator &cell,
                             const Point<dim>                                 &p) const
{
				   // first a Newton iteration based
				   // on a Q1 mapping
  Point<dim> p_unit = MappingQ1<dim>::transform_real_to_unit_cell(cell, p);
  
                                   // then a Newton iteration based on
                                   // the full MappingQ if we need
                                   // this
  if (cell->has_boundary_lines() || use_mapping_q_on_all_cells)
    {
      const Quadrature<dim> point_quadrature(p_unit);
      std::auto_ptr<InternalData>
        mdata (dynamic_cast<InternalData *> (
                 get_data(update_transformation_values |
                          update_transformation_gradients,
                          point_quadrature)));
      
      mdata->use_mapping_q1_on_current_cell = false;

      std::vector<Point<dim> > &points = mdata->mapping_support_points;
      compute_mapping_support_points (cell, points);

      this->transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
    }
  
  return p_unit;
}



template <int dim>
unsigned int
MappingQ<dim>::get_degree() const
{
  return degree;
}

  
// explicit instantiation
template class MappingQ<deal_II_dimension>;

invivomapping_q__ZNK8MappingQILi3EE20apply_laplace_vectorERK5TableILi2EdERSt6vectorI5PointILi3EESaIS7_EE_792

Real time: 4.953406e+08 -- Predicted time: 5.913699e+08

Callcount: 210830

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
151044	### 1	187932.892258208	9.580000e+02	6.400000e+02	33.1
50040	### 2	42563.1369266758	9.664000e+03	8.800000e+03	8.9

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q1_include_base_tensor_baseZN6TensorILi1ELi3EEpLERKS0448

Real time: 6.682630e+08 -- Predicted time: 5.569616e+09

Callcount: 27825616

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5470403	### 1	27848078	2.000000e+02	2.400000e+01	88.0

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EE5clearEv_534

Real time: 1.557901e+09 -- Predicted time: 2.300684e+10

Callcount: 64872848

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
12753738	### 1	64721732.8249564	3.540000e+02	2.400000e+01	93.2
56154385	### 2	198644.435751622	4.800000e+02	2.800000e+01	94.1

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivofull_matrix_double_include_lac_full_matrix_templates__ZNK10FullMatrixIdE6TmmultIdEEvRS_IT_ERKS3_b_573

Real time: 2.545289e+09 -- Predicted time: 2.486162e+09

Callcount: 448

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	221.745526838966	3.140000e+03	2.012000e+03	35.9
264	### 2	433.306805254859	5.736042e+06	5.847540e+06	1.9

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470

Real time: 3.261866e+09 -- Predicted time: 4.744976e+10

Callcount: 135858852

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
26709298	### 1	135959206.5	3.490000e+02	2.400000e+01	93.1

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q1_include_base_tensor_baseZN6TensorILi1ELi3EEC2ERKS0350

Real time: 1.348959e+09 -- Predicted time: 1.529519e+10

Callcount: 56176922

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
11044155	### 1	56232332.5	2.720000e+02	2.400000e+01	91.1

                //----------------------------------------------------------------
//    $Id: fe_poly.cc,v 1.1 2004/09/14 00:51:24 wolf Exp $
//    Version: $Name:  $
//
//    Copyright (C) 2004 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//----------------------------------------------------------------

#include <base/tensor_product_polynomials.h>
#include <base/polynomials_p.h>
#include <fe/fe_poly.h>
#include <fe/fe_values.h>


template <class POLY, int dim>
FE_Poly<POLY,dim>::FE_Poly (unsigned int degree,
			    const POLY& poly_space,
			    const FiniteElementData<dim> &fe_data,
			    const std::vector<bool> &restriction_is_additive_flags,
			    const std::vector<std::vector<bool> > &nonzero_components):
		FiniteElement<dim> (fe_data,
				    restriction_is_additive_flags,
				    nonzero_components),
                degree(degree),
                poly_space(poly_space)
{}


template <class POLY, int dim>
unsigned int
FE_Poly<POLY,dim>::get_degree () const
{
  return degree;
}


template <class POLY, int dim>
double
FE_Poly<POLY,dim>::shape_value (const unsigned int i,
				const Point<dim> &p) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  return poly_space.compute_value(i, p);
}


template <class POLY, int dim>
double
FE_Poly<POLY,dim>::shape_value_component (const unsigned int i,
					   const Point<dim> &p,
					   const unsigned int component) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  Assert (component == 0, ExcIndexRange (component, 0, 1));
  return poly_space.compute_value(i, p);
}



template <class POLY, int dim>
Tensor<1,dim>
FE_Poly<POLY,dim>::shape_grad (const unsigned int i,
				const Point<dim> &p) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  return poly_space.compute_grad(i, p);
}



template <class POLY, int dim>
Tensor<1,dim>
FE_Poly<POLY,dim>::shape_grad_component (const unsigned int i,
					  const Point<dim> &p,
					  const unsigned int component) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  Assert (component == 0, ExcIndexRange (component, 0, 1));
  return poly_space.compute_grad(i, p);
}



template <class POLY, int dim>
Tensor<2,dim>
FE_Poly<POLY,dim>::shape_grad_grad (const unsigned int i,
				     const Point<dim> &p) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  return poly_space.compute_grad_grad(i, p);
}



template <class POLY, int dim>
Tensor<2,dim>
FE_Poly<POLY,dim>::shape_grad_grad_component (const unsigned int i,
					       const Point<dim> &p,
					       const unsigned int component) const
{
  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
  Assert (component == 0, ExcIndexRange (component, 0, 1));
  return poly_space.compute_grad_grad(i, p);
}



//----------------------------------------------------------------------
// Auxiliary functions
//----------------------------------------------------------------------




template <class POLY, int dim>
UpdateFlags
FE_Poly<POLY,dim>::update_once (const UpdateFlags flags) const
{
				   // for this kind of elements, only
				   // the values can be precomputed
				   // once and for all. set this flag
				   // if the values are requested at
				   // all
  return (update_default | (flags & update_values));
}



template <class POLY, int dim>
UpdateFlags
FE_Poly<POLY,dim>::update_each (const UpdateFlags flags) const
{
  UpdateFlags out = update_default;

  if (flags & update_gradients)
    out |= update_gradients | update_covariant_transformation;
  if (flags & update_second_derivatives)
    out |= update_second_derivatives | update_covariant_transformation;

  return out;
}



//----------------------------------------------------------------------
// Data field initialization
//----------------------------------------------------------------------

template <class POLY, int dim>
typename Mapping<dim>::InternalDataBase *
FE_Poly<POLY,dim>::get_data (const UpdateFlags      update_flags,
			     const Mapping<dim>    &mapping,
			     const Quadrature<dim> &quadrature) const
{
				   // generate a new data object and
				   // initialize some fields
  InternalData* data = new InternalData;

				   // check what needs to be
				   // initialized only once and what
				   // on every cell/face/subface we
				   // visit
  data->update_once = update_once(update_flags);
  data->update_each = update_each(update_flags);
  data->update_flags = data->update_once | data->update_each;

  const UpdateFlags flags(data->update_flags);
  const unsigned int n_q_points = quadrature.n_quadrature_points;

				   // some scratch arrays
  std::vector<double> values(0);
  std::vector<Tensor<1,dim> > grads(0);
  std::vector<Tensor<2,dim> > grad_grads(0);

				   // initialize fields only if really
				   // necessary. otherwise, don't
				   // allocate memory
  if (flags & update_values)
    {
      values.resize (this->dofs_per_cell);
      data->shape_values.reinit (this->dofs_per_cell,
				 n_q_points);
    }

  if (flags & update_gradients)
    {
      grads.resize (this->dofs_per_cell);
      data->shape_gradients.reinit (this->dofs_per_cell,
				    n_q_points);
    }

				   // if second derivatives through
				   // finite differencing is required,
				   // then initialize some objects for
				   // that
  if (flags & update_second_derivatives)
    data->initialize_2nd (this, mapping, quadrature);

				   // next already fill those fields
				   // of which we have information by
				   // now. note that the shape
				   // gradients are only those on the
				   // unit cell, and need to be
				   // transformed when visiting an
				   // actual cell
  if (flags & (update_values | update_gradients))
    for (unsigned int i=0; i<n_q_points; ++i)
      {
	poly_space.compute(quadrature.point(i),
			   values, grads, grad_grads);
	
	if (flags & update_values)
	  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
	    data->shape_values[k][i] = values[k];
	
	if (flags & update_gradients)
	  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
	    data->shape_gradients[k][i] = grads[k];
      }
  return data;
}




//----------------------------------------------------------------------
// Fill data of FEValues
//----------------------------------------------------------------------

template <class POLY, int dim>
void
FE_Poly<POLY,dim>::fill_fe_values (const Mapping<dim>                   &mapping,
			   const typename DoFHandler<dim>::cell_iterator &cell,
			   const Quadrature<dim>                &quadrature,
			   typename Mapping<dim>::InternalDataBase &mapping_data,
			   typename Mapping<dim>::InternalDataBase &fedata,
			   FEValuesData<dim>                    &data) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
  
  const UpdateFlags flags(fe_data.current_update_flags());

  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
    {
      if (flags & update_values)
	for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
	  data.shape_values(k,i) = fe_data.shape_values[k][i];
      
      if (flags & update_gradients)
	{
	  Assert (data.shape_gradients[k].size() <=
		  fe_data.shape_gradients[k].size(),
		  ExcInternalError());	  
	  mapping.transform_covariant(data.shape_gradients[k].begin(),
				      data.shape_gradients[k].end(),
				      fe_data.shape_gradients[k].begin(),
				      mapping_data);
	};
    }

  const typename QProjector<dim>::DataSetDescriptor dsd;
  if (flags & update_second_derivatives)
    this->compute_2nd (mapping, cell, dsd.cell(),
		       mapping_data, fe_data, data);
}



template <class POLY, int dim>
void
FE_Poly<POLY,dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
				const typename DoFHandler<dim>::cell_iterator &cell,
				const unsigned int                    face,
				const Quadrature<dim-1>              &quadrature,
				typename Mapping<dim>::InternalDataBase       &mapping_data,
				typename Mapping<dim>::InternalDataBase       &fedata,
				FEValuesData<dim>                    &data) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);

				   // offset determines which data set
				   // to take (all data sets for all
				   // faces are stored contiguously)
  
  const typename QProjector<dim>::DataSetDescriptor dsd;
  const typename QProjector<dim>::DataSetDescriptor offset
    = dsd.face (face, cell->face_orientation(face),
		quadrature.n_quadrature_points);
  
  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);

  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
    {
      if (flags & update_values)
        for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
	  data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
      
      if (flags & update_gradients)
	{
	  Assert (data.shape_gradients[k].size() + offset <=
		  fe_data.shape_gradients[k].size(),
		  ExcInternalError());
	  mapping.transform_covariant(data.shape_gradients[k].begin(),
				      data.shape_gradients[k].end(),
				      fe_data.shape_gradients[k].begin()+offset,
				      mapping_data);
	};
    }

  if (flags & update_second_derivatives)
    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
}



template <class POLY, int dim>
void
FE_Poly<POLY,dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
				   const typename DoFHandler<dim>::cell_iterator &cell,
				   const unsigned int                    face,
				   const unsigned int                    subface,
				   const Quadrature<dim-1>              &quadrature,
				   typename Mapping<dim>::InternalDataBase       &mapping_data,
				   typename Mapping<dim>::InternalDataBase       &fedata,
				   FEValuesData<dim>                    &data) const
{
				   // convert data object to internal
				   // data for this class. fails with
				   // an exception if that is not
				   // possible
  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);

				   // offset determines which data set
				   // to take (all data sets for all
				 // sub-faces are stored contiguously)

  const typename QProjector<dim>::DataSetDescriptor dsd;
  const typename QProjector<dim>::DataSetDescriptor offset
    = dsd.sub_face (face, subface, cell->face_orientation(face),
		    quadrature.n_quadrature_points);

  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);

  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
    {
      if (flags & update_values)
        for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
	  data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
      
      if (flags & update_gradients)
	{
	  Assert (data.shape_gradients[k].size() + offset <=
		  fe_data.shape_gradients[k].size(),
		  ExcInternalError());
	  mapping.transform_covariant(data.shape_gradients[k].begin(),
				      data.shape_gradients[k].end(),
				      fe_data.shape_gradients[k].begin()+offset,
				      mapping_data);
	};
    }
  
  if (flags & update_second_derivatives)
    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
}



template <class POLY, int dim>
unsigned int
FE_Poly<POLY,dim>::n_base_elements () const
{
  return 1;
}



template <class POLY, int dim>
const FiniteElement<dim> &
FE_Poly<POLY,dim>::base_element (const unsigned int index) const
{
  Assert (index==0, ExcIndexRange(index, 0, 1));
  return *this;
}



template <class POLY, int dim>
unsigned int
FE_Poly<POLY,dim>::element_multiplicity (const unsigned int index) const
{
  Assert (index==0, ExcIndexRange(index, 0, 1));
  return 1;
}



template class FE_Poly<TensorProductPolynomials<deal_II_dimension>, deal_II_dimension>;
template class FE_Poly<PolynomialSpace<deal_II_dimension>, deal_II_dimension>;
template class FE_Poly<PolynomialsP<deal_II_dimension>, deal_II_dimension>;

invivofe_poly__ZNK7FE_PolyI24TensorProductPolynomialsILi3EELi3EE19fill_fe_face_valuesERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjRK10QuadratureILi2EERNS4_16InternalDataBaseESI_R12FEValuesDataILi3EE_305

Real time: 4.585134e+08 -- Predicted time: 0.000000e+00

Callcount: 77250

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
37445	### 1	36970.5076976696	0.000000e+00	1.013200e+04	100.0
61667	### 2	38480.007770491	0.000000e+00	2.060000e+03	100.0

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivomapping_q1_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470

Real time: 6.681448e+08 -- Predicted time: 3.953498e+09

Callcount: 27825616

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5470403	### 1	27841534.6666667	1.420000e+02	2.400000e+01	83.0

447-dealII

Exec Time : 3.343264e+10 cycles

Matching : 30.6%

__invivo__fe_values_include_base_tensor_base__ZN6TensorILi1ELi3EEpLERKS0__448

Real time: 6.815835e+08 -- Predicted time: 1.007903e+10

Callcount: 28375652

__invivo__mapping_q1_include_base_tensor_base__ZN6TensorILi1ELi3EEaSERKS0__398

Real time: 4.319898e+08 -- Predicted time: 6.303630e+09

Callcount: 16561775

__invivo__sparse_matrix_double_include_lac_sparse_matrix_templates__ZNK12SparseMatrixIdE5vmultI6VectorIdES3_EEvRT_RKT0__348

Real time: 8.875573e+08 -- Predicted time: 9.106381e+08

Callcount: 668

__invivo__polynomial__ZNK11Polynomials10PolynomialIdE5valueEd_61

Real time: 1.230907e+09 -- Predicted time: 1.202738e+10

Callcount: 43046364

__invivo__fe_values_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470

Real time: 6.813907e+08 -- Predicted time: 1.391051e+10

Callcount: 28375652

__invivo__mapping_q1__ZNK9MappingQ1ILi3EE12compute_fillERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjN10QProjectorILi3EE17DataSetDescriptorERNS0_12InternalDataERSt6vectorI5PointILi3EESaISE_EE_511

Real time: 6.839042e+09 -- Predicted time: 8.004070e+09

Callcount: 200996

__invivo__mapping_q_include_fe_mapping_q__ZNK8MappingQILi3EE19transform_covariantEP6TensorILi1ELi3EES3_PKS2_RKN7MappingILi3EE16InternalDataBaseE_1159

Real time: 1.061856e+09 -- Predicted time: 2.817437e+09

Callcount: 3736156

__invivo__fe_include_base_tensor_base__ZN6TensorILi1ELi3EEC2Eb_330

Real time: 7.223488e+08 -- Predicted time: 4.667876e+09

Callcount: 30017177

__invivo__compressed_sparsity_pattern__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_iterator_base_funcs__ZSt9__advanceISt23_Rb_tree_const_iteratorIjElEvRT_T0_St26bidirectional_iterator_tag_141

Real time: 4.813426e+09 -- Predicted time: 2.089787e+10

Callcount: 14910185

__invivo__mapping_q_include_base_tensor__Z8contractILi3EEvR6TensorILi1EXT_EERKS1_RKS0_ILi2EXT_EE_502

Real time: 1.722230e+09 -- Predicted time: 1.163764e+10

Callcount: 63821468

__invivo__step_14__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_construct__ZNSt12_Destroy_auxILb0EE9__destroyIPSt3setIjSt4lessIjESaIjEEEEvT_S8__102

Real time: 3.677175e+08 -- Predicted time: 0.000000e+00

Callcount: 18

__invivo__compressed_sparsity_pattern__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_tree__ZNSt8_Rb_treeIjjSt9_IdentityIjESt4lessIjESaIjEE24_M_get_insert_unique_posERKj_1321

Real time: 1.393543e+09 -- Predicted time: 5.988079e+09

Callcount: 13395685

__invivo__step_14_include_base_tensor_base__ZNK6TensorILi1ELi3EEmlERKS0__493

Real time: 3.439522e+09 -- Predicted time: 3.286173e+10

Callcount: 143234298

__invivo__mapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EEC2ERKS0__350

Real time: 6.715702e+09 -- Predicted time: 1.143546e+11

Callcount: 279315570

__invivo__mapping_q_include_base_tensor__Z8contractILi3EEdRK6TensorILi1EXT_EES3__460

Real time: 5.469959e+08 -- Predicted time: 4.856730e+09

Callcount: 22781250

__invivo__fe__ZNK17FiniteElementBaseILi3EE11compute_2ndERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjRNS2_16InternalDataBaseERNS0_16InternalDataBaseER12FEValuesDataILi3EE_489

Real time: 3.816822e+09 -- Predicted time: 0.000000e+00

Callcount: 6664

__invivo__fe_include_base_tensor_base__ZN6TensorILi1ELi3EEaSERKS0__398

Real time: 7.204908e+08 -- Predicted time: 5.554267e+09

Callcount: 30004893

__invivo__fe_values_include_base_tensor_base__ZN6TensorILi1ELi3EEC2ERKS0__350

Real time: 6.820756e+08 -- Predicted time: 1.105129e+10

Callcount: 28405766

__invivo__tensor_product_polynomials__ZNK24TensorProductPolynomialsILi3EE13compute_valueEjRK5PointILi3EE_108

Real time: 1.509820e+09 -- Predicted time: 9.303152e+09

Callcount: 14347366

__invivo__fe_poly__ZNK7FE_PolyI24TensorProductPolynomialsILi3EELi3EE14fill_fe_valuesERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEERK10QuadratureILi3EERNS4_16InternalDataBaseESI_R12FEValuesDataILi3EE_252

Real time: 7.792083e+08 -- Predicted time: 3.152913e+07

Callcount: 167658

__invivo__mapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EEpLERKS0__448

Real time: 3.261883e+09 -- Predicted time: 1.046865e+10

Callcount: 135858852

__invivo__mapping_q__ZNK8MappingQILi3EE20apply_laplace_vectorERK5TableILi2EdERSt6vectorI5PointILi3EESaIS7_EE_792

Real time: 4.953406e+08 -- Predicted time: 5.913699e+08

Callcount: 210830

__invivo__mapping_q1_include_base_tensor_base__ZN6TensorILi1ELi3EEpLERKS0__448

Real time: 6.682630e+08 -- Predicted time: 5.569616e+09

Callcount: 27825616

__invivo__mapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EE5clearEv_534

Real time: 1.557901e+09 -- Predicted time: 2.300684e+10

Callcount: 64872848

__invivo__full_matrix_double_include_lac_full_matrix_templates__ZNK10FullMatrixIdE6TmmultIdEEvRS_IT_ERKS3_b_573

Real time: 2.545289e+09 -- Predicted time: 2.486162e+09

Callcount: 448

__invivo__mapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470

Real time: 3.261866e+09 -- Predicted time: 4.744976e+10

invivofe_values_include_base_tensor_baseZN6TensorILi1ELi3EEpLERKS0448

invivomapping_q1_include_base_tensor_baseZN6TensorILi1ELi3EEaSERKS0398

invivosparse_matrix_double_include_lac_sparse_matrix_templates__ZNK12SparseMatrixIdE5vmultI6VectorIdES3_EEvRT_RKT0__348

invivopolynomial__ZNK11Polynomials10PolynomialIdE5valueEd_61

invivofe_values_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470

invivomapping_q1__ZNK9MappingQ1ILi3EE12compute_fillERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjN10QProjectorILi3EE17DataSetDescriptorERNS0_12InternalDataERSt6vectorI5PointILi3EESaISE_EE_511

invivomapping_q_include_fe_mapping_q__ZNK8MappingQILi3EE19transform_covariantEP6TensorILi1ELi3EES3_PKS2_RKN7MappingILi3EE16InternalDataBaseE_1159

invivofe_include_base_tensor_base__ZN6TensorILi1ELi3EEC2Eb_330

invivocompressed_sparsity_pattern__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_iterator_base_funcsZSt9advanceISt23_Rb_tree_const_iteratorIjElEvRT_T0_St26bidirectional_iterator_tag_141

invivomapping_q_include_base_tensor__Z8contractILi3EEvR6TensorILi1EXT_EERKS1_RKS0_ILi2EXT_EE_502

invivostep_14__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_construct__ZNSt12_Destroy_auxILb0EE9__destroyIPSt3setIjSt4lessIjESaIjEEEEvT_S8__102

invivocompressed_sparsity_pattern__usr_lib_gcc_x86_64_linux_gnu_4_8_____________include_c___4_8_bits_stl_tree__ZNSt8_Rb_treeIjjSt9_IdentityIjESt4lessIjESaIjEE24_M_get_insert_unique_posERKj_1321

invivostep_14_include_base_tensor_baseZNK6TensorILi1ELi3EEmlERKS0493

invivomapping_q_include_base_tensor_baseZN6TensorILi1ELi3EEC2ERKS0350

invivomapping_q_include_base_tensor__Z8contractILi3EEdRK6TensorILi1EXT_EES3__460

invivofe__ZNK17FiniteElementBaseILi3EE11compute_2ndERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjRNS2_16InternalDataBaseERNS0_16InternalDataBaseER12FEValuesDataILi3EE_489

invivofe_include_base_tensor_baseZN6TensorILi1ELi3EEaSERKS0398

invivofe_values_include_base_tensor_baseZN6TensorILi1ELi3EEC2ERKS0350

invivotensor_product_polynomials__ZNK24TensorProductPolynomialsILi3EE13compute_valueEjRK5PointILi3EE_108

invivofe_poly__ZNK7FE_PolyI24TensorProductPolynomialsILi3EELi3EE14fill_fe_valuesERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEERK10QuadratureILi3EERNS4_16InternalDataBaseESI_R12FEValuesDataILi3EE_252

invivomapping_q_include_base_tensor_baseZN6TensorILi1ELi3EEpLERKS0448

invivomapping_q__ZNK8MappingQILi3EE20apply_laplace_vectorERK5TableILi2EdERSt6vectorI5PointILi3EESaIS7_EE_792

invivomapping_q1_include_base_tensor_baseZN6TensorILi1ELi3EEpLERKS0448

invivomapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EE5clearEv_534

invivofull_matrix_double_include_lac_full_matrix_templates__ZNK10FullMatrixIdE6TmmultIdEEvRS_IT_ERKS3_b_573

invivomapping_q_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470

invivomapping_q1_include_base_tensor_baseZN6TensorILi1ELi3EEC2ERKS0350

invivofe_poly__ZNK7FE_PolyI24TensorProductPolynomialsILi3EELi3EE19fill_fe_face_valuesERK7MappingILi3EERK12TriaIteratorILi3E15DoFCellAccessorILi3EEEjRK10QuadratureILi2EERNS4_16InternalDataBaseESI_R12FEValuesDataILi3EE_305

invivomapping_q1_include_base_tensor_base__ZN6TensorILi1ELi3EEmLERKd_470