libMesh
adaptivity_ex5.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 5 - Periodic Boundary Conditions with Adaptive Mesh Refinement</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This example uses the same simple, linear transient
25 // system as in example 10; however in this case periodic boundary
26 // conditions are applied at the sides of the domain.
27 //
28 // This code also contains an example use of ParsedFunction, to
29 // allow users to specify an exact solution on the command line.
30 
31 // C++ include files that we need
32 #include <iostream>
33 #include <algorithm>
34 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
35 #include <cmath>
36 #include <memory>
37 
38 // libMesh includes
39 #include "libmesh/libmesh.h"
40 #include "libmesh/replicated_mesh.h"
41 #include "libmesh/mesh_refinement.h"
42 #include "libmesh/gmv_io.h"
43 #include "libmesh/exodusII_io.h"
44 #include "libmesh/equation_systems.h"
45 #include "libmesh/fe.h"
46 #include "libmesh/quadrature_gauss.h"
47 #include "libmesh/dof_map.h"
48 #include "libmesh/sparse_matrix.h"
49 #include "libmesh/numeric_vector.h"
50 #include "libmesh/dense_matrix.h"
51 #include "libmesh/dense_vector.h"
52 #include "libmesh/periodic_boundaries.h"
53 #include "libmesh/periodic_boundary.h"
54 #include "libmesh/mesh_generation.h"
55 #include "libmesh/parsed_function.h"
56 #include "libmesh/getpot.h"
57 #include "libmesh/enum_solver_package.h"
58 #include "libmesh/enum_xdr_mode.h"
59 #include "libmesh/enum_norm_type.h"
60 
61 // This example will solve a linear transient system,
62 // so we need to include the TransientLinearImplicitSystem definition.
63 #include "libmesh/transient_system.h"
64 #include "libmesh/linear_implicit_system.h"
65 #include "libmesh/vector_value.h"
66 
67 // To refine the mesh we need an ErrorEstimator
68 // object to figure out which elements to refine.
69 #include "libmesh/error_vector.h"
70 #include "libmesh/kelly_error_estimator.h"
71 
72 // The definition of a geometric element
73 #include "libmesh/elem.h"
74 
75 // Bring in everything from the libMesh namespace
76 using namespace libMesh;
77 
78 // Function prototype. This function will assemble the system
79 // matrix and right-hand-side at each time step. Note that
80 // since the system is linear we technically do not need to
81 // assemble the matrix at each time step, but we will anyway.
82 // In subsequent examples we will employ adaptive mesh refinement,
83 // and with a changing mesh it will be necessary to rebuild the
84 // system matrix.
85 #ifdef LIBMESH_ENABLE_AMR
86 void assemble_cd (EquationSystems & es,
87  const std::string & system_name);
88 #endif
89 
90 // Function prototype. This function will initialize the system.
91 // Initialization functions are optional for systems. They allow
92 // you to specify the initial values of the solution. If an
93 // initialization function is not provided then the default (0)
94 // solution is provided.
95 void init_cd (EquationSystems & es,
96  const std::string & system_name);
97 
98 // Exact solution function prototype. This gives the exact
99 // solution as a function of space and time. In this case the
100 // initial condition will be taken as the exact solution at time 0,
101 // as will the Dirichlet boundary conditions at time t.
102 Real exact_solution (const Real x,
103  const Real y,
104  const Real t);
105 
107  const Parameters & parameters,
108  const std::string &,
109  const std::string &)
110 {
111  return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
112 }
113 
114 // With --enable-fparser, the user can also optionally set their own
115 // exact solution equations.
116 std::unique_ptr<FunctionBase<Number>> parsed_solution;
117 
118 
119 // Returns a string with 'number' formatted and placed directly
120 // into the string in some way
121 std::string exodus_filename(unsigned number);
122 
123 
124 // Begin the main program. Note that the first
125 // statement in the program throws an error if
126 // you are in complex number mode, since this
127 // example is only intended to work with real
128 // numbers.
129 int main (int argc, char ** argv)
130 {
131  // Initialize libMesh.
132  LibMeshInit init (argc, argv);
133 
134  // This example requires a linear solver package.
135  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
136  "--enable-petsc, --enable-trilinos, or --enable-eigen");
137 
138  // Skip this 2D example if libMesh was compiled as 1D-only.
139  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
140 
141 #if !defined(LIBMESH_ENABLE_AMR)
142  libmesh_example_requires(false, "--enable-amr");
143 #elif !defined(LIBMESH_HAVE_XDR)
144  // We use XDR support in our output here
145  libmesh_example_requires(false, "--enable-xdr");
146 #elif !defined(LIBMESH_ENABLE_PERIODIC)
147  libmesh_example_requires(false, "--enable-periodic");
148 #elif (LIBMESH_DOF_ID_BYTES == 8)
149  libmesh_example_requires(false, "--with-dof-id-bytes=4");
150 #else
151 
152  // Our Trilinos interface does not yet support adaptive transient
153  // problems
154  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc");
155 
156  // Brief message to the user regarding the program name
157  // and command line arguments.
158 
159  // Use commandline parameter to specify if we are to
160  // read in an initial solution or generate it ourself
161  libMesh::out << "Usage:\n"
162  <<"\t " << argv[0] << " -init_timestep 0\n"
163  << "OR\n"
164  <<"\t " << argv[0] << " -read_solution -init_timestep 26\n"
165  << std::endl;
166 
167  libMesh::out << "Running: " << argv[0];
168 
169  for (int i=1; i<argc; i++)
170  libMesh::out << " " << argv[i];
171 
172  libMesh::out << std::endl << std::endl;
173 
174  // This boolean value is obtained from the command line, it is true
175  // if the flag "-read_solution" is present, false otherwise.
176  // It indicates whether we are going to read in
177  // the mesh and solution files "saved_mesh.xda" and "saved_solution.xda"
178  // or whether we are going to start from scratch by just reading
179  // "mesh.xda"
180  const bool read_solution = libMesh::on_command_line("-read_solution");
181 
182  // This value is also obtained from the commandline and it specifies the
183  // initial value for the t_step looping variable. We must
184  // distinguish between the two cases here, whether we read in the
185  // solution or we started from scratch, so that we do not overwrite the
186  // gmv output files.
187  const unsigned int init_timestep =
188  libMesh::command_line_next("-init_timestep",
190 
191  if (init_timestep == libMesh::invalid_uint)
192  {
193  // This handy function will print the file name, line number,
194  // specified message, and then throw an exception.
195  libmesh_error_msg("ERROR: Initial timestep not specified!");
196  }
197 
198  // This value is also obtained from the command line, and specifies
199  // the number of time steps to take.
200  const unsigned int n_timesteps =
201  libMesh::command_line_next("-n_timesteps",
203 
204  if (n_timesteps == libMesh::invalid_uint)
205  libmesh_error_msg("ERROR: Number of timesteps not specified");
206 
207  // The user can specify a different exact solution on the command
208  // line, if we have an expression parser compiled in
209 #ifdef LIBMESH_HAVE_FPARSER
210  const bool have_expression = libMesh::on_command_line("-exact_solution");
211 #else
212  const bool have_expression = false;
213 #endif
214  if (have_expression)
215  parsed_solution = std::make_unique<ParsedFunction<Number>>
216  (libMesh::command_line_next("-exact_solution", std::string()));
217 
218  // Skip this 2D example if libMesh was compiled as 1D-only.
219  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
220 
221  // Create a new mesh on the default MPI communicator.
222  // There seems to be a DistributedMesh bug triggered here.
223  ReplicatedMesh mesh(init.comm());
224 
225  // Create an equation systems object.
226  EquationSystems equation_systems (mesh);
227  MeshRefinement mesh_refinement (mesh);
228 
229  // Declare the system and its variables.
230  // Begin by creating a transient system
231  // named "Convection-Diffusion".
233  equation_systems.add_system<TransientLinearImplicitSystem>
234  ("Convection-Diffusion");
235 
236  // Give the system a pointer to the assembly function.
237  system.attach_assemble_function (assemble_cd);
238 
239  // Creating and attaching Periodic Boundaries
240  DofMap & dof_map = system.get_dof_map();
241 
242  // Create a boundary periodic with one displaced 2.0 in the x
243  // direction
244  PeriodicBoundary horz(RealVectorValue(2.0, 0., 0.));
245 
246  // Connect boundary ids 3 and 1 with it
247  horz.myboundary = 3;
248  horz.pairedboundary = 1;
249 
250  // Add it to the PeriodicBoundaries
251  dof_map.add_periodic_boundary(horz);
252 
253  // Create a boundary periodic with one displaced 2.0 in the y
254  // direction
255  PeriodicBoundary vert(RealVectorValue(0., 2.0, 0.));
256 
257  // Connect boundary ids 0 and 2 with it
258  vert.myboundary = 0;
259  vert.pairedboundary = 2;
260 
261  // Add it to the PeriodicBoundaries
262  dof_map.add_periodic_boundary(vert);
263 
264  // Next build or read the mesh. We do this only *after* generating
265  // periodic boundaries; otherwise a DistributedMesh won't know to
266  // retain periodic neighbor elements.
267 
268  // First we process the case where we do not read in the solution
269  if (!read_solution)
270  {
271  MeshTools::Generation::build_square(mesh, 2, 2, 0., 2., 0., 2., QUAD4);
272 
273  // Again do a search on the command line for an argument
274  const unsigned int n_refinements =
275  libMesh::command_line_next("-n_refinements", 5);
276 
277  // Uniformly refine the mesh 5 times
278  if (!read_solution)
279  mesh_refinement.uniformly_refine (n_refinements);
280 
281  // Print information about the mesh to the screen.
282  mesh.print_info();
283 
284  // Adds the variable "u" to "Convection-Diffusion". "u"
285  // will be approximated using first-order approximation.
286  system.add_variable ("u", FIRST);
287 
288  // Give the system a pointer to the initialization function.
289  system.attach_init_function (init_cd);
290  }
291  // Otherwise we read in the solution and mesh
292  else
293  {
294  // Read in the mesh stored in "saved_mesh.xda"
295  mesh.read("saved_mesh.xdr");
296 
297  // Print information about the mesh to the screen.
298  mesh.print_info();
299 
300  // Read in the solution stored in "saved_solution.xda"
301  equation_systems.read("saved_solution.xdr", DECODE);
302  }
303 
304  // Initialize the data structures for the equation system.
305  if (!read_solution)
306  equation_systems.init ();
307  else
308  equation_systems.reinit ();
309 
310  // Print out the H1 norm of the initialized or saved solution, for
311  // verification purposes:
312  Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
313 
314  libMesh::out << "Initial H1 norm = " << H1norm << std::endl << std::endl;
315 
316  // Prints information about the system to the screen.
317  equation_systems.print_info();
318 
319  equation_systems.parameters.set<unsigned int>
320  ("linear solver maximum iterations") = 250;
321  equation_systems.parameters.set<Real>
322  ("linear solver tolerance") = TOLERANCE;
323 
324  if (!read_solution)
325  {
326  // Write out the initial condition
327 #ifdef LIBMESH_HAVE_GMV
328  GMVIO(mesh).write_equation_systems ("out.gmv.000",
329  equation_systems);
330 #endif
331 #ifdef LIBMESH_HAVE_EXODUS_API
333  equation_systems);
334 #endif
335  }
336  else
337  {
338  // Write out the solution that was read in
339 #ifdef LIBMESH_HAVE_GMV
340  GMVIO(mesh).write_equation_systems ("solution_read_in.gmv",
341  equation_systems);
342 #endif
343 #ifdef LIBMESH_HAVE_EXODUS_API
344  ExodusII_IO(mesh).write_equation_systems ("solution_read_in.e",
345  equation_systems);
346 #endif
347  }
348 
349 
350  // The Convection-Diffusion system requires that we specify
351  // the flow velocity. We will specify it as a RealVectorValue
352  // data type and then use the Parameters object to pass it to
353  // the assemble function.
354  equation_systems.parameters.set<RealVectorValue>("velocity") =
355  RealVectorValue (0.8, 0.8);
356 
357  // The Convection-Diffusion system also requires a specified
358  // diffusivity. We use an isotropic (hence Real) value.
359  equation_systems.parameters.set<Real>("diffusivity") = 0.01;
360 
361  // Solve the system "Convection-Diffusion". This will be done by
362  // looping over the specified time interval and calling the
363  // solve() member at each time step. This will assemble the
364  // system and call the linear solver.
365  const Real dt = 0.025;
366  system.time = init_timestep*dt;
367 
368  // Tell the MeshRefinement object about the periodic boundaries so
369  // that it can get heuristics like level-one conformity and
370  // unrefined island elimination right.
371  mesh_refinement.set_periodic_boundaries_ptr(dof_map.get_periodic_boundaries());
372 
373  // We do 25 timesteps both before and after writing out the
374  // intermediate solution
375  for (unsigned int t_step=init_timestep; t_step<(init_timestep+n_timesteps); t_step++)
376  {
377  // Increment the time counter, set the time and the
378  // time step size as parameters in the EquationSystem.
379  system.time += dt;
380 
381  equation_systems.parameters.set<Real> ("time") = system.time;
382  equation_systems.parameters.set<Real> ("dt") = dt;
383 
384  // A pretty update message
385  libMesh::out << " Solving time step ";
386 
387  {
388  // Save flags to avoid polluting cout with custom precision values, etc.
389  std::ios_base::fmtflags os_flags = libMesh::out.flags();
390 
391  libMesh::out << t_step
392  << ", time="
393  << std::setw(6)
394  << std::setprecision(3)
395  << std::setfill('0')
396  << std::left
397  << system.time
398  << "..."
399  << std::endl;
400 
401  // Restore flags
402  libMesh::out.flags(os_flags);
403  }
404 
405  // At this point we need to update the old
406  // solution vector. The old solution vector
407  // will be the current solution vector from the
408  // previous time step. We will do this by extracting the
409  // system from the EquationSystems object and using
410  // vector assignment. Since only TransientLinearImplicitSystems
411  // (and systems derived from them) contain old solutions
412  // we need to specify the system type when we ask for it.
413  *system.old_local_solution = *system.current_local_solution;
414 
415  // The number of refinement steps per time step.
416  const unsigned int max_r_steps =
417  libMesh::command_line_next("-max_r_steps", 1);
418 
419  // A refinement loop.
420  for (unsigned int r_step=0; r_step<max_r_steps+1; r_step++)
421  {
422  // Assemble & solve the linear system
423  system.solve();
424 
425  // Print out the H1 norm, for verification purposes:
426  H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
427  libMesh::out << "H1 norm = " << H1norm << std::endl;
428 
429  // Possibly refine the mesh
430  if (r_step+1 <= max_r_steps)
431  {
432  libMesh::out << " Refining the mesh..." << std::endl;
433 
434  // The ErrorVector is a particular StatisticsVector
435  // for computing error information on a finite element mesh.
436  ErrorVector error;
437 
438  // The ErrorEstimator class interrogates a finite element
439  // solution and assigns to each element a positive error value.
440  // This value is used for deciding which elements to refine
441  // and which to coarsen.
442  KellyErrorEstimator error_estimator;
443 
444  // This is a subclass of JumpErrorEstimator, based
445  // on measuring discontinuities across sides between
446  // elements, and we can tell it to use a cheaper
447  // "unweighted" quadrature rule when numerically
448  // integrating those discontinuities.
449  error_estimator.use_unweighted_quadrature_rules = true;
450 
451  // Compute the error for each active element using the provided
452  // flux_jump indicator. Note in general you will need to
453  // provide an error estimator specifically designed for your
454  // application.
455  error_estimator.estimate_error (system,
456  error);
457 
458  // This takes the error in error and decides which elements
459  // will be coarsened or refined. Any element within 20% of the
460  // maximum error on any element will be refined, and any
461  // element within 7% of the minimum error on any element might
462  // be coarsened. Note that the elements flagged for refinement
463  // will be refined, but those flagged for coarsening _might_ be
464  // coarsened.
465  mesh_refinement.refine_fraction() = 0.80;
466  mesh_refinement.coarsen_fraction() = 0.07;
467  mesh_refinement.max_h_level() = 5;
468  mesh_refinement.flag_elements_by_error_fraction (error);
469 
470  // This call actually refines and coarsens the flagged
471  // elements.
472  mesh_refinement.refine_and_coarsen_elements();
473 
474  // This call reinitializes the EquationSystems object for
475  // the newly refined mesh. One of the steps in the
476  // reinitialization is projecting the solution,
477  // old_solution, etc... vectors from the old mesh to
478  // the current one.
479  equation_systems.reinit ();
480  }
481  }
482 
483  // Again do a search on the command line for an argument
484  const unsigned int output_freq =
485  libMesh::command_line_next("-output_freq", 10);
486 
487  // Output every 10 timesteps to file.
488  if ((t_step+1)%output_freq == 0)
489  {
490 #ifdef LIBMESH_HAVE_GMV
491  // std::ostringstream file_name;
492  // out << "out.gmv."
493  // << std::setw(3)
494  // << std::setfill('0')
495  // << std::right
496  // << t_step+1;
497  // GMVIO(mesh).write_equation_systems (file_name.str(),
498  // equation_systems);
499 #endif
500 #ifdef LIBMESH_HAVE_EXODUS_API
501  // So... if paraview is told to open a file called out.e.{N}, it automatically tries to
502  // open out.e.{N-1}, out.e.{N-2}, etc. If we name the file something else, we can work
503  // around that issue, but the right thing to do (for adaptive meshes) is to write a filename
504  // with the adaptation step into a separate file.
506  equation_systems);
507 #endif
508  }
509  }
510 
511  if (!read_solution)
512  {
513  // Print out the H1 norm of the saved solution, for verification purposes:
514  H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
515 
516  libMesh::out << "Final H1 norm = " << H1norm << std::endl << std::endl;
517 
518  mesh.write("saved_mesh.xdr");
519  equation_systems.write("saved_solution.xdr", ENCODE);
520 #ifdef LIBMESH_HAVE_GMV
521  GMVIO(mesh).write_equation_systems ("saved_solution.gmv",
522  equation_systems);
523 #endif
524 #ifdef LIBMESH_HAVE_EXODUS_API
525  ExodusII_IO(mesh).write_equation_systems ("saved_solution.e",
526  equation_systems);
527 #endif
528  }
529 #endif // #ifndef LIBMESH_ENABLE_AMR
530 
531  return 0;
532 }
533 
534 // Here we define the initialization routine for the
535 // Convection-Diffusion system. This routine is
536 // responsible for applying the initial conditions to
537 // the system.
539  const std::string & libmesh_dbg_var(system_name))
540 {
541  // It is a good idea to make sure we are initializing
542  // the proper system.
543  libmesh_assert_equal_to (system_name, "Convection-Diffusion");
544 
545  // Get a reference to the Convection-Diffusion system object.
547  es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
548 
549  // Project initial conditions at time 0
550  es.parameters.set<Real> ("time") = system.time = 0;
551 
552  if (parsed_solution.get())
553  system.project_solution(parsed_solution.get(), nullptr);
554  else
555  system.project_solution(exact_value, nullptr, es.parameters);
556 }
557 
558 
559 
560 // This function defines the assembly routine which
561 // will be called at each time step. It is responsible
562 // for computing the proper matrix entries for the
563 // element stiffness matrices and right-hand sides.
564 #ifdef LIBMESH_ENABLE_AMR
566  const std::string & libmesh_dbg_var(system_name))
567 {
568  // It is a good idea to make sure we are assembling
569  // the proper system.
570  libmesh_assert_equal_to (system_name, "Convection-Diffusion");
571 
572  // Get a constant reference to the mesh object.
573  const MeshBase & mesh = es.get_mesh();
574 
575  // The dimension that we are running
576  const unsigned int dim = mesh.mesh_dimension();
577 
578  // Get a reference to the Convection-Diffusion system object.
580  es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
581 
582  // Get the Finite Element type for the first (and only)
583  // variable in the system.
584  FEType fe_type = system.variable_type(0);
585 
586  // Build a Finite Element object of the specified type. Since the
587  // FEBase::build() member dynamically creates memory we will
588  // store the object as a std::unique_ptr<FEBase>. This can be thought
589  // of as a pointer that will clean up after itself.
590  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
591  std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
592 
593  // A Gauss quadrature rule for numerical integration.
594  // Let the FEType object decide what order rule is appropriate.
595  QGauss qrule (dim, fe_type.default_quadrature_order());
596  QGauss qface (dim-1, fe_type.default_quadrature_order());
597 
598  // Tell the finite element object to use our quadrature rule.
599  fe->attach_quadrature_rule (&qrule);
600  fe_face->attach_quadrature_rule (&qface);
601 
602  // Here we define some references to cell-specific data that
603  // will be used to assemble the linear system. We will start
604  // with the element Jacobian * quadrature weight at each integration point.
605  const std::vector<Real> & JxW = fe->get_JxW();
606 
607  // The element shape functions evaluated at the quadrature points.
608  const std::vector<std::vector<Real>> & phi = fe->get_phi();
609 
610  // The element shape function gradients evaluated at the quadrature
611  // points.
612  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
613 
614  // A reference to the DofMap object for this system. The DofMap
615  // object handles the index translation from node and element numbers
616  // to degree of freedom numbers. We will talk more about the DofMap
617  // in future examples.
618  const DofMap & dof_map = system.get_dof_map();
619 
620  // Define data structures to contain the element matrix
621  // and right-hand-side vector contribution. Following
622  // basic finite element terminology we will denote these
623  // "Ke" and "Fe".
626 
627  // This vector will hold the degree of freedom indices for
628  // the element. These define where in the global system
629  // the element degrees of freedom get mapped.
630  std::vector<dof_id_type> dof_indices;
631 
632  // Here we extract the velocity & parameters that we put in the
633  // EquationSystems object.
634  const RealVectorValue velocity =
635  es.parameters.get<RealVectorValue> ("velocity");
636 
637  const Real diffusivity =
638  es.parameters.get<Real> ("diffusivity");
639 
640  const Real dt = es.parameters.get<Real> ("dt");
641 
642  // The global system matrix
643  SparseMatrix<Number> & matrix = system.get_system_matrix();
644 
645  // Now we will loop over all the elements in the mesh that
646  // live on the local processor. We will compute the element
647  // matrix and right-hand-side contribution. Since the mesh
648  // will be refined we want to only consider the ACTIVE elements,
649  // hence we use a variant of the active_elem_iterator.
650  for (const auto & elem : mesh.active_local_element_ptr_range())
651  {
652  // Get the degree of freedom indices for the
653  // current element. These define where in the global
654  // matrix and right-hand-side this element will
655  // contribute to.
656  dof_map.dof_indices (elem, dof_indices);
657 
658  // Compute the element-specific data for the current
659  // element. This involves computing the location of the
660  // quadrature points (q_point) and the shape functions
661  // (phi, dphi) for the current element.
662  fe->reinit (elem);
663 
664  const unsigned int n_dofs =
665  cast_int<unsigned int>(dof_indices.size());
666  libmesh_assert_equal_to (n_dofs, phi.size());
667 
668  // Zero the element matrix and right-hand side before
669  // summing them. We use the resize member here because
670  // the number of degrees of freedom might have changed from
671  // the last element. Note that this will be the case if the
672  // element type is different (i.e. the last element was a
673  // triangle, now we are on a quadrilateral).
674  Ke.resize (n_dofs, n_dofs);
675 
676  Fe.resize (n_dofs);
677 
678  // Now we will build the element matrix and right-hand-side.
679  // Constructing the RHS requires the solution and its
680  // gradient from the previous timestep. This myst be
681  // calculated at each quadrature point by summing the
682  // solution degree-of-freedom values by the appropriate
683  // weight functions.
684  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
685  {
686  // Values to hold the old solution & its gradient.
687  Number u_old = 0.;
688  Gradient grad_u_old;
689 
690  // Compute the old solution & its gradient.
691  for (unsigned int l=0; l != n_dofs; l++)
692  {
693  u_old += phi[l][qp]*system.old_solution (dof_indices[l]);
694 
695  // This will work,
696  // grad_u_old += dphi[l][qp]*system.old_solution (dof_indices[l]);
697  // but we can do it without creating a temporary like this:
698  grad_u_old.add_scaled (dphi[l][qp], system.old_solution (dof_indices[l]));
699  }
700 
701  // Now compute the element matrix and RHS contributions.
702  for (unsigned int i=0; i != n_dofs; i++)
703  {
704  // The RHS contribution
705  Fe(i) += JxW[qp]*(
706  // Mass matrix term
707  u_old*phi[i][qp] +
708  -.5*dt*(
709  // Convection term
710  // (grad_u_old may be complex, so the
711  // order here is important!)
712  (grad_u_old*velocity)*phi[i][qp] +
713 
714  // Diffusion term
715  diffusivity*(grad_u_old*dphi[i][qp]))
716  );
717 
718  for (unsigned int j=0; j != n_dofs; j++)
719  {
720  // The matrix contribution
721  Ke(i,j) += JxW[qp]*(
722  // Mass-matrix
723  phi[i][qp]*phi[j][qp] +
724  .5*dt*(
725  // Convection term
726  (velocity*dphi[j][qp])*phi[i][qp] +
727  // Diffusion term
728  diffusivity*(dphi[i][qp]*dphi[j][qp]))
729  );
730  }
731  }
732  }
733 
734  // We have now built the element matrix and RHS vector in terms
735  // of the element degrees of freedom. However, it is possible
736  // that some of the element DOFs are constrained to enforce
737  // solution continuity, i.e. they are not really "free". We need
738  // to constrain those DOFs in terms of non-constrained DOFs to
739  // ensure a continuous solution. The
740  // DofMap::constrain_element_matrix_and_vector() method does
741  // just that.
742  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
743 
744  // The element matrix and right-hand-side are now built
745  // for this element. Add them to the global matrix and
746  // right-hand-side vector. The SparseMatrix::add_matrix()
747  // and NumericVector::add_vector() members do this for us.
748  matrix.add_matrix (Ke, dof_indices);
749  system.rhs->add_vector (Fe, dof_indices);
750 
751  }
752  // Finished computing the system matrix and right-hand side.
753 }
754 #endif // #ifdef LIBMESH_ENABLE_AMR
755 
756 
757 
758 
759 std::string exodus_filename(unsigned number)
760 {
761  std::ostringstream oss;
762 
763  oss << "out_";
764  oss << std::setw(3) << std::setfill('0') << number;
765  oss << ".e";
766 
767  return oss.str();
768 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
void assemble_cd(EquationSystems &es, const std::string &system_name)
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1078
Real exact_solution(const Real x, const Real y, const Real t)
This is the exact solution that we are trying to obtain.
This is the EquationSystems class.
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
void add_periodic_boundary(const PeriodicBoundaryBase &periodic_boundary)
Adds a copy of the specified periodic boundary to the system.
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:67
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:2330
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=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 dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
VectorValue< Real > RealVectorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
static constexpr Real TOLERANCE
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
int main(int argc, char **argv)
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
Real & coarsen_fraction()
The coarsen_fraction sets either a desired target or a desired maximum number of elements to flag for...
MeshBase & mesh
Manages storage and variables for transient systems.
Real & refine_fraction()
The refine_fraction sets either a desired target or a desired maximum number of elements to flag for ...
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
Definition: system_norm.h:49
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.
void init_cd(EquationSystems &es, const std::string &system_name)
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:90
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
The definition of a periodic boundary.
const T_sys & get_system(std::string_view name) const
unsigned int & max_h_level()
The max_h_level is the greatest refinement level an element should reach.
NumericVector< Number > * old_local_solution
All the values I need to compute my contribution to the simulation at hand.
This is the MeshBase class.
Definition: mesh_base.h:75
SolverPackage default_solver_package()
Definition: libmesh.C:1117
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
boundary_id_type myboundary
The boundary ID of this boundary and its counterpart.
const T & get(std::string_view) const
Definition: parameters.h:426
Implements (adaptive) mesh refinement algorithms for a MeshBase.
std::string exodus_filename(unsigned number)
PeriodicBoundaries * get_periodic_boundaries()
Definition: dof_map.h:1455
virtual void write_equation_systems(const std::string &fname, const EquationSystems &es, const std::set< std::string > *system_names=nullptr) override
Writes out the solution for no specific time or timestep.
Definition: exodusII_io.C:2033
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1562
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
std::unique_ptr< FunctionBase< Number > > parsed_solution
Number old_solution(const dof_id_type global_dof_number) const
virtual void write(const std::string &name) const =0
This class implements the Kelly error indicator which is based on the flux jumps between elements...
std::ios_base::fmtflags flags() const
Get the associated format flags.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
bool refine_and_coarsen_elements()
Refines and coarsens user-requested elements.
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
bool use_unweighted_quadrature_rules
This boolean flag allows you to use "unweighted" quadrature rules (sized to exactly integrate unweigh...
This class implements specific orders of Gauss quadrature.
Number exact_value(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
Parameters parameters
Data structure holding arbitrary parameters.
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...
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override
This function uses the derived class&#39;s jump error estimate formula to estimate the error on each cell...
bool on_command_line(std::string arg)
Definition: libmesh.C:987
void set_periodic_boundaries_ptr(PeriodicBoundaries *pb_ptr)
Sets the PeriodicBoundaries pointer.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.