libMesh
adaptivity_ex4.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 // <h1>Adaptivity Example 4 - Biharmonic Equation</h1>
19 // \author Benjamin S. Kirk
20 // \date 2004
21 //
22 // This example solves the Biharmonic equation on a square or cube,
23 // using a Galerkin formulation with C1 elements approximating the
24 // H^2_0 function space.
25 // The initial mesh contains two TRI6, one QUAD9 or one HEX27
26 // An input file named "ex15.in"
27 // is provided which allows the user to set several parameters for
28 // the solution so that the problem can be re-run without a
29 // re-compile. The solution technique employed is to have a
30 // refinement loop with a linear solve inside followed by a
31 // refinement of the grid and projection of the solution to the new grid
32 // In the final loop iteration, there is no additional
33 // refinement after the solve. In the input file "ex15.in", the variable
34 // "max_r_steps" controls the number of refinement steps, and
35 // "max_r_level" controls the maximum element refinement level.
36 
37 // LibMesh include files.
38 #include "libmesh/mesh.h"
39 #include "libmesh/equation_systems.h"
40 #include "libmesh/linear_implicit_system.h"
41 #include "libmesh/exodusII_io.h"
42 #include "libmesh/fe.h"
43 #include "libmesh/quadrature.h"
44 #include "libmesh/dense_matrix.h"
45 #include "libmesh/dense_vector.h"
46 #include "libmesh/sparse_matrix.h"
47 #include "libmesh/mesh_generation.h"
48 #include "libmesh/mesh_modification.h"
49 #include "libmesh/mesh_refinement.h"
50 #include "libmesh/error_vector.h"
51 #include "libmesh/fourth_error_estimators.h"
52 #include "libmesh/getpot.h"
53 #include "libmesh/exact_solution.h"
54 #include "libmesh/dof_map.h"
55 #include "libmesh/numeric_vector.h"
56 #include "libmesh/elem.h"
57 #include "libmesh/tensor_value.h"
58 #include "libmesh/perf_log.h"
59 #include "libmesh/string_to_enum.h"
60 
61 // Bring in everything from the libMesh namespace
62 using namespace libMesh;
63 
64 // Function prototype. This is the function that will assemble
65 // the linear system for our Biharmonic problem. Note that the
66 // function will take the EquationSystems object and the
67 // name of the system we are assembling as input. From the
68 // EquationSystems object we have access to the Mesh and
69 // other objects we might need.
71  const std::string & system_name);
72 
73 
74 // Prototypes for calculation of the exact solution. Necessary
75 // for setting boundary conditions.
77  const Parameters &,
78  const std::string &,
79  const std::string &);
80 
82  const Parameters &, // parameters, not needed
83  const std::string &, // sys_name, not needed
84  const std::string &); // unk_name, not needed);
85 
87  const Parameters &,
88  const std::string &,
89  const std::string &);
90 
91 // Prototypes for calculation of the gradient of the exact solution.
92 // Necessary for setting boundary conditions in H^2_0 and testing
93 // H^1 convergence of the solution
95  const Parameters &,
96  const std::string &,
97  const std::string &);
98 
100  const Parameters &,
101  const std::string &,
102  const std::string &);
103 
105  const Parameters &,
106  const std::string &,
107  const std::string &);
108 
109 Tensor exact_1D_hessian(const Point & p,
110  const Parameters &,
111  const std::string &,
112  const std::string &);
113 
114 Tensor exact_2D_hessian(const Point & p,
115  const Parameters &,
116  const std::string &,
117  const std::string &);
118 
119 Tensor exact_3D_hessian(const Point & p,
120  const Parameters &,
121  const std::string &,
122  const std::string &);
123 
124 Number forcing_function_1D(const Point & p);
125 
126 Number forcing_function_2D(const Point & p);
127 
128 Number forcing_function_3D(const Point & p);
129 
130 // Pointers to dimension-independent functions
132  const Parameters &,
133  const std::string &,
134  const std::string &);
135 
137  const Parameters &,
138  const std::string &,
139  const std::string &);
140 
142  const Parameters &,
143  const std::string &,
144  const std::string &);
145 
147 
148 
149 
150 int main(int argc, char ** argv)
151 {
152  // Initialize libMesh.
153  LibMeshInit init (argc, argv);
154 
155  // This example requires a linear solver package.
156  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
157  "--enable-petsc, --enable-trilinos, or --enable-eigen");
158 
159  // Adaptive constraint calculations for fine Hermite elements seems
160  // to require half-decent precision
161 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
162  libmesh_example_requires(false, "double precision");
163 #endif
164 
165  // This example requires Adaptive Mesh Refinement support
166 #ifndef LIBMESH_ENABLE_AMR
167  libmesh_example_requires(false, "--enable-amr");
168 #else
169 
170  // This example requires second derivative calculation support
171 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
172  libmesh_example_requires(false, "--enable-second");
173 #else
174 
175  // Parse the input file
176  GetPot input_file("adaptivity_ex4.in");
177 
178  // But allow the command line to override it.
179  input_file.parse_command_line(argc, argv);
180 
181  // Read in parameters from the input file
182  const unsigned int max_r_level = input_file("max_r_level", 10);
183  const unsigned int max_r_steps = input_file("max_r_steps", 4);
184  const std::string approx_type = input_file("approx_type", "HERMITE");
185  const std::string approx_order_string = input_file("approx_order", "THIRD");
186  const unsigned int uniform_refine = input_file("uniform_refine", 0);
187  const Real refine_percentage = input_file("refine_percentage", 0.5);
188  const Real coarsen_percentage = input_file("coarsen_percentage", 0.5);
189  const unsigned int dim = input_file("dimension", 2);
190  const unsigned int max_linear_iterations = input_file("max_linear_iterations", 10000);
191 
192  // Skip higher-dimensional examples on a lower-dimensional libMesh build
193  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
194 
195  // Currently only the Hermite cubics give a 1D or 3D C^1 basis
196  libmesh_assert (dim == 2 || approx_type == "HERMITE");
197 
198  // Create a mesh, with dimension to be overridden later, on the
199  // default MPI communicator.
200  Mesh mesh(init.comm());
201 
202  // Output file for plotting the error
203  std::string output_file = "";
204 
205  if (dim == 1)
206  output_file += "1D_";
207  else if (dim == 2)
208  output_file += "2D_";
209  else if (dim == 3)
210  output_file += "3D_";
211 
212  if (approx_type == "HERMITE")
213  output_file += "hermite_";
214  else if (approx_type == "SECOND")
215  output_file += "reducedclough_";
216  else
217  output_file += "clough_";
218 
219  if (uniform_refine == 0)
220  output_file += "adaptive";
221  else
222  output_file += "uniform";
223 
224 #ifdef LIBMESH_HAVE_EXODUS_API
225  // If we have Exodus, use the same base output filename
226  std::string exd_file = output_file;
227  exd_file += ".e";
228 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
229 
230  output_file += ".m";
231 
232  std::ofstream out (output_file.c_str());
233  out << "% dofs L2-error H1-error H2-error\n"
234  << "e = [\n";
235 
236  // Set up the dimension-dependent coarse mesh and solution
237  if (dim == 1)
238  {
244  }
245 
246  if (dim == 2)
247  {
253  }
254  else if (dim == 3)
255  {
261  }
262 
263  // Convert the mesh to second order: necessary for computing with
264  // Clough-Tocher elements, useful for getting slightly less
265  // broken visualization output with Hermite elements
267 
268  // Convert it to triangles if necessary
269  if (approx_type != "HERMITE")
271 
272  // Mesh Refinement object
273  MeshRefinement mesh_refinement(mesh);
274  mesh_refinement.refine_fraction() = refine_percentage;
275  mesh_refinement.coarsen_fraction() = coarsen_percentage;
276  mesh_refinement.max_h_level() = max_r_level;
277 
278  // Create an equation systems object.
279  EquationSystems equation_systems (mesh);
280 
281  // Declare the system and its variables.
282  // Create a system named "Biharmonic"
283  LinearImplicitSystem & system =
284  equation_systems.add_system<LinearImplicitSystem> ("Biharmonic");
285 
286  Order approx_order = approx_type == "SECOND" ? SECOND :
287  Utility::string_to_enum<Order>(approx_order_string);
288 
289  // Adds the variable "u" to "Biharmonic". "u" will be approximated
290  // using Hermite tensor product squares or Clough-Tocher triangles
291 
292  if (approx_type == "HERMITE")
293  system.add_variable("u", approx_order, HERMITE);
294  else if (approx_type == "SECOND")
295  system.add_variable("u", SECOND, CLOUGH);
296  else if (approx_type == "CLOUGH")
297  system.add_variable("u", approx_order, CLOUGH);
298  else
299  libmesh_error_msg("Invalid approx_type = " << approx_type);
300 
301  // Give the system a pointer to the matrix assembly
302  // function.
304 
305  // Initialize the data structures for the equation system.
306  equation_systems.init();
307 
308  // Set linear solver max iterations
309  equation_systems.parameters.set<unsigned int>
310  ("linear solver maximum iterations") = max_linear_iterations;
311 
312  // Linear solver tolerance.
313  equation_systems.parameters.set<Real>
314  ("linear solver tolerance") = TOLERANCE * TOLERANCE;
315 
316  // Prints information about the system to the screen.
317  equation_systems.print_info();
318 
319  // Construct ExactSolution object and attach function to compute exact solution
320  ExactSolution exact_sol(equation_systems);
324 
325  // Construct zero solution object, useful for computing solution norms
326  // Attaching "zero_solution" functions is unnecessary
327  ExactSolution zero_sol(equation_systems);
328 
329  // A refinement loop.
330  for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
331  {
332  mesh.print_info();
333  equation_systems.print_info();
334 
335  libMesh::out << "Beginning Solve " << r_step << std::endl;
336 
337  // Solve the system "Biharmonic", just like example 2.
338  system.solve();
339 
340  libMesh::out << "Linear solver converged at step: "
341  << system.n_linear_iterations()
342  << ", final residual: "
343  << system.final_linear_residual()
344  << std::endl;
345 
346  // Compute the error.
347  exact_sol.compute_error("Biharmonic", "u");
348  // Compute the norm.
349  zero_sol.compute_error("Biharmonic", "u");
350 
351  // Print out the error values
352  libMesh::out << "L2-Norm is: "
353  << zero_sol.l2_error("Biharmonic", "u")
354  << std::endl;
355  libMesh::out << "H1-Norm is: "
356  << zero_sol.h1_error("Biharmonic", "u")
357  << std::endl;
358  libMesh::out << "H2-Norm is: "
359  << zero_sol.h2_error("Biharmonic", "u")
360  << std::endl
361  << std::endl;
362  libMesh::out << "L2-Error is: "
363  << exact_sol.l2_error("Biharmonic", "u")
364  << std::endl;
365  libMesh::out << "H1-Error is: "
366  << exact_sol.h1_error("Biharmonic", "u")
367  << std::endl;
368  libMesh::out << "H2-Error is: "
369  << exact_sol.h2_error("Biharmonic", "u")
370  << std::endl
371  << std::endl;
372 
373  // Print to output file
374  out << equation_systems.n_active_dofs() << " "
375  << exact_sol.l2_error("Biharmonic", "u") << " "
376  << exact_sol.h1_error("Biharmonic", "u") << " "
377  << exact_sol.h2_error("Biharmonic", "u") << std::endl;
378 
379  // Possibly refine the mesh
380  if (r_step+1 != max_r_steps)
381  {
382  libMesh::out << " Refining the mesh..." << std::endl;
383 
384  if (uniform_refine == 0)
385  {
386  ErrorVector error;
387  LaplacianErrorEstimator error_estimator;
388 
389  error_estimator.estimate_error(system, error);
390  mesh_refinement.flag_elements_by_elem_fraction (error);
391 
392  libMesh::out << "Mean Error: " << error.mean() << std::endl;
393  libMesh::out << "Error Variance: " << error.variance() << std::endl;
394 
395  mesh_refinement.refine_and_coarsen_elements();
396  }
397  else
398  {
399  mesh_refinement.uniformly_refine(1);
400  }
401 
402  // This call reinitializes the EquationSystems object for
403  // the newly refined mesh. One of the steps in the
404  // reinitialization is projecting the solution,
405  // old_solution, etc... vectors from the old mesh to
406  // the current one.
407  equation_systems.reinit ();
408  }
409  }
410 
411 #ifdef LIBMESH_HAVE_EXODUS_API
412  // After solving the system write the solution
413  // to a ExodusII-formatted plot file.
415  equation_systems);
416 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
417 
418  // Close up the output file.
419  out << "];\n"
420  << "hold on\n"
421  << "loglog(e(:,1), e(:,2), 'bo-');\n"
422  << "loglog(e(:,1), e(:,3), 'ro-');\n"
423  << "loglog(e(:,1), e(:,4), 'go-');\n"
424  << "xlabel('log(dofs)');\n"
425  << "ylabel('log(error)');\n"
426  << "title('C1 " << approx_type << " elements');\n"
427  << "legend('L2-error', 'H1-error', 'H2-error');\n";
428 
429  // All done.
430  return 0;
431 #endif // #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
432 #endif // #ifndef LIBMESH_ENABLE_AMR
433 }
434 
435 
436 
438  const Parameters &, // parameters, not needed
439  const std::string &, // sys_name, not needed
440  const std::string &) // unk_name, not needed
441 {
442  // x coordinate in space
443  const Real x = p(0);
444 
445  // analytic solution value
446  return 256.*(x-x*x)*(x-x*x);
447 }
448 
449 
450 // We now define the gradient of the exact solution
452  const Parameters &, // parameters, not needed
453  const std::string &, // sys_name, not needed
454  const std::string &) // unk_name, not needed
455 {
456  // x coordinate in space
457  const Real x = p(0);
458 
459  // First derivatives to be returned.
460  Gradient gradu;
461 
462  gradu(0) = 256.*2.*(x-x*x)*(1-2*x);
463 
464  return gradu;
465 }
466 
467 
468 // We now define the hessian of the exact solution
470  const Parameters &, // parameters, not needed
471  const std::string &, // sys_name, not needed
472  const std::string &) // unk_name, not needed
473 {
474  // Second derivatives to be returned.
475  Tensor hessu;
476 
477  // x coordinate in space
478  const Real x = p(0);
479 
480  hessu(0,0) = 256.*2.*(1-6.*x+6.*x*x);
481 
482  return hessu;
483 }
484 
485 
486 
488 {
489  // Equals laplacian(laplacian(u)), u'''' in 1D
490  return 256. * 2. * 12.;
491 }
492 
493 
495  const Parameters &, // parameters, not needed
496  const std::string &, // sys_name, not needed
497  const std::string &) // unk_name, not needed
498 {
499  // x and y coordinates in space
500  const Real x = p(0);
501  const Real y = p(1);
502 
503  // analytic solution value
504  return 256.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y);
505 }
506 
507 
508 // We now define the gradient of the exact solution
510  const Parameters &, // parameters, not needed
511  const std::string &, // sys_name, not needed
512  const std::string &) // unk_name, not needed
513 {
514  // x and y coordinates in space
515  const Real x = p(0);
516  const Real y = p(1);
517 
518  // First derivatives to be returned.
519  Gradient gradu;
520 
521  gradu(0) = 256.*2.*(x-x*x)*(1-2*x)*(y-y*y)*(y-y*y);
522  gradu(1) = 256.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1-2*y);
523 
524  return gradu;
525 }
526 
527 
528 // We now define the hessian of the exact solution
530  const Parameters &, // parameters, not needed
531  const std::string &, // sys_name, not needed
532  const std::string &) // unk_name, not needed
533 {
534  // Second derivatives to be returned.
535  Tensor hessu;
536 
537  // x and y coordinates in space
538  const Real x = p(0);
539  const Real y = p(1);
540 
541  hessu(0,0) = 256.*2.*(1-6.*x+6.*x*x)*(y-y*y)*(y-y*y);
542  hessu(0,1) = 256.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y);
543  hessu(1,1) = 256.*2.*(x-x*x)*(x-x*x)*(1.-6.*y+6.*y*y);
544 
545  // Hessians are always symmetric
546  hessu(1,0) = hessu(0,1);
547  return hessu;
548 }
549 
550 
551 
553 {
554  // x and y coordinates in space
555  const Real x = p(0);
556  const Real y = p(1);
557 
558  // Equals laplacian(laplacian(u))
559  return 256. * 8. * (3.*((y-y*y)*(y-y*y)+(x-x*x)*(x-x*x))
560  + (1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y));
561 }
562 
563 
564 
566  const Parameters &, // parameters, not needed
567  const std::string &, // sys_name, not needed
568  const std::string &) // unk_name, not needed
569 {
570  // xyz coordinates in space
571  const Real x = p(0);
572  const Real y = p(1);
573  const Real z = p(2);
574 
575  // analytic solution value
576  return 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
577 }
578 
579 
581  const Parameters &, // parameters, not needed
582  const std::string &, // sys_name, not needed
583  const std::string &) // unk_name, not needed
584 {
585  // First derivatives to be returned.
586  Gradient gradu;
587 
588  // xyz coordinates in space
589  const Real x = p(0);
590  const Real y = p(1);
591  const Real z = p(2);
592 
593  gradu(0) = 4096.*2.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
594  gradu(1) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
595  gradu(2) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);
596 
597  return gradu;
598 }
599 
600 
601 // We now define the hessian of the exact solution
603  const Parameters &, // parameters, not needed
604  const std::string &, // sys_name, not needed
605  const std::string &) // unk_name, not needed
606 {
607  // Second derivatives to be returned.
608  Tensor hessu;
609 
610  // xyz coordinates in space
611  const Real x = p(0);
612  const Real y = p(1);
613  const Real z = p(2);
614 
615  hessu(0,0) = 4096.*(2.-12.*x+12.*x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
616  hessu(0,1) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
617  hessu(0,2) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);
618  hessu(1,1) = 4096.*(x-x*x)*(x-x*x)*(2.-12.*y+12.*y*y)*(z-z*z)*(z-z*z);
619  hessu(1,2) = 4096.*4.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(1.-2.*z);
620  hessu(2,2) = 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(2.-12.*z+12.*z*z);
621 
622  // Hessians are always symmetric
623  hessu(1,0) = hessu(0,1);
624  hessu(2,0) = hessu(0,2);
625  hessu(2,1) = hessu(1,2);
626 
627  return hessu;
628 }
629 
630 
631 
633 {
634  // xyz coordinates in space
635  const Real x = p(0);
636  const Real y = p(1);
637  const Real z = p(2);
638 
639  // Equals laplacian(laplacian(u))
640  return 4096. * 8. * (3.*((y-y*y)*(y-y*y)*(x-x*x)*(x-x*x) +
641  (z-z*z)*(z-z*z)*(x-x*x)*(x-x*x) +
642  (z-z*z)*(z-z*z)*(y-y*y)*(y-y*y)) +
643  (1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y)*(z-z*z)*(z-z*z) +
644  (1.-6.*x+6.*x*x)*(1.-6.*z+6.*z*z)*(y-y*y)*(y-y*y) +
645  (1.-6.*y+6.*y*y)*(1.-6.*z+6.*z*z)*(x-x*x)*(x-x*x));
646 }
647 
648 
649 
650 // We now define the matrix assembly function for the
651 // Biharmonic system. We need to first compute element
652 // matrices and right-hand sides, and then take into
653 // account the boundary conditions, which will be handled
654 // via a penalty method.
656  const std::string & system_name)
657 {
658  // Ignore unused parameter warnings when libmesh is configured without certain options.
659  libmesh_ignore(es);
660  libmesh_ignore(system_name);
661 
662 #ifdef LIBMESH_ENABLE_AMR
663 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
664 
665  // It is a good idea to make sure we are assembling
666  // the proper system.
667  libmesh_assert_equal_to (system_name, "Biharmonic");
668 
669  // Declare a performance log. Give it a descriptive
670  // string to identify what part of the code we are
671  // logging, since there may be many PerfLogs in an
672  // application.
673  PerfLog perf_log ("Matrix Assembly", false);
674 
675  // Get a constant reference to the mesh object.
676  const MeshBase & mesh = es.get_mesh();
677 
678  // The dimension that we are running
679  const unsigned int dim = mesh.mesh_dimension();
680 
681  // Get a reference to the LinearImplicitSystem we are solving
682  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Biharmonic");
683 
684  // A reference to the DofMap object for this system. The DofMap
685  // object handles the index translation from node and element numbers
686  // to degree of freedom numbers. We will talk more about the DofMap
687  // in future examples.
688  const DofMap & dof_map = system.get_dof_map();
689 
690  // Get a constant reference to the Finite Element type
691  // for the first (and only) variable in the system.
692  FEType fe_type = dof_map.variable_type(0);
693 
694  // Build a Finite Element object of the specified type. Since the
695  // FEBase::build() member dynamically creates memory we will
696  // store the object as a UniquePtr<FEBase>. This can be thought
697  // of as a pointer that will clean up after itself.
698  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
699 
700  // Quadrature rule for numerical integration.
701  // With 2D triangles, the Clough quadrature rule puts a Gaussian
702  // quadrature rule on each of the 3 subelements
703  UniquePtr<QBase> qrule(fe_type.default_quadrature_rule(dim));
704 
705  // Tell the finite element object to use our quadrature rule.
706  fe->attach_quadrature_rule (qrule.get());
707 
708  // Declare a special finite element object for
709  // boundary integration.
710  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
711 
712  // Boundary integration requires another quadrature rule,
713  // with dimensionality one less than the dimensionality
714  // of the element.
715  // In 1D, the Clough and Gauss quadrature rules are identical.
716  UniquePtr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
717 
718  // Tell the finite element object to use our
719  // quadrature rule.
720  fe_face->attach_quadrature_rule (qface.get());
721 
722  // Here we define some references to cell-specific data that
723  // will be used to assemble the linear system.
724  // We begin with the element Jacobian * quadrature weight at each
725  // integration point.
726  const std::vector<Real> & JxW = fe->get_JxW();
727 
728  // The physical XY locations of the quadrature points on the element.
729  // These might be useful for evaluating spatially varying material
730  // properties at the quadrature points.
731  const std::vector<Point> & q_point = fe->get_xyz();
732 
733  // The element shape functions evaluated at the quadrature points.
734  const std::vector<std::vector<Real>> & phi = fe->get_phi();
735 
736  // The element shape function second derivatives evaluated at the
737  // quadrature points. Note that for the simple biharmonic, shape
738  // function first derivatives are unnecessary.
739  const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
740 
741  // For efficiency we will compute shape function laplacians n times,
742  // not n^2
743  std::vector<Real> shape_laplacian;
744 
745  // Define data structures to contain the element matrix
746  // and right-hand-side vector contribution. Following
747  // basic finite element terminology we will denote these
748  // "Ke" and "Fe". More detail is in example 3.
751 
752  // This vector will hold the degree of freedom indices for
753  // the element. These define where in the global system
754  // the element degrees of freedom get mapped.
755  std::vector<dof_id_type> dof_indices;
756 
757  // Now we will loop over all the elements in the mesh. We will
758  // compute the element matrix and right-hand-side contribution. See
759  // example 3 for a discussion of the element iterators.
760  for (const auto & elem : mesh.active_local_element_ptr_range())
761  {
762  // Start logging the shape function initialization.
763  // This is done through a simple function call with
764  // the name of the event to log.
765  perf_log.push("elem init");
766 
767  // Get the degree of freedom indices for the
768  // current element. These define where in the global
769  // matrix and right-hand-side this element will
770  // contribute to.
771  dof_map.dof_indices (elem, dof_indices);
772 
773  // Compute the element-specific data for the current
774  // element. This involves computing the location of the
775  // quadrature points (q_point) and the shape functions
776  // (phi, dphi) for the current element.
777  fe->reinit (elem);
778 
779  // Zero the element matrix and right-hand side before
780  // summing them.
781  Ke.resize (dof_indices.size(),
782  dof_indices.size());
783 
784  Fe.resize (dof_indices.size());
785 
786  // Make sure there is enough room in this cache
787  shape_laplacian.resize(dof_indices.size());
788 
789  // Stop logging the shape function initialization.
790  // If you forget to stop logging an event the PerfLog
791  // object will probably catch the error and abort.
792  perf_log.pop("elem init");
793 
794  // Now we will build the element matrix. This involves
795  // a double loop to integrate laplacians of the test functions
796  // (i) against laplacians of the trial functions (j).
797  //
798  // This step is why we need the Clough-Tocher elements -
799  // these C1 differentiable elements have square-integrable
800  // second derivatives.
801  //
802  // Now start logging the element matrix computation
803  perf_log.push ("Ke");
804 
805  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
806  {
807  for (std::size_t i=0; i<phi.size(); i++)
808  {
809  shape_laplacian[i] = d2phi[i][qp](0,0);
810  if (dim > 1)
811  shape_laplacian[i] += d2phi[i][qp](1,1);
812  if (dim == 3)
813  shape_laplacian[i] += d2phi[i][qp](2,2);
814  }
815  for (std::size_t i=0; i<phi.size(); i++)
816  for (std::size_t j=0; j<phi.size(); j++)
817  Ke(i,j) += JxW[qp]*
818  shape_laplacian[i]*shape_laplacian[j];
819  }
820 
821  // Stop logging the matrix computation
822  perf_log.pop ("Ke");
823 
824 
825  // At this point the interior element integration has
826  // been completed. However, we have not yet addressed
827  // boundary conditions. For this example we will only
828  // consider simple Dirichlet boundary conditions imposed
829  // via the penalty method. Note that this is a fourth-order
830  // problem: Dirichlet boundary conditions include *both*
831  // boundary values and boundary normal fluxes.
832  {
833  // Start logging the boundary condition computation
834  perf_log.push ("BCs");
835 
836  // The penalty values, for solution boundary trace and flux.
837  const Real penalty = 1e10;
838  const Real penalty2 = 1e10;
839 
840  // The following loops over the sides of the element.
841  // If the element has no neighbor on a side then that
842  // side MUST live on a boundary of the domain.
843  for (auto s : elem->side_index_range())
844  if (elem->neighbor_ptr(s) == libmesh_nullptr)
845  {
846  // The value of the shape functions at the quadrature
847  // points.
848  const std::vector<std::vector<Real>> & phi_face =
849  fe_face->get_phi();
850 
851  // The value of the shape function derivatives at the
852  // quadrature points.
853  const std::vector<std::vector<RealGradient>> & dphi_face =
854  fe_face->get_dphi();
855 
856  // The Jacobian * Quadrature Weight at the quadrature
857  // points on the face.
858  const std::vector<Real> & JxW_face = fe_face->get_JxW();
859 
860  // The XYZ locations (in physical space) of the
861  // quadrature points on the face. This is where
862  // we will interpolate the boundary value function.
863  const std::vector<Point> & qface_point = fe_face->get_xyz();
864  const std::vector<Point> & face_normals = fe_face->get_normals();
865 
866  // Compute the shape function values on the element
867  // face.
868  fe_face->reinit(elem, s);
869 
870  // Loop over the face quadrature points for integration.
871  for (unsigned int qp=0; qp<qface->n_points(); qp++)
872  {
873  // The boundary value.
874  Number value = exact_solution(qface_point[qp],
875  es.parameters, "null",
876  "void");
877  Gradient flux = exact_2D_derivative(qface_point[qp],
878  es.parameters,
879  "null", "void");
880 
881  // Matrix contribution of the L2 projection.
882  // Note that the basis function values are
883  // integrated against test function values while
884  // basis fluxes are integrated against test function
885  // fluxes.
886  for (std::size_t i=0; i<phi_face.size(); i++)
887  for (std::size_t j=0; j<phi_face.size(); j++)
888  Ke(i,j) += JxW_face[qp] *
889  (penalty * phi_face[i][qp] *
890  phi_face[j][qp] + penalty2
891  * (dphi_face[i][qp] *
892  face_normals[qp]) *
893  (dphi_face[j][qp] *
894  face_normals[qp]));
895 
896  // Right-hand-side contribution of the L2
897  // projection.
898  for (std::size_t i=0; i<phi_face.size(); i++)
899  Fe(i) += JxW_face[qp] *
900  (penalty * value * phi_face[i][qp]
901  + penalty2 *
902  (flux * face_normals[qp])
903  * (dphi_face[i][qp]
904  * face_normals[qp]));
905 
906  }
907  }
908 
909  // Stop logging the boundary condition computation
910  perf_log.pop ("BCs");
911  }
912 
913  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
914  for (std::size_t i=0; i<phi.size(); i++)
915  Fe(i) += JxW[qp]*phi[i][qp]*forcing_function(q_point[qp]);
916 
917  // The element matrix and right-hand-side are now built
918  // for this element. Add them to the global matrix and
919  // right-hand-side vector. The SparseMatrix::add_matrix()
920  // and NumericVector::add_vector() members do this for us.
921  // Start logging the insertion of the local (element)
922  // matrix and vector into the global matrix and vector
923  perf_log.push ("matrix insertion");
924 
925  dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
926  system.matrix->add_matrix (Ke, dof_indices);
927  system.rhs->add_vector (Fe, dof_indices);
928 
929  // Stop logging the insertion of the local (element)
930  // matrix and vector into the global matrix and vector
931  perf_log.pop ("matrix insertion");
932  }
933 
934  // That's it. We don't need to do anything else to the
935  // PerfLog. When it goes out of scope (at this function return)
936  // it will print its log to the screen. Pretty easy, huh?
937 
938 #else
939 
940 #endif // #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
941 #endif // #ifdef LIBMESH_ENABLE_AMR
942 }
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...
virtual Real variance() const libmesh_override
Definition: error_vector.h:115
void assemble_biharmonic(EquationSystems &es, const std::string &system_name)
This is the EquationSystems class.
void pop(const char *label, const char *header="")
Pop the event label off the stack, resuming any lower event.
Definition: perf_log.C:162
Gradient(* exact_derivative)(const Point &p, const Parameters &, const std::string &, const std::string &)
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
Number(* exact_solution)(const Point &p, const Parameters &, const std::string &, const std::string &)
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.
Gradient exact_1D_derivative(const Point &p, const Parameters &, const std::string &, const std::string &)
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 h2_error(const std::string &sys_name, const std::string &unknown_name)
Gradient exact_2D_derivative(const Point &p, const Parameters &, const std::string &, const std::string &)
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.
Tensor exact_1D_hessian(const Point &p, const Parameters &, const std::string &, const std::string &)
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
static const Real TOLERANCE
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
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.
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 is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
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
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
Number exact_1D_solution(const Point &p, const Parameters &, const std::string &, const std::string &)
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.
Tensor exact_3D_hessian(const Point &p, const Parameters &, const std::string &, const std::string &)
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 print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
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.
NumberVectorValue Gradient
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.
Tensor(* exact_hessian)(const Point &p, const Parameters &, const std::string &, const std::string &)
void push(const char *label, const char *header="")
Push the event label onto the stack, pausing any active event.
Definition: perf_log.C:132
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
Tensor exact_2D_hessian(const Point &p, const Parameters &, const std::string &, const std::string &)
Gradient exact_3D_derivative(const Point &p, const Parameters &, const std::string &, const std::string &)
Number forcing_function_2D(const Point &p)
Number(* forcing_function)(const Point &p)
void attach_exact_hessian(unsigned int sys_num, FunctionBase< Tensor > *h)
Clone and attach an arbitrary functor which computes the exact second derivatives of the system sys_n...
void flag_elements_by_elem_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 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
void libmesh_ignore(const T &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
NumberTensorValue Tensor
Real h1_error(const std::string &sys_name, const std::string &unknown_name)
T & set(const std::string &)
Definition: parameters.h:469
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
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
Number forcing_function_1D(const Point &p)
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
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.
int main(int argc, char **argv)
Number exact_2D_solution(const Point &p, const Parameters &, const std::string &, const std::string &)
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...
Number exact_3D_solution(const Point &p, const Parameters &, const std::string &, const std::string &)
Number forcing_function_3D(const Point &p)
virtual void init()
Initialize all the systems.
virtual Real mean() const libmesh_override
Definition: error_vector.C:67
This class is an error indicator based on laplacian jumps between elements.
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
std::size_t n_active_dofs() const
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
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
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
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.