source: XIOS/dev/dev_olga/extern/remap/src/elt.hpp @ 1158

Last change on this file since 1158 was 1158, checked in by oabramkina, 4 years ago

Two server levels: merging with trunk r1137.
There are bugs.

File size: 4.1 KB
Line 
1#ifndef __ELT_H__
2#define __ELT_H__
3#include <list>
4#include "triple.hpp"
5
6#define NMAX 10 /**< maximum number of vertices for polygons */
7
8#define NOT_FOUND -1
9
10namespace sphereRemap {
11
12using namespace std;
13
14Coord barycentre(const Coord *x, int n);
15
16/** Two great or small circles (or mixed) have two or zero intersections.
17    The coordinates of the intersections are stored in `pt` and `nb` holds the number of intersections (0 or 2).
18    */
19struct Ipt
20{
21        int nb;
22        Coord pt[2];
23};
24
25struct Sgm // edge
26{
27        Coord n; // normal
28        Coord xt[2]; // endpoints
29        double d; // (see Elt)
30};
31
32struct GloId
33{
34        int rank;
35        int ind; /* local id */
36  long globalId ;
37
38        bool operator<(const GloId& other) const {
39                return (rank == other.rank) ? (ind < other.ind) : (rank < other.rank);
40        }
41};
42
43struct Polyg
44{ 
45        /* Note:  for the target grid elements the id (rank and local id) depends on the order of the target grid elements as read from the nc-file whereas for source grid elements it depends on the SS-tree (i.e. super mesh distribution, not the order in the nc-file) */
46        struct GloId id;
47        struct GloId src_id;
48        int n; /* number of vertices */
49        double area;
50        Coord x; /* barycentre */
51};
52
53struct Elt : Polyg
54{
55        Elt() {}
56        Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert)
57        {
58                int k = 0;
59                vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]);
60                for (int i = 1; i < max_num_vert; i++)
61                {
62                        vertex[k] = xyz(bounds_lon[i], bounds_lat[i]);
63                        /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */
64                        if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS) 
65                                break;
66                        /* eliminate zero edges: move to next vertex only if it is different */
67                        if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS)
68                                k++;
69                        else
70                                /* cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl */ ;
71                }
72                n = k;
73                x = barycentre(vertex, n);
74        }
75
76        Elt& operator=(const Elt& rhs)
77        {
78                id    = rhs.id;
79                src_id    = rhs.src_id;
80                n    = rhs.n;
81                area = rhs.area;
82                x    = rhs.x;
83                val  = rhs.val;
84                grad = rhs.grad;
85                is   = rhs.is;
86
87                for(int i = 0; i < NMAX; i++)
88                {
89                        neighbour[i] = rhs.neighbour[i];
90                        d[i]         = rhs.d[i];
91                        edge[i]      = rhs.edge[i];
92                        vertex[i]    = rhs.vertex[i];
93                        gradNeigh[i] = rhs.gradNeigh[i];
94                }
95                return *this;
96        }
97
98        void delete_intersections()
99        {
100                for (list<Polyg*>::iterator it = this->is.begin(); it != this->is.end(); it++)
101                {
102                        Polyg* poly = *it;
103                        delete poly;
104                }
105        }
106
107  void insert_vertex(int i, const Coord& v)
108  {
109    for(int j=n; j > i ; j--)
110    {
111      vertex[j]=vertex[j-1] ;
112      edge[j]=edge[j-1] ;
113      d[j]=d[j-1] ;
114      neighbour[j]=neighbour[j-1] ;
115    }
116    vertex[i+1]=v ;
117    n++ ;
118  }
119 
120        int neighbour[NMAX];
121        double d[NMAX]; /**< distance of centre of small circle to origin, zero if great circle */
122        double val;     /**< value (sample if src element, interpolated if dest element) */
123        Coord vertex[NMAX];
124        Coord edge[NMAX]; /**< edge normals: if great circle tangential to sphere, if small circle parallel to pole */
125        Coord grad;    /**< gradient of the reconstructed linear function over this element */
126        Coord gradNeigh[NMAX]; /**< for weight computation: gradients for val=1 at individual neighbours */
127        struct GloId neighId[NMAX]; /**< weight computation needs to know global IDs for all sources with "link" */
128        std::list<Polyg*> is; /**< intersections */
129};
130
131static double normals(Elt &elt, const Coord &pole)
132{
133        double nmin = 17.;
134        for (int i = 0; i < elt.n; i++)  // supposed oriented
135        {
136                int j = (i+1) % elt.n;
137                elt.edge[i] = crossprod(elt.vertex[j], elt.vertex[i]); 
138                Coord t = elt.vertex[j] - elt.vertex[i];
139                /* polygonal grid || vertices not on same latitude */
140                if (pole == ORIGIN || fabs(scalarprod(t, pole)) > EPS)  // great circle
141                {
142                        double n = norm(elt.edge[i]);
143                        //assert(n > 0);
144                        if (n < nmin) nmin = n;
145                        elt.edge[i] = proj(elt.edge[i]);
146                        elt.d[i] = 0.0;
147                }
148                else /* lan lot grid && vertices on same latitude => small circle */
149                {
150                        int north = (scalarprod(elt.edge[i], pole) < 0) ? -1 : 1;
151                        elt.edge[i] = pole * north;
152                        elt.d[i] = scalarprod(elt.vertex[i], elt.edge[i]);
153                }
154        }
155        return nmin;
156}
157
158}
159
160#endif
Note: See TracBrowser for help on using the repository browser.