libMesh
introduction_ex4.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>Introduction Example 4 - Solving a 1D, 2D or 3D Poisson Problem in Parallel</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This is the fourth example program. It builds on
25 // the third example program by showing how to formulate
26 // the code in a dimension-independent way. Very minor
27 // changes to the example will allow the problem to be
28 // solved in one, two or three dimensions.
29 //
30 // This example will also introduce the PerfLog class
31 // as a way to monitor your code's performance. We will
32 // use it to instrument the matrix assembly code and look
33 // for bottlenecks where we should focus optimization efforts.
34 //
35 // This example also shows how to extend example 3 to run in
36 // parallel. Notice how little has changed!
37 
38 
39 // C++ include files that we need
40 #include <iostream>
41 #include <algorithm>
42 #include <math.h>
43 #include <set>
44 
45 // Basic include file needed for the mesh functionality.
46 #include "libmesh/libmesh.h"
47 #include "libmesh/mesh.h"
48 #include "libmesh/mesh_generation.h"
49 #include "libmesh/exodusII_io.h"
50 #include "libmesh/gnuplot_io.h"
51 #include "libmesh/linear_implicit_system.h"
52 #include "libmesh/equation_systems.h"
53 
54 // Define the Finite Element object.
55 #include "libmesh/fe.h"
56 
57 // Define Gauss quadrature rules.
58 #include "libmesh/quadrature_gauss.h"
59 
60 // Define the DofMap, which handles degree of freedom
61 // indexing.
62 #include "libmesh/dof_map.h"
63 
64 // Define useful datatypes for finite element
65 // matrix and vector components.
66 #include "libmesh/sparse_matrix.h"
67 #include "libmesh/numeric_vector.h"
68 #include "libmesh/dense_matrix.h"
69 #include "libmesh/dense_vector.h"
70 
71 // Define the PerfLog, a performance logging utility.
72 // It is useful for timing events in a code and giving
73 // you an idea where bottlenecks lie.
74 #include "libmesh/perf_log.h"
75 
76 // The definition of a geometric element
77 #include "libmesh/elem.h"
78 
79 // To impose Dirichlet boundary conditions
80 #include "libmesh/dirichlet_boundaries.h"
81 #include "libmesh/analytic_function.h"
82 
83 #include "libmesh/string_to_enum.h"
84 #include "libmesh/getpot.h"
85 #include "libmesh/enum_solver_package.h"
86 
87 // Bring in everything from the libMesh namespace
88 using namespace libMesh;
89 
90 
91 
92 // Function prototype. This is the function that will assemble
93 // the linear system for our Poisson problem. Note that the
94 // function will take the EquationSystems object and the
95 // name of the system we are assembling as input. From the
96 // EquationSystems object we have access to the Mesh and
97 // other objects we might need.
99  const std::string & system_name);
100 
101 // Exact solution function prototype.
102 Real exact_solution (const Real x,
103  const Real y,
104  const Real z = 0.);
105 
106 // Define a wrapper for exact_solution that will be needed below
108  const Point & p,
109  const Real)
110 {
111  output(0) = exact_solution(p(0),
112  (LIBMESH_DIM>1)?p(1):0,
113  (LIBMESH_DIM>2)?p(2):0);
114 }
115 
116 // Begin the main program.
117 int main (int argc, char ** argv)
118 {
119  // Initialize libMesh and any dependent libraries, like in example 2.
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  // Declare a performance log for the main program
127  // PerfLog perf_main("Main Program");
128 
129  // Check for proper calling arguments.
130  libmesh_error_msg_if(argc < 3, "Usage:\n" << "\t " << argv[0] << " -d 2(3)" << " -n 15");
131 
132  // Brief message to the user regarding the program name
133  // and command line arguments.
134  libMesh::out << "Running " << argv[0];
135 
136  for (int i=1; i<argc; i++)
137  libMesh::out << " " << argv[i];
138 
139  libMesh::out << std::endl << std::endl;
140 
141  // Read problem dimension from command line. Use int
142  // instead of unsigned since the GetPot overload is ambiguous
143  // otherwise.
144  const int dim = libMesh::command_line_next("-d", 2);
145 
146  // Skip higher-dimensional examples on a lower-dimensional libMesh build
147  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
148 
149  // We use Dirichlet boundary conditions here
150 #ifndef LIBMESH_ENABLE_DIRICHLET
151  libmesh_example_requires(false, "--enable-dirichlet");
152 #endif
153 
154  // Create a mesh with user-defined dimension.
155  // Read number of elements from command line
156  const int ps = libMesh::command_line_next("-n", 15);
157 
158  // Read FE order from command line
159  std::string order = "SECOND";
160  order = libMesh::command_line_next("-o", order);
161  order = libMesh::command_line_next("-Order", order);
162 
163  // Read FE Family from command line
164  std::string family = "LAGRANGE";
165  family = libMesh::command_line_next("-f", family);
166  family = libMesh::command_line_next("-FEFamily", family);
167 
168  // Cannot use discontinuous basis.
169  libmesh_error_msg_if((family == "MONOMIAL") || (family == "XYZ"),
170  "ex4 currently requires a C^0 (or higher) FE basis.");
171 
172  // Create a mesh, with dimension to be overridden later, distributed
173  // across the default MPI communicator.
174  Mesh mesh(init.comm());
175 
176  // Use the MeshTools::Generation mesh generator to create a uniform
177  // grid on the square [-1,1]^D. We instruct the mesh generator
178  // to build a mesh of 8x8 Quad9 elements in 2D, or Hex27
179  // elements in 3D. Building these higher-order elements allows
180  // us to use higher-order approximation, as in example 3.
181 
182  Real halfwidth = dim > 1 ? 1. : 0.;
183  Real halfheight = dim > 2 ? 1. : 0.;
184 
185  if ((family == "LAGRANGE") && (order == "FIRST"))
186  {
187  // No reason to use high-order geometric elements if we are
188  // solving with low-order finite elements.
190  ps,
191  (dim>1) ? ps : 0,
192  (dim>2) ? ps : 0,
193  -1., 1.,
194  -halfwidth, halfwidth,
195  -halfheight, halfheight,
196  (dim==1) ? EDGE2 :
197  ((dim == 2) ? QUAD4 : HEX8));
198  }
199 
200  else
201  {
203  ps,
204  (dim>1) ? ps : 0,
205  (dim>2) ? ps : 0,
206  -1., 1.,
207  -halfwidth, halfwidth,
208  -halfheight, halfheight,
209  (dim==1) ? EDGE3 :
210  ((dim == 2) ? QUAD9 : HEX27));
211  }
212 
213 
214  // Print information about the mesh to the screen.
215  mesh.print_info();
216 
217 
218  // Create an equation systems object.
219  EquationSystems equation_systems (mesh);
220 
221  // Declare the system and its variables.
222  // Create a system named "Poisson"
223  LinearImplicitSystem & system =
224  equation_systems.add_system<LinearImplicitSystem> ("Poisson");
225 
226 
227  // Add the variable "u" to "Poisson". "u"
228  // will be approximated using second-order approximation by default
229  unsigned int u_var = system.add_variable("u",
230  Utility::string_to_enum<Order> (order),
231  Utility::string_to_enum<FEFamily>(family));
232 
233  // Give the system a pointer to the matrix assembly
234  // function.
236 
237  // Construct a Dirichlet boundary condition object
238 
239  // Indicate which boundary IDs we impose the BC on
240  // We either build a line, a square or a cube, and
241  // here we indicate the boundaries IDs in each case
242  std::set<boundary_id_type> boundary_ids;
243  // the dim==1 mesh has two boundaries with IDs 0 and 1
244  boundary_ids.insert(0);
245  boundary_ids.insert(1);
246  // the dim==2 mesh has four boundaries with IDs 0, 1, 2 and 3
247  if (dim>=2)
248  {
249  boundary_ids.insert(2);
250  boundary_ids.insert(3);
251  }
252  // the dim==3 mesh has four boundaries with IDs 0, 1, 2, 3, 4 and 5
253  if (dim==3)
254  {
255  boundary_ids.insert(4);
256  boundary_ids.insert(5);
257  }
258 
259  // Create an AnalyticFunction object that we use to project the BC
260  // This function just calls the function exact_solution via exact_solution_wrapper
261  AnalyticFunction<> exact_solution_object(exact_solution_wrapper);
262 
263 #ifdef LIBMESH_ENABLE_DIRICHLET
264  // In general, when reusing a system-indexed exact solution, we want
265  // to use the default system-ordering constructor for
266  // DirichletBoundary, so we demonstrate that here. In this case,
267  // though, we have only one variable, so system- and local-
268  // orderings are the same.
269  DirichletBoundary dirichlet_bc
270  (boundary_ids, {u_var}, exact_solution_object);
271 
272  // We must add the Dirichlet boundary condition _before_
273  // we call equation_systems.init()
274  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
275 #endif
276 
277  // Initialize the data structures for the equation system.
278  equation_systems.init();
279 
280  // Print information about the system to the screen.
281  equation_systems.print_info();
282  mesh.print_info();
283 
284  // Solve the system "Poisson", just like example 2.
285  system.solve();
286 
287  // After solving the system write the solution
288  // to a GMV-formatted plot file.
289  if (dim == 1)
290  {
291  GnuPlotIO plot(mesh, "Introduction Example 4, 1D", GnuPlotIO::GRID_ON);
292  plot.write_equation_systems("gnuplot_script", equation_systems);
293  }
294 #ifdef LIBMESH_HAVE_EXODUS_API
295  else
296  {
298  "out_3.e" : "out_2.e", equation_systems);
299  }
300 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
301 
302  // All done.
303  return 0;
304 }
305 
306 
307 
308 
309 // We now define the matrix assembly function for the
310 // Poisson system. We need to first compute element
311 // matrices and right-hand sides, and then take into
312 // account the boundary conditions.
314  const std::string & libmesh_dbg_var(system_name))
315 {
316  // It is a good idea to make sure we are assembling
317  // the proper system.
318  libmesh_assert_equal_to (system_name, "Poisson");
319 
320  // Declare a performance log. Give it a descriptive
321  // string to identify what part of the code we are
322  // logging, since there may be many PerfLogs in an
323  // application.
324  PerfLog perf_log ("Matrix Assembly");
325 
326  // Get a constant reference to the mesh object.
327  const MeshBase & mesh = es.get_mesh();
328 
329  // The dimension that we are running
330  const unsigned int dim = mesh.mesh_dimension();
331 
332  // Get a reference to the LinearImplicitSystem we are solving
333  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Poisson");
334 
335  // A reference to the DofMap object for this system. The DofMap
336  // object handles the index translation from node and element numbers
337  // to degree of freedom numbers. We will talk more about the DofMap
338  // in future examples.
339  const DofMap & dof_map = system.get_dof_map();
340 
341  // Get a constant reference to the Finite Element type
342  // for the first (and only) variable in the system.
343  FEType fe_type = dof_map.variable_type(0);
344 
345  // Build a Finite Element object of the specified type. Since the
346  // FEBase::build() member dynamically creates memory we will
347  // store the object as a std::unique_ptr<FEBase>. This can be thought
348  // of as a pointer that will clean up after itself.
349  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
350 
351  // A 5th order Gauss quadrature rule for numerical integration.
352  QGauss qrule (dim, FIFTH);
353 
354  // Tell the finite element object to use our quadrature rule.
355  fe->attach_quadrature_rule (&qrule);
356 
357  // Declare a special finite element object for
358  // boundary integration.
359  std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
360 
361  // Boundary integration requires one quadrature rule,
362  // with dimensionality one less than the dimensionality
363  // of the element.
364  QGauss qface(dim-1, FIFTH);
365 
366  // Tell the finite element object to use our
367  // quadrature rule.
368  fe_face->attach_quadrature_rule (&qface);
369 
370  // Here we define some references to cell-specific data that
371  // will be used to assemble the linear system.
372  // We begin with the element Jacobian * quadrature weight at each
373  // integration point.
374  const std::vector<Real> & JxW = fe->get_JxW();
375 
376  // The physical XY locations of the quadrature points on the element.
377  // These might be useful for evaluating spatially varying material
378  // properties at the quadrature points.
379  const std::vector<Point> & q_point = fe->get_xyz();
380 
381  // The element shape functions evaluated at the quadrature points.
382  const std::vector<std::vector<Real>> & phi = fe->get_phi();
383 
384  // The element shape function gradients evaluated at the quadrature
385  // points.
386  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
387 
388  // Define data structures to contain the element matrix
389  // and right-hand-side vector contribution. Following
390  // basic finite element terminology we will denote these
391  // "Ke" and "Fe". More detail is in example 3.
394 
395  // This vector will hold the degree of freedom indices for
396  // the element. These define where in the global system
397  // the element degrees of freedom get mapped.
398  std::vector<dof_id_type> dof_indices;
399 
400  // The global system matrix
401  SparseMatrix<Number> & matrix = system.get_system_matrix();
402 
403  // Now we will loop over all the elements in the mesh.
404  // We will compute the element matrix and right-hand-side
405  // contribution. See example 3 for a discussion of the
406  // element iterators.
407  for (const auto & elem : mesh.active_local_element_ptr_range())
408  {
409  // Start logging the shape function initialization.
410  // This is done through a simple function call with
411  // the name of the event to log.
412  perf_log.push("elem init");
413 
414  // Get the degree of freedom indices for the
415  // current element. These define where in the global
416  // matrix and right-hand-side this element will
417  // contribute to.
418  dof_map.dof_indices (elem, dof_indices);
419 
420  // Cache the number of degrees of freedom on this element, for
421  // use as a loop bound later. We use cast_int to explicitly
422  // convert from size() (which may be 64-bit) to unsigned int
423  // (which may be 32-bit but which is definitely enough to count
424  // *local* degrees of freedom.
425  const unsigned int n_dofs =
426  cast_int<unsigned int>(dof_indices.size());
427 
428  // Compute the element-specific data for the current
429  // element. This involves computing the location of the
430  // quadrature points (q_point) and the shape functions
431  // (phi, dphi) for the current element.
432  fe->reinit (elem);
433 
434  // With one variable, we should have the same number of degrees
435  // of freedom as shape functions.
436  libmesh_assert_equal_to (n_dofs, phi.size());
437 
438  // Zero the element matrix and right-hand side before
439  // summing them. We use the resize member here because
440  // the number of degrees of freedom might have changed from
441  // the last element. Note that this will be the case if the
442  // element type is different (i.e. the last element was a
443  // triangle, now we are on a quadrilateral).
444  Ke.resize (n_dofs, n_dofs);
445 
446  Fe.resize (n_dofs);
447 
448  // Stop logging the shape function initialization.
449  // If you forget to stop logging an event the PerfLog
450  // object will probably catch the error and abort.
451  perf_log.pop("elem init");
452 
453  // Now we will build the element matrix. This involves
454  // a double loop to integrate the test functions (i) against
455  // the trial functions (j).
456  //
457  // We have split the numeric integration into two loops
458  // so that we can log the matrix and right-hand-side
459  // computation separately.
460  //
461  // Now start logging the element matrix computation
462  perf_log.push ("Ke");
463 
464  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
465  for (unsigned int i=0; i != n_dofs; i++)
466  for (unsigned int j=0; j != n_dofs; j++)
467  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
468 
469 
470  // Stop logging the matrix computation
471  perf_log.pop ("Ke");
472 
473  // Now we build the element right-hand-side contribution.
474  // This involves a single loop in which we integrate the
475  // "forcing function" in the PDE against the test functions.
476  //
477  // Start logging the right-hand-side computation
478  perf_log.push ("Fe");
479 
480  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
481  {
482  // fxy is the forcing function for the Poisson equation.
483  // In this case we set fxy to be a finite difference
484  // Laplacian approximation to the (known) exact solution.
485  //
486  // We will use the second-order accurate FD Laplacian
487  // approximation, which in 2D on a structured grid is
488  //
489  // u_xx + u_yy = (u(i-1,j) + u(i+1,j) +
490  // u(i,j-1) + u(i,j+1) +
491  // -4*u(i,j))/h^2
492  //
493  // Since the value of the forcing function depends only
494  // on the location of the quadrature point (q_point[qp])
495  // we will compute it here, outside of the i-loop
496  const Real x = q_point[qp](0);
497 #if LIBMESH_DIM > 1
498  const Real y = q_point[qp](1);
499 #else
500  const Real y = 0.;
501 #endif
502 #if LIBMESH_DIM > 2
503  const Real z = q_point[qp](2);
504 #else
505  const Real z = 0.;
506 #endif
507  const Real eps = 1.e-3;
508 
509  const Real uxx = (exact_solution(x-eps, y, z) +
510  exact_solution(x+eps, y, z) +
511  -2.*exact_solution(x, y, z))/eps/eps;
512 
513  const Real uyy = (exact_solution(x, y-eps, z) +
514  exact_solution(x, y+eps, z) +
515  -2.*exact_solution(x, y, z))/eps/eps;
516 
517  const Real uzz = (exact_solution(x, y, z-eps) +
518  exact_solution(x, y, z+eps) +
519  -2.*exact_solution(x, y, z))/eps/eps;
520 
521  Real fxy;
522  if (dim==1)
523  {
524  // In 1D, compute the rhs by differentiating the
525  // exact solution twice.
526  const Real pi = libMesh::pi;
527  fxy = (0.25*pi*pi)*sin(.5*pi*x);
528  }
529  else
530  {
531  fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));
532  }
533 
534  // Add the RHS contribution
535  for (unsigned int i=0; i != n_dofs; i++)
536  Fe(i) += JxW[qp]*fxy*phi[i][qp];
537  }
538 
539  // Stop logging the right-hand-side computation
540  perf_log.pop ("Fe");
541 
542  // If this assembly program were to be used on an adaptive mesh,
543  // we would have to apply any hanging node constraint equations
544  // Also, note that here we call heterogenously_constrain_element_matrix_and_vector
545  // to impose a inhomogeneous Dirichlet boundary conditions.
546  dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
547 
548  // The element matrix and right-hand-side are now built
549  // for this element. Add them to the global matrix and
550  // right-hand-side vector. The SparseMatrix::add_matrix()
551  // and NumericVector::add_vector() members do this for us.
552  // Start logging the insertion of the local (element)
553  // matrix and vector into the global matrix and vector
554  LOG_SCOPE_WITH("matrix insertion", "", perf_log);
555 
556  matrix.add_matrix (Ke, dof_indices);
557  system.rhs->add_vector (Fe, dof_indices);
558  }
559 
560  // That's it. We don't need to do anything else to the
561  // PerfLog. When it goes out of scope (at this function return)
562  // it will print its log to the screen. Pretty easy, huh?
563 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
int main(int argc, char **argv)
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
This is the EquationSystems class.
void pop(const char *label, const char *header="")
Pop the event label off the stack, resuming any lower event.
Definition: perf_log.C:168
void exact_solution_wrapper(DenseVector< Number > &output, const Point &p, const Real)
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
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
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.
virtual void solve() override
Assembles & solves the linear system A*x=b.
const T_sys & get_system(std::string_view name) const
Wraps a function pointer into a FunctionBase object.
This is the MeshBase class.
Definition: mesh_base.h:75
SolverPackage default_solver_package()
Definition: libmesh.C:1117
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:145
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
This class implements writing meshes using GNUplot, designed for use only with 1D meshes...
Definition: gnuplot_io.h:43
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
void assemble_poisson(EquationSystems &es, const std::string &system_name)
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.
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1357
void push(const char *label, const char *header="")
Push the event label onto the stack, pausing any active event.
Definition: perf_log.C:138
unsigned int n_points() const
Definition: quadrature.h:131
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:2161
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
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
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
virtual void init()
Initialize all the systems.
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
const DofMap & get_dof_map() const
Definition: system.h:2374
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const SparseMatrix< Number > & get_system_matrix() const
Real exact_solution(const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
const Real pi
.
Definition: libmesh.h:299