mesh.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#ifndef MFEM_MESH
#define MFEM_MESH

#include "../config/config.hpp"
#include "../general/stable3d.hpp"
#include "triangle.hpp"
#include "tetrahedron.hpp"
#include "vertex.hpp"
#include "ncmesh.hpp"
#include "../fem/eltrans.hpp"
#include "../fem/coefficient.hpp"
#include <iostream>

namespace mfem
{

// Data type mesh

class KnotVector;
class NURBSExtension;
class FiniteElementSpace;
class GridFunction;
struct Refinement;

#ifdef MFEM_USE_MPI
class ParMesh;
class ParNCMesh;
#endif

class Mesh
{
#ifdef MFEM_USE_MPI
   friend class ParMesh;
   friend class ParNCMesh;
#endif
   friend class NURBSExtension;

protected:
   int Dim;
   int spaceDim;

   int NumOfVertices, NumOfElements, NumOfBdrElements;
   int NumOfEdges, NumOfFaces;

   int State, WantTwoLevelState;

   // 0 = Empty, 1 = Standard (NetGen), 2 = TrueGrid
   int meshgen;

   int c_NumOfVertices, c_NumOfElements, c_NumOfBdrElements;
   int f_NumOfVertices, f_NumOfElements, f_NumOfBdrElements;
   int c_NumOfEdges, c_NumOfFaces;
   int f_NumOfEdges, f_NumOfFaces;

   Array<Element *> elements;
   // Vertices are only at the corners of elements, where you would expect them
   // in the lowest-order mesh.
   Array<Vertex> vertices;
   Array<Element *> boundary;
   Array<Element *> faces;

   struct FaceInfo
   {
      int Elem1No, Elem2No, Elem1Inf, Elem2Inf;
      int NCFace; /* -1 if this is a regular conforming/boundary face;
                     index into 'nc_faces_info' if >= 0. */
   };
   // NOTE: in NC meshes, master faces have Elem2No == -1. Slave faces on the
   // other hand have Elem2No and Elem2Inf set to the master face's element and
   // its local face number.

   struct NCFaceInfo
   {
      bool Slave; // true if this is a slave face, false if master face
      int MasterFace; // if Slave, this is the index of the master face
      const DenseMatrix* PointMatrix; // if Slave, position within master face
      // (NOTE: PointMatrix points to a matrix owned by NCMesh.)

      NCFaceInfo(bool slave, int master, const DenseMatrix* pm)
         : Slave(slave), MasterFace(master), PointMatrix(pm) {}
   };

   Array<FaceInfo> faces_info;
   Array<NCFaceInfo> nc_faces_info;

   Table *el_to_edge;
   Table *el_to_face;
   Table *el_to_el;
   Array<int> be_to_edge;  // for 2D
   Table *bel_to_edge;     // for 3D
   Array<int> be_to_face;
   mutable Table *face_edge;
   mutable Table *edge_vertex;

   Table *c_el_to_edge, *f_el_to_edge, *c_bel_to_edge, *f_bel_to_edge;
   Array<int> fc_be_to_edge; // swapped with be_to_edge when switching state
   Table *c_el_to_face, *f_el_to_face; // for 3D two-level state
   Array<FaceInfo> fc_faces_info;      // for 3D two-level state

   IsoparametricTransformation Transformation, Transformation2;
   IsoparametricTransformation FaceTransformation, EdgeTransformation;
   FaceElementTransformations FaceElemTr;

   // Nodes are only active for higher order meshes, and share locations with
   // the vertices, plus all the higher- order control points within the
   // element and along the edges and on the faces.
   GridFunction *Nodes;
   int own_nodes;

   // Backup of the coarse mesh. Only used if WantTwoLevelState == 1 and
   // nonconforming refinements are used.
   Mesh* nc_coarse_level;

   static const int tet_faces[4][3];
   static const int hex_faces[6][4]; // same as Hexahedron::faces

   static const int tri_orientations[6][3];
   static const int quad_orientations[8][4];

#ifdef MFEM_USE_MEMALLOC
   friend class Tetrahedron;

   MemAlloc <Tetrahedron, 1024> TetMemory;
   MemAlloc <BisectedElement, 1024> BEMemory;
#endif

public:
   Array<int> attributes;
   Array<int> bdr_attributes;

   NURBSExtension *NURBSext;
   NCMesh *ncmesh;

protected:
   void Init();

   void InitTables();

   void DeleteTables();

   void DeleteCoarseTables();

   Element *ReadElementWithoutAttr(std::istream &);
   static void PrintElementWithoutAttr(const Element *, std::ostream &);

   Element *ReadElement(std::istream &);
   static void PrintElement(const Element *, std::ostream &);

   void SetMeshGen(); // set 'meshgen'

   double GetLength(int i, int j) const;

   void GetElementJacobian(int i, DenseMatrix &J);

   void MarkForRefinement();
   void MarkTriMeshForRefinement();
   void GetEdgeOrdering(DSTable &v_to_v, Array<int> &order);
   void MarkTetMeshForRefinement();

   void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert);
   void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert);

   STable3D *GetFacesTable();
   STable3D *GetElementToFaceTable(int ret_ftbl = 0);

   void RedRefinement(int i, const DSTable &v_to_v,
                      int *edge1, int *edge2, int *middle)
   { UniformRefinement(i, v_to_v, edge1, edge2, middle); }

   void GreenRefinement(int i, const DSTable &v_to_v,
                        int *edge1, int *edge2, int *middle)
   { Bisection(i, v_to_v, edge1, edge2, middle); }

   void Bisection(int i, const DSTable &, int *, int *, int *);

   void Bisection(int i, const DSTable &, int *);

   void UniformRefinement(int i, const DSTable &, int *, int *, int *);

   void AverageVertices (int * indexes, int n, int result);

   void UpdateNodes();

   virtual void QuadUniformRefinement();

   virtual void HexUniformRefinement();

   virtual void NURBSUniformRefinement();

   virtual void LocalRefinement(const Array<int> &marked_el, int type = 3);

   virtual void NonconformingRefinement(const Array<Refinement> &refinements,
                                        int nc_limit = 0);

   void LoadPatchTopo(std::istream &input, Array<int> &edge_to_knot);

   void UpdateNURBS();

   void PrintTopo(std::ostream &out, const Array<int> &e_to_k) const;

   void BisectTriTrans (DenseMatrix &pointmat, Triangle *tri,
                        int child);

   void BisectTetTrans (DenseMatrix &pointmat, Tetrahedron *tet,
                        int child);

   int GetFineElemPath (int i, int j);

   int GetBisectionHierarchy (Element *E);

   void GetLocalPtToSegTransformation(IsoparametricTransformation &, int);
   void GetLocalSegToTriTransformation (IsoparametricTransformation &loc,
                                        int i);
   void GetLocalSegToQuadTransformation (IsoparametricTransformation &loc,
                                         int i);
   void GetLocalTriToTetTransformation (IsoparametricTransformation &loc,
                                        int i);
   void GetLocalQuadToHexTransformation (IsoparametricTransformation &loc,
                                         int i);
   void ApplySlaveTransformation(IsoparametricTransformation &transf,
                                 const FaceInfo &fi);
   bool IsSlaveFace(const FaceInfo &fi);

   static int GetTriOrientation (const int * base, const int * test);
   static int GetQuadOrientation (const int * base, const int * test);

   static void GetElementArrayEdgeTable(const Array<Element*> &elem_array,
                                        const DSTable &v_to_v,
                                        Table &el_to_edge);

   void GetVertexToVertexTable(DSTable &) const;

   int GetElementToEdgeTable(Table &, Array<int> &);

   void AddPointFaceElement(int lf, int gf, int el);

   void AddSegmentFaceElement (int lf, int gf, int el, int v0, int v1);

   void AddTriangleFaceElement (int lf, int gf, int el,
                                int v0, int v1, int v2);

   void AddQuadFaceElement (int lf, int gf, int el,
                            int v0, int v1, int v2, int v3);
   bool FaceIsTrueInterior(int FaceNo) const
   {
      return FaceIsInterior(FaceNo) || (faces_info[FaceNo].Elem2Inf >= 0);
   }

   // shift cyclically 3 integers left-to-right
   inline static void ShiftL2R(int &, int &, int &);
   // shift cyclically 3 integers so that the smallest is first
   inline static void Rotate3(int &, int &, int &);

   void FreeElement (Element *E);

   void GenerateFaces();
   void GenerateNCFaceInfo();

   void InitMesh(int _Dim, int _spaceDim, int NVert, int NElem, int NBdrElem);

   void Make3D(int nx, int ny, int nz, Element::Type type, int generate_edges,
               double sx, double sy, double sz);

   void Make2D(int nx, int ny, Element::Type type, int generate_edges,
               double sx, double sy);

   void Make1D(int n, double sx = 1.0);

   void InitFromNCMesh(const NCMesh &ncmesh);

   Mesh(const NCMesh &ncmesh);

   void Swap(Mesh& other, bool non_geometry = false);

public:

   enum { NORMAL, TWO_LEVEL_COARSE, TWO_LEVEL_FINE };

   Mesh() { Init(); InitTables(); meshgen = 0; Dim = 0; }

   explicit Mesh(const Mesh &mesh, bool copy_nodes = true);

   Mesh(int _Dim, int NVert, int NElem, int NBdrElem = 0, int _spaceDim= -1)
   {
      if (_spaceDim == -1)
      {
         _spaceDim = _Dim;
      }
      InitMesh(_Dim, _spaceDim, NVert, NElem, NBdrElem);
   }

   Element *NewElement(int geom);

   void AddVertex(const double *);
   void AddTri(const int *vi, int attr = 1);
   void AddTriangle(const int *vi, int attr = 1);
   void AddQuad(const int *vi, int attr = 1);
   void AddTet(const int *vi, int attr = 1);
   void AddHex(const int *vi, int attr = 1);
   void AddHexAsTets(const int *vi, int attr = 1);
   // 'elem' should be allocated using the NewElement method
   void AddElement(Element *elem)     { elements[NumOfElements++] = elem; }
   void AddBdrElement(Element *elem)  { boundary[NumOfBdrElements++] = elem; }
   void AddBdrSegment(const int *vi, int attr = 1);
   void AddBdrTriangle(const int *vi, int attr = 1);
   void AddBdrQuad(const int *vi, int attr = 1);
   void AddBdrQuadAsTriangles(const int *vi, int attr = 1);
   void GenerateBoundaryElements();
   void FinalizeTriMesh(int generate_edges = 0, int refine = 0,
                        bool fix_orientation = true);
   void FinalizeQuadMesh(int generate_edges = 0, int refine = 0,
                         bool fix_orientation = true);
   void FinalizeTetMesh(int generate_edges = 0, int refine = 0,
                        bool fix_orientation = true);
   void FinalizeHexMesh(int generate_edges = 0, int refine = 0,
                        bool fix_orientation = true);

   void SetAttributes();

#ifdef MFEM_USE_GECKO

   void GetGeckoElementReordering(Array<int> &ordering);
#endif

   void ReorderElements(const Array<int> &ordering, bool reorder_vertices = true);

   Mesh(int nx, int ny, int nz, Element::Type type, int generate_edges = 0,
        double sx = 1.0, double sy = 1.0, double sz = 1.0)
   {
      Make3D(nx, ny, nz, type, generate_edges, sx, sy, sz);
   }

   Mesh(int nx, int ny, Element::Type type, int generate_edges = 0,
        double sx = 1.0, double sy = 1.0)
   {
      Make2D(nx, ny, type, generate_edges, sx, sy);
   }

   explicit Mesh(int n, double sx = 1.0)
   {
      Make1D(n, sx);
   }

   Mesh(std::istream &input, int generate_edges = 0, int refine = 1,
        bool fix_orientation = true);

   Mesh(Mesh *mesh_array[], int num_pieces);

   /* This is similar to the above mesh constructor, but here the current
      mesh is destroyed and another one created based on the data stream
      again given in MFEM, netgen, or VTK format. If generate_edges = 0
      (default) edges are not generated, if 1 edges are generated. */
   void Load(std::istream &input, int generate_edges = 0, int refine = 1,
             bool fix_orientation = true);

   inline int MeshGenerator() { return meshgen; }

   inline int GetNV() const { return NumOfVertices; }

   inline int GetNE() const { return NumOfElements; }

   inline int GetNBE() const { return NumOfBdrElements; }

   inline int GetNEdges() const { return NumOfEdges; }

   inline int GetNFaces() const { return NumOfFaces; }

   int GetNumFaces() const;

   inline int EulerNumber() const
   { return NumOfVertices - NumOfEdges + NumOfFaces - NumOfElements; }
   inline int EulerNumber2D() const
   { return NumOfVertices - NumOfEdges + NumOfElements; }

   int Dimension() const { return Dim; }
   int SpaceDimension() const { return spaceDim; }

   const double *GetVertex(int i) const { return vertices[i](); }
   double *GetVertex(int i) { return vertices[i](); }

   const Element *GetElement(int i) const { return elements[i]; }

   Element *GetElement(int i) { return elements[i]; }

   const Element *GetBdrElement(int i) const { return boundary[i]; }

   Element *GetBdrElement(int i) { return boundary[i]; }

   const Element *GetFace(int i) const { return faces[i]; }

   int GetFaceBaseGeometry(int i) const;

   int GetElementBaseGeometry(int i) const
   { return elements[i]->GetGeometryType(); }

   int GetBdrElementBaseGeometry(int i) const
   { return boundary[i]->GetGeometryType(); }

   void GetElementVertices(int i, Array<int> &dofs) const
   { elements[i]->GetVertices(dofs); }

   void GetBdrElementVertices(int i, Array<int> &dofs) const
   { boundary[i]->GetVertices(dofs); }

   void GetElementEdges(int i, Array<int> &edges, Array<int> &cor) const;

   void GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const;

   void GetFaceEdges(int i, Array<int> &, Array<int> &) const;

   void GetFaceVertices(int i, Array<int> &vert) const
   {
      if (Dim == 1)
      {
         vert.SetSize(1); vert[0] = i;
      }
      else
      {
         faces[i]->GetVertices(vert);
      }
   }

   void GetEdgeVertices(int i, Array<int> &vert) const;

   Table *GetFaceEdgeTable() const;

   Table *GetEdgeVertexTable() const;

   void GetElementFaces(int i, Array<int> &, Array<int> &) const;

   void GetBdrElementFace(int i, int *, int *) const;

   int GetBdrElementEdgeIndex(int i) const;

   int GetElementType(int i) const;

   int GetBdrElementType(int i) const;

   /* Return point matrix of element i of dimension Dim X #dofs, where for
      every degree of freedom we give its coordinates in space of dimension
      Dim. */
   void GetPointMatrix(int i, DenseMatrix &pointmat) const;

   /* Return point matrix of boundary element i of dimension Dim X #dofs,
      where for every degree of freedom we give its coordinates in space
      of dimension Dim. */
   void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const;

   static FiniteElement *GetTransformationFEforElementType(int);

   void GetElementTransformation(int i, IsoparametricTransformation *ElTr);

   ElementTransformation *GetElementTransformation(int i);

   void GetElementTransformation(int i, const Vector &nodes,
                                 IsoparametricTransformation *ElTr);

   ElementTransformation * GetBdrElementTransformation(int i);
   void GetBdrElementTransformation(int i, IsoparametricTransformation *ElTr);

   void GetFaceTransformation(int i, IsoparametricTransformation *FTr);

   ElementTransformation *GetFaceTransformation(int FaceNo);

   void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr);

   ElementTransformation *GetEdgeTransformation(int EdgeNo);

   FaceElementTransformations *GetFaceElementTransformations(int FaceNo,
                                                             int mask = 31);

   FaceElementTransformations *GetInteriorFaceTransformations (int FaceNo)
   {
      if (faces_info[FaceNo].Elem2No < 0) { return NULL; }
      return GetFaceElementTransformations (FaceNo);
   }

   FaceElementTransformations *GetBdrFaceTransformations (int BdrElemNo);

   bool FaceIsInterior(int FaceNo) const
   {
      return (faces_info[FaceNo].Elem2No >= 0);
   }
   void GetFaceElements (int Face, int *Elem1, int *Elem2);
   void GetFaceInfos (int Face, int *Inf1, int *Inf2);

   void CheckElementOrientation(bool fix_it = true);
   void CheckBdrElementOrientation(bool fix_it = true);

   int GetAttribute(int i) const { return elements[i]->GetAttribute();}

   int GetBdrAttribute(int i) const { return boundary[i]->GetAttribute(); }

   const Table &ElementToElementTable();

   const Table &ElementToFaceTable() const;

   const Table &ElementToEdgeTable() const;

   Table *GetVertexToElementTable();

   Table *GetFaceToElementTable() const;

   virtual void ReorientTetMesh();

   int *CartesianPartitioning(int nxyz[]);
   int *GeneratePartitioning(int nparts, int part_method = 1);
   void CheckPartitioning(int *partitioning);

   void CheckDisplacements(const Vector &displacements, double &tmax);

   // Vertices are only at the corners of elements, where you would expect them
   // in the lowest-order mesh.
   void MoveVertices(const Vector &displacements);
   void GetVertices(Vector &vert_coord) const;
   void SetVertices(const Vector &vert_coord);

   // Nodes are only active for higher order meshes, and share locations with
   // the vertices, plus all the higher- order control points within the element
   // and along the edges and on the faces.
   void GetNode(int i, double *coord);
   void SetNode(int i, const double *coord);

   // Node operations for curved mesh.
   // They call the corresponding '...Vertices' method if the
   // mesh is not curved (i.e. Nodes == NULL).
   void MoveNodes(const Vector &displacements);
   void GetNodes(Vector &node_coord) const;
   void SetNodes(const Vector &node_coord);

   GridFunction *GetNodes() { return Nodes; }
   void NewNodes(GridFunction &nodes, bool make_owner = false);
   void SwapNodes(GridFunction *&nodes, int &own_nodes_);

   void GetNodes(GridFunction &nodes) const;
   void SetNodalFESpace(FiniteElementSpace *nfes);
   void SetNodalGridFunction(GridFunction *nodes, bool make_owner = false);
   const FiniteElementSpace *GetNodalFESpace();

   void SetCurvature(int order, bool discont = false, int space_dim = -1,
                     int ordering = 1);

   void UniformRefinement();

   void GeneralRefinement(const Array<Refinement> &refinements,
                          int nonconforming = -1, int nc_limit = 0);

   void GeneralRefinement(const Array<int> &el_to_refine,
                          int nonconforming = -1, int nc_limit = 0);

   void EnsureNCMesh();

   void RandomRefinement(int levels, int frac = 2, bool aniso = false,
                         int nonconforming = -1, int nc_limit = -1,
                         int seed = 0 /* should be the same on all CPUs */);

   void RefineAtVertex(const Vertex& vert, int levels,
                       double eps = 0.0, int nonconforming = -1);

   // NURBS mesh refinement methods
   void KnotInsert(Array<KnotVector *> &kv);
   void DegreeElevate(int t);

   void UseTwoLevelState (int use)
   {
      if (!use && State != Mesh::NORMAL)
      {
         SetState (Mesh::NORMAL);
      }
      WantTwoLevelState = use;
   }

   void SetState (int s);

   int GetState() const { return State; }

   int GetNumFineElems (int i);

   int GetRefinementType (int i);

   int GetFineElem (int i, int j);

   ElementTransformation * GetFineElemTrans (int i, int j);

   virtual void PrintXG(std::ostream &out = std::cout) const;

   virtual void Print(std::ostream &out = std::cout) const;

   void PrintVTK(std::ostream &out);

   void PrintVTK(std::ostream &out, int ref, int field_data=0);

   void GetElementColoring(Array<int> &colors, int el0 = 0);

   void PrintWithPartitioning (int *partitioning,
                               std::ostream &out, int elem_attr = 0) const;

   void PrintElementsWithPartitioning (int *partitioning,
                                       std::ostream &out,
                                       int interior_faces = 0);


   void PrintSurfaces(const Table &Aface_face, std::ostream &out) const;

   void ScaleSubdomains (double sf);
   void ScaleElements (double sf);

   void Transform(void (*f)(const Vector&, Vector&));
   void Transform(VectorCoefficient &deformation);

   void RemoveUnusedVertices();

   void RemoveInternalBoundaries();

   double GetElementSize(int i, int type = 0);

   double GetElementSize(int i, const Vector &dir);

   double GetElementVolume(int i);

   void PrintCharacteristics(Vector *Vh = NULL, Vector *Vk = NULL,
                             std::ostream &out = std::cout);

   virtual void PrintInfo(std::ostream &out = std::cout)
   {
      PrintCharacteristics(NULL, NULL, out);
   }

   void MesquiteSmooth(const int mesquite_option = 0);

   virtual ~Mesh();
};

std::ostream &operator<<(std::ostream &out, const Mesh &mesh);


class NodeExtrudeCoefficient : public VectorCoefficient
{
private:
   int n, layer;
   double p[2], s;
   Vector tip;
public:
   NodeExtrudeCoefficient(const int dim, const int _n, const double _s);
   void SetLayer(const int l) { layer = l; }
   using VectorCoefficient::Eval;
   virtual void Eval(Vector &V, ElementTransformation &T,
                     const IntegrationPoint &ip);
   virtual ~NodeExtrudeCoefficient() { }
};


Mesh *Extrude1D(Mesh *mesh, const int ny, const double sy,
                const bool closed = false);


// inline functions
inline void Mesh::ShiftL2R(int &a, int &b, int &c)
{
   int t = a;
   a = c;  c = b;  b = t;
}

inline void Mesh::Rotate3(int &a, int &b, int &c)
{
   if (a < b)
   {
      if (a > c)
      {
         ShiftL2R(a, b, c);
      }
   }
   else
   {
      if (b < c)
      {
         ShiftL2R(c, b, a);
      }
      else
      {
         ShiftL2R(a, b, c);
      }
   }
}

}

#endif