41 #include "libmesh/libmesh.h" 42 #include "libmesh/mesh.h" 43 #include "libmesh/mesh_generation.h" 44 #include "libmesh/exodusII_io.h" 45 #include "libmesh/linear_implicit_system.h" 46 #include "libmesh/equation_systems.h" 47 #include "libmesh/fe.h" 48 #include "libmesh/quadrature_gauss.h" 49 #include "libmesh/dof_map.h" 50 #include "libmesh/sparse_matrix.h" 51 #include "libmesh/numeric_vector.h" 52 #include "libmesh/dense_matrix.h" 53 #include "libmesh/dense_submatrix.h" 54 #include "libmesh/dense_vector.h" 55 #include "libmesh/dense_subvector.h" 56 #include "libmesh/perf_log.h" 57 #include "libmesh/elem.h" 58 #include "libmesh/boundary_info.h" 59 #include "libmesh/zero_function.h" 60 #include "libmesh/dirichlet_boundaries.h" 61 #include "libmesh/string_to_enum.h" 62 #include "libmesh/getpot.h" 63 #include "libmesh/enum_solver_package.h" 70 const std::string & system_name);
80 int main (
int argc,
char ** argv)
87 "--enable-petsc, --enable-trilinos, or --enable-eigen");
90 const unsigned int dim = 2;
93 libmesh_example_requires(
dim <= LIBMESH_DIM,
"2D support");
96 #ifndef LIBMESH_ENABLE_DIRICHLET 97 libmesh_example_requires(
false,
"--enable-dirichlet");
101 GetPot command_line (argc, argv);
125 #ifdef LIBMESH_ENABLE_DIRICHLET 150 #endif // LIBMESH_ENABLE_DIRICHLET 153 equation_systems.
init();
162 #ifdef LIBMESH_HAVE_EXODUS_API 164 #endif // #ifdef LIBMESH_HAVE_EXODUS_API 172 const std::string & libmesh_dbg_var(system_name))
174 libmesh_assert_equal_to (system_name,
"Elasticity");
190 fe->attach_quadrature_rule (&qrule);
194 fe_face->attach_quadrature_rule (&qface);
196 const std::vector<Real> & JxW = fe->get_JxW();
197 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
211 std::vector<dof_id_type> dof_indices;
212 std::vector<dof_id_type> dof_indices_u;
213 std::vector<dof_id_type> dof_indices_v;
214 std::vector<dof_id_type> dof_indices_lambda;
218 for (
const auto & elem :
mesh.active_local_element_ptr_range())
223 dof_map.
dof_indices (elem, dof_indices_lambda, lambda_var);
225 const unsigned int n_dofs = dof_indices.size();
226 const unsigned int n_u_dofs = dof_indices_u.size();
227 const unsigned int n_v_dofs = dof_indices_v.size();
228 const unsigned int n_lambda_dofs = dof_indices_lambda.size();
232 Ke.
resize (n_dofs, n_dofs);
235 Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
236 Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
238 Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
239 Kvv.
reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
242 Kv_lambda.
reposition (v_var*n_u_dofs, v_var*n_u_dofs+n_v_dofs, n_v_dofs, 1);
243 Klambda_v.reposition (v_var*n_v_dofs+n_v_dofs, v_var*n_v_dofs, 1, n_v_dofs);
245 Fu.reposition (u_var*n_u_dofs, n_u_dofs);
248 for (
unsigned int qp=0; qp<qrule.n_points(); qp++)
250 for (
unsigned int i=0; i<n_u_dofs; i++)
251 for (
unsigned int j=0; j<n_u_dofs; j++)
254 unsigned int C_i, C_j, C_k, C_l;
270 for (
unsigned int i=0; i<n_u_dofs; i++)
271 for (
unsigned int j=0; j<n_v_dofs; j++)
274 unsigned int C_i, C_j, C_k, C_l;
290 for (
unsigned int i=0; i<n_v_dofs; i++)
291 for (
unsigned int j=0; j<n_u_dofs; j++)
294 unsigned int C_i, C_j, C_k, C_l;
310 for (
unsigned int i=0; i<n_v_dofs; i++)
311 for (
unsigned int j=0; j<n_v_dofs; j++)
314 unsigned int C_i, C_j, C_k, C_l;
332 std::vector<boundary_id_type> bc_ids;
333 for (
auto side : elem->side_index_range())
334 if (elem->neighbor_ptr(side) ==
nullptr)
338 const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
339 const std::vector<Real> & JxW_face = fe_face->get_JxW();
341 fe_face->reinit(elem, side);
343 for (std::vector<boundary_id_type>::const_iterator b =
344 bc_ids.begin(); b != bc_ids.end(); ++b)
347 for (
unsigned int qp=0; qp<qface.n_points(); qp++)
351 for (
unsigned int i=0; i<n_v_dofs; i++)
352 Fv(i) += JxW_face[qp] * (-1.) * phi_face[i][qp];
357 for (
unsigned int i=0; i<n_v_dofs; i++)
358 for (
unsigned int j=0; j<n_lambda_dofs; j++)
359 Kv_lambda(i,j) += JxW_face[qp] * (-1.) * phi_face[i][qp];
361 for (
unsigned int i=0; i<n_lambda_dofs; i++)
362 for (
unsigned int j=0; j<n_v_dofs; j++)
363 Klambda_v(i,j) += JxW_face[qp] * (-1.) * phi_face[j][qp];
386 const Real lambda_1 = nu / ((1. + nu) * (1. - 2.*nu));
387 const Real lambda_2 = 0.5 / (1 + nu);
390 Real delta_ij = (i == j) ? 1. : 0.;
391 Real delta_il = (i == l) ? 1. : 0.;
392 Real delta_ik = (i == k) ? 1. : 0.;
393 Real delta_jl = (j == l) ? 1. : 0.;
394 Real delta_jk = (j == k) ? 1. : 0.;
395 Real delta_kl = (k == l) ? 1. : 0.;
397 return lambda_1 * delta_ij * delta_kl + lambda_2 * (delta_ik * delta_jl + delta_il * delta_jk);
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
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.
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.
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 ...
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]...
const FEType & variable_type(const unsigned int c) const
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
Order default_quadrature_order() const
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
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
This is the MeshBase class.
unsigned int variable_number(std::string_view var) const
Defines a dense subvector for use in finite element computations.
SolverPackage default_solver_package()
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.
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.
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.
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.
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.
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
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
Real eval_elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
void reposition(const unsigned int ioff, const unsigned int joff, const unsigned int new_m, const unsigned int new_n)
Changes the location of the submatrix in the parent matrix.
virtual void init()
Initialize all the systems.
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.
const DofMap & get_dof_map() const
const SparseMatrix< Number > & get_system_matrix() const
void assemble_elasticity(EquationSystems &es, const std::string &system_name)
int main(int argc, char **argv)