Commit 6718826e authored by Santiago Ospina's avatar Santiago Ospina

comment out debug info

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 95d0a970
......@@ -187,13 +187,12 @@ public:
{
lfsv_volume.bind( eg.entity() );
std::cout << "lfsv_volume.size(): " << lfsv_volume.size() << std::endl;
std::cout << "lfsv_volume.debug():" << std::endl; lfsv_volume.debug();
// std::cout << "lfsv_volume.size(): " << lfsv_volume.size() << std::endl;
// std::cout << "lfsv_volume.debug():" << std::endl; lfsv_volume.debug();
global_pl_view.bind(lfsw_cache);
pl.resize(lfsw_cache.size());
rl.assign(lfsv_volume.size(),0.0);
}
......@@ -203,16 +202,16 @@ public:
lfsv_skeleton.bind( ig.intersection() );
// inside part
std::cout << "lfsu_cache.size(): " << lfsu_cache.size() << std::endl;
std::cout << "lfsu.debug():" << std::endl; lfsu_cache.localFunctionSpace().debug();
// std::cout << "lfsu_cache.size(): " << lfsu_cache.size() << std::endl;
// std::cout << "lfsu.debug():" << std::endl; lfsu_cache.localFunctionSpace().debug();
std::cout << "lfsw_cache.size(): " << lfsw_cache.size() << std::endl;
std::cout << "lfsw.debug():" << std::endl; lfsw_cache.localFunctionSpace().debug();
// std::cout << "lfsw_cache.size(): " << lfsw_cache.size() << std::endl;
// std::cout << "lfsw.debug():" << std::endl; lfsw_cache.localFunctionSpace().debug();
std::cout << "lfsv_skeleton.size(): " << lfsv_skeleton.size() << std::endl;
std::cout << "lfsv_skeleton.debug():" << std::endl; lfsv_skeleton.debug();
// std::cout << "lfsv_skeleton.size(): " << lfsv_skeleton.size() << std::endl;
// std::cout << "lfsv_skeleton.debug():" << std::endl; lfsv_skeleton.debug();
std::cout << "ig.entity().subEntities(1): " << ig.inside().subEntities(1) << std::endl;
// std::cout << "ig.entity().subEntities(1): " << ig.inside().subEntities(1) << std::endl;
}
template<typename IG, typename LFSUC, typename LFSWC>
......@@ -256,16 +255,16 @@ public:
const LFSUC & lfsu_cache,
const LFSVC & lfsv_cache)
{
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
// std::cout << "r_vec" << std::endl;
// std::cout << r_vec << std::endl;
// std::cout << "m_matrix" << std::endl;
// std::cout << m_matrix << std::endl;
x_vec.resize(lfsu_cache.size(),0.0);
m_matrix.solve(x_vec,r_vec);
std::cout << "x_vec" << std::endl;
std::cout << x_vec << std::endl;
// std::cout << "x_vec" << std::endl;
// std::cout << x_vec << std::endl;
auto& lfsu = lfsu_cache.localFunctionSpace();
......@@ -389,11 +388,11 @@ public:
m_matrix[i][j] += (phiu[j]*phiv[i])/detB * factor/detB;
}
std::cout << "assembleUVVolume" << std::endl;
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
// std::cout << "assembleUVVolume" << std::endl;
// std::cout << "r_vec" << std::endl;
// std::cout << r_vec << std::endl;
// std::cout << "m_matrix" << std::endl;
// std::cout << m_matrix << std::endl;
}
......@@ -408,7 +407,7 @@ public:
// if (not ig.intersection().conforming())
// std::cout << WARNING: Non-conforming entities produce non-continuos fluxes in normal direction of the entity << std::endl;
std::cout << "assembleUVSkeleton" << std::endl;
// std::cout << "assembleUVSkeleton" << std::endl;
rl_view.setWeight(local_assembler.weight());
rn_view.setWeight(local_assembler.weight());
......@@ -429,7 +428,7 @@ public:
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << " , " << rn(lfsv_skeleton_n,i) << std::endl;
// std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << " , " << rn(lfsv_skeleton_n,i) << std::endl;
r_vec[i+offset] += weight*rl(lfsv_skeleton,i);
}
......@@ -493,10 +492,10 @@ public:
}
}
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
// std::cout << "r_vec" << std::endl;
// std::cout << r_vec << std::endl;
// std::cout << "m_matrix" << std::endl;
// std::cout << m_matrix << std::endl;
}
template<typename IG, typename LFSWC>
......@@ -510,13 +509,13 @@ public:
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
alpha_boundary(lop,ig,lfsw_s_cache.localFunctionSpace(),pl,lfsv_skeleton,rl_view);
std::cout << "assembleUVBoundary" << std::endl;
// std::cout << "assembleUVBoundary" << std::endl;
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersection().indexInInside();
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << std::endl;
// std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << std::endl;
r_vec[i+offset] += rl(lfsv_skeleton,i);
}
......@@ -567,10 +566,10 @@ public:
m_matrix[i+offset][j] += (y*phiu_s[j])/detB * phiv_s[i] * factor;
}
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
// std::cout << "r_vec" << std::endl;
// std::cout << r_vec << std::endl;
// std::cout << "m_matrix" << std::endl;
// std::cout << m_matrix << std::endl;
}
template<typename IG, typename LFSWC>
......
// -*- tab-width: 2; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_RAVIART_THOMAS_BASIS_FACTORY_HH
#define DUNE_DORIE_RAVIART_THOMAS_BASIS_FACTORY_HH
// #include "raviartthomas2cube3d.hh"
#include <dune/localfunctions/lagrange/pk.hh>
#include <dune/localfunctions/lagrange/qk.hh>
#include <dune/localfunctions/raviartthomas.hh>
#include <dune/common/power.hh>
//======================================================================
// provide size and types depending on element type, dimension and order
//======================================================================
namespace Dune{
namespace Dorie{
template<class D, class R, int k, int d, Dune::GeometryType::BasicType gt>
struct RTBasisFactory
{
};
template<class D, class R, int k, int d>
struct RTBasisFactory<D,R,k,d,Dune::GeometryType::simplex>
{
static constexpr int order = k;
static constexpr int dim = d;
static constexpr int faces = d+1;
static constexpr int facesize = Dune::PB::PkSize<order,dim-1>::value;
static constexpr int size = dim == 2 ? (order+1)*(order+3) : (order+1)*(order+2)*(order+4)/2;
static constexpr int elementsize = size - faces*facesize;
typedef Dune::RaviartThomasSimplexLocalFiniteElement<dim,D,R> LFE;
static LFE get ()
{
return {Dune::GeometryType(Dune::GeometryType::simplex,dim),order};
}
typedef Dune::PkLocalFiniteElement<D,R,dim-1,order> TESTLFEFACE;
};
template<class D, class R, int k, int d>
struct RTBasisFactory<D,R,k,d,Dune::GeometryType::cube>
{
static constexpr int order = k;
static constexpr int dim = d;
static constexpr int faces = 2*dim;
static constexpr int facesize = Dune::StaticPower<order+1,dim-1>::power;
static constexpr int size = dim*facesize*(order+2);
static constexpr int elementsize = size - faces*facesize;
typedef Dune::RaviartThomasCubeLocalFiniteElement<D,R,dim,order> LFE;
static LFE get()
{
return {};
}
typedef Dune::QkLocalFiniteElement<D,R,dim-1,order> TESTLFEFACE;
};
//======================================================================
// the RT element test spaces
//======================================================================
template<class D, class R, int k, int d, Dune::GeometryType::BasicType gt>
struct RTPsi
{
};
template<class D, class R>
struct RTPsi<D,R,0,2,Dune::GeometryType::simplex>
{
static constexpr int order = 0;
static constexpr int dim = 2;
static constexpr int elementsize = dim == 2 ? dim*order*(order+1)/2 : dim*order*(order+1)*(order+2)/6;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{}
};
template<class D, class R>
struct RTPsi<D,R,0,3,Dune::GeometryType::simplex>
{
static constexpr int order = 0;
static constexpr int dim = 3;
static constexpr int elementsize = dim == 2 ? dim*order*(order+1)/2 : dim*order*(order+1)*(order+2)/6;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{}
};
template<class D, class R>
struct RTPsi<D,R,1,2,Dune::GeometryType::simplex>
{
static constexpr int order = 1;
static constexpr int dim = 2;
static constexpr int elementsize = dim == 2 ? dim*order*(order+1)/2 : dim*order*(order+1)*(order+2)/6;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
psi[0][0] = 1.0; psi[1][0] = 0.0;
psi[0][1] = 0.0; psi[1][1] = 1.0;
}
};
template<class D, class R>
struct RTPsi<D,R,1,3,Dune::GeometryType::simplex>
{
static constexpr int order = 1;
static constexpr int dim = 3;
static constexpr int elementsize = dim == 2 ? dim*order*(order+1)/2 : dim*order*(order+1)*(order+2)/6;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
psi[0][0] = 1.0; psi[1][0] = 0.0; psi[2][0] = 0.0;
psi[0][1] = 0.0; psi[1][1] = 1.0; psi[2][1] = 0.0;
psi[0][2] = 0.0; psi[1][2] = 0.0; psi[2][2] = 1.0;
}
};
template<class D, class R>
struct RTPsi<D,R,2,2,Dune::GeometryType::simplex>
{
static constexpr int order = 2;
static constexpr int dim = 2;
static constexpr int elementsize = dim == 2 ? dim*order*(order+1)/2 : dim*order*(order+1)*(order+2)/6;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
psi[0][0] = 1.0; psi[1][0] = x[0]; psi[2][0] = x[1];
psi[0][1] = 0.0; psi[1][1] = 0.0; psi[2][1] = 0.0;
psi[3][0] = 0.0; psi[4][0] = 0.0; psi[5][0] = 0.0;
psi[3][1] = 1.0; psi[4][1] = x[0]; psi[5][1] = x[1];
}
};
template<class D, class R>
struct RTPsi<D,R,2,3,Dune::GeometryType::simplex>
{
static constexpr int order = 2;
static constexpr int dim = 3;
static constexpr int elementsize = dim == 2 ? dim*order*(order+1)/2 : dim*order*(order+1)*(order+2)/6;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
psi[0][0] = 1.0; psi[1][0] = x[0]; psi[2][0] = x[1]; psi[3][0] = x[2];
psi[0][1] = 0.0; psi[1][1] = 0.0; psi[2][1] = 0.0; psi[3][1] = 0.0;
psi[0][2] = 0.0; psi[1][2] = 0.0; psi[2][2] = 0.0; psi[3][2] = 0.0;
psi[4][0] = 0.0; psi[5][0] = 0.0; psi[6][0] = 0.0; psi[7][0] = 0.0;
psi[4][1] = 1.0; psi[5][1] = x[0]; psi[6][1] = x[1]; psi[7][1] = x[2];
psi[4][2] = 0.0; psi[5][2] = 0.0; psi[6][2] = 0.0; psi[7][2] = 0.0;
psi[8][0] = 0.0; psi[9][0] = 0.0; psi[10][0] = 0.0; psi[11][0] = 0.0;
psi[8][1] = 0.0; psi[9][1] = 0.0; psi[10][1] = 0.0; psi[11][1] = 0.0;
psi[8][2] = 1.0; psi[9][2] = x[0]; psi[10][2] = x[1]; psi[11][2] = x[2];
}
};
//********************** now cubes **********************//
template<class D, class R>
struct RTPsi<D,R,0,2,Dune::GeometryType::cube>
{
static constexpr int order = 0;
static constexpr int dim = 2;
static constexpr int elementsize = dim*Dune::StaticPower<order+1,dim-1>::power*order;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{}
};
template<class D, class R>
struct RTPsi<D,R,0,3,Dune::GeometryType::cube>
{
static constexpr int order = 0;
static constexpr int dim = 3;
static constexpr int elementsize = dim*Dune::StaticPower<order+1,dim-1>::power*order;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{}
};
template<class D, class R>
struct RTPsi<D,R,1,2,Dune::GeometryType::cube>
{
static constexpr int order = 1;
static constexpr int dim = 2;
static constexpr int elementsize = dim*Dune::StaticPower<order+1,dim-1>::power*order;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
psi[0][0] = 1.0; psi[1][0] = x[1]; psi[2][0] = 0.0; psi[3][0] = 0.0;
psi[0][1] = 0.0; psi[1][1] = 0.0; psi[2][1] = 1.0; psi[3][1] = x[0];
}
};
template<class D, class R>
struct RTPsi<D,R,1,3,Dune::GeometryType::cube>
{
static constexpr int order = 1;
static constexpr int dim = 3;
static constexpr int elementsize = dim*Dune::StaticPower<order+1,dim-1>::power*order;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
psi[0] = {1.0, 0.0, 0.0};
psi[1] = {x[1], 0.0, 0.0};
psi[2] = {x[2], 0.0, 0.0};
psi[3] = {x[1]*x[2], 0.0, 0.0};
psi[4] = {0.0, 1.0, 0.0};
psi[5] = {0.0, x[0], 0.0};
psi[6] = {0.0, x[2], 0.0};
psi[7] = {0.0, x[0]*x[2], 0.0};
psi[8] = {0.0, 0.0, 1.0};
psi[9] = {0.0, 0.0, x[0]};
psi[10] = {0.0, 0.0, x[1]};
psi[11] = {0.0, 0.0, x[0]*x[1]};
#if 0
psi[0][0] = 1.0; psi[1][0] = x[1]; psi[2][0] = x[2]; psi[3][0] = x[1]*x[2];
psi[0][1] = 0.0; psi[1][1] = 0.0; psi[2][1] = 0.0; psi[3][1] = 0.0;
psi[0][2] = 0.0; psi[1][2] = 0.0; psi[2][2] = 0.0; psi[3][2] = 0.0;
psi[4][0] = 0.0; psi[5][0] = 0.0; psi[6][0] = 0.0; psi[7][0] = 0.0;
psi[4][1] = 1.0; psi[5][1] = x[0]; psi[6][1] = x[2]; psi[7][1] = x[0]*x[2];
psi[4][2] = 0.0; psi[5][2] = 0.0; psi[6][2] = 0.0; psi[7][2] = 0.0;
psi[8][0] = 0.0; psi[9][0] = 0.0; psi[10][0] = 0.0; psi[11][0] = 0.0;
psi[8][1] = 0.0; psi[9][1] = 0.0; psi[10][1] = 0.0; psi[11][1] = 0.0;
psi[8][2] = 1.0; psi[9][2] = x[0]; psi[10][2] = x[1]; psi[11][2] = x[0]*x[1];
#endif
}
};
template<class D, class R>
struct RTPsi<D,R,2,2,Dune::GeometryType::cube>
{
static constexpr int order = 2;
static constexpr int dim = 2;
static constexpr int elementsize = dim*Dune::StaticPower<order+1,dim-1>::power*order;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
psi[0][0] = 1.0; psi[1][0] = x[1]; psi[2][0] = x[1]*x[1];
psi[0][1] = 0.0; psi[1][1] = 0.0; psi[2][1] = 0.0;
psi[3][0] = 0.0; psi[4][0] = 0.0; psi[5][0] = 0.0;
psi[3][1] = 1.0; psi[4][1] = x[0]; psi[5][1] = x[0]*x[0];
psi[6][0] = x[0]; psi[7][0] = x[0]*x[1]; psi[8][0] = x[0]*x[1]*x[1];
psi[6][1] = 0.0; psi[7][1] = 0.0; psi[8][1] = 0.0;
psi[9][0] = 0.0; psi[10][0] = 0.0; psi[11][0] = 0.0;
psi[9][1] = x[1]; psi[10][1] = x[1]*x[0]; psi[11][1] = x[1]*x[0]*x[0];
}
};
template<class D, class R>
struct RTPsi<D,R,2,3,Dune::GeometryType::cube>
{
static constexpr int order = 2;
static constexpr int dim = 3;
static constexpr int elementsize = dim*Dune::StaticPower<order+1,dim-1>::power*order;
typedef Dune::FieldVector<D,dim> DomainType;
typedef Dune::FieldVector<R,dim> RangeType;
Dune::QkStuff::GaussLobattoLagrangePolynomials<D,R,order> poly_k;
Dune::QkStuff::GaussLobattoLagrangePolynomials<D,R,order-1> poly_km1;
constexpr unsigned int size() const
{
return elementsize;
}
inline void evaluateFunction (const DomainType& x, std::vector<RangeType>& psi) const
{
psi.resize(elementsize);
// x-direction fastest, y-direction middle, z-direction slowest
std::size_t global_iter = 0.0;
for(std::size_t i2=0; i2<order+1; ++i2)
for(std::size_t i1=0; i1<order+1; ++i1)
for(std::size_t i0=0; i0<order; ++i0) // reduced
{
psi[global_iter] = {poly_km1.p(i0,x[0])*poly_k.p(i1,x[1])*poly_k.p(i2,x[2]), 0.0, 0.0};
global_iter++;
}
for(std::size_t i2=0; i2<order+1; ++i2)
for(std::size_t i1=0; i1<order; ++i1) // reduced
for(std::size_t i0=0; i0<order+1; ++i0)
{
psi[global_iter] = {0.0, poly_k.p(i0,x[0])*poly_km1.p(i1,x[1])*poly_k.p(i2,x[2]), 0.0};
global_iter++;
}
for(std::size_t i2=0; i2<order; ++i2) // reduced
for(std::size_t i1=0; i1<order+1; ++i1)
for(std::size_t i0=0; i0<order+1; ++i0)
{
psi[global_iter] = {0.0, 0.0, poly_k.p(i0,x[0])*poly_k.p(i1,x[1])*poly_km1.p(i2,x[2])};
global_iter++;
}
#if 0
// x-direction fastest, y-direction middle, z-direction slowest
psi[0] = {1.0, 0.0, 0.0};
psi[1] = {x[0], 0.0, 0.0};
psi[2] = {x[1], 0.0, 0.0};
psi[3] = {x[0]*x[1], 0.0, 0.0};
psi[4] = {x[1]*x[1], 0.0, 0.0};
psi[5] = {x[0]*x[1]*x[1], 0.0, 0.0};
//--------------------------------------------
psi[6] = {x[2], 0.0, 0.0};
psi[7] = {x[0]*x[2], 0.0, 0.0};
psi[8] = {x[1]*x[2], 0.0, 0.0};
psi[9] = {x[0]*x[1]*x[2], 0.0, 0.0};
psi[10] = {x[1]*x[1]*x[2], 0.0, 0.0};
psi[11] = {x[0]*x[1]*x[1]*x[2], 0.0, 0.0};
//--------------------------------------------
psi[12] = {x[2]*x[2], 0.0, 0.0};
psi[13] = {x[0]*x[2]*x[2], 0.0, 0.0};
psi[14] = {x[1]*x[2]*x[2], 0.0, 0.0};
psi[15] = {x[0]*x[1]*x[2]*x[2], 0.0, 0.0};
psi[16] = {x[1]*x[1]*x[2]*x[2], 0.0, 0.0};
psi[17] = {x[0]*x[1]*x[1]*x[2]*x[2], 0.0, 0.0};
psi[18] = {0.0, 1.0, 0.0};
psi[19] = {0.0, x[0], 0.0};
psi[20] = {0.0, x[0]*x[0], 0.0};
psi[21] = {0.0, x[1], 0.0};
psi[22] = {0.0, x[0]*x[1], 0.0};
psi[23] = {0.0, x[0]*x[0]*x[1], 0.0};
//--------------------------------------------
psi[24] = {0.0, x[2], 0.0};
psi[25] = {0.0, x[0]*x[2], 0.0};
psi[26] = {0.0, x[0]*x[0]*x[2], 0.0};
psi[27] = {0.0, x[1]*x[2], 0.0};
psi[28] = {0.0, x[0]*x[1]*x[2], 0.0};
psi[29] = {0.0, x[0]*x[0]*x[1]*x[2], 0.0};
//--------------------------------------------
psi[30] = {0.0, x[2]*x[2], 0.0};
psi[31] = {0.0, x[0]*x[2]*x[2], 0.0};
psi[32] = {0.0, x[0]*x[0]*x[2]*x[2], 0.0};
psi[33] = {0.0, x[1]*x[2]*x[2], 0.0};
psi[34] = {0.0, x[0]*x[1]*x[2]*x[2], 0.0};
psi[35] = {0.0, x[0]*x[0]*x[1]*x[2]*x[2], 0.0};
psi[36] = {0.0, 0.0, 1.0};
psi[37] = {0.0, 0.0, x[0]};
psi[38] = {0.0, 0.0, x[0]*x[0]};
psi[39] = {0.0, 0.0, x[1]};
psi[40] = {0.0, 0.0, x[0]*x[1]};
psi[41] = {0.0, 0.0, x[0]*x[0]*x[1]};
psi[42] = {0.0, 0.0, x[1]*x[1]};
psi[43] = {0.0, 0.0, x[0]*x[1]*x[1]};
psi[44] = {0.0, 0.0, x[0]*x[0]*x[1]*x[1]};
//--------------------------------------------
psi[45] = {0.0, 0.0, x[2]};
psi[46] = {0.0, 0.0, x[0]*x[2]};
psi[47] = {0.0, 0.0, x[0]*x[0]*x[2]};
psi[48] = {0.0, 0.0, x[1]*x[2]};
psi[49] = {0.0, 0.0, x[0]*x[1]*x[2]};
psi[50] = {0.0, 0.0, x[0]*x[0]*x[1]*x[2]};
psi[51] = {0.0, 0.0, x[1]*x[1]*x[2]};
psi[52] = {0.0, 0.0, x[0]*x[1]*x[1]*x[2]};
psi[53] = {0.0, 0.0, x[0]*x[0]*x[1]*x[1]*x[2]};
#endif
#if 0
psi[0][0] = 1.0; psi[1][0] = x[1]; psi[2][0] = x[1]*x[1]; psi[3][0] = x[2]; psi[4][0] = x[2]*x[2]; psi[5][0] = x[1]*x[2]; psi[6][0] = x[1]*x[2]*x[2]; psi[7][0] = x[1]*x[1]*x[2]; psi[8][0] = x[1]*x[1]*x[2]*x[2];
psi[0][1] = 0.0; psi[1][1] = 0.0; psi[2][1] = 0.0; psi[3][1] = 0.0; psi[4][1] = 0.0; psi[5][1] = 0.0; psi[6][1] = 0.0; psi[7][1] = 0.0; psi[8][1] = 0.0;
psi[0][2] = 0.0; psi[1][2] = 0.0; psi[2][2] = 0.0; psi[3][2] = 0.0; psi[4][2] = 0.0; psi[5][2] = 0.0; psi[6][2] = 0.0; psi[7][2] = 0.0; psi[8][2] = 0.0;
psi[9][0] = 0.0; psi[10][0] = 0.0; psi[11][0] = 0.0; psi[12][0] = 0.0; psi[13][0] = 0.0; psi[14][0] = 0.0; psi[15][0] = 0.0; psi[16][0] = 0.0; psi[17][0] = 0.0;
psi[9][1] = 1.0; psi[10][1] = x[0]; psi[11][1] = x[0]*x[0]; psi[12][1] = x[2]; psi[13][1] = x[2]*x[2]; psi[14][1] = x[0]*x[2]; psi[15][1] = x[0]*x[2]*x[2]; psi[16][1] = x[0]*x[0]*x[2]; psi[17][1] = x[0]*x[0]*x[2]*x[2];
psi[9][2] = 0.0; psi[10][2] = 0.0; psi[11][2] = 0.0; psi[12][2] = 0.0; psi[13][2] = 0.0; psi[14][2] = 0.0; psi[15][2] = 0.0; psi[16][2] = 0.0; psi[17][2] = 0.0;
psi[18][0] = 0.0; psi[19][0] = 0.0; psi[20][0] = 0.0; psi[21][0] = 0.0; psi[22][0] = 0.0; psi[23][0] = 0.0; psi[24][0] = 0.0; psi[25][0] = 0.0; psi[26][0] = 0.0;
psi[18][1] = 0.0; psi[19][1] = 0.0; psi[20][1] = 0.0; psi[21][1] = 0.0; psi[22][1] = 0.0; psi[23][1] = 0.0; psi[24][1] = 0.0; psi[25][1] = 0.0; psi[26][1] = 0.0;
psi[18][2] = 1.0; psi[19][2] = x[0]; psi[20][2] = x[0]*x[0]; psi[21][2] = x[1]; psi[22][2] = x[1]*x[1]; psi[23][2] = x[0]*x[1]; psi[24][2] = x[0]*x[1]*x[1]; psi[25][2] = x[0]*x[0]*x[1]; psi[26][2] = x[0]*x[0]*x[1]*x[1];
psi[27][0] = x[0]; psi[28][0] = x[0]*x[1]; psi[29][0] = x[0]*x[1]*x[1]; psi[30][0] = x[0]*x[2]; psi[31][0] = x[0]*x[2]*x[2]; psi[32][0] = x[0]*x[1]*x[2]; psi[33][0] = x[0]*x[1]*x[2]*x[2]; psi[34][0] = x[0]*x[1]*x[1]*x[2]; psi[35][0] = x[0]*x[1]*x[1]*x[2]*x[2];
psi[27][1] = 0.0; psi[28][1] = 0.0; psi[29][1] = 0.0; psi[30][1] = 0.0; psi[31][1] = 0.0; psi[32][1] = 0.0; psi[33][1] = 0.0; psi[34][1] = 0.0; psi[35][1] = 0.0;
psi[27][2] = 0.0; psi[28][2] = 0.0; psi[29][2] = 0.0; psi[30][2] = 0.0; psi[31][2] = 0.0; psi[32][2] = 0.0; psi[33][2] = 0.0; psi[34][2] = 0.0; psi[35][2] = 0.0;