libMesh
introduction_ex4.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
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.
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
86 // Bring in everything from the libMesh namespace
87 using namespace libMesh;
88
89
90
91 // Function prototype. This is the function that will assemble
92 // the linear system for our Poisson problem. Note that the
93 // function will take the EquationSystems object and the
94 // name of the system we are assembling as input. From the
96 // other objects we might need.
98  const std::string & system_name);
99
100 // Exact solution function prototype.
101 Real exact_solution (const Real x,
102  const Real y,
103  const Real z = 0.);
104
105 // Define a wrapper for exact_solution that will be needed below
107  const Point & p,
108  const Real)
109 {
110  output(0) = exact_solution(p(0),
111  (LIBMESH_DIM>1)?p(1):0,
112  (LIBMESH_DIM>2)?p(2):0);
113 }
114
115 // Begin the main program.
116 int main (int argc, char ** argv)
117 {
118  // Initialize libMesh and any dependent libraries, like in example 2.
119  LibMeshInit init (argc, argv);
120
121  // This example requires a linear solver package.
122  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
123  "--enable-petsc, --enable-trilinos, or --enable-eigen");
124
125  // Declare a performance log for the main program
126  // PerfLog perf_main("Main Program");
127
128  // Create a GetPot object to parse the command line
129  GetPot command_line (argc, argv);
130
131  // Check for proper calling arguments.
132  if (argc < 3)
133  {
134  // This handy function will print the file name, line number,
135  // specified message, and then throw an exception.
136  libmesh_error_msg("Usage:\n" << "\t " << argv[0] << " -d 2(3)" << " -n 15");
137  }
138
139  // Brief message to the user regarding the program name
140  // and command line arguments.
141  else
142  {
143  libMesh::out << "Running " << argv[0];
144
145  for (int i=1; i<argc; i++)
146  libMesh::out << " " << argv[i];
147
148  libMesh::out << std::endl << std::endl;
149  }
150
151  // Read problem dimension from command line. Use int
153  // otherwise.
154  int dim = 2;
155  if (command_line.search(1, "-d"))
156  dim = command_line.next(dim);
157
158  // Skip higher-dimensional examples on a lower-dimensional libMesh build
159  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
160
161  // Create a mesh with user-defined dimension.
162  // Read number of elements from command line
163  int ps = 15;
164  if (command_line.search(1, "-n"))
165  ps = command_line.next(ps);
166
167  // Read FE order from command line
168  std::string order = "SECOND";
169  if (command_line.search(2, "-Order", "-o"))
170  order = command_line.next(order);
171
172  // Read FE Family from command line
173  std::string family = "LAGRANGE";
174  if (command_line.search(2, "-FEFamily", "-f"))
175  family = command_line.next(family);
176
177  // Cannot use discontinuous basis.
178  if ((family == "MONOMIAL") || (family == "XYZ"))
179  libmesh_error_msg("ex4 currently requires a C^0 (or higher) FE basis.");
180
181  // Create a mesh, with dimension to be overridden later, distributed
182  // across the default MPI communicator.
183  Mesh mesh(init.comm());
184
185  // Use the MeshTools::Generation mesh generator to create a uniform
186  // grid on the square [-1,1]^D. We instruct the mesh generator
187  // to build a mesh of 8x8 Quad9 elements in 2D, or Hex27
188  // elements in 3D. Building these higher-order elements allows
189  // us to use higher-order approximation, as in example 3.
190
191  Real halfwidth = dim > 1 ? 1. : 0.;
192  Real halfheight = dim > 2 ? 1. : 0.;
193
194  if ((family == "LAGRANGE") && (order == "FIRST"))
195  {
196  // No reason to use high-order geometric elements if we are
197  // solving with low-order finite elements.
199  ps,
200  (dim>1) ? ps : 0,
201  (dim>2) ? ps : 0,
202  -1., 1.,
203  -halfwidth, halfwidth,
204  -halfheight, halfheight,
205  (dim==1) ? EDGE2 :
206  ((dim == 2) ? QUAD4 : HEX8));
207  }
208
209  else
210  {
212  ps,
213  (dim>1) ? ps : 0,
214  (dim>2) ? ps : 0,
215  -1., 1.,
216  -halfwidth, halfwidth,
217  -halfheight, halfheight,
218  (dim==1) ? EDGE3 :
219  ((dim == 2) ? QUAD9 : HEX27));
220  }
221
222
223  // Print information about the mesh to the screen.
224  mesh.print_info();
225
226
227  // Create an equation systems object.
228  EquationSystems equation_systems (mesh);
229
230  // Declare the system and its variables.
231  // Create a system named "Poisson"
232  LinearImplicitSystem & system =
234
235
236  // Add the variable "u" to "Poisson". "u"
237  // will be approximated using second-order approximation.
238  unsigned int u_var = system.add_variable("u",
239  Utility::string_to_enum<Order> (order),
240  Utility::string_to_enum<FEFamily>(family));
241
242  // Give the system a pointer to the matrix assembly
243  // function.
245
246  // Construct a Dirichlet boundary condition object
247
248  // Indicate which boundary IDs we impose the BC on
249  // We either build a line, a square or a cube, and
250  // here we indicate the boundaries IDs in each case
251  std::set<boundary_id_type> boundary_ids;
252  // the dim==1 mesh has two boundaries with IDs 0 and 1
253  boundary_ids.insert(0);
254  boundary_ids.insert(1);
255  // the dim==2 mesh has four boundaries with IDs 0, 1, 2 and 3
256  if (dim>=2)
257  {
258  boundary_ids.insert(2);
259  boundary_ids.insert(3);
260  }
261  // the dim==3 mesh has four boundaries with IDs 0, 1, 2, 3, 4 and 5
262  if (dim==3)
263  {
264  boundary_ids.insert(4);
265  boundary_ids.insert(5);
266  }
267
268  // Create a vector storing the variable numbers which the BC applies to
269  std::vector<unsigned int> variables(1);
270  variables[0] = u_var;
271
272  // Create an AnalyticFunction object that we use to project the BC
273  // This function just calls the function exact_solution via exact_solution_wrapper
274  AnalyticFunction<> exact_solution_object(exact_solution_wrapper);
275
276  // In general, when reusing a system-indexed exact solution, we want
277  // to use the default system-ordering constructor for
278  // DirichletBoundary, so we demonstrate that here. In this case,
279  // though, we have only one variable, so system- and local-
280  // orderings are the same.
281  DirichletBoundary dirichlet_bc
282  (boundary_ids, variables, exact_solution_object);
283
284  // We must add the Dirichlet boundary condition _before_
285  // we call equation_systems.init()
287
288
289  // Initialize the data structures for the equation system.
290  equation_systems.init();
291
292  // Print information about the system to the screen.
293  equation_systems.print_info();
294  mesh.print_info();
295
296  // Solve the system "Poisson", just like example 2.
297  system.solve();
298
299  // After solving the system write the solution
300  // to a GMV-formatted plot file.
301  if (dim == 1)
302  {
303  GnuPlotIO plot(mesh, "Introduction Example 4, 1D", GnuPlotIO::GRID_ON);
304  plot.write_equation_systems("gnuplot_script", equation_systems);
305  }
306 #ifdef LIBMESH_HAVE_EXODUS_API
307  else
308  {
309  ExodusII_IO (mesh).write_equation_systems ((dim == 3) ?
310  "out_3.e" : "out_2.e", equation_systems);
311  }
312 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
313
314  // All done.
315  return 0;
316 }
317
318
319
320
321 // We now define the matrix assembly function for the
322 // Poisson system. We need to first compute element
323 // matrices and right-hand sides, and then take into
324 // account the boundary conditions.
326  const std::string & libmesh_dbg_var(system_name))
327 {
328  // It is a good idea to make sure we are assembling
329  // the proper system.
330  libmesh_assert_equal_to (system_name, "Poisson");
331
332  // Declare a performance log. Give it a descriptive
333  // string to identify what part of the code we are
334  // logging, since there may be many PerfLogs in an
335  // application.
336  PerfLog perf_log ("Matrix Assembly");
337
338  // Get a constant reference to the mesh object.
339  const MeshBase & mesh = es.get_mesh();
340
341  // The dimension that we are running
342  const unsigned int dim = mesh.mesh_dimension();
343
344  // Get a reference to the LinearImplicitSystem we are solving
345  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Poisson");
346
347  // A reference to the DofMap object for this system. The DofMap
348  // object handles the index translation from node and element numbers
349  // to degree of freedom numbers. We will talk more about the DofMap
350  // in future examples.
351  const DofMap & dof_map = system.get_dof_map();
352
353  // Get a constant reference to the Finite Element type
354  // for the first (and only) variable in the system.
355  FEType fe_type = dof_map.variable_type(0);
356
357  // Build a Finite Element object of the specified type. Since the
358  // FEBase::build() member dynamically creates memory we will
359  // store the object as a UniquePtr<FEBase>. This can be thought
360  // of as a pointer that will clean up after itself.
361  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
362
363  // A 5th order Gauss quadrature rule for numerical integration.
364  QGauss qrule (dim, FIFTH);
365
366  // Tell the finite element object to use our quadrature rule.
368
369  // Declare a special finite element object for
370  // boundary integration.
371  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
372
373  // Boundary integration requires one quadrature rule,
374  // with dimensionality one less than the dimensionality
375  // of the element.
376  QGauss qface(dim-1, FIFTH);
377
378  // Tell the finite element object to use our
381
382  // Here we define some references to cell-specific data that
383  // will be used to assemble the linear system.
384  // We begin with the element Jacobian * quadrature weight at each
385  // integration point.
386  const std::vector<Real> & JxW = fe->get_JxW();
387
388  // The physical XY locations of the quadrature points on the element.
389  // These might be useful for evaluating spatially varying material
390  // properties at the quadrature points.
391  const std::vector<Point> & q_point = fe->get_xyz();
392
393  // The element shape functions evaluated at the quadrature points.
394  const std::vector<std::vector<Real>> & phi = fe->get_phi();
395
397  // points.
398  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
399
400  // Define data structures to contain the element matrix
401  // and right-hand-side vector contribution. Following
402  // basic finite element terminology we will denote these
403  // "Ke" and "Fe". More detail is in example 3.
406
407  // This vector will hold the degree of freedom indices for
408  // the element. These define where in the global system
409  // the element degrees of freedom get mapped.
410  std::vector<dof_id_type> dof_indices;
411
412  // Now we will loop over all the elements in the mesh.
413  // We will compute the element matrix and right-hand-side
414  // contribution. See example 3 for a discussion of the
415  // element iterators.
416  for (const auto & elem : mesh.active_local_element_ptr_range())
417  {
418  // Start logging the shape function initialization.
419  // This is done through a simple function call with
420  // the name of the event to log.
421  perf_log.push("elem init");
422
423  // Get the degree of freedom indices for the
424  // current element. These define where in the global
425  // matrix and right-hand-side this element will
426  // contribute to.
427  dof_map.dof_indices (elem, dof_indices);
428
429  // Compute the element-specific data for the current
430  // element. This involves computing the location of the
431  // quadrature points (q_point) and the shape functions
432  // (phi, dphi) for the current element.
433  fe->reinit (elem);
434
435  // Zero the element matrix and right-hand side before
436  // summing them. We use the resize member here because
437  // the number of degrees of freedom might have changed from
438  // the last element. Note that this will be the case if the
439  // element type is different (i.e. the last element was a
440  // triangle, now we are on a quadrilateral).
441  Ke.resize (dof_indices.size(),
442  dof_indices.size());
443
444  Fe.resize (dof_indices.size());
445
446  // Stop logging the shape function initialization.
447  // If you forget to stop logging an event the PerfLog
448  // object will probably catch the error and abort.
449  perf_log.pop("elem init");
450
451  // Now we will build the element matrix. This involves
452  // a double loop to integrate the test functions (i) against
453  // the trial functions (j).
454  //
455  // We have split the numeric integration into two loops
456  // so that we can log the matrix and right-hand-side
457  // computation separately.
458  //
459  // Now start logging the element matrix computation
460  perf_log.push ("Ke");
461
462  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
463  for (std::size_t i=0; i<phi.size(); i++)
464  for (std::size_t j=0; j<phi.size(); j++)
465  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
466
467
468  // Stop logging the matrix computation
469  perf_log.pop ("Ke");
470
471  // Now we build the element right-hand-side contribution.
472  // This involves a single loop in which we integrate the
473  // "forcing function" in the PDE against the test functions.
474  //
475  // Start logging the right-hand-side computation
476  perf_log.push ("Fe");
477
478  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
479  {
480  // fxy is the forcing function for the Poisson equation.
481  // In this case we set fxy to be a finite difference
482  // Laplacian approximation to the (known) exact solution.
483  //
484  // We will use the second-order accurate FD Laplacian
485  // approximation, which in 2D on a structured grid is
486  //
487  // u_xx + u_yy = (u(i-1,j) + u(i+1,j) +
488  // u(i,j-1) + u(i,j+1) +
489  // -4*u(i,j))/h^2
490  //
491  // Since the value of the forcing function depends only
492  // on the location of the quadrature point (q_point[qp])
493  // we will compute it here, outside of the i-loop
494  const Real x = q_point[qp](0);
495 #if LIBMESH_DIM > 1
496  const Real y = q_point[qp](1);
497 #else
498  const Real y = 0.;
499 #endif
500 #if LIBMESH_DIM > 2
501  const Real z = q_point[qp](2);
502 #else
503  const Real z = 0.;
504 #endif
505  const Real eps = 1.e-3;
506
507  const Real uxx = (exact_solution(x-eps, y, z) +
508  exact_solution(x+eps, y, z) +
509  -2.*exact_solution(x, y, z))/eps/eps;
510
511  const Real uyy = (exact_solution(x, y-eps, z) +
512  exact_solution(x, y+eps, z) +
513  -2.*exact_solution(x, y, z))/eps/eps;
514
515  const Real uzz = (exact_solution(x, y, z-eps) +
516  exact_solution(x, y, z+eps) +
517  -2.*exact_solution(x, y, z))/eps/eps;
518
519  Real fxy;
520  if (dim==1)
521  {
522  // In 1D, compute the rhs by differentiating the
523  // exact solution twice.
524  const Real pi = libMesh::pi;
525  fxy = (0.25*pi*pi)*sin(.5*pi*x);
526  }
527  else
528  {
529  fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));
530  }
531
532  // Add the RHS contribution
533  for (std::size_t i=0; i<phi.size(); i++)
534  Fe(i) += JxW[qp]*fxy*phi[i][qp];
535  }
536
537  // Stop logging the right-hand-side computation
538  perf_log.pop ("Fe");
539
540  // If this assembly program were to be used on an adaptive mesh,
541  // we would have to apply any hanging node constraint equations
542  // Also, note that here we call heterogenously_constrain_element_matrix_and_vector
543  // to impose a inhomogeneous Dirichlet boundary conditions.
544  dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
545
546  // The element matrix and right-hand-side are now built
547  // for this element. Add them to the global matrix and
548  // right-hand-side vector. The SparseMatrix::add_matrix()
549  // and NumericVector::add_vector() members do this for us.
550  // Start logging the insertion of the local (element)
551  // matrix and vector into the global matrix and vector
552  perf_log.push ("matrix insertion");
553
556
557  // Start logging the insertion of the local (element)
558  // matrix and vector into the global matrix and vector
559  perf_log.pop ("matrix insertion");
560  }
561
562  // That's it. We don't need to do anything else to the
563  // PerfLog. When it goes out of scope (at this function return)
564  // it will print its log to the screen. Pretty easy, huh?
565 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
int main(int argc, char **argv)
This is the EquationSystems class.
void pop(const char *label, const char *header="")
Pop the event label off the stack, resuming any lower event.
Definition: perf_log.C:162
void exact_solution_wrapper(DenseVector< Number > &output, const Point &p, const Real)
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
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:62
The libMesh namespace provides an interface to certain functionality in the library.
This class provides function-like objects for which an analytical expression can be provided...
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
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:125
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
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)
PetscErrorCode Vec x
const DofMap & get_dof_map() const
Definition: system.h:2030
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:1795
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
void push(const char *label, const char *header="")
Push the event label onto the stack, pausing any active event.
Definition: perf_log.C:132
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
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.
const MeshBase & get_mesh() const
Adds a copy of the specified Dirichlet boundary to the system.
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
unsigned int n_points() const