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?

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

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