libMesh
adaptivity_ex3.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 3 - Laplace Equation in the L-Shaped Domain</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This example solves the Laplace equation on the classic "L-shaped"
25 // domain with adaptive mesh refinement. In this case, the exact
26 // solution is u(r,\theta) = r^{2/3} * \sin ( (2/3) * \theta), but
27 // the standard Kelly error indicator is used to estimate the error.
28 // The initial mesh contains three QUAD9 elements which represent the
29 // standard quadrants I, II, and III of the domain [-1,1]x[-1,1],
30 // i.e.
31 // Element 0: [-1,0]x[ 0,1]
32 // Element 1: [ 0,1]x[ 0,1]
33 // Element 2: [-1,0]x[-1,0]
34 // The mesh is provided in the standard libMesh ASCII format file
35 // named "lshaped.xda". In addition, an input file named "adaptivity_ex3.in"
36 // is provided which allows the user to set several parameters for
37 // the solution so that the problem can be re-run without a
38 // re-compile. The solution technique employed is to have a
39 // refinement loop with a linear solve inside followed by a
40 // refinement of the grid and projection of the solution to the new grid
41 // In the final loop iteration, there is no additional
42 // refinement after the solve. In the input file "adaptivity_ex3.in",
43 // the variable "max_r_steps" controls the number of refinement steps,
44 // "max_r_level" controls the maximum element refinement level, and
45 // "refine_percentage" / "coarsen_percentage" determine the number of
46 // elements which will be refined / coarsened at each step.
47 
48 // LibMesh include files.
49 #include "libmesh/mesh.h"
50 #include "libmesh/equation_systems.h"
51 #include "libmesh/linear_implicit_system.h"
52 #include "libmesh/exodusII_io.h"
53 #include "libmesh/tecplot_io.h"
54 #include "libmesh/fe.h"
55 #include "libmesh/quadrature_gauss.h"
56 #include "libmesh/dense_matrix.h"
57 #include "libmesh/dense_vector.h"
58 #include "libmesh/sparse_matrix.h"
59 #include "libmesh/mesh_refinement.h"
60 #include "libmesh/error_vector.h"
61 #include "libmesh/discontinuity_measure.h"
62 #include "libmesh/exact_error_estimator.h"
63 #include "libmesh/kelly_error_estimator.h"
64 #include "libmesh/patch_recovery_error_estimator.h"
65 #include "libmesh/uniform_refinement_estimator.h"
66 #include "libmesh/hp_coarsentest.h"
67 #include "libmesh/hp_singular.h"
68 #include "libmesh/sibling_coupling.h"
69 #include "libmesh/mesh_generation.h"
70 #include "libmesh/mesh_modification.h"
71 #include "libmesh/perf_log.h"
72 #include "libmesh/getpot.h"
73 #include "libmesh/exact_solution.h"
74 #include "libmesh/dof_map.h"
75 #include "libmesh/numeric_vector.h"
76 #include "libmesh/elem.h"
77 #include "libmesh/string_to_enum.h"
78 
79 // Bring in everything from the libMesh namespace
80 using namespace libMesh;
81 
82 // Function prototype. This is the function that will assemble
83 // the linear system for our Laplace problem. Note that the
84 // function will take the EquationSystems object and the
85 // name of the system we are assembling as input. From the
86 // EquationSystems object we have access to the Mesh and
87 // other objects we might need.
89  const std::string & system_name);
90 
91 
92 // Prototype for calculation of the exact solution. Useful
93 // for setting boundary conditions.
94 Number exact_solution(const Point & p,
95  const Parameters &, // EquationSystem parameters, not needed
96  const std::string &, // sys_name, not needed
97  const std::string &); // unk_name, not needed);
98 
99 // Prototype for calculation of the gradient of the exact solution.
101  const Parameters &, // EquationSystems parameters, not needed
102  const std::string &, // sys_name, not needed
103  const std::string &); // unk_name, not needed);
104 
105 
106 // These are non-const because the input file may change it,
107 // It is global because our exact_* functions use it.
108 
109 // Set the dimensionality of the mesh
110 unsigned int dim = 2;
111 
112 // Choose whether or not to use the singular solution
113 bool singularity = true;
114 
115 
116 int main(int argc, char ** argv)
117 {
118  // Initialize libMesh.
119  LibMeshInit init (argc, argv);
120 
121  // This example requires a linear solver package.
122  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
123  "--enable-petsc, --enable-trilinos, or --enable-eigen");
124 
125  // Single precision is inadequate for p refinement
126  libmesh_example_requires(sizeof(Real) > 4, "--disable-singleprecision");
127 
128  // Skip adaptive examples on a non-adaptive libMesh build
129 #ifndef LIBMESH_ENABLE_AMR
130  libmesh_example_requires(false, "--enable-amr");
131 #else
132 
133  // Parse the input file
134  GetPot input_file("adaptivity_ex3.in");
135 
136  // But allow the command line to override it.
137  input_file.parse_command_line(argc, argv);
138 
139  // Read in parameters from the input file
140  const unsigned int max_r_steps = input_file("max_r_steps", 3);
141  const unsigned int max_r_level = input_file("max_r_level", 3);
142  const Real refine_percentage = input_file("refine_percentage", 0.5);
143  const Real coarsen_percentage = input_file("coarsen_percentage", 0.5);
144  const unsigned int uniform_refine = input_file("uniform_refine",0);
145  const std::string refine_type = input_file("refinement_type", "h");
146  const std::string approx_type = input_file("approx_type", "LAGRANGE");
147  const unsigned int approx_order = input_file("approx_order", 1);
148  const std::string element_type = input_file("element_type", "tensor");
149  const int extra_error_quadrature = input_file("extra_error_quadrature", 0);
150  const int max_linear_iterations = input_file("max_linear_iterations", 5000);
151 
152 #ifdef LIBMESH_HAVE_EXODUS_API
153  const bool output_intermediate = input_file("output_intermediate", false);
154 #endif
155 
156  // If libmesh is configured without second derivative support, we
157  // can't run this example with Hermite elements and will therefore
158  // fail gracefully.
159 #if !defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
160  libmesh_example_requires(approx_type != "HERMITE", "--enable-second");
161 #endif
162 
163  dim = input_file("dimension", 2);
164  const std::string indicator_type = input_file("indicator_type", "kelly");
165  singularity = input_file("singularity", true);
166 
167  // Skip higher-dimensional examples on a lower-dimensional libMesh build
168  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
169 
170  // Output file for plotting the error as a function of
171  // the number of degrees of freedom.
172  std::string approx_name = "";
173  if (element_type == "tensor")
174  approx_name += "bi";
175  if (approx_order == 1)
176  approx_name += "linear";
177  else if (approx_order == 2)
178  approx_name += "quadratic";
179  else if (approx_order == 3)
180  approx_name += "cubic";
181  else if (approx_order == 4)
182  approx_name += "quartic";
183 
184  std::string output_file = approx_name;
185  output_file += "_";
186  output_file += refine_type;
187  if (uniform_refine == 0)
188  output_file += "_adaptive.m";
189  else
190  output_file += "_uniform.m";
191 
192  std::ofstream out (output_file.c_str());
193  out << "% dofs L2-error H1-error" << std::endl;
194  out << "e = [" << std::endl;
195 
196  // Create a mesh, with dimension to be overridden later, on the
197  // default MPI communicator.
198  Mesh mesh(init.comm());
199 
200  // Create an equation systems object.
201  EquationSystems equation_systems (mesh);
202 
203  // Declare the system we will solve, named Laplace
204  LinearImplicitSystem & system =
205  equation_systems.add_system<LinearImplicitSystem> ("Laplace");
206 
207  // If we're doing HP refinement, then we'll need to be able to
208  // evaluate data on elements' siblings in HPCoarsenTest, which means
209  // we should instruct our system's DofMap to distribute that data.
210  // Create and add (and keep around! the DofMap will only be holding
211  // a pointer to it!) a SiblingCoupling functor to provide that
212  // instruction.
213  SiblingCoupling sibling_coupling;
214  if (refine_type == "hp")
215  system.get_dof_map().add_algebraic_ghosting_functor(sibling_coupling);
216 
217  // Read in the mesh
218  if (dim == 1)
220  else if (dim == 2)
221  mesh.read("lshaped.xda");
222  else
223  mesh.read("lshaped3D.xda");
224 
225  // Use triangles if the config file says so
226  if (element_type == "simplex")
228 
229  // We used first order elements to describe the geometry,
230  // but we may need second order elements to hold the degrees
231  // of freedom
232  if (approx_order > 1 || refine_type != "h")
234 
235  // Mesh Refinement object
236  MeshRefinement mesh_refinement(mesh);
237  mesh_refinement.refine_fraction() = refine_percentage;
238  mesh_refinement.coarsen_fraction() = coarsen_percentage;
239  mesh_refinement.max_h_level() = max_r_level;
240 
241  // Adds the variable "u" to "Laplace", using
242  // the finite element type and order specified
243  // in the config file
244  system.add_variable("u", static_cast<Order>(approx_order),
245  Utility::string_to_enum<FEFamily>(approx_type));
246 
247  // Give the system a pointer to the matrix assembly
248  // function.
250 
251  // Initialize the data structures for the equation system.
252  equation_systems.init();
253 
254  // Set linear solver max iterations
255  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations")
256  = max_linear_iterations;
257 
258  // Linear solver tolerance.
259  equation_systems.parameters.set<Real>("linear solver tolerance") =
260  std::pow(TOLERANCE, 2.5);
261 
262  // Prints information about the system to the screen.
263  equation_systems.print_info();
264 
265  // Construct ExactSolution object and attach solution functions
266  ExactSolution exact_sol(equation_systems);
269 
270  // Use higher quadrature order for more accurate error results
271  exact_sol.extra_quadrature_order(extra_error_quadrature);
272 
273  // A refinement loop.
274  for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
275  {
276  libMesh::out << "Beginning Solve " << r_step << std::endl;
277 
278  // Solve the system "Laplace", just like example 2.
279  system.solve();
280 
281  libMesh::out << "System has: "
282  << equation_systems.n_active_dofs()
283  << " degrees of freedom."
284  << std::endl;
285 
286  libMesh::out << "Linear solver converged at step: "
287  << system.n_linear_iterations()
288  << ", final residual: "
289  << system.final_linear_residual()
290  << std::endl;
291 
292 #ifdef LIBMESH_HAVE_EXODUS_API
293  // After solving the system write the solution
294  // to a ExodusII-formatted plot file.
295  if (output_intermediate)
296  {
297  std::ostringstream outfile;
298  outfile << "lshaped_" << r_step << ".e";
299  ExodusII_IO (mesh).write_equation_systems (outfile.str(),
300  equation_systems);
301  }
302 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
303 
304  // Compute the error.
305  exact_sol.compute_error("Laplace", "u");
306 
307  // Print out the error values
308  libMesh::out << "L2-Error is: "
309  << exact_sol.l2_error("Laplace", "u")
310  << std::endl;
311  libMesh::out << "H1-Error is: "
312  << exact_sol.h1_error("Laplace", "u")
313  << std::endl;
314 
315  // Compute any discontinuity. There should be none.
316  {
318  ErrorVector disc_error;
319  disc.estimate_error(system, disc_error);
320 
321  Real mean_disc_error = disc_error.mean();
322 
323  libMesh::out << "Mean discontinuity error = " << mean_disc_error << std::endl;
324 
325  // FIXME - this test fails when solving with Eigen?
326 #ifdef LIBMESH_ENABLE_PETSC
327  libmesh_assert_less (mean_disc_error, 1e-14);
328 #endif
329  }
330 
331  // Print to output file
332  out << equation_systems.n_active_dofs() << " "
333  << exact_sol.l2_error("Laplace", "u") << " "
334  << exact_sol.h1_error("Laplace", "u") << std::endl;
335 
336  // Possibly refine the mesh
337  if (r_step+1 != max_r_steps)
338  {
339  libMesh::out << " Refining the mesh..." << std::endl;
340 
341  if (uniform_refine == 0)
342  {
343 
344  // The ErrorVector is a particular StatisticsVector
345  // for computing error information on a finite element mesh.
346  ErrorVector error;
347 
348  if (indicator_type == "exact")
349  {
350  // The ErrorEstimator class interrogates a
351  // finite element solution and assigns to each
352  // element a positive error value.
353  // This value is used for deciding which elements to
354  // refine and which to coarsen.
355  // For these simple test problems, we can use
356  // numerical quadrature of the exact error between
357  // the approximate and analytic solutions.
358  // However, for real problems, we would need an error
359  // indicator which only relies on the approximate
360  // solution.
361  ExactErrorEstimator error_estimator;
362 
363  error_estimator.attach_exact_value(exact_solution);
364  error_estimator.attach_exact_deriv(exact_derivative);
365 
366  // We optimize in H1 norm, the default
367  // error_estimator.error_norm = H1;
368 
369  // Compute the error for each active element using
370  // the provided indicator. Note in general you
371  // will need to provide an error estimator
372  // specifically designed for your application.
373  error_estimator.estimate_error (system, error);
374  }
375  else if (indicator_type == "patch")
376  {
377  // The patch recovery estimator should give a
378  // good estimate of the solution interpolation
379  // error.
380  PatchRecoveryErrorEstimator error_estimator;
381 
382  error_estimator.estimate_error (system, error);
383  }
384  else if (indicator_type == "uniform")
385  {
386  // Error indication based on uniform refinement
387  // is reliable, but very expensive.
388  UniformRefinementEstimator error_estimator;
389 
390  error_estimator.estimate_error (system, error);
391  }
392  else
393  {
394  libmesh_assert_equal_to (indicator_type, "kelly");
395 
396  // The Kelly error estimator is based on
397  // an error bound for the Poisson problem
398  // on linear elements, but is useful for
399  // driving adaptive refinement in many problems
400  KellyErrorEstimator error_estimator;
401 
402  error_estimator.estimate_error (system, error);
403  }
404 
405  // Write out the error distribution
406  std::ostringstream ss;
407  ss << r_step;
408 #ifdef LIBMESH_HAVE_EXODUS_API
409  std::string error_output = "error_" + ss.str() + ".e";
410 #else
411  std::string error_output = "error_" + ss.str() + ".gmv";
412 #endif
413  error.plot_error(error_output, mesh);
414 
415  // This takes the error in error and decides which elements
416  // will be coarsened or refined. Any element within 20% of the
417  // maximum error on any element will be refined, and any
418  // element within 10% of the minimum error on any element might
419  // be coarsened. Note that the elements flagged for refinement
420  // will be refined, but those flagged for coarsening _might_ be
421  // coarsened.
422  mesh_refinement.flag_elements_by_error_fraction (error);
423 
424  // If we are doing adaptive p refinement, we want
425  // elements flagged for that instead.
426  if (refine_type == "p")
427  mesh_refinement.switch_h_to_p_refinement();
428  // If we are doing "matched hp" refinement, we
429  // flag elements for both h and p
430  else if (refine_type == "matchedhp")
431  mesh_refinement.add_p_to_h_refinement();
432  // If we are doing hp refinement, we
433  // try switching some elements from h to p
434  else if (refine_type == "hp")
435  {
436  HPCoarsenTest hpselector;
437  hpselector.select_refinement(system);
438  }
439  // If we are doing "singular hp" refinement, we
440  // try switching most elements from h to p
441  else if (refine_type == "singularhp")
442  {
443  // This only differs from p refinement for
444  // the singular problem
446  HPSingularity hpselector;
447  // Our only singular point is at the origin
448  hpselector.singular_points.push_back(Point());
449  hpselector.select_refinement(system);
450  }
451  else if (refine_type != "h")
452  libmesh_error_msg("Unknown refinement_type = " << refine_type);
453 
454  // This call actually refines and coarsens the flagged
455  // elements.
456  mesh_refinement.refine_and_coarsen_elements();
457  }
458 
459  else if (uniform_refine == 1)
460  {
461  if (refine_type == "h" || refine_type == "hp" ||
462  refine_type == "matchedhp")
463  mesh_refinement.uniformly_refine(1);
464  if (refine_type == "p" || refine_type == "hp" ||
465  refine_type == "matchedhp")
466  mesh_refinement.uniformly_p_refine(1);
467  }
468 
469  // This call reinitializes the EquationSystems object for
470  // the newly refined mesh. One of the steps in the
471  // reinitialization is projecting the solution,
472  // old_solution, etc... vectors from the old mesh to
473  // the current one.
474  equation_systems.reinit ();
475  }
476  }
477 
478 #ifdef LIBMESH_HAVE_EXODUS_API
479  // Write out the solution
480  // After solving the system write the solution
481  // to a ExodusII-formatted plot file.
482  ExodusII_IO (mesh).write_equation_systems ("lshaped.e",
483  equation_systems);
484 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
485 
486  // Close up the output file.
487  out << "];" << std::endl;
488  out << "hold on" << std::endl;
489  out << "plot(e(:,1), e(:,2), 'bo-');" << std::endl;
490  out << "plot(e(:,1), e(:,3), 'ro-');" << std::endl;
491  // out << "set(gca,'XScale', 'Log');" << std::endl;
492  // out << "set(gca,'YScale', 'Log');" << std::endl;
493  out << "xlabel('dofs');" << std::endl;
494  out << "title('" << approx_name << " elements');" << std::endl;
495  out << "legend('L2-error', 'H1-error');" << std::endl;
496  // out << "disp('L2-error linear fit');" << std::endl;
497  // out << "polyfit(log10(e(:,1)), log10(e(:,2)), 1)" << std::endl;
498  // out << "disp('H1-error linear fit');" << std::endl;
499  // out << "polyfit(log10(e(:,1)), log10(e(:,3)), 1)" << std::endl;
500 #endif // #ifndef LIBMESH_ENABLE_AMR
501 
502  // All done.
503  return 0;
504 }
505 
506 
507 
508 
509 // We now define the exact solution, being careful
510 // to obtain an angle from atan2 in the correct
511 // quadrant.
513  const Parameters &, // parameters, not needed
514  const std::string &, // sys_name, not needed
515  const std::string &) // unk_name, not needed
516 {
517  const Real x = p(0);
518  const Real y = (dim > 1) ? p(1) : 0.;
519 
520  if (singularity)
521  {
522  // The exact solution to the singular problem,
523  // u_exact = r^(2/3)*sin(2*theta/3).
524  Real theta = atan2(y, x);
525 
526  // Make sure 0 <= theta <= 2*pi
527  if (theta < 0)
528  theta += 2. * libMesh::pi;
529 
530  // Make the 3D solution similar
531  const Real z = (dim > 2) ? p(2) : 0;
532 
533  return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
534  }
535  else
536  {
537  // The exact solution to a nonsingular problem,
538  // good for testing ideal convergence rates
539  const Real z = (dim > 2) ? p(2) : 0;
540 
541  return cos(x) * exp(y) * (1. - z);
542  }
543 }
544 
545 
546 
547 
548 
549 // We now define the gradient of the exact solution, again being careful
550 // to obtain an angle from atan2 in the correct
551 // quadrant.
553  const Parameters &, // parameters, not needed
554  const std::string &, // sys_name, not needed
555  const std::string &) // unk_name, not needed
556 {
557  // Gradient value to be returned.
558  Gradient gradu;
559 
560  // x and y coordinates in space
561  const Real x = p(0);
562  const Real y = dim > 1 ? p(1) : 0.;
563 
564  if (singularity)
565  {
566  // We can't compute the gradient at x=0, it is not defined.
567  libmesh_assert_not_equal_to (x, 0.);
568 
569  // For convenience...
570  const Real tt = 2./3.;
571  const Real ot = 1./3.;
572 
573  // The value of the radius, squared
574  const Real r2 = x*x + y*y;
575 
576  // The boundary value, given by the exact solution,
577  // u_exact = r^(2/3)*sin(2*theta/3).
578  Real theta = atan2(y, x);
579 
580  // Make sure 0 <= theta <= 2*pi
581  if (theta < 0)
582  theta += 2. * libMesh::pi;
583 
584  // du/dx
585  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;
586 
587  // du/dy
588  if (dim > 1)
589  gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
590 
591  if (dim > 2)
592  gradu(2) = 1.;
593  }
594  else
595  {
596  const Real z = (dim > 2) ? p(2) : 0;
597 
598  gradu(0) = -sin(x) * exp(y) * (1. - z);
599  if (dim > 1)
600  gradu(1) = cos(x) * exp(y) * (1. - z);
601  if (dim > 2)
602  gradu(2) = -cos(x) * exp(y);
603  }
604 
605  return gradu;
606 }
607 
608 
609 
610 
611 
612 
613 // We now define the matrix assembly function for the
614 // Laplace system. We need to first compute element
615 // matrices and right-hand sides, and then take into
616 // account the boundary conditions, which will be handled
617 // via a penalty method.
619  const std::string & system_name)
620 {
621  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
622  libmesh_ignore(es);
623  libmesh_ignore(system_name);
624 
625 #ifdef LIBMESH_ENABLE_AMR
626  // It is a good idea to make sure we are assembling
627  // the proper system.
628  libmesh_assert_equal_to (system_name, "Laplace");
629 
630  // Declare a performance log. Give it a descriptive
631  // string to identify what part of the code we are
632  // logging, since there may be many PerfLogs in an
633  // application.
634  PerfLog perf_log ("Matrix Assembly", false);
635 
636  // Get a constant reference to the mesh object.
637  const MeshBase & mesh = es.get_mesh();
638 
639  // The dimension that we are running
640  const unsigned int mesh_dim = mesh.mesh_dimension();
641 
642  // Get a reference to the LinearImplicitSystem we are solving
643  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Laplace");
644 
645  // A reference to the DofMap object for this system. The DofMap
646  // object handles the index translation from node and element numbers
647  // to degree of freedom numbers. We will talk more about the DofMap
648  // in future examples.
649  const DofMap & dof_map = system.get_dof_map();
650 
651  // Get a constant reference to the Finite Element type
652  // for the first (and only) variable in the system.
653  FEType fe_type = dof_map.variable_type(0);
654 
655  // Build a Finite Element object of the specified type. Since the
656  // FEBase::build() member dynamically creates memory we will
657  // store the object as a UniquePtr<FEBase>. This can be thought
658  // of as a pointer that will clean up after itself.
659  UniquePtr<FEBase> fe (FEBase::build(mesh_dim, fe_type));
660  UniquePtr<FEBase> fe_face (FEBase::build(mesh_dim, fe_type));
661 
662  // Quadrature rules for numerical integration.
663  UniquePtr<QBase> qrule(fe_type.default_quadrature_rule(mesh_dim));
664  UniquePtr<QBase> qface(fe_type.default_quadrature_rule(mesh_dim-1));
665 
666  // Tell the finite element object to use our quadrature rule.
667  fe->attach_quadrature_rule (qrule.get());
668  fe_face->attach_quadrature_rule (qface.get());
669 
670  // Here we define some references to cell-specific data that
671  // will be used to assemble the linear system.
672  // We begin with the element Jacobian * quadrature weight at each
673  // integration point.
674  const std::vector<Real> & JxW = fe->get_JxW();
675  const std::vector<Real> & JxW_face = fe_face->get_JxW();
676 
677  // The physical XY locations of the quadrature points on the element.
678  // These might be useful for evaluating spatially varying material
679  // properties or forcing functions at the quadrature points.
680  const std::vector<Point> & q_point = fe->get_xyz();
681 
682  // The element shape functions evaluated at the quadrature points.
683  // For this simple problem we usually only need them on element
684  // boundaries.
685  const std::vector<std::vector<Real>> & phi = fe->get_phi();
686  const std::vector<std::vector<Real>> & psi = fe_face->get_phi();
687 
688  // The element shape function gradients evaluated at the quadrature
689  // points.
690  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
691 
692  // The XY locations of the quadrature points used for face integration
693  const std::vector<Point> & qface_points = fe_face->get_xyz();
694 
695  // Define data structures to contain the element matrix
696  // and right-hand-side vector contribution. Following
697  // basic finite element terminology we will denote these
698  // "Ke" and "Fe". More detail is in example 3.
701 
702  // This vector will hold the degree of freedom indices for
703  // the element. These define where in the global system
704  // the element degrees of freedom get mapped.
705  std::vector<dof_id_type> dof_indices;
706 
707  // Now we will loop over all the elements in the mesh. We will
708  // compute the element matrix and right-hand-side contribution. See
709  // example 3 for a discussion of the element iterators. Here we use
710  // the const_active_local_elem_iterator to indicate we only want
711  // to loop over elements that are assigned to the local processor
712  // which are "active" in the sense of AMR. This allows each
713  // processor to compute its components of the global matrix for
714  // active elements while ignoring parent elements which have been
715  // refined.
716  for (const auto & elem : mesh.active_local_element_ptr_range())
717  {
718  // Start logging the shape function initialization.
719  // This is done through a simple function call with
720  // the name of the event to log.
721  perf_log.push("elem init");
722 
723  // Get the degree of freedom indices for the
724  // current element. These define where in the global
725  // matrix and right-hand-side this element will
726  // contribute to.
727  dof_map.dof_indices (elem, dof_indices);
728 
729  // Compute the element-specific data for the current
730  // element. This involves computing the location of the
731  // quadrature points (q_point) and the shape functions
732  // (phi, dphi) for the current element.
733  fe->reinit (elem);
734 
735  // Zero the element matrix and right-hand side before
736  // summing them. We use the resize member here because
737  // the number of degrees of freedom might have changed from
738  // the last element. Note that this will be the case if the
739  // element type is different (i.e. the last element was a
740  // triangle, now we are on a quadrilateral).
741  Ke.resize (dof_indices.size(),
742  dof_indices.size());
743 
744  Fe.resize (dof_indices.size());
745 
746  // Stop logging the shape function initialization.
747  // If you forget to stop logging an event the PerfLog
748  // object will probably catch the error and abort.
749  perf_log.pop("elem init");
750 
751  // Now we will build the element matrix. This involves
752  // a double loop to integrate the test functions (i) against
753  // the trial functions (j).
754  //
755  // Now start logging the element matrix computation
756  perf_log.push ("Ke");
757 
758  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
759  for (std::size_t i=0; i<dphi.size(); i++)
760  for (std::size_t j=0; j<dphi.size(); j++)
761  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
762 
763  // We need a forcing function to make the 1D case interesting
764  if (mesh_dim == 1)
765  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
766  {
767  Real x = q_point[qp](0);
768  Real f = singularity ? sqrt(3.)/9.*pow(-x, -4./3.) :
769  cos(x);
770  for (std::size_t i=0; i<dphi.size(); ++i)
771  Fe(i) += JxW[qp]*phi[i][qp]*f;
772  }
773 
774  // Stop logging the matrix computation
775  perf_log.pop ("Ke");
776 
777 
778  // At this point the interior element integration has
779  // been completed. However, we have not yet addressed
780  // boundary conditions. For this example we will only
781  // consider simple Dirichlet boundary conditions imposed
782  // via the penalty method.
783  //
784  // This approach adds the L2 projection of the boundary
785  // data in penalty form to the weak statement. This is
786  // a more generic approach for applying Dirichlet BCs
787  // which is applicable to non-Lagrange finite element
788  // discretizations.
789  {
790  // Start logging the boundary condition computation
791  perf_log.push ("BCs");
792 
793  // The penalty value.
794  const Real penalty = 1.e10;
795 
796  // The following loops over the sides of the element.
797  // If the element has no neighbor on a side then that
798  // side MUST live on a boundary of the domain.
799  for (auto s : elem->side_index_range())
800  if (elem->neighbor_ptr(s) == libmesh_nullptr)
801  {
802  fe_face->reinit(elem, s);
803 
804  for (unsigned int qp=0; qp<qface->n_points(); qp++)
805  {
806  const Number value = exact_solution (qface_points[qp],
807  es.parameters,
808  "null",
809  "void");
810 
811  // RHS contribution
812  for (std::size_t i=0; i<psi.size(); i++)
813  Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
814 
815  // Matrix contribution
816  for (std::size_t i=0; i<psi.size(); i++)
817  for (std::size_t j=0; j<psi.size(); j++)
818  Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
819  }
820  }
821 
822  // Stop logging the boundary condition computation
823  perf_log.pop ("BCs");
824  }
825 
826 
827  // The element matrix and right-hand-side are now built
828  // for this element. Add them to the global matrix and
829  // right-hand-side vector. The SparseMatrix::add_matrix()
830  // and NumericVector::add_vector() members do this for us.
831  // Start logging the insertion of the local (element)
832  // matrix and vector into the global matrix and vector
833  perf_log.push ("matrix insertion");
834 
835  dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
836  system.matrix->add_matrix (Ke, dof_indices);
837  system.rhs->add_vector (Fe, dof_indices);
838 
839  // Start logging the insertion of the local (element)
840  // matrix and vector into the global matrix and vector
841  perf_log.pop ("matrix insertion");
842  }
843 
844  // That's it. We don't need to do anything else to the
845  // PerfLog. When it goes out of scope (at this function return)
846  // it will print its log to the screen. Pretty easy, huh?
847 #endif // #ifdef LIBMESH_ENABLE_AMR
848 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
unsigned int n_linear_iterations() const
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
This is the EquationSystems class.
void assemble_laplace(EquationSystems &es, const std::string &system_name)
void pop(const char *label, const char *header="")
Pop the event label off the stack, resuming any lower event.
Definition: perf_log.C:162
virtual void select_refinement(System &system) libmesh_override
This pure virtual function must be redefined in derived classes to take a mesh flagged for h refineme...
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 does uniform refinements and a solve to get an improved solution on each cell...
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 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
void extra_quadrature_order(const int extraorder)
Increases or decreases the order of the quadrature rule used for numerical integration.
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 ...
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.
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...
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
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
std::list< Point > singular_points
This list, to be filled by the user, should include all singular points in the solution.
Definition: hp_singular.h:79
void add_algebraic_ghosting_functor(GhostingFunctor &ghosting_functor)
Adds a functor which can specify algebraic ghosting requirements for use with distributed vectors...
Definition: dof_map.C:1813
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:125
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.
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.
int main(int argc, char **argv)
virtual void all_second_order(const bool full_ordered=true)=0
Converts a (conforming, non-refined) mesh with linear elements into a mesh with second-order elements...
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
This class implements a ``brute force&#39;&#39; error estimator which integrates differences between the curr...
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
void push(const char *label, const char *header="")
Push the event label onto the stack, pausing any active event.
Definition: perf_log.C:132
Number exact_solution(const Point &p, const Parameters &, const std::string &, const std::string &)
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
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 Patch Recovery error estimate to estimate the error on each cell...
void plot_error(const std::string &filename, const MeshBase &mesh) const
Plots a data file, of a type determined by looking at the file extension in filename, of the error values on the active elements of mesh.
Definition: error_vector.C:208
double pow(double a, int b)
virtual void select_refinement(System &system)
This pure virtual function must be redefined in derived classes to take a mesh flagged for h refineme...
Definition: hp_singular.C:36
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...
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
Gradient exact_derivative(const Point &p, const Parameters &, const std::string &, const std::string &)
void libmesh_ignore(const T &)
This class implements the Kelly error indicator which is based on the flux jumps between elements...
This class uses the error estimate given by different types of derefinement (h coarsening or p reduct...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This class adds coupling (for use in send_list construction) between active elements and all descenda...
This class implements an "error estimator" based on the difference between the approximate and exact ...
Real h1_error(const std::string &sys_name, const std::string &unknown_name)
This class implements the Patch Recovery error indicator.
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
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.
static const bool value
Definition: xdr_io.C:108
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
Parameters parameters
Data structure holding arbitrary parameters.
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 exact solution function to estimate the error on each cell.
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 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 ...
This class uses a user-provided list of singularity locations to choose between h refining and p elev...
Definition: hp_singular.h:49
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...
virtual Real mean() const libmesh_override
Definition: error_vector.C:67
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
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.
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
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.