libMesh
Functions | Variables
introduction_ex5.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))
 

Variables

QuadratureType quad_type =INVALID_Q_RULE
 

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 217 of file introduction_ex5.C.

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

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.
243  UniquePtr<QBase> qrule(QBase::build(quad_type, dim, THIRD));
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.
275  UniquePtr<QBase> qface (QBase::build(quad_type,
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
unsigned int dim
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 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
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
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 84 of file introduction_ex5.C.

References exact_solution().

Referenced by main().

87 {
88  output(0) = exact_solution(p(0), p(1), p(2));
89 }
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 98 of file introduction_ex5.C.

References libMesh::EquationSystems::add_system(), assemble_poisson(), libMesh::MeshTools::Generation::build_cube(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), exact_solution_wrapper(), libMesh::FIRST, libMesh::EquationSystems::get_system(), libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), quad_type, and libMesh::MeshOutput< MT >::write_equation_systems().

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 }
QuadratureType
Defines an enum for currently available quadrature rules.
This is the EquationSystems class.
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
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
void exact_solution_wrapper(DenseVector< Number > &output, const Point &p, const Real)
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
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
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)

Variable Documentation

QuadratureType quad_type =INVALID_Q_RULE

Definition at line 93 of file introduction_ex5.C.

Referenced by assemble_poisson(), and main().