First cut of GradateMesh, moved SetReferenceMesh into fem_adapt_algs
[charm.git] / src / libs / ck-libs / fem / fem_adapt_algs.h
1 /* File: fem_adapt_algs.h
2  * Authors: Terry Wilmarth, Nilesh Choudhury
3  * 
4  */
5
6 // This module implements high level mesh adaptivity algorithms that make use 
7 // of the primitive mesh adaptivity operations provided by fem_adapt(_new).
8 // Ask: TLW
9 #ifndef __CHARM_FEM_ADAPT_ALGS_H
10 #define __CHARM_FEM_ADAPT_ALGS_H
11
12 #include "charm++.h"
13 #include "tcharm.h"
14 #include "charm-api.h"
15 #include "ckvector3d.h"
16 #include "fem.h"
17 #include "fem_mesh.h"
18 #include "vector2d.h"  // **CW** tentative
19
20 #define SLIVERAREA 1.0e-18
21 #define REFINE_TOL 1.5
22 #define COARSEN_TOL 0.85
23 #define QUALITY_MIN 0.6
24
25 class FEM_Adapt_Algs;
26 CtvExtern(FEM_Adapt_Algs *, _adaptAlgs);
27
28 class femMeshModify;
29 class FEM_Adapt;
30 class FEM_AdaptL;
31
32 class FEM_Adapt_Algs {
33   friend class FEM_AdaptL;
34   friend class FEM_Adapt;
35   friend class femMeshModify;
36   friend class FEM_Interpolate;
37   friend class FEM_MUtil;
38
39  protected: 
40   FEM_Mesh *theMesh;
41   femMeshModify *theMod;
42   //FEM_Adapt *theAdaptor;
43   FEM_AdaptL *theAdaptor;
44   int numNodes, numElements, dim;
45   int coord_attr;
46   int bc_attr;
47   // These are for element sorting
48   typedef struct {
49     int elID;
50     double len;
51   } elemHeap;
52   elemHeap *coarsenElements;
53   elemHeap *refineElements;
54   elemHeap *refineStack;
55   int refineTop, refineHeapSize, coarsenHeapSize;
56
57  public:
58   FEM_Adapt_Algs() {
59     theMesh = NULL; theMod = NULL; theAdaptor = NULL;
60   }
61   /// Initialize FEM_Adapt_Algs with a chunk of the mesh
62   FEM_Adapt_Algs(FEM_Mesh *m, femMeshModify *fm, int dimension);
63   void FEM_Adapt_Algs_Init(int coord_at, int bc_at) {
64     coord_attr = coord_at;
65     bc_attr = bc_at;
66   }
67   /// Perform refinements on a mesh
68   /** Perform refinements on a mesh.  Tries to maintain/improve element quality
69       as specified by a quality measure qm;
70       if method = 0, refine areas with size larger than factor down to factor
71       if method = 1, refine elements down to sizes specified in sizes array
72       Negative entries in size array indicate no refinement. */
73   void FEM_Refine(int qm, int method, double factor, double *sizes);
74   /// Perform coarsening on a mesh
75   /** Perform coarsening on a mesh.  Tries to maintain/improve element quality
76       as specified by a quality measure qm;
77       if method = 0, coarsen areas with size smaller than factor up to factor
78       if method = 1, coarsen elements up to sizes specified in sizes array
79       Negative entries in size array indicate no coarsening. */
80   void FEM_Coarsen(int qm, int method, double factor, double *sizes);
81   /// Smooth the mesh using method according to some quality measure qm
82   void FEM_Smooth(int qm, int method);
83   /// Repair the mesh according to some quality measure qm
84
85   // FEM_Mesh_mooth
86   //    Inputs  : meshP - a pointer to the FEM_Mesh object to smooth
87   //            : nodes - an array of local node numbers to be smoothed.  Send
88   //                      NULL pointer to smooth all nodes.
89   //            : nNodes - the size of the nodes array
90   //            : attrNo - the attribute number where the coords are registered
91   //    Shifts nodes around to improve mesh quality.  FEM_BOUNDARY attribute
92   //    and interpolator function must be registered by user to maintain 
93   //    boundary information.
94   void FEM_mesh_smooth(FEM_Mesh *meshP, int *nodes, int nNodes, int attrNo);
95
96   void FEM_Repair(int qm);
97   /// Remesh entire mesh
98   /** Remesh entire mesh according to quality measure qm
99       if method = 0, set entire mesh size to factor
100       if method = 1, keep regional mesh sizes, and scale by factor
101       if method = 2, uses sizes to size mesh by regions */
102   void FEM_Remesh(int qm, int method, double factor, double *sizes);
103   
104   /// Set sizes on mesh elements based on their average edge length
105   void SetReferenceMesh();
106   /// Adjust sizes on mesh elements to avoid sharp discontinuities
107   void GradateMesh(double smoothness);
108  private:
109   // Helper methods
110   /// Performs refinement; returns number of modifications
111   int Refine(int qm, int method, double factor, double *sizes);
112   /// Performs coarsening; returns number of modifications
113   int Coarsen(int qm, int method, double factor, double *sizes);
114   /// Set sizes on elements throughout the mesh; note: size is edge length
115   void SetMeshSize(int method, double factor, double *sizes);
116   /// Insert element to be refined/coarsened
117   void Insert(int elID, double len, int cFlag);
118   /// Get next element to be refined/coarsened
119   int Delete_Min(int cFlag);
120  public:
121   /// Initiate instance of Longest Edge Bisection on an element
122   /** Initiate instance of Longest Edge Bisection on element e.  Propagates
123       throughout the mesh to maintain the requirement that only longest edges
124       are bisected; returns 1 if successful, 0 if not **/
125   virtual int refine_element_leb(int e);
126   virtual void refine_flip_element_leb(int e, int p, int n1, int n2, 
127                                        double le);
128
129   int simple_refine(double targetA, double xmin=0.0, double ymin=0.0, double xmax=1.0, double ymax=1.0);
130   int simple_coarsen(double targetA, double xmin=0.0, double ymin=0.0, double xmax=1.0, double ymax=1.0);
131   double length(int n1, int n2);
132   double getArea(int n1, int n2, int n3);
133   double length(double *n1_coord, double *n2_coord);
134   double getArea(double *n1_coord, double *n2_coord, double *n3_coord);
135   int getCoord(int n1, double *crds);
136   int getShortestEdge(int n1, int n2, int n3, int* shortestEdge);
137   double getAreaQuality(int elem);
138   bool didItFlip(int n1, int n2, int n3, double *n4_coord);
139   bool didItFlip(double *n1_coord, double *n2_coord, double *n3_coord, double *n4_coord);
140   double getSignedArea(double *n1_coord, double *n2_coord, double *n3_coord);
141   void tests(void);
142 };
143
144
145 #endif