cloasest_point_projection.cpp File Reference

#include <MoFEM.hpp>
#include <boost/program_options.hpp>
#include <ADOLCPlasticity.hpp>
#include <ADOLCPlasticityMaterialModels.hpp>

Go to the source code of this file.

Functions
int	main (int argc, char *argv[])

Variables
static char	help []

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 37 of file cloasest_point_projection.cpp.

                                 {
 
  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;
 
    auto cp_ptr = createMaterial<ParaboloidalPlasticity>();
 
    VectorDouble active_variables(12 + cp_ptr->nbInternalVariables);
    VectorDouble grad_w(12 + cp_ptr->nbInternalVariables);
 
    auto tmp_active = getVectorAdaptor(&(active_variables[0]),
                                       12 + cp_ptr->nbInternalVariables);
    cp_ptr->activeVariablesW.swap(tmp_active);
    auto tmp_gradient =
        getVectorAdaptor(&(grad_w[0]), 12 + cp_ptr->nbInternalVariables);
    cp_ptr->gradientW.swap(tmp_gradient);
 
    auto strain = cp_ptr->getTotalStrain();
    auto internal_variables = cp_ptr->getInternalVariables();
    auto stress = cp_ptr->getStress();
    auto plastic_strain = cp_ptr->getPlasticStrain();
 
    cp_ptr->createMatAVecR();
    cp_ptr->snesCreate();
 
    const int idx = 3;
    double alpha = 0.1;
    int N = 40;
    for (int ii = 0; ii < N; ii++) {
      cout << "Step " << ii << endl;
 
      if (ii <= 14) {
        strain[idx] += +alpha;
        strain[0] += +alpha;
      } else if (ii <= 100) {
        strain[idx] += -alpha;
        strain[0] += -alpha;
      } else {
        strain[idx] += +alpha;
      }
 
      cp_ptr->plasticStrain0 = plastic_strain;
      cp_ptr->internalVariables0 = internal_variables;
      CHKERR VecSetValue(cp_ptr->Chi, 6 + cp_ptr->nbInternalVariables, 0,
                         INSERT_VALUES);
      CHKERR VecAssemblyBegin(cp_ptr->Chi);
      CHKERR VecAssemblyEnd(cp_ptr->Chi);
      CHKERR ADOLCPlasticityRes(cp_ptr->sNes, cp_ptr->Chi, cp_ptr->R,
                                       cp_ptr.get());
 
      double tol = 1e-8;
      bool nonlinear = false;
      if (cp_ptr->y > tol) {
        nonlinear = true;
        CHKERR cp_ptr->solveClosetProjection();
      }
 
      cout << "plot " << strain[idx] << " " << stress[idx] << " "
           << internal_variables[0] << " " << plastic_strain[idx] << " "
           << cp_ptr->deltaGamma << endl;
 
      if (nonlinear) {
        CHKERR cp_ptr->consistentTangent();
        cerr << endl;
        cerr << "Cep " << cp_ptr->Cep << endl;
        MatrixDouble e = cp_ptr->Cep;
        VectorDouble strain0 = strain;
        const double eps = 1e-6;
        for (int ii = 0; ii < 6; ii++) {
          strain = strain0;
          strain[ii] -= eps;
          CHKERR cp_ptr->solveClosetProjection();
          for (int dd = 0; dd < 6; dd++) {
            e(dd, ii) += stress[dd] / (2 * eps);
          }
          strain = strain0;
          strain[ii] += eps;
          CHKERR cp_ptr->solveClosetProjection();
          for (int dd = 0; dd < 6; dd++) {
            e(dd, ii) -= stress[dd] / (2 * eps);
          }
        }
        strain = strain0;
        cerr << "e   " << e << endl;
      }
 
      // std::string wait;
      // std::cin >> wait;
    }
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

help

char help[]

static

Initial value:

= "...\n"
                     "\n"

Definition at line 34 of file cloasest_point_projection.cpp.