30 #include "libmesh/libmesh_config.h" 38 #include "libmesh/libmesh.h" 39 #include "libmesh/mesh.h" 40 #include "libmesh/mesh_generation.h" 41 #include "libmesh/linear_implicit_system.h" 42 #include "libmesh/equation_systems.h" 45 #include "libmesh/fe.h" 48 #include "libmesh/quadrature_gauss.h" 52 #include "libmesh/sparse_matrix.h" 53 #include "libmesh/numeric_vector.h" 54 #include "libmesh/dense_matrix.h" 55 #include "libmesh/dense_vector.h" 56 #include "libmesh/elem.h" 57 #include "libmesh/enum_solver_package.h" 58 #include "libmesh/static_condensation.h" 62 #include "libmesh/dof_map.h" 65 #include "libmesh/mesh_refinement.h" 66 #include "libmesh/error_vector.h" 67 #include "libmesh/kelly_error_estimator.h" 70 #include "libmesh/getpot.h" 71 #include "libmesh/exodusII_io.h" 74 #include "libmesh/petsc_linear_solver.h" 79 #if defined(LIBMESH_HAVE_EIGEN_DENSE) && defined(LIBMESH_HAVE_PETSC) 93 main(
int argc,
char ** argv)
98 #if !defined(LIBMESH_HAVE_EIGEN_DENSE) 99 libmesh_example_requires(
false,
"--enable-eigen");
100 #elif !defined(LIBMESH_HAVE_PETSC) 101 libmesh_example_requires(
false,
"--enable-petsc");
108 for (
int i = 1; i < argc; i++)
113 libmesh_example_requires(2 <= LIBMESH_DIM,
"2D support");
116 GetPot infile(
"miscellaneous_ex16.in");
119 infile.parse_command_line(argc, argv);
121 const unsigned int grid_size = infile(
"grid_size", 5);
159 sc_sys.create_static_condensation();
161 #ifdef LIBMESH_ENABLE_AMR 177 equation_systems.
init();
179 #ifdef LIBMESH_ENABLE_AMR 181 const unsigned int max_r_steps = 2;
183 for (
const auto r_step :
make_range(max_r_steps + 1))
193 KSP sc_ksp = sc_solver->ksp();
194 LibmeshPetscCall2(sc_solver->comm(), KSPSetType(sc_ksp, KSPPREONLY));
195 LibmeshPetscCall2(sc_solver->comm(), KSPSetInitialGuessNonzero(sc_ksp, PETSC_FALSE));
199 "mismatching solution");
200 libMesh::out <<
"Static condensation reduced problem size to " 201 << sc_sys.get_static_condensation().get_condensed_mat().m() << std::endl << std::endl;
203 #if defined(LIBMESH_HAVE_EXODUS_API) && !defined(LIBMESH_ENABLE_PARMESH) 207 const std::string file_name =
208 #ifdef LIBMESH_ENABLE_AMR 209 "out_" + std::to_string(r_step) +
".e";
216 #ifdef LIBMESH_ENABLE_AMR 220 if (r_step != max_r_steps)
232 <<
"\nmaximum = " << error.
maximum() << std::endl << std::endl;
243 equation_systems.
reinit();
246 #endif // LIBMESH_ENABLE_AMR 247 #endif // defined(LIBMESH_HAVE_EIGEN_DENSE) && defined(LIBMESH_HAVE_PETSC) 253 #if defined(LIBMESH_HAVE_EIGEN_DENSE) && defined(LIBMESH_HAVE_PETSC) 299 fe->attach_quadrature_rule(&qrule);
312 fe_face->attach_quadrature_rule(&qface);
318 const std::vector<Real> & JxW = fe->get_JxW();
323 const std::vector<Point> & q_point = fe->get_xyz();
326 const std::vector<std::vector<Real>> & phi = fe->get_phi();
330 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
344 std::vector<dof_id_type> dof_indices;
362 for (
const auto & elem :
mesh.active_local_element_ptr_range())
375 const unsigned int n_dofs = cast_int<unsigned int>(dof_indices.size());
385 libmesh_assert_equal_to(n_dofs, phi.size());
396 Ke.
resize(n_dofs, n_dofs);
402 for (
unsigned int qp = 0; qp < qrule.
n_points(); qp++)
408 for (
unsigned int i = 0; i != n_dofs; i++)
409 for (
unsigned int j = 0; j != n_dofs; j++)
411 Ke(i, j) += JxW[qp] * (dphi[i][qp] * dphi[j][qp]);
419 const Real x = q_point[qp](0);
420 const Real y = q_point[qp](1);
421 const Real eps = 1.e-3;
442 for (
unsigned int i = 0; i != n_dofs; i++)
443 Fe(i) += JxW[qp] * fxy * phi[i][qp];
483 for (
auto side : elem->side_index_range())
484 if (elem->neighbor_ptr(side) ==
nullptr)
488 const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
492 const std::vector<Real> & JxW_face = fe_face->get_JxW();
497 const std::vector<Point> & qface_point = fe_face->get_xyz();
501 fe_face->reinit(elem, side);
506 libmesh_assert_equal_to(n_dofs, phi_face.size());
509 for (
unsigned int qp = 0; qp < qface.
n_points(); qp++)
513 const Real xf = qface_point[qp](0);
514 const Real yf = qface_point[qp](1);
518 const Real penalty = 1.e10;
524 for (
unsigned int i = 0; i != n_dofs; i++)
525 for (
unsigned int j = 0; j != n_dofs; j++)
526 Ke(i, j) += JxW_face[qp] * penalty * phi_face[i][qp] * phi_face[j][qp];
530 for (
unsigned int i = 0; i != n_dofs; i++)
531 Fe(i) += JxW_face[qp] * penalty *
value * phi_face[i][qp];
544 sc->set_current_elem(*elem);
557 #endif // defined(LIBMESH_HAVE_EIGEN_DENSE) && defined(LIBMESH_HAVE_PETSC) virtual T maximum() const
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
This is the EquationSystems class.
Real exact_solution(const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.
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
virtual Real l2_norm() const
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
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
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
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.
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
bool has_static_condensation() 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 T_sys & get_system(std::string_view name) const
unsigned int & max_h_level()
The max_h_level is the greatest refinement level an element should reach.
This is the MeshBase class.
This class handles the numbering of degrees of freedom on a mesh.
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.
Implements (adaptive) mesh refinement algorithms for a MeshBase.
int main(int argc, char **argv)
virtual void write_equation_systems(const std::string &fname, const EquationSystems &es, const std::set< std::string > *system_names=nullptr) override
Writes out the solution for no specific time or timestep.
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
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.
StaticCondensation & get_static_condensation()
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
void assemble_poisson(EquationSystems &es, const std::string &system_name)
This class implements the Kelly error indicator which is based on the flux jumps between elements...
virtual void close()=0
Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across p...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool refine_and_coarsen_elements()
Refines and coarsens user-requested elements.
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().
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
bool use_unweighted_quadrature_rules
This boolean flag allows you to use "unweighted" quadrature rules (sized to exactly integrate unweigh...
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
void flag_elements_by_error_fraction(const ErrorVector &error_per_cell, const Real refine_fraction=0.3, const Real coarsen_fraction=0.0, const unsigned int max_level=libMesh::invalid_uint)
Flags elements for coarsening and refinement based on the computed error passed in error_per_cell...
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override
This function uses the derived class's jump error estimate formula to estimate the error on each cell...
bool relative_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
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