Ignore:
Timestamp:
06/06/17 17:58:16 (3 years ago)
Author:
oabramkina
Message:

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

Location:
XIOS/dev/dev_olga/extern/remap/src
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • XIOS/dev/dev_olga/extern/remap/src/elt.hpp

    r923 r1158  
    105105        } 
    106106 
     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   
    107120        int neighbour[NMAX]; 
    108121        double d[NMAX]; /**< distance of centre of small circle to origin, zero if great circle */ 
  • XIOS/dev/dev_olga/extern/remap/src/intersect.cpp

    r688 r1158  
    3737} 
    3838 
     39/** New methods to find an insert a neighbour in a cell of the source mesh. 
     40 *  return true/false if cell b is a neighbour of a. if "insert" is true, then b will be inserted as a neighbour 
     41 * in cell a . This is needed for 2 order interpolation that need neighboround for gradient computing. 
     42 * A cell is a neighbour if : 
     43 *  - it shares 2 countiguous vertex (ie an edge) with a 
     44 *  - A vertex of b is located on one of an edge of a. 
     45 **/ 
     46bool insertNeighbour( Elt& a, const Elt& b, bool insert ) 
     47{ 
     48  // for now suppose pole -> Oz 
     49  Coord pole(0,0,1) ; 
     50  Coord O, Oa1, Oa2,Ob1,Ob2,V1,V2 ; 
     51  double da,db,alpha,alpha1,alpha2,delta ; 
     52   
     53     
     54  for (int i = 0; i < a.n; i++) 
     55  { 
     56    for (int j = 0; j < b.n; j++) 
     57    { 
     58// share a full edge ? be carefull at the orientation 
     59      assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > 1e-10*1e-10 || 
     60             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > 1e-10*1e-10); 
     61      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-10*1e-10 && 
     62             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-10*1e-10) 
     63      { 
     64        if (insert) a.neighbour[i] = b.id.ind ; 
     65        return true; 
     66      } 
     67       
     68 
     69// 1 or 2 vertices of b is located on an edge of a 
     70       da=a.d[i] ; 
     71       if (scalarprod(a.edge[i], pole) < 0) da=-da ; 
     72       db=b.d[(j+b.n-1)%b.n] ; 
     73       if (scalarprod(b.edge[(j+b.n-1)%b.n], pole) < 0) db=-db ; 
     74       
     75      if ( fabs(da-db)<1e-10 )  
     76      { 
     77        O=pole*da ; 
     78        Oa1=a.vertex[i]-O ; 
     79        Oa2=a.vertex[(i+1)%a.n]-O ;  
     80        Ob1=b.vertex[j]-O ; 
     81        Ob2=b.vertex[(j+b.n-1)%b.n]-O ; 
     82        V1=crossprod(Oa1,Oa2) ; 
     83        V2=crossprod(Ob1,Ob2) ; 
     84        if (norm(crossprod(V1,V2))/(norm(V1)*norm(V2)) < 1e-10) 
     85        { 
     86          alpha = vectAngle(Oa1,Oa2,V1) ; 
     87          alpha1= vectAngle(Oa1,Ob1,V1) ; 
     88          alpha2= vectAngle(Oa1,Ob2,V1) ; 
     89          delta= alpha2-alpha1 ; 
     90          if (delta >= M_PI) delta=2*M_PI-delta ; 
     91          else if (delta <= -M_PI) delta=2*M_PI+delta ; 
     92           
     93          if (alpha >= 0) 
     94          { 
     95            if (alpha1 > 1e-10 && alpha1 < alpha-1e-10) 
     96            { 
     97              if (alpha2 > 1e-10 && alpha2 < alpha-1e-10) 
     98              { 
     99                assert(delta > 0) ; 
     100                if (insert) 
     101                { 
     102                // insert both 
     103                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     104                  a.insert_vertex(i,b.vertex[j]); 
     105                  a.neighbour[i+1] = b.id.ind ; 
     106                } 
     107                return true ; 
     108              } 
     109              else 
     110              { 
     111                assert( delta > 0 ) ; 
     112                if (insert) 
     113                { 
     114                //insert alpha1 
     115                  a.insert_vertex(i,b.vertex[j]); 
     116                  a.neighbour[i+1] = b.id.ind ; 
     117                } 
     118                return true ; 
     119              } 
     120            } 
     121            else if (alpha2 > 1e-10 && alpha2 < alpha-1e-10) 
     122            { 
     123              assert( delta > 0 ) ; 
     124              if (insert) 
     125              { 
     126              // insert alpha2 
     127                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     128                a.neighbour[i] = b.id.ind ; 
     129              } 
     130              return true ; 
     131            } 
     132            else 
     133            { 
     134              // nothing to do 
     135            }  
     136 
     137          } 
     138          else  // alpha < 0 
     139          { 
     140            if (alpha1 < -1e-10 && alpha1 > alpha+1e-10) 
     141            { 
     142              if (alpha2 < -1e-10 && alpha2 > alpha+1e-10) 
     143              { 
     144                assert(delta < 0) ; 
     145                if (insert) 
     146                { 
     147                // insert both 
     148                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     149                  a.insert_vertex(i,b.vertex[j]); 
     150                  a.neighbour[i+1] = b.id.ind ; 
     151                } 
     152                return true ; 
     153              } 
     154              else 
     155              { 
     156                assert(delta < 0) ; 
     157                if (insert) 
     158                { 
     159                //insert alpha1 
     160                  a.insert_vertex(i,b.vertex[j]); 
     161                  a.neighbour[i+1] = b.id.ind ; 
     162                } 
     163                return true ; 
     164              } 
     165            } 
     166            else if (alpha2 < -1e-10 && alpha2 > alpha+1e-10) 
     167            { 
     168              assert(delta < 0) ; 
     169              if (insert) 
     170              { 
     171              // insert alpha2 
     172                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     173                a.neighbour[i] = b.id.ind ; 
     174              } 
     175              return true ; 
     176            } 
     177            else 
     178            { 
     179               // nothing to do 
     180            } 
     181          } 
     182        } 
     183      } 
     184    } 
     185  } 
     186  return false; 
     187} 
    39188 
    40189/** 
     
    44193void set_neighbour(Elt& a, const Elt& b) 
    45194{ 
    46         if (b.id.ind == a.id.ind) return; 
    47         int idx = neighbour_idx(a, b); 
    48         if (idx != NOT_FOUND) 
    49                 a.neighbour[idx] = b.id.ind; 
     195  if (b.id.ind == a.id.ind) return; 
     196/* 
     197  int idx = neighbour_idx(a, b); 
     198  if (idx != NOT_FOUND) 
     199  a.neighbour[idx] = b.id.ind; 
     200*/ 
     201  insertNeighbour(a,b,true) ;  
    50202} 
    51203 
    52204/** return true if `a` and `b` share an edge */ 
    53 bool isNeighbour(const Elt& a, const Elt& b) 
    54 { 
    55         return neighbour_idx(a, b) != NOT_FOUND; 
     205bool isNeighbour(Elt& a, const Elt& b) 
     206{ 
     207        // return neighbour_idx(a, b) != NOT_FOUND; 
     208  return insertNeighbour(a,b,false) ; 
    56209} 
    57210 
  • XIOS/dev/dev_olga/extern/remap/src/intersect.hpp

    r688 r1158  
    22 
    33void set_neighbour(Elt& elt, const Elt& elt2); 
    4 bool isNeighbour(const Elt& elt, const Elt& elt2); 
     4bool isNeighbour(Elt& elt, const Elt& elt2); 
    55 
    66void intersect(Elt *a, Elt *b); 
  • XIOS/dev/dev_olga/extern/remap/src/mapper.cpp

    r923 r1158  
    109109} 
    110110 
    111 vector<double> Mapper::computeWeights(int interpOrder, bool renormalize) 
     111vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 
    112112{ 
    113113        vector<double> timings; 
     
    152152        if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
    153153        tic = cputime(); 
    154         nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize); 
     154        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
    155155        timings.push_back(cputime() - tic); 
    156156 
     
    166166   @param order is the order of interpolaton (must be 1 or 2). 
    167167*/ 
    168 int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize) 
     168int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 
    169169{ 
    170170        int mpiSize, mpiRank; 
     
    184184        int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
    185185        double **recvValue = new double*[mpiSize]; 
     186        double **recvArea = new double*[mpiSize]; 
    186187        Coord **recvGrad = new Coord*[mpiSize]; 
    187188        GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
     
    205206                        sendElement[rank] = new int[nbSendElement[rank]]; 
    206207                        recvValue[rank]   = new double[nbSendElement[rank]]; 
     208                        recvArea[rank]    = new double[nbSendElement[rank]]; 
    207209                        if (order == 2) 
    208210                        { 
     
    234236        int **recvElement = new int*[mpiSize]; 
    235237        double **sendValue = new double*[mpiSize]; 
     238        double **sendArea = new double*[mpiSize]; 
    236239        Coord **sendGrad = new Coord*[mpiSize]; 
    237240        GloId **sendNeighIds = new GloId*[mpiSize]; 
    238         MPI_Request *sendRequest = new MPI_Request[3*mpiSize]; 
    239         MPI_Request *recvRequest = new MPI_Request[3*mpiSize]; 
     241        MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
     242        MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
    240243        for (int rank = 0; rank < mpiSize; rank++) 
    241244        { 
     
    250253                        recvElement[rank] = new int[nbRecvElement[rank]]; 
    251254                        sendValue[rank]   = new double[nbRecvElement[rank]]; 
     255                        sendArea[rank]   = new double[nbRecvElement[rank]]; 
    252256                        if (order == 2) 
    253257                        { 
     
    263267                } 
    264268        } 
    265         MPI_Status *status = new MPI_Status[3*mpiSize]; 
    266         MPI_Waitall(nbRecvRequest, recvRequest, status); 
     269        MPI_Status *status = new MPI_Status[4*mpiSize]; 
     270         
    267271        MPI_Waitall(nbSendRequest, sendRequest, status); 
     272        MPI_Waitall(nbRecvRequest, recvRequest, status); 
    268273 
    269274        /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
     
    278283                        { 
    279284                                sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
     285                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
    280286                                if (order == 2) 
    281287                                { 
     
    297303                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    298304                        nbSendRequest++; 
     305                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     306                        nbSendRequest++; 
    299307                        if (order == 2) 
    300308                        { 
     
    317325                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    318326                        nbRecvRequest++; 
     327                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     328                        nbRecvRequest++; 
    319329                        if (order == 2) 
    320330                        { 
     
    334344                } 
    335345        } 
     346         
     347        MPI_Waitall(nbSendRequest, sendRequest, status); 
    336348        MPI_Waitall(nbRecvRequest, recvRequest, status); 
    337         MPI_Waitall(nbSendRequest, sendRequest, status); 
     349         
    338350 
    339351        /* now that all values and gradients are available use them to computed interpolated values on target 
     
    357369                        int rank = (*it)->id.rank; 
    358370                        double fk = recvValue[rank][n1]; 
     371                        double srcArea = recvArea[rank][n1]; 
    359372                        double w = (*it)->area; 
     373      if (quantity) w/=srcArea ; 
    360374 
    361375                        /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
    362376                        int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
    363377                        GloId neighID = recvNeighIds[rank][kk]; 
    364                         wgt_map[neighID] += (*it)->area; 
     378                        wgt_map[neighID] += w; 
    365379 
    366380                        if (order == 2) 
     
    383397    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    384398                { 
    385                         this->remapMatrix[i] = (it->second / e.area) / renorm; 
     399      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
     400                        else this->remapMatrix[i] = (it->second / e.area) / renorm; 
    386401                        this->srcAddress[i] = it->first.ind; 
    387402                        this->srcRank[i] = it->first.rank; 
     
    400415                        delete[] sendElement[rank]; 
    401416                        delete[] recvValue[rank]; 
     417                        delete[] recvArea[rank]; 
    402418                        if (order == 2) 
    403419                        { 
     
    410426                        delete[] recvElement[rank]; 
    411427                        delete[] sendValue[rank]; 
     428                        delete[] sendArea[rank]; 
    412429                        if (order == 2) 
    413430                                delete[] sendGrad[rank]; 
  • XIOS/dev/dev_olga/extern/remap/src/mapper.hpp

    r844 r1158  
    3434       /** @param trgElts are the elements of the unstructured target grid 
    3535           Returns the timings for substeps: */ 
    36        vector<double> computeWeights(int interpOrder, bool renormalize=false); 
     36       vector<double> computeWeights(int interpOrder, bool renormalize=false, bool quantity=false); 
    3737       int getNbWeights(void) { return nWeights ; } 
    3838/* 
     
    5151private: 
    5252       /** @return number of weights (local to cpu) */ 
    53        int remap(Elt* elements, int nbElements, int order, bool renormalize=false); 
     53       int remap(Elt* elements, int nbElements, int order, bool renormalize=false, bool quantity=false); 
    5454 
    5555       void buildMeshTopology(); 
  • XIOS/dev/dev_olga/extern/remap/src/tree.cpp

    r951 r1158  
    142142        root->parent = 0; 
    143143        root->leafCount = 0; 
    144         root->centre = ORIGIN; 
     144// initialize root node on the sphere 
     145  root->centre.x=1 ; root->centre.y=0 ; root->centre.z=0 ;  
    145146        root->radius = 0.; 
    146147        root->reinserted = false; 
  • XIOS/dev/dev_olga/extern/remap/src/triple.cpp

    r849 r1158  
    124124} 
    125125 
     126// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b 
     127double vectAngle(const Coord &a, const Coord &b, const Coord &pole) 
     128{  
     129  double nab = 1./(norm(a)*norm(b)) ; 
     130   
     131  Coord a_cross_b=crossprod(a, b)*nab ; 
     132  double sinVect ; 
     133  if (scalarprod(a_cross_b, pole) >= 0) sinVect=norm(a_cross_b) ; 
     134  else sinVect=-norm(a_cross_b) ; 
     135  double cosVect=scalarprod(a,b)*nab ; 
     136 
     137  return atan2(sinVect,cosVect) ; 
    126138} 
     139 
     140} 
  • XIOS/dev/dev_olga/extern/remap/src/triple.hpp

    r688 r1158  
    121121 
    122122double angle(const Coord &a, const Coord &b, const Coord &pole); 
     123 
     124// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b 
     125double vectAngle(const Coord &a, const Coord &b, const Coord &pole) ; 
     126 
    123127void print_Coord(Coord &p); 
    124128 
Note: See TracChangeset for help on using the changeset viewer.