cloasest_point_projection.cpp

Go to the documentation of this file.

/** \file small_strain_plasticity.cpp
 * \ingroup small_strain_plasticity
 * \brief Small strain plasticity example
 *
 */
 
/* This file is part of MoFEM.
 * MoFEM 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 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM 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 MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <boost/program_options.hpp>
using namespace std;
namespace po = boost::program_options;
 
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
 
#include <adolc/adolc.h>
 
#include <SmallStrainPlasticity.hpp>
#include <SmallStrainPlasticityMaterialModels.hpp>
 
using namespace boost::numeric;
 
static char help[] = "...\n"
                     "\n";
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    SmallStrainJ2Plasticity cp;
 
    cp.tAgs.resize(3);
    cp.tAgs[0] = 0;
    cp.tAgs[1] = 1;
    cp.tAgs[2] = 2;
 
    const double young_modulus = 1;
    const double poisson_ratio = 0.0;
    cp.mu = MU(young_modulus, poisson_ratio);
    cp.lambda = LAMBDA(young_modulus, poisson_ratio);
    cp.sIgma_y = 1;
    cp.H = 0.1;
    cp.K = 0;
    cp.phi = 1;
 
    cp.sTrain.resize(6, false);
    cp.sTrain.clear();
    cp.plasticStrain.resize(6, false);
    cp.plasticStrain.clear();
    cp.internalVariables.resize(7, false);
    cp.internalVariables.clear();
 
    cp.createMatAVecR();
    cp.snesCreate();
    CHKERR SNESSetFromOptions(cp.sNes);
 
    const int idx = 3;
    double alpha = 0.1;
    int N = 40;
    for (int ii = 0; ii < N; ii++) {
      cout << "Step " << ii << endl;
 
      cp.gG = ii;
 
      if (ii <= 14) {
        cp.sTrain[idx] += +alpha;
        cp.sTrain[0] += +alpha;
      } else if (ii <= 100) {
        cp.sTrain[idx] += -alpha;
        cp.sTrain[0] += -alpha;
      } else {
        cp.sTrain[idx] += +alpha;
      }
 
      cp.plasticStrain0 = cp.plasticStrain;
      cp.internalVariables0 = cp.internalVariables;
      CHKERR VecSetValue(cp.Chi, 6 + cp.internalVariables.size(), 0,
                         INSERT_VALUES);
      CHKERR VecAssemblyBegin(cp.Chi);
      CHKERR VecAssemblyEnd(cp.Chi);
      CHKERR SmallStrainPlasticityfRes(cp.sNes, cp.Chi, cp.R, &cp);
 
      double tol = 1e-8;
      bool nonlinear = false;
      if (cp.y > tol) {
        nonlinear = true;
        CHKERR cp.solveColasetProjection();
      }
 
      cout << "plot " << cp.sTrain[idx] << " " << cp.sTress[idx] << " "
           << cp.internalVariables[0] << " " << cp.plasticStrain[idx] << " "
           << cp.deltaGamma << endl;
 
      if (nonlinear) {
        CHKERR cp.consistentTangent();
        cerr << endl;
        cerr << "Cep " << cp.Cep << endl;
        MatrixDouble e = cp.Cep;
        VectorDouble strain0 = cp.sTrain;
        const double eps = 1e-6;
        for (int ii = 0; ii < 6; ii++) {
          cp.sTrain = strain0;
          cp.sTrain[ii] -= eps;
          CHKERR cp.solveColasetProjection();
          for (int dd = 0; dd < 6; dd++) {
            e(dd, ii) += cp.sTress[dd] / (2 * eps);
          }
          cp.sTrain = strain0;
          cp.sTrain[ii] += eps;
          CHKERR cp.solveColasetProjection();
          for (int dd = 0; dd < 6; dd++) {
            e(dd, ii) -= cp.sTress[dd] / (2 * eps);
          }
        }
        cp.sTrain = strain0;
        cerr << "e   " << e << endl;
      }
 
      // std::string wait;
      // std::cin >> wait;
    }
 
    cp.destroyMatAVecR();
    cp.snesDestroy();
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
 
  return 0;
}