error_estimator.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_ERROR_ESTIMATOR_NAMESPACE_HEADER
#define OOMPH_ERROR_ESTIMATOR_NAMESPACE_HEADER
 
#include "mesh.h"
#include "quadtree.h"
#include "nodes.h"
#include "algebraic_elements.h"
 
namespace oomph
{
  //========================================================================
  //========================================================================
  class ErrorEstimator
  {
  public:
    ErrorEstimator() {}
 
    ErrorEstimator(const ErrorEstimator&) = delete;
 
    void operator=(const ErrorEstimator&) = delete;
 
    virtual ~ErrorEstimator() {}
 
    void get_element_errors(Mesh*& mesh_pt, Vector<double>& elemental_error)
    {
      // Create dummy doc info object and switch off output
      DocInfo doc_info;
      doc_info.disable_doc();
      // Forward call to version with doc.
      get_element_errors(mesh_pt, elemental_error, doc_info);
    }
 
 
    virtual void get_element_errors(Mesh*& mesh_pt,
                                    Vector<double>& elemental_error,
                                    DocInfo& doc_info) = 0;
  };
 
 
  //========================================================================
  //========================================================================
  class ElementWithZ2ErrorEstimator : public virtual FiniteElement
  {
  public:
    ElementWithZ2ErrorEstimator() {}
 
    ElementWithZ2ErrorEstimator(const ElementWithZ2ErrorEstimator&) = delete;
 
    void operator=(const ElementWithZ2ErrorEstimator&) = delete;
 
    virtual unsigned num_Z2_flux_terms() = 0;
 
    virtual unsigned ncompound_fluxes()
    {
      return 1;
    }
 
    virtual void get_Z2_flux(const Vector<double>& s, Vector<double>& flux) = 0;
 
    virtual void compute_exact_Z2_error(
      std::ostream& outfile,
      FiniteElement::SteadyExactSolutionFctPt exact_flux_pt,
      double& error,
      double& norm)
    {
      std::string error_message =
        "compute_exact_Z2_error undefined for this element \n";
      outfile << error_message;
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void get_Z2_compound_flux_indices(Vector<unsigned>& flux_index) {}
 
    virtual unsigned nvertex_node() const = 0;
 
    virtual Node* vertex_node_pt(const unsigned& j) const = 0;
 
    virtual unsigned nrecovery_order() = 0;
 
    virtual double geometric_jacobian(const Vector<double>& x)
    {
      return 1.0;
    }
  };
 
 
  //========================================================================
  //========================================================================
  class Z2ErrorEstimator : public virtual ErrorEstimator
  {
  public:
    typedef double (*CombinedErrorEstimateFctPt)(const Vector<double>& errors);
 
    Z2ErrorEstimator(const unsigned& recovery_order)
      : Recovery_order(recovery_order),
        Recovery_order_from_first_element(false),
        Reference_flux_norm(0.0),
        Combined_error_fct_pt(0)
    {
    }
 
 
    Z2ErrorEstimator()
      : Recovery_order(0),
        Recovery_order_from_first_element(true),
        Reference_flux_norm(0.0),
        Combined_error_fct_pt(0)
    {
    }
 
    Z2ErrorEstimator(const Z2ErrorEstimator&) = delete;
 
    void operator=(const Z2ErrorEstimator&) = delete;
 
    virtual ~Z2ErrorEstimator() {}
 
    void get_element_errors(Mesh*& mesh_pt, Vector<double>& elemental_error)
    {
      // Create dummy doc info object and switch off output
      DocInfo doc_info;
      doc_info.disable_doc();
      // Forward call to version with doc.
      get_element_errors(mesh_pt, elemental_error, doc_info);
    }
 
    void get_element_errors(Mesh*& mesh_pt,
                            Vector<double>& elemental_error,
                            DocInfo& doc_info);
 
    unsigned& recovery_order()
    {
      return Recovery_order;
    }
 
    unsigned recovery_order() const
    {
      return Recovery_order;
    }
 
    CombinedErrorEstimateFctPt& combined_error_fct_pt()
    {
      return Combined_error_fct_pt;
    }
 
    CombinedErrorEstimateFctPt combined_error_fct_pt() const
    {
      return Combined_error_fct_pt;
    }
 
    void setup_patches(Mesh*& mesh_pt,
                       std::map<Node*, Vector<ElementWithZ2ErrorEstimator*>*>&
                         adjacent_elements_pt,
                       Vector<Node*>& vertex_node_pt);
 
    double& reference_flux_norm()
    {
      return Reference_flux_norm;
    }
 
    double reference_flux_norm() const
    {
      return Reference_flux_norm;
    }
 
    double get_combined_error_estimate(const Vector<double>& compound_error);
 
  private:
    void get_recovered_flux_in_patch(
      const Vector<ElementWithZ2ErrorEstimator*>& patch_el_pt,
      const unsigned& num_recovery_terms,
      const unsigned& num_flux_terms,
      const unsigned& dim,
      DenseMatrix<double>*& recovered_flux_coefficient_pt);
 
 
    unsigned nrecovery_terms(const unsigned& dim);
 
 
    void shape_rec(const Vector<double>& x,
                   const unsigned& dim,
                   Vector<double>& psi_r);
 
    Integral* integral_rec(const unsigned& dim, const bool& is_q_mesh);
 
    unsigned Recovery_order;
 
    bool Recovery_order_from_first_element;
 
    void doc_flux(Mesh* mesh_pt,
                  const unsigned& num_flux_terms,
                  MapMatrixMixed<Node*, int, double>& rec_flux_map,
                  const Vector<double>& elemental_error,
                  DocInfo& doc_info);
 
    double Reference_flux_norm;
 
    CombinedErrorEstimateFctPt Combined_error_fct_pt;
  };
 
 
 
 
  //========================================================================
  //========================================================================
  class DummyErrorEstimator : public virtual ErrorEstimator
  {
  public:
    DummyErrorEstimator(Mesh* mesh_pt,
                        const Vector<unsigned>& elements_to_refine,
                        const unsigned& central_node_number,
                        const bool& use_lagrangian_coordinates = false)
      : Use_lagrangian_coordinates(use_lagrangian_coordinates),
        Central_node_number(central_node_number)
    {
#ifdef PARANOID
#ifdef OOMPH_HAS_MPI
      if (mesh_pt->is_mesh_distributed())
      {
        throw OomphLibError(
          "Can't use this error estimator on distributed meshes!",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
#endif
 
#ifdef PARANOID
      if (mesh_pt->nelement() == 0)
      {
        throw OomphLibError(
          "Can't build error estimator if there are no elements in mesh\n",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      unsigned dim = mesh_pt->finite_element_pt(0)->node_pt(0)->ndim();
      if (use_lagrangian_coordinates)
      {
        SolidNode* solid_nod_pt =
          dynamic_cast<SolidNode*>(mesh_pt->finite_element_pt(0)->node_pt(0));
        if (solid_nod_pt != 0)
        {
          dim = solid_nod_pt->nlagrangian();
        }
      }
      unsigned nregion = elements_to_refine.size();
      Region_low_bound.resize(nregion);
      Region_upp_bound.resize(nregion);
      for (unsigned e = 0; e < nregion; e++)
      {
        Region_low_bound[e].resize(dim, 1.0e20);
        Region_upp_bound[e].resize(dim, -1.0e20);
        FiniteElement* el_pt =
          mesh_pt->finite_element_pt(elements_to_refine[e]);
        unsigned nnod = el_pt->nnode();
        for (unsigned j = 0; j < nnod; j++)
        {
          Node* nod_pt = el_pt->node_pt(j);
          for (unsigned i = 0; i < dim; i++)
          {
            double x = nod_pt->x(i);
            if (use_lagrangian_coordinates)
            {
              SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
              if (solid_nod_pt != 0)
              {
                x = solid_nod_pt->xi(i);
              }
            }
            if (x < Region_low_bound[e][i])
            {
              Region_low_bound[e][i] = x;
            }
            if (x > Region_upp_bound[e][i])
            {
              Region_upp_bound[e][i] = x;
            }
          }
        }
      }
    }
 
 
    DummyErrorEstimator(Mesh* mesh_pt,
                        const Vector<double>& lower_left,
                        const Vector<double>& upper_right,
                        const unsigned& central_node_number,
                        const bool& use_lagrangian_coordinates = false)
      : Use_lagrangian_coordinates(use_lagrangian_coordinates),
        Central_node_number(central_node_number)
    {
#ifdef PARANOID
      if (mesh_pt->nelement() == 0)
      {
        throw OomphLibError(
          "Can't build error estimator if there are no elements in mesh\n",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      unsigned nregion = 1;
      Region_low_bound.resize(nregion);
      Region_upp_bound.resize(nregion);
      Region_low_bound[0] = lower_left;
      Region_upp_bound[0] = upper_right;
    }
 
    DummyErrorEstimator(const DummyErrorEstimator&) = delete;
 
    void operator=(const DummyErrorEstimator&) = delete;
 
    virtual ~DummyErrorEstimator() {}
 
 
    virtual void get_element_errors(Mesh*& mesh_pt,
                                    Vector<double>& elemental_error,
                                    DocInfo& doc_info)
    {
#ifdef PARANOID
      if (doc_info.is_doc_enabled())
      {
        std::ostringstream warning_stream;
        warning_stream
          << "No output defined in DummyErrorEstimator::get_element_errors()\n"
          << "Ignoring doc_info flag.\n";
        OomphLibWarning(warning_stream.str(),
                        "DummyErrorEstimator::get_element_errors()",
                        OOMPH_EXCEPTION_LOCATION);
      }
#endif
      unsigned nregion = Region_low_bound.size();
      unsigned nelem = mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        elemental_error[e] = 0.0;
 
        // Check if element is in the regions to be refined
        // (based on coords of its central node)
        Node* nod_pt =
          mesh_pt->finite_element_pt(e)->node_pt(Central_node_number);
        for (unsigned r = 0; r < nregion; r++)
        {
          bool is_inside = true;
          unsigned dim = Region_low_bound[r].size();
          for (unsigned i = 0; i < dim; i++)
          {
            double x = nod_pt->x(i);
            if (Use_lagrangian_coordinates)
            {
              SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
              if (solid_nod_pt != 0)
              {
                x = solid_nod_pt->xi(i);
              }
            }
            if (x < Region_low_bound[r][i])
            {
              is_inside = false;
              break;
            }
            if (x > Region_upp_bound[r][i])
            {
              is_inside = false;
              break;
            }
          }
          if (is_inside)
          {
            elemental_error[e] = 1.0;
            break;
          }
        }
      }
    }
 
  private:
    bool Use_lagrangian_coordinates;
 
    unsigned Central_node_number;
 
    Vector<Vector<double>> Region_upp_bound;
 
    Vector<Vector<double>> Region_low_bound;
  };
 
 
} // namespace oomph
 
#endif