#ifndef _Hex8_hpp_

#define _Hex8_hpp_

//@HEADER

// ************************************************************************

//

// MiniFE: Simple Finite Element Assembly and Solve

// Copyright (2006-2013) Sandia Corporation

//

// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive

// license for use of this work by or on behalf of the U.S. Government.

//

// This library is free software; you can redistribute it and/or modify

// it under the terms of the GNU Lesser General Public License as

// published by the Free Software Foundation; either version 2.1 of the

// License, or (at your option) any later version.

//

// This library 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

// Lesser General Public License for more details.

//

// You should have received a copy of the GNU Lesser General Public

// License along with this library; if not, write to the Free Software

// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307

// USA

//

// ************************************************************************

//@HEADER

#ifndef KERNEL_PREFIX

#define KERNEL_PREFIX

#endif

#include <gauss_pts.hpp>

#include <matrix_algebra_3x3.hpp>

#include <Hex8_enums.hpp>

namespace miniFE {

namespace Hex8 {

template<typename Scalar>

KERNEL_PREFIX void shape_fns(const Scalar* x, Scalar* values_at_nodes)

{

//assumptions: values_at_nodes has length numNodesPerElem

// x has length 3 (hard-coded spatialDim)

const Scalar u = 1.0 - x[0];

const Scalar v = 1.0 - x[1];

const Scalar w = 1.0 - x[2];

const Scalar up1 = 1.0 + x[0];

const Scalar vp1 = 1.0 + x[1];

const Scalar wp1 = 1.0 + x[2];

values_at_nodes[0] = 0.125 * u * v * w;//(1-x)*(1-y)*(1-z)

values_at_nodes[1] = 0.125 * up1 * v * w;//(1+x)*(1-y)*(1-z)

values_at_nodes[2] = 0.125 * up1 * vp1 * w;//(1+x)*(1+y)*(1-z)

values_at_nodes[3] = 0.125 * u * vp1 * w;//(1-x)*(1+y)*(1-z)

values_at_nodes[4] = 0.125 * u * v * wp1;//(1-x)*(1-y)*(1+z)

values_at_nodes[5] = 0.125 * up1 * v * wp1;//(1+x)*(1-y)*(1+z)

values_at_nodes[6] = 0.125 * up1 * vp1 * wp1;//(1+x)*(1+y)*(1+z)

values_at_nodes[7] = 0.125 * u * vp1 * wp1;//(1-x)*(1+y)*(1+z)

}

template<typename Scalar>

KERNEL_PREFIX void gradients(const Scalar* x, Scalar* values_per_fn)

{

//assumptions values_per_fn has length 24 (numNodesPerElem*spatialDim)

// spatialDim == 3

const Scalar u = 1.0 - x[0];

const Scalar v = 1.0 - x[1];

const Scalar w = 1.0 - x[2];

const Scalar up1 = 1.0 + x[0];

const Scalar vp1 = 1.0 + x[1];

const Scalar wp1 = 1.0 + x[2];

//fn 0

values_per_fn[0] = -0.125 * v * w;

values_per_fn[1] = -0.125 * u * w;

values_per_fn[2] = -0.125 * u * v;

//fn 1

values_per_fn[3] = 0.125 * v * w;

values_per_fn[4] = -0.125 * up1 * w;

values_per_fn[5] = -0.125 * up1 * v;

//fn 2

values_per_fn[6] = 0.125 * vp1 * w;

values_per_fn[7] = 0.125 * up1 * w;

values_per_fn[8] = -0.125 * up1 * vp1;

//fn 3

values_per_fn[9] = -0.125 * vp1 * w;

values_per_fn[10] = 0.125 * u * w;

values_per_fn[11] = -0.125 * u * vp1;

//fn 4

values_per_fn[12] = -0.125 * v * wp1;

values_per_fn[13] = -0.125 * u * wp1;

values_per_fn[14] = 0.125 * u * v;

//fn 5

values_per_fn[15] = 0.125 * v * wp1;

values_per_fn[16] = -0.125 * up1 * wp1;

values_per_fn[17] = 0.125 * up1 * v;

//fn 6

values_per_fn[18] = 0.125 * vp1 * wp1;

values_per_fn[19] = 0.125 * up1 * wp1;

values_per_fn[20] = 0.125 * up1 * vp1;

//fn 7

values_per_fn[21] = -0.125 * vp1 * wp1;

values_per_fn[22] = 0.125 * u * wp1;

values_per_fn[23] = 0.125 * u * vp1;

}

template<typename Scalar>

KERNEL_PREFIX void gradients_and_detJ(const Scalar* elemNodeCoords,

const Scalar* grad_vals,

Scalar& detJ)

{

/**

pt is the point at which the jacobian is to be computed.

*/

//assumptions on the lengths of input arguments:

//elemNodeCoords has length numNodesPerElem*spatialDim,

//grad_vals has length numNodesPerElem*spatialDim

const Scalar zero = 0;

Scalar J00 = zero;

Scalar J01 = zero;

Scalar J02 = zero;

Scalar J10 = zero;

Scalar J11 = zero;

Scalar J12 = zero;

Scalar J20 = zero;

Scalar J21 = zero;

Scalar J22 = zero;

size_t i_X_spatialDim = 0;

for(size_t i=0; i<numNodesPerElem; ++i) {

// size_t offset = 0;

// for(size_t gd=0; gd<spatialDim; ++gd) {

//

// Scalar gval = grad_vals[i_X_spatialDim+gd];

//

// for(size_t jd=0; jd<spatialDim; ++jd) {

// J[offset++] += gval*elemNodeCoords[i_X_spatialDim+jd];

// }

// }

//for optimization, unroll the above double-loop over spatialDim:

//(hard-coded assumption that spatialDim == 3)

J00 += grad_vals[i_X_spatialDim+0]*elemNodeCoords[i_X_spatialDim+0];

J01 += grad_vals[i_X_spatialDim+0]*elemNodeCoords[i_X_spatialDim+1];

J02 += grad_vals[i_X_spatialDim+0]*elemNodeCoords[i_X_spatialDim+2];

J10 += grad_vals[i_X_spatialDim+1]*elemNodeCoords[i_X_spatialDim+0];

J11 += grad_vals[i_X_spatialDim+1]*elemNodeCoords[i_X_spatialDim+1];

J12 += grad_vals[i_X_spatialDim+1]*elemNodeCoords[i_X_spatialDim+2];

J20 += grad_vals[i_X_spatialDim+2]*elemNodeCoords[i_X_spatialDim+0];

J21 += grad_vals[i_X_spatialDim+2]*elemNodeCoords[i_X_spatialDim+1];

J22 += grad_vals[i_X_spatialDim+2]*elemNodeCoords[i_X_spatialDim+2];

i_X_spatialDim += spatialDim;

}

Scalar term0 = J22*J11 - J21*J12;

Scalar term1 = J22*J01 - J21*J02;

Scalar term2 = J12*J01 - J11*J02;

detJ = J00*term0 - J10*term1 + J20*term2;

}

template<typename Scalar>

KERNEL_PREFIX void gradients_and_invJ_and_detJ(const Scalar* elemNodeCoords,

const Scalar* grad_vals,

Scalar* invJ,

Scalar& detJ)

{

/**

pt is the point at which the jacobian is to be computed.

*/

//assumptions on the lengths of input arguments:

//pt has length spatialDim,

//elemNodeCoords has length numNodesPerElem*spatialDim,

//grad_vals has length numNodesPerElem*spatialDim, and

//J has length spatialDim*spatialDim

const Scalar zero = 0;

//

//First we compute the jacobian J:

//

Scalar J00 = zero;

Scalar J01 = zero;

Scalar J02 = zero;

Scalar J10 = zero;

Scalar J11 = zero;

Scalar J12 = zero;

Scalar J20 = zero;

Scalar J21 = zero;

Scalar J22 = zero;

size_t i_X_spatialDim = 0;

for(size_t i=0; i<numNodesPerElem; ++i) {

// size_t offset = 0;

// for(size_t gd=0; gd<spatialDim; ++gd) {

//

// Scalar gval = grad_vals[i_X_spatialDim+gd];

//

// for(size_t jd=0; jd<spatialDim; ++jd) {

// J[offset++] += gval*elemNodeCoords[i_X_spatialDim+jd];

// }

// }

//for optimization, unroll the above double-loop over spatialDim:

//(a hard-coded assumption that spatialDim == 3)

J00 += grad_vals[i_X_spatialDim+0]*elemNodeCoords[i_X_spatialDim+0];

J01 += grad_vals[i_X_spatialDim+0]*elemNodeCoords[i_X_spatialDim+1];

J02 += grad_vals[i_X_spatialDim+0]*elemNodeCoords[i_X_spatialDim+2];

J10 += grad_vals[i_X_spatialDim+1]*elemNodeCoords[i_X_spatialDim+0];

J11 += grad_vals[i_X_spatialDim+1]*elemNodeCoords[i_X_spatialDim+1];

J12 += grad_vals[i_X_spatialDim+1]*elemNodeCoords[i_X_spatialDim+2];

J20 += grad_vals[i_X_spatialDim+2]*elemNodeCoords[i_X_spatialDim+0];

J21 += grad_vals[i_X_spatialDim+2]*elemNodeCoords[i_X_spatialDim+1];

J22 += grad_vals[i_X_spatialDim+2]*elemNodeCoords[i_X_spatialDim+2];

i_X_spatialDim += spatialDim;

}

Scalar term0 = J22*J11 - J21*J12;

Scalar term1 = J22*J01 - J21*J02;

Scalar term2 = J12*J01 - J11*J02;

detJ = J00*term0 - J10*term1 + J20*term2;

Scalar inv_detJ = 1.0/detJ;

invJ[0] = term0*inv_detJ;

invJ[1] = -term1*inv_detJ;

invJ[2] = term2*inv_detJ;

invJ[3] = -(J22*J10 - J20*J12)*inv_detJ;

invJ[4] = (J22*J00 - J20*J02)*inv_detJ;

invJ[5] = -(J12*J00 - J10*J02)*inv_detJ;

invJ[6] = (J21*J10 - J20*J11)*inv_detJ;

invJ[7] = -(J21*J00 - J20*J01)*inv_detJ;

invJ[8] = (J11*J00 - J10*J01)*inv_detJ;

}

template<typename Scalar>

KERNEL_PREFIX void diffusionMatrix_symm(const Scalar* elemNodeCoords,

const Scalar* grad_vals,

Scalar* elem_mat)

{

int len = (numNodesPerElem * (numNodesPerElem+1))/2;

const Scalar zero = 0;

miniFE::fill(elem_mat, elem_mat+len, zero);

Scalar gpts[numGaussPointsPerDim];

Scalar gwts[numGaussPointsPerDim];

gauss_pts(numGaussPointsPerDim, gpts, gwts);

const Scalar k = 1.0;

Scalar detJ = 0.0;

Scalar dpsidx[numNodesPerElem], dpsidy[numNodesPerElem], dpsidz[numNodesPerElem];

Scalar invJ[spatialDim*spatialDim];

//The following nested loop implements equations 3.4.5 and 3.4.7 on page 88

//of Reddy & Gartling, "The Finite Element Method in Heat Transfer and Fluid

//Dynamics", 2nd edition,

//to compute the element diffusion matrix for the steady conduction equation.

Scalar pt[spatialDim];

#ifdef MINIFE_DEBUG

Scalar volume = zero;

#endif

size_t gv_offset = 0;

for(size_t ig=0; ig<numGaussPointsPerDim; ++ig) {

Scalar wi = gwts[ig];

for(size_t jg=0; jg<numGaussPointsPerDim; ++jg) {

Scalar wi_wj = wi*gwts[jg];

for(size_t kg=0; kg<numGaussPointsPerDim; ++kg) {

Scalar wi_wj_wk = wi_wj*gwts[kg];

const Scalar* grad_vals_ptr = &grad_vals[gv_offset];

gv_offset += numNodesPerElem*spatialDim;

gradients_and_invJ_and_detJ(elemNodeCoords, grad_vals_ptr, invJ, detJ);

#ifdef MINIFE_DEBUG

volume += detJ;

#endif

Scalar k_detJ_wi_wj_wk = k*detJ*wi_wj_wk;

const Scalar* gv = grad_vals_ptr;

for(int i=0; i<numNodesPerElem; ++i) {

Scalar gv0 = gv[0], gv1 = gv[1], gv2 = gv[2];

dpsidx[i] = gv0 * invJ[0] +

gv1 * invJ[1] +

gv2 * invJ[2];

dpsidy[i] = gv0 * invJ[3] +

gv1 * invJ[4] +

gv2 * invJ[5];

dpsidz[i] = gv0 * invJ[6] +

gv1 * invJ[7] +

gv2 * invJ[8];

gv += spatialDim;

}

int offset = 0;

for(int m=0; m<numNodesPerElem; ++m) {

const Scalar dpsidx_m = dpsidx[m];

const Scalar dpsidy_m = dpsidy[m];

const Scalar dpsidz_m = dpsidz[m];

elem_mat[offset++] += k_detJ_wi_wj_wk *

((dpsidx_m*dpsidx_m) +

(dpsidy_m*dpsidy_m) +

(dpsidz_m*dpsidz_m));

for(int n=m+1; n<numNodesPerElem; ++n) {

elem_mat[offset++] += k_detJ_wi_wj_wk *

((dpsidx_m * dpsidx[n]) +

(dpsidy_m * dpsidy[n]) +

(dpsidz_m * dpsidz[n]));

}

}

}//for kg

}//for jg

}//for ig

//int offset = 0;

//std::cout.precision(16);

//for(int m=0; m<numNodesPerElem; ++m) {

// for(int n=m; n<numNodesPerElem; ++n) {

//std::cout<<"elem_mat["<<offset<<"] = "<<elem_mat[offset]<<";"<<std::endl;

// ++offset;

// }

//}

#ifdef MINIFE_DEBUG

// std::cout << "element volume: " << volume << std::endl;

// if (std::abs(volume - 1) > 1.e-7) {

// std::cout << "element volume is "<<volume<<", expected 1.0."<<std::endl;

// }

#endif

}

template<typename Scalar>

KERNEL_PREFIX void sourceVector(const Scalar* elemNodeCoords,

const Scalar* grad_vals,

Scalar* elem_vec)

{

int len = numNodesPerElem;

const Scalar zero = 0;

miniFE::fill(elem_vec, elem_vec+len, zero);

Scalar gpts[numGaussPointsPerDim];

Scalar gwts[numGaussPointsPerDim];

Scalar psi[numNodesPerElem];

gauss_pts(numGaussPointsPerDim, gpts, gwts);

Scalar Q = 1.0;

Scalar pt[spatialDim];

size_t gv_offset = 0;

for(size_t ig=0; ig<numGaussPointsPerDim; ++ig) {

pt[0] = gpts[ig];

Scalar wi = gwts[ig];

for(size_t jg=0; jg<numGaussPointsPerDim; ++jg) {

pt[1] = gpts[jg];

Scalar wj = gwts[jg];

for(size_t kg=0; kg<numGaussPointsPerDim; ++kg) {

pt[2] = gpts[kg];

Scalar wk = gwts[kg];

shape_fns(pt, psi);

const Scalar* grad_vals_ptr = &grad_vals[gv_offset];

gv_offset += numNodesPerElem*spatialDim;

Scalar detJ;

gradients_and_detJ(elemNodeCoords, grad_vals_ptr, detJ);

Scalar term = Q*detJ*wi*wj*wk;

for(int i=0; i<numNodesPerElem; ++i) {

elem_vec[i] += psi[i]*term;

}

}

}

}

}

}//namespace Hex8

}//namespace miniFE

#endif