The source file nonlinear_neohooke_cc.h with comments:

        #ifndef NONLINEAR_NEOHOOKE_CC_H_
        #define NONLINEAR_NEOHOOKE_CC_H_
        
        #include "libmesh/dense_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/vector_value.h"
        #include "libmesh/tensor_value.h"
        #include "libmesh/getpot.h"
        
        using namespace libMesh;
        
        /**
         * This class implements a constitutive formulation for an Neo-Hookean elastic solid
         * in terms of the current configuration. This implementation is suitable for computing
         * large deformations. See e.g. "Nonlinear finite element methods" (P. Wriggers, 2008,
         * Springer) for details.
         */
        class NonlinearNeoHookeCurrentConfig {
        public:
          NonlinearNeoHookeCurrentConfig(
              const std::vector<std::vector<RealGradient> >& dphi_in, GetPot& args) :
              dphi(dphi_in) {
            E = args("material/neohooke/e_modulus", 10000.0);
            nu = args("material/neohooke/nu", 0.3);
          }
        
          /**
           * Initialize the class for the given displacement gradient at the
           * specified quadrature point.
           */
          void init_for_qp(VectorValue<Gradient> & grad_u, unsigned int qp);
        
          /**
           * Return the residual vector for the current state.
           */
          void get_residual(DenseVector<Real> & residuum, unsigned int & i);
        
          /**
           * Return the stiffness matrix for the current state.
           */
          void get_linearized_stiffness(DenseMatrix<Real> & stiffness,
              unsigned int & i, unsigned int & j);
        
          /**
           * Flag to indicate if it is necessary to calculate values for stiffness
           * matrix during initialization.
           */
          bool calculate_linearized_stiffness;
        private:
          void build_b_0_mat(int i, DenseMatrix<Real>& b_l_mat);
          void calculate_stress();
          void calculate_tangent();
          static void tensor_to_voigt(const RealTensor &tensor, DenseVector<Real> &vec);
        
          unsigned int current_qp;
          const std::vector<std::vector<RealGradient> >& dphi;
        
          DenseMatrix<Real> C_mat;
          Real E;
          Real nu;
          RealTensor F, S, tau, sigma;
          DenseMatrix<Real> B_L;
          DenseMatrix<Real> B_K;
        };
        
        #endif /* NONLINEAR_NEOHOOKE_CC_H_ */



The source file solid_system.h with comments:

        #ifndef SOLID_SYSTEM_H_
        #define SOLID_SYSTEM_H_
        
        #include "libmesh/fem_system.h"
        #include "libmesh/auto_ptr.h"
        
        using namespace libMesh;
        
        class SolidSystem: public FEMSystem {
        public:
Constructor
          SolidSystem(EquationSystems& es, const std::string& name,
              const unsigned int number);
        
System initialization
          virtual void init_data();
        
Context initialization
          virtual void init_context(DiffContext &context);
        
Element residual and jacobian calculations
          virtual bool element_time_derivative(bool request_jacobian,
              DiffContext& context);
        
Contributions for adding boundary conditions
          virtual bool side_time_derivative(bool request_jacobian,
              DiffContext& context);
        
          virtual bool eulerian_residual(bool, DiffContext &) {
            return false;
          }
        
Simulation parameters
          GetPot args;
        
Custom Identifier
          virtual std::string system_type() const {
            return "Solid";
          }
        
override method to update mesh also
          virtual void update();
        
save the undeformed mesh to an auxiliary system
          void save_initial_mesh();
        
variable numbers of primary variables in the current system
          unsigned int var[3];
        
variable numbers of primary variables in auxiliary system (for accessing the undeformed configuration)
          unsigned int undefo_var[3];
        };
        
        #endif /* SOLID_SYSTEM_H_ */



The source file fem_system_ex2.C with comments:

        #include "libmesh/equation_systems.h"
        #include "libmesh/getpot.h"
        #include "libmesh/libmesh.h"
        #include "libmesh/libmesh_logging.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/string_to_enum.h"
        #include "libmesh/time_solver.h"
        #include "libmesh/transient_system.h"
        #include "libmesh/vtk_io.h"
        
        #include <cstdio>
        #include <ctime>
        #include <fstream>
        #include <iomanip>
        #include <iostream>
        #include <sstream>
        #include <string>
        #include <unistd.h>
        
        using namespace libMesh;
        
        #include "solid_system.h"
        
        void setup(EquationSystems& systems, Mesh& mesh, GetPot& args)
        {
          const unsigned int dim = mesh.mesh_dimension();
We currently invert tensors with the assumption that they're 3x3
          libmesh_assert (dim == 3);
        
Generating Mesh
          ElemType eltype = Utility::string_to_enum<ElemType>(args("mesh/generation/element_type", "hex8"));
          int nx = args("mesh/generation/num_elem", 4, 0);
          int ny = args("mesh/generation/num_elem", 4, 1);
          int nz = dim > 2 ? args("mesh/generation/num_elem", 4, 2) : 0;
          double origx = args("mesh/generation/origin", -1.0, 0);
          double origy = args("mesh/generation/origin", -1.0, 1);
          double origz = args("mesh/generation/origin", 0.0, 2);
          double sizex = args("mesh/generation/size", 2.0, 0);
          double sizey = args("mesh/generation/size", 2.0, 1);
          double sizez = args("mesh/generation/size", 2.0, 2);
          MeshTools::Generation::build_cube(mesh, nx, ny, nz,
              origx, origx+sizex, origy, origy+sizey, origz, origz+sizez, eltype);
        
Creating Systems
          SolidSystem& imms = systems.add_system<SolidSystem> ("solid");
          imms.args = args;
        
Build up auxiliary system
          ExplicitSystem& aux_sys = systems.add_system<TransientExplicitSystem>("auxiliary");
        
Initialize the system
          systems.parameters.set<unsigned int>("phase") = 0;
          systems.init();
        
          imms.save_initial_mesh();
        
Fill global solution vector from local ones
          aux_sys.reinit();
        }
        
        
        
        void run_timestepping(EquationSystems& systems, GetPot& args)
        {
          TransientExplicitSystem& aux_system = systems.get_system<TransientExplicitSystem>("auxiliary");
        
          SolidSystem& solid_system = systems.get_system<SolidSystem>("solid");
        
          AutoPtr<VTKIO> io = AutoPtr<VTKIO>(new VTKIO(systems.get_mesh()));
        
          Real duration = args("duration", 1.0);
        
          for (unsigned int t_step = 0; t_step < duration/solid_system.deltat; t_step++) {
Progress in current phase [0..1]
            Real progress = t_step * solid_system.deltat / duration;
            systems.parameters.set<Real>("progress") = progress;
            systems.parameters.set<unsigned int>("step") = t_step;
        
Update message

            out << "===== Time Step " << std::setw(4) << t_step;
            out << " (" << std::fixed << std::setprecision(2) << std::setw(6) << (progress*100.) << "%)";
            out << ", time = " << std::setw(7) << solid_system.time;
            out << " =====" << std::endl;
        
Advance timestep in auxiliary system
            aux_system.current_local_solution->close();
            aux_system.old_local_solution->close();
            *aux_system.older_local_solution = *aux_system.old_local_solution;
            *aux_system.old_local_solution = *aux_system.current_local_solution;
        
            out << "Solving Solid" << std::endl;
            solid_system.solve();
            aux_system.reinit();
        
Carry out the adaptive mesh refinement/coarsening
            out << "Doing a reinit of the equation systems" << std::endl;
            systems.reinit();
        
            if (t_step % args("output/frequency", 1) == 0) {
              std::string result;
              std::stringstream file_name;
              file_name << args("results_directory", "./") << "fem_";
              file_name << std::setw(6) << std::setfill('0') << t_step;
              file_name << ".pvtu";
        
        
              io->write_equation_systems(file_name.str(), systems);
            }
Advance to the next timestep in a transient problem
            out << "Advancing to next step" << std::endl;
            solid_system.time_solver->advance_timestep();
          }
        }
        
        
        
        int main(int argc, char** argv)
        {
Initialize libMesh and any dependent libraries
          LibMeshInit init(argc, argv);
        
Skip this example if we do not meet certain requirements
        #ifndef LIBMESH_HAVE_VTK
          libmesh_example_assert(false, "--enable-vtk");
        #endif
        
Trilinos gives us an inverted element on this one...
          libmesh_example_assert(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc");
        
Threaded assembly doesn't currently work with the moving mesh code. We'll skip this example for now.
          if (libMesh::n_threads() > 1)
            {
              std::cout << "We skip fem_system_ex2 when using threads.\n"
                        << std::endl;
              return 0;
            }
        
read simulation parameters from file
          GetPot args = GetPot("solid.in");
        
Create System and Mesh
          int dim = args("mesh/generation/dimension", 3);
          libmesh_example_assert(dim <= LIBMESH_DIM, "3D support");
        
Create a mesh distributed across the default MPI communicator.
          Mesh mesh(init.comm(), dim);
        
          EquationSystems systems(mesh);
        
Create and set systems up
          setup(systems, mesh, args);
        
run the systems
          run_timestepping(systems, args);
        
          out << "Finished calculations" << std::endl;
          return 0;
        }



The source file nonlinear_neohooke_cc.C with comments:

        #include "nonlinear_neohooke_cc.h"
        
        /**
         * Return the inverse of the given TypeTensor. Algorithm taken from the tensor classes
         * of Ton van den Boogaard (http://tmku209.ctw.utwente.nl/~ton/tensor.html)
         */
        template <typename T> TypeTensor<T> inv(const TypeTensor<T> &A ) {
          double Sub11, Sub12, Sub13;
          Sub11 = A._coords[4]*A._coords[8] - A._coords[5]*A._coords[7];
          Sub12 = A._coords[3]*A._coords[8] - A._coords[6]*A._coords[5];
          Sub13 = A._coords[3]*A._coords[7] - A._coords[6]*A._coords[4];
          double detA = A._coords[0]*Sub11 - A._coords[1]*Sub12 + A._coords[2]*Sub13;
          libmesh_assert( std::fabs(detA)>1.e-15 );
        
          TypeTensor<T> Ainv(A);
        
          Ainv._coords[0] =  Sub11/detA;
          Ainv._coords[1] = (-A._coords[1]*A._coords[8]+A._coords[2]*A._coords[7])/detA;
          Ainv._coords[2] = ( A._coords[1]*A._coords[5]-A._coords[2]*A._coords[4])/detA;
          Ainv._coords[3] = -Sub12/detA;
          Ainv._coords[4] = ( A._coords[0]*A._coords[8]-A._coords[2]*A._coords[6])/detA;
          Ainv._coords[5] = (-A._coords[0]*A._coords[5]+A._coords[2]*A._coords[3])/detA;
          Ainv._coords[6] =  Sub13/detA;
          Ainv._coords[7] = (-A._coords[0]*A._coords[7]+A._coords[1]*A._coords[6])/detA;
          Ainv._coords[8] = ( A._coords[0]*A._coords[4]-A._coords[1]*A._coords[3])/detA;
        
          return Ainv;
        }
        
        void NonlinearNeoHookeCurrentConfig::init_for_qp(VectorValue<Gradient> & grad_u, unsigned int qp) {
        	this->current_qp = qp;
        	F.zero();
        	S.zero();
        
        	{
        	  RealTensor invF;
        	  invF.zero();
        	  for (unsigned int i = 0; i < 3; ++i)
        	    for (unsigned int j = 0; j < 3; ++j) {
        	      invF(i, j) += libmesh_real(grad_u(i)(j));
        	    }
        	  F.add(inv(invF));
        
        	  libmesh_assert_greater (F.det(), -TOLERANCE);
        	}
        
        	if (this->calculate_linearized_stiffness) {
        		this->calculate_tangent();
        	}
        
        	this->calculate_stress();
        }
        
        
        
        void NonlinearNeoHookeCurrentConfig::calculate_tangent() {
        	Real mu = E / (2 * (1 + nu));
        	Real lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
        
        	Real detF = F.det();
        
        	C_mat.resize(6, 6);
        	for (unsigned int i = 0; i < 3; ++i) {
        		for (unsigned int j = 0; j < 3; ++j) {
        			if (i == j) {
        				C_mat(i, j) = 2 * mu + lambda;
        				C_mat(i + 3, j + 3) = mu - 0.5 * lambda * (detF * detF - 1);
        			} else {
        				C_mat(i, j) = lambda * detF * detF;
        			}
        		}
        	}
        }
        
        
        void NonlinearNeoHookeCurrentConfig::calculate_stress() {
        
          double mu = E / (2.0 * (1.0 + nu));
        	double lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
        
        	Real detF = F.det();
        	RealTensor Ft = F.transpose();
        
        	RealTensor C = Ft * F;
        	RealTensor b = F * Ft;
        	RealTensor identity;
        	identity(0, 0) = 1.0; identity(1, 1) = 1.0; identity(2, 2) = 1.0;
        	RealTensor invC = inv(C);
        
        	S = 0.5 * lambda * (detF * detF - 1) * invC + mu * (identity - invC);
        
        	tau = (F * S) * Ft;
        	sigma = 1/detF * tau;
        }
        
        void NonlinearNeoHookeCurrentConfig::get_residual(DenseVector<Real> & residuum, unsigned int & i) {
        	B_L.resize(3, 6);
        	DenseVector<Real> sigma_voigt(6);
        
        	this->build_b_0_mat(i, B_L);
        
        	tensor_to_voigt(sigma, sigma_voigt);
        
        	B_L.vector_mult(residuum, sigma_voigt);
        }
        
        void NonlinearNeoHookeCurrentConfig::tensor_to_voigt(const RealTensor &tensor, DenseVector<Real> &vec) {
          vec(0) = tensor(0, 0);
          vec(1) = tensor(1, 1);
          vec(2) = tensor(2, 2);
          vec(3) = tensor(0, 1);
          vec(4) = tensor(1, 2);
          vec(5) = tensor(0, 2);
        
        }
        
        void NonlinearNeoHookeCurrentConfig::get_linearized_stiffness(DenseMatrix<Real> & stiffness, unsigned int & i, unsigned int & j) {
        	stiffness.resize(3, 3);
        
        	double G_IK = (sigma * dphi[i][current_qp]) * dphi[j][current_qp];
        	stiffness(0, 0) += G_IK;
        	stiffness(1, 1) += G_IK;
        	stiffness(2, 2) += G_IK;
        
        	B_L.resize(3, 6);
        	this->build_b_0_mat(i, B_L);
        	B_K.resize(3, 6);
        	this->build_b_0_mat(j, B_K);
        
        	B_L.right_multiply(C_mat);
        	B_L.right_multiply_transpose(B_K);
        	B_L *= 1/F.det();
        
        	stiffness += B_L;
        }
        
        void NonlinearNeoHookeCurrentConfig::build_b_0_mat(int i, DenseMatrix<Real>& b_0_mat) {
        	for (unsigned int ii = 0; ii < 3; ++ii) {
        		b_0_mat(ii, ii) = dphi[i][current_qp](ii);
        	}
        	b_0_mat(0, 3) = dphi[i][current_qp](1);
        	b_0_mat(1, 3) = dphi[i][current_qp](0);
        	b_0_mat(1, 4) = dphi[i][current_qp](2);
        	b_0_mat(2, 4) = dphi[i][current_qp](1);
        	b_0_mat(0, 5) = dphi[i][current_qp](2);
        	b_0_mat(2, 5) = dphi[i][current_qp](0);
        }



The source file solid_system.C with comments:

        #include "libmesh/boundary_info.h"
        #include "libmesh/diff_solver.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/getpot.h"
        #include "libmesh/mesh.h"
        #include "libmesh/newton_solver.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/quadrature.h"
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/steady_solver.h"
        #include "libmesh/transient_system.h"
        
        #include "nonlinear_neohooke_cc.h"
        #include "solid_system.h"
        
Solaris Studio has no NAN
        #ifdef __SUNPRO_CC
          #define NAN (1.0/0.0)
        #endif
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        SolidSystem::SolidSystem(EquationSystems& es, const std::string& name_in,
            const unsigned int number_in) :
            FEMSystem(es, name_in, number_in) {
        
Add a time solver. We are just looking at a steady state problem here.
          this->time_solver = AutoPtr<TimeSolver>(new SteadySolver(*this));
        }
        
        void SolidSystem::save_initial_mesh() {
          System & aux_sys = this->get_equation_systems().get_system("auxiliary");
        
          const unsigned int dim = this->get_mesh().mesh_dimension();
        
Loop over all nodes and copy the location from the current system to the auxiliary system.
          const MeshBase::const_node_iterator nd_end =
              this->get_mesh().local_nodes_end();
          for (MeshBase::const_node_iterator nd = this->get_mesh().local_nodes_begin();
              nd != nd_end; ++nd) {
            const Node *node = *nd;
            for (unsigned int d = 0; d < dim; ++d) {
              unsigned int source_dof = node->dof_number(this->number(), var[d], 0);
              unsigned int dest_dof = node->dof_number(aux_sys.number(), undefo_var[d],
                  0);
              Number value = this->current_local_solution->el(source_dof);
              aux_sys.current_local_solution->set(dest_dof, value);
            }
          }
        }
        
        void SolidSystem::init_data() {
          const unsigned int dim = this->get_mesh().mesh_dimension();
        
Get the default order of the used elements. Assumption: Just one type of elements in the mesh.
          Order order = (*(this->get_mesh().elements_begin()))->default_order();
        
Add the node positions as primary variables.
          var[0] = this->add_variable("x", order);
          var[1] = this->add_variable("y", order);
          if (dim == 3)
            var[2] = this->add_variable("z", order);
          else
            var[2] = var[1];
        
Add variables for storing the initial mesh to an auxiliary system.
          System& aux_sys = this->get_equation_systems().get_system("auxiliary");
          undefo_var[0] = aux_sys.add_variable("undefo_x", order);
          undefo_var[1] = aux_sys.add_variable("undefo_y", order);
          undefo_var[2] = aux_sys.add_variable("undefo_z", order);
        
Set the time stepping options
          this->deltat = args("schedule/dt", 0.2);
        
Do the parent's initialization after variables are defined
          FEMSystem::init_data();
        
// Tell the system to march velocity forward in time, but // leave p as a constraint only this->time_evolving(var[0]); this->time_evolving(var[1]); if (dim == 3) this->time_evolving(var[2]);

Tell the system which variables are containing the node positions
          set_mesh_system(this);
        
          this->set_mesh_x_var(var[0]);
          this->set_mesh_y_var(var[1]);
          if (dim == 3)
            this->set_mesh_z_var(var[2]);
        
Fill the variables with the position of the nodes
          this->mesh_position_get();
        
          System::reinit();
        
Set some options for the DiffSolver
          DiffSolver &solver = *(this->time_solver->diff_solver().get());
          solver.quiet = args("solver/quiet", false);
          solver.max_nonlinear_iterations = args(
              "solver/nonlinear/max_nonlinear_iterations", 100);
          solver.relative_step_tolerance = args(
              "solver/nonlinear/relative_step_tolerance", 1.e-3);
          solver.relative_residual_tolerance = args(
              "solver/nonlinear/relative_residual_tolerance", 1.e-8);
          solver.absolute_residual_tolerance = args(
              "solver/nonlinear/absolute_residual_tolerance", 1.e-8);
          solver.verbose = !args("solver/quiet", false);
        
          ((NewtonSolver&) solver).require_residual_reduction = args(
              "solver/nonlinear/require_reduction", false);
        
And the linear solver options
          solver.max_linear_iterations = args("max_linear_iterations", 50000);
          solver.initial_linear_tolerance = args("initial_linear_tolerance", 1.e-3);
        }
        
        void SolidSystem::update() {
          System::update();
          this->mesh_position_set();
        }
        
        void SolidSystem::init_context(DiffContext &context) {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
Pre-request all the data needed
          c.element_fe_var[var[0]]->get_JxW();
          c.element_fe_var[var[0]]->get_phi();
          c.element_fe_var[var[0]]->get_dphi();
          c.element_fe_var[var[0]]->get_xyz();
        
          c.side_fe_var[var[0]]->get_JxW();
          c.side_fe_var[var[0]]->get_phi();
          c.side_fe_var[var[0]]->get_xyz();
        }
        
        /**
         * Compute contribution from internal forces in elem_residual and contribution from
         * linearized internal forces (stiffness matrix) in elem_jacobian.
         */
        bool SolidSystem::element_time_derivative(bool request_jacobian,
            DiffContext &context) {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
First we get some references to cell-specific data that will be used to assemble the linear system.

Element Jacobian * quadrature weights for interior integration
          const std::vector<Real> &JxW = c.element_fe_var[var[0]]->get_JxW();
        
The gradients of the shape functions at interior quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi =
              c.element_fe_var[var[0]]->get_dphi();
        
Dimension of the mesh
          const unsigned int dim = this->get_mesh().mesh_dimension();
        
The number of local degrees of freedom in each variable
          const unsigned int n_u_dofs = c.dof_indices_var[var[0]].size();
          libmesh_assert(n_u_dofs == c.dof_indices_var[var[1]].size());
          if (dim == 3) {
            libmesh_assert(n_u_dofs == c.dof_indices_var[var[2]].size());
          }
        
          unsigned int n_qpoints = c.element_qrule->n_points();
        
Some matrices and vectors for storing the results of the constitutive law
          DenseMatrix<Real> stiff;
          DenseVector<Real> res;
          VectorValue<Gradient> grad_u;
        
Instantiate the constitutive law
          NonlinearNeoHookeCurrentConfig material(dphi, args);
        
Just calculate jacobian contribution when we need to
          material.calculate_linearized_stiffness = request_jacobian;
        
Get a reference to the auxiliary system
          TransientExplicitSystem& aux_system = this->get_equation_systems().get_system<
              TransientExplicitSystem>("auxiliary");
          std::vector<dof_id_type> undefo_index;
        
Assume symmetry of local stiffness matrices
          bool use_symmetry = args("assembly/use_symmetry", false);
        
Now we will build the element Jacobian and residual. This must be calculated at each quadrature point by summing the solution degree-of-freedom values by the appropriate weight functions. This class just takes care of the assembly. The matrix of the jacobian and the residual vector are provided by the constitutive formulation.

          for (unsigned int qp = 0; qp != n_qpoints; qp++) {
Compute the displacement gradient
            grad_u(0) = grad_u(1) = grad_u(2) = 0;
            for (unsigned int d = 0; d < dim; ++d) {
              std::vector<Number> u_undefo;
              aux_system.get_dof_map().dof_indices(c.elem, undefo_index, undefo_var[d]);
              aux_system.current_local_solution->get(undefo_index, u_undefo);
              for (unsigned int l = 0; l != n_u_dofs; l++)
                grad_u(d).add_scaled(dphi[l][qp], u_undefo[l]); // u_current(l)); // -
            }
        
initialize the constitutive formulation with the current displacement gradient
            material.init_for_qp(grad_u, qp);
        
Aquire, scale and assemble residual and stiffness
            for (unsigned int i = 0; i < n_u_dofs; i++) {
              res.resize(dim);
              material.get_residual(res, i);
              res.scale(JxW[qp]);
              for (unsigned int ii = 0; ii < dim; ++ii) {
                c.elem_subresiduals[ii]->operator ()(i) += res(ii);
              }
        
              if (request_jacobian && c.elem_solution_derivative) {
                libmesh_assert(c.elem_solution_derivative == 1.0);
                for (unsigned int j = (use_symmetry ? i : 0); j < n_u_dofs; j++) {
                  material.get_linearized_stiffness(stiff, i, j);
                  stiff.scale(JxW[qp]);
                  for (unsigned int ii = 0; ii < dim; ++ii) {
                    for (unsigned int jj = 0; jj < dim; ++jj) {
                      c.elem_subjacobians[ii][jj]->operator ()(i, j) += stiff(ii, jj);
                      if (use_symmetry && i != j) {
                        c.elem_subjacobians[ii][jj]->operator ()(j, i) += stiff(jj, ii);
                      }
                    }
                  }
                }
              }
            }
          } // end of the quadrature point qp-loop
        
          return request_jacobian;
        }
        
        bool SolidSystem::side_time_derivative(bool request_jacobian,
            DiffContext &context) {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
Apply displacement boundary conditions with penalty method

Get the current load step
          Real ratio = this->get_equation_systems().parameters.get<Real>("progress")
              + 0.001;
        
The BC are stored in the simulation parameters as array containing sequences of four numbers: Id of the side for the displacements and three values describing the displacement. E.g.: bc/displacement = '5 nan nan -1.0'. This will move all nodes of side 5 about 1.0 units down the z-axis while leaving all other directions unrestricted

Get number of BCs to enforce
          unsigned int num_bc = args.vector_variable_size("bc/displacement");
          if (num_bc % 4 != 0) {
            libMesh::err
                << "ERROR, Odd number of values in displacement boundary condition.\n"
                << std::endl;
            libmesh_error();
          }
          num_bc /= 4;
        
Loop over all BCs
          for (unsigned int nbc = 0; nbc < num_bc; nbc++) {
Get IDs of the side for this BC
            short int positive_boundary_id = args("bc/displacement", 1, nbc * 4);
        
The current side may not be on the boundary to be restricted
            if (!this->get_mesh().boundary_info->has_boundary_id
        	  (c.elem,c.side,positive_boundary_id))
              continue;
        
Read values from configuration file
            Point diff_value;
            for (unsigned int d = 0; d < c.dim; ++d) {
              diff_value(d) = args("bc/displacement", NAN, nbc * 4 + 1 + d);
            }
Scale according to current load step
            diff_value *= ratio;
        
            Real penalty_number = args("bc/displacement_penalty", 1e7);
        
            FEBase * fe = c.side_fe_var[var[0]];
            const std::vector<std::vector<Real> > & phi = fe->get_phi();
            const std::vector<Real>& JxW = fe->get_JxW();
            const std::vector<Point>& coords = fe->get_xyz();
        
            unsigned int n_x_dofs = c.dof_indices_var[this->var[0]].size();
        
get mappings for dofs for auxiliary system for original mesh positions
            const System & auxsys = this->get_equation_systems().get_system(
                "auxiliary");
            const DofMap & auxmap = auxsys.get_dof_map();
            std::vector<dof_id_type> undefo_dofs[3];
            for (unsigned int d = 0; d < c.dim; ++d) {
              auxmap.dof_indices(c.elem, undefo_dofs[d], undefo_var[d]);
            }
        
            for (unsigned int qp = 0; qp < c.side_qrule->n_points(); ++qp) {
calculate coordinates of qp on undeformed mesh
              Point orig_point;
              for (unsigned int i = 0; i < n_x_dofs; ++i) {
                for (unsigned int d = 0; d < c.dim; ++d) {
                  Number orig_val = auxsys.current_solution(undefo_dofs[d][i]);
        
        #if LIBMESH_USE_COMPLEX_NUMBERS
                  orig_point(d) += phi[i][qp] * orig_val.real();
        #else
                  orig_point(d) += phi[i][qp] * orig_val;
        #endif
                }
              }
        
Calculate displacement to be enforced.
              Point diff = coords[qp] - orig_point - diff_value;
        
Assemble
              for (unsigned int i = 0; i < n_x_dofs; ++i) {
                for (unsigned int d1 = 0; d1 < c.dim; ++d1) {
                  if (libmesh_isnan(diff(d1)))
                    continue;
                  Real val = JxW[qp] * phi[i][qp] * diff(d1) * penalty_number;
                  c.elem_subresiduals[var[d1]]->operator ()(i) += val;
                }
                if (request_jacobian) {
                  for (unsigned int j = 0; j < n_x_dofs; ++j) {
                    for (unsigned int d1 = 0; d1 < c.dim; ++d1) {
                      if (libmesh_isnan(diff(d1)))
                        continue;
                      Real val = JxW[qp] * phi[i][qp] * phi[j][qp] * penalty_number;
                      c.elem_subjacobians[var[d1]][var[d1]]->operator ()(i, j) += val;
                    }
                  }
                }
              }
            }
          }
        
          return request_jacobian;
        }



The source file nonlinear_neohooke_cc.h without comments:

 
  #ifndef NONLINEAR_NEOHOOKE_CC_H_
  #define NONLINEAR_NEOHOOKE_CC_H_
  
  #include "libmesh/dense_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/vector_value.h"
  #include "libmesh/tensor_value.h"
  #include "libmesh/getpot.h"
  
  using namespace libMesh;
  
  /**
   * This class implements a constitutive formulation for an Neo-Hookean elastic solid
   * in terms of the current configuration. This implementation is suitable for computing
   * large deformations. See e.g. "Nonlinear finite element methods" (P. Wriggers, 2008,
   * Springer) for details.
   */
  class NonlinearNeoHookeCurrentConfig {
  public:
    NonlinearNeoHookeCurrentConfig(
        const std::vector<std::vector<RealGradient> >& dphi_in, GetPot& args) :
        dphi(dphi_in) {
      E = args("material/neohooke/e_modulus", 10000.0);
      nu = args("material/neohooke/nu", 0.3);
    }
  
    /**
     * Initialize the class for the given displacement gradient at the
     * specified quadrature point.
     */
    void init_for_qp(VectorValue<Gradient> & grad_u, unsigned int qp);
  
    /**
     * Return the residual vector for the current state.
     */
    void get_residual(DenseVector<Real> & residuum, unsigned int & i);
  
    /**
     * Return the stiffness matrix for the current state.
     */
    void get_linearized_stiffness(DenseMatrix<Real> & stiffness,
        unsigned int & i, unsigned int & j);
  
    /**
     * Flag to indicate if it is necessary to calculate values for stiffness
     * matrix during initialization.
     */
    bool calculate_linearized_stiffness;
  private:
    void build_b_0_mat(int i, DenseMatrix<Real>& b_l_mat);
    void calculate_stress();
    void calculate_tangent();
    static void tensor_to_voigt(const RealTensor &tensor, DenseVector<Real> &vec);
  
    unsigned int current_qp;
    const std::vector<std::vector<RealGradient> >& dphi;
  
    DenseMatrix<Real> C_mat;
    Real E;
    Real nu;
    RealTensor F, S, tau, sigma;
    DenseMatrix<Real> B_L;
    DenseMatrix<Real> B_K;
  };
  
  #endif /* NONLINEAR_NEOHOOKE_CC_H_ */



The source file solid_system.h without comments:

 
  #ifndef SOLID_SYSTEM_H_
  #define SOLID_SYSTEM_H_
  
  #include "libmesh/fem_system.h"
  #include "libmesh/auto_ptr.h"
  
  using namespace libMesh;
  
  class SolidSystem: public FEMSystem {
  public:
    SolidSystem(EquationSystems& es, const std::string& name,
        const unsigned int number);
  
    virtual void init_data();
  
    virtual void init_context(DiffContext &context);
  
    virtual bool element_time_derivative(bool request_jacobian,
        DiffContext& context);
  
    virtual bool side_time_derivative(bool request_jacobian,
        DiffContext& context);
  
    virtual bool eulerian_residual(bool, DiffContext &) {
      return false;
    }
  
    GetPot args;
  
    virtual std::string system_type() const {
      return "Solid";
    }
  
    virtual void update();
  
    void save_initial_mesh();
  
    unsigned int var[3];
  
    unsigned int undefo_var[3];
  };
  
  #endif /* SOLID_SYSTEM_H_ */



The source file fem_system_ex2.C without comments:

 
  #include "libmesh/equation_systems.h"
  #include "libmesh/getpot.h"
  #include "libmesh/libmesh.h"
  #include "libmesh/libmesh_logging.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/string_to_enum.h"
  #include "libmesh/time_solver.h"
  #include "libmesh/transient_system.h"
  #include "libmesh/vtk_io.h"
  
  #include <cstdio>
  #include <ctime>
  #include <fstream>
  #include <iomanip>
  #include <iostream>
  #include <sstream>
  #include <string>
  #include <unistd.h>
  
  using namespace libMesh;
  
  #include "solid_system.h"
  
  void setup(EquationSystems& systems, Mesh& mesh, GetPot& args)
  {
    const unsigned int dim = mesh.mesh_dimension();
    libmesh_assert (dim == 3);
  
    ElemType eltype = Utility::string_to_enum<ElemType>(args("mesh/generation/element_type", "hex8"));
    int nx = args("mesh/generation/num_elem", 4, 0);
    int ny = args("mesh/generation/num_elem", 4, 1);
    int nz = dim > 2 ? args("mesh/generation/num_elem", 4, 2) : 0;
    double origx = args("mesh/generation/origin", -1.0, 0);
    double origy = args("mesh/generation/origin", -1.0, 1);
    double origz = args("mesh/generation/origin", 0.0, 2);
    double sizex = args("mesh/generation/size", 2.0, 0);
    double sizey = args("mesh/generation/size", 2.0, 1);
    double sizez = args("mesh/generation/size", 2.0, 2);
    MeshTools::Generation::build_cube(mesh, nx, ny, nz,
        origx, origx+sizex, origy, origy+sizey, origz, origz+sizez, eltype);
  
    SolidSystem& imms = systems.add_system<SolidSystem> ("solid");
    imms.args = args;
  
    ExplicitSystem& aux_sys = systems.add_system<TransientExplicitSystem>("auxiliary");
  
    systems.parameters.set<unsigned int>("phase") = 0;
    systems.init();
  
    imms.save_initial_mesh();
  
    aux_sys.reinit();
  }
  
  
  
  void run_timestepping(EquationSystems& systems, GetPot& args)
  {
    TransientExplicitSystem& aux_system = systems.get_system<TransientExplicitSystem>("auxiliary");
  
    SolidSystem& solid_system = systems.get_system<SolidSystem>("solid");
  
    AutoPtr<VTKIO> io = AutoPtr<VTKIO>(new VTKIO(systems.get_mesh()));
  
    Real duration = args("duration", 1.0);
  
    for (unsigned int t_step = 0; t_step < duration/solid_system.deltat; t_step++) {
      Real progress = t_step * solid_system.deltat / duration;
      systems.parameters.set<Real>("progress") = progress;
      systems.parameters.set<unsigned int>("step") = t_step;
  
  
      out << "===== Time Step " << std::setw(4) << t_step;
      out << " (" << std::fixed << std::setprecision(2) << std::setw(6) << (progress*100.) << "%)";
      out << ", time = " << std::setw(7) << solid_system.time;
      out << " =====" << std::endl;
  
      aux_system.current_local_solution->close();
      aux_system.old_local_solution->close();
      *aux_system.older_local_solution = *aux_system.old_local_solution;
      *aux_system.old_local_solution = *aux_system.current_local_solution;
  
      out << "Solving Solid" << std::endl;
      solid_system.solve();
      aux_system.reinit();
  
      out << "Doing a reinit of the equation systems" << std::endl;
      systems.reinit();
  
      if (t_step % args("output/frequency", 1) == 0) {
        std::string result;
        std::stringstream file_name;
        file_name << args("results_directory", "./") << "fem_";
        file_name << std::setw(6) << std::setfill('0') << t_step;
        file_name << ".pvtu";
  
  
        io->write_equation_systems(file_name.str(), systems);
      }
      out << "Advancing to next step" << std::endl;
      solid_system.time_solver->advance_timestep();
    }
  }
  
  
  
  int main(int argc, char** argv)
  {
    LibMeshInit init(argc, argv);
  
  #ifndef LIBMESH_HAVE_VTK
    libmesh_example_assert(false, "--enable-vtk");
  #endif
  
    libmesh_example_assert(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc");
  
    if (libMesh::n_threads() > 1)
      {
        std::cout << "We skip fem_system_ex2 when using threads.\n"
                  << std::endl;
        return 0;
      }
  
    GetPot args = GetPot("solid.in");
  
    int dim = args("mesh/generation/dimension", 3);
    libmesh_example_assert(dim <= LIBMESH_DIM, "3D support");
  
    Mesh mesh(init.comm(), dim);
  
    EquationSystems systems(mesh);
  
    setup(systems, mesh, args);
  
    run_timestepping(systems, args);
  
    out << "Finished calculations" << std::endl;
    return 0;
  }



The source file nonlinear_neohooke_cc.C without comments:

 
  #include "nonlinear_neohooke_cc.h"
  
  /**
   * Return the inverse of the given TypeTensor. Algorithm taken from the tensor classes
   * of Ton van den Boogaard (http://tmku209.ctw.utwente.nl/~ton/tensor.html)
   */
  template <typename T> TypeTensor<T> inv(const TypeTensor<T> &A ) {
    double Sub11, Sub12, Sub13;
    Sub11 = A._coords[4]*A._coords[8] - A._coords[5]*A._coords[7];
    Sub12 = A._coords[3]*A._coords[8] - A._coords[6]*A._coords[5];
    Sub13 = A._coords[3]*A._coords[7] - A._coords[6]*A._coords[4];
    double detA = A._coords[0]*Sub11 - A._coords[1]*Sub12 + A._coords[2]*Sub13;
    libmesh_assert( std::fabs(detA)>1.e-15 );
  
    TypeTensor<T> Ainv(A);
  
    Ainv._coords[0] =  Sub11/detA;
    Ainv._coords[1] = (-A._coords[1]*A._coords[8]+A._coords[2]*A._coords[7])/detA;
    Ainv._coords[2] = ( A._coords[1]*A._coords[5]-A._coords[2]*A._coords[4])/detA;
    Ainv._coords[3] = -Sub12/detA;
    Ainv._coords[4] = ( A._coords[0]*A._coords[8]-A._coords[2]*A._coords[6])/detA;
    Ainv._coords[5] = (-A._coords[0]*A._coords[5]+A._coords[2]*A._coords[3])/detA;
    Ainv._coords[6] =  Sub13/detA;
    Ainv._coords[7] = (-A._coords[0]*A._coords[7]+A._coords[1]*A._coords[6])/detA;
    Ainv._coords[8] = ( A._coords[0]*A._coords[4]-A._coords[1]*A._coords[3])/detA;
  
    return Ainv;
  }
  
  void NonlinearNeoHookeCurrentConfig::init_for_qp(VectorValue<Gradient> & grad_u, unsigned int qp) {
  	this->current_qp = qp;
  	F.zero();
  	S.zero();
  
  	{
  	  RealTensor invF;
  	  invF.zero();
  	  for (unsigned int i = 0; i < 3; ++i)
  	    for (unsigned int j = 0; j < 3; ++j) {
  	      invF(i, j) += libmesh_real(grad_u(i)(j));
  	    }
  	  F.add(inv(invF));
  
  	  libmesh_assert_greater (F.det(), -TOLERANCE);
  	}
  
  	if (this->calculate_linearized_stiffness) {
  		this->calculate_tangent();
  	}
  
  	this->calculate_stress();
  }
  
  
  
  void NonlinearNeoHookeCurrentConfig::calculate_tangent() {
  	Real mu = E / (2 * (1 + nu));
  	Real lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
  
  	Real detF = F.det();
  
  	C_mat.resize(6, 6);
  	for (unsigned int i = 0; i < 3; ++i) {
  		for (unsigned int j = 0; j < 3; ++j) {
  			if (i == j) {
  				C_mat(i, j) = 2 * mu + lambda;
  				C_mat(i + 3, j + 3) = mu - 0.5 * lambda * (detF * detF - 1);
  			} else {
  				C_mat(i, j) = lambda * detF * detF;
  			}
  		}
  	}
  }
  
  
  void NonlinearNeoHookeCurrentConfig::calculate_stress() {
  
    double mu = E / (2.0 * (1.0 + nu));
  	double lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
  
  	Real detF = F.det();
  	RealTensor Ft = F.transpose();
  
  	RealTensor C = Ft * F;
  	RealTensor b = F * Ft;
  	RealTensor identity;
  	identity(0, 0) = 1.0; identity(1, 1) = 1.0; identity(2, 2) = 1.0;
  	RealTensor invC = inv(C);
  
  	S = 0.5 * lambda * (detF * detF - 1) * invC + mu * (identity - invC);
  
  	tau = (F * S) * Ft;
  	sigma = 1/detF * tau;
  }
  
  void NonlinearNeoHookeCurrentConfig::get_residual(DenseVector<Real> & residuum, unsigned int & i) {
  	B_L.resize(3, 6);
  	DenseVector<Real> sigma_voigt(6);
  
  	this->build_b_0_mat(i, B_L);
  
  	tensor_to_voigt(sigma, sigma_voigt);
  
  	B_L.vector_mult(residuum, sigma_voigt);
  }
  
  void NonlinearNeoHookeCurrentConfig::tensor_to_voigt(const RealTensor &tensor, DenseVector<Real> &vec) {
    vec(0) = tensor(0, 0);
    vec(1) = tensor(1, 1);
    vec(2) = tensor(2, 2);
    vec(3) = tensor(0, 1);
    vec(4) = tensor(1, 2);
    vec(5) = tensor(0, 2);
  
  }
  
  void NonlinearNeoHookeCurrentConfig::get_linearized_stiffness(DenseMatrix<Real> & stiffness, unsigned int & i, unsigned int & j) {
  	stiffness.resize(3, 3);
  
  	double G_IK = (sigma * dphi[i][current_qp]) * dphi[j][current_qp];
  	stiffness(0, 0) += G_IK;
  	stiffness(1, 1) += G_IK;
  	stiffness(2, 2) += G_IK;
  
  	B_L.resize(3, 6);
  	this->build_b_0_mat(i, B_L);
  	B_K.resize(3, 6);
  	this->build_b_0_mat(j, B_K);
  
  	B_L.right_multiply(C_mat);
  	B_L.right_multiply_transpose(B_K);
  	B_L *= 1/F.det();
  
  	stiffness += B_L;
  }
  
  void NonlinearNeoHookeCurrentConfig::build_b_0_mat(int i, DenseMatrix<Real>& b_0_mat) {
  	for (unsigned int ii = 0; ii < 3; ++ii) {
  		b_0_mat(ii, ii) = dphi[i][current_qp](ii);
  	}
  	b_0_mat(0, 3) = dphi[i][current_qp](1);
  	b_0_mat(1, 3) = dphi[i][current_qp](0);
  	b_0_mat(1, 4) = dphi[i][current_qp](2);
  	b_0_mat(2, 4) = dphi[i][current_qp](1);
  	b_0_mat(0, 5) = dphi[i][current_qp](2);
  	b_0_mat(2, 5) = dphi[i][current_qp](0);
  }



The source file solid_system.C without comments:

 
  #include "libmesh/boundary_info.h"
  #include "libmesh/diff_solver.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/fe_base.h"
  #include "libmesh/fem_context.h"
  #include "libmesh/getpot.h"
  #include "libmesh/mesh.h"
  #include "libmesh/newton_solver.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/quadrature.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/steady_solver.h"
  #include "libmesh/transient_system.h"
  
  #include "nonlinear_neohooke_cc.h"
  #include "solid_system.h"
  
  #ifdef __SUNPRO_CC
    #define NAN (1.0/0.0)
  #endif
  
  using namespace libMesh;
  
  SolidSystem::SolidSystem(EquationSystems& es, const std::string& name_in,
      const unsigned int number_in) :
      FEMSystem(es, name_in, number_in) {
  
    this->time_solver = AutoPtr<TimeSolver>(new SteadySolver(*this));
  }
  
  void SolidSystem::save_initial_mesh() {
    System & aux_sys = this->get_equation_systems().get_system("auxiliary");
  
    const unsigned int dim = this->get_mesh().mesh_dimension();
  
    const MeshBase::const_node_iterator nd_end =
        this->get_mesh().local_nodes_end();
    for (MeshBase::const_node_iterator nd = this->get_mesh().local_nodes_begin();
        nd != nd_end; ++nd) {
      const Node *node = *nd;
      for (unsigned int d = 0; d < dim; ++d) {
        unsigned int source_dof = node->dof_number(this->number(), var[d], 0);
        unsigned int dest_dof = node->dof_number(aux_sys.number(), undefo_var[d],
            0);
        Number value = this->current_local_solution->el(source_dof);
        aux_sys.current_local_solution->set(dest_dof, value);
      }
    }
  }
  
  void SolidSystem::init_data() {
    const unsigned int dim = this->get_mesh().mesh_dimension();
  
    Order order = (*(this->get_mesh().elements_begin()))->default_order();
  
    var[0] = this->add_variable("x", order);
    var[1] = this->add_variable("y", order);
    if (dim == 3)
      var[2] = this->add_variable("z", order);
    else
      var[2] = var[1];
  
    System& aux_sys = this->get_equation_systems().get_system("auxiliary");
    undefo_var[0] = aux_sys.add_variable("undefo_x", order);
    undefo_var[1] = aux_sys.add_variable("undefo_y", order);
    undefo_var[2] = aux_sys.add_variable("undefo_z", order);
  
    this->deltat = args("schedule/dt", 0.2);
  
    FEMSystem::init_data();
  
  
    set_mesh_system(this);
  
    this->set_mesh_x_var(var[0]);
    this->set_mesh_y_var(var[1]);
    if (dim == 3)
      this->set_mesh_z_var(var[2]);
  
    this->mesh_position_get();
  
    System::reinit();
  
    DiffSolver &solver = *(this->time_solver->diff_solver().get());
    solver.quiet = args("solver/quiet", false);
    solver.max_nonlinear_iterations = args(
        "solver/nonlinear/max_nonlinear_iterations", 100);
    solver.relative_step_tolerance = args(
        "solver/nonlinear/relative_step_tolerance", 1.e-3);
    solver.relative_residual_tolerance = args(
        "solver/nonlinear/relative_residual_tolerance", 1.e-8);
    solver.absolute_residual_tolerance = args(
        "solver/nonlinear/absolute_residual_tolerance", 1.e-8);
    solver.verbose = !args("solver/quiet", false);
  
    ((NewtonSolver&) solver).require_residual_reduction = args(
        "solver/nonlinear/require_reduction", false);
  
    solver.max_linear_iterations = args("max_linear_iterations", 50000);
    solver.initial_linear_tolerance = args("initial_linear_tolerance", 1.e-3);
  }
  
  void SolidSystem::update() {
    System::update();
    this->mesh_position_set();
  }
  
  void SolidSystem::init_context(DiffContext &context) {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    c.element_fe_var[var[0]]->get_JxW();
    c.element_fe_var[var[0]]->get_phi();
    c.element_fe_var[var[0]]->get_dphi();
    c.element_fe_var[var[0]]->get_xyz();
  
    c.side_fe_var[var[0]]->get_JxW();
    c.side_fe_var[var[0]]->get_phi();
    c.side_fe_var[var[0]]->get_xyz();
  }
  
  /**
   * Compute contribution from internal forces in elem_residual and contribution from
   * linearized internal forces (stiffness matrix) in elem_jacobian.
   */
  bool SolidSystem::element_time_derivative(bool request_jacobian,
      DiffContext &context) {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW = c.element_fe_var[var[0]]->get_JxW();
  
    const std::vector<std::vector<RealGradient> >& dphi =
        c.element_fe_var[var[0]]->get_dphi();
  
    const unsigned int dim = this->get_mesh().mesh_dimension();
  
    const unsigned int n_u_dofs = c.dof_indices_var[var[0]].size();
    libmesh_assert(n_u_dofs == c.dof_indices_var[var[1]].size());
    if (dim == 3) {
      libmesh_assert(n_u_dofs == c.dof_indices_var[var[2]].size());
    }
  
    unsigned int n_qpoints = c.element_qrule->n_points();
  
    DenseMatrix<Real> stiff;
    DenseVector<Real> res;
    VectorValue<Gradient> grad_u;
  
    NonlinearNeoHookeCurrentConfig material(dphi, args);
  
    material.calculate_linearized_stiffness = request_jacobian;
  
    TransientExplicitSystem& aux_system = this->get_equation_systems().get_system<
        TransientExplicitSystem>("auxiliary");
    std::vector<dof_id_type> undefo_index;
  
    bool use_symmetry = args("assembly/use_symmetry", false);
  
  
    for (unsigned int qp = 0; qp != n_qpoints; qp++) {
      grad_u(0) = grad_u(1) = grad_u(2) = 0;
      for (unsigned int d = 0; d < dim; ++d) {
        std::vector<Number> u_undefo;
        aux_system.get_dof_map().dof_indices(c.elem, undefo_index, undefo_var[d]);
        aux_system.current_local_solution->get(undefo_index, u_undefo);
        for (unsigned int l = 0; l != n_u_dofs; l++)
          grad_u(d).add_scaled(dphi[l][qp], u_undefo[l]); // u_current(l)); // -
      }
  
      material.init_for_qp(grad_u, qp);
  
      for (unsigned int i = 0; i < n_u_dofs; i++) {
        res.resize(dim);
        material.get_residual(res, i);
        res.scale(JxW[qp]);
        for (unsigned int ii = 0; ii < dim; ++ii) {
          c.elem_subresiduals[ii]->operator ()(i) += res(ii);
        }
  
        if (request_jacobian && c.elem_solution_derivative) {
          libmesh_assert(c.elem_solution_derivative == 1.0);
          for (unsigned int j = (use_symmetry ? i : 0); j < n_u_dofs; j++) {
            material.get_linearized_stiffness(stiff, i, j);
            stiff.scale(JxW[qp]);
            for (unsigned int ii = 0; ii < dim; ++ii) {
              for (unsigned int jj = 0; jj < dim; ++jj) {
                c.elem_subjacobians[ii][jj]->operator ()(i, j) += stiff(ii, jj);
                if (use_symmetry && i != j) {
                  c.elem_subjacobians[ii][jj]->operator ()(j, i) += stiff(jj, ii);
                }
              }
            }
          }
        }
      }
    } // end of the quadrature point qp-loop
  
    return request_jacobian;
  }
  
  bool SolidSystem::side_time_derivative(bool request_jacobian,
      DiffContext &context) {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    Real ratio = this->get_equation_systems().parameters.get<Real>("progress")
        + 0.001;
  
  
    unsigned int num_bc = args.vector_variable_size("bc/displacement");
    if (num_bc % 4 != 0) {
      libMesh::err
          << "ERROR, Odd number of values in displacement boundary condition.\n"
          << std::endl;
      libmesh_error();
    }
    num_bc /= 4;
  
    for (unsigned int nbc = 0; nbc < num_bc; nbc++) {
      short int positive_boundary_id = args("bc/displacement", 1, nbc * 4);
  
      if (!this->get_mesh().boundary_info->has_boundary_id
  	  (c.elem,c.side,positive_boundary_id))
        continue;
  
      Point diff_value;
      for (unsigned int d = 0; d < c.dim; ++d) {
        diff_value(d) = args("bc/displacement", NAN, nbc * 4 + 1 + d);
      }
      diff_value *= ratio;
  
      Real penalty_number = args("bc/displacement_penalty", 1e7);
  
      FEBase * fe = c.side_fe_var[var[0]];
      const std::vector<std::vector<Real> > & phi = fe->get_phi();
      const std::vector<Real>& JxW = fe->get_JxW();
      const std::vector<Point>& coords = fe->get_xyz();
  
      unsigned int n_x_dofs = c.dof_indices_var[this->var[0]].size();
  
      const System & auxsys = this->get_equation_systems().get_system(
          "auxiliary");
      const DofMap & auxmap = auxsys.get_dof_map();
      std::vector<dof_id_type> undefo_dofs[3];
      for (unsigned int d = 0; d < c.dim; ++d) {
        auxmap.dof_indices(c.elem, undefo_dofs[d], undefo_var[d]);
      }
  
      for (unsigned int qp = 0; qp < c.side_qrule->n_points(); ++qp) {
        Point orig_point;
        for (unsigned int i = 0; i < n_x_dofs; ++i) {
          for (unsigned int d = 0; d < c.dim; ++d) {
            Number orig_val = auxsys.current_solution(undefo_dofs[d][i]);
  
  #if LIBMESH_USE_COMPLEX_NUMBERS
            orig_point(d) += phi[i][qp] * orig_val.real();
  #else
            orig_point(d) += phi[i][qp] * orig_val;
  #endif
          }
        }
  
        Point diff = coords[qp] - orig_point - diff_value;
  
        for (unsigned int i = 0; i < n_x_dofs; ++i) {
          for (unsigned int d1 = 0; d1 < c.dim; ++d1) {
            if (libmesh_isnan(diff(d1)))
              continue;
            Real val = JxW[qp] * phi[i][qp] * diff(d1) * penalty_number;
            c.elem_subresiduals[var[d1]]->operator ()(i) += val;
          }
          if (request_jacobian) {
            for (unsigned int j = 0; j < n_x_dofs; ++j) {
              for (unsigned int d1 = 0; d1 < c.dim; ++d1) {
                if (libmesh_isnan(diff(d1)))
                  continue;
                Real val = JxW[qp] * phi[i][qp] * phi[j][qp] * penalty_number;
                c.elem_subjacobians[var[d1]][var[d1]]->operator ()(i, j) += val;
              }
            }
          }
        }
      }
    }
  
    return request_jacobian;
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/fem_system/fem_system_ex2'
***************************************************************
* Running Example fem_system_ex2:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
*** Warning, This code is untested, experimental, or likely to see future API changes: ../../../include/libmesh/diff_physics.h, line 382, compiled Apr 19 2013 at 11:50:09 ***
*** Warning, This code is untested, experimental, or likely to see future API changes: ../src/mesh/vtk_io.C, line 358, compiled Apr 19 2013 at 11:42:41 ***
===== Time Step    0 (  0.00%), time =    0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 3691.41
Linear solve starting, tolerance 0.00
Linear solve finished, step 12, residual 0.00
Trying full Newton step
  Current Residual: 0.16
  Nonlinear solver current_residual 0.16 > 0.00
  Nonlinear solver relative residual 0.00 > 0.00
  Nonlinear solver converged, step 0, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step    1 ( 20.00%), time =    0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 738516.15
Linear solve starting, tolerance 0.00
Linear solve finished, step 12, residual 0.00
Trying full Newton step
  Current Residual: 52.19
  Nonlinear step: |du|/|u| = 0.06, |du| = 1.04
Assembling the System
Nonlinear Residual: 52.19
Linear solve starting, tolerance 0.00
Linear solve finished, step 15, residual 0.00
Trying full Newton step
  Current Residual: 1.56
  Nonlinear solver current_residual 1.56 > 0.00
  Nonlinear solver relative residual 0.00 > 0.00
  Nonlinear solver converged, step 1, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step    2 ( 40.00%), time =    0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 788391.92
Linear solve starting, tolerance 0.00
Linear solve finished, step 12, residual 0.00
Trying full Newton step
  Current Residual: 58.83
  Nonlinear step: |du|/|u| = 0.06, |du| = 1.05
Assembling the System
Nonlinear Residual: 58.83
Linear solve starting, tolerance 0.00
Linear solve finished, step 16, residual 0.00
Trying full Newton step
  Current Residual: 3.93
  Nonlinear step: |du|/|u| = 0.00, |du| = 0.02
Assembling the System
Nonlinear Residual: 3.93
Linear solve starting, tolerance 0.00
Linear solve finished, step 15, residual 0.00
Trying full Newton step
  Current Residual: 0.01
  Nonlinear solver current_residual 0.01 > 0.00
  Nonlinear solver relative residual 0.00 > 0.00
  Nonlinear solver converged, step 2, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step    3 ( 60.00%), time =    0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 844767.72
Linear solve starting, tolerance 0.00
Linear solve finished, step 12, residual 0.00
Trying full Newton step
  Current Residual: 84.75
  Nonlinear step: |du|/|u| = 0.06, |du| = 1.06
Assembling the System
Nonlinear Residual: 84.75
Linear solve starting, tolerance 0.00
Linear solve finished, step 16, residual 0.00
Trying full Newton step
  Current Residual: 8.60
  Nonlinear step: |du|/|u| = 0.00, |du| = 0.02
Assembling the System
Nonlinear Residual: 8.60
Linear solve starting, tolerance 0.00
Linear solve finished, step 15, residual 0.00
Trying full Newton step
  Current Residual: 0.05
  Nonlinear solver current_residual 0.05 > 0.00
  Nonlinear solver relative residual 0.00 > 0.00
  Nonlinear solver converged, step 2, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step    4 ( 80.00%), time =    0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 909688.38
Linear solve starting, tolerance 0.00
Linear solve finished, step 13, residual 0.00
Trying full Newton step
  Current Residual: 165.58
  Nonlinear step: |du|/|u| = 0.07, |du| = 1.07
Assembling the System
Nonlinear Residual: 165.58
Linear solve starting, tolerance 0.00
Linear solve finished, step 17, residual 0.00
Trying full Newton step
  Current Residual: 20.40
  Nonlinear step: |du|/|u| = 0.00, |du| = 0.03
Assembling the System
Nonlinear Residual: 20.40
Linear solve starting, tolerance 0.00
Linear solve finished, step 20, residual 0.00
Trying full Newton step
  Current Residual: 0.14
  Nonlinear solver current_residual 0.14 > 0.00
  Nonlinear solver relative residual 0.00 > 0.00
  Nonlinear solver converged, step 2, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
Finished calculations

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:50:22 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-10.3.12.361-mpich2-1.4.1p1-cxx-opt'                                                |
|  'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt'                                     |
|  'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.862794, Active time=0.741512                                                 |
 ----------------------------------------------------------------------------------------------------------------
| Event                              nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                              w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|----------------------------------------------------------------------------------------------------------------|
|                                                                                                                |
|                                                                                                                |
| DofMap                                                                                                         |
|   add_neighbors_to_send_list()     12        0.0027      0.000227    0.0095      0.000795    0.37     1.29     |
|   build_sparsity()                 6         0.0020      0.000332    0.0068      0.001128    0.27     0.91     |
|   create_dof_constraints()         12        0.0005      0.000041    0.0005      0.000041    0.07     0.07     |
|   distribute_dofs()                12        0.0023      0.000195    0.0117      0.000974    0.32     1.58     |
|   dof_indices()                    16316     0.0858      0.000005    0.0858      0.000005    11.58    11.58    |
|   old_dof_indices()                1280      0.0079      0.000006    0.0079      0.000006    1.07     1.07     |
|   prepare_send_list()              12        0.0003      0.000027    0.0003      0.000027    0.04     0.04     |
|   reinit()                         12        0.0028      0.000233    0.0028      0.000233    0.38     0.38     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          5         0.0010      0.000191    0.0065      0.001297    0.13     0.87     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        3200      0.0318      0.000010    0.0318      0.000010    4.29     4.29     |
|   init_shape_functions()           1100      0.0098      0.000009    0.0098      0.000009    1.32     1.32     |
|                                                                                                                |
| FEMSystem                                                                                                      |
|   assembly()                       12        0.1370      0.011416    0.2158      0.017980    18.47    29.10    |
|   assembly(get_residual)           12        0.0288      0.002402    0.1080      0.009000    3.89     14.56    |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             138       0.0006      0.000005    0.0006      0.000005    0.09     0.09     |
|   compute_face_map()               960       0.0034      0.000004    0.0034      0.000004    0.46     0.46     |
|   compute_map()                    3062      0.0255      0.000008    0.0255      0.000008    3.44     3.44     |
|   init_face_shape_functions()      48        0.0004      0.000007    0.0004      0.000007    0.05     0.05     |
|   init_reference_to_physical_map() 1100      0.0121      0.000011    0.0121      0.000011    1.63     1.63     |
|                                                                                                                |
| LocationMap                                                                                                    |
|   init()                           5         0.0002      0.000046    0.0002      0.000046    0.03     0.03     |
|                                                                                                                |
| Mesh                                                                                                           |
|   contract()                       5         0.0000      0.000010    0.0001      0.000023    0.01     0.02     |
|   find_neighbors()                 1         0.0003      0.000280    0.0005      0.000496    0.04     0.07     |
|   renumber_nodes_and_elem()        7         0.0001      0.000014    0.0001      0.000014    0.01     0.01     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        2         0.0005      0.000231    0.0005      0.000231    0.06     0.06     |
|   find_global_indices()            2         0.0001      0.000075    0.0015      0.000737    0.02     0.20     |
|   parallel_sort()                  2         0.0003      0.000163    0.0006      0.000280    0.04     0.08     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         5         0.0073      0.001460    0.0141      0.002826    0.98     1.91     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _coarsen_elements()              5         0.0000      0.000006    0.0001      0.000015    0.00     0.01     |
|   _refine_elements()               5         0.0001      0.000017    0.0003      0.000051    0.01     0.03     |
|   make_coarsening_compatible()     10        0.0001      0.000008    0.0001      0.000008    0.01     0.01     |
|   make_flags_parallel_consistent() 10        0.0005      0.000052    0.0034      0.000338    0.07     0.46     |
|   make_refinement_compatible()     10        0.0000      0.000001    0.0001      0.000005    0.00     0.01     |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0002      0.000160    0.0002      0.000160    0.02     0.02     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0010      0.000959    0.0016      0.001582    0.13     0.21     |
|                                                                                                                |
| NewtonSolver                                                                                                   |
|   solve()                          5         0.2210      0.044191    0.6325      0.126505    29.80    85.30    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      40        0.0026      0.000064    0.0026      0.000065    0.35     0.35     |
|   max(bool)                        21        0.0003      0.000015    0.0003      0.000015    0.04     0.04     |
|   max(scalar)                      1325      0.0048      0.000004    0.0048      0.000004    0.65     0.65     |
|   max(vector)                      323       0.0018      0.000006    0.0054      0.000017    0.24     0.72     |
|   min(bool)                        1663      0.0059      0.000004    0.0059      0.000004    0.80     0.80     |
|   min(scalar)                      1312      0.0611      0.000047    0.0611      0.000047    8.24     8.24     |
|   min(vector)                      323       0.0020      0.000006    0.0057      0.000018    0.28     0.76     |
|   probe()                          594       0.0040      0.000007    0.0040      0.000007    0.54     0.54     |
|   receive()                        594       0.0013      0.000002    0.0054      0.000009    0.18     0.73     |
|   send()                           594       0.0007      0.000001    0.0007      0.000001    0.10     0.10     |
|   send_receive()                   598       0.0019      0.000003    0.0086      0.000014    0.26     1.16     |
|   sum()                            68        0.0006      0.000009    0.0018      0.000027    0.08     0.24     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           594       0.0004      0.000001    0.0004      0.000001    0.06     0.06     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0001      0.000088    0.0006      0.000621    0.01     0.08     |
|   set_parent_processor_ids()       1         0.0000      0.000022    0.0000      0.000022    0.00     0.00     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          12        0.0475      0.003957    0.0475      0.003957    6.40     6.40     |
|                                                                                                                |
| ProjectVector                                                                                                  |
|   operator()                       20        0.0023      0.000113    0.0108      0.000538    0.30     1.45     |
|                                                                                                                |
| System                                                                                                         |
|   project_vector()                 20        0.0176      0.000881    0.0347      0.001736    2.38     4.68     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            35478     0.7415                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example fem_system_ex2:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/fem_system/fem_system_ex2'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo