macro_element.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_MACROELEMENT_HEADER
#define OOMPH_MACROELEMENT_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
// oomph-lib headers
#include "Vector.h"
#include "oomph_utilities.h"
#include "quadtree.h"
#include "octree.h"
 
namespace oomph
{
  class Domain;
 
  //================================================================
  //=================================================================
  class MacroElement
  {
  public:
    MacroElement(Domain* domain_pt, const unsigned& macro_element_number)
      : Domain_pt(domain_pt), Macro_element_number(macro_element_number)
    {
#ifdef LEAK_CHECK
      LeakCheckNames::MacroElement_build += 1;
#endif
    }
 
    MacroElement()
    {
      throw OomphLibError("Don't call empty constructor for MacroElement!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    MacroElement(const MacroElement& dummy) = delete;
 
    void operator=(const MacroElement&) = delete;
 
    virtual ~MacroElement()
    {
#ifdef LEAK_CHECK
      LeakCheckNames::MacroElement_build -= 1;
#endif
    }
 
 
    void output(std::ostream& outfile, const int& nplot)
    {
      unsigned t = 0;
      output(t, outfile, nplot);
    }
 
 
    virtual void output(const unsigned& t,
                        std::ostream& outfile,
                        const unsigned& nplot) = 0;
 
 
    void macro_map(const Vector<double>& s, Vector<double>& r)
    {
      // Evaluate at current timestep
      unsigned t = 0;
      macro_map(t, s, r);
    }
 
 
    virtual void macro_map(const unsigned& t,
                           const Vector<double>& s,
                           Vector<double>& r) = 0;
 
 
    virtual void macro_map(const double& t,
                           const Vector<double>& s,
                           Vector<double>& r)
    {
      // Create an output stream
      std::ostringstream error_message_stream;
 
      // Create an error message
      error_message_stream << "The function macro_map(...) is broken virtual\n"
                           << "If you need it, please implement it!"
                           << std::endl;
 
      // Throw an error
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    } // End of macro_map
 
 
    virtual void output_macro_element_boundaries(std::ostream& outfile,
                                                 const unsigned& nplot) = 0;
 
 
    virtual void assemble_macro_to_eulerian_jacobian(
      const unsigned& t, const Vector<double>& s, DenseMatrix<double>& jacobian)
    {
      // error message stream
      std::ostringstream error_message;
      error_message << "assemble_macro_to_eulerian_jacobian(...) not \n"
                    << "implemented for this element\n"
                    << std::endl;
      // throw error
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void assemble_macro_to_eulerian_jacobian2(
      const unsigned& t,
      const Vector<double>& s,
      DenseMatrix<double>& jacobian2)
    {
      // error message stream
      std::ostringstream error_message;
      error_message << "assemble_macro_to_eulerian_jacobian2(...) not \n"
                    << "implemented for this element\n"
                    << std::endl;
      // throw error
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void assemble_macro_to_eulerian_jacobian(const Vector<double>& s,
                                             DenseMatrix<double>& jacobian)
    {
      unsigned t = 0;
      assemble_macro_to_eulerian_jacobian(t, s, jacobian);
    }
 
 
    void assemble_macro_to_eulerian_jacobian2(const Vector<double>& s,
                                              DenseMatrix<double>& jacobian2)
    {
      unsigned t = 0;
      assemble_macro_to_eulerian_jacobian2(t, s, jacobian2);
    }
 
    unsigned& macro_element_number()
    {
      return Macro_element_number;
    }
 
    Domain*& domain_pt()
    {
      return Domain_pt;
    }
 
  protected:
    Domain* Domain_pt;
 
    unsigned Macro_element_number;
  };
 
 
 
 
  //================================================================
  //=================================================================
  template<int DIM>
  class QMacroElement : public MacroElement
  {
  };
 
 
  //================================================================
  //=================================================================
  template<>
  class QMacroElement<2> : public MacroElement
  {
  public:
    QMacroElement(Domain* domain_pt, const unsigned& macro_element_number)
      : MacroElement(domain_pt, macro_element_number){};
 
 
    QMacroElement()
    {
      throw OomphLibError("Don't call empty constructor for QMacroElement!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    QMacroElement(const QMacroElement& dummy) = delete;
 
    void operator=(const QMacroElement&) = delete;
 
    virtual ~QMacroElement(){};
 
    void output(const unsigned& t, std::ostream& outfile, const unsigned& nplot)
    {
      Vector<double> x(2), f(2);
      outfile << "ZONE I=" << nplot << ", J=" << nplot << std::endl;
      for (unsigned i = 0; i < nplot; i++)
      {
        x[1] = -1.0 + 2.0 * double(i) / double(nplot - 1);
        for (unsigned j = 0; j < nplot; j++)
        {
          x[0] = -1.0 + 2.0 * double(j) / double(nplot - 1);
          macro_map(t, x, f);
          outfile << f[0] << " " << f[1] << std::endl;
        }
      }
    }
 
 
    void output_macro_element_boundaries(std::ostream& outfile,
                                         const unsigned& nplot);
 
    void macro_map(const unsigned& t,
                   const Vector<double>& S,
                   Vector<double>& r);
 
 
    void macro_map(const double& t, const Vector<double>& s, Vector<double>& r);
 
 
    virtual void assemble_macro_to_eulerian_jacobian(
      const unsigned& t,
      const Vector<double>& s,
      DenseMatrix<double>& jacobian);
 
 
    virtual void assemble_macro_to_eulerian_jacobian2(
      const unsigned& t,
      const Vector<double>& s,
      DenseMatrix<double>& jacobian2);
  };
 
 
  //================================================================
  //=================================================================
  template<>
  class QMacroElement<3> : public MacroElement
  {
  public:
    QMacroElement(Domain* domain_pt, const unsigned& macro_element_number)
      : MacroElement(domain_pt, macro_element_number){};
 
    QMacroElement()
    {
      throw OomphLibError("Don't call empty constructor for QMacroElement!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    QMacroElement(const QMacroElement& dummy) = delete;
 
    void operator=(const QMacroElement&) = delete;
 
    virtual ~QMacroElement(){};
 
 
    void output(const unsigned& t, std::ostream& outfile, const unsigned& nplot)
    {
      Vector<double> x(3), f(3);
 
      outfile << "ZONE I=" << nplot << ", J=" << nplot << ", k=" << nplot
              << std::endl;
      for (unsigned i = 0; i < nplot; i++)
      {
        x[2] = -1.0 + 2.0 * double(i) / double(nplot - 1);
 
        for (unsigned j = 0; j < nplot; j++)
        {
          x[1] = -1.0 + 2.0 * double(j) / double(nplot - 1);
 
          for (unsigned k = 0; k < nplot; k++)
          {
            x[0] = -1.0 + 2.0 * double(k) / double(nplot - 1);
 
            macro_map(t, x, f);
 
            outfile << f[0] << " " << f[1] << " " << f[2] << std::endl;
          }
        }
      }
    }
 
 
    void output_macro_element_boundaries(std::ostream& outfile,
                                         const unsigned& nplot);
 
    void macro_map(const unsigned& t,
                   const Vector<double>& S,
                   Vector<double>& r);
  };
 
} // namespace oomph
 
#endif