Ignore:
Timestamp:
Jun 3, 2009, 1:16:10 PM (16 years ago)
Author:
ole
Message:

Added another memset in quantity_ext.c and fixed some indentation issues to
comply with the Emacs c-style.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r6654 r7154  
    3030                        double* b){
    3131
    32         int i, k, k0, k1, k2, index3;
    33         double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
    34 
    35 
    36         for (k=0; k<N; k++) {
    37                 index3 = 3*k;
    38 
    39                 if (number_of_boundaries[k] < 2) {
    40                         // Two or three true neighbours
    41 
    42                         // Get indices of neighbours (or self when used as surrogate)
    43                         // k0, k1, k2 = surrogate_neighbours[k,:]
    44 
    45                         k0 = surrogate_neighbours[index3 + 0];
    46                         k1 = surrogate_neighbours[index3 + 1];
    47                         k2 = surrogate_neighbours[index3 + 2];
    48 
    49 
    50                         if (k0 == k1 || k1 == k2) return -1;
    51 
    52                         // Get data
    53                         q0 = centroid_values[k0];
    54                         q1 = centroid_values[k1];
    55                         q2 = centroid_values[k2];
    56 
    57                         x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    58                         x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    59                         x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
    60 
    61                         // Gradient
    62                         _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
    63 
    64                 } else if (number_of_boundaries[k] == 2) {
    65                         // One true neighbour
    66 
    67                         // Get index of the one neighbour
    68                         i=0; k0 = k;
    69                         while (i<3 && k0==k) {
    70                                 k0 = surrogate_neighbours[index3 + i];
    71                                 i++;
    72                         }
    73                         if (k0 == k) return -1;
    74 
    75                         k1 = k; //self
    76 
    77                         // Get data
    78                         q0 = centroid_values[k0];
    79                         q1 = centroid_values[k1];
    80 
    81                         x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    82                         x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    83 
    84                         // Two point gradient
    85                         _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
    86 
    87                 }
    88                 //    else:
    89                 //        #No true neighbours -
    90                 //        #Fall back to first order scheme
    91         }
    92         return 0;
     32  int i, k, k0, k1, k2, index3;
     33  double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
     34
     35
     36  for (k=0; k<N; k++) {
     37    index3 = 3*k;
     38   
     39    if (number_of_boundaries[k] < 2) {
     40      // Two or three true neighbours
     41     
     42      // Get indices of neighbours (or self when used as surrogate)
     43      // k0, k1, k2 = surrogate_neighbours[k,:]
     44     
     45      k0 = surrogate_neighbours[index3 + 0];
     46      k1 = surrogate_neighbours[index3 + 1];
     47      k2 = surrogate_neighbours[index3 + 2];
     48
     49     
     50      if (k0 == k1 || k1 == k2) return -1;
     51     
     52      // Get data
     53      q0 = centroid_values[k0];
     54      q1 = centroid_values[k1];
     55      q2 = centroid_values[k2];
     56
     57      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     58      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     59      x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
     60
     61      // Gradient
     62      _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
     63
     64    } else if (number_of_boundaries[k] == 2) {
     65      // One true neighbour
     66
     67      // Get index of the one neighbour
     68      i=0; k0 = k;
     69      while (i<3 && k0==k) {
     70        k0 = surrogate_neighbours[index3 + i];
     71        i++;
     72      }
     73      if (k0 == k) return -1;
     74     
     75      k1 = k; //self
     76
     77      // Get data
     78      q0 = centroid_values[k0];
     79      q1 = centroid_values[k1];
     80     
     81      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     82      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     83     
     84      // Two point gradient
     85      _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
     86     
     87    }
     88    //    else:
     89    //        #No true neighbours -
     90    //        #Fall back to first order scheme
     91  }
     92  return 0;
    9393}
    9494
    9595
    9696int _extrapolate_from_gradient(int N,
    97                 double* centroids,
    98                 double* centroid_values,
    99                 double* vertex_coordinates,
    100                 double* vertex_values,
    101                 double* edge_values,
    102                 double* a,
    103                 double* b) {
    104 
    105         int k, k2, k3, k6;
    106         double x, y, x0, y0, x1, y1, x2, y2;
    107 
    108         for (k=0; k<N; k++){
    109                 k6 = 6*k;
    110                 k3 = 3*k;
    111                 k2 = 2*k;
    112 
    113                 // Centroid coordinates
    114                 x = centroids[k2]; y = centroids[k2+1];
    115 
    116                 // vertex coordinates
    117                 // x0, y0, x1, y1, x2, y2 = X[k,:]
    118                 x0 = vertex_coordinates[k6 + 0];
    119                 y0 = vertex_coordinates[k6 + 1];
    120                 x1 = vertex_coordinates[k6 + 2];
    121                 y1 = vertex_coordinates[k6 + 3];
    122                 x2 = vertex_coordinates[k6 + 4];
    123                 y2 = vertex_coordinates[k6 + 5];
    124 
    125                 // Extrapolate to Vertices
    126                 vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
    127                 vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
    128                 vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
    129 
    130                 // Extrapolate to Edges (midpoints)
    131                 edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
    132                 edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
    133                 edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
    134 
    135         }
    136         return 0;
     97                              double* centroids,
     98                              double* centroid_values,
     99                              double* vertex_coordinates,
     100                              double* vertex_values,
     101                              double* edge_values,
     102                              double* a,
     103                              double* b) {
     104
     105  int k, k2, k3, k6;
     106  double x, y, x0, y0, x1, y1, x2, y2;
     107 
     108  for (k=0; k<N; k++){
     109    k6 = 6*k;
     110    k3 = 3*k;
     111    k2 = 2*k;
     112   
     113    // Centroid coordinates
     114    x = centroids[k2]; y = centroids[k2+1];
     115   
     116    // vertex coordinates
     117    // x0, y0, x1, y1, x2, y2 = X[k,:]
     118    x0 = vertex_coordinates[k6 + 0];
     119    y0 = vertex_coordinates[k6 + 1];
     120    x1 = vertex_coordinates[k6 + 2];
     121    y1 = vertex_coordinates[k6 + 3];
     122    x2 = vertex_coordinates[k6 + 4];
     123    y2 = vertex_coordinates[k6 + 5];
     124   
     125    // Extrapolate to Vertices
     126    vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
     127    vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
     128    vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
     129   
     130    // Extrapolate to Edges (midpoints)
     131    edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
     132    edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
     133    edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
     134   
     135  }
     136  return 0;
    137137}
    138138
    … …  
    149149                                         double* y_gradient) {
    150150
    151         int i, k, k2, k3, k6;
    152         double x, y, x0, y0, x1, y1, x2, y2;
    153         long n;
    154         double qmin, qmax, qc;
    155         double qn[3];
    156         double dq, dqa[3], r;
    157 
    158         for (k=0; k<N; k++){
    159                 k6 = 6*k;
    160                 k3 = 3*k;
    161                 k2 = 2*k;
    162 
    163                 // Centroid coordinates
    164                 x = centroids[k2+0];
    165                 y = centroids[k2+1];
    166 
    167                 // vertex coordinates
    168                 // x0, y0, x1, y1, x2, y2 = X[k,:]
    169                 x0 = vertex_coordinates[k6 + 0];
    170                 y0 = vertex_coordinates[k6 + 1];
    171                 x1 = vertex_coordinates[k6 + 2];
    172                 y1 = vertex_coordinates[k6 + 3];
    173                 x2 = vertex_coordinates[k6 + 4];
    174                 y2 = vertex_coordinates[k6 + 5];
    175 
    176                 // Extrapolate to Vertices
    177                 vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
    178                 vertex_values[k3+1] = centroid_values[k] + x_gradient[k]*(x1-x) + y_gradient[k]*(y1-y);
    179                 vertex_values[k3+2] = centroid_values[k] + x_gradient[k]*(x2-x) + y_gradient[k]*(y2-y);
    180 
    181                 // Extrapolate to Edges (midpoints)
    182                 edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
    183                 edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
    184                 edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
    185         }
    186 
    187 
    188 
    189         for (k=0; k<N; k++){
    190                 k6 = 6*k;
    191                 k3 = 3*k;
    192                 k2 = 2*k;
    193 
    194 
    195                 qc = centroid_values[k];
    196                
    197                 qmin = qc;
    198                 qmax = qc;
    199                
    200                 for (i=0; i<3; i++) {
    201                     n = neighbours[k3+i];
    202                     if (n < 0) {
    203                         qn[i] = qc;
    204                     } else {
    205                         qn[i] = centroid_values[n];
    206                     }
    207 
    208                     qmin = min(qmin, qn[i]);
    209                     qmax = max(qmax, qn[i]);
    210                 }
    211 
    212                 //qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
    213                 //qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
    214 
    215 /*              for (i=0; i<3; i++) { */
    216 /*                  n = neighbours[k3+i]; */
    217 /*                  if (n < 0) { */
    218 /*                      qn[i] = qc; */
    219 /*                      qmin[i] = qtmin; */
    220 /*                      qmax[i] = qtmax; */
    221 /*                  }  */
    222 /*              } */
    223                        
    224                 phi[k] = 1.0;
    225 
    226                 for (i=0; i<3; i++) {
    227                     dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
    228                     dqa[i] = dq;                      //Save dq for use in updating vertex values
    229 
    230                     r = 1.0;
    231      
    232                     if (dq > 0.0) r = (qmax - qc)/dq;
    233                     if (dq < 0.0) r = (qmin - qc)/dq;
    234                    
    235                     phi[k] = min( min(r*beta, 1.0), phi[k]);
    236                    
    237                 }
    238 
    239 
    240 
    241                 //Update gradient, edge and vertex values using phi limiter
    242                 x_gradient[k] = x_gradient[k]*phi[k];
    243                 y_gradient[k] = y_gradient[k]*phi[k];
    244 
    245                 edge_values[k3+0] = qc + phi[k]*dqa[0];
    246                 edge_values[k3+1] = qc + phi[k]*dqa[1];
    247                 edge_values[k3+2] = qc + phi[k]*dqa[2];
    248 
    249 
    250                 vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
    251                 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
    252                 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
    253                
    254 
    255         }
    256 
    257         return 0;
    258 
     151  int i, k, k2, k3, k6;
     152  double x, y, x0, y0, x1, y1, x2, y2;
     153  long n;
     154  double qmin, qmax, qc;
     155  double qn[3];
     156  double dq, dqa[3], r;
     157
     158  for (k=0; k<N; k++){
     159    k6 = 6*k;
     160    k3 = 3*k;
     161    k2 = 2*k;
     162
     163    // Centroid coordinates
     164    x = centroids[k2+0];
     165    y = centroids[k2+1];
     166
     167    // vertex coordinates
     168    // x0, y0, x1, y1, x2, y2 = X[k,:]
     169    x0 = vertex_coordinates[k6 + 0];
     170    y0 = vertex_coordinates[k6 + 1];
     171    x1 = vertex_coordinates[k6 + 2];
     172    y1 = vertex_coordinates[k6 + 3];
     173    x2 = vertex_coordinates[k6 + 4];
     174    y2 = vertex_coordinates[k6 + 5];
     175   
     176    // Extrapolate to Vertices
     177    vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
     178    vertex_values[k3+1] = centroid_values[k] + x_gradient[k]*(x1-x) + y_gradient[k]*(y1-y);
     179    vertex_values[k3+2] = centroid_values[k] + x_gradient[k]*(x2-x) + y_gradient[k]*(y2-y);
     180   
     181    // Extrapolate to Edges (midpoints)
     182    edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
     183    edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
     184    edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
     185  }
     186 
     187 
     188 
     189  for (k=0; k<N; k++){
     190    k6 = 6*k;
     191    k3 = 3*k;
     192    k2 = 2*k;
     193   
     194   
     195    qc = centroid_values[k];
     196   
     197    qmin = qc;
     198    qmax = qc;
     199   
     200    for (i=0; i<3; i++) {
     201      n = neighbours[k3+i];
     202      if (n < 0) {
     203        qn[i] = qc;
     204      } else {
     205        qn[i] = centroid_values[n];
     206      }
     207     
     208      qmin = min(qmin, qn[i]);
     209      qmax = max(qmax, qn[i]);
     210    }
     211   
     212    //qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
     213    //qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
     214   
     215    /*          for (i=0; i<3; i++) { */
     216    /*              n = neighbours[k3+i]; */
     217    /*              if (n < 0) { */
     218    /*                  qn[i] = qc; */
     219    /*                  qmin[i] = qtmin; */
     220    /*                  qmax[i] = qtmax; */
     221    /*              }  */
     222    /*          } */
     223   
     224    phi[k] = 1.0;
     225   
     226    for (i=0; i<3; i++) {
     227      dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
     228      dqa[i] = dq;                      //Save dq for use in updating vertex values
     229     
     230      r = 1.0;
     231     
     232      if (dq > 0.0) r = (qmax - qc)/dq;
     233      if (dq < 0.0) r = (qmin - qc)/dq;
     234     
     235      phi[k] = min( min(r*beta, 1.0), phi[k]);
     236     
     237    }
     238   
     239   
     240   
     241    //Update gradient, edge and vertex values using phi limiter
     242    x_gradient[k] = x_gradient[k]*phi[k];
     243    y_gradient[k] = y_gradient[k]*phi[k];
     244   
     245    edge_values[k3+0] = qc + phi[k]*dqa[0];
     246    edge_values[k3+1] = qc + phi[k]*dqa[1];
     247    edge_values[k3+2] = qc + phi[k]*dqa[2];
     248   
     249   
     250    vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
     251    vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
     252    vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     253   
     254   
     255  }
     256 
     257  return 0;
     258 
    259259}
    260260
    … …  
    263263
    264264int _limit_vertices_by_all_neighbours(int N, double beta,
    265                      double* centroid_values,
    266                      double* vertex_values,
    267                      double* edge_values,
    268                      long*   neighbours,
    269                      double* x_gradient,
    270                      double* y_gradient) {
    271 
    272 
    273         int i, k, k2, k3, k6;
    274         long n;
    275         double qmin, qmax, qn, qc;
    276         double dq, dqa[3], phi, r;
    277 
    278         for (k=0; k<N; k++){
    279                 k6 = 6*k;
    280                 k3 = 3*k;
    281                 k2 = 2*k;
    282 
    283                 qc = centroid_values[k];
    284                 qmin = qc;
    285                 qmax = qc;
    286 
    287                 for (i=0; i<3; i++) {
    288                     n = neighbours[k3+i];
    289                     if (n >= 0) {
    290                         qn = centroid_values[n]; //Neighbour's centroid value
    291 
    292                         qmin = min(qmin, qn);
    293                         qmax = max(qmax, qn);
    294                     }
    295                 }
    296 
    297                 phi = 1.0;
    298                 for (i=0; i<3; i++) {   
    299                     r = 1.0;
    300      
    301                     dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
    302                     dqa[i] = dq;                      //Save dq for use in updating vertex values
    303      
    304                     if (dq > 0.0) r = (qmax - qc)/dq;
    305                     if (dq < 0.0) r = (qmin - qc)/dq;     
     265                                      double* centroid_values,
     266                                      double* vertex_values,
     267                                      double* edge_values,
     268                                      long*   neighbours,
     269                                      double* x_gradient,
     270                                      double* y_gradient) {
    306271 
    307                    
    308                     phi = min( min(r*beta, 1.0), phi);   
    309                 }
    310 
    311                 //Update gradient, vertex and edge values using phi limiter
    312                 x_gradient[k] = x_gradient[k]*phi;
    313                 y_gradient[k] = y_gradient[k]*phi;
    314    
    315                 vertex_values[k3+0] = qc + phi*dqa[0];
    316                 vertex_values[k3+1] = qc + phi*dqa[1];
    317                 vertex_values[k3+2] = qc + phi*dqa[2];
    318                
    319                 edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
    320                 edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
    321                 edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
    322 
    323         }
    324 
    325         return 0;
     272 
     273  int i, k, k2, k3, k6;
     274  long n;
     275  double qmin, qmax, qn, qc;
     276  double dq, dqa[3], phi, r;
     277 
     278  for (k=0; k<N; k++){
     279    k6 = 6*k;
     280    k3 = 3*k;
     281    k2 = 2*k;
     282   
     283    qc = centroid_values[k];
     284    qmin = qc;
     285    qmax = qc;
     286   
     287    for (i=0; i<3; i++) {
     288      n = neighbours[k3+i];
     289      if (n >= 0) {
     290        qn = centroid_values[n]; //Neighbour's centroid value
     291       
     292        qmin = min(qmin, qn);
     293        qmax = max(qmax, qn);
     294      }
     295    }
     296   
     297    phi = 1.0;
     298    for (i=0; i<3; i++) {   
     299      r = 1.0;
     300     
     301      dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
     302      dqa[i] = dq;                      //Save dq for use in updating vertex values
     303     
     304      if (dq > 0.0) r = (qmax - qc)/dq;
     305      if (dq < 0.0) r = (qmin - qc)/dq;     
     306     
     307     
     308      phi = min( min(r*beta, 1.0), phi);   
     309    }
     310   
     311    //Update gradient, vertex and edge values using phi limiter
     312    x_gradient[k] = x_gradient[k]*phi;
     313    y_gradient[k] = y_gradient[k]*phi;
     314   
     315    vertex_values[k3+0] = qc + phi*dqa[0];
     316    vertex_values[k3+1] = qc + phi*dqa[1];
     317    vertex_values[k3+2] = qc + phi*dqa[2];
     318   
     319    edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
     320    edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
     321    edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
     322   
     323  }
     324 
     325  return 0;
    326326}
    327327
    … …  
    337337                                   double* y_gradient) {
    338338
    339         int i, k, k2, k3, k6;
    340         long n;
    341         double qmin, qmax, qn, qc;
    342         double dq, dqa[3], phi, r;
    343 
    344         for (k=0; k<N; k++){
    345                 k6 = 6*k;
    346                 k3 = 3*k;
    347                 k2 = 2*k;
    348 
    349                 qc = centroid_values[k];
    350                 qmin = qc;
    351                 qmax = qc;
    352 
    353                 for (i=0; i<3; i++) {
    354                     n = neighbours[k3+i];
    355                     if (n >= 0) {
    356                         qn = centroid_values[n]; //Neighbour's centroid value
    357 
    358                         qmin = min(qmin, qn);
    359                         qmax = max(qmax, qn);
    360                     }
    361                 }
    362 
    363                 phi = 1.0;
    364                 for (i=0; i<3; i++) {   
    365                     r = 1.0;
    366      
    367                     dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
    368                     dqa[i] = dq;                      //Save dq for use in updating vertex values
    369      
    370                     if (dq > 0.0) r = (qmax - qc)/dq;
    371                     if (dq < 0.0) r = (qmin - qc)/dq;     
     339  int i, k, k2, k3, k6;
     340  long n;
     341  double qmin, qmax, qn, qc;
     342  double dq, dqa[3], phi, r;
    372343 
    373                    
    374                     phi = min( min(r*beta, 1.0), phi);   
    375                 }
    376    
    377                 //Update gradient, vertex and edge values using phi limiter
    378                 x_gradient[k] = x_gradient[k]*phi;
    379                 y_gradient[k] = y_gradient[k]*phi;
    380 
    381                 edge_values[k3+0] = qc + phi*dqa[0];
    382                 edge_values[k3+1] = qc + phi*dqa[1];
    383                 edge_values[k3+2] = qc + phi*dqa[2];
    384                
    385                 vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
    386                 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
    387                 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
    388 
    389         }
    390 
    391         return 0;
     344  for (k=0; k<N; k++){
     345    k6 = 6*k;
     346    k3 = 3*k;
     347    k2 = 2*k;
     348   
     349    qc = centroid_values[k];
     350    qmin = qc;
     351    qmax = qc;
     352   
     353    for (i=0; i<3; i++) {
     354      n = neighbours[k3+i];
     355      if (n >= 0) {
     356        qn = centroid_values[n]; //Neighbour's centroid value
     357       
     358        qmin = min(qmin, qn);
     359        qmax = max(qmax, qn);
     360      }
     361    }
     362   
     363    phi = 1.0;
     364    for (i=0; i<3; i++) {   
     365      r = 1.0;
     366     
     367      dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
     368      dqa[i] = dq;                      //Save dq for use in updating vertex values
     369     
     370      if (dq > 0.0) r = (qmax - qc)/dq;
     371      if (dq < 0.0) r = (qmin - qc)/dq;     
     372     
     373     
     374      phi = min( min(r*beta, 1.0), phi);   
     375    }
     376   
     377    //Update gradient, vertex and edge values using phi limiter
     378    x_gradient[k] = x_gradient[k]*phi;
     379    y_gradient[k] = y_gradient[k]*phi;
     380   
     381    edge_values[k3+0] = qc + phi*dqa[0];
     382    edge_values[k3+1] = qc + phi*dqa[1];
     383    edge_values[k3+2] = qc + phi*dqa[2];
     384   
     385    vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
     386    vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
     387    vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     388   
     389  }
     390 
     391  return 0;
    392392}
    393393
    … …  
    730730
    731731
    732         // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
    733         for (k=0;k<N;k++){
    734                 semi_implicit_update[k]=0.0;
    735         }
     732        // Reset semi_implicit_update here ready for next time step
     733        memset(semi_implicit_update, 0, N*sizeof(double));
    736734
    737735        return 0;
Note: See TracChangeset for help on using the changeset viewer.