Changeset 844 for XIOS/trunk/extern

Timestamp:

04/27/16 17:23:48 (8 years ago)

Author:

ymipsl

Message:

Management of masked cell for remapping algorithm.

YM

Location:

XIOS/trunk/extern/remap/src

Files:

• mapper.cpp (modified) (8 diffs)
• mapper.hpp (modified) (2 diffs)
• meshutil.cpp (modified) (2 diffs)

XIOS/trunk/extern/remap/src/mapper.cpp

@@ -109 +109 @@
 }
 
-vector<double> Mapper::computeWeights(int interpOrder)
+vector<double> Mapper::computeWeights(int interpOrder, bool renormalize)
 {
         vector<double> timings;
@@ -152 +152 @@
         if (mpiRank == 0 && verbose) cout << "Remapping..." << endl;
         tic = cputime();
-        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder);
+        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize);
         timings.push_back(cputime() - tic);
 
@@ -166 +166 @@
    @param order is the order of interpolaton (must be 1 or 2).
 */
-int Mapper::remap(Elt *elements, int nbElements, int order)
+int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize)
 {
         int mpiSize, mpiRank;
@@ -281 +281 @@
                                 {
                                         sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad;
+//          cout<<"grad  "<<jj<<"  "<<recvElement[rank][j]<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
                                         sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id;
                                         jj++;
@@ -286 +287 @@
                                         {
                                                 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i];
-                                                sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i];
+//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
+            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i];
                                                 jj++;
                                         }
@@ -368 +370 @@
                                         int kk = n1 * (NMAX + 1) + k;
                                         GloId neighID = recvNeighIds[rank][kk];
-                                        if (neighID.ind == -1) break;
-                                        wgt_map[neighID] +=
-                                                w * scalarprod(recvGrad[rank][kk], (*it)->x);
+                                        if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x);
                                 }
 
                         }
                 }
-                for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++)
-                {
-                        this->remapMatrix[i] = it->second / e.area;
+
+    double renorm=0;
+    if (renormalize) 
+      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area;
+    else renorm=1. ;
+
+    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++)
+                {
+                        this->remapMatrix[i] = (it->second / e.area) / renorm;
                         this->srcAddress[i] = it->first.ind;
                         this->srcRank[i] = it->first.rank;
@@ -669 +675 @@
                 }
 
+/*
                 for (int i = 0; i < elt->n; i++)
                 {
@@ -674 +681 @@
                                 error_exit("neighbour not found");
                 }
+*/
         }
 }

XIOS/trunk/extern/remap/src/mapper.hpp

@@ -34 +34 @@
        /** @param trgElts are the elements of the unstructured target grid
            Returns the timings for substeps: */
-       vector<double> computeWeights(int interpOrder);
+       vector<double> computeWeights(int interpOrder, bool renormalize=false);
        int getNbWeights(void) { return nWeights ; }
 /*
@@ -51 +51 @@
 private:
        /** @return number of weights (local to cpu) */
-       int remap(Elt* elements, int nbElements, int order);
+       int remap(Elt* elements, int nbElements, int order, bool renormalize=false);
 
        void buildMeshTopology();

XIOS/trunk/extern/remap/src/meshutil.cpp

@@ -9 +9 @@
 void cptEltGeom(Elt& elt, const Coord &pole)
 {
-        orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x);
-        normals(elt, pole);
-        Coord gg;
-        elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
-        elt.x = gg;
+  orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x);
+  normals(elt, pole);
+  Coord gg;
+  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
+  elt.x = gg;
 }
 
 void cptAllEltsGeom(Elt *elt, int N, const Coord &pole)
 {
-        for (int ne=0; ne<N; ne++)
-                cptEltGeom(elt[ne], pole);
+  for (int ne=0; ne<N; ne++)
+    cptEltGeom(elt[ne], pole);
 }
 
@@ -26 +26 @@
 void update_baryc(Elt *elt, int N)
 {
-        for (int ne=0; ne<N; ne++)
-        {
-                Elt &e = elt[ne];
-                int ns = e.is.size();  // sous-elements
-                Coord *sx = new Coord[ns];
-                int i=0;
-                for (list<Polyg*>::iterator it = e.is.begin(); it != e.is.end(); i++, it++)
-                {
-                        sx[i] = (*it)->x * (*it)->area;
-                }
-                e.x = barycentre(sx, ns);
-        }
+  for (int ne=0; ne<N; ne++)
+  {
+    Elt &e = elt[ne];
+    int ns = e.is.size();  // sous-elements
+    Coord *sx = new Coord[ns];
+    int i=0;
+    for (list<Polyg*>::iterator it = e.is.begin(); it != e.is.end(); i++, it++)
+    {
+      sx[i] = (*it)->x * (*it)->area;
+    }
+    e.x = barycentre(sx, ns);
+  }
 }
+
+
+Coord gradient_old(Elt& elt, Elt **neighElts)
+{
+  Coord grad = ORIGIN;
+  Coord *neighBaryc = new Coord[elt.n];
+  for (int j = 0; j < elt.n; j++)
+  {
+    int k = (j + 1) % elt.n;
+    neighBaryc[j] = neighElts[j]->x;
+    Coord edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x);
+
+    // use nomenclauture form paper 
+    double f_i = elt.val;
+    double f_j = neighElts[j]->val;
+    double f_k = neighElts[k]->val;
+    grad = grad + edgeNormal * (0.5*(f_j + f_k) - f_i);
+  }
+  // area of the polygon whoes vertices are the barycentres the neighbours 
+  grad = grad * (1./polygonarea(neighBaryc, elt.n));
+  delete[] neighBaryc;
+
+  return grad - elt.x * scalarprod(elt.x, grad);
+}
+
+
 
 Coord gradient(Elt& elt, Elt **neighElts)
 {
-        Coord grad = ORIGIN;
-        Coord *neighBaryc = new Coord[elt.n];
-        for (int j = 0; j < elt.n; j++)
-        {
-                int k = (j + 1) % elt.n;
-                neighBaryc[j] = neighElts[j]->x;
-                Coord edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x);
+    
+  Coord grad = ORIGIN;
+  Coord neighBaryc[3] ;
 
-                /* use nomenclauture form paper */
-                double f_i = elt.val;
-                double f_j = neighElts[j]->val;
-                double f_k = neighElts[k]->val;
-                grad = grad + edgeNormal * (0.5*(f_j + f_k) - f_i);
-        }
-        /* area of the polygon whoes vertices are the barycentres the neighbours */
-        grad = grad * (1./polygonarea(neighBaryc, elt.n));
-        delete[] neighBaryc;
+  double f_i ;
+  double f_j ;
+  double f_k ;
+ 
+  Coord edgeNormal ;
+  double area=0 ;
+  int k ;
 
-        return grad - elt.x * scalarprod(elt.x, grad);
+  for (int j = 0; j < elt.n; j++)
+  {
+    f_i = elt.val;
+    k = (j + 1) % elt.n;
+    if (neighElts[j]==NULL || neighElts[k]==NULL) continue ;
+
+    // use nomenclauture form paper 
+    f_j = neighElts[j]->val;
+    f_k = neighElts[k]->val;
+
+    
+   
+    neighBaryc[0] = elt.x;
+    neighBaryc[1] = neighElts[j]->x;
+    neighBaryc[2] = neighElts[k]->x;
+
+    edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x);
+    grad = grad + edgeNormal * (0.5*(f_k + f_j) - f_i);
+
+    edgeNormal = crossprod(neighElts[j]->x, elt.x);
+    grad = grad + edgeNormal * (0.5*(f_j + f_i) - f_i);
+
+    edgeNormal = crossprod(elt.x, neighElts[k]->x);
+    grad = grad + edgeNormal * (0.5*(f_i + f_k) - f_i);
+
+  // area of the polygon whoes vertices are the barycentres the neighbours 
+
+    area+=polygonarea(neighBaryc, 3) ;
+
+  }
+  grad=grad*(1./area) ;
+  return grad - elt.x * scalarprod(elt.x, grad);
 }
+
+
+
 
 void computeGradients(Elt **elts, int N)
 {
-        
-        for (int j = 0; j < N; j++)
-                elts[j]->val = 0;
+  
+  for (int j = 0; j < N; j++)
+    elts[j]->val = 0;
 
-        Elt *neighbours[NMAX];
-        for (int j = 0; j < N; j++)
-        {
-                for (int i = 0; i < elts[j]->n; i++)
-                        neighbours[i] = elts[elts[j]->neighbour[i]];
+  Elt *neighbours[NMAX];
+  for (int j = 0; j < N; j++)
+  {
+    for (int i = 0; i < elts[j]->n; i++) 
+      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ;
+      else  neighbours[i] = elts[elts[j]->neighbour[i]];
 
-                for (int i = 0; i < elts[j]->n; i++)
-                        neighbours[i]->val = 0;
-                for (int i = 0; i < elts[j]->n; i++)
-                {
-                        elts[j]->neighId[i] = neighbours[i]->src_id;
-                        /* for weight computation all values are always kept zero and only set to one when used .. */
-                        neighbours[i]->val = 1;
-                        elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours);
-                        /* .. and right after zeroed again */
-                        neighbours[i]->val = 0;
-                }
-                elts[j]->neighId[elts[j]->n].rank = -1; // mark end
-                elts[j]->neighId[elts[j]->n].ind = -1; // mark end
-                /* For the most naive algorithm the case where the element itself is one must also be considered.
-                   Thomas says this can later be optimized out. */
-                elts[j]->val = 1;
-                elts[j]->grad = gradient(*(elts[j]), neighbours);
-                elts[j]->val = 0;
-        }
+    for (int i = 0; i < elts[j]->n; i++)
+      if (neighbours[i]!=NULL) neighbours[i]->val = 0;
+      
+    for (int i = 0; i < elts[j]->n; i++)
+    {
+      if (neighbours[i]!=NULL)
+      {
+        elts[j]->neighId[i] = neighbours[i]->src_id;
+        /* for weight computation all values are always kept zero and only set to one when used .. */
+        neighbours[i]->val = 1;
+        elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours);
+        /* .. and right after zeroed again */
+        neighbours[i]->val = 0;
+      }
+      else
+      {
+        elts[j]->neighId[i].rank = -1; // mark end
+        elts[j]->neighId[i].ind = -1; // mark end
+      }
+    }
+
+    for(int i = elts[j]->n ; i < NMAX; i++)
+    {
+      elts[j]->neighId[i].rank = -1; // mark end
+      elts[j]->neighId[i].ind = -1; // mark end
+    }
+    /* For the most naive algorithm the case where the element itself is one must also be considered.
+       Thomas says this can later be optimized out. */
+    elts[j]->val = 1;
+    elts[j]->grad = gradient(*(elts[j]), neighbours);
+    elts[j]->val = 0;
+  }
 }