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