libMesh
Functions
introduction_ex4.C File Reference

Go to the source code of this file.

Functions

void assemble_poisson (EquationSystems &es, const std::string &system_name)
 
Real exact_solution (const Real x, const Real y, const Real z=0.)
 This is the exact solution that we are trying to obtain. More...
 
void exact_solution_wrapper (DenseVector< Number > &output, const Point &p, const Real)
 
int main (int argc, char **argv)
 
void assemble_poisson (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
 

Function Documentation

void assemble_poisson ( EquationSystems es,
const std::string &  system_name 
)

Referenced by main().

void assemble_poisson ( EquationSystems es,
const std::string &  libmesh_dbg_varsystem_name 
)

Definition at line 325 of file introduction_ex4.C.

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), dim, exact_solution(), libMesh::FIFTH, libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::pi, libMesh::PerfLog::pop(), libMesh::PerfLog::push(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, and libMesh::x.

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.
367  fe->attach_quadrature_rule (&qrule);
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
379  // quadrature rule.
380  fe_face->attach_quadrature_rule (&qface);
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 
396  // The element shape function gradients evaluated at the quadrature
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 
554  system.matrix->add_matrix (Ke, dof_indices);
555  system.rhs->add_vector (Fe, dof_indices);
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
unsigned int dim
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]...
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
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.
PetscErrorCode Vec x
const DofMap & get_dof_map() const
Definition: system.h:2030
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
const T_sys & get_system(const std::string &name) 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.
const Real pi
.
Definition: libmesh.h:172
Real exact_solution ( const Real  x,
const Real  y,
const Real  t 
)

This is the exact solution that we are trying to obtain.

We will solve

  • (u_xx + u_yy) = f

and take a finite difference approximation using this function to get f. This is the well-known "method of manufactured solutions".

Definition at line 43 of file exact_solution.C.

Referenced by assemble_poisson(), and exact_solution_wrapper().

46 {
47  static const Real pi = acos(-1.);
48 
49  return cos(.5*pi*x)*sin(.5*pi*y)*cos(.5*pi*z);
50 }
PetscErrorCode Vec x
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Real pi
.
Definition: libmesh.h:172
void exact_solution_wrapper ( DenseVector< Number > &  output,
const Point p,
const Real   
)

Definition at line 106 of file introduction_ex4.C.

References exact_solution().

Referenced by main().

109 {
110  output(0) = exact_solution(p(0),
111  (LIBMESH_DIM>1)?p(1):0,
112  (LIBMESH_DIM>2)?p(2):0);
113 }
Real exact_solution(const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.
int main ( int  argc,
char **  argv 
)

Definition at line 116 of file introduction_ex4.C.

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_poisson(), libMesh::System::attach_assemble_function(), libMesh::MeshTools::Generation::build_cube(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), dim, libMesh::EDGE2, libMesh::EDGE3, exact_solution_wrapper(), libMesh::System::get_dof_map(), libMesh::GnuPlotIO::GRID_ON, libMesh::HEX27, libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD4, libMesh::QUAD9, libMesh::Real, libMesh::LinearImplicitSystem::solve(), and libMesh::MeshOutput< MT >::write_equation_systems().

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
152  // instead of unsigned since the GetPot overload is ambiguous
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 =
233  equation_systems.add_system<LinearImplicitSystem> ("Poisson");
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()
286  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
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 }
This is the EquationSystems class.
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.
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
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
This class provides function-like objects for which an analytical expression can be provided...
SolverPackage default_solver_package()
Definition: libmesh.C:995
This class implements writing meshes using GNUplot, designed for use only with 1D meshes...
Definition: gnuplot_io.h:43
void assemble_poisson(EquationSystems &es, const std::string &system_name)
const DofMap & get_dof_map() const
Definition: system.h:2030
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 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
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
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.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
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.