Skip to content
Snippets Groups Projects
Select Git revision
  • 759f7cfc73ef61a371fcfbe856d3681298df1576
  • master default protected
  • relaying
  • overlaps_tags_and_distributed_export
  • alphashapes
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

Options.h

Blame
  • MElement.h 21.77 KiB
    // 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 MELEMENT_H
    #define MELEMENT_H
    
    #include <stdio.h>
    #include <algorithm>
    #include <string>
    #include <fstream>
    
    #include "GmshMessage.h"
    #include "ElementType.h"
    #include "MVertex.h"
    #include "MEdge.h"
    #include "MFace.h"
    #include "nodalBasis.h"
    #include "polynomialBasis.h"
    #include "GaussIntegration.h"
    #include "FuncSpaceData.h"
    
    class GModel;
    class JacobianBasis;
    
    // A mesh element.
    class MElement {
    private:
      // the id number of the element (this number is unique and is guaranteed never
      // to change once a mesh has been generated, unless the mesh is explicitly
      // renumbered)
      std::size_t _num;
      // the number of the mesh partition the element belongs to
      short _partition;
      // a visibility flag
      char _visible;
    
    protected:
      // the tolerance used to determine if a point is inside an element, in
      // parametric coordinates
      static double _isInsideTolerance;
    
    protected:
      void _getEdgeRep(MVertex *v0, MVertex *v1, double *x, double *y, double *z,
                       SVector3 *n, int faceIndex = -1);
      void _getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2, double *x, double *y,
                       double *z, SVector3 *n);
    #if defined(HAVE_VISUDEV)
      void _getFaceRepQuad(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
                           double *x, double *y, double *z, SVector3 *n);
    #endif
    
      static bool _getFaceInfo(const MFace &face, const MFace &other, int &sign,
                               int &rot);
    
    public:
      MElement(std::size_t num = 0, int part = 0);
      virtual ~MElement() {}
      // set/get the tolerance for isInside() test
      static void setTolerance(const double tol);
      static double getTolerance();
    
      // return the tag of the element
      virtual std::size_t getNum() const { return _num; }
    
      // force the immutable number (this should never be used, except when
      // explicitly renumbering the mesh)
      void forceNum(std::size_t num);
    
      // return the geometrical dimension of the element
      virtual int getDim() const = 0;
    
      // return the polynomial order the element
      virtual int getPolynomialOrder() const { return 1; }
    
      // return true if the element can be considered as a serendipity element
      virtual bool getIsAssimilatedSerendipity() const
      {
        return ElementType::getSerendipity(getTypeForMSH()) > 0;
      }
      // return true if the element has to be considered as a serendipity element
      virtual bool getIsOnlySerendipity() const
      {
        return ElementType::getSerendipity(getTypeForMSH()) > 1;
      }
    
      // get/set the partition to which the element belongs
      virtual int getPartition() const { return _partition; }
      virtual void setPartition(int num) { _partition = (short)num; }
    
      // get/set the visibility flag
      virtual char getVisibility() const;
      virtual void setVisibility(char val) { _visible = val; }
    
      // get & set the vertices
      virtual std::size_t getNumVertices() const = 0;
      virtual const MVertex *getVertex(int num) const = 0;
      virtual MVertex *getVertex(int num) = 0;
      void getVertices(std::vector<MVertex *> &verts)
      {
        int N = getNumVertices();
        verts.resize(N);
        for(int i = 0; i < N; i++) verts[i] = getVertex(i);
      }
      virtual void setVertex(int num, MVertex *v)
      {
        Msg::Error("Vertex set not supported for this element");
      }
    
      // give an MVertex as input and get its local number
      virtual void getVertexInfo(const MVertex *vertex, int &ithVertex) const
      {
        Msg::Error("Vertex information not available for this element");
      }
      // get the face  using the local orientation defined by Solin
      virtual MFace getFaceSolin(int numFace) = 0;
      // get the global vertex num of a edge using the local orientation defined by
      // Solin
      virtual int getVertexSolin(int numEdge, int numVertex) = 0;
      // get the vertex using the I-deas UNV ordering
      virtual MVertex *getVertexUNV(int num) { return getVertex(num); }
    
      // get the vertex using the VTK ordering
      virtual MVertex *getVertexVTK(int num) { return getVertex(num); }
    
      // get the vertex using the MATLAB ordering
      virtual MVertex *getVertexMATLAB(int num) { return getVertex(num); }
    
      // get the vertex using the Tochnog ordering
      virtual MVertex *getVertexTOCHNOG(int num) { return getVertex(num); }
    
      // get the vertex using the Nastran BDF ordering
      virtual MVertex *getVertexBDF(int num) { return getVertex(num); }
    
      // get the vertex using DIFF ordering
      virtual MVertex *getVertexDIFF(int num) { return getVertex(num); }
    
      // get the vertex using INP ordering
      virtual MVertex *getVertexINP(int num) { return getVertex(num); }
    
      // get the vertex using KEY ordering
      virtual MVertex *getVertexKEY(int num) { return getVertex(num); }
    
      // get the number of vertices associated with edges, faces and volumes (i.e.
      // internal vertices on edges, faces and volume; nonzero only for higher order
      // elements, polygons or polyhedra) - These functions should be renamed to
      // getNumInternal*Vertices(), to avoid confusion with getEdgeVertices() and
      // getFaceVertices()
      virtual int getNumEdgeVertices() const { return 0; }
      virtual int getNumFaceVertices() const { return 0; }
      virtual int getNumVolumeVertices() const { return 0; }
    
      // get the number of primary vertices (first-order element)
      int getNumPrimaryVertices() const
      {
        return getNumVertices() - getNumEdgeVertices() - getNumFaceVertices() -
               getNumVolumeVertices();
      }
    
      // get the edges
      virtual int getNumEdges() const = 0;
      virtual MEdge getEdge(int num) const = 0;
      virtual MEdgeN getHighOrderEdge(int num, int sign);
      virtual MEdgeN getHighOrderEdge(const MEdge &edge)
      {
        int num, sign;
        if(!getEdgeInfo(edge, num, sign)) return MEdgeN();
        return getHighOrderEdge(num, sign);
      }
    
      // give an MEdge as input and get its local number and sign
      virtual bool getEdgeInfo(const MEdge &edge, int &ithEdge, int &sign) const;
      virtual int numEdge2numVertex(int numEdge, int numVert) const
      {
        Msg::Error("Edge information not available for this element");
        return -1;
      }
    
      // get an edge representation for drawing
      virtual int getNumEdgesRep(bool curved) = 0;
      virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z,
                              SVector3 *n) = 0;
    
      // get the faces
      virtual int getNumFaces() = 0;
      virtual MFace getFace(int num) const = 0;
      virtual MFaceN getHighOrderFace(int num, int sign, int rot);
      virtual MFaceN getHighOrderFace(const MFace &face)
      {
        int num, sign, rot;
        if(!getFaceInfo(face, num, sign, rot)) return MFaceN();
        return this->getHighOrderFace(num, sign, rot);
      }
    
      // give an MFace as input and get its local number, sign and rotation
      virtual bool getFaceInfo(const MFace &face, int &ithFace, int &sign,
                               int &rot) const
      {
        Msg::Error("Face information not available for this element");
        return false;
      }
    
      // get a face representation for drawing
      virtual int getNumFacesRep(bool curved) = 0;
      virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z,
                              SVector3 *n) = 0;
    
      // get all the vertices on a edge or a face
      virtual void getEdgeVertices(const int num, std::vector<MVertex *> &v) const
      {
        v.resize(0);
      }
      virtual void getFaceVertices(const int num, std::vector<MVertex *> &v) const
      {
        v.resize(0);
      }
    
      // get and set parent and children for hierarchial grids
      virtual MElement *getParent() const { return NULL; }
      virtual void setParent(MElement *p, bool owner = false) {}
      virtual void updateParent(GModel *gm) {} // used only for reading msh files
      virtual int getNumChildren() const { return 0; }
      virtual MElement *getChild(int i) const { return NULL; }
      virtual bool ownsParent() const { return false; }
    
      // get base element in case of MSubElement
      virtual const MElement *getBaseElement() const { return this; }
      virtual MElement *getBaseElement() { return this; }
    
      // get and set domain for borders
      virtual MElement *getDomain(int i) const { return NULL; }
      virtual void setDomain(MElement *e, int i) {}
    
      // get the type of the element
      virtual int getType() const = 0;
    
      // get the max/min edge length
      virtual double maxEdge();
      virtual double minEdge();
    
      // max. distance between curved and straight element among all high-order
      // nodes
      double maxDistToStraight() const;
    
      // get the quality measures
      double skewness();
      virtual double gammaShapeMeasure() { return 0.; }
      virtual double etaShapeMeasure() { return 0.; }
      double minSICNShapeMeasure()
      {
        double sICNMin, sICNMax;
        signedInvCondNumRange(sICNMin, sICNMax);
        return sICNMin;
      }
      double minSIGEShapeMeasure()
      {
        double minSIGE, maxSIGE;
        signedInvGradErrorRange(minSIGE, maxSIGE);
        return minSIGE;
      }
      double distoShapeMeasure()
      {
        double jmin, jmax;
        scaledJacRange(jmin, jmax);
        return jmin;
      }
      double minIsotropyMeasure(bool knownValid = false, bool reversedOk = false);
      double minScaledJacobian(bool knownValid = false, bool reversedOk = false);
      double specialQuality();
      double specialQuality2();
      virtual double angleShapeMeasure() { return 1.0; }
      virtual void scaledJacRange(double &jmin, double &jmax,
                                  GEntity *ge = 0) const;
      virtual void idealJacRange(double &jmin, double &jmax, GEntity *ge = 0);
      virtual void signedInvCondNumRange(double &iCNMin, double &iCNMax,
                                         GEntity *ge = 0);
      virtual void signedInvGradErrorRange(double &minSIGE, double &maxSIGE);
    
      // get the radius of the inscribed circle/sphere if it exists, otherwise get
      // the minimum radius of all the circles/spheres tangent to the most
      // boundaries of the element.
      virtual double getInnerRadius() { return 0.; }
    
      // get the radius of the circumscribed circle/sphere if it exists, otherwise
      // get the maximum radius of all the circles/spheres tangent to the most
      // boundaries of the element.
      virtual double getOuterRadius() { return 0.; }
    
      // compute the barycenter
      virtual SPoint3 barycenter(bool primary = false) const;
      // compute the barycenter without divided by the number of nodes
      virtual SPoint3 fastBarycenter(bool primary = false) const;
      virtual SPoint3 barycenterUVW() const;
      // compute the barycenter in infinity norm
      virtual SPoint3 barycenter_infty() const;
    
      // reverse the orientation of the element
      virtual void reverse() {}
    
      // get volume of element
      virtual double getVolume();
    
      // return sign of volume (+1 or -1) for 3D elements (or 0 if element has zero
      // volume)
      virtual int getVolumeSign();
    
      // compute and change the orientation of 3D elements to get positive volume
      // (return false if element has zero volume)
      virtual bool setVolumePositive();
    
      // compute the extrema of the Jacobian determinant return 1 if the element is
      // valid, 0 if the element is invalid, -1 if the element is reversed
      int getValidity();
    
      // return an information string for the element
      virtual std::string getInfoString(bool multline);
    
      // get the function space for the element
      virtual const nodalBasis *getFunctionSpace(int order = -1,
                                                 bool serendip = false) const;
      virtual const FuncSpaceData getFuncSpaceData(int order = -1,
                                                   bool serendip = false) const;
    
      // get the function space for the jacobian of the element
      virtual const JacobianBasis *
      getJacobianFuncSpace(int orderElement = -1) const;
      virtual const FuncSpaceData
      getJacobianFuncSpaceData(int orderElement = -1) const;
    
      // return parametric coordinates (u,v,w) of a vertex
      virtual void getNode(int num, double &u, double &v, double &w) const;
    
      // return the interpolating nodal shape functions evaluated at point (u,v,w)
      // in parametric coordinates (if order == -1, use the polynomial order of the
      // element)
      virtual void getShapeFunctions(double u, double v, double w, double s[],
                                     int order = -1) const;
    
      // return the gradient of the nodal shape functions evaluated at point (u,v,w)
      // in parametric coordinates (if order == -1, use the polynomial order of the
      // element)
      virtual void getGradShapeFunctions(double u, double v, double w,
                                         double s[][3], int order = -1) const;
      virtual void getHessShapeFunctions(double u, double v, double w,
                                         double s[][3][3], int order = -1) const;
      virtual void getThirdDerivativeShapeFunctions(double u, double v, double w,
                                                    double s[][3][3][3],
                                                    int order = -1) const;
      // return the Jacobian of the element evaluated at point (u,v,w) in parametric
      // coordinates: jac[i][j] = \partial x_j / \partial u_i (beware the
      // transposition compared to the usual definition of the Jacobian)
      virtual double getJacobian(const fullMatrix<double> &gsf,
                                 double jac[3][3]) const;
      // To be compatible with _vgrads of functionSpace without having to put under
      // fullMatrix form
      virtual double getJacobian(const std::vector<SVector3> &gsf,
                                 double jac[3][3]) const;
      // jac is an row-major order array
      virtual double getJacobian(const std::vector<SVector3> &gsf,
                                 double *jac) const;
      virtual double getJacobian(double u, double v, double w,
                                 double jac[3][3]) const;
      double getJacobian(double u, double v, double w, fullMatrix<double> &j) const
      {
        double JAC[3][3];
        const double detJ = getJacobian(u, v, w, JAC);
        for(int i = 0; i < 3; i++) {
          j(i, 0) = JAC[i][0];
          j(i, 1) = JAC[i][1];
          j(i, 2) = JAC[i][2];
        }
        return detJ;
      }
      virtual double getPrimaryJacobian(double u, double v, double w,
                                        double jac[3][3]) const;
      double getJacobianDeterminant(double u, double v, double w) const
      {
        double jac[3][3];
        return getJacobian(u, v, w, jac);
      }
      void getSignedJacobian(fullVector<double> &jacobian, int o = -1) const;
      void getNodesCoord(fullMatrix<double> &nodesXYZ) const;
      virtual int getNumShapeFunctions() const { return getNumVertices(); }
      virtual int getNumPrimaryShapeFunctions() const
      {
        return getNumPrimaryVertices();
      }
      virtual const MVertex *getShapeFunctionNode(int i) const
      {
        return getVertex(i);
      }
      virtual MVertex *getShapeFunctionNode(int i) { return getVertex(i); }
    
      // return the eigenvalues of the metric evaluated at point (u,v,w) in
      // parametric coordinates
      virtual double getEigenvaluesMetric(double u, double v, double w,
                                          double values[3]) const;
    
      // get the point in cartesian coordinates corresponding to the point (u,v,w)
      // in parametric coordinates
      virtual void pnt(double u, double v, double w, SPoint3 &p) const;
      virtual void pnt(double u, double v, double w, double *p) const;
      // To be compatible with functionSpace without changing form
      virtual void pnt(const std::vector<double> &sf, SPoint3 &p) const;
      virtual void primaryPnt(double u, double v, double w, SPoint3 &p);
    
      // invert the parametrisation
      virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
    
      // move point between parent and element parametric spaces
      virtual void movePointFromParentSpaceToElementSpace(double &u, double &v,
                                                          double &w) const;
      virtual void movePointFromElementSpaceToParentSpace(double &u, double &v,
                                                          double &w) const;
    
      // test if a point, given in parametric coordinates, belongs to the element
      virtual bool isInside(double u, double v, double w) const = 0;
    
      // interpolate the given nodal data (resp. its gradient, curl and divergence)
      // at point (u,v,w) in parametric coordinates
      double interpolate(double val[], double u, double v, double w, int stride = 1,
                         int order = -1);
      void interpolateGrad(double val[], double u, double v, double w, double f[],
                           int stride = 1, double invjac[3][3] = 0, int order = -1);
      void interpolateCurl(double val[], double u, double v, double w, double f[],
                           int stride = 3, int order = -1);
      double interpolateDiv(double val[], double u, double v, double w,
                            int stride = 3, int order = -1);
    
      // integration routines
      virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
      {
        Msg::Error("No integration points defined for this type of element: %d",
                   this->getType());
        *npts = 0;
        *pts = 0;
      }
      double integrate(double val[], int pOrder, int stride = 1, int order = -1);
      // val[] must contain interpolation data for face/edge vertices of given
      // edge/face
      double integrateCirc(double val[], int edge, int pOrder, int order = -1);
      double integrateFlux(double val[], int face, int pOrder, int order = -1);
    
      // IO routines
      virtual void writeMSH2(FILE *fp, double version = 1.0, bool binary = false,
                             int num = 0, int elementary = 1, int physical = 1,
                             int parentNum = 0, int dom1Num = 0, int dom2Num = 0,
                             std::vector<short> *ghosts = 0);
      virtual void writeMSH3(FILE *fp, bool binary = false, int elementary = 1,
                             std::vector<short> *ghosts = 0);
      virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber,
                            bool printSICN, bool printSIGE, bool printGamma,
                            bool printDisto, double scalingFactor = 1.0,
                            int elementary = 1);
      virtual void writeSTL(FILE *fp, bool binary = false,
                            double scalingFactor = 1.0);
      virtual void writeVRML(FILE *fp);
      virtual void writePLY2(FILE *fp);
      virtual void writeUNV(FILE *fp, int num = 0, int elementary = 1,
                            int physical = 1);
      virtual void writeVTK(FILE *fp, bool binary = false, bool bigEndian = false);
      virtual void writeMATLAB(FILE *fp, int filetype, int elementary = 0,
                               int physical = 0, bool binary = false);
      virtual void writeTOCHNOG(FILE *fp, int num);
      virtual void writeMESH(FILE *fp, int elementTagType = 1, int elementary = 1,
                             int physical = 0);
      virtual void writeNEU(FILE *fp, unsigned gambitType, int adjust,
                            int phys = 0);
      virtual void writeIR3(FILE *fp, int elementTagType, int num, int elementary,
                            int physical);
      virtual void writeBDF(FILE *fp, int format = 0, int elementTagType = 1,
                            int elementary = 1, int physical = 0);
      virtual void writeDIFF(FILE *fp, int num, bool binary = false,
                             int physical_property = 1);
      virtual void writeINP(FILE *fp, int num);
      virtual void writeKEY(FILE *fp, int pid, int num);
      virtual void writeSU2(FILE *fp, int num);
    
      // info for specific IO formats (returning 0 means that the element is not
      // implemented in that format)
      virtual int getTypeForMSH() const { return 0; }
      virtual int getTypeForUNV() const { return 0; }
      virtual int getTypeForVTK() const { return 0; }
      virtual const char *getStringForTOCHNOG() const { return 0; }
      virtual const char *getStringForPOS() const { return 0; }
      virtual const char *getStringForBDF() const { return 0; }
      virtual const char *getStringForDIFF() const { return 0; }
      virtual const char *getStringForINP() const { return 0; }
      virtual const char *getStringForKEY() const { return 0; }
    
      // return the number of vertices, as well as the element name if 'name' != 0
      static unsigned int getInfoMSH(const int typeMSH,
                                     const char **const name = 0);
      virtual int getNumVerticesForMSH() { return getNumVertices(); }
      virtual void getVerticesIdForMSH(std::vector<int> &verts);
    
      // copy element and parent if any, vertexMap contains the new vertices
      virtual MElement *copy(std::map<int, MVertex *> &vertexMap,
                             std::map<MElement *, MElement *> &newParents,
                             std::map<MElement *, MElement *> &newDomains);
    
      // Return the number of nodes that this element must have with the other in
      // order to put an edge between them in the dual graph used during the
      // partitioning.
      virtual int numCommonNodesInDualGraph(const MElement *const other) const = 0;
    };
    
    class MElementFactory {
    public:
      MElement *create(int type, std::vector<MVertex *> &v, std::size_t num = 0,
                       int part = 0, bool owner = false, int parent = 0,
                       MElement *parent_ptr = NULL, MElement *d1 = 0,
                       MElement *d2 = 0);
      MElement *create(int num, int type, const std::vector<int> &data,
                       GModel *model);
    };
    
    // Traits of various elements based on the dimension.  These generally define
    // the faces of 2-D elements as MEdge and 3-D elements as MFace.
    
    template <unsigned DIM> struct DimTr;
    template <> struct DimTr<2> {
      typedef MEdge FaceT;
      static int getNumFace(MElement *const element)
      {
        return element->getNumEdges();
      }
      static MEdge getFace(MElement *const element, const int iFace)
      {
        return element->getEdge(iFace);
      }
      static void getAllFaceVertices(MElement *const element, const int iFace,
                                     std::vector<MVertex *> &v)
      {
        element->getEdgeVertices(iFace, v);
      }
    };
    
    template <> struct DimTr<3> {
      typedef MFace FaceT;
      static int getNumFace(MElement *const element)
      {
        return element->getNumFaces();
      }
      static MFace getFace(MElement *const element, const int iFace)
      {
        return element->getFace(iFace);
      }
      static void getAllFaceVertices(MElement *const element, const int iFace,
                                     std::vector<MVertex *> &v)
      {
        element->getFaceVertices(iFace, v);
      }
    };
    
    struct Less_ElementPtr {
      bool operator()(const MElement *e1, const MElement *e2) const
      {
        return e1->getNum() < e2->getNum();
      }
    };
    
    #endif