TMR: Use standard C++ headers
[charm.git] / src / libs / ck-libs / tmr / tri.h
1 // Triangular Mesh Refinement Framework - 2D (TMR)
2 // Created by: Terry L. Wilmarth
3
4 #include <math.h>
5 #include <vector>
6 #include "charm++.h"
7 #include "tcharm.h"
8 #include "charm-api.h"
9
10 // Constants to tell FEM interface whether node is on a boudary between chunks
11 // and if it is the first of two split operations
12 #define LOCAL_FIRST 0x2
13 #define LOCAL_SECOND 0x0
14 #define BOUND_FIRST 0x3
15 #define BOUND_SECOND 0x1
16
17 class node;
18 class chunk;
19 class elemRef;
20
21 // --------------------------- Helper Classes -------------------------------
22 // objRef: References to mesh data require a chunk ID (cid) and an
23 // index on that chunk (idx). Subclasses for nodeRefs, edgeRefs and
24 // elemRefs define the same operators as the data they reference.
25 class objRef {
26  public:
27   int cid, idx;
28   objRef() { cid = -1;  idx = -1; }
29   objRef(int chunkId, int objIdx) { cid = chunkId; idx = objIdx; }
30   void init() { cid = -1; idx = -1; }
31   void init(int chunkId, int objIdx) { cid = chunkId; idx = objIdx; }
32   int operator==(const objRef& o) const { return((cid == o.cid) && (idx == o.idx)); }
33   int operator!=(const objRef& o) const { return !( (*this)==o ); }
34   int isNull(void) const {return cid==-1;}
35   void sanityCheck(chunk *C);
36   void pup(PUP::er &p) { p(cid); p(idx); }
37 };
38
39 // edge and element References: see actual classes below for
40 // method descriptions
41 class edgeRef : public objRef {
42  public:
43   edgeRef() :objRef() {}
44   edgeRef(int c,int i) :objRef(c,i) {}
45
46   void updateElement(chunk *C, elemRef oldval, elemRef newval);
47   int lock(chunk *C);
48   void unlock(chunk *C);
49   int locked(chunk *C) const;
50 };
51
52 class elemRef : public objRef {
53  public:
54   elemRef() : objRef() {}
55   elemRef(int c, int i) : objRef(c,i) {}
56
57   int checkIfLongEdge(chunk *C, edgeRef e);
58   double getArea(chunk *C);
59   void setTargetArea(chunk *C, double ta);
60   void updateEdges(chunk *C, edgeRef e0, edgeRef e1, edgeRef e2);
61   void unsetDependency(chunk *C);
62   void setDependent(chunk *C, int anIdx, int aCid);
63   int hasDependent(chunk *C);
64 };
65
66 #include "refine.decl.h"
67 // ------------------------ Global Read-only Data ---------------------------
68 extern CProxy_chunk mesh;
69 CtvExtern(chunk *, _refineChunk);
70
71 // ------------------------------ Messages ----------------------------------
72 // chunkMsg: information needed at startup
73 class chunkMsg : public CMessage_chunkMsg {
74  public:
75   int nChunks;
76   CProxy_TCharm myThreads;
77 };
78
79 // refMsg: generic message for sending/receiving a reference to/from an element
80 class refMsg : public CMessage_refMsg {
81  public:
82   objRef aRef;
83   int idx;
84 };
85
86 // doubleMsg: used to send a double to an element idx
87 class doubleMsg : public CMessage_doubleMsg {
88  public:
89   double aDouble;
90 };
91
92 // intMsg: used to return integers from sync methods
93 class intMsg : public CMessage_intMsg {
94  public:
95   int anInt;
96 };
97
98 // ---------------------------- Data Classes --------------------------------
99
100 class node {
101   // node: a coordinate
102
103   // Each node has coordinates (x, y) and a local pointer (C).
104   double x, y;  
105   chunk *C;
106
107  public:
108   // Basic node constructor
109   node() { C = NULL; }
110   node(double a, double b) { init(a,b); }
111
112   // Initializers for a node: can set coordinate separately, basic
113   // chunk and index info separately, or everything at once. 
114   void init() { x = -1.0; y = -1.0; C = NULL; }
115   void init(double a, double b) { x = a; y = b; }
116   void init(chunk *cPtr) { C = cPtr; }
117   void init(double a, double b, chunk *cPtr) { x = a; y = b; C = cPtr; }
118
119   // Assignment
120   node& operator=(const node& n) { x = n.x; y = n.y; return *this; }
121   // Equality
122   int operator==(const node& n) const { return ((x == n.x) && (y == n.y)); }
123
124   // Access X and Y coordinates
125   double X() const { return x; }
126   double Y() const { return y; }
127
128   // distance: compute distance between this and another node
129   double distance(const node& n) const {
130     double dx = n.x - x, dy = n.y - y;
131     return (sqrt ((dx * dx) + (dy * dy)));
132   }
133
134   // midpoint: compute midpoint between this and another node
135   void midpoint(const node& n, node *result) const {
136     result->x = (x + n.x) / 2.0;  result->y = (y + n.y) / 2.0;
137   }
138 };
139
140
141 class edge {
142   // edge: acts as a passageway between elements; thus, we can lock
143   // passageways to prevent multiple refinement paths from accessing a
144   // single element simultaneously.  
145
146   // * Only ONE edge exists between two REAL nodes. 
147   // * An edge has indices of its two endpoints on the chunk on which the 
148   // edge exists. 
149   // * An edge has references to the two elements that on either side of it. 
150   // * Also, it has a reference to itself (myRef).
151   // in addition, it has a pointer to the local chunk (C) on which it resides.
152   // * Finally, it has a lock which can cut off access to the elements on 
153   // either side of the edge. 
154   edgeRef myRef;
155   chunk *C;
156   int theLock;  // 0 if open; 1 if locked
157
158  public:
159   elemRef elements[2];
160
161   // Basic edge contructor: unlocked by default
162   edge() { 
163     for (int i=0; i<2; i++)
164       elements[i].init();
165     myRef.init();  C = NULL;  theLock = 0; 
166   }
167   
168   void sanityCheck(chunk *C,edgeRef ref);
169   // Initializers for an edge: either set all fields, or just the
170   // index and chunk info.  Element references can be set or
171   // modified later with updateElement
172   void init() { theLock = 0; } // equivalent to constructor initializations
173   void init(int i, chunk *cPtr);
174   void init(elemRef *e, int i, chunk *cPtr);
175
176   // updateElement: Set or modify the element references of the edge.  
177   // If the edge has not been initialized, passing
178   // a -1 or the null reference as oldval will initialize any uninitialized
179   // field in the edge
180   void updateElement(elemRef oldval, elemRef newval);
181   const edgeRef &getRef() const { return myRef; } // get reference to this edge
182
183   // getNbrRef: called by an element on one of it's edges to get a
184   // reference to the neighbor on the other side of the edge; the
185   // element passes in its own reference for comparison
186   const elemRef &getNbrRef(const elemRef &er) { 
187     if (er == elements[0])
188       return elements[1];
189     else if (!(er == elements[1]))
190       CkAbort("ERROR: edge::getNbrRef: input edgeRef not on edge\n");
191     return elements[0];
192   }
193   // lock, unlock, and locked control access to the edge's lock
194   void lock() { theLock = 1; }
195   void unlock() { theLock = 0; }
196   int locked() { return theLock; }
197 };
198
199 class element {
200   // element: triangular elements defined by three node references,
201   // and having three edge references to control refinement access to
202   // this element, and provide connectivity to adjacent elements
203
204   // targetArea is the area to attempt to refine this element
205   // below. It's unset state is -1.0, and this will be detected as a
206   // 'no refinement desired' state. currentArea is the actual area of
207   // the triangle that was cached most recently.  This gets updated by
208   // calls to either getArea or calculateArea.
209   double targetArea, currentArea;
210
211   // When refining this element, it may discover that an adjacent
212   // element needs to be refined first. This element then sets its
213   // depend field to indicate that it is dependent on another element's
214   // refinement.  It also tells the adjacent element that it is
215   // dependent, and sends its reference to the adjacent element.  The
216   // adjacent element sets its dependent field to that reference. THUS:
217   // depend=1 if this element depends on another to refine first, 0 otherwise
218   int depend;
219   // dependent is a valid reference if some adjacent element depends
220   // on this one to refine first; dependent = {-1, -1} (nullRef) otherwise
221   elemRef dependent;
222
223   // elements on different chunks use the following data fields to
224   // negotiate and communicate refinements; see below for more details
225   int specialRequest, pendingRequest, requestResponse;
226   elemRef specialRequester;
227   node newNode, otherNode;
228   edgeRef newLongEdgeRef;
229   
230   // the reference for this edge and its chunk
231   elemRef myRef;
232   chunk *C;
233
234  public:
235   /* node and edge numberings follow this convention
236                          0
237                         / \
238                       2/   \0
239                       /     \
240                      2_______1
241                          1
242   */
243   int nodes[3];
244   edgeRef edges[3];
245   
246   element(); // Basic element constructor
247   
248   void sanityCheck(chunk *C,elemRef ref);
249     
250   // Initializers: specifying edges is optional; they can be added
251   // with updateEdges later on
252   void init(); // equivalent to constructor initializations
253   void init(int *n, edgeRef *e, int index, chunk *chk);
254   void init(int *n, int index, chunk *chk);
255   void updateEdge(int idx, edgeRef e) { edges[idx] = e; }
256   void updateEdges(edgeRef e0, edgeRef e1, edgeRef e2);
257   
258   // Access a particular node or edge reference
259   node getNode(int i) const;
260   const edgeRef &getEdge(int i) const { return edges[i]; }
261   // getOpNode returns the node index for node opposite edge[e]
262   int getOpNode(int e) { return (e+2)%3; }
263   // getOtherNode returns the node index for one endpoint on edge[e]
264   int getOtherNode(int e) { return e; }
265   // getNeighbor returns the neighboring element reference along the edge[e].
266   elemRef getNeighbor(int e) const;
267   
268   // getArea & calculateArea both set currentArea; getArea returns it
269   // setTargetArea initializes or minimizes targetArea
270   // getTargetArea & getCachedArea provide access to targetArea & currentArea
271   double getArea();
272   void calculateArea();
273   void setTargetArea(double area) {
274     if (((targetArea > area) || (targetArea < 0.0))  &&  (area >= 0.0))
275       targetArea = area;
276   }
277   double getTargetArea() { return targetArea; }
278   double getCachedArea() { return currentArea; }
279
280   // These methods manipulate this element's dependence on another
281   void setDependency() { depend = 1; }
282   void unsetDependency() { depend = 0; }
283   int hasDependency() { return (depend); }
284
285   // These methods manipulate what element is dependent on this one
286   void setDependent(elemRef e) { dependent = e; }
287   void setDependent(int cId, int i) { dependent.cid = cId; dependent.idx = i; }
288   void unsetDependent() { dependent.idx = dependent.cid = -1; }
289   int hasDependent() { return ((dependent.idx!=-1) && (dependent.cid!=-1)); }
290   void tellDepend() {
291     if (hasDependent()) {
292       dependent.unsetDependency(C);
293       unsetDependent();
294     }
295   }
296         
297   // These methods handle the negotiation of refinement between two
298   // elements on different chunks.  In addition, once the refinement
299   // relationship has been established, they control how information
300   // on that refinement is passed between the two elements.  Elements
301   // A needs to refine along with element B, but they are on different
302   // chunks. Element A makes a special request of B to refine first. A
303   // labels itself as 'pendingRequest' to avoid resending the request.
304   // B receives the request and performs a test based on local id and
305   // chunk id.  The greater index gets to refine first.  If B fails,
306   // it sends a special request to A, otherwise it dos half of the
307   // refinement sends the refinement info to A in the form of a
308   // request response.  When A receives the response, it completes its
309   // half of the refinement. Note that both A & B could have sent
310   // special requests simultaneously.  The test resolves which one
311   // gets to proceed first, and when one element is preparing a
312   // request response, it ignores incoming special requests.
313   void setSpecialRequest(elemRef r) { specialRequest=1; specialRequester=r; }
314   int isSpecialRequest() { return (specialRequest == 1); }
315   int isPendingRequest() { return (pendingRequest == 1); }
316   int isRequestResponse() { return (requestResponse == 1); }
317   void setRequestResponse(node n, node o, edgeRef e) { 
318     requestResponse = 1; 
319     newNode = n;
320     otherNode = o;
321     newLongEdgeRef = e;
322   }
323
324   // These methods handle various types and stages of refinements.
325   //
326   // refine is the method by which a refinement on an element is
327   // initiated. refine first checks if there are any requests arrived
328   // or pending, and handles these appropriately.  If there are none,
329   // the longest edge is determined, and some tests are performed to
330   // determine:
331   // 1. if the edge is on a border, splitBorder is called
332   // 2. if the edge is also the neighbor's longest, splitNeighbors is called
333   // 3. if the edge is not the neighbor's longest, refineNeighbor is called
334   void refine();
335   // refineNeighbor has to tell this element's neighbor to refine, and
336   // tell the neighbor that this element depends on it for its own
337   // refinement.  In addition, this element notes that it is dependent
338   // on a neighbor, and does not attempt further refinement until it
339   // hears from the neighbor
340   void refineNeighbor(int longEdge);
341   // splitBorder and splitNeighbors set up a locked perimeter of edges
342   // around the area of refinement and call splitBorderLocal and
343   // splitNeighborsLocal respectively.  splitNeighbors handles the
344   // sending of special requests for refinement should the neighbor
345   // element be located on a different chunk.  In this case, locking
346   // the perimeter is put off until the refinement actually happens.
347   void splitBorder(int longEdge);
348   void splitNeighbors(int longEdge);
349   // These methods handle local refinement -- they create a new node,
350   // new edge(s), new element(s), and make sure everything is properly
351   // connected.
352   void splitBorderLocal(int longEdge, int opnode, int othernode, int modEdge);
353   void splitNeighborsLocal(int longEdge, int opnode, int othernode, 
354                            int modEdge, int nbrLongEdge, int nbrOpnode,
355                            int nbrOthernode, int nbrModEdge, const elemRef &nbr);
356   
357   // These methods handle the two-phase split of neighboring elements
358   // on differing chunks.  splitResponse is called by the element that
359   // has precedence, while splitHelp is called by the secondary
360   // element to complete the refinement.  Note that an element that
361   // receives a special request may choose to simply accept it and
362   // continue with the splitResponse regardless of precedence when it
363   // is not scheduled for any refinements.  This is because the status
364   // of the elements are only checked if the elements need to be
365   // refined.
366   void splitHelp(int longEdge);
367   void splitResponse(int longEdge);
368
369   // These are helper methods.
370   //
371   // Finds the element's longes edge and returns its index in edges.
372   int findLongestEdge();
373   // Checks the status of neighbor on longEdge: returns -1 if there is
374   // no neighbor (border case), 1 if there is a neighbor and that
375   // neighbor has the same longEdge, and 0 otherwise (neighbor does
376   // have longEdge as its longest edge)
377   int checkNeighbor(int longEdge);
378   // Checks if e is the edge reference for the longest edge of this element
379   int checkIfLongEdge(edgeRef e);
380 };
381
382 /**
383  * The user inherits from this class to receive "split" calls,
384  * and to be informed when the refinement is complete.
385  */
386 class refineClient {
387 public:
388   virtual ~refineClient() {}
389
390   /**
391    * This triangle of our chunk is being split along this edge.
392    *
393    * For our purposes, edges are numbered 0 (connecting nodes 0 and 1), 
394    * 1 (connecting 1 and 2), and 2 (connecting 2 and 0).
395    * 
396    * Taking as A and B the (triangle-order) nodes of the splitting edge:
397    *
398    *                     C                      C                 
399    *                    / \                    /|\                  
400    *                   /   \                  / | \                 
401    *                  /     \      =>        /  |  \                
402    *                 /       \              /   |   \               
403    *                /         \            /old | new\            
404    *               B --------- A          B --- D --- A         
405    *
406    *   The original triangle's node A should be replaced by D;
407    * while a new triangle should be inserted with nodes CAD.
408    *
409    *   The new node D's location should equal A*(1-frac)+B*frac.
410    * For a simple splitter, frac will always be 0.5.
411    *
412    *   If nodes A and B are shared with some other processor,
413    * that processor will also receive a "split" call for the
414    * same edge.  If nodes A and B are shared by some other local
415    * triangle, that triangle will immediately receive a "split" call
416    * for the same edge.  
417          *
418          * flag denotes the properties of the new node added by the split       
419    * 0x1 - node is on the chunk boundary
420          * 0x2 - since split will be called twice for each new node,
421          *       this bit shows whether it is the first time or not
422          *
423    * Client's responsibilities:
424    *   -Add the new node D.  Since both sides of a shared local edge
425    *      will receive a "split" call, you must ensure the node is
426    *      not added twice.
427    *   -Update connectivity for source triangle
428    *   -Add new triangle.
429    */
430   virtual void split(int triNo,int edgeOfTri,int movingNode,double frac) =0;
431   virtual void split(int triNo,int edgeOfTri,int movingNode,double frac,int flags) =0;
432
433 };
434
435 class refineResults; //Used by refinement API to store intermediate results
436
437 // ---------------------------- Chare Arrays -------------------------------
438 class chunk : public TCharmClient1D {
439   // Data fields for this chunk's array index, and counts of elements,
440   // edges, and nodes located on this chunk; numGhosts is numElements
441   // plus number of ghost elements surrounding this chunk
442
443   // current sizes of arrays allocated for the mesh
444   int sizeElements, sizeEdges, sizeNodes;
445   
446   // debug_counter is used to print successive snapshots of the chunk
447   // and match them up to other chunk snapshots; refineInProgress
448   // flags that the refinement loop is active; modified flags that a
449   // target area for some element on this chunk has been modified
450   int debug_counter, refineInProgress, modified;
451
452   // meshLock is used to lock the mesh for expansion; if meshlock is
453   // zero, the mesh can be either accessed or locked; accesses to the
454   // mesh (by a chunk method) decrement the lock, and when the
455   // accesses are complete, the lock is incremented; when an expansion
456   // of the mesh is required, the meshExpandFlag is set, indicating
457   // that no more accesses will be allowed to the mesh until the
458   // adjuster gets control and completes the expansion; when the
459   // adjuster gets control, it sets meshLock to 1 and when it is
460   // finished, it resets both variables to zero.  See methods below.
461   int meshLock, meshExpandFlag;
462
463   // private helper methods used by FEM interface functions
464   void deriveNodes();
465   int edgeLocal(elemRef e1, elemRef e2);
466   int findEdge(int n1, int n2);
467   int addNewEdge();
468   int getNbrRefOnEdge(int n1, int n2, int *conn, int nGhost, int *gid, 
469                       int idx, elemRef *er);
470   int hasEdge(int n1, int n2, int *conn, int idx);
471   
472  public:
473   refineResults *refineResultsStorage;
474  
475   // the chunk's components, left public for sanity's sake
476   int cid;
477   int numElements, numEdges, numNodes, numGhosts, numChunks;
478   std::vector<element> theElements;
479   std::vector<edge> theEdges;
480   std::vector<node> theNodes;
481
482   // client to report refinement split information to
483   refineClient *theClient;
484
485   // Basic constructor
486   chunk(chunkMsg *);
487   chunk(CkMigrateMessage *m) : TCharmClient1D(m) { };
488   
489   void sanityCheck(void);
490   
491   void setupThreadPrivate(CthThread forThread) {
492     CtvAccessOther(forThread, _refineChunk) = this;
493   }
494   // entry methods
495   
496   // Initiates a refinement for a single element
497   void refineElement(int i, double area);
498   // Loops through all elements performing refinements as needed
499   void refiningElements();
500
501   // The following methods simply provide remote access to local data
502   // See above for details of each
503   void updateElement(int i, objRef oldval, objRef newval);
504   void specialRequest(int reqestee, elemRef requester);
505   void specialRequestResponse(int i, double newNodeX, double newNodeY, 
506                               double otherNodeX, double otherNodeY, 
507                               edgeRef newLongEdgeRef);
508   doubleMsg *getArea(int i);
509   intMsg *lock(int i);
510   void unlock(int i);
511   intMsg *locked(int i);
512   intMsg *checkElement(objRef oR, int i);
513   refMsg *getNeighbor(objRef oR, int i);
514   void setTargetArea(int i, double area);
515   void updateEdges(int i, edgeRef e0, edgeRef e1, edgeRef e2);
516   void unsetDependency(int i);
517   void setDependent(objRef oR, int i);
518   intMsg *hasDependent(int i);
519
520   // meshLock methods
521   void accessLock();  // waits until meshExpandFlag not set, then decs meshLock
522   void releaseLock(); // incs meshLock
523   void adjustFlag();  // sets meshExpandFlag
524   void adjustLock();  // waits until meshLock is 0, then sets it to 1
525   void adjustRelease();  // resets meshLock and meshExpandFlag to 0
526
527   // used to print snapshots of all chunks at once (more or less)
528   void print();
529
530   // local methods
531
532   // These methods are part of the interface with the FEM framework:
533
534   // Sets the node coordinates and recalculates element areas
535   void updateNodeCoords(int nNode, double *coord, int nEl);
536   // multipleRefine sets target areas of elements as specified by
537   // desiredArea, and starts refining.  Each split that occurs is
538   // reported to client, and the methods returns only when all
539   // refinement is completed.
540   void multipleRefine(double *desiredArea, refineClient *client);
541   void newMesh(int nEl, int nGhost,const int *conn_,const int *gid_, int idxOffset);
542   void addRemoteEdge(int elem, int localEdge, edgeRef er);
543
544   
545   // These access and set local flags
546   void setModified() { modified = 1; }
547   int isModified() { return modified; }
548   void setRefining() { refineInProgress = 1; }
549   int isRefining() { return refineInProgress; }
550
551   // these methods allow for run-time additions/modifications to the chunk
552   void allocMesh(int nEl);
553   void adjustMesh();
554   int addNode(node n);
555   edgeRef addEdge();
556   elemRef addElement(int n1, int n2, int n3);
557   elemRef addElement(int n1, int n2, int n3,
558                      edgeRef er1, edgeRef er2, edgeRef er3);
559
560   void debug_print(int c);  // prints a snapshot of the chunk to file
561   void out_print();  // prints a readable meshfile that can be used as input
562 };