libMesh::ErrorVector Class Reference

#include <error_vector.h>

Inheritance diagram for libMesh::ErrorVector:

Public Member Functions
	ErrorVector (dof_id_type i=0, MeshBase *mesh=nullptr)
	ErrorVector (dof_id_type i, ErrorVectorReal val)
virtual ErrorVectorReal	minimum () const override
virtual Real	mean () const override
virtual Real	median () override
virtual Real	median () const override
virtual Real	variance () const override
virtual Real	variance (const Real mean) const override
virtual std::vector< dof_id_type >	cut_below (Real cut) const override
virtual std::vector< dof_id_type >	cut_above (Real cut) const override
void	plot_error (const std::string &filename, const MeshBase &mesh) const
virtual Real	l2_norm () const
virtual ErrorVectorReal	maximum () const
virtual Real	stddev () const
virtual Real	stddev (const Real known_mean) const
void	normalize ()
virtual void	histogram (std::vector< dof_id_type > &bin_members, unsigned int n_bins=10)
virtual void	histogram (std::vector< dof_id_type > &bin_members, unsigned int n_bins=10) const
void	plot_histogram (const processor_id_type my_procid, const std::string &filename, unsigned int n_bins)

Protected Member Functions
bool	is_active_elem (dof_id_type i) const

Protected Attributes
MeshBase *	_mesh

Detailed Description

The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite element mesh. In general, when computing the error on a mesh only the active elements are considered, but the ErrorVector is sized according to the total number of elements in the mesh. The ErrorVector is thus padded with zeros for all the inactive elements, and this must be taken into account when calculating the statistics. Since the error is a positive quantity this class assumes it contains positive data (i.e. min_val >= 0.).

Author

Benjamin S. Kirk

Date

2003

Definition at line 50 of file error_vector.h.

Constructor & Destructor Documentation

ErrorVector() [1/2]

libMesh::ErrorVector::ErrorVector ( dof_id_type  i = 0, MeshBase *  mesh = nullptr  )

inline

ErrorVector constructor; sets initial length to i.

If mesh is not null, MeshBase::elem() and Elem::is_active() will be used to distinguish active and inactive elements. If mesh is null, ErrorVector will assume that all 0.0 error values correspond to inactive elements and all non-zero error values correspond to active elements.

Definition at line 62 of file error_vector.h.

                                                          :
    StatisticsVector<ErrorVectorReal> (i),
    _mesh(mesh)
  {}

ErrorVector() [2/2]

libMesh::ErrorVector::ErrorVector ( dof_id_type  i, ErrorVectorReal  val  )

inline

ErrorVector constructor; sets initial length to i and initial values to val.

If mesh is not null, MeshBase::elem() and Elem::is_active() will be used to distinguish active and inactive elements. If mesh is null, ErrorVector will assume that all 0.0 error values correspond to inactive elements and all non-zero error values correspond to active elements.

Definition at line 75 of file error_vector.h.

                                                  :
    StatisticsVector<ErrorVectorReal> (i,val) {}

Member Function Documentation

cut_above()

std::vector< dof_id_type > libMesh::ErrorVector::cut_above ( Real  cut ) const

overridevirtual

Returns

A vector of dof_id_types which correspond to the indices of every member of the data set above the cutoff value cut ignoring inactive elements.

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 174 of file error_vector.C.

References is_active_elem().

{
  LOG_SCOPE ("cut_above()", "ErrorVector");

  const dof_id_type n = cast_int<dof_id_type>(this->size());

  std::vector<dof_id_type> cut_indices;
  cut_indices.reserve(n/2);  // Arbitrary

  for (dof_id_type i=0; i<n; i++)
    if (this->is_active_elem(i))
      {
        if ((*this)[i] > cut)
          {
            cut_indices.push_back(i);
          }
      }

  return cut_indices;
}

cut_below()

std::vector< dof_id_type > libMesh::ErrorVector::cut_below ( Real  cut ) const

overridevirtual

Returns

A vector of dof_id_types which correspond to the indices of every member of the data set below the cutoff value cut ignoring inactive elements.

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 150 of file error_vector.C.

References is_active_elem().

{
  LOG_SCOPE ("cut_below()", "ErrorVector");

  const dof_id_type n = cast_int<dof_id_type>(this->size());

  std::vector<dof_id_type> cut_indices;
  cut_indices.reserve(n/2);  // Arbitrary

  for (dof_id_type i=0; i<n; i++)
    if (this->is_active_elem(i))
      {
        if ((*this)[i] < cut)
          {
            cut_indices.push_back(i);
          }
      }

  return cut_indices;
}

histogram() [1/2]

void libMesh::StatisticsVector< ErrorVectorReal >::histogram ( std::vector< dof_id_type > &  bin_members, unsigned int  n_bins = 10  )

virtualinherited

Returns

A histogram with n_bins bins for the data set.

For simplicity, the bins are assumed to be of uniform size. Upon return, the bin_members vector will contain unsigned integers which give the number of members in each bin. WARNING: This non-const function sorts the vector, changing its order. Source: GNU Scientific Library.

Definition at line 178 of file statistics.C.

References end, std::max(), std::min(), libMesh::out, and libMesh::Real.

{
  // Must have at least 1 bin
  libmesh_assert (n_bins>0);

  const dof_id_type n = cast_int<dof_id_type>(this->size());

  std::sort(this->begin(), this->end());

  // The StatisticsVector can hold both integer and float types.
  // We will define all the bins, etc. using Reals.
  Real min      = static_cast<Real>(this->minimum());
  Real max      = static_cast<Real>(this->maximum());
  Real bin_size = (max - min) / static_cast<Real>(n_bins);

  LOG_SCOPE ("histogram()", "StatisticsVector");

  std::vector<Real> bin_bounds(n_bins+1);
  for (std::size_t i=0; i<bin_bounds.size(); i++)
    bin_bounds[i] = min + Real(i) * bin_size;

  // Give the last bin boundary a little wiggle room: we don't want
  // it to be just barely less than the max, otherwise our bin test below
  // may fail.
  bin_bounds.back() += 1.e-6 * bin_size;

  // This vector will store the number of members each bin has.
  bin_members.resize(n_bins);

  dof_id_type data_index = 0;
  for (std::size_t j=0; j<bin_members.size(); j++) // bin vector indexing
    {
      // libMesh::out << "(debug) Filling bin " << j << std::endl;

      for (dof_id_type i=data_index; i<n; i++) // data vector indexing
        {
          //libMesh::out << "(debug) Processing index=" << i << std::endl;
          Real current_val = static_cast<Real>( (*this)[i] );

          // There may be entries in the vector smaller than the value
          // reported by this->minimum().  (e.g. inactive elements in an
          // ErrorVector.)  We just skip entries like that.
          if (current_val < min)
            {
              //     libMesh::out << "(debug) Skipping entry v[" << i << "]="
              //       << (*this)[i]
              //       << " which is less than the min value: min="
              //       << min << std::endl;
              continue;
            }

          if (current_val > bin_bounds[j+1]) // if outside the current bin (bin[j] is bounded
            // by bin_bounds[j] and bin_bounds[j+1])
            {
              // libMesh::out.precision(16);
              //     libMesh::out.setf(std::ios_base::fixed);
              //     libMesh::out << "(debug) (*this)[i]= " << (*this)[i]
              //       << " is greater than bin_bounds[j+1]="
              //      << bin_bounds[j+1] << std::endl;
              data_index = i; // start searching here for next bin
              break; // go to next bin
            }

          // Otherwise, increment current bin's count
          bin_members[j]++;
          // libMesh::out << "(debug) Binned index=" << i << std::endl;
        }
    }

#ifdef DEBUG
  // Check the number of binned entries
  const dof_id_type n_binned = std::accumulate(bin_members.begin(),
                                               bin_members.end(),
                                               static_cast<dof_id_type>(0),
                                               std::plus<dof_id_type>());

  if (n != n_binned)
    {
      libMesh::out << "Warning: The number of binned entries, n_binned="
                   << n_binned
                   << ", did not match the total number of entries, n="
                   << n << "." << std::endl;
    }
#endif
}

histogram() [2/2]

void libMesh::StatisticsVector< ErrorVectorReal >::histogram ( std::vector< dof_id_type > &  bin_members, unsigned int  n_bins = 10  ) const

virtualinherited

A const version of the histogram function.

Definition at line 313 of file statistics.C.

References libMesh::StatisticsVector< T >::histogram().

{
  StatisticsVector<T> sv = (*this);

  return sv.histogram(bin_members, n_bins);
}

is_active_elem()

bool libMesh::ErrorVector::is_active_elem ( dof_id_type  i ) const

protected

Utility function to decide whether element i is active

Definition at line 197 of file error_vector.C.

References _mesh, libMesh::Elem::active(), and libMesh::MeshBase::elem_ptr().

Referenced by cut_above(), cut_below(), mean(), median(), minimum(), and variance().

{
  libmesh_assert_less (i, this->size());

  if (_mesh)
    {
      return _mesh->elem_ptr(i)->active();
    }
  else
    return ((*this)[i] != 0.);
}

l2_norm()

Real libMesh::StatisticsVector< ErrorVectorReal >::l2_norm ( ) const

virtualinherited

Returns

The l2 norm of the data set.

Definition at line 36 of file statistics.C.

References libMesh::Real.

{
  Real normsq = 0.;
  const dof_id_type n = cast_int<dof_id_type>(this->size());
  for (dof_id_type i = 0; i != n; ++i)
    normsq += ((*this)[i] * (*this)[i]);

  return std::sqrt(normsq);
}

maximum()

ErrorVectorReal libMesh::StatisticsVector< ErrorVectorReal >::maximum ( ) const

virtualinherited

Returns

The maximum value in the data set.

Definition at line 61 of file statistics.C.

References end, and std::max().

{
  LOG_SCOPE ("maximum()", "StatisticsVector");

  const T max = *(std::max_element(this->begin(), this->end()));

  return max;
}

mean()

Real libMesh::ErrorVector::mean ( ) const

overridevirtual

Returns

The mean value of the data set. Ignores zero values.

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 69 of file error_vector.C.

References is_active_elem(), and libMesh::Real.

Referenced by libMesh::MeshRefinement::flag_elements_by_mean_stddev(), and variance().

{
  LOG_SCOPE ("mean()", "ErrorVector");

  const dof_id_type n = cast_int<dof_id_type>(this->size());

  Real the_mean  = 0;
  dof_id_type nnz = 0;

  for (dof_id_type i=0; i<n; i++)
    if (this->is_active_elem(i))
      {
        the_mean += ( static_cast<Real>((*this)[i]) - the_mean ) / (nnz + 1);

        nnz++;
      }

  return the_mean;
}

median() [1/2]

Real libMesh::ErrorVector::median ( )

overridevirtual

Returns

The median (e.g. the middle) value of the data set, ignoring inactive elements.

This function modifies the original data by sorting, so it can't be called on const objects. Source: GNU Scientific Library

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 92 of file error_vector.C.

References is_active_elem(), and libMesh::StatisticsVector< T >::median().

Referenced by median().

{
  const dof_id_type n = cast_int<dof_id_type>(this->size());

  if (n == 0)
    return 0.;


  // Build a StatisticsVector<ErrorVectorReal> containing
  // only our active entries and take its mean
  StatisticsVector<ErrorVectorReal> sv;

  sv.reserve (n);

  for (dof_id_type i=0; i<n; i++)
    if (this->is_active_elem(i))
      sv.push_back((*this)[i]);

  return sv.median();
}

median() [2/2]

Real libMesh::ErrorVector::median ( ) const

overridevirtual

A const version of the median function. Requires twice the memory of original data set but does not change the original.

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 116 of file error_vector.C.

References median().

{
  ErrorVector ev = (*this);

  return ev.median();
}

minimum()

ErrorVectorReal libMesh::ErrorVector::minimum ( ) const

overridevirtual

Returns

The minimum nonzero value in the data set.

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 46 of file error_vector.C.

References libMesh::ErrorVectorReal, is_active_elem(), std::max(), and std::min().

{
  LOG_SCOPE ("minimum()", "ErrorVector");

  const dof_id_type n = cast_int<dof_id_type>(this->size());
  ErrorVectorReal min = std::numeric_limits<ErrorVectorReal>::max();

  for (dof_id_type i=0; i<n; i++)
    {
      // Only positive (or zero) values in the error vector
      libmesh_assert_greater_equal ((*this)[i], 0.);
      if (this->is_active_elem(i))
        min = std::min (min, (*this)[i]);
    }

  // ErrorVectors are for positive values
  libmesh_assert_greater_equal (min, 0.);

  return min;
}

normalize()

void libMesh::StatisticsVector< ErrorVectorReal >::normalize ( )

inherited

Divides all entries by the largest entry and stores the result.

Definition at line 164 of file statistics.C.

References std::max(), and libMesh::Real.

{
  const dof_id_type n = cast_int<dof_id_type>(this->size());
  const Real max = this->maximum();

  for (dof_id_type i=0; i<n; i++)
    (*this)[i] = static_cast<T>((*this)[i] / max);
}

plot_error()

void libMesh::ErrorVector::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 at line 210 of file error_vector.C.

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::MeshBase::clone(), libMesh::CONSTANT, libMesh::Elem::DO_NOTHING, libMesh::DofMap::dof_indices(), libMesh::ENCODE, libMesh::err, libMesh::System::get_dof_map(), libMesh::EquationSystems::init(), mesh, libMesh::MONOMIAL, libMesh::System::solution, libMesh::WRITE, libMesh::Nemesis_IO::write(), libMesh::ExodusII_IO::write(), libMesh::XdrIO::write(), libMesh::EquationSystems::write(), libMesh::EquationSystems::WRITE_ADDITIONAL_DATA, libMesh::EquationSystems::WRITE_DATA, libMesh::GMVIO::write_discontinuous_gmv(), libMesh::Nemesis_IO::write_element_data(), libMesh::ExodusII_IO::write_element_data(), and libMesh::MeshOutput< MT >::write_equation_systems().

Referenced by libMesh::AdjointResidualErrorEstimator::estimate_error().

{
  std::unique_ptr<MeshBase> meshptr = oldmesh.clone();
  MeshBase & mesh = *meshptr;

  // The all_first_order routine will prepare_for_use(), which would
  // break our ordering if elements get changed.
  mesh.allow_renumbering(false);
  mesh.all_first_order();

#ifdef LIBMESH_ENABLE_AMR
  // We don't want p elevation when plotting a single constant value
  // per element
  for (auto & elem : mesh.element_ptr_range())
    {
      elem->set_p_refinement_flag(Elem::DO_NOTHING);
      elem->set_p_level(0);
    }
#endif // LIBMESH_ENABLE_AMR

  EquationSystems temp_es (mesh);
  ExplicitSystem & error_system
    = temp_es.add_system<ExplicitSystem> ("Error");
  error_system.add_variable("error", CONSTANT, MONOMIAL);
  temp_es.init();

  const DofMap & error_dof_map = error_system.get_dof_map();
  std::vector<dof_id_type> dof_indices;

  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      error_dof_map.dof_indices(elem, dof_indices);

      const dof_id_type elem_id = elem->id();

      //0 for the monomial basis
      const dof_id_type solution_index = dof_indices[0];

      // libMesh::out << "elem_number=" << elem_number << std::endl;
      libmesh_assert_less (elem_id, (*this).size());

      // We may have zero error values in special circumstances
      // libmesh_assert_greater ((*this)[elem_id], 0.);
      error_system.solution->set(solution_index, (*this)[elem_id]);
    }

  error_system.solution->close();

  // We may have to renumber if the original numbering was not
  // contiguous.  Since this is just a temporary mesh, that's probably
  // fine.
  if (mesh.max_elem_id() != mesh.n_elem() ||
      mesh.max_node_id() != mesh.n_nodes())
    {
      mesh.allow_renumbering(true);
      mesh.renumber_nodes_and_elements();
    }

  if (filename.rfind(".gmv") < filename.size())
    {
      GMVIO(mesh).write_discontinuous_gmv(filename,
                                          temp_es, false);
    }
  else if (filename.rfind(".plt") < filename.size())
    {
      TecplotIO (mesh).write_equation_systems
        (filename, temp_es);
    }
#if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
  else if ((filename.rfind(".nem") < filename.size()) ||
           (filename.rfind(".n") < filename.size()))
    {
      Nemesis_IO io(mesh);
      io.write(filename);
      io.write_element_data(temp_es);
    }
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
  else if ((filename.rfind(".exo") < filename.size()) ||
           (filename.rfind(".e") < filename.size()))
    {
      ExodusII_IO io(mesh);
      io.write(filename);
      io.write_element_data(temp_es);
    }
#endif
  else if (filename.rfind(".xda") < filename.size())
    {
      XdrIO(mesh).write("mesh-"+filename);
      temp_es.write("soln-"+filename,WRITE,
                    EquationSystems::WRITE_DATA |
                    EquationSystems::WRITE_ADDITIONAL_DATA);
    }
  else if (filename.rfind(".xdr") < filename.size())
    {
      XdrIO(mesh,true).write("mesh-"+filename);
      temp_es.write("soln-"+filename,ENCODE,
                    EquationSystems::WRITE_DATA |
                    EquationSystems::WRITE_ADDITIONAL_DATA);
    }
  else
    {
      libmesh_here();
      libMesh::err << "Warning: ErrorVector::plot_error currently only"
                   << " supports .gmv, .plt, .xdr/.xda, and .exo/.e (if enabled) output;" << std::endl;
      libMesh::err << "Could not recognize filename: " << filename
                   << std::endl;
    }
}

plot_histogram()

void libMesh::StatisticsVector< ErrorVectorReal >::plot_histogram ( const processor_id_type  my_procid, const std::string &  filename, unsigned int  n_bins  )

inherited

Generates a Matlab/Octave style file which can be used to make a plot of the histogram having the desired number of bins. Uses the histogram(...) function in this class WARNING: The histogram(...) function is non-const, and changes the order of the vector.

Definition at line 270 of file statistics.C.

References std::max(), std::min(), and libMesh::Real.

{
  // First generate the histogram with the desired number of bins
  std::vector<dof_id_type> bin_members;
  this->histogram(bin_members, n_bins);

  // The max, min and bin size are used to generate x-axis values.
  T min      = this->minimum();
  T max      = this->maximum();
  T bin_size = (max - min) / static_cast<T>(n_bins);

  // On processor 0: Write histogram to file
  if (my_procid==0)
    {
      std::ofstream out_stream (filename.c_str());

      out_stream << "clear all\n";
      out_stream << "clf\n";
      //out_stream << "x=linspace(" << min << "," << max << "," << n_bins+1 << ");\n";

      // abscissa values are located at the center of each bin.
      out_stream << "x=[";
      for (std::size_t i=0; i<bin_members.size(); ++i)
        {
          out_stream << min + (Real(i)+0.5)*bin_size << " ";
        }
      out_stream << "];\n";

      out_stream << "y=[";
      for (std::size_t i=0; i<bin_members.size(); ++i)
        {
          out_stream << bin_members[i] << " ";
        }
      out_stream << "];\n";
      out_stream << "bar(x,y);\n";
    }
}

stddev() [1/2]

virtual Real libMesh::StatisticsVector< ErrorVectorReal >::stddev ( ) const

inlinevirtualinherited

Returns

The standard deviation of the data set, which is simply the square-root of the variance.

Definition at line 154 of file statistics.h.

References libMesh::StatisticsVector< T >::variance().

  { return std::sqrt(this->variance()); }

stddev() [2/2]

virtual Real libMesh::StatisticsVector< ErrorVectorReal >::stddev ( const Real  known_mean ) const

inlinevirtualinherited

Returns

Computes the standard deviation of the data set, which is simply the square-root of the variance.

This method can be used for efficiency when the mean has already been computed.

Definition at line 164 of file statistics.h.

References libMesh::StatisticsVector< T >::variance().

  { return std::sqrt(this->variance(known_mean)); }

variance() [1/2]

virtual Real libMesh::ErrorVector::variance ( ) const

inlineoverridevirtual

Returns

The variance of the data set ignoring inactive elements.

Uses a recurrence relation to prevent data overflow for large sums.

Note

The variance is equal to the standard deviation squared. The variance is normalized by N in this case. Source: GNU Scientific Library.

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 115 of file error_vector.h.

References mean().

Referenced by libMesh::MeshRefinement::flag_elements_by_mean_stddev().

  { return this->variance(this->mean()); }

variance() [2/2]

Real libMesh::ErrorVector::variance ( const Real  mean ) const

overridevirtual

Returns

The variance of the data set ignoring inactive elements and given the mean.

This is useful for efficiency when you have already calculated the mean. Uses a recurrence relation to prevent data overflow for large sums.

Note

The variance is equal to the standard deviation squared. Source: GNU Scientific Library.

Reimplemented from libMesh::StatisticsVector< ErrorVectorReal >.

Definition at line 126 of file error_vector.C.

References is_active_elem(), and libMesh::Real.

{
  const dof_id_type n = cast_int<dof_id_type>(this->size());

  LOG_SCOPE ("variance()", "ErrorVector");

  Real the_variance = 0;
  dof_id_type nnz = 0;

  for (dof_id_type i=0; i<n; i++)
    if (this->is_active_elem(i))
      {
        const Real delta = ( static_cast<Real>((*this)[i]) - mean_in );
        the_variance += (delta * delta - the_variance) / (nnz + 1);

        nnz++;
      }

  return the_variance;
}

Member Data Documentation

_mesh

MeshBase* libMesh::ErrorVector::_mesh

protected

Pointer to the mesh, which may be used to decide which elements are active

Definition at line 163 of file error_vector.h.

Referenced by is_active_elem().

---

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

• error_vector.h
• error_vector.C