Changeset 844 for XIOS/trunk/extern/remap/src/meshutil.cpp
- Timestamp:
- 04/27/16 17:23:48 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/trunk/extern/remap/src/meshutil.cpp
r688 r844 9 9 void cptEltGeom(Elt& elt, const Coord &pole) 10 10 { 11 12 13 14 15 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; 16 16 } 17 17 18 18 void cptAllEltsGeom(Elt *elt, int N, const Coord &pole) 19 19 { 20 21 20 for (int ne=0; ne<N; ne++) 21 cptEltGeom(elt[ne], pole); 22 22 } 23 23 … … 26 26 void update_baryc(Elt *elt, int N) 27 27 { 28 29 30 31 32 33 34 35 36 37 38 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 } 40 40 } 41 42 43 Coord 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 41 67 42 68 Coord gradient(Elt& elt, Elt **neighElts) 43 69 { 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] ; 51 73 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 ; 61 81 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); 63 114 } 115 116 117 64 118 65 119 void computeGradients(Elt **elts, int N) 66 120 { 67 68 69 121 122 for (int j = 0; j < N; j++) 123 elts[j]->val = 0; 70 124 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]]; 76 131 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 } 96 164 } 97 165
Note: See TracChangeset
for help on using the changeset viewer.