330 #if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES) 333 libmesh_assert_equal_to (system_name,
"Shell");
347 const bool distributed_load = es.
parameters.
get<
bool> (
"distributed load");
354 0., 0., 0.5 * (1-nu);
355 Hm *= h * E/(1-nu*nu);
362 0., 0., 0.5 * (1-nu);
363 Hf *= h*h*h/12 * E/(1-nu*nu);
367 Hc0 *= h * 5./6*E/(2*(1+nu));
370 Hc1 *= h*h*h/12 * 5./6*E/(2*(1+nu));
376 std::unique_ptr<FEBase> fe (FEBase::build(
dim, fe_type));
378 fe->attach_quadrature_rule (&qrule);
381 const std::vector<Real> & JxW = fe->get_JxW();
385 const std::vector<RealGradient> & dxyzdxi = fe->get_dxyzdxi();
386 const std::vector<RealGradient> & dxyzdeta = fe->get_dxyzdeta();
387 const std::vector<RealGradient> & d2xyzdxi2 = fe->get_d2xyzdxi2();
388 const std::vector<RealGradient> & d2xyzdeta2 = fe->get_d2xyzdeta2();
389 const std::vector<RealGradient> & d2xyzdxideta = fe->get_d2xyzdxideta();
390 const std::vector<std::vector<Real>> & dphidxi = fe->get_dphidxi();
391 const std::vector<std::vector<Real>> & dphideta = fe->get_dphideta();
392 const std::vector<std::vector<Real>> & phi = fe->get_phi();
424 std::vector<dof_id_type> dof_indices;
425 std::vector<std::vector<dof_id_type>> dof_indices_var(6);
429 for (
const auto & elem :
mesh.active_local_element_ptr_range())
432 for (
unsigned int var=0; var<6; var++)
433 dof_map.
dof_indices (elem, dof_indices_var[var], var);
435 const unsigned int n_dofs = dof_indices.size();
436 const unsigned int n_var_dofs = dof_indices_var[0].size();
439 std::vector<Point> nodes;
440 for (
unsigned int i=0; i<elem->n_nodes(); ++i)
441 nodes.push_back(elem->reference_elem()->node_ref(i));
442 fe->reinit (elem, &nodes);
445 std::vector<MyMatrix3d> Qnode;
446 for (
unsigned int i=0; i<elem->n_nodes(); ++i)
449 a1 << dxyzdxi[i](0), dxyzdxi[i](1), dxyzdxi[i](2);
451 a2 << dxyzdeta[i](0), dxyzdeta[i](1), dxyzdeta[i](2);
472 C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
473 -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
479 Ke.
resize (n_dofs, n_dofs);
480 for (
unsigned int var_i=0; var_i<6; var_i++)
481 for (
unsigned int var_j=0; var_j<6; var_j++)
482 Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
485 Fe_w.reposition(2*n_var_dofs,n_var_dofs);
491 for (
unsigned int qp=0; qp<qrule.n_points(); ++qp)
496 a1 << dxyzdxi[qp](0), dxyzdxi[qp](1), dxyzdxi[qp](2);
498 a2 << dxyzdeta[qp](0), dxyzdeta[qp](1), dxyzdeta[qp](2);
510 F0it =
F0.inverse().transpose();
527 C+1./(1+C)*ny*ny, -1./(1+C)*nx*ny, nx,
528 -1./(1+C)*nx*ny, C+1./(1+C)*nx*nx, ny,
533 C0 = F0it.block<3,2>(0,0).transpose()*Q.block<3,2>(0,0);
536 MyVector3d d2Xdxi2(d2xyzdxi2[qp](0), d2xyzdxi2[qp](1), d2xyzdxi2[qp](2));
537 MyVector3d d2Xdeta2(d2xyzdeta2[qp](0), d2xyzdeta2[qp](1), d2xyzdeta2[qp](2));
538 MyVector3d d2Xdxideta(d2xyzdxideta[qp](0), d2xyzdxideta[qp](1), d2xyzdxideta[qp](2));
542 n.dot(d2Xdxi2), n.dot(d2Xdxideta),
543 n.dot(d2Xdxideta), n.dot(d2Xdeta2);
545 MyVector3d dndxi = -b(0,0)*F0it.col(0) - b(0,1)*F0it.col(1);
546 MyVector3d dndeta = -b(1,0)*F0it.col(0) - b(1,1)*F0it.col(1);
550 F0it.col(1).dot(dndeta), -F0it.col(0).dot(dndeta),
551 -F0it.col(1).dot(dndxi), F0it.col(0).dot(dndxi);
557 Real H = 0.5*(dndxi.dot(F0it.col(0))+dndeta.dot(F0it.col(1)));
560 for (
unsigned int i=0; i<n_var_dofs; ++i)
563 Real C1i = dphidxi[i][qp]*C0(0,0) + dphideta[i][qp]*C0(1,0);
564 Real C2i = dphidxi[i][qp]*C0(0,1) + dphideta[i][qp]*C0(1,1);
567 B0I = MyMatrixXd::Zero(3, 5);
568 B0I.block<1,3>(0,0) = C1i*Q.col(0).transpose();
569 B0I.block<1,3>(1,0) = C2i*Q.col(1).transpose();
570 B0I.block<1,3>(2,0) = C2i*Q.col(0).transpose()+C1i*Q.col(1).transpose();
573 Real bc1i = dphidxi[i][qp]*bc(0,0) + dphideta[i][qp]*bc(1,0);
574 Real bc2i = dphidxi[i][qp]*bc(0,1) + dphideta[i][qp]*bc(1,1);
576 MyVector2d V1i(-Q.col(0).dot(Qnode[i].col(1)),
577 Q.col(0).dot(Qnode[i].col(0)));
579 MyVector2d V2i(-Q.col(1).dot(Qnode[i].col(1)),
580 Q.col(1).dot(Qnode[i].col(0)));
583 B1I = MyMatrixXd::Zero(3,5);
584 B1I.block<1,3>(0,0) = bc1i*Q.col(0).transpose();
585 B1I.block<1,3>(1,0) = bc2i*Q.col(1).transpose();
586 B1I.block<1,3>(2,0) = bc2i*Q.col(0).transpose()+bc1i*Q.col(1).transpose();
588 B1I.block<1,2>(0,3) = C1i*V1i.transpose();
589 B1I.block<1,2>(1,3) = C2i*V2i.transpose();
590 B1I.block<1,2>(2,3) = C2i*V1i.transpose()+C1i*V2i.transpose();
594 B2I = MyMatrixXd::Zero(3,5);
596 B2I.block<1,2>(0,3) = bc1i*V1i.transpose();
597 B2I.block<1,2>(1,3) = bc2i*V2i.transpose();
598 B2I.block<1,2>(2,3) = bc2i*V1i.transpose()+bc1i*V2i.transpose();
602 Bc0I = MyMatrixXd::Zero(2,5);
603 Bc0I.block<1,3>(0,0) = C1i*Q.col(2).transpose();
604 Bc0I.block<1,3>(1,0) = C2i*Q.col(2).transpose();
605 Bc0I.block<1,2>(0,3) = phi[i][qp]*V1i.transpose();
606 Bc0I.block<1,2>(1,3) = phi[i][qp]*V2i.transpose();
610 Bc1I = MyMatrixXd::Zero(2,5);
611 Bc1I.block<1,3>(0,0) = bc1i*Q.col(2).transpose();
612 Bc1I.block<1,3>(1,0) = bc2i*Q.col(2).transpose();
615 MyVector2d BdxiI(dphidxi[i][qp],dphideta[i][qp]);
618 for (
unsigned int j=0; j<n_var_dofs; ++j)
622 Real C1j = dphidxi[j][qp]*C0(0,0) + dphideta[j][qp]*C0(1,0);
623 Real C2j = dphidxi[j][qp]*C0(0,1) + dphideta[j][qp]*C0(1,1);
626 B0J = MyMatrixXd::Zero(3,5);
627 B0J.block<1,3>(0,0) = C1j*Q.col(0).transpose();
628 B0J.block<1,3>(1,0) = C2j*Q.col(1).transpose();
629 B0J.block<1,3>(2,0) = C2j*Q.col(0).transpose()+C1j*Q.col(1).transpose();
632 Real bc1j = dphidxi[j][qp]*bc(0,0) + dphideta[j][qp]*bc(1,0);
633 Real bc2j = dphidxi[j][qp]*bc(0,1) + dphideta[j][qp]*bc(1,1);
635 MyVector2d V1j(-Q.col(0).dot(Qnode[j].col(1)),
636 Q.col(0).dot(Qnode[j].col(0)));
638 MyVector2d V2j(-Q.col(1).dot(Qnode[j].col(1)),
639 Q.col(1).dot(Qnode[j].col(0)));
642 B1J = MyMatrixXd::Zero(3,5);
643 B1J.block<1,3>(0,0) = bc1j*Q.col(0).transpose();
644 B1J.block<1,3>(1,0) = bc2j*Q.col(1).transpose();
645 B1J.block<1,3>(2,0) = bc2j*Q.col(0).transpose()+bc1j*Q.col(1).transpose();
647 B1J.block<1,2>(0,3) = C1j*V1j.transpose();
648 B1J.block<1,2>(1,3) = C2j*V2j.transpose();
649 B1J.block<1,2>(2,3) = C2j*V1j.transpose()+C1j*V2j.transpose();
653 B2J = MyMatrixXd::Zero(3,5);
655 B2J.block<1,2>(0,3) = bc1j*V1j.transpose();
656 B2J.block<1,2>(1,3) = bc2j*V2j.transpose();
657 B2J.block<1,2>(2,3) = bc2j*V1j.transpose()+bc1j*V2j.transpose();
661 Bc0J = MyMatrixXd::Zero(2,5);
662 Bc0J.block<1,3>(0,0) = C1j*Q.col(2).transpose();
663 Bc0J.block<1,3>(1,0) = C2j*Q.col(2).transpose();
664 Bc0J.block<1,2>(0,3) = phi[j][qp]*V1j.transpose();
665 Bc0J.block<1,2>(1,3) = phi[j][qp]*V2j.transpose();
669 Bc1J = MyMatrixXd::Zero(2,5);
670 Bc1J.block<1,3>(0,0) = bc1j*Q.col(2).transpose();
671 Bc1J.block<1,3>(1,0) = bc2j*Q.col(2).transpose();
674 MyVector2d BdxiJ(dphidxi[j][qp], dphideta[j][qp]);
680 local_KIJ = JxW[qp] * (
681 B0I.transpose() * Hm * B0J
682 + B2I.transpose() * Hf * B0J
683 + B0I.transpose() * Hf * B2J
684 + B1I.transpose() * Hf * B1J
685 + 2*H * B0I.transpose() * Hf * B1J
686 + 2*H * B1I.transpose() * Hf * B0J
687 + Bc0I.transpose() * Hc0 * Bc0J
688 + Bc1I.transpose() * Hc1 * Bc1J
689 + 2*H * Bc0I.transpose() * Hc1 * Bc1J
690 + 2*H * Bc1I.transpose() * Hc1 * Bc0J
695 full_local_KIJ = MyMatrixXd::Zero(6, 6);
696 full_local_KIJ.block<5,5>(0,0)=local_KIJ;
702 full_local_KIJ(5,5) =
Real(Hf(0,0)*JxW[qp]*BdI.transpose()*BdJ);
707 TI = MyMatrixXd::Identity(6,6);
708 TI.block<3,3>(3,3) = Qnode[i].transpose();
710 TJ = MyMatrixXd::Identity(6,6);
711 TJ.block<3,3>(3,3) = Qnode[j].transpose();
712 global_KIJ = TI.transpose()*full_local_KIJ*TJ;
717 for (
unsigned int k=0;k<6;k++)
718 for (
unsigned int l=0;l<6;l++)
719 Ke_var[k][l](i,j) += global_KIJ(k,l);
725 if (distributed_load)
728 for (
unsigned int shellface=0; shellface<2; shellface++)
730 std::vector<boundary_id_type> bids;
733 for (std::size_t k=0; k<bids.size(); k++)
735 for (
unsigned int qp=0; qp<qrule.n_points(); ++qp)
736 for (
unsigned int i=0; i<n_var_dofs; ++i)
737 Fe_w(i) -= JxW[qp] * phi[i][qp];
751 if (!distributed_load)
761 for (
const auto & node :
mesh.node_ptr_range())
762 if (((*node) - C).norm() < 1e-3)
763 system.
rhs->
set(node->dof_number(0, 2, 0), -q/4);
769 #endif // defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_ENABLE_SECOND_DERIVATIVES) class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
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.
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.
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
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]...
Eigen::Matrix< libMesh::Real, Eigen::Dynamic, Eigen::Dynamic > MyMatrixXd
NumericVector< Number > * rhs
The system matrix.
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
const T_sys & get_system(std::string_view name) const
Eigen::Matrix< libMesh::Real, 3, 3 > MyMatrix3d
This is the MeshBase class.
Defines a dense subvector for use in finite element computations.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
This class handles the numbering of degrees of freedom on a mesh.
void libmesh_ignore(const Args &...)
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 T & get(std::string_view) const
Defines a dense submatrix for use in Finite Element-type computations.
Eigen::Matrix< libMesh::Real, 2, 2 > MyMatrix2d
virtual void close()=0
Calls the NumericVector's internal assembly routines, ensuring that the values are consistent across ...
const FEType & variable_type(const unsigned int i) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Parameters parameters
Data structure holding arbitrary parameters.
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
Eigen::Matrix< libMesh::Real, 2, 1 > MyVector2d
const DofMap & get_dof_map() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
const SparseMatrix< Number > & get_system_matrix() const
Eigen::Matrix< libMesh::Real, 3, 1 > MyVector3d