Changeset 844


Ignore:
Timestamp:
04/27/16 17:23:48 (6 years ago)
Author:
ymipsl
Message:

Management of masked cell for remapping algorithm.

YM

Location:
XIOS/trunk/extern/remap/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/extern/remap/src/mapper.cpp

    r694 r844  
    109109} 
    110110 
    111 vector<double> Mapper::computeWeights(int interpOrder) 
     111vector<double> Mapper::computeWeights(int interpOrder, bool renormalize) 
    112112{ 
    113113        vector<double> timings; 
     
    152152        if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
    153153        tic = cputime(); 
    154         nWeights = remap(&targetElements[0], targetElements.size(), interpOrder); 
     154        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize); 
    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) 
     168int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize) 
    169169{ 
    170170        int mpiSize, mpiRank; 
     
    281281                                { 
    282282                                        sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
     283//          cout<<"grad  "<<jj<<"  "<<recvElement[rank][j]<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 
    283284                                        sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
    284285                                        jj++; 
     
    286287                                        { 
    287288                                                sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
    288                                                 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
     289//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 
     290            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
    289291                                                jj++; 
    290292                                        } 
     
    368370                                        int kk = n1 * (NMAX + 1) + k; 
    369371                                        GloId neighID = recvNeighIds[rank][kk]; 
    370                                         if (neighID.ind == -1) break; 
    371                                         wgt_map[neighID] += 
    372                                                 w * scalarprod(recvGrad[rank][kk], (*it)->x); 
     372                                        if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
    373373                                } 
    374374 
    375375                        } 
    376376                } 
    377                 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    378                 { 
    379                         this->remapMatrix[i] = it->second / e.area; 
     377 
     378    double renorm=0; 
     379    if (renormalize)  
     380      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 
     381    else renorm=1. ; 
     382 
     383    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
     384                { 
     385                        this->remapMatrix[i] = (it->second / e.area) / renorm; 
    380386                        this->srcAddress[i] = it->first.ind; 
    381387                        this->srcRank[i] = it->first.rank; 
     
    669675                } 
    670676 
     677/* 
    671678                for (int i = 0; i < elt->n; i++) 
    672679                { 
     
    674681                                error_exit("neighbour not found"); 
    675682                } 
     683*/ 
    676684        } 
    677685} 
  • XIOS/trunk/extern/remap/src/mapper.hpp

    r694 r844  
    3434       /** @param trgElts are the elements of the unstructured target grid 
    3535           Returns the timings for substeps: */ 
    36        vector<double> computeWeights(int interpOrder); 
     36       vector<double> computeWeights(int interpOrder, bool renormalize=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); 
     53       int remap(Elt* elements, int nbElements, int order, bool renormalize=false); 
    5454 
    5555       void buildMeshTopology(); 
  • XIOS/trunk/extern/remap/src/meshutil.cpp

    r688 r844  
    99void cptEltGeom(Elt& elt, const Coord &pole) 
    1010{ 
    11         orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x); 
    12         normals(elt, pole); 
    13         Coord gg; 
    14         elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg); 
    15         elt.x = gg; 
     11  orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x); 
     12  normals(elt, pole); 
     13  Coord gg; 
     14  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg); 
     15  elt.x = gg; 
    1616} 
    1717 
    1818void cptAllEltsGeom(Elt *elt, int N, const Coord &pole) 
    1919{ 
    20         for (int ne=0; ne<N; ne++) 
    21                 cptEltGeom(elt[ne], pole); 
     20  for (int ne=0; ne<N; ne++) 
     21    cptEltGeom(elt[ne], pole); 
    2222} 
    2323 
     
    2626void update_baryc(Elt *elt, int N) 
    2727{ 
    28         for (int ne=0; ne<N; ne++) 
    29         { 
    30                 Elt &e = elt[ne]; 
    31                 int ns = e.is.size();  // sous-elements 
    32                 Coord *sx = new Coord[ns]; 
    33                 int i=0; 
    34                 for (list<Polyg*>::iterator it = e.is.begin(); it != e.is.end(); i++, it++) 
    35                 { 
    36                         sx[i] = (*it)->x * (*it)->area; 
    37                 } 
    38                 e.x = barycentre(sx, ns); 
    39         } 
     28  for (int ne=0; ne<N; ne++) 
     29  { 
     30    Elt &e = elt[ne]; 
     31    int ns = e.is.size();  // sous-elements 
     32    Coord *sx = new Coord[ns]; 
     33    int i=0; 
     34    for (list<Polyg*>::iterator it = e.is.begin(); it != e.is.end(); i++, it++) 
     35    { 
     36      sx[i] = (*it)->x * (*it)->area; 
     37    } 
     38    e.x = barycentre(sx, ns); 
     39  } 
    4040} 
     41 
     42 
     43Coord gradient_old(Elt& elt, Elt **neighElts) 
     44{ 
     45  Coord grad = ORIGIN; 
     46  Coord *neighBaryc = new Coord[elt.n]; 
     47  for (int j = 0; j < elt.n; j++) 
     48  { 
     49    int k = (j + 1) % elt.n; 
     50    neighBaryc[j] = neighElts[j]->x; 
     51    Coord edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x); 
     52 
     53    // use nomenclauture form paper  
     54    double f_i = elt.val; 
     55    double f_j = neighElts[j]->val; 
     56    double f_k = neighElts[k]->val; 
     57    grad = grad + edgeNormal * (0.5*(f_j + f_k) - f_i); 
     58  } 
     59  // area of the polygon whoes vertices are the barycentres the neighbours  
     60  grad = grad * (1./polygonarea(neighBaryc, elt.n)); 
     61  delete[] neighBaryc; 
     62 
     63  return grad - elt.x * scalarprod(elt.x, grad); 
     64} 
     65 
     66 
    4167 
    4268Coord gradient(Elt& elt, Elt **neighElts) 
    4369{ 
    44         Coord grad = ORIGIN; 
    45         Coord *neighBaryc = new Coord[elt.n]; 
    46         for (int j = 0; j < elt.n; j++) 
    47         { 
    48                 int k = (j + 1) % elt.n; 
    49                 neighBaryc[j] = neighElts[j]->x; 
    50                 Coord edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x); 
     70     
     71  Coord grad = ORIGIN; 
     72  Coord neighBaryc[3] ; 
    5173 
    52                 /* use nomenclauture form paper */ 
    53                 double f_i = elt.val; 
    54                 double f_j = neighElts[j]->val; 
    55                 double f_k = neighElts[k]->val; 
    56                 grad = grad + edgeNormal * (0.5*(f_j + f_k) - f_i); 
    57         } 
    58         /* area of the polygon whoes vertices are the barycentres the neighbours */ 
    59         grad = grad * (1./polygonarea(neighBaryc, elt.n)); 
    60         delete[] neighBaryc; 
     74  double f_i ; 
     75  double f_j ; 
     76  double f_k ; 
     77  
     78  Coord edgeNormal ; 
     79  double area=0 ; 
     80  int k ; 
    6181 
    62         return grad - elt.x * scalarprod(elt.x, grad); 
     82  for (int j = 0; j < elt.n; j++) 
     83  { 
     84    f_i = elt.val; 
     85    k = (j + 1) % elt.n; 
     86    if (neighElts[j]==NULL || neighElts[k]==NULL) continue ; 
     87 
     88    // use nomenclauture form paper  
     89    f_j = neighElts[j]->val; 
     90    f_k = neighElts[k]->val; 
     91 
     92     
     93    
     94    neighBaryc[0] = elt.x; 
     95    neighBaryc[1] = neighElts[j]->x; 
     96    neighBaryc[2] = neighElts[k]->x; 
     97 
     98    edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x); 
     99    grad = grad + edgeNormal * (0.5*(f_k + f_j) - f_i); 
     100 
     101    edgeNormal = crossprod(neighElts[j]->x, elt.x); 
     102    grad = grad + edgeNormal * (0.5*(f_j + f_i) - f_i); 
     103 
     104    edgeNormal = crossprod(elt.x, neighElts[k]->x); 
     105    grad = grad + edgeNormal * (0.5*(f_i + f_k) - f_i); 
     106 
     107  // area of the polygon whoes vertices are the barycentres the neighbours  
     108 
     109    area+=polygonarea(neighBaryc, 3) ; 
     110 
     111  } 
     112  grad=grad*(1./area) ; 
     113  return grad - elt.x * scalarprod(elt.x, grad); 
    63114} 
     115 
     116 
     117 
    64118 
    65119void computeGradients(Elt **elts, int N) 
    66120{ 
    67          
    68         for (int j = 0; j < N; j++) 
    69                 elts[j]->val = 0; 
     121   
     122  for (int j = 0; j < N; j++) 
     123    elts[j]->val = 0; 
    70124 
    71         Elt *neighbours[NMAX]; 
    72         for (int j = 0; j < N; j++) 
    73         { 
    74                 for (int i = 0; i < elts[j]->n; i++) 
    75                         neighbours[i] = elts[elts[j]->neighbour[i]]; 
     125  Elt *neighbours[NMAX]; 
     126  for (int j = 0; j < N; j++) 
     127  { 
     128    for (int i = 0; i < elts[j]->n; i++)  
     129      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; 
     130      else  neighbours[i] = elts[elts[j]->neighbour[i]]; 
    76131 
    77                 for (int i = 0; i < elts[j]->n; i++) 
    78                         neighbours[i]->val = 0; 
    79                 for (int i = 0; i < elts[j]->n; i++) 
    80                 { 
    81                         elts[j]->neighId[i] = neighbours[i]->src_id; 
    82                         /* for weight computation all values are always kept zero and only set to one when used .. */ 
    83                         neighbours[i]->val = 1; 
    84                         elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours); 
    85                         /* .. and right after zeroed again */ 
    86                         neighbours[i]->val = 0; 
    87                 } 
    88                 elts[j]->neighId[elts[j]->n].rank = -1; // mark end 
    89                 elts[j]->neighId[elts[j]->n].ind = -1; // mark end 
    90                 /* For the most naive algorithm the case where the element itself is one must also be considered. 
    91                    Thomas says this can later be optimized out. */ 
    92                 elts[j]->val = 1; 
    93                 elts[j]->grad = gradient(*(elts[j]), neighbours); 
    94                 elts[j]->val = 0; 
    95         } 
     132    for (int i = 0; i < elts[j]->n; i++) 
     133      if (neighbours[i]!=NULL) neighbours[i]->val = 0; 
     134       
     135    for (int i = 0; i < elts[j]->n; i++) 
     136    { 
     137      if (neighbours[i]!=NULL) 
     138      { 
     139        elts[j]->neighId[i] = neighbours[i]->src_id; 
     140        /* for weight computation all values are always kept zero and only set to one when used .. */ 
     141        neighbours[i]->val = 1; 
     142        elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours); 
     143        /* .. and right after zeroed again */ 
     144        neighbours[i]->val = 0; 
     145      } 
     146      else 
     147      { 
     148        elts[j]->neighId[i].rank = -1; // mark end 
     149        elts[j]->neighId[i].ind = -1; // mark end 
     150      } 
     151    } 
     152 
     153    for(int i = elts[j]->n ; i < NMAX; i++) 
     154    { 
     155      elts[j]->neighId[i].rank = -1; // mark end 
     156      elts[j]->neighId[i].ind = -1; // mark end 
     157    } 
     158    /* For the most naive algorithm the case where the element itself is one must also be considered. 
     159       Thomas says this can later be optimized out. */ 
     160    elts[j]->val = 1; 
     161    elts[j]->grad = gradient(*(elts[j]), neighbours); 
     162    elts[j]->val = 0; 
     163  } 
    96164} 
    97165 
Note: See TracChangeset for help on using the changeset viewer.