//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: Q2HexahedronFiniteElement.cpp,v 1.2 2006/04/08 10:43:20 delpinux Exp $

#include <Q2HexahedronFiniteElement.hpp>

TinyVector<3, real_t> Q2HexahedronFiniteElement::__massCenter(0.5, 0.5, 0.5);

const size_t
Q2HexahedronFiniteElement::
facesDOF[Hexahedron::NumberOfFaces][Q2HexahedronFiniteElement::numberOfFaceLivingDegreesOfFreedom]
= {{ 0, 1, 2, 3, 8, 9,10,11,20},
   { 0, 1, 4, 5, 8,12,13,16,21},
   { 1, 2, 5, 6, 9,13,14,17,22},
   { 2, 3, 6, 7,10,14,15,18,23},
   { 0, 3, 4, 7,11,12,15,19,24},
   { 4, 5, 6, 7,16,17,18,19,25}};


real_t
Q2HexahedronFiniteElement::W(const size_t& i, const TinyVector<3, real_t>& X) const
{
  const real_t& x = X[0];
  const real_t& y = X[1];
  const real_t& z = X[2];

  switch (i) {
    // Vertices basis functions 
  case 0: {
    return __w1(x)*__w1(y)*__w1(z);
  }
  case 1: {
    return __w3(x)*__w1(y)*__w1(z);
  }
  case 2: {
    return __w3(x)*__w3(y)*__w1(z);
  }
  case 3: {
    return __w1(x)*__w3(y)*__w1(z);
  }
  case 4: {
    return __w1(x)*__w1(y)*__w3(z);
  }
  case 5: {
    return __w3(x)*__w1(y)*__w3(z);
  }
  case 6: {
    return __w3(x)*__w3(y)*__w3(z);
  }
  case 7: {
    return __w1(x)*__w3(y)*__w3(z);
  }
    // Edges basis functions
  case 8: {
    return __w2(x)*__w1(y)*__w1(z);
  }
  case 9: {
    return __w3(x)*__w2(y)*__w1(z);
  }
  case 10: {
    return __w2(x)*__w3(y)*__w1(z);
  }
  case 11: {
    return __w1(x)*__w2(y)*__w1(z);
  }
  case 12: {
    return __w1(x)*__w1(y)*__w2(z);
  }
  case 13: {
    return __w3(x)*__w1(y)*__w2(z);
  }
  case 14: {
    return __w3(x)*__w3(y)*__w2(z);
  }
  case 15: {
    return __w1(x)*__w3(y)*__w2(z);
  }
  case 16: {
    return __w2(x)*__w1(y)*__w3(z);
  }
  case 17: {
    return __w3(x)*__w2(y)*__w3(z);
  }
  case 18: {
    return __w2(x)*__w3(y)*__w3(z);
  }
  case 19: {
    return __w1(x)*__w2(y)*__w3(z);
  }
    // Faces basis functions
  case 20: {
    return __w2(x)*__w2(y)*__w1(z);
  }
  case 21: {
    return __w2(x)*__w1(y)*__w2(z);
  }
  case 22: {
    return __w3(x)*__w2(y)*__w2(z);
  }
  case 23: {
    return __w2(x)*__w3(y)*__w2(z);
  }
  case 24: {
    return __w1(x)*__w2(y)*__w2(z);
  }
  case 25: {
    return __w2(x)*__w2(y)*__w3(z);
  }
    // Volume basis functions
  case 26: {
    return __w2(x)*__w2(y)*__w2(z);
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected basis function number",
		       ErrorHandler::unexpected);
    return 0.;
  }
  }
}

real_t
Q2HexahedronFiniteElement::dxW(const size_t& i, const TinyVector<3, real_t>& X) const
{
  const real_t& x = X[0];
  const real_t& y = X[1];
  const real_t& z = X[2];

  switch (i) {
    // Vertices basis functions 
  case 0: {
    return __dw1(x)*__w1(y)*__w1(z);
  }
  case 1: {
    return __dw3(x)*__w1(y)*__w1(z);
  }
  case 2: {
    return __dw3(x)*__w3(y)*__w1(z);
  }
  case 3: {
    return __dw1(x)*__w3(y)*__w1(z);
  }
  case 4: {
    return __dw1(x)*__w1(y)*__w3(z);
  }
  case 5: {
    return __dw3(x)*__w1(y)*__w3(z);
  }
  case 6: {
    return __dw3(x)*__w3(y)*__w3(z);
  }
  case 7: {
    return __dw1(x)*__w3(y)*__w3(z);
  }
    // Edges basis functions
  case 8: {
    return __dw2(x)*__w1(y)*__w1(z);
  }
  case 9: {
    return __dw3(x)*__w2(y)*__w1(z);
  }
  case 10: {
    return __dw2(x)*__w3(y)*__w1(z);
  }
  case 11: {
    return __dw1(x)*__w2(y)*__w1(z);
  }
  case 12: {
    return __dw1(x)*__w1(y)*__w2(z);
  }
  case 13: {
    return __dw3(x)*__w1(y)*__w2(z);
  }
  case 14: {
    return __dw3(x)*__w3(y)*__w2(z);
  }
  case 15: {
    return __dw1(x)*__w3(y)*__w2(z);
  }
  case 16: {
    return __dw2(x)*__w1(y)*__w3(z);
  }
  case 17: {
    return __dw3(x)*__w2(y)*__w3(z);
  }
  case 18: {
    return __dw2(x)*__w3(y)*__w3(z);
  }
  case 19: {
    return __dw1(x)*__w2(y)*__w3(z);
  }
    // Faces basis functions
  case 20: {
    return __dw2(x)*__w2(y)*__w1(z);
  }
  case 21: {
    return __dw2(x)*__w1(y)*__w2(z);
  }
  case 22: {
    return __dw3(x)*__w2(y)*__w2(z);
  }
  case 23: {
    return __dw2(x)*__w3(y)*__w2(z);
  }
  case 24: {
    return __dw1(x)*__w2(y)*__w2(z);
  }
  case 25: {
    return __dw2(x)*__w2(y)*__w3(z);
  }
    // Volume basis functions
  case 26: {
    return __dw2(x)*__w2(y)*__w2(z);
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected basis function number",
		       ErrorHandler::unexpected);
    return 0.;
  }
  }
}

real_t
Q2HexahedronFiniteElement::dyW(const size_t& i, const TinyVector<3, real_t>& X) const
{
  const real_t& x = X[0];
  const real_t& y = X[1];
  const real_t& z = X[2];

  switch (i) {
    // Vertices basis functions 
  case 0: {
    return __w1(x)*__dw1(y)*__w1(z);
  }
  case 1: {
    return __w3(x)*__dw1(y)*__w1(z);
  }
  case 2: {
    return __w3(x)*__dw3(y)*__w1(z);
  }
  case 3: {
    return __w1(x)*__dw3(y)*__w1(z);
  }
  case 4: {
    return __w1(x)*__dw1(y)*__w3(z);
  }
  case 5: {
    return __w3(x)*__dw1(y)*__w3(z);
  }
  case 6: {
    return __w3(x)*__dw3(y)*__w3(z);
  }
  case 7: {
    return __w1(x)*__dw3(y)*__w3(z);
  }
    // Edges basis functions
  case 8: {
    return __w2(x)*__dw1(y)*__w1(z);
  }
  case 9: {
    return __w3(x)*__dw2(y)*__w1(z);
  }
  case 10: {
    return __w2(x)*__dw3(y)*__w1(z);
  }
  case 11: {
    return __w1(x)*__dw2(y)*__w1(z);
  }
  case 12: {
    return __w1(x)*__dw1(y)*__w2(z);
  }
  case 13: {
    return __w3(x)*__dw1(y)*__w2(z);
  }
  case 14: {
    return __w3(x)*__dw3(y)*__w2(z);
  }
  case 15: {
    return __w1(x)*__dw3(y)*__w2(z);
  }
  case 16: {
    return __w2(x)*__dw1(y)*__w3(z);
  }
  case 17: {
    return __w3(x)*__dw2(y)*__w3(z);
  }
  case 18: {
    return __w2(x)*__dw3(y)*__w3(z);
  }
  case 19: {
    return __w1(x)*__dw2(y)*__w3(z);
  }
    // Faces basis functions
  case 20: {
    return __w2(x)*__dw2(y)*__w1(z);
  }
  case 21: {
    return __w2(x)*__dw1(y)*__w2(z);
  }
  case 22: {
    return __w3(x)*__dw2(y)*__w2(z);
  }
  case 23: {
    return __w2(x)*__dw3(y)*__w2(z);
  }
  case 24: {
    return __w1(x)*__dw2(y)*__w2(z);
  }
  case 25: {
    return __w2(x)*__dw2(y)*__w3(z);
  }
    // Volume basis functions
  case 26: {
    return __w2(x)*__dw2(y)*__w2(z);
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected basis function number",
		       ErrorHandler::unexpected);
    return 0.;
  }
  }
}

real_t 
Q2HexahedronFiniteElement::dzW(const size_t& i, const TinyVector<3, real_t>& X) const
{
  const real_t& x = X[0];
  const real_t& y = X[1];
  const real_t& z = X[2];

  switch (i) {
    // Vertices basis functions 
  case 0: {
    return __w1(x)*__w1(y)*__dw1(z);
  }
  case 1: {
    return __w3(x)*__w1(y)*__dw1(z);
  }
  case 2: {
    return __w3(x)*__w3(y)*__dw1(z);
  }
  case 3: {
    return __w1(x)*__w3(y)*__dw1(z);
  }
  case 4: {
    return __w1(x)*__w1(y)*__dw3(z);
  }
  case 5: {
    return __w3(x)*__w1(y)*__dw3(z);
  }
  case 6: {
    return __w3(x)*__w3(y)*__dw3(z);
  }
  case 7: {
    return __w1(x)*__w3(y)*__dw3(z);
  }
    // Edges basis functions
  case 8: {
    return __w2(x)*__w1(y)*__dw1(z);
  }
  case 9: {
    return __w3(x)*__w2(y)*__dw1(z);
  }
  case 10: {
    return __w2(x)*__w3(y)*__dw1(z);
  }
  case 11: {
    return __w1(x)*__w2(y)*__dw1(z);
  }
  case 12: {
    return __w1(x)*__w1(y)*__dw2(z);
  }
  case 13: {
    return __w3(x)*__w1(y)*__dw2(z);
  }
  case 14: {
    return __w3(x)*__w3(y)*__dw2(z);
  }
  case 15: {
    return __w1(x)*__w3(y)*__dw2(z);
  }
  case 16: {
    return __w2(x)*__w1(y)*__dw3(z);
  }
  case 17: {
    return __w3(x)*__w2(y)*__dw3(z);
  }
  case 18: {
    return __w2(x)*__w3(y)*__dw3(z);
  }
  case 19: {
    return __w1(x)*__w2(y)*__dw3(z);
  }
    // Faces basis functions
  case 20: {
    return __w2(x)*__w2(y)*__dw1(z);
  }
  case 21: {
    return __w2(x)*__w1(y)*__dw2(z);
  }
  case 22: {
    return __w3(x)*__w2(y)*__dw2(z);
  }
  case 23: {
    return __w2(x)*__w3(y)*__dw2(z);
  }
  case 24: {
    return __w1(x)*__w2(y)*__dw2(z);
  }
  case 25: {
    return __w2(x)*__w2(y)*__dw3(z);
  }
    // Volume basis functions
  case 26: {
    return __w2(x)*__w2(y)*__dw2(z);
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected basis function number",
		       ErrorHandler::unexpected);
    return 0.;
  }
  }
}