ncmesh.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_NCMESH
#define MFEM_NCMESH

#include <vector>
#include <map>
#include <iostream>

#include "../config/config.hpp"
#include "../general/hash.hpp"
#include "../linalg/densemat.hpp"
#include "element.hpp"
#include "vertex.hpp"
#include "../fem/geom.hpp"

namespace mfem
{

struct Refinement
{
   int index; 
   char ref_type; 

   Refinement(int index, int type = 7) : index(index), ref_type(type) {}
};

struct Embedding
{
   int parent; 
   int matrix; 

   Embedding(int elem, int matrix = 0) : parent(elem), matrix(matrix) {}
};

struct CoarseFineTransformations
{
   DenseTensor point_matrices;  
   Array<Embedding> embeddings; 
};


class NCMesh
{
protected:
   struct Element; // forward

public:

   NCMesh(const Mesh *mesh, std::istream *vertex_parents = NULL);

   NCMesh(const NCMesh &other);

   virtual ~NCMesh();

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

   virtual void Refine(const Array<Refinement> &refinements);

   virtual void LimitNCLevel(int max_nc_level);

   const Table &GetDerefinementTable();

   virtual void CheckDerefinementNCLevel(const Table &deref_table,
                                         Array<int> &level_ok, int max_nc_level);

   virtual void Derefine(const Array<int> &derefs);


   // master/slave lists

   struct MeshId
   {
      int index; 
      int local; 
      Element* element; 

      MeshId(int index = -1, Element* element = NULL, int local = -1)
         : index(index), local(local), element(element) {}
   };

   struct Master : public MeshId
   {
      int slaves_begin, slaves_end; 

      Master(int index, Element* element, int local, int sb, int se)
         : MeshId(index, element, local), slaves_begin(sb), slaves_end(se) {}
   };

   struct Slave : public MeshId
   {
      int master; 
      int edge_flags; 
      DenseMatrix point_matrix; 

      Slave(int index, Element* element, int local)
         : MeshId(index, element, local), master(-1), edge_flags(0) {}

      void OrientedPointMatrix(DenseMatrix &oriented_matrix) const;
   };

   struct NCList
   {
      std::vector<MeshId> conforming;
      std::vector<Master> masters;
      std::vector<Slave> slaves;
      // TODO: switch to Arrays when fixed for non-POD types
      // TODO: make a list of unique slave matrices to save memory (+ time later)

      void Clear() { conforming.clear(); masters.clear(); slaves.clear(); }
      bool Empty() const { return !conforming.size() && !masters.size(); }
      long MemoryUsage() const;
   };

   const NCList& GetFaceList()
   {
      if (face_list.Empty()) { BuildFaceList(); }
      return face_list;
   }

   const NCList& GetEdgeList()
   {
      if (edge_list.Empty()) { BuildEdgeList(); }
      return edge_list;
   }


   // coarse/fine transforms

   void MarkCoarseLevel();

   const CoarseFineTransformations& GetRefinementTransforms();

   const CoarseFineTransformations& GetDerefinementTransforms();

   void ClearTransforms();


   // utility

   int GetEdgeMaster(int v1, int v2) const;

   virtual void GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
                                   Array<int> &bdr_vertices,
                                   Array<int> &bdr_edges);

   int GetElementGeometry() const { return root_elements[0]->geom; }

   int GetElementDepth(int i) const;

   void PrintVertexParents(std::ostream &out) const;

   void PrintCoarseElements(std::ostream &out) const;

   void LoadVertexParents(std::istream &input);

   void LoadCoarseElements(std::istream &input);

   void SetVertexPositions(const Array<mfem::Vertex> &vertices);

   long MemoryUsage() const;

   void PrintMemoryDetail() const;


protected: // interface for Mesh to be able to construct itself from NCMesh

   friend class Mesh;

   void GetMeshComponents(Array<mfem::Vertex>& vertices,
                          Array<mfem::Element*>& elements,
                          Array<mfem::Element*>& boundary) const;

   virtual void OnMeshUpdated(Mesh *mesh);


protected: // implementation

   int Dim, spaceDim; 
   bool Iso; 

   Element* CopyHierarchy(Element* elem);
   void DeleteHierarchy(Element* elem);


   // primary data

   struct RefCount
   {
      int ref_count;

      RefCount() : ref_count(0) {}

      int Ref()
      {
         return ++ref_count;
      }
      int Unref()
      {
         int ret = --ref_count;
         if (!ret) { delete this; }
         return ret;
      }
   };

   struct Vertex : public RefCount
   {
      int index;     
      double pos[3]; 

      Vertex() {}
      Vertex(double x, double y, double z) : index(-1)
      { pos[0] = x, pos[1] = y, pos[2] = z; }
      Vertex(const Vertex& other) { std::memcpy(this, &other, sizeof(*this)); }
   };

   struct Edge : public RefCount
   {
      int attribute; 
      int index;     

      Edge() : attribute(-1), index(-1) {}
      Edge(const Edge &other) { std::memcpy(this, &other, sizeof(*this)); }

      bool Boundary() const { return attribute >= 0; }
   };

   struct Node : public Hashed2<Node>
   {
      Vertex* vertex;
      Edge* edge;

      Node(int id) : Hashed2<Node>(id), vertex(NULL), edge(NULL) {}
      Node(const Node &other);
      ~Node();

      // Bump ref count on a vertex or an edge, or create them. Used when an
      // element starts using a vertex or an edge.
      void RefVertex();
      void RefEdge();

      // Decrement ref on vertex or edge when an element is not using them
      // anymore. The vertex, edge or the whole Node can autodestruct.
      // (The hash-table pointer needs to be known then to remove the node.)
      void UnrefVertex(HashTable<Node>& nodes);
      // Implement as a static method since we check if node == NULL.
      static void UnrefEdge(Node *node, HashTable<Node>& nodes);
   };

   struct Face : public RefCount, public Hashed4<Face>
   {
      int attribute;    
      int index;        
      Element* elem[2]; 

      Face(int id);
      Face(const Face& other);

      bool Boundary() const { return attribute >= 0; }

      // add or remove an element from the 'elem[2]' array
      void RegisterElement(Element* e);
      void ForgetElement(Element* e);

      Element* GetSingleElement() const;

      // overloaded Unref without auto-destruction
      int Unref() { return --ref_count; }
   };

   struct Element
   {
      char geom;     
      char ref_type; 
      char flag;     
      int index;     
      int rank;      
      int attribute;
      union
      {
         Node* node[8];  
         Element* child[8]; 
      };
      Element* parent; 

      Element(int geom, int attr);
      Element(const Element& other) { std::memcpy(this, &other, sizeof(*this)); }
   };

   // TODO: use MemAlloc or similar to store Element, Face, Node, etc., replace
   // most pointers with indices (ints) into the allocators to save memory;
   // ref_count and attributes could be chars

   Array<Element*> root_elements; // coarsest mesh, initialized by constructor

   HashTable<Node> nodes; // associative container holding all Nodes
   HashTable<Face> faces; // associative container holding all Faces


   // secondary data

   virtual void Update();

   Array<Element*> leaf_elements; // finest level, updated by UpdateLeafElements

   Array<int> vertex_nodeId; // vertex-index to node-id map, see UpdateVertices

   NCList face_list; 
   NCList edge_list; 

   Array<Face*> boundary_faces; 
   Array<Node*> boundary_edges; 

   Table element_vertex; 
   int num_vertices;     

   virtual void UpdateVertices(); 

   void CollectLeafElements(Element* elem, int state);
   void UpdateLeafElements();

   virtual void AssignLeafIndices();

   virtual bool IsGhost(const Element* elem) const { return false; }
   virtual int GetNumGhosts() const { return 0; }


   // refinement/derefinement

   struct ElemRefType
   {
      Element* elem;
      int ref_type;

      ElemRefType(Element* elem, int type)
         : elem(elem), ref_type(type) {}
   };

   Array<ElemRefType> ref_stack; 

   Table derefinements; 

   void RefineElement(Element* elem, char ref_type);
   void DerefineElement(Element* elem);

   Element* NewHexahedron(Node* n0, Node* n1, Node* n2, Node* n3,
                          Node* n4, Node* n5, Node* n6, Node* n7,
                          int attr,
                          int fattr0, int fattr1, int fattr2,
                          int fattr3, int fattr4, int fattr5);

   Element* NewQuadrilateral(Node* n0, Node* n1, Node* n2, Node* n3,
                             int attr,
                             int eattr0, int eattr1, int eattr2, int eattr3);

   Element* NewTriangle(Node* n0, Node* n1, Node* n2,
                        int attr, int eattr0, int eattr1, int eattr2);

   Vertex* NewVertex(Node* v1, Node* v2);

   mfem::Element* NewMeshElement(int geom) const;

   Node* GetMidEdgeVertex(Node* v1, Node* v2);
   Node* GetMidEdgeVertexSimple(Node* v1, Node* v2);
   Node* GetMidFaceVertex(Node* e1, Node* e2, Node* e3, Node* e4);

   int FaceSplitType(Node* v1, Node* v2, Node* v3, Node* v4,
                     Node* mid[4] = NULL /* optional output of mid-edge nodes*/)
   const;

   void ForceRefinement(Node* v1, Node* v2, Node* v3, Node* v4);

   void CheckAnisoFace(Node* v1, Node* v2, Node* v3, Node* v4,
                       Node* mid12, Node* mid34, int level = 0);

   void CheckIsoFace(Node* v1, Node* v2, Node* v3, Node* v4,
                     Node* e1, Node* e2, Node* e3, Node* e4, Node* midf);

   void RefElementNodes(Element *elem);
   void UnrefElementNodes(Element *elem);
   Face* GetFace(Element* elem, int face_no);
   void RegisterFaces(Element* elem, int *fattr = NULL);

   Node* PeekAltParents(Node* v1, Node* v2);

   bool NodeSetX1(Node* node, Node** n);
   bool NodeSetX2(Node* node, Node** n);
   bool NodeSetY1(Node* node, Node** n);
   bool NodeSetY2(Node* node, Node** n);
   bool NodeSetZ1(Node* node, Node** n);
   bool NodeSetZ2(Node* node, Node** n);

   void CollectDerefinements(Element* elem, Array<Connection> &list);

   // face/edge lists

   static int find_node(Element* elem, Node* node);
   static int find_node(Element* elem, int node_id);
   static int find_element_edge(Element* elem, int v0, int v1);
   static int find_hex_face(int a, int b, int c);

   int ReorderFacePointMat(Node* v0, Node* v1, Node* v2, Node* v3,
                           Element* elem, DenseMatrix& mat) const;
   struct PointMatrix;

   void TraverseFace(Node* v0, Node* v1, Node* v2, Node* v3,
                     const PointMatrix& pm, int level);

   void TraverseEdge(Node* v0, Node* v1, double t0, double t1, int flags,
                     int level);

   virtual void BuildFaceList();
   virtual void BuildEdgeList();

   virtual void ElementSharesEdge(Element* elem, Edge* edge) {} // ParNCMesh
   virtual void ElementSharesFace(Element* elem, Face* face) {} // ParNCMesh


   // neighbors / element_vertex table

   void FindSetNeighbors(const Array<char> &elem_set,
                         Array<Element*> *neighbors, /* append */
                         Array<char> *neighbor_set = NULL);

   void FindNeighbors(const Element* elem,
                      Array<Element*> &neighbors, /* append */
                      const Array<Element*> *search_set = NULL);

   void NeighborExpand(const Array<Element*> &elements,
                       Array<Element*> &expanded,
                       const Array<Element*> *search_set = NULL);


   void CollectEdgeVertices(Node *v0, Node *v1, Array<int> &indices);
   void CollectFaceVertices(Node* v0, Node* v1, Node* v2, Node* v3,
                            Array<int> &indices);
   void BuildElementToVertexTable();

   void UpdateElementToVertexTable()
   {
      if (element_vertex.Size() < 0) { BuildElementToVertexTable(); }
   }


   // coarse/fine transformations

   struct Point
   {
      int dim;
      double coord[3];

      Point() { dim = 0; }

      Point(double x, double y)
      { dim = 2; coord[0] = x; coord[1] = y; }

      Point(double x, double y, double z)
      { dim = 3; coord[0] = x; coord[1] = y; coord[2] = z; }

      Point(const Point& p0, const Point& p1)
      {
         dim = p0.dim;
         for (int i = 0; i < dim; i++)
         {
            coord[i] = (p0.coord[i] + p1.coord[i]) * 0.5;
         }
      }

      Point(const Point& p0, const Point& p1, const Point& p2, const Point& p3)
      {
         dim = p0.dim;
         MFEM_ASSERT(p1.dim == dim && p2.dim == dim && p3.dim == dim, "");
         for (int i = 0; i < dim; i++)
         {
            coord[i] = (p0.coord[i] + p1.coord[i] + p2.coord[i] + p3.coord[i])
                       * 0.25;
         }
      }

      Point& operator=(const Point& src)
      {
         dim = src.dim;
         for (int i = 0; i < dim; i++) { coord[i] = src.coord[i]; }
         return *this;
      }
   };

   struct PointMatrix
   {
      int np;
      Point points[8];

      PointMatrix(const Point& p0, const Point& p1, const Point& p2)
      { np = 3; points[0] = p0; points[1] = p1; points[2] = p2; }

      PointMatrix(const Point& p0, const Point& p1, const Point& p2, const Point& p3)
      { np = 4; points[0] = p0; points[1] = p1; points[2] = p2; points[3] = p3; }

      PointMatrix(const Point& p0, const Point& p1, const Point& p2,
                  const Point& p3, const Point& p4, const Point& p5,
                  const Point& p6, const Point& p7)
      {
         np = 8;
         points[0] = p0; points[1] = p1; points[2] = p2; points[3] = p3;
         points[4] = p4; points[5] = p5; points[6] = p6; points[7] = p7;
      }

      Point& operator()(int i) { return points[i]; }
      const Point& operator()(int i) const { return points[i]; }

      void GetMatrix(DenseMatrix& point_matrix) const;
   };

   static PointMatrix pm_tri_identity;
   static PointMatrix pm_quad_identity;
   static PointMatrix pm_hex_identity;

   static const PointMatrix& GetGeomIdentity(int geom);

   void GetPointMatrix(int geom, const char* ref_path, DenseMatrix& matrix);

   typedef std::map<std::string, int> RefPathMap;

   void TraverseRefinements(Element* elem, int coarse_index,
                            std::string &ref_path, RefPathMap &map);

   CoarseFineTransformations transforms;

   Array<Element*> coarse_elements;

   void InitDerefTransforms();
   void SetDerefMatrixCodes(Element* parent, Array<Element*> &coarse);


   // utility

   Node* GetEdgeMaster(Node* node) const;

   static void find_face_nodes(const Face *face, Node* node[4]);

   int  EdgeSplitLevel(Node* v1, Node* v2) const;
   void FaceSplitLevel(Node* v1, Node* v2, Node* v3, Node* v4,
                       int& h_level, int& v_level) const;

   void CountSplits(Element* elem, int splits[3]) const;
   void GetLimitRefinements(Array<Refinement> &refinements, int max_level);

   int CountElements(Element* elem) const;

   int PrintElements(std::ostream &out, Element* elem, int &coarse_id) const;

   void CountObjects(int &nelem, int &nvert, int &nedges) const;


public: // TODO: maybe make this part of mfem::Geometry?

   struct GeomInfo
   {
      int nv, ne, nf, nfv; // number of: vertices, edges, faces, face vertices
      int edges[12][2];    // edge vertices (up to 12 edges)
      int faces[6][4];     // face vertices (up to 6 faces)

      bool initialized;
      GeomInfo() : initialized(false) {}
      void Initialize(const mfem::Element* elem);
   };

   static GeomInfo GI[Geometry::NumGeom];

#ifdef MFEM_DEBUG
public:
   void DebugNeighbors(Array<char> &elem_set);

   void DebugLeafOrder() const;
#endif

   friend class ParNCMesh/*::ElementSet*/;
};

}

#endif