36 #include "libmesh/libmesh.h" 37 #include "libmesh/mesh.h" 38 #include "libmesh/linear_implicit_system.h" 39 #include "libmesh/equation_systems.h" 40 #include "libmesh/fe.h" 41 #include "libmesh/quadrature.h" 42 #include "libmesh/node.h" 43 #include "libmesh/elem.h" 44 #include "libmesh/dof_map.h" 45 #include "libmesh/quadrature_gauss.h" 46 #include "libmesh/vector_value.h" 47 #include "libmesh/tensor_value.h" 48 #include "libmesh/dense_matrix.h" 49 #include "libmesh/dense_submatrix.h" 50 #include "libmesh/dense_vector.h" 51 #include "libmesh/dense_subvector.h" 52 #include "libmesh/sparse_matrix.h" 53 #include "libmesh/numeric_vector.h" 54 #include "libmesh/exodusII_io.h" 55 #include "libmesh/dirichlet_boundaries.h" 56 #include "libmesh/zero_function.h" 57 #include "libmesh/linear_solver.h" 58 #include "libmesh/getpot.h" 59 #include "libmesh/enum_solver_package.h" 60 #include "libmesh/enum_solver_type.h" 61 #include "libmesh/parallel.h" 64 #ifdef LIBMESH_HAVE_EIGEN 65 #include "libmesh/ignore_warnings.h" 66 # include <Eigen/Dense> 67 #include "libmesh/restore_warnings.h" 74 typedef Eigen::Matrix<libMesh::Real, Eigen::Dynamic, Eigen::Dynamic>
MyMatrixXd;
84 const std::string & system_name);
87 int main (
int argc,
char ** argv)
93 libmesh_example_requires (3 == LIBMESH_DIM,
"3D support");
96 #ifndef LIBMESH_ENABLE_DIRICHLET 97 libmesh_example_requires(
false,
"--enable-dirichlet");
101 #ifndef LIBMESH_HAVE_EXODUS_API 102 libmesh_example_requires (
false,
"ExodusII support");
105 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES 106 libmesh_example_requires (
false,
"second derivatives enabled");
111 #ifndef LIBMESH_HAVE_EIGEN 112 libmesh_example_requires(
false,
"--enable-eigen");
117 #ifndef LIBMESH_HAVE_XDR 118 libmesh_example_requires (
false,
"XDR support");
123 #ifdef LIBMESH_DEFAULT_QUADRUPLE_PRECISION 124 libmesh_example_requires (
false,
"--disable-quadruple-precision");
178 equation_systems.
parameters.
set<
bool>(
"distributed load") = (distributed_load != 0);
188 #ifdef LIBMESH_ENABLE_DIRICHLET 227 #endif // LIBMESH_ENABLE_DIRICHLET 230 equation_systems.
init();
244 if (distributed_load==0)
247 Node * node_C =
nullptr;
248 Point point_C(0, 3, 3);
250 Real nearest_dist_sq = std::numeric_limits<Real>::max();
255 for (
const auto & node :
mesh.local_node_ptr_range())
258 if (dist_sq < nearest_dist_sq)
260 nearest_dist_sq = dist_sq;
266 unsigned int minrank = 0;
267 system.
comm().
minloc(nearest_dist_sq, minrank);
273 nearest_node_id = node_C->
id();
294 Number w_C_bar = -E*h*w/q;
295 const Real w_C_bar_analytic = 164.24;
299 libMesh::out <<
"z-displacement of the point C: " << w_C_bar << std::endl;
300 libMesh::out <<
"Analytic solution: " << w_C_bar_analytic << std::endl;
304 Point point_D(0, 0, 3);
309 const Real v_D_bar_analytic = 4.114;
313 libMesh::out <<
"y-displacement of the point D: " << v_D_bar << std::endl;
314 libMesh::out <<
"Analytic solution: " << v_D_bar_analytic << std::endl;
326 const std::string & system_name)
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));
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);
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...
int main(int argc, char **argv)
T command_line_next(std::string name, T default_value)
Use GetPot's search()/next() functions to get following arguments from the command line...
This is the EquationSystems class.
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
A Node is like a Point, but with more information.
ConstFunction that simply returns 0.
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.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
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.
void minloc(T &r, unsigned int &min_id) const
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
Eigen::Matrix< libMesh::Real, 3, 3 > MyMatrix3d
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
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Number point_value(unsigned int var, const Point &p, const bool insist_on_success=true, const NumericVector< Number > *sol=nullptr) const
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
const Parallel::Communicator & comm() const
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
virtual void solve() override
Assembles & solves the linear system A*x=b.
const T_sys & get_system(std::string_view name) const
void assemble_shell(EquationSystems &es, const std::string &system_name)
Eigen::Matrix< libMesh::Real, 3, 3 > MyMatrix3d
This is the MeshBase class.
Eigen::Matrix< libMesh::Real, 2, 1 > MyVector2d
Number current_solution(const dof_id_type global_dof_number) const
unsigned int variable_number(std::string_view var) const
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 reposition(const unsigned int ioff, const unsigned int n)
Changes the location of the subvector in the parent vector.
void libmesh_ignore(const Args &...)
unsigned int number() const
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
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Defines a dense submatrix for use in Finite Element-type computations.
virtual const Node * query_node_ptr(const dof_id_type i) const =0
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
unsigned int n_points() const
virtual void write(const std::string &name)=0
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.
Eigen::Matrix< libMesh::Real, 2, 2 > MyMatrix2d
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
virtual void close()=0
Calls the NumericVector's internal assembly routines, ensuring that the values are consistent across ...
Eigen::Matrix< libMesh::Real, Eigen::Dynamic, Eigen::Dynamic > MyMatrixXd
Eigen::Matrix< libMesh::Real, 3, 1 > MyVector3d
const FEType & variable_type(const unsigned int i) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Eigen::Matrix< libMesh::Real, 2, 2 > MyMatrix2d
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.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
virtual void init()
Initialize all the systems.
dof_id_type first_dof(const processor_id_type proc) const
dof_id_type end_dof(const processor_id_type proc) const
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
processor_id_type processor_id() const
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