proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
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*/
12void 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 ***********************************************************************/
49double msign(double a)
50{
51 return (a < 0.0 ? -1.0 : 1.0);
52}
53double 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}
59double musclmod(double dU, double dU1, double dU2)
60{
61 double tmp = minmod2(dU1,dU2);
62 return minmod2(dU,2.0*tmp);
63}
64double minmod3(double dU, double dU1, double dU2)
65{
66 double tmp = minmod2(dU1,dU2);
67 return minmod2(dU,tmp);
68}
69double 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*/
80int 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 */
118void 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 */
128void 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
164void 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}
214void 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
310void 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
409void 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
649void 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
736void 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}
839void 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*/
862void 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}
879void 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
1107void 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
1397void 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
1643
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
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
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
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
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
musclmod
double musclmod(double dU, double dU1, double dU2)
Definition timeIntegration.c:59
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
mminmod2
double mminmod2(double dU, double dU1, double M, int *tag)
Definition timeIntegration.c:69
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
minmod3
double minmod3(double dU, double dU1, double dU2)
Definition timeIntegration.c:64
computeElementAveragesP1Lagrange
void computeElementAveragesP1Lagrange(int nElements_global, int nDOF_element, const int *l2g, const double *Uin, double *elementAverages)
Definition timeIntegration.c:839
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
shortSortDescending
void shortSortDescending(double *A, int *order, int len)
Definition timeIntegration.c:862
minmod2
double minmod2(double a, double b)
Definition timeIntegration.c:53
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
msign
double msign(double a)
Definition timeIntegration.c:49
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
invertLocal3d
int invertLocal3d(double A[3][3], double AI[3][3])
Definition timeIntegration.c:80
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
Generated on for proteus by doxygen 1.18.0