tfuns.cpp

/*************************************************************************************
 * 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                                           *
 *************************************************************************************/

#include <iostream>

#include "tensors/tensors.h"
#include "tensors/functions.h"

using namespace std;

int main(int argc, char **argv)
{
    //Tensors::Tensor2 a("1 2 3 4 5 6");
    //Tensors::Tensor2 a; a=2,2,2,0,0,0.0001;
    Tensors::Tensor2 a; a(0)=10; a(1)=20; a(2)=30; a(3)=4*sqrt(2.0); a(4)=6*sqrt(2.0);  a(5)=5*sqrt(2.0);
    Tensors::Tensor2 b; b(0)=20; b(1)=40; b(2)=60; b(3)=8*sqrt(2.0); b(4)=12*sqrt(2.0); b(5)=10*sqrt(2.0); // with this tensor, the matrix product a*b will be symmetric
    //Tensors::Tensor2 b(" 6  5  4  3  2  1"); // with this tensor, the matrix product a*b will NOT be symmetric

    cout << "a = " << a << endl;
    cout << "b = " << b << endl;

    Tensors::Tensor4 T;
    for (int i=0; i<6; ++i)
    for (int j=0; j<6; ++j)
        T(i,j) = 10*(i+1) + (j+1);

    Tensors::Tensor2 ab;
    Tensors::Tensor2 inv_a;
    REAL             only_ev_a[3];
    REAL             ev_a[3];
    Tensors::Tensor2 ep_a[3];
    Tensors::Tensor2 sqrt_a;

    Tensors::Mult(a,b,ab);

    if (!Tensors::Inv(a,inv_a))
    {
        cout << "ERROR: Inv(a) failed\n";
    }
    if (!Tensors::Eigenvals(a,only_ev_a))
    {
        cout << "ERROR: Eigenvals(a) failed\n";
    }
    if (!Tensors::Eigenvp(a,ev_a,ep_a))
    {
        cout << "ERROR: Eigenvp(a) failed\n";
    }
    if (!Tensors::Sqrt(a,sqrt_a))
    {
        cout << "ERROR: sqrt(a) failed\n";
    }

    cout << "ab =     \n" << ab << endl;
    cout << "Norm(a) = " << Tensors::Norm(a) << endl;
    cout << "inv_a =  \n" << inv_a << endl;
    cout << "only_ev_a =   \n"; for (int i=0; i<3; ++i) cout << only_ev_a[i] << endl;
    cout << "ev_a =   \n"; for (int i=0; i<3; ++i) cout << ev_a[i] << endl;
    cout << "ep_a =   \n"; for (int i=0; i<3; ++i) cout << ep_a[i] << endl;
    cout << "sqrt_a = \n" << sqrt_a << endl;

    REAL Ia[3];
    REAL Ev,Ed,p,q,tn,ts,P,Q;

    Tensors::Tensor2 I; I(0)=1.0; I(1)=1.0; I(2)=1.0; I(3)=0.0; I(4)=0.0; I(5)=0.0; 
    //Tensors::Tensor2 isoA; isoA = 
    //Tensors::Tensor2 devA;
    Tensors::CharInvs(a,Ia);
    Tensors::Strain_Ev_Ed(a,Ev,Ed);
    Tensors::Stress_p_q(a,p,q);
    Tensors::Stress_tn_ts(a,tn,ts);
    Tensors::Stress_P_Q(a,P,Q);

    cout << "a = " << a << endl;
    cout << "Ia[0]=" << Ia[0] << ", Ia[1]=" << Ia[1] << ", Ia[2]=" << Ia[2] << endl;
    cout << "Ev=" << Ev << ", Ed=" << Ed << ", p=" << p << ", q=" << q << endl;
    cout << "tn=" << tn << ", ts=" << ts << endl;
    cout << "P=" << P << ", Q=" << Q << endl;

    cout << "I=" << Tensors::I << endl;
    cout << "IIsym=\n" << Tensors::IIsym << endl;
    cout << "IdyI=\n"  << Tensors::IdyI << endl;

    REAL sin3th;
    Tensors::Tensor2 S;
    Tensors::Stress_p_q_S_sin3th(a, p,q,S,sin3th);
    cout << "p=" << p << ", q=" << q << ", S=" << S << ", sin3th=" << sin3th << endl;

    Tensors::Stress_P_Q_S_sin3th(a, P,Q,S,sin3th);
    cout << "P=" << P << ", Q=" << Q << ", sin3th=" << sin3th << endl;

    return 0;
}

---