miscellaneous_ex11.C File Reference

Go to the source code of this file.

Functions
void	assemble_shell (EquationSystems &es, const std::string &system_name)
int	main (int argc, char **argv)
void	assemble_shell (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))

Function Documentation

assemble_shell() [1/2]

void assemble_shell	(	EquationSystems &	es,
		const std::string &	system_name
	)

Definition at line 328 of file miscellaneous_ex12.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::NumericVector< T >::close(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::Parameters::get(), libMesh::MeshBase::get_boundary_info(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), libMesh::libmesh_ignore(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseSubVector< T >::reposition(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, libMesh::SECOND, libMesh::NumericVector< T >::set(), libMesh::BoundaryInfo::shellface_boundary_ids(), and libMesh::System::variable_type().

Referenced by main().

{
  // This example requires Eigen to actually work, but we should still
  // let it compile and throw a runtime error if you don't.

  // The same holds for second derivatives,
  // since they are class-members only depending on the config.
#if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "Shell");

  // Get a constant reference to the mesh object.
  const MeshBase & mesh = es.get_mesh();
  const unsigned int dim = mesh.mesh_dimension();

  // Get a reference to the shell system object.
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> (system_name);

  // Get the shell parameters that we need during assembly.
  const Real h  = es.parameters.get<Real> ("thickness");
  const Real E  = es.parameters.get<Real> ("young's modulus");
  const Real nu = es.parameters.get<Real> ("poisson ratio");
  const Real q  = es.parameters.get<Real> ("point load");
  const bool distributed_load  = es.parameters.get<bool> ("distributed load");

  // The membrane elastic matrix.
  MyMatrix3d Hm;
  Hm <<
    1., nu, 0.,
    nu, 1., 0.,
    0., 0., 0.5 * (1-nu);
  Hm *= h * E/(1-nu*nu);

  // The bending elastic matrix.
  MyMatrix3d Hf;
  Hf <<
    1., nu, 0.,
    nu, 1., 0.,
    0., 0., 0.5 * (1-nu);
  Hf *= h*h*h/12 * E/(1-nu*nu);

  // The shear elastic matrices.
  MyMatrix2d Hc0 = MyMatrix2d::Identity();
  Hc0 *= h * 5./6*E/(2*(1+nu));

  MyMatrix2d Hc1 = MyMatrix2d::Identity();
  Hc1 *= h*h*h/12 * 5./6*E/(2*(1+nu));

  // Get the Finite Element type, this will be
  // the same for all variables.
  FEType fe_type = system.variable_type (0);

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

  // The element Jacobian * quadrature weight at each integration point.
  const std::vector<Real> & JxW = fe->get_JxW();

  // The element shape function and its derivatives evaluated at the
  // quadrature points.
  const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
  const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();

  const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
  const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
  const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
  const std::vector<std::vector<Real>> & dphidxi = fe->get_dphidxi();
  const std::vector<std::vector<Real>> & dphideta = fe->get_dphideta();
  const std::vector<std::vector<Real>> & phi = fe->get_phi();

  // A reference to the DofMap object for this system.  The DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.
  const DofMap & dof_map = system.get_dof_map();

  // The global system matrix
  SparseMatrix<Number> & matrix = system.get_system_matrix();

  // Define data structures to contain the element stiffness matrix.
  DenseMatrix<Number> Ke;
  DenseSubMatrix<Number> Ke_var[6][6] =
    {
      {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),
       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)},
      {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),
       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)}
    };

  // Define data structures to contain the element rhs vector.
  DenseVector<Number> Fe;
  DenseSubVector<Number> Fe_w(Fe);

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

  // Now we will loop over all the elements in the mesh.  We will
  // compute the element matrix and right-hand-side contribution.
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      dof_map.dof_indices (elem, dof_indices);
      for (unsigned int var=0; var<6; 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();

      // First compute element data at the nodes
      std::vector<Point> nodes;
      for (auto i : elem->node_index_range())
        nodes.push_back(elem->reference_elem()->node_ref(i));
      fe->reinit (elem, &nodes);

      // Convenient notation for the element node positions
      MyVector3d X1(elem->node_ref(0)(0), elem->node_ref(0)(1), elem->node_ref(0)(2));
      MyVector3d X2(elem->node_ref(1)(0), elem->node_ref(1)(1), elem->node_ref(1)(2));
      MyVector3d X3(elem->node_ref(2)(0), elem->node_ref(2)(1), elem->node_ref(2)(2));
      MyVector3d X4(elem->node_ref(3)(0), elem->node_ref(3)(1), elem->node_ref(3)(2));

      //Store covariant basis and local orthonormal basis at the nodes
      std::vector<MyMatrix3d> F0node;
      std::vector<MyMatrix3d> Qnode;
      for (auto i : elem->node_index_range())
        {
          MyVector3d a1;
          a1 << dxyzdxi[i](0), dxyzdxi[i](1), dxyzdxi[i](2);
          MyVector3d a2;
          a2 << dxyzdeta[i](0), dxyzdeta[i](1), dxyzdeta[i](2);
          MyVector3d n;
          n = a1.cross(a2);
          n /= n.norm();
          MyMatrix3d F0;
          F0 <<
            a1(0), a2(0), n(0),
            a1(1), a2(1), n(1),
            a1(2), a2(2), n(2);
          F0node.push_back(F0);

          Real nx = n(0);
          Real ny = n(1);
          Real C  = n(2);
          if (std::abs(1.+C)<1e-6)
            {
              MyMatrix3d Q;
              Q <<
                1, 0, 0,
                0, -1, 0,
                0, 0, -1;
              Qnode.push_back(Q);
            }
          else
            {
              MyMatrix3d Q;
              Q <<
                C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
                -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
                -nx,             -ny,              C;
              Qnode.push_back(Q);
            }
        }

      Ke.resize (n_dofs, n_dofs);
      for (unsigned int var_i=0; var_i<6; var_i++)
        for (unsigned int var_j=0; var_j<6; 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);

      Fe.resize(n_dofs);
      Fe_w.reposition(2*n_var_dofs,n_var_dofs);

      // Reinit element data at the regular Gauss quadrature points
      fe->reinit (elem);

      // Now we will build the element matrix and right-hand-side.
      for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
        {

          //Covariant basis at the quadrature point
          MyVector3d a1;
          a1 << dxyzdxi[qp](0), dxyzdxi[qp](1), dxyzdxi[qp](2);
          MyVector3d a2;
          a2 << dxyzdeta[qp](0), dxyzdeta[qp](1), dxyzdeta[qp](2);
          MyVector3d n;
          n = a1.cross(a2);
          n /= n.norm();
          MyMatrix3d F0;
          F0 <<
            a1(0), a2(0), n(0),
            a1(1), a2(1), n(1),
            a1(2), a2(2), n(2);

          //Contravariant basis
          MyMatrix3d F0it;
          F0it = F0.inverse().transpose();

          //Local orthonormal basis at the quadrature point
          Real nx = n(0);
          Real ny = n(1);
          Real C  = n(2);
          MyMatrix3d Q;
          if (std::abs(1.+C) < 1e-6)
            {
              Q <<
                1, 0, 0,
                0, -1, 0,
                0, 0, -1;
            }
          else
            {
              Q <<
                C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
                -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
                -nx,             -ny,              C;
            }

          MyMatrix2d C0;
          C0 = F0it.block<3,2>(0,0).transpose()*Q.block<3,2>(0,0);

          // Normal derivatives in reference coordinates
          MyVector3d d2Xdxi2(d2xyzdxi2[qp](0), d2xyzdxi2[qp](1), d2xyzdxi2[qp](2));
          MyVector3d d2Xdeta2(d2xyzdeta2[qp](0), d2xyzdeta2[qp](1), d2xyzdeta2[qp](2));
          MyVector3d d2Xdxideta(d2xyzdxideta[qp](0), d2xyzdxideta[qp](1), d2xyzdxideta[qp](2));


          MyMatrix2d b;
          b <<
            n.dot(d2Xdxi2), n.dot(d2Xdxideta),
            n.dot(d2Xdxideta), n.dot(d2Xdeta2);

          MyVector3d dndxi = -b(0,0)*F0it.col(0) - b(0,1)*F0it.col(1);
          MyVector3d dndeta = -b(1,0)*F0it.col(0) - b(1,1)*F0it.col(1);

          MyMatrix2d bhat;
          bhat <<
            F0it.col(1).dot(dndeta), -F0it.col(0).dot(dndeta),
            -F0it.col(1).dot(dndxi), F0it.col(0).dot(dndxi);

          MyMatrix2d bc;
          bc = bhat*C0;

          // Mean curvature
          Real H = 0.5*(dndxi.dot(F0it.col(0))+dndeta.dot(F0it.col(1)));

          // Quadrature point reference coordinates
          Real xi = qrule.qp(qp)(0);
          Real eta = qrule.qp(qp)(1);

          // Preassemble the MITC4 shear strain matrix for all nodes as they involve
          // cross references to midside nodes.
          // The QUAD4 element has nodes X1,X2,X3,X4 with coordinates (xi,eta)
          // in the reference element: (-1,-1),(1,-1),(1,1),(-1,1).
          // The midside nodes are denoted A1=(X1+X2)/2, B2=(X2+X3)/2, A2=(X3+X4)/2, B1=(X4+X1)/2.

          // Normals at the midside nodes (average of normals at the edge corners).
          // Multiplication by the assumed shear strain shape function.
          MyVector3d nA1 = 0.5*(Qnode[0].col(2)+Qnode[1].col(2));
          nA1 /= nA1.norm();
          nA1 *= (1-eta)/4;
          MyVector3d nB2 = 0.5*(Qnode[1].col(2)+Qnode[2].col(2));
          nB2 /= nB2.norm();
          nB2 *= (1+xi)/4;
          MyVector3d nA2 = 0.5*(Qnode[2].col(2)+Qnode[3].col(2));
          nA2 /= nA2.norm();
          nA2 *= (1+eta)/4;
          MyVector3d nB1 = 0.5*(Qnode[3].col(2)+Qnode[0].col(2));
          nB1 /= nB1.norm();
          nB1 *= (1-xi)/4;

          // Edge tangents
          MyVector3d aA1 = 0.5*(X2-X1);
          MyVector3d aA2 = 0.5*(X3-X4);
          MyVector3d aB1 = 0.5*(X4-X1);
          MyVector3d aB2 = 0.5*(X3-X2);

          // Contribution of the rotational dofs to the shear strain
          MyVector2d AS1A1(-aA1.dot(Qnode[0].col(1)), aA1.dot(Qnode[0].col(0)));
          MyVector2d AS2A1(-aA1.dot(Qnode[1].col(1)), aA1.dot(Qnode[1].col(0)));
          AS1A1 *= (1-eta)/4;
          AS2A1 *= (1-eta)/4;

          MyVector2d AS1A2(-aA2.dot(Qnode[3].col(1)), aA2.dot(Qnode[3].col(0)));
          MyVector2d AS2A2(-aA2.dot(Qnode[2].col(1)), aA2.dot(Qnode[2].col(0)));
          AS1A2 *= (1+eta)/4;
          AS2A2 *= (1+eta)/4;

          MyVector2d AS1B1(-aB1.dot(Qnode[0].col(1)), aB1.dot(Qnode[0].col(0)));
          MyVector2d AS2B1(-aB1.dot(Qnode[3].col(1)), aB1.dot(Qnode[3].col(0)));
          AS1B1 *= (1-xi)/4;
          AS2B1 *= (1-xi)/4;

          MyVector2d AS1B2(-aB2.dot(Qnode[1].col(1)), aB2.dot(Qnode[1].col(0)));
          MyVector2d AS2B2(-aB2.dot(Qnode[2].col(1)), aB2.dot(Qnode[2].col(0)));
          AS1B2 *= (1+xi)/4;
          AS2B2 *= (1+xi)/4;

          // Store previous quantities in the shear strain matrices for each node
          std::vector<MyMatrixXd> Bcnode;
          MyMatrixXd Bc(2, 5);
          // Node 1
          Bc.block<1,3>(0,0) = -nA1.transpose();
          Bc.block<1,2>(0,3) = AS1A1.transpose();
          Bc.block<1,3>(1,0) = -nB1.transpose();
          Bc.block<1,2>(1,3) = AS1B1.transpose();
          Bcnode.push_back(Bc);
          // Node 2
          Bc.block<1,3>(0,0) = nA1.transpose();
          Bc.block<1,2>(0,3) = AS2A1.transpose();
          Bc.block<1,3>(1,0) = -nB2.transpose();
          Bc.block<1,2>(1,3) = AS1B2.transpose();
          Bcnode.push_back(Bc);
          // Node 3
          Bc.block<1,3>(0,0) = nA2.transpose();
          Bc.block<1,2>(0,3) = AS2A2.transpose();
          Bc.block<1,3>(1,0) = nB2.transpose();
          Bc.block<1,2>(1,3) = AS2B2.transpose();
          Bcnode.push_back(Bc);
          // Node 4
          Bc.block<1,3>(0,0) = -nA2.transpose();
          Bc.block<1,2>(0,3) = AS1A2.transpose();
          Bc.block<1,3>(1,0) = nB1.transpose();
          Bc.block<1,2>(1,3) = AS2B1.transpose();
          Bcnode.push_back(Bc);

          // Loop over all pairs of nodes I,J.
          for (unsigned int i=0; i<n_var_dofs; ++i)
            {
              // Matrix B0, zeroth order (through thickness) membrane-bending strain
              Real C1i = dphidxi[i][qp]*C0(0,0) + dphideta[i][qp]*C0(1,0);
              Real C2i = dphidxi[i][qp]*C0(0,1) + dphideta[i][qp]*C0(1,1);

              MyMatrixXd B0I(3, 5);
              B0I = MyMatrixXd::Zero(3, 5);
              B0I.block<1,3>(0,0) = C1i*Q.col(0).transpose();
              B0I.block<1,3>(1,0) = C2i*Q.col(1).transpose();
              B0I.block<1,3>(2,0) = C2i*Q.col(0).transpose()+C1i*Q.col(1).transpose();

              // Matrix B1, first order membrane-bending strain
              Real bc1i = dphidxi[i][qp]*bc(0,0) + dphideta[i][qp]*bc(1,0);
              Real bc2i = dphidxi[i][qp]*bc(0,1) + dphideta[i][qp]*bc(1,1);

              MyVector2d V1i(-Q.col(0).dot(Qnode[i].col(1)),
                             Q.col(0).dot(Qnode[i].col(0)));

              MyVector2d V2i(-Q.col(1).dot(Qnode[i].col(1)),
                             Q.col(1).dot(Qnode[i].col(0)));

              MyMatrixXd B1I(3,5);
              B1I = MyMatrixXd::Zero(3,5);
              B1I.block<1,3>(0,0) = bc1i*Q.col(0).transpose();
              B1I.block<1,3>(1,0) = bc2i*Q.col(1).transpose();
              B1I.block<1,3>(2,0) = bc2i*Q.col(0).transpose()+bc1i*Q.col(1).transpose();

              B1I.block<1,2>(0,3) = C1i*V1i.transpose();
              B1I.block<1,2>(1,3) = C2i*V2i.transpose();
              B1I.block<1,2>(2,3) = C2i*V1i.transpose()+C1i*V2i.transpose();

              // Matrix B2, second order membrane-bending strain
              MyMatrixXd B2I(3,5);
              B2I = MyMatrixXd::Zero(3,5);

              B2I.block<1,2>(0,3) = bc1i*V1i.transpose();
              B2I.block<1,2>(1,3) = bc2i*V2i.transpose();
              B2I.block<1,2>(2,3) = bc2i*V1i.transpose()+bc1i*V2i.transpose();

              // Matrix Bc0, zeroth order shear strain
              MyMatrixXd Bc0I(2,5);
              Bc0I = C0.transpose()*Bcnode[i];

              // Matrix Bc1, first order shear strain
              MyMatrixXd Bc1I(2,5);
              Bc1I = bc.transpose()*Bcnode[i];

              // Drilling dof (in-plane rotation)
              MyVector2d BdxiI(dphidxi[i][qp],dphideta[i][qp]);
              MyVector2d BdI = C0.transpose()*BdxiI;

              for (unsigned int j=0; j<n_var_dofs; ++j)
                {

                  // Matrix B0, zeroth order membrane-bending strain
                  Real C1j = dphidxi[j][qp]*C0(0,0) + dphideta[j][qp]*C0(1,0);
                  Real C2j = dphidxi[j][qp]*C0(0,1) + dphideta[j][qp]*C0(1,1);

                  MyMatrixXd B0J(3,5);
                  B0J = MyMatrixXd::Zero(3,5);
                  B0J.block<1,3>(0,0) = C1j*Q.col(0).transpose();
                  B0J.block<1,3>(1,0) = C2j*Q.col(1).transpose();
                  B0J.block<1,3>(2,0) = C2j*Q.col(0).transpose()+C1j*Q.col(1).transpose();

                  // Matrix B1, first order membrane-bending strain
                  Real bc1j = dphidxi[j][qp]*bc(0,0) + dphideta[j][qp]*bc(1,0);
                  Real bc2j = dphidxi[j][qp]*bc(0,1) + dphideta[j][qp]*bc(1,1);

                  MyVector2d V1j(-Q.col(0).dot(Qnode[j].col(1)),
                                 Q.col(0).dot(Qnode[j].col(0)));

                  MyVector2d V2j(-Q.col(1).dot(Qnode[j].col(1)),
                                 Q.col(1).dot(Qnode[j].col(0)));

                  MyMatrixXd B1J(3,5);
                  B1J = MyMatrixXd::Zero(3,5);
                  B1J.block<1,3>(0,0) = bc1j*Q.col(0).transpose();
                  B1J.block<1,3>(1,0) = bc2j*Q.col(1).transpose();
                  B1J.block<1,3>(2,0) = bc2j*Q.col(0).transpose()+bc1j*Q.col(1).transpose();

                  B1J.block<1,2>(0,3) = C1j*V1j.transpose();
                  B1J.block<1,2>(1,3) = C2j*V2j.transpose();
                  B1J.block<1,2>(2,3) = C2j*V1j.transpose()+C1j*V2j.transpose();

                  // Matrix B2, second order membrane-bending strain
                  MyMatrixXd B2J(3,5);
                  B2J = MyMatrixXd::Zero(3,5);

                  B2J.block<1,2>(0,3) = bc1j*V1j.transpose();
                  B2J.block<1,2>(1,3) = bc2j*V2j.transpose();
                  B2J.block<1,2>(2,3) = bc2j*V1j.transpose()+bc1j*V2j.transpose();

                  // Matrix Bc0, zeroth order shear strain
                  MyMatrixXd Bc0J(2, 5);
                  Bc0J = C0.transpose()*Bcnode[j];

                  // Matrix Bc1, first order shear strain
                  MyMatrixXd Bc1J(2, 5);
                  Bc1J = bc.transpose()*Bcnode[j];

                  // Drilling dof
                  MyVector2d BdxiJ(dphidxi[j][qp], dphideta[j][qp]);
                  MyVector2d BdJ = C0.transpose()*BdxiJ;

                  // The total stiffness matrix coupling the nodes
                  // I and J is a sum of membrane, bending and shear contributions.
                  MyMatrixXd local_KIJ(5, 5);
                  local_KIJ = JxW[qp] * (
                                         B0I.transpose() * Hm * B0J
                                         +  B2I.transpose() * Hf * B0J
                                         +  B0I.transpose() * Hf * B2J
                                         +  B1I.transpose() * Hf * B1J
                                         +  2*H * B0I.transpose() * Hf * B1J
                                         +  2*H * B1I.transpose() * Hf * B0J
                                         +  Bc0I.transpose() * Hc0 * Bc0J
                                         +  Bc1I.transpose() * Hc1 * Bc1J
                                         +  2*H * Bc0I.transpose() * Hc1 * Bc1J
                                         +  2*H * Bc1I.transpose() * Hc1 * Bc0J
                                         );

                  // Going from 5 to 6 dofs to add drilling dof
                  MyMatrixXd full_local_KIJ(6, 6);
                  full_local_KIJ = MyMatrixXd::Zero(6, 6);
                  full_local_KIJ.block<5,5>(0,0)=local_KIJ;

                  // Drilling dof stiffness contribution
                  // Note that in the original book, there is a coefficient of
                  // alpha between 1e-4 and 1e-7 to make the fictitious
                  // drilling stiffness small while preventing the stiffness
                  // matrix from being singular. For this problem, we can use
                  // alpha = 1 to also get a good result.
                  //
                  // The explicit conversion to Real here is to work
                  // around an Eigen+boost::float128 incompatibility.
                  full_local_KIJ(5,5) = Real(Hf(0,0)*JxW[qp]*BdI.transpose()*BdJ);

                  // Transform the stiffness matrix to global coordinates
                  MyMatrixXd global_KIJ(6,6);
                  MyMatrixXd TI(6,6);
                  TI = MyMatrixXd::Identity(6,6);
                  TI.block<3,3>(3,3) = Qnode[i].transpose();
                  MyMatrixXd TJ(6,6);
                  TJ = MyMatrixXd::Identity(6,6);
                  TJ.block<3,3>(3,3) = Qnode[j].transpose();
                  global_KIJ = TI.transpose()*full_local_KIJ*TJ;

                  // Insert the components of the coupling stiffness
                  // matrix KIJ into the corresponding directional
                  // submatrices.
                  for (unsigned int k=0;k<6;k++)
                    for (unsigned int l=0;l<6;l++)
                      Ke_var[k][l](i,j) += global_KIJ(k,l);
                }
            }

        } // end of the quadrature point qp-loop

      if (distributed_load)
        {
          // Loop on shell faces
          for (unsigned int shellface=0; shellface<2; shellface++)
            {
              std::vector<boundary_id_type> bids;
              mesh.get_boundary_info().shellface_boundary_ids(elem, shellface, bids);

              for (std::size_t k=0; k<bids.size(); k++)
                if (bids[k]==11) // sideset id for surface load
                  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
                    for (unsigned int i=0; i<n_var_dofs; ++i)
                      Fe_w(i) -= JxW[qp] * phi[i][qp];
            }
        }

      // The element matrix is now built for this element.
      // Add it to the global matrix.

      dof_map.constrain_element_matrix_and_vector (Ke,Fe,dof_indices);

      matrix.add_matrix         (Ke, dof_indices);
      system.rhs->add_vector    (Fe, dof_indices);

    }

  if (!distributed_load)
    {
      //Adding point load to the RHS

      //Pinch position
      Point C(0, 3, 3);

      //Finish assembling rhs so we can set one value
      system.rhs->close();

      for (const auto & node : mesh.node_ptr_range())
        if (((*node) - C).norm() < 1e-3)
          system.rhs->set(node->dof_number(0, 2, 0), -q/4);
    }

#else
  // Avoid compiler warnings
  libmesh_ignore(es, system_name);
#endif // defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
}

assemble_shell() [2/2]

void assemble_shell	(	EquationSystems &	es,
		const std::string &	libmesh_dbg_varsystem_name
	)

Definition at line 264 of file miscellaneous_ex11.C.

References libMesh::SparseMatrix< T >::add(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::TypeVector< T >::cross(), libMesh::FEType::default_quadrature_rule(), libMesh::DofMap::dof_indices(), libMesh::DofObject::dof_number(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), libMesh::Tri3Subdivision::is_ghost(), mesh, libMesh::MeshTools::Subdivision::next, libMesh::Elem::node_ptr(), libMesh::TypeVector< T >::norm(), libMesh::System::number(), libMesh::EquationSystems::parameters, libMesh::MeshTools::Subdivision::prev, libMesh::Real, libMesh::DenseSubVector< T >::reposition(), libMesh::DenseSubMatrix< T >::reposition(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, libMesh::TypeTensor< T >::transpose(), libMesh::TRI3SUBDIVISION, libMesh::System::variable_number(), and libMesh::System::variable_type().

{
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "Shell");

  // Get a constant reference to the mesh object.
  const MeshBase & mesh = es.get_mesh();

  // Get a reference to the shell system object.
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> ("Shell");

  // Get the shell parameters that we need during assembly.
  const Real h  = es.parameters.get<Real> ("thickness");
  const Real E  = es.parameters.get<Real> ("young's modulus");
  const Real nu = es.parameters.get<Real> ("poisson ratio");
  const Real q  = es.parameters.get<Real> ("uniform load");

  // Compute the membrane stiffness K and the bending
  // rigidity D from these parameters.
  const Real K = E * h     /     (1-nu*nu);
  const Real D = E * h*h*h / (12*(1-nu*nu));

  // Numeric ids corresponding to each variable in the system.
  const unsigned int u_var = system.variable_number ("u");
  const unsigned int v_var = system.variable_number ("v");
  const unsigned int w_var = system.variable_number ("w");

  // Get the Finite Element type for "u".  Note this will be
  // the same as the type for "v" and "w".
  FEType fe_type = system.variable_type (u_var);

  // Build a Finite Element object of the specified type.
  std::unique_ptr<FEBase> fe (FEBase::build(2, fe_type));

  // A Gauss quadrature rule for numerical integration.
  // For subdivision shell elements, a single Gauss point per
  // element is sufficient, hence we use extraorder = 0.
  const int extraorder = 0;
  std::unique_ptr<QBase> qrule (fe_type.default_quadrature_rule (2, extraorder));

  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule (qrule.get());

  // The element Jacobian * quadrature weight at each integration point.
  const std::vector<Real> & JxW = fe->get_JxW();

  // The surface tangents in both directions at the quadrature points.
  const std::vector<RealGradient> & dxyzdxi  = fe->get_dxyzdxi();
  const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();

  // The second partial derivatives at the quadrature points.
  const std::vector<RealGradient> & d2xyzdxi2    = fe->get_d2xyzdxi2();
  const std::vector<RealGradient> & d2xyzdeta2   = fe->get_d2xyzdeta2();
  const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();

  // The element shape function and its derivatives evaluated at the
  // quadrature points.
  const std::vector<std::vector<Real>> &          phi = fe->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
  const std::vector<std::vector<RealTensor>> &  d2phi = fe->get_d2phi();

  // A reference to the DofMap object for this system.  The DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.
  const DofMap & dof_map = system.get_dof_map();

  // The global system matrix
  SparseMatrix<Number> & matrix = system.get_system_matrix();

  // Define data structures to contain the element stiffness matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke" and "Fe".
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  DenseSubMatrix<Number>
    Kuu(Ke), Kuv(Ke), Kuw(Ke),
    Kvu(Ke), Kvv(Ke), Kvw(Ke),
    Kwu(Ke), Kwv(Ke), Kww(Ke);

  DenseSubVector<Number>
    Fu(Fe),
    Fv(Fe),
    Fw(Fe);

  // This vector will hold the degree of freedom indices for
  // the element.  These define where in the global system
  // the element degrees of freedom get mapped.
  std::vector<dof_id_type> dof_indices;
  std::vector<dof_id_type> dof_indices_u;
  std::vector<dof_id_type> dof_indices_v;
  std::vector<dof_id_type> dof_indices_w;

  // Now we will loop over all the elements in the mesh.  We will
  // compute the element matrix and right-hand-side contribution.
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // The ghost elements at the boundaries need to be excluded
      // here, as they don't belong to the physical shell,
      // but serve for a proper boundary treatment only.
      libmesh_assert_equal_to (elem->type(), TRI3SUBDIVISION);
      const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *> (elem);
      if (sd_elem->is_ghost())
        continue;

      // Get the degree of freedom indices for the
      // current element.  These define where in the global
      // matrix and right-hand-side this element will
      // contribute to.
      dof_map.dof_indices (elem, dof_indices);
      dof_map.dof_indices (elem, dof_indices_u, u_var);
      dof_map.dof_indices (elem, dof_indices_v, v_var);
      dof_map.dof_indices (elem, dof_indices_w, w_var);

      const std::size_t n_dofs   = dof_indices.size();
      const std::size_t n_u_dofs = dof_indices_u.size();
      const std::size_t n_v_dofs = dof_indices_v.size();
      const std::size_t n_w_dofs = dof_indices_w.size();

      // Compute the element-specific data for the current
      // element.  This involves computing the location of the
      // quadrature points and the shape functions
      // (phi, dphi, d2phi) for the current element.
      fe->reinit (elem);

      // Zero the element matrix and right-hand side before
      // summing them.  We use the resize member here because
      // the number of degrees of freedom might have changed from
      // the last element.
      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);

      // Reposition the submatrices...  The idea is this:
      //
      //         -           -          -  -
      //        | Kuu Kuv Kuw |        | Fu |
      //   Ke = | Kvu Kvv Kvw |;  Fe = | Fv |
      //        | Kwu Kwv Kww |        | Fw |
      //         -           -          -  -
      //
      // The DenseSubMatrix.reposition () member takes the
      // (row_offset, column_offset, row_size, column_size).
      //
      // Similarly, the DenseSubVector.reposition () member
      // takes the (row_offset, row_size)
      Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
      Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
      Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);

      Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
      Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
      Kvw.reposition (v_var*n_v_dofs, w_var*n_v_dofs, n_v_dofs, n_w_dofs);

      Kwu.reposition (w_var*n_w_dofs, u_var*n_w_dofs, n_w_dofs, n_u_dofs);
      Kwv.reposition (w_var*n_w_dofs, v_var*n_w_dofs, n_w_dofs, n_v_dofs);
      Kww.reposition (w_var*n_w_dofs, w_var*n_w_dofs, n_w_dofs, n_w_dofs);

      Fu.reposition (u_var*n_u_dofs, n_u_dofs);
      Fv.reposition (v_var*n_u_dofs, n_v_dofs);
      Fw.reposition (w_var*n_u_dofs, n_w_dofs);

      // Now we will build the element matrix and right-hand-side.
      for (unsigned int qp=0; qp<qrule->n_points(); ++qp)
        {
          // First, we compute the external force resulting
          // from a load q distributed uniformly across the plate.
          // Since the load is supposed to be transverse to the plate,
          // it affects the z-direction, i.e. the "w" variable.
          for (unsigned int i=0; i<n_u_dofs; ++i)
            Fw(i) += JxW[qp] * phi[i][qp] * q;

          // Next, we assemble the stiffness matrix.  This is only valid
          // for the linear theory, i.e., for small deformations, where
          // reference and deformed surface metrics are indistinguishable.

          // Get the three surface basis vectors.
          const RealVectorValue & a1 = dxyzdxi[qp];
          const RealVectorValue & a2 = dxyzdeta[qp];
          RealVectorValue   a3 = a1.cross(a2);
          const Real jac = a3.norm(); // the surface Jacobian
          libmesh_assert_greater (jac, 0);
          a3 /= jac; // the shell director a3 is normalized to unit length

          // Get the derivatives of the surface tangents.
          const RealVectorValue & a11 = d2xyzdxi2[qp];
          const RealVectorValue & a22 = d2xyzdeta2[qp];
          const RealVectorValue & a12 = d2xyzdxideta[qp];

          // Compute the three covariant components of the first
          // fundamental form of the surface.
          const RealVectorValue a(a1*a1, a2*a2, a1*a2);

          // The elastic H matrix in Voigt's notation, computed from the
          // covariant components of the first fundamental form rather
          // than the contravariant components, exploiting that the
          // contravariant first fundamental form is the inverse of the
          // covariant first fundamental form (hence the determinant etc.).
          RealTensorValue H;
          H(0,0) = a(1) * a(1);
          H(0,1) = H(1,0) = nu * a(1) * a(0) + (1-nu) * a(2) * a(2);
          H(0,2) = H(2,0) = -a(1) * a(2);
          H(1,1) = a(0) * a(0);
          H(1,2) = H(2,1) = -a(0) * a(2);
          H(2,2) = 0.5 * ((1-nu) * a(1) * a(0) + (1+nu) * a(2) * a(2));
          const Real det = a(0) * a(1) - a(2) * a(2);
          libmesh_assert_not_equal_to (det * det, 0);
          H /= det * det;

          // Precompute come cross products for the bending part below.
          const RealVectorValue a11xa2 = a11.cross(a2);
          const RealVectorValue a22xa2 = a22.cross(a2);
          const RealVectorValue a12xa2 = a12.cross(a2);
          const RealVectorValue a1xa11 =  a1.cross(a11);
          const RealVectorValue a1xa22 =  a1.cross(a22);
          const RealVectorValue a1xa12 =  a1.cross(a12);
          const RealVectorValue a2xa3  =  a2.cross(a3);
          const RealVectorValue a3xa1  =  a3.cross(a1);

          // Loop over all pairs of nodes I,J.
          for (unsigned int i=0; i<n_u_dofs; ++i)
            {
              for (unsigned int j=0; j<n_u_dofs; ++j)
                {
                  // The membrane strain matrices in Voigt's notation.
                  RealTensorValue MI, MJ;
                  for (unsigned int k=0; k<3; ++k)
                    {
                      MI(0,k) = dphi[i][qp](0) * a1(k);
                      MI(1,k) = dphi[i][qp](1) * a2(k);
                      MI(2,k) = dphi[i][qp](1) * a1(k)
                        + dphi[i][qp](0) * a2(k);

                      MJ(0,k) = dphi[j][qp](0) * a1(k);
                      MJ(1,k) = dphi[j][qp](1) * a2(k);
                      MJ(2,k) = dphi[j][qp](1) * a1(k)
                        + dphi[j][qp](0) * a2(k);
                    }

                  // The bending strain matrices in Voigt's notation.
                  RealTensorValue BI, BJ;
                  for (unsigned int k=0; k<3; ++k)
                    {
                      const Real term_ik = dphi[i][qp](0) * a2xa3(k)
                        + dphi[i][qp](1) * a3xa1(k);
                      BI(0,k) = -d2phi[i][qp](0,0) * a3(k)
                        +(dphi[i][qp](0) * a11xa2(k)
                          + dphi[i][qp](1) * a1xa11(k)
                          + (a3*a11) * term_ik) / jac;
                      BI(1,k) = -d2phi[i][qp](1,1) * a3(k)
                        +(dphi[i][qp](0) * a22xa2(k)
                          + dphi[i][qp](1) * a1xa22(k)
                          + (a3*a22) * term_ik) / jac;
                      BI(2,k) = 2 * (-d2phi[i][qp](0,1) * a3(k)
                                     +(dphi[i][qp](0) * a12xa2(k)
                                       + dphi[i][qp](1) * a1xa12(k)
                                       + (a3*a12) * term_ik) / jac);

                      const Real term_jk = dphi[j][qp](0) * a2xa3(k)
                        + dphi[j][qp](1) * a3xa1(k);
                      BJ(0,k) = -d2phi[j][qp](0,0) * a3(k)
                        +(dphi[j][qp](0) * a11xa2(k)
                          + dphi[j][qp](1) * a1xa11(k)
                          + (a3*a11) * term_jk) / jac;
                      BJ(1,k) = -d2phi[j][qp](1,1) * a3(k)
                        +(dphi[j][qp](0) * a22xa2(k)
                          + dphi[j][qp](1) * a1xa22(k)
                          + (a3*a22) * term_jk) / jac;
                      BJ(2,k) = 2 * (-d2phi[j][qp](0,1) * a3(k)
                                     +(dphi[j][qp](0) * a12xa2(k)
                                       + dphi[j][qp](1) * a1xa12(k)
                                       + (a3*a12) * term_jk) / jac);
                    }

                  // The total stiffness matrix coupling the nodes
                  // I and J is a sum of membrane and bending
                  // contributions according to the following formula.
                  const RealTensorValue KIJ = JxW[qp] * K * MI.transpose() * H * MJ
                    + JxW[qp] * D * BI.transpose() * H * BJ;

                  // Insert the components of the coupling stiffness
                  // matrix KIJ into the corresponding directional
                  // submatrices.
                  Kuu(i,j) += KIJ(0,0);
                  Kuv(i,j) += KIJ(0,1);
                  Kuw(i,j) += KIJ(0,2);

                  Kvu(i,j) += KIJ(1,0);
                  Kvv(i,j) += KIJ(1,1);
                  Kvw(i,j) += KIJ(1,2);

                  Kwu(i,j) += KIJ(2,0);
                  Kwv(i,j) += KIJ(2,1);
                  Kww(i,j) += KIJ(2,2);
                }
            }

        } // end of the quadrature point qp-loop

      // The element matrix and right-hand-side are now built
      // for this element.  Add them to the global matrix and
      // right-hand-side vector.  The NumericMatrix::add_matrix()
      // and NumericVector::add_vector() members do this for us.
      matrix.add_matrix         (Ke, dof_indices);
      system.rhs->add_vector    (Fe, dof_indices);
    } // end of non-ghost element loop

  // Next, we apply the boundary conditions.  In this case,
  // all boundaries are clamped by the penalty method, using
  // the special "ghost" nodes along the boundaries.  Note
  // that there are better ways to implement boundary conditions
  // for subdivision shells.  We use the simplest way here,
  // which is known to be overly restrictive and will lead to
  // a slightly too small deformation of the plate.
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // For the boundary conditions, we only need to loop over
      // the ghost elements.
      libmesh_assert_equal_to (elem->type(), TRI3SUBDIVISION);
      const Tri3Subdivision * gh_elem = static_cast<const Tri3Subdivision *> (elem);
      if (!gh_elem->is_ghost())
        continue;

      // Find the side which is part of the physical plate boundary,
      // that is, the boundary of the original mesh without ghosts.
      for (auto s : elem->side_index_range())
        {
          const Tri3Subdivision * nb_elem = static_cast<const Tri3Subdivision *> (elem->neighbor_ptr(s));
          if (nb_elem == nullptr || nb_elem->is_ghost())
            continue;

          /*
           * Determine the four nodes involved in the boundary
           * condition treatment of this side.  The MeshTools::Subdiv
           * namespace provides lookup tables next and prev
           * for an efficient determination of the next and previous
           * nodes of an element, respectively.
           *
           *      n4
           *     /  \
           *    / gh \
           *  n2 ---- n3
           *    \ nb /
           *     \  /
           *      n1
           */
          const Node * nodes [4]; // n1, n2, n3, n4
          nodes[1] = gh_elem->node_ptr(s); // n2
          nodes[2] = gh_elem->node_ptr(MeshTools::Subdivision::next[s]); // n3
          nodes[3] = gh_elem->node_ptr(MeshTools::Subdivision::prev[s]); // n4

          // The node in the interior of the domain, n1, is the
          // hardest to find.  Walk along the edges of element nb until
          // we have identified it.
          unsigned int n_int = 0;
          nodes[0] = nb_elem->node_ptr(0);
          while (nodes[0]->id() == nodes[1]->id() || nodes[0]->id() == nodes[2]->id())
            nodes[0] = nb_elem->node_ptr(++n_int);

          // The penalty value.  \f$ \frac{1}{\epsilon} \f$
          const Real penalty = 1.e10;

          // With this simple method, clamped boundary conditions are
          // obtained by penalizing the displacements of all four nodes.
          // This ensures that the displacement field vanishes on the
          // boundary side s.
          for (unsigned int n=0; n<4; ++n)
            {
              const dof_id_type u_dof = nodes[n]->dof_number (system.number(), u_var, 0);
              const dof_id_type v_dof = nodes[n]->dof_number (system.number(), v_var, 0);
              const dof_id_type w_dof = nodes[n]->dof_number (system.number(), w_var, 0);
              matrix.add (u_dof, u_dof, penalty);
              matrix.add (v_dof, v_dof, penalty);
              matrix.add (w_dof, w_dof, penalty);
            }
        }
    } // end of ghost element loop
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 80 of file miscellaneous_ex11.C.

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_shell(), libMesh::System::attach_assemble_function(), libMesh::ParallelObject::comm(), libMesh::command_line_next(), libMesh::System::current_solution(), libMesh::default_solver_package(), libMesh::DofObject::dof_number(), libMesh::DofMapBase::end_dof(), libMesh::DofMapBase::first_dof(), libMesh::MeshTools::Modification::flatten(), libMesh::FOURTH, libMesh::System::get_dof_map(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), libMesh::TypeVector< T >::norm_sq(), libMesh::System::number(), libMesh::out, libMesh::EquationSystems::parameters, libMesh::MeshBase::point(), libMesh::MeshTools::Subdivision::prepare_subdivision_mesh(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::MeshBase::read(), libMesh::Real, libMesh::MeshTools::Modification::scale(), libMesh::Parameters::set(), libMesh::LinearImplicitSystem::solve(), libMesh::SUBDIVISION, TIMPI::Communicator::sum(), libMesh::MeshRefinement::uniformly_refine(), libMesh::System::variable_number(), libMesh::VTKIO::write(), libMesh::ExodusII_IO::write(), libMesh::MeshOutput< MT >::write_equation_systems(), and libMesh::ExodusII_IO::write_equation_systems().

{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);

  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
  libmesh_example_requires (3 == LIBMESH_DIM, "3D support");

  // Skip this example without --enable-node-valence
#ifndef LIBMESH_ENABLE_NODE_VALENCE
  libmesh_example_requires (false, "--enable-node-valence");
#endif

  // Skip this example without --enable-amr; requires MeshRefinement
#ifndef LIBMESH_ENABLE_AMR
  libmesh_example_requires(false, "--enable-amr");
#else

  // Skip this example without --enable-second; requires d2phi
#ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
  libmesh_example_requires(false, "--enable-second");
#elif LIBMESH_DIM > 2

  // Create a 2D mesh distributed across the default MPI communicator.
  // Subdivision surfaces do not appear to work with DistributedMesh yet.
  ReplicatedMesh mesh (init.comm(), 2);

  // Read the coarse square mesh.
  mesh.read ("square_mesh.off");

  // Resize the square plate to edge length L.
  const Real L = 100.;
  MeshTools::Modification::scale(mesh, L, L, L);

  // Get the number of mesh refinements from the command line
  const int n_refinements =
    libMesh::command_line_next("-n_refinements", 3);

  // Quadrisect the mesh triangles a few times to obtain a
  // finer mesh.  Subdivision surface elements require the
  // refinement data to be removed afterward.
  MeshRefinement mesh_refinement (mesh);
  mesh_refinement.uniformly_refine (n_refinements);
  MeshTools::Modification::flatten (mesh);

  // Write the mesh before the ghost elements are added.
#if defined(LIBMESH_HAVE_VTK)
  VTKIO(mesh).write ("without_ghosts.pvtu");
#endif
#if defined(LIBMESH_HAVE_EXODUS_API)
  ExodusII_IO(mesh).write ("without_ghosts.e");
#endif

  // Print information about the triangulated mesh to the screen.
  mesh.print_info();

  // Turn the triangulated mesh into a subdivision mesh
  // and add an additional row of "ghost" elements around
  // it in order to complete the extended local support of
  // the triangles at the boundaries.  If the second
  // argument is set to true, the outermost existing
  // elements are converted into ghost elements, and the
  // actual physical mesh is thus getting smaller.
  MeshTools::Subdivision::prepare_subdivision_mesh (mesh, false);

  // Print information about the subdivision mesh to the screen.
  mesh.print_info();

  // Write the mesh with the ghost elements added.
  // Compare this to the original mesh to see the difference.
#if defined(LIBMESH_HAVE_VTK)
  VTKIO(mesh).write ("with_ghosts.pvtu");
#endif
#if defined(LIBMESH_HAVE_EXODUS_API)
  ExodusII_IO(mesh).write ("with_ghosts.e");
#endif

  // Create an equation systems object.
  EquationSystems equation_systems (mesh);

  // Declare the system and its variables.
  // Create a linear implicit system named "Shell".
  LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem> ("Shell");

  // Add the three translational deformation variables
  // "u", "v", "w" to "Shell".  Since subdivision shell
  // elements meet the C1-continuity requirement, no
  // rotational or other auxiliary variables are needed.
  // Loop Subdivision Elements are always interpolated
  // by quartic box splines, hence the order must always
  // be FOURTH.
  system.add_variable ("u", FOURTH, SUBDIVISION);
  system.add_variable ("v", FOURTH, SUBDIVISION);
  system.add_variable ("w", FOURTH, SUBDIVISION);

  // Give the system a pointer to the matrix and rhs assembly
  // function.
  system.attach_assemble_function (assemble_shell);

  // Use the parameters of the equation systems object to
  // tell the shell system about the material properties, the
  // shell thickness, and the external load.
  const Real h  = 1.;
  const Real E  = 1.e7;
  const Real nu = 0.;
  const Real q  = 1.;
  equation_systems.parameters.set<Real> ("thickness")       = h;
  equation_systems.parameters.set<Real> ("young's modulus") = E;
  equation_systems.parameters.set<Real> ("poisson ratio")   = nu;
  equation_systems.parameters.set<Real> ("uniform load")    = q;

  // Initialize the data structures for the equation system.
  equation_systems.init();

  // Print information about the system to the screen.
  equation_systems.print_info();

  // Solve the linear system.
  system.solve();

  // After solving the system, write the solution to a VTK
  // or ExodusII output file ready for import in, e.g.,
  // Paraview.
#if defined(LIBMESH_HAVE_VTK)
  VTKIO(mesh).write_equation_systems ("out.pvtu", equation_systems);
#endif
#if defined(LIBMESH_HAVE_EXODUS_API)
  ExodusII_IO(mesh).write_equation_systems ("out.e", equation_systems);
#endif

  // Find the center node to measure the maximum deformation of the plate.
  Node * center_node = 0;
  Real nearest_dist_sq = mesh.point(0).norm_sq();
  for (unsigned int nid=1; nid<mesh.n_nodes(); ++nid)
    {
      const Real dist_sq = mesh.point(nid).norm_sq();
      if (dist_sq < nearest_dist_sq)
        {
          nearest_dist_sq = dist_sq;
          center_node = mesh.node_ptr(nid);
        }
    }

  // Finally, we evaluate the z-displacement "w" at the center node.
  const unsigned int w_var = system.variable_number ("w");
  dof_id_type w_dof = center_node->dof_number (system.number(), w_var, 0);
  Number w = 0;
  if (w_dof >= system.get_dof_map().first_dof() &&
      w_dof <  system.get_dof_map().end_dof())
    w = system.current_solution(w_dof);
  system.comm().sum(w);


  // The analytic solution for the maximum displacement of
  // a clamped square plate in pure bending, from Taylor,
  // Govindjee, Commun. Numer. Meth. Eng. 20, 757-765, 2004.
  const Real D = E * h*h*h / (12*(1-nu*nu));
  const Real w_analytic = 0.001265319 * L*L*L*L * q / D;

  // Print the finite element solution and the analytic
  // prediction of the maximum displacement of the clamped
  // square plate to the screen.
  libMesh::out << "z-displacement of the center point: " << w << std::endl;
  libMesh::out << "Analytic solution for pure bending: " << w_analytic << std::endl;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES, LIBMESH_DIM > 2

#endif // #ifdef LIBMESH_ENABLE_AMR

  // All done.
  return 0;
}