libMesh
Functions
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

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

Definition at line 326 of file miscellaneous_ex12.C.

References std::abs(), libMesh::MeshBase::active_local_element_ptr_range(), 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::DofObject::dof_number(), libMesh::Parameters::get(), libMesh::MeshBase::get_boundary_info(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::libmesh_ignore(), libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), 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().

328 {
329  // This example requires Eigen to actually work, but we should still
330  // let it compile and throw a runtime error if you don't.
331 #ifdef LIBMESH_HAVE_EIGEN
332  // It is a good idea to make sure we are assembling
333  // the proper system.
334  libmesh_assert_equal_to (system_name, "Shell");
335 
336  // Get a constant reference to the mesh object.
337  const MeshBase & mesh = es.get_mesh();
338  const unsigned int dim = mesh.mesh_dimension();
339 
340  // Get a reference to the shell system object.
341  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> (system_name);
342 
343  // Get the shell parameters that we need during assembly.
344  const Real h = es.parameters.get<Real> ("thickness");
345  const Real E = es.parameters.get<Real> ("young's modulus");
346  const Real nu = es.parameters.get<Real> ("poisson ratio");
347  const Real q = es.parameters.get<Real> ("point load");
348  const bool distributed_load = es.parameters.get<bool> ("distributed load");
349 
350  // The membrane elastic matrix.
351  Eigen::Matrix3d Hm;
352  Hm <<
353  1., nu, 0.,
354  nu, 1., 0.,
355  0., 0., 0.5 * (1-nu);
356  Hm *= h * E/(1-nu*nu);
357 
358  // The bending elastic matrix.
359  Eigen::Matrix3d Hf;
360  Hf <<
361  1., nu, 0.,
362  nu, 1., 0.,
363  0., 0., 0.5 * (1-nu);
364  Hf *= h*h*h/12 * E/(1-nu*nu);
365 
366  // The shear elastic matrices.
367  Eigen::Matrix2d Hc0 = Eigen::Matrix2d::Identity();
368  Hc0 *= h * 5./6*E/(2*(1+nu));
369 
370  Eigen::Matrix2d Hc1 = Eigen::Matrix2d::Identity();
371  Hc1 *= h*h*h/12 * 5./6*E/(2*(1+nu));
372 
373  // Get the Finite Element type, this will be
374  // the same for all variables.
375  FEType fe_type = system.variable_type (0);
376 
377  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
378  QGauss qrule (dim, fe_type.default_quadrature_order());
379  fe->attach_quadrature_rule (&qrule);
380 
381  // The element Jacobian * quadrature weight at each integration point.
382  const std::vector<Real> & JxW = fe->get_JxW();
383 
384  // The element shape function and its derivatives evaluated at the
385  // quadrature points.
386  const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
387  const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();
388  const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
389  const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
390  const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
391  const std::vector<std::vector<Real>> & dphidxi = fe->get_dphidxi();
392  const std::vector<std::vector<Real>> & dphideta = fe->get_dphideta();
393  const std::vector<std::vector<Real>> & phi = fe->get_phi();
394 
395  // A reference to the DofMap object for this system. The DofMap
396  // object handles the index translation from node and element numbers
397  // to degree of freedom numbers.
398  const DofMap & dof_map = system.get_dof_map();
399 
400  // Define data structures to contain the element stiffness matrix.
402  DenseSubMatrix<Number> Ke_var[6][6] =
403  {
416  };
417 
418  // Define data structures to contain the element rhs vector.
420  DenseSubVector<Number> Fe_w(Fe);
421 
422  std::vector<dof_id_type> dof_indices;
423  std::vector<std::vector<dof_id_type>> dof_indices_var(6);
424 
425  // Now we will loop over all the elements in the mesh. We will
426  // compute the element matrix and right-hand-side contribution.
427  for (const auto & elem : mesh.active_local_element_ptr_range())
428  {
429  dof_map.dof_indices (elem, dof_indices);
430  for (unsigned int var=0; var<6; var++)
431  dof_map.dof_indices (elem, dof_indices_var[var], var);
432 
433  const unsigned int n_dofs = dof_indices.size();
434  const unsigned int n_var_dofs = dof_indices_var[0].size();
435 
436  // First compute element data at the nodes
437  std::vector<Point> nodes;
438  for (auto i : elem->node_index_range())
439  nodes.push_back(elem->reference_elem()->node_ref(i));
440  fe->reinit (elem, &nodes);
441 
442  // Convenient notation for the element node positions
443  Eigen::Vector3d X1(elem->node_ref(0)(0), elem->node_ref(0)(1), elem->node_ref(0)(2));
444  Eigen::Vector3d X2(elem->node_ref(1)(0), elem->node_ref(1)(1), elem->node_ref(1)(2));
445  Eigen::Vector3d X3(elem->node_ref(2)(0), elem->node_ref(2)(1), elem->node_ref(2)(2));
446  Eigen::Vector3d X4(elem->node_ref(3)(0), elem->node_ref(3)(1), elem->node_ref(3)(2));
447 
448  //Store covariant basis and local orthonormal basis at the nodes
449  std::vector<Eigen::Matrix3d> F0node;
450  std::vector<Eigen::Matrix3d> Qnode;
451  for (auto i : elem->node_index_range())
452  {
453  Eigen::Vector3d a1;
454  a1 << dxyzdxi[i](0), dxyzdxi[i](1), dxyzdxi[i](2);
455  Eigen::Vector3d a2;
456  a2 << dxyzdeta[i](0), dxyzdeta[i](1), dxyzdeta[i](2);
457  Eigen::Vector3d n;
458  n = a1.cross(a2);
459  n /= n.norm();
460  Eigen::Matrix3d F0;
461  F0 <<
462  a1(0), a2(0), n(0),
463  a1(1), a2(1), n(1),
464  a1(2), a2(2), n(2);
465  F0node.push_back(F0);
466 
467  Real nx = n(0);
468  Real ny = n(1);
469  Real C = n(2);
470  if (std::abs(1.+C)<1e-6)
471  {
472  Eigen::Matrix3d Q;
473  Q <<
474  1, 0, 0,
475  0, -1, 0,
476  0, 0, -1;
477  Qnode.push_back(Q);
478  }
479  else
480  {
481  Eigen::Matrix3d Q;
482  Q <<
483  C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
484  -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
485  -nx, -ny, C;
486  Qnode.push_back(Q);
487  }
488  }
489 
490  Ke.resize (n_dofs, n_dofs);
491  for (unsigned int var_i=0; var_i<6; var_i++)
492  for (unsigned int var_j=0; var_j<6; var_j++)
493  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
494 
495  Fe.resize(n_dofs);
496  Fe_w.reposition(2*n_var_dofs,n_var_dofs);
497 
498  // Reinit element data at the regular Gauss quadrature points
499  fe->reinit (elem);
500 
501  // Now we will build the element matrix and right-hand-side.
502  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
503  {
504 
505  //Covariant basis at the quadrature point
506  Eigen::Vector3d a1;
507  a1 << dxyzdxi[qp](0), dxyzdxi[qp](1), dxyzdxi[qp](2);
508  Eigen::Vector3d a2;
509  a2 << dxyzdeta[qp](0), dxyzdeta[qp](1), dxyzdeta[qp](2);
510  Eigen::Vector3d n;
511  n = a1.cross(a2);
512  n /= n.norm();
513  Eigen::Matrix3d F0;
514  F0 <<
515  a1(0), a2(0), n(0),
516  a1(1), a2(1), n(1),
517  a1(2), a2(2), n(2);
518 
519  //Contravariant basis
520  Eigen::Matrix3d F0it;
521  F0it = F0.inverse().transpose();
522 
523  //Local orthonormal basis at the quadrature point
524  Real nx = n(0);
525  Real ny = n(1);
526  Real C = n(2);
527  Eigen::Matrix3d Q;
528  if (std::abs(1.+C) < 1e-6)
529  {
530  Q <<
531  1, 0, 0,
532  0, -1, 0,
533  0, 0, -1;
534  }
535  else
536  {
537  Q <<
538  C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
539  -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
540  -nx, -ny, C;
541  }
542 
543  Eigen::Matrix2d C0;
544  C0 = F0it.block<3,2>(0,0).transpose()*Q.block<3,2>(0,0);
545 
546  // Normal derivatives in reference coordinates
547  Eigen::Vector3d d2Xdxi2(d2xyzdxi2[qp](0), d2xyzdxi2[qp](1), d2xyzdxi2[qp](2));
548  Eigen::Vector3d d2Xdeta2(d2xyzdeta2[qp](0), d2xyzdeta2[qp](1), d2xyzdeta2[qp](2));
549  Eigen::Vector3d d2Xdxideta(d2xyzdxideta[qp](0), d2xyzdxideta[qp](1), d2xyzdxideta[qp](2));
550 
551  Eigen::Matrix2d b;
552  b <<
553  n.dot(d2Xdxi2), n.dot(d2Xdxideta),
554  n.dot(d2Xdxideta), n.dot(d2Xdeta2);
555 
556  Eigen::Vector3d dndxi = -b(0,0)*F0it.col(0) - b(0,1)*F0it.col(1);
557  Eigen::Vector3d dndeta = -b(1,0)*F0it.col(0) - b(1,1)*F0it.col(1);
558 
559  Eigen::Matrix2d bhat;
560  bhat <<
561  F0it.col(1).dot(dndeta), -F0it.col(0).dot(dndeta),
562  -F0it.col(1).dot(dndxi), F0it.col(0).dot(dndxi);
563 
564  Eigen::Matrix2d bc;
565  bc = bhat*C0;
566 
567  // Mean curvature
568  Real H = 0.5*(dndxi.dot(F0it.col(0))+dndeta.dot(F0it.col(1)));
569 
570  // Quadrature point reference coordinates
571  Real xi = qrule.qp(qp)(0);
572  Real eta = qrule.qp(qp)(1);
573 
574  // Preassemble the MITC4 shear strain matrix for all nodes as they involve
575  // cross references to midside nodes.
576  // The QUAD4 element has nodes X1,X2,X3,X4 with coordinates (xi,eta)
577  // in the reference element: (-1,-1),(1,-1),(1,1),(-1,1).
578  // The midside nodes are denoted A1=(X1+X2)/2, B2=(X2+X3)/2, A2=(X3+X4)/2, B1=(X4+X1)/2.
579 
580  // Normals at the midside nodes (average of normals at the edge corners).
581  // Multiplication by the assumed shear strain shape function.
582  Eigen::Vector3d nA1 = 0.5*(Qnode[0].col(2)+Qnode[1].col(2));
583  nA1 /= nA1.norm();
584  nA1 *= (1-eta)/4;
585  Eigen::Vector3d nB2 = 0.5*(Qnode[1].col(2)+Qnode[2].col(2));
586  nB2 /= nB2.norm();
587  nB2 *= (1+xi)/4;
588  Eigen::Vector3d nA2 = 0.5*(Qnode[2].col(2)+Qnode[3].col(2));
589  nA2 /= nA2.norm();
590  nA2 *= (1+eta)/4;
591  Eigen::Vector3d nB1 = 0.5*(Qnode[3].col(2)+Qnode[0].col(2));
592  nB1 /= nB1.norm();
593  nB1 *= (1-xi)/4;
594 
595  // Edge tangents
596  Eigen::Vector3d aA1 = 0.5*(X2-X1);
597  Eigen::Vector3d aA2 = 0.5*(X3-X4);
598  Eigen::Vector3d aB1 = 0.5*(X4-X1);
599  Eigen::Vector3d aB2 = 0.5*(X3-X2);
600 
601  // Contribution of the rotational dofs to the shear strain
602  Eigen::Vector2d AS1A1(-aA1.dot(Qnode[0].col(1)), aA1.dot(Qnode[0].col(0)));
603  Eigen::Vector2d AS2A1(-aA1.dot(Qnode[1].col(1)), aA1.dot(Qnode[1].col(0)));
604  AS1A1 *= (1-eta)/4;
605  AS2A1 *= (1-eta)/4;
606 
607  Eigen::Vector2d AS1A2(-aA2.dot(Qnode[3].col(1)), aA2.dot(Qnode[3].col(0)));
608  Eigen::Vector2d AS2A2(-aA2.dot(Qnode[2].col(1)), aA2.dot(Qnode[2].col(0)));
609  AS1A2 *= (1+eta)/4;
610  AS2A2 *= (1+eta)/4;
611 
612  Eigen::Vector2d AS1B1(-aB1.dot(Qnode[0].col(1)), aB1.dot(Qnode[0].col(0)));
613  Eigen::Vector2d AS2B1(-aB1.dot(Qnode[3].col(1)), aB1.dot(Qnode[3].col(0)));
614  AS1B1 *= (1-xi)/4;
615  AS2B1 *= (1-xi)/4;
616 
617  Eigen::Vector2d AS1B2(-aB2.dot(Qnode[1].col(1)), aB2.dot(Qnode[1].col(0)));
618  Eigen::Vector2d AS2B2(-aB2.dot(Qnode[2].col(1)), aB2.dot(Qnode[2].col(0)));
619  AS1B2 *= (1+xi)/4;
620  AS2B2 *= (1+xi)/4;
621 
622  // Store previous quantities in the shear strain matrices for each node
623  std::vector<Eigen::MatrixXd> Bcnode;
624  Eigen::MatrixXd Bc(2, 5);
625  // Node 1
626  Bc.block<1,3>(0,0) = -nA1.transpose();
627  Bc.block<1,2>(0,3) = AS1A1.transpose();
628  Bc.block<1,3>(1,0) = -nB1.transpose();
629  Bc.block<1,2>(1,3) = AS1B1.transpose();
630  Bcnode.push_back(Bc);
631  // Node 2
632  Bc.block<1,3>(0,0) = nA1.transpose();
633  Bc.block<1,2>(0,3) = AS2A1.transpose();
634  Bc.block<1,3>(1,0) = -nB2.transpose();
635  Bc.block<1,2>(1,3) = AS1B2.transpose();
636  Bcnode.push_back(Bc);
637  // Node 3
638  Bc.block<1,3>(0,0) = nA2.transpose();
639  Bc.block<1,2>(0,3) = AS2A2.transpose();
640  Bc.block<1,3>(1,0) = nB2.transpose();
641  Bc.block<1,2>(1,3) = AS2B2.transpose();
642  Bcnode.push_back(Bc);
643  // Node 4
644  Bc.block<1,3>(0,0) = -nA2.transpose();
645  Bc.block<1,2>(0,3) = AS1A2.transpose();
646  Bc.block<1,3>(1,0) = nB1.transpose();
647  Bc.block<1,2>(1,3) = AS2B1.transpose();
648  Bcnode.push_back(Bc);
649 
650  // Loop over all pairs of nodes I,J.
651  for (unsigned int i=0; i<n_var_dofs; ++i)
652  {
653  // Matrix B0, zeroth order (through thickness) membrane-bending strain
654  Real C1i = dphidxi[i][qp]*C0(0,0) + dphideta[i][qp]*C0(1,0);
655  Real C2i = dphidxi[i][qp]*C0(0,1) + dphideta[i][qp]*C0(1,1);
656 
657  Eigen::MatrixXd B0I(3, 5);
658  B0I = Eigen::MatrixXd::Zero(3, 5);
659  B0I.block<1,3>(0,0) = C1i*Q.col(0).transpose();
660  B0I.block<1,3>(1,0) = C2i*Q.col(1).transpose();
661  B0I.block<1,3>(2,0) = C2i*Q.col(0).transpose()+C1i*Q.col(1).transpose();
662 
663  // Matrix B1, first order membrane-bending strain
664  Real bc1i = dphidxi[i][qp]*bc(0,0) + dphideta[i][qp]*bc(1,0);
665  Real bc2i = dphidxi[i][qp]*bc(0,1) + dphideta[i][qp]*bc(1,1);
666 
667  Eigen::Vector2d V1i(-Q.col(0).dot(Qnode[i].col(1)),
668  Q.col(0).dot(Qnode[i].col(0)));
669 
670  Eigen::Vector2d V2i(-Q.col(1).dot(Qnode[i].col(1)),
671  Q.col(1).dot(Qnode[i].col(0)));
672 
673  Eigen::MatrixXd B1I(3,5);
674  B1I = Eigen::MatrixXd::Zero(3,5);
675  B1I.block<1,3>(0,0) = bc1i*Q.col(0).transpose();
676  B1I.block<1,3>(1,0) = bc2i*Q.col(1).transpose();
677  B1I.block<1,3>(2,0) = bc2i*Q.col(0).transpose()+bc1i*Q.col(1).transpose();
678 
679  B1I.block<1,2>(0,3) = C1i*V1i.transpose();
680  B1I.block<1,2>(1,3) = C2i*V2i.transpose();
681  B1I.block<1,2>(2,3) = C2i*V1i.transpose()+C1i*V2i.transpose();
682 
683  // Matrix B2, second order membrane-bending strain
684  Eigen::MatrixXd B2I(3,5);
685  B2I = Eigen::MatrixXd::Zero(3,5);
686 
687  B2I.block<1,2>(0,3) = bc1i*V1i.transpose();
688  B2I.block<1,2>(1,3) = bc2i*V2i.transpose();
689  B2I.block<1,2>(2,3) = bc2i*V1i.transpose()+bc1i*V2i.transpose();
690 
691  // Matrix Bc0, zeroth order shear strain
692  Eigen::MatrixXd Bc0I(2,5);
693  Bc0I = C0.transpose()*Bcnode[i];
694 
695  // Matrix Bc1, first order shear strain
696  Eigen::MatrixXd Bc1I(2,5);
697  Bc1I = bc.transpose()*Bcnode[i];
698 
699  // Drilling dof (in-plane rotation)
700  Eigen::Vector2d BdxiI(dphidxi[i][qp],dphideta[i][qp]);
701  Eigen::Vector2d BdI = C0.transpose()*BdxiI;
702 
703  for (unsigned int j=0; j<n_var_dofs; ++j)
704  {
705 
706  // Matrix B0, zeroth order membrane-bending strain
707  Real C1j = dphidxi[j][qp]*C0(0,0) + dphideta[j][qp]*C0(1,0);
708  Real C2j = dphidxi[j][qp]*C0(0,1) + dphideta[j][qp]*C0(1,1);
709 
710  Eigen::MatrixXd B0J(3,5);
711  B0J = Eigen::MatrixXd::Zero(3,5);
712  B0J.block<1,3>(0,0) = C1j*Q.col(0).transpose();
713  B0J.block<1,3>(1,0) = C2j*Q.col(1).transpose();
714  B0J.block<1,3>(2,0) = C2j*Q.col(0).transpose()+C1j*Q.col(1).transpose();
715 
716  // Matrix B1, first order membrane-bending strain
717  Real bc1j = dphidxi[j][qp]*bc(0,0) + dphideta[j][qp]*bc(1,0);
718  Real bc2j = dphidxi[j][qp]*bc(0,1) + dphideta[j][qp]*bc(1,1);
719 
720  Eigen::Vector2d V1j(-Q.col(0).dot(Qnode[j].col(1)),
721  Q.col(0).dot(Qnode[j].col(0)));
722 
723  Eigen::Vector2d V2j(-Q.col(1).dot(Qnode[j].col(1)),
724  Q.col(1).dot(Qnode[j].col(0)));
725 
726  Eigen::MatrixXd B1J(3,5);
727  B1J = Eigen::MatrixXd::Zero(3,5);
728  B1J.block<1,3>(0,0) = bc1j*Q.col(0).transpose();
729  B1J.block<1,3>(1,0) = bc2j*Q.col(1).transpose();
730  B1J.block<1,3>(2,0) = bc2j*Q.col(0).transpose()+bc1j*Q.col(1).transpose();
731 
732  B1J.block<1,2>(0,3) = C1j*V1j.transpose();
733  B1J.block<1,2>(1,3) = C2j*V2j.transpose();
734  B1J.block<1,2>(2,3) = C2j*V1j.transpose()+C1j*V2j.transpose();
735 
736  // Matrix B2, second order membrane-bending strain
737  Eigen::MatrixXd B2J(3,5);
738  B2J = Eigen::MatrixXd::Zero(3,5);
739 
740  B2J.block<1,2>(0,3) = bc1j*V1j.transpose();
741  B2J.block<1,2>(1,3) = bc2j*V2j.transpose();
742  B2J.block<1,2>(2,3) = bc2j*V1j.transpose()+bc1j*V2j.transpose();
743 
744  // Matrix Bc0, zeroth order shear strain
745  Eigen::MatrixXd Bc0J(2, 5);
746  Bc0J = C0.transpose()*Bcnode[j];
747 
748  // Matrix Bc1, first order shear strain
749  Eigen::MatrixXd Bc1J(2, 5);
750  Bc1J = bc.transpose()*Bcnode[j];
751 
752  // Drilling dof
753  Eigen::Vector2d BdxiJ(dphidxi[j][qp], dphideta[j][qp]);
754  Eigen::Vector2d BdJ = C0.transpose()*BdxiJ;
755 
756  // The total stiffness matrix coupling the nodes
757  // I and J is a sum of membrane, bending and shear contributions.
758  Eigen::MatrixXd local_KIJ(5, 5);
759  local_KIJ = JxW[qp] * (
760  B0I.transpose() * Hm * B0J
761  + B2I.transpose() * Hf * B0J
762  + B0I.transpose() * Hf * B2J
763  + B1I.transpose() * Hf * B1J
764  + 2*H * B0I.transpose() * Hf * B1J
765  + 2*H * B1I.transpose() * Hf * B0J
766  + Bc0I.transpose() * Hc0 * Bc0J
767  + Bc1I.transpose() * Hc1 * Bc1J
768  + 2*H * Bc0I.transpose() * Hc1 * Bc1J
769  + 2*H * Bc1I.transpose() * Hc1 * Bc0J
770  );
771 
772  // Going from 5 to 6 dofs to add drilling dof
773  Eigen::MatrixXd full_local_KIJ(6, 6);
774  full_local_KIJ = Eigen::MatrixXd::Zero(6, 6);
775  full_local_KIJ.block<5,5>(0,0)=local_KIJ;
776 
777  // Drilling dof stiffness contribution
778  full_local_KIJ(5,5) = Hf(0,0)*JxW[qp]*BdI.transpose()*BdJ;
779 
780  // Transform the stiffness matrix to global coordinates
781  Eigen::MatrixXd global_KIJ(6,6);
782  Eigen::MatrixXd TI(6,6);
783  TI = Eigen::MatrixXd::Identity(6,6);
784  TI.block<3,3>(3,3) = Qnode[i].transpose();
785  Eigen::MatrixXd TJ(6,6);
786  TJ = Eigen::MatrixXd::Identity(6,6);
787  TJ.block<3,3>(3,3) = Qnode[j].transpose();
788  global_KIJ = TI.transpose()*full_local_KIJ*TJ;
789 
790  // Insert the components of the coupling stiffness
791  // matrix KIJ into the corresponding directional
792  // submatrices.
793  for (unsigned int k=0;k<6;k++)
794  for (unsigned int l=0;l<6;l++)
795  Ke_var[k][l](i,j) += global_KIJ(k,l);
796  }
797  }
798 
799  } // end of the quadrature point qp-loop
800 
801  if (distributed_load)
802  {
803  // Loop on shell faces
804  for (unsigned int shellface=0; shellface<2; shellface++)
805  {
806  std::vector<boundary_id_type> bids;
807  mesh.get_boundary_info().shellface_boundary_ids(elem, shellface, bids);
808 
809  for (std::size_t k=0; k<bids.size(); k++)
810  if (bids[k]==11) // sideset id for surface load
811  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
812  for (unsigned int i=0; i<n_var_dofs; ++i)
813  Fe_w(i) -= JxW[qp] * phi[i][qp];
814  }
815  }
816 
817  // The element matrix is now built for this element.
818  // Add it to the global matrix.
819 
820  dof_map.constrain_element_matrix_and_vector (Ke,Fe,dof_indices);
821 
822  system.matrix->add_matrix (Ke, dof_indices);
823  system.rhs->add_vector (Fe, dof_indices);
824 
825  }
826 
827  if (!distributed_load)
828  {
829  //Adding point load to the RHS
830 
831  //Pinch position
832  Point C(0, 3, 3);
833 
834  //Finish assembling rhs so we can set one value
835  system.rhs->close();
836 
838  const MeshBase::const_node_iterator node_end = mesh.nodes_end();
839 
840  for ( ; nodeit!=node_end; ++nodeit)
841  {
842  Node & node = **nodeit;
843  if ((node-C).norm() < 1e-3)
844  system.rhs->set(node.dof_number(0, 2, 0), -q/4);
845  }
846  }
847 
848 #else
849  // Avoid compiler warnings
850  libmesh_ignore(es);
851  libmesh_ignore(system_name);
852 #endif // LIBMESH_HAVE_EIGEN
853 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
double abs(double a)
A Node is like a Point, but with more information.
Definition: node.h:52
unsigned int dim
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
const T & get(const std::string &) const
Definition: parameters.h:430
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
virtual node_iterator nodes_begin()=0
Iterate over all the nodes in the Mesh.
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const DofMap & get_dof_map() const
Definition: system.h:2030
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
Order default_quadrature_order() const
Definition: fe_type.h:332
Defines a dense submatrix for use in Finite Element-type computations.
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
void libmesh_ignore(const T &)
Definition: assembly.h:127
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
The definition of the const_node_iterator struct.
Definition: mesh_base.h:1548
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
void assemble_shell ( EquationSystems es,
const std::string &  libmesh_dbg_varsystem_name 
)

Definition at line 257 of file miscellaneous_ex11.C.

References libMesh::MeshBase::active_local_element_ptr_range(), 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_nullptr, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshTools::Subdivision::next, 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().

259 {
260  // It is a good idea to make sure we are assembling
261  // the proper system.
262  libmesh_assert_equal_to (system_name, "Shell");
263 
264  // Get a constant reference to the mesh object.
265  const MeshBase & mesh = es.get_mesh();
266 
267  // Get a reference to the shell system object.
268  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> ("Shell");
269 
270  // Get the shell parameters that we need during assembly.
271  const Real h = es.parameters.get<Real> ("thickness");
272  const Real E = es.parameters.get<Real> ("young's modulus");
273  const Real nu = es.parameters.get<Real> ("poisson ratio");
274  const Real q = es.parameters.get<Real> ("uniform load");
275 
276  // Compute the membrane stiffness K and the bending
277  // rigidity D from these parameters.
278  const Real K = E * h / (1-nu*nu);
279  const Real D = E * h*h*h / (12*(1-nu*nu));
280 
281  // Numeric ids corresponding to each variable in the system.
282  const unsigned int u_var = system.variable_number ("u");
283  const unsigned int v_var = system.variable_number ("v");
284  const unsigned int w_var = system.variable_number ("w");
285 
286  // Get the Finite Element type for "u". Note this will be
287  // the same as the type for "v" and "w".
288  FEType fe_type = system.variable_type (u_var);
289 
290  // Build a Finite Element object of the specified type.
291  UniquePtr<FEBase> fe (FEBase::build(2, fe_type));
292 
293  // A Gauss quadrature rule for numerical integration.
294  // For subdivision shell elements, a single Gauss point per
295  // element is sufficient, hence we use extraorder = 0.
296  const int extraorder = 0;
297  UniquePtr<QBase> qrule (fe_type.default_quadrature_rule (2, extraorder));
298 
299  // Tell the finite element object to use our quadrature rule.
300  fe->attach_quadrature_rule (qrule.get());
301 
302  // The element Jacobian * quadrature weight at each integration point.
303  const std::vector<Real> & JxW = fe->get_JxW();
304 
305  // The surface tangents in both directions at the quadrature points.
306  const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
307  const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();
308 
309  // The second partial derivatives at the quadrature points.
310  const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
311  const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
312  const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
313 
314  // The element shape function and its derivatives evaluated at the
315  // quadrature points.
316  const std::vector<std::vector<Real>> & phi = fe->get_phi();
317  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
318  const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
319 
320  // A reference to the DofMap object for this system. The DofMap
321  // object handles the index translation from node and element numbers
322  // to degree of freedom numbers.
323  const DofMap & dof_map = system.get_dof_map();
324 
325  // Define data structures to contain the element stiffness matrix
326  // and right-hand-side vector contribution. Following
327  // basic finite element terminology we will denote these
328  // "Ke" and "Fe".
331 
333  Kuu(Ke), Kuv(Ke), Kuw(Ke),
334  Kvu(Ke), Kvv(Ke), Kvw(Ke),
335  Kwu(Ke), Kwv(Ke), Kww(Ke);
336 
338  Fu(Fe),
339  Fv(Fe),
340  Fw(Fe);
341 
342  // This vector will hold the degree of freedom indices for
343  // the element. These define where in the global system
344  // the element degrees of freedom get mapped.
345  std::vector<dof_id_type> dof_indices;
346  std::vector<dof_id_type> dof_indices_u;
347  std::vector<dof_id_type> dof_indices_v;
348  std::vector<dof_id_type> dof_indices_w;
349 
350  // Now we will loop over all the elements in the mesh. We will
351  // compute the element matrix and right-hand-side contribution.
352  for (const auto & elem : mesh.active_local_element_ptr_range())
353  {
354  // The ghost elements at the boundaries need to be excluded
355  // here, as they don't belong to the physical shell,
356  // but serve for a proper boundary treatment only.
357  libmesh_assert_equal_to (elem->type(), TRI3SUBDIVISION);
358  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *> (elem);
359  if (sd_elem->is_ghost())
360  continue;
361 
362  // Get the degree of freedom indices for the
363  // current element. These define where in the global
364  // matrix and right-hand-side this element will
365  // contribute to.
366  dof_map.dof_indices (elem, dof_indices);
367  dof_map.dof_indices (elem, dof_indices_u, u_var);
368  dof_map.dof_indices (elem, dof_indices_v, v_var);
369  dof_map.dof_indices (elem, dof_indices_w, w_var);
370 
371  const std::size_t n_dofs = dof_indices.size();
372  const std::size_t n_u_dofs = dof_indices_u.size();
373  const std::size_t n_v_dofs = dof_indices_v.size();
374  const std::size_t n_w_dofs = dof_indices_w.size();
375 
376  // Compute the element-specific data for the current
377  // element. This involves computing the location of the
378  // quadrature points and the shape functions
379  // (phi, dphi, d2phi) for the current element.
380  fe->reinit (elem);
381 
382  // Zero the element matrix and right-hand side before
383  // summing them. We use the resize member here because
384  // the number of degrees of freedom might have changed from
385  // the last element.
386  Ke.resize (n_dofs, n_dofs);
387  Fe.resize (n_dofs);
388 
389  // Reposition the submatrices... The idea is this:
390  //
391  // - - - -
392  // | Kuu Kuv Kuw | | Fu |
393  // Ke = | Kvu Kvv Kvw |; Fe = | Fv |
394  // | Kwu Kwv Kww | | Fw |
395  // - - - -
396  //
397  // The DenseSubMatrix.reposition () member takes the
398  // (row_offset, column_offset, row_size, column_size).
399  //
400  // Similarly, the DenseSubVector.reposition () member
401  // takes the (row_offset, row_size)
402  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
403  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
404  Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
405 
406  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
407  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
408  Kvw.reposition (v_var*n_v_dofs, w_var*n_v_dofs, n_v_dofs, n_w_dofs);
409 
410  Kwu.reposition (w_var*n_w_dofs, u_var*n_w_dofs, n_w_dofs, n_u_dofs);
411  Kwv.reposition (w_var*n_w_dofs, v_var*n_w_dofs, n_w_dofs, n_v_dofs);
412  Kww.reposition (w_var*n_w_dofs, w_var*n_w_dofs, n_w_dofs, n_w_dofs);
413 
414  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
415  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
416  Fw.reposition (w_var*n_u_dofs, n_w_dofs);
417 
418  // Now we will build the element matrix and right-hand-side.
419  for (unsigned int qp=0; qp<qrule->n_points(); ++qp)
420  {
421  // First, we compute the external force resulting
422  // from a load q distributed uniformly across the plate.
423  // Since the load is supposed to be transverse to the plate,
424  // it affects the z-direction, i.e. the "w" variable.
425  for (unsigned int i=0; i<n_u_dofs; ++i)
426  Fw(i) += JxW[qp] * phi[i][qp] * q;
427 
428  // Next, we assemble the stiffness matrix. This is only valid
429  // for the linear theory, i.e., for small deformations, where
430  // reference and deformed surface metrics are indistinguishable.
431 
432  // Get the three surface basis vectors.
433  const RealVectorValue & a1 = dxyzdxi[qp];
434  const RealVectorValue & a2 = dxyzdeta[qp];
435  RealVectorValue a3 = a1.cross(a2);
436  const Real jac = a3.norm(); // the surface Jacobian
437  libmesh_assert_greater (jac, 0);
438  a3 /= jac; // the shell director a3 is normalized to unit length
439 
440  // Get the derivatives of the surface tangents.
441  const RealVectorValue & a11 = d2xyzdxi2[qp];
442  const RealVectorValue & a22 = d2xyzdeta2[qp];
443  const RealVectorValue & a12 = d2xyzdxideta[qp];
444 
445  // Compute the three covariant components of the first
446  // fundamental form of the surface.
447  const RealVectorValue a(a1*a1, a2*a2, a1*a2);
448 
449  // The elastic H matrix in Voigt's notation, computed from the
450  // covariant components of the first fundamental form rather
451  // than the contravariant components, exploiting that the
452  // contravariant first fundamental form is the inverse of the
453  // covariant first fundamental form (hence the determinant etc.).
454  RealTensorValue H;
455  H(0,0) = a(1) * a(1);
456  H(0,1) = H(1,0) = nu * a(1) * a(0) + (1-nu) * a(2) * a(2);
457  H(0,2) = H(2,0) = -a(1) * a(2);
458  H(1,1) = a(0) * a(0);
459  H(1,2) = H(2,1) = -a(0) * a(2);
460  H(2,2) = 0.5 * ((1-nu) * a(1) * a(0) + (1+nu) * a(2) * a(2));
461  const Real det = a(0) * a(1) - a(2) * a(2);
462  libmesh_assert_not_equal_to (det * det, 0);
463  H /= det * det;
464 
465  // Precompute come cross products for the bending part below.
466  const RealVectorValue a11xa2 = a11.cross(a2);
467  const RealVectorValue a22xa2 = a22.cross(a2);
468  const RealVectorValue a12xa2 = a12.cross(a2);
469  const RealVectorValue a1xa11 = a1.cross(a11);
470  const RealVectorValue a1xa22 = a1.cross(a22);
471  const RealVectorValue a1xa12 = a1.cross(a12);
472  const RealVectorValue a2xa3 = a2.cross(a3);
473  const RealVectorValue a3xa1 = a3.cross(a1);
474 
475  // Loop over all pairs of nodes I,J.
476  for (unsigned int i=0; i<n_u_dofs; ++i)
477  {
478  for (unsigned int j=0; j<n_u_dofs; ++j)
479  {
480  // The membrane strain matrices in Voigt's notation.
481  RealTensorValue MI, MJ;
482  for (unsigned int k=0; k<3; ++k)
483  {
484  MI(0,k) = dphi[i][qp](0) * a1(k);
485  MI(1,k) = dphi[i][qp](1) * a2(k);
486  MI(2,k) = dphi[i][qp](1) * a1(k)
487  + dphi[i][qp](0) * a2(k);
488 
489  MJ(0,k) = dphi[j][qp](0) * a1(k);
490  MJ(1,k) = dphi[j][qp](1) * a2(k);
491  MJ(2,k) = dphi[j][qp](1) * a1(k)
492  + dphi[j][qp](0) * a2(k);
493  }
494 
495  // The bending strain matrices in Voigt's notation.
496  RealTensorValue BI, BJ;
497  for (unsigned int k=0; k<3; ++k)
498  {
499  const Real term_ik = dphi[i][qp](0) * a2xa3(k)
500  + dphi[i][qp](1) * a3xa1(k);
501  BI(0,k) = -d2phi[i][qp](0,0) * a3(k)
502  +(dphi[i][qp](0) * a11xa2(k)
503  + dphi[i][qp](1) * a1xa11(k)
504  + (a3*a11) * term_ik) / jac;
505  BI(1,k) = -d2phi[i][qp](1,1) * a3(k)
506  +(dphi[i][qp](0) * a22xa2(k)
507  + dphi[i][qp](1) * a1xa22(k)
508  + (a3*a22) * term_ik) / jac;
509  BI(2,k) = 2 * (-d2phi[i][qp](0,1) * a3(k)
510  +(dphi[i][qp](0) * a12xa2(k)
511  + dphi[i][qp](1) * a1xa12(k)
512  + (a3*a12) * term_ik) / jac);
513 
514  const Real term_jk = dphi[j][qp](0) * a2xa3(k)
515  + dphi[j][qp](1) * a3xa1(k);
516  BJ(0,k) = -d2phi[j][qp](0,0) * a3(k)
517  +(dphi[j][qp](0) * a11xa2(k)
518  + dphi[j][qp](1) * a1xa11(k)
519  + (a3*a11) * term_jk) / jac;
520  BJ(1,k) = -d2phi[j][qp](1,1) * a3(k)
521  +(dphi[j][qp](0) * a22xa2(k)
522  + dphi[j][qp](1) * a1xa22(k)
523  + (a3*a22) * term_jk) / jac;
524  BJ(2,k) = 2 * (-d2phi[j][qp](0,1) * a3(k)
525  +(dphi[j][qp](0) * a12xa2(k)
526  + dphi[j][qp](1) * a1xa12(k)
527  + (a3*a12) * term_jk) / jac);
528  }
529 
530  // The total stiffness matrix coupling the nodes
531  // I and J is a sum of membrane and bending
532  // contributions according to the following formula.
533  const RealTensorValue KIJ = JxW[qp] * K * MI.transpose() * H * MJ
534  + JxW[qp] * D * BI.transpose() * H * BJ;
535 
536  // Insert the components of the coupling stiffness
537  // matrix KIJ into the corresponding directional
538  // submatrices.
539  Kuu(i,j) += KIJ(0,0);
540  Kuv(i,j) += KIJ(0,1);
541  Kuw(i,j) += KIJ(0,2);
542 
543  Kvu(i,j) += KIJ(1,0);
544  Kvv(i,j) += KIJ(1,1);
545  Kvw(i,j) += KIJ(1,2);
546 
547  Kwu(i,j) += KIJ(2,0);
548  Kwv(i,j) += KIJ(2,1);
549  Kww(i,j) += KIJ(2,2);
550  }
551  }
552 
553  } // end of the quadrature point qp-loop
554 
555  // The element matrix and right-hand-side are now built
556  // for this element. Add them to the global matrix and
557  // right-hand-side vector. The NumericMatrix::add_matrix()
558  // and NumericVector::add_vector() members do this for us.
559  system.matrix->add_matrix (Ke, dof_indices);
560  system.rhs->add_vector (Fe, dof_indices);
561  } // end of non-ghost element loop
562 
563  // Next, we apply the boundary conditions. In this case,
564  // all boundaries are clamped by the penalty method, using
565  // the special "ghost" nodes along the boundaries. Note
566  // that there are better ways to implement boundary conditions
567  // for subdivision shells. We use the simplest way here,
568  // which is known to be overly restrictive and will lead to
569  // a slightly too small deformation of the plate.
570  for (const auto & elem : mesh.active_local_element_ptr_range())
571  {
572  // For the boundary conditions, we only need to loop over
573  // the ghost elements.
574  libmesh_assert_equal_to (elem->type(), TRI3SUBDIVISION);
575  const Tri3Subdivision * gh_elem = static_cast<const Tri3Subdivision *> (elem);
576  if (!gh_elem->is_ghost())
577  continue;
578 
579  // Find the side which is part of the physical plate boundary,
580  // that is, the boundary of the original mesh without ghosts.
581  for (auto s : elem->side_index_range())
582  {
583  const Tri3Subdivision * nb_elem = static_cast<const Tri3Subdivision *> (elem->neighbor_ptr(s));
584  if (nb_elem == libmesh_nullptr || nb_elem->is_ghost())
585  continue;
586 
587  /*
588  * Determine the four nodes involved in the boundary
589  * condition treatment of this side. The MeshTools::Subdiv
590  * namespace provides lookup tables next and prev
591  * for an efficient determination of the next and previous
592  * nodes of an element, respectively.
593  *
594  * n4
595  * / \
596  * / gh \
597  * n2 ---- n3
598  * \ nb /
599  * \ /
600  * n1
601  */
602  const Node * nodes [4]; // n1, n2, n3, n4
603  nodes[1] = gh_elem->node_ptr(s); // n2
604  nodes[2] = gh_elem->node_ptr(MeshTools::Subdivision::next[s]); // n3
605  nodes[3] = gh_elem->node_ptr(MeshTools::Subdivision::prev[s]); // n4
606 
607  // The node in the interior of the domain, n1, is the
608  // hardest to find. Walk along the edges of element nb until
609  // we have identified it.
610  unsigned int n_int = 0;
611  nodes[0] = nb_elem->node_ptr(0);
612  while (nodes[0]->id() == nodes[1]->id() || nodes[0]->id() == nodes[2]->id())
613  nodes[0] = nb_elem->node_ptr(++n_int);
614 
615  // The penalty value. \f$ \frac{1}{\epsilon} \f$
616  const Real penalty = 1.e10;
617 
618  // With this simple method, clamped boundary conditions are
619  // obtained by penalizing the displacements of all four nodes.
620  // This ensures that the displacement field vanishes on the
621  // boundary side s.
622  for (unsigned int n=0; n<4; ++n)
623  {
624  const dof_id_type u_dof = nodes[n]->dof_number (system.number(), u_var, 0);
625  const dof_id_type v_dof = nodes[n]->dof_number (system.number(), v_var, 0);
626  const dof_id_type w_dof = nodes[n]->dof_number (system.number(), w_var, 0);
627  system.matrix->add (u_dof, u_dof, penalty);
628  system.matrix->add (v_dof, v_dof, penalty);
629  system.matrix->add (w_dof, w_dof, penalty);
630  }
631  }
632  } // end of ghost element loop
633 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
A Node is like a Point, but with more information.
Definition: node.h:52
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
static const unsigned int prev[3]
A lookup table for the decrement modulo 3 operation, for iterating through the three nodes per elemen...
Real norm() const
Definition: type_vector.h:909
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
const T & get(const std::string &) const
Definition: parameters.h:430
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
TypeTensor< T > transpose() const
Definition: type_tensor.h:992
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const DofMap & get_dof_map() const
Definition: system.h:2030
Defines a dense submatrix for use in Finite Element-type computations.
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:879
static const unsigned int next[3]
A lookup table for the increment modulo 3 operation, for iterating through the three nodes per elemen...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int number() const
Definition: system.h:2006
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
SparseMatrix< Number > * matrix
The system matrix.
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
int main ( int  argc,
char **  argv 
)

Definition at line 77 of file miscellaneous_ex11.C.

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_shell(), libMesh::System::attach_assemble_function(), libMesh::LibMeshInit::comm(), libMesh::ParallelObject::comm(), libMesh::System::current_solution(), libMesh::default_solver_package(), libMesh::DofObject::dof_number(), libMesh::DofMap::end_dof(), libMesh::DofMap::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, libMesh::Parallel::Communicator::sum(), libMesh::MeshRefinement::uniformly_refine(), libMesh::System::variable_number(), libMesh::ExodusII_IO::write(), libMesh::VTKIO::write(), and libMesh::MeshOutput< MT >::write_equation_systems().

78 {
79  // Initialize libMesh.
80  LibMeshInit init (argc, argv);
81 
82  // This example requires a linear solver package.
83  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
84  "--enable-petsc, --enable-trilinos, or --enable-eigen");
85 
86  // Skip this 3D example if libMesh was compiled as 1D/2D-only.
87  libmesh_example_requires (3 == LIBMESH_DIM, "3D support");
88 
89  // Skip this example without --enable-node-valence
90 #ifndef LIBMESH_ENABLE_NODE_VALENCE
91  libmesh_example_requires (false, "--enable-node-valence");
92 #endif
93 
94  // Skip this example without --enable-amr; requires MeshRefinement
95 #ifndef LIBMESH_ENABLE_AMR
96  libmesh_example_requires(false, "--enable-amr");
97 #else
98 
99  // Skip this example without --enable-second; requires d2phi
100 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
101  libmesh_example_requires(false, "--enable-second");
102 #else
103 
104  // Create a 2D mesh distributed across the default MPI communicator.
105  // Subdivision surfaces do not appear to work with DistributedMesh yet.
106  ReplicatedMesh mesh (init.comm(), 2);
107 
108  // Read the coarse square mesh.
109  mesh.read ("square_mesh.off");
110 
111  // Resize the square plate to edge length L.
112  const Real L = 100.;
114 
115  // Quadrisect the mesh triangles a few times to obtain a
116  // finer mesh. Subdivision surface elements require the
117  // refinement data to be removed afterward.
118  MeshRefinement mesh_refinement (mesh);
119  mesh_refinement.uniformly_refine (3);
121 
122  // Write the mesh before the ghost elements are added.
123 #if defined(LIBMESH_HAVE_VTK)
124  VTKIO(mesh).write ("without_ghosts.pvtu");
125 #endif
126 #if defined(LIBMESH_HAVE_EXODUS_API)
127  ExodusII_IO(mesh).write ("without_ghosts.e");
128 #endif
129 
130  // Print information about the triangulated mesh to the screen.
131  mesh.print_info();
132 
133  // Turn the triangulated mesh into a subdivision mesh
134  // and add an additional row of "ghost" elements around
135  // it in order to complete the extended local support of
136  // the triangles at the boundaries. If the second
137  // argument is set to true, the outermost existing
138  // elements are converted into ghost elements, and the
139  // actual physical mesh is thus getting smaller.
141 
142  // Print information about the subdivision mesh to the screen.
143  mesh.print_info();
144 
145  // Write the mesh with the ghost elements added.
146  // Compare this to the original mesh to see the difference.
147 #if defined(LIBMESH_HAVE_VTK)
148  VTKIO(mesh).write ("with_ghosts.pvtu");
149 #endif
150 #if defined(LIBMESH_HAVE_EXODUS_API)
151  ExodusII_IO(mesh).write ("with_ghosts.e");
152 #endif
153 
154  // Create an equation systems object.
155  EquationSystems equation_systems (mesh);
156 
157  // Declare the system and its variables.
158  // Create a linear implicit system named "Shell".
159  LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem> ("Shell");
160 
161  // Add the three translational deformation variables
162  // "u", "v", "w" to "Shell". Since subdivision shell
163  // elements meet the C1-continuity requirement, no
164  // rotational or other auxiliary variables are needed.
165  // Loop Subdivision Elements are always interpolated
166  // by quartic box splines, hence the order must always
167  // be FOURTH.
168  system.add_variable ("u", FOURTH, SUBDIVISION);
169  system.add_variable ("v", FOURTH, SUBDIVISION);
170  system.add_variable ("w", FOURTH, SUBDIVISION);
171 
172  // Give the system a pointer to the matrix and rhs assembly
173  // function.
175 
176  // Use the parameters of the equation systems object to
177  // tell the shell system about the material properties, the
178  // shell thickness, and the external load.
179  const Real h = 1.;
180  const Real E = 1.e7;
181  const Real nu = 0.;
182  const Real q = 1.;
183  equation_systems.parameters.set<Real> ("thickness") = h;
184  equation_systems.parameters.set<Real> ("young's modulus") = E;
185  equation_systems.parameters.set<Real> ("poisson ratio") = nu;
186  equation_systems.parameters.set<Real> ("uniform load") = q;
187 
188  // Initialize the data structures for the equation system.
189  equation_systems.init();
190 
191  // Print information about the system to the screen.
192  equation_systems.print_info();
193 
194  // Solve the linear system.
195  system.solve();
196 
197  // After solving the system, write the solution to a VTK
198  // or ExodusII output file ready for import in, e.g.,
199  // Paraview.
200 #if defined(LIBMESH_HAVE_VTK)
201  VTKIO(mesh).write_equation_systems ("out.pvtu", equation_systems);
202 #endif
203 #if defined(LIBMESH_HAVE_EXODUS_API)
204  ExodusII_IO(mesh).write_equation_systems ("out.e", equation_systems);
205 #endif
206 
207  // Find the center node to measure the maximum deformation of the plate.
208  Node * center_node = 0;
209  Real nearest_dist_sq = mesh.point(0).norm_sq();
210  for (unsigned int nid=1; nid<mesh.n_nodes(); ++nid)
211  {
212  const Real dist_sq = mesh.point(nid).norm_sq();
213  if (dist_sq < nearest_dist_sq)
214  {
215  nearest_dist_sq = dist_sq;
216  center_node = mesh.node_ptr(nid);
217  }
218  }
219 
220  // Finally, we evaluate the z-displacement "w" at the center node.
221  const unsigned int w_var = system.variable_number ("w");
222  dof_id_type w_dof = center_node->dof_number (system.number(), w_var, 0);
223  Number w = 0;
224  if (w_dof >= system.get_dof_map().first_dof() &&
225  w_dof < system.get_dof_map().end_dof())
226  w = system.current_solution(w_dof);
227  system.comm().sum(w);
228 
229 
230  // The analytic solution for the maximum displacement of
231  // a clamped square plate in pure bending, from Taylor,
232  // Govindjee, Commun. Numer. Meth. Eng. 20, 757-765, 2004.
233  const Real D = E * h*h*h / (12*(1-nu*nu));
234  const Real w_analytic = 0.001265319 * L*L*L*L * q / D;
235 
236  // Print the finite element solution and the analytic
237  // prediction of the maximum displacement of the clamped
238  // square plate to the screen.
239  libMesh::out << "z-displacement of the center point: " << w << std::endl;
240  libMesh::out << "Analytic solution for pure bending: " << w_analytic << std::endl;
241 
242 #endif // #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
243 
244 #endif // #ifdef LIBMESH_ENABLE_AMR
245 
246  // All done.
247  return 0;
248 }
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
A Node is like a Point, but with more information.
Definition: node.h:52
virtual const Point & point(const dof_id_type i) const =0
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
This class provides a specific system class.
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
MeshBase & mesh
void assemble_shell(EquationSystems &es, const std::string &system_name)
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
virtual const Node * node_ptr(const dof_id_type i) const =0
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
This class implements reading and writing meshes in the VTK format.
Definition: vtk_io.h:59
SolverPackage default_solver_package()
Definition: libmesh.C:995
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:1795
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
Real norm_sq() const
Definition: type_vector.h:940
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
virtual void write(const std::string &) libmesh_override
Output the mesh without solutions to a .pvtu file.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
unsigned int number() const
Definition: system.h:2006
virtual void write(const std::string &fname) libmesh_override
This method implements writing a mesh to a specified file.
Definition: exodusII_io.C:789
const Parallel::Communicator & comm() const
OStreamProxy out
void prepare_subdivision_mesh(MeshBase &mesh, bool ghosted=false)
Prepares the mesh for use with subdivision elements.
virtual dof_id_type n_nodes() const =0
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
void flatten(MeshBase &mesh)
Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elem...
uint8_t dof_id_type
Definition: id_types.h:64