LargeDeformationElasticity Class Reference

Inheritance diagram for LargeDeformationElasticity:

Public Member Functions
	LargeDeformationElasticity (EquationSystems &es_in)
Real	kronecker_delta (unsigned int i, unsigned int j)
	Kronecker delta function. More...
Real	elasticity_tensor (Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
	Evaluate the fourth order tensor (C_ijkl) that relates stress to strain. More...
virtual void	jacobian (const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, NonlinearImplicitSystem &)
	Evaluate the Jacobian of the nonlinear system. More...
virtual void	residual (const NumericVector< Number > &soln, NumericVector< Number > &residual, NonlinearImplicitSystem &)
	Evaluate the residual of the nonlinear system. More...
void	compute_stresses ()
	Compute the Cauchy stress for the current solution. More...
virtual void	residual (const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)=0
	Residual function. More...
virtual void	jacobian (const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)=0
	Jacobian function. More...

Private Attributes
EquationSystems &	es

Detailed Description

Definition at line 87 of file systems_of_equations_ex7.C.

Constructor & Destructor Documentation

LargeDeformationElasticity()

LargeDeformationElasticity::LargeDeformationElasticity ( EquationSystems &  es_in )

inline

Definition at line 95 of file systems_of_equations_ex7.C.

                                                       :
    es(es_in)
  {}

Member Function Documentation

compute_stresses()

void LargeDeformationElasticity::compute_stresses ( )

inline

Compute the Cauchy stress for the current solution.

Definition at line 381 of file systems_of_equations_ex7.C.

References libMesh::TypeTensor< T >::add_scaled(), libMesh::FEGenericBase< OutputType >::build(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_order(), libMesh::TypeTensor< T >::det(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::n_vars(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::System::solution, libMesh::TypeTensor< T >::transpose(), libMesh::System::update(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), and libMesh::TypeTensor< T >::zero().

Referenced by main().

  {
    const Real young_modulus = es.parameters.get<Real>("young_modulus");
    const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");

    const MeshBase & mesh = es.get_mesh();
    const unsigned int dim = mesh.mesh_dimension();

    NonlinearImplicitSystem & system =
      es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");

    unsigned int displacement_vars[] = {
      system.variable_number ("u"),
      system.variable_number ("v"),
      system.variable_number ("w")};
    const unsigned int u_var = system.variable_number ("u");

    const DofMap & dof_map = system.get_dof_map();
    FEType fe_type = dof_map.variable_type(u_var);
    std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);

    const std::vector<Real> & JxW = fe->get_JxW();
    const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

    // Also, get a reference to the ExplicitSystem
    ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
    const DofMap & stress_dof_map = stress_system.get_dof_map();
    unsigned int sigma_vars[] = {
      stress_system.variable_number ("sigma_00"),
      stress_system.variable_number ("sigma_01"),
      stress_system.variable_number ("sigma_02"),
      stress_system.variable_number ("sigma_11"),
      stress_system.variable_number ("sigma_12"),
      stress_system.variable_number ("sigma_22")};

    // Storage for the stress dof indices on each element
    std::vector<std::vector<dof_id_type>> dof_indices_var(system.n_vars());
    std::vector<dof_id_type> stress_dof_indices_var;

    // To store the stress tensor on each element
    TensorValue<Number> elem_avg_stress_tensor;

    for (const auto & elem : mesh.active_local_element_ptr_range())
      {
        for (unsigned int var=0; var<3; var++)
          dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);

        const unsigned int n_var_dofs = dof_indices_var[0].size();

        fe->reinit (elem);

        // clear the stress tensor
        elem_avg_stress_tensor.zero();

        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            TensorValue<Number> grad_u;
            // Row is variable u1, u2, or u3, column is x, y, or z
            for (unsigned int var_i=0; var_i<3; var_i++)
              for (unsigned int var_j=0; var_j<3; var_j++)
                for (unsigned int j=0; j<n_var_dofs; j++)
                  grad_u(var_i,var_j) += dphi[j][qp](var_j) * system.current_solution(dof_indices_var[var_i][j]);

            TensorValue<Number> strain_tensor;
            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                {
                  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));

                  for (unsigned int k=0; k<3; k++)
                    strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
                }

            // Define the deformation gradient
            auto F = grad_u;
            for (unsigned int var=0; var<3; var++)
              F(var, var) += 1.;

            TensorValue<Number> stress_tensor;
            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                for (unsigned int k=0; k<3; k++)
                  for (unsigned int l=0; l<3; l++)
                    stress_tensor(i,j) +=
                      elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k, l);

            // stress_tensor now holds the second Piola-Kirchoff stress (PK2) at point qp.
            // However, in this example we want to compute the Cauchy stress which is given by
            // 1/det(F) * F * PK2 * F^T, hence we now apply this transformation.
            stress_tensor = 1. / F.det() * F * stress_tensor * F.transpose();

            // We want to plot the average Cauchy stress on each element, hence
            // we integrate stress_tensor
            elem_avg_stress_tensor.add_scaled(stress_tensor, JxW[qp]);
          }

        // Get the average stress per element by dividing by volume
        elem_avg_stress_tensor /= elem->volume();

        // load elem_sigma data into stress_system
        unsigned int stress_var_index = 0;
        for (unsigned int i=0; i<3; i++)
          for (unsigned int j=i; j<3; j++)
            {
              stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);

              // We are using CONSTANT MONOMIAL basis functions, hence we only need to get
              // one dof index per variable
              dof_id_type dof_index = stress_dof_indices_var[0];

              if ((stress_system.solution->first_local_index() <= dof_index) &&
                  (dof_index < stress_system.solution->last_local_index()))
                stress_system.solution->set(dof_index, elem_avg_stress_tensor(i,j));

              stress_var_index++;
            }
      }

    // Should call close and update when we set vector entries directly
    stress_system.solution->close();
    stress_system.update();
  }

elasticity_tensor()

Real LargeDeformationElasticity::elasticity_tensor ( Real  young_modulus, Real  poisson_ratio, unsigned int  i, unsigned int  j, unsigned int  k, unsigned int  l  )

inline

Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.

Definition at line 111 of file systems_of_equations_ex7.C.

References kronecker_delta(), and libMesh::Real.

  {
    // Define the Lame constants
    const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
    const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));

    return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l) +
      lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
  }

jacobian() [1/2]

virtual void libMesh::NonlinearImplicitSystem::ComputeJacobian::jacobian ( const NumericVector< Number > &  X, SparseMatrix< Number > &  J, sys_type &  S  )

pure virtualinherited

Jacobian function.

This function will be called to compute the jacobian and must be implemented by the user in a derived class.

Referenced by libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), and libMesh::libmesh_petsc_snes_jacobian().

jacobian() [2/2]

virtual void LargeDeformationElasticity::jacobian ( const NumericVector< Number > &  soln, SparseMatrix< Number > &  jacobian, NonlinearImplicitSystem &    )

inlinevirtual

Evaluate the Jacobian of the nonlinear system.

Definition at line 130 of file systems_of_equations_ex7.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_matrix(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseMatrix< T >::resize(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), and libMesh::SparseMatrix< T >::zero().

  {
    const Real young_modulus = es.parameters.get<Real>("young_modulus");
    const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");

    const MeshBase & mesh = es.get_mesh();
    const unsigned int dim = mesh.mesh_dimension();

    NonlinearImplicitSystem & system =
      es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");

    const unsigned int u_var = system.variable_number ("u");

    const DofMap & dof_map = system.get_dof_map();

    FEType fe_type = dof_map.variable_type(u_var);
    std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);

    std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
    QGauss qface (dim-1, fe_type.default_quadrature_order());
    fe_face->attach_quadrature_rule (&qface);

    const std::vector<Real> & JxW = fe->get_JxW();
    const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

    DenseMatrix<Number> Ke;
    DenseSubMatrix<Number> Ke_var[3][3] =
      {
        {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)},
        {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)},
        {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke)}
      };

    std::vector<dof_id_type> dof_indices;
    std::vector<std::vector<dof_id_type>> dof_indices_var(3);

    jacobian.zero();

    for (const auto & elem : mesh.active_local_element_ptr_range())
      {
        dof_map.dof_indices (elem, dof_indices);
        for (unsigned int var=0; var<3; var++)
          dof_map.dof_indices (elem, dof_indices_var[var], var);

        const unsigned int n_dofs = dof_indices.size();
        const unsigned int n_var_dofs = dof_indices_var[0].size();

        fe->reinit (elem);

        Ke.resize (n_dofs, n_dofs);
        for (unsigned int var_i=0; var_i<3; var_i++)
          for (unsigned int var_j=0; var_j<3; var_j++)
            Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);

        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            TensorValue<Number> grad_u;
            for (unsigned int var_i=0; var_i<3; var_i++)
              {
                // Row is variable u1, u2, or u3, column is x, y, or z
                for (unsigned int var_j=0; var_j<3; var_j++)
                  for (unsigned int j=0; j<n_var_dofs; j++)
                    grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
              }

            TensorValue<Number> strain_tensor;
            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                {
                  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));

                  for (unsigned int k=0; k<3; k++)
                    strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
                }

            // Define the deformation gradient
            auto F = grad_u;
            for (unsigned int var=0; var<3; var++)
              F(var, var) += 1.;

            TensorValue<Number> stress_tensor;

            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                for (unsigned int k=0; k<3; k++)
                  for (unsigned int l=0; l<3; l++)
                    stress_tensor(i,j) +=
                      elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k, l);

            for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
              for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
                {
                  for (unsigned int i=0; i<3; i++)
                    for (unsigned int j=0; j<3; j++)
                      for (unsigned int m=0; m<3; m++)
                        Ke_var[i][i](dof_i,dof_j) += JxW[qp] *
                          (-dphi[dof_j][qp](m) * stress_tensor(m,j) * dphi[dof_i][qp](j));

                  for (unsigned int i=0; i<3; i++)
                    for (unsigned int j=0; j<3; j++)
                      for (unsigned int k=0; k<3; k++)
                        for (unsigned int l=0; l<3; l++)
                          {
                            Number FxC_ijkl = 0.;
                            for (unsigned int m=0; m<3; m++)
                              FxC_ijkl += F(i,m) * elasticity_tensor(young_modulus, poisson_ratio, m, j, k, l);

                            Ke_var[i][k](dof_i,dof_j) += JxW[qp] *
                              (-0.5 * FxC_ijkl * dphi[dof_j][qp](l) * dphi[dof_i][qp](j));

                            Ke_var[i][l](dof_i,dof_j) += JxW[qp] *
                              (-0.5 * FxC_ijkl * dphi[dof_j][qp](k) * dphi[dof_i][qp](j));

                            for (unsigned int n=0; n<3; n++)
                              Ke_var[i][n](dof_i,dof_j) += JxW[qp] *
                                (-0.5 * FxC_ijkl * (dphi[dof_j][qp](k) * grad_u(n,l) + dphi[dof_j][qp](l) * grad_u(n,k)) * dphi[dof_i][qp](j));
                          }
                }
          }

        dof_map.constrain_element_matrix (Ke, dof_indices);
        jacobian.add_matrix (Ke, dof_indices);
      }
  }

kronecker_delta()

Real LargeDeformationElasticity::kronecker_delta ( unsigned int  i, unsigned int  j  )

inline

Kronecker delta function.

Definition at line 102 of file systems_of_equations_ex7.C.

  {
    return i == j ? 1. : 0.;
  }

residual() [1/2]

virtual void libMesh::NonlinearImplicitSystem::ComputeResidual::residual ( const NumericVector< Number > &  X, NumericVector< Number > &  R, sys_type &  S  )

pure virtualinherited

Residual function.

This function will be called to compute the residual and must be implemented by the user in a derived class.

Referenced by libMesh::Problem_Interface::computeF(), libMesh::libmesh_petsc_snes_fd_residual(), libMesh::libmesh_petsc_snes_mffd_residual(), and libMesh::libmesh_petsc_snes_residual().

residual() [2/2]

virtual void LargeDeformationElasticity::residual ( const NumericVector< Number > &  soln, NumericVector< Number > &  residual, NonlinearImplicitSystem &    )

inlinevirtual

Evaluate the residual of the nonlinear system.

Definition at line 262 of file systems_of_equations_ex7.C.

References libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_vector(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), and libMesh::NumericVector< T >::zero().

  {
    const Real young_modulus = es.parameters.get<Real>("young_modulus");
    const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
    const Real forcing_magnitude = es.parameters.get<Real>("forcing_magnitude");

    const MeshBase & mesh = es.get_mesh();
    const unsigned int dim = mesh.mesh_dimension();

    NonlinearImplicitSystem & system =
      es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");

    const unsigned int u_var = system.variable_number ("u");

    const DofMap & dof_map = system.get_dof_map();

    FEType fe_type = dof_map.variable_type(u_var);
    std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);

    std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
    QGauss qface (dim-1, fe_type.default_quadrature_order());
    fe_face->attach_quadrature_rule (&qface);

    const std::vector<Real> & JxW = fe->get_JxW();
    const std::vector<std::vector<Real>> & phi = fe->get_phi();
    const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

    DenseVector<Number> Re;

    DenseSubVector<Number> Re_var[3] =
      {DenseSubVector<Number>(Re),
       DenseSubVector<Number>(Re),
       DenseSubVector<Number>(Re)};

    std::vector<dof_id_type> dof_indices;
    std::vector<std::vector<dof_id_type>> dof_indices_var(3);

    residual.zero();

    for (const auto & elem : mesh.active_local_element_ptr_range())
      {
        dof_map.dof_indices (elem, dof_indices);
        for (unsigned int var=0; var<3; var++)
          dof_map.dof_indices (elem, dof_indices_var[var], var);

        const unsigned int n_dofs = dof_indices.size();
        const unsigned int n_var_dofs = dof_indices_var[0].size();

        fe->reinit (elem);

        Re.resize (n_dofs);
        for (unsigned int var=0; var<3; var++)
          Re_var[var].reposition (var*n_var_dofs, n_var_dofs);

        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            TensorValue<Number> grad_u;
            for (unsigned int var_i=0; var_i<3; var_i++)
              {
                // Row is variable u, v, or w column is x, y, or z
                for (unsigned int var_j=0; var_j<3; var_j++)
                  for (unsigned int j=0; j<n_var_dofs; j++)
                    grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
              }

            TensorValue<Number> strain_tensor;
            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                {
                  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));

                  for (unsigned int k=0; k<3; k++)
                    strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
                }

            // Define the deformation gradient
            auto F = grad_u;
            for (unsigned int var=0; var<3; var++)
              F(var, var) += 1.;

            TensorValue<Number> stress_tensor;

            for (unsigned int i=0; i<3; i++)
              for (unsigned int j=0; j<3; j++)
                for (unsigned int k=0; k<3; k++)
                  for (unsigned int l=0; l<3; l++)
                    stress_tensor(i,j) +=
                      elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k,l);

            VectorValue<Number> f_vec(0., 0., -forcing_magnitude);

            for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
              for (unsigned int i=0; i<3; i++)
                {
                  for (unsigned int j=0; j<3; j++)
                    {
                      Number FxStress_ij = 0.;
                      for (unsigned int m=0; m<3; m++)
                        FxStress_ij += F(i,m) * stress_tensor(m,j);

                      Re_var[i](dof_i) += JxW[qp] * (-FxStress_ij * dphi[dof_i][qp](j));
                    }

                  Re_var[i](dof_i) += JxW[qp] * (f_vec(i) * phi[dof_i][qp]);
                }
          }

        dof_map.constrain_element_vector (Re, dof_indices);
        residual.add_vector (Re, dof_indices);
      }
  }

Member Data Documentation

es

EquationSystems& LargeDeformationElasticity::es

private

Definition at line 91 of file systems_of_equations_ex7.C.

---

The documentation for this class was generated from the following file:

• examples/systems_of_equations/systems_of_equations_ex7/systems_of_equations_ex7.C