libMesh
Functions
introduction_ex3.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...
 
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 178 of file introduction_ex3.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_nullptr, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, side, value, and libMesh::x.

180 {
181 
182  // It is a good idea to make sure we are assembling
183  // the proper system.
184  libmesh_assert_equal_to (system_name, "Poisson");
185 
186  // Get a constant reference to the mesh object.
187  const MeshBase & mesh = es.get_mesh();
188 
189  // The dimension that we are running
190  const unsigned int dim = mesh.mesh_dimension();
191 
192  // Get a reference to the LinearImplicitSystem we are solving
193  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> ("Poisson");
194 
195  // A reference to the DofMap object for this system. The DofMap
196  // object handles the index translation from node and element numbers
197  // to degree of freedom numbers. We will talk more about the DofMap
198  // in future examples.
199  const DofMap & dof_map = system.get_dof_map();
200 
201  // Get a constant reference to the Finite Element type
202  // for the first (and only) variable in the system.
203  FEType fe_type = dof_map.variable_type(0);
204 
205  // Build a Finite Element object of the specified type. Since the
206  // FEBase::build() member dynamically creates memory we will
207  // store the object as a UniquePtr<FEBase>. This can be thought
208  // of as a pointer that will clean up after itself. Introduction Example 4
209  // describes some advantages of UniquePtr's in the context of
210  // quadrature rules.
211  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
212 
213  // A 5th order Gauss quadrature rule for numerical integration.
214  QGauss qrule (dim, FIFTH);
215 
216  // Tell the finite element object to use our quadrature rule.
217  fe->attach_quadrature_rule (&qrule);
218 
219  // Declare a special finite element object for
220  // boundary integration.
221  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
222 
223  // Boundary integration requires one quadrature rule,
224  // with dimensionality one less than the dimensionality
225  // of the element.
226  QGauss qface(dim-1, FIFTH);
227 
228  // Tell the finite element object to use our
229  // quadrature rule.
230  fe_face->attach_quadrature_rule (&qface);
231 
232  // Here we define some references to cell-specific data that
233  // will be used to assemble the linear system.
234  //
235  // The element Jacobian * quadrature weight at each integration point.
236  const std::vector<Real> & JxW = fe->get_JxW();
237 
238  // The physical XY locations of the quadrature points on the element.
239  // These might be useful for evaluating spatially varying material
240  // properties at the quadrature points.
241  const std::vector<Point> & q_point = fe->get_xyz();
242 
243  // The element shape functions evaluated at the quadrature points.
244  const std::vector<std::vector<Real>> & phi = fe->get_phi();
245 
246  // The element shape function gradients evaluated at the quadrature
247  // points.
248  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
249 
250  // Define data structures to contain the element matrix
251  // and right-hand-side vector contribution. Following
252  // basic finite element terminology we will denote these
253  // "Ke" and "Fe". These datatypes are templated on
254  // Number, which allows the same code to work for real
255  // or complex numbers.
258 
259  // This vector will hold the degree of freedom indices for
260  // the element. These define where in the global system
261  // the element degrees of freedom get mapped.
262  std::vector<dof_id_type> dof_indices;
263 
264  // Now we will loop over all the elements in the mesh.
265  // We will compute the element matrix and right-hand-side
266  // contribution.
267  //
268  // Element iterators are a nice way to iterate through all the
269  // elements, or all the elements that have some property. The
270  // iterator el will iterate from the first to the last element on
271  // the local processor. The iterator end_el tells us when to stop.
272  // It is smart to make this one const so that we don't accidentally
273  // mess it up! In case users later modify this program to include
274  // refinement, we will be safe and will only consider the active
275  // elements; hence we use a variant of the active_elem_iterator.
276  for (const auto & elem : mesh.active_local_element_ptr_range())
277  {
278  // Get the degree of freedom indices for the
279  // current element. These define where in the global
280  // matrix and right-hand-side this element will
281  // contribute to.
282  dof_map.dof_indices (elem, dof_indices);
283 
284  // Compute the element-specific data for the current
285  // element. This involves computing the location of the
286  // quadrature points (q_point) and the shape functions
287  // (phi, dphi) for the current element.
288  fe->reinit (elem);
289 
290  // Zero the element matrix and right-hand side before
291  // summing them. We use the resize member here because
292  // the number of degrees of freedom might have changed from
293  // the last element. Note that this will be the case if the
294  // element type is different (i.e. the last element was a
295  // triangle, now we are on a quadrilateral).
296 
297  // The DenseMatrix::resize() and the DenseVector::resize()
298  // members will automatically zero out the matrix and vector.
299  Ke.resize (dof_indices.size(),
300  dof_indices.size());
301 
302  Fe.resize (dof_indices.size());
303 
304  // Now loop over the quadrature points. This handles
305  // the numeric integration.
306  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
307  {
308 
309  // Now we will build the element matrix. This involves
310  // a double loop to integrate the test functions (i) against
311  // the trial functions (j).
312  for (std::size_t i=0; i<phi.size(); i++)
313  for (std::size_t j=0; j<phi.size(); j++)
314  {
315  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
316  }
317 
318  // This is the end of the matrix summation loop
319  // Now we build the element right-hand-side contribution.
320  // This involves a single loop in which we integrate the
321  // "forcing function" in the PDE against the test functions.
322  {
323  const Real x = q_point[qp](0);
324  const Real y = q_point[qp](1);
325  const Real eps = 1.e-3;
326 
327 
328  // "fxy" is the forcing function for the Poisson equation.
329  // In this case we set fxy to be a finite difference
330  // Laplacian approximation to the (known) exact solution.
331  //
332  // We will use the second-order accurate FD Laplacian
333  // approximation, which in 2D is
334  //
335  // u_xx + u_yy = (u(i,j-1) + u(i,j+1) +
336  // u(i-1,j) + u(i+1,j) +
337  // -4*u(i,j))/h^2
338  //
339  // Since the value of the forcing function depends only
340  // on the location of the quadrature point (q_point[qp])
341  // we will compute it here, outside of the i-loop
342  const Real fxy = -(exact_solution(x, y-eps) +
343  exact_solution(x, y+eps) +
344  exact_solution(x-eps, y) +
345  exact_solution(x+eps, y) -
346  4.*exact_solution(x, y))/eps/eps;
347 
348  for (std::size_t i=0; i<phi.size(); i++)
349  Fe(i) += JxW[qp]*fxy*phi[i][qp];
350  }
351  }
352 
353  // We have now reached the end of the RHS summation,
354  // and the end of quadrature point loop, so
355  // the interior element integration has
356  // been completed. However, we have not yet addressed
357  // boundary conditions. For this example we will only
358  // consider simple Dirichlet boundary conditions.
359  //
360  // There are several ways Dirichlet boundary conditions
361  // can be imposed. A simple approach, which works for
362  // interpolary bases like the standard Lagrange polynomials,
363  // is to assign function values to the
364  // degrees of freedom living on the domain boundary. This
365  // works well for interpolary bases, but is more difficult
366  // when non-interpolary (e.g Legendre or Hierarchic) bases
367  // are used.
368  //
369  // Dirichlet boundary conditions can also be imposed with a
370  // "penalty" method. In this case essentially the L2 projection
371  // of the boundary values are added to the matrix. The
372  // projection is multiplied by some large factor so that, in
373  // floating point arithmetic, the existing (smaller) entries
374  // in the matrix and right-hand-side are effectively ignored.
375  //
376  // This amounts to adding a term of the form (in latex notation)
377  //
378  // \frac{1}{\epsilon} \int_{\delta \Omega} \phi_i \phi_j = \frac{1}{\epsilon} \int_{\delta \Omega} u \phi_i
379  //
380  // where
381  //
382  // \frac{1}{\epsilon} is the penalty parameter, defined such that \epsilon << 1
383  {
384 
385  // The following loop is over the sides of the element.
386  // If the element has no neighbor on a side then that
387  // side MUST live on a boundary of the domain.
388  for (auto side : elem->side_index_range())
389  if (elem->neighbor_ptr(side) == libmesh_nullptr)
390  {
391  // The value of the shape functions at the quadrature
392  // points.
393  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
394 
395  // The Jacobian * Quadrature Weight at the quadrature
396  // points on the face.
397  const std::vector<Real> & JxW_face = fe_face->get_JxW();
398 
399  // The XYZ locations (in physical space) of the
400  // quadrature points on the face. This is where
401  // we will interpolate the boundary value function.
402  const std::vector<Point> & qface_point = fe_face->get_xyz();
403 
404  // Compute the shape function values on the element
405  // face.
406  fe_face->reinit(elem, side);
407 
408  // Loop over the face quadrature points for integration.
409  for (unsigned int qp=0; qp<qface.n_points(); qp++)
410  {
411  // The location on the boundary of the current
412  // face quadrature point.
413  const Real xf = qface_point[qp](0);
414  const Real yf = qface_point[qp](1);
415 
416  // The penalty value. \frac{1}{\epsilon}
417  // in the discussion above.
418  const Real penalty = 1.e10;
419 
420  // The boundary value.
421  const Real value = exact_solution(xf, yf);
422 
423  // Matrix contribution of the L2 projection.
424  for (std::size_t i=0; i<phi_face.size(); i++)
425  for (std::size_t j=0; j<phi_face.size(); j++)
426  Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
427 
428  // Right-hand-side contribution of the L2
429  // projection.
430  for (std::size_t i=0; i<phi_face.size(); i++)
431  Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
432  }
433  }
434  }
435 
436  // We have now finished the quadrature point loop,
437  // and have therefore applied all the boundary conditions.
438 
439  // If this assembly program were to be used on an adaptive mesh,
440  // we would have to apply any hanging node constraint equations
441  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
442 
443  // The element matrix and right-hand-side are now built
444  // for this element. Add them to the global matrix and
445  // right-hand-side vector. The SparseMatrix::add_matrix()
446  // and NumericVector::add_vector() members do this for us.
447  system.matrix->add_matrix (Ke, dof_indices);
448  system.rhs->add_vector (Fe, dof_indices);
449  }
450 
451  // All done!
452 }
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]...
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
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
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
Real exact_solution(const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.
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
static const bool value
Definition: xdr_io.C:108
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  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().

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
int main ( int  argc,
char **  argv 
)

Definition at line 79 of file introduction_ex3.C.

References libMesh::EquationSystems::add_system(), assemble_poisson(), libMesh::MeshTools::Generation::build_square(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), libMesh::EquationSystems::get_system(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::SECOND, and libMesh::MeshOutput< MT >::write_equation_systems().

80 {
81  // Initialize libraries, like in example 2.
82  LibMeshInit init (argc, argv);
83 
84  // This example requires a linear solver package.
85  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
86  "--enable-petsc, --enable-trilinos, or --enable-eigen");
87 
88  // Brief message to the user regarding the program name
89  // and command line arguments.
90  libMesh::out << "Running " << argv[0];
91 
92  for (int i=1; i<argc; i++)
93  libMesh::out << " " << argv[i];
94 
95  libMesh::out << std::endl << std::endl;
96 
97  // Skip this 2D example if libMesh was compiled as 1D-only.
98  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
99 
100  // Create a mesh, with dimension to be overridden later, distributed
101  // across the default MPI communicator.
102  Mesh mesh(init.comm());
103 
104  // Use the MeshTools::Generation mesh generator to create a uniform
105  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
106  // to build a mesh of 15x15 QUAD9 elements. Building QUAD9
107  // elements instead of the default QUAD4's we used in example 2
108  // allow us to use higher-order approximation.
110  15, 15,
111  -1., 1.,
112  -1., 1.,
113  QUAD9);
114 
115  // Print information about the mesh to the screen.
116  // Note that 5x5 QUAD9 elements actually has 11x11 nodes,
117  // so this mesh is significantly larger than the one in example 2.
118  mesh.print_info();
119 
120  // Create an equation systems object.
121  EquationSystems equation_systems (mesh);
122 
123  // Declare the Poisson system and its variables.
124  // The Poisson system is another example of a steady system.
125  equation_systems.add_system<LinearImplicitSystem> ("Poisson");
126 
127  // Adds the variable "u" to "Poisson". "u"
128  // will be approximated using second-order approximation.
129  equation_systems.get_system("Poisson").add_variable("u", SECOND);
130 
131  // Give the system a pointer to the matrix assembly
132  // function. This will be called when needed by the
133  // library.
134  equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
135 
136  // Initialize the data structures for the equation system.
137  equation_systems.init();
138 
139  // Prints information about the system to the screen.
140  equation_systems.print_info();
141 
142  // Solve the system "Poisson". Note that calling this
143  // member will assemble the linear system and invoke
144  // the default numerical solver. With PETSc the solver can be
145  // controlled from the command line. For example,
146  // you can invoke conjugate gradient with:
147  //
148  // ./introduction_ex3 -ksp_type cg
149  //
150  // You can also get a nice X-window that monitors the solver
151  // convergence with:
152  //
153  // ./introduction-ex3 -ksp_xmonitor
154  //
155  // if you linked against the appropriate X libraries when you
156  // built PETSc.
157  equation_systems.get_system("Poisson").solve();
158 
159 #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
160 
161  // After solving the system write the solution
162  // to a VTK-formatted plot file.
163  VTKIO (mesh).write_equation_systems ("out.pvtu", equation_systems);
164 
165 #endif // #ifdef LIBMESH_HAVE_VTK
166 
167  // All done.
168  return 0;
169 }
This is the EquationSystems class.
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
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
This class implements reading and writing meshes in the VTK format.
Definition: vtk_io.h:59
SolverPackage default_solver_package()
Definition: libmesh.C:995
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
void assemble_poisson(EquationSystems &es, const std::string &system_name)
OStreamProxy out
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448