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.googlecode.com.
//
// 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

// Data type mesh

class NURBSExtension;
class FiniteElementSpace;
class GridFunction;

#ifdef MFEM_USE_MPI
class ParMesh;
#endif

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

protected:
   int Dim;

   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;
   Array<Vertex> vertices;
   Array<Element *> boundary;
   Array<Element *> faces;

   class FaceInfo { public: int Elem1No, Elem2No, Elem1Inf, Elem2Inf; };
   Array<FaceInfo> 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;
   FaceElementTransformations FaceElemTr;

   GridFunction *Nodes;
   int own_nodes;

#ifdef MFEM_USE_MEMALLOC
   friend class Tetrahedron;

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

   Element *NewElement(int geom);

   void Init();

   void InitTables();

   void DeleteTables();

   void DeleteCoarseTables();

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

   double GetLength(int i, int j) const;

   void GetElementJacobian(int i, DenseMatrix &J);

   void MarkForRefinement();
   void MarkTriMeshForRefinement();
   void MarkTetMeshForRefinement();

   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();

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

   void UpdateNURBS();

   void PrintTopo(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);

   FiniteElement *GetTransformationFEforElementType (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);

   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 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);

   // 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();

public:

   enum { NORMAL, TWO_LEVEL_COARSE, TWO_LEVEL_FINE };

   Array<int> attributes;
   Array<int> bdr_attributes;

   NURBSExtension *NURBSext;

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

   Mesh (int _Dim, int NVert, int NElem, int NBdrElem = 0);
   void AddVertex (double *);
   void AddElement (Element *elem)  { elements[NumOfElements++] = elem; }
   void AddTri (int *vi, int attr = 1);
   void AddTriangle (int *vi, int attr = 1);
   void AddQuad (int *vi, int attr = 1);
   void AddTet (int *vi, int attr = 1);
   void AddHex (int *vi, int attr = 1);
   void AddBdrSegment (int *vi, int attr = 1);
   void AddBdrTriangle (int *vi, int attr = 1);
   void AddBdrQuad (int *vi, int attr = 1);
   void GenerateBoundaryElements();
   void FinalizeTriMesh (int generate_edges = 0, int refine = 0);
   void FinalizeQuadMesh (int generate_edges = 0, int refine = 0);
   void FinalizeTetMesh (int generate_edges = 0, int refine = 0);
   void FinalizeHexMesh (int generate_edges = 0, int refine = 0);

   void SetAttributes();

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

   explicit Mesh(int n);

   Mesh ( istream &input, int generate_edges = 0, int refine = 1);

   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 netgen format. If generate_edges = 0 (default) edges
      are not generated, if 1 edges are generated. */
   void Load ( istream &input, int generate_edges = 0, int refine = 1);

   void SetNodalFESpace(FiniteElementSpace *nfes);
   void SetNodalGridFunction(GridFunction *nodes);
   const FiniteElementSpace *GetNodalFESpace();

   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; }

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

   int Dimension() const { return Dim; }

   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> &, Array<int> &) const;

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

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

   void GetFaceVertices(int i, Array<int> &vert) const
   { 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;

   void GetElementTransformation(int i, IsoparametricTransformation *ElTr);

   ElementTransformation *GetElementTransformation(int i);

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

   ElementTransformation * GetBdrElementTransformation(int i);

   void GetFaceTransformation(int i, IsoparametricTransformation *FTr);

   ElementTransformation *GetFaceTransformation(int FaceNo);

   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);

   void GetFaceElements (int Face, int *Elem1, int *Elem2);
   void GetFaceInfos (int Face, int *Inf1, int *Inf2);

   void CheckElementOrientation ();
   void CheckBdrElementOrientation ();

   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();

   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);
   void MoveVertices(const Vector &displacements);
   void GetVertices(Vector &vert_coord) const;
   void SetVertices(const Vector &vert_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; }
   // use the provided GridFunction as Nodes
   void NewNodes(GridFunction &nodes);

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

   void UniformRefinement();

   // 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 GetNumFineElems (int i);

   int GetRefinementType (int i);

   int GetFineElem (int i, int j);

   ElementTransformation * GetFineElemTrans (int i, int j);

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

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

   void PrintVTK(ostream &out);

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

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

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

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

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

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

   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);

   virtual ~Mesh();
};


// 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