libMesh
miscellaneous_ex5.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>Miscellaneous Example 5 - Interior Penalty Discontinuous Galerkin</h1>
21 // \author Lorenzo Botti
22 // \date 2010
23 //
24 // This example is based on Adaptivity Example 3, but uses an
25 // Interior Penalty Discontinuous Galerkin formulation.
26 
27 #include <iostream>
28 
29 // LibMesh include files.
30 #include "libmesh/libmesh.h"
31 #include "libmesh/mesh.h"
32 #include "libmesh/equation_systems.h"
33 #include "libmesh/mesh_generation.h"
34 #include "libmesh/mesh_modification.h"
35 #include "libmesh/elem.h"
36 #include "libmesh/transient_system.h"
37 #include "libmesh/fe.h"
38 #include "libmesh/quadrature_gauss.h"
39 #include "libmesh/dof_map.h"
40 #include "libmesh/sparse_matrix.h"
41 #include "libmesh/dense_matrix.h"
42 #include "libmesh/dense_vector.h"
43 #include "libmesh/dense_submatrix.h"
44 #include "libmesh/dense_subvector.h"
45 #include "libmesh/numeric_vector.h"
46 #include "libmesh/linear_implicit_system.h"
47 #include "libmesh/exodusII_io.h"
48 #include "libmesh/fe_interface.h"
49 #include "libmesh/getpot.h"
50 #include "libmesh/mesh_refinement.h"
51 #include "libmesh/error_vector.h"
52 #include "libmesh/kelly_error_estimator.h"
53 #include "libmesh/discontinuity_measure.h"
54 #include "libmesh/string_to_enum.h"
55 
56 #include "libmesh/exact_solution.h"
57 //#define QORDER TWENTYSIXTH
58 
59 // Bring in everything from the libMesh namespace
60 using namespace libMesh;
61 
63  const Parameters & parameters,
64  const std::string &,
65  const std::string &)
66 {
67  const Real x = p(0);
68  const Real y = p(1);
69  const Real z = p(2);
70 
71  if (parameters.get<bool>("singularity"))
72  {
73  Real theta = atan2(y, x);
74 
75  if (theta < 0)
76  theta += 2. * libMesh::pi;
77 
78  return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
79  }
80  else
81  {
82  return cos(x) * exp(y) * (1. - z);
83  }
84 }
85 
86 // We now define the gradient of the exact solution, again being careful
87 // to obtain an angle from atan2 in the correct
88 // quadrant.
90  const Parameters & parameters, // es parameters
91  const std::string &, // sys_name, not needed
92  const std::string &) // unk_name, not needed
93 {
94  // Gradient value to be returned.
95  Gradient gradu;
96 
97  // x and y coordinates in space
98  const Real x = p(0);
99  const Real y = p(1);
100  const Real z = p(2);
101 
102  if (parameters.get<bool>("singularity"))
103  {
104  libmesh_assert_not_equal_to (x, 0.);
105 
106  // For convenience...
107  const Real tt = 2./3.;
108  const Real ot = 1./3.;
109 
110  // The value of the radius, squared
111  const Real r2 = x*x + y*y;
112 
113  // The boundary value, given by the exact solution,
114  // u_exact = r^(2/3)*sin(2*theta/3).
115  Real theta = atan2(y, x);
116 
117  // Make sure 0 <= theta <= 2*pi
118  if (theta < 0)
119  theta += 2. * libMesh::pi;
120 
121  // du/dx
122  gradu(0) = tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x;
123  gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
124  gradu(2) = 1.;
125  }
126  else
127  {
128  gradu(0) = -sin(x) * exp(y) * (1. - z);
129  gradu(1) = cos(x) * exp(y) * (1. - z);
130  gradu(2) = -cos(x) * exp(y);
131  }
132  return gradu;
133 }
134 
135 // We now define the matrix assembly function for the
136 // Laplace system. We need to first compute element volume
137 // matrices, and then take into account the boundary
138 // conditions and the flux integrals, which will be handled
139 // via an interior penalty method.
141  const std::string & libmesh_dbg_var(system_name))
142 {
143  libMesh::out << " assembling elliptic dg system... ";
145 
146  // It is a good idea to make sure we are assembling
147  // the proper system.
148  libmesh_assert_equal_to (system_name, "EllipticDG");
149 
150  // Get a constant reference to the mesh object.
151  const MeshBase & mesh = es.get_mesh();
152  // The dimension that we are running
153  const unsigned int dim = mesh.mesh_dimension();
154 
155  // Get a reference to the LinearImplicitSystem we are solving
156  LinearImplicitSystem & ellipticdg_system = es.get_system<LinearImplicitSystem> ("EllipticDG");
157  // Get some parameters that we need during assembly
158  const Real penalty = es.parameters.get<Real> ("penalty");
159  std::string refinement_type = es.parameters.get<std::string> ("refinement");
160 
161  // A reference to the DofMap object for this system. The DofMap
162  // object handles the index translation from node and element numbers
163  // to degree of freedom numbers. We will talk more about the DofMap
164  const DofMap & dof_map = ellipticdg_system.get_dof_map();
165 
166  // Get a constant reference to the Finite Element type
167  // for the first (and only) variable in the system.
168  FEType fe_type = ellipticdg_system.variable_type(0);
169 
170  // Build a Finite Element object of the specified type. Since the
171  // FEBase::build() member dynamically creates memory we will
172  // store the object as a UniquePtr<FEBase>. This can be thought
173  // of as a pointer that will clean up after itself.
174  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
175  UniquePtr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
176  UniquePtr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));
177 
178  // Quadrature rules for numerical integration.
179 #ifdef QORDER
180  QGauss qrule (dim, QORDER);
181 #else
182  QGauss qrule (dim, fe_type.default_quadrature_order());
183 #endif
184  fe->attach_quadrature_rule (&qrule);
185 
186 #ifdef QORDER
187  QGauss qface(dim-1, QORDER);
188 #else
189  QGauss qface(dim-1, fe_type.default_quadrature_order());
190 #endif
191 
192  // Tell the finite element object to use our quadrature rule.
193  fe_elem_face->attach_quadrature_rule(&qface);
194  fe_neighbor_face->attach_quadrature_rule(&qface);
195 
196  // Here we define some references to cell-specific data that
197  // will be used to assemble the linear system.
198  // Data for interior volume integrals
199  const std::vector<Real> & JxW = fe->get_JxW();
200  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
201 
202  // Data for surface integrals on the element boundary
203  const std::vector<std::vector<Real>> & phi_face = fe_elem_face->get_phi();
204  const std::vector<std::vector<RealGradient>> & dphi_face = fe_elem_face->get_dphi();
205  const std::vector<Real> & JxW_face = fe_elem_face->get_JxW();
206  const std::vector<Point> & qface_normals = fe_elem_face->get_normals();
207  const std::vector<Point> & qface_points = fe_elem_face->get_xyz();
208 
209  // Data for surface integrals on the neighbor boundary
210  const std::vector<std::vector<Real>> & phi_neighbor_face = fe_neighbor_face->get_phi();
211  const std::vector<std::vector<RealGradient>> & dphi_neighbor_face = fe_neighbor_face->get_dphi();
212 
213  // Define data structures to contain the element interior matrix
214  // and right-hand-side vector contribution. Following
215  // basic finite element terminology we will denote these
216  // "Ke" and "Fe".
219 
220  // Data structures to contain the element and neighbor boundary matrix
221  // contribution. This matrices will do the coupling between the dofs of
222  // the element and those of his neighbors.
223  // Ken: matrix coupling elem and neighbor dofs
228 
229  // This vector will hold the degree of freedom indices for
230  // the element. These define where in the global system
231  // the element degrees of freedom get mapped.
232  std::vector<dof_id_type> dof_indices;
233 
234  // Now we will loop over all the elements in the mesh. We will
235  // compute first the element interior matrix and right-hand-side contribution
236  // and then the element and neighbors boundary matrix contributions.
237  for (const auto & elem : mesh.active_local_element_ptr_range())
238  {
239  // Get the degree of freedom indices for the
240  // current element. These define where in the global
241  // matrix and right-hand-side this element will
242  // contribute to.
243  dof_map.dof_indices (elem, dof_indices);
244  const unsigned int n_dofs = dof_indices.size();
245 
246  // Compute the element-specific data for the current
247  // element. This involves computing the location of the
248  // quadrature points (q_point) and the shape functions
249  // (phi, dphi) for the current element.
250  fe->reinit (elem);
251 
252  // Zero the element matrix and right-hand side before
253  // summing them. We use the resize member here because
254  // the number of degrees of freedom might have changed from
255  // the last element.
256  Ke.resize (n_dofs, n_dofs);
257  Fe.resize (n_dofs);
258 
259  // Now we will build the element interior matrix. This involves
260  // a double loop to integrate the test functions (i) against
261  // the trial functions (j).
262  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
263  for (unsigned int i=0; i<n_dofs; i++)
264  for (unsigned int j=0; j<n_dofs; j++)
265  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
266 
267  // Now we address boundary conditions.
268  // We consider Dirichlet bc imposed via the interior penalty method
269  // The following loops over the sides of the element.
270  // If the element has no neighbor on a side then that
271  // side MUST live on a boundary of the domain.
272  for (auto side : elem->side_index_range())
273  {
274  if (elem->neighbor_ptr(side) == libmesh_nullptr)
275  {
276  // Pointer to the element face
277  fe_elem_face->reinit(elem, side);
278 
279  UniquePtr<const Elem> elem_side (elem->build_side_ptr(side));
280  // h element dimension to compute the interior penalty penalty parameter
281  const unsigned int elem_b_order = static_cast<unsigned int> (fe_elem_face->get_order());
282  const double h_elem = elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);
283 
284  for (unsigned int qp=0; qp<qface.n_points(); qp++)
285  {
286  Number bc_value = exact_solution(qface_points[qp], es.parameters, "null", "void");
287  for (unsigned int i=0; i<n_dofs; i++)
288  {
289  // Matrix contribution
290  for (unsigned int j=0; j<n_dofs; j++)
291  {
292  // stability
293  Ke(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[i][qp] * phi_face[j][qp];
294 
295  // consistency
296  Ke(i,j) -=
297  JxW_face[qp] *
298  (phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) +
299  phi_face[j][qp] * (dphi_face[i][qp]*qface_normals[qp]));
300  }
301 
302  // RHS contributions
303 
304  // stability
305  Fe(i) += JxW_face[qp] * bc_value * penalty/h_elem * phi_face[i][qp];
306 
307  // consistency
308  Fe(i) -= JxW_face[qp] * dphi_face[i][qp] * (bc_value*qface_normals[qp]);
309  }
310  }
311  }
312 
313  // If the element is not on a boundary of the domain
314  // we loop over his neighbors to compute the element
315  // and neighbor boundary matrix contributions
316  else
317  {
318  // Store a pointer to the neighbor we are currently
319  // working on.
320  const Elem * neighbor = elem->neighbor_ptr(side);
321 
322  // Get the global id of the element and the neighbor
323  const unsigned int elem_id = elem->id();
324  const unsigned int neighbor_id = neighbor->id();
325 
326  // If the neighbor has the same h level and is active
327  // perform integration only if our global id is bigger than our neighbor id.
328  // We don't want to compute twice the same contributions.
329  // If the neighbor has a different h level perform integration
330  // only if the neighbor is at a lower level.
331  if ((neighbor->active() &&
332  (neighbor->level() == elem->level()) &&
333  (elem_id < neighbor_id)) ||
334  (neighbor->level() < elem->level()))
335  {
336  // Pointer to the element side
337  UniquePtr<const Elem> elem_side (elem->build_side_ptr(side));
338 
339  // h dimension to compute the interior penalty penalty parameter
340  const unsigned int elem_b_order = static_cast<unsigned int>(fe_elem_face->get_order());
341  const unsigned int neighbor_b_order = static_cast<unsigned int>(fe_neighbor_face->get_order());
342  const double side_order = (elem_b_order + neighbor_b_order)/2.;
343  const double h_elem = (elem->volume()/elem_side->volume()) * 1./pow(side_order,2.);
344 
345  // The quadrature point locations on the neighbor side
346  std::vector<Point> qface_neighbor_point;
347 
348  // The quadrature point locations on the element side
349  std::vector<Point > qface_point;
350 
351  // Reinitialize shape functions on the element side
352  fe_elem_face->reinit(elem, side);
353 
354  // Get the physical locations of the element quadrature points
355  qface_point = fe_elem_face->get_xyz();
356 
357  // Find their locations on the neighbor
358  unsigned int side_neighbor = neighbor->which_neighbor_am_i(elem);
359  if (refinement_type == "p")
360  fe_neighbor_face->side_map (neighbor,
361  elem_side.get(),
362  side_neighbor,
363  qface.get_points(),
364  qface_neighbor_point);
365  else
366  FEInterface::inverse_map (elem->dim(),
367  fe->get_fe_type(),
368  neighbor,
369  qface_point,
370  qface_neighbor_point);
371 
372  // Calculate the neighbor element shape functions at those locations
373  fe_neighbor_face->reinit(neighbor, &qface_neighbor_point);
374 
375  // Get the degree of freedom indices for the
376  // neighbor. These define where in the global
377  // matrix this neighbor will contribute to.
378  std::vector<dof_id_type> neighbor_dof_indices;
379  dof_map.dof_indices (neighbor, neighbor_dof_indices);
380  const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
381 
382  // Zero the element and neighbor side matrix before
383  // summing them. We use the resize member here because
384  // the number of degrees of freedom might have changed from
385  // the last element or neighbor.
386  // Note that Kne and Ken are not square matrices if neighbor
387  // and element have a different p level
388  Kne.resize (n_neighbor_dofs, n_dofs);
389  Ken.resize (n_dofs, n_neighbor_dofs);
390  Kee.resize (n_dofs, n_dofs);
391  Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
392 
393  // Now we will build the element and neighbor
394  // boundary matrices. This involves
395  // a double loop to integrate the test functions
396  // (i) against the trial functions (j).
397  for (unsigned int qp=0; qp<qface.n_points(); qp++)
398  {
399  // Kee Matrix. Integrate the element test function i
400  // against the element test function j
401  for (unsigned int i=0; i<n_dofs; i++)
402  {
403  for (unsigned int j=0; j<n_dofs; j++)
404  {
405  // consistency
406  Kee(i,j) -=
407  0.5 * JxW_face[qp] *
408  (phi_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) +
409  phi_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]));
410 
411  // stability
412  Kee(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_face[i][qp];
413  }
414  }
415 
416  // Knn Matrix. Integrate the neighbor test function i
417  // against the neighbor test function j
418  for (unsigned int i=0; i<n_neighbor_dofs; i++)
419  {
420  for (unsigned int j=0; j<n_neighbor_dofs; j++)
421  {
422  // consistency
423  Knn(i,j) +=
424  0.5 * JxW_face[qp] *
425  (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]) +
426  phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp]));
427 
428  // stability
429  Knn(i,j) +=
430  JxW_face[qp] * penalty/h_elem * phi_neighbor_face[j][qp]*phi_neighbor_face[i][qp];
431  }
432  }
433 
434  // Kne Matrix. Integrate the neighbor test function i
435  // against the element test function j
436  for (unsigned int i=0; i<n_neighbor_dofs; i++)
437  {
438  for (unsigned int j=0; j<n_dofs; j++)
439  {
440  // consistency
441  Kne(i,j) +=
442  0.5 * JxW_face[qp] *
443  (phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]) -
444  phi_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]));
445 
446  // stability
447  Kne(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_neighbor_face[i][qp];
448  }
449  }
450 
451  // Ken Matrix. Integrate the element test function i
452  // against the neighbor test function j
453  for (unsigned int i=0; i<n_dofs; i++)
454  {
455  for (unsigned int j=0; j<n_neighbor_dofs; j++)
456  {
457  // consistency
458  Ken(i,j) +=
459  0.5 * JxW_face[qp] *
460  (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) -
461  phi_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp]));
462 
463  // stability
464  Ken(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[i][qp]*phi_neighbor_face[j][qp];
465  }
466  }
467  }
468 
469  // The element and neighbor boundary matrix are now built
470  // for this side. Add them to the global matrix
471  // The SparseMatrix::add_matrix() members do this for us.
472  ellipticdg_system.matrix->add_matrix(Kne, neighbor_dof_indices, dof_indices);
473  ellipticdg_system.matrix->add_matrix(Ken, dof_indices, neighbor_dof_indices);
474  ellipticdg_system.matrix->add_matrix(Kee, dof_indices);
475  ellipticdg_system.matrix->add_matrix(Knn, neighbor_dof_indices);
476  }
477  }
478  }
479  // The element interior matrix and right-hand-side are now built
480  // for this element. Add them to the global matrix and
481  // right-hand-side vector. The SparseMatrix::add_matrix()
482  // and NumericVector::add_vector() members do this for us.
483  ellipticdg_system.matrix->add_matrix(Ke, dof_indices);
484  ellipticdg_system.rhs->add_vector(Fe, dof_indices);
485  }
486 
487  libMesh::out << "done" << std::endl;
488 }
489 
490 
491 
492 int main (int argc, char** argv)
493 {
494  LibMeshInit init(argc, argv);
495 
496  // This example requires a linear solver package.
497  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
498  "--enable-petsc, --enable-trilinos, or --enable-eigen");
499 
500  // Skip adaptive examples on a non-adaptive libMesh build
501 #ifndef LIBMESH_ENABLE_AMR
502  libmesh_example_requires(false, "--enable-amr");
503 #else
504 
505  //Parse the input file
506  GetPot input_file("miscellaneous_ex5.in");
507 
508  //Read in parameters from the input file
509  const unsigned int adaptive_refinement_steps = input_file("max_adaptive_r_steps", 3);
510  const unsigned int uniform_refinement_steps = input_file("uniform_h_r_steps", 3);
511  const Real refine_fraction = input_file("refine_fraction", 0.5);
512  const Real coarsen_fraction = input_file("coarsen_fraction", 0.);
513  const unsigned int max_h_level = input_file("max_h_level", 10);
514  const std::string refinement_type = input_file("refinement_type","p");
515  Order p_order = static_cast<Order>(input_file("p_order", 1));
516  const std::string element_type = input_file("element_type", "tensor");
517  const Real penalty = input_file("ip_penalty", 10.);
518  const bool singularity = input_file("singularity", true);
519  const unsigned int dim = input_file("dimension", 3);
520 
521  // Skip higher-dimensional examples on a lower-dimensional libMesh build
522  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
523 
524 
525  // Create a mesh, with dimension to be overridden later, distributed
526  // across the default MPI communicator.
527  Mesh mesh(init.comm());
528 
529  if (dim == 1)
531  else if (dim == 2)
532  mesh.read("lshaped.xda");
533  else
534  mesh.read ("lshaped3D.xda");
535 
536  // Use triangles if the config file says so
537  if (element_type == "simplex")
539 
540  // Mesh Refinement object
541  MeshRefinement mesh_refinement(mesh);
542  mesh_refinement.refine_fraction() = refine_fraction;
543  mesh_refinement.coarsen_fraction() = coarsen_fraction;
544  mesh_refinement.max_h_level() = max_h_level;
545 
546  // Do uniform refinement
547  for (unsigned int rstep=0; rstep<uniform_refinement_steps; rstep++)
548  mesh_refinement.uniformly_refine(1);
549 
550  // Crate an equation system object
551  EquationSystems equation_system (mesh);
552 
553  // Set parameters for the equation system and the solver
554  equation_system.parameters.set<Real>("linear solver tolerance") = TOLERANCE * TOLERANCE;
555  equation_system.parameters.set<unsigned int>("linear solver maximum iterations") = 1000;
556  equation_system.parameters.set<Real>("penalty") = penalty;
557  equation_system.parameters.set<bool>("singularity") = singularity;
558  equation_system.parameters.set<std::string>("refinement") = refinement_type;
559 
560  // Create a system named ellipticdg
561  LinearImplicitSystem & ellipticdg_system = equation_system.add_system<LinearImplicitSystem> ("EllipticDG");
562 
563  // Add a variable "u" to "ellipticdg" using the p_order specified in the config file
564  if (on_command_line("element_type"))
565  {
566  std::string fe_str =
567  command_line_value(std::string("element_type"),
568  std::string("MONOMIAL"));
569 
570  if (fe_str != "MONOMIAL" || fe_str != "XYZ")
571  libmesh_error_msg("Error: This example must be run with MONOMIAL or XYZ element types.");
572 
573  ellipticdg_system.add_variable ("u", p_order, Utility::string_to_enum<FEFamily>(fe_str));
574  }
575  else
576  ellipticdg_system.add_variable ("u", p_order, MONOMIAL);
577 
578  // Give the system a pointer to the matrix assembly function
579  ellipticdg_system.attach_assemble_function (assemble_ellipticdg);
580 
581  // Initialize the data structures for the equation system
582  equation_system.init();
583 
584  // Construct ExactSolution object and attach solution functions
585  ExactSolution exact_sol(equation_system);
588 
589  // A refinement loop.
590  for (unsigned int rstep=0; rstep<adaptive_refinement_steps; ++rstep)
591  {
592  libMesh::out << " Beginning Solve " << rstep << std::endl;
593  libMesh::out << "Number of elements: " << mesh.n_elem() << std::endl;
594 
595  // Solve the system
596  ellipticdg_system.solve();
597 
598  libMesh::out << "System has: "
599  << equation_system.n_active_dofs()
600  << " degrees of freedom."
601  << std::endl;
602 
603  libMesh::out << "Linear solver converged at step: "
604  << ellipticdg_system.n_linear_iterations()
605  << ", final residual: "
606  << ellipticdg_system.final_linear_residual()
607  << std::endl;
608 
609  // Compute the error
610  exact_sol.compute_error("EllipticDG", "u");
611 
612  // Print out the error values
613  libMesh::out << "L2-Error is: "
614  << exact_sol.l2_error("EllipticDG", "u")
615  << std::endl;
616 
617  // Possibly refine the mesh
618  if (rstep+1 < adaptive_refinement_steps)
619  {
620  // The ErrorVector is a particular StatisticsVector
621  // for computing error information on a finite element mesh.
622  ErrorVector error;
623 
624  // The discontinuity error estimator
625  // evaluate the jump of the solution
626  // on elements faces
627  DiscontinuityMeasure error_estimator;
628  error_estimator.estimate_error(ellipticdg_system,error);
629 
630  // Take the error in error and decide which elements will be coarsened or refined
631  mesh_refinement.flag_elements_by_error_fraction(error);
632  if (refinement_type == "p")
633  mesh_refinement.switch_h_to_p_refinement();
634  if (refinement_type == "hp")
635  mesh_refinement.add_p_to_h_refinement();
636 
637  // Refine and coarsen the flagged elements
638  mesh_refinement.refine_and_coarsen_elements();
639  equation_system.reinit();
640  }
641  }
642 
643  // Write out the solution
644  // After solving the system write the solution
645  // to a ExodusII-formatted plot file.
646 #ifdef LIBMESH_HAVE_EXODUS_API
647  ExodusII_IO (mesh).write_discontinuous_exodusII("lshaped_dg.e", equation_system);
648 #endif
649 
650 #endif // #ifndef LIBMESH_ENABLE_AMR
651 
652  // All done.
653  return 0;
654 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
This is the EquationSystems class.
bool active() const
Definition: elem.h:2257
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
unsigned int dim
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
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]...
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition: elem.h:2181
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
MeshBase & mesh
int main(int argc, char **argv)
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 ...
const T & get(const std::string &) const
Definition: parameters.h:430
static const Real TOLERANCE
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.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
bool singularity
unsigned int & max_h_level()
The max_h_level is the greatest refinement level an element should reach.
void compute_error(const std::string &sys_name, const std::string &unknown_name)
Computes and stores the error in the solution value e = u-u_h, the gradient grad(e) = grad(u) - grad(...
This class measures discontinuities between elements for debugging purposes.
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
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
void switch_h_to_p_refinement()
Takes a mesh whose elements are flagged for h refinement and coarsening, and switches those flags to ...
void all_tri(MeshBase &mesh)
Converts the 2D quadrilateral elements of a Mesh into triangular elements.
Gradient exact_derivative(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
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.
PetscErrorCode Vec x
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
T command_line_value(const std::string &, T)
Definition: libmesh.C:932
Order default_quadrature_order() const
Definition: fe_type.h:332
virtual void reinit()
Reinitialize all the systems.
Number exact_solution(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
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.
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
double pow(double a, int b)
void attach_exact_deriv(unsigned int sys_num, FunctionBase< Gradient > *g)
Clone and attach an arbitrary functor which computes the exact gradient of the system sys_num solutio...
BasicOStreamProxy & flush()
Flush the associated stream buffer.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
T & set(const std::string &)
Definition: parameters.h:469
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
const std::vector< Point > & get_points() const
Definition: quadrature.h:128
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
unsigned int level() const
Definition: elem.h:2388
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
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.
void assemble_ellipticdg(EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
const T_sys & get_system(const std::string &name) const
void attach_exact_value(unsigned int sys_num, FunctionBase< Number > *f)
Clone and attach an arbitrary functor which computes the exact value of the system sys_num solution a...
void add_p_to_h_refinement()
Takes a mesh whose elements are flagged for h refinement and coarsening, and adds flags to request p ...
Real size() const
Definition: type_vector.h:898
virtual void init()
Initialize all the systems.
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
std::size_t n_active_dofs() const
dof_id_type id() const
Definition: dof_object.h:632
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual dof_id_type n_elem() const =0
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113
void write_discontinuous_exodusII(const std::string &name, const EquationSystems &es, const std::set< std::string > *system_names=libmesh_nullptr)
Writes a exodusII file with discontinuous data.
Definition: exodusII_io.C:89
const Real pi
.
Definition: libmesh.h:172
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.