elplast solver and mesh

miniapps/CMakeLists.txt

@@ -31,4 +31,5 @@ add_subdirectory(shifted)
 add_subdirectory(mtop)
 add_subdirectory(autodiff)
 add_subdirectory(parelag)
+add_subdirectory(elplast)
 add_subdirectory(hooke)

miniapps/elplast/CMakeLists.txt

@@ -0,0 +1,49 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+if (MFEM_USE_MPI)
+  list(APPEND DIST_COMMON_SOURCES
+	  coefficients.cpp
+	  elplast.cpp
+      )
+  list(APPEND DIST_COMMON_HEADERS
+	  coefficients.hpp
+	  elplast.hpp
+      )
+
+  convert_filenames_to_full_paths(DIST_COMMON_SOURCES)
+  convert_filenames_to_full_paths(DIST_COMMON_HEADERS)
+
+  set(DIST_COMMON_FILES
+      EXTRA_SOURCES ${DIST_COMMON_SOURCES}
+      EXTRA_HEADERS ${DIST_COMMON_HEADERS})
+
+  add_mfem_miniapp(adapt
+    MAIN adaptive_el.cpp
+    ${DIST_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(test_elplinteg
+    MAIN test_elinteg.cpp
+    ${DIST_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(elpl_elem
+    MAIN elpl_test.cpp
+    ${DIST_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(holep
+    MAIN holep.cpp
+    ${DIST_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+
+endif ()

miniapps/elplast/adaptive_el.cpp

@@ -0,0 +1,401 @@
+//                       MFEM Example 21 - Parallel Version
+//
+// Compile with: make ex21p
+//
+// Sample runs:  mpirun -np 4 ex21p
+//               mpirun -np 4 ex21p -o 3
+//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh
+//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh -o 2
+//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh -o 2
+//
+// Description:  This is a version of Example 2p with a simple adaptive mesh
+//               refinement loop. The problem being solved is again the linear
+//               elasticity describing a multi-material cantilever beam.
+//               The problem is solved on a sequence of meshes which
+//               are locally refined in a conforming (triangles, tetrahedrons)
+//               or non-conforming (quadrilaterals, hexahedra) manner according
+//               to a simple ZZ error estimator.
+//
+//               The example demonstrates MFEM's capability to work with both
+//               conforming and nonconforming refinements, in 2D and 3D, on
+//               linear and curved meshes. Interpolation of functions from
+//               coarse to fine meshes, as well as persistent GLVis
+//               visualization are also illustrated.
+//
+//               We recommend viewing Examples 2p and 6p before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "coefficients.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/beam-tri.mesh";
+   int serial_ref_levels = 0;
+   int order = 2;
+   //bool static_cond = false;
+   bool static_cond = true;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&serial_ref_levels, "-rs", "--refine-serial",
+                  "Number of uniform serial refinements (before parallel"
+                  " partitioning)");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, and hexahedral meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
+
+   // 3. Refine the mesh before parallel partitioning. Since a NURBS mesh can
+   //    currently only be refined uniformly, we need to convert it to a
+   //    piecewise-polynomial curved mesh. First we refine the NURBS mesh a bit
+   //    more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext && serial_ref_levels == 0)
+   {
+      serial_ref_levels = 2;
+   }
+   for (int i = 0; i < serial_ref_levels; i++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(2);
+   }
+   mesh.EnsureNCMesh();
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // 4. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec, dim);
+
+   // 5. As in Example 2, we set up the linear form b(.) which corresponds to
+   //    the right-hand side of the FEM linear system. In this case, b_i equals
+   //    the boundary integral of f*phi_i where f represents a "pull down"
+   //    force on the Neumann part of the boundary and phi_i are the basis
+   //    functions in the finite element fespace. The force is defined by the
+   //    VectorArrayCoefficient object f, which is a vector of Coefficient
+   //    objects. The fact that f is non-zero on boundary attribute 2 is
+   //    indicated by the use of piece-wise constants coefficient for its last
+   //    component. We don't assemble the discrete problem yet, this will be
+   //    done in the main loop.
+   VectorArrayCoefficient f(dim);
+   f.Set(0, new ConstantCoefficient(0.0));
+   f.Set(1, new ConstantCoefficient(1.0));
+   if(dim==3){
+       f.Set(2, new ConstantCoefficient(0.0));
+   }
+
+   ParLinearForm b(&fespace);
+
+   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f));
+
+
+   // 6. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the linear elasticity integrator with piece-wise
+   //    constants coefficient lambda and mu.
+   /*
+   Vector lambda(pmesh.attributes.Max());
+   lambda = 1.0;
+   lambda(0) = lambda(1)*50;
+   PWConstCoefficient lambda_func(lambda);
+   Vector mu(pmesh.attributes.Max());
+   mu = 1.0;
+   mu(0) = mu(1)*50;
+   PWConstCoefficient mu_func(mu);
+   */
+
+   //ConstantCoefficient lambda_func(1.0);
+   ConstantCoefficient lambda_func(0.576923E+7);
+   //ConstantCoefficient mu_func(1.0);
+   ConstantCoefficient mu_func(0.384615E+7);
+
+   ParBilinearForm a(&fespace);
+   BilinearFormIntegrator *integ =
+      new ElasticityIntegrator(lambda_func,mu_func);
+   a.AddDomainIntegrator(integ);
+   if (static_cond) { a.EnableStaticCondensation(); }
+
+   // 7. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   Vector zero_vec(dim);
+   zero_vec = 0.0;
+   VectorConstantCoefficient zero_vec_coeff(zero_vec);
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking only
+   //    boundary attribute 1 from the mesh as essential and converting it to a
+   //    list of true dofs.  The conversion to true dofs will be done in the
+   //    main loop.
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 0;
+   ess_bdr[0] = 1;
+
+   // 9. GLVis visualization.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+
+   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     that uses the ComputeElementFlux method of the ElasticityIntegrator to
+   //     recover a smoothed flux (stress) that is subtracted from the element
+   //     flux to get an error indicator. We need to supply the space for the
+   //     smoothed flux: an (H1)^tdim (i.e., vector-valued) space is used here.
+   //     Here, tdim represents the number of components for a symmetric (dim x
+   //     dim) tensor.
+   const int tdim = dim*(dim+1)/2;
+   L2_FECollection flux_fec(order, dim);
+   ParFiniteElementSpace flux_fespace(&pmesh, &flux_fec, tdim);
+   ParFiniteElementSpace smooth_flux_fespace(&pmesh, &fec, tdim);
+   L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace,
+                                       smooth_flux_fespace);
+
+
+   //strain fields
+   L2_FECollection strain_fec(order-1,dim);
+   ParFiniteElementSpace strain_fespace(&pmesh,&strain_fec, dim*(dim+1)/2);
+   ParGridFunction estrains; estrains.SetSpace(&strain_fespace);
+   ParGridFunction sstrains; sstrains.SetSpace(&strain_fespace);
+
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaView", &pmesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.SetDataFormat(VTKFormat::ASCII);
+   dacol.RegisterField("disp", &x);
+   dacol.RegisterField("estr",&estrains);
+   dacol.RegisterField("sstr",&sstrains);
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+   // 11. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.7);
+
+   // 12. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   const int max_dofs = 100000;
+   const int max_amr_itr = 20;
+   for (int it = 0; it <= max_amr_itr; it++)
+   {
+      HYPRE_BigInt global_dofs = fespace.GlobalTrueVSize();
+      if (myid == 0)
+      {
+         cout << "\nAMR iteration " << it << endl;
+         cout << "Number of unknowns: " << global_dofs << endl;
+      }
+
+      // 13. Assemble the stiffness matrix and the right-hand side.
+      a.Assemble();
+      b.Assemble();
+
+      // 14. Set Dirichlet boundary values in the GridFunction x.
+      //     Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_tdof_list;
+      x.ProjectBdrCoefficient(zero_vec_coeff, ess_bdr);
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      // 15. Create the linear system: eliminate boundary conditions, constrain
+      //     hanging nodes and possibly apply other transformations. The system
+      //     will be solved for true (unconstrained) DOFs only.
+
+      HypreParMatrix A;
+      Vector B, X;
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+
+      // 16. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
+      //     preconditioner from hypre.
+      HypreBoomerAMG amg;
+      amg.SetPrintLevel(0);
+      // amg.SetSystemsOptions(dim); // optional
+      CGSolver pcg(A.GetComm());
+      pcg.SetPreconditioner(amg);
+      pcg.SetOperator(A);
+      pcg.SetRelTol(1e-6);
+      pcg.SetMaxIter(500);
+      pcg.SetPrintLevel(3); // print the first and the last iterations only
+      pcg.Mult(B, X);
+
+      // 17. After solving the linear system, reconstruct the solution as a
+      //     finite element GridFunction. Constrained nodes are interpolated
+      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
+      a.RecoverFEMSolution(X, b, x);
+
+      //compute the starins
+      EngStrainCoefficient estr(x);
+      StrainCoefficient sstr(x);
+
+      estrains.ProjectCoefficient(estr);
+      sstrains.ProjectCoefficient(sstr);
+
+      //ParaView output
+      dacol.SetTime(double(it+1));
+      dacol.SetCycle(it+1);
+      dacol.Save();
+
+
+      // 18. Send solution by socket to the GLVis server.
+      if (visualization && it == 0)
+      {
+         sol_sock.open(vishost, visport);
+         sol_sock.precision(8);
+      }
+      if (visualization && sol_sock.good())
+      {
+         GridFunction nodes(&fespace), *nodes_p = &nodes;
+         pmesh.GetNodes(nodes);
+         nodes += x;
+         int own_nodes = 0;
+         pmesh.SwapNodes(nodes_p, own_nodes);
+         x.Neg(); // visualize the backward displacement
+         sol_sock << "parallel " << num_procs << ' ' << myid << '\n';
+         sol_sock << "solution\n" << pmesh << x << flush;
+         x.Neg();
+         pmesh.SwapNodes(nodes_p, own_nodes);
+         if (it == 0)
+         {
+            sol_sock << "keys '" << ((dim == 2) ? "Rjl" : "") << "m'" << endl;
+         }
+         sol_sock << "window_title 'AMR iteration: " << it << "'\n"
+                  << "pause" << endl;
+         if (myid == 0)
+         {
+            cout << "Visualization paused. "
+                 "Press <space> in the GLVis window to continue." << endl;
+         }
+      }
+
+      if (global_dofs > max_dofs)
+      {
+         if (myid == 0)
+         {
+            cout << "Reached the maximum number of dofs. Stop." << endl;
+         }
+         break;
+      }
+
+      // 19. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+      refiner.Apply(pmesh);
+      if (refiner.Stop())
+      {
+         if (myid == 0)
+         {
+            cout << "Stopping criterion satisfied. Stop." << endl;
+         }
+         break;
+      }
+
+      // 20. Update the space to reflect the new state of the mesh. Also,
+      //     interpolate the solution x so that it lies in the new space but
+      //     represents the same function. This saves solver iterations later
+      //     since we'll have a good initial guess of x in the next step.
+      //     Internally, FiniteElementSpace::Update() calculates an
+      //     interpolation matrix which is then used by GridFunction::Update().
+      fespace.Update();
+      x.Update();
+
+      strain_fespace.Update();
+      estrains.Update();
+      sstrains.Update();
+
+      // 21. Load balance the mesh, and update the space and solution. Currently
+      //     available only for nonconforming meshes.
+      if (pmesh.Nonconforming())
+      {
+         pmesh.Rebalance();
+
+         // Update the space and the GridFunction. This time the update matrix
+         // redistributes the GridFunction among the processors.
+         fespace.Update();
+         x.Update();
+      }
+
+      // 21. Inform also the bilinear and linear forms that the space has
+      //     changed.
+      a.Update();
+      b.Update();
+   }
+
+   dacol.SetTime(double(max_amr_itr));
+   dacol.SetCycle(max_amr_itr);
+   dacol.Save();
+
+   {
+      ostringstream mref_name, mesh_name, sol_name;
+      mref_name << "ex21p_reference_mesh." << setfill('0') << setw(6) << myid;
+      mesh_name << "ex21p_deformed_mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "ex21p_displacement." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ref_out(mref_name.str().c_str());
+      mesh_ref_out.precision(16);
+      pmesh.Print(mesh_ref_out);
+
+      ofstream mesh_out(mesh_name.str().c_str());
+      mesh_out.precision(16);
+      GridFunction nodes(&fespace), *nodes_p = &nodes;
+      pmesh.GetNodes(nodes);
+      nodes += x;
+      int own_nodes = 0;
+      pmesh.SwapNodes(nodes_p, own_nodes);
+      pmesh.Print(mesh_out);
+      pmesh.SwapNodes(nodes_p, own_nodes);
+
+      ofstream x_out(sol_name.str().c_str());
+      x_out.precision(16);
+      x.Save(x_out);
+   }
+
+   return 0;
+}

miniapps/elplast/coefficients.hpp

@@ -0,0 +1,733 @@
+#ifndef COEFFICIENTS_HPP
+#define COEFFICIENTS_HPP
+
+
+#include "mfem.hpp"
+#include "linalg/dual.hpp"
+
+namespace mfem {
+
+
+class StrainCoefficient:public VectorCoefficient{
+public:
+    StrainCoefficient(GridFunction& disp_):VectorCoefficient(disp_.VectorDim()*(disp_.VectorDim()+1)/2)
+    {
+        disp=&disp_;
+        g.SetSize(disp->VectorDim());
+        e.SetSize(disp->VectorDim());
+    }
+
+    virtual
+    ~StrainCoefficient(){}
+
+    virtual
+    void Eval(Vector &V, ElementTransformation &T,
+              const IntegrationPoint &ip)
+    {
+        T.SetIntPoint(&ip);
+        disp->GetVectorGradient(T,g);
+        for(int i=0;i<disp->VectorDim();i++){
+        for(int j=i;j<disp->VectorDim();j++){
+            e(i,j)=0.5*(g(i,j)+g(j,i));
+        }}
+
+        V.SetSize(GetVDim());
+
+        if(disp->VectorDim()==3){
+            V(0)=e(0,0);
+            V(1)=e(1,1);
+            V(2)=e(2,2);
+            V(3)=e(1,2);
+            V(4)=e(0,2);
+            V(5)=e(0,1);
+        }else{
+            V(0)=e(0,0);
+            V(1)=e(1,1);
+            V(2)=e(0,1);
+        }
+    }
+
+private:
+    GridFunction* disp;
+    DenseMatrix g;
+    DenseMatrix e;
+
+};
+
+
+class EngStrainCoefficient: public mfem::VectorCoefficient{
+public:
+
+    EngStrainCoefficient(GridFunction& disp_):VectorCoefficient(disp_.VectorDim()*(disp_.VectorDim()+1)/2)
+    {
+        disp=&disp_;
+        g.SetSize(disp->VectorDim());
+        e.SetSize(disp->VectorDim());
+    }
+
+    virtual
+    ~EngStrainCoefficient(){}
+
+    virtual
+    void Eval(Vector &V, ElementTransformation &T,
+              const IntegrationPoint &ip)
+    {
+        T.SetIntPoint(&ip);
+        disp->GetVectorGradient(T,g);
+        for(int i=0;i<disp->VectorDim();i++){
+        for(int j=i;j<disp->VectorDim();j++){
+            e(i,j)=0.5*(g(i,j)+g(j,i));
+        }}
+
+        V.SetSize(GetVDim());
+
+        if(disp->VectorDim()==3){
+            V(0)=e(0,0);
+            V(1)=e(1,1);
+            V(2)=e(2,2);
+            V(3)=2.0*e(1,2);
+            V(4)=2.0*e(0,2);
+            V(5)=2.0*e(0,1);
+        }else{
+            V(0)=e(0,0);
+            V(1)=e(1,1);
+            V(2)=2.0*e(0,1);
+        }
+    }
+
+
+private:
+    GridFunction* disp;
+    DenseMatrix g;
+    DenseMatrix e;
+
+
+};
+
+template<typename ElastMaterial>
+class StressCoefficient:public VectorCoefficient
+{
+public:
+    StressCoefficient(ElastMaterial& mat_, MatrixCoefficient& str_):
+        VectorCoefficient(str_.GetWidth()*(str_.GetWidth()+1)/2)
+    {
+        mat=&mat_;
+        str=&str_;
+        ss.SetSize(str_.GetWidth());
+
+    }
+
+    virtual
+    void Eval(Vector &V, ElementTransformation &T,
+              const IntegrationPoint &ip)
+    {
+        str->Eval(ss,T,ip);
+        mat->EvalStress(st,ss,T,ip);
+        if(str->GetWidth()==3){
+            V(0)=st(0,0);
+            V(1)=st(1,1);
+            V(2)=st(2,2);
+            V(3)=st(1,2);
+            V(4)=st(0,2);
+            V(5)=st(0,1);
+        }else{
+            V(0)=st(0,0);
+            V(1)=st(1,1);
+            V(2)=st(0,1);
+        }
+    }
+
+private:
+    ElastMaterial* mat;
+    MatrixCoefficient* str;
+
+    DenseMatrix ss; //strain
+    DenseMatrix st; //stress
+};
+
+template<typename ElastMaterial>
+class VonMisesStressCoefficient:public Coefficient
+{
+public:
+    VonMisesStressCoefficient(ElastMaterial& mat_, MatrixCoefficient& str_){
+        mat=&mat_;
+        str=&str_;//stress
+        ss.SetSize(str->GetWidth());
+        st.SetSize(str->GetWidth());
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+                 const IntegrationPoint &ip)
+    {
+        str->Eval(ss,T,ip);
+        mat->EvalStress(st,ss,T,ip);
+
+        double res=0.0;
+
+        if(str->GetWidth()==3){
+            res=(st(0,0)-st(1,1))*(st(0,0)-st(1,1));
+            res=res+(st(1,1)-st(2,2))*(st(1,1)-st(2,2));
+            res=res+(st(0,0)-st(2,2))*(st(0,0)-st(2,2));
+            res=res+6*st(1,2)*st(1,2);
+            res=res+6*st(0,2)*st(0,2);
+            res=res+6*st(0,1)*st(0,1);
+        }else{
+            res=(st(0,0)-st(1,1))*(st(0,0)-st(1,1));
+            res=res+st(1,1)*st(1,1);
+            res=res+st(0,0)*st(0,0);
+            res=res+6*st(0,1)*st(0,1);
+        }
+
+        return sqrt(res);
+
+    }
+
+private:
+   ElastMaterial* mat;
+   MatrixCoefficient* str;
+   DenseMatrix st; //stress
+   DenseMatrix ss; //strain
+
+};
+
+class IsoElastMat
+{
+public:
+    IsoElastMat()
+    {
+        E=1.0;
+        nu=0.2;
+    }
+
+    void SetLameParam(double lam, double mu)
+    {
+        E=mu*(3.0*lam+2.0*mu)/(lam+mu);
+        nu=lam/(2.0*(lam+mu));
+    }
+
+    void SetElastParam(double E_,double nu_)
+    {
+        E=E_;
+        nu=nu_;
+    }
+
+    void SetE(double E_){ E=E_;}
+
+    void SetPoisson(double nu_){ nu=nu_;}
+
+
+    void EvalStress(DenseMatrix st, DenseMatrix ss,
+                    ElementTransformation &T,
+                    const IntegrationPoint &ip)
+    {
+        double cc=E/((1.0+nu)*(1.0-2.0*nu));
+        if(ss.Width()==3){
+            st(0,0)=cc*((1.0-nu)*ss(0,0)+nu*ss(1,1)+nu*ss(2,2));
+            st(1,1)=cc*(nu*ss(0,0)+(1.0-nu)*ss(1,1)+nu*ss(2,2));
+            st(2,2)=cc*(nu*ss(0,0)+nu*ss(1,1)+(1.0-nu)*ss(2,2));
+            st(1,2)=cc*(1.0-2.0*nu)*ss(1,2); st(2,1)=cc*(1.0-2.0*nu)*ss(2,1);
+            st(0,2)=cc*(1.0-2.0*nu)*ss(0,2); st(2,0)=cc*(1.0-2.0*nu)*ss(2,0);
+            st(0,1)=cc*(1.0-2.0*nu)*ss(0,1); st(1,0)=cc*(1.0-2.0*nu)*ss(1,0);
+        }else{
+            st(0,0)=cc*((1.0-nu)*ss(0,0)+nu*ss(1,1));
+            st(1,1)=cc*(nu*ss(0,0)+(1.0-nu)*ss(1,1));
+            //st(2,2)=cc*(nu*ss(0,0)+nu*ss(1,1));
+            st(0,1)=cc*(1.0-2.0*nu)*ss(0,1); st(1,0)=cc*(1.0-2.0*nu)*ss(1,0);
+        }
+    }
+
+    template<typename dtype>
+    void EvalStress(dtype* st, dtype* ss,
+                    ElementTransformation &T,
+                    const IntegrationPoint &ip)
+    {
+        double cc=E/((1.0+nu)*(1.0-2.0*nu));
+        st[0+3*0]=cc*((1.0-nu)*ss[0+3*0]+nu*ss[1+3*1]+nu*ss[2+3*2]);
+        st[1+3*1]=cc*(nu*ss[0+3*0]+(1.0-nu)*ss[1+3*1]+nu*ss[2+3*2]);
+        st[2+3*2]=cc*(nu*ss[0+3*0]+nu*ss[1+3*1]+(1.0-nu)*ss[2+3*2]);
+        st[1+3*2]=cc*(1.0-2.0*nu)*ss[1+3*2]; st[2+3*1]=cc*(1.0-2.0*nu)*ss[2+3*1];
+        st[0+3*2]=cc*(1.0-2.0*nu)*ss[0+3*2]; st[2+3*0]=cc*(1.0-2.0*nu)*ss[2+3*0];
+        st[0+3*1]=cc*(1.0-2.0*nu)*ss[0+3*1]; st[1+3*0]=cc*(1.0-2.0*nu)*ss[1+3*0];
+
+    }
+
+    // mat[9x9] ss[9]
+    template<typename dtype>
+    void EvalGrad(dtype* mat, dtype* ss,
+                  ElementTransformation &T,
+                  const IntegrationPoint &ip)
+    {
+        double cc=E/((1.0+nu)*(1.0-2.0*nu));
+        for(int i=0;i<81;i++){ mat[i]=dtype(0.0);}
+
+        mat[0*9 + 0]=cc*(1.0-nu);
+        mat[0*9 + 4]=cc*nu;
+        mat[0*9 + 8]=cc*nu;
+
+        mat[4*9 + 0]=cc*nu;
+        mat[4*9 + 4]=cc*(1.0-nu);
+        mat[4*9 + 8]=cc*nu;
+
+        mat[8*9+  0]=cc*nu;
+        mat[8*9+  4]=cc*nu;
+        mat[8*9+  8]=cc*(1.0-nu);
+
+        mat[1+3*2+(1+3*2)*9]=cc*(1.0-2.0*nu);
+        mat[0+3*2+(0+3*2)*9]=cc*(1.0-2.0*nu);
+        mat[0+3*1+(0+3*1)*9]=cc*(1.0-2.0*nu);
+
+        mat[2+3*1+(2+3*1)*9]=cc*(1.0-2.0*nu);
+        mat[2+3*0+(2+3*0)*9]=cc*(1.0-2.0*nu);
+        mat[1+3*0+(1+3*0)*9]=cc*(1.0-2.0*nu);
+    }
+
+
+    void EvalGrad(DenseMatrix& mat,
+                  Vector& ee, // strains
+                  ElementTransformation &T,
+                  const IntegrationPoint &ip)
+    {
+        mat.SetSize(9);
+        EvalGrad(mat.GetData(), ee.GetData(), T, ip);
+    }
+
+private:
+    double E;
+    double nu;
+};
+
+
+class J2YieldFunction
+{
+public:
+
+    J2YieldFunction()
+    {
+        sigma_0=1.0;
+        H=0.0;
+        beta=0.0;
+    }
+
+    J2YieldFunction(double ssy_,double H_, double beta_)
+    {
+        sigma_0=ssy_;
+        H=H_;
+        beta=beta_;
+    }
+
+    void Set(double ssy_,double H_, double beta_)
+    {
+        sigma_0=ssy_;
+        H=H_;
+        beta=beta_;
+    }
+
+    void SetJ2Param(double ssy_, double H_, double beta_)
+    {
+        sigma_0=ssy_;
+        H=H_;
+        beta=beta_;
+    }
+
+    //str[9] - stress tensor
+    //ip[2]:
+    //      ip[0] accumulated plastic strain (possitive)
+    //      ip[1] filtered accumulated plastic strain (possitive)
+    template<typename dtype0,typename dtype1,typename dtype2 >
+    dtype0 EvalI(dtype1* str,dtype2* ip)
+    {
+        dtype0 p=(str[0]+str[4]+str[8])/3.0;
+        dtype0 dstr[9];
+        for(int i=0;i<9;i++){ dstr[i]=str[i];}
+        dstr[0]=dstr[0]-p;
+        dstr[4]=dstr[4]-p;
+        dstr[8]=dstr[8]-p;
+        dtype0 se=dstr[0]*dstr[0];
+        for(int i=1;i<9;i++){ se=se+dstr[i]*dstr[i];}
+        se=sqrt(3.0*se/2.0);
+
+        dtype0 sy=sigma_0+H*ip[0];
+        dtype0 tv=exp(-beta*ip[1]);
+        //return se-tv*tv*sy*sy;
+        return se-tv*sy;
+    }
+
+    //return gradients of the yield function with respect to str and ip
+    template<typename dtype>
+    dtype EvalG(dtype* gstr, dtype* gip, dtype* str, dtype* ip)
+    {
+        dtype t1,t7,t9,t11,t13,t17,t19,t21,t23,t25,t26,t27,t29,t30,t40;
+        t1 = str[0]*str[0];
+        t7 = str[1]*str[1];
+        t9 = str[2]*str[2];
+        t11 = str[3]*str[3];
+        t13 = str[4]*str[4];
+        t17 = str[5]*str[5];
+        t19 = str[6]*str[6];
+        t21 = str[7]*str[7];
+        t23 = str[8]*str[8];
+        t25 = -4.0*str[0]*str[4]-4.0*str[0]*str[8]-4.0*str[4]*str[8]+4.0*t1+6.0*
+                t11+4.0*t13+6.0*t17+6.0*t19+6.0*t21+4.0*t23+6.0*t7+6.0*t9;
+        t26 = sqrt(t25);
+        t27 = 1.0/t26;
+        t29 = 4.0*str[4];
+        t30 = 4.0*str[8];
+        t40 = 4.0*str[0];
+        gstr[0] = t27*(8.0*str[0]-t29-t30)/4.0;
+        gstr[1] = 3.0*t27*str[1];
+        gstr[2] = 3.0*t27*str[2];
+        gstr[3] = 3.0*t27*str[3];
+        gstr[4] = t27*(-t40+8.0*str[4]-t30)/4.0;
+        gstr[5] = 3.0*t27*str[5];
+        gstr[6] = 3.0*t27*str[6];
+        gstr[7] = 3.0*t27*str[7];
+        gstr[8] = t27*(-t40-t29+8.0*str[8])/4.0;
+
+        gip[0] = -exp(-beta*ip[1])*H;
+        gip[1] = beta*exp(-beta*ip[1])*(H*ip[0]+sigma_0);
+
+        dtype sy=sigma_0+H*ip[0];
+        dtype tv=exp(-beta*ip[1]);
+        return t26/2.0-tv*sy;
+    }
+
+    double Eval(const Vector& str,const Vector& ivar)
+    {
+        return EvalI<double,double,double>(str.GetData(),ivar.GetData());
+    }
+
+    double Eval(double* str, double* ivar)
+    {
+        return EvalI<double,double,double>(str,ivar);
+    }
+
+    double EvalFGrad(Vector& drstr, Vector& divar,
+                    Vector& str, Vector& ivar)
+    {
+        //drstr.SetSize(str.Size());
+        //divar.SetSize(ivar.Size());
+        typedef internal::dual<double, double> ADFloatType;
+        ADFloatType tst[str.Size()];
+        ADFloatType var[ivar.Size()];
+
+        for(int i=0;i<str.Size();i++){
+            tst[i].value=str[i];
+            tst[i].gradient=0.0;
+        }
+
+        for(int i=0;i<ivar.Size();i++){
+            var[i].value=ivar[i];
+            var[i].gradient=0.0;
+        }
+
+        ADFloatType rez;
+        for(int i=0;i<str.Size();i++){
+            tst[i].gradient=1.0;
+            rez=EvalI<ADFloatType,ADFloatType,ADFloatType>(tst,var);
+            drstr[i]=rez.gradient;
+            tst[i].gradient=0.0;
+        }
+
+        for(int i=0;i<ivar.Size();i++){
+            var[i].gradient=1.0;
+            rez=EvalI<ADFloatType,ADFloatType,ADFloatType>(tst,var);
+            divar[i]=rez.gradient;
+            var[i].gradient=0.0;
+        }
+
+        return rez.value;
+    }
+
+private:
+    double H;
+    double beta;
+    double sigma_0;
+
+};
+
+
+class DruckerPragerYieldFunction
+{
+public:
+    DruckerPragerYieldFunction(double ssy_,double alpha_)
+    {
+        sigma_0=ssy_;
+        alpha=alpha_;
+    }
+
+    template<typename dtype0, typename dtype1, typename dtype2>
+    dtype0 Eval(dtype1* str,dtype2* ip)
+    {
+        dtype0 p=(str[0]+str[4]+str[8])/3.0;
+        dtype0 dstr[9];
+        for(int i=0;i<9;i++){ dstr[i]=str[i];}
+        dstr[0]=dstr[0]-p;
+        dstr[4]=dstr[1]-p;
+        dstr[8]=dstr[8]-p;
+        dtype0 J2=dstr[0]*dstr[0];
+        for(int i=1;i<9;i++){ J2=J2+dstr[i]*dstr[i];}
+        J2=J2/2.0;
+        return sqrt(J2)-3.0*alpha*p-sigma_0/sqrt(3.0);
+    }
+
+    double operator()(Vector& str,Vector& ivar)
+    {
+        return Eval<double>(str.GetData(),ivar.GetData());
+    }
+
+private:
+
+    double sigma_0;
+    double alpha;
+};
+
+class MatsuokaNakaiYieldFunction
+{
+public:
+    MatsuokaNakaiYieldFunction(double phi_)
+    {
+        phi=phi_;
+    }
+
+    template<typename dtype0, typename dtype1, typename dtype2>
+    dtype0 Eval(dtype1* str,dtype2* ip)
+    {
+        dtype0 I1=(str[0]+str[4]+str[8]);
+        dtype0 I2=str[0]*str[4]+str[4]*str[8]+str[8]*str[0]
+                              -str[3]*str[3]-str[6]*str[6]-str[7]*str[7];
+        dtype0 I3= str[0]*str[4]*str[8]+str[1]*str[5]*str[6]+str[3]*str[7]*str[2]
+                 -str[2]*str[4]*str[6]-str[1]*str[3]*str[8]-str[5]*str[7]*str[0];
+
+        double tp=tan(phi);
+        return I1*I2 - (9.0+8.0*tp*tp)*I3;
+    }
+
+    double operator()(Vector& str,Vector& ivar)
+    {
+        return Eval<double>(str.GetData(),ivar.GetData());
+    }
+private:
+    double phi;
+};
+
+
+
+
+
+template<typename elmat, typename yfunc>
+class StressEval
+{
+public:
+
+    StressEval(elmat* mat_, yfunc* yf_, double lerr_=1e-8){
+        mat=mat_;
+        yf=yf_;
+        lerr=lerr_;
+    }
+
+    void SetPlasticStrain(const Vector& ep_)
+    {
+        for(int i=0;i<9;i++){
+            epn[i]=ep_[i];
+        }
+    }
+
+    void SetInternalParameters(const Vector& ip_)
+    {
+        ipn[0]=ip_[0];
+        ipn[1]=ip_[1];
+    }
+
+    void SetStrain(const Vector &ee_)
+    {
+        for(int i=0;i<9;i++){
+            ee[i]=ee_[i];
+        }
+    }
+
+    void EvalTangent(DenseMatrix& Cep,
+                     ElementTransformation &T,
+                     const IntegrationPoint &ip)
+    {
+        Vector ees(9);
+        Vector css(9);
+        for(int i=0;i<9;i++){
+            ees[i]=ee[i]-epn[i];
+        }
+
+        mat->EvalStress(css.GetData(),ees.GetData(),T,ip);
+        mat->EvalGrad(Cep,ees,T,ip);
+        double f=yf->Eval(css.GetData(),ipn);
+        if(f<0.0){
+            return;}
+        Vector dip(2);
+        //elasto-plastic behaviour - modify the elastic tensor
+        Vector r(9);
+        //yf->EvalFGrad(r,dip,css,ipn);
+        yf->EvalG(r.GetData(), dip.GetData(), css.GetData(),ipn);
+        double H;
+
+        mat->EvalStress(css.GetData(),r.GetData(),T,ip);
+        H=(css*r)-dip[0]*sqrt(2.0*(r*r)/3.0);
+        for(int i=0;i<9;i++){
+        for(int j=0;j<9;j++){
+            Cep(i,j)=Cep(i,j)-css(i)*css(j)/H;
+        }}
+
+    }
+
+    /// vin[stress[9], strainp[9], ip[2], multiplier[1]]
+    template<typename dtype>
+    void EvalResidual(dtype* rr, dtype* vin,
+                      ElementTransformation &T,
+                      const IntegrationPoint &ip)
+    {
+        dtype* css=vin;   //current stress
+        dtype* cep=vin+9; //current plastic strain
+        dtype* cip=vin+18;//current internal variables
+        dtype* lam=vin+20;//current Largrange multiplier
+
+        dtype* rr1=rr;
+        dtype* rr2=rr+9;
+        dtype* rr3=rr+18;
+        dtype* rr4=rr+20;
+
+        dtype tv[9];
+        dtype hh[2];
+        dtype vv;
+
+        //stress residual
+        for(int i=0;i<9;i++){
+            tv[i]=ee[i]-cep[i];
+        }
+        mat->EvalStress(rr1,tv,T,ip);
+        for(int i=0;i<9;i++){
+            rr1[i]=css[i]-rr1[i];
+        }
+
+        rr4[0]=yf->EvalG(tv, hh, css, cip);
+        for(int i=0;i<9;i++){
+            rr2[i]=cep[i]-epn[i]-(*lam)*tv[i];
+        }
+
+        //internal parameters
+        vv=tv[0]*tv[0];
+        for(int i=1;i<9;i++){ vv=vv+tv[i]*tv[i];}
+        rr3[0]=cip[0]-ipn[0]-(*lam)*sqrt((2.0/3.0)*vv);
+        rr3[1]=cip[1]-ipn[1];
+    }
+
+    //E
+    void EvalNewton(DenseMatrix& tmat, Vector& rr,
+                    Vector& vv,
+                    ElementTransformation &T,
+                    const IntegrationPoint &ip){
+        typedef internal::dual<double, double> ADFloatType;
+        ADFloatType tv[21];
+        ADFloatType cv[21];
+        for(int i=0;i<21;i++){
+            cv[i].value=vv[i];
+            cv[i].gradient=0.0;
+        }
+
+        for(int i=0;i<21;i++){
+            cv[i].gradient=1.0;
+            EvalResidual(tv, cv,T,ip);
+            for(int j=0;j<21;j++){
+                tmat(j,i)=tv[j].gradient;
+            }
+            cv[i].gradient=0.0;
+        }
+
+        for(int i=0;i<21;i++){
+            rr[i]=tv[i].value;
+        }
+    }
+
+    void Solve(Vector& css, Vector& cep, Vector& cip,
+               ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        Vector ees(9);
+        for(int i=0;i<9;i++){
+            ees[i]=ee[i]-epn[i];
+            cep[i]=epn[i];
+        }
+
+        cip[0]=ipn[0];
+        cip[1]=ipn[1];
+
+        mat->EvalStress(css.GetData(),ees.GetData(),T,ip);
+
+        double f=yf->Eval(css,cip);
+        if(f<0.0){return;}
+        SolveEP(css,cep,cip,T,ip);
+    }
+
+    void SolveEP(Vector& css, Vector& cep, Vector& cip,
+               ElementTransformation &T,
+               const IntegrationPoint &ip, double err=1e-8){
+        Vector vv(21);
+        Vector rr(21);
+        Vector tv(21);
+        DenseMatrix mm(21);
+        for(int i=0;i<9;i++){
+            vv[i]=css[i];
+            vv[i+9]=cep[i];
+        }
+        vv[18]=cip[0];
+        vv[19]=cip[1];
+        vv[20]=0.0;
+
+        EvalNewton(mm,rr,vv,T,ip);
+
+        /*
+        std::cout<<std::endl;
+        mm.PrintMatlab(std::cout);
+        std::cout<<"rr="<<std::endl;
+        rr.Print(std::cout,21);
+        */
+
+        double cerr=rr.Norml2();
+        int it=0;
+        while(cerr>err){
+            DenseMatrixInverse im(mm);
+            im.Mult(rr,tv);
+            vv.Add(-1.0,tv);
+            EvalNewton(mm,rr,vv,T,ip);
+            cerr=rr.Norml2();
+            //std::cout<<"it:"<<it<<" err="<<cerr<<std::endl;
+            it++;
+            if(it>100){
+                mfem_error("Maximum number (100) of iterations have been"
+                           "reached in the elasto-plastic material update! \n");
+                break;}
+        }
+
+        for(int i=0;i<9;i++){
+            css[i]=vv[i];
+            cep[i]=vv[i+9];
+        }
+        cip[0]=vv[18];
+        cip[1]=vv[19];
+    }
+
+
+
+private:
+    elmat* mat;
+    yfunc* yf;
+
+    double ee[9]; //total strain
+    double epn[9]; //plastic strain
+    double ipn[2]; //internal variables
+
+    double lerr;
+
+};
+
+
+}
+
+#endif

miniapps/elplast/elpl_test.cpp

@@ -0,0 +1,331 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "coefficients.hpp"
+#include "elplast.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+class MyVectorCoeff:public VectorCoefficient
+{
+public:
+    MyVectorCoeff():VectorCoefficient(3)
+    {
+        SetTime(1.0);
+    }
+
+    virtual
+    ~MyVectorCoeff(){}
+
+
+    void Eval(Vector &V, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        Vector xx(3);
+        T.Transform(ip,xx);
+        V=xx; V*=0.0;
+        V[0]=xx[0]*GetTime();
+        //V[1]=xx[0]*GetTime();
+        //V[2]=xx[2]*GetTime();
+    }
+
+private:
+};
+
+class MyForceCoeff:public VectorCoefficient
+{
+public:
+    MyForceCoeff():VectorCoefficient(3)
+    {
+        SetTime(1.0);
+    }
+
+    virtual
+    ~MyForceCoeff(){}
+
+
+    void Eval(Vector &V, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        Vector xx(3);
+        T.Transform(ip,xx);
+        V=xx; V*=0.0;
+        if(xx[0]>1.9){V[0]=GetTime();}
+        //V[1]=xx[1]*GetTime();
+        //V[2]=xx[2]*GetTime();
+    }
+
+private:
+};
+
+
+class MyBlockPrec:public BlockDiagonalPreconditioner
+{
+public:
+    MyBlockPrec(const Array<int> & offsets, ParFiniteElementSpace* fesel):BlockDiagonalPreconditioner(offsets)
+    {
+        pr1=new HypreBoomerAMG();
+        pr1->SetPrintLevel(0);
+        pr1->SetElasticityOptions(fesel);
+        pr2=new HypreBoomerAMG();
+        pr2->SetPrintLevel(0);
+
+    }
+
+    virtual void SetOperator(const Operator &op)
+    {
+        const BlockOperator& blo=static_cast<const BlockOperator&>(op);
+        pr1->SetOperator(blo.GetBlock(0,0));
+        pr2->SetOperator(blo.GetBlock(1,1));
+        SetDiagonalBlock(0,pr1);
+        SetDiagonalBlock(1,pr2);
+    }
+
+    virtual ~MyBlockPrec()
+    {
+        delete pr1;
+        delete pr2;
+    }
+
+private:
+    HypreBoomerAMG* pr1;
+    HypreBoomerAMG* pr2;
+
+};
+
+
+
+
+int main(int argc, char *argv[])
+{
+    // 0. Initialize MPI and HYPRE.
+    Mpi::Init(argc, argv);
+    int num_procs = Mpi::WorldSize();
+    int myid = Mpi::WorldRank();
+    Hypre::Init();
+
+    // 1. Parse command-line options.
+    const char *mesh_file = "../data/beam-tri.mesh";
+    int serial_ref_levels = 0;
+    int order = 2;
+    //bool static_cond = false;
+    bool static_cond = true;
+    bool visualization = 1;
+
+    OptionsParser args(argc, argv);
+    args.AddOption(&mesh_file, "-m", "--mesh",
+                   "Mesh file to use.");
+    args.AddOption(&serial_ref_levels, "-rs", "--refine-serial",
+                   "Number of uniform serial refinements (before parallel"
+                   " partitioning)");
+    args.AddOption(&order, "-o", "--order",
+                   "Finite element order (polynomial degree).");
+    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                   "--no-static-condensation", "Enable static condensation.");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                   "--no-visualization",
+                   "Enable or disable GLVis visualization.");
+    args.Parse();
+    if (!args.Good())
+    {
+       if (myid == 0)
+       {
+          args.PrintUsage(cout);
+       }
+       return 1;
+    }
+    if (myid == 0)
+    {
+       args.PrintOptions(cout);
+    }
+
+    // 2. Read the mesh from the given mesh file. We can handle triangular,
+    //    quadrilateral, tetrahedral, and hexahedral meshes with the same code.
+    Mesh mesh(mesh_file, 1, 1);
+    int dim = mesh.Dimension();
+    MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
+
+    // 3. Refine the mesh before parallel partitioning. Since a NURBS mesh can
+    //    currently only be refined uniformly, we need to convert it to a
+    //    piecewise-polynomial curved mesh. First we refine the NURBS mesh a bit
+    //    more and then project the curvature to quadratic Nodes.
+    if (mesh.NURBSext && serial_ref_levels == 0)
+    {
+       serial_ref_levels = 2;
+    }
+    for (int i = 0; i < serial_ref_levels; i++)
+    {
+       mesh.UniformRefinement();
+    }
+    if (mesh.NURBSext)
+    {
+       mesh.SetCurvature(2);
+    }
+    mesh.EnsureNCMesh();
+
+    ParMesh pmesh(MPI_COMM_WORLD, mesh);
+    mesh.Clear();
+
+    // 4. Define a finite element space on the mesh. The polynomial order is
+    //    one (linear) by default, but this can be changed on the command line.
+    H1_FECollection fec(order, dim);
+    ParFiniteElementSpace ufespace(&pmesh, &fec, dim);
+    ParFiniteElementSpace efespace(&pmesh, &fec, 1);
+
+    QuadratureSpace  qfes(&pmesh,3*order);
+    QuadratureFunction kappa;
+    QuadratureFunction eep;
+
+    kappa.SetSpace(&qfes,1); kappa=0.0;
+    eep.SetSpace(&qfes,6);   eep=0.0;
+
+    ParGridFunction u(&ufespace); u=0.0;
+    ParGridFunction ep(&efespace); ep=0.0;
+    MyVectorCoeff co; co.SetTime(0.5);
+    MyForceCoeff  fo; fo.SetTime(0.5);
+
+    u.ProjectCoefficient(co);
+
+    Array<ParFiniteElementSpace *> pfes(2);
+    pfes[0]=&ufespace;
+    pfes[1]=&efespace;
+    ParBlockNonlinearForm* bnl=new ParBlockNonlinearForm();
+    bnl->SetParSpaces(pfes);
+
+    //add the integrator
+    NLElPlastIntegrator* itgr=new NLElPlastIntegrator(1,0.2,1.0,0.01);
+    itgr->SetPlasticStrains(eep,kappa);
+    //itgr->SetForce(fo);
+    bnl->AddDomainIntegrator(itgr);
+
+    Array<int> offset;
+    BlockVector bv; bv.Update(bnl->GetBlockTrueOffsets()); bv=0.0;
+    BlockVector rv; rv.Update(bnl->GetBlockTrueOffsets()); rv=0.0;
+
+    u.GetTrueDofs(bv.GetBlock(0));
+
+    bnl->Mult(bv,rv);
+
+    BlockOperator& op=bnl->GetGradient(bv);
+
+
+    // Define the essential boundary attributes
+    Array<int> ess_bdr_u(pmesh.bdr_attributes.Max());
+    Array<int> ess_bdr_e(pmesh.bdr_attributes.Max());
+    Array<Array<int>*> ess_bdr(2);
+    ess_bdr[0]=&ess_bdr_u;
+    ess_bdr[1]=&ess_bdr_e;
+
+    ess_bdr_u = 1; ess_bdr_u[0] = 1;
+    ess_bdr_e = 0;
+
+    Array<Vector*> nrhs(2);
+    nrhs[0]=nullptr;
+    nrhs[1]=nullptr;
+
+    bnl->SetEssentialBC(ess_bdr,nrhs);
+
+
+    ParGridFunction du(u);
+    du.SetFromTrueDofs(bv.GetBlock(0));
+    ep.SetFromTrueDofs(bv.GetBlock(1));
+
+    // ParaView output.
+    ParaViewDataCollection dacol("ParaView", &pmesh);
+    dacol.SetLevelsOfDetail(order);
+    dacol.SetDataFormat(VTKFormat::ASCII);
+    dacol.RegisterField("odisp", &u);
+    dacol.RegisterField("ddisp", &du);
+    dacol.RegisterField("ep",&ep);
+
+    /*
+    dacol.SetTime(1.0);
+    dacol.SetCycle(1);
+    dacol.Save();
+    */
+
+
+
+
+    /*
+    BlockDiagonalPreconditioner* blprec=new BlockDiagonalPreconditioner(bnl->GetBlockTrueOffsets());
+    HypreBoomerAMG* pr1=new HypreBoomerAMG();
+    pr1->SetPrintLevel(1);
+    HypreBoomerAMG* pr2=new HypreBoomerAMG();
+    pr2->SetPrintLevel(1);
+    */
+
+
+    MyBlockPrec* blprec=new MyBlockPrec(bnl->GetBlockTrueOffsets(), &ufespace);
+
+    /*
+    GMRESSolver* gmres = new GMRESSolver(MPI_COMM_WORLD);
+    gmres->SetAbsTol(1e-12);
+    gmres->SetRelTol(1e-9);
+    gmres->SetMaxIter(500);
+    gmres->SetPrintLevel(2);
+    gmres->SetPreconditioner(*blprec);
+    */
+
+    CGSolver* cg = new CGSolver(MPI_COMM_WORLD);
+    cg->SetAbsTol(1e-12);
+    cg->SetRelTol(1e-9);
+    cg->SetMaxIter(500);
+    cg->SetPrintLevel(2);
+    cg->SetPreconditioner(*blprec);
+
+
+    NewtonSolver* ns=new  NewtonSolver(MPI_COMM_WORLD);
+
+    ns->iterative_mode = true;
+    ns->SetSolver(*cg);
+    ns->SetOperator(*bnl);
+    ns->SetPrintLevel(1);
+    ns->SetRelTol(1e-5);
+    ns->SetAbsTol(1e-12);
+    ns->SetMaxIter(8);
+
+    u.GetTrueDofs(bv.GetBlock(0));
+    //bv=0.0;
+    Vector b;
+    //ns->Mult(b,bv);
+
+    for(int bi=0;bi<10;bi++){
+        co.SetTime(0.4+bi*0.1);
+        u.ProjectCoefficient(co);
+        u.GetTrueDofs(bv.GetBlock(0));
+        ns->Mult(b,bv);
+
+        //update the plasticity vars
+        itgr->SetUpdateFlag(true);
+        bnl->Mult(bv,rv);
+        itgr->SetUpdateFlag(false);
+
+        u.SetFromTrueDofs(bv.GetBlock(0));
+        ep.SetFromTrueDofs(bv.GetBlock(1));
+
+        dacol.SetTime(1.0+bi);
+        dacol.SetCycle(1+bi);
+        dacol.Save();
+
+    }
+
+
+    delete ns;
+    delete cg;
+    delete blprec;
+
+    MPI_Barrier(MPI_COMM_WORLD);
+
+
+
+    delete bnl;
+
+
+
+
+    return 0;
+}

miniapps/elplast/elplast.cpp

@@ -0,0 +1,904 @@
+#include "elplast.hpp"
+
+#include <fstream>
+#include <iostream>
+
+
+namespace mfem{
+
+
+ElPlastSolver::ElPlastSolver(mfem::ParMesh* mesh_,int vorder, int forder)
+{
+
+    pmesh=mesh_;
+
+    int dim=pmesh->Dimension();
+    vfec=new H1_FECollection(vorder,dim);
+    ffec=new H1_FECollection(forder,dim);
+
+    vfes=new ParFiniteElementSpace(pmesh,vfec,dim,Ordering::byVDIM);
+    ffes=new ParFiniteElementSpace(pmesh,ffec);
+
+    fdisp.SetSpace(vfes); fdisp=0.0;
+    adisp.SetSpace(vfes); adisp=0.0;
+
+
+    qfes=new QuadratureSpace(pmesh,3*std::max(vorder,forder));
+    kappa.SetSpace(qfes,1); kappa=0.0;
+    eep.SetSpace(qfes,6);   eep=0.0;
+    eee.SetSpace(qfes,6); eee=0.0;
+
+    nf=nullptr;
+    Array<mfem::ParFiniteElementSpace*> pf;
+    pf.Append(vfes);
+    pf.Append(ffes);
+
+    nf=new ParBlockNonlinearForm(pf);
+    rhs.Update(nf->GetBlockTrueOffsets()); rhs=0.0;
+    sol.Update(nf->GetBlockTrueOffsets()); sol=0.0;
+    adj.Update(nf->GetBlockTrueOffsets()); adj=0.0;
+
+
+    SetNewtonSolver();
+    SetLinearSolver();
+
+}
+
+
+ElPlastSolver::~ElPlastSolver()
+{
+
+    delete nf;
+
+    delete vfes;
+    delete ffes;
+    delete ffec;
+    delete vfec;
+
+    //delete the local forces
+    for(auto it=lvforce.begin();it!=lvforce.end();it++)
+    {
+        delete it->second;
+    }
+}
+
+void ElPlastSolver::SetNewtonSolver(double rtol, double atol,int miter, int prt_level)
+{
+    rel_tol=rtol;
+    abs_tol=atol;
+    max_iter=miter;
+    print_level=prt_level;
+}
+
+void ElPlastSolver::SetLinearSolver(double rtol, double atol, int miter)
+{
+    linear_rtol=rtol;
+    linear_atol=atol;
+    linear_iter=miter;
+}
+
+void ElPlastSolver::AddDispBC(int id, int dir, double val)
+{
+    if(dir==0){
+        bcx[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,dir,bcx[id]);
+    }
+    if(dir==1){
+        bcy[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,dir,bcy[id]);
+
+    }
+    if(dir==2){
+        bcz[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,dir,bcz[id]);
+    }
+    if(dir==4){
+        bcx[id]=mfem::ConstantCoefficient(val);
+        bcy[id]=mfem::ConstantCoefficient(val);
+        bcz[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,0,bcx[id]);
+        AddDispBC(id,1,bcy[id]);
+        AddDispBC(id,2,bcz[id]);
+    }
+}
+
+void ElPlastSolver::AddDispBC(int id, int dir, Coefficient &val)
+{
+    if(dir==0){ bccx[id]=&val; }
+    if(dir==1){ bccy[id]=&val; }
+    if(dir==2){ bccz[id]=&val; }
+    if(dir==4){ bccx[id]=&val; bccy[id]=&val; bccz[id]=&val;}
+    if(pmesh->Dimension()==2)
+    {
+        bccz.clear();
+    }
+}
+
+void ElPlastSolver::AddDispBC(int id, mfem::VectorCoefficient& val)
+{
+    bcca[id]=&val;
+}
+
+void ElPlastSolver::AddVolForce(int id, double fx, double fy, double fz)
+{
+    //check if the force id is already in the database
+    if(lvforce.find(id)!=lvforce.end()){
+        delete lvforce[id];
+    }
+    Vector tv(3); tv[0]=fx; tv[1]=fy; tv[2]=fz;
+    lvforce[id]=new mfem::VectorConstantCoefficient(tv);
+    volforce[id]=lvforce[id];
+}
+
+/// Adds vol force
+void ElPlastSolver::AddVolForce(int id, mfem::VectorCoefficient& ff)
+{
+    volforce[id]=&ff;
+}
+
+void ElPlastSolver::FSolve()
+{
+    // Set the BC
+    ess_tdofv.DeleteAll();
+    Array<int> ess_tdofx;
+    Array<int> ess_tdofy;
+    Array<int> ess_tdofz;
+
+    int dim=pmesh->Dimension();
+    {
+        for(auto it=bccx.begin();it!=bccx.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,0);
+            ess_tdofx.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(0, it->second, false);
+            fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+
+        //copy tdofsx from velocity grid function
+        {
+            fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
+            for(int ii=0;ii<ess_tdofx.Size();ii++)
+            {
+                sol.GetBlock(0)[ess_tdofx[ii]]=rhs.GetBlock(0)[ess_tdofx[ii]];
+            }
+        }
+        ess_tdofv.Append(ess_tdofx);
+
+        for(auto it=bccy.begin();it!=bccy.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,1);
+            ess_tdofy.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(1, it->second, false);
+            fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsy from velocity grid function
+        {
+            fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
+            for(int ii=0;ii<ess_tdofy.Size();ii++)
+            {
+                sol.GetBlock(0)[ess_tdofy[ii]]=rhs.GetBlock(0)[ess_tdofy[ii]];
+            }
+        }
+        ess_tdofv.Append(ess_tdofy);
+
+        if(dim==3){
+            for(auto it=bccz.begin();it!=bccz.end();it++)
+            {
+                mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+                ess_bdr=0;
+                ess_bdr[it->first -1]=1;
+                mfem::Array<int> ess_tdof_list;
+                vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,2);
+                ess_tdofz.Append(ess_tdof_list);
+
+                mfem::VectorArrayCoefficient pcoeff(dim);
+                pcoeff.Set(2, it->second, false);
+                fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
+            }
+
+            //copy tdofsz from velocity grid function
+            {
+                fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
+                for(int ii=0;ii<ess_tdofz.Size();ii++)
+                {
+                    sol[ess_tdofz[ii]]=rhs[ess_tdofz[ii]];
+                }
+            }
+            ess_tdofv.Append(ess_tdofz);
+        }
+
+        //set vector coefficients
+        for(auto it=bcca.begin();it!=bcca.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+            fdisp.ProjectBdrCoefficient(*(it->second), ess_bdr);
+            //copy tdofs from velocity grid function
+            fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
+            for(int ii=0;ii<ess_tdof_list.Size();ii++)
+            {
+                sol[ess_tdof_list[ii]]=rhs[ess_tdof_list[ii]];
+            }
+            ess_tdofv.Append(ess_tdof_list);
+        }
+    }
+
+
+
+}
+
+
+
+void NLElPlastIntegrator::AssembleElementVector(const Array<const FiniteElement *> &el,
+                                                ElementTransformation &Tr,
+                                                const Array<const Vector *> &elfun,
+                                                const Array<Vector *> &elvec)
+{
+    //the integrator works only for 3 dimensional problems
+    int dof_u = el[0]->GetDof();
+    int dof_e = el[1]->GetDof();
+
+    int dim =  Tr.GetSpaceDim();
+
+    elvec[0]->SetSize(dim*dof_u);
+    elvec[1]->SetSize(dof_e);
+
+    if (dim != 3)
+    {
+        mfem::mfem_error("NLElPlastIntegrator::AssembleElementVector"
+                                " is not defined on manifold meshes");
+    }
+
+    Vector uu(elfun[0]->GetData()+0*dof_u, dof_u);
+    Vector vv(elfun[0]->GetData()+1*dof_u, dof_u);
+    Vector ww(elfun[0]->GetData()+2*dof_u, dof_u);
+
+    Vector ru(elvec[0]->GetData()+0*dof_u, dof_u); ru=0.0;
+    Vector rv(elvec[0]->GetData()+1*dof_u, dof_u); rv=0.0;
+    Vector rw(elvec[0]->GetData()+2*dof_u, dof_u); rw=0.0;
+
+    Vector ep(elfun[1]->GetData(), dof_e);
+    Vector rp(elvec[1]->GetData(), dof_e); rp=0.0;
+
+    // temp storages for vectors and matrices
+    Vector su(dof_u); //shape functions for displacements
+    DenseMatrix du(dof_u,dim); //gradients of the shape functions
+    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
+    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
+    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
+    Vector se(dof_e); //shape functions for plastic strains
+    DenseMatrix de(dof_e,dim); //gradients of the shape functions
+    Vector grade(dim);
+
+    DenseMatrix vpss; //plastic strain at the integration points
+    DenseMatrix vkap; //accumulated plastic strain at the integration points
+    eep->GetElementValues(Tr.ElementNo,vpss);
+    kappa->GetElementValues(Tr.ElementNo,vkap);
+
+
+    DenseMatrix estrain(3,3);
+    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
+    DenseMatrix pstrain(3,3);
+    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
+    DenseMatrix gradu(3,3); gradu=0.0;
+    Vector tv;
+
+    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
+
+    double H=0.001;
+    double beta=0.0;
+    Vector vip(2); vip=0.0;
+
+    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
+    Vector stress(9);stress=0.0;
+
+    double w;
+    for (int i = 0; i < ir.GetNPoints(); i++)
+    {
+        const IntegrationPoint &ip = ir.IntPoint(i);
+        Tr.SetIntPoint(&ip);
+        w = Tr.Weight();
+        w = ip.weight * w;
+
+
+        //calculate the total strain
+        el[0]->CalcPhysShape(Tr,su);
+        el[0]->CalcPhysDShape(Tr,du);
+        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
+        du.MultTranspose(uu,tv);
+        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
+        du.MultTranspose(vv,tv);
+        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
+        du.MultTranspose(ww,tv);
+        //compute the strain tensor
+        for(int ii=0;ii<dim;ii++){
+            for(int jj=ii+1;jj<dim;jj++){
+                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
+                estrain(jj,ii)=estrain(ii,jj);
+            }
+            estrain(ii,ii)=gradu(ii,ii);
+        }
+
+        //estrain.PrintMatlab(std::cout);
+
+        //set current plastic strain Voight indexing
+        {
+            pstrain(0,0)=vpss(0,i);
+            pstrain(1,1)=vpss(1,i);
+            pstrain(2,2)=vpss(2,i);
+            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
+            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
+            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
+        }
+
+        //calculate the filtered plastic strain
+        el[1]->CalcPhysShape(Tr,se);
+        el[1]->CalcPhysDShape(Tr,de);
+        double epf=se*ep;//filtered accumulated plastic strain
+        double rry=ll->Eval(Tr,ip); //filter radius
+
+        //evaluate the material behaviour at the integration point
+        {
+            double Em=E->Eval(Tr,ip);
+            double nm=nu->Eval(Tr,ip);
+            double ssy=ss_y->Eval(Tr,ip);
+            mat.SetE(Em);
+            mat.SetPoisson(nm);
+            yf.Set(ssy,H,beta);
+
+            vip[0]=vkap(0,i); //set the accumulated plastic starin
+            vip[1]=epf;   //set the filtered plastic strain
+            seval.SetStrain(vestrain);
+            seval.SetPlasticStrain(vpstrain);
+            seval.SetInternalParameters(vip);
+            seval.Solve(stress,vpstrain,vip,Tr,ip);
+
+            if(flag_update){
+                vpss(0,i)=pstrain(0,0);
+                vpss(1,i)=pstrain(1,1);
+                vpss(2,i)=pstrain(2,2);
+                vpss(3,i)=pstrain(1,2);
+                vpss(4,i)=pstrain(0,2);
+                vpss(5,i)=pstrain(0,1);
+
+                vkap(0,i)=vip[0];
+            }
+
+            //stress.Print(std::cout,3);
+            //std::cout<<std::endl;
+        }
+
+        //assemble the RHS for displacements
+        ru.Add(stress[0]*w,dux);
+        ru.Add(0.5*stress[1]*w,duy); rv.Add(0.5*stress[1]*w,dux);
+        ru.Add(0.5*stress[2]*w,duz); rw.Add(0.5*stress[2]*w,dux);
+        ru.Add(0.5*stress[3]*w,duy); rv.Add(0.5*stress[3]*w,dux);
+        rv.Add(stress[4]*w,duy);
+        rv.Add(0.5*stress[5]*w,duz); rw.Add(0.5*stress[5]*w,duy);
+        ru.Add(0.5*stress[6]*w,duz); rw.Add(0.5*stress[6]*w,dux);
+        rv.Add(0.5*stress[7]*w,duz); rw.Add(0.5*stress[7]*w,duy);
+        rw.Add(stress[8]*w,duz);
+
+
+        //assemble the RHS for the accumulated plastic strain
+        rp.Add(-vip(0)*w,se);
+        rp.Add(epf*w,se);
+        de.MultTranspose(ep,grade);
+        tv.SetDataAndSize(de.GetData()+0*dof_e,dof_e);
+        rp.Add(rry*rry*grade[0]*w,tv);
+        tv.SetDataAndSize(de.GetData()+1*dof_e,dof_e);
+        rp.Add(rry*rry*grade[1]*w,tv);
+        tv.SetDataAndSize(de.GetData()+2*dof_e,dof_e);
+        rp.Add(rry*rry*grade[2]*w,tv);
+    }
+
+    //add the force
+    /*
+    if(force)
+    {
+        Vector lv(3*dof_u);
+        VectorDomainLFIntegrator li(*force);
+        li.AssembleRHSElementVect(*el[0],Tr,lv);
+        elvec[0]->Add(-1.0,lv);
+    }*/
+
+    //check the linear elasticity
+
+    /*
+    {
+        const IntegrationPoint &ip = ir.IntPoint(0);
+        double Em=E->Eval(Tr,ip);
+        double num=nu->Eval(Tr,ip);
+        double la=Em*num/((1.0+num)*(1.0-2.0*num));
+        double mu=Em/(2.0*(1+num));
+        ConstantCoefficient lc(la);
+        ConstantCoefficient mc(mu);
+        ElasticityIntegrator eli(lc,mc);
+        DenseMatrix K(dim*dof_u);
+        eli.AssembleElementMatrix(*el[0],Tr,K);
+        Vector rr(dim*dof_u);
+        //K.Mult(*(elfun[0]),rr);
+        //rr.Add(-1.0,*(elvec[0]));
+        *(elvec[0])=0.0;
+        K.Mult(*(elfun[0]),*(elvec[0]));
+        //std::cout<<"|rr|="<<rr.Norml2()<<" "<<elvec[0]->Norml2()<<std::endl;
+    }
+    */
+
+
+
+}
+
+void NLElPlastIntegrator::AssembleElementGrad(const Array<const FiniteElement *> &el,
+                                              ElementTransformation &Tr,
+                                              const Array<const Vector *> &elfun,
+                                              const Array2D<DenseMatrix *> &elmat)
+{
+    //the integrator works only for 3 dimensional problems
+    int dof_u = el[0]->GetDof();
+    int dof_e = el[1]->GetDof();
+
+    elmat(0,0)->SetSize(3*dof_u,3*dof_u); (*elmat(0,0))=0.0;
+    elmat(0,1)->SetSize(3*dof_u,dof_e); (*elmat(0,1))=0.0;
+    elmat(1,0)->SetSize(dof_e,3*dof_u); (*elmat(1,0))=0.0;
+    elmat(1,1)->SetSize(dof_e,dof_e); (*elmat(1,1))=0.0;
+
+
+    int dim =  Tr.GetSpaceDim();
+    if (dim != 3)
+    {
+        mfem::mfem_error("NLElPlastIntegrator::AssembleElementVector"
+                                " is not defined on manifold meshes");
+    }
+
+    Vector uu(elfun[0]->GetData()+0*dof_u, dof_u);
+    Vector vv(elfun[0]->GetData()+1*dof_u, dof_u);
+    Vector ww(elfun[0]->GetData()+2*dof_u, dof_u);
+
+    Vector ep(elfun[1]->GetData(), dof_e);
+
+    // temp storages for vectors and matrices
+    Vector su(dof_u); //shape functions for displacements
+    DenseMatrix du(dof_u,dim); //gradients of the shape functions
+    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
+    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
+    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
+    Vector se(dof_e); //shape functions for plastic strains
+    DenseMatrix de(dof_e,dim); //gradients of the shape functions
+    Vector grade(dim);
+
+    DenseMatrix vpss; //plastic strain at the integration points
+    DenseMatrix vkap; //accumulated plastic strain at the integration points
+    eep->GetElementValues(Tr.ElementNo,vpss);
+    kappa->GetElementValues(Tr.ElementNo,vkap);
+
+
+    DenseMatrix estrain(3,3);
+    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
+    DenseMatrix pstrain(3,3);
+    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
+    DenseMatrix gradu(3,3); gradu=0.0;
+    Vector tv;
+
+    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
+
+    double H=0.001;
+    double beta=0.0;
+    Vector vip(2); vip=0.0;
+
+    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
+    Vector stress(9);stress=0.0;
+
+    DenseMatrix Cep(9,9);
+    DenseMatrix tm(3*dof_u,3*dof_u);
+
+    DenseMatrix B(3*dof_u,9); B=0.0;
+    DenseMatrix Bm(3*dof_u,9);
+
+
+    double w;
+    for (int i = 0; i < ir.GetNPoints(); i++)
+    {
+        const IntegrationPoint &ip = ir.IntPoint(i);
+        Tr.SetIntPoint(&ip);
+        w = Tr.Weight();
+        w = ip.weight * w;
+
+        //calculate the total strain
+        el[0]->CalcPhysShape(Tr,su);
+        el[0]->CalcPhysDShape(Tr,du);
+        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
+        du.MultTranspose(uu,tv);
+        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
+        du.MultTranspose(vv,tv);
+        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
+        du.MultTranspose(ww,tv);
+        //compute the strain tensor
+        for(int ii=0;ii<dim;ii++){
+            for(int jj=ii+1;jj<dim;jj++){
+                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
+                estrain(jj,ii)=estrain(ii,jj);
+            }
+            estrain(ii,ii)=gradu(ii,ii);
+        }
+
+        //estrain.PrintMatlab(std::cout);
+
+        //set current plastic strain Voight indexing
+        {
+            pstrain(0,0)=vpss(0,i);
+            pstrain(1,1)=vpss(1,i);
+            pstrain(2,2)=vpss(2,i);
+            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
+            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
+            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
+        }
+
+        //calculate the filtered plastic strain
+        el[1]->CalcPhysShape(Tr,se);
+        el[1]->CalcPhysDShape(Tr,de);
+        double epf=se*ep;//filtered accumulated plastic strain
+        double rry=ll->Eval(Tr,ip); //filter radius
+
+        //evaluate the material behaviour at the integration point
+        {
+            double Em=E->Eval(Tr,ip);
+            double nm=nu->Eval(Tr,ip);
+            double ssy=ss_y->Eval(Tr,ip);
+            mat.SetE(Em);
+            mat.SetPoisson(nm);
+            yf.Set(ssy,H,beta);
+
+            vip[0]=vkap(0,i); //set the accumulated plastic starin
+            vip[1]=epf;   //set the filtered plastic strain
+            seval.SetStrain(vestrain);
+            seval.SetPlasticStrain(vpstrain);
+            seval.SetInternalParameters(vip);
+            //seval.Solve(stress,vpstrain,vip,Tr,ip);
+            seval.EvalTangent(Cep,Tr,ip);
+
+            //stress.Print(std::cout,3);
+            //Cep.PrintMatlab(std::cut);
+            //std::cout<<std::endl;
+
+        }
+
+        {
+            //form B
+            for(int i=0;i<dof_u;i++)
+            {
+                B(i+0*dof_u,0)=dux[i];
+                B(i+0*dof_u,1)=0.5*duy[i]; B(i+1*dof_u,1)=0.5*dux[i];
+                B(i+0*dof_u,2)=0.5*duz[i]; B(i+2*dof_u,2)=0.5*dux[i];
+                B(i+0*dof_u,3)=0.5*duy[i]; B(i+1*dof_u,3)=0.5*dux[i];
+                B(i+1*dof_u,4)=duy[i];
+                B(i+1*dof_u,5)=0.5*duz[i]; B(i+2*dof_u,5)=0.5*duy[i];
+                B(i+0*dof_u,6)=0.5*duz[i]; B(i+2*dof_u,6)=0.5*dux[i];
+                B(i+1*dof_u,7)=0.5*duz[i]; B(i+2*dof_u,7)=0.5*duy[i];
+                B(i+2*dof_u,8)=duz[i];
+            }
+
+            MultABt(B,Cep,Bm);
+            MultABt(B,Bm,tm);
+            elmat(0,0)->Add(w,tm);
+        }
+
+
+    }
+
+
+    //check the linear elasticity
+    /*
+    {
+        const IntegrationPoint &ip = ir.IntPoint(0);
+        double Em=E->Eval(Tr,ip);
+        double num=nu->Eval(Tr,ip);
+        double la=Em*num/((1.0+num)*(1.0-2.0*num));
+        double mu=Em/(2.0*(1+num));
+        ConstantCoefficient lc(la);
+        ConstantCoefficient mc(mu);
+        ElasticityIntegrator eli(lc,mc);
+        DenseMatrix K(dim*dof_u);
+        eli.AssembleElementMatrix(*el[0],Tr,K);
+        //std::fstream osi("k_mat.dat",std::ios::out);
+        //K.PrintMatlab(osi);
+
+        //K.Add(-1.0,*(elmat(0,0)));
+        //std::cout<<"|K|"<<K.FNorm()<<std::endl;
+        (*elmat(0,0))=K;
+
+    }
+    */
+
+
+    {
+        DenseMatrix K(dof_e,dof_e); K=0.0;
+        PowerCoefficient pc(*ll,2.0);
+        ConstantCoefficient one(1.0);
+        MassIntegrator mi(one);
+        mi.AssembleElementMatrix(*el[1],Tr,*elmat(1,1));
+        DiffusionIntegrator di(pc);
+        di.AssembleElementMatrix(*el[1],Tr,K);
+        elmat(1,1)->Add(1.0,K);
+    }
+
+}
+
+
+void NLElPlastIntegratorS::AssembleElementVector(const FiniteElement &el,
+                                                 ElementTransformation &Tr,
+                                                 const Vector &elfun, Vector &elvect)
+{
+    //the integrator works only for 3 dimensional problems
+    int dof_u = el.GetDof();
+    int dim =  Tr.GetSpaceDim();
+    elvect.SetSize(dim*dof_u); elvect=0.0;
+
+    Vector uu(elfun.GetData()+0*dof_u, dof_u);
+    Vector vv(elfun.GetData()+1*dof_u, dof_u);
+    Vector ww(elfun.GetData()+2*dof_u, dof_u);
+
+    Vector ru(elvect.GetData()+0*dof_u, dof_u); ru=0.0;
+    Vector rv(elvect.GetData()+1*dof_u, dof_u); rv=0.0;
+    Vector rw(elvect.GetData()+2*dof_u, dof_u); rw=0.0;
+
+    // temp storages for vectors and matrices
+    Vector su(dof_u); //shape functions for displacements
+    DenseMatrix du(dof_u,dim); //gradients of the shape functions
+    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
+    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
+    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
+
+    DenseMatrix vpss; //plastic strain at the integration points
+    DenseMatrix vkap; //accumulated plastic strain at the integration points
+    eep->GetElementValues(Tr.ElementNo,vpss);
+    kappa->GetElementValues(Tr.ElementNo,vkap);
+
+    DenseMatrix estrain(3,3);
+    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
+    DenseMatrix pstrain(3,3);
+    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
+    DenseMatrix gradu(3,3); gradu=0.0;
+    Vector tv;
+
+    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
+
+    double H=0.001;
+    double beta=0.0;
+    Vector vip(2); vip=0.0;
+
+    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
+    Vector stress(9);stress=0.0;
+
+    double w;
+    for (int i = 0; i < ir.GetNPoints(); i++)
+    {
+        const IntegrationPoint &ip = ir.IntPoint(i);
+        Tr.SetIntPoint(&ip);
+        w = Tr.Weight();
+        w = ip.weight * w;
+
+
+        //calculate the total strain
+        el.CalcPhysShape(Tr,su);
+        el.CalcPhysDShape(Tr,du);
+        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
+        du.MultTranspose(uu,tv);
+        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
+        du.MultTranspose(vv,tv);
+        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
+        du.MultTranspose(ww,tv);
+        //compute the strain tensor
+        for(int ii=0;ii<dim;ii++){
+            for(int jj=ii+1;jj<dim;jj++){
+                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
+                estrain(jj,ii)=estrain(ii,jj);
+            }
+            estrain(ii,ii)=gradu(ii,ii);
+        }
+
+        //set current plastic strain Voight indexing
+        {
+            pstrain(0,0)=vpss(0,i);
+            pstrain(1,1)=vpss(1,i);
+            pstrain(2,2)=vpss(2,i);
+            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
+            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
+            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
+        }
+
+        double epf=eef->GetValue(Tr,ip);
+
+        //evaluate the material behaviour at the integration point
+        {
+            double Em=E->Eval(Tr,ip);
+            double nm=nu->Eval(Tr,ip);
+            double ssy=ss_y->Eval(Tr,ip);
+            mat.SetE(Em);
+            mat.SetPoisson(nm);
+            yf.Set(ssy,H,beta);
+
+            vip[0]=vkap(0,i); //set the accumulated plastic starin
+            vip[1]=epf;   //set the filtered plastic strain
+            seval.SetStrain(vestrain);
+            seval.SetPlasticStrain(vpstrain);
+            seval.SetInternalParameters(vip);
+            seval.Solve(stress,vpstrain,vip,Tr,ip);
+            //stress.Print(std::cout,3);
+            //std::cout<<std::endl;
+        }
+
+        //assemble the RHS for displacements
+        ru.Add(stress[0]*w,dux);
+        ru.Add(0.5*stress[1]*w,duy); rv.Add(0.5*stress[1]*w,dux);
+        ru.Add(0.5*stress[2]*w,duz); rw.Add(0.5*stress[2]*w,dux);
+        ru.Add(0.5*stress[3]*w,duy); rv.Add(0.5*stress[3]*w,dux);
+        rv.Add(stress[4]*w,duy);
+        rv.Add(0.5*stress[5]*w,duz); rw.Add(0.5*stress[5]*w,duy);
+        ru.Add(0.5*stress[6]*w,duz); rw.Add(0.5*stress[6]*w,dux);
+        rv.Add(0.5*stress[7]*w,duz); rw.Add(0.5*stress[7]*w,duy);
+        rw.Add(stress[8]*w,duz);
+    }
+}
+
+void NLElPlastIntegratorS::AssembleElementGrad (const FiniteElement &el,
+                                                ElementTransformation &Tr,
+                                                const Vector &elfun, DenseMatrix &elmat)
+{
+    //the integrator works only for 3 dimensional problems
+    int dof_u = el.GetDof();
+    int dim =  Tr.GetSpaceDim();
+
+    elmat.SetSize(dof_u,dof_u); elmat=0.0;
+
+    if (dim != 3)
+    {
+        mfem::mfem_error("NLElPlastIntegrator::AssembleElementVector"
+                                " is not defined on manifold meshes");
+    }
+
+    Vector uu(elfun.GetData()+0*dof_u, dof_u);
+    Vector vv(elfun.GetData()+1*dof_u, dof_u);
+    Vector ww(elfun.GetData()+2*dof_u, dof_u);
+
+    // temp storages for vectors and matrices
+    Vector su(dof_u); //shape functions for displacements
+    DenseMatrix du(dof_u,dim); //gradients of the shape functions
+    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
+    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
+    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
+
+    DenseMatrix vpss; //plastic strain at the integration points
+    DenseMatrix vkap; //accumulated plastic strain at the integration points
+    eep->GetElementValues(Tr.ElementNo,vpss);
+    kappa->GetElementValues(Tr.ElementNo,vkap);
+
+
+    DenseMatrix estrain(3,3);
+    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
+    DenseMatrix pstrain(3,3);
+    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
+    DenseMatrix gradu(3,3); gradu=0.0;
+    Vector tv;
+
+
+
+
+    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
+
+
+    double H=0.001;
+    double beta=0.0;
+    Vector vip(2); vip=0.0;
+
+
+    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
+    Vector stress(9);stress=0.0;
+
+    DenseMatrix Cep(9,9);
+    DenseMatrix tm(3*dof_u,3*dof_u);
+
+    DenseMatrix B(3*dof_u,9); B=0.0;
+    DenseMatrix Bm(3*dof_u,9);
+
+    double w;
+    for (int i = 0; i < ir.GetNPoints(); i++)
+    {
+        const IntegrationPoint &ip = ir.IntPoint(i);
+        Tr.SetIntPoint(&ip);
+        w = Tr.Weight();
+        w = ip.weight * w;
+
+        //calculate the total strain
+        el.CalcPhysShape(Tr,su);
+        el.CalcPhysDShape(Tr,du);
+        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
+        du.MultTranspose(uu,tv);
+        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
+        du.MultTranspose(vv,tv);
+        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
+        du.MultTranspose(ww,tv);
+        //compute the strain tensor
+        for(int ii=0;ii<dim;ii++){
+            for(int jj=ii+1;jj<dim;jj++){
+                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
+                estrain(jj,ii)=estrain(ii,jj);
+            }
+            estrain(ii,ii)=gradu(ii,ii);
+        }
+
+        //set current plastic strain Voight indexing
+        {
+            pstrain(0,0)=vpss(0,i);
+            pstrain(1,1)=vpss(1,i);
+            pstrain(2,2)=vpss(2,i);
+            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
+            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
+            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
+        }
+
+        double epf=eef->GetValue(Tr,ip);
+        //evaluate the material behaviour at the integration point
+        {
+            double Em=E->Eval(Tr,ip);
+            double nm=nu->Eval(Tr,ip);
+            double ssy=ss_y->Eval(Tr,ip);
+            mat.SetE(Em);
+            mat.SetPoisson(nm);
+            yf.Set(ssy,H,beta);
+
+            vip[0]=vkap(0,i); //set the accumulated plastic starin
+            vip[1]=epf;   //set the filtered plastic strain
+            seval.SetStrain(vestrain);
+            seval.SetPlasticStrain(vpstrain);
+            seval.SetInternalParameters(vip);
+            //seval.Solve(stress,vpstrain,vip,Tr,ip);
+            seval.EvalTangent(Cep,Tr,ip);
+
+            //stress.Print(std::cout,3);
+            //Cep.PrintMatlab(std::cut);
+            //std::cout<<std::endl;
+
+        }
+
+        {
+            //form B
+            for(int i=0;i<dof_u;i++)
+            {
+                B(i+0*dof_u,0)=dux[i];
+                B(i+0*dof_u,1)=0.5*duy[i]; B(i+1*dof_u,1)=0.5*dux[i];
+                B(i+0*dof_u,2)=0.5*duz[i]; B(i+2*dof_u,2)=0.5*dux[i];
+                B(i+0*dof_u,3)=0.5*duy[i]; B(i+1*dof_u,3)=0.5*dux[i];
+                B(i+1*dof_u,4)=duy[i];
+                B(i+1*dof_u,5)=0.5*duz[i]; B(i+2*dof_u,5)=0.5*duy[i];
+                B(i+0*dof_u,6)=0.5*duz[i]; B(i+2*dof_u,6)=0.5*dux[i];
+                B(i+1*dof_u,7)=0.5*duz[i]; B(i+2*dof_u,7)=0.5*duy[i];
+                B(i+2*dof_u,8)=duz[i];
+            }
+
+            MultABt(B,Cep,Bm);
+            MultABt(B,Bm,tm);
+            elmat.Add(w,tm);
+        }
+
+    }
+
+}
+
+
+}

miniapps/elplast/elplast.hpp

@@ -0,0 +1,310 @@
+#ifndef ELPLAST_HPP
+#define ELPLAST_HPP
+
+#include "mfem.hpp"
+#include "coefficients.hpp"
+
+namespace mfem{
+
+
+
+
+class NLElPlastIntegrator:public BlockNonlinearFormIntegrator
+{
+public:
+    NLElPlastIntegrator(double E_=1, double nu_=0.2, double ss_y_=1.0, double ll_=0.1)
+    {
+        lE=new ConstantCoefficient(E_);
+        lnu=new ConstantCoefficient(nu_);
+        lss_y=new ConstantCoefficient(ss_y_);
+        lll=new ConstantCoefficient(ll_);
+
+        E=lE;
+        nu=lnu;
+        ss_y=lss_y;
+        ll=lll;
+
+        eep=nullptr;
+        kappa=nullptr;
+
+        force=nullptr;
+        flag_update=false;
+    }
+
+    virtual ~NLElPlastIntegrator()
+    {
+        delete lE;
+        delete lnu;
+        delete lss_y;
+        delete lll;
+    }
+
+    void SetPlasticStrains(mfem::QuadratureFunction& eep_,
+                           mfem::QuadratureFunction& kappa_)
+    {
+        eep=&eep_;
+        kappa=&kappa_;
+    }
+
+    void SetForce(VectorCoefficient& vc)
+    {
+        force=&vc;
+    }
+
+    void SetUpdateFlag(bool fl)
+    {
+        flag_update=fl;
+    }
+
+    virtual
+    double GetElementEnergy(const Array<const FiniteElement *> &el,
+                            ElementTransformation &Tr,
+                            const Array<const Vector *> &elfun)
+    {
+        return 0.0;
+    }
+
+
+    virtual
+    void AssembleElementVector(const Array<const FiniteElement *> &el,
+                               ElementTransformation &Tr,
+                               const Array<const Vector *> &elfun,
+                               const Array<Vector *> &elvec);
+
+    virtual
+    void AssembleElementGrad(const Array<const FiniteElement *> &el,
+                             ElementTransformation &Tr,
+                             const Array<const Vector *> &elfun,
+                             const Array2D<DenseMatrix *> &elmats);
+
+private:
+
+    mfem::Coefficient* lE;
+    mfem::Coefficient* lnu;
+    mfem::Coefficient* lss_y;
+    mfem::Coefficient* lll;
+
+
+    mfem::Coefficient* E;
+    mfem::Coefficient* nu;
+    mfem::Coefficient* ss_y;
+    mfem::Coefficient* ll;
+
+    mfem::VectorCoefficient* force;
+
+    mfem::QuadratureFunction* eep;
+    mfem::QuadratureFunction* kappa;
+
+    mfem::IsoElastMat mat;
+    mfem::J2YieldFunction yf;
+
+    bool flag_update;
+
+};
+
+class NLElPlastIntegratorS:public NonlinearFormIntegrator
+{
+public:
+
+    NLElPlastIntegratorS(double E_=1, double nu_=0.2, double ss_y_=1.0, double ll_=0.1)
+    {
+        lE=new ConstantCoefficient(E_);
+        lnu=new ConstantCoefficient(nu_);
+        lss_y=new ConstantCoefficient(ss_y_);
+        lll=new ConstantCoefficient(ll_);
+
+        E=lE;
+        nu=lnu;
+        ss_y=lss_y;
+        ll=lll;
+
+        eep=nullptr;
+        kappa=nullptr;
+    }
+
+    virtual ~NLElPlastIntegratorS()
+    {
+        delete lE;
+        delete lnu;
+        delete lss_y;
+        delete lll;
+    }
+
+    void SetPlasticStrains(mfem::QuadratureFunction& eep_,
+                           mfem::QuadratureFunction& kappa_)
+    {
+        eep=&eep_;
+        kappa=&kappa_;
+    }
+
+    void SetFilteredPlasticStrain(GridFunction& epf)
+    {
+        eef=&epf;
+    }
+
+    virtual
+    double GetElementEnrgy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
+    {
+        return 0.0;
+    }
+
+    virtual
+    void AssembleElementVector (const FiniteElement &el,
+                                ElementTransformation &Tr,
+                                const Vector &elfun, Vector &elvect);
+
+    virtual
+    void AssembleElementGrad (const FiniteElement &el,
+                              ElementTransformation &Tr,
+                              const Vector &elfun, DenseMatrix &elmat);
+
+
+private:
+    mfem::Coefficient* lE;
+    mfem::Coefficient* lnu;
+    mfem::Coefficient* lss_y;
+    mfem::Coefficient* lll;
+
+
+    mfem::Coefficient* E;
+    mfem::Coefficient* nu;
+    mfem::Coefficient* ss_y;
+    mfem::Coefficient* ll;
+
+    mfem::QuadratureFunction* eep;
+    mfem::QuadratureFunction* kappa;
+    mfem::GridFunction* eef;
+
+    mfem::IsoElastMat mat;
+    mfem::J2YieldFunction yf;
+};
+
+
+
+class ElPlastSolver
+{
+public:
+    ElPlastSolver(mfem::ParMesh* mesh_,int vorder=1, int forder=2);
+
+    ~ElPlastSolver();
+
+    /// Set the Newton Solver
+    void SetNewtonSolver(double rtol=1e-7, double atol=1e-12,int miter=1000, int prt_level=1);
+
+    /// Set the Linear Solver
+    void SetLinearSolver(double rtol=1e-8, double atol=1e-12, int miter=1000);
+
+    /// Solves the forward problem.
+    void FSolve();
+
+    /// Solves the adjoint with the provided rhs.
+    void ASolve(mfem::Vector& rhs);
+
+    /// Adds displacement BC in direction 0(x),1(y),2(z), or 4(all).
+    void AddDispBC(int id, int dir, double val);
+
+    /// Adds displacement BC in direction 0(x),1(y),2(z), or 4(all).
+    void AddDispBC(int id, int dir, mfem::Coefficient& val);
+
+    /// Adds displacement BC specified by the vector coefficient val.
+    void AddDispBC(int id, mfem::VectorCoefficient& val);
+
+    /// Adds vol force
+    void AddVolForce(int id, double fx, double fy, double fz);
+
+    /// Adds vol force
+    void AddVolForce(int id, mfem::VectorCoefficient& ff);
+
+    /// Returns the displacements.
+    mfem::ParGridFunction& GetDisplacements()
+    {
+        fdisp.SetFromTrueDofs(sol.GetBlock(0));
+        return fdisp;
+    }
+
+    /// Returns the adjoint displacements.
+    mfem::ParGridFunction& GetADisplacements()
+    {
+        adisp.SetFromTrueDofs(adj.GetBlock(0));
+        return adisp;
+    }
+
+    /// Returns the solution vector.
+    mfem::Vector& GetSol(){return sol;}
+
+    /// Returns the adjoint solution vector.
+    mfem::Vector& GetAdj(){return adj;}
+
+    void GetSol(ParGridFunction& sgf){
+        sgf.SetSpace(vfes); sgf.SetFromTrueDofs(sol.GetBlock(0));}
+
+    void GetAdj(ParGridFunction& agf){
+        agf.SetSpace(vfes); agf.SetFromTrueDofs(adj.GetBlock(0));}
+private:
+
+    double current_time;
+
+    mfem::ParMesh* pmesh;
+
+    //solution vector
+    mfem::BlockVector sol;
+    //adjoint vector
+    mfem::BlockVector adj;
+    //RHS
+    mfem::BlockVector rhs;
+
+    // localy defined volumetric forces
+    std::map<int, mfem::VectorConstantCoefficient*> lvforce;
+    // globaly defined volumetric forces
+    std::map<int, mfem::VectorCoefficient*> volforce;
+
+    // boundary conditions for x,y, and z directions
+    std::map<int, mfem::ConstantCoefficient> bcx;
+    std::map<int, mfem::ConstantCoefficient> bcy;
+    std::map<int, mfem::ConstantCoefficient> bcz;
+
+    // holds BC in coefficient form
+    std::map<int, mfem::Coefficient*> bccx;
+    std::map<int, mfem::Coefficient*> bccy;
+    std::map<int, mfem::Coefficient*> bccz;
+    std::map<int, mfem::VectorCoefficient*> bcca;
+
+    // holds the displacement contrained DOFs
+    mfem::Array<int> ess_tdofv;
+
+    //forward solution
+    mfem::ParGridFunction fdisp;
+    //adjoint solution
+    mfem::ParGridFunction adisp;
+
+    //total strains
+    mfem::QuadratureFunction eee;
+    //plastic strains
+    mfem::QuadratureFunction eep;
+    //accumulated plastic strain
+    mfem::QuadratureFunction kappa;
+
+    //Newton solver parameters
+    double abs_tol;
+    double rel_tol;
+    int print_level;
+    int max_iter;
+
+    //Linear solver parameters
+    double linear_rtol;
+    double linear_atol;
+    int linear_iter;
+
+    mfem::ParBlockNonlinearForm *nf;
+    mfem::ParFiniteElementSpace* vfes; //displacements fes
+    mfem::ParFiniteElementSpace* ffes; //filter fes
+
+    mfem::FiniteElementCollection* vfec;
+    mfem::FiniteElementCollection* ffec;
+
+    mfem::QuadratureSpace* qfes;
+};
+
+}
+
+#endif

miniapps/elplast/holep.cpp

@@ -0,0 +1,401 @@
+//                       MFEM Example 21 - Parallel Version
+//
+// Compile with: make ex21p
+//
+// Sample runs:  mpirun -np 4 ex21p
+//               mpirun -np 4 ex21p -o 3
+//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh
+//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh -o 2
+//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh -o 2
+//
+// Description:  This is a version of Example 2p with a simple adaptive mesh
+//               refinement loop. The problem being solved is again the linear
+//               elasticity describing a multi-material cantilever beam.
+//               The problem is solved on a sequence of meshes which
+//               are locally refined in a conforming (triangles, tetrahedrons)
+//               or non-conforming (quadrilaterals, hexahedra) manner according
+//               to a simple ZZ error estimator.
+//
+//               The example demonstrates MFEM's capability to work with both
+//               conforming and nonconforming refinements, in 2D and 3D, on
+//               linear and curved meshes. Interpolation of functions from
+//               coarse to fine meshes, as well as persistent GLVis
+//               visualization are also illustrated.
+//
+//               We recommend viewing Examples 2p and 6p before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "hole-abaqus.mesh";
+   int serial_ref_levels = 0;
+   int order = 2;
+   //bool static_cond = false;
+   bool static_cond = true;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&serial_ref_levels, "-rs", "--refine-serial",
+                  "Number of uniform serial refinements (before parallel"
+                  " partitioning)");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, and hexahedral meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
+
+   std::cout<<"NE="<<mesh.GetNE()<<std::endl;
+   std::cout<<"NN="<<mesh.GetNV()<<std::endl;
+
+   double* v1=mesh.GetVertex(0);
+   std::cout<<v1[0]<<" "<<v1[1]<<" "<<v1[2]<<std::endl;
+
+
+   // 3. Refine the mesh before parallel partitioning. Since a NURBS mesh can
+   //    currently only be refined uniformly, we need to convert it to a
+   //    piecewise-polynomial curved mesh. First we refine the NURBS mesh a bit
+   //    more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext && serial_ref_levels == 0)
+   {
+      serial_ref_levels = 2;
+   }
+   for (int i = 0; i < serial_ref_levels; i++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(2);
+   }
+   mesh.EnsureNCMesh();
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // 4. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec, dim);
+
+   // 5. As in Example 2, we set up the linear form b(.) which corresponds to
+   //    the right-hand side of the FEM linear system. In this case, b_i equals
+   //    the boundary integral of f*phi_i where f represents a "pull down"
+   //    force on the Neumann part of the boundary and phi_i are the basis
+   //    functions in the finite element fespace. The force is defined by the
+   //    VectorArrayCoefficient object f, which is a vector of Coefficient
+   //    objects. The fact that f is non-zero on boundary attribute 2 is
+   //    indicated by the use of piece-wise constants coefficient for its last
+   //    component. We don't assemble the discrete problem yet, this will be
+   //    done in the main loop.
+   VectorArrayCoefficient f(dim);
+   for (int i = 0; i < dim-1; i++)
+   {
+      f.Set(i, new ConstantCoefficient(0.0));
+   }
+
+      Vector pull_force(pmesh.attributes.Max());
+      pull_force = 0.0;
+   // pull_force(1) = -1.0e-2;
+      pull_force(1) = 0.757576E-2;
+      pull_force(2) = 0.757576E-2;
+
+      pull_force*=100000;
+
+      f.Set(0, new PWConstCoefficient(pull_force));
+
+   ParLinearForm b(&fespace);
+
+   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f));
+
+
+   // 6. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the linear elasticity integrator with piece-wise
+   //    constants coefficient lambda and mu.
+
+   std::cout<<"pmesh atr="<<pmesh.attributes.Max()<<std::endl;
+
+   Vector lambda(pmesh.attributes.Max());
+   lambda = 1.0;
+   //lambda(0) = lambda(1)*50;
+   lambda(0) = 0.138651E+8;
+   lambda(1) = 0.576923E+7;
+   lambda(2) = 0.288462E+6;
+
+   PWConstCoefficient lambda_func(lambda);
+   Vector mu(pmesh.attributes.Max());
+   mu = 1.0;
+   //mu(0) = mu(1)*50;
+   mu(0) = 0.11811E+8;
+   mu(1) = 0.384615E+7;
+   mu(2) = 0.192308E+6;
+
+   PWConstCoefficient mu_func(mu);
+
+   ParBilinearForm a(&fespace);
+   BilinearFormIntegrator *integ =
+      new ElasticityIntegrator(lambda_func,mu_func);
+   a.AddDomainIntegrator(integ);
+   if (static_cond) { a.EnableStaticCondensation(); }
+
+   // 7. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   Vector zero_vec(dim);
+   zero_vec = 0.0;
+   VectorConstantCoefficient zero_vec_coeff(zero_vec);
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking only
+   //    boundary attribute 1 from the mesh as essential and converting it to a
+   //    list of true dofs.  The conversion to true dofs will be done in the
+   //    main loop.
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 0;
+   ess_bdr[0] = 1;
+
+   // 9. GLVis visualization.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+
+   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     that uses the ComputeElementFlux method of the ElasticityIntegrator to
+   //     recover a smoothed flux (stress) that is subtracted from the element
+   //     flux to get an error indicator. We need to supply the space for the
+   //     smoothed flux: an (H1)^tdim (i.e., vector-valued) space is used here.
+   //     Here, tdim represents the number of components for a symmetric (dim x
+   //     dim) tensor.
+   const int tdim = dim*(dim+1)/2;
+   L2_FECollection flux_fec(order, dim);
+   ParFiniteElementSpace flux_fespace(&pmesh, &flux_fec, tdim);
+   ParFiniteElementSpace smooth_flux_fespace(&pmesh, &fec, tdim);
+   L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace,
+                                       smooth_flux_fespace);
+
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaView", &pmesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.SetDataFormat(VTKFormat::ASCII);
+   dacol.RegisterField("disp", &x);
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+
+   // 11. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.7);
+
+   // 12. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   const int max_dofs = 500000;
+   const int max_amr_itr = 20;
+   for (int it = 0; it <= max_amr_itr; it++)
+   {
+      HYPRE_BigInt global_dofs = fespace.GlobalTrueVSize();
+      if (myid == 0)
+      {
+         cout << "\nAMR iteration " << it << endl;
+         cout << "Number of unknowns: " << global_dofs << endl;
+      }
+
+      // 13. Assemble the stiffness matrix and the right-hand side.
+      a.Assemble();
+      b.Assemble();
+
+      // 14. Set Dirichlet boundary values in the GridFunction x.
+      //     Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_tdof_list;
+      x.ProjectBdrCoefficient(zero_vec_coeff, ess_bdr);
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      // 15. Create the linear system: eliminate boundary conditions, constrain
+      //     hanging nodes and possibly apply other transformations. The system
+      //     will be solved for true (unconstrained) DOFs only.
+
+      HypreParMatrix A;
+      Vector B, X;
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+
+      // 16. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
+      //     preconditioner from hypre.
+      HypreBoomerAMG amg;
+      amg.SetPrintLevel(0);
+      // amg.SetSystemsOptions(dim); // optional
+      CGSolver pcg(A.GetComm());
+      pcg.SetPreconditioner(amg);
+      pcg.SetOperator(A);
+      pcg.SetRelTol(1e-6);
+      pcg.SetMaxIter(5000);
+      pcg.SetPrintLevel(3); // print the first and the last iterations only
+      pcg.Mult(B, X);
+
+      // 17. After solving the linear system, reconstruct the solution as a
+      //     finite element GridFunction. Constrained nodes are interpolated
+      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
+      a.RecoverFEMSolution(X, b, x);
+
+      //ParaView output
+      dacol.SetTime(double(it+1));
+      dacol.SetCycle(it+1);
+      dacol.Save();
+
+
+      // 18. Send solution by socket to the GLVis server.
+      if (visualization && it == 0)
+      {
+         sol_sock.open(vishost, visport);
+         sol_sock.precision(8);
+      }
+      if (visualization && sol_sock.good())
+      {
+         GridFunction nodes(&fespace), *nodes_p = &nodes;
+         pmesh.GetNodes(nodes);
+         nodes += x;
+         int own_nodes = 0;
+         pmesh.SwapNodes(nodes_p, own_nodes);
+         x.Neg(); // visualize the backward displacement
+         sol_sock << "parallel " << num_procs << ' ' << myid << '\n';
+         sol_sock << "solution\n" << pmesh << x << flush;
+         x.Neg();
+         pmesh.SwapNodes(nodes_p, own_nodes);
+         if (it == 0)
+         {
+            sol_sock << "keys '" << ((dim == 2) ? "Rjl" : "") << "m'" << endl;
+         }
+         sol_sock << "window_title 'AMR iteration: " << it << "'\n"
+                  << "pause" << endl;
+         if (myid == 0)
+         {
+            cout << "Visualization paused. "
+                 "Press <space> in the GLVis window to continue." << endl;
+         }
+      }
+
+      if (global_dofs > max_dofs)
+      {
+         if (myid == 0)
+         {
+            cout << "Reached the maximum number of dofs. Stop." << endl;
+         }
+         break;
+      }
+
+      // 19. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+      refiner.Apply(pmesh);
+      if (refiner.Stop())
+      {
+         if (myid == 0)
+         {
+            cout << "Stopping criterion satisfied. Stop." << endl;
+         }
+         break;
+      }
+
+      // 20. Update the space to reflect the new state of the mesh. Also,
+      //     interpolate the solution x so that it lies in the new space but
+      //     represents the same function. This saves solver iterations later
+      //     since we'll have a good initial guess of x in the next step.
+      //     Internally, FiniteElementSpace::Update() calculates an
+      //     interpolation matrix which is then used by GridFunction::Update().
+      fespace.Update();
+      x.Update();
+
+      // 21. Load balance the mesh, and update the space and solution. Currently
+      //     available only for nonconforming meshes.
+      if (pmesh.Nonconforming())
+      {
+         pmesh.Rebalance();
+
+         // Update the space and the GridFunction. This time the update matrix
+         // redistributes the GridFunction among the processors.
+         fespace.Update();
+         x.Update();
+      }
+
+      // 21. Inform also the bilinear and linear forms that the space has
+      //     changed.
+      a.Update();
+      b.Update();
+   }
+
+   dacol.SetTime(double(max_amr_itr));
+   dacol.SetCycle(max_amr_itr);
+   dacol.Save();
+
+   {
+      ostringstream mref_name, mesh_name, sol_name;
+      mref_name << "ex21p_reference_mesh." << setfill('0') << setw(6) << myid;
+      mesh_name << "ex21p_deformed_mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "ex21p_displacement." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ref_out(mref_name.str().c_str());
+      mesh_ref_out.precision(16);
+      pmesh.Print(mesh_ref_out);
+
+      ofstream mesh_out(mesh_name.str().c_str());
+      mesh_out.precision(16);
+      GridFunction nodes(&fespace), *nodes_p = &nodes;
+      pmesh.GetNodes(nodes);
+      nodes += x;
+      int own_nodes = 0;
+      pmesh.SwapNodes(nodes_p, own_nodes);
+      pmesh.Print(mesh_out);
+      pmesh.SwapNodes(nodes_p, own_nodes);
+
+      ofstream x_out(sol_name.str().c_str());
+      x_out.precision(16);
+      x.Save(x_out);
+   }
+
+   return 0;
+}

miniapps/elplast/makefile

@@ -0,0 +1,89 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/elplast/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+ADAPT_SRC = adaptive_el.cpp
+ADAPT_OBJ = $(ADAPT_SRC:.cpp=.o)
+
+PAR_MINIAPPS = adaptel
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = 
+else
+   MINIAPPS = $(PAR_MINIAPPS)
+endif
+
+COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
+
+# If MFEM_SHARED is set, add the ../common rpath
+COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all lib-common clean clean-build clean-exec
+
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+%: %.o $(NAVIER_COMMON_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $^ -o $@ $(MFEM_LIBS)
+
+%.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+
+
+
+all: $(MINIAPPS)
+
+adaptel: $(ADAPT_OBJ) lib-common
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(ADAPT_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
+# Rule for building lib-common
+lib-common:
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ adaptel
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -f diffusion.mesh diffusion.gf
+	@rm -rf ParaViewDistance ParaViewDiffusion ParaViewExtrapolate ParaViewLSF

miniapps/elplast/test_elinteg.cpp

@@ -0,0 +1,148 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "coefficients.hpp"
+
+
+int main(int argc, char *argv[])
+{
+
+    mfem::ElementTransformation* Tr;
+    mfem::IntegrationPoint* ip;
+
+
+    mfem::IsoElastMat mat; mat.SetElastParam(1,0.0);
+    mfem::J2YieldFunction yf(1.0,0.10,1.0);
+
+    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> se(&mat,&yf);
+
+    mfem::Vector ss(9); ss=0.0;
+    mfem::Vector ep(9); ep=0.0;
+    mfem::Vector iv(2); iv=0.0;
+
+    mfem::Vector dee(9); dee=0.0;
+    mfem::Vector ee(9); ee=0.0;
+    mfem::Vector epn(9); epn=0.0;
+    mfem::Vector ivn(2); ivn=0.0;
+
+    int bb=8;
+
+    for(int i=0;i<20;i++){
+
+
+        std::cout<<i<<" " ;//std::endl;
+
+        dee(bb)=0.11;
+        ee.Add(1.0,dee);
+        se.SetStrain(ee);
+        se.SetPlasticStrain(epn);
+        se.SetInternalParameters(ivn);
+        se.Solve(ss,epn,ivn,*Tr,*ip);
+        //epn=ep;
+        //ivn=iv;
+
+        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
+
+        /*
+        std::cout<<"strain="<<std::endl;
+        ee.Print(std::cout,3);
+
+        std::cout<<"stress="<<std::endl;
+        std::cout.precision(6);
+        ss.Print(std::cout,3);
+
+        iv.Print(std::cout);
+        */
+
+    }
+
+    for(int i=0;i<40;i++){
+
+
+        std::cout<<i<<" " ;//std::endl;
+
+        dee(bb)=-0.11;
+        ee.Add(1.0,dee);
+        se.SetStrain(ee);
+        se.SetPlasticStrain(epn);
+        se.SetInternalParameters(ivn);
+        se.Solve(ss,ep,iv,*Tr,*ip);
+        epn=ep;
+        ivn=iv;
+
+        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
+
+        /*
+        std::cout<<"strain="<<std::endl;
+        ee.Print(std::cout,3);
+
+        std::cout<<"stress="<<std::endl;
+        std::cout.precision(6);
+        ss.Print(std::cout,3);
+
+        iv.Print(std::cout);
+        */
+
+    }
+
+    for(int i=0;i<25;i++){
+
+
+        std::cout<<i<<" " ;//std::endl;
+
+        dee(bb)=0.11;
+        ee.Add(1.0,dee);
+        se.SetStrain(ee);
+        se.SetPlasticStrain(epn);
+        se.SetInternalParameters(ivn);
+        se.Solve(ss,ep,iv,*Tr,*ip);
+        epn=ep;
+        ivn=iv;
+
+        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
+
+        /*
+        std::cout<<"strain="<<std::endl;
+        ee.Print(std::cout,3);
+
+        std::cout<<"stress="<<std::endl;
+        std::cout.precision(6);
+        ss.Print(std::cout,3);
+
+        iv.Print(std::cout);
+        */
+
+    }
+
+    for(int i=0;i<15;i++){
+
+
+        std::cout<<i<<" " ;//std::endl;
+
+        dee(bb)=-0.11;
+        ee.Add(1.0,dee);
+        se.SetStrain(ee);
+        se.SetPlasticStrain(epn);
+        se.SetInternalParameters(ivn);
+        se.Solve(ss,ep,iv,*Tr,*ip);
+        epn=ep;
+        ivn=iv;
+
+        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
+
+        /*
+        std::cout<<"strain="<<std::endl;
+        ee.Print(std::cout,3);
+
+        std::cout<<"stress="<<std::endl;
+        std::cout.precision(6);
+        ss.Print(std::cout,3);
+
+        iv.Print(std::cout);
+        */
+
+    }
+
+
+    return 0;
+}