tensortest.cc

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=8 sw=2 et sts=2:
#include <config.h>

#include <dune/common/parallel/mpihelper.hh>

#include <dune/fufem/symmetrictensor.hh>
#include <dune/fufem/mechanics/isotropictensor.hh>
#include <dune/fufem/mechanics/cubictensor.hh>
#include <dune/fufem/mechanics/tetratensor.hh>

#include <cmath>
#include <iostream>

int main (int argc, char *argv[])
{
  bool passed = true;
  Dune::MPIHelper::instance(argc, argv);

  SymmetricTensor<3> e1;
  e1(0,0) = 1;
  e1(1,1) = 2;
  e1(2,2) = 3;
  e1(1,0) = 4;
  e1(2,0) = 5;
  e1(2,1) = 6;
  passed = passed and std::abs(e1.trace() - 6.0) < 1e-10;

  SymmetricTensor<3> e2;
  e2 = e1;

  passed = passed and   e1 == e2;
  passed = passed and !(e1 != e2);

  passed = passed and std::abs(e2(0,0) - 1) < 1e-10;
  passed = passed and std::abs(e2(1,1) - 2) < 1e-10;
  passed = passed and std::abs(e2(2,2) - 3) < 1e-10;
  passed = passed and std::abs(e2(0,1) - 4) < 1e-10;
  passed = passed and std::abs(e2(0,2) - 5) < 1e-10;
  passed = passed and std::abs(e2(1,2) - 6) < 1e-10;

  e2.addToDiag(1);
  passed = passed and std::abs(e2(0,0) - 2) < 1e-10;
  passed = passed and std::abs(e2(1,1) - 3) < 1e-10;
  passed = passed and std::abs(e2(2,2) - 4) < 1e-10;
  passed = passed and std::abs(e2(0,1) - 4) < 1e-10;
  passed = passed and std::abs(e2(0,2) - 5) < 1e-10;
  passed = passed and std::abs(e2(1,2) - 6) < 1e-10;

  e2.setDiag(1);
  e2(1,1) += 1.0;
  e2(2,2) += 2.0;
  e2(0,1) += 4.0;
  e2(0,2) += 5.0;
  e2(1,2) += 6.0;
  e2 += e1;
  e2 -= e1;

  e2 *= 2;
  passed = passed and std::abs(e1 * e2 - 2 * (e1 * e1)) < 1e-10;

  passed = passed and !(e1 == e2);
  passed = passed and   e1 != e2;

  e2 /= 2;
  passed = passed and std::abs(e1 * e1 - e1 * e2) < 1e-10;

  IsotropicTensor<3> const C(.5,.25);
  SymmetricTensor<3> h;
  C.mv(e1, h);

  std::cout << e1 << std::endl;
  std::cout << h << std::endl;
  std::cout << e1 * h << std::endl;
  passed = passed and std::abs(e1 * h - 74.4) < 1e-10; // e : C : e

  SymmetricTensor<3> const e3(e2);
  std::cout << e3 << std::endl;
  std::cout << e3(1,1) << std::endl;

  TetraTensor<2> const tt2(1,2,3,4);
  TetraTensor<3> const tt3(1,2,3,4,5,6);
  CubicTensor<2> const ct2(1,2,3);
  CubicTensor<3> const ct3(1,2,3);

  return passed ? 0 : 1;
}