geom_obj_with_boundary.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_GEOM_OBJ_WITH_BOUNDARY_HEADER
#define OOMPH_GEOM_OBJ_WITH_BOUNDARY_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
// oomph-lib headers
#include "geom_objects.h"
 
 
namespace oomph
{
  //===========================================================================
  //===========================================================================
  class DiskLikeGeomObjectWithBoundaries : public virtual GeomObject
  {
  public:
    DiskLikeGeomObjectWithBoundaries() : GeomObject(2, 3) {}
 
    unsigned nboundary() const
    {
      return Boundary_parametrising_geom_object_pt.size();
    }
 
    void position_on_boundary(const unsigned& b,
                              const double& zeta_bound,
                              Vector<double>& r) const
    {
      Vector<double> zeta(2);
      zeta_on_boundary(b, zeta_bound, zeta);
      position(zeta, r);
    }
 
    void zeta_on_boundary(const unsigned& b,
                          const double& zeta_bound,
                          Vector<double>& zeta) const
    {
#ifdef PARANOID
      if (Boundary_parametrising_geom_object_pt[b] == 0)
      {
        std::ostringstream error_message;
        error_message << "You must create a <1,2> Geom Object that provides a\n"
                      << "mapping from the 1D boundary coordinate to the 2D\n"
                      << "intrinsic Lagrangian coordinate of disk-like object\n"
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (Boundary_parametrising_geom_object_pt[b]->nlagrangian() != 1)
      {
        std::ostringstream error_message;
        error_message << "Boundary_parametrising_geom_object_pt must point to\n"
                      << "GeomObject with one Lagrangian coordinate. Yours has "
                      << Boundary_parametrising_geom_object_pt[b]->nlagrangian()
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      if (Boundary_parametrising_geom_object_pt[b]->ndim() != 2)
      {
        std::ostringstream error_message;
        error_message << "Boundary_parametrising_geom_object_pt must point to\n"
                      << "GeomObject with two Eulerian coordinates. Yours has "
                      << Boundary_parametrising_geom_object_pt[b]->ndim()
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      Vector<double> zeta_bound_vector(1, zeta_bound);
      Boundary_parametrising_geom_object_pt[b]->position(zeta_bound_vector,
                                                         zeta);
    }
 
    GeomObject* boundary_parametrising_geom_object_pt(const unsigned& b) const
    {
      return Boundary_parametrising_geom_object_pt[b];
    }
 
    double zeta_boundary_start(const unsigned& b) const
    {
      return Zeta_boundary_start[b];
    }
 
    double zeta_boundary_end(const unsigned& b) const
    {
      return Zeta_boundary_end[b];
    }
 
 
    virtual void boundary_triad(const unsigned& b,
                                const double& zeta_bound,
                                Vector<double>& r,
                                Vector<double>& tangent,
                                Vector<double>& normal,
                                Vector<double>& binormal)
    {
      std::ostringstream error_message;
      error_message << "Broken virtual function; please implement for your\n"
                    << "derived version of this class!" << std::endl;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    void output_boundaries(const unsigned& nplot,
                           std::ofstream& two_d_boundaries_file,
                           std::ofstream& three_d_boundaries_file)
    {
      std::ofstream boundaries_tangent_file;
      std::ofstream boundaries_normal_file;
      std::ofstream boundaries_binormal_file;
      output_boundaries_and_triads(nplot,
                                   two_d_boundaries_file,
                                   three_d_boundaries_file,
                                   boundaries_tangent_file,
                                   boundaries_normal_file,
                                   boundaries_binormal_file);
    }
 
 
    void output_boundaries_and_triads(const unsigned& nplot,
                                      std::ofstream& two_d_boundaries_file,
                                      std::ofstream& three_d_boundaries_file,
                                      std::ofstream& boundaries_tangent_file,
                                      std::ofstream& boundaries_normal_file,
                                      std::ofstream& boundaries_binormal_file)
    {
      Vector<double> r(3);
      Vector<double> zeta(2);
      double zeta_bound = 0.0;
      Vector<double> tangent(3);
      Vector<double> normal(3);
      Vector<double> binormal(3);
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        two_d_boundaries_file << "ZONE" << std::endl;
        three_d_boundaries_file << "ZONE" << std::endl;
        boundaries_tangent_file << "ZONE" << std::endl;
        boundaries_normal_file << "ZONE" << std::endl;
        boundaries_binormal_file << "ZONE" << std::endl;
 
        double zeta_min = zeta_boundary_start(b);
        double zeta_max = zeta_boundary_end(b);
        unsigned n = 100;
        for (unsigned i = 0; i < n; i++)
        {
          zeta_bound =
            zeta_min + (zeta_max - zeta_min) * double(i) / double(n - 1);
          position_on_boundary(b, zeta_bound, r);
          zeta_on_boundary(b, zeta_bound, zeta);
          boundary_triad(b, zeta_bound, r, tangent, normal, binormal);
 
          two_d_boundaries_file << zeta[0] << " " << zeta[1] << " "
                                << zeta_bound << " " << std::endl;
 
          three_d_boundaries_file << r[0] << " " << r[1] << " " << r[2] << " "
                                  << zeta[0] << " " << zeta[1] << " "
                                  << zeta_bound << " " << std::endl;
 
          boundaries_tangent_file << r[0] << " " << r[1] << " " << r[2] << " "
                                  << tangent[0] << " " << tangent[1] << " "
                                  << tangent[2] << " " << std::endl;
 
          boundaries_normal_file << r[0] << " " << r[1] << " " << r[2] << " "
                                 << normal[0] << " " << normal[1] << " "
                                 << normal[2] << " " << std::endl;
 
          boundaries_binormal_file << r[0] << " " << r[1] << " " << r[2] << " "
                                   << binormal[0] << " " << binormal[1] << " "
                                   << binormal[2] << " " << std::endl;
        }
      }
    }
 
 
    void add_region_coordinates(const unsigned& r,
                                Vector<double>& zeta_in_region)
    {
      // Verify if not using the default region number (zero)
      if (r == 0)
      {
        std::ostringstream error_message;
        error_message
          << "Please use another region id different from zero.\n"
          << "It is internally used as the default region number.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      // Need two coordinates
      unsigned n = zeta_in_region.size();
      if (n != 2)
      {
        std::ostringstream error_message;
        error_message << "Vector specifying zeta coordinates of point in\n"
                      << "region has be length 2; yours has length " << n
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
      // First check if the region with the specified id does not already exist
      std::map<unsigned, Vector<double>>::iterator it;
      it = Zeta_in_region.find(r);
 
      // If it is already a region defined with that id throw an error
      if (it != Zeta_in_region.end())
      {
        std::ostringstream error_message;
        error_message << "The region id (" << r << ") that you are using for"
                      << "defining\n"
                      << "your region is already in use. Use another\n"
                      << "region id and verify that you are not re-using\n"
                      << " previously defined regions ids.\n"
                      << std::endl;
        OomphLibWarning(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      // If it does not exist then create the map
      Zeta_in_region[r] = zeta_in_region;
    }
 
    std::map<unsigned, Vector<double>> zeta_in_region() const
    {
      return Zeta_in_region;
    }
 
  protected:
    Vector<double> Zeta_boundary_start;
 
    Vector<double> Zeta_boundary_end;
 
    Vector<GeomObject*> Boundary_parametrising_geom_object_pt;
 
    std::map<unsigned, Vector<double>> Zeta_in_region;
  };
 
 
 
 
  //=========================================================================
  //=========================================================================
  class WarpedCircularDisk : public DiskLikeGeomObjectWithBoundaries
  {
  public:
    WarpedCircularDisk(const double& epsilon,
                       const unsigned& n,
                       const double& z_offset = 0.0)
      : Epsilon(epsilon), N(n), Z_offset(z_offset)
    {
      // How many boundaries do we have?
      unsigned nb = 2;
      Boundary_parametrising_geom_object_pt.resize(nb);
      Zeta_boundary_start.resize(nb, 0.0);
      Zeta_boundary_end.resize(nb, 0.0);
 
      // GeomObject<1,2> representing the first boundary
      Boundary_parametrising_geom_object_pt[0] = new Ellipse(1.0, 1.0);
      Zeta_boundary_start[0] = 0.0;
      Zeta_boundary_end[0] = MathematicalConstants::Pi;
 
      // GeomObject<1,2> representing the second boundary
      Boundary_parametrising_geom_object_pt[1] = new Ellipse(1.0, 1.0);
      Zeta_boundary_start[1] = MathematicalConstants::Pi;
      Zeta_boundary_end[1] = 2.0 * MathematicalConstants::Pi;
    }
 
    WarpedCircularDisk()
    {
      throw OomphLibError("Don't call default constructor!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    WarpedCircularDisk(const WarpedCircularDisk& dummy) = delete;
 
    void operator=(const WarpedCircularDisk&) = delete;
 
    virtual ~WarpedCircularDisk()
    {
      unsigned n = nboundary();
      for (unsigned b = 0; b < n; b++)
      {
        delete Boundary_parametrising_geom_object_pt[b];
        Boundary_parametrising_geom_object_pt[b] = 0;
      }
    }
 
    double& epsilon()
    {
      return Epsilon;
    }
 
    void position(const Vector<double>& zeta, Vector<double>& r) const
    {
      // Position Vector
      r[0] = zeta[0];
      r[1] = zeta[1];
      double radius = sqrt(r[0] * r[0] + r[1] * r[1]);
      double phi = atan2(r[1], r[0]);
      r[2] = Z_offset + w(radius, phi);
    }
 
 
    void position(const unsigned& t,
                  const Vector<double>& zeta,
                  Vector<double>& r) const
    {
      position(zeta, r);
    }
 
 
    void boundary_triad(const unsigned& b,
                        const double& zeta_bound,
                        Vector<double>& r,
                        Vector<double>& tangent,
                        Vector<double>& normal,
                        Vector<double>& binormal)
    {
      double phi = zeta_bound;
 
      // Position Vector
      r[0] = cos(phi);
      r[1] = sin(phi);
      r[2] = Z_offset + w(1.0, phi);
 
      Vector<double> dr_dr(3);
      dr_dr[0] = cos(phi);
      dr_dr[1] = sin(phi);
      dr_dr[2] = dwdr(1.0, phi);
      double inv_norm = 1.0 / sqrt(dr_dr[0] * dr_dr[0] + dr_dr[1] * dr_dr[1] +
                                   dr_dr[2] * dr_dr[2]);
 
      normal[0] = dr_dr[0] * inv_norm;
      normal[1] = dr_dr[1] * inv_norm;
      normal[2] = dr_dr[2] * inv_norm;
 
      Vector<double> dr_dphi(3);
      dr_dphi[0] = -sin(phi);
      dr_dphi[1] = cos(phi);
      dr_dphi[2] = dwdphi(1.0, phi);
 
      inv_norm = 1.0 / sqrt(dr_dphi[0] * dr_dphi[0] + dr_dphi[1] * dr_dphi[1] +
                            dr_dphi[2] * dr_dphi[2]);
 
      tangent[0] = dr_dphi[0] * inv_norm;
      tangent[1] = dr_dphi[1] * inv_norm;
      tangent[2] = dr_dphi[2] * inv_norm;
 
      binormal[0] = tangent[1] * normal[2] - tangent[2] * normal[1];
      binormal[1] = tangent[2] * normal[0] - tangent[0] * normal[2];
      binormal[2] = tangent[0] * normal[1] - tangent[1] * normal[0];
    }
 
  private:
    double w(const double& r, const double& phi) const
    {
      return Epsilon * cos(double(N) * phi) * pow(r, 2);
    }
 
    double dwdr(const double& r, const double& phi) const
    {
      return Epsilon * cos(double(N) * phi) * 2.0 * r;
    }
 
    double dwdphi(const double& r, const double& phi) const
    {
      return -Epsilon * double(N) * sin(double(N) * phi) * pow(r, 2);
    }
 
    double Epsilon;
 
    unsigned N;
 
    double Z_offset;
  };
 
 
 
  //=========================================================================
  //=========================================================================
  class WarpedCircularDiskWithAnnularInternalBoundary
    : public virtual WarpedCircularDisk
  {
  public:
    WarpedCircularDiskWithAnnularInternalBoundary(const double& h_annulus,
                                                  const double& epsilon,
                                                  const unsigned& n,
                                                  const double& z_offset = 0.0)
      : WarpedCircularDisk(epsilon, n, z_offset), H_annulus(h_annulus)
    {
      // We have two more boundaries!
      Boundary_parametrising_geom_object_pt.resize(4);
      Zeta_boundary_start.resize(4);
      Zeta_boundary_end.resize(4);
 
      // Radius of the internal annular boundary
      double r_annulus = 1.0 - h_annulus;
 
      // GeomObject<1,2> representing the third boundary
      Boundary_parametrising_geom_object_pt[2] =
        new Ellipse(r_annulus, r_annulus);
      Zeta_boundary_start[2] = 0.0;
      Zeta_boundary_end[2] = MathematicalConstants::Pi;
 
      // GeomObject<1,2> representing the fourth boundary
      Boundary_parametrising_geom_object_pt[3] =
        new Ellipse(r_annulus, r_annulus);
      Zeta_boundary_start[3] = MathematicalConstants::Pi;
      Zeta_boundary_end[3] = 2.0 * MathematicalConstants::Pi;
 
      // Region 1 is the annular region; identify it by a point in
      // this region.
      unsigned r = 1;
      Vector<double> zeta_in_region(2);
      zeta_in_region[0] = 0.0;
      zeta_in_region[0] = 1.0 - 0.5 * h_annulus;
      add_region_coordinates(r, zeta_in_region);
    }
 
    WarpedCircularDiskWithAnnularInternalBoundary(
      const WarpedCircularDiskWithAnnularInternalBoundary& dummy) = delete;
 
    void operator=(const WarpedCircularDiskWithAnnularInternalBoundary&) =
      delete;
 
    virtual ~WarpedCircularDiskWithAnnularInternalBoundary() {}
 
 
    double h_annulus() const
    {
      return H_annulus;
    }
 
  protected:
    double H_annulus;
  };
 
 
} // namespace oomph
 
#endif