39 #include "libmesh/libmesh.h" 40 #include "libmesh/replicated_mesh.h" 41 #include "libmesh/mesh_refinement.h" 42 #include "libmesh/gmv_io.h" 43 #include "libmesh/exodusII_io.h" 44 #include "libmesh/equation_systems.h" 45 #include "libmesh/fe.h" 46 #include "libmesh/quadrature_gauss.h" 47 #include "libmesh/dof_map.h" 48 #include "libmesh/sparse_matrix.h" 49 #include "libmesh/numeric_vector.h" 50 #include "libmesh/dense_matrix.h" 51 #include "libmesh/dense_vector.h" 52 #include "libmesh/periodic_boundaries.h" 53 #include "libmesh/periodic_boundary.h" 54 #include "libmesh/mesh_generation.h" 55 #include "libmesh/parsed_function.h" 56 #include "libmesh/getpot.h" 57 #include "libmesh/enum_solver_package.h" 58 #include "libmesh/enum_xdr_mode.h" 59 #include "libmesh/enum_norm_type.h" 63 #include "libmesh/transient_system.h" 64 #include "libmesh/linear_implicit_system.h" 65 #include "libmesh/vector_value.h" 69 #include "libmesh/error_vector.h" 70 #include "libmesh/kelly_error_estimator.h" 73 #include "libmesh/elem.h" 85 #ifdef LIBMESH_ENABLE_AMR 87 const std::string & system_name);
96 const std::string & system_name);
129 int main (
int argc,
char ** argv)
136 "--enable-petsc, --enable-trilinos, or --enable-eigen");
139 libmesh_example_requires(2 <= LIBMESH_DIM,
"2D support");
141 #if !defined(LIBMESH_ENABLE_AMR) 142 libmesh_example_requires(
false,
"--enable-amr");
143 #elif !defined(LIBMESH_HAVE_XDR) 145 libmesh_example_requires(
false,
"--enable-xdr");
146 #elif !defined(LIBMESH_ENABLE_PERIODIC) 147 libmesh_example_requires(
false,
"--enable-periodic");
148 #elif (LIBMESH_DOF_ID_BYTES == 8) 149 libmesh_example_requires(
false,
"--with-dof-id-bytes=4");
162 <<
"\t " << argv[0] <<
" -init_timestep 0\n" 164 <<
"\t " << argv[0] <<
" -read_solution -init_timestep 26\n" 169 for (
int i=1; i<argc; i++)
187 const unsigned int init_timestep =
195 libmesh_error_msg(
"ERROR: Initial timestep not specified!");
200 const unsigned int n_timesteps =
205 libmesh_error_msg(
"ERROR: Number of timesteps not specified");
209 #ifdef LIBMESH_HAVE_FPARSER 212 const bool have_expression =
false;
219 libmesh_example_requires(2 <= LIBMESH_DIM,
"2D support");
234 (
"Convection-Diffusion");
240 DofMap & dof_map = system.get_dof_map();
274 const unsigned int n_refinements =
286 system.add_variable (
"u",
FIRST);
289 system.attach_init_function (
init_cd);
301 equation_systems.read(
"saved_solution.xdr",
DECODE);
306 equation_systems.init ();
308 equation_systems.reinit ();
314 libMesh::out <<
"Initial H1 norm = " << H1norm << std::endl << std::endl;
317 equation_systems.print_info();
319 equation_systems.parameters.set<
unsigned int>
320 (
"linear solver maximum iterations") = 250;
321 equation_systems.parameters.set<
Real>
327 #ifdef LIBMESH_HAVE_GMV 331 #ifdef LIBMESH_HAVE_EXODUS_API 339 #ifdef LIBMESH_HAVE_GMV 343 #ifdef LIBMESH_HAVE_EXODUS_API 359 equation_systems.parameters.set<
Real>(
"diffusivity") = 0.01;
365 const Real dt = 0.025;
366 system.time = init_timestep*dt;
375 for (
unsigned int t_step=init_timestep; t_step<(init_timestep+n_timesteps); t_step++)
381 equation_systems.parameters.set<
Real> (
"time") = system.time;
382 equation_systems.parameters.set<
Real> (
"dt") = dt;
394 << std::setprecision(3)
416 const unsigned int max_r_steps =
420 for (
unsigned int r_step=0; r_step<max_r_steps+1; r_step++)
426 H1norm = system.calculate_norm(*system.solution,
SystemNorm(
H1));
430 if (r_step+1 <= max_r_steps)
479 equation_systems.reinit ();
484 const unsigned int output_freq =
488 if ((t_step+1)%output_freq == 0)
490 #ifdef LIBMESH_HAVE_GMV 500 #ifdef LIBMESH_HAVE_EXODUS_API 514 H1norm = system.calculate_norm(*system.solution,
SystemNorm(
H1));
516 libMesh::out <<
"Final H1 norm = " << H1norm << std::endl << std::endl;
519 equation_systems.write(
"saved_solution.xdr",
ENCODE);
520 #ifdef LIBMESH_HAVE_GMV 524 #ifdef LIBMESH_HAVE_EXODUS_API 529 #endif // #ifndef LIBMESH_ENABLE_AMR 539 const std::string & libmesh_dbg_var(system_name))
543 libmesh_assert_equal_to (system_name,
"Convection-Diffusion");
564 #ifdef LIBMESH_ENABLE_AMR 566 const std::string & libmesh_dbg_var(system_name))
570 libmesh_assert_equal_to (system_name,
"Convection-Diffusion");
584 FEType fe_type = system.variable_type(0);
595 QGauss qrule (
dim, fe_type.default_quadrature_order());
596 QGauss qface (
dim-1, fe_type.default_quadrature_order());
599 fe->attach_quadrature_rule (&qrule);
600 fe_face->attach_quadrature_rule (&qface);
605 const std::vector<Real> & JxW = fe->get_JxW();
608 const std::vector<std::vector<Real>> & phi = fe->get_phi();
612 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
618 const DofMap & dof_map = system.get_dof_map();
630 std::vector<dof_id_type> dof_indices;
637 const Real diffusivity =
650 for (
const auto & elem :
mesh.active_local_element_ptr_range())
664 const unsigned int n_dofs =
665 cast_int<unsigned int>(dof_indices.size());
666 libmesh_assert_equal_to (n_dofs, phi.size());
674 Ke.
resize (n_dofs, n_dofs);
684 for (
unsigned int qp=0; qp<qrule.n_points(); qp++)
691 for (
unsigned int l=0; l != n_dofs; l++)
693 u_old += phi[l][qp]*system.
old_solution (dof_indices[l]);
702 for (
unsigned int i=0; i != n_dofs; i++)
712 (grad_u_old*velocity)*phi[i][qp] +
715 diffusivity*(grad_u_old*dphi[i][qp]))
718 for (
unsigned int j=0; j != n_dofs; j++)
723 phi[i][qp]*phi[j][qp] +
726 (velocity*dphi[j][qp])*phi[i][qp] +
728 diffusivity*(dphi[i][qp]*dphi[j][qp]))
748 matrix.add_matrix (Ke, dof_indices);
749 system.rhs->add_vector (Fe, dof_indices);
754 #endif // #ifdef LIBMESH_ENABLE_AMR 761 std::ostringstream oss;
764 oss << std::setw(3) << std::setfill(
'0') << number;
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
void assemble_cd(EquationSystems &es, const std::string &system_name)
T command_line_next(std::string name, T default_value)
Use GetPot's search()/next() functions to get following arguments from the command line...
Real exact_solution(const Real x, const Real y, const Real t)
This is the exact solution that we are trying to obtain.
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.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
void add_periodic_boundary(const PeriodicBoundaryBase &periodic_boundary)
Adds a copy of the specified periodic boundary to the system.
This class provides the ability to map between arbitrary, user-defined strings and several data types...
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
VectorValue< Real > RealVectorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
static constexpr Real TOLERANCE
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...
int main(int argc, char **argv)
void resize(const unsigned int n)
Resize the vector.
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
boundary_id_type pairedboundary
Manages storage and variables for transient systems.
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
void init_cd(EquationSystems &es, const std::string &system_name)
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
The definition of a periodic boundary.
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.
NumericVector< Number > * old_local_solution
All the values I need to compute my contribution to the simulation at hand.
This is the MeshBase class.
SolverPackage default_solver_package()
This class handles the numbering of degrees of freedom on a mesh.
boundary_id_type myboundary
The boundary ID of this boundary and its counterpart.
const T & get(std::string_view) const
Implements (adaptive) mesh refinement algorithms for a MeshBase.
std::string exodus_filename(unsigned number)
PeriodicBoundaries * get_periodic_boundaries()
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.
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.
std::unique_ptr< FunctionBase< Number > > parsed_solution
Number old_solution(const dof_id_type global_dof_number) const
virtual void write(const std::string &name) const =0
This class implements the Kelly error indicator which is based on the flux jumps between elements...
std::ios_base::fmtflags flags() const
Get the associated format flags.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
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().
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.
Number exact_value(const Point &p, const Parameters ¶meters, const std::string &, const std::string &)
unsigned int mesh_dimension() const
Parameters parameters
Data structure holding arbitrary parameters.
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 on_command_line(std::string arg)
void set_periodic_boundaries_ptr(PeriodicBoundaries *pb_ptr)
Sets the PeriodicBoundaries pointer.
A Point defines a location in LIBMESH_DIM dimensional Real space.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.