libMesh
Functions | Variables
adaptivity_ex3.C File Reference

Go to the source code of this file.

Functions

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

Variables

unsigned int dim = 2
 
bool singularity = true
 

Function Documentation

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

Definition at line 618 of file adaptivity_ex3.C.

References 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_rule(), exact_solution(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::libmesh_ignore(), libmesh_nullptr, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::EquationSystems::parameters, libMesh::PerfLog::pop(), std::pow(), libMesh::PerfLog::push(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, singularity, value, and libMesh::x.

Referenced by main().

620 {
621  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
622  libmesh_ignore(es);
623  libmesh_ignore(system_name);
624 
625 #ifdef LIBMESH_ENABLE_AMR
626  // It is a good idea to make sure we are assembling
627  // the proper system.
628  libmesh_assert_equal_to (system_name, "Laplace");
629 
630  // Declare a performance log. Give it a descriptive
631  // string to identify what part of the code we are
632  // logging, since there may be many PerfLogs in an
633  // application.
634  PerfLog perf_log ("Matrix Assembly", false);
635 
636  // Get a constant reference to the mesh object.
637  const MeshBase & mesh = es.get_mesh();
638 
639  // The dimension that we are running
640  const unsigned int mesh_dim = mesh.mesh_dimension();
641 
642  // Get a reference to the LinearImplicitSystem we are solving
643  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Laplace");
644 
645  // A reference to the DofMap object for this system. The DofMap
646  // object handles the index translation from node and element numbers
647  // to degree of freedom numbers. We will talk more about the DofMap
648  // in future examples.
649  const DofMap & dof_map = system.get_dof_map();
650 
651  // Get a constant reference to the Finite Element type
652  // for the first (and only) variable in the system.
653  FEType fe_type = dof_map.variable_type(0);
654 
655  // Build a Finite Element object of the specified type. Since the
656  // FEBase::build() member dynamically creates memory we will
657  // store the object as a UniquePtr<FEBase>. This can be thought
658  // of as a pointer that will clean up after itself.
659  UniquePtr<FEBase> fe (FEBase::build(mesh_dim, fe_type));
660  UniquePtr<FEBase> fe_face (FEBase::build(mesh_dim, fe_type));
661 
662  // Quadrature rules for numerical integration.
663  UniquePtr<QBase> qrule(fe_type.default_quadrature_rule(mesh_dim));
664  UniquePtr<QBase> qface(fe_type.default_quadrature_rule(mesh_dim-1));
665 
666  // Tell the finite element object to use our quadrature rule.
667  fe->attach_quadrature_rule (qrule.get());
668  fe_face->attach_quadrature_rule (qface.get());
669 
670  // Here we define some references to cell-specific data that
671  // will be used to assemble the linear system.
672  // We begin with the element Jacobian * quadrature weight at each
673  // integration point.
674  const std::vector<Real> & JxW = fe->get_JxW();
675  const std::vector<Real> & JxW_face = fe_face->get_JxW();
676 
677  // The physical XY locations of the quadrature points on the element.
678  // These might be useful for evaluating spatially varying material
679  // properties or forcing functions at the quadrature points.
680  const std::vector<Point> & q_point = fe->get_xyz();
681 
682  // The element shape functions evaluated at the quadrature points.
683  // For this simple problem we usually only need them on element
684  // boundaries.
685  const std::vector<std::vector<Real>> & phi = fe->get_phi();
686  const std::vector<std::vector<Real>> & psi = fe_face->get_phi();
687 
688  // The element shape function gradients evaluated at the quadrature
689  // points.
690  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
691 
692  // The XY locations of the quadrature points used for face integration
693  const std::vector<Point> & qface_points = fe_face->get_xyz();
694 
695  // Define data structures to contain the element matrix
696  // and right-hand-side vector contribution. Following
697  // basic finite element terminology we will denote these
698  // "Ke" and "Fe". More detail is in example 3.
701 
702  // This vector will hold the degree of freedom indices for
703  // the element. These define where in the global system
704  // the element degrees of freedom get mapped.
705  std::vector<dof_id_type> dof_indices;
706 
707  // Now we will loop over all the elements in the mesh. We will
708  // compute the element matrix and right-hand-side contribution. See
709  // example 3 for a discussion of the element iterators. Here we use
710  // the const_active_local_elem_iterator to indicate we only want
711  // to loop over elements that are assigned to the local processor
712  // which are "active" in the sense of AMR. This allows each
713  // processor to compute its components of the global matrix for
714  // active elements while ignoring parent elements which have been
715  // refined.
716  for (const auto & elem : mesh.active_local_element_ptr_range())
717  {
718  // Start logging the shape function initialization.
719  // This is done through a simple function call with
720  // the name of the event to log.
721  perf_log.push("elem init");
722 
723  // Get the degree of freedom indices for the
724  // current element. These define where in the global
725  // matrix and right-hand-side this element will
726  // contribute to.
727  dof_map.dof_indices (elem, dof_indices);
728 
729  // Compute the element-specific data for the current
730  // element. This involves computing the location of the
731  // quadrature points (q_point) and the shape functions
732  // (phi, dphi) for the current element.
733  fe->reinit (elem);
734 
735  // Zero the element matrix and right-hand side before
736  // summing them. We use the resize member here because
737  // the number of degrees of freedom might have changed from
738  // the last element. Note that this will be the case if the
739  // element type is different (i.e. the last element was a
740  // triangle, now we are on a quadrilateral).
741  Ke.resize (dof_indices.size(),
742  dof_indices.size());
743 
744  Fe.resize (dof_indices.size());
745 
746  // Stop logging the shape function initialization.
747  // If you forget to stop logging an event the PerfLog
748  // object will probably catch the error and abort.
749  perf_log.pop("elem init");
750 
751  // Now we will build the element matrix. This involves
752  // a double loop to integrate the test functions (i) against
753  // the trial functions (j).
754  //
755  // Now start logging the element matrix computation
756  perf_log.push ("Ke");
757 
758  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
759  for (std::size_t i=0; i<dphi.size(); i++)
760  for (std::size_t j=0; j<dphi.size(); j++)
761  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
762 
763  // We need a forcing function to make the 1D case interesting
764  if (mesh_dim == 1)
765  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
766  {
767  Real x = q_point[qp](0);
768  Real f = singularity ? sqrt(3.)/9.*pow(-x, -4./3.) :
769  cos(x);
770  for (std::size_t i=0; i<dphi.size(); ++i)
771  Fe(i) += JxW[qp]*phi[i][qp]*f;
772  }
773 
774  // Stop logging the matrix computation
775  perf_log.pop ("Ke");
776 
777 
778  // At this point the interior element integration has
779  // been completed. However, we have not yet addressed
780  // boundary conditions. For this example we will only
781  // consider simple Dirichlet boundary conditions imposed
782  // via the penalty method.
783  //
784  // This approach adds the L2 projection of the boundary
785  // data in penalty form to the weak statement. This is
786  // a more generic approach for applying Dirichlet BCs
787  // which is applicable to non-Lagrange finite element
788  // discretizations.
789  {
790  // Start logging the boundary condition computation
791  perf_log.push ("BCs");
792 
793  // The penalty value.
794  const Real penalty = 1.e10;
795 
796  // The following loops over the sides of the element.
797  // If the element has no neighbor on a side then that
798  // side MUST live on a boundary of the domain.
799  for (auto s : elem->side_index_range())
800  if (elem->neighbor_ptr(s) == libmesh_nullptr)
801  {
802  fe_face->reinit(elem, s);
803 
804  for (unsigned int qp=0; qp<qface->n_points(); qp++)
805  {
806  const Number value = exact_solution (qface_points[qp],
807  es.parameters,
808  "null",
809  "void");
810 
811  // RHS contribution
812  for (std::size_t i=0; i<psi.size(); i++)
813  Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
814 
815  // Matrix contribution
816  for (std::size_t i=0; i<psi.size(); i++)
817  for (std::size_t j=0; j<psi.size(); j++)
818  Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
819  }
820  }
821 
822  // Stop logging the boundary condition computation
823  perf_log.pop ("BCs");
824  }
825 
826 
827  // The element matrix and right-hand-side are now built
828  // for this element. Add them to the global matrix and
829  // right-hand-side vector. The SparseMatrix::add_matrix()
830  // and NumericVector::add_vector() members do this for us.
831  // Start logging the insertion of the local (element)
832  // matrix and vector into the global matrix and vector
833  perf_log.push ("matrix insertion");
834 
835  dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
836  system.matrix->add_matrix (Ke, dof_indices);
837  system.rhs->add_vector (Fe, dof_indices);
838 
839  // Start logging the insertion of the local (element)
840  // matrix and vector into the global matrix and vector
841  perf_log.pop ("matrix insertion");
842  }
843 
844  // That's it. We don't need to do anything else to the
845  // PerfLog. When it goes out of scope (at this function return)
846  // it will print its log to the screen. Pretty easy, huh?
847 #endif // #ifdef LIBMESH_ENABLE_AMR
848 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
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
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
bool singularity
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
Number exact_solution(const Point &p, const Parameters &, const std::string &, const std::string &)
double pow(double a, int b)
void libmesh_ignore(const T &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
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
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
Gradient exact_derivative ( const Point p,
const Parameters ,
const std::string &  ,
const std::string &   
)

Definition at line 552 of file adaptivity_ex3.C.

References dim, libMesh::pi, std::pow(), libMesh::Real, singularity, and libMesh::x.

Referenced by main().

556 {
557  // Gradient value to be returned.
558  Gradient gradu;
559 
560  // x and y coordinates in space
561  const Real x = p(0);
562  const Real y = dim > 1 ? p(1) : 0.;
563 
564  if (singularity)
565  {
566  // We can't compute the gradient at x=0, it is not defined.
567  libmesh_assert_not_equal_to (x, 0.);
568 
569  // For convenience...
570  const Real tt = 2./3.;
571  const Real ot = 1./3.;
572 
573  // The value of the radius, squared
574  const Real r2 = x*x + y*y;
575 
576  // The boundary value, given by the exact solution,
577  // u_exact = r^(2/3)*sin(2*theta/3).
578  Real theta = atan2(y, x);
579 
580  // Make sure 0 <= theta <= 2*pi
581  if (theta < 0)
582  theta += 2. * libMesh::pi;
583 
584  // du/dx
585  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;
586 
587  // du/dy
588  if (dim > 1)
589  gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
590 
591  if (dim > 2)
592  gradu(2) = 1.;
593  }
594  else
595  {
596  const Real z = (dim > 2) ? p(2) : 0;
597 
598  gradu(0) = -sin(x) * exp(y) * (1. - z);
599  if (dim > 1)
600  gradu(1) = cos(x) * exp(y) * (1. - z);
601  if (dim > 2)
602  gradu(2) = -cos(x) * exp(y);
603  }
604 
605  return gradu;
606 }
unsigned int dim
bool singularity
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 ,
const std::string &  ,
const std::string &   
)

Definition at line 512 of file adaptivity_ex3.C.

References dim, libMesh::pi, std::pow(), libMesh::Real, singularity, and libMesh::x.

Referenced by assemble_laplace(), and main().

516 {
517  const Real x = p(0);
518  const Real y = (dim > 1) ? p(1) : 0.;
519 
520  if (singularity)
521  {
522  // The exact solution to the singular problem,
523  // u_exact = r^(2/3)*sin(2*theta/3).
524  Real theta = atan2(y, x);
525 
526  // Make sure 0 <= theta <= 2*pi
527  if (theta < 0)
528  theta += 2. * libMesh::pi;
529 
530  // Make the 3D solution similar
531  const Real z = (dim > 2) ? p(2) : 0;
532 
533  return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
534  }
535  else
536  {
537  // The exact solution to a nonsingular problem,
538  // good for testing ideal convergence rates
539  const Real z = (dim > 2) ? p(2) : 0;
540 
541  return cos(x) * exp(y) * (1. - z);
542  }
543 }
unsigned int dim
bool singularity
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 116 of file adaptivity_ex3.C.

References libMesh::DofMap::add_algebraic_ghosting_functor(), libMesh::MeshRefinement::add_p_to_h_refinement(), libMesh::System::add_variable(), libMesh::MeshBase::all_second_order(), libMesh::MeshTools::Modification::all_tri(), assemble_laplace(), libMesh::System::attach_assemble_function(), libMesh::ExactSolution::attach_exact_deriv(), libMesh::ExactErrorEstimator::attach_exact_deriv(), libMesh::ExactSolution::attach_exact_value(), libMesh::ExactErrorEstimator::attach_exact_value(), libMesh::MeshTools::Generation::build_line(), libMesh::MeshRefinement::coarsen_fraction(), libMesh::LibMeshInit::comm(), libMesh::ExactSolution::compute_error(), libMesh::default_solver_package(), dim, libMesh::JumpErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::UniformRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), exact_derivative(), exact_solution(), libMesh::ExactSolution::extra_quadrature_order(), libMesh::LinearImplicitSystem::final_linear_residual(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::System::get_dof_map(), libMesh::ExactSolution::h1_error(), libMesh::TriangleWrapper::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::ExactSolution::l2_error(), libMesh::libmesh_assert(), libMesh::MeshRefinement::max_h_level(), libMesh::ErrorVector::mean(), mesh, libMesh::LinearImplicitSystem::n_linear_iterations(), libMesh::out, libMesh::ErrorVector::plot_error(), std::pow(), libMesh::MeshBase::read(), libMesh::Real, libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::MeshRefinement::refine_fraction(), libMesh::HPSingularity::select_refinement(), libMesh::HPCoarsenTest::select_refinement(), libMesh::HPSingularity::singular_points, singularity, libMesh::LinearImplicitSystem::solve(), libMesh::MeshRefinement::switch_h_to_p_refinement(), libMesh::TOLERANCE, libMesh::MeshRefinement::uniformly_p_refine(), libMesh::MeshRefinement::uniformly_refine(), and libMesh::MeshOutput< MT >::write_equation_systems().

117 {
118  // Initialize libMesh.
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  // Single precision is inadequate for p refinement
126  libmesh_example_requires(sizeof(Real) > 4, "--disable-singleprecision");
127 
128  // Skip adaptive examples on a non-adaptive libMesh build
129 #ifndef LIBMESH_ENABLE_AMR
130  libmesh_example_requires(false, "--enable-amr");
131 #else
132 
133  // Parse the input file
134  GetPot input_file("adaptivity_ex3.in");
135 
136  // But allow the command line to override it.
137  input_file.parse_command_line(argc, argv);
138 
139  // Read in parameters from the input file
140  const unsigned int max_r_steps = input_file("max_r_steps", 3);
141  const unsigned int max_r_level = input_file("max_r_level", 3);
142  const Real refine_percentage = input_file("refine_percentage", 0.5);
143  const Real coarsen_percentage = input_file("coarsen_percentage", 0.5);
144  const unsigned int uniform_refine = input_file("uniform_refine",0);
145  const std::string refine_type = input_file("refinement_type", "h");
146  const std::string approx_type = input_file("approx_type", "LAGRANGE");
147  const unsigned int approx_order = input_file("approx_order", 1);
148  const std::string element_type = input_file("element_type", "tensor");
149  const int extra_error_quadrature = input_file("extra_error_quadrature", 0);
150  const int max_linear_iterations = input_file("max_linear_iterations", 5000);
151 
152 #ifdef LIBMESH_HAVE_EXODUS_API
153  const bool output_intermediate = input_file("output_intermediate", false);
154 #endif
155 
156  // If libmesh is configured without second derivative support, we
157  // can't run this example with Hermite elements and will therefore
158  // fail gracefully.
159 #if !defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
160  libmesh_example_requires(approx_type != "HERMITE", "--enable-second");
161 #endif
162 
163  dim = input_file("dimension", 2);
164  const std::string indicator_type = input_file("indicator_type", "kelly");
165  singularity = input_file("singularity", true);
166 
167  // Skip higher-dimensional examples on a lower-dimensional libMesh build
168  libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
169 
170  // Output file for plotting the error as a function of
171  // the number of degrees of freedom.
172  std::string approx_name = "";
173  if (element_type == "tensor")
174  approx_name += "bi";
175  if (approx_order == 1)
176  approx_name += "linear";
177  else if (approx_order == 2)
178  approx_name += "quadratic";
179  else if (approx_order == 3)
180  approx_name += "cubic";
181  else if (approx_order == 4)
182  approx_name += "quartic";
183 
184  std::string output_file = approx_name;
185  output_file += "_";
186  output_file += refine_type;
187  if (uniform_refine == 0)
188  output_file += "_adaptive.m";
189  else
190  output_file += "_uniform.m";
191 
192  std::ofstream out (output_file.c_str());
193  out << "% dofs L2-error H1-error" << std::endl;
194  out << "e = [" << std::endl;
195 
196  // Create a mesh, with dimension to be overridden later, on the
197  // default MPI communicator.
198  Mesh mesh(init.comm());
199 
200  // Create an equation systems object.
201  EquationSystems equation_systems (mesh);
202 
203  // Declare the system we will solve, named Laplace
204  LinearImplicitSystem & system =
205  equation_systems.add_system<LinearImplicitSystem> ("Laplace");
206 
207  // If we're doing HP refinement, then we'll need to be able to
208  // evaluate data on elements' siblings in HPCoarsenTest, which means
209  // we should instruct our system's DofMap to distribute that data.
210  // Create and add (and keep around! the DofMap will only be holding
211  // a pointer to it!) a SiblingCoupling functor to provide that
212  // instruction.
213  SiblingCoupling sibling_coupling;
214  if (refine_type == "hp")
215  system.get_dof_map().add_algebraic_ghosting_functor(sibling_coupling);
216 
217  // Read in the mesh
218  if (dim == 1)
220  else if (dim == 2)
221  mesh.read("lshaped.xda");
222  else
223  mesh.read("lshaped3D.xda");
224 
225  // Use triangles if the config file says so
226  if (element_type == "simplex")
228 
229  // We used first order elements to describe the geometry,
230  // but we may need second order elements to hold the degrees
231  // of freedom
232  if (approx_order > 1 || refine_type != "h")
234 
235  // Mesh Refinement object
236  MeshRefinement mesh_refinement(mesh);
237  mesh_refinement.refine_fraction() = refine_percentage;
238  mesh_refinement.coarsen_fraction() = coarsen_percentage;
239  mesh_refinement.max_h_level() = max_r_level;
240 
241  // Adds the variable "u" to "Laplace", using
242  // the finite element type and order specified
243  // in the config file
244  system.add_variable("u", static_cast<Order>(approx_order),
245  Utility::string_to_enum<FEFamily>(approx_type));
246 
247  // Give the system a pointer to the matrix assembly
248  // function.
250 
251  // Initialize the data structures for the equation system.
252  equation_systems.init();
253 
254  // Set linear solver max iterations
255  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations")
256  = max_linear_iterations;
257 
258  // Linear solver tolerance.
259  equation_systems.parameters.set<Real>("linear solver tolerance") =
260  std::pow(TOLERANCE, 2.5);
261 
262  // Prints information about the system to the screen.
263  equation_systems.print_info();
264 
265  // Construct ExactSolution object and attach solution functions
266  ExactSolution exact_sol(equation_systems);
267  exact_sol.attach_exact_value(exact_solution);
268  exact_sol.attach_exact_deriv(exact_derivative);
269 
270  // Use higher quadrature order for more accurate error results
271  exact_sol.extra_quadrature_order(extra_error_quadrature);
272 
273  // A refinement loop.
274  for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
275  {
276  libMesh::out << "Beginning Solve " << r_step << std::endl;
277 
278  // Solve the system "Laplace", just like example 2.
279  system.solve();
280 
281  libMesh::out << "System has: "
282  << equation_systems.n_active_dofs()
283  << " degrees of freedom."
284  << std::endl;
285 
286  libMesh::out << "Linear solver converged at step: "
287  << system.n_linear_iterations()
288  << ", final residual: "
289  << system.final_linear_residual()
290  << std::endl;
291 
292 #ifdef LIBMESH_HAVE_EXODUS_API
293  // After solving the system write the solution
294  // to a ExodusII-formatted plot file.
295  if (output_intermediate)
296  {
297  std::ostringstream outfile;
298  outfile << "lshaped_" << r_step << ".e";
299  ExodusII_IO (mesh).write_equation_systems (outfile.str(),
300  equation_systems);
301  }
302 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
303 
304  // Compute the error.
305  exact_sol.compute_error("Laplace", "u");
306 
307  // Print out the error values
308  libMesh::out << "L2-Error is: "
309  << exact_sol.l2_error("Laplace", "u")
310  << std::endl;
311  libMesh::out << "H1-Error is: "
312  << exact_sol.h1_error("Laplace", "u")
313  << std::endl;
314 
315  // Compute any discontinuity. There should be none.
316  {
318  ErrorVector disc_error;
319  disc.estimate_error(system, disc_error);
320 
321  Real mean_disc_error = disc_error.mean();
322 
323  libMesh::out << "Mean discontinuity error = " << mean_disc_error << std::endl;
324 
325  // FIXME - this test fails when solving with Eigen?
326 #ifdef LIBMESH_ENABLE_PETSC
327  libmesh_assert_less (mean_disc_error, 1e-14);
328 #endif
329  }
330 
331  // Print to output file
332  out << equation_systems.n_active_dofs() << " "
333  << exact_sol.l2_error("Laplace", "u") << " "
334  << exact_sol.h1_error("Laplace", "u") << std::endl;
335 
336  // Possibly refine the mesh
337  if (r_step+1 != max_r_steps)
338  {
339  libMesh::out << " Refining the mesh..." << std::endl;
340 
341  if (uniform_refine == 0)
342  {
343 
344  // The ErrorVector is a particular StatisticsVector
345  // for computing error information on a finite element mesh.
346  ErrorVector error;
347 
348  if (indicator_type == "exact")
349  {
350  // The ErrorEstimator class interrogates a
351  // finite element solution and assigns to each
352  // element a positive error value.
353  // This value is used for deciding which elements to
354  // refine and which to coarsen.
355  // For these simple test problems, we can use
356  // numerical quadrature of the exact error between
357  // the approximate and analytic solutions.
358  // However, for real problems, we would need an error
359  // indicator which only relies on the approximate
360  // solution.
361  ExactErrorEstimator error_estimator;
362 
363  error_estimator.attach_exact_value(exact_solution);
364  error_estimator.attach_exact_deriv(exact_derivative);
365 
366  // We optimize in H1 norm, the default
367  // error_estimator.error_norm = H1;
368 
369  // Compute the error for each active element using
370  // the provided indicator. Note in general you
371  // will need to provide an error estimator
372  // specifically designed for your application.
373  error_estimator.estimate_error (system, error);
374  }
375  else if (indicator_type == "patch")
376  {
377  // The patch recovery estimator should give a
378  // good estimate of the solution interpolation
379  // error.
380  PatchRecoveryErrorEstimator error_estimator;
381 
382  error_estimator.estimate_error (system, error);
383  }
384  else if (indicator_type == "uniform")
385  {
386  // Error indication based on uniform refinement
387  // is reliable, but very expensive.
388  UniformRefinementEstimator error_estimator;
389 
390  error_estimator.estimate_error (system, error);
391  }
392  else
393  {
394  libmesh_assert_equal_to (indicator_type, "kelly");
395 
396  // The Kelly error estimator is based on
397  // an error bound for the Poisson problem
398  // on linear elements, but is useful for
399  // driving adaptive refinement in many problems
400  KellyErrorEstimator error_estimator;
401 
402  error_estimator.estimate_error (system, error);
403  }
404 
405  // Write out the error distribution
406  std::ostringstream ss;
407  ss << r_step;
408 #ifdef LIBMESH_HAVE_EXODUS_API
409  std::string error_output = "error_" + ss.str() + ".e";
410 #else
411  std::string error_output = "error_" + ss.str() + ".gmv";
412 #endif
413  error.plot_error(error_output, mesh);
414 
415  // This takes the error in error and decides which elements
416  // will be coarsened or refined. Any element within 20% of the
417  // maximum error on any element will be refined, and any
418  // element within 10% of the minimum error on any element might
419  // be coarsened. Note that the elements flagged for refinement
420  // will be refined, but those flagged for coarsening _might_ be
421  // coarsened.
422  mesh_refinement.flag_elements_by_error_fraction (error);
423 
424  // If we are doing adaptive p refinement, we want
425  // elements flagged for that instead.
426  if (refine_type == "p")
427  mesh_refinement.switch_h_to_p_refinement();
428  // If we are doing "matched hp" refinement, we
429  // flag elements for both h and p
430  else if (refine_type == "matchedhp")
431  mesh_refinement.add_p_to_h_refinement();
432  // If we are doing hp refinement, we
433  // try switching some elements from h to p
434  else if (refine_type == "hp")
435  {
436  HPCoarsenTest hpselector;
437  hpselector.select_refinement(system);
438  }
439  // If we are doing "singular hp" refinement, we
440  // try switching most elements from h to p
441  else if (refine_type == "singularhp")
442  {
443  // This only differs from p refinement for
444  // the singular problem
446  HPSingularity hpselector;
447  // Our only singular point is at the origin
448  hpselector.singular_points.push_back(Point());
449  hpselector.select_refinement(system);
450  }
451  else if (refine_type != "h")
452  libmesh_error_msg("Unknown refinement_type = " << refine_type);
453 
454  // This call actually refines and coarsens the flagged
455  // elements.
456  mesh_refinement.refine_and_coarsen_elements();
457  }
458 
459  else if (uniform_refine == 1)
460  {
461  if (refine_type == "h" || refine_type == "hp" ||
462  refine_type == "matchedhp")
463  mesh_refinement.uniformly_refine(1);
464  if (refine_type == "p" || refine_type == "hp" ||
465  refine_type == "matchedhp")
466  mesh_refinement.uniformly_p_refine(1);
467  }
468 
469  // This call reinitializes the EquationSystems object for
470  // the newly refined mesh. One of the steps in the
471  // reinitialization is projecting the solution,
472  // old_solution, etc... vectors from the old mesh to
473  // the current one.
474  equation_systems.reinit ();
475  }
476  }
477 
478 #ifdef LIBMESH_HAVE_EXODUS_API
479  // Write out the solution
480  // After solving the system write the solution
481  // to a ExodusII-formatted plot file.
482  ExodusII_IO (mesh).write_equation_systems ("lshaped.e",
483  equation_systems);
484 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
485 
486  // Close up the output file.
487  out << "];" << std::endl;
488  out << "hold on" << std::endl;
489  out << "plot(e(:,1), e(:,2), 'bo-');" << std::endl;
490  out << "plot(e(:,1), e(:,3), 'ro-');" << std::endl;
491  // out << "set(gca,'XScale', 'Log');" << std::endl;
492  // out << "set(gca,'YScale', 'Log');" << std::endl;
493  out << "xlabel('dofs');" << std::endl;
494  out << "title('" << approx_name << " elements');" << std::endl;
495  out << "legend('L2-error', 'H1-error');" << std::endl;
496  // out << "disp('L2-error linear fit');" << std::endl;
497  // out << "polyfit(log10(e(:,1)), log10(e(:,2)), 1)" << std::endl;
498  // out << "disp('H1-error linear fit');" << std::endl;
499  // out << "polyfit(log10(e(:,1)), log10(e(:,3)), 1)" << std::endl;
500 #endif // #ifndef LIBMESH_ENABLE_AMR
501 
502  // All done.
503  return 0;
504 }
unsigned int n_linear_iterations() const
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
This is the EquationSystems class.
void assemble_laplace(EquationSystems &es, const std::string &system_name)
virtual void select_refinement(System &system) libmesh_override
This pure virtual function must be redefined in derived classes to take a mesh flagged for h refineme...
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 does uniform refinements and a solve to get an improved solution on each cell...
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.
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
static const Real TOLERANCE
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
void attach_exact_value(unsigned int sys_num, FunctionBase< Number > *f)
Clone and attach an arbitrary functor which computes the exact value of the system sys_num solution a...
bool singularity
This class measures discontinuities between elements for debugging purposes.
libmesh_assert(j)
std::list< Point > singular_points
This list, to be filled by the user, should include all singular points in the solution.
Definition: hp_singular.h:79
void add_algebraic_ghosting_functor(GhostingFunctor &ghosting_functor)
Adds a functor which can specify algebraic ghosting requirements for use with distributed vectors...
Definition: dof_map.C:1813
SolverPackage default_solver_package()
Definition: libmesh.C:995
void all_tri(MeshBase &mesh)
Converts the 2D quadrilateral elements of a Mesh into triangular elements.
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
virtual void all_second_order(const bool full_ordered=true)=0
Converts a (conforming, non-refined) mesh with linear elements into a mesh with second-order elements...
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
This class implements a ``brute force&#39;&#39; error estimator which integrates differences between the curr...
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
Number exact_solution(const Point &p, const Parameters &, const std::string &, const std::string &)
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 Patch Recovery error estimate to estimate the error on each cell...
void plot_error(const std::string &filename, const MeshBase &mesh) const
Plots a data file, of a type determined by looking at the file extension in filename, of the error values on the active elements of mesh.
Definition: error_vector.C:208
double pow(double a, int b)
virtual void select_refinement(System &system)
This pure virtual function must be redefined in derived classes to take a mesh flagged for h refineme...
Definition: hp_singular.C:36
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
Gradient exact_derivative(const Point &p, const Parameters &, const std::string &, const std::string &)
This class implements the Kelly error indicator which is based on the flux jumps between elements...
This class uses the error estimate given by different types of derefinement (h coarsening or p reduct...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This class adds coupling (for use in send_list construction) between active elements and all descenda...
This class implements an "error estimator" based on the difference between the approximate and exact ...
This class implements the Patch Recovery error indicator.
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...
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 exact solution function to estimate the error on each cell.
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.
This class uses a user-provided list of singularity locations to choose between h refining and p elev...
Definition: hp_singular.h:49
void attach_exact_deriv(unsigned int sys_num, FunctionBase< Gradient > *g)
Clone and attach an arbitrary functor which computes the exact gradient of the system sys_num solutio...
virtual Real mean() const libmesh_override
Definition: error_vector.C:67
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
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.
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38

Variable Documentation

unsigned int dim = 2

Definition at line 110 of file adaptivity_ex3.C.

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::FEMContext::_do_elem_position_set(), libMesh::DofMap::_dof_indices(), libMesh::UniformRefinementEstimator::_estimate_error(), assemble(), LinearElasticity::assemble(), assemble_1D(), AssembleOptimization::assemble_A_and_F(), assemble_biharmonic(), assemble_cd(), assemble_elasticity(), assemble_ellipticdg(), assemble_helmholtz(), assemble_mass(), assemble_matrices(), assemble_poisson(), assemble_shell(), assemble_stokes(), assemble_wave(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::EquationSystems::build_variable_names(), libMesh::System::calculate_norm(), libMesh::JumpErrorEstimator::coarse_n_flux_faces_increment(), libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::FEXYZ< Dim >::compute_face_values(), libMesh::FEGenericBase< OutputType >::compute_shape_functions(), libMesh::FEXYZ< Dim >::compute_shape_functions(), compute_stresses(), LinearElasticityWithContact::compute_stresses(), LinearElasticity::compute_stresses(), LargeDeformationElasticity::compute_stresses(), libMesh::MeshFunction::discontinuous_gradient(), libMesh::MeshFunction::discontinuous_value(), libMesh::FEMContext::elem_fe_reinit(), libMesh::FEMContext::elem_position_get(), NavierSystem::element_constraint(), NavierSystem::element_time_derivative(), SolidSystem::element_time_derivative(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), exact_derivative(), exact_solution(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), NavierSystem::forcing(), libMesh::FEMContext::get_element_fe(), libMesh::FEMContext::get_element_qrule(), libMesh::FEMContext::get_side_fe(), libMesh::FEMContext::get_side_qrule(), libMesh::MeshFunction::gradient(), libMesh::MeshFunction::hessian(), libMesh::LaplacianErrorEstimator::init_context(), libMesh::DiscontinuityMeasure::init_context(), libMesh::KellyErrorEstimator::init_context(), libMesh::OldSolutionValue< Output, point_output >::init_context(), NavierSystem::init_data(), SolidSystem::init_data(), libMesh::FEMContext::init_internal_data(), libMesh::RBEIMAssembly::initialize_fe(), libMesh::TreeNode< N >::insert(), libMesh::LaplacianErrorEstimator::internal_side_integration(), LaplaceYoung::jacobian(), LargeDeformationElasticity::jacobian(), main(), NavierSystem::mass_residual(), libMesh::DofMap::old_dof_indices(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::MeshFunction::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::FEType::operator<(), NavierSystem::postprocess(), libMesh::FE< Dim, T >::reinit(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), libMesh::DofMap::reinit(), libMesh::JumpErrorEstimator::reinit_sides(), LaplaceYoung::residual(), LargeDeformationElasticity::residual(), LinearElasticityWithContact::residual_and_jacobian(), SolidSystem::save_initial_mesh(), libMesh::HPCoarsenTest::select_refinement(), ElemTest< elem_type >::setUp(), setup(), libMesh::Side< SideType, ParentType >::Side(), NavierSystem::side_constraint(), libMesh::FEMContext::side_fe_reinit(), libMesh::SideEdge< EdgeType, ParentType >::SideEdge(), SlitMeshRefinedSystemTest::testRestart(), SlitMeshRefinedSystemTest::testSystem(), libMesh::PatchRecoveryErrorEstimator::type(), libMesh::libmesh_final< T >::type(), libMesh::ExodusII_IO_Helper::write_as_dimension(), libMesh::EnsightIO::write_scalar_ascii(), libMesh::EnsightIO::write_vector_ascii(), libMesh::FEInterface::~FEInterface(), libMesh::FEMap::~FEMap(), libMesh::FETransformationBase< OutputShape >::~FETransformationBase(), libMesh::FEXYZMap::~FEXYZMap(), libMesh::H1FETransformation< OutputShape >::~H1FETransformation(), libMesh::HCurlFETransformation< OutputShape >::~HCurlFETransformation(), and libMesh::QBase::~QBase().

bool singularity = true

Definition at line 113 of file adaptivity_ex3.C.

Referenced by assemble_laplace(), exact_derivative(), exact_solution(), and main().