DiscontinuityMeasure Class Reference

#include <discontinuity_measure.h>

Inheritance diagram for DiscontinuityMeasure:

[legend]

List of all members.

Public Types
typedef std::map< std::pair < const System *, unsigned int > , ErrorVector * >	ErrorMap
Public Member Functions
	DiscontinuityMeasure ()
	~DiscontinuityMeasure ()
void	attach_essential_bc_function (std::pair< bool, Real > fptr(const System &system, const Point &p, const std::string &var_name))
virtual void	estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=NULL, bool estimate_parent_error=false)
virtual void	estimate_errors (const EquationSystems &equation_systems, ErrorVector &error_per_cell, const std::map< const System *, SystemNorm > &error_norms, const std::map< const System *, const NumericVector< Number > * > *solution_vectors=NULL, bool estimate_parent_error=false)
virtual void	estimate_errors (const EquationSystems &equation_systems, ErrorMap &errors_per_cell, const std::map< const System *, const NumericVector< Number > * > *solution_vectors=NULL, bool estimate_parent_error=false)
Public Attributes
bool	scale_by_n_flux_faces
SystemNorm	error_norm
Protected Member Functions
virtual void	initialize (const System &system, ErrorVector &error_per_cell, bool estimate_parent_error)
virtual void	internal_side_integration ()
virtual bool	boundary_side_integration ()
void	reinit_sides ()
float	coarse_n_flux_faces_increment ()
void	reduce_error (std::vector< float > &error_per_cell) const
Protected Attributes
const System *	my_system
std::pair< bool, Real >(*	_bc_function )(const System &system, const Point &p, const std::string &var_name)
bool	integrate_boundary_sides
const Elem *	fine_elem
const Elem *	coarse_elem
Real	fine_error
Real	coarse_error
unsigned int	fine_side
unsigned int	var
DenseVector< Number >	Ufine
DenseVector< Number >	Ucoarse
AutoPtr< FEBase >	fe_fine
AutoPtr< FEBase >	fe_coarse

---

Detailed Description

This class measures discontinuities between elements for debugging purposes. It derives from ErrorEstimator just in case someone finds it useful in a DG framework.

Author:

Roy H. Stogner, 2006.

Definition at line 46 of file discontinuity_measure.h.

---

Member Typedef Documentation

typedef std::map<std::pair<const System*, unsigned int>, ErrorVector*> ErrorEstimator::ErrorMap [inherited]

When calculating many error vectors at once, we need a data structure to hold them all

Definition at line 105 of file error_estimator.h.

---

Constructor & Destructor Documentation

DiscontinuityMeasure::DiscontinuityMeasure

(

)

[inline]

Constructor. Responsible for initializing the _bc_function function pointer to NULL. Defaults to L2 norm; changes to system norm are ignored.

Definition at line 55 of file discontinuity_measure.h.

References ErrorEstimator::error_norm, and libMeshEnums::L2.

: _bc_function(NULL) { error_norm = L2; }

DiscontinuityMeasure::~DiscontinuityMeasure

(

)

[inline]

Destructor.

Definition at line 60 of file discontinuity_measure.h.

{}

---

Member Function Documentation

void DiscontinuityMeasure::attach_essential_bc_function

(

std::pair< bool, Real >

fptrconst System &system,const Point &p,const std::string &var_name

)

Register a user function to use in computing the essential BCs. The return value is std::pair<bool, Real>

Definition at line 163 of file discontinuity_measure.C.

{
  _bc_function = fptr;

// We may be turning boundary side integration on or off
  if (fptr)
    integrate_boundary_sides = true;
  else
    integrate_boundary_sides = false;
}

bool DiscontinuityMeasure::boundary_side_integration

(

)

[protected, virtual]

The function which calculates a normal derivative jump based error term on a boundary side. Returns true if the flux bc function is in fact defined on the current side.

Reimplemented from JumpErrorEstimator.

Definition at line 95 of file discontinuity_measure.C.

References _bc_function, ErrorEstimator::error_norm, JumpErrorEstimator::fe_fine, JumpErrorEstimator::fine_elem, JumpErrorEstimator::fine_error, Elem::hmax(), libmesh_norm(), my_system, JumpErrorEstimator::Ufine, JumpErrorEstimator::var, System::variable_name(), and SystemNorm::weight().

{
  const std::string &var_name = my_system->variable_name(var);

  std::vector<std::vector<Real> > phi_fine = fe_fine->get_phi();
  std::vector<Real> JxW_face = fe_fine->get_JxW();
  std::vector<Point> qface_point = fe_fine->get_xyz();

  // The reinitialization also recomputes the locations of
  // the quadrature points on the side.  By checking if the
  // first quadrature point on the side is on an essential boundary
  // for a particular variable, we will determine if the whole
  // element is on an essential boundary (assuming quadrature points
  // are strictly contained in the side).
  if (this->_bc_function(*my_system, qface_point[0], var_name).first)
    {
      const Real h = fine_elem->hmax();
                    
      // The number of quadrature points
      const unsigned int n_qp = fe_fine->n_quadrature_points();

      // The error contribution from this face
      Real error = 1.e-30;
                    
      // loop over the integration points on the face.
      for (unsigned int qp=0; qp<n_qp; qp++)
        {
          // Value of the imposed essential BC at this quadrature point.
          const std::pair<bool,Real> essential_bc =
            this->_bc_function(*my_system, qface_point[qp], var_name);

          // Be sure the BC function still thinks we're on the 
          // essential boundary.
          libmesh_assert (essential_bc.first == true);

          // The solution gradient from each element
          Number u_fine = 0.;

          // Compute the solution gradient on element e
          for (unsigned int i=0; i != Ufine.size(); i++)
            u_fine += phi_fine[i][qp] * Ufine(i);

          // The difference between the desired BC and the approximate solution. 
          const Number jump = essential_bc.second - u_fine;

          // The flux jump squared.  If using complex numbers,
          // libmesh_norm(z) returns |z|^2, where |z| is the modulus of z.
          const Real jump2 = libmesh_norm(jump);

          // Integrate the error on the face.  The error is
          // scaled by an additional power of h, where h is
          // the maximum side length for the element.  This
          // arises in the definition of the indicator.
          error += JxW_face[qp]*jump2;                  
                        
        } // End quadrature point loop

      fine_error = error*h*error_norm.weight(var);

      return true;
    } // end if side on flux boundary
  return false;
}

float JumpErrorEstimator::coarse_n_flux_faces_increment

(

)

[protected, inherited]

A utility function to correctly increase n_flux_faces for the coarse element

Definition at line 433 of file jump_error_estimator.C.

References JumpErrorEstimator::coarse_elem, Elem::dim(), dim, JumpErrorEstimator::fine_elem, and Elem::level().

Referenced by JumpErrorEstimator::estimate_error().

{
  // Keep track of the number of internal flux sides found on each
  // element
  unsigned int dim = coarse_elem->dim();

  const unsigned int divisor =
    1 << (dim-1)*(fine_elem->level() - coarse_elem->level());

  // With a difference of n levels between fine and coarse elements,
  // we compute a fractional flux face for the coarse element by adding:
  // 1/2^n in 2D
  // 1/4^n in 3D
  // each time.  This code will get hit 2^n times in 2D and 4^n
  // times in 3D so that the final flux face count for the coarse
  // element will be an integer value.

  return 1.0 / static_cast<Real>(divisor);
}

void JumpErrorEstimator::estimate_error	(	const System &	system,
		ErrorVector &	error_per_cell,
		const NumericVector< Number > *	solution_vector = NULL,
		bool	estimate_parent_error = false
	)			[virtual, inherited]

This function uses the derived class's jump error estimate formula to estimate the error on each cell. The estimated error is output in the vector error_per_cell

Implements ErrorEstimator.

Definition at line 52 of file jump_error_estimator.C.

References Elem::active(), MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), JumpErrorEstimator::boundary_side_integration(), FEBase::build(), Elem::child(), JumpErrorEstimator::coarse_elem, JumpErrorEstimator::coarse_error, JumpErrorEstimator::coarse_n_flux_faces_increment(), FEBase::coarsened_dof_values(), System::current_solution(), FEType::default_quadrature_order(), dim, DofMap::dof_indices(), ErrorEstimator::error_norm, JumpErrorEstimator::fe_coarse, JumpErrorEstimator::fe_fine, FEBase::fe_type, JumpErrorEstimator::fine_elem, JumpErrorEstimator::fine_error, JumpErrorEstimator::fine_side, System::get_dof_map(), System::get_mesh(), DofObject::id(), JumpErrorEstimator::initialize(), JumpErrorEstimator::integrate_boundary_sides, JumpErrorEstimator::internal_side_integration(), Elem::level(), MeshBase::max_elem_id(), mesh, MeshBase::mesh_dimension(), Elem::n_children(), Elem::n_neighbors(), System::n_vars(), n_vars, Elem::neighbor(), Elem::parent(), FEBase::qrule, ErrorEstimator::reduce_error(), JumpErrorEstimator::reinit_sides(), JumpErrorEstimator::scale_by_n_flux_faces, System::solution, NumericVector< T >::swap(), JumpErrorEstimator::Ucoarse, JumpErrorEstimator::Ufine, JumpErrorEstimator::var, DofMap::variable_type(), and SystemNorm::weight().

Referenced by Adaptive< T >::solve().

{
  START_LOG("estimate_error()", "JumpErrorEstimator");
  /*

  Conventions for assigning the direction of the normal:
  
  - e & f are global element ids
  
  Case (1.) Elements are at the same level, e<f
            Compute the flux jump on the face and
            add it as a contribution to error_per_cell[e]
            and error_per_cell[f]
  
                   ----------------------
                  |           |          |
                  |           |    f     |
                  |           |          |
                  |    e      |---> n    | 
                  |           |          |
                  |           |          |
                   ----------------------


   Case (2.) The neighbor is at a higher level.
             Compute the flux jump on e's face and
             add it as a contribution to error_per_cell[e]
             and error_per_cell[f]

                   ----------------------
                  |     |     |          |
                  |     |  e  |---> n    |
                  |     |     |          |
                  |-----------|    f     | 
                  |     |     |          |
                  |     |     |          |
                  |     |     |          |
                   ----------------------
  */
   
  // The current mesh
  const MeshBase& mesh = system.get_mesh();

  // The dimensionality of the mesh
  const unsigned int dim = mesh.mesh_dimension();
  
  // The number of variables in the system
  const unsigned int n_vars = system.n_vars();
  
  // The DofMap for this system
  const DofMap& dof_map = system.get_dof_map();

  // Resize the error_per_cell vector to be
  // the number of elements, initialize it to 0.
  error_per_cell.resize (mesh.max_elem_id());
  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);

  // Declare a vector of floats which is as long as
  // error_per_cell above, and fill with zeros.  This vector will be
  // used to keep track of the number of edges (faces) on each active
  // element which are either:
  // 1) an internal edge
  // 2) an edge on a Neumann boundary for which a boundary condition
  //    function has been specified.
  // The error estimator can be scaled by the number of flux edges (faces)
  // which the element actually has to obtain a more uniform measure
  // of the error.  Use floats instead of ints since in case 2 (above)
  // f gets 1/2 of a flux face contribution from each of his
  // neighbors
  std::vector<float> n_flux_faces (error_per_cell.size());
  
  // Prepare current_local_solution to localize a non-standard
  // solution vector if necessary
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number>* newsol =
        const_cast<NumericVector<Number>*>(solution_vector);
      System &sys = const_cast<System&>(system);
      newsol->swap(*sys.solution);
      sys.update();
    }
  
  // Loop over all the variables in the system
  for (var=0; var<n_vars; var++)
    {
      // Possibly skip this variable
      if (error_norm.weight(var) == 0.0) continue;

      // The type of finite element to use for this variable
      const FEType& fe_type = dof_map.variable_type (var);
      
      // Finite element objects for the same face from
      // different sides
      fe_fine = FEBase::build (dim, fe_type);
      fe_coarse = FEBase::build (dim, fe_type);

      // Build an appropriate Gaussian quadrature rule
      QGauss qrule (dim-1, fe_type.default_quadrature_order());

      // Tell the finite element for the fine element about the quadrature
      // rule.  The finite element for the coarse element need not know about it
      fe_fine->attach_quadrature_rule (&qrule);
      
      // By convention we will always do the integration
      // on the face of element e.  We'll need its Jacobian values and
      // physical point locations, at least
      fe_fine->get_JxW();
      fe_fine->get_xyz();

      // Our derived classes may want to do some initialization here
      this->initialize(system, error_per_cell, estimate_parent_error);
      
      // The global DOF indices for elements e & f
      std::vector<unsigned int> dof_indices_fine;
      std::vector<unsigned int> dof_indices_coarse;


      
      // Iterate over all the active elements in the mesh
      // that live on this processor.
      MeshBase::const_element_iterator       elem_it  = mesh.active_local_elements_begin();
      const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end(); 

      for (; elem_it != elem_end; ++elem_it)
        {
          // e is necessarily an active element on the local processor
          const Elem* e = *elem_it;
          const unsigned int e_id = e->id();

#ifdef LIBMESH_ENABLE_AMR
          // See if the parent of element e has been examined yet;
          // if not, we may want to compute the estimator on it
          const Elem* parent = e->parent();

          // We only can compute and only need to compute on
          // parents with all active children
          bool compute_on_parent = true;
          if (!parent || !estimate_parent_error)
            compute_on_parent = false;
          else
            for (unsigned int c=0; c != parent->n_children(); ++c)
              if (!parent->child(c)->active())
                compute_on_parent = false;
             
          if (compute_on_parent &&
              !error_per_cell[parent->id()])
            {
              // Compute a projection onto the parent
              DenseVector<Number> Uparent;
              FEBase::coarsened_dof_values(*(system.solution),
                                           dof_map, parent, Uparent,
                                           var, false);

              // Loop over the neighbors of the parent
              for (unsigned int n_p=0; n_p<parent->n_neighbors(); n_p++)
                {
                  if (parent->neighbor(n_p) != NULL) // parent has a neighbor here
                    {
                      // Find the active neighbors in this direction
                      std::vector<const Elem*> active_neighbors;
                      parent->neighbor(n_p)->
                        active_family_tree_by_neighbor(active_neighbors,
                                                       parent);
                      // Compute the flux to each active neighbor
                      for (unsigned int a=0; 
                           a != active_neighbors.size(); ++a)
                        {
                          const Elem *f = active_neighbors[a];
                      // FIXME - what about when f->level <
                      // parent->level()??
                          if (f->level() >= parent->level())
                            {
                              fine_elem = f;
                              coarse_elem = parent;
                              Ucoarse = Uparent;

                              dof_map.dof_indices (fine_elem, dof_indices_fine, var);
                              const unsigned int n_dofs_fine = dof_indices_fine.size();
                              Ufine.resize(n_dofs_fine);
                        
                              for (unsigned int i=0; i<n_dofs_fine; i++)
                                Ufine(i) = system.current_solution(dof_indices_fine[i]);
                              this->reinit_sides();
                              this->internal_side_integration();

                              error_per_cell[fine_elem->id()] += fine_error;
                              error_per_cell[coarse_elem->id()] += coarse_error;

                              // Keep track of the number of internal flux
                              // sides found on each element
                              n_flux_faces[fine_elem->id()]++;
                              n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
                            }
                        }
                    }
                  else if (integrate_boundary_sides)
                    {
                      fine_elem = parent;
                      Ufine = Uparent;

                      // Reinitialize shape functions on the fine element side
                      fe_fine->reinit (fine_elem, fine_side);

                      if (this->boundary_side_integration())
                        {
                          error_per_cell[fine_elem->id()] += fine_error;
                          n_flux_faces[fine_elem->id()]++;
                        }
                    } 
                }
            }
#endif // #ifdef LIBMESH_ENABLE_AMR

          // If we do any more flux integration, e will be the fine element
          fine_elem = e;

          // Loop over the neighbors of element e
          for (unsigned int n_e=0; n_e<e->n_neighbors(); n_e++)
            {
              fine_side = n_e;

              if (e->neighbor(n_e) != NULL) // e is not on the boundary
                {
                  const Elem* f           = e->neighbor(n_e);
                  const unsigned int f_id = f->id();

                  // Compute flux jumps if we are in case 1 or case 2.
                  if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
                      || (f->level() < e->level()))
                    {               
                      // f is now the coarse element
                      coarse_elem = f;

                      // Get the DOF indices for the two elements
                      dof_map.dof_indices (fine_elem, dof_indices_fine, var);
                      dof_map.dof_indices (coarse_elem, dof_indices_coarse, var);

                      // The number of DOFS on each element
                      const unsigned int n_dofs_fine = dof_indices_fine.size();
                      const unsigned int n_dofs_coarse = dof_indices_coarse.size();
                      Ufine.resize(n_dofs_fine);
                      Ucoarse.resize(n_dofs_coarse);

                      // The local solutions on each element
                      for (unsigned int i=0; i<n_dofs_fine; i++)
                        Ufine(i) = system.current_solution(dof_indices_fine[i]);
                      for (unsigned int i=0; i<n_dofs_coarse; i++)
                        Ucoarse(i) = system.current_solution(dof_indices_coarse[i]);
                        
                      this->reinit_sides();
                      this->internal_side_integration();

                      error_per_cell[fine_elem->id()] += fine_error;
                      error_per_cell[coarse_elem->id()] += coarse_error;

                      // Keep track of the number of internal flux
                      // sides found on each element
                      n_flux_faces[fine_elem->id()]++;
                      n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
                    } // end if (case1 || case2)
                } // if (e->neigbor(n_e) != NULL)

              // Otherwise, e is on the boundary.  If it happens to
              // be on a Dirichlet boundary, we need not do anything.
              // On the other hand, if e is on a Neumann (flux) boundary
              // with grad(u).n = g, we need to compute the additional residual
              // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
              // We can only do this with some knowledge of the boundary
              // conditions, i.e. the user must have attached an appropriate
              // BC function.
              else
                {
                  if (integrate_boundary_sides)
                    {
                      // Reinitialize shape functions on the fine element side
                      fe_fine->reinit (fine_elem, fine_side);

                      // Get the DOF indices 
                      dof_map.dof_indices (fine_elem, dof_indices_fine, var);

                      // The number of DOFS on each element
                      const unsigned int n_dofs_fine = dof_indices_fine.size();
                      Ufine.resize(n_dofs_fine);

                      for (unsigned int i=0; i<n_dofs_fine; i++)
                        Ufine(i) = system.current_solution(dof_indices_fine[i]);

                      if (this->boundary_side_integration())
                        {
                          error_per_cell[fine_elem->id()] += fine_error;
                          n_flux_faces[fine_elem->id()]++;
                        }
                    } // end if _bc_function != NULL
                } // end if (e->neighbor(n_e) == NULL)
            } // end loop over neighbors
        } // End loop over active local elements
    } // End loop over variables


  
  // Each processor has now computed the error contribuions
  // for its local elements.  We need to sum the vector
  // and then take the square-root of each component.  Note
  // that we only need to sum if we are running on multiple
  // processors, and we only need to take the square-root
  // if the value is nonzero.  There will in general be many
  // zeros for the inactive elements.

  // First sum the vector of estimated error values
  this->reduce_error(error_per_cell);

  // Compute the square-root of each component.
  for (unsigned int i=0; i<error_per_cell.size(); i++)
    if (error_per_cell[i] != 0.)
      error_per_cell[i] = std::sqrt(error_per_cell[i]);


  if (this->scale_by_n_flux_faces)
    {
      // Sum the vector of flux face counts
      this->reduce_error(n_flux_faces);

      // Sanity check: Make sure the number of flux faces is
      // always an integer value
#ifdef DEBUG
      for (unsigned int i=0; i<n_flux_faces.size(); ++i)
        libmesh_assert (n_flux_faces[i] == static_cast<float>(static_cast<unsigned int>(n_flux_faces[i])) );
#endif
  
      // Scale the error by the number of flux faces for each element
      for (unsigned int i=0; i<n_flux_faces.size(); ++i)
        {
          if (n_flux_faces[i] == 0.0) // inactive or non-local element
            continue;
      
          //std::cout << "Element " << i << " has " << n_flux_faces[i] << " flux faces." << std::endl;
          error_per_cell[i] /= static_cast<Real>(n_flux_faces[i]); 
        }
    }
  
  // If we used a non-standard solution before, now is the time to fix
  // the current_local_solution
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number>* newsol =
        const_cast<NumericVector<Number>*>(solution_vector);
      System &sys = const_cast<System&>(system);
      newsol->swap(*sys.solution);
      sys.update();
    }

  STOP_LOG("estimate_error()", "JumpErrorEstimator");
}

void ErrorEstimator::estimate_errors	(	const EquationSystems &	equation_systems,
		ErrorMap &	errors_per_cell,
		const std::map< const System *, const NumericVector< Number > * > *	solution_vectors = NULL,
		bool	estimate_parent_error = false
	)			[virtual, inherited]

This virtual function can be redefined in derived classes, but by default it calls estimate_error repeatedly to calculate the requested error vectors.

Currently this function ignores the error_norm.weight() values because it calculates each variable's error individually, unscaled.

The user selects which errors get computed by filling a map with error vectors: If errors_per_cell[&system][v] exists, it will be filled with the error values in variable v of system

FIXME: This is a default implementation - derived classes should reimplement it for efficiency.

Definition at line 92 of file error_estimator.C.

References ErrorEstimator::error_norm, ErrorEstimator::estimate_error(), EquationSystems::get_system(), EquationSystems::n_systems(), System::n_vars(), n_vars, and SystemNorm::type().

{
  SystemNorm old_error_norm = this->error_norm;

  // Find the requested error values from each system
  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
    {
      const System &sys = equation_systems.get_system(s);

      unsigned int n_vars = sys.n_vars();

      for (unsigned int v = 0; v != n_vars; ++v)
        {
          // Only fill in ErrorVectors the user asks for
          if (errors_per_cell.find(std::make_pair(&sys, v)) ==
              errors_per_cell.end())
            continue;

          // Calculate error in only one variable
          std::vector<Real> weights(n_vars, 0.0);
          weights[v] = 1.0;
          this->error_norm =
            SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(0)),
                       weights);

          const NumericVector<Number>* solution_vector = NULL;
          if (solution_vectors &&
              solution_vectors->find(&sys) != solution_vectors->end())
            solution_vector = solution_vectors->find(&sys)->second;

          this->estimate_error
            (sys, *errors_per_cell[std::make_pair(&sys, v)],
             solution_vector, estimate_parent_error);
        }
    }

  // Restore our old state before returning
  this->error_norm = old_error_norm;
}

void ErrorEstimator::estimate_errors	(	const EquationSystems &	equation_systems,
		ErrorVector &	error_per_cell,
		const std::map< const System *, SystemNorm > &	error_norms,
		const std::map< const System *, const NumericVector< Number > * > *	solution_vectors = NULL,
		bool	estimate_parent_error = false
	)			[virtual, inherited]

This virtual function can be redefined in derived classes, but by default computes the sum of the error_per_cell for each system in the equation_systems.

Currently this function ignores the error_norm member variable, and uses the function argument error_norms instead.

This function is named estimate_errors instead of estimate_error because otherwise C++ can get confused.

Reimplemented in UniformRefinementEstimator.

Definition at line 46 of file error_estimator.C.

References ErrorEstimator::error_norm, ErrorEstimator::estimate_error(), EquationSystems::get_system(), and EquationSystems::n_systems().

{
  SystemNorm old_error_norm = this->error_norm;

  // Sum the error values from each system
  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
    {
      ErrorVector system_error_per_cell;
      const System &sys = equation_systems.get_system(s);
      if (error_norms.find(&sys) == error_norms.end())
        this->error_norm = old_error_norm;
      else
        this->error_norm = error_norms.find(&sys)->second;

      const NumericVector<Number>* solution_vector = NULL;
      if (solution_vectors &&
          solution_vectors->find(&sys) != solution_vectors->end())
        solution_vector = solution_vectors->find(&sys)->second;

      this->estimate_error(sys, system_error_per_cell,
                           solution_vector, estimate_parent_error);

      if (s)
        {
          libmesh_assert(error_per_cell.size() == system_error_per_cell.size());
          for (unsigned int i=0; i != error_per_cell.size(); ++i)
            error_per_cell[i] += system_error_per_cell[i];
        }
      else
        error_per_cell = system_error_per_cell;
    }

  // Restore our old state before returning
  this->error_norm = old_error_norm;
}

void DiscontinuityMeasure::initialize	(	const System &	system,
		ErrorVector &	error_per_cell,
		bool	estimate_parent_error
	)			[protected, virtual]

An initialization function, for requesting specific data from the FE objects

Reimplemented from JumpErrorEstimator.

Definition at line 41 of file discontinuity_measure.C.

References JumpErrorEstimator::fe_coarse, JumpErrorEstimator::fe_fine, and my_system.

{
  // Hang onto the system - we may need it for variable names later.
  my_system = &system;

  // We'll need values for jump computation
  fe_fine->get_phi();
  fe_coarse->get_phi();
}

void DiscontinuityMeasure::internal_side_integration

(

)

[protected, virtual]

The function which calculates a normal derivative jump based error term on an internal side

Implements JumpErrorEstimator.

Definition at line 56 of file discontinuity_measure.C.

References JumpErrorEstimator::coarse_elem, JumpErrorEstimator::coarse_error, ErrorEstimator::error_norm, JumpErrorEstimator::fe_coarse, JumpErrorEstimator::fe_fine, JumpErrorEstimator::fine_elem, JumpErrorEstimator::fine_error, Elem::hmax(), libmesh_norm(), JumpErrorEstimator::Ucoarse, JumpErrorEstimator::Ufine, JumpErrorEstimator::var, and SystemNorm::weight().

{
  Real error = 1.e-30;
  unsigned int n_qp = fe_fine->n_quadrature_points();
  unsigned int n_fine_dofs = Ufine.size();
  unsigned int n_coarse_dofs = Ucoarse.size();

  std::vector<std::vector<Real> > phi_coarse = fe_coarse->get_phi();
  std::vector<std::vector<Real> > phi_fine = fe_fine->get_phi();
  std::vector<Real> JxW_face = fe_fine->get_JxW();

  for (unsigned int qp=0; qp != n_qp; ++qp)
    {
      // Calculate solution values on fine and coarse elements
      // at this quadrature point
      Number u_fine = 0., u_coarse = 0.;
      for (unsigned int i=0; i != n_coarse_dofs; ++i)
        u_coarse += phi_coarse[i][qp] * Ucoarse(i);

      for (unsigned int i=0; i != n_fine_dofs; ++i)
        u_fine += phi_fine[i][qp] * Ufine(i);
                                
      // Find the jump in the value
      // at this quadrature point
      const Number jump = u_fine - u_coarse;
      const Real jump2 = libmesh_norm(jump);
      // Accumulate the jump integral
      error += JxW_face[qp] * jump2;
    }

  // Add the h-weighted jump integral to each error term
  fine_error =
    error * fine_elem->hmax() * error_norm.weight(var);
  coarse_error =
    error * coarse_elem->hmax() * error_norm.weight(var);
}

void ErrorEstimator::reduce_error

(

std::vector< float > &

error_per_cell

)

const [protected, inherited]

This method takes the local error contributions in error_per_cell from each processor and combines them to get the global error vector.

Referenced by PatchRecoveryErrorEstimator::estimate_error(), and JumpErrorEstimator::estimate_error().

void JumpErrorEstimator::reinit_sides

(

)

[protected, inherited]

A utility function to reinit the finite element data on elements sharing a side

Definition at line 412 of file jump_error_estimator.C.

References JumpErrorEstimator::coarse_elem, Elem::dim(), JumpErrorEstimator::fe_coarse, JumpErrorEstimator::fe_fine, JumpErrorEstimator::fine_elem, JumpErrorEstimator::fine_side, and InfFE< Dim, T_radial, T_map >::inverse_map().

Referenced by JumpErrorEstimator::estimate_error().

{
  // The master quadrature point locations on the coarse element
  std::vector<Point> qp_coarse;

  // Reinitialize shape functions on the fine element side
  fe_fine->reinit (fine_elem, fine_side);

  // Get the physical locations of the fine element quadrature points
  std::vector<Point> qface_point = fe_fine->get_xyz();

  // Find their locations on the coarse element
  FEInterface::inverse_map (coarse_elem->dim(), fe_coarse->get_fe_type(),
                            coarse_elem, qface_point, qp_coarse);

  // Calculate the coarse element shape functions at those locations
  fe_coarse->reinit (coarse_elem, &qp_coarse);
}

---

Member Data Documentation

std::pair<bool,Real>(* DiscontinuityMeasure::_bc_function)(const System &system, const Point &p, const std::string &var_name) [protected]

Pointer to function that returns BC information.

Referenced by boundary_side_integration().

const Elem * JumpErrorEstimator::coarse_elem [protected, inherited]

Definition at line 130 of file jump_error_estimator.h.

Referenced by JumpErrorEstimator::coarse_n_flux_faces_increment(), JumpErrorEstimator::estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), internal_side_integration(), and JumpErrorEstimator::reinit_sides().

Real JumpErrorEstimator::coarse_error [protected, inherited]

Definition at line 135 of file jump_error_estimator.h.

Referenced by JumpErrorEstimator::estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and internal_side_integration().

SystemNorm ErrorEstimator::error_norm [inherited]

When estimating the error in a single system, the error_norm is used to control the scaling and norm choice for each variable. Not all estimators will support all norm choices. The default scaling is for all variables to be weighted equally. The default norm choice depends on the error estimator.

Part of this functionality was supported via component_scale and sobolev_order in older libMesh versions, and a small part was supported via component_mask in even older versions. Hopefully the encapsulation here will allow us to avoid changing this API again.

Definition at line 137 of file error_estimator.h.

Referenced by KellyErrorEstimator::boundary_side_integration(), boundary_side_integration(), DiscontinuityMeasure(), JumpErrorEstimator::estimate_error(), ErrorEstimator::estimate_errors(), ExactErrorEstimator::ExactErrorEstimator(), ExactErrorEstimator::find_squared_element_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), internal_side_integration(), KellyErrorEstimator::KellyErrorEstimator(), LaplacianErrorEstimator::LaplacianErrorEstimator(), PatchRecoveryErrorEstimator::EstimateError::operator()(), PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), and UniformRefinementEstimator::UniformRefinementEstimator().

AutoPtr<FEBase> JumpErrorEstimator::fe_coarse [protected, inherited]

Definition at line 155 of file jump_error_estimator.h.

Referenced by JumpErrorEstimator::estimate_error(), KellyErrorEstimator::initialize(), LaplacianErrorEstimator::initialize(), initialize(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), internal_side_integration(), and JumpErrorEstimator::reinit_sides().

AutoPtr<FEBase> JumpErrorEstimator::fe_fine [protected, inherited]

The finite element objects for fine and coarse elements

Definition at line 155 of file jump_error_estimator.h.

Referenced by KellyErrorEstimator::boundary_side_integration(), boundary_side_integration(), JumpErrorEstimator::estimate_error(), KellyErrorEstimator::initialize(), LaplacianErrorEstimator::initialize(), initialize(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), internal_side_integration(), and JumpErrorEstimator::reinit_sides().

const Elem* JumpErrorEstimator::fine_elem [protected, inherited]

The fine and coarse elements sharing a face

Definition at line 130 of file jump_error_estimator.h.

Referenced by KellyErrorEstimator::boundary_side_integration(), boundary_side_integration(), JumpErrorEstimator::coarse_n_flux_faces_increment(), JumpErrorEstimator::estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), internal_side_integration(), and JumpErrorEstimator::reinit_sides().

Real JumpErrorEstimator::fine_error [protected, inherited]

The fine and coarse error values to be set by each side_integration();

Definition at line 135 of file jump_error_estimator.h.

Referenced by KellyErrorEstimator::boundary_side_integration(), boundary_side_integration(), JumpErrorEstimator::estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and internal_side_integration().

unsigned int JumpErrorEstimator::fine_side [protected, inherited]

Which side of the fine element is this?

Definition at line 140 of file jump_error_estimator.h.

Referenced by JumpErrorEstimator::estimate_error(), and JumpErrorEstimator::reinit_sides().

bool JumpErrorEstimator::integrate_boundary_sides [protected, inherited]

A boolean flag, by default false, to be set to true if integrations with boundary_side_integration() should be performed

Definition at line 125 of file jump_error_estimator.h.

Referenced by KellyErrorEstimator::attach_flux_bc_function(), and JumpErrorEstimator::estimate_error().

const System* DiscontinuityMeasure::my_system [protected]

A pointer to the current System

Definition at line 96 of file discontinuity_measure.h.

Referenced by boundary_side_integration(), and initialize().

bool JumpErrorEstimator::scale_by_n_flux_faces [inherited]

This boolean flag allows you to scale the error indicator result for each element by the number of "flux faces" the element actually has. This tends to weight more evenly cells which are on the boundaries and thus have fewer contributions to their flux. The value is initialized to false, simply set it to true if you want to use the feature.

Definition at line 86 of file jump_error_estimator.h.

Referenced by JumpErrorEstimator::estimate_error().

DenseVector<Number> JumpErrorEstimator::Ucoarse [protected, inherited]

Definition at line 150 of file jump_error_estimator.h.

Referenced by JumpErrorEstimator::estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and internal_side_integration().

DenseVector<Number> JumpErrorEstimator::Ufine [protected, inherited]

The local degree of freedom values on fine and coarse elements

Definition at line 150 of file jump_error_estimator.h.

Referenced by KellyErrorEstimator::boundary_side_integration(), boundary_side_integration(), JumpErrorEstimator::estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and internal_side_integration().

unsigned int JumpErrorEstimator::var [protected, inherited]

The variable number currently being evaluated

Definition at line 145 of file jump_error_estimator.h.

Referenced by KellyErrorEstimator::boundary_side_integration(), boundary_side_integration(), JumpErrorEstimator::estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and internal_side_integration().

---

The documentation for this class was generated from the following files:

• discontinuity_measure.h
• discontinuity_measure.C