libMesh
introduction_ex5.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Introduction Example 5 - Run-Time Quadrature Rule Selection</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This is the fifth example program. It builds on
25 // the previous two examples, and extends the use
26 // of the UniquePtr as a convenient build method to
27 // determine the quadrature rule at run time.
28 
29 
30 // C++ include files that we need
31 #include <iostream>
32 #include <sstream>
33 #include <algorithm>
34 #include <math.h>
35 
36 // Basic include file needed for the mesh functionality.
37 #include "libmesh/libmesh.h"
38 #include "libmesh/mesh.h"
39 #include "libmesh/mesh_generation.h"
40 #include "libmesh/exodusII_io.h"
41 #include "libmesh/linear_implicit_system.h"
42 #include "libmesh/equation_systems.h"
43 
44 // Define the Finite Element object.
45 #include "libmesh/fe.h"
46 
47 // Define the base quadrature class, with which
48 // specialized quadrature rules will be built.
49 #include "libmesh/quadrature.h"
50 
51 // Define useful datatypes for finite element
52 // matrix and vector components.
53 #include "libmesh/sparse_matrix.h"
54 #include "libmesh/numeric_vector.h"
55 #include "libmesh/dense_matrix.h"
56 #include "libmesh/dense_vector.h"
57 
58 // Define the DofMap, which handles degree of freedom
59 // indexing.
60 #include "libmesh/dof_map.h"
61 
62 // To impose Dirichlet boundary conditions
63 #include "libmesh/dirichlet_boundaries.h"
64 #include "libmesh/analytic_function.h"
65 
66 // The definition of a geometric element
67 #include "libmesh/elem.h"
68 
69 // Bring in everything from the libMesh namespace
70 using namespace libMesh;
71 
72 
73 
74 // Function prototype, as before.
76  const std::string & system_name);
77 
78 // Exact solution function prototype, as before.
79 Real exact_solution (const Real x,
80  const Real y,
81  const Real z = 0.);
82 
83 // Define a wrapper for exact_solution that will be needed below
85  const Point & p,
86  const Real)
87 {
88  output(0) = exact_solution(p(0), p(1), p(2));
89 }
90 
91 
92 // The quadrature type the user requests.
94 
95 
96 
97 // Begin the main program.
98 int main (int argc, char ** argv)
99 {
100  // Initialize libMesh and any dependent libraries, like in example 2.
101  LibMeshInit init (argc, argv);
102 
103  // This example requires a linear solver package.
104  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
105  "--enable-petsc, --enable-trilinos, or --enable-eigen");
106 
107  // Check for proper usage. The quadrature rule
108  // must be given at run time.
109  if (argc < 3)
110  {
111  libmesh_error_msg("Usage: " << argv[0] << " -q <rule>\n" \
112  << " where <rule> is one of QGAUSS, QSIMPSON, or QTRAP.");
113  }
114 
115 
116  // Tell the user what we are doing.
117  else
118  {
119  libMesh::out << "Running " << argv[0];
120 
121  for (int i=1; i<argc; i++)
122  libMesh::out << " " << argv[i];
123 
124  libMesh::out << std::endl << std::endl;
125  }
126 
127 
128  // Set the quadrature rule type that the user wants from argv[2]
129  quad_type = static_cast<QuadratureType>(std::atoi(argv[2]));
130 
131  // Skip this 3D example if libMesh was compiled as 1D-only.
132  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");
133 
134  // The following is identical to example 4, and therefore
135  // not commented. Differences are mentioned when present.
136  Mesh mesh(init.comm());
137 
138  // We will use a linear approximation space in this example,
139  // hence 8-noded hexahedral elements are sufficient. This
140  // is different than example 4 where we used 27-noded
141  // hexahedral elements to support a second-order approximation
142  // space.
144  16, 16, 16,
145  -1., 1.,
146  -1., 1.,
147  -1., 1.,
148  HEX8);
149 
150  mesh.print_info();
151 
152  EquationSystems equation_systems (mesh);
153 
154  equation_systems.add_system<LinearImplicitSystem> ("Poisson");
155 
156  unsigned int u_var = equation_systems.get_system("Poisson").add_variable("u", FIRST);
157 
158  equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
159 
160  // Construct a Dirichlet boundary condition object
161 
162  // Indicate which boundary IDs we impose the BC on
163  // We either build a line, a square or a cube, and
164  // here we indicate the boundaries IDs in each case
165  std::set<boundary_id_type> boundary_ids;
166  // the dim==1 mesh has two boundaries with IDs 0 and 1
167  boundary_ids.insert(0);
168  boundary_ids.insert(1);
169  boundary_ids.insert(2);
170  boundary_ids.insert(3);
171  boundary_ids.insert(4);
172  boundary_ids.insert(5);
173 
174  // Create a vector storing the variable numbers which the BC applies to
175  std::vector<unsigned int> variables(1);
176  variables[0] = u_var;
177 
178  // Create an AnalyticFunction object that we use to project the BC
179  // This function just calls the function exact_solution via exact_solution_wrapper
180  AnalyticFunction<> exact_solution_object(exact_solution_wrapper);
181 
182  // In general, when reusing a system-indexed exact solution, we want
183  // to use the default system-ordering constructor for
184  // DirichletBoundary, so we demonstrate that here. In this case,
185  // though, we have only one variable, so system- and local-
186  // orderings are the same.
187  DirichletBoundary dirichlet_bc
188  (boundary_ids, variables, exact_solution_object);
189 
190  // We must add the Dirichlet boundary condition _before_
191  // we call equation_systems.init()
192  equation_systems.get_system("Poisson").get_dof_map().add_dirichlet_boundary(dirichlet_bc);
193 
194  equation_systems.init();
195 
196  equation_systems.print_info();
197 
198  equation_systems.get_system("Poisson").solve();
199 
200  // "Personalize" the output, with the
201  // number of the quadrature rule appended.
202  std::ostringstream f_name;
203  f_name << "out_" << quad_type << ".e";
204 
205 #ifdef LIBMESH_HAVE_EXODUS_API
206  ExodusII_IO(mesh).write_equation_systems (f_name.str(),
207  equation_systems);
208 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
209 
210  // All done.
211  return 0;
212 }
213 
214 
215 
216 
218  const std::string & libmesh_dbg_var(system_name))
219 {
220  libmesh_assert_equal_to (system_name, "Poisson");
221 
222  const MeshBase & mesh = es.get_mesh();
223 
224  const unsigned int dim = mesh.mesh_dimension();
225 
226  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Poisson");
227 
228  const DofMap & dof_map = system.get_dof_map();
229 
230  FEType fe_type = dof_map.variable_type(0);
231 
232  // Build a Finite Element object of the specified type. Since the
233  // FEBase::build() member dynamically creates memory we will
234  // store the object as a UniquePtr<FEBase>. Below, the
235  // functionality of UniquePtr's is described more detailed in
236  // the context of building quadrature rules.
237  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
238 
239  // Now this deviates from example 4. we create a
240  // 5th order quadrature rule of user-specified type
241  // for numerical integration. Note that not all
242  // quadrature rules support this order.
244 
245  // Tell the finite element object to use our
246  // quadrature rule. Note that a UniquePtr<QBase> returns
247  // a QBase* pointer to the object it handles with get().
248  // However, using get(), the UniquePtr<QBase> qrule is
249  // still in charge of this pointer. I.e., when qrule goes
250  // out of scope, it will safely delete the QBase object it
251  // points to. This behavior may be overridden using
252  // UniquePtr<Xyz>::release(), but is currently not
253  // recommended.
254  fe->attach_quadrature_rule (qrule.get());
255 
256  // Declare a special finite element object for
257  // boundary integration.
258  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
259 
260  // As already seen in example 3, boundary integration
261  // requires a quadrature rule. Here, however,
262  // we use the more convenient way of building this
263  // rule at run-time using quad_type. Note that one
264  // could also have initialized the face quadrature rules
265  // with the type directly determined from qrule, namely
266  // through:
267  // \verbatim
268  // UniquePtr<QBase> qface (QBase::build(qrule->type(),
269  // dim-1,
270  // THIRD));
271  // \endverbatim
272  // And again: using the UniquePtr<QBase> relaxes
273  // the need to delete the object afterward,
274  // they clean up themselves.
276  dim-1,
277  THIRD));
278 
279  // Tell the finite element object to use our
280  // quadrature rule. Note that a UniquePtr<QBase> returns
281  // a QBase* pointer to the object it handles with get().
282  // However, using get(), the UniquePtr<QBase> qface is
283  // still in charge of this pointer. I.e., when qface goes
284  // out of scope, it will safely delete the QBase object it
285  // points to. This behavior may be overridden using
286  // UniquePtr<Xyz>::release(), but is not recommended.
287  fe_face->attach_quadrature_rule (qface.get());
288 
289  // This is again identical to example 4, and not commented.
290  const std::vector<Real> & JxW = fe->get_JxW();
291 
292  const std::vector<Point> & q_point = fe->get_xyz();
293 
294  const std::vector<std::vector<Real>> & phi = fe->get_phi();
295 
296  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
297 
300  std::vector<dof_id_type> dof_indices;
301 
302  // Now we will loop over all the elements in the mesh.
303  // See example 3 for details.
304  for (const auto & elem : mesh.active_local_element_ptr_range())
305  {
306  dof_map.dof_indices (elem, dof_indices);
307 
308  fe->reinit (elem);
309 
310  Ke.resize (dof_indices.size(),
311  dof_indices.size());
312 
313  Fe.resize (dof_indices.size());
314 
315  // Now loop over the quadrature points. This handles
316  // the numeric integration. Note the slightly different
317  // access to the QBase members!
318  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
319  {
320  // Add the matrix contribution
321  for (std::size_t i=0; i<phi.size(); i++)
322  for (std::size_t j=0; j<phi.size(); j++)
323  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
324 
325  // fxy is the forcing function for the Poisson equation.
326  // In this case we set fxy to be a finite difference
327  // Laplacian approximation to the (known) exact solution.
328  //
329  // We will use the second-order accurate FD Laplacian
330  // approximation, which in 2D on a structured grid is
331  //
332  // u_xx + u_yy = (u(i-1,j) + u(i+1,j) +
333  // u(i,j-1) + u(i,j+1) +
334  // -4*u(i,j))/h^2
335  //
336  // Since the value of the forcing function depends only
337  // on the location of the quadrature point (q_point[qp])
338  // we will compute it here, outside of the i-loop
339  const Real x = q_point[qp](0);
340  const Real y = q_point[qp](1);
341  const Real z = q_point[qp](2);
342  const Real eps = 1.e-3;
343 
344  const Real uxx = (exact_solution(x-eps, y, z) +
345  exact_solution(x+eps, y, z) +
346  -2.*exact_solution(x, y, z))/eps/eps;
347 
348  const Real uyy = (exact_solution(x, y-eps, z) +
349  exact_solution(x, y+eps, z) +
350  -2.*exact_solution(x, y, z))/eps/eps;
351 
352  const Real uzz = (exact_solution(x, y, z-eps) +
353  exact_solution(x, y, z+eps) +
354  -2.*exact_solution(x, y, z))/eps/eps;
355 
356  const Real fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));
357 
358 
359  // Add the RHS contribution
360  for (std::size_t i=0; i<phi.size(); i++)
361  Fe(i) += JxW[qp]*fxy*phi[i][qp];
362  }
363 
364  // If this assembly program were to be used on an adaptive mesh,
365  // we would have to apply any hanging node constraint equations
366  // Call heterogenously_constrain_element_matrix_and_vector to impose
367  // non-homogeneous Dirichlet BCs
368  dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
369 
370  // The element matrix and right-hand-side are now built
371  // for this element. Add them to the global matrix and
372  // right-hand-side vector. The SparseMatrix::add_matrix()
373  // and NumericVector::add_vector() members do this for us.
374  system.matrix->add_matrix (Ke, dof_indices);
375  system.rhs->add_vector (Fe, dof_indices);
376 
377  } // end of element loop
378 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
QuadratureType
Defines an enum for currently available quadrature rules.
This is the EquationSystems class.
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.
QuadratureType quad_type
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 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
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
void exact_solution_wrapper(DenseVector< Number > &output, const Point &p, const Real)
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
int main(int argc, char **argv)
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
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 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.
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
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
static UniquePtr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
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 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.
void assemble_poisson(EquationSystems &es, const std::string &system_name)
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.