operators.h

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys 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 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys 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      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#ifndef MECHSYS_TENSORS_OPERATORS_H
#define MECHSYS_TENSORS_OPERATORS_H

#ifdef HAVE_CONFIG_H
  #include "config.h"
#else
  #ifndef REAL
    #define REAL double
  #endif
#endif

#include "tensors/tensors.h"

/* Operators: {{{
 * sc: scalar
 * T2: 2nd order tensors;  6 components (Mandel representation)
 * T4: 4th order tensors; 36 components (Mandel representation)
 *
 * + add:
 *    T2 = T2 + T2
 *    T4 = T4 + T4
 * - sub:
 *    T2 = T2 - T2
 *    T4 = T4 - T4
 * * dot:
 *    sc = T2 * T2;     T2 = T2 * T4;
 *    T2 = T4 * T2;     T4 = T4 * T4;
 * & dyadic:
 *    T4 = (T2 & T2);
 * 
 }}} */

namespace Tensors
{

// 0) v = A * u
void Dot (Tensor2 const & A,
          Tensor1 const & u,
          Tensor1       & v);

// 1) sc = x * y
REAL Dot (Tensor2 const & x,
          Tensor2 const & y);

// 2) y = x * A
void Dot (Tensor2 const & x,
          Tensor4 const & A,
          Tensor2       & y);

// 3) y = A * x
void Dot (Tensor4 const & A,
          Tensor2 const & x,
          Tensor2       & y);

// 4) C = A * B
void Dot (Tensor4 const & A,
          Tensor4 const & B,
          Tensor4       & C);

// 5) A = x dyadic y
void Dyad (Tensor2 const & x,
           Tensor2 const & y,
           Tensor4       & A);

// 6) Z = a*X + b*Y
void AddScaled (REAL    const & a,
                Tensor4 const & X,
                REAL    const & b,
                Tensor4 const & Y,
                Tensor4       & Z);

// 7) A = a * (x dyadic y) + A
void Ger (REAL    const & a,
          Tensor2 const & x,
          Tensor2 const & y,
          Tensor4       & A);

// 8) D = a * (A*x) dyadic (y*B) + C
void GerX (REAL    const & a,
           Tensor4 const & A,
           Tensor2 const & x,
           Tensor2 const & y,
           Tensor4 const & B,
           Tensor4 const & C, 
           Tensor4       & D);

// 9) B = a * B
void Scale (REAL    const & a,
            Tensor4       & B);

// 10) B = a * A
void CopyScale (REAL    const & a,
                Tensor4 const & A,
                Tensor4       & B);

// 11) s = xt * A * y
REAL Reduce(Tensor2 const & x,
            Tensor4 const & A,
            Tensor2 const & y);
 // end of group TensorOperators



// 0) v = A * u
inline void Dot(Tensor2 const & A, // {{{
                Tensor1 const & u,
                Tensor1       & v)
{
    v(0) = A(0)*u(0) + A(3)*u(1) + A(5)*u(2);
    v(1) = A(3)*u(0) + A(1)*u(1) + A(4)*u(2);
    v(2) = A(5)*u(0) + A(4)*u(1) + A(2)*u(2);
} // }}}

// 1) sc = x * y
inline REAL Dot (Tensor2 const & x, Tensor2 const & y) // {{{
{
    return blitz::dot(x,y);
} // }}}

// 2) y = x * A
inline void Dot (Tensor2 const & x, Tensor4 const & A, Tensor2  & y) // {{{
{
    y(0) = x(0)*A(0,0) + x(1)*A(1,0) + x(2)*A(2,0) + x(3)*A(3,0) + x(4)*A(4,0) + x(5)*A(5,0);
    y(1) = x(0)*A(0,1) + x(1)*A(1,1) + x(2)*A(2,1) + x(3)*A(3,1) + x(4)*A(4,1) + x(5)*A(5,1);
    y(2) = x(0)*A(0,2) + x(1)*A(1,2) + x(2)*A(2,2) + x(3)*A(3,2) + x(4)*A(4,2) + x(5)*A(5,2);
    y(3) = x(0)*A(0,3) + x(1)*A(1,3) + x(2)*A(2,3) + x(3)*A(3,3) + x(4)*A(4,3) + x(5)*A(5,3);
    y(4) = x(0)*A(0,4) + x(1)*A(1,4) + x(2)*A(2,4) + x(3)*A(3,4) + x(4)*A(4,4) + x(5)*A(5,4);
    y(5) = x(0)*A(0,5) + x(1)*A(1,5) + x(2)*A(2,5) + x(3)*A(3,5) + x(4)*A(4,5) + x(5)*A(5,5);
} // }}}

// 3) y = A * x
inline void Dot (Tensor4 const & A, Tensor2 const & x, Tensor2 & y) // {{{
{
    y(0) = A(0,0)*x(0) + A(0,1)*x(1) + A(0,2)*x(2) + A(0,3)*x(3) + A(0,4)*x(4) + A(0,5)*x(5);
    y(1) = A(1,0)*x(0) + A(1,1)*x(1) + A(1,2)*x(2) + A(1,3)*x(3) + A(1,4)*x(4) + A(1,5)*x(5);
    y(2) = A(2,0)*x(0) + A(2,1)*x(1) + A(2,2)*x(2) + A(2,3)*x(3) + A(2,4)*x(4) + A(2,5)*x(5);
    y(3) = A(3,0)*x(0) + A(3,1)*x(1) + A(3,2)*x(2) + A(3,3)*x(3) + A(3,4)*x(4) + A(3,5)*x(5);
    y(4) = A(4,0)*x(0) + A(4,1)*x(1) + A(4,2)*x(2) + A(4,3)*x(3) + A(4,4)*x(4) + A(4,5)*x(5);
    y(5) = A(5,0)*x(0) + A(5,1)*x(1) + A(5,2)*x(2) + A(5,3)*x(3) + A(5,4)*x(4) + A(5,5)*x(5);
} // }}}

// 4) C = A * B
inline void Dot (Tensor4 const & A, Tensor4 const & B, Tensor4 & C) // {{{
{
    C = blitz::product(A,B);
} // }}}

// 5) A = x dyadic y
inline void Dyad (Tensor2 const & x, Tensor2 const & y, Tensor4  & A) // {{{
{
    A(0,0)=x(0)*y(0); A(0,1)=x(0)*y(1); A(0,2)=x(0)*y(2); A(0,3)=x(0)*y(3); A(0,4)=x(0)*y(4); A(0,5)=x(0)*y(5);
    A(1,0)=x(1)*y(0); A(1,1)=x(1)*y(1); A(1,2)=x(1)*y(2); A(1,3)=x(1)*y(3); A(1,4)=x(1)*y(4); A(1,5)=x(1)*y(5);
    A(2,0)=x(2)*y(0); A(2,1)=x(2)*y(1); A(2,2)=x(2)*y(2); A(2,3)=x(2)*y(3); A(2,4)=x(2)*y(4); A(2,5)=x(2)*y(5);
    A(3,0)=x(3)*y(0); A(3,1)=x(3)*y(1); A(3,2)=x(3)*y(2); A(3,3)=x(3)*y(3); A(3,4)=x(3)*y(4); A(3,5)=x(3)*y(5);
    A(4,0)=x(4)*y(0); A(4,1)=x(4)*y(1); A(4,2)=x(4)*y(2); A(4,3)=x(4)*y(3); A(4,4)=x(4)*y(4); A(4,5)=x(4)*y(5);
    A(5,0)=x(5)*y(0); A(5,1)=x(5)*y(1); A(5,2)=x(5)*y(2); A(5,3)=x(5)*y(3); A(5,4)=x(5)*y(4); A(5,5)=x(5)*y(5);
} // }}}

// 6) Z = a*X + b*Y 
inline void AddScaled (REAL const & a, Tensor4 const & X, REAL const & b, Tensor4 const & Y, Tensor4 & Z) // {{{
{
    Z(0,0)=a*X(0,0)+b*Y(0,0); Z(0,1)=a*X(0,1)+b*Y(0,1); Z(0,2)=a*X(0,2)+b*Y(0,2); Z(0,3)=a*X(0,3)+b*Y(0,3); Z(0,4)=a*X(0,4)+b*Y(0,4); Z(0,5)=a*X(0,5)+b*Y(0,5);
    Z(1,0)=a*X(1,0)+b*Y(1,0); Z(1,1)=a*X(1,1)+b*Y(1,1); Z(1,2)=a*X(1,2)+b*Y(1,2); Z(1,3)=a*X(1,3)+b*Y(1,3); Z(1,4)=a*X(1,4)+b*Y(1,4); Z(1,5)=a*X(1,5)+b*Y(1,5);
    Z(2,0)=a*X(2,0)+b*Y(2,0); Z(2,1)=a*X(2,1)+b*Y(2,1); Z(2,2)=a*X(2,2)+b*Y(2,2); Z(2,3)=a*X(2,3)+b*Y(2,3); Z(2,4)=a*X(2,4)+b*Y(2,4); Z(2,5)=a*X(2,5)+b*Y(2,5);
    Z(3,0)=a*X(3,0)+b*Y(3,0); Z(3,1)=a*X(3,1)+b*Y(3,1); Z(3,2)=a*X(3,2)+b*Y(3,2); Z(3,3)=a*X(3,3)+b*Y(3,3); Z(3,4)=a*X(3,4)+b*Y(3,4); Z(3,5)=a*X(3,5)+b*Y(3,5);
    Z(4,0)=a*X(4,0)+b*Y(4,0); Z(4,1)=a*X(4,1)+b*Y(4,1); Z(4,2)=a*X(4,2)+b*Y(4,2); Z(4,3)=a*X(4,3)+b*Y(4,3); Z(4,4)=a*X(4,4)+b*Y(4,4); Z(4,5)=a*X(4,5)+b*Y(4,5);
    Z(5,0)=a*X(5,0)+b*Y(5,0); Z(5,1)=a*X(5,1)+b*Y(5,1); Z(5,2)=a*X(5,2)+b*Y(5,2); Z(5,3)=a*X(5,3)+b*Y(5,3); Z(5,4)=a*X(5,4)+b*Y(5,4); Z(5,5)=a*X(5,5)+b*Y(5,5);
} // }}}

// 7) A = a * (x dyadic y) + A
inline void Ger (REAL const & a, Tensor2 const & x, Tensor2 const & y, Tensor4  & A) // {{{
{
    A(0,0)=a*x(0)*y(0)+A(0,0); A(0,1)=a*x(0)*y(1)+A(0,1); A(0,2)=a*x(0)*y(2)+A(0,2); A(0,3)=a*x(0)*y(3)+A(0,3); A(0,4)=a*x(0)*y(4)+A(0,4); A(0,5)=a*x(0)*y(5)+A(0,5);
    A(1,0)=a*x(1)*y(0)+A(1,0); A(1,1)=a*x(1)*y(1)+A(1,1); A(1,2)=a*x(1)*y(2)+A(1,2); A(1,3)=a*x(1)*y(3)+A(1,3); A(1,4)=a*x(1)*y(4)+A(1,4); A(1,5)=a*x(1)*y(5)+A(1,5);
    A(2,0)=a*x(2)*y(0)+A(2,0); A(2,1)=a*x(2)*y(1)+A(2,1); A(2,2)=a*x(2)*y(2)+A(2,2); A(2,3)=a*x(2)*y(3)+A(2,3); A(2,4)=a*x(2)*y(4)+A(2,4); A(2,5)=a*x(2)*y(5)+A(2,5);
    A(3,0)=a*x(3)*y(0)+A(3,0); A(3,1)=a*x(3)*y(1)+A(3,1); A(3,2)=a*x(3)*y(2)+A(3,2); A(3,3)=a*x(3)*y(3)+A(3,3); A(3,4)=a*x(3)*y(4)+A(3,4); A(3,5)=a*x(3)*y(5)+A(3,5);
    A(4,0)=a*x(4)*y(0)+A(4,0); A(4,1)=a*x(4)*y(1)+A(4,1); A(4,2)=a*x(4)*y(2)+A(4,2); A(4,3)=a*x(4)*y(3)+A(4,3); A(4,4)=a*x(4)*y(4)+A(4,4); A(4,5)=a*x(4)*y(5)+A(4,5);
    A(5,0)=a*x(5)*y(0)+A(5,0); A(5,1)=a*x(5)*y(1)+A(5,1); A(5,2)=a*x(5)*y(2)+A(5,2); A(5,3)=a*x(5)*y(3)+A(5,3); A(5,4)=a*x(5)*y(4)+A(5,4); A(5,5)=a*x(5)*y(5)+A(5,5);
} // }}}

// 8) D = a * (A*x) dyadic (y*B) + C
inline void GerX (REAL const & a, Tensor4 const & A, Tensor2 const & x, Tensor2 const & y, Tensor4 const & B, Tensor4 const & C, Tensor4 & D) // {{{
{
    //C(i,j)=a* (A(i,k)*x(k)) * (y(l)*B(l,j)) +D(i,j);

    D(0,0)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(0,0);
    D(0,1)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(0,1);
    D(0,2)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(0,2);
    D(0,3)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(0,3);
    D(0,4)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(0,4);
    D(0,5)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(0,5);
    
    D(1,0)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(1,0);
    D(1,1)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(1,1);
    D(1,2)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(1,2);
    D(1,3)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(1,3);
    D(1,4)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(1,4);
    D(1,5)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(1,5);
    
    D(2,0)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(2,0);
    D(2,1)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(2,1);
    D(2,2)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(2,2);
    D(2,3)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(2,3);
    D(2,4)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(2,4);
    D(2,5)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(2,5);
    
    D(3,0)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(3,0);
    D(3,1)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(3,1);
    D(3,2)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(3,2);
    D(3,3)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(3,3);
    D(3,4)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(3,4);
    D(3,5)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(3,5);
    
    D(4,0)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(4,0);
    D(4,1)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(4,1);
    D(4,2)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(4,2);
    D(4,3)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(4,3);
    D(4,4)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(4,4);
    D(4,5)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(4,5);
    
    D(5,0)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(5,0);
    D(5,1)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(5,1);
    D(5,2)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(5,2);
    D(5,3)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(5,3);
    D(5,4)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(5,4);
    D(5,5)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(5,5);
} // }}}

// 9) B = a * B
inline void Scale (REAL const & a, Tensor4 & B) // {{{
{
    B(0,0)=a*B(0,0); B(0,1)=a*B(0,1); B(0,2)=a*B(0,2); B(0,3)=a*B(0,3); B(0,4)=a*B(0,4); B(0,5)=a*B(0,5);
    B(1,0)=a*B(1,0); B(1,1)=a*B(1,1); B(1,2)=a*B(1,2); B(1,3)=a*B(1,3); B(1,4)=a*B(1,4); B(1,5)=a*B(1,5);
    B(2,0)=a*B(2,0); B(2,1)=a*B(2,1); B(2,2)=a*B(2,2); B(2,3)=a*B(2,3); B(2,4)=a*B(2,4); B(2,5)=a*B(2,5);
    B(3,0)=a*B(3,0); B(3,1)=a*B(3,1); B(3,2)=a*B(3,2); B(3,3)=a*B(3,3); B(3,4)=a*B(3,4); B(3,5)=a*B(3,5);
    B(4,0)=a*B(4,0); B(4,1)=a*B(4,1); B(4,2)=a*B(4,2); B(4,3)=a*B(4,3); B(4,4)=a*B(4,4); B(4,5)=a*B(4,5);
    B(5,0)=a*B(5,0); B(5,1)=a*B(5,1); B(5,2)=a*B(5,2); B(5,3)=a*B(5,3); B(5,4)=a*B(5,4); B(5,5)=a*B(5,5);
} // }}}

// 10) B = a * A
inline void CopyScale (REAL const & a, Tensor4 const & A, Tensor4 & B) // {{{
{
    B(0,0)=a*A(0,0); B(0,1)=a*A(0,1); B(0,2)=a*A(0,2); B(0,3)=a*A(0,3); B(0,4)=a*A(0,4); B(0,5)=a*A(0,5);
    B(1,0)=a*A(1,0); B(1,1)=a*A(1,1); B(1,2)=a*A(1,2); B(1,3)=a*A(1,3); B(1,4)=a*A(1,4); B(1,5)=a*A(1,5);
    B(2,0)=a*A(2,0); B(2,1)=a*A(2,1); B(2,2)=a*A(2,2); B(2,3)=a*A(2,3); B(2,4)=a*A(2,4); B(2,5)=a*A(2,5);
    B(3,0)=a*A(3,0); B(3,1)=a*A(3,1); B(3,2)=a*A(3,2); B(3,3)=a*A(3,3); B(3,4)=a*A(3,4); B(3,5)=a*A(3,5);
    B(4,0)=a*A(4,0); B(4,1)=a*A(4,1); B(4,2)=a*A(4,2); B(4,3)=a*A(4,3); B(4,4)=a*A(4,4); B(4,5)=a*A(4,5);
    B(5,0)=a*A(5,0); B(5,1)=a*A(5,1); B(5,2)=a*A(5,2); B(5,3)=a*A(5,3); B(5,4)=a*A(5,4); B(5,5)=a*A(5,5);
} // }}}

// 11) s = xt * A * y
inline REAL Reduce(Tensor2 const & x, Tensor4 const & A, Tensor2 const & y) // {{{
{
    return x(0)*(A(0,0)*y(0) + A(0,1)*y(1) + A(0,2)*y(2) + A(0,3)*y(3) + A(0,4)*y(4) + A(0,5)*y(5)) +
           x(1)*(A(1,0)*y(0) + A(1,1)*y(1) + A(1,2)*y(2) + A(1,3)*y(3) + A(1,4)*y(4) + A(1,5)*y(5)) +
           x(2)*(A(2,0)*y(0) + A(2,1)*y(1) + A(2,2)*y(2) + A(2,3)*y(3) + A(2,4)*y(4) + A(2,5)*y(5)) +
           x(3)*(A(3,0)*y(0) + A(3,1)*y(1) + A(3,2)*y(2) + A(3,3)*y(3) + A(3,4)*y(4) + A(3,5)*y(5)) +
           x(4)*(A(4,0)*y(0) + A(4,1)*y(1) + A(4,2)*y(2) + A(4,3)*y(3) + A(4,4)*y(4) + A(4,5)*y(5)) +
           x(5)*(A(5,0)*y(0) + A(5,1)*y(1) + A(5,2)*y(2) + A(5,3)*y(3) + A(5,4)*y(4) + A(5,5)*y(5));
} // }}}

}; // namespace Tensors

#endif // MECHSYS_TENSORS_OPERATORS_H

// vim:fdm=marker

---