GFace.h

// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.

#ifndef _GFACE_H_
#define _GFACE_H_

#include <list>
#include <string>
#include <vector>
#include <map>
#include "GEntity.h"
#include "GPoint.h"
#include "GEdgeLoop.h"
#include "SPoint2.h"
#include "SVector3.h"
#include "Pair.h"
#include "Numeric.h"
#include "boundaryLayersData.h"

class MElement;
class MTriangle;
class MQuadrangle;
class MPolygon;
class ExtrudeParams;

class GRegion;

// A model face.
class GFace : public GEntity {
protected:
  // edge loops might replace what follows (list of all the edges of
  // the face + directions)
  std::vector<GEdge *> l_edges;
  std::vector<int> l_dirs;
  GRegion *r1, *r2;
  mean_plane meanPlane;
  std::vector<GEdge *> embedded_edges;
  std::set<GVertex *, GEntityLessThan> embedded_vertices;

  BoundaryLayerColumns _columns;

public: // this will become protected or private
  std::list<GEdgeLoop> edgeLoops;

  // periodic counterparts of edges
  std::map<GEdge *, std::pair<GEdge *, int> > edgeCounterparts;

  // specify mesh master with transformation, deduce edgeCounterparts
  void setMeshMaster(GFace *master, const std::vector<double> &);

  // specify mesh master and edgeCounterparts, deduce transformation
  void setMeshMaster(GFace *master, const std::map<int, int> &);

  // align elements with mesh master
  void alignElementsWithMaster();

public:
  GFace(GModel *model, int tag);
  virtual ~GFace();

  // delete mesh data
  virtual void deleteMesh(bool onlyDeleteElements = false);

  // add/delete regions that are bounded by the face
  void addRegion(GRegion *r)
  {
    if(r == r1 || r == r2) return;
    r1 ? r2 = r : r1 = r;
  }
  void delRegion(GRegion *r)
  {
    if(r1 == r) r1 = r2;
    r2 = NULL;
  }
  GRegion *getRegion(int const num) const { return num == 0 ? r1 : r2; }

  // get number of regions
  int numRegions() const { return (r1 != NULL) + (r2 != NULL); }

  std::list<GRegion *> regions() const
  {
    std::list<GRegion *> r;
    for(int i = 0; i < numRegions(); i++) r.push_back(getRegion(i));
    return r;
  }

  // add embedded vertices/edges
  void addEmbeddedVertex(GVertex *v) { embedded_vertices.insert(v); }
  void addEmbeddedEdge(GEdge *e) { embedded_edges.push_back(e); }

  // delete the edge from the face (the edge is supposed to be a free
  // edge in the face, not part of any edge loops--use with caution!)
  void delFreeEdge(GEdge *e);

  // edge orientations
  virtual std::vector<int> const &orientations() const { return l_dirs; }

  // edges that bound the face
  int delEdge(GEdge *edge);
  virtual std::vector<GEdge *> const &edges() const { return l_edges; }
  void set(const std::vector<GEdge *> &f) { l_edges = f; }
  void setOrientations(const std::vector<int> &f) { l_dirs = f; }
  void setEdge(GEdge *const f, int const orientation)
  {
    l_edges.push_back(f);
    l_dirs.push_back(orientation);
  }

  virtual std::vector<int> const &edgeOrientations() const { return l_dirs; }

  bool containsEdge(int const iEdge) const
  {
    for(std::vector<GEdge *>::const_iterator it = l_edges.begin();
        it != l_edges.end(); ++it)
      if((*it)->tag() == iEdge) return true;
    return false;
  }

  // edges that are embedded in the face
  virtual std::vector<GEdge *> embeddedEdges() const { return embedded_edges; }

  // edges that are embedded in the face
  virtual std::set<GVertex *, GEntityLessThan> embeddedVertices() const
  {
    return embedded_vertices;
  }

  std::vector<MVertex *> getEmbeddedMeshVertices() const;

  // vertices that bound the face
  virtual std::vector<GVertex *> vertices() const;

  // dimension of the face (2)
  virtual int dim() const { return 2; }

  // returns the parent entity for partitioned entities
  virtual GEntity *getParentEntity() { return 0; }

  // set visibility flag
  virtual void setVisibility(char val, bool recursive = false);

  // set color
  virtual void setColor(unsigned int val, bool recursive = false);

  // compute the parameters UV from a point XYZ
  void XYZtoUV(double X, double Y, double Z, double &U, double &V, double relax,
               bool onSurface = true) const;

  // get the bounding box
  virtual SBoundingBox3d bounds(bool fast = false) const;

  // get the oriented bounding box
  virtual SOrientedBoundingBox getOBB();

  // compute the genus G of the surface
  virtual int genusGeom() const;
  virtual bool checkTopology() const { return true; }

  // return the point on the face corresponding to the given parameter
  virtual GPoint point(double par1, double par2) const = 0;
  virtual GPoint point(const SPoint2 &pt) const
  {
    return point(pt.x(), pt.y());
  }

  // if the mapping is a conforming mapping, i.e. a mapping that
  // conserves angles, this function returns the eigenvalue of the
  // metric at a given point this is a special feature for
  // stereographic mappings of the sphere that is used in 2D mesh
  // generation !
  virtual double getMetricEigenvalue(const SPoint2 &);

  // eigen values are absolute values and sorted from min to max of absolute
  // values eigen vectors are the COLUMNS of eigVec
  virtual void getMetricEigenVectors(const SPoint2 &param, double eigVal[2],
                                     double eigVec[4]) const;

  // return the parameter location on the face given a point in space
  // that is on the face
  virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface = true) const;

  // true if the parameter value is interior to the face
  virtual bool containsParam(const SPoint2 &pt);

  // return the point on the face closest to the given point
  virtual GPoint closestPoint(const SPoint3 &queryPoint,
                              const double initialGuess[2]) const;

  // return the normal to the face at the given parameter location
  virtual SVector3 normal(const SPoint2 &param) const;

  // return the first derivate of the face at the parameter location
  virtual Pair<SVector3, SVector3> firstDer(const SPoint2 &param) const = 0;

  // compute the second derivates of the face at the parameter location
  virtual void secondDer(const SPoint2 &param, SVector3 &dudu, SVector3 &dvdv,
                         SVector3 &dudv) const = 0;

  // return the curvature computed as the divergence of the normal
  inline double curvature(const SPoint2 &param) const
  {
    return curvatureMax(param);
  }
  virtual double curvatureDiv(const SPoint2 &param) const;

  // return the maximum curvature at a point
  virtual double curvatureMax(const SPoint2 &param) const;

  // compute the min and max curvatures and the corresponding directions
  // return the max curvature
  virtual double curvatures(const SPoint2 &param, SVector3 &dirMax,
                            SVector3 &dirMin, double &curvMax,
                            double &curvMin) const;

  // return a type-specific additional information string
  virtual std::string getAdditionalInfoString(bool multline = false);

  // export in GEO format
  virtual void writeGEO(FILE *fp);

  // fill the crude representation cross
  virtual bool buildRepresentationCross(bool force = false);

  // build an STL triangulation (or do nothing if it already exists,
  // unless force=true)
  virtual bool buildSTLTriangulation(bool force = false);

  // fill the vertex array using an STL triangulation
  bool fillVertexArray(bool force = false);

  // recompute the mean plane of the surface from a list of points
  void computeMeanPlane(const std::vector<MVertex *> &points);
  void computeMeanPlane(const std::vector<SPoint3> &points);

  // recompute the mean plane of the surface from its bounding vertices
  void computeMeanPlane();

  // get the mean plane info
  void getMeanPlaneData(double VX[3], double VY[3], double &x, double &y,
                        double &z) const;
  void getMeanPlaneData(double plan[3][3]) const;

  // number of types of elements
  int getNumElementTypes() const { return 3; }

  // get total/by-type number of elements in the mesh
  size_type getNumMeshElements() const;
  unsigned int getNumMeshElementsByType(const int familyType) const;
  unsigned int getNumMeshParentElements();
  void getNumMeshElements(unsigned *const c) const;

  // get the start of the array of a type of element
  MElement *const *getStartElementType(int type) const;

  // get the element at the given index
  MElement *getMeshElement(unsigned int index) const;
  // get the element at the given index for a given familyType
  MElement *getMeshElementByType(const int familyType,
                                 const unsigned int index) const;

  // reset the mesh attributes to default values
  virtual void resetMeshAttributes();

  // for periodic faces, move parameters into the range chosen
  // for that face
  void moveToValidRange(SPoint2 &pt) const;

  // compute mesh statistics
  void computeMeshSizeFieldAccuracy(double &avg, double &max_e, double &min_e,
                                    int &nE, int &GS);

  // add points (and optionally normals) in vectors so that two
  // points are at most maxDist apart
  bool fillPointCloud(double maxDist, std::vector<SPoint3> *points,
                      std::vector<SPoint2> *uvpoints = 0,
                      std::vector<SVector3> *normals = 0);

  // tells if it's a sphere, and if it is, returns parameters
  virtual bool isSphere(double &radius, SPoint3 &center) const { return false; }

  // new interface for meshing
  virtual void mesh(bool verbose);

  struct {
    // do we recombine the triangles of the mesh?
    int recombine;
    // what is the treshold angle for recombination
    double recombineAngle;
    // is this surface meshed using a transfinite interpolation
    char method;
    // corners of the transfinite interpolation
    std::vector<GVertex *> corners;
    // all diagonals of the triangulation are left (-1), right (1) or
    // alternated starting at right (2) or left (-2)
    int transfiniteArrangement;
    // do we smooth (transfinite) mesh? (<0 to use default smoothing)
    int transfiniteSmoothing;
    // the extrusion parameters (if any)
    ExtrudeParams *extrude;
    // reverse mesh orientation
    bool reverseMesh;
    // global mesh size constraint for the surface
    double meshSize;
  } meshAttributes;

  int getMeshingAlgo() const;
  void setMeshingAlgo(int) const;
  void unsetMeshingAlgo() const;
  int getCurvatureControlParameter() const;
  void setCurvatureControlParameter(int);
  virtual double getMeshSize() const { return meshAttributes.meshSize; }

  struct {
    mutable GEntity::MeshGenerationStatus status;
    bool refineAllEdges;
    double worst_element_shape, best_element_shape, average_element_shape;
    double smallest_edge_length, longest_edge_length, efficiency_index;
    int nbEdge, nbTriangle;
    int nbGoodQuality, nbGoodLength;
  } meshStatistics;

  // a crude graphical representation using a "cross" represented by points
  // along u and v parameter lines
  std::vector<std::vector<SPoint3> > cross[2];

  // the STL mesh
  std::vector<SPoint2> stl_vertices_uv;
  std::vector<SPoint3> stl_vertices_xyz;
  std::vector<SVector3> stl_normals;
  std::vector<int> stl_triangles;

  // a vertex array containing a geometrical representation of the
  // surface
  VertexArray *va_geom_triangles;

  // a array for accessing the transfinite vertices using a pair of
  // indices
  std::vector<std::vector<MVertex *> > transfinite_vertices;

  // FIXME: testing for constrained packing of parallelogram algo
  std::set<MVertex *> constr_vertices;

  // relocate mesh vertices using parametric coordinates
  void relocateMeshVertices();

  std::vector<MTriangle *> triangles;
  std::vector<MQuadrangle *> quadrangles;
  std::vector<MPolygon *> polygons;

  void addTriangle(MTriangle *t) { triangles.push_back(t); }
  void addQuadrangle(MQuadrangle *q) { quadrangles.push_back(q); }
  void addPolygon(MPolygon *p) { polygons.push_back(p); }
  void addElement(int type, MElement *e);
  void removeElement(int type, MElement *e);

  // get the boundary layer columns
  BoundaryLayerColumns *getColumns() { return &_columns; }

  std::vector<SPoint3> storage1; // sizes and directions storage
  std::vector<SVector3> storage2; // sizes and directions storage
  std::vector<SVector3> storage3; // sizes and directions storage
  std::vector<double> storage4; // sizes and directions storage

  virtual bool reorder(const int elementType, const std::vector<int> &ordering);
};

#endif