text: /**
   * \copyright
   * Copyright (c) 2015, OpenGeoSys Community (http://www.opengeosys.org)
   *            Distributed under a Modified BSD License.
   *              See accompanying file LICENSE.txt or
   *              http://www.opengeosys.org/project/license
   *
   */
  
  /**
  * \file msh_mesh.h
  */
  
  /**************************************************************************
     MSHLib - Object:
     Task:
     Programing:
     08/2005 WW/OK Encapsulation from rf_ele_msh
     last modified
  **************************************************************************/
  #ifndef msh_mesh_INC
  #define msh_mesh_INC
  
  // C++
  #include <string>
  
  /** depreciated includes of GEOLib */
  #include "geo_lib.h"
  //#include "geo_pnt.h"
  #include "geo_sfc.h"
  #include "geo_vol.h"
  
  // GEOLIB
  #include "GEOObjects.h"
  #include "Point.h"
  #include "Polyline.h"
  #include "Surface.h"
  #include "Grid.h"
  
  // MSHLib
  #include "MSHEnums.h" // KR 2010/11/15
  #include "MeshNodesAlongPolyline.h"
  
  // FileIO
  #include "MeshIO/OGSMeshIO.h"
  
  #include "msh_elem.h"
  
  class RandomWalk;
  class CFluidMomentum;
  
  class Problem;
  
  #ifdef NEW_EQS // 1.11.2007 WW
  namespace Math_Group
  {
  class SparseTable;
  }
  using Math_Group::SparseTable;
  #endif
  
  //------------------------------------------------------------------------
  namespace MeshLib
  {
  /*!
     Class to handle topologic relationship among grids
     Designed by WW
     First version:    13.05.2009
   */
  class GridsTopo
  {
  private:
  	long* local_indices; // of border nodes in local node array of a grid
  	// double *comm_data;   // data for communication
  	std::string neighbor_name;
  	long bnodes;
  	friend class CFEMesh;
  	friend class ::CRFProcess;
  
  public:
  	GridsTopo(std::istream& in, std::string sec_name);
  	std::string getNeighbor_Name() const { return neighbor_name; }
  	long* getBorderNodeIndicies() const { return local_indices; }
  	long getBorderNodeNumber() const { return bnodes; }
  	// void Write(std::ostream &os=cout);
  	~GridsTopo();
  };
  
  /// For parallel computing. 03.2012. WW
  #if defined(USE_PETSC) // || defined(using other parallel scheme)
  typedef long MyInt;
  struct MeshNodes
  {
  	MyInt index;
  	double x;
  	double y;
  	double z;
  };
  
  #endif
  
  //------------------------------------------------------------------------
  // Class definition
  class CFEMesh
  {
  public:
  	/// Constructor using geometric information.
  	CFEMesh(GEOLIB::GEOObjects* geo_obj = NULL, std::string* unique_name = NULL);
  
  	/// Copy-Constructor.
  	/// Note that this is not a real copy-constructor. It copies only nodes and elements and calls ConstructGrid()
  	/// afterwards.
  	CFEMesh(CFEMesh const& mesh);
  
  	/// Destructor
  	~CFEMesh();
  
  	GEOLIB::GEOObjects* getGEOObjects() const { return _geo_obj; }
  	std::string* getProjectName() const { return _geo_name; }
  	/**
  	 * sets the value for element type
  	 * @param ele_type
  	 */
  	void setElementType(MshElemType::type ele_type);
  
  	/**
  	 * set the number of mesh layer
  	 * @param n_msh_layer
  	 */
  	void setNumberOfMeshLayers(size_t n_msh_layer);
  
  	/**
  	 * returns the number of mesh layers
  	 * @return the number of mesh layers
  	 */
  	size_t getNumberOfMeshLayers() const; // TF
  
  	/**
  	 *
  	 * @return
  	 */
  	bool hasCrossSection() const;
  
  	size_t getNumberOfLines() const;
  	size_t getNumberOfQuads() const;
  	size_t getNumberOfHexs() const;
  	size_t getNumberOfTris() const;
  	size_t getNumberOfTets() const;
  	size_t getNumberOfPrisms() const;
  	size_t getNumberOfPyramids() const;
  	double getMinEdgeLength() const;
  	CNode* const* getNodes() const // WW 05.2012.
  	{
  		return &nod_vector[0];
  	}
  	/**
  	 * do not use this method REMOVE CANDIDATE
  	 * @param val
  	 */
  	void setMinEdgeLength(double val) // TF temporary
  	{
  		_min_edge_length = val;
  	}
  
  	/**
  	 * Access method to the search length for geometric search algorithms computed
  	 * by method \c computeSearchLength().
  	 * @return the search length
  	 */
  	double getSearchLength() const;
  
  	/**
  	 * set search length to half of minimum edge length.
  	 * @set the search length
  	 */
  	void setSearchLength(double len);
  
  	/**
  	 * @brief Compute the search length for geometric search algorithms.
  	 *
  	 * Let \f$\mu\f$ the mean value of all edge length and \f$s\f$ the
  	 * standard deviation. The search length \f$\ell\f$ is computed by the formula
  	 * \f$\ell = \mu - c \cdot s.\f$
  	 * @param c (input) scaling constant, default value = 2
  	 */
  	void computeSearchLength(double c = 2);
  
  	/*!
  	   \brief Read mesh data.
  	   \param fem_file Stream for input, which has to read a line before calling this function
  	   \return Return true if there are other mesh to be read.
  	*/
  	bool Read(std::ifstream* fem_file);
  
  	friend class FileIO::OGSMeshIO;
  	std::ios::pos_type GMSReadTIN(std::ifstream*);
  	//
  	void ConstructGrid();
  	void GenerateHighOrderNodes();
  /// For parallel computing. 03.2012. WW
  #if defined(USE_PETSC) // || defined(other parallel solver libs)
  	void ConfigHighOrderElements();
  
  	/*!
  	   Fill data for subdomain mesh
  	       @param header  : mesh header
  	       @param s_nodes : mesh nodes
  	*/
  	void setSubdomainNodes(MyInt* header, const MeshNodes* s_nodes);
  	/*!
  	   Fill data for subdomain mesh
  	       @param header    : mesh header
  	       @param elem_info : element information
  	       @param inside    : indicator for elements that are inside the subdomain
  	*/
  	void setSubdomainElements(MyInt* header, const MyInt* elem_info, const bool inside);
  	int calMaximumConnectedNodes();
  	/// Get number of nodes of the entire mesh
  	int getNumNodesGlobal() const { return glb_NodesNumber_Linear; }
  	/// Get number of nodes of the entire mesh of quadratic elements
  	int getNumNodesGlobal_Q() const { return glb_NodesNumber_Quadratic; }
  	/// Get number of nodes of the subdomain mesh
  	int getNumNodesLocal() const { return loc_NodesNumber_Linear; }
  	/// Get number of nodes of the subdomain mesh of quadratic elements
  	int getNumNodesLocal_Q() const { return loc_NodesNumber_Quadratic; }
  	/// Get the largest ID of active nodes for higher order interpolation
  	size_t getLargestActiveNodeID_Quadratic() const
  	{
  		return static_cast<size_t>(NodesNumber_Linear + loc_NodesNumber_Quadratic - loc_NodesNumber_Linear);
  	}
  #endif
  
  	//
  	//         void RenumberNodesForGlobalAssembly();
  	// For number of nodes
  	int GetMaxElementDim() const { return max_ele_dim; }
  	void SwitchOnQuadraticNodes(bool quad) { useQuadratic = quad; }
  	bool getOrder() const { return useQuadratic; }
  	bool isAxisymmetry() const { return _axisymmetry; }
  	// Get number of nodes
  	// CMCD int to long
  	size_t GetNodesNumber(const bool quadr) const
  	{
  		if (quadr)
  			return NodesNumber_Quadratic;
  		else
  			return NodesNumber_Linear;
  	}
  
  	size_t NodesInUsage() const
  	{
  		if (useQuadratic)
  			return NodesNumber_Quadratic;
  		else
  			return NodesNumber_Linear;
  	}
  
  	void InitialNodesNumber() // WW
  	{
  		NodesNumber_Quadratic = NodesNumber_Linear = nod_vector.size();
  	}
  
  	void setNumberOfNodesFromNodesVectorSize() { NodesNumber_Linear = nod_vector.size(); }
  	/// Free the memory occupied by edges
  	void FreeEdgeMemory(); // 09.2012. WW
  
  	/**
  	 * @{
  	 * */
  	/*!
  	    brief Find the element by a point
  	*/
  	size_t FindElementByPoint(const double* xyz);
  
  	/**
  	 * \brief depreciated method - uses old surface class
  	 */
  	void GetNODOnPLY_XY(CGLPolyline* m_ply, std::vector<long>& msh_nod_vector);
  	//	/**
  	//	 * \brief depreciated method - uses old surface class
  	//	 */
  	void GetNODOnSFC_Vertical(Surface* m_sfc, std::vector<long>& msh_nod_vector);
  
  	/**
  	 * \brief depreciated method
  	 */
  	void CreateLineELEFromPLY(CGLPolyline*);
  	/**
  	 * \brief depreciated method
  	 */
  	void CreateLayerPolylines(CGLPolyline*); // OK
  
  	// GEO-SFC
  	/**
  	 * \brief depreciated method
  	 */
  	void GetNODOnSFC(Surface* m_sfc, std::vector<long>& msh_nod_vector, const bool for_s_term = false);
  	/**
  	 * \brief depreciated method
  	 */
  	void GetNODOnSFC_PLY(Surface const* m_sfc, std::vector<long>& msh_nod_vector, const bool for_s_term = false) const;
  	/**
  	 * \brief depreciated method
  	 */
  	void GetNODOnSFC_PLY_XY(Surface* m_sfc, std::vector<long>& msh_nod_vector,
  	                        bool givenNodesOnSurface = false); // givenNodeOnSurface by WW
  	/**
  	 * \brief depreciated method
  	 */
  	// 02.2009/OK
  	void GetNODOnSFC_PLY_Z(Surface*, std::vector<long>&);
  
  	/**
  	 * \brief depreciated method
  	 */
  	void GetNODOnSFC_TIN(Surface*, std::vector<long>&);
  	/**
  	 * \brief deprecated method
  	 */
  	void GetNodesOnCylindricalSurface(Surface* m_sfc, std::vector<long>& NodesS);
  	/**
  	 * \brief depreciated method
  	 */
  	void CreateQuadELEFromSFC(Surface*);
  
  	/**
  	 * GetNODOnPNT searchs the nearest node to the geometric point
  	 * */
  	long GetNODOnPNT(const GEOLIB::Point* const pnt) const;
  	/**
  	 * GetNearestELEOnPNT searchs the nearest element (gravity center)
  	 * to the geometric point
  	 * */
  	long GetNearestELEOnPNT(const GEOLIB::Point* const pnt) const;
  
  	/**
  	 * GetNODOnPLY search the nearest nodes along the Polyline object
  	 * @param ply constant pointer to a constant Polyline object
  	 * @param msh_nod_vector the mesh node indices are saved in this vector
  	 * @param automatic use the computed search length
  	 * @param eps if automatic is false use eps as search length
  	 * */
  	void GetNODOnPLY(const GEOLIB::Polyline* const ply,
  	                 std::vector<size_t>& msh_nod_vector,
  	                 bool automatic = true,
  	                 double eps = std::numeric_limits<double>::epsilon());
  
  	/**
  	 *
  	 * @param ply
  	 * @return
  	 */
  	const MeshNodesAlongPolyline& GetMeshNodesAlongPolyline(const GEOLIB::Polyline* const ply);
  
  	/**
  	 *
  	 * @param ply
  	 * @param points
  	 */
  	void getPointsForInterpolationAlongPolyline(const GEOLIB::Polyline* const ply, std::vector<double>& points);
  
  	/**
  	 * GetNODOnPLY search the nearest nodes to the Polyline
  	 * */
  	void GetNODOnPLY(const GEOLIB::Polyline* const ply,
  	                 std::vector<long>& msh_nod_vector,
  	                 const bool for_s_term = false,
  	                 bool automatic = true,
  	                 double search_radius = std::numeric_limits<double>::epsilon());
  
  	/**
  	 * \brief gives the indices of CElement elements, which have an edge
  	 * in common with the polyline.
  	 */
  	void GetELEOnPLY(const GEOLIB::Polyline*, std::vector<size_t>&, bool With1DElements); // 11/2011 BG
  
  	/**
  	 * \brief gives the indices of nodes, which are contained in the surface
  	 */
  	void GetNODOnSFC(const GEOLIB::Surface* sfc,
  	                 std::vector<size_t>& msh_nod_vector,
  	                 const bool for_s_term = false) const;
  
  	/** @} */ // close doxygen group
  
  //....................................................................
  #ifdef ObsoleteGUI // WW 03.2012
  	// QUAD->HEX
  	void CreateHexELEFromQuad(int, double);
  	// QUAD->LINE
  	void CreateLineELEFromQuad(int, double, int);
  	void SetActiveElements(std::vector<long>&);
  	// MB
  	void SetMSHPart(std::vector<long>&, long);
  	bool NodeExists(size_t node);
  	// LINE->LINE
  	void AppendLineELE();
  
  	// TRI->PRIS
  	void CreatePriELEFromTri(int, double);
  
  #endif
  
  	// Coordinate system
  	int GetCoordinateFlag() const { return coordinate_system; }
  	void FillTransformMatrix();
  
  	void RenumberNodesForGlobalAssembly();
  	/**
  	 * returns the vector storing pointers to all nodes (class CNode) of the mesh
  	 * @return
  	 */
  	const std::vector<MeshLib::CNode*>& getNodeVector() const { return nod_vector; }
  	// All nodes - should be private!!!
  	std::vector<MeshLib::CNode*> nod_vector;
  
  	// All edges
  	std::vector<MeshLib::CEdge*> edge_vector;
  	// All surface faces
  	std::vector<MeshLib::CElem*> face_vector;
  	// All surface normal
  	std::vector<double*> face_normal; // YD
  
  	const std::vector<MeshLib::CElem*>& getElementVector() const { return ele_vector; }
  	/**
  	 * all elements stored in this vector
  	 * */
  	std::vector<MeshLib::CElem*> ele_vector;
  	// Nodes in usage
  	// To record eqs_index->global node index
  	std::vector<long> Eqs2Global_NodeIndex;
  
  	// OK
  	void PrismRefine(int Layer, int subdivision);
  
  	void ConnectedNodes(bool quadratic) const;
  	// WW
  	void ConnectedElements2Node(bool quadratic = false);
  	// OK
  	std::vector<std::string> mat_names_vector;
  	void DefineMobileNodes(CRFProcess*); // OK
  	void FaceNormal(); // YD
  	void SetNODPatchAreas(); // OK4310
  	void SetNetworkIntersectionNodes(); // OK4319->PCH
  
  	void GetConnectedElements(std::vector<long>& nodes_on_sfc, std::vector<long>& vec_elements);
  
  #ifdef NEW_EQS // 1.11.2007 WW
  	// Compute the graph of the sparse matrix related to this mesh. 1.11.2007 WW
  	void CreateSparseTable();
  	// Get the sparse graph   1.11.2007 WW
  	SparseTable* GetSparseTable(bool quad = false) const
  	{
  		if (!quad)
  			return sparse_graph;
  		else
  			return sparse_graph_H;
  	}
  #endif
  
  	std::string pcs_name;
  	std::string geo_name; // MB
  	std::string geo_type_name; // OK10_4310
  
  	size_t max_mmp_groups; // OKCC
  	int msh_max_dim;
  
  	RandomWalk* PT; // PCH
  
  	CFluidMomentum* fm_pcs; // by PCH
  
  	std::vector<size_t> sorted_nodes;
  	std::vector<size_t> xy_change;
  	bool nodes_are_sorted;
  
  	/// Import MODFlow grid. 10.2009 WW
  	void ImportMODFlowGrid(std::string const& fname);
  	/// Convert raster cells into grid. 12.2009 WW
  	void ConvertShapeCells(std::string const& fname);
  #ifdef USE_HydSysMshGen
  	// Be activated if it is still needed.
  	// WWmesh_header
  	/// Generate Column-surface grid system for the modeling of surafce-subsuface coupled processes
  	// 15.05.2009. WW // removed useless const TF
  	void HydroSysMeshGenerator(std::string fname, int nlayers, double thickness, int mapping);
  
  #endif
  	void MarkInterface_mHM_Hydro_3D(); // 07.06.2010. WW
  	void mHM2NeumannBC();
  	/// Comptute int {f} a dA on top surface.
  	void TopSurfaceIntegration();
  
  #ifndef NDEBUG
  	/**
  	 * This is a getter method to access the private attribute _mesh_grid
  	 * that is an instance of class Grid.
  	 * @return
  	 */
  	GEOLIB::Grid<MeshLib::CNode> const* getGrid() const;
  #endif
  
  private:
  	// private attributes
  	/**
  	 * reference to object of class GEOObject, that manages the geometry data
  	 */
  	GEOLIB::GEOObjects* _geo_obj;
  	/**
  	 * identifier for geometry
  	 */
  	std::string* _geo_name;
  	std::vector<MeshLib::MeshNodesAlongPolyline> _mesh_nodes_along_polylines;
  
  	MshElemType::type _ele_type;
  	size_t _n_msh_layer; // OK
  	bool _cross_section;
  	size_t _msh_n_lines;
  	size_t _msh_n_quads;
  	size_t _msh_n_hexs;
  	size_t _msh_n_tris;
  	size_t _msh_n_tets;
  	size_t _msh_n_prisms;
  	size_t _msh_n_pyras;
  
  	/**
  	 * method initializes the minimal edge length that is used in search algorithms
  	 */
  	void computeMinEdgeLength();
  	double _min_edge_length; // TK
  
  	/**
  	 * length for geometrical search algorithms, calculated in method
  	 * \c computeSearchLength() and returned by method \c getSearchLength()
  	 */
  	double _search_length;
  
  	// Process friends
  	friend class ::CRFProcess;
  
  	size_t NodesNumber_Linear;
  	size_t NodesNumber_Quadratic;
  /// For parallel computing. 03.2012. WW
  #if defined(USE_PETSC) // || defined(using other parallel scheme)
  	// int n_sub_elements;
  	int glb_NodesNumber_Linear;
  	int glb_NodesNumber_Quadratic;
  	int loc_NodesNumber_Linear; // index of shadow nodes starts from this number
  	int loc_NodesNumber_Quadratic;
  
  #endif
  	bool useQuadratic;
  	bool _axisymmetry;
  	bool top_surface_checked; // 07.06.2010.  WW
  
  	// Coordinate indicator
  	// 10:  X component only
  	// 11: Y component only
  	// 12: Z component only
  	// 21:  X, Y component
  	// 22:  X, Z component
  	// 23:  Y, Z component
  	// 32:  X, Y, Z component
  	int coordinate_system;
  	bool has_multi_dim_ele;
  	int max_ele_dim;
  	int map_counter; // 21.01.2009 WW
  
  	bool mapping_check; // 23.01.2009 WW
  	/// Import shape file. 16.03.2026. WW
  	size_t ncols, nrows;
  	/// (x_0, y_0): coordinate of the left down corner
  	double x0, y0, csize, ndata_v;
  	std::vector<double> zz; // Elevation
  	inline void ReadShapeFile(std::string const& fname);
  	// 03.2010. WW
  	inline void Precipitation2NeumannBC(std::string const& fname, std::string const& ofname, double ratio = 0.8);
  
  	/// Store border nodes among different grids.
  	std::vector<GridsTopo*> grid_neighbors;
  	friend class ::Problem;
  
  //
  // Sparse graph of this mesh. 1.11.2007 WW
  #ifdef NEW_EQS
  	SparseTable* sparse_graph;
  	SparseTable* sparse_graph_H; // For high order interpolation
  #endif
  
  	void CreateLineElementsFromMarkedEdges(CFEMesh* m_msh_ply, std::vector<long>& ele_vector_at_ply); // NW
  
  	/// Find nodes in convex polygon
  	void findNodesInPolygon(const double area_orig, const double tol, const size_t start_id, const size_t end_id,
  	                        const CGLPolyline* ply, std::vector<long>& node_id_vector) const;
  
  public:
  	void constructMeshGrid();
  
  private:
  	GEOLIB::Grid<MeshLib::CNode>* _mesh_grid;
  };
  
  } // namespace MeshLib
  
  #endif
  