libMesh
adaptivity_ex2.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 2 - Solving a Transient System with Adaptive Mesh Refinement</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This example shows how a simple, linear transient
25 // system can be solved in parallel. The system is simple
26 // scalar convection-diffusion with a specified external
27 // velocity. The initial condition is given, and the
28 // solution is advanced in time with a standard Crank-Nicolson
29 // time-stepping strategy.
30 //
31 // Also, we use this example to demonstrate writing out and reading
32 // in of solutions. We do 25 time steps, then save the solution
33 // and do another 25 time steps starting from the saved solution.
34 
35 // C++ include files that we need
36 #include <iostream>
37 #include <algorithm>
38 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
39 #include <cmath>
40 
41 // Basic include file needed for the mesh functionality.
42 #include "libmesh/libmesh.h"
43 #include "libmesh/replicated_mesh.h"
44 #include "libmesh/mesh_refinement.h"
45 #include "libmesh/gmv_io.h"
46 #include "libmesh/equation_systems.h"
47 #include "libmesh/fe.h"
48 #include "libmesh/quadrature_gauss.h"
49 #include "libmesh/dof_map.h"
50 #include "libmesh/sparse_matrix.h"
51 #include "libmesh/numeric_vector.h"
52 #include "libmesh/dense_matrix.h"
53 #include "libmesh/dense_vector.h"
54 
55 #include "libmesh/getpot.h"
56 
57 // This example will solve a linear transient system,
58 // so we need to include the TransientLinearImplicitSystem definition.
59 #include "libmesh/transient_system.h"
60 #include "libmesh/linear_implicit_system.h"
61 #include "libmesh/vector_value.h"
62 
63 // To refine the mesh we need an ErrorEstimator
64 // object to figure out which elements to refine.
65 #include "libmesh/error_vector.h"
66 #include "libmesh/kelly_error_estimator.h"
67 
68 // The definition of a geometric element
69 #include "libmesh/elem.h"
70 
71 // Bring in everything from the libMesh namespace
72 using namespace libMesh;
73 
74 // Function prototype. This function will assemble the system
75 // matrix and right-hand-side at each time step. Note that
76 // since the system is linear we technically do not need to
77 // assemble the matrix at each time step, but we will anyway.
78 // In subsequent examples we will employ adaptive mesh refinement,
79 // and with a changing mesh it will be necessary to rebuild the
80 // system matrix.
81 #ifdef LIBMESH_ENABLE_AMR
82 void assemble_cd (EquationSystems & es,
83  const std::string & system_name);
84 #endif
85 
86 // Function prototype. This function will initialize the system.
87 // Initialization functions are optional for systems. They allow
88 // you to specify the initial values of the solution. If an
89 // initialization function is not provided then the default (0)
90 // solution is provided.
91 void init_cd (EquationSystems & es,
92  const std::string & system_name);
93 
94 // Exact solution function prototype. This gives the exact
95 // solution as a function of space and time. In this case the
96 // initial condition will be taken as the exact solution at time 0,
97 // as will the Dirichlet boundary conditions at time t.
98 Real exact_solution (const Real x,
99  const Real y,
100  const Real t);
101 
103  const Parameters & parameters,
104  const std::string &,
105  const std::string &)
106 {
107  return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
108 }
109 
110 
111 
112 // Begin the main program. Note that the first
113 // statement in the program throws an error if
114 // you are in complex number mode, since this
115 // example is only intended to work with real
116 // numbers.
117 int main (int argc, char ** argv)
118 {
119  // Initialize libMesh.
120  LibMeshInit init (argc, argv);
121 
122  // This example requires a linear solver package.
123  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
124  "--enable-petsc, --enable-trilinos, or --enable-eigen");
125 
126 #ifndef LIBMESH_ENABLE_AMR
127  libmesh_example_requires(false, "--enable-amr");
128 #else
129 
130  // Our Trilinos interface does not yet support adaptive transient
131  // problems
132  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc");
133 
134  // Brief message to the user regarding the program name
135  // and command line arguments.
136 
137  // Use commandline parameter to specify if we are to
138  // read in an initial solution or generate it ourself
139  libMesh::out << "Usage:\n"
140  <<"\t " << argv[0] << " -init_timestep 0 -n_timesteps 25 [-n_refinements 5]\n"
141  << "OR\n"
142  <<"\t " << argv[0] << " -read_solution -init_timestep 26 -n_timesteps 25\n"
143  << std::endl;
144 
145  libMesh::out << "Running: " << argv[0];
146 
147  for (int i=1; i<argc; i++)
148  libMesh::out << " " << argv[i];
149 
150  libMesh::out << std::endl << std::endl;
151 
152  // Create a GetPot object to parse the command line
153  GetPot command_line (argc, argv);
154 
155 
156  // This boolean value is obtained from the command line, it is true
157  // if the flag "-read_solution" is present, false otherwise.
158  // It indicates whether we are going to read in
159  // the mesh and solution files "saved_mesh.xda" and "saved_solution.xda"
160  // or whether we are going to start from scratch by just reading
161  // "mesh.xda"
162  const bool read_solution = command_line.search("-read_solution");
163 
164  // This value is also obtained from the commandline and it specifies the
165  // initial value for the t_step looping variable. We must
166  // distinguish between the two cases here, whether we read in the
167  // solution or we started from scratch, so that we do not overwrite the
168  // gmv output files.
169  unsigned int init_timestep = 0;
170 
171  // Search the command line for the "init_timestep" flag and if it is
172  // present, set init_timestep accordingly.
173  if (command_line.search("-init_timestep"))
174  init_timestep = command_line.next(0);
175  else
176  {
177  // This handy function will print the file name, line number,
178  // specified message, and then throw an exception.
179  libmesh_error_msg("ERROR: Initial timestep not specified!");
180  }
181 
182  // This value is also obtained from the command line, and specifies
183  // the number of time steps to take.
184  unsigned int n_timesteps = 0;
185 
186  // Again do a search on the command line for the argument
187  if (command_line.search("-n_timesteps"))
188  n_timesteps = command_line.next(0);
189  else
190  libmesh_error_msg("ERROR: Number of timesteps not specified");
191 
192 
193  // Skip this 2D example if libMesh was compiled as 1D-only.
194  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
195 
196  // Create a new mesh on the default MPI communicator.
197  // We still need some work on automatic parallel restarts with
198  // DistributedMesh
199  ReplicatedMesh mesh(init.comm());
200 
201  // Create an equation systems object.
202  EquationSystems equation_systems (mesh);
203  MeshRefinement mesh_refinement (mesh);
204 
205  // First we process the case where we do not read in the solution
206  if (!read_solution)
207  {
208  // Read the mesh from file.
209  mesh.read ("mesh.xda");
210 
211  // Again do a search on the command line for an argument
212  unsigned int n_refinements = 5;
213  if (command_line.search("-n_refinements"))
214  n_refinements = command_line.next(0);
215 
216  // Uniformly refine the mesh 5 times
217  if (!read_solution)
218  mesh_refinement.uniformly_refine (n_refinements);
219 
220  // Print information about the mesh to the screen.
221  mesh.print_info();
222 
223  // Declare the system and its variables.
224  // Begin by creating a transient system
225  // named "Convection-Diffusion".
227  equation_systems.add_system<TransientLinearImplicitSystem>("Convection-Diffusion");
228 
229  // Adds the variable "u" to "Convection-Diffusion". "u"
230  // will be approximated using first-order approximation.
231  system.add_variable ("u", FIRST);
232 
233  // Give the system a pointer to the matrix assembly
234  // and initialization functions.
235  system.attach_assemble_function (assemble_cd);
236  system.attach_init_function (init_cd);
237 
238  // Initialize the data structures for the equation system.
239  equation_systems.init ();
240  }
241  // Otherwise we read in the solution and mesh
242  else
243  {
244  // Read in the mesh stored in "saved_mesh.xda"
245  mesh.read("saved_mesh.xda");
246 
247  // Print information about the mesh to the screen.
248  mesh.print_info();
249 
250  // Read in the solution stored in "saved_solution.xda"
251  equation_systems.read("saved_solution.xda", READ);
252 
253  // Get a reference to the system so that we can call update() on it
255  equation_systems.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
256 
257  // We need to call update to put system in a consistent state
258  // with the solution that was read in
259  system.update();
260 
261  // Attach the same matrix assembly function as above. Note, we do not
262  // have to attach an init() function since we are initializing the
263  // system by reading in "saved_solution.xda"
264  system.attach_assemble_function (assemble_cd);
265 
266  // Print out the H1 norm of the saved solution, for verification purposes:
267  Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
268 
269  libMesh::out << "Initial H1 norm = " << H1norm << std::endl << std::endl;
270  }
271 
272  // Prints information about the system to the screen.
273  equation_systems.print_info();
274 
275  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;
276  equation_systems.parameters.set<Real>("linear solver tolerance") = TOLERANCE;
277 
278  if (!read_solution)
279  // Write out the initial condition
280  GMVIO(mesh).write_equation_systems ("out.gmv.000", equation_systems);
281  else
282  // Write out the solution that was read in
283  GMVIO(mesh).write_equation_systems ("solution_read_in.gmv", equation_systems);
284 
285  // The Convection-Diffusion system requires that we specify
286  // the flow velocity. We will specify it as a RealVectorValue
287  // data type and then use the Parameters object to pass it to
288  // the assemble function.
289  equation_systems.parameters.set<RealVectorValue>("velocity") =
290  RealVectorValue (0.8, 0.8);
291 
292  // The Convection-Diffusion system also requires a specified
293  // diffusivity. We use an isotropic (hence Real) value.
294  equation_systems.parameters.set<Real>("diffusivity") = 0.01;
295 
296  // Solve the system "Convection-Diffusion". This will be done by
297  // looping over the specified time interval and calling the
298  // solve() member at each time step. This will assemble the
299  // system and call the linear solver.
300 
301  // Since only TransientLinearImplicitSystems (and systems
302  // derived from them) contain old solutions, to use the
303  // old_local_solution later we now need to specify the system
304  // type when we ask for it.
306  equation_systems.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
307 
308  const Real dt = 0.025;
309  system.time = init_timestep*dt;
310 
311  // We do 25 timesteps both before and after writing out the
312  // intermediate solution
313  for (unsigned int t_step=init_timestep; t_step<(init_timestep+n_timesteps); t_step++)
314  {
315  // Increment the time counter, set the time and the
316  // time step size as parameters in the EquationSystem.
317  system.time += dt;
318 
319  equation_systems.parameters.set<Real> ("time") = system.time;
320  equation_systems.parameters.set<Real> ("dt") = dt;
321 
322  // A pretty update message
323  libMesh::out << " Solving time step ";
324 
325  // Add a set of scope braces to enforce data locality.
326  {
327  std::ostringstream out;
328 
329  out << std::setw(2)
330  << std::right
331  << t_step
332  << ", time="
333  << std::fixed
334  << std::setw(6)
335  << std::setprecision(3)
336  << std::setfill('0')
337  << std::left
338  << system.time
339  << "...";
340 
341  libMesh::out << out.str() << std::endl;
342  }
343 
344  // At this point we need to update the old
345  // solution vector. The old solution vector
346  // will be the current solution vector from the
347  // previous time step.
348 
349  *system.old_local_solution = *system.current_local_solution;
350 
351  // The number of refinement steps per time step.
352  const unsigned int max_r_steps = 2;
353 
354  // A refinement loop.
355  for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
356  {
357  // Assemble & solve the linear system
358  system.solve();
359 
360  // Print out the H1 norm, for verification purposes:
361  Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
362 
363  libMesh::out << "H1 norm = " << H1norm << std::endl;
364 
365  // Possibly refine the mesh
366  if (r_step+1 != max_r_steps)
367  {
368  libMesh::out << " Refining the mesh..." << std::endl;
369 
370  // The ErrorVector is a particular StatisticsVector
371  // for computing error information on a finite element mesh.
372  ErrorVector error;
373 
374  // The ErrorEstimator class interrogates a finite element
375  // solution and assigns to each element a positive error value.
376  // This value is used for deciding which elements to refine
377  // and which to coarsen.
378  //ErrorEstimator* error_estimator = new KellyErrorEstimator;
379  KellyErrorEstimator error_estimator;
380 
381  // Compute the error for each active element using the provided
382  // flux_jump indicator. Note in general you will need to
383  // provide an error estimator specifically designed for your
384  // application.
385  error_estimator.estimate_error (system,
386  error);
387 
388  // This takes the error in error and decides which elements
389  // will be coarsened or refined. Any element within 20% of the
390  // maximum error on any element will be refined, and any
391  // element within 7% of the minimum error on any element might
392  // be coarsened. Note that the elements flagged for refinement
393  // will be refined, but those flagged for coarsening _might_ be
394  // coarsened.
395  mesh_refinement.refine_fraction() = 0.80;
396  mesh_refinement.coarsen_fraction() = 0.07;
397  mesh_refinement.max_h_level() = 5;
398  mesh_refinement.flag_elements_by_error_fraction (error);
399 
400  // This call actually refines and coarsens the flagged
401  // elements.
402  mesh_refinement.refine_and_coarsen_elements();
403 
404  // This call reinitializes the EquationSystems object for
405  // the newly refined mesh. One of the steps in the
406  // reinitialization is projecting the solution,
407  // old_solution, etc... vectors from the old mesh to
408  // the current one.
409  equation_systems.reinit ();
410  }
411  }
412 
413  // Again do a search on the command line for an argument
414  unsigned int output_freq = 10;
415  if (command_line.search("-output_freq"))
416  output_freq = command_line.next(0);
417 
418  // Output every 10 timesteps to file.
419  if ((t_step+1)%output_freq == 0)
420  {
421  std::ostringstream file_name;
422 
423  file_name << "out.gmv."
424  << std::setw(3)
425  << std::setfill('0')
426  << std::right
427  << t_step+1;
428 
429  GMVIO(mesh).write_equation_systems (file_name.str(),
430  equation_systems);
431  }
432  }
433 
434  if (!read_solution)
435  {
436  // Print out the H1 norm of the saved solution, for verification purposes:
437  Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
438 
439  libMesh::out << "Final H1 norm = " << H1norm << std::endl << std::endl;
440 
441  mesh.write("saved_mesh.xda");
442  equation_systems.write("saved_solution.xda", WRITE);
443  GMVIO(mesh).write_equation_systems ("saved_solution.gmv",
444  equation_systems);
445  }
446 #endif // #ifndef LIBMESH_ENABLE_AMR
447 
448  return 0;
449 }
450 
451 // Here we define the initialization routine for the
452 // Convection-Diffusion system. This routine is
453 // responsible for applying the initial conditions to
454 // the system.
456  const std::string & libmesh_dbg_var(system_name))
457 {
458  // It is a good idea to make sure we are initializing
459  // the proper system.
460  libmesh_assert_equal_to (system_name, "Convection-Diffusion");
461 
462  // Get a reference to the Convection-Diffusion system object.
464  es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
465 
466  // Project initial conditions at time 0
467  es.parameters.set<Real> ("time") = system.time = 0;
468 
469  system.project_solution(exact_value, libmesh_nullptr, es.parameters);
470 }
471 
472 
473 
474 // This function defines the assembly routine which
475 // will be called at each time step. It is responsible
476 // for computing the proper matrix entries for the
477 // element stiffness matrices and right-hand sides.
478 #ifdef LIBMESH_ENABLE_AMR
480  const std::string & libmesh_dbg_var(system_name))
481 {
482  // It is a good idea to make sure we are assembling
483  // the proper system.
484  libmesh_assert_equal_to (system_name, "Convection-Diffusion");
485 
486  // Get a constant reference to the mesh object.
487  const MeshBase & mesh = es.get_mesh();
488 
489  // The dimension that we are running
490  const unsigned int dim = mesh.mesh_dimension();
491 
492  // Get a reference to the Convection-Diffusion system object.
494  es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
495 
496  // Get the Finite Element type for the first (and only)
497  // variable in the system.
498  FEType fe_type = system.variable_type(0);
499 
500  // Build a Finite Element object of the specified type. Since the
501  // FEBase::build() member dynamically creates memory we will
502  // store the object as a UniquePtr<FEBase>. This can be thought
503  // of as a pointer that will clean up after itself.
504  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
505  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
506 
507  // A Gauss quadrature rule for numerical integration.
508  // Let the FEType object decide what order rule is appropriate.
509  QGauss qrule (dim, fe_type.default_quadrature_order());
510  QGauss qface (dim-1, fe_type.default_quadrature_order());
511 
512  // Tell the finite element object to use our quadrature rule.
513  fe->attach_quadrature_rule (&qrule);
514  fe_face->attach_quadrature_rule (&qface);
515 
516  // Here we define some references to cell-specific data that
517  // will be used to assemble the linear system. We will start
518  // with the element Jacobian * quadrature weight at each integration point.
519  const std::vector<Real> & JxW = fe->get_JxW();
520  const std::vector<Real> & JxW_face = fe_face->get_JxW();
521 
522  // The element shape functions evaluated at the quadrature points.
523  const std::vector<std::vector<Real>> & phi = fe->get_phi();
524  const std::vector<std::vector<Real>> & psi = fe_face->get_phi();
525 
526  // The element shape function gradients evaluated at the quadrature
527  // points.
528  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
529 
530  // The XY locations of the quadrature points used for face integration
531  const std::vector<Point> & qface_points = fe_face->get_xyz();
532 
533  // A reference to the DofMap object for this system. The DofMap
534  // object handles the index translation from node and element numbers
535  // to degree of freedom numbers. We will talk more about the DofMap
536  // in future examples.
537  const DofMap & dof_map = system.get_dof_map();
538 
539  // Define data structures to contain the element matrix
540  // and right-hand-side vector contribution. Following
541  // basic finite element terminology we will denote these
542  // "Ke" and "Fe".
545 
546  // This vector will hold the degree of freedom indices for
547  // the element. These define where in the global system
548  // the element degrees of freedom get mapped.
549  std::vector<dof_id_type> dof_indices;
550 
551  // Here we extract the velocity & parameters that we put in the
552  // EquationSystems object.
553  const RealVectorValue velocity =
554  es.parameters.get<RealVectorValue> ("velocity");
555 
556  const Real diffusivity =
557  es.parameters.get<Real> ("diffusivity");
558 
559  const Real dt = es.parameters.get<Real> ("dt");
560 
561  // Now we will loop over all the elements in the mesh that
562  // live on the local processor. We will compute the element
563  // matrix and right-hand-side contribution. Since the mesh
564  // will be refined we want to only consider the ACTIVE elements,
565  // hence we use a variant of the active_elem_iterator.
566  for (const auto & elem : mesh.active_local_element_ptr_range())
567  {
568  // Get the degree of freedom indices for the
569  // current element. These define where in the global
570  // matrix and right-hand-side this element will
571  // contribute to.
572  dof_map.dof_indices (elem, dof_indices);
573 
574  // Compute the element-specific data for the current
575  // element. This involves computing the location of the
576  // quadrature points (q_point) and the shape functions
577  // (phi, dphi) for the current element.
578  fe->reinit (elem);
579 
580  // Zero the element matrix and right-hand side before
581  // summing them. We use the resize member here because
582  // the number of degrees of freedom might have changed from
583  // the last element. Note that this will be the case if the
584  // element type is different (i.e. the last element was a
585  // triangle, now we are on a quadrilateral).
586  Ke.resize (dof_indices.size(),
587  dof_indices.size());
588 
589  Fe.resize (dof_indices.size());
590 
591  // Now we will build the element matrix and right-hand-side.
592  // Constructing the RHS requires the solution and its
593  // gradient from the previous timestep. This myst be
594  // calculated at each quadrature point by summing the
595  // solution degree-of-freedom values by the appropriate
596  // weight functions.
597  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
598  {
599  // Values to hold the old solution & its gradient.
600  Number u_old = 0.;
601  Gradient grad_u_old;
602 
603  // Compute the old solution & its gradient.
604  for (std::size_t l=0; l<phi.size(); l++)
605  {
606  u_old += phi[l][qp]*system.old_solution (dof_indices[l]);
607 
608  // This will work,
609  // grad_u_old += dphi[l][qp]*system.old_solution (dof_indices[l]);
610  // but we can do it without creating a temporary like this:
611  grad_u_old.add_scaled (dphi[l][qp], system.old_solution (dof_indices[l]));
612  }
613 
614  // Now compute the element matrix and RHS contributions.
615  for (std::size_t i=0; i<phi.size(); i++)
616  {
617  // The RHS contribution
618  Fe(i) += JxW[qp]*(
619  // Mass matrix term
620  u_old*phi[i][qp] +
621  -.5*dt*(
622  // Convection term
623  // (grad_u_old may be complex, so the
624  // order here is important!)
625  (grad_u_old*velocity)*phi[i][qp] +
626 
627  // Diffusion term
628  diffusivity*(grad_u_old*dphi[i][qp]))
629  );
630 
631  for (std::size_t j=0; j<phi.size(); j++)
632  {
633  // The matrix contribution
634  Ke(i,j) += JxW[qp]*(
635  // Mass-matrix
636  phi[i][qp]*phi[j][qp] +
637  .5*dt*(
638  // Convection term
639  (velocity*dphi[j][qp])*phi[i][qp] +
640  // Diffusion term
641  diffusivity*(dphi[i][qp]*dphi[j][qp]))
642  );
643  }
644  }
645  }
646 
647  // At this point the interior element integration has
648  // been completed. However, we have not yet addressed
649  // boundary conditions. For this example we will only
650  // consider simple Dirichlet boundary conditions imposed
651  // via the penalty method.
652  //
653  // The following loops over the sides of the element.
654  // If the element has no neighbor on a side then that
655  // side MUST live on a boundary of the domain.
656  {
657  // The penalty value.
658  const Real penalty = 1.e10;
659 
660  // The following loops over the sides of the element.
661  // If the element has no neighbor on a side then that
662  // side MUST live on a boundary of the domain.
663  for (auto s : elem->side_index_range())
664  if (elem->neighbor_ptr(s) == libmesh_nullptr)
665  {
666  fe_face->reinit(elem, s);
667 
668  for (unsigned int qp=0; qp<qface.n_points(); qp++)
669  {
670  const Number value = exact_solution (qface_points[qp](0),
671  qface_points[qp](1),
672  system.time);
673 
674  // RHS contribution
675  for (std::size_t i=0; i<psi.size(); i++)
676  Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
677 
678  // Matrix contribution
679  for (std::size_t i=0; i<psi.size(); i++)
680  for (std::size_t j=0; j<psi.size(); j++)
681  Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
682  }
683  }
684  }
685 
686 
687  // We have now built the element matrix and RHS vector in terms
688  // of the element degrees of freedom. However, it is possible
689  // that some of the element DOFs are constrained to enforce
690  // solution continuity, i.e. they are not really "free". We need
691  // to constrain those DOFs in terms of non-constrained DOFs to
692  // ensure a continuous solution. The
693  // DofMap::constrain_element_matrix_and_vector() method does
694  // just that.
695  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
696 
697  // The element matrix and right-hand-side are now built
698  // for this element. Add them to the global matrix and
699  // right-hand-side vector. The SparseMatrix::add_matrix()
700  // and NumericVector::add_vector() members do this for us.
701  system.matrix->add_matrix (Ke, dof_indices);
702  system.rhs->add_vector (Fe, dof_indices);
703 
704  }
705  // Finished computing the system matrix and right-hand side.
706 }
707 #endif // #ifdef LIBMESH_ENABLE_AMR
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
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
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
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
static const Real TOLERANCE
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
unsigned int & max_h_level()
The max_h_level is the greatest refinement level an element should reach.
Number exact_value(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
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
PetscErrorCode Vec x
This is the MeshRefinement class.
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.
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
int main(int argc, char **argv)
void init_cd(EquationSystems &es, const std::string &system_name)
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.
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
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
void flag_elements_by_error_fraction(const ErrorVector &error_per_cell, const Real refine_fraction=0.3, const Real coarsen_fraction=0.0, const unsigned int max_level=libMesh::invalid_uint)
Flags elements for coarsening and refinement based on the computed error passed in error_per_cell...
Number old_solution(const dof_id_type global_dof_number) const
const T_sys & get_system(const std::string &name) const
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
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.