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

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

Definition at line 183 of file adaptivity_ex1.C.

References libMesh::MeshBase::active_local_element_ptr_range(), 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::libmesh_ignore(), libmesh_nullptr, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, and libMesh::DenseVector< T >::size().

Referenced by main().

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
unsigned int dim
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]...
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
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
void libmesh_ignore(const T &)
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
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
const T_sys & get_system(const std::string &name) const
int main ( int  argc,
char **  argv 
)

Definition at line 59 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::LibMeshInit::comm(), 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, and libMesh::MeshOutput< MT >::write_equation_systems().

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 }
This is the EquationSystems class.
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.
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
MeshBase & mesh
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
SolverPackage default_solver_package()
Definition: libmesh.C:995
This class implements writing meshes using GNUplot, designed for use only with 1D meshes...
Definition: gnuplot_io.h:43
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
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
This class implements the Kelly error indicator which is based on the flux jumps between elements...
OStreamProxy out
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 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.
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