- Timestamp:
- 11/09/18 09:44:17 (5 years ago)
- Location:
- XIOS/dev/dev_olga
- Files:
-
- 2 added
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/dev_olga/extern/remap/src/intersection_ym.cpp
r849 r1590 8 8 #include <stdlib.h> 9 9 #include <limits> 10 #include <array> 11 #include <cstdint> 12 #include "earcut.hpp" 13 #include <fstream> 14 10 15 11 16 #define epsilon 1e-3 // epsilon distance ratio over side lenght for approximate small circle by great circle … … 16 21 using namespace std; 17 22 using namespace ClipperLib ; 23 18 24 19 25 double intersect_ym(Elt *a, Elt *b) 20 26 { 27 28 using N = uint32_t; 29 using Point = array<double, 2>; 30 vector<Point> vect_points; 31 vector< vector<Point> > polyline; 21 32 22 33 // transform small circle into piece of great circle if necessary … … 45 56 double cos_alpha; 46 57 58 /// vector<p2t::Point*> polyline; 47 59 for(int n=0; n<na;n++) 48 60 { … … 51 63 a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 52 64 a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 53 } 65 66 vect_points.push_back( array<double, 2>() ); 67 vect_points[n][0] = a_gno[n].x; 68 vect_points[n][1] = a_gno[n].y; 69 70 } 71 72 polyline.push_back(vect_points); 73 vector<N> indices_a_gno = mapbox::earcut<N>(polyline); 74 75 double area_a_gno=0 ; 76 for(int i=0;i<indices_a_gno.size()/3;++i) 77 { 78 Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ; 79 Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ; 80 Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ; 81 area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 82 } 83 84 vect_points.clear(); 85 polyline.clear(); 86 indices_a_gno.clear(); 87 88 54 89 55 90 for(int n=0; n<nb;n++) … … 59 94 b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ; 60 95 b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1 61 } 62 96 97 vect_points.push_back( array<double, 2>() ); 98 vect_points[n][0] = b_gno[n].x; 99 vect_points[n][1] = b_gno[n].y; 100 } 101 102 103 polyline.push_back(vect_points); 104 vector<N> indices_b_gno = mapbox::earcut<N>(polyline); 105 106 double area_b_gno=0 ; 107 for(int i=0;i<indices_b_gno.size()/3;++i) 108 { 109 Coord x0 = Ox * polyline[0][indices_b_gno[3*i]][0] + Oy* polyline[0][indices_b_gno[3*i]][1] + Oz ; 110 Coord x1 = Ox * polyline[0][indices_b_gno[3*i+1]][0] + Oy* polyline[0][indices_b_gno[3*i+1]][1] + Oz ; 111 Coord x2 = Ox * polyline[0][indices_b_gno[3*i+2]][0] + Oy* polyline[0][indices_b_gno[3*i+2]][1] + Oz ; 112 area_b_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 113 } 114 115 vect_points.clear(); 116 polyline.clear(); 117 indices_b_gno.clear(); 63 118 64 119 … … 109 164 clip.AddPaths(dst, ptClip, true); 110 165 clip.Execute(ctIntersection, intersection); 111 166 112 167 double area=0 ; 113 168 114 169 for(int ni=0;ni<intersection.size(); ni++) 115 170 { 116 // go back into real coordinate on the sphere117 171 Coord* intersectPolygon=new Coord[intersection[ni].size()] ; 118 172 for(int n=0; n < intersection[ni].size(); n++) 119 173 { 120 double x=intersection[ni][n].X/xscale+xoffset ; 121 double y=intersection[ni][n].Y/yscale+yoffset ; 122 123 intersectPolygon[n]=Ox*x+Oy*y+Oz ; 124 intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 125 } 126 127 // remove redondants vertex 128 int nv=0 ; 174 intersectPolygon[n].x=intersection[ni][n].X/xscale+xoffset ; 175 intersectPolygon[n].y=intersection[ni][n].Y/yscale+yoffset ; 176 } 177 178 179 int nv=0; 180 129 181 for(int n=0; n < intersection[ni].size(); n++) 130 182 { 131 if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex) 183 double dx=intersectPolygon[n].x-intersectPolygon[(n+1)%intersection[ni].size()].x ; 184 double dy=intersectPolygon[n].y-intersectPolygon[(n+1)%intersection[ni].size()].y ; 185 186 if (dx*dx+dy*dy>fusion_vertex*fusion_vertex) 187 { 188 intersectPolygon[nv]=intersectPolygon[n] ; 189 190 vect_points.push_back( array<double, 2>() ); 191 vect_points[nv][0] = intersectPolygon[n].x; 192 vect_points[nv][1] = intersectPolygon[n].y; 193 194 nv++ ; 195 } 196 197 198 } 199 200 polyline.push_back(vect_points); 201 vect_points.clear(); 202 203 if (nv>2) 204 { 205 206 vector<N> indices = mapbox::earcut<N>(polyline); 207 208 double area2=0 ; 209 for(int i=0;i<indices.size()/3;++i) 132 210 { 133 intersectPolygon[nv]=intersectPolygon[n] ; 134 nv++ ; 211 Coord x0 = Ox * polyline[0][indices[3*i]][0] + Oy* polyline[0][indices[3*i]][1] + Oz ; 212 Coord x1 = Ox * polyline[0][indices[3*i+1]][0] + Oy* polyline[0][indices[3*i+1]][1] + Oz ; 213 Coord x2 = Ox * polyline[0][indices[3*i+2]][0] + Oy* polyline[0][indices[3*i+2]][1] + Oz ; 214 area2+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 135 215 } 136 } 137 138 139 if (nv>2) 140 { 216 217 polyline.clear(); 218 219 for(int n=0; n < nv; n++) 220 { 221 intersectPolygon[n] = Ox*intersectPolygon[n].x+Oy*intersectPolygon[n].y+Oz; 222 intersectPolygon[n] = intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 223 } 224 225 141 226 // assign intersection to source and destination polygons 142 227 Polyg *is = new Polyg; 143 228 is->x = exact_barycentre(intersectPolygon,nv); 144 is->area = polygonarea(intersectPolygon,nv) ; 229 // is->area = polygonarea(intersectPolygon,nv) ; 230 is->area = area2 ; 231 145 232 // if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ; 146 233 if (is->area==0.) delete is ; … … 161 248 delete[] b_gno ; 162 249 return area ; 250 163 251 } 252 253 164 254 165 255 void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates) -
XIOS/dev/dev_olga/extern/remap/src/meshutil.cpp
r877 r1590 133 133 { 134 134 for (int i = 0; i < elts[j]->n; i++) 135 if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; 135 if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; // no neighbour 136 else if (elts[elts[j]->neighbour[i]]->is.size() == 0) neighbours[i]=NULL ; // neighbour has none supermesh cell 136 137 else neighbours[i] = elts[elts[j]->neighbour[i]]; 137 138 -
XIOS/dev/dev_olga/extern/remap/src/polyg.cpp
r950 r1590 20 20 void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g) 21 21 { 22 23 24 25 26 27 28 29 30 31 32 33 34 35 22 Coord ga = vertex[0] - g; 23 Coord gb = vertex[1] - g; 24 Coord vertical = crossprod(ga, gb); 25 if (N > 2 && scalarprod(g, vertical) < 0) // (GAxGB).G 26 { 27 for (int i = 0; i < N/2; i++) 28 swap(vertex[i], vertex[N-1-i]); 29 30 for (int i = 0; i < (N-1)/2; i++) 31 { 32 swap(edge[N-2-i], edge[i]); 33 swap(d[i], d[N-2-i]); 34 } 35 } 36 36 } 37 37 … … 39 39 void normals(Coord *x, int n, Coord *a) 40 40 { 41 42 41 for (int i = 0; i<n; i++) 42 a[i] = crossprod(x[(i+1)%n], x[i]); 43 43 } 44 44 45 45 Coord barycentre(const Coord *x, int n) 46 46 { 47 48 49 50 51 52 53 54 55 56 57 47 if (n == 0) return ORIGIN; 48 Coord bc = ORIGIN; 49 for (int i = 0; i < n; i++) 50 bc = bc + x[i]; 51 /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 52 which can occur when weighted with tiny area */ 53 54 assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 55 //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 56 57 return proj(bc); 58 58 } 59 59 … … 62 62 static Coord tetrah_side_diff_centre(Coord a, Coord b) 63 63 { 64 64 Coord n = crossprod(a,b); 65 65 double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z; 66 66 assert(sinc2 < 1.0 + EPS); 67 67 68 68 // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab 69 69 // approx: 70 70 // double u = sinc2/6. + (3./40.)*sinc2*sinc2; 71 71 72 73 74 75 72 // exact 73 if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */ 74 return n * (M_PI_2 - 1); 75 double sinc = sqrt(sinc2); 76 76 double u = asin(sinc)/sinc - 1; 77 77 … … 82 82 Coord gc_normalintegral(const Coord *x, int n) 83 83 { 84 85 86 87 88 84 Coord m = barycentre(x, n); 85 Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]); 86 for (int i = 1; i < n; i++) 87 bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]); 88 return bc*0.5; 89 89 } 90 90 91 91 Coord exact_barycentre(const Coord *x, int n) 92 92 { 93 94 95 96 97 98 99 93 if (n >= 3) 94 { 95 return proj(gc_normalintegral(x, n)); 96 } 97 else if (n == 0) return ORIGIN; 98 else if (n == 2) return midpoint(x[0], x[1]); 99 else if (n == 1) return x[0]; 100 100 } 101 101 102 102 Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole) 103 103 { 104 double hemisphere = (a.z > 0) ? 1: -1; 105 106 double lat = hemisphere * (M_PI_2 - acos(a.z)); 107 double lon1 = atan2(a.y, a.x); 108 double lon2 = atan2(b.y, b.x); 109 double lon_diff = lon2 - lon1; 110 111 // wraparound at lon=-pi=pi 112 if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 113 else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 114 115 // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 116 Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 117 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 118 hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 119 Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 120 Coord t[] = {a, b, p}; 121 if (hemisphere < 0) swap(t[0], t[1]); 122 return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 123 } 124 125 126 double triarea(Coord& A, Coord& B, Coord& C) 127 { 128 double a = ds(B, C); 129 double b = ds(C, A); 130 double c = ds(A, B); 131 double s = 0.5 * (a + b + c); 132 double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c)); 133 // assert(t >= 0); 134 if (t<1e-20) return 0. ; 135 return 4 * atan(sqrt(t)); 104 double hemisphere = (a.z > 0) ? 1: -1; 105 106 double lat = hemisphere * (M_PI_2 - acos(a.z)); 107 double lon1 = atan2(a.y, a.x); 108 double lon2 = atan2(b.y, b.x); 109 double lon_diff = lon2 - lon1; 110 111 // wraparound at lon=-pi=pi 112 if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 113 else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 114 115 // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 116 Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 117 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 118 hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 119 Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 120 Coord t[] = {a, b, p}; 121 if (hemisphere < 0) swap(t[0], t[1]); 122 return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 123 } 124 125 126 double triarea(const Coord& A, const Coord& B, const Coord& C) 127 { 128 double a = ds(B, C); 129 double b = ds(C, A); 130 double c = ds(A, B); 131 double tmp ; 132 133 if (a<b) { tmp=a ; a=b ; b=tmp ; } 134 if (c > a) { tmp=a ; a=c ; c=b, b=tmp; } 135 else if (c > b) { tmp=c ; c=b ; b=tmp ; } 136 137 double s = 0.5 * (a + b + c); 138 double t = tan(0.25*(a+(b+c))) * tan(0.25*(c-(a-b))) * tan(0.25*( c + (a-b) )) * tan(0.25*( a + (b - c))); 139 if (t>0) return 4 * atan(sqrt(t)); 140 else 141 { 142 std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ; 143 return 0 ; 144 } 136 145 } 137 146 … … 142 151 double alun(double b, double d) 143 152 { 144 145 146 147 148 149 150 153 double a = acos(d); 154 assert(b <= 2 * a); 155 double s = a + 0.5 * b; 156 double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); 157 double r = sqrt(1 - d*d); 158 double p = 2 * asin(sin(0.5*b) / r); 159 return p*(1 - d) - 4*atan(sqrt(t)); 151 160 } 152 161 … … 160 169 double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg) 161 170 { 162 163 164 165 166 167 168 169 170 171 172 173 174 171 if (N < 3) 172 return 0; /* polygons with less then three vertices have zero area */ 173 Coord t[3]; 174 t[0] = barycentre(x, N); 175 Coord *g = new Coord[N]; 176 double area = 0; 177 Coord gg_exact = gc_normalintegral(x, N); 178 for (int i = 0; i < N; i++) 179 { 180 /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ 181 int ii = (i + 1) % N; 182 t[1] = x[i]; 183 t[2] = x[ii]; 175 184 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 185 assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 186 double area_gc = triarea(t[0], t[1], t[2]); 187 double area_sc_gc_moon = 0; 188 if (d[i]) /* handle small circle case */ 189 { 190 Coord m = midpoint(t[1], t[2]); 191 double mext = scalarprod(m, c[i]) - d[i]; 192 char sgl = (mext > 0) ? -1 : 1; 193 area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 194 gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 195 } 196 area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 197 g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 198 } 199 gg = barycentre(g, N); 200 gg_exact = proj(gg_exact); 201 delete[] g; 202 gg = gg_exact; 203 return area; 195 204 } 196 205 197 206 double polygonarea(Coord *vertices, int N) 198 207 { 199 200 201 202 203 204 205 206 208 assert(N >= 3); 209 210 /* compute polygon area as sum of triangles */ 211 Coord centre = barycentre(vertices, N); 212 double area = 0; 213 for (int i = 0; i < N; i++) 214 area += triarea(centre, vertices[i], vertices[(i+1)%N]); 215 return area; 207 216 } 208 217 209 218 int packedPolygonSize(const Elt& e) 210 219 { 211 212 220 return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + 221 sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); 213 222 } 214 223 215 224 void packPolygon(const Elt& e, char *buffer, int& pos) 216 225 { 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 226 *((GloId *) &(buffer[pos])) = e.id; 227 pos += sizeof(e.id); 228 *((GloId *) &(buffer[pos])) = e.src_id; 229 pos += sizeof(e.src_id); 230 231 *((Coord *) &(buffer[pos])) = e.x; 232 pos += sizeof(e.x); 233 234 *((double*) &(buffer[pos])) = e.val; 235 pos += sizeof(e.val); 236 237 *((int *) & (buffer[pos])) = e.n; 238 pos += sizeof(e.n); 239 240 for (int i = 0; i < e.n; i++) 241 { 242 *((double *) & (buffer[pos])) = e.d[i]; 243 pos += sizeof(e.d[i]); 244 245 *((Coord *) &(buffer[pos])) = e.vertex[i]; 246 pos += sizeof(e.vertex[i]); 247 } 239 248 240 249 } … … 242 251 void unpackPolygon(Elt& e, const char *buffer, int& pos) 243 252 { 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 253 e.id = *((GloId *) &(buffer[pos])); 254 pos += sizeof(e.id); 255 e.src_id = *((GloId *) &(buffer[pos])); 256 pos += sizeof(e.src_id); 257 258 e.x = *((Coord *) & (buffer[pos]) ); 259 pos += sizeof(e.x); 260 261 e.val = *((double *) & (buffer[pos])); 262 pos += sizeof(double); 263 264 e.n = *((int *) & (buffer[pos])); 265 pos += sizeof(int); 266 267 for (int i = 0; i < e.n; i++) 268 { 269 e.d[i] = *((double *) & (buffer[pos])); 270 pos += sizeof(double); 271 272 e.vertex[i] = *((Coord *) & (buffer[pos])); 273 pos += sizeof(Coord); 274 } 266 275 } 267 276 268 277 int packIntersectionSize(const Elt& elt) 269 278 { 270 279 return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); 271 280 } 272 281 … … 274 283 { 275 284 for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) 276 277 278 279 280 285 { 286 *((int *) &(buffer[0])) += 1; 287 288 *((int *) &(buffer[pos])) = e.id.ind; 289 pos += sizeof(int); 281 290 282 291 *((double *) &(buffer[pos])) = e.area; 283 292 pos += sizeof(double); 284 293 285 286 294 *((GloId *) &(buffer[pos])) = (*it)->id; 295 pos += sizeof(GloId); 287 296 288 289 290 291 292 293 294 295 297 *((int *) &(buffer[pos])) = (*it)->n; 298 pos += sizeof(int); 299 *((double *) &(buffer[pos])) = (*it)->area; 300 pos += sizeof(double); 301 302 *((Coord *) &(buffer[pos])) = (*it)->x; 303 pos += sizeof(Coord); 304 } 296 305 } 297 306 298 307 void unpackIntersection(Elt* e, const char* buffer) 299 308 { 300 301 309 int ind; 310 int pos = 0; 302 311 303 304 305 306 307 308 309 310 312 int n = *((int *) & (buffer[pos])); 313 pos += sizeof(int); 314 for (int i = 0; i < n; i++) 315 { 316 ind = *((int*) &(buffer[pos])); 317 pos+=sizeof(int); 318 319 Elt& elt= e[ind]; 311 320 312 321 elt.area=*((double *) & (buffer[pos])); 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 } 332 333 } 322 pos += sizeof(double); 323 324 Polyg *polygon = new Polyg; 325 326 polygon->id = *((GloId *) & (buffer[pos])); 327 pos += sizeof(GloId); 328 329 polygon->n = *((int *) & (buffer[pos])); 330 pos += sizeof(int); 331 332 polygon->area = *((double *) & (buffer[pos])); 333 pos += sizeof(double); 334 335 polygon->x = *((Coord *) & (buffer[pos])); 336 pos += sizeof(Coord); 337 338 elt.is.push_back(polygon); 339 } 340 } 341 342 } -
XIOS/dev/dev_olga/extern/remap/src/polyg.hpp
r688 r1590 23 23 void unpackIntersection(Elt *e, const char *buffer); 24 24 int packIntersectionSize(const Elt& e); 25 double triarea( const Coord& A, const Coord& B, const Coord& C) ; 25 26 26 27 } -
XIOS/dev/dev_olga/src/client.cpp
r1542 r1590 11 11 #include "timer.hpp" 12 12 #include "buffer_client.hpp" 13 #include "string_tools.hpp" 13 14 14 15 namespace xios … … 23 24 StdOFStream CClient::m_infoStream; 24 25 StdOFStream CClient::m_errorStream; 25 26 26 MPI_Comm& CClient::getInterComm(void) { return (interComm); } 27 27 28 28 ///--------------------------------------------------------------- 29 29 /*! … … 230 230 } 231 231 } 232 233 /*! 234 * \fn void CClient::callOasisEnddef(void) 235 * \brief Send the order to the servers to call "oasis_enddef". It must be done by each compound of models before calling oasis_enddef on client side 236 * Function is only called by client. 237 */ 238 void CClient::callOasisEnddef(void) 239 { 240 bool oasisEnddef=CXios::getin<bool>("call_oasis_enddef",true) ; 241 if (!oasisEnddef) ERROR("void CClient::callOasisEnddef(void)", <<"Function xios_oasis_enddef called but variable <call_oasis_enddef> is set to false."<<endl 242 <<"Variable <call_oasis_enddef> must be set to true"<<endl) ; 243 if (CXios::isServer) 244 // Attached mode 245 { 246 // nothing to do 247 } 248 else 249 { 250 int rank ; 251 int msg=0 ; 252 253 MPI_Comm_rank(intraComm,&rank) ; 254 if (rank==0) 255 { 256 MPI_Send(&msg,1,MPI_INT,0,5,interComm) ; // tags oasis_endded = 5 257 } 258 259 } 260 } 261 232 262 233 263 void CClient::finalize(void) -
XIOS/dev/dev_olga/src/client.hpp
r1243 r1590 13 13 static void finalize(void); 14 14 static void registerContext(const string& id, MPI_Comm contextComm); 15 static void callOasisEnddef(void) ; 15 16 16 17 static MPI_Comm intraComm; -
XIOS/dev/dev_olga/src/interface/c/icdata.cpp
r1542 r1590 81 81 CClient::registerContext(str, comm); 82 82 CTimer::get("XIOS init context").suspend(); 83 CTimer::get("XIOS").suspend(); 84 } 85 86 void cxios_oasis_enddef() 87 { 88 CTimer::get("XIOS").resume(); 89 CClient::callOasisEnddef(); 83 90 CTimer::get("XIOS").suspend(); 84 91 } -
XIOS/dev/dev_olga/src/interface/fortran/idata.F90
r1158 r1590 42 42 END SUBROUTINE cxios_context_finalize 43 43 44 44 SUBROUTINE cxios_oasis_enddef() BIND(C) 45 USE ISO_C_BINDING 46 END SUBROUTINE cxios_oasis_enddef 47 45 48 SUBROUTINE cxios_finalize() BIND(C) 46 49 END SUBROUTINE cxios_finalize … … 507 510 END SUBROUTINE xios(finalize) 508 511 512 SUBROUTINE xios(oasis_enddef) 513 IMPLICIT NONE 514 515 CALL cxios_oasis_enddef 516 517 END SUBROUTINE xios(oasis_enddef) 509 518 510 519 SUBROUTINE xios(close_context_definition)() -
XIOS/dev/dev_olga/src/interface/fortran/ixios.F90
r981 r1590 13 13 14 14 USE idata, ONLY : xios(initialize), xios(init_server), xios(finalize), xios(context_initialize), xios(context_is_initialized), & 15 xios(close_context_definition), xios(context_finalize), xios(solve_inheritance) 15 xios(close_context_definition), xios(context_finalize), xios(solve_inheritance), xios(oasis_enddef) 16 16 17 17 USE idomain, ONLY : txios(domain), txios(domaingroup), xios(is_valid_domain), xios(is_valid_domaingroup) -
XIOS/dev/dev_olga/src/io/nc4_data_output.cpp
r1589 r1590 56 56 { 57 57 StdString lonName,latName ; 58 58 59 59 domain->computeWrittenIndex(); 60 60 domain->computeWrittenCompressedIndex(comm_file); -
XIOS/dev/dev_olga/src/server.cpp
r1542 r1590 14 14 #include "timer.hpp" 15 15 #include "event_scheduler.hpp" 16 #include "string_tools.hpp" 16 17 17 18 namespace xios … … 343 344 344 345 string codesId=CXios::getin<string>("oasis_codes_id") ; 345 vector<string> splitted;346 boost::split( splitted, codesId, boost::is_any_of(","), boost::token_compress_on ) ;346 vector<string> oasisCodeId=splitRegex(codesId,"\\s*,\\s*") ; 347 347 348 vector<string>::iterator it ; 348 349 … … 352 353 353 354 // (2) Create interComms with models 354 for(it= splitted.begin();it!=splitted.end();it++)355 for(it=oasisCodeId.begin();it!=oasisCodeId.end();it++) 355 356 { 356 357 oasis_get_intercomm(newComm,*it) ; … … 386 387 } 387 388 if (CXios::usingServer2) delete [] srvGlobalRanks ; 388 oasis_enddef() ; 389 390 bool oasisEnddef=CXios::getin<bool>("call_oasis_enddef",true) ; 391 if (!oasisEnddef) oasis_enddef() ; 389 392 } 390 393 … … 442 445 listenContext(); 443 446 listenRootContext(); 447 listenOasisEnddef() ; 448 listenRootOasisEnddef() ; 444 449 if (!finished) listenFinalize() ; 445 450 } … … 447 452 { 448 453 listenRootContext(); 454 listenRootOasisEnddef() ; 449 455 if (!finished) listenRootFinalize() ; 450 456 } … … 516 522 } 517 523 } 524 525 526 /*! 527 * Root process is listening for an order sent by client to call "oasis_enddef". 528 * The root client of a compound send the order (tag 5). It is probed and received. 529 * When the order has been received from each coumpound, the server root process ping the order to the root processes of the secondary levels of servers (if any). 530 * After, it also inform (asynchronous call) other processes of the communicator that the oasis_enddef call must be done 531 */ 532 533 void CServer::listenOasisEnddef(void) 534 { 535 int flag ; 536 MPI_Status status ; 537 list<MPI_Comm>::iterator it; 538 int msg ; 539 static int nbCompound=0 ; 540 int size ; 541 static bool sent=false ; 542 static MPI_Request* allRequests ; 543 static MPI_Status* allStatus ; 544 545 546 if (sent) 547 { 548 MPI_Comm_size(intraComm,&size) ; 549 MPI_Testall(size,allRequests, &flag, allStatus) ; 550 if (flag==true) 551 { 552 delete [] allRequests ; 553 delete [] allStatus ; 554 sent=false ; 555 } 556 } 557 558 559 for(it=interCommLeft.begin();it!=interCommLeft.end();it++) 560 { 561 MPI_Status status ; 562 traceOff() ; 563 MPI_Iprobe(0,5,*it,&flag,&status) ; // tags oasis_endded = 5 564 traceOn() ; 565 if (flag==true) 566 { 567 MPI_Recv(&msg,1,MPI_INT,0,5,*it,&status) ; // tags oasis_endded = 5 568 nbCompound++ ; 569 if (nbCompound==interCommLeft.size()) 570 { 571 for (std::list<MPI_Comm>::iterator it = interCommRight.begin(); it != interCommRight.end(); it++) 572 { 573 MPI_Send(&msg,1,MPI_INT,0,5,*it) ; // tags oasis_endded = 5 574 } 575 MPI_Comm_size(intraComm,&size) ; 576 allRequests= new MPI_Request[size] ; 577 allStatus= new MPI_Status[size] ; 578 for(int i=0;i<size;i++) MPI_Isend(&msg,1,MPI_INT,i,5,intraComm,&allRequests[i]) ; // tags oasis_endded = 5 579 sent=true ; 580 } 581 } 582 } 583 } 584 585 /*! 586 * Processes probes message from root process if oasis_enddef call must be done. 587 * When the order is received it is scheduled to be treated in a synchronized way by all server processes of the communicator 588 */ 589 void CServer::listenRootOasisEnddef(void) 590 { 591 int flag ; 592 MPI_Status status ; 593 const int root=0 ; 594 int msg ; 595 static bool eventSent=false ; 596 597 if (eventSent) 598 { 599 boost::hash<string> hashString; 600 size_t hashId = hashString("oasis_enddef"); 601 if (eventScheduler->queryEvent(0,hashId)) 602 { 603 oasis_enddef() ; 604 eventSent=false ; 605 } 606 } 607 608 traceOff() ; 609 MPI_Iprobe(root,5,intraComm, &flag, &status) ; 610 traceOn() ; 611 if (flag==true) 612 { 613 MPI_Recv(&msg,1,MPI_INT,root,5,intraComm,&status) ; // tags oasis_endded = 5 614 boost::hash<string> hashString; 615 size_t hashId = hashString("oasis_enddef"); 616 eventScheduler->registerEvent(0,hashId); 617 eventSent=true ; 618 } 619 } 620 621 622 623 518 624 519 625 void CServer::listenContext(void) -
XIOS/dev/dev_olga/src/server.hpp
r1378 r1590 22 22 static void listenRootContext(void); 23 23 static void listenRootFinalize(void); 24 static void listenRootOasisEnddef(void); 25 static void listenOasisEnddef(void); 24 26 static void registerContext(void* buff,int count, int leaderRank=0); 25 27 … … 72 74 static StdOFStream m_infoStream; 73 75 static StdOFStream m_errorStream; 74 75 76 static void openStream(const StdString& fileName, const StdString& ext, std::filebuf* fb); 76 77 };
Note: See TracChangeset
for help on using the changeset viewer.