proteus  1.8.1
C/C++/Fortran libraries
timeIntegration.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <strings.h>
5 #include <assert.h>
6 #include "timeIntegration.h"
11 /*Pseudo-Transient Continuation TTE formula choice for dt*/
12 void psiTCtteDT(int nPoints,
13  double tau,
14  double dtn,
15  double dtnm1,
16  double * yn,
17  double * ypn,
18  double * ypnm1,
19  double * dtnp1)
20 {
21  int k;
22  double dy2,dtavgInv,maxdy2,dtOut,safety,eps;
23 
24  dtavgInv = 2.0/(dtn+dtnm1);
25  maxdy2 = 0.0;
26  safety = 0.9;
27  eps = 1.0e-8;
28  /*mwf debug
29  printf("psiTCtte nPoints=%d tau= %g dtn=%g dtnm1=%g \n",nPoints,tau,dtn,dtnm1);
30  */
31  for (k=0; k < nPoints; k++)
32  {
33  dy2 = ypn[k]-ypnm1[k];
34  maxdy2 = fmax(maxdy2,fabs(dy2)*0.5/(1.0+fabs(yn[k])));
35  /*mwf debug
36  printf("k=%d ypn= %g ypnm1= %g maxdy2= %g\n",k,ypn[k],ypnm1[k],maxdy2);
37  */
38  }
39  /*mwf debug
40  printf("psiTCtte maxdy2= %g\n",maxdy2);
41  */
42  dtOut = sqrt(tau/(maxdy2+eps))*safety;
43  *dtnp1 = dtOut;
44 }
45 
46 /***********************************************************************
47  RKDG limiting procedures, probably need a better spot
48  ***********************************************************************/
49 double msign(double a)
50 {
51  return (a < 0.0 ? -1.0 : 1.0);
52 }
53 double minmod2(double a, double b)
54 {
55  double aa = fabs(a), bb = fabs(b);
56  if (a*b < 0.0) return 0.0;
57  return aa <= bb ? a : b;
58 }
59 double musclmod(double dU, double dU1, double dU2)
60 {
61  double tmp = minmod2(dU1,dU2);
62  return minmod2(dU,2.0*tmp);
63 }
64 double minmod3(double dU, double dU1, double dU2)
65 {
66  double tmp = minmod2(dU1,dU2);
67  return minmod2(dU,tmp);
68 }
69 double mminmod2(double dU, double dU1, double M, int* tag)
70 {
71  if (fabs(dU) <= M)
72  {
73  *tag = 1;
74  return dU;
75  }
76  *tag = 0;
77  return minmod2(dU,dU1);
78 }
79 /*make sure doesn't conflict with postprocessing's version*/
80 int invertLocal3d(double A[3][3], double AI[3][3])
81 {
82  int i,j;
83  double detA,detAinv;
84  detA =
85  A[0][0]*(A[1][1]*A[2][2]-A[2][1]*A[1][2])-
86  A[0][1]*(A[1][0]*A[2][2]-A[2][0]*A[1][2])+
87  A[0][2]*(A[1][0]*A[2][1]-A[2][0]*A[1][1]);
88  /*mwf debug*/
89  if (fabs(detA) <= 0.0)
90  {
91  printf("WARNING invertLocal3d detA= %g A= \n",detA);
92  for (i=0; i < 3; i++)
93  {
94  for (j=0; j < 3; j++)
95  printf("%g ",A[i][j]);
96  printf("\n");
97  }
98  return -1;
99  }
100  assert(fabs(detA) > 0.0);
101  detAinv = 1.0/detA;
102  AI[0][0] = detAinv*(A[1][1]*A[2][2]-A[1][2]*A[2][1]);
103  AI[0][1] = detAinv*(A[0][2]*A[2][1]-A[0][1]*A[2][2]);
104  AI[0][2] = detAinv*(A[0][1]*A[1][2]-A[0][2]*A[1][1]);
105 
106  AI[1][0] = detAinv*(A[1][2]*A[2][0]-A[1][0]*A[2][2]);
107  AI[1][1] = detAinv*(A[0][0]*A[2][2]-A[0][2]*A[2][0]);
108  AI[1][2] = detAinv*(A[0][2]*A[1][0]-A[0][0]*A[1][2]);
109 
110  AI[2][0] = detAinv*(A[1][0]*A[2][1]-A[1][1]*A[2][0]);
111  AI[2][1] = detAinv*(A[0][1]*A[2][0]-A[0][0]*A[2][1]);
112  AI[2][2] = detAinv*(A[0][0]*A[1][1]-A[0][1]*A[1][0]);
113  return 0;
114 }
115 
116 
117 /* local computation ignores orientation */
118 void computeLocalGradient1d(double * x0, double * x1, double * dN0, double * dN1,
119  double u0, double u1, double * dU)
120 {
121  register double volume;
122  volume = fabs(x1[0]-x0[0]);
123  dN0[0]= -(x1[0]-x0[0]);/* -n0*area/volume/nSpace */
124  dN1[0]= -(x0[0]-x1[0]);/* -n1*area/volume/nSpace */
125  dU[0] = dN0[0]*u0 + dN1[0]*u1;
126 }
127 /* local computation ignores orientation */
128 void computeLocalGradient2d(double * x0, double * x1, double * x2, double * dN0,
129  double * dN1, double * dN2,
130  double u0, double u1, double u2, double * dU)
131 {
132  double detJ,dx10[2],dx20[2],JinvT[2][2];
133  /***************************************************
134  \hat{N}_0 = 1-\hat{x}-\hat{y}, \hat{N}_1 = \hat{x}
135  \hat{N}_2 = \hat{y}
136  ***************************************************/
137  const double dN0h[2] = {-1.,-1.};
138  const double dN1h[2] = { 1., 0.};
139  const double dN2h[2] = { 0., 1.};
140  dx10[0] = x1[0]-x0[0]; dx10[1] = x1[1]-x0[1];
141  dx20[0] = x2[0]-x0[0]; dx20[1] = x2[1]-x0[1];
142  detJ = dx10[0]*dx20[1]-dx10[1]*dx20[0];
143 
144  /***************************************************
145  J = |x1-x0 x2-x0| Jinv = |y2-y0 x0-x2|*1/det JinvT = |y2-y0 y0-y1|*1/det
146  |y1-y0 y2-y0| |y0-y1 x1-x0| |x0-x2 x1-x0|
147  *******************************/
148  JinvT[0][0] = dx20[1]/detJ; JinvT[0][1] =-dx10[1]/detJ;
149  JinvT[1][0] =-dx20[0]/detJ; JinvT[1][1] = dx10[0]/detJ;
150 
151  dN0[0] = JinvT[0][0]*dN0h[0]+JinvT[0][1]*dN0h[1];
152  dN0[1] = JinvT[1][0]*dN0h[0]+JinvT[1][1]*dN0h[1];
153 
154  dN1[0] = JinvT[0][0]*dN1h[0]+JinvT[0][1]*dN1h[1];
155  dN1[1] = JinvT[1][0]*dN1h[0]+JinvT[1][1]*dN1h[1];
156 
157  dN2[0] = JinvT[0][0]*dN2h[0]+JinvT[0][1]*dN2h[1];
158  dN2[1] = JinvT[1][0]*dN2h[0]+JinvT[1][1]*dN2h[1];
159 
160  dU[0] = u0*dN0[0] + u1*dN1[0] + u2*dN2[0];
161  dU[1] = u0*dN0[1] + u1*dN1[1] + u2*dN2[1];
162 }
163 
164 void computeLocalGradient3d(double * x0, double * x1, double * x2, double * x3,
165  double * dN0, double * dN1, double * dN2, double * dN3,
166  double u0, double u1, double u2, double u3, double * dU)
167 {
168  double dx[3][3],Jinv[3][3],JinvT[3][3];
169  int gradientFailed = 0;
170  /***************************************************
171  \hat{N}_0 = 1-\hat{x}-\hat{y}-\hat{z}, \hat{N}_1 = \hat{x}
172  \hat{N}_2 = \hat{y}, \hat{N}_3 = \hat{z}
173  ***************************************************/
174  const double dN0h[3] = {-1.,-1.,-1.};
175  const double dN1h[3] = { 1., 0.,0.};
176  const double dN2h[3] = { 0., 1.,0.};
177  const double dN3h[3] = { 0., 0.,1.};
178  dx[0][0] = x1[0]-x0[0]; dx[1][0] = x1[1]-x0[1]; dx[2][0] = x1[2]-x0[2];
179  dx[0][1] = x2[0]-x0[0]; dx[1][1] = x2[1]-x0[1]; dx[2][1] = x2[2]-x0[2];
180  dx[0][2] = x3[0]-x0[0]; dx[1][2] = x3[1]-x0[1]; dx[2][2] = x3[2]-x0[2];
181 
182  /*find routine for 3x3 inverse and jacobian*/
183  gradientFailed = invertLocal3d(dx,Jinv);
184  if (gradientFailed != 0)
185  {
186  dU[0]=0.0; dU[1]=0.0; dU[2]=0.0;
187  }
188  JinvT[0][0] = Jinv[0][0]; JinvT[0][1] = Jinv[1][0]; JinvT[0][2] = Jinv[2][0];
189  JinvT[1][0] = Jinv[0][1]; JinvT[1][1] = Jinv[1][1]; JinvT[1][2] = Jinv[2][1];
190  JinvT[2][0] = Jinv[0][2]; JinvT[2][1] = Jinv[1][2]; JinvT[2][2] = Jinv[2][2];
191 
192 
193 
194  dN0[0] = JinvT[0][0]*dN0h[0]+JinvT[0][1]*dN0h[1]+JinvT[0][2]*dN0h[2];
195  dN0[1] = JinvT[1][0]*dN0h[0]+JinvT[1][1]*dN0h[1]+JinvT[1][2]*dN0h[2];
196  dN0[2] = JinvT[2][0]*dN0h[0]+JinvT[2][1]*dN0h[1]+JinvT[2][2]*dN0h[2];
197 
198  dN1[0] = JinvT[0][0]*dN1h[0]+JinvT[0][1]*dN1h[1]+JinvT[0][2]*dN1h[2];
199  dN1[1] = JinvT[1][0]*dN1h[0]+JinvT[1][1]*dN1h[1]+JinvT[1][2]*dN1h[2];
200  dN1[2] = JinvT[2][0]*dN1h[0]+JinvT[2][1]*dN1h[1]+JinvT[2][2]*dN1h[2];
201 
202  dN2[0] = JinvT[0][0]*dN2h[0]+JinvT[0][1]*dN2h[1]+JinvT[0][2]*dN2h[2];
203  dN2[1] = JinvT[1][0]*dN2h[0]+JinvT[1][1]*dN2h[1]+JinvT[1][2]*dN2h[2];
204  dN2[2] = JinvT[2][0]*dN2h[0]+JinvT[2][1]*dN2h[1]+JinvT[2][2]*dN2h[2];
205 
206  dN3[0] = JinvT[0][0]*dN3h[0]+JinvT[0][1]*dN3h[1]+JinvT[0][2]*dN3h[2];
207  dN3[1] = JinvT[1][0]*dN3h[0]+JinvT[1][1]*dN3h[1]+JinvT[1][2]*dN3h[2];
208  dN3[2] = JinvT[2][0]*dN3h[0]+JinvT[2][1]*dN3h[1]+JinvT[2][2]*dN3h[2];
209 
210  dU[0] = u0*dN0[0] + u1*dN1[0] + u2*dN2[0] + u3*dN3[0];
211  dU[1] = u0*dN0[1] + u1*dN1[1] + u2*dN2[1] + u3*dN3[1];
212  dU[2] = u0*dN0[2] + u1*dN1[2] + u2*dN2[2] + u3*dN3[2];
213 }
214 void applyDGlimitingP1Lagrange1d(int limiterFlag,
215  int nElements_global,
216  int nNodes_element,
217  int nElementBoundaries_element,
218  int nDOF_element,
219  int * elementNodesArray,
220  int * elementNeighborsArray,
221  double * nodeArray,
222  double * elementBarycentersArray,
223  int * l2g,
224  int * tag,
225  double * Uin,
226  double * Uout)
227 {
228  /***********************************************************************
229  Apply basic DG limiting procedure in 1d assuming input/output fem functions
230  have a P1 Lagrange basis representation over each element
231 
232  Assumes that local node numbering agrees with local dof numbering
233  Assumes that local elementBoundary i is across from node i
234 
235  Input:
236  limitingFlag --- which limiter to use on slopes
237  0 minmod3(\pd{u}{x}_i,\Delta_x u_i+1/2,\Delta_x u_i-1/2)
238  1 muscl(\pd{u}{x}_i,\Delta_x u_i+1/2,\Delta_x u_i-1/2) (doubles slopes)
239  2 minmod2(\Delta_x u_i+1/2,\Delta_x u_i-1/2)
240  elementNodesArray[eN,:] --- global node numbers on element eN
241  elementBoundariesArray[eN,:] --- global element boundaries on element eN
242  elementNeighborsArray[eN,i] --- global element number across local elementBoundary i
243  -1 ---> exterior boundary
244  nodeArray[nN,:] --- physical coordinates of node nN (always 3)
245  elementBarycentersArray[eN,:]--- physical coordiantes of barycenter for element eN (always 3)
246  l2g[eN,:] --- global degree of freedom numbers for local unknowns
247  Uin --- input P1 degrees of freedom values
248 
249 
250  Output:
251  Uout --- output P1 degrees of freedom values
252  tag[eN] --- 1 if original P1 representation on element eN ok (not limited)
253 
254  ***********************************************************************/
255 
256  int eN;
257  int eN0,eN1,nN0,nN1,ndof0,ndof1;
258 
259  double dUlim,dU,dU0,dU1,Ubar,Ubar0,Ubar1;
260  double xbar,xbar0,xbar1,dx;
261 
262  for (eN=0; eN < nElements_global; eN++)
263  {
264  ndof0 = l2g[eN*nDOF_element + 0]; ndof1 = l2g[eN*nDOF_element + 1];
265  nN0 = elementNodesArray[eN*nNodes_element + 0];
266  nN1 = elementNodesArray[eN*nNodes_element + 1];
267  dx = nodeArray[nN1*3 + 0] - nodeArray[nN0*3 + 0];
268  xbar = elementBarycentersArray[eN*3 + 0];
269 
270  dU = (Uin[ndof1]-Uin[ndof0])/dx;
271  Ubar = 0.5*(Uin[ndof1]+Uin[ndof0]);
272 
273  /*** neighbor across from node 0 and adjacent to elementBoundary 0***/
274  /* in case on physical boundary*/
275  xbar0 = nodeArray[nN1*3 + 0]; Ubar0 = Uin[ndof1];
276  eN0 = elementNeighborsArray[eN*nElementBoundaries_element + 0];
277  if (eN0 >= 0)
278  {
279  xbar0 = elementBarycentersArray[eN0*3 + 0];
280  Ubar0 = 0.5*(Uin[l2g[eN0*nDOF_element + 0]] + Uin[l2g[eN0*nDOF_element + 1]]);
281  }
282  dU0 = (Ubar0-Ubar)/(xbar0-xbar);
283 
284  /*** neighbor across from node 1 and adjacent to elementBoundary 1***/
285  /* in case on physical boundary*/
286  xbar1 = nodeArray[nN0*3 + 0]; Ubar1 = Uin[ndof0];
287  eN1 = elementNeighborsArray[eN*nElementBoundaries_element + 1];
288  if (eN1 >= 0)
289  {
290  xbar1 = elementBarycentersArray[eN1*3 + 0];
291  Ubar1 = 0.5*(Uin[l2g[eN1*nDOF_element + 0]] + Uin[l2g[eN1*nDOF_element + 1]]);
292  }
293  dU1 = (Ubar1-Ubar)/(xbar1-xbar);
294  /*dU's have dx's */
295  if (limiterFlag == 2)
296  dUlim = minmod2(dU0,dU1);
297  else if (limiterFlag == 1)
298  dUlim = musclmod(dU,dU0,dU1);
299  else
300  dUlim = minmod3(dU,dU0,dU1);
301  tag[eN] = fabs(dUlim-dU) < 1.0e-6; /*need better tolerance*/
302 
303  /*back out nodal dofs from average and slope*/
304  Uout[ndof0] = Ubar + (nodeArray[nN0*3 + 0]-xbar)*dUlim;
305  Uout[ndof1] = Ubar + (nodeArray[nN1*3 + 0]-xbar)*dUlim;
306  }/*eN*/
307 
308 }
309 
310 void applyDGlimitingP1Lagrange1d_withVacuumTol(int enforcePositivity,
311  double vacuumTol,
312  int nElements_global,
313  int nNodes_element,
314  int nElementBoundaries_element,
315  int nDOF_element,
316  int * elementNodesArray,
317  int * elementNeighborsArray,
318  double * nodeArray,
319  double * elementBarycentersArray,
320  int * l2g,
321  int * tag,
322  double * Uin,
323  double * Uout)
324 {
325  /***********************************************************************
326  Apply basic DG limiting procedure in 1d assuming input/output fem functions
327  have a P1 Lagrange basis representation over each element
328 
329  Assumes that local node numbering agrees with local dof numbering
330  Assumes that local elementBoundary i is across from node i
331 
332  Input:
333  vacuumTol --- if |\bar{u}| < tol use minmod2 else muscl
334  minmod3(\pd{u}{x}_i,\Delta_x u_i+1/2,\Delta_x u_i-1/2)
335  muscl(\pd{u}{x}_i,\Delta_x u_i+1/2,\Delta_x u_i-1/2) (doubles slopes)
336  minmod2(\Delta_x u_i+1/2,\Delta_x u_i-1/2)
337  elementNodesArray[eN,:] --- global node numbers on element eN
338  elementBoundariesArray[eN,:] --- global element boundaries on element eN
339  elementNeighborsArray[eN,i] --- global element number across local elementBoundary i
340  -1 ---> exterior boundary
341  nodeArray[nN,:] --- physical coordinates of node nN (always 3)
342  elementBarycentersArray[eN,:]--- physical coordiantes of barycenter for element eN (always 3)
343  l2g[eN,:] --- global degree of freedom numbers for local unknowns
344  Uin --- input P1 degrees of freedom values
345 
346 
347  Output:
348  Uout --- output P1 degrees of freedom values
349  tag[eN] --- 1 if original P1 representation on element eN ok (not limited)
350 
351  ***********************************************************************/
352 
353  int eN;
354  int eN0,eN1,nN0,nN1,ndof0,ndof1;
355 
356  double dUlim,dU,dU0,dU1,Ubar,Ubar0,Ubar1,Umin,UminA;
357  double xbar,xbar0,xbar1,dx;
358 
359  for (eN=0; eN < nElements_global; eN++)
360  {
361  ndof0 = l2g[eN*nDOF_element + 0]; ndof1 = l2g[eN*nDOF_element + 1];
362  nN0 = elementNodesArray[eN*nNodes_element + 0];
363  nN1 = elementNodesArray[eN*nNodes_element + 1];
364  dx = nodeArray[nN1*3 + 0] - nodeArray[nN0*3 + 0];
365  xbar = elementBarycentersArray[eN*3 + 0];
366 
367  dU = (Uin[ndof1]-Uin[ndof0])/dx;
368  Ubar = 0.5*(Uin[ndof1]+Uin[ndof0]);
369  Umin = fmin(Uin[ndof1],Uin[ndof0]);/*or Ubar?*/
370  UminA = fmin(fabs(Uin[ndof1]),fabs(Uin[ndof0])); /*or |Ubar| ?*/
371 
372  /*** neighbor across from node 0 and adjacent to elementBoundary 0***/
373  /* in case on physical boundary*/
374  xbar0 = nodeArray[nN1*3 + 0]; Ubar0 = Uin[ndof1];
375  eN0 = elementNeighborsArray[eN*nElementBoundaries_element + 0];
376  if (eN0 >= 0)
377  {
378  xbar0 = elementBarycentersArray[eN0*3 + 0];
379  Ubar0 = 0.5*(Uin[l2g[eN0*nDOF_element + 0]] + Uin[l2g[eN0*nDOF_element + 1]]);
380  }
381  dU0 = (Ubar0-Ubar)/(xbar0-xbar);
382 
383  /*** neighbor across from node 1 and adjacent to elementBoundary 1***/
384  /* in case on physical boundary*/
385  xbar1 = nodeArray[nN0*3 + 0]; Ubar1 = Uin[ndof0];
386  eN1 = elementNeighborsArray[eN*nElementBoundaries_element + 1];
387  if (eN1 >= 0)
388  {
389  xbar1 = elementBarycentersArray[eN1*3 + 0];
390  Ubar1 = 0.5*(Uin[l2g[eN1*nDOF_element + 0]] + Uin[l2g[eN1*nDOF_element + 1]]);
391  }
392  dU1 = (Ubar1-Ubar)/(xbar1-xbar);
393  /*dU's have dx's */
394  if (enforcePositivity && Umin < vacuumTol)
395  dUlim = minmod3(dU,dU0,dU1);
396  else if (UminA < vacuumTol)
397  dUlim = minmod3(dU,dU0,dU1);
398  else
399  dUlim = musclmod(dU,dU0,dU1);
400  tag[eN] = fabs(dUlim-dU) < 1.0e-6; /*need better tolerance*/
401 
402  /*back out nodal dofs from average and slope*/
403  Uout[ndof0] = Ubar + (nodeArray[nN0*3 + 0]-xbar)*dUlim;
404  Uout[ndof1] = Ubar + (nodeArray[nN1*3 + 0]-xbar)*dUlim;
405  }/*eN*/
406 
407 }
408 
409 void computeElementNeighborShapeGradients(int nElements_global,
410  int nElementBoundaries_element,
411  int nSpace,
412  const int * elementBoundariesArray,
413  const int * elementNeighborsArray,
414  double * elementBarycentersArray,
415  double * elementBoundaryBarycentersArray,
416  double * elementNeighborShapeGradients)
417 {
418  /***********************************************************************
419  compute gradient of linear shape functions (aka barycentric coordinates)
420  for the simplices formed by connecting element neighbor barycenters
421 
422  store these in
423  elementNeighborShapeGradients[eN,i,j,k] --->
424  kth coordinate of jth barycentric coordinate for ith simplex
425  formed from element eN and its nd+1 neighbors
426 
427  ith neighbor consists of element eN, neighbor i and neighbor i+1 % (nSpace+1)
428 
429 
430  ***********************************************************************/
431  int eN,ebN,eN_opposite,ebN_global,ebNplus1, ebNplus2,I;
432  double u[4] = {1.0,1.0,1.0,1.0};
433  double dU[3]= {0.0,0.0,0.0};
434  /*local shape gradients (nSpace+1)*/
435  double dN[4][3] = {{0.0,0.0,0.0},
436  {0.0,0.0,0.0},
437  {0.0,0.0,0.0},
438  {0.0,0.0,0.0}};
439  /*neighboring barycenters (nSpace+2)*/
440  double xn[5][3] = {{0.0,0.0,0.0},
441  {0.0,0.0,0.0},
442  {0.0,0.0,0.0},
443  {0.0,0.0,0.0},
444  {0.0,0.0,0.0}};
445  if (nSpace == 1)
446  {
447  for (eN = 0; eN < nElements_global; eN++)
448  {
449  xn[0][0] = elementBarycentersArray[eN*3 + 0];
450  xn[0][1] = elementBarycentersArray[eN*3 + 1];
451  xn[0][2] = elementBarycentersArray[eN*3 + 2];
452 
453  /*go through boundaries on this element and form simpleces (x0,x1), (x0,x2) etc*/
454  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
455  {
456  eN_opposite = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
457  if (eN_opposite < 0)/*outside domain*/
458  {
459  ebN_global = elementBoundariesArray[eN*nElementBoundaries_element + ebN];
460  xn[1][0] = elementBoundaryBarycentersArray[ebN_global*3 +0];
461  xn[1][1] = elementBoundaryBarycentersArray[ebN_global*3 +1];
462  xn[1][2] = elementBoundaryBarycentersArray[ebN_global*3 +2];
463  }
464  else
465  {
466  xn[1][0] = elementBarycentersArray[eN_opposite*3 +0];
467  xn[1][1] = elementBarycentersArray[eN_opposite*3 +1];
468  xn[1][2] = elementBarycentersArray[eN_opposite*3 +2];
469  }
470 
471  computeLocalGradient1d(xn[0],xn[1],dN[0],dN[1],u[0],u[1],dU);
472  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
473  ebN*nElementBoundaries_element*nSpace +
474  0*nSpace +
475  0] = dN[0][0];
476  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
477  ebN*nElementBoundaries_element*nSpace +
478  1*nSpace +
479  0] = dN[1][0];
480  }/*ebN*/
481  }/*eN*/
482  }/*1d*/
483  else if (nSpace == 2)
484  {
485  for (eN = 0; eN < nElements_global; eN++)
486  {
487  xn[0][0] = elementBarycentersArray[eN*3 + 0];
488  xn[0][1] = elementBarycentersArray[eN*3 + 1];
489  xn[0][2] = elementBarycentersArray[eN*3 + 2];
490 
491  /*first collect all the relevant barycenters*/
492  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
493  {
494  eN_opposite = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
495 
496  if (eN_opposite < 0)/*outside domain*/
497  {
498  ebN_global = elementBoundariesArray[eN*nElementBoundaries_element + ebN];
499  xn[ebN+1][0] = elementBoundaryBarycentersArray[ebN_global*3 +0];
500  xn[ebN+1][1] = elementBoundaryBarycentersArray[ebN_global*3 +1];
501  xn[ebN+1][2] = elementBoundaryBarycentersArray[ebN_global*3 +2];
502  }
503  else
504  {
505  xn[ebN+1][0] = elementBarycentersArray[eN_opposite*3 +0];
506  xn[ebN+1][1] = elementBarycentersArray[eN_opposite*3 +1];
507  xn[ebN+1][2] = elementBarycentersArray[eN_opposite*3 +2];
508  }
509  }
510 
511  /*
512  now form local approximations, local neighbor simplex i consists of
513  (\bar{x}_{eN},\bar{x}^i_{eN},\bar{x}^{i+1}_{eN}
514  */
515  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
516  {
517  ebNplus1 = (ebN+1) % nElementBoundaries_element;
518  /*offset of 1 is because xn[0] is this element*/
519  computeLocalGradient2d(xn[0],xn[ebN+1],xn[ebNplus1+1],
520  dN[0],dN[1],dN[2],
521  u[0],u[1],u[2],dU);
522  /*mwf debug
523  printf("eN=%d ebN=%d ebNplus1=%d xn[0]=[%g,%g], xn[ebN+1]=[%g,%g],xn[ebNplus1+1]=[%g,%g]\n",
524  eN,ebN,ebNplus1,xn[0][0],xn[0][1],xn[ebN+1][0],xn[ebN+1][1],
525  xn[ebNplus1+1][0],xn[ebNplus1+1][1]);
526  printf("\t dN[0]=[%g,%g],dN[1]=[%g,%g],dN[2]=[%g,%g]\n",
527  dN[0][0],dN[0][1],dN[1][0],dN[1][1],dN[2][0],dN[2][1]);
528  */
529  for (I = 0; I < nSpace; I++)
530  {
531  /*this element*/
532  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
533  ebN*nElementBoundaries_element*nSpace +
534  0*nSpace +
535  I] = dN[0][I];
536  /*neighbor ebN*/
537  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
538  ebN*nElementBoundaries_element*nSpace +
539  1*nSpace +
540  I] = dN[1][I];
541  /*neighbor ebN+1*/
542  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
543  ebN*nElementBoundaries_element*nSpace +
544  2*nSpace +
545  I] = dN[2][I];
546 
547  }/*I*/
548  }/*ebN*/
549  }/*eN*/
550  }/*2d*/
551  else /*3d*/
552  {
553  for (eN = 0; eN < nElements_global; eN++)
554  {
555  xn[0][0] = elementBarycentersArray[eN*3 + 0];
556  xn[0][1] = elementBarycentersArray[eN*3 + 1];
557  xn[0][2] = elementBarycentersArray[eN*3 + 2];
558 
559  /*first collect all the relevant barycenters*/
560  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
561  {
562  eN_opposite = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
563 
564  if (eN_opposite < 0)/*outside domain*/
565  {
566  ebN_global = elementBoundariesArray[eN*nElementBoundaries_element + ebN];
567  xn[ebN+1][0] = elementBoundaryBarycentersArray[ebN_global*3 +0];
568  xn[ebN+1][1] = elementBoundaryBarycentersArray[ebN_global*3 +1];
569  xn[ebN+1][2] = elementBoundaryBarycentersArray[ebN_global*3 +2];
570  /*mwf debug*/
571 /* printf("Durlofsky3d calling computeLocalGradient3d eN=%d ebN=%d eN_opposite= %d\n", */
572 /* eN,ebN,eN_opposite); */
573 
574 /* for (I=0; I < nSpace; I++) */
575 /* { */
576 /* printf("xn[%d][%d]= %g \n",ebN+1,I,xn[ebN+1][I]); */
577 /* } */
578 
579  }
580  else
581  {
582  xn[ebN+1][0] = elementBarycentersArray[eN_opposite*3 +0];
583  xn[ebN+1][1] = elementBarycentersArray[eN_opposite*3 +1];
584  xn[ebN+1][2] = elementBarycentersArray[eN_opposite*3 +2];
585  /*mwf debug*/
586 /* printf("Durlofsky3d calling computeLocalGradient3d eN=%d ebN=%d eN_opposite= %d\n", */
587 /* eN,ebN,eN_opposite); */
588 
589 /* for (I=0; I < nSpace; I++) */
590 /* { */
591 /* printf("xn[%d][%d]= %g \n",ebN+1,I,xn[ebN+1][I]); */
592 /* } */
593 
594  }
595  }
596 
597  /*
598  now form local approximations, local neighbor simplex i consists of
599  (\bar{x}_{eN},\bar{x}^i_{eN},\bar{x}^{i+1}_{eN},\bar{x}^{i+2}
600  */
601  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
602  {
603  ebNplus1 = (ebN+1) % nElementBoundaries_element;
604  ebNplus2 = (ebN+2) % nElementBoundaries_element;
605  /*offset of 1 is because xn[0] is this element*/
606  /*mwf debug*/
607 /* printf("Durlofsky3d calling computeLocalGradient3d eN=%d ebN=%d nElementBoundaries_element=%d \n", */
608 /* eN,ebN,nElementBoundaries_element); */
609 /* for (I=0; I < nSpace; I++) */
610 /* { */
611 /* printf("xn[%d][%d]= %g xn[%d][%d]= %g xn[%d][%d]= %g xn[%d][%d]= %g \n", */
612 /* 0,I,xn[0][I], */
613 /* ebN+1,I,xn[ebN+1][I], */
614 /* ebNplus1+1,I,xn[ebNplus1+1][I], */
615 /* ebNplus2+1,I,xn[ebNplus2+1][I]); */
616 /* } */
617  computeLocalGradient3d(xn[0],xn[ebN+1],xn[ebNplus1+1],xn[ebNplus2+1],
618  dN[0],dN[1],dN[2],dN[3],
619  u[0],u[1],u[2],u[3],dU);
620  for (I = 0; I < nSpace; I++)
621  {
622  /*this element*/
623  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
624  ebN*nElementBoundaries_element*nSpace +
625  0*nSpace +
626  I] = dN[0][I];
627  /*neighbor ebN*/
628  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
629  ebN*nElementBoundaries_element*nSpace +
630  1*nSpace +
631  I] = dN[1][I];
632  /*neighbor ebN+1*/
633  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
634  ebN*nElementBoundaries_element*nSpace +
635  2*nSpace +
636  I] = dN[2][I];
637  /*neighbor ebN+2*/
638  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
639  ebN*nElementBoundaries_element*nSpace +
640  3*nSpace +
641  I] = dN[3][I];
642 
643  }/*I*/
644  }/*ebN*/
645  }/*eN*/
646  }/*3d*/
647 }
648 
649 void computeCockburnDGlimiterArrays2d(int nElements_global,
650  int nElementBoundaries_element,
651  int nSpace,
652  const int * elementBoundariesArray,
653  const int * elementNeighborsArray,
654  const double * elementBarycentersArray,
655  const double * elementBoundaryBarycentersArray,
656  const double * elementNeighborShapeGradients,
657  double * alphas,
658  int * alphaNeighbors)
659 {
660  int eN,ebN,ebN_global,ebN_I,ebN_Iplus1,eN_ebN_I,eN_ebN_Iplus1,lamn[2],
661  onBoundary,positiveFound;
662  double lam[2];
663 
664  for (eN=0; eN < nElements_global; eN++)
665  {
666  onBoundary =
667  elementNeighborsArray[eN*nElementBoundaries_element + 0] < 0 ||
668  elementNeighborsArray[eN*nElementBoundaries_element + 1] < 0 ||
669  elementNeighborsArray[eN*nElementBoundaries_element + 2] < 0;
670 
671  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
672  {
673  ebN_global = elementBoundariesArray[eN*nElementBoundaries_element+ebN];
674  if (onBoundary > 0)
675  {
676  lam[0] = 0.0; lam[1] = 0.0;
677  lamn[0]=-1; lamn[1]=-1;
678  alphas[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 0] = lam[0];
679  alphas[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 1] = lam[1];
680  alphaNeighbors[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 0] = lamn[0];
681  alphaNeighbors[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 1] = lamn[1];
682  }
683  else
684  {
685  positiveFound = 0; ebN_I=0;
686  while (ebN_I < nElementBoundaries_element && positiveFound == 0)
687  {
688  lam[0] = 0.0; lam[1] = 0.0; lamn[0] = -1; lamn[1] = -1;
689  /***********
690  in local neighbor simplex ordering
691  0 is always this element
692  1 is element across from current face (say i)
693  2 is i+1
694  ***********/
695  ebN_Iplus1 = (ebN_I + 1) % nElementBoundaries_element;
696  eN_ebN_I = elementNeighborsArray[eN*nElementBoundaries_element+ebN_I];
697  eN_ebN_Iplus1 = elementNeighborsArray[eN*nElementBoundaries_element+ebN_Iplus1];
698  lamn[0] = eN_ebN_I; lamn[1] = eN_ebN_Iplus1;
699 
700  lam[0] = 1.0 +
701  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
702  ebN_I*nElementBoundaries_element*nSpace +
703  1*nSpace + 0]
704  *(elementBoundaryBarycentersArray[ebN_global*3 + 0] - elementBarycentersArray[eN_ebN_I*3 + 0])
705  +
706  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
707  ebN_I*nElementBoundaries_element*nSpace +
708  1*nSpace + 1]
709  *(elementBoundaryBarycentersArray[ebN_global*3 + 1] - elementBarycentersArray[eN_ebN_I*3 + 1]);
710  lam[1] = 1.0 +
711  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
712  ebN_I*nElementBoundaries_element*nSpace +
713  2*nSpace + 0]
714  *(elementBoundaryBarycentersArray[ebN_global*3 + 0] - elementBarycentersArray[eN_ebN_Iplus1*3 + 0])
715  +
716  elementNeighborShapeGradients[eN*nElementBoundaries_element*nElementBoundaries_element*nSpace +
717  ebN_I*nElementBoundaries_element*nSpace +
718  2*nSpace + 1]
719  *(elementBoundaryBarycentersArray[ebN_global*3 + 1] - elementBarycentersArray[eN_ebN_Iplus1*3 + 1]);
720  if (lam[0] >= -1.0e-10 && lam[1] >= -1.0e-10)
721  {
722  positiveFound = 1;
723  alphas[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 0] = lam[0];
724  alphas[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 1] = lam[1];
725  alphaNeighbors[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 0] = lamn[0];
726  alphaNeighbors[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 1] = lamn[1];
727  }
728  ebN_I += 1;
729  }/*while*/
730  assert(positiveFound); /*asserts not turned on right now stupid gcc python compile options*/
731  }/*not on boundary*/
732  }/*ebn*/
733  }/*eN*/
734 }
735 
736 void applyCockburnDGlimiterP1Lagrange2d(double nu,
737  double Mh2,
738  int nElements_global,
739  int nElementBoundaries_element,
740  int nSpace,
741  int nDOF_element,
742  int * elementNeighborsArray,
743  int * l2g,
744  int * tag,
745  double * alphas,
746  int * alphaNeighbors,
747  double * Uin,
748  double * Uout)
749 {
750  /***********************************************************************
751  try to implement DG limiter for triangles according to Cockburn notes
752  assumes P1 lagrange representation coming in
753  ***********************************************************************/
754  double dU,deltaU[3],tdeltaU[3],pos[3],neg[3],pos_eN,neg_eN,thp,thm,sumDU,
755  uM[3],uMlim[3],uBar,uBar_0,uBar_1;
756  int eN,ebN,ebNplus1,ebNplus2,eN_0,eN_1,onBoundary,tag_eN[3];
757  /*for debugging*/
758  double uavgIn,uavgOut;
759  double oneThird=1.0/3.0;
760  /*coulud check dimensions all agree but asserts turned off*/
761  for (eN = 0; eN < nElements_global; eN++)
762  {
763  uBar = (Uin[l2g[eN*nDOF_element+0]]+Uin[l2g[eN*nDOF_element+1]]+Uin[l2g[eN*nDOF_element+2]])*oneThird;
764  uavgIn = uBar;
765  onBoundary =
766  elementNeighborsArray[eN*nElementBoundaries_element + 0] < 0 ||
767  elementNeighborsArray[eN*nElementBoundaries_element + 1] < 0 ||
768  elementNeighborsArray[eN*nElementBoundaries_element + 2] < 0;
769  if (onBoundary > 0)
770  {
771  tdeltaU[0] = 0.0; tdeltaU[1] = 0.0; tdeltaU[2] = 0.0;
772  }
773  else
774  {
775  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
776  {
777  ebNplus1 = (ebN+1) % nElementBoundaries_element;
778  ebNplus2 = (ebN+2) % nElementBoundaries_element;
779  /*
780  midpoint value for this face assuming correspondence
781  between local node numbers and faces (across)
782  and nodal dofs (same point)
783  */
784  uM[ebN] = 0.5*(Uin[l2g[eN*nDOF_element+ebNplus1]]+Uin[l2g[eN*nDOF_element+ebNplus2]]);
785  eN_0 = alphaNeighbors[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 0];
786  eN_1 = alphaNeighbors[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 1];
787  uBar_0 = (Uin[l2g[eN_0*nDOF_element+0]]+Uin[l2g[eN_0*nDOF_element+1]]+Uin[l2g[eN_0*nDOF_element+2]])*oneThird;
788  uBar_1 = (Uin[l2g[eN_1*nDOF_element+0]]+Uin[l2g[eN_1*nDOF_element+1]]+Uin[l2g[eN_1*nDOF_element+2]])*oneThird;
789  dU = alphas[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 0]*(uBar_0-uBar)+
790  alphas[eN*nElementBoundaries_element*nSpace + ebN*nSpace + 1]*(uBar_1-uBar);
791  deltaU[ebN] = mminmod2(uM[ebN]-uBar,nu*dU,Mh2,&tag_eN[ebN]);
792  tdeltaU[ebN]= deltaU[ebN];
793  }
794  /*check that slopes sum to zero for mass conservation*/
795  sumDU = tdeltaU[0] + tdeltaU[1] + tdeltaU[2];
796  if (fabs(sumDU) > 1.0e-6) /*need tolerance*/
797  {
798  /*should unroll this stuff too*/
799  pos_eN = 0.0; neg_eN = 0.0;
800  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
801  {
802  pos[ebN] = ( deltaU[ebN] > 0.0) ? deltaU[ebN] : 0.0;
803  neg[ebN] = (-deltaU[ebN] > 0.0) ? -deltaU[ebN] : 0.0;
804  pos_eN += pos[ebN];
805  neg_eN += neg[ebN];
806  }
807  thp = (1.0 < neg_eN/pos_eN) ? 1.0 : neg_eN/pos_eN;
808  thm = (1.0 < pos_eN/neg_eN) ? 1.0 : pos_eN/neg_eN;
809  for (ebN=0; ebN < nElementBoundaries_element; ebN++)
810  {
811  tdeltaU[ebN] = thp*pos[ebN] - thm*neg[ebN];
812  tag_eN[ebN] = 0;
813  }
814  }/*extra limit for mass bal*/
815  }/*not at boundary*/
816  /*get limited values at midpoints of element boundaries*/
817  /*based on Crouzeix-Raviart type representation with limited slopes
818  so phi_{i}(\bar{x^{j}) = \delta_ij */
819  uMlim[0] = uBar + tdeltaU[0];
820  uMlim[1] = uBar + tdeltaU[1];
821  uMlim[2] = uBar + tdeltaU[2];
822 
823  /*convert back to Lagrange representation for output*/
824  Uout[l2g[eN*nDOF_element + 0]] = uMlim[1]+uMlim[2]-uMlim[0];
825  Uout[l2g[eN*nDOF_element + 1]] = uMlim[2]+uMlim[0]-uMlim[1];
826  Uout[l2g[eN*nDOF_element + 2]] = uMlim[0]+uMlim[1]-uMlim[2];
827  tag[eN] = tag_eN[0] + tag_eN[1] + tag_eN[2] == 3;
828  /*mwf debug*/
829  uavgOut = (Uout[l2g[eN*nDOF_element + 0]]+Uout[l2g[eN*nDOF_element + 1]]+Uout[l2g[eN*nDOF_element + 2]])*oneThird;
830  /*
831  if (fabs(uavgOut-uavgIn) > 1.0e-7)
832  printf("problem CockburnLimit eN=%d uavgIn=%g uavgOut=%g Uin=[%g,%g,%g] Uout=[%g,%g,%g]\n",
833  eN,uavgIn,uavgOut,
834  Uin[l2g[eN*nDOF_element + 0]],Uin[l2g[eN*nDOF_element + 1]],Uin[l2g[eN*nDOF_element + 2]],
835  Uout[l2g[eN*nDOF_element + 0]],Uout[l2g[eN*nDOF_element + 1]],Uout[l2g[eN*nDOF_element + 2]]);
836  */
837  }/*eN*/
838 }
839 void computeElementAveragesP1Lagrange(int nElements_global,
840  int nDOF_element,
841  const int * l2g,
842  const double * Uin,
843  double * elementAverages)
844 {
845  int eN,i,I;
846  const double denom = 1.0/nDOF_element;
847  for (eN = 0; eN < nElements_global; eN++)
848  {
849  elementAverages[eN] = 0.0;
850  for (i=0; i < nDOF_element; i++)
851  {
852  I = l2g[eN*nDOF_element+i];
853  elementAverages[eN] += Uin[I];
854  }
855  elementAverages[eN] *= denom;
856  }
857 }
858 /*
859  try simple sorting of A, but also keep track of sorted order
860  for other quantity based on A
861 */
862 void shortSortDescending(double *A, int * order, int len)
863 {
864  register int i,j,itmp;
865  register double tmp;
866  for (i=0; i < len; i++)
867  {
868  /*order[i] = i;*/
869  for (j = i+1; j < len; j++)
870  {
871  if (A[j] > A[i])
872  {
873  tmp = A[i]; A[i]= A[j]; A[j]= tmp;
874  itmp= order[i]; order[i] = order[j]; order[j] = itmp;
875  }
876  }
877  }
878 }
879 void applyDurlofskyDGlimiterP1Lagrange2d(int killExtrema,
880  int allowMinWithUndershoot,
881  int nElements_global,
882  int nElementBoundaries_element,
883  int nNodes_element,
884  int nSpace,
885  int nDOF_element,
886  const int * elementNeighborsArray,
887  const int * elementBoundariesArray,
888  const int * elementNodesArray,
889  const double * nodeArray,
890  const double * elementBarycentersArray,
891  const double * elementBoundaryBarycentersArray,
892  const double * elementNeighborShapeGradients,
893  const int * l2g,
894  const double * grad_v0,
895  double * elementAverages,
896  int * tag,
897  double * Uin,
898  double * Uout)
899 {
900  /*this should be Durlofsky's original version with Phi limiter (not the min one)*/
901  int eN,nN_global,ebN,ebN_global,ebN_1,eN_ebN,eN_ebN_1,onBoundary,isExtremum;
902  int dUdescendingOrder[4];
903  double normDU[4],dU[4][2],max_ebN,min_ebN,uM;
904  /*mwf debug*/
905  double normDUsave[4];
906  int maxFound,minFound,islot,i,itmp,okGradient;
907  int nElementBoundaries_element2 = nElementBoundaries_element*nElementBoundaries_element;
908  register int DOF0,DOF1,DOF2;
909  register double dUlim0,dUlim1;
910  register double uBar,uBar_ebN,uBar_ebN_1;
911 
912  const double otol = 1.0e-5;
913 
914  computeElementAveragesP1Lagrange(nElements_global,nDOF_element,l2g,Uin,elementAverages);
915 
916  for (eN = 0; eN < nElements_global; eN++)
917  {
918  /*current element*/
919  uBar = elementAverages[eN];
920  DOF0 = l2g[eN*nDOF_element+0]; DOF1 = l2g[eN*nDOF_element+1]; DOF2 = l2g[eN*nDOF_element+2];
921  onBoundary =
922  elementNeighborsArray[eN*nElementBoundaries_element + 0] < 0 ||
923  elementNeighborsArray[eN*nElementBoundaries_element + 1] < 0 ||
924  elementNeighborsArray[eN*nElementBoundaries_element + 2] < 0;
925  if (onBoundary > 0)
926  {
927  dUlim0 = 0.0; dUlim1 = 0.0;
928  tag[eN] = 0;
929  }
930  else
931  {
932  /*check to see if this is a local extremum*/
933  maxFound = 0; minFound = 0;
934  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
935  {
936  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
937  uBar_ebN = elementAverages[eN_ebN];
938  maxFound = uBar_ebN > uBar ? 1 : maxFound;
939  minFound = uBar_ebN < uBar ? 1 : minFound;
940  }
941  isExtremum = (maxFound*minFound == 0);
942  if (isExtremum && killExtrema > 0)
943  {
944  /*mwf debug
945  printf("Durlofsky extremum found eN=%d uBar= %g uBar_ebN[%g,%g,%g]\n",
946  eN,uBar,uBar_ebN[0],uBar_ebN[1],uBar_ebN[2]);
947  */
948  dUlim0 = 0.0; dUlim1 = 0.0;
949  tag[eN] = 0;
950  }
951  else
952  {
953  /*by default take zero slope*/
954  dUlim0 = 0.0; dUlim1 = 0.0;
955  dU[0][0] =
956  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 0*nSpace + 0]*Uin[DOF0]+
957  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 1*nSpace + 0]*Uin[DOF1]+
958  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 2*nSpace + 0]*Uin[DOF2];
959  dU[0][1] =
960  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 0*nSpace + 1]*Uin[DOF0]+
961  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 1*nSpace + 1]*Uin[DOF1]+
962  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 2*nSpace + 1]*Uin[DOF2];
963  normDU[0] = sqrt(dU[0][0]*dU[0][0] + dU[0][1]*dU[0][1]);
964  normDU[1] = 0.0; normDU[2] = 0.0; normDU[3] = 0.0;
965  /*loop through neighbor simplexes and compute local interpolant gradients*/
966  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
967  {
968  ebN_1 = (ebN+1) % nElementBoundaries_element;
969  /*two neighbors for this simplex
970  local numbering is
971  0 <--> this element
972  1 <--> ebN neighbor
973  2 <--> ebN+1 neighbor
974  */
975  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
976  eN_ebN_1 = elementNeighborsArray[eN*nElementBoundaries_element + ebN_1];
977  uBar_ebN = elementAverages[eN_ebN];
978  uBar_ebN_1= elementAverages[eN_ebN_1];
979 
980  /*local number zero is always this element*/
981  dU[ebN+1][0] = 0.0; dU[ebN+1][1] = 0.0;
982  dU[ebN+1][0]=
983  uBar*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
984  ebN*nElementBoundaries_element*nSpace +
985  0*nSpace + 0]
986  +
987  uBar_ebN*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
988  ebN*nElementBoundaries_element*nSpace +
989  1*nSpace + 0]
990  +
991  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
992  ebN*nElementBoundaries_element*nSpace +
993  2*nSpace + 0];
994 
995  dU[ebN+1][1]=
996  uBar*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
997  ebN*nElementBoundaries_element*nSpace +
998  0*nSpace + 1]
999  +
1000  uBar_ebN*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1001  ebN*nElementBoundaries_element*nSpace +
1002  1*nSpace + 1]
1003  +
1004  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1005  ebN*nElementBoundaries_element*nSpace +
1006  2*nSpace + 1];
1007  normDU[ebN+1] = sqrt(dU[ebN+1][0]*dU[ebN+1][0] + dU[ebN+1][1]*dU[ebN+1][1]);
1008  }/*ebN*/
1009  /*now sort the gradients from largest to smallest*/
1010  dUdescendingOrder[0] = 0; dUdescendingOrder[1] = 1; dUdescendingOrder[2]=2; dUdescendingOrder[3]=3;
1011  normDUsave[0] = normDU[0]; normDUsave[1]=normDU[1]; normDUsave[2]=normDU[2]; normDUsave[3]=normDU[3];
1012  shortSortDescending(normDU,dUdescendingOrder,4);
1013  /*mwf debug check ordering
1014  for (i = 0; i < nElementBoundaries_element; i++)
1015  {
1016  if (normDU[i+1] > normDU[i])
1017  {
1018  printf("\nPROBLEM Durlofsky out of order eN=%d normDU[%d]=%g > normDU[%d] =%g \n",
1019  eN,i+1,normDU[i+1],
1020  i,normDU[i]);
1021  printf("normDU=[%g,%g,%g,%g] dUdescendingOrder=[%d,%d,%d,%d]\n",
1022  normDU[0],normDU[1],normDU[2],normDU[3],
1023  dUdescendingOrder[0],dUdescendingOrder[1],dUdescendingOrder[2],dUdescendingOrder[3]);
1024  exit(1);
1025  }
1026  if (fabs(normDU[i]-normDUsave[dUdescendingOrder[i]]) > 1.0e-8)
1027  {
1028  printf("\nPROBLEM Durlofsky order wrong eN=%d normDU[%d]=%g normDUsave[%d] =%g \n",
1029  eN,i,normDU[i],dUdescendingOrder[i],normDUsave[dUdescendingOrder[i]]);
1030  printf("normDU=[%g,%g,%g,%g] normDUsave=[%g,%g,%g,%g] dUdescendingOrder=[%d,%d,%d,%d]\n",
1031  normDU[0],normDU[1],normDU[2],normDU[3],
1032  normDUsave[0],normDUsave[1],normDUsave[2],normDUsave[3],
1033  dUdescendingOrder[0],dUdescendingOrder[1],dUdescendingOrder[2],dUdescendingOrder[3]);
1034  exit(1);
1035 
1036  }
1037 
1038  }
1039  mwf debug end */
1040  /*now start checking for overshoot, starting with largest du*/
1041  okGradient = 0; islot = 0;
1042  /*use minimum slope if undershoot/overshoot detected for others*/
1043  while (okGradient == 0 && islot < nElementBoundaries_element)
1044  {
1045  itmp = dUdescendingOrder[islot];
1046  /*start with ok*/
1047  okGradient = 1;
1048  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1049  {
1050  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
1051  uBar_ebN = elementAverages[eN_ebN];
1052  ebN_global = elementBoundariesArray[eN*nElementBoundaries_element + ebN];
1053  uM = uBar
1054  +
1055  dU[itmp][0]*(elementBoundaryBarycentersArray[ebN_global*3 + 0]-
1056  elementBarycentersArray[eN*3 + 0])
1057  +
1058  dU[itmp][1]*(elementBoundaryBarycentersArray[ebN_global*3 + 1]-
1059  elementBarycentersArray[eN*3 + 1]);
1060  max_ebN = uBar > uBar_ebN ? uBar : uBar_ebN;
1061  min_ebN = uBar < uBar_ebN ? uBar : uBar_ebN;
1062 
1063  okGradient = (okGradient > 0) && min_ebN <= uM -otol && uM <= max_ebN + otol;
1064 
1065  }
1066  islot++;/*undo this later if okGradient==1*/
1067  }
1068  if (okGradient == 1)
1069  islot--;
1070  if (okGradient || allowMinWithUndershoot > 0)
1071  {
1072  itmp = dUdescendingOrder[islot];
1073  dUlim0 = dU[itmp][0];
1074  dUlim1 = dU[itmp][1];
1075  tag[eN] = islot == 0 ? 1 : 0;
1076  /*mwf debug
1077  printf("Durlofsky eN=%d okGradient islot=%d dUlim[0]=%g dUlim[1]=%g \n",
1078  eN,islot,dUlim[0],dUlim[1]);
1079  */
1080  }
1081 
1082  }/*not an extremum*/
1083  }/*in interior*/
1084  /*interpolate values to nodes assuming correspondence with local dofs*/
1085  nN_global = elementNodesArray[eN*nNodes_element+0];
1086  Uout[DOF0] = uBar
1087  +
1088  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1089  +
1090  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1]);
1091  nN_global = elementNodesArray[eN*nNodes_element+1];
1092  Uout[DOF1] = uBar
1093  +
1094  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1095  +
1096  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1]);
1097  nN_global = elementNodesArray[eN*nNodes_element+2];
1098  Uout[DOF2] = uBar
1099  +
1100  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1101  +
1102  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1]);
1103 
1104  }/*eN*/
1105 }
1106 
1107 void applyDurlofskyDGlimiterP1Lagrange3d(int killExtrema,
1108  int allowMinWithUndershoot,
1109  int nElements_global,
1110  int nElementBoundaries_element,
1111  int nNodes_element,
1112  int nSpace,
1113  int nDOF_element,
1114  const int * elementNeighborsArray,
1115  const int * elementBoundariesArray,
1116  const int * elementNodesArray,
1117  const double * nodeArray,
1118  const double * elementBarycentersArray,
1119  const double * elementBoundaryBarycentersArray,
1120  const double * elementNeighborShapeGradients,
1121  const int * l2g,
1122  const double * grad_v0,
1123  double * elementAverages,
1124  int * tag,
1125  double * Uin,
1126  double * Uout)
1127 {
1128  /*try 3d version of Durlofsky's original version with Phi limiter (not the min one)*/
1129  int eN,nN_global,ebN,ebN_global,ebN_1,ebN_2,eN_ebN,eN_ebN_1,eN_ebN_2,onBoundary,isExtremum;
1130  int dUdescendingOrder[5];
1131  double normDU[5],dU[5][3],max_ebN,min_ebN,uM;
1132  /*mwf for debugging*/
1133  double normDUsave[5];
1134  int maxFound,minFound,islot,i,itmp,okGradient;
1135  int nElementBoundaries_element2 = nElementBoundaries_element*nElementBoundaries_element;
1136  register int DOF0,DOF1,DOF2,DOF3;
1137  register double dUlim0,dUlim1,dUlim2;
1138  register double uBar,uBar_ebN,uBar_ebN_1,uBar_ebN_2;
1139 
1140  const double otol = 1.0e-5;
1141  computeElementAveragesP1Lagrange(nElements_global,nDOF_element,l2g,Uin,elementAverages);
1142 
1143  for (eN = 0; eN < nElements_global; eN++)
1144  {
1145  /*current element*/
1146  uBar = elementAverages[eN];
1147  DOF0 = l2g[eN*nDOF_element+0]; DOF1 = l2g[eN*nDOF_element+1];
1148  DOF2 = l2g[eN*nDOF_element+2]; DOF3 = l2g[eN*nDOF_element+3];
1149  onBoundary =
1150  elementNeighborsArray[eN*nElementBoundaries_element + 0] < 0 ||
1151  elementNeighborsArray[eN*nElementBoundaries_element + 1] < 0 ||
1152  elementNeighborsArray[eN*nElementBoundaries_element + 2] < 0 ||
1153  elementNeighborsArray[eN*nElementBoundaries_element + 3] < 0;
1154  if (onBoundary > 0)
1155  {
1156  dUlim0 = 0.0; dUlim1 = 0.0; dUlim2 = 0.0;
1157  tag[eN] = 0;
1158  }
1159  else
1160  {
1161  /*check to see if this is a local extremum*/
1162  maxFound = 0; minFound = 0;
1163  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1164  {
1165  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
1166  uBar_ebN = elementAverages[eN_ebN];
1167  maxFound = uBar_ebN > uBar ? 1 : maxFound;
1168  minFound = uBar_ebN < uBar ? 1 : minFound;
1169  }
1170  isExtremum = (maxFound*minFound == 0);
1171  if (isExtremum && killExtrema > 0)
1172  {
1173  /*mwf debug
1174  printf("Durlofsky extremum found eN=%d uBar= %g uBar_ebN[%g,%g,%g]\n",
1175  eN,uBar,uBar_ebN[0],uBar_ebN[1],uBar_ebN[2]);
1176  */
1177  dUlim0 = 0.0; dUlim1 = 0.0; dUlim2 = 0.0;
1178  tag[eN] = 0;
1179  }
1180  else
1181  {
1182  /*by default take zero slope*/
1183  dUlim0 = 0.0; dUlim1 = 0.0; dUlim2 = 0.0;
1184  dU[0][0] =
1185  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 0*nSpace + 0]*Uin[DOF0]+
1186  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 1*nSpace + 0]*Uin[DOF1]+
1187  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 2*nSpace + 0]*Uin[DOF2]+
1188  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 3*nSpace + 0]*Uin[DOF3];
1189  dU[0][1] =
1190  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 0*nSpace + 1]*Uin[DOF0]+
1191  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 1*nSpace + 1]*Uin[DOF1]+
1192  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 2*nSpace + 1]*Uin[DOF2]+
1193  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 3*nSpace + 1]*Uin[DOF3];
1194  dU[0][2] =
1195  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 0*nSpace + 2]*Uin[DOF0]+
1196  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 1*nSpace + 2]*Uin[DOF1]+
1197  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 2*nSpace + 2]*Uin[DOF2]+
1198  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 3*nSpace + 2]*Uin[DOF3];
1199  normDU[0] = sqrt(dU[0][0]*dU[0][0] + dU[0][1]*dU[0][1] + dU[0][2]*dU[0][2]);
1200  normDU[1] = 0.0; normDU[2] = 0.0; normDU[3] = 0.0; normDU[4] = 0.0;
1201  /*loop through neighbor simplexes and compute local interpolant gradients*/
1202  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1203  {
1204  ebN_1 = (ebN+1) % nElementBoundaries_element;
1205  ebN_2 = (ebN+2) % nElementBoundaries_element;
1206  /*two neighbors for this simplex
1207  local numbering is
1208  0 <--> this element
1209  1 <--> ebN neighbor
1210  2 <--> ebN+1 neighbor
1211  3 <--> ebN+2 neighbor
1212  */
1213  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
1214  eN_ebN_1 = elementNeighborsArray[eN*nElementBoundaries_element + ebN_1];
1215  eN_ebN_2 = elementNeighborsArray[eN*nElementBoundaries_element + ebN_2];
1216  uBar_ebN = elementAverages[eN_ebN];
1217  uBar_ebN_1= elementAverages[eN_ebN_1];
1218  uBar_ebN_2= elementAverages[eN_ebN_2];
1219 
1220  /*local number zero is always this element*/
1221  dU[ebN+1][0] = 0.0; dU[ebN+1][1] = 0.0; dU[ebN+1][2] = 0.0;
1222  dU[ebN+1][0]=
1223  uBar*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1224  ebN*nElementBoundaries_element*nSpace +
1225  0*nSpace + 0]
1226  +
1227  uBar_ebN*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1228  ebN*nElementBoundaries_element*nSpace +
1229  1*nSpace + 0]
1230  +
1231  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1232  ebN*nElementBoundaries_element*nSpace +
1233  2*nSpace + 0]
1234  +
1235  uBar_ebN_2*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1236  ebN*nElementBoundaries_element*nSpace +
1237  3*nSpace + 0];
1238 
1239  dU[ebN+1][1]=
1240  uBar*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1241  ebN*nElementBoundaries_element*nSpace +
1242  0*nSpace + 1]
1243  +
1244  uBar_ebN*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1245  ebN*nElementBoundaries_element*nSpace +
1246  1*nSpace + 1]
1247  +
1248  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1249  ebN*nElementBoundaries_element*nSpace +
1250  2*nSpace + 1]
1251  +
1252  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1253  ebN*nElementBoundaries_element*nSpace +
1254  3*nSpace + 1];
1255 
1256  dU[ebN+1][2]=
1257  uBar*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1258  ebN*nElementBoundaries_element*nSpace +
1259  0*nSpace + 2]
1260  +
1261  uBar_ebN*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1262  ebN*nElementBoundaries_element*nSpace +
1263  1*nSpace + 2]
1264  +
1265  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1266  ebN*nElementBoundaries_element*nSpace +
1267  2*nSpace + 2]
1268  +
1269  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1270  ebN*nElementBoundaries_element*nSpace +
1271  3*nSpace + 2];
1272 
1273  normDU[ebN+1] = sqrt(dU[ebN+1][0]*dU[ebN+1][0] + dU[ebN+1][1]*dU[ebN+1][1] + dU[ebN+1][2]*dU[ebN+1][2]);
1274  }/*ebN*/
1275  /*now sort the gradients from largest to smallest*/
1276  dUdescendingOrder[0] = 0; dUdescendingOrder[1] = 1; dUdescendingOrder[2]=2;
1277  dUdescendingOrder[3]=3; dUdescendingOrder[4] = 4;
1278  /*save for debugging*/
1279  normDUsave[0] = normDU[0]; normDUsave[1]=normDU[1]; normDUsave[2]=normDU[2]; normDUsave[3]=normDU[3];
1280  normDUsave[4] = normDU[4];
1281  shortSortDescending(normDU,dUdescendingOrder,5);
1282  /*mwf debug check ordering
1283  for (i = 0; i < nElementBoundaries_element; i++)
1284  {
1285  if (normDU[i+1] > normDU[i])
1286  {
1287  printf("\nPROBLEM Durlofsky 3d out of order eN=%d normDU[%d]=%g > normDU[%d] =%g \n",
1288  eN,i+1,normDU[i+1],
1289  i,normDU[i]);
1290  printf("normDU=[%g,%g,%g,%g,%g] dUdescendingOrder=[%d,%d,%d,%d,%d]\n",
1291  normDU[0],normDU[1],normDU[2],normDU[3],normDU[4],
1292  dUdescendingOrder[0],dUdescendingOrder[1],dUdescendingOrder[2],
1293  dUdescendingOrder[3],dUdescendingOrder[4]);
1294  exit(1);
1295  }
1296  if (fabs(normDU[i]-normDUsave[dUdescendingOrder[i]]) > 1.0e-8)
1297  {
1298  printf("\nPROBLEM Durlofsky order wrong eN=%d normDU[%d]=%g normDUsave[%d] =%g \n",
1299  eN,i,normDU[i],dUdescendingOrder[i],normDUsave[dUdescendingOrder[i]]);
1300  printf("normDU=[%g,%g,%g,%g,%g] normDUsave=[%g,%g,%g,%g,%g] dUdescendingOrder=[%d,%d,%d,%d,%d]\n",
1301  normDU[0],normDU[1],normDU[2],normDU[3],normDU[4],
1302  normDUsave[0],normDUsave[1],normDUsave[2],normDUsave[3],normDUsave[4],
1303  dUdescendingOrder[0],dUdescendingOrder[1],dUdescendingOrder[2],dUdescendingOrder[3],
1304  dUdescendingOrder[4]);
1305  exit(1);
1306 
1307  }
1308  }
1309  mwf debug end */
1310  /*now start checking for overshoot, starting with largest du*/
1311  okGradient = 0; islot = 0;
1312  /*use minimum slope if undershoot/overshoot detected for others*/
1313  while (okGradient == 0 && islot < nElementBoundaries_element)
1314  {
1315  itmp = dUdescendingOrder[islot];
1316  /*start with ok*/
1317  okGradient = 1;
1318  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1319  {
1320  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
1321  uBar_ebN = elementAverages[eN_ebN];
1322  ebN_global = elementBoundariesArray[eN*nElementBoundaries_element + ebN];
1323  uM = uBar
1324  +
1325  dU[itmp][0]*(elementBoundaryBarycentersArray[ebN_global*3 + 0]-
1326  elementBarycentersArray[eN*3 + 0])
1327  +
1328  dU[itmp][1]*(elementBoundaryBarycentersArray[ebN_global*3 + 1]-
1329  elementBarycentersArray[eN*3 + 1])
1330  +
1331  dU[itmp][2]*(elementBoundaryBarycentersArray[ebN_global*3 + 2]-
1332  elementBarycentersArray[eN*3 + 2]);
1333 
1334  max_ebN = uBar > uBar_ebN ? uBar : uBar_ebN;
1335  min_ebN = uBar < uBar_ebN ? uBar : uBar_ebN;
1336 
1337  okGradient = (okGradient > 0) && min_ebN <= uM -otol && uM <= max_ebN + otol;
1338 
1339  }
1340  islot++;/*undo this later if okGradient==1*/
1341  }
1342  if (okGradient == 1)
1343  islot--;
1344  if (okGradient || allowMinWithUndershoot > 0)
1345  {
1346  itmp = dUdescendingOrder[islot];
1347  dUlim0 = dU[itmp][0];
1348  dUlim1 = dU[itmp][1];
1349  dUlim2 = dU[itmp][2];
1350  tag[eN] = islot == 0 ? 1 : 0;
1351  /*mwf debug
1352  printf("Durlofsky eN=%d okGradient islot=%d dUlim[0]=%g dUlim[1]=%g dUlim[2]=%g\n",
1353  eN,islot,dUlim0,dUlim1,dUlim2);
1354  */
1355  }
1356 
1357 
1358  }/*not an extremum*/
1359  }/*in interior*/
1360  /*interpolate values to nodes assuming correspondence with local dofs*/
1361  nN_global = elementNodesArray[eN*nNodes_element+0];
1362  Uout[DOF0] = uBar
1363  +
1364  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1365  +
1366  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1])
1367  +
1368  dUlim2*(nodeArray[nN_global*3 + 2] - elementBarycentersArray[eN*3 + 2]);
1369  nN_global = elementNodesArray[eN*nNodes_element+1];
1370  Uout[DOF1] = uBar
1371  +
1372  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1373  +
1374  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1])
1375  +
1376  dUlim2*(nodeArray[nN_global*3 + 2] - elementBarycentersArray[eN*3 + 2]);
1377  nN_global = elementNodesArray[eN*nNodes_element+2];
1378  Uout[DOF2] = uBar
1379  +
1380  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1381  +
1382  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1])
1383  +
1384  dUlim2*(nodeArray[nN_global*3 + 2] - elementBarycentersArray[eN*3 + 2]);
1385  nN_global = elementNodesArray[eN*nNodes_element+3];
1386  Uout[DOF3] = uBar
1387  +
1388  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1389  +
1390  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1])
1391  +
1392  dUlim2*(nodeArray[nN_global*3 + 2] - elementBarycentersArray[eN*3 + 2]);
1393 
1394  }/*eN*/
1395 }
1396 
1397 void applyDurlofskyDGlimiterP1Lagrange2d_withVacuumTol(int killExtrema,
1398  int allowMinWithUndershoot,
1399  int enforcePositivity,
1400  double vacuumTol,
1401  int nElements_global,
1402  int nElementBoundaries_element,
1403  int nNodes_element,
1404  int nSpace,
1405  int nDOF_element,
1406  const int * elementNeighborsArray,
1407  const int * elementBoundariesArray,
1408  const int * elementNodesArray,
1409  const double * nodeArray,
1410  const double * elementBarycentersArray,
1411  const double * elementBoundaryBarycentersArray,
1412  const double * elementNeighborShapeGradients,
1413  const int * l2g,
1414  const double * grad_v0,
1415  double * elementAverages,
1416  int * tag,
1417  double * Uin,
1418  double * Uout)
1419 {
1420  /*this should be Durlofsky's original version with Phi limiter (not the min one)*/
1421  int eN,nN_global,ebN,ebN_global,ebN_1,eN_ebN,eN_ebN_1,onBoundary,isExtremum;
1422  int dUdescendingOrder[4];
1423  double normDU[4],dU[4][2],max_ebN,min_ebN,uM;
1424  /*mwf debug*/
1425  double normDUsave[4];
1426  int maxFound,minFound,islot,i,itmp,okGradient;
1427  int nElementBoundaries_element2 = nElementBoundaries_element*nElementBoundaries_element;
1428  register int DOF0,DOF1,DOF2;
1429  register double dUlim0,dUlim1;
1430  register double uBar,uBar_ebN,uBar_ebN_1;
1431 
1432  const double otol = 1.0e-5;
1433 
1434  computeElementAveragesP1Lagrange(nElements_global,nDOF_element,l2g,Uin,elementAverages);
1435 
1436  for (eN = 0; eN < nElements_global; eN++)
1437  {
1438  /*current element*/
1439  uBar = elementAverages[eN];
1440  DOF0 = l2g[eN*nDOF_element+0]; DOF1 = l2g[eN*nDOF_element+1]; DOF2 = l2g[eN*nDOF_element+2];
1441  onBoundary =
1442  elementNeighborsArray[eN*nElementBoundaries_element + 0] < 0 ||
1443  elementNeighborsArray[eN*nElementBoundaries_element + 1] < 0 ||
1444  elementNeighborsArray[eN*nElementBoundaries_element + 2] < 0;
1445  if (onBoundary > 0)
1446  {
1447  dUlim0 = 0.0; dUlim1 = 0.0;
1448  tag[eN] = 0;
1449  }
1450  else
1451  {
1452  /*check to see if this is a local extremum*/
1453  maxFound = 0; minFound = 0;
1454  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1455  {
1456  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
1457  uBar_ebN = elementAverages[eN_ebN];
1458  maxFound = uBar_ebN > uBar ? 1 : maxFound;
1459  minFound = uBar_ebN < uBar ? 1 : minFound;
1460  }
1461  isExtremum = (maxFound*minFound == 0);
1462  if (isExtremum && killExtrema > 0)
1463  {
1464  /*mwf debug
1465  printf("Durlofsky extremum found eN=%d uBar= %g uBar_ebN[%g,%g,%g]\n",
1466  eN,uBar,uBar_ebN[0],uBar_ebN[1],uBar_ebN[2]);
1467  */
1468  dUlim0 = 0.0; dUlim1 = 0.0;
1469  tag[eN] = 0;
1470  }
1471  else if (fabs(uBar) < vacuumTol)
1472  {
1473  dUlim0 = 0.0; dUlim1 = 0.0;
1474  tag[eN] = 0;
1475  }
1476  else if (enforcePositivity && uBar < vacuumTol)
1477  {
1478  dUlim0 = 0.0; dUlim1 = 0.0;
1479  tag[eN] = 0;
1480  }
1481  else
1482  {
1483  /*by default take zero slope*/
1484  dUlim0 = 0.0; dUlim1 = 0.0;
1485  dU[0][0] =
1486  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 0*nSpace + 0]*Uin[DOF0]+
1487  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 1*nSpace + 0]*Uin[DOF1]+
1488  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 2*nSpace + 0]*Uin[DOF2];
1489  dU[0][1] =
1490  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 0*nSpace + 1]*Uin[DOF0]+
1491  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 1*nSpace + 1]*Uin[DOF1]+
1492  grad_v0[eN*1*nDOF_element*nSpace + 0*nDOF_element*nSpace + 2*nSpace + 1]*Uin[DOF2];
1493  normDU[0] = sqrt(dU[0][0]*dU[0][0] + dU[0][1]*dU[0][1]);
1494  normDU[1] = 0.0; normDU[2] = 0.0; normDU[3] = 0.0;
1495  /*loop through neighbor simplexes and compute local interpolant gradients*/
1496  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1497  {
1498  ebN_1 = (ebN+1) % nElementBoundaries_element;
1499  /*two neighbors for this simplex
1500  local numbering is
1501  0 <--> this element
1502  1 <--> ebN neighbor
1503  2 <--> ebN+1 neighbor
1504  */
1505  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
1506  eN_ebN_1 = elementNeighborsArray[eN*nElementBoundaries_element + ebN_1];
1507  uBar_ebN = elementAverages[eN_ebN];
1508  uBar_ebN_1= elementAverages[eN_ebN_1];
1509 
1510  /*local number zero is always this element*/
1511  dU[ebN+1][0] = 0.0; dU[ebN+1][1] = 0.0;
1512  dU[ebN+1][0]=
1513  uBar*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1514  ebN*nElementBoundaries_element*nSpace +
1515  0*nSpace + 0]
1516  +
1517  uBar_ebN*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1518  ebN*nElementBoundaries_element*nSpace +
1519  1*nSpace + 0]
1520  +
1521  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1522  ebN*nElementBoundaries_element*nSpace +
1523  2*nSpace + 0];
1524 
1525  dU[ebN+1][1]=
1526  uBar*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1527  ebN*nElementBoundaries_element*nSpace +
1528  0*nSpace + 1]
1529  +
1530  uBar_ebN*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1531  ebN*nElementBoundaries_element*nSpace +
1532  1*nSpace + 1]
1533  +
1534  uBar_ebN_1*elementNeighborShapeGradients[eN*nElementBoundaries_element2*nSpace +
1535  ebN*nElementBoundaries_element*nSpace +
1536  2*nSpace + 1];
1537  normDU[ebN+1] = sqrt(dU[ebN+1][0]*dU[ebN+1][0] + dU[ebN+1][1]*dU[ebN+1][1]);
1538  }/*ebN*/
1539  /*now sort the gradients from largest to smallest*/
1540  dUdescendingOrder[0] = 0; dUdescendingOrder[1] = 1; dUdescendingOrder[2]=2; dUdescendingOrder[3]=3;
1541  normDUsave[0] = normDU[0]; normDUsave[1]=normDU[1]; normDUsave[2]=normDU[2]; normDUsave[3]=normDU[3];
1542  shortSortDescending(normDU,dUdescendingOrder,4);
1543  /*mwf debug check ordering
1544  for (i = 0; i < nElementBoundaries_element; i++)
1545  {
1546  if (normDU[i+1] > normDU[i])
1547  {
1548  printf("\nPROBLEM Durlofsky out of order eN=%d normDU[%d]=%g > normDU[%d] =%g \n",
1549  eN,i+1,normDU[i+1],
1550  i,normDU[i]);
1551  printf("normDU=[%g,%g,%g,%g] dUdescendingOrder=[%d,%d,%d,%d]\n",
1552  normDU[0],normDU[1],normDU[2],normDU[3],
1553  dUdescendingOrder[0],dUdescendingOrder[1],dUdescendingOrder[2],dUdescendingOrder[3]);
1554  exit(1);
1555  }
1556  if (fabs(normDU[i]-normDUsave[dUdescendingOrder[i]]) > 1.0e-8)
1557  {
1558  printf("\nPROBLEM Durlofsky order wrong eN=%d normDU[%d]=%g normDUsave[%d] =%g \n",
1559  eN,i,normDU[i],dUdescendingOrder[i],normDUsave[dUdescendingOrder[i]]);
1560  printf("normDU=[%g,%g,%g,%g] normDUsave=[%g,%g,%g,%g] dUdescendingOrder=[%d,%d,%d,%d]\n",
1561  normDU[0],normDU[1],normDU[2],normDU[3],
1562  normDUsave[0],normDUsave[1],normDUsave[2],normDUsave[3],
1563  dUdescendingOrder[0],dUdescendingOrder[1],dUdescendingOrder[2],dUdescendingOrder[3]);
1564  exit(1);
1565 
1566  }
1567 
1568  }
1569  mwf debug end */
1570  /*now start checking for overshoot, starting with largest du*/
1571  okGradient = 0; islot = 0;
1572  /*use minimum slope if undershoot/overshoot detected for others*/
1573  while (okGradient == 0 && islot < nElementBoundaries_element)
1574  {
1575  itmp = dUdescendingOrder[islot];
1576  /*start with ok*/
1577  okGradient = 1;
1578  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1579  {
1580  eN_ebN = elementNeighborsArray[eN*nElementBoundaries_element + ebN];
1581  uBar_ebN = elementAverages[eN_ebN];
1582  ebN_global = elementBoundariesArray[eN*nElementBoundaries_element + ebN];
1583  uM = uBar
1584  +
1585  dU[itmp][0]*(elementBoundaryBarycentersArray[ebN_global*3 + 0]-
1586  elementBarycentersArray[eN*3 + 0])
1587  +
1588  dU[itmp][1]*(elementBoundaryBarycentersArray[ebN_global*3 + 1]-
1589  elementBarycentersArray[eN*3 + 1]);
1590  max_ebN = uBar > uBar_ebN ? uBar : uBar_ebN;
1591  min_ebN = uBar < uBar_ebN ? uBar : uBar_ebN;
1592 
1593  okGradient = (okGradient > 0) && min_ebN <= uM -otol && uM <= max_ebN + otol &&
1594  fabs(min_ebN) >= vacuumTol && fabs(max_ebN) >= vacuumTol;
1595  if (enforcePositivity)
1596  {
1597  okGradient = (okGradient > 0) && min_ebN >= vacuumTol && max_ebN >= vacuumTol;
1598  }
1599 
1600  }
1601  islot++;/*undo this later if okGradient==1*/
1602  }
1603  if (okGradient == 1)
1604  islot--;
1605  if ((okGradient || allowMinWithUndershoot > 0) && uBar > vacuumTol)
1606  {
1607  itmp = dUdescendingOrder[islot];
1608  dUlim0 = dU[itmp][0];
1609  dUlim1 = dU[itmp][1];
1610  tag[eN] = islot == 0 ? 1 : 0;
1611  /*mwf debug
1612  printf("Durlofsky eN=%d okGradient islot=%d dUlim[0]=%g dUlim[1]=%g \n",
1613  eN,islot,dUlim[0],dUlim[1]);
1614  */
1615  }
1616 
1617  }/*not an extremum*/
1618  }/*in interior*/
1619  /*interpolate values to nodes assuming correspondence with local dofs*/
1620  nN_global = elementNodesArray[eN*nNodes_element+0];
1621  Uout[DOF0] = uBar
1622  +
1623  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1624  +
1625  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1]);
1626  nN_global = elementNodesArray[eN*nNodes_element+1];
1627  Uout[DOF1] = uBar
1628  +
1629  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1630  +
1631  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1]);
1632  nN_global = elementNodesArray[eN*nNodes_element+2];
1633  Uout[DOF2] = uBar
1634  +
1635  dUlim0*(nodeArray[nN_global*3 + 0] - elementBarycentersArray[eN*3 + 0])
1636  +
1637  dUlim1*(nodeArray[nN_global*3 + 1] - elementBarycentersArray[eN*3 + 1]);
1638 
1639  }/*eN*/
1640 }
1641 
computeElementAveragesP1Lagrange
void computeElementAveragesP1Lagrange(int nElements_global, int nDOF_element, const int *l2g, const double *Uin, double *elementAverages)
Definition: timeIntegration.c:839
mminmod2
double mminmod2(double dU, double dU1, double M, int *tag)
Definition: timeIntegration.c:69
applyDurlofskyDGlimiterP1Lagrange3d
void applyDurlofskyDGlimiterP1Lagrange3d(int killExtrema, int allowMinWithUndershoot, int nElements_global, int nElementBoundaries_element, int nNodes_element, int nSpace, int nDOF_element, const int *elementNeighborsArray, const int *elementBoundariesArray, const int *elementNodesArray, const double *nodeArray, const double *elementBarycentersArray, const double *elementBoundaryBarycentersArray, const double *elementNeighborShapeGradients, const int *l2g, const double *grad_v0, double *elementAverages, int *tag, double *Uin, double *Uout)
Definition: timeIntegration.c:1107
neg
double neg(double a)
Definition: testFMMandFSW.cpp:10
computeLocalGradient3d
void computeLocalGradient3d(double *x0, double *x1, double *x2, double *x3, double *dN0, double *dN1, double *dN2, double *dN3, double u0, double u1, double u2, double u3, double *dU)
Definition: timeIntegration.c:164
computeCockburnDGlimiterArrays2d
void computeCockburnDGlimiterArrays2d(int nElements_global, int nElementBoundaries_element, int nSpace, const int *elementBoundariesArray, const int *elementNeighborsArray, const double *elementBarycentersArray, const double *elementBoundaryBarycentersArray, const double *elementNeighborShapeGradients, double *alphas, int *alphaNeighbors)
Definition: timeIntegration.c:649
applyDGlimitingP1Lagrange1d_withVacuumTol
void applyDGlimitingP1Lagrange1d_withVacuumTol(int enforcePositivity, double vacuumTol, int nElements_global, int nNodes_element, int nElementBoundaries_element, int nDOF_element, int *elementNodesArray, int *elementNeighborsArray, double *nodeArray, double *elementBarycentersArray, int *l2g, int *tag, double *Uin, double *Uout)
Definition: timeIntegration.c:310
psiTCtteDT
void psiTCtteDT(int nPoints, double tau, double dtn, double dtnm1, double *yn, double *ypn, double *ypnm1, double *dtnp1)
Definition: timeIntegration.c:12
msign
double msign(double a)
Definition: timeIntegration.c:49
invertLocal3d
int invertLocal3d(double A[3][3], double AI[3][3])
Definition: timeIntegration.c:80
minmod3
double minmod3(double dU, double dU1, double dU2)
Definition: timeIntegration.c:64
applyDurlofskyDGlimiterP1Lagrange2d_withVacuumTol
void applyDurlofskyDGlimiterP1Lagrange2d_withVacuumTol(int killExtrema, int allowMinWithUndershoot, int enforcePositivity, double vacuumTol, int nElements_global, int nElementBoundaries_element, int nNodes_element, int nSpace, int nDOF_element, const int *elementNeighborsArray, const int *elementBoundariesArray, const int *elementNodesArray, const double *nodeArray, const double *elementBarycentersArray, const double *elementBoundaryBarycentersArray, const double *elementNeighborShapeGradients, const int *l2g, const double *grad_v0, double *elementAverages, int *tag, double *Uin, double *Uout)
Definition: timeIntegration.c:1397
musclmod
double musclmod(double dU, double dU1, double dU2)
Definition: timeIntegration.c:59
shortSortDescending
void shortSortDescending(double *A, int *order, int len)
Definition: timeIntegration.c:862
applyDGlimitingP1Lagrange1d
void applyDGlimitingP1Lagrange1d(int limiterFlag, int nElements_global, int nNodes_element, int nElementBoundaries_element, int nDOF_element, int *elementNodesArray, int *elementNeighborsArray, double *nodeArray, double *elementBarycentersArray, int *l2g, int *tag, double *Uin, double *Uout)
Definition: timeIntegration.c:214
u
Double u
Definition: Headers.h:89
computeLocalGradient1d
void computeLocalGradient1d(double *x0, double *x1, double *dN0, double *dN1, double u0, double u1, double *dU)
Definition: timeIntegration.c:118
applyCockburnDGlimiterP1Lagrange2d
void applyCockburnDGlimiterP1Lagrange2d(double nu, double Mh2, int nElements_global, int nElementBoundaries_element, int nSpace, int nDOF_element, int *elementNeighborsArray, int *l2g, int *tag, double *alphas, int *alphaNeighbors, double *Uin, double *Uout)
Definition: timeIntegration.c:736
computeLocalGradient2d
void computeLocalGradient2d(double *x0, double *x1, double *x2, double *dN0, double *dN1, double *dN2, double u0, double u1, double u2, double *dU)
Definition: timeIntegration.c:128
pos
double pos(double a)
Definition: testFMMandFSW.cpp:8
computeElementNeighborShapeGradients
void computeElementNeighborShapeGradients(int nElements_global, int nElementBoundaries_element, int nSpace, const int *elementBoundariesArray, const int *elementNeighborsArray, double *elementBarycentersArray, double *elementBoundaryBarycentersArray, double *elementNeighborShapeGradients)
Definition: timeIntegration.c:409
timeIntegration.h
c interface for implementation of time discretization algorithms
applyDurlofskyDGlimiterP1Lagrange2d
void applyDurlofskyDGlimiterP1Lagrange2d(int killExtrema, int allowMinWithUndershoot, int nElements_global, int nElementBoundaries_element, int nNodes_element, int nSpace, int nDOF_element, const int *elementNeighborsArray, const int *elementBoundariesArray, const int *elementNodesArray, const double *nodeArray, const double *elementBarycentersArray, const double *elementBoundaryBarycentersArray, const double *elementNeighborShapeGradients, const int *l2g, const double *grad_v0, double *elementAverages, int *tag, double *Uin, double *Uout)
Definition: timeIntegration.c:879
minmod2
double minmod2(double a, double b)
Definition: timeIntegration.c:53
Generated on Fri Jul 1 2022 11:17:46 for proteus by doxygen 1.8.20