libMesh
Functions
adaptivity_ex1.C File Reference

Go to the source code of this file.

Functions

void assemble_1D (EquationSystems &es, const std::string &system_name)
 
int main (int argc, char **argv)
 

Function Documentation

◆ assemble_1D()

void assemble_1D ( EquationSystems es,
const std::string &  system_name 
)

Definition at line 181 of file adaptivity_ex1.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), dim, libMesh::FIFTH, libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), libMesh::libmesh_ignore(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), and libMesh::ExplicitSystem::rhs.

Referenced by main().

183 {
184  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
185  libmesh_ignore(es, system_name);
186 
187 #ifdef LIBMESH_ENABLE_AMR
188 
189  // It is a good idea to check we are solving the correct system
190  libmesh_assert_equal_to (system_name, "1D");
191 
192  // Get a reference to the mesh object
193  const MeshBase & mesh = es.get_mesh();
194 
195  // The dimension we are using, i.e. dim==1
196  const unsigned int dim = mesh.mesh_dimension();
197 
198  // Get a reference to the system we are solving
200 
201  // Get a reference to the DofMap object for this system. The DofMap object
202  // handles the index translation from node and element numbers to degree of
203  // freedom numbers. DofMap's are discussed in more detail in future examples.
204  const DofMap & dof_map = system.get_dof_map();
205 
206  // Get a constant reference to the Finite Element type for the first
207  // (and only) variable in the system.
208  FEType fe_type = dof_map.variable_type(0);
209 
210  // Build a finite element object of the specified type. The build
211  // function dynamically allocates memory so we use a std::unique_ptr in this case.
212  // A std::unique_ptr is a pointer that cleans up after itself. See examples 3 and 4
213  // for more details on std::unique_ptr.
214  std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
215 
216  // Tell the finite element object to use fifth order Gaussian quadrature
217  QGauss qrule(dim, FIFTH);
218  fe->attach_quadrature_rule(&qrule);
219 
220  // Here we define some references to cell-specific data that will be used to
221  // assemble the linear system.
222 
223  // The element Jacobian * quadrature weight at each integration point.
224  const std::vector<Real> & JxW = fe->get_JxW();
225 
226  // The element shape functions evaluated at the quadrature points.
227  const std::vector<std::vector<Real>> & phi = fe->get_phi();
228 
229  // The element shape function gradients evaluated at the quadrature points.
230  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
231 
232  // Declare a dense matrix and dense vector to hold the element matrix
233  // and right-hand-side contribution
236 
237  // This vector will hold the degree of freedom indices for the element.
238  // These define where in the global system the element degrees of freedom
239  // get mapped.
240  std::vector<dof_id_type> dof_indices;
241 
242  // The global system matrix
243  SparseMatrix<Number> & matrix = system.get_system_matrix();
244 
245  // We now loop over all the active elements in the mesh in order to calculate
246  // the matrix and right-hand-side contribution from each element. Use a
247  // const_element_iterator to loop over the elements. We make
248  // el_end const as it is used only for the stopping condition of the loop.
249  for (const auto & elem : mesh.active_local_element_ptr_range())
250  {
251  // Get the degree of freedom indices for the current element.
252  // These define where in the global matrix and right-hand-side this
253  // element will contribute to.
254  dof_map.dof_indices(elem, dof_indices);
255 
256  // Compute the element-specific data for the current element. This
257  // involves computing the location of the quadrature points (q_point)
258  // and the shape functions (phi, dphi) for the current element.
259  fe->reinit(elem);
260 
261  // Store the number of local degrees of freedom contained in this element
262  const unsigned int n_dofs =
263  cast_int<unsigned int>(dof_indices.size());
264  libmesh_assert_equal_to (n_dofs, phi.size());
265 
266  // We resize and zero out Ke and Fe (resize() also clears the matrix and
267  // vector). In this example, all elements in the mesh are EDGE3's, so
268  // Ke will always be 3x3, and Fe will always be 3x1. If the mesh contained
269  // different element types, then the size of Ke and Fe would change.
270  Ke.resize(n_dofs, n_dofs);
271  Fe.resize(n_dofs);
272 
273 
274  // Now loop over quadrature points to handle numerical integration
275  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
276  {
277  // Now build the element matrix and right-hand-side using loops to
278  // integrate the test functions (i) against the trial functions (j).
279  for (unsigned int i=0; i != n_dofs; i++)
280  {
281  Fe(i) += JxW[qp]*phi[i][qp];
282 
283  for (unsigned int j=0; j != n_dofs; j++)
284  {
285  Ke(i,j) += JxW[qp]*(1.e-3*dphi[i][qp]*dphi[j][qp] +
286  phi[i][qp]*phi[j][qp]);
287  }
288  }
289  }
290 
291 
292  // At this point we have completed the matrix and RHS summation. The
293  // final step is to apply boundary conditions, which in this case are
294  // simple Dirichlet conditions with u(0) = u(1) = 0.
295 
296  // Define the penalty parameter used to enforce the BC's
297  double penalty = 1.e10;
298 
299  // Loop over the sides of this element. For a 1D element, the "sides"
300  // are defined as the nodes on each edge of the element, i.e. 1D elements
301  // have 2 sides.
302  for (auto s : elem->side_index_range())
303  {
304  // If this element has a nullptr neighbor, then it is on the edge of the
305  // mesh and we need to enforce a boundary condition using the penalty
306  // method.
307  if (elem->neighbor_ptr(s) == nullptr)
308  {
309  Ke(s,s) += penalty;
310  Fe(s) += 0*penalty;
311  }
312  }
313 
314  // This is a function call that is necessary when using adaptive
315  // mesh refinement. See Adaptivity Example 2 for more details.
316  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
317 
318  // Add Ke and Fe to the global matrix and right-hand-side.
319  matrix.add_matrix(Ke, dof_indices);
320  system.rhs->add_vector(Fe, dof_indices);
321  }
322 #endif // #ifdef LIBMESH_ENABLE_AMR
323 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
unsigned int dim
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
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]...
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Definition: mesh_base.h:74
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
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 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().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:324
const DofMap & get_dof_map() const
Definition: system.h:2293
const SparseMatrix< Number > & get_system_matrix() const

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 60 of file adaptivity_ex1.C.

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_1D(), libMesh::System::attach_assemble_function(), libMesh::MeshTools::Generation::build_line(), libMesh::MeshRefinement::coarsen_fraction(), libMesh::command_line_next(), libMesh::default_solver_package(), libMesh::EDGE3, libMesh::JumpErrorEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::EquationSystems::get_system(), libMesh::GnuPlotIO::GRID_ON, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::StatisticsVector< T >::l2_norm(), libMesh::MeshRefinement::max_h_level(), libMesh::StatisticsVector< T >::maximum(), mesh, libMesh::out, libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::MeshRefinement::refine_fraction(), libMesh::EquationSystems::reinit(), libMesh::SECOND, libMesh::JumpErrorEstimator::use_unweighted_quadrature_rules, and libMesh::MeshOutput< MT >::write_equation_systems().

61 {
62  // Initialize the library. This is necessary because the library
63  // may depend on a number of other libraries (i.e. MPI and PETSc)
64  // that require initialization before use. When the LibMeshInit
65  // object goes out of scope, other libraries and resources are
66  // finalized.
67  LibMeshInit init (argc, argv);
68 
69  // This example requires a linear solver package.
70  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
71  "--enable-petsc, --enable-trilinos, or --enable-eigen");
72 
73  // Skip adaptive examples on a non-adaptive libMesh build
74 #ifndef LIBMESH_ENABLE_AMR
75  libmesh_example_requires(false, "--enable-amr");
76 #else
77 
78  // Create a mesh, with dimension to be overridden later, on the
79  // default MPI communicator.
80  Mesh mesh(init.comm());
81 
82  const int n = libMesh::command_line_next("-n", 4);
83 
84  // Build a 1D mesh with 4 elements from x=0 to x=1, using
85  // EDGE3 (i.e. quadratic) 1D elements. They are called EDGE3 elements
86  // because a quadratic element contains 3 nodes.
88 
89  // Define the equation systems object and the system we are going
90  // to solve. See Introduction Example 2 for more details.
91  EquationSystems equation_systems(mesh);
92  LinearImplicitSystem & system = equation_systems.add_system
93  <LinearImplicitSystem>("1D");
94 
95  // Add a variable "u" to the system, using second-order approximation
96  system.add_variable("u", SECOND);
97 
98  // Give the system a pointer to the matrix assembly function. This
99  // will be called when needed by the library.
101 
102  // Define the mesh refinement object that takes care of adaptively
103  // refining the mesh.
104  MeshRefinement mesh_refinement(mesh);
105 
106  // These parameters determine the proportion of elements that will
107  // be refined and coarsened. Any element within 30% of the maximum
108  // error on any element will be refined, and any element within 30%
109  // of the minimum error on any element might be coarsened
110  mesh_refinement.refine_fraction() = 0.7;
111  mesh_refinement.coarsen_fraction() = 0.3;
112  // We won't refine any element more than 5 times in total
113  mesh_refinement.max_h_level() = 5;
114 
115  // Initialize the data structures for the equation system.
116  equation_systems.init();
117 
118  // Refinement parameters
119  const unsigned int max_r_steps = 5; // Refine the mesh 5 times
120 
121  // Define the refinement loop
122  for (unsigned int r_step=0; r_step<=max_r_steps; r_step++)
123  {
124  // Solve the equation system
125  equation_systems.get_system("1D").solve();
126 
127  // We need to ensure that the mesh is not refined on the last iteration
128  // of this loop, since we do not want to refine the mesh unless we are
129  // going to solve the equation system for that refined mesh.
130  if (r_step != max_r_steps)
131  {
132  // Error estimation objects, see Adaptivity Example 2 for details
133  ErrorVector error;
134  KellyErrorEstimator error_estimator;
135  error_estimator.use_unweighted_quadrature_rules = true;
136 
137  // Compute the error for each active element
138  error_estimator.estimate_error(system, error);
139 
140  // Output error estimate magnitude
141  libMesh::out << "Error estimate\nl2 norm = "
142  << error.l2_norm()
143  << "\nmaximum = "
144  << error.maximum()
145  << std::endl;
146 
147  // Flag elements to be refined and coarsened
148  mesh_refinement.flag_elements_by_error_fraction (error);
149 
150  // Perform refinement and coarsening
151  mesh_refinement.refine_and_coarsen_elements();
152 
153  // Reinitialize the equation_systems object for the newly refined
154  // mesh. One of the steps in this is project the solution onto the
155  // new mesh
156  equation_systems.reinit();
157  }
158  }
159 
160  // Construct gnuplot plotting object, pass in mesh, title of plot
161  // and boolean to indicate use of grid in plot. The grid is used to
162  // show the edges of each element in the mesh.
163  GnuPlotIO plot(mesh, "Adaptivity Example 1", GnuPlotIO::GRID_ON);
164 
165  // Write out script to be called from within gnuplot:
166  // Load gnuplot, then type "call 'gnuplot_script'" from gnuplot prompt
167  plot.write_equation_systems("gnuplot_script", equation_systems);
168 #endif // #ifndef LIBMESH_ENABLE_AMR
169 
170  // All done. libMesh objects are destroyed here. Because the
171  // LibMeshInit object was created first, its destruction occurs
172  // last, and it's destructor finalizes any external libraries and
173  // checks for leaked memory.
174  return 0;
175 }
virtual T maximum() const
Definition: statistics.C:62
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1011
This is the EquationSystems class.
virtual Real l2_norm() const
Definition: statistics.C:37
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
MeshBase & mesh
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
SolverPackage default_solver_package()
Definition: libmesh.C:1050
This class implements writing meshes using GNUplot, designed for use only with 1D meshes...
Definition: gnuplot_io.h:43
Implements (adaptive) mesh refinement algorithms for a MeshBase.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
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.
Definition: system.C:1305
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.
Definition: system.C:2109
This class implements the Kelly error indicator which is based on the flux jumps between elements...
OStreamProxy out
bool use_unweighted_quadrature_rules
This boolean flag allows you to use "unweighted" quadrature rules (sized to exactly integrate unweigh...
void build_line(UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 1D meshes.
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&#39;s jump error estimate formula to estimate the error on each cell...
void assemble_1D(EquationSystems &es, const std::string &system_name)
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50