libMesh
adaptivity_ex1.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Adaptivity Example 1 - Solving 1D PDE Using Adaptive Mesh Refinement</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This example demonstrates how to solve a simple 1D problem
25 // using adaptive mesh refinement. The PDE that is solved is:
26 // -epsilon*u''(x) + u(x) = 1, on the domain [0,1] with boundary conditions
27 // u(0) = u(1) = 0 and where epsilon << 1.
28 //
29 // The approach used to solve 1D problems in libMesh is virtually identical to
30 // solving 2D or 3D problems, so in this sense this example represents a good
31 // starting point for new users. Note that many concepts are used in this
32 // example which are explained more fully in subsequent examples.
33 
34 // Libmesh includes
35 #include "libmesh/mesh.h"
36 #include "libmesh/mesh_generation.h"
37 #include "libmesh/edge_edge3.h"
38 #include "libmesh/gnuplot_io.h"
39 #include "libmesh/equation_systems.h"
40 #include "libmesh/linear_implicit_system.h"
41 #include "libmesh/fe.h"
42 #include "libmesh/getpot.h"
43 #include "libmesh/quadrature_gauss.h"
44 #include "libmesh/sparse_matrix.h"
45 #include "libmesh/dof_map.h"
46 #include "libmesh/numeric_vector.h"
47 #include "libmesh/dense_matrix.h"
48 #include "libmesh/dense_vector.h"
49 #include "libmesh/error_vector.h"
50 #include "libmesh/kelly_error_estimator.h"
51 #include "libmesh/mesh_refinement.h"
52 
53 // Bring in everything from the libMesh namespace
54 using namespace libMesh;
55 
57  const std::string & system_name);
58 
59 int main(int argc, char ** argv)
60 {
61  // Initialize the library. This is necessary because the library
62  // may depend on a number of other libraries (i.e. MPI and PETSc)
63  // that require initialization before use. When the LibMeshInit
64  // object goes out of scope, other libraries and resources are
65  // finalized.
66  LibMeshInit init (argc, argv);
67 
68  // This example requires a linear solver package.
69  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
70  "--enable-petsc, --enable-trilinos, or --enable-eigen");
71 
72  // Skip adaptive examples on a non-adaptive libMesh build
73 #ifndef LIBMESH_ENABLE_AMR
74  libmesh_example_requires(false, "--enable-amr");
75 #else
76 
77  // Create a mesh, with dimension to be overridden later, on the
78  // default MPI communicator.
79  Mesh mesh(init.comm());
80 
81  GetPot command_line (argc, argv);
82 
83  int n = 4;
84  if (command_line.search(1, "-n"))
85  n = command_line.next(n);
86 
87  // Build a 1D mesh with 4 elements from x=0 to x=1, using
88  // EDGE3 (i.e. quadratic) 1D elements. They are called EDGE3 elements
89  // because a quadratic element contains 3 nodes.
91 
92  // Define the equation systems object and the system we are going
93  // to solve. See Introduction Example 2 for more details.
94  EquationSystems equation_systems(mesh);
95  LinearImplicitSystem & system = equation_systems.add_system
96  <LinearImplicitSystem>("1D");
97 
98  // Add a variable "u" to the system, using second-order approximation
99  system.add_variable("u", SECOND);
100 
101  // Give the system a pointer to the matrix assembly function. This
102  // will be called when needed by the library.
104 
105  // Define the mesh refinement object that takes care of adaptively
106  // refining the mesh.
107  MeshRefinement mesh_refinement(mesh);
108 
109  // These parameters determine the proportion of elements that will
110  // be refined and coarsened. Any element within 30% of the maximum
111  // error on any element will be refined, and any element within 30%
112  // of the minimum error on any element might be coarsened
113  mesh_refinement.refine_fraction() = 0.7;
114  mesh_refinement.coarsen_fraction() = 0.3;
115  // We won't refine any element more than 5 times in total
116  mesh_refinement.max_h_level() = 5;
117 
118  // Initialize the data structures for the equation system.
119  equation_systems.init();
120 
121  // Refinement parameters
122  const unsigned int max_r_steps = 5; // Refine the mesh 5 times
123 
124  // Define the refinement loop
125  for (unsigned int r_step=0; r_step<=max_r_steps; r_step++)
126  {
127  // Solve the equation system
128  equation_systems.get_system("1D").solve();
129 
130  // We need to ensure that the mesh is not refined on the last iteration
131  // of this loop, since we do not want to refine the mesh unless we are
132  // going to solve the equation system for that refined mesh.
133  if (r_step != max_r_steps)
134  {
135  // Error estimation objects, see Adaptivity Example 2 for details
136  ErrorVector error;
137  KellyErrorEstimator error_estimator;
138 
139  // Compute the error for each active element
140  error_estimator.estimate_error(system, error);
141 
142  // Output error estimate magnitude
143  libMesh::out << "Error estimate\nl2 norm = "
144  << error.l2_norm()
145  << "\nmaximum = "
146  << error.maximum()
147  << std::endl;
148 
149  // Flag elements to be refined and coarsened
150  mesh_refinement.flag_elements_by_error_fraction (error);
151 
152  // Perform refinement and coarsening
153  mesh_refinement.refine_and_coarsen_elements();
154 
155  // Reinitialize the equation_systems object for the newly refined
156  // mesh. One of the steps in this is project the solution onto the
157  // new mesh
158  equation_systems.reinit();
159  }
160  }
161 
162  // Construct gnuplot plotting object, pass in mesh, title of plot
163  // and boolean to indicate use of grid in plot. The grid is used to
164  // show the edges of each element in the mesh.
165  GnuPlotIO plot(mesh, "Adaptivity Example 1", GnuPlotIO::GRID_ON);
166 
167  // Write out script to be called from within gnuplot:
168  // Load gnuplot, then type "call 'gnuplot_script'" from gnuplot prompt
169  plot.write_equation_systems("gnuplot_script", equation_systems);
170 #endif // #ifndef LIBMESH_ENABLE_AMR
171 
172  // All done. libMesh objects are destroyed here. Because the
173  // LibMeshInit object was created first, its destruction occurs
174  // last, and it's destructor finalizes any external libraries and
175  // checks for leaked memory.
176  return 0;
177 }
178 
179 
180 
181 
182 // Define the matrix assembly function for the 1D PDE we are solving
184  const std::string & system_name)
185 {
186  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
187  libmesh_ignore(es);
188  libmesh_ignore(system_name);
189 
190 #ifdef LIBMESH_ENABLE_AMR
191 
192  // It is a good idea to check we are solving the correct system
193  libmesh_assert_equal_to (system_name, "1D");
194 
195  // Get a reference to the mesh object
196  const MeshBase & mesh = es.get_mesh();
197 
198  // The dimension we are using, i.e. dim==1
199  const unsigned int dim = mesh.mesh_dimension();
200 
201  // Get a reference to the system we are solving
203 
204  // Get a reference to the DofMap object for this system. The DofMap object
205  // handles the index translation from node and element numbers to degree of
206  // freedom numbers. DofMap's are discussed in more detail in future examples.
207  const DofMap & dof_map = system.get_dof_map();
208 
209  // Get a constant reference to the Finite Element type for the first
210  // (and only) variable in the system.
211  FEType fe_type = dof_map.variable_type(0);
212 
213  // Build a finite element object of the specified type. The build
214  // function dynamically allocates memory so we use a UniquePtr in this case.
215  // A UniquePtr is a pointer that cleans up after itself. See examples 3 and 4
216  // for more details on UniquePtr.
217  UniquePtr<FEBase> fe(FEBase::build(dim, fe_type));
218 
219  // Tell the finite element object to use fifth order Gaussian quadrature
220  QGauss qrule(dim, FIFTH);
221  fe->attach_quadrature_rule(&qrule);
222 
223  // Here we define some references to cell-specific data that will be used to
224  // assemble the linear system.
225 
226  // The element Jacobian * quadrature weight at each integration point.
227  const std::vector<Real> & JxW = fe->get_JxW();
228 
229  // The element shape functions evaluated at the quadrature points.
230  const std::vector<std::vector<Real>> & phi = fe->get_phi();
231 
232  // The element shape function gradients evaluated at the quadrature points.
233  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
234 
235  // Declare a dense matrix and dense vector to hold the element matrix
236  // and right-hand-side contribution
239 
240  // This vector will hold the degree of freedom indices for the element.
241  // These define where in the global system the element degrees of freedom
242  // get mapped.
243  std::vector<dof_id_type> dof_indices;
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 int n_dofs = dof_indices.size();
263 
264  // We resize and zero out Ke and Fe (resize() also clears the matrix and
265  // vector). In this example, all elements in the mesh are EDGE3's, so
266  // Ke will always be 3x3, and Fe will always be 3x1. If the mesh contained
267  // different element types, then the size of Ke and Fe would change.
268  Ke.resize(n_dofs, n_dofs);
269  Fe.resize(n_dofs);
270 
271 
272  // Now loop over quadrature points to handle numerical integration
273  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
274  {
275  // Now build the element matrix and right-hand-side using loops to
276  // integrate the test functions (i) against the trial functions (j).
277  for (std::size_t i=0; i<phi.size(); i++)
278  {
279  Fe(i) += JxW[qp]*phi[i][qp];
280 
281  for (std::size_t j=0; j<phi.size(); j++)
282  {
283  Ke(i,j) += JxW[qp]*(1.e-3*dphi[i][qp]*dphi[j][qp] +
284  phi[i][qp]*phi[j][qp]);
285  }
286  }
287  }
288 
289 
290  // At this point we have completed the matrix and RHS summation. The
291  // final step is to apply boundary conditions, which in this case are
292  // simple Dirichlet conditions with u(0) = u(1) = 0.
293 
294  // Define the penalty parameter used to enforce the BC's
295  double penalty = 1.e10;
296 
297  // Loop over the sides of this element. For a 1D element, the "sides"
298  // are defined as the nodes on each edge of the element, i.e. 1D elements
299  // have 2 sides.
300  for (auto s : elem->side_index_range())
301  {
302  // If this element has a NULL neighbor, then it is on the edge of the
303  // mesh and we need to enforce a boundary condition using the penalty
304  // method.
305  if (elem->neighbor_ptr(s) == libmesh_nullptr)
306  {
307  Ke(s,s) += penalty;
308  Fe(s) += 0*penalty;
309  }
310  }
311 
312  // This is a function call that is necessary when using adaptive
313  // mesh refinement. See Adaptivity Example 2 for more details.
314  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
315 
316  // Add Ke and Fe to the global matrix and right-hand-side.
317  system.matrix->add_matrix(Ke, dof_indices);
318  system.rhs->add_vector(Fe, dof_indices);
319  }
320 #endif // #ifdef LIBMESH_ENABLE_AMR
321 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
unsigned int dim
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
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]...
virtual Real l2_norm() const
Definition: statistics.C:36
virtual T maximum() const
Definition: statistics.C:61
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
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 ...
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
unsigned int & max_h_level()
The max_h_level is the greatest refinement level an element should reach.
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
This class implements writing meshes using GNUplot, designed for use only with 1D meshes...
Definition: gnuplot_io.h:43
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 DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
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:1795
virtual void reinit()
Reinitialize all the systems.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
void libmesh_ignore(const T &)
This class implements the Kelly error indicator which is based on the flux jumps between elements...
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
bool refine_and_coarsen_elements()
Refines and coarsens user-requested elements.
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
This function uses the derived class&#39;s jump error estimate formula to estimate the error on each cell...
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
const MeshBase & get_mesh() 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...
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.
const T_sys & get_system(const std::string &name) const
void assemble_1D(EquationSystems &es, const std::string &system_name)
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
unsigned int n_points() const
Definition: quadrature.h:113
int main(int argc, char **argv)
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.