CD_ComputationalGeometry.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_COMPUTATIONALGEOMETRY_H
#define CD_COMPUTATIONALGEOMETRY_H
 
// Std includes
#include <memory>
#include <string>
 
// EBGeometry includes
#include <EBGeometry.hpp>
 
// Chombo includes
#include <BaseIF.H>
#include <GeometryService.H>
#include <MFIndexSpace.H>
#include <Box.H>
#include <RealVect.H>
#include <ProblemDomain.H>
 
// Our includes
#include <CD_Dielectric.H>
#include <CD_Electrode.H>
#include <CD_MultiFluidIndexSpace.H>
#include <CD_RealBox.H>
#include <CD_NamespaceHeader.H>
 
class ComputationalGeometry
{
public:
  static constexpr int K = 4;
 
  using Vec3 = EBGeometry::Vec3T<Real>;
 
  using BV = EBGeometry::BoundingVolumes::AABBT<Real>;
 
  using BoxTree = EBGeometry::BVH::PackedBVH<Real, int, K, EBGeometry::BVH::ValueStorage<int>>;
 
  ComputationalGeometry();
 
  virtual ~ComputationalGeometry();
 
  const Vector<Dielectric>&
  getDielectrics() const;
 
  const Vector<Electrode>&
  getElectrodes() const;
 
  Real
  getGasPermittivity() const;
 
  void
  useScanShop(const ProblemDomain& a_beginDomain);
 
  void
  useChomboShop();
 
  void
  usePolyhedralShop(const ProblemDomain& a_beginDomain);
 
  void
  setDielectrics(const Vector<Dielectric>& a_dielectrics);
 
  void
  setElectrodes(const Vector<Electrode>& a_electrodes);
 
  void
  setGasPermittivity(const Real a_eps0);
 
  const RefCountedPtr<MultiFluidIndexSpace>&
  getMfIndexSpace() const;
 
  const RefCountedPtr<BaseIF>&
  getGasImplicitFunction() const;
 
  const RefCountedPtr<BaseIF>&
  getSolidImplicitFunction() const;
 
  const RefCountedPtr<BaseIF>&
  getImplicitFunction(const phase::which_phase a_phase) const;
 
  virtual void
  makeGrids(const ProblemDomain& a_startDomain,
            const ProblemDomain& a_stopDomain,
            const RealVect&      a_probLo,
            const Real           a_startDx,
            const Real           a_refineAngle,
            const int            a_maxGhostEB,
            const int            a_minBlockSize,
            const int            a_maxBlockSize);
 
  int
  getNumGridLevels() const noexcept;
 
  GeometryService::InOut
  classify(const Box& a_box, const int a_level, const phase::which_phase a_phase) const;
 
  GeometryService::InOut
  classify(const Box& a_box, const ProblemDomain& a_domain, const phase::which_phase a_phase) const;
 
  int
  getStartLevel() const noexcept;
 
  int
  getLevel(const ProblemDomain& a_domain) const noexcept;
 
  Vector<Box>
  getBoxes(const phase::which_phase a_phase, const int a_level, const GeometryService::InOut a_type) const noexcept;
 
  const Vector<Box>&
  getBoxes(const int a_level) const noexcept;
 
  const Vector<Box>&
  getCutTiles(const int a_level) const noexcept;
 
  const Vector<Box>&
  getSplitBoxes(const int a_level, Vector<int>& a_reasons) const noexcept;
 
  const Vector<GeometryService::InOut>&
  getTypes(const phase::which_phase a_phase, const int a_level) const noexcept;
 
  const ProblemDomain&
  getDomain(const int a_level) const noexcept;
 
  Real
  getDx(const int a_level) const noexcept;
 
  virtual void
  buildGeometries(const ProblemDomain& a_finestDomain,
                  const RealVect&      a_probLo,
                  const Real           a_finestDx,
                  const int            a_nCellMax,
                  const int            a_maxGhostEB,
                  const int            a_maxCoarsen = -1);
 
protected:
  constexpr static Real s_thresh = 1.E-15;
 
  constexpr static bool s_strictGeometry = true;
 
  enum class Generator
  {
    GeometryShop,
    ScanShop,
    PolyhedralShop
  };
 
  Generator m_generator;
 
  RefCountedPtr<MultiFluidIndexSpace> m_multifluidIndexSpace;
 
  RefCountedPtr<BaseIF> m_implicitFunctionGas;
 
  RefCountedPtr<BaseIF> m_implicitFunctionSolid;
 
  Vector<Dielectric> m_dielectrics;
 
  Vector<Electrode> m_electrodes;
 
  ProblemDomain m_scanDomain;
 
  RealVect m_probLo;
 
  Real m_eps0;
 
  Real m_refineAngle;
 
  bool m_refineSaddles;
 
  int m_maxGhostEB;
 
  int m_startLevel;
 
  int m_stopLevel;
 
  int m_minBlockSize;
 
  int m_maxBlockSize;
 
  bool m_profile;
 
  bool m_verbose;
 
  Vector<ProblemDomain> m_domains;
 
  Vector<Real> m_dx;
 
  Vector<Vector<Box>> m_cutTiles;
 
  Vector<Vector<Box>> m_boxes;
 
  Vector<std::shared_ptr<BoxTree>> m_boxTrees;
 
  Vector<std::shared_ptr<BoxTree>> m_tileTrees;
 
  Vector<Vector<int>> m_splitCounts;
 
  Vector<Vector<Box>> m_splitBoxes;
 
  Vector<Vector<int>> m_splitReasons;
 
  Vector<Vector<GeometryService::InOut>> m_gasTypes;
 
  Vector<Vector<GeometryService::InOut>> m_solidTypes;
 
  void
  reportGrids() const;
 
  void
  buildImplicitFunctions();
 
  Vector<Vector<GeometryService::InOut>>&
  types(const phase::which_phase a_phase) noexcept;
 
  const Vector<Vector<GeometryService::InOut>>&
  types(const phase::which_phase a_phase) const noexcept;
 
  void
  buildStartLevel();
 
  void
  buildFinerLevels(Vector<Vector<int>>& a_firstChild, Vector<Vector<int>>& a_numChildren);
 
  void
  classifyBoxes(const Vector<Box>&              a_boxes,
                const int                       a_level,
                Vector<GeometryService::InOut>& a_gasTypes,
                Vector<GeometryService::InOut>& a_solidTypes) const;
 
  Vector<int>
  splitFlags(const Vector<Box>&                    a_boxes,
             const int                             a_level,
             const Vector<GeometryService::InOut>& a_gasTypes,
             const Vector<GeometryService::InOut>& a_solidTypes) const;
 
  GeometryService::InOut
  classifyBox(const Box& a_box, const int a_level, const phase::which_phase a_phase) const;
 
  enum class SplitReason
  {
    None,
    Interior,
    Ring,
    Medial,
    DoubleCrossing,
    Twist
  };
 
  bool
  hasTwistedPatch(const Box& a_box, const int a_level, const phase::which_phase a_phase) const;
 
  SplitReason
  exceedsCurvature(const Box& a_box, const int a_level, const phase::which_phase a_phase) const;
 
  bool
  doublyCrossedEdge(const Box& a_box, const int a_level, const phase::which_phase a_phase) const;
 
  void
  makeTiles();
 
  void
  classifyTiles(const Vector<Vector<int>>&              a_firstChild,
                const Vector<Vector<int>>&              a_numChildren,
                Vector<Vector<GeometryService::InOut>>& a_gasTileTypes,
                Vector<Vector<GeometryService::InOut>>& a_solidTileTypes,
                Vector<Vector<int>>&                    a_tileHosts) const;
 
  static constexpr int s_treeLeafSize = 8;
 
  static constexpr int s_maxTilePasses = 8;
 
  static BV
  boundingVolume(const Box& a_box) noexcept;
 
  void
  buildBoxTrees();
 
  std::shared_ptr<BoxTree>
  buildTree(const Vector<Box>& a_boxes) const;
 
  static Vector<int>
  meeting(const std::shared_ptr<BoxTree>& a_tree, const Vector<Box>& a_boxes, const Box& a_box);
 
  bool
  tagUnresolvedSeams(Vector<IntVectSet>& a_tags) const;
 
  void
  buildTileTrees();
 
  void
  buildBoxTree(const int a_level);
 
  Vector<int>
  boxesMeeting(const int a_level, const Box& a_box) const;
 
  Vector<int>
  tilesMeeting(const int a_level, const Box& a_box) const;
 
  int
  containingBox(const int a_level, const IntVect& a_cell) const;
 
  void
  decimateBoxes(const Vector<Vector<GeometryService::InOut>>& a_gasTileTypes,
                const Vector<Vector<GeometryService::InOut>>& a_solidTileTypes,
                const Vector<Vector<int>>&                    a_tileHosts);
 
  void
  buildCoarserLevels();
 
  void
  buildGasGeometry(GeometryService*&    a_geoserver,
                   const ProblemDomain& a_finestDomain,
                   const RealVect&      a_probLo,
                   const Real           a_finestDx);
 
  void
  buildSolidGeometry(GeometryService*&    a_geoserver,
                     const ProblemDomain& a_finestDomain,
                     const RealVect&      a_probLo,
                     const Real           a_finestDx);
};
 
#include <CD_NamespaceFooter.H>
 
#endif