libMesh
Functions
miscellaneous_ex5.C File Reference

Go to the source code of this file.

Functions

Number exact_solution (const Point &p, const Parameters &parameters, const std::string &, const std::string &)
 
Gradient exact_derivative (const Point &p, const Parameters &parameters, const std::string &, const std::string &)
 
void assemble_ellipticdg (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
 
int main (int argc, char **argv)
 

Function Documentation

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

Definition at line 140 of file miscellaneous_ex5.C.

References libMesh::Elem::active(), libMesh::MeshBase::active_local_element_ptr_range(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, exact_solution(), libMesh::BasicOStreamProxy< charT, traits >::flush(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::QBase::get_points(), libMesh::EquationSystems::get_system(), libMesh::DofObject::id(), libMesh::FEInterface::inverse_map(), libMesh::Elem::level(), libmesh_nullptr, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::Elem::neighbor_ptr(), libMesh::out, libMesh::EquationSystems::parameters, std::pow(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, side, libMesh::TypeVector< T >::size(), libMesh::System::variable_type(), and libMesh::Elem::which_neighbor_am_i().

Referenced by main().

142 {
143  libMesh::out << " assembling elliptic dg system... ";
145 
146  // It is a good idea to make sure we are assembling
147  // the proper system.
148  libmesh_assert_equal_to (system_name, "EllipticDG");
149 
150  // Get a constant reference to the mesh object.
151  const MeshBase & mesh = es.get_mesh();
152  // The dimension that we are running
153  const unsigned int dim = mesh.mesh_dimension();
154 
155  // Get a reference to the LinearImplicitSystem we are solving
156  LinearImplicitSystem & ellipticdg_system = es.get_system<LinearImplicitSystem> ("EllipticDG");
157  // Get some parameters that we need during assembly
158  const Real penalty = es.parameters.get<Real> ("penalty");
159  std::string refinement_type = es.parameters.get<std::string> ("refinement");
160 
161  // A reference to the DofMap object for this system. The DofMap
162  // object handles the index translation from node and element numbers
163  // to degree of freedom numbers. We will talk more about the DofMap
164  const DofMap & dof_map = ellipticdg_system.get_dof_map();
165 
166  // Get a constant reference to the Finite Element type
167  // for the first (and only) variable in the system.
168  FEType fe_type = ellipticdg_system.variable_type(0);
169 
170  // Build a Finite Element object of the specified type. Since the
171  // FEBase::build() member dynamically creates memory we will
172  // store the object as a UniquePtr<FEBase>. This can be thought
173  // of as a pointer that will clean up after itself.
174  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
175  UniquePtr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
176  UniquePtr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));
177 
178  // Quadrature rules for numerical integration.
179 #ifdef QORDER
180  QGauss qrule (dim, QORDER);
181 #else
182  QGauss qrule (dim, fe_type.default_quadrature_order());
183 #endif
184  fe->attach_quadrature_rule (&qrule);
185 
186 #ifdef QORDER
187  QGauss qface(dim-1, QORDER);
188 #else
189  QGauss qface(dim-1, fe_type.default_quadrature_order());
190 #endif
191 
192  // Tell the finite element object to use our quadrature rule.
193  fe_elem_face->attach_quadrature_rule(&qface);
194  fe_neighbor_face->attach_quadrature_rule(&qface);
195 
196  // Here we define some references to cell-specific data that
197  // will be used to assemble the linear system.
198  // Data for interior volume integrals
199  const std::vector<Real> & JxW = fe->get_JxW();
200  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
201 
202  // Data for surface integrals on the element boundary
203  const std::vector<std::vector<Real>> & phi_face = fe_elem_face->get_phi();
204  const std::vector<std::vector<RealGradient>> & dphi_face = fe_elem_face->get_dphi();
205  const std::vector<Real> & JxW_face = fe_elem_face->get_JxW();
206  const std::vector<Point> & qface_normals = fe_elem_face->get_normals();
207  const std::vector<Point> & qface_points = fe_elem_face->get_xyz();
208 
209  // Data for surface integrals on the neighbor boundary
210  const std::vector<std::vector<Real>> & phi_neighbor_face = fe_neighbor_face->get_phi();
211  const std::vector<std::vector<RealGradient>> & dphi_neighbor_face = fe_neighbor_face->get_dphi();
212 
213  // Define data structures to contain the element interior matrix
214  // and right-hand-side vector contribution. Following
215  // basic finite element terminology we will denote these
216  // "Ke" and "Fe".
219 
220  // Data structures to contain the element and neighbor boundary matrix
221  // contribution. This matrices will do the coupling between the dofs of
222  // the element and those of his neighbors.
223  // Ken: matrix coupling elem and neighbor dofs
228 
229  // This vector will hold the degree of freedom indices for
230  // the element. These define where in the global system
231  // the element degrees of freedom get mapped.
232  std::vector<dof_id_type> dof_indices;
233 
234  // Now we will loop over all the elements in the mesh. We will
235  // compute first the element interior matrix and right-hand-side contribution
236  // and then the element and neighbors boundary matrix contributions.
237  for (const auto & elem : mesh.active_local_element_ptr_range())
238  {
239  // Get the degree of freedom indices for the
240  // current element. These define where in the global
241  // matrix and right-hand-side this element will
242  // contribute to.
243  dof_map.dof_indices (elem, dof_indices);
244  const unsigned int n_dofs = dof_indices.size();
245 
246  // Compute the element-specific data for the current
247  // element. This involves computing the location of the
248  // quadrature points (q_point) and the shape functions
249  // (phi, dphi) for the current element.
250  fe->reinit (elem);
251 
252  // Zero the element matrix and right-hand side before
253  // summing them. We use the resize member here because
254  // the number of degrees of freedom might have changed from
255  // the last element.
256  Ke.resize (n_dofs, n_dofs);
257  Fe.resize (n_dofs);
258 
259  // Now we will build the element interior matrix. This involves
260  // a double loop to integrate the test functions (i) against
261  // the trial functions (j).
262  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
263  for (unsigned int i=0; i<n_dofs; i++)
264  for (unsigned int j=0; j<n_dofs; j++)
265  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
266 
267  // Now we address boundary conditions.
268  // We consider Dirichlet bc imposed via the interior penalty method
269  // The following loops over the sides of the element.
270  // If the element has no neighbor on a side then that
271  // side MUST live on a boundary of the domain.
272  for (auto side : elem->side_index_range())
273  {
274  if (elem->neighbor_ptr(side) == libmesh_nullptr)
275  {
276  // Pointer to the element face
277  fe_elem_face->reinit(elem, side);
278 
279  UniquePtr<const Elem> elem_side (elem->build_side_ptr(side));
280  // h element dimension to compute the interior penalty penalty parameter
281  const unsigned int elem_b_order = static_cast<unsigned int> (fe_elem_face->get_order());
282  const double h_elem = elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);
283 
284  for (unsigned int qp=0; qp<qface.n_points(); qp++)
285  {
286  Number bc_value = exact_solution(qface_points[qp], es.parameters, "null", "void");
287  for (unsigned int i=0; i<n_dofs; i++)
288  {
289  // Matrix contribution
290  for (unsigned int j=0; j<n_dofs; j++)
291  {
292  // stability
293  Ke(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[i][qp] * phi_face[j][qp];
294 
295  // consistency
296  Ke(i,j) -=
297  JxW_face[qp] *
298  (phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) +
299  phi_face[j][qp] * (dphi_face[i][qp]*qface_normals[qp]));
300  }
301 
302  // RHS contributions
303 
304  // stability
305  Fe(i) += JxW_face[qp] * bc_value * penalty/h_elem * phi_face[i][qp];
306 
307  // consistency
308  Fe(i) -= JxW_face[qp] * dphi_face[i][qp] * (bc_value*qface_normals[qp]);
309  }
310  }
311  }
312 
313  // If the element is not on a boundary of the domain
314  // we loop over his neighbors to compute the element
315  // and neighbor boundary matrix contributions
316  else
317  {
318  // Store a pointer to the neighbor we are currently
319  // working on.
320  const Elem * neighbor = elem->neighbor_ptr(side);
321 
322  // Get the global id of the element and the neighbor
323  const unsigned int elem_id = elem->id();
324  const unsigned int neighbor_id = neighbor->id();
325 
326  // If the neighbor has the same h level and is active
327  // perform integration only if our global id is bigger than our neighbor id.
328  // We don't want to compute twice the same contributions.
329  // If the neighbor has a different h level perform integration
330  // only if the neighbor is at a lower level.
331  if ((neighbor->active() &&
332  (neighbor->level() == elem->level()) &&
333  (elem_id < neighbor_id)) ||
334  (neighbor->level() < elem->level()))
335  {
336  // Pointer to the element side
337  UniquePtr<const Elem> elem_side (elem->build_side_ptr(side));
338 
339  // h dimension to compute the interior penalty penalty parameter
340  const unsigned int elem_b_order = static_cast<unsigned int>(fe_elem_face->get_order());
341  const unsigned int neighbor_b_order = static_cast<unsigned int>(fe_neighbor_face->get_order());
342  const double side_order = (elem_b_order + neighbor_b_order)/2.;
343  const double h_elem = (elem->volume()/elem_side->volume()) * 1./pow(side_order,2.);
344 
345  // The quadrature point locations on the neighbor side
346  std::vector<Point> qface_neighbor_point;
347 
348  // The quadrature point locations on the element side
349  std::vector<Point > qface_point;
350 
351  // Reinitialize shape functions on the element side
352  fe_elem_face->reinit(elem, side);
353 
354  // Get the physical locations of the element quadrature points
355  qface_point = fe_elem_face->get_xyz();
356 
357  // Find their locations on the neighbor
358  unsigned int side_neighbor = neighbor->which_neighbor_am_i(elem);
359  if (refinement_type == "p")
360  fe_neighbor_face->side_map (neighbor,
361  elem_side.get(),
362  side_neighbor,
363  qface.get_points(),
364  qface_neighbor_point);
365  else
366  FEInterface::inverse_map (elem->dim(),
367  fe->get_fe_type(),
368  neighbor,
369  qface_point,
370  qface_neighbor_point);
371 
372  // Calculate the neighbor element shape functions at those locations
373  fe_neighbor_face->reinit(neighbor, &qface_neighbor_point);
374 
375  // Get the degree of freedom indices for the
376  // neighbor. These define where in the global
377  // matrix this neighbor will contribute to.
378  std::vector<dof_id_type> neighbor_dof_indices;
379  dof_map.dof_indices (neighbor, neighbor_dof_indices);
380  const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
381 
382  // Zero the element and neighbor side matrix before
383  // summing them. We use the resize member here because
384  // the number of degrees of freedom might have changed from
385  // the last element or neighbor.
386  // Note that Kne and Ken are not square matrices if neighbor
387  // and element have a different p level
388  Kne.resize (n_neighbor_dofs, n_dofs);
389  Ken.resize (n_dofs, n_neighbor_dofs);
390  Kee.resize (n_dofs, n_dofs);
391  Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
392 
393  // Now we will build the element and neighbor
394  // boundary matrices. This involves
395  // a double loop to integrate the test functions
396  // (i) against the trial functions (j).
397  for (unsigned int qp=0; qp<qface.n_points(); qp++)
398  {
399  // Kee Matrix. Integrate the element test function i
400  // against the element test function j
401  for (unsigned int i=0; i<n_dofs; i++)
402  {
403  for (unsigned int j=0; j<n_dofs; j++)
404  {
405  // consistency
406  Kee(i,j) -=
407  0.5 * JxW_face[qp] *
408  (phi_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) +
409  phi_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]));
410 
411  // stability
412  Kee(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_face[i][qp];
413  }
414  }
415 
416  // Knn Matrix. Integrate the neighbor test function i
417  // against the neighbor test function j
418  for (unsigned int i=0; i<n_neighbor_dofs; i++)
419  {
420  for (unsigned int j=0; j<n_neighbor_dofs; j++)
421  {
422  // consistency
423  Knn(i,j) +=
424  0.5 * JxW_face[qp] *
425  (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]) +
426  phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp]));
427 
428  // stability
429  Knn(i,j) +=
430  JxW_face[qp] * penalty/h_elem * phi_neighbor_face[j][qp]*phi_neighbor_face[i][qp];
431  }
432  }
433 
434  // Kne Matrix. Integrate the neighbor test function i
435  // against the element test function j
436  for (unsigned int i=0; i<n_neighbor_dofs; i++)
437  {
438  for (unsigned int j=0; j<n_dofs; j++)
439  {
440  // consistency
441  Kne(i,j) +=
442  0.5 * JxW_face[qp] *
443  (phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]) -
444  phi_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]));
445 
446  // stability
447  Kne(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_neighbor_face[i][qp];
448  }
449  }
450 
451  // Ken Matrix. Integrate the element test function i
452  // against the neighbor test function j
453  for (unsigned int i=0; i<n_dofs; i++)
454  {
455  for (unsigned int j=0; j<n_neighbor_dofs; j++)
456  {
457  // consistency
458  Ken(i,j) +=
459  0.5 * JxW_face[qp] *
460  (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) -
461  phi_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp]));
462 
463  // stability
464  Ken(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[i][qp]*phi_neighbor_face[j][qp];
465  }
466  }
467  }
468 
469  // The element and neighbor boundary matrix are now built
470  // for this side. Add them to the global matrix
471  // The SparseMatrix::add_matrix() members do this for us.
472  ellipticdg_system.matrix->add_matrix(Kne, neighbor_dof_indices, dof_indices);
473  ellipticdg_system.matrix->add_matrix(Ken, dof_indices, neighbor_dof_indices);
474  ellipticdg_system.matrix->add_matrix(Kee, dof_indices);
475  ellipticdg_system.matrix->add_matrix(Knn, neighbor_dof_indices);
476  }
477  }
478  }
479  // The element interior matrix and right-hand-side are now built
480  // for this element. Add them to the global matrix and
481  // right-hand-side vector. The SparseMatrix::add_matrix()
482  // and NumericVector::add_vector() members do this for us.
483  ellipticdg_system.matrix->add_matrix(Ke, dof_indices);
484  ellipticdg_system.rhs->add_vector(Fe, dof_indices);
485  }
486 
487  libMesh::out << "done" << std::endl;
488 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
bool active() const
Definition: elem.h:2257
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 int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition: elem.h:2181
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
const T & get(const std::string &) const
Definition: parameters.h:430
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
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.
const DofMap & get_dof_map() const
Definition: system.h:2030
Order default_quadrature_order() const
Definition: fe_type.h:332
Number exact_solution(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
double pow(double a, int b)
BasicOStreamProxy & flush()
Flush the associated stream buffer.
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
unsigned int level() const
Definition: elem.h:2388
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
dof_id_type id() const
Definition: dof_object.h:632
Gradient exact_derivative ( const Point p,
const Parameters parameters,
const std::string &  ,
const std::string &   
)

Definition at line 89 of file miscellaneous_ex5.C.

References libMesh::Parameters::get(), libMesh::pi, std::pow(), libMesh::Real, and libMesh::x.

Referenced by main().

93 {
94  // Gradient value to be returned.
95  Gradient gradu;
96 
97  // x and y coordinates in space
98  const Real x = p(0);
99  const Real y = p(1);
100  const Real z = p(2);
101 
102  if (parameters.get<bool>("singularity"))
103  {
104  libmesh_assert_not_equal_to (x, 0.);
105 
106  // For convenience...
107  const Real tt = 2./3.;
108  const Real ot = 1./3.;
109 
110  // The value of the radius, squared
111  const Real r2 = x*x + y*y;
112 
113  // The boundary value, given by the exact solution,
114  // u_exact = r^(2/3)*sin(2*theta/3).
115  Real theta = atan2(y, x);
116 
117  // Make sure 0 <= theta <= 2*pi
118  if (theta < 0)
119  theta += 2. * libMesh::pi;
120 
121  // du/dx
122  gradu(0) = tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x;
123  gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
124  gradu(2) = 1.;
125  }
126  else
127  {
128  gradu(0) = -sin(x) * exp(y) * (1. - z);
129  gradu(1) = cos(x) * exp(y) * (1. - z);
130  gradu(2) = -cos(x) * exp(y);
131  }
132  return gradu;
133 }
const T & get(const std::string &) const
Definition: parameters.h:430
PetscErrorCode Vec x
double pow(double a, int b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Real pi
.
Definition: libmesh.h:172
Number exact_solution ( const Point p,
const Parameters parameters,
const std::string &  ,
const std::string &   
)

Definition at line 62 of file miscellaneous_ex5.C.

References libMesh::Parameters::get(), libMesh::pi, std::pow(), libMesh::Real, and libMesh::x.

Referenced by assemble_ellipticdg(), and main().

66 {
67  const Real x = p(0);
68  const Real y = p(1);
69  const Real z = p(2);
70 
71  if (parameters.get<bool>("singularity"))
72  {
73  Real theta = atan2(y, x);
74 
75  if (theta < 0)
76  theta += 2. * libMesh::pi;
77 
78  return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
79  }
80  else
81  {
82  return cos(x) * exp(y) * (1. - z);
83  }
84 }
const T & get(const std::string &) const
Definition: parameters.h:430
PetscErrorCode Vec x
double pow(double a, int b)
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 492 of file miscellaneous_ex5.C.

References libMesh::MeshRefinement::add_p_to_h_refinement(), libMesh::EquationSystems::add_system(), libMesh::MeshTools::Modification::all_tri(), assemble_ellipticdg(), libMesh::ExactSolution::attach_exact_deriv(), libMesh::ExactSolution::attach_exact_value(), libMesh::MeshTools::Generation::build_line(), libMesh::MeshRefinement::coarsen_fraction(), libMesh::LibMeshInit::comm(), libMesh::command_line_value(), libMesh::ExactSolution::compute_error(), libMesh::default_solver_package(), dim, libMesh::JumpErrorEstimator::estimate_error(), exact_derivative(), exact_solution(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::ExactSolution::l2_error(), libMesh::MeshRefinement::max_h_level(), mesh, libMesh::MONOMIAL, libMesh::EquationSystems::n_active_dofs(), libMesh::MeshBase::n_elem(), libMesh::on_command_line(), libMesh::out, libMesh::EquationSystems::parameters, libMesh::MeshBase::read(), libMesh::Real, libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::MeshRefinement::refine_fraction(), libMesh::EquationSystems::reinit(), libMesh::Parameters::set(), singularity, libMesh::MeshRefinement::switch_h_to_p_refinement(), libMesh::TOLERANCE, libMesh::MeshRefinement::uniformly_refine(), and libMesh::ExodusII_IO::write_discontinuous_exodusII().

493 {
494  LibMeshInit init(argc, argv);
495 
496  // This example requires a linear solver package.
497  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
498  "--enable-petsc, --enable-trilinos, or --enable-eigen");
499 
500  // Skip adaptive examples on a non-adaptive libMesh build
501 #ifndef LIBMESH_ENABLE_AMR
502  libmesh_example_requires(false, "--enable-amr");
503 #else
504 
505  //Parse the input file
506  GetPot input_file("miscellaneous_ex5.in");
507 
508  //Read in parameters from the input file
509  const unsigned int adaptive_refinement_steps = input_file("max_adaptive_r_steps", 3);
510  const unsigned int uniform_refinement_steps = input_file("uniform_h_r_steps", 3);
511  const Real refine_fraction = input_file("refine_fraction", 0.5);
512  const Real coarsen_fraction = input_file("coarsen_fraction", 0.);
513  const unsigned int max_h_level = input_file("max_h_level", 10);
514  const std::string refinement_type = input_file("refinement_type","p");
515  Order p_order = static_cast<Order>(input_file("p_order", 1));
516  const std::string element_type = input_file("element_type", "tensor");
517  const Real penalty = input_file("ip_penalty", 10.);
518  const bool singularity = input_file("singularity", true);
519  const unsigned int dim = input_file("dimension", 3);
520 
521  // Skip higher-dimensional examples on a lower-dimensional libMesh build
522  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
523 
524 
525  // Create a mesh, with dimension to be overridden later, distributed
526  // across the default MPI communicator.
527  Mesh mesh(init.comm());
528 
529  if (dim == 1)
531  else if (dim == 2)
532  mesh.read("lshaped.xda");
533  else
534  mesh.read ("lshaped3D.xda");
535 
536  // Use triangles if the config file says so
537  if (element_type == "simplex")
539 
540  // Mesh Refinement object
541  MeshRefinement mesh_refinement(mesh);
542  mesh_refinement.refine_fraction() = refine_fraction;
543  mesh_refinement.coarsen_fraction() = coarsen_fraction;
544  mesh_refinement.max_h_level() = max_h_level;
545 
546  // Do uniform refinement
547  for (unsigned int rstep=0; rstep<uniform_refinement_steps; rstep++)
548  mesh_refinement.uniformly_refine(1);
549 
550  // Crate an equation system object
551  EquationSystems equation_system (mesh);
552 
553  // Set parameters for the equation system and the solver
554  equation_system.parameters.set<Real>("linear solver tolerance") = TOLERANCE * TOLERANCE;
555  equation_system.parameters.set<unsigned int>("linear solver maximum iterations") = 1000;
556  equation_system.parameters.set<Real>("penalty") = penalty;
557  equation_system.parameters.set<bool>("singularity") = singularity;
558  equation_system.parameters.set<std::string>("refinement") = refinement_type;
559 
560  // Create a system named ellipticdg
561  LinearImplicitSystem & ellipticdg_system = equation_system.add_system<LinearImplicitSystem> ("EllipticDG");
562 
563  // Add a variable "u" to "ellipticdg" using the p_order specified in the config file
564  if (on_command_line("element_type"))
565  {
566  std::string fe_str =
567  command_line_value(std::string("element_type"),
568  std::string("MONOMIAL"));
569 
570  if (fe_str != "MONOMIAL" || fe_str != "XYZ")
571  libmesh_error_msg("Error: This example must be run with MONOMIAL or XYZ element types.");
572 
573  ellipticdg_system.add_variable ("u", p_order, Utility::string_to_enum<FEFamily>(fe_str));
574  }
575  else
576  ellipticdg_system.add_variable ("u", p_order, MONOMIAL);
577 
578  // Give the system a pointer to the matrix assembly function
579  ellipticdg_system.attach_assemble_function (assemble_ellipticdg);
580 
581  // Initialize the data structures for the equation system
582  equation_system.init();
583 
584  // Construct ExactSolution object and attach solution functions
585  ExactSolution exact_sol(equation_system);
586  exact_sol.attach_exact_value(exact_solution);
587  exact_sol.attach_exact_deriv(exact_derivative);
588 
589  // A refinement loop.
590  for (unsigned int rstep=0; rstep<adaptive_refinement_steps; ++rstep)
591  {
592  libMesh::out << " Beginning Solve " << rstep << std::endl;
593  libMesh::out << "Number of elements: " << mesh.n_elem() << std::endl;
594 
595  // Solve the system
596  ellipticdg_system.solve();
597 
598  libMesh::out << "System has: "
599  << equation_system.n_active_dofs()
600  << " degrees of freedom."
601  << std::endl;
602 
603  libMesh::out << "Linear solver converged at step: "
604  << ellipticdg_system.n_linear_iterations()
605  << ", final residual: "
606  << ellipticdg_system.final_linear_residual()
607  << std::endl;
608 
609  // Compute the error
610  exact_sol.compute_error("EllipticDG", "u");
611 
612  // Print out the error values
613  libMesh::out << "L2-Error is: "
614  << exact_sol.l2_error("EllipticDG", "u")
615  << std::endl;
616 
617  // Possibly refine the mesh
618  if (rstep+1 < adaptive_refinement_steps)
619  {
620  // The ErrorVector is a particular StatisticsVector
621  // for computing error information on a finite element mesh.
622  ErrorVector error;
623 
624  // The discontinuity error estimator
625  // evaluate the jump of the solution
626  // on elements faces
627  DiscontinuityMeasure error_estimator;
628  error_estimator.estimate_error(ellipticdg_system,error);
629 
630  // Take the error in error and decide which elements will be coarsened or refined
631  mesh_refinement.flag_elements_by_error_fraction(error);
632  if (refinement_type == "p")
633  mesh_refinement.switch_h_to_p_refinement();
634  if (refinement_type == "hp")
635  mesh_refinement.add_p_to_h_refinement();
636 
637  // Refine and coarsen the flagged elements
638  mesh_refinement.refine_and_coarsen_elements();
639  equation_system.reinit();
640  }
641  }
642 
643  // Write out the solution
644  // After solving the system write the solution
645  // to a ExodusII-formatted plot file.
646 #ifdef LIBMESH_HAVE_EXODUS_API
647  ExodusII_IO (mesh).write_discontinuous_exodusII("lshaped_dg.e", equation_system);
648 #endif
649 
650 #endif // #ifndef LIBMESH_ENABLE_AMR
651 
652  // All done.
653  return 0;
654 }
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
This is the EquationSystems class.
unsigned int dim
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
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.
MeshBase & mesh
static const Real TOLERANCE
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
bool singularity
This class measures discontinuities between elements for debugging purposes.
SolverPackage default_solver_package()
Definition: libmesh.C:995
void all_tri(MeshBase &mesh)
Converts the 2D quadrilateral elements of a Mesh into triangular elements.
Gradient exact_derivative(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
This is the MeshRefinement class.
T command_line_value(const std::string &, T)
Definition: libmesh.C:932
Number exact_solution(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
OStreamProxy out
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
This function uses the derived class&#39;s jump error estimate formula to estimate the error on each cell...
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
void build_line(UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 1D meshes.
void assemble_ellipticdg(EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual dof_id_type n_elem() const =0
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
void write_discontinuous_exodusII(const std::string &name, const EquationSystems &es, const std::set< std::string > *system_names=libmesh_nullptr)
Writes a exodusII file with discontinuous data.
Definition: exodusII_io.C:89