proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
postprocessing.c
Go to the documentation of this file.
1#include <stdlib.h>
2#include <stdio.h>
3#include <math.h>
4#include <string.h>
5#include <strings.h>
6#include <assert.h>
7#include "postprocessing.h"
8#include PROTEUS_LAPACK_H
13void invertLocal(int nSpace,double A[3][3], double AI[3][3])
14{
15 double detA,detAinv;
16 if (nSpace == 1)
17 {
18 assert(fabs(A[0][0]) > 0.0);
19 AI[0][0] = 1.0/A[0][0];
20 }
21 else if (nSpace == 2)
22 {
23 detA = A[0][0]*A[1][1]-A[0][1]*A[1][0];
24 assert(fabs(detA) > 0.0);
25 detAinv = 1.0/detA;
26 AI[0][0] = detAinv*A[1][1]; AI[0][1] =-detAinv*A[0][1];
27 AI[1][0] =-detAinv*A[1][0]; AI[1][1] = detAinv*A[0][0];
28 }
29 else
30 {
31 detA =
32 A[0][0]*(A[1][1]*A[2][2]-A[2][1]*A[1][2])-
33 A[0][1]*(A[1][0]*A[2][2]-A[2][0]*A[1][2])+
34 A[0][2]*(A[1][0]*A[2][1]-A[2][0]*A[1][1]);
35 assert(fabs(detA) > 0.0);
36 detAinv = 1.0/detA;
37 AI[0][0] = detAinv*(A[1][1]*A[2][2]-A[1][2]*A[2][1]);
38 AI[0][1] = detAinv*(A[0][2]*A[2][1]-A[0][1]*A[2][2]);
39 AI[0][2] = detAinv*(A[0][1]*A[1][2]-A[0][2]*A[1][1]);
40
41 AI[1][0] = detAinv*(A[1][2]*A[2][0]-A[1][0]*A[2][2]);
42 AI[1][1] = detAinv*(A[0][0]*A[2][2]-A[0][2]*A[2][0]);
43 AI[1][2] = detAinv*(A[0][2]*A[1][0]-A[0][0]*A[1][2]);
44
45 AI[2][0] = detAinv*(A[1][0]*A[2][1]-A[1][1]*A[2][0]);
46 AI[2][1] = detAinv*(A[0][1]*A[2][0]-A[0][0]*A[2][1]);
47 AI[2][2] = detAinv*(A[0][0]*A[1][1]-A[0][1]*A[1][0]);
48 }
49}
50
51
52void postProcessRT0velocityFromP1nc(int nElements_global,
53 int nQuadraturePoints_element,
54 int nDOF_test_element,
55 int nElementBoundaries_element,
56 int nQuadraturePoints_elementBoundary,
57 int nSpace,
58 int * nFreeDOF_element,
59 int * freeLocal_element,
60 double * detJ,
61 double * sqrt_det_g,
62 double * n,
63 double * elementBarycenters,
64 double * quad_a,
65 double * quad_f,
66 double * w_dV_r,
67 double * w_dV_m,
68 double * u,
69 double * gradu,
70 double * a,
71 double * f,
72 double * r,
73 double * mt,
74 double * rt0vdofs)
75{
76 /***********************************************************************
77 combine example in Chen '94 paper and Chou and Tang '00 paper to
78 take scalar finite element function, assuming it's from NC_Affine
79 ... space and compute an RT0 velocity.
80
81 This version includes a piecewise constant 'gravity'/advection
82 term, f as well as mass and source terms that may not be computed
83 using L_2 projection. It also likely only holds for diagonal tensor A.
84
85 This means that we lose the correspondence with RT_0 mixed finite
86 element approximation, but should still be able to get a velocity
87 in RT_0 that's locally conservative (and has continuous normal
88 fluxes). Chou and Tang show that their approach should still
89 converge to the correct solution O(h) I believe.
90
91 Velocity (Darcy flux) definition is
92
93 q = -A\grad u + \vec f
94
95 where \vec f wouldy be \rho A\vec g, A is the (diagonal)
96 conductivity \vec g is gravity and \rho is a density.
97
98 A and f need to be evaluated as the 'averages' aka L_2 projections
99 on element-wise constants. Note that in the P^1_nc weak formulation, we get
100
101 (A_h\grad u_h,\grad w_h)_E = (\bar{A}\grad u_h,\grad w_h)_E
102 (f_h,\grad w_h)_E = (\bar{f},grad w_h)_E
103
104 for a diagonal A, since \grad u_h and grad w_h are constant on E.
105
106 Weak formulation for P^1_nc should be
107
108 (m_{h,t},w_h) + \sum_{E}(A_h \grad u_h,\grad w_h) - \sum_{E}(f_h,\grad w_h)_E +
109 (r_h,w_h) = (s_h,w_h) - <q^b,w_h>_{\Gamma_N}
110
111
112 Where, q^b = total flux on Neumann boundary
113
114 To represent the RT0 velocity, we use
115
116 \vec q_h = \vec a_E + b_E\vec x
117
118 which turns out to be on triangle E
119
120 \vec q_h = -\bar{A}_h\grad u_h + \bar{\vec f}_h +
121 (\bar{s}_E-\bar{r}_{E}-\bar{m}_{t})/d(\vec x - \vec x_{E}) + \vec c_{E}
122
123 The scalar b_E is determined exclusively by the mass-conservation constraint:
124
125 b_E = (\bar{s}_E-\bar{r}_{E}-\bar{m}_{t})/d
126
127 and the constant \vec a_{E} comes about through normal flux continuity condition
128
129 \vec a_T = -\bar{A}_{h,E}\grad u_h + \bar{\vec f} - b_E \vec x_{E} + \vec c_E
130
131 To get the correspondence with the RT_0 solution, we have to use L^2 projections for
132 the mass, reaction, and source terms. Otherwise, we can still choose \vec c_E to enforce
133 continuity of the normal fluxes.
134
135 To solve for \vec c_E, we solve a d x d system on each element
136
137 \mat{G}_{E}\vec c_E = \vec j_E
138
139 G_{E,i:} = |e_i|\vec n_i, the unit outer normal for face i times its area/length
140
141 j_{E,i} = (s_h - r_h - m_{h,t},N_i)_{E} - |E|(\bar{s}-\bar{r}-bar{m}_t)/(d+1)
142
143 where the integral (,)_{E} and averages should be computed with
144 the same quadrature used for the P^1_nc solution, N_i is the test
145 function that's 1 on e_i. For Dirichlet boundaries, we set j_{E,i}
146 = 0
147
148
149 In the notation
150 \bar{s}_{E} is the average source over element E
151 \bar{r}_{E} is the average reaction term
152 \bar{m}_{t,E} is the average accumulation term
153 \bar{\vec f}_{E} is the constant advection (gravity) term
154 \bar{A}_{h} is the average diagonal conductivity over E
155 \vec x_{E} is the element barycenter
156 d is the spatial dimension
157
158 assumes that solution holds all degrees of freedom and boundary
159 values.
160
161 Note, r array holds -s + r terms.
162
163 returns array rt0vdofs that's nelements by nd+1, that stores
164 rt0vdofs[eN,:] = [\vec a_E,b_E]
165
166 ***********************************************************************/
167 int eN,ebN,I,J,i,k,ebN_free,ebN_dir;
168 int nSpace2 = nSpace*nSpace;
169 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
170 double volume,area,volFact,areaFact;
171 double ah[3][3] = {{0.0,0.0,0.0},
172 {0.0,0.0,0.0},
173 {0.0,0.0,0.0}};
174 double vecf[3] = {0.0,0.0,0.0};
175 double rh[4] = {0.0,0.0,0.0,0.0}; /*int_{E}r N_i\dx */
176 double mth[4] = {0.0,0.0,0.0,0.0}; /*int_{E}m_t N_i\dx */
177 double rbar,mtbar;
178 double G_E[3][3] = {{0.0,0.0,0.0},
179 {0.0,0.0,0.0},
180 {0.0,0.0,0.0}};
181 double G_Ei[3][3]={{0.0,0.0,0.0},
182 {0.0,0.0,0.0},
183 {0.0,0.0,0.0}};
184 double rhs_c[3] ={0.0,0.0,0.0};
185 double c_E[3] = {0.0,0.0,0.0};
186
187 double dDim = (double) nSpace;
188 double b_E = 0.;
189
190 assert(nDOF_test_element == nSpace+1);
191 volFact = 1.0; areaFact = 1.0;
192 if (nSpace == 2)
193 volFact = 0.5;
194 if (nSpace == 3)
195 {
196 volFact = 1.0/6.0; areaFact = 0.5;
197 }
198 /*mwf debug
199 printf("P1ncV2 MASS CALLED\n");
200 */
201 /*
202 compute average for r,mt, f and A from integration point values
203 */
204 for (eN = 0; eN < nElements_global; eN++)
205 {
206 b_E = 0.0;
207 rh[0] = 0.0; rh[1] = 0.0; rh[2] = 0.0; rh[3] = 0.0; rbar = 0.0;
208 mth[0]= 0.0; mth[1]= 0.0; mth[2] = 0.0; mth[3] = 0.0; mtbar= 0.0;
209 vecf[0]=0.0; vecf[1] = 0.0; vecf[2] = 0.0;
210 ah[0][0] = 0.0; ah[0][1]=0.0; ah[0][2] = 0.0;
211 ah[1][0] = 0.0; ah[1][1]=0.0; ah[1][2] = 0.0;
212 ah[2][0] = 0.0; ah[2][1]=0.0; ah[2][2] = 0.0;
213
214 /*assume affine*/
215 volume = fabs(detJ[eN*nQuadraturePoints_element + 0])*volFact;
216
217 for (k = 0; k < nQuadraturePoints_element; k++)
218 {
219
220 for (I = 0; I < nSpace; I++)
221 {
222 vecf[I] +=
223 f[eN*nQuadraturePoints_element*nSpace +
224 k*nSpace +
225 I]
226 *
227 quad_f[k]
228 *
229 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
230
231 for (J = 0; J < nSpace; J++)
232 {
233 ah[I][J] +=
234 a[eN*nQuadraturePoints_element*nSpace2 +
235 k*nSpace2 +
236 I*nSpace +
237 J]
238 *
239 quad_a[k]
240 *
241 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
242 }/*J*/
243 }/*I*/
244
245 for (i = 0; i < nDOF_test_element; i++)
246 {
247 rh[i] +=
248 r[eN*nQuadraturePoints_element + k]
249 *
250 w_dV_r[eN*nQuadraturePoints_element*nDOF_test_element +
251 k*nDOF_test_element +
252 i];
253
254 mth[i] +=
255 mt[eN*nQuadraturePoints_element + k]
256 *
257 w_dV_m[eN*nQuadraturePoints_element*nDOF_test_element +
258 k*nDOF_test_element +
259 i];
260 }/*i*/
261 }/*k*/
262 /*sum over test functions to get avgs, can't do this on the
263 fly because of more points than just thoughs used in r*/
264 for (i = 0; i < nDOF_test_element; i++)
265 {
266 rbar += rh[i]/volume; mtbar += mth[i]/volume;
267 }
268 b_E = (-rbar - mtbar)/dDim; /*r holds -s + r in notes*/
269 /*b_E*/
270 rt0vdofs[eN*nDOF_RT0V_element + nSpace] = b_E;
271 /*compute base part of constant term as before*/
272 for (I=0; I < nSpace; I++)
273 {
274 rt0vdofs[eN*nDOF_RT0V_element + I] = 0.0;
275 for (J=0; J < nSpace; J++)
276 {
277 /*know that gradu is constant over element*/
278 rt0vdofs[eN*nDOF_RT0V_element + I] -=
279 ah[I][J]
280 *
281 gradu[eN*nQuadraturePoints_element*nSpace +
282 0*nSpace +
283 I];
284 }/*J*/
285 /*advection term*/
286 /*mwf debug
287 printf("rt0 pp eN=%d vecc[%d]=%g \n",eN,I,vecc[I]);
288 */
289 rt0vdofs[eN*nDOF_RT0V_element + I] += vecf[I];
290
291 rt0vdofs[eN*nDOF_RT0V_element + I] -=
292 b_E
293 *
294 elementBarycenters[eN*3 + I]; /*points always 3d*/
295 }/*I*/
296 /*mwf debug
297 printf("v2pp eN=%d b_E=%g b4 corr, a_E= ",eN,b_E);
298 for (I=0; I < nSpace; I++)
299 printf("%g ",rt0vdofs[eN*nDOF_RT0V_element + I]);
300 printf("\n");
301 */
302 /*now compute correction due to not using L_2 projection for
303 mass, reaction, and source terms*/
304 rhs_c[0] = 0.0; rhs_c[1] = 0.0; rhs_c[2] = 0.0;
305 G_E[0][0] = 1.0; G_E[0][1] = 0.0; G_E[0][2] = 0.0;
306 G_E[1][0] = 0.0; G_E[1][1] = 1.0; G_E[1][2] = 0.0;
307 G_E[2][0] = 0.0; G_E[2][1] = 0.0; G_E[2][2] = 1.0;
308
309 ebN_free = nSpace;
310 if (nFreeDOF_element[eN] < ebN_free)
311 ebN_free = nFreeDOF_element[eN];
312 for (i = 0; i < ebN_free; i++)
313 {
314 ebN = freeLocal_element[eN*nDOF_test_element + i];/*should be _trial_element*/
315 /*assumed affine*/
316 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
317 ebN*nQuadraturePoints_elementBoundary + 0]
318 *
319 areaFact;
320 for (I = 0; I < nSpace; I++)
321 G_E[i][I] =
322 area
323 *
324 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
325 ebN*nQuadraturePoints_elementBoundary*nSpace +
326 0*nSpace +
327 I];
328 /*recall r holds -s + r from notes*/
329 rhs_c[i] = -rh[ebN] + volume*rbar/(dDim+1.) -mth[ebN] + mtbar*volume/(dDim+1.);
330 /*mwf debug
331 if (fabs(rhs_c[i]) > 1.0e-16)
332 {
333 printf("v2pp eN=%d rhs_c[%d]= %g ebN=%d rh=%g rbar=%g vol*rbar/(d+1)= %g mth=%g mtbar=%g \n",
334 eN,i,rhs_c[i],ebN,rh[ebN],rbar,volume*rbar/(dDim+1.),mth[ebN],mtbar);
335
336 }
337 */
338 }
339
340 /*mwf debug
341 printf("v2pp eN=%d after nonDir, G_E= \n",eN);
342 for (i=0; i < ebN_free; i++)
343 {
344 for (I = 0; I < nSpace; I++)
345 printf("%g ",G_E[i][I]);
346 printf("\n");
347 }
348 */
349 ebN = 0;
350 while (ebN_free < nSpace && ebN < nElementBoundaries_element)
351 {
352 /*
353 at most d-1 non-Dirichlet boundaries, so have to include
354 dirichlet boundaries with rhs = 0
355 */
356 ebN_dir = 1;
357 for (k = 0; k < ebN_free; k++)
358 {
359 if (ebN == freeLocal_element[eN*nDOF_test_element + k])
360 ebN_dir = 0;
361 }
362 if (ebN_dir > 0)
363 {
364 /*assumed affine*/
365 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
366 ebN*nQuadraturePoints_elementBoundary + 0]
367 *
368 areaFact;
369 for (I = 0; I < nSpace; I++)
370 G_E[ebN_free][I] =
371 area
372 *
373 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
374 ebN*nQuadraturePoints_elementBoundary*nSpace +
375 0*nSpace +
376 I];
377 rhs_c[ebN_free] = 0.0;
378
379 ebN_free++;
380 }
381 ebN++;
382 }/*ebN_free*/
383 assert (ebN_free >= nSpace);
384 /*mwf debug
385 printf("v2pp eN=%d after dir, G_E= \n",eN);
386 for (i=0; i < nSpace; i++)
387 {
388 for (I = 0; I < nSpace; I++)
389 printf("%g ",G_E[i][I]);
390 printf("\n");
391 }
392 */
393 invertLocal(nSpace,G_E,G_Ei);
394 for (I = 0; I < nSpace; I++)
395 {
396 c_E[I] = 0.0;
397 for (J = 0; J < nSpace; J++)
398 c_E[I] += G_Ei[I][J]*rhs_c[J];
399 }
400 /*mwf debug
401 printf("v2pp eN=%d after solve, \n\t rhs_c= ",eN);
402 for (I=0; I < nSpace; I++)
403 printf("%g ",rhs_c[I]);
404 printf("\n\t c_E= ",eN);
405 for (I=0; I < nSpace; I++)
406 printf("%g ",c_E[I]);
407 printf("\n");
408 */
409 /*now correct RT_0 approximation*/
410 for (I = 0; I < nSpace; I++)
411 rt0vdofs[eN*nDOF_RT0V_element + I] += c_E[I];
412 }/*eN*/
413}
414
415void postProcessRT0velocityFromP1nc_sd(int nElements_global,
416 int nQuadraturePoints_element,
417 int nDOF_test_element,
418 int nElementBoundaries_element,
419 int nQuadraturePoints_elementBoundary,
420 int nSpace,
421 int* rowptr,
422 int* colind,
423 int * nFreeDOF_element,
424 int * freeLocal_element,
425 double * detJ,
426 double * sqrt_det_g,
427 double * n,
428 double * elementBarycenters,
429 double * quad_a,
430 double * quad_f,
431 double * w_dV_r,
432 double * w_dV_m,
433 double * u,
434 double * gradu,
435 double * a,
436 double * f,
437 double * r,
438 double * mt,
439 double * rt0vdofs)
440{
441 /***********************************************************************
442 combine example in Chen '94 paper and Chou and Tang '00 paper to
443 take scalar finite element function, assuming it's from NC_Affine
444 ... space and compute an RT0 velocity.
445
446 This version includes a piecewise constant 'gravity'/advection
447 term, f as well as mass and source terms that may not be computed
448 using L_2 projection. It also likely only holds for diagonal tensor A.
449
450 This means that we lose the correspondence with RT_0 mixed finite
451 element approximation, but should still be able to get a velocity
452 in RT_0 that's locally conservative (and has continuous normal
453 fluxes). Chou and Tang show that their approach should still
454 converge to the correct solution O(h) I believe.
455
456 Velocity (Darcy flux) definition is
457
458 q = -A\grad u + \vec f
459
460 where \vec f wouldy be \rho A\vec g, A is the (diagonal)
461 conductivity \vec g is gravity and \rho is a density.
462
463 A and f need to be evaluated as the 'averages' aka L_2 projections
464 on element-wise constants. Note that in the P^1_nc weak formulation, we get
465
466 (A_h\grad u_h,\grad w_h)_E = (\bar{A}\grad u_h,\grad w_h)_E
467 (f_h,\grad w_h)_E = (\bar{f},grad w_h)_E
468
469 for a diagonal A, since \grad u_h and grad w_h are constant on E.
470
471 Weak formulation for P^1_nc should be
472
473 (m_{h,t},w_h) + \sum_{E}(A_h \grad u_h,\grad w_h) - \sum_{E}(f_h,\grad w_h)_E +
474 (r_h,w_h) = (s_h,w_h) - <q^b,w_h>_{\Gamma_N}
475
476
477 Where, q^b = total flux on Neumann boundary
478
479 To represent the RT0 velocity, we use
480
481 \vec q_h = \vec a_E + b_E\vec x
482
483 which turns out to be on triangle E
484
485 \vec q_h = -\bar{A}_h\grad u_h + \bar{\vec f}_h +
486 (\bar{s}_E-\bar{r}_{E}-\bar{m}_{t})/d(\vec x - \vec x_{E}) + \vec c_{E}
487
488 The scalar b_E is determined exclusively by the mass-conservation constraint:
489
490 b_E = (\bar{s}_E-\bar{r}_{E}-\bar{m}_{t})/d
491
492 and the constant \vec a_{E} comes about through normal flux continuity condition
493
494 \vec a_T = -\bar{A}_{h,E}\grad u_h + \bar{\vec f} - b_E \vec x_{E} + \vec c_E
495
496 To get the correspondence with the RT_0 solution, we have to use L^2 projections for
497 the mass, reaction, and source terms. Otherwise, we can still choose \vec c_E to enforce
498 continuity of the normal fluxes.
499
500 To solve for \vec c_E, we solve a d x d system on each element
501
502 \mat{G}_{E}\vec c_E = \vec j_E
503
504 G_{E,i:} = |e_i|\vec n_i, the unit outer normal for face i times its area/length
505
506 j_{E,i} = (s_h - r_h - m_{h,t},N_i)_{E} - |E|(\bar{s}-\bar{r}-bar{m}_t)/(d+1)
507
508 where the integral (,)_{E} and averages should be computed with
509 the same quadrature used for the P^1_nc solution, N_i is the test
510 function that's 1 on e_i. For Dirichlet boundaries, we set j_{E,i}
511 = 0
512
513
514 In the notation
515 \bar{s}_{E} is the average source over element E
516 \bar{r}_{E} is the average reaction term
517 \bar{m}_{t,E} is the average accumulation term
518 \bar{\vec f}_{E} is the constant advection (gravity) term
519 \bar{A}_{h} is the average diagonal conductivity over E
520 \vec x_{E} is the element barycenter
521 d is the spatial dimension
522
523 assumes that solution holds all degrees of freedom and boundary
524 values.
525
526 Note, r array holds -s + r terms.
527
528 returns array rt0vdofs that's nelements by nd+1, that stores
529 rt0vdofs[eN,:] = [\vec a_E,b_E]
530
531 ***********************************************************************/
532 int eN,ebN,I,J,i,k,ebN_free,ebN_dir;
533 int m,nnz=rowptr[nSpace];
534 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
535 double volume,area,volFact,areaFact;
536 double ah[3][3] = {{0.0,0.0,0.0},
537 {0.0,0.0,0.0},
538 {0.0,0.0,0.0}};
539 double vecf[3] = {0.0,0.0,0.0};
540 double rh[4] = {0.0,0.0,0.0,0.0}; /*int_{E}r N_i\dx */
541 double mth[4] = {0.0,0.0,0.0,0.0}; /*int_{E}m_t N_i\dx */
542 double rbar,mtbar;
543 double G_E[3][3] = {{0.0,0.0,0.0},
544 {0.0,0.0,0.0},
545 {0.0,0.0,0.0}};
546 double G_Ei[3][3]={{0.0,0.0,0.0},
547 {0.0,0.0,0.0},
548 {0.0,0.0,0.0}};
549 double rhs_c[3] ={0.0,0.0,0.0};
550 double c_E[3] = {0.0,0.0,0.0};
551
552 double dDim = (double) nSpace;
553 double b_E = 0.;
554
555 assert(nDOF_test_element == nSpace+1);
556 volFact = 1.0; areaFact = 1.0;
557 if (nSpace == 2)
558 volFact = 0.5;
559 if (nSpace == 3)
560 {
561 volFact = 1.0/6.0; areaFact = 0.5;
562 }
563 /*mwf debug
564 printf("P1ncV2 MASS CALLED\n");
565 */
566 /*
567 compute average for r,mt, f and A from integration point values
568 */
569 for (eN = 0; eN < nElements_global; eN++)
570 {
571 b_E = 0.0;
572 rh[0] = 0.0; rh[1] = 0.0; rh[2] = 0.0; rh[3] = 0.0; rbar = 0.0;
573 mth[0]= 0.0; mth[1]= 0.0; mth[2] = 0.0; mth[3] = 0.0; mtbar= 0.0;
574 vecf[0]=0.0; vecf[1] = 0.0; vecf[2] = 0.0;
575 ah[0][0] = 0.0; ah[0][1]=0.0; ah[0][2] = 0.0;
576 ah[1][0] = 0.0; ah[1][1]=0.0; ah[1][2] = 0.0;
577 ah[2][0] = 0.0; ah[2][1]=0.0; ah[2][2] = 0.0;
578
579 /*assume affine*/
580 volume = fabs(detJ[eN*nQuadraturePoints_element + 0])*volFact;
581
582 for (k = 0; k < nQuadraturePoints_element; k++)
583 {
584
585 for (I = 0; I < nSpace; I++)
586 {
587 vecf[I] +=
588 f[eN*nQuadraturePoints_element*nSpace +
589 k*nSpace +
590 I]
591 *
592 quad_f[k]
593 *
594 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
595
596 for (m=rowptr[I];m<rowptr[I+1];m++)
597 {
598 ah[I][colind[m]] +=
599 a[eN*nQuadraturePoints_element*nnz+
600 k*nnz+
601 m]
602 *
603 quad_a[k]
604 *
605 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
606 }/*J*/
607 }/*I*/
608
609 for (i = 0; i < nDOF_test_element; i++)
610 {
611 rh[i] +=
612 r[eN*nQuadraturePoints_element + k]
613 *
614 w_dV_r[eN*nQuadraturePoints_element*nDOF_test_element +
615 k*nDOF_test_element +
616 i];
617
618 mth[i] +=
619 mt[eN*nQuadraturePoints_element + k]
620 *
621 w_dV_m[eN*nQuadraturePoints_element*nDOF_test_element +
622 k*nDOF_test_element +
623 i];
624 }/*i*/
625 }/*k*/
626 /*sum over test functions to get avgs, can't do this on the
627 fly because of more points than just thoughs used in r*/
628 for (i = 0; i < nDOF_test_element; i++)
629 {
630 rbar += rh[i]/volume; mtbar += mth[i]/volume;
631 }
632 b_E = (-rbar - mtbar)/dDim; /*r holds -s + r in notes*/
633 /*b_E*/
634 rt0vdofs[eN*nDOF_RT0V_element + nSpace] = b_E;
635 /*compute base part of constant term as before*/
636 for (I=0; I < nSpace; I++)
637 {
638 rt0vdofs[eN*nDOF_RT0V_element + I] = 0.0;
639 for (J=0; J < nSpace; J++)
640 {
641 /*know that gradu is constant over element*/
642 rt0vdofs[eN*nDOF_RT0V_element + I] -=
643 ah[I][J]
644 *
645 gradu[eN*nQuadraturePoints_element*nSpace +
646 0*nSpace +
647 I];
648 }/*J*/
649 /*advection term*/
650 /*mwf debug
651 printf("rt0 pp eN=%d vecc[%d]=%g \n",eN,I,vecc[I]);
652 */
653 rt0vdofs[eN*nDOF_RT0V_element + I] += vecf[I];
654
655 rt0vdofs[eN*nDOF_RT0V_element + I] -=
656 b_E
657 *
658 elementBarycenters[eN*3 + I]; /*points always 3d*/
659 }/*I*/
660 /*mwf debug
661 printf("v2pp eN=%d b_E=%g b4 corr, a_E= ",eN,b_E);
662 for (I=0; I < nSpace; I++)
663 printf("%g ",rt0vdofs[eN*nDOF_RT0V_element + I]);
664 printf("\n");
665 */
666 /*now compute correction due to not using L_2 projection for
667 mass, reaction, and source terms*/
668 rhs_c[0] = 0.0; rhs_c[1] = 0.0; rhs_c[2] = 0.0;
669 G_E[0][0] = 1.0; G_E[0][1] = 0.0; G_E[0][2] = 0.0;
670 G_E[1][0] = 0.0; G_E[1][1] = 1.0; G_E[1][2] = 0.0;
671 G_E[2][0] = 0.0; G_E[2][1] = 0.0; G_E[2][2] = 1.0;
672
673 ebN_free = nSpace;
674 if (nFreeDOF_element[eN] < ebN_free)
675 ebN_free = nFreeDOF_element[eN];
676 for (i = 0; i < ebN_free; i++)
677 {
678 ebN = freeLocal_element[eN*nDOF_test_element + i];/*should be _trial_element*/
679 /*assumed affine*/
680 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
681 ebN*nQuadraturePoints_elementBoundary + 0]
682 *
683 areaFact;
684 for (I = 0; I < nSpace; I++)
685 G_E[i][I] =
686 area
687 *
688 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
689 ebN*nQuadraturePoints_elementBoundary*nSpace +
690 0*nSpace +
691 I];
692 /*recall r holds -s + r from notes*/
693 rhs_c[i] = -rh[ebN] + volume*rbar/(dDim+1.) -mth[ebN] + mtbar*volume/(dDim+1.);
694 /*mwf debug
695 if (fabs(rhs_c[i]) > 1.0e-16)
696 {
697 printf("v2pp eN=%d rhs_c[%d]= %g ebN=%d rh=%g rbar=%g vol*rbar/(d+1)= %g mth=%g mtbar=%g \n",
698 eN,i,rhs_c[i],ebN,rh[ebN],rbar,volume*rbar/(dDim+1.),mth[ebN],mtbar);
699
700 }
701 */
702 }
703
704 /*mwf debug
705 printf("v2pp eN=%d after nonDir, G_E= \n",eN);
706 for (i=0; i < ebN_free; i++)
707 {
708 for (I = 0; I < nSpace; I++)
709 printf("%g ",G_E[i][I]);
710 printf("\n");
711 }
712 */
713 ebN = 0;
714 while (ebN_free < nSpace && ebN < nElementBoundaries_element)
715 {
716 /*
717 at most d-1 non-Dirichlet boundaries, so have to include
718 dirichlet boundaries with rhs = 0
719 */
720 ebN_dir = 1;
721 for (k = 0; k < ebN_free; k++)
722 {
723 if (ebN == freeLocal_element[eN*nDOF_test_element + k])
724 ebN_dir = 0;
725 }
726 if (ebN_dir > 0)
727 {
728 /*assumed affine*/
729 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
730 ebN*nQuadraturePoints_elementBoundary + 0]
731 *
732 areaFact;
733 for (I = 0; I < nSpace; I++)
734 G_E[ebN_free][I] =
735 area
736 *
737 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
738 ebN*nQuadraturePoints_elementBoundary*nSpace +
739 0*nSpace +
740 I];
741 rhs_c[ebN_free] = 0.0;
742
743 ebN_free++;
744 }
745 ebN++;
746 }/*ebN_free*/
747 assert (ebN_free >= nSpace);
748 /*mwf debug
749 printf("v2pp eN=%d after dir, G_E= \n",eN);
750 for (i=0; i < nSpace; i++)
751 {
752 for (I = 0; I < nSpace; I++)
753 printf("%g ",G_E[i][I]);
754 printf("\n");
755 }
756 */
757 invertLocal(nSpace,G_E,G_Ei);
758 for (I = 0; I < nSpace; I++)
759 {
760 c_E[I] = 0.0;
761 for (J = 0; J < nSpace; J++)
762 c_E[I] += G_Ei[I][J]*rhs_c[J];
763 }
764 /*mwf debug
765 printf("v2pp eN=%d after solve, \n\t rhs_c= ",eN);
766 for (I=0; I < nSpace; I++)
767 printf("%g ",rhs_c[I]);
768 printf("\n\t c_E= ",eN);
769 for (I=0; I < nSpace; I++)
770 printf("%g ",c_E[I]);
771 printf("\n");
772 */
773 /*now correct RT_0 approximation*/
774 for (I = 0; I < nSpace; I++)
775 rt0vdofs[eN*nDOF_RT0V_element + I] += c_E[I];
776 }/*eN*/
777}
778
779void postProcessRT0velocityFromP1ncNoMass(int nElements_global,
780 int nQuadraturePoints_element,
781 int nDOF_test_element,
782 int nElementBoundaries_element,
783 int nQuadraturePoints_elementBoundary,
784 int nSpace,
785 int * nFreeDOF_element,
786 int * freeLocal_element,
787 double * detJ,
788 double * sqrt_det_g,
789 double * n,
790 double * elementBarycenters,
791 double * quad_a,
792 double * quad_f,
793 double * w_dV_r,
794 double * u,
795 double * gradu,
796 double * a,
797 double * f,
798 double * r,
799 double * rt0vdofs)
800{
801 /***********************************************************************
802 version of postProcessRT0velocityFromP1nc without a mass variable
803 ***********************************************************************/
804 int eN,ebN,I,J,i,k,ebN_free,ebN_dir;
805 int nSpace2 = nSpace*nSpace;
806 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
807 double volume,area,volFact,areaFact;
808 double ah[3][3] = {{0.0,0.0,0.0},
809 {0.0,0.0,0.0},
810 {0.0,0.0,0.0}};
811 double vecf[3] = {0.0,0.0,0.0};
812 double rh[4] = {0.0,0.0,0.0,0.0}; /*int_{E}r N_i\dx */
813 double rbar;
814 double G_E[3][3] = {{0.0,0.0,0.0},
815 {0.0,0.0,0.0},
816 {0.0,0.0,0.0}};
817 double G_Ei[3][3]={{0.0,0.0,0.0},
818 {0.0,0.0,0.0},
819 {0.0,0.0,0.0}};
820 double rhs_c[3] ={0.0,0.0,0.0};
821 double c_E[3] = {0.0,0.0,0.0};
822
823 double dDim = (double) nSpace;
824 double b_E = 0.;
825
826 assert(nDOF_test_element == nSpace+1);
827 volFact = 1.0; areaFact = 1.0;
828 if (nSpace == 2)
829 volFact = 0.5;
830 if (nSpace == 3)
831 {
832 volFact = 1.0/6.0; areaFact = 0.5;
833 }
834
835 /*mwf debug
836 printf("P1ncV2 NO MASS CALLED\n");
837 */
838 /*
839 compute average for r, f and A from integration point values
840 */
841
842
843 for (eN = 0; eN < nElements_global; eN++)
844 {
845 b_E = 0.0;
846 rh[0] = 0.0; rh[1] = 0.0; rh[2] = 0.0; rh[3] = 0.0; rbar = 0.0;
847 vecf[0]=0.0; vecf[1] = 0.0; vecf[2] = 0.0;
848 ah[0][0] = 0.0; ah[0][1]=0.0; ah[0][2] = 0.0;
849 ah[1][0] = 0.0; ah[1][1]=0.0; ah[1][2] = 0.0;
850 ah[2][0] = 0.0; ah[2][1]=0.0; ah[2][2] = 0.0;
851
852 /*assume affine*/
853 volume = fabs(detJ[eN*nQuadraturePoints_element + 0])*volFact;
854
855 for (k = 0; k < nQuadraturePoints_element; k++)
856 {
857
858 for (I = 0; I < nSpace; I++)
859 {
860 vecf[I] +=
861 f[eN*nQuadraturePoints_element*nSpace +
862 k*nSpace +
863 I]
864 *
865 quad_f[k]
866 *
867 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
868
869 for (J = 0; J < nSpace; J++)
870 {
871 ah[I][J] +=
872 a[eN*nQuadraturePoints_element*nSpace2 +
873 k*nSpace2 +
874 I*nSpace +
875 J]
876 *
877 quad_a[k]
878 *
879 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
880 }/*J*/
881 }/*I*/
882
883 for (i = 0; i < nDOF_test_element; i++)
884 {
885 rh[i] +=
886 r[eN*nQuadraturePoints_element + k]
887 *
888 w_dV_r[eN*nQuadraturePoints_element*nDOF_test_element +
889 k*nDOF_test_element +
890 i];
891
892 }/*i*/
893 }/*k*/
894 /*sum over test functions to get avgs, can't do this on the
895 fly because of more points than just thoughs used in r*/
896 for (i = 0; i < nDOF_test_element; i++)
897 {
898 rbar += rh[i]/volume;
899 }
900 b_E = -rbar/dDim; /*r holds -s + r in notes*/
901 /*b_E*/
902 rt0vdofs[eN*nDOF_RT0V_element + nSpace] = b_E;
903 /*compute base part of constant term as before*/
904 for (I=0; I < nSpace; I++)
905 {
906 rt0vdofs[eN*nDOF_RT0V_element + I] = 0.0;
907 for (J=0; J < nSpace; J++)
908 {
909 /*know that gradu is constant over element*/
910 rt0vdofs[eN*nDOF_RT0V_element + I] -=
911 ah[I][J]
912 *
913 gradu[eN*nQuadraturePoints_element*nSpace +
914 0*nSpace +
915 I];
916 }/*J*/
917 /*advection term*/
918 /*mwf debug
919 printf("rt0 pp eN=%d vecc[%d]=%g \n",eN,I,vecc[I]);
920 */
921 rt0vdofs[eN*nDOF_RT0V_element + I] += vecf[I];
922
923 rt0vdofs[eN*nDOF_RT0V_element + I] -=
924 b_E
925 *
926 elementBarycenters[eN*3 + I]; /*points always 3d*/
927 }/*I*/
928 /*mwf debug
929 printf("v2pp eN=%d b_E=%g b4 corr, a_E= ",eN,b_E);
930 for (I=0; I < nSpace; I++)
931 printf("%g ",rt0vdofs[eN*nDOF_RT0V_element + I]);
932 printf("\n");
933 */
934 /*now compute correction due to not using L_2 projection for
935 mass, reaction, and source terms*/
936 rhs_c[0] = 0.0; rhs_c[1] = 0.0; rhs_c[2] = 0.0;
937 G_E[0][0] = 1.0; G_E[0][1] = 0.0; G_E[0][2] = 0.0;
938 G_E[1][0] = 0.0; G_E[1][1] = 1.0; G_E[1][2] = 0.0;
939 G_E[2][0] = 0.0; G_E[2][1] = 0.0; G_E[2][2] = 1.0;
940
941 ebN_free = nSpace;
942 if (nFreeDOF_element[eN] < ebN_free)
943 ebN_free = nFreeDOF_element[eN];
944 for (i = 0; i < ebN_free; i++)
945 {
946 ebN = freeLocal_element[eN*nDOF_test_element + i];/*should be _trial_element*/
947 /*assumed affine*/
948 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
949 ebN*nQuadraturePoints_elementBoundary + 0]
950 *
951 areaFact;
952 for (I = 0; I < nSpace; I++)
953 G_E[i][I] =
954 area
955 *
956 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
957 ebN*nQuadraturePoints_elementBoundary*nSpace +
958 0*nSpace +
959 I];
960 /*recall r holds -s + r from notes*/
961 rhs_c[i] = -rh[ebN] + volume*rbar/(dDim+1.);
962 /*mwf debug
963 if (fabs(rhs_c[i]) > 1.0e-16)
964 {
965 printf("v2pp eN=%d rhs_c[%d]= %g ebN=%d rh=%g rbar=%g vol*rbar/(d+1)= %g mth=%g mtbar=%g \n",
966 eN,i,rhs_c[i],ebN,rh[ebN],rbar,volume*rbar/(dDim+1.),mth[ebN],mtbar);
967
968 }
969 */
970 }
971
972 /*mwf debug
973 printf("v2pp eN=%d after nonDir, G_E= \n",eN);
974 for (i=0; i < ebN_free; i++)
975 {
976 for (I = 0; I < nSpace; I++)
977 printf("%g ",G_E[i][I]);
978 printf("\n");
979 }
980 */
981 ebN = 0;
982 while (ebN_free < nSpace && ebN < nElementBoundaries_element)
983 {
984 /*
985 at most d-1 non-Dirichlet boundaries, so have to include
986 dirichlet boundaries with rhs = 0
987 */
988 ebN_dir = 1;
989 for (k = 0; k < ebN_free; k++)
990 {
991 if (ebN == freeLocal_element[eN*nDOF_test_element + k])
992 ebN_dir = 0;
993 }
994 if (ebN_dir > 0)
995 {
996 /*assumed affine*/
997 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
998 ebN*nQuadraturePoints_elementBoundary + 0]
999 *
1000 areaFact;
1001 for (I = 0; I < nSpace; I++)
1002 G_E[ebN_free][I] =
1003 area
1004 *
1005 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
1006 ebN*nQuadraturePoints_elementBoundary*nSpace +
1007 0*nSpace +
1008 I];
1009 rhs_c[ebN_free] = 0.0;
1010
1011 ebN_free++;
1012 }
1013 ebN++;
1014 }/*ebN_free*/
1015 assert (ebN_free >= nSpace);
1016 /*mwf debug
1017 printf("v2pp eN=%d after dir, G_E= \n",eN);
1018 for (i=0; i < nSpace; i++)
1019 {
1020 for (I = 0; I < nSpace; I++)
1021 printf("%g ",G_E[i][I]);
1022 printf("\n");
1023 }
1024 */
1025 invertLocal(nSpace,G_E,G_Ei);
1026 for (I = 0; I < nSpace; I++)
1027 {
1028 c_E[I] = 0.0;
1029 for (J = 0; J < nSpace; J++)
1030 c_E[I] += G_Ei[I][J]*rhs_c[J];
1031 }
1032 /*mwf debug
1033 printf("v2pp eN=%d after solve, \n\t rhs_c= ",eN);
1034 for (I=0; I < nSpace; I++)
1035 printf("%g ",rhs_c[I]);
1036 printf("\n\t c_E= ",eN);
1037 for (I=0; I < nSpace; I++)
1038 printf("%g ",c_E[I]);
1039 printf("\n");
1040 */
1041 /*now correct RT_0 approximation*/
1042 for (I = 0; I < nSpace; I++)
1043 rt0vdofs[eN*nDOF_RT0V_element + I] += c_E[I];
1044 }/*eN*/
1045}
1046
1047void postProcessRT0velocityFromP1ncNoMass_sd(int nElements_global,
1048 int nQuadraturePoints_element,
1049 int nDOF_test_element,
1050 int nElementBoundaries_element,
1051 int nQuadraturePoints_elementBoundary,
1052 int nSpace,
1053 int* rowptr,
1054 int* colind,
1055 int * nFreeDOF_element,
1056 int * freeLocal_element,
1057 double * detJ,
1058 double * sqrt_det_g,
1059 double * n,
1060 double * elementBarycenters,
1061 double * quad_a,
1062 double * quad_f,
1063 double * w_dV_r,
1064 double * u,
1065 double * gradu,
1066 double * a,
1067 double * f,
1068 double * r,
1069 double * rt0vdofs)
1070{
1071 /***********************************************************************
1072 version of postProcessRT0velocityFromP1nc without a mass variable
1073 ***********************************************************************/
1074 int eN,ebN,I,J,i,k,ebN_free,ebN_dir;
1075 int m,nnz=rowptr[nSpace];
1076 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1077 double volume,area,volFact,areaFact;
1078 double ah[3][3] = {{0.0,0.0,0.0},
1079 {0.0,0.0,0.0},
1080 {0.0,0.0,0.0}};
1081 double vecf[3] = {0.0,0.0,0.0};
1082 double rh[4] = {0.0,0.0,0.0,0.0}; /*int_{E}r N_i\dx */
1083 double rbar;
1084 double G_E[3][3] = {{0.0,0.0,0.0},
1085 {0.0,0.0,0.0},
1086 {0.0,0.0,0.0}};
1087 double G_Ei[3][3]={{0.0,0.0,0.0},
1088 {0.0,0.0,0.0},
1089 {0.0,0.0,0.0}};
1090 double rhs_c[3] ={0.0,0.0,0.0};
1091 double c_E[3] = {0.0,0.0,0.0};
1092
1093 double dDim = (double) nSpace;
1094 double b_E = 0.;
1095
1096 assert(nDOF_test_element == nSpace+1);
1097 volFact = 1.0; areaFact = 1.0;
1098 if (nSpace == 2)
1099 volFact = 0.5;
1100 if (nSpace == 3)
1101 {
1102 volFact = 1.0/6.0; areaFact = 0.5;
1103 }
1104
1105 /*mwf debug
1106 printf("P1ncV2 NO MASS CALLED\n");
1107 */
1108 /*
1109 compute average for r, f and A from integration point values
1110 */
1111
1112
1113 for (eN = 0; eN < nElements_global; eN++)
1114 {
1115 b_E = 0.0;
1116 rh[0] = 0.0; rh[1] = 0.0; rh[2] = 0.0; rh[3] = 0.0; rbar = 0.0;
1117 vecf[0]=0.0; vecf[1] = 0.0; vecf[2] = 0.0;
1118 ah[0][0] = 0.0; ah[0][1]=0.0; ah[0][2] = 0.0;
1119 ah[1][0] = 0.0; ah[1][1]=0.0; ah[1][2] = 0.0;
1120 ah[2][0] = 0.0; ah[2][1]=0.0; ah[2][2] = 0.0;
1121
1122 /*assume affine*/
1123 volume = fabs(detJ[eN*nQuadraturePoints_element + 0])*volFact;
1124
1125 for (k = 0; k < nQuadraturePoints_element; k++)
1126 {
1127
1128 for (I = 0; I < nSpace; I++)
1129 {
1130 vecf[I] +=
1131 f[eN*nQuadraturePoints_element*nSpace +
1132 k*nSpace +
1133 I]
1134 *
1135 quad_f[k]
1136 *
1137 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
1138
1139 for (m=rowptr[I];m<rowptr[I+1];m++)
1140 {
1141 ah[I][colind[m]] +=
1142 a[eN*nQuadraturePoints_element*nnz+
1143 k*nnz+
1144 m]
1145 *
1146 quad_a[k]
1147 *
1148 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
1149 }/*J*/
1150 }/*I*/
1151
1152 for (i = 0; i < nDOF_test_element; i++)
1153 {
1154 rh[i] +=
1155 r[eN*nQuadraturePoints_element + k]
1156 *
1157 w_dV_r[eN*nQuadraturePoints_element*nDOF_test_element +
1158 k*nDOF_test_element +
1159 i];
1160
1161 }/*i*/
1162 }/*k*/
1163 /*sum over test functions to get avgs, can't do this on the
1164 fly because of more points than just thoughs used in r*/
1165 for (i = 0; i < nDOF_test_element; i++)
1166 {
1167 rbar += rh[i]/volume;
1168 }
1169 b_E = -rbar/dDim; /*r holds -s + r in notes*/
1170 /*b_E*/
1171 rt0vdofs[eN*nDOF_RT0V_element + nSpace] = b_E;
1172 /*compute base part of constant term as before*/
1173 for (I=0; I < nSpace; I++)
1174 {
1175 rt0vdofs[eN*nDOF_RT0V_element + I] = 0.0;
1176 for (J=0; J < nSpace; J++)
1177 {
1178 /*know that gradu is constant over element*/
1179 rt0vdofs[eN*nDOF_RT0V_element + I] -=
1180 ah[I][J]
1181 *
1182 gradu[eN*nQuadraturePoints_element*nSpace +
1183 0*nSpace +
1184 I];
1185 }/*J*/
1186 /*advection term*/
1187 /*mwf debug
1188 printf("rt0 pp eN=%d vecc[%d]=%g \n",eN,I,vecc[I]);
1189 */
1190 rt0vdofs[eN*nDOF_RT0V_element + I] += vecf[I];
1191
1192 rt0vdofs[eN*nDOF_RT0V_element + I] -=
1193 b_E
1194 *
1195 elementBarycenters[eN*3 + I]; /*points always 3d*/
1196 }/*I*/
1197 /*mwf debug
1198 printf("v2pp eN=%d b_E=%g b4 corr, a_E= ",eN,b_E);
1199 for (I=0; I < nSpace; I++)
1200 printf("%g ",rt0vdofs[eN*nDOF_RT0V_element + I]);
1201 printf("\n");
1202 */
1203 /*now compute correction due to not using L_2 projection for
1204 mass, reaction, and source terms*/
1205 rhs_c[0] = 0.0; rhs_c[1] = 0.0; rhs_c[2] = 0.0;
1206 G_E[0][0] = 1.0; G_E[0][1] = 0.0; G_E[0][2] = 0.0;
1207 G_E[1][0] = 0.0; G_E[1][1] = 1.0; G_E[1][2] = 0.0;
1208 G_E[2][0] = 0.0; G_E[2][1] = 0.0; G_E[2][2] = 1.0;
1209
1210 ebN_free = nSpace;
1211 if (nFreeDOF_element[eN] < ebN_free)
1212 ebN_free = nFreeDOF_element[eN];
1213 for (i = 0; i < ebN_free; i++)
1214 {
1215 ebN = freeLocal_element[eN*nDOF_test_element + i];/*should be _trial_element*/
1216 /*assumed affine*/
1217 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
1218 ebN*nQuadraturePoints_elementBoundary + 0]
1219 *
1220 areaFact;
1221 for (I = 0; I < nSpace; I++)
1222 G_E[i][I] =
1223 area
1224 *
1225 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
1226 ebN*nQuadraturePoints_elementBoundary*nSpace +
1227 0*nSpace +
1228 I];
1229 /*recall r holds -s + r from notes*/
1230 rhs_c[i] = -rh[ebN] + volume*rbar/(dDim+1.);
1231 /*mwf debug
1232 if (fabs(rhs_c[i]) > 1.0e-16)
1233 {
1234 printf("v2pp eN=%d rhs_c[%d]= %g ebN=%d rh=%g rbar=%g vol*rbar/(d+1)= %g mth=%g mtbar=%g \n",
1235 eN,i,rhs_c[i],ebN,rh[ebN],rbar,volume*rbar/(dDim+1.),mth[ebN],mtbar);
1236
1237 }
1238 */
1239 }
1240
1241 /*mwf debug
1242 printf("v2pp eN=%d after nonDir, G_E= \n",eN);
1243 for (i=0; i < ebN_free; i++)
1244 {
1245 for (I = 0; I < nSpace; I++)
1246 printf("%g ",G_E[i][I]);
1247 printf("\n");
1248 }
1249 */
1250 ebN = 0;
1251 while (ebN_free < nSpace && ebN < nElementBoundaries_element)
1252 {
1253 /*
1254 at most d-1 non-Dirichlet boundaries, so have to include
1255 dirichlet boundaries with rhs = 0
1256 */
1257 ebN_dir = 1;
1258 for (k = 0; k < ebN_free; k++)
1259 {
1260 if (ebN == freeLocal_element[eN*nDOF_test_element + k])
1261 ebN_dir = 0;
1262 }
1263 if (ebN_dir > 0)
1264 {
1265 /*assumed affine*/
1266 area = sqrt_det_g[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
1267 ebN*nQuadraturePoints_elementBoundary + 0]
1268 *
1269 areaFact;
1270 for (I = 0; I < nSpace; I++)
1271 G_E[ebN_free][I] =
1272 area
1273 *
1274 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
1275 ebN*nQuadraturePoints_elementBoundary*nSpace +
1276 0*nSpace +
1277 I];
1278 rhs_c[ebN_free] = 0.0;
1279
1280 ebN_free++;
1281 }
1282 ebN++;
1283 }/*ebN_free*/
1284 assert (ebN_free >= nSpace);
1285 /*mwf debug
1286 printf("v2pp eN=%d after dir, G_E= \n",eN);
1287 for (i=0; i < nSpace; i++)
1288 {
1289 for (I = 0; I < nSpace; I++)
1290 printf("%g ",G_E[i][I]);
1291 printf("\n");
1292 }
1293 */
1294 invertLocal(nSpace,G_E,G_Ei);
1295 for (I = 0; I < nSpace; I++)
1296 {
1297 c_E[I] = 0.0;
1298 for (J = 0; J < nSpace; J++)
1299 c_E[I] += G_Ei[I][J]*rhs_c[J];
1300 }
1301 /*mwf debug
1302 printf("v2pp eN=%d after solve, \n\t rhs_c= ",eN);
1303 for (I=0; I < nSpace; I++)
1304 printf("%g ",rhs_c[I]);
1305 printf("\n\t c_E= ",eN);
1306 for (I=0; I < nSpace; I++)
1307 printf("%g ",c_E[I]);
1308 printf("\n");
1309 */
1310 /*now correct RT_0 approximation*/
1311 for (I = 0; I < nSpace; I++)
1312 rt0vdofs[eN*nDOF_RT0V_element + I] += c_E[I];
1313 }/*eN*/
1314}
1315
1316void updateRT0velocityWithAveragedPotentialP1nc(int nElements_global,
1317 int nQuadraturePoints_element,
1318 int nSpace,
1319 double * detJ,
1320 double * quad_a,
1321 double * phi,
1322 double * gradphi,
1323 double * a,
1324 double * rt0vdofs)
1325{
1326 /***********************************************************************
1327 In case have multiple potentials for conservation equation:
1328 -\ten{a}_{j}\grad \phi_j
1329
1330 correct constant part of p1nc RT0 flux with the corresponding
1331 product. This routine computes average for \ten_{a}_j and product
1332 that then goes into constant term for velocity (\vec a_E)
1333
1334 returns array rt0vdofs that's nelements by nd+1, that stores
1335 rt0vdofs[eN,:] = [\vec a_E,b_E]
1336
1337 ***********************************************************************/
1338 int eN,I,J,k;
1339 int nSpace2 = nSpace*nSpace;
1340 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1341 double volume,volFact;
1342 double ah[3][3] = {{0.0,0.0,0.0},
1343 {0.0,0.0,0.0},
1344 {0.0,0.0,0.0}};
1345
1346 volFact = 1.0;
1347 if (nSpace == 2)
1348 volFact = 0.5;
1349 if (nSpace == 3)
1350 volFact = 1.0/6.0;
1351
1352 /*mwf debug
1353 printf("P1ncV2 potential correction CALLED\n");
1354 */
1355 /*
1356 compute average for A from integration point values
1357 */
1358 for (eN = 0; eN < nElements_global; eN++)
1359 {
1360 ah[0][0] = 0.0; ah[0][1]=0.0; ah[0][2] = 0.0;
1361 ah[1][0] = 0.0; ah[1][1]=0.0; ah[1][2] = 0.0;
1362 ah[2][0] = 0.0; ah[2][1]=0.0; ah[2][2] = 0.0;
1363
1364 /*assume affine*/
1365 volume = fabs(detJ[eN*nQuadraturePoints_element + 0])*volFact;
1366
1367 for (k = 0; k < nQuadraturePoints_element; k++)
1368 {
1369
1370 for (I = 0; I < nSpace; I++)
1371 {
1372 for (J = 0; J < nSpace; J++)
1373 {
1374 ah[I][J] +=
1375 a[eN*nQuadraturePoints_element*nSpace2 +
1376 k*nSpace2 +
1377 I*nSpace +
1378 J]
1379 *
1380 quad_a[k]
1381 *
1382 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
1383 }/*J*/
1384 }/*I*/
1385
1386 }/*k*/
1387 /*compute base part of constant term as before*/
1388 for (I=0; I < nSpace; I++)
1389 {
1390 for (J=0; J < nSpace; J++)
1391 {
1392 /*know that gradu is constant over element*/
1393 rt0vdofs[eN*nDOF_RT0V_element + I] -=
1394 ah[I][J]
1395 *
1396 gradphi[eN*nQuadraturePoints_element*nSpace +
1397 0*nSpace +
1398 I];
1399 }/*J*/
1400 }/*I*/
1401 /*mwf debug
1402 printf("v2pp eN=%d b_E=%g b4 corr, a_E= ",eN,b_E);
1403 for (I=0; I < nSpace; I++)
1404 printf("%g ",rt0vdofs[eN*nDOF_RT0V_element + I]);
1405 printf("\n");
1406 */
1407 }/*eN*/
1408}
1409
1410void updateRT0velocityWithAveragedPotentialP1nc_sd(int nElements_global,
1411 int nQuadraturePoints_element,
1412 int nSpace,
1413 int* rowptr,
1414 int* colind,
1415 double * detJ,
1416 double * quad_a,
1417 double * phi,
1418 double * gradphi,
1419 double * a,
1420 double * rt0vdofs)
1421{
1422 /***********************************************************************
1423 In case have multiple potentials for conservation equation:
1424 -\ten{a}_{j}\grad \phi_j
1425
1426 correct constant part of p1nc RT0 flux with the corresponding
1427 product. This routine computes average for \ten_{a}_j and product
1428 that then goes into constant term for velocity (\vec a_E)
1429
1430 returns array rt0vdofs that's nelements by nd+1, that stores
1431 rt0vdofs[eN,:] = [\vec a_E,b_E]
1432
1433 ***********************************************************************/
1434 int eN,I,J,k;
1435 int m,nnz=rowptr[nSpace];
1436 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1437 double volume,volFact;
1438 double ah[3][3] = {{0.0,0.0,0.0},
1439 {0.0,0.0,0.0},
1440 {0.0,0.0,0.0}};
1441
1442 volFact = 1.0;
1443 if (nSpace == 2)
1444 volFact = 0.5;
1445 if (nSpace == 3)
1446 volFact = 1.0/6.0;
1447
1448 /*mwf debug
1449 printf("P1ncV2 potential correction CALLED\n");
1450 */
1451 /*
1452 compute average for A from integration point values
1453 */
1454 for (eN = 0; eN < nElements_global; eN++)
1455 {
1456 ah[0][0] = 0.0; ah[0][1]=0.0; ah[0][2] = 0.0;
1457 ah[1][0] = 0.0; ah[1][1]=0.0; ah[1][2] = 0.0;
1458 ah[2][0] = 0.0; ah[2][1]=0.0; ah[2][2] = 0.0;
1459
1460 /*assume affine*/
1461 volume = fabs(detJ[eN*nQuadraturePoints_element + 0])*volFact;
1462
1463 for (k = 0; k < nQuadraturePoints_element; k++)
1464 {
1465
1466 for (I = 0; I < nSpace; I++)
1467 {
1468 for(m=rowptr[I];m<rowptr[I+1];m++)
1469 {
1470 ah[I][colind[m]] +=
1471 a[eN*nQuadraturePoints_element*nnz+
1472 k*nnz+
1473 m]
1474 *
1475 quad_a[k]
1476 *
1477 fabs(detJ[eN*nQuadraturePoints_element + k])/volume;
1478 }/*J*/
1479 }/*I*/
1480
1481 }/*k*/
1482 /*compute base part of constant term as before*/
1483 for (I=0; I < nSpace; I++)
1484 {
1485 for (J=0; J < nSpace; J++)
1486 {
1487 /*know that gradu is constant over element*/
1488 rt0vdofs[eN*nDOF_RT0V_element + I] -=
1489 ah[I][J]
1490 *
1491 gradphi[eN*nQuadraturePoints_element*nSpace +
1492 0*nSpace +
1493 I];
1494 }/*J*/
1495 }/*I*/
1496 /*mwf debug
1497 printf("v2pp eN=%d b_E=%g b4 corr, a_E= ",eN,b_E);
1498 for (I=0; I < nSpace; I++)
1499 printf("%g ",rt0vdofs[eN*nDOF_RT0V_element + I]);
1500 printf("\n");
1501 */
1502 }/*eN*/
1503}
1504
1505void getElementRT0velocityValues(int nElements_global,
1506 int nPoints_element,
1507 int nSpace,
1508 double * x_element,
1509 double * rt0vdofs_element,
1510 double * v_element)
1511{
1512 /***********************************************************************
1513 Compute \vec q_h at physical points stored on each element
1514
1515 On element T:
1516
1517 \vec q_h = \vec a_T + b_t\vec x
1518
1519 where rt0vdofs_element is logically nElements_global x nSpace+1
1520 assumes
1521 x stored as nElements_global \times nPoints \times 3
1522 v_element stored as nElements_global \times nPoints \times nSpace
1523
1524 ***********************************************************************/
1525 int eN,I,k;
1526 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1527 double b_T = 0.;
1528 for (eN = 0; eN < nElements_global; eN++)
1529 {
1530
1531 b_T = rt0vdofs_element[eN*nDOF_RT0V_element + nSpace];
1532 for (k = 0; k < nPoints_element; k++)
1533 {
1534 for (I = 0; I < nSpace; I++)
1535 {
1536 v_element[eN*nPoints_element*nSpace + k*nSpace + I] =
1537 rt0vdofs_element[eN*nDOF_RT0V_element + I] +
1538 b_T * x_element[eN*nPoints_element*3 + k*3 + I];
1539 }
1540 }/*k*/
1541 }/*eN*/
1542
1543}
1544void getElementBoundaryRT0velocityValues(int nElements_global,
1545 int nElementBoundaries_element,
1546 int nPoints_elementBoundary,
1547 int nSpace,
1548 double * x_elementBoundary,
1549 double * rt0vdofs_element,
1550 double * v_elementBoundary)
1551{
1552 /***********************************************************************
1553 Compute \vec q_h at physical points stored on each element boundary
1554
1555 On element T:
1556
1557 \vec q_h = \vec a_T + b_t\vec x
1558
1559 where rt0vdofs_element is logically nElements_global x nSpace+1
1560 assumes
1561 x stored as nElements_global \times nElementBoundaries_element
1562 \times nPoints \times 3
1563 v_element stored as nElements_global \times nElementBoundaries_element
1564 \times nPoints \times nSpace
1565
1566 ***********************************************************************/
1567 int eN,ebN,I,k;
1568 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1569 double b_T = 0.;
1570 for (eN = 0; eN < nElements_global; eN++)
1571 {
1572
1573 b_T = rt0vdofs_element[eN*nDOF_RT0V_element + nSpace];
1574 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1575 {
1576 for (k = 0; k < nPoints_elementBoundary; k++)
1577 {
1578 for (I = 0; I < nSpace; I++)
1579 {
1580 v_elementBoundary[eN*nElementBoundaries_element*nPoints_elementBoundary*nSpace +
1581 ebN*nPoints_elementBoundary*nSpace +
1582 k*nSpace + I]
1583 =
1584 rt0vdofs_element[eN*nDOF_RT0V_element + I] +
1585 b_T * x_elementBoundary[eN*nElementBoundaries_element*nPoints_elementBoundary*3 +
1586 ebN*nPoints_elementBoundary*3 +
1587 k*3 + I];
1588 }/*I*/
1589 }/*k*/
1590 }/*ebN*/
1591 }/*eN*/
1592
1593}
1594void getGlobalElementBoundaryRT0velocityValues(int nElementBoundaries_global,
1595 int nPoints_elementBoundary,
1596 int nSpace,
1597 int * elementBoundaryElementsArray,
1598 double * x_elementBoundary_global,
1599 double * rt0vdofs_element,
1600 double * v_elementBoundary_global)
1601{
1602 /***********************************************************************
1603 Compute \vec q_h at physical points stored on each global
1604 elementBoundary Just use the left neighbor because it always
1605 exists, and the flux is "known" to be continuous. A good check
1606 would be to compute the values from the right too.
1607
1608 Recall on element T:
1609
1610 \vec q_h = \vec a_T + b_t\vec x
1611
1612 rt0vdofs_element is logically nElements_global x nSpace+1
1613 assumes
1614 x stored as nElementBoundaries_global \times nPoints \times 3
1615 v_elementBoundary_global stored as nElementBoundaries_global \times nPoints \times nSpace
1616
1617 ***********************************************************************/
1618 int ebN,eN,I,k;
1619 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1620 double b_T = 0.;
1621 for (ebN = 0; ebN < nElementBoundaries_global; ebN++)
1622 {
1623 eN = elementBoundaryElementsArray[ebN*2 + 0]; /*left neighbor*/
1624 b_T = rt0vdofs_element[eN*nDOF_RT0V_element + nSpace];
1625 for (k = 0; k < nPoints_elementBoundary; k++)
1626 {
1627 for (I = 0; I < nSpace; I++)
1628 {
1629 v_elementBoundary_global[ebN*nPoints_elementBoundary*nSpace + k*nSpace + I] =
1630 rt0vdofs_element[eN*nDOF_RT0V_element + I] +
1631 b_T * x_elementBoundary_global[ebN*nPoints_elementBoundary*3 + k*3 + I];
1632 }/*I*/
1633 }/*k*/
1634 }/*ebN*/
1635
1636}
1637void getGlobalExteriorElementBoundaryRT0velocityValues(int nExteriorElementBoundaries_global,
1638 int nPoints_elementBoundary,
1639 int nSpace,
1640 int * elementBoundaryElementsArray,
1641 int * exteriorElementBoundariesArray,
1642 double * x_elementBoundary_global,
1643 double * rt0vdofs_element,
1644 double * v_elementBoundary_global)
1645{
1646 /***********************************************************************
1647 Compute \vec q_h at physical points stored on each global exterior
1648 elementBoundary Just use the left neighbor because it always
1649 exists, and the flux is "known" to be continuous. A good check
1650 would be to compute the values from the right too.
1651
1652 Recall on element T:
1653
1654 \vec q_h = \vec a_T + b_t\vec x
1655
1656 rt0vdofs_element is logically nElements_global x nSpace+1
1657 assumes
1658 x stored as nElementBoundaries_global \times nPoints \times 3
1659 v_elementBoundary_global stored as nElementBoundaries_global \times nPoints \times nSpace
1660
1661 ***********************************************************************/
1662 int ebNE,ebN,eN,I,k;
1663 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1664 double b_T = 0.;
1665 for (ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1666 {
1667 ebN = exteriorElementBoundariesArray[ebNE];
1668 eN = elementBoundaryElementsArray[ebN*2 + 0]; /*left neighbor*/
1669 b_T = rt0vdofs_element[eN*nDOF_RT0V_element + nSpace];
1670 for (k = 0; k < nPoints_elementBoundary; k++)
1671 {
1672 for (I = 0; I < nSpace; I++)
1673 {
1674 v_elementBoundary_global[ebNE*nPoints_elementBoundary*nSpace + k*nSpace + I] =
1675 rt0vdofs_element[eN*nDOF_RT0V_element + I] +
1676 b_T * x_elementBoundary_global[ebNE*nPoints_elementBoundary*3 + k*3 + I];
1677 }/*I*/
1678 }/*k*/
1679 }/*ebN*/
1680
1681}
1682
1683void postProcessRT0potentialFromP1nc(int nElements_global,
1684 int nQuadraturePoints_element,
1685 int nElementBoundaries_element,
1686 int nQuadraturePoints_elementBoundary,
1687 int nSpace,
1688 double * uQuadratureWeights_element,
1689 double * elementBarycenters,
1690 double * aElementQuadratureWeights,
1691 double * detJ,
1692 double * uQuadratureWeights_elementBoundary,
1693 double * x,
1694 double * u,
1695 double * gradu,
1696 double * x_elementBoundary,
1697 double * u_elementBoundary,
1698 double * n,
1699 double * a,
1700 double * f,
1701 double * r,
1702 double * rt0vdofs,
1703 double * rt0potential)
1704{
1705 /***********************************************************************
1706 follow example in Chen '94 paper to get RTO potential value.
1707 Assumes RT0 flux has already been calculated!
1708
1709 here, we basically have to plug in the P^1_{nc} solution,
1710 and the postprocessed flux into the local Darcy's law expression
1711 from the mixed hybrid formulation using the test function
1712 \vec v_h= \vec x.
1713 Then we solve for the pressure unknown.
1714
1715 (\mat{A}_h^{-1}\vec q_h,\vec v_h)_T - (p_h,\div \vec v_h)_T
1716 + (\lambda_h,\vec v_h\cdot n_{T})_{\partial T} - (\bar{\vec c},\vec v_h)_T = 0
1717
1718 Recall, that \lambda^{j}_h = p^j_h, the P^1_{nc} solution for edge j
1719
1720 To evaluate the weak integrals, I'll use the quadrature formula
1721
1722 \int_{T} f \dx \approx |T|/3 \sum_{j}f(\vec \bar{x}^j)
1723
1724 for the mass matrix term, and basically exact integration (midpoint rule)
1725 for the boundary integrals of the Lagrange multipler term since it's
1726 linear on each edge, and midpoint rule for the "advection" term since it's linear
1727
1728 ***********************************************************************/
1729 int eN,ebN,I,J,k;
1730 int nSpace2 = nSpace*nSpace;
1731 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1732 double volume,vmass,bndsum,adot,ahqInvI,xdotn,gravsum;
1733 double ah[3][3] = {{0.0,0.0,0.0},
1734 {0.0,0.0,0.0},
1735 {0.0,0.0,0.0}};
1736 double ahInv[3][3] = {{0.0,0.0,0.0},
1737 {0.0,0.0,0.0},
1738 {0.0,0.0,0.0}};
1739 double vecc[3] = {0.0,0.0,0.0};
1740
1741 double dDim = nSpace;
1742 for (eN = 0; eN < nElements_global; eN++)
1743 {
1744 volume = 0.0;
1745 for (k = 0; k < nQuadraturePoints_element; k++)
1746 {
1747 volume += uQuadratureWeights_element[eN*nQuadraturePoints_element + k];
1748 }
1749 /*mwf debug
1750 printf("RT0pot. volume[%d]= %g \n",eN,volume);
1751 */
1752 for (I = 0; I < nSpace; I++)
1753 {
1754 for (J = 0; J < nSpace; J++)
1755 {
1756 /*need to stick to a quad points and weights to be consistent with
1757 original calculation*/
1758 ah[I][J] = 0.0;
1759 for (k = 0; k < nQuadraturePoints_element; k++)
1760 {
1761 ah[I][J] +=
1762 a[eN*nQuadraturePoints_element*nSpace2 +
1763 k*nSpace2 +
1764 I*nSpace +
1765 J]
1766 *
1767 aElementQuadratureWeights[k]*fabs(detJ[eN*nQuadraturePoints_element+k]);
1768 }/*k*/
1769 ah[I][J] = ah[I][J]/volume;
1770 }/*J*/
1771 vecc[I] = 0.0;
1772 for (k = 0; k < nQuadraturePoints_element; k++)
1773 {
1774 vecc[I] +=
1775 f[eN*nQuadraturePoints_element*nSpace +
1776 k*nSpace + I]
1777 *
1778 aElementQuadratureWeights[k]*fabs(detJ[eN*nQuadraturePoints_element+k]);
1779 }/*k*/
1780 vecc[I] = vecc[I]/volume;
1781 }/*I*/
1782 /*mwf debug
1783 printf("RT0pot. before invertLocal ah= \n %g %g \n %g %g \n",ah[0][0],ah[0][1],
1784 ah[1][0],ah[1][1]);
1785 */
1786 invertLocal(nSpace,ah,ahInv);
1787 /*mwf debug
1788 printf("RT0pot. after invertLocal AhInv= \n %g %g \n %g %g \n",ahInv[0][0],ahInv[0][1],
1789 ahInv[1][0],ahInv[1][1]);
1790 */
1791 vmass = 0.0;
1792 for (k = 0; k < nQuadraturePoints_element; k++)
1793 {
1794 adot = 0.0;
1795 for (I = 0; I < nSpace; I++)
1796 {
1797 ahqInvI = 0.0;
1798 for (J = 0; J < nSpace; J++)
1799 ahqInvI+= ahInv[I][J]*(rt0vdofs[eN*nDOF_RT0V_element + J]
1800 +
1801 rt0vdofs[eN*nDOF_RT0V_element + nSpace]
1802 *
1803 x[eN*nQuadraturePoints_element*3 +
1804 k*3 + J]);
1805
1806 adot += ahqInvI*x[eN*nQuadraturePoints_element*3 +
1807 k*3 + I];
1808 }/*I*/
1809 vmass += adot*uQuadratureWeights_element[eN*nQuadraturePoints_element + k];
1810 }
1811
1812 bndsum = 0.0;
1813 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1814 {
1815 for (k = 0; k < nQuadraturePoints_elementBoundary; k++)
1816 {
1817 xdotn = 0.0;
1818 for (I = 0; I < nSpace; I++)
1819 {
1820 xdotn +=
1821 x_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*3 +
1822 ebN*nQuadraturePoints_elementBoundary*3 +
1823 k*3 + I]
1824 *
1825 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
1826 ebN*nQuadraturePoints_elementBoundary*nSpace +
1827 k*nSpace + I];
1828 }
1829 bndsum += u_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
1830 ebN*nQuadraturePoints_elementBoundary + k]
1831 * xdotn
1832 * uQuadratureWeights_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
1833 ebN*nQuadraturePoints_elementBoundary + k];
1834 }
1835 }/*ebN*/
1836 gravsum = 0.0;
1837 for (I=0; I < nSpace; I++)
1838 gravsum -= vecc[I]*elementBarycenters[eN*3 + I]*volume;
1839 rt0potential[eN] = (bndsum + vmass + gravsum)/(dDim*volume);
1840 }/*eN*/
1841}
1842
1843void postProcessRT0potentialFromP1nc_sd(int nElements_global,
1844 int nQuadraturePoints_element,
1845 int nElementBoundaries_element,
1846 int nQuadraturePoints_elementBoundary,
1847 int nSpace,
1848 int* rowptr,
1849 int* colind,
1850 double * uQuadratureWeights_element,
1851 double * elementBarycenters,
1852 double * aElementQuadratureWeights,
1853 double * detJ,
1854 double * uQuadratureWeights_elementBoundary,
1855 double * x,
1856 double * u,
1857 double * gradu,
1858 double * x_elementBoundary,
1859 double * u_elementBoundary,
1860 double * n,
1861 double * a,
1862 double * f,
1863 double * r,
1864 double * rt0vdofs,
1865 double * rt0potential)
1866{
1867 /***********************************************************************
1868 follow example in Chen '94 paper to get RTO potential value.
1869 Assumes RT0 flux has already been calculated!
1870
1871 here, we basically have to plug in the P^1_{nc} solution,
1872 and the postprocessed flux into the local Darcy's law expression
1873 from the mixed hybrid formulation using the test function
1874 \vec v_h= \vec x.
1875 Then we solve for the pressure unknown.
1876
1877 (\mat{A}_h^{-1}\vec q_h,\vec v_h)_T - (p_h,\div \vec v_h)_T
1878 + (\lambda_h,\vec v_h\cdot n_{T})_{\partial T} - (\bar{\vec c},\vec v_h)_T = 0
1879
1880 Recall, that \lambda^{j}_h = p^j_h, the P^1_{nc} solution for edge j
1881
1882 To evaluate the weak integrals, I'll use the quadrature formula
1883
1884 \int_{T} f \dx \approx |T|/3 \sum_{j}f(\vec \bar{x}^j)
1885
1886 for the mass matrix term, and basically exact integration (midpoint rule)
1887 for the boundary integrals of the Lagrange multipler term since it's
1888 linear on each edge, and midpoint rule for the "advection" term since it's linear
1889
1890 ***********************************************************************/
1891 int eN,ebN,I,J,k;
1892 int m,nnz=rowptr[nSpace];
1893 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
1894 double volume,vmass,bndsum,adot,ahqInvI,xdotn,gravsum;
1895 double ah[3][3] = {{0.0,0.0,0.0},
1896 {0.0,0.0,0.0},
1897 {0.0,0.0,0.0}};
1898 double ahInv[3][3] = {{0.0,0.0,0.0},
1899 {0.0,0.0,0.0},
1900 {0.0,0.0,0.0}};
1901 double vecc[3] = {0.0,0.0,0.0};
1902
1903 double dDim = nSpace;
1904 for (eN = 0; eN < nElements_global; eN++)
1905 {
1906 volume = 0.0;
1907 for (k = 0; k < nQuadraturePoints_element; k++)
1908 {
1909 volume += uQuadratureWeights_element[eN*nQuadraturePoints_element + k];
1910 }
1911 /*mwf debug
1912 printf("RT0pot. volume[%d]= %g \n",eN,volume);
1913 */
1914 for (I = 0; I < nSpace; I++)
1915 {
1916 for(m=rowptr[I];m<rowptr[I+1];m++)
1917 {
1918 /*need to stick to a quad points and weights to be consistent with
1919 original calculation*/
1920 ah[I][colind[m]] = 0.0;
1921 for (k = 0; k < nQuadraturePoints_element; k++)
1922 {
1923 ah[I][colind[m]] +=
1924 a[eN*nQuadraturePoints_element*nnz+
1925 k*nnz+
1926 m]
1927 *
1928 aElementQuadratureWeights[k]*fabs(detJ[eN*nQuadraturePoints_element+k]);
1929 }/*k*/
1930 ah[I][colind[m]] = ah[I][colind[m]]/volume;
1931 }/*J*/
1932 vecc[I] = 0.0;
1933 for (k = 0; k < nQuadraturePoints_element; k++)
1934 {
1935 vecc[I] +=
1936 f[eN*nQuadraturePoints_element*nSpace +
1937 k*nSpace + I]
1938 *
1939 aElementQuadratureWeights[k]*fabs(detJ[eN*nQuadraturePoints_element+k]);
1940 }/*k*/
1941 vecc[I] = vecc[I]/volume;
1942 }/*I*/
1943 /*mwf debug
1944 printf("RT0pot. before invertLocal ah= \n %g %g \n %g %g \n",ah[0][0],ah[0][1],
1945 ah[1][0],ah[1][1]);
1946 */
1947 invertLocal(nSpace,ah,ahInv);
1948 /*mwf debug
1949 printf("RT0pot. after invertLocal AhInv= \n %g %g \n %g %g \n",ahInv[0][0],ahInv[0][1],
1950 ahInv[1][0],ahInv[1][1]);
1951 */
1952 vmass = 0.0;
1953 for (k = 0; k < nQuadraturePoints_element; k++)
1954 {
1955 adot = 0.0;
1956 for (I = 0; I < nSpace; I++)
1957 {
1958 ahqInvI = 0.0;
1959 for (J = 0; J < nSpace; J++)
1960 ahqInvI+= ahInv[I][J]*(rt0vdofs[eN*nDOF_RT0V_element + J]
1961 +
1962 rt0vdofs[eN*nDOF_RT0V_element + nSpace]
1963 *
1964 x[eN*nQuadraturePoints_element*3 +
1965 k*3 + J]);
1966
1967 adot += ahqInvI*x[eN*nQuadraturePoints_element*3 +
1968 k*3 + I];
1969 }/*I*/
1970 vmass += adot*uQuadratureWeights_element[eN*nQuadraturePoints_element + k];
1971 }
1972
1973 bndsum = 0.0;
1974 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
1975 {
1976 for (k = 0; k < nQuadraturePoints_elementBoundary; k++)
1977 {
1978 xdotn = 0.0;
1979 for (I = 0; I < nSpace; I++)
1980 {
1981 xdotn +=
1982 x_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*3 +
1983 ebN*nQuadraturePoints_elementBoundary*3 +
1984 k*3 + I]
1985 *
1986 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
1987 ebN*nQuadraturePoints_elementBoundary*nSpace +
1988 k*nSpace + I];
1989 }
1990 bndsum += u_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
1991 ebN*nQuadraturePoints_elementBoundary + k]
1992 * xdotn
1993 * uQuadratureWeights_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
1994 ebN*nQuadraturePoints_elementBoundary + k];
1995 }
1996 }/*ebN*/
1997 gravsum = 0.0;
1998 for (I=0; I < nSpace; I++)
1999 gravsum -= vecc[I]*elementBarycenters[eN*3 + I]*volume;
2000 rt0potential[eN] = (bndsum + vmass + gravsum)/(dDim*volume);
2001 }/*eN*/
2002}
2003
2004void projectElementBoundaryVelocityToRT0fluxRep(int nElements_global,
2005 int nElementBoundaries_element,
2006 int nQuadraturePoints_elementBoundary,
2007 int nSpace,
2008 double * elementBoundaryQuadratureWeights,
2009 double * n,
2010 double * v_elementBoundary,
2011 double * rt0vdofs_element)
2012{
2013 /***********************************************************************
2014 Compute local projection of velocity to RT_0, where the local basis
2015 representation is
2016
2017 \vec N_i = \frac{1}{d|E|}(\vec x - p_i), i=0,...,d
2018
2019 where p_i is the vertex across from face i, |E| is the volume of the element,
2020 and d is the space dimension.
2021
2022 The degrees of freedom are
2023 V^i = \int_{e_i}\vec v\dot n_{i}\ds
2024
2025 Assumes velocity is already consistent so that the normal fluxes are
2026 the same for neighboring elements and that the velocity is stored
2027 in an ebq array of size
2028 nElements x nElementBoundaries_element
2029 x nQuadraturePoints_elementBoundary x nSpace
2030
2031 Uses physical quadrature points on element boundary to calculate flux integral
2032 ***********************************************************************/
2033
2034 int eN,ebN,I,k;
2035 double fluxsum,dotk;
2036 int nDOF_RT0V_element = nSpace+1;
2037 /*mwf debug*/
2038 double area;
2039 for (eN = 0; eN < nElements_global; eN++)
2040 {
2041 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
2042 {
2043 fluxsum = 0.0;
2044 /*mwf debug*/
2045 area = 0.0;
2046 for (k = 0; k < nQuadraturePoints_elementBoundary; k++)
2047 {
2048 dotk = 0.0;
2049 for (I=0; I < nSpace; I++)
2050 {
2051 dotk +=
2052 v_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2053 ebN*nQuadraturePoints_elementBoundary*nSpace+
2054 k*nSpace + I]
2055 *
2056 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2057 ebN*nQuadraturePoints_elementBoundary*nSpace+
2058 k*nSpace + I];
2059 /*mwf debug
2060 printf("getRT0 flux rep dofs v_eb[%d,%d,%d,%d]=%g ; n_eb=%g \n",eN,ebN,k,I,
2061 v_elementBoundary[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2062 ebN*nQuadraturePoints_elementBoundary*nSpace+
2063 k*nSpace + I],
2064 n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2065 ebN*nQuadraturePoints_elementBoundary*nSpace+
2066 k*nSpace + I]);
2067 */
2068
2069 }/*I*/
2070 fluxsum += dotk*elementBoundaryQuadratureWeights[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
2071 ebN*nQuadraturePoints_elementBoundary+
2072 k];
2073 /*mwf debug*/
2074 area += elementBoundaryQuadratureWeights[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
2075 ebN*nQuadraturePoints_elementBoundary+
2076 k];
2077 }/*k*/
2078 rt0vdofs_element[eN*nDOF_RT0V_element + ebN] = fluxsum;
2079 /*mwf debug
2080 printf("rt0vdofs[%d,%d]=%g ; area= %g \n",eN,ebN,fluxsum,area);
2081 */
2082 }/*ebN*/
2083 }/*eN*/
2084}
2085void projectElementBoundaryFluxToRT0fluxRep(int nElements_global,
2086 int nElementBoundaries_element,
2087 int nQuadraturePoints_elementBoundary,
2088 int nDOF_RT0V_element,
2089 int* elementBoundaryElementsArray,
2090 int* elementBoundariesArray,
2091 double * elementBoundaryQuadratureWeights,
2092 double * flux_elementBoundary,
2093 double * rt0vdofs_element)
2094{
2095 /***********************************************************************
2096 Compute local projection of normal flux to RT_0, where the local basis
2097 representation is
2098
2099 \vec N_i = \frac{1}{d|E|}(\vec x - p_i), i=0,...,d
2100
2101 where p_i is the vertex across from face i, |E| is the volume of the element,
2102 and d is the space dimension.
2103
2104 The degrees of freedom are
2105 V^i = \int_{e_i}\vec v\dot n_{i}\ds
2106
2107 Assumes velocity is already consistent so that the normal fluxes are
2108 the same for neighboring elements and that the velocity is stored
2109 in an ebq array of size
2110 nElements x nElementBoundaries_element
2111 x nQuadraturePoints_elementBoundary x nSpace
2112
2113 Uses physical quadrature points on element boundary to calculate flux integral
2114 ***********************************************************************/
2115
2116 int eN,ebN,ebN_global,k;
2117 double fluxsum,sign;
2118 for (eN = 0; eN < nElements_global; eN++)
2119 {
2120 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
2121 {
2122 ebN_global = elementBoundariesArray[eN*nElementBoundaries_element+ebN];
2123 fluxsum = 0.0;
2124 sign=1.0;
2125 if(elementBoundaryElementsArray[2*ebN_global+1] == eN)
2126 sign=-1.0;
2127 for (k = 0; k < nQuadraturePoints_elementBoundary; k++)
2128 {
2129 fluxsum += sign*flux_elementBoundary[ebN_global*nQuadraturePoints_elementBoundary+
2130 k]
2131 *
2132 elementBoundaryQuadratureWeights[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
2133 ebN*nQuadraturePoints_elementBoundary+
2134 k];
2135 }
2136 rt0vdofs_element[eN*nDOF_RT0V_element + ebN] = fluxsum;
2137 /*mwf debug
2138 printf("rt0vdofs[%d,%d]=%g \n",eN,ebN,fluxsum);
2139 */
2140 }/*ebN*/
2141 }/*eN*/
2142}
2143
2144void getElementRT0velocityValuesFluxRep(int nElements_global,
2145 int nElementBoundaries_element,
2146 int nPoints_element,
2147 int nSpace,
2148 int nDetVals_element,
2149 double * nodeArray,
2150 int * elementNodesArray,
2151 double * abs_det_J,
2152 double * x_element,
2153 double * rt0vdofs_element,
2154 double * v_element)
2155{
2156 /***********************************************************************
2157 Compute \vec q_h at physical points stored on each element using the
2158 standard local basis representation
2159
2160 On element T:
2161 \vec q_h = \sum^d_{i=0}V^i\vec N_{T,i}
2162 for
2163 \vec N_{T,i} = \frac{1}{d|E|}(\vec x - p_i), i=0,...,d
2164
2165
2166 where rt0vdofs_element is logically nElements_global x nSpace+1
2167 assumes
2168 x stored as nElements_global \times nPoints \times 3
2169 v_element stored as nElements_global \times nPoints \times nSpace
2170
2171 ***********************************************************************/
2172 int eN,I,k,j,jg;
2173 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
2174 double volFact,dvolInv;
2175 double ddim = nSpace;
2176 double volume = 0.0;
2177 volFact = 1.0;
2178 if (nSpace > 1) volFact = 0.5;
2179 if (nSpace > 2) volFact = 1.0/6.0;
2180
2181 for (eN = 0; eN < nElements_global; eN++)
2182 {
2183 volume = volFact*abs_det_J[eN*nDetVals_element + 0];/*assumed affine*/
2184 assert(volume > 0.0);
2185 dvolInv = 1.0/(ddim*volume);
2186
2187 for (k = 0; k < nPoints_element; k++)
2188 {
2189 for (I = 0; I < nSpace; I++)
2190 {
2191 v_element[eN*nPoints_element*nSpace + k*nSpace + I] = 0.0;
2192 for (j = 0; j < nElementBoundaries_element; j++)
2193 {
2194 jg = elementNodesArray[eN*nElementBoundaries_element + j];
2195 v_element[eN*nPoints_element*nSpace + k*nSpace + I] +=
2196 rt0vdofs_element[eN*nDOF_RT0V_element + j]
2197 *
2198 dvolInv
2199 *(x_element[eN*nPoints_element*3 + k*3 + I]- nodeArray[jg*3 + I]);
2200 /*mwf debug
2201 printf("getRT0 flux rep v[%d,%d,%d]=%g \n",eN,k,I,
2202 v_element[eN*nPoints_element*nSpace + k*nSpace + I]);
2203 */
2204 }/*j*/
2205 }/*I*/
2206 }/*k*/
2207 }/*eN*/
2208}
2209
2210void getElementBoundaryRT0velocityValuesFluxRep(int nElements_global,
2211 int nElementBoundaries_element,
2212 int nPoints_elementBoundary,
2213 int nSpace,
2214 int nDetVals_element,
2215 double * nodeArray,
2216 int * elementNodesArray,
2217 double * abs_det_J,
2218 double * x_elementBoundary,
2219 double * rt0vdofs_element,
2220 double * v_elementBoundary)
2221{
2222 /***********************************************************************
2223 Compute \vec q_h at physical points stored on each local element boundary using the
2224 standard local basis representation
2225
2226 On element T:
2227 \vec q_h = \sum^d_{i=0}V^i\vec N_{T,i}
2228 for
2229 \vec N_{T,i} = \frac{1}{d|E|}(\vec x - p_i), i=0,...,d
2230
2231
2232 where rt0vdofs_element is logically nElements_global x nSpace+1
2233 assumes
2234 x stored as nElements_global \times nPoints \times 3
2235 v_element stored as nElements_global \times nPoints \times nSpace
2236
2237 ***********************************************************************/
2238 int eN,ebN,I,k,j,jg;
2239 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
2240 double volFact,dvolInv;
2241 double ddim = nSpace;
2242 double volume = 0.0;
2243 volFact = 1.0;
2244 if (nSpace > 1) volFact = 0.5;
2245 if (nSpace > 2) volFact = 1.0/6.0;
2246
2247 for (eN = 0; eN < nElements_global; eN++)
2248 {
2249 volume = volFact*abs_det_J[eN*nDetVals_element + 0];/*assumed affine*/
2250 assert(volume > 0.0);
2251 dvolInv = 1.0/(ddim*volume);
2252 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
2253 {
2254 for (k = 0; k < nPoints_elementBoundary; k++)
2255 {
2256 for (I = 0; I < nSpace; I++)
2257 {
2258 v_elementBoundary[eN*nElementBoundaries_element*nPoints_elementBoundary*nSpace +
2259 ebN*nPoints_elementBoundary*nSpace+
2260 k*nSpace + I] = 0.0;
2261 for (j = 0; j < nElementBoundaries_element; j++)
2262 {
2263 jg = elementNodesArray[eN*nElementBoundaries_element + j];
2264 v_elementBoundary[eN*nElementBoundaries_element*nPoints_elementBoundary*nSpace +
2265 ebN*nPoints_elementBoundary*nSpace+
2266 k*nSpace + I] +=
2267 rt0vdofs_element[eN*nDOF_RT0V_element + j]
2268 *
2269 dvolInv
2270 *(x_elementBoundary[eN*nElementBoundaries_element*nPoints_elementBoundary*3 + ebN*nPoints_elementBoundary*3+k*3 + I]- nodeArray[jg*3 + I]);
2271 }/*j*/
2272 }/*I*/
2273 }/*k*/
2274 }/*ebN*/
2275 }/*eN*/
2276}
2277
2278void getGlobalElementBoundaryRT0velocityValuesFluxRep(int nElementBoundaries_global,
2279 int nPoints_elementBoundary_global,
2280 int nSpace,
2281 int nDetVals_element,
2282 double * nodeArray,
2283 int *elementNodesArray,
2284 int *elementBoundaryElementsArray,
2285 double * abs_det_J,
2286 double * x_elementBoundary_global,
2287 double * rt0vdofs_element,
2288 double * v_elementBoundary_global)
2289{
2290 /***********************************************************************
2291 Compute \vec q_h at physical points stored on each local element boundary using the
2292 standard local basis representation
2293
2294 On element T:
2295 \vec q_h = \sum^d_{i=0}V^i\vec N_{T,i}
2296 for
2297 \vec N_{T,i} = \frac{1}{d|E|}(\vec x - p_i), i=0,...,d
2298
2299
2300 where rt0vdofs_element is logically nElements_global x nSpace+1
2301 assumes
2302 x stored as nElements_global \times nPoints \times 3
2303 v_element stored as nElements_global \times nPoints \times nSpace
2304
2305 ***********************************************************************/
2306 int eN,ebN,I,k,j,jg;
2307 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
2308 double volFact,dvolInv;
2309 double ddim = nSpace;
2310 double volume = 0.0;
2311 volFact = 1.0;
2312 if (nSpace > 1) volFact = 0.5;
2313 if (nSpace > 2) volFact = 1.0/6.0;
2314
2315 for (ebN = 0; ebN < nElementBoundaries_global; ebN++)
2316 {
2317 eN = elementBoundaryElementsArray[ebN*2 + 0]; /*left neighbor*/
2318 volume = volFact*abs_det_J[eN*nDetVals_element + 0];/*assumed affine*/
2319 assert(volume > 0.0);
2320 dvolInv = 1.0/(ddim*volume);
2321 for (k = 0; k < nPoints_elementBoundary_global; k++)
2322 {
2323 for (I = 0; I < nSpace; I++)
2324 {
2325 v_elementBoundary_global[ebN*nPoints_elementBoundary_global*nSpace +
2326 k*nSpace +
2327 I] = 0.0;
2328 for (j = 0; j < nDOF_RT0V_element; j++)
2329 {
2330 jg = elementNodesArray[eN*nDOF_RT0V_element + j];
2331 v_elementBoundary_global[ebN*nPoints_elementBoundary_global*nSpace+k*nSpace + I] +=
2332 rt0vdofs_element[eN*nDOF_RT0V_element + j]
2333 *
2334 dvolInv
2335 *(x_elementBoundary_global[ebN*nPoints_elementBoundary_global*3+k*3 + I]- nodeArray[jg*3 + I]);
2336 }/*j*/
2337 }/*I*/
2338 }/*k*/
2339 }/*ebN*/
2340}
2341void getGlobalExteriorElementBoundaryRT0velocityValuesFluxRep(int nExteriorElementBoundaries_global,
2342 int nPoints_elementBoundary_global,
2343 int nSpace,
2344 int nDetVals_element,
2345 double * nodeArray,
2346 int *elementNodesArray,
2347 int *elementBoundaryElementsArray,
2348 int* exteriorElementBoundariesArray,
2349 double * abs_det_J,
2350 double * x_ebqe,
2351 double * rt0vdofs_element,
2352 double * v_ebqe)
2353{
2354 /***********************************************************************
2355 Compute \vec q_h at physical points stored on each local element boundary using the
2356 standard local basis representation
2357
2358 On element T:
2359 \vec q_h = \sum^d_{i=0}V^i\vec N_{T,i}
2360 for
2361 \vec N_{T,i} = \frac{1}{d|E|}(\vec x - p_i), i=0,...,d
2362
2363
2364 where rt0vdofs_element is logically nElements_global x nSpace+1
2365 assumes
2366 x stored as nElements_global \times nPoints \times 3
2367 v_element stored as nElements_global \times nPoints \times nSpace
2368
2369 ***********************************************************************/
2370 int eN,ebN,ebNE,I,k,j,jg;
2371 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
2372 double volFact,dvolInv;
2373 double ddim = nSpace;
2374 double volume = 0.0;
2375 volFact = 1.0;
2376 if (nSpace > 1) volFact = 0.5;
2377 if (nSpace > 2) volFact = 1.0/6.0;
2378 for (ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2379 {
2380 ebN = exteriorElementBoundariesArray[ebNE];
2381 eN = elementBoundaryElementsArray[ebN*2 + 0]; /*left neighbor*/
2382 volume = volFact*abs_det_J[eN*nDetVals_element + 0];/*assumed affine*/
2383 assert(volume > 0.0);
2384 dvolInv = 1.0/(ddim*volume);
2385 for (k = 0; k < nPoints_elementBoundary_global; k++)
2386 {
2387 for (I = 0; I < nSpace; I++)
2388 {
2389 v_ebqe[ebNE*nPoints_elementBoundary_global*nSpace +
2390 k*nSpace +
2391 I] = 0.0;
2392 for (j = 0; j < nDOF_RT0V_element; j++)
2393 {
2394 jg = elementNodesArray[eN*nDOF_RT0V_element + j];
2395 v_ebqe[ebNE*nPoints_elementBoundary_global*nSpace+k*nSpace + I] +=
2396 rt0vdofs_element[eN*nDOF_RT0V_element + j]
2397 *
2398 dvolInv
2399 *(x_ebqe[ebNE*nPoints_elementBoundary_global*3+k*3 + I]- nodeArray[jg*3 + I]);
2400 }/*j*/
2401 }/*I*/
2402 }/*k*/
2403 }/*ebN*/
2404}
2405
2406void getRT0velocityValuesFluxRep_arbitraryElementMembership(int nElements_global,
2407 int nElementBoundaries_element,
2408 int nPoints,
2409 int nSpace,
2410 int nDetVals_element,
2411 const double * nodeArray,
2412 const int * elementNodesArray,
2413 const double * abs_det_J,
2414 const double * x,
2415 const int * element_locations,
2416 const double * rt0vdofs_element,
2417 double * v_element)
2418{
2419 /***********************************************************************
2420 Compute \vec q_h at physical points stored on in x belonging to element_location[x] using the
2421 standard local basis representation
2422
2423 On element T:
2424 \vec q_h = \sum^d_{i=0}V^i\vec N_{T,i}
2425 for
2426 \vec N_{T,i} = \frac{1}{d|E|}(\vec x - p_i), i=0,...,d
2427
2428
2429 where rt0vdofs_element is logically nElements_global x nSpace+1
2430 assumes
2431 x stored as nElements_global \times nPoints \times 3
2432 v_element stored as nElements_global \times nPoints \times nSpace
2433
2434 ***********************************************************************/
2435 int eN,I,k,j,jg;
2436 int nDOF_RT0V_element = nSpace+1; /*number of dofs for vector part of RT0*/
2437 double volFact,dvolInv;
2438 double ddim = nSpace;
2439 double volume = 0.0;
2440 volFact = 1.0;
2441 if (nSpace > 1) volFact = 0.5;
2442 if (nSpace > 2) volFact = 1.0/6.0;
2443 for (k = 0; k < nPoints; k++)
2444 {
2445 eN = element_locations[k];
2446 assert(0 <= eN && eN < nElements_global);
2447 volume = volFact*abs_det_J[eN*nDetVals_element + 0];/*assumed affine*/
2448 assert(volume > 0.0);
2449 dvolInv = 1.0/(ddim*volume);
2450 for (I = 0; I < nSpace; I++)
2451 {
2452 v_element[k*nSpace + I] = 0.0;
2453 for (j = 0; j < nElementBoundaries_element; j++)
2454 {
2455 jg = elementNodesArray[eN*nElementBoundaries_element + j];
2456 v_element[k*nSpace + I] +=
2457 rt0vdofs_element[eN*nDOF_RT0V_element + j]
2458 *
2459 dvolInv
2460 *(x[k*3 + I]- nodeArray[jg*3 + I]);
2461 /*mwf debug
2462 printf("getRT0 flux rep v[%d,%d,%d]=%g \n",eN,k,I,
2463 v_element[eN*nPoints_element*nSpace + k*nSpace + I]);
2464 */
2465 }/*j*/
2466 }/*I*/
2467 }/*k*/
2468}
2469
2470void buildLocalBDM1projectionMatrices_orig(int nElements_global,
2471 int nElementBoundaries_element,
2472 int nQuadraturePoints_elementBoundary,
2473 int nSpace,
2474 int nDOFs_test_element,
2475 int nDOFs_trial_element,
2476 int nVDOFs_element,
2477 double * w_dS_f,
2478 double * ebq_n,
2479 double * ebq_v,
2480 double * BDMprojectionMat_element)
2481{
2482 /***********************************************************************
2483 loop through and build local \f$BDM_1\f$ projection representation for
2484 each element for a simplicial mesh in 2d or 3d. Involves
2485 integration over each face, e, of normal flux weighted by basis
2486 for \f$P^1(e)\f$. Local velocity space is \f$P^1(E)\f$
2487 \f{eqnarray}
2488 P_{ij} &=& \int_{ebN} w_{s} \vec N_j \cdot \vec n_{ebN}\ds \\
2489
2490 i &=& ebN*nd + s (\mbox{local test function index}) \\
2491 ebN &=& 0,\ldots,nd (\mbox{local element boundaries}) \\
2492 s &=& 0,\ldots,nd (\mbox{index into local basis for} P^1(e)) \\
2493 j &=& 0,\ldots,nd*(nd+1)-1 \mbox{index for local basis functions} \\
2494\f{eqnarray}
2495
2496 Assumes local basis functions are defined as
2497\f]
2498 \vec N_j = \lambda_k \vec e_l where k = j % nd+1, l = j/(nd+1)
2499\f]
2500 \f$\lambda_k\f$ is barycentric coordinate associated with node \f$k\f$
2501 \f$\vec e_l\f$ is coordinate vector for axis \f$l\f$.
2502
2503 Assumes nodes \f$k\f$ is corresponds to face \f$k \f$across from it.
2504 ***********************************************************************/
2505
2506 int eN,ebN,s,j,k,l,kp,irow,ibq,nVDOFs_element2,nSimplex;
2507 int TRANSPOSE_FOR_LAPACK=1;
2508 double pval;
2509 nSimplex = nSpace+1;
2510 assert(nVDOFs_element == nSpace*(nSpace+1));
2511 assert(nSimplex == nDOFs_trial_element);
2512 assert(nSimplex == nDOFs_test_element);
2513 nVDOFs_element2 = nVDOFs_element*nVDOFs_element;
2514 /*mwf debug
2515 printf("build local BDM nE= %d nEb=%d nBq=%d nvd=%d nd=%d\n",
2516 nElements_global,nElementBoundaries_element,nQuadraturePoints_elementBoundary,
2517 nVDOFs_element,nSpace);
2518 */
2519
2520 for (eN=0; eN < nElements_global; eN++)
2521 {
2522 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
2523 {
2524 for (s = 0; s < nSpace; s++)
2525 {
2526 irow = ebN*nSpace + s;
2527 for (j = 0; j < nVDOFs_element; j++)
2528 {
2529 k = j % nSimplex;
2530 l = j/nSimplex;
2531 kp= (ebN+s+1) % nSimplex; /*neighbor (s) of node k*/
2532 /*mwf debug
2533 printf("eN=%d ebN=%d s=%d irow=%d j=%d k=%d l=%d kp=%d\n",
2534 eN,ebN,s,irow,j,k,l,kp);
2535 */
2536 if (TRANSPOSE_FOR_LAPACK > 0)
2537 BDMprojectionMat_element[eN*nVDOFs_element2 + irow + j*nVDOFs_element] = 0.0;
2538 else
2539 BDMprojectionMat_element[eN*nVDOFs_element2 + irow*nVDOFs_element + j] = 0.0;
2540 for (ibq = 0; ibq < nQuadraturePoints_elementBoundary; ibq++)
2541 {
2542
2543 pval =
2544 ebq_n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
2545 ebN*nQuadraturePoints_elementBoundary*nSpace+
2546 ibq*nSpace+
2547 l]
2548 *
2549 w_dS_f[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_test_element+
2550 ebN*nQuadraturePoints_elementBoundary*nDOFs_test_element+
2551 ibq*nDOFs_test_element+
2552 kp]
2553 *
2554 ebq_v[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_trial_element+
2555 ebN*nQuadraturePoints_elementBoundary*nDOFs_trial_element+
2556 ibq*nDOFs_trial_element+
2557 k];
2558 if (TRANSPOSE_FOR_LAPACK > 0)
2559 BDMprojectionMat_element[eN*nVDOFs_element2 + irow + j*nVDOFs_element] += pval;
2560 else
2561 BDMprojectionMat_element[eN*nVDOFs_element2 + irow*nVDOFs_element + j] += pval;
2562
2563 }/*ibq*/
2564 }/*j*/
2565 }/*s*/
2566 }/*ebN*/
2567
2568 }/*eN*/
2569
2570}
2571
2572void buildLocalBDM1projectionMatrices(int nElements_global,
2573 int nElementBoundaries_element,
2574 int nQuadraturePoints_elementBoundary,
2575 int nSpace,
2576 int nDOFs_test_element,
2577 int nDOFs_trial_element,
2578 int nVDOFs_element,
2579 double * w_dS_f,
2580 double * ebq_n,
2581 double * ebq_v,
2582 double * BDMprojectionMat_element)
2583{
2584 /***********************************************************************
2585
2586Input Variables
2587
2588 nElements_global - Number of elements in triangulation
2589
2590 nElementBoundaries_element - Number of boundaries per element for
2591 2D triangles this is 3, quarilateral - 4 etc.
2592
2593 nQuadraturePoints_elementBoundary - This is the number of quadrature
2594 points taken along the boundary. This value is typically set in
2595 the numerics file with a flag like quad_order.
2596
2597 nSpace - dimension of the problem (typically 2 or 3)
2598
2599 nDOFs_test_element -
2600
2601 nDOFs_trial_element -
2602
2603 nVDOFs_element - number of velocity DoF per element
2604
2605 loop through and build local \f$BDM_1\f$ projection representation for
2606 each element for a simplicial mesh in 2d or 3d. Involves
2607 integration over each face, e, of normal flux weighted by basis
2608 for \f$P^1(e)\f$. Local velocity space is \f$P^1(E)\f$
2609 \f{eqnarray}
2610 P_{ij} &=& \int_{ebN} w_{s} \vec N_j \cdot \vec n_{ebN}\ds \\
2611
2612 i &=& ebN*nd + s (\mbox{local test function index}) \\
2613 ebN &=& 0,\ldots,nd (\mbox{local element boundaries}) \\
2614 s &=& 0,\ldots,nd-1 (\mbox{index into local basis for} P^1(e)) \\
2615 j &=& 0,\ldots,nd*(nd+1)-1 \mbox{index for local basis functions} \\
2616\f{eqnarray}
2617
2618 Assumes local basis functions are defined as
2619\f]
2620 \vec N_j = \lambda_k \vec e_l where k = j / nd, l = j % nd
2621\f]
2622 \f$\lambda_k\f$ is barycentric coordinate associated with node \f$k\f$
2623 \f$\vec e_l\f$ is coordinate vector for axis \f$l\f$.
2624
2625 Assumes nodes \f$k\f$ is corresponds to face \f$k \f$across from it.
2626 ***********************************************************************/
2627
2628 int eN,ebN,s,j,k,l,kp,irow,ibq,nVDOFs_element2,nSimplex;
2629 int TRANSPOSE_FOR_LAPACK=1;
2630 double pval;
2631 nSimplex = nSpace+1;
2632 assert(nVDOFs_element == nSpace*(nSpace+1));
2633 assert(nSimplex == nDOFs_trial_element);
2634 assert(nSimplex == nDOFs_test_element);
2635 nVDOFs_element2 = nVDOFs_element*nVDOFs_element;
2636 /*mwf debug
2637 printf("build local BDM nE= %d nEb=%d nBq=%d nvd=%d nd=%d\n",
2638 nElements_global,nElementBoundaries_element,nQuadraturePoints_elementBoundary,
2639 nVDOFs_element,nSpace);
2640 */
2641
2642 /* printf("BDM1 test information: \n"); */
2643 /* printf("nVDOFs_element: %d\n",nVDOFs_element); */
2644
2645 for (eN=0; eN < nElements_global; eN++)
2646 {
2647 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
2648 {
2649 for (s = 0; s < nSpace; s++)
2650 {
2651 irow = ebN*nSpace + s;
2652 for (j = 0; j < nVDOFs_element; j++)
2653 {
2654 k = j / nSpace;
2655 l = j % nSpace;
2656 kp= (ebN+s+1) % nSimplex; /*neighbor (s) of node k*/
2657 /*mwf debug
2658 printf("BDM1 new eN=%d ebN=%d s=%d irow=%d j=%d k=%d l=%d kp=%d\n",
2659 eN,ebN,s,irow,j,k,l,kp);
2660 printf("BDMprojectionMat_element: %d\n", eN*nVDOFs_element2 + irow + j*nVDOFs_element);
2661 printf("nQuadraturePoints_elementBoundary: %d \n", nQuadraturePoints_elementBoundary);
2662 */
2663
2664 if (TRANSPOSE_FOR_LAPACK > 0)
2665 BDMprojectionMat_element[eN*nVDOFs_element2 + irow + j*nVDOFs_element] = 0.0;
2666 else
2667 BDMprojectionMat_element[eN*nVDOFs_element2 + irow*nVDOFs_element + j] = 0.0;
2668 for (ibq = 0; ibq < nQuadraturePoints_elementBoundary; ibq++)
2669 {
2670
2671 pval =
2672 ebq_n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
2673 ebN*nQuadraturePoints_elementBoundary*nSpace+
2674 ibq*nSpace+
2675 l]
2676 *
2677 w_dS_f[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_test_element+
2678 ebN*nQuadraturePoints_elementBoundary*nDOFs_test_element+
2679 ibq*nDOFs_test_element+
2680 kp]
2681 *
2682 ebq_v[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_trial_element+
2683 ebN*nQuadraturePoints_elementBoundary*nDOFs_trial_element+
2684 ibq*nDOFs_trial_element+
2685 k];
2686 if (TRANSPOSE_FOR_LAPACK > 0){
2687 BDMprojectionMat_element[eN*nVDOFs_element2 + irow + j*nVDOFs_element] += pval;
2688 }
2689 else
2690 BDMprojectionMat_element[eN*nVDOFs_element2 + irow*nVDOFs_element + j] += pval;
2691
2692 }/*ibq*/
2693 /*
2694 if (TRANSPOSE_FOR_LAPACK > 0)
2695 printf("BDM1 new eN=%d B(%d,%d)=%g\n",
2696 eN,irow,j,BDMprojectionMat_element[eN*nVDOFs_element2 + irow + j*nVDOFs_element]);
2697 else
2698 printf("BDM1 new eN=%d B(%d,%d)=%g\n",
2699 eN,irow,j,BDMprojectionMat_element[eN*nVDOFs_element2 + irow*nVDOFs_element + j]);
2700 */
2701 }/*j*/
2702 }/*s*/
2703 }/*ebN*/
2704
2705 }/*eN*/
2706
2707}
2708
2709
2710void buildLocalBDM2projectionMatrices(int degree,
2711 int nElements_global,
2712 int nElementBoundaries_element,
2713 int nQuadraturePoints_elementBoundary,
2714 int nQuadraturePoints_elementInterior,
2715 int nSpace,
2716 int nDOFs_test_element,
2717 int nDOFs_trial_boundary_element,
2718 int nDOFs_trial_interior_element,
2719 int nVDOFs_element,
2720 int *edgeFlags,
2721 double * w_dS_f,
2722 double * ebq_n,
2723 double * ebq_v,
2724 double * BDMprojectionMat_element,
2725 double * q_basis_vals,
2726 double * w_int_test_grads,
2727 double * w_int_div_free,
2728 double * piola_trial_fun)
2729{
2730 /***********************************************************************
2731 This function builds the LocalBDM2projectionMatrices. This includes
2732 three boundary integrals per edge and three interior degrees of freedom.
2733 NOTE - this function has only been tested for TRANSPOSE_FOR_LAPACK=1.
2734 ***********************************************************************/
2735
2736 int eN,ebN,s,i,j,k,l,kp,irow,ibq,nSimplex;
2737 int dof, dof_edge, boundary_dof;
2738 int num_div_free;
2739 int TRANSPOSE_FOR_LAPACK=1;
2740 double pval, pvalx, pvaly;
2741 nSimplex = nSpace+1;
2742 assert(degree == 2);
2743
2744 int interiorPspace = nDOFs_trial_interior_element;
2745
2746 if (nSpace == 2){
2747 dof = (degree+1)*(degree+2);
2748 dof_edge = degree + 1;
2749 boundary_dof = nElementBoundaries_element*dof_edge;
2750 num_div_free = 1;
2751 }
2752 else if (nSpace == 3){
2753 dof = (degree+1)*(degree+2)*(degree+3) / 2;
2754 dof_edge = degree*(degree+1);
2755 boundary_dof = nElementBoundaries_element*dof_edge;
2756 num_div_free = 3;
2757 }
2758 else {
2759 assert(1 == 0);
2760 }
2761
2762 int interior_dof = dof - boundary_dof;
2763 // Begin populating the projection matrix
2764
2765 // Loop over elements of the triangulation
2766 for (eN=0; eN < nElements_global; eN++)
2767 {
2768 // Boundary Integrals
2769
2770 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
2771 {
2772 for (s = 0; s < dof_edge; s++)
2773 {
2774 irow = ebN*dof_edge + s;
2775
2776 for (j = 0; j < dof; j++)
2777 {
2778
2779 k = j / nSpace;
2780 l = j % nSpace;
2781 kp = edgeFlags[ebN*dof_edge+s];
2782
2783 if (TRANSPOSE_FOR_LAPACK > 0)
2784 BDMprojectionMat_element[eN*dof*dof + irow + j*nVDOFs_element] = 0.;
2785 else
2786 BDMprojectionMat_element[eN*dof*dof + irow*nVDOFs_element + j] = 0.;
2787
2788 for (ibq = 0; ibq < nQuadraturePoints_elementBoundary; ibq++)
2789 {
2790
2791 pval =
2792 ebq_n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace +
2793 ebN*nQuadraturePoints_elementBoundary*nSpace+
2794 ibq*nSpace+
2795 l]
2796 *
2797 w_dS_f[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary* nDOFs_test_element+
2798 ebN*nQuadraturePoints_elementBoundary*nDOFs_test_element+
2799 ibq*nDOFs_test_element+
2800 kp]
2801 *
2802 ebq_v[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_trial_boundary_element+
2803 ebN*nQuadraturePoints_elementBoundary*nDOFs_trial_boundary_element+
2804 ibq*nDOFs_trial_boundary_element+
2805 k];
2806
2807
2808 if (TRANSPOSE_FOR_LAPACK > 0)
2809 BDMprojectionMat_element[eN*dof*dof + irow + j*nVDOFs_element] += pval;
2810 else
2811 BDMprojectionMat_element[eN*dof*dof + irow*nVDOFs_element + j] += pval;
2812
2813 }/*ibq*/
2814 }/*j*/
2815 }/*s*/
2816 }/*ebN*/
2817
2818 // **** Interior Integrals (two interior integrals come from gradients, one comes from div-free element) ****
2819
2820 // **** Gradient Integrals ****
2821 for (s = 0; s < interiorPspace-1; s++){
2822 // Iterate over interior polynomial test space
2823 irow = boundary_dof + s;
2824 for (j=0; j < (dof/nSpace); j++){
2825
2826 // Iterate over trial functions
2827 if (TRANSPOSE_FOR_LAPACK > 0){
2828 for (i = 0; i < nSpace; i++){
2829 BDMprojectionMat_element[eN*dof*dof + irow + j*dof*nSpace + i*dof] = 0.0;
2830 }
2831
2832 }
2833 else {
2834 for (i = 0; i < nSpace; i++){
2835 BDMprojectionMat_element[eN*dof*dof + irow*nVDOFs_element + j + i] = 0.0;
2836 }
2837 }
2838
2839 for (ibq=0; ibq < nQuadraturePoints_elementInterior; ibq++){
2840 // Iterate over quadrature points
2841
2842 for (i = 0; i < nSpace; i++){
2843
2844 if (TRANSPOSE_FOR_LAPACK > 0){
2845
2846 BDMprojectionMat_element[eN*dof*dof + irow + j*dof*nSpace + i*dof] +=
2847
2848 q_basis_vals[eN* (dof/nSpace) *nQuadraturePoints_elementInterior +
2849 ibq*nDOFs_test_element +
2850 j]
2851 * w_int_test_grads[eN*nSpace*interiorPspace*nQuadraturePoints_elementInterior+
2852 s*nSpace +
2853 ibq* nSpace *interiorPspace + i];
2854 }
2855
2856 else {
2857 BDMprojectionMat_element[eN*dof*dof + irow*nVDOFs_element + j] += pval;
2858 BDMprojectionMat_element[eN*dof*dof + irow*nVDOFs_element + j + 1] += pval;
2859 }
2860
2861 }
2862
2863 } /* end ibq */
2864 } /* end j*/
2865 } /* end s */
2866
2867 /* // **** DIV-FREE Integrals **** */
2868 irow = boundary_dof + interiorPspace - 1;
2869 for (j=0; j < dof; j++){
2870
2871 if (TRANSPOSE_FOR_LAPACK > 0){
2872 for (s = 0; s < num_div_free; s++){
2873 BDMprojectionMat_element[eN*dof*dof + (irow + s) + j*nVDOFs_element] = 0.0;
2874 }
2875 }
2876 else
2877 BDMprojectionMat_element[eN*dof*dof + irow*nVDOFs_element + j] = 0.0;
2878
2879 for (ibq=0; ibq<nQuadraturePoints_elementInterior; ibq++)
2880 {
2881
2882 if (TRANSPOSE_FOR_LAPACK > 0){
2883 for(s = 0; s < num_div_free; s++){
2884 for (i=0 ; i < nSpace; i++){
2885 BDMprojectionMat_element[eN*dof*dof + (irow + s) + j*nVDOFs_element] +=
2886
2887 w_int_div_free[eN * nQuadraturePoints_elementInterior * nSpace * num_div_free +
2888 ibq * nSpace * num_div_free +
2889 i * num_div_free +
2890 s]
2891
2892 * piola_trial_fun[eN * nQuadraturePoints_elementInterior * dof * nSpace +
2893 ibq * dof * nSpace +
2894 j * nSpace +
2895 i] ;
2896 }
2897 }
2898 }
2899 else
2900 BDMprojectionMat_element[eN*dof*dof + irow*nVDOFs_element + j] = 0.0;
2901 }
2902 }
2903 }/*eN*/
2904}
2905
2906
2907void factorLocalBDM1projectionMatrices(int nElements_global,
2908 int nVDOFs_element,
2909 double *BDMprojectionMat_element,
2910 int *BDMprojectionMatPivots_element)
2911{
2912 PROTEUS_LAPACK_INTEGER INFO=0;
2913 int eN,i,nVDOFs_element2;
2914 PROTEUS_LAPACK_INTEGER pivots_element[12]; /*maximum size for local space is 3*(3+1)*/
2915 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER) nVDOFs_element);
2916 nVDOFs_element2 = nVDOFs_element*nVDOFs_element;
2917 for (eN = 0; eN < nElements_global; eN++)
2918 {
2919 dgetrf_(&nE_n,
2920 &nE_n,
2921 &BDMprojectionMat_element[eN*nVDOFs_element2],
2922 &nE_n,
2923 &pivots_element[0],
2924 &INFO);
2925 // /*mwf debug
2926 // printf("factor local BDM eN=%d INFO=%d",eN,INFO);
2927 // */
2928 for (i = 0; i < nVDOFs_element; i++)
2929 BDMprojectionMatPivots_element[eN*nVDOFs_element+i] = (int) pivots_element[i];
2930
2931 }
2932
2933}
2934
2935
2936void factorLocalBDM2projectionMatrices(int nElements_global,
2937 int nVDOFs_element,
2938 double *BDMprojectionMat_element,
2939 int *BDMprojectionMatPivots_element)
2940{
2941 PROTEUS_LAPACK_INTEGER INFO=0;
2942 int eN,i,nVDOFs_element2;
2943 PROTEUS_LAPACK_INTEGER pivots_element[30]; /*maximum size for local space is 12*/
2944 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER) nVDOFs_element);
2945 nVDOFs_element2 = nVDOFs_element*nVDOFs_element;
2946 for (eN = 0; eN < nElements_global; eN++)
2947 {
2948 dgetrf_(&nE_n,
2949 &nE_n,
2950 &BDMprojectionMat_element[eN*nVDOFs_element2],
2951 &nE_n,
2952 &pivots_element[0],
2953 &INFO);
2954
2955 // printf("factor local BDM eN=%d INFO=%d\n",eN,INFO);
2956
2957 for (i = 0; i < nVDOFs_element; i++)
2958 BDMprojectionMatPivots_element[eN*nVDOFs_element+i] = (int) pivots_element[i];
2959
2960 }
2961
2962}
2963
2964
2965void solveLocalBDM1projection(int nElements_global,
2966 int nElementBoundaries_element,
2967 int nQuadraturePoints_elementBoundary,
2968 int nSpace,
2969 int nDOFs_test_element,
2970 int nVDOFs_element,
2971 double * BDMprojectionMatFact_element,
2972 int* BDMprojectionMatPivots_element,
2973 double * w_dS_f,
2974 double * ebq_n,
2975 double * ebq_velocity,
2976 double * p1_velocity_dofs)
2977{
2978 /***********************************************************************
2979 build right hand side for projection to BDM1 space and then solve
2980 the local projection system.
2981
2982 Assumes ebq_velocity holds the velocity values at element boundaries
2983 that will be used in the projection. \f$w_dS_f\f$ holds test function values
2984 times surface quadrature weights. The test functions are technically
2985 defined for \f$P^1(E)\f$ but take advantage of the fact that \f$ k+1,k+2 (\mod d+1)\f$
2986 for basis for \f$P^1(e_k)\f$ where \f$e_k\f$ is the face across from node \f$k\f$.
2987
2988
2989 Also assumes local projection has been factored and
2990 BDMprojectionMatPivots holds the pivots.
2991
2992
2993 ***********************************************************************/
2994
2995 PROTEUS_LAPACK_INTEGER INFO=0,NRHS=1;
2996 char TRANS='N';
2997 int eN,ebN,s,irow,kp,ibq,J,nSimplex,nVDOFs_element2;
2998 double btmp;
2999 PROTEUS_LAPACK_INTEGER pivots_element[12]; /*maximum size for local space is 3*(3+1)*/
3000 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER) nVDOFs_element);
3001
3002 nSimplex = nSpace+1;
3003 assert(nVDOFs_element == nSpace*(nSpace+1));
3004 assert(nSimplex == nDOFs_test_element);
3005 assert(BDMprojectionMatPivots_element);
3006 nVDOFs_element2 = nVDOFs_element*nVDOFs_element;
3007
3008 for (eN = 0; eN < nElements_global; eN++)
3009 {
3010 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
3011 {
3012 for (s = 0; s < nSpace; s++)
3013 {
3014 irow = ebN*nSpace + s;
3015 kp = (ebN+s+1) % nSimplex;
3016 btmp = 0.0;
3017 for (ibq = 0; ibq < nQuadraturePoints_elementBoundary; ibq++)
3018 {
3019 for (J = 0; J < nSpace; J++)
3020 {
3021 btmp +=
3022 ebq_n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3023 ebN*nQuadraturePoints_elementBoundary*nSpace+
3024 ibq*nSpace+
3025 J]
3026 *
3027 ebq_velocity[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3028 ebN*nQuadraturePoints_elementBoundary*nSpace+
3029 ibq*nSpace+
3030 J]
3031 * w_dS_f[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_test_element+
3032 ebN*nQuadraturePoints_elementBoundary*nDOFs_test_element+
3033 ibq*nDOFs_test_element+
3034 kp];
3035 }/*J*/
3036 }/*ibq*/
3037 p1_velocity_dofs[eN*nVDOFs_element+irow] = btmp;
3038 }
3039 }/*ebN*/
3040 for (irow = 0; irow < nVDOFs_element; irow++)
3041 pivots_element[irow] = BDMprojectionMatPivots_element[eN*nVDOFs_element+irow];
3042 dgetrs_(&TRANS,
3043 &nE_n,
3044 &NRHS,
3045 &BDMprojectionMatFact_element[eN*nVDOFs_element2],
3046 &nE_n,
3047 &pivots_element[0],
3048 &p1_velocity_dofs[eN*nVDOFs_element],
3049 &nE_n,
3050 &INFO);
3051
3052 }/*eN*/
3053}
3054
3055void buildBDM2rhs(int nElements_global,
3056 int nElementBoundaries_element,
3057 int nQuadraturePoints_elementBoundary,
3058 int nQuadraturePoints_elementInterior,
3059 int nSpace,
3060 int nDOFs_test_element,
3061 int nVDOFs_element,
3062 int nDOFs_trial_interior_element,
3063 double * BDMprojectionMatFact_element,
3064 int* BDMprojectionMatPivots_element,
3065 int *edgeFlags,
3066 double * w_dS_f,
3067 double * ebq_n,
3068 double * w_interior_grads,
3069 double * w_interior_divfree,
3070 double * ebq_velocity,
3071 double * q_velocity,
3072 double * p1_velocity_dofs)
3073{
3074 /***********************************************************************
3075 NOTE - *b represents the constructed righthand side
3076
3077 build right hand side for projection to BDM2 space and then solve
3078 the local projection system.
3079
3080 Assumes ebq_velocity holds the velocity values at element boundaries
3081 that will be used in the projection. \f$w_dS_f\f$ holds test function values
3082 times surface quadrature weights. The test functions are technically
3083 defined for \f$P^1(E)\f$ but take advantage of the fact that \f$ k+1,k+2 (\mod d+1)\f$
3084 for basis for \f$P^1(e_k)\f$ where \f$e_k\f$ is the face across from node \f$k\f$.
3085
3086
3087 Also assumes local projection has been factored and
3088 BDMprojectionMatPivots holds the pivots.
3089
3090
3091 ***********************************************************************/
3092
3093 PROTEUS_LAPACK_INTEGER INFO=0,NRHS=1;
3094 char TRANS='N';
3095 int eN,ebN,s,irow,kp,ibq,j,dof_edge,num_div_free,boundary_dof,dof;
3096 double btmp,pvalx,pvaly;
3097 PROTEUS_LAPACK_INTEGER pivots_element[30]; /*maximum size for local space is ???*/
3098 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER) nVDOFs_element);
3099
3100 // temporary variables...these will be added as inputs
3101 int degree = 2;
3102 int interiorPspace = nDOFs_trial_interior_element;
3103
3104 if (nSpace == 2){
3105 dof = (degree+1)*(degree+2);
3106 dof_edge = degree + 1;
3107 boundary_dof = nElementBoundaries_element*dof_edge;
3108 num_div_free = 1;
3109 }
3110 else if (nSpace == 3){
3111 dof = (degree+1)*(degree+2)*(degree+3) / 2;
3112 dof_edge = degree*(degree+1);
3113 boundary_dof = nElementBoundaries_element*dof_edge;
3114 num_div_free = 3;
3115 }
3116 else {
3117 assert(1 == 0);
3118 }
3119
3120 assert(nVDOFs_element == dof);
3121 assert(BDMprojectionMatPivots_element);
3122
3123 for (eN = 0; eN < nElements_global; eN++)
3124 {
3125 // Boundary Integrals
3126 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
3127 {
3128 for (s = 0; s < dof_edge; s++)
3129 {
3130 irow = ebN*dof_edge + s;
3131 kp = edgeFlags[ebN*dof_edge + s];
3132 btmp = 0.0;
3133
3134 for (ibq = 0; ibq < nQuadraturePoints_elementBoundary; ibq++)
3135 {
3136 for (j = 0; j < nSpace; j++) {
3137
3138 btmp +=
3139 ebq_n[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3140 ebN*nQuadraturePoints_elementBoundary*nSpace+
3141 ibq*nSpace+
3142 j]
3143 *
3144 ebq_velocity[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3145 ebN*nQuadraturePoints_elementBoundary*nSpace+
3146 ibq*nSpace+
3147 j]
3148 * w_dS_f[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_test_element+
3149 ebN*nQuadraturePoints_elementBoundary*nDOFs_test_element+
3150 ibq*nDOFs_test_element+
3151 kp];
3152
3153 }/*J*/
3154 }/*ibq*/
3155 p1_velocity_dofs[eN*nVDOFs_element+irow] = btmp;
3156 }
3157 }/*ebN*/
3158
3159 // interior gradient integrals
3160
3161 for (s = 0; s < interiorPspace-1; s++){
3162 // Iterate over interior polynomial test space
3163 irow = boundary_dof + s;
3164 p1_velocity_dofs[eN*nVDOFs_element+irow] = 0. ;
3165
3166 for (ibq=0; ibq < nQuadraturePoints_elementInterior; ibq++){
3167 // Iterate over quadrature points
3168
3169 for (j = 0 ; j < nSpace; j++){
3170
3171 p1_velocity_dofs[eN*nVDOFs_element+irow] +=
3172
3173 q_velocity[eN*nQuadraturePoints_elementInterior*nSpace +
3174 ibq*nSpace +
3175 j] *
3176 w_interior_grads[eN*nSpace*interiorPspace*nQuadraturePoints_elementInterior +
3177 s* nSpace +
3178 ibq * nSpace * interiorPspace +
3179 j];
3180 }
3181
3182 } /* end ibq */
3183 } /* end s */
3184
3185 // div free elements
3186
3187 irow = boundary_dof + interiorPspace - 1;
3188 for (s = 0 ; s < num_div_free ; s++){
3189 p1_velocity_dofs[eN*nVDOFs_element + (irow + s) ] = 0.0 ;
3190 }
3191
3192 for (ibq = 0; ibq < nQuadraturePoints_elementInterior; ibq++){
3193 for (s = 0 ; s < num_div_free ; s++){
3194 for (j = 0 ; j < nSpace ; j++){
3195
3196 p1_velocity_dofs[eN*nVDOFs_element + (irow + s) ] +=
3197
3198 q_velocity[eN * nQuadraturePoints_elementInterior * nSpace +
3199 ibq * nSpace +
3200 j ] *
3201
3202 w_interior_divfree[eN * nQuadraturePoints_elementInterior * num_div_free * nSpace+
3203 ibq * num_div_free * nSpace +
3204 j * num_div_free +
3205 s ];
3206 }
3207 }
3208 }
3209
3210 }/*eN*/
3211
3212}
3213
3214void solveLocalBDM2projection(int nElements_global,
3215 int nElementBoundaries_element,
3216 int nQuadraturePoints_elementBoundary,
3217 int nSpace,
3218 int nDOFs_test_element,
3219 int nVDOFs_element,
3220 double * BDMprojectionMatFact_element,
3221 int* BDMprojectionMatPivots_element,
3222 double * w_dS_f,
3223 double * ebq_n,
3224 double * w_interior_gradients,
3225 double * q_velocity,
3226 double * ebq_velocity,
3227 double * p1_velocity_dofs)
3228{
3229 /***********************************************************************
3230 Solve the BDM2 projection and save answer in p1_velocity_dofs vector.
3231
3232 ***********************************************************************/
3233
3234 PROTEUS_LAPACK_INTEGER INFO=0,NRHS=1;
3235 char TRANS='N';
3236 int eN,ebN,s,irow,kp,ibq,J,nSimplex,nVDOFs_element2;
3237 double btmp;
3238 PROTEUS_LAPACK_INTEGER pivots_element[30];
3239 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER) nVDOFs_element);
3240
3241 for (eN = 0; eN < nElements_global; eN++)
3242 {
3243
3244 for (irow = 0; irow < nVDOFs_element; irow++)
3245 pivots_element[irow] = BDMprojectionMatPivots_element[eN*nVDOFs_element+irow];
3246
3247 dgetrs_(&TRANS,
3248 &nE_n,
3249 &NRHS,
3250 &BDMprojectionMatFact_element[eN*nVDOFs_element*nVDOFs_element],
3251 &nE_n,
3252 &pivots_element[0],
3253 &p1_velocity_dofs[eN*nVDOFs_element],
3254 &nE_n,
3255 &INFO);
3256
3257 }/*eN*/
3258
3259}
3260
3261
3262void solveLocalBDM1projectionFromFlux(int nElements_global,
3263 int nElementBoundaries_element,
3264 int nQuadraturePoints_elementBoundary,
3265 int nDOFs_test_element,
3266 int nVDOFs_element,
3267 double * BDMprojectionMatFact_element,
3268 int* BDMprojectionMatPivots_element,
3269 int* elementBoundaryElementsArray,
3270 int* elementBoundariesArray,
3271 double * w_dS_f,
3272 double * ebq_global_flux,
3273 double * p1_velocity_dofs)
3274{
3275 /***********************************************************************
3276 build right hand side for projection to BDM1 space and then solve
3277 the local projection system.
3278
3279 Assumes ebq_velocity holds the velocity values at element boundaries
3280 that will be used in the projection. \f$w_dS_f\f$ holds test function values
3281 times surface quadrature weights. The test functions are technically
3282 defined for \f$P^1(E)\f$ but take advantage of the fact that \f$ k+1,k+2 (\mod d+1)\f$
3283 for basis for \f$P^1(e_k)\f$ where \f$e_k\f$ is the face across from node \f$k\f$.
3284
3285
3286 Also assumes local projection has been factored and
3287 BDMprojectionMatPivots holds the pivots.
3288
3289
3290 ***********************************************************************/
3291
3292 PROTEUS_LAPACK_INTEGER INFO=0,NRHS=1;
3293 char TRANS='N';
3294 int eN,ebN,ebN_global,nSpace,s,irow,kp,ibq,nSimplex,nVDOFs_element2;
3295 double btmp,sign;
3296 PROTEUS_LAPACK_INTEGER pivots_element[12]; /*maximum size for local space is 3*(3+1)*/
3297 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER) nVDOFs_element);
3298
3299 nSimplex = nDOFs_test_element;
3300 nSpace = nSimplex - 1;
3301 assert(nVDOFs_element == nSpace*(nSpace+1));
3302 assert(nSimplex == nDOFs_test_element);
3303 assert(BDMprojectionMatPivots_element);
3304 nVDOFs_element2 = nVDOFs_element*nVDOFs_element;
3305
3306 for (eN = 0; eN < nElements_global; eN++)
3307 {
3308 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
3309 {
3310 ebN_global = elementBoundariesArray[eN*nElementBoundaries_element+ebN];
3311 sign=1.0;
3312 if(elementBoundaryElementsArray[2*ebN_global+1] == eN)
3313 sign=-1.0;
3314 for (s = 0; s < nSimplex; s++)
3315 {
3316 irow = ebN*nSpace + s;
3317 kp = (ebN+s+1) % nSimplex;
3318 btmp = 0.0;
3319 for (ibq = 0; ibq < nQuadraturePoints_elementBoundary; ibq++)
3320 {
3321 btmp += sign*ebq_global_flux[ebN_global*nQuadraturePoints_elementBoundary+ibq]
3322 *
3323 w_dS_f[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOFs_test_element+
3324 ebN*nQuadraturePoints_elementBoundary*nDOFs_test_element+
3325 ibq*nDOFs_test_element+
3326 kp];
3327 }/*ibq*/
3328 p1_velocity_dofs[eN*nVDOFs_element+irow] = btmp;
3329 }
3330 }/*ebN*/
3331 for (irow = 0; irow < nVDOFs_element; irow++)
3332 pivots_element[irow] = BDMprojectionMatPivots_element[eN*nVDOFs_element+irow];
3333 dgetrs_(&TRANS,
3334 &nE_n,
3335 &NRHS,
3336 &BDMprojectionMatFact_element[eN*nVDOFs_element2],
3337 &nE_n,
3338 &pivots_element[0],
3339 &p1_velocity_dofs[eN*nVDOFs_element],
3340 &nE_n,
3341 &INFO);
3342
3343 }/*eN*/
3344}
3345
3346void getElementBDM1velocityValuesLagrangeRep_orig(int nElements_global,
3347 int nQuadraturePoints_element,
3348 int nSpace,
3349 int nDOF_trial_element,
3350 int nVDOF_element,
3351 double * q_v, /*scalar P^1 shape fncts*/
3352 double * p1_velocity_dofs,
3353 double * q_velocity)
3354{
3355 /***********************************************************************
3356 Assumes local representation for \f$[P^1(E)]^d\f$ is
3357\f[
3358 \vec N_j= \lambda_k \vec e_{id}
3359\f]
3360\f[
3361 k = j % (nd+1), id = j/(nd+1)
3362\f]
3363 **********************************************************************/
3364 int eN,iq,id,k,j;
3365
3366 for (eN = 0; eN < nElements_global; eN++)
3367 {
3368 for (iq = 0; iq < nQuadraturePoints_element; iq++)
3369 {
3370 for (id = 0; id < nSpace; id++)
3371 {
3372 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] = 0.0;
3373 for (k = 0; k < nSpace+1; k++)
3374 {
3375 j = id*(nSpace+1) + k;
3376 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] +=
3377 q_v[eN*nQuadraturePoints_element*nDOF_trial_element + iq*nDOF_trial_element + k]
3378 *
3379 p1_velocity_dofs[eN*nVDOF_element + j];
3380 }/*k*/
3381 }/*id*/
3382 }/*iq*/
3383 }/*eN*/
3384
3385}
3386
3387void getElementBDM1velocityValuesLagrangeRep(int nElements_global,
3388 int nQuadraturePoints_element,
3389 int nSpace,
3390 int nDOF_trial_element,
3391 int nVDOF_element,
3392 double * q_v, /*scalar P^1 shape fncts*/
3393 double * p1_velocity_dofs,
3394 double * q_velocity)
3395{
3396 /***********************************************************************
3397 Assumes local representation for \f$[P^1(E)]^d\f$ is
3398\f[
3399 \vec N_j= \lambda_k \vec e_{id}
3400\f]
3401\f[
3402 k = j / nd, id = j % nd
3403\f]
3404 **********************************************************************/
3405 int eN,iq,id,k,j;
3406
3407 for (eN = 0; eN < nElements_global; eN++)
3408 {
3409 for (iq = 0; iq < nQuadraturePoints_element; iq++)
3410 {
3411 for (id = 0; id < nSpace; id++)
3412 {
3413 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] = 0.0;
3414 for (k = 0; k < nSpace+1; k++)
3415 {
3416 j = k*nSpace+ id; /*id*(nSpace+1) + k;*/
3417 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] +=
3418 q_v[eN*nQuadraturePoints_element*nDOF_trial_element + iq*nDOF_trial_element + k]
3419 *
3420 p1_velocity_dofs[eN*nVDOF_element + j];
3421 }/*k*/
3422 }/*id*/
3423 }/*iq*/
3424 }/*eN*/
3425
3426}
3427
3428void getElementBDM2velocityValuesLagrangeRep(int nElements_global,
3429 int nQuadraturePoints_element,
3430 int nSpace,
3431 int nDOF_trial_element,
3432 int nVDOF_element,
3433 double * q_v, /*scalar P^1 shape fncts*/
3434 double * p1_velocity_dofs,
3435 double * q_velocity)
3436{
3437 /***********************************************************************
3438 Assumes local representation for \f$[P^1(E)]^d\f$ is
3439\f[
3440 \vec N_j= \lambda_k \vec e_{id}
3441\f]
3442\f[
3443 k = j / nd, id = j % nd
3444\f]
3445 **********************************************************************/
3446 int eN,iq,id,k,j;
3447 for (eN = 0; eN < nElements_global; eN++)
3448 {
3449 for (iq = 0; iq < nQuadraturePoints_element; iq++)
3450 {
3451 for (id = 0; id < nSpace; id++)
3452 {
3453 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] = 0.0;
3454 for (k = 0; k < nDOF_trial_element; k++)
3455 {
3456 j = k*nSpace+ id; /*id*(nSpace+1) + k;*/
3457 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] +=
3458 q_v[eN*nQuadraturePoints_element*nDOF_trial_element + iq*nDOF_trial_element + k]
3459 *
3460 p1_velocity_dofs[eN*nVDOF_element + j];
3461
3462 }/*k*/
3463 }/*id*/
3464 }/*iq*/
3465 }/*eN*/
3466
3467}
3468
3469void getElementLDGvelocityValuesLagrangeRep(int nElements_global,
3470 int nQuadraturePoints_element,
3471 int nSpace,
3472 int nDOF_trial_element,
3473 int nVDOF_element,
3474 double * q_v, /*scalar shape fncts*/
3475 double * velocity_dofs,
3476 double * q_velocity)
3477{
3478 /***********************************************************************
3479 Assumes local representation for \f$[P^k(E)]^d\f$ where k=1 or 2
3480
3481 **********************************************************************/
3482 int eN,iq,id,k,j;
3483
3484 for (eN = 0; eN < nElements_global; eN++)
3485 {
3486 for (iq = 0; iq < nQuadraturePoints_element; iq++)
3487 {
3488 for (id = 0; id < nSpace; id++)
3489 {
3490 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] = 0.0;
3491
3492 for (k=0; k < nDOF_trial_element; k++)
3493 /*for (k = 0; k < nSpace+1; k++)*/
3494 {
3495 j = k*nSpace+ id; /*id*(nSpace+1) + k;*/
3496 q_velocity[eN*nQuadraturePoints_element*nSpace + iq*nSpace + id] +=
3497 q_v[eN*nQuadraturePoints_element*nDOF_trial_element + iq*nDOF_trial_element + k]
3498 *
3499 velocity_dofs[eN*nVDOF_element + j];
3500 }/*k*/
3501 }/*id*/
3502 }/*iq*/
3503 }/*eN*/
3504
3505}
3506
3507void getGlobalExteriorElementBoundaryBDM1velocityValuesLagrangeRep(int nExteriorElementBoundaries_global,
3508 int nQuadraturePoints_elementBoundary,
3509 int nSpace,
3510 int nDOF_trial_element,
3511 int nVDOF_element,
3512 int *elementBoundaryElementsArray,
3513 int *exteriorElementBoundariesArray,
3514 double * ebqe_v, /*scalar P^1 shape fncts*/
3515 double * p1_velocity_dofs,
3516 double * ebqe_velocity)
3517{
3518 /***********************************************************************
3519 Assumes local representation for \f$[P^1(E)]^d\f$ is
3520\f[
3521 \vec N_j= \lambda_k \vec e_{id}
3522\f]
3523\f[
3524 k = j / nd, id = j % nd
3525\f]
3526 **********************************************************************/
3527 int ebN,ebNE,eN,iq,id,k,j;
3528
3529 for (ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3530 {
3531 ebN = exteriorElementBoundariesArray[ebNE];
3532 eN = elementBoundaryElementsArray[ebN*2 + 0];
3533
3534 for (iq = 0; iq < nQuadraturePoints_elementBoundary; iq++)
3535 {
3536 for (id = 0; id < nSpace; id++)
3537 {
3538 ebqe_velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace + iq*nSpace + id] = 0.0;
3539 for (k = 0; k < nSpace+1; k++)
3540 {
3541 j = k*nSpace+ id; /*id*(nSpace+1) + k;*/
3542 ebqe_velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace + iq*nSpace + id] +=
3543 ebqe_v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element + iq*nDOF_trial_element + k]
3544 *
3545 p1_velocity_dofs[eN*nVDOF_element + j];
3546 }/*k*/
3547 }/*id*/
3548 }/*iq*/
3549 }/*eN*/
3550
3551}
3552void getGlobalElementBoundaryBDM1velocityValuesLagrangeRep(int nExteriorElementBoundaries_global,
3553 int nQuadraturePoints_elementBoundary,
3554 int nSpace,
3555 int nDOF_trial_element,
3556 int nVDOF_element,
3557 int *elementBoundaryElementsArray,
3558 int *exteriorElementBoundariesArray,
3559 double * ebqe_v, /*scalar P^1 shape fncts*/
3560 double * p1_velocity_dofs,
3561 double * ebq_global_velocity)
3562{
3563 /***********************************************************************
3564 Assumes local representation for \f$[P^1(E)]^d\f$ is
3565\f[
3566 \vec N_j= \lambda_k \vec e_{id}
3567\f]
3568\f[
3569 k = j / nd, id = j % nd
3570\f]
3571 **********************************************************************/
3572 int ebN,ebNE,eN,iq,id,k,j;
3573
3574 for (ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3575 {
3576 ebN = exteriorElementBoundariesArray[ebNE]; /* global boundary number */
3577 eN = elementBoundaryElementsArray[ebN*2 + 0];
3578
3579 for (iq = 0; iq < nQuadraturePoints_elementBoundary; iq++)
3580 {
3581 for (id = 0; id < nSpace; id++)
3582 {
3583 ebq_global_velocity[ebN*nQuadraturePoints_elementBoundary*nSpace + iq*nSpace + id] = 0.0;
3584 for (k = 0; k < nSpace+1; k++) /* looping over P1 degrees of freedom */
3585 {
3586 j = k*nSpace+ id; /*id*(nSpace+1) + k;*/
3587 ebq_global_velocity[ebN*nQuadraturePoints_elementBoundary*nSpace + iq*nSpace + id] +=
3588 ebqe_v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element + iq*nDOF_trial_element + k]
3589 *
3590 p1_velocity_dofs[eN*nVDOF_element + j];
3591 }/*k*/
3592 }/*id*/
3593 }/*iq*/
3594 }/*eN*/
3595
3596}
3597
3598void getElementBoundaryBDM1velocityValuesLagrangeRep(int nElements_global,
3599 int nBoundaries_Element,
3600 int nQuadraturePoints_elementBoundary,
3601 int nSpace,
3602 int nDOF_trial_element,
3603 int nVDOF_element,
3604 int *elementBoundaryElementsArray, /* not used */
3605 int *exteriorElementBoundariesArray, /*not used */
3606 double * ebq_v, /*scalar P^1 shape fncts*/
3607 double * p1_velocity_dofs,
3608 double * ebq_velocity)
3609{
3610 /***********************************************************************
3611 Assumes local representation for \f$[P^1(E)]^d\f$ is
3612\f[
3613 \vec N_j= \lambda_k \vec e_{id}
3614\f]
3615\f[
3616 k = j / nd, id = j % nd
3617\f]
3618 **********************************************************************/
3619 int ebN,eN,iq,id,k,j;
3620
3621 for (eN = 0; eN < nElements_global; eN++)
3622 {
3623
3624 for (ebN=0; ebN < nBoundaries_Element; ebN ++)
3625 {
3626 for (iq = 0; iq < nQuadraturePoints_elementBoundary; iq++)
3627 {
3628 for (id = 0; id < nSpace; id++)
3629 {
3630 ebq_velocity[eN*nBoundaries_Element*nQuadraturePoints_elementBoundary*nSpace
3631 + ebN*nQuadraturePoints_elementBoundary*nSpace + iq*nSpace + id] = 0.0;
3632 for (k = 0; k < nSpace+1; k++)
3633 {
3634 j = k*nSpace+ id; /*id*(nSpace+1) + k;*/
3635 ebq_velocity[eN*ebN*nQuadraturePoints_elementBoundary*nSpace
3636 + nBoundaries_Element*nQuadraturePoints_elementBoundary*nSpace + iq*nSpace + id] +=
3637 ebq_v[eN*nBoundaries_Element*nQuadraturePoints_elementBoundary*nDOF_trial_element
3638 + ebN*nQuadraturePoints_elementBoundary*iq*nDOF_trial_element + k]
3639 *
3640 p1_velocity_dofs[eN*nVDOF_element + j];
3641 }/*k*/
3642 }/*id*/
3643 }/*iq*/
3644 }/*ebN*/
3645 }/*eN*/
3646
3647}
3648
3649/***********************************************************************
3650 try implementing Sun-Wheeler Gauss-Seidel velocity postprocessing
3651 ***********************************************************************/
3652void calculateConservationResidualGlobalBoundaries(int nElements_global,
3653 int nInteriorElementBoundaries_global,
3654 int nExteriorElementBoundaries_global,
3655 int nElementBoundaries_element,
3656 int nQuadraturePoints_elementBoundary,
3657 int nNodes_element,
3658 int nSpace,
3659 int* interiorElementBoundaries,
3660 int* exteriorElementBoundaries,
3661 int* elementBoundaryElements,
3662 int* elementBoundaryLocalElementBoundaries,
3663 int* exteriorElementBoundariesToSkip,
3664 double* dS,
3665 double* normal,
3666 double* elementResidual,
3667 double* velocity,
3668 double* conservationResidual)
3669{
3670 int ebNI,ebNE,ebN,eN,nN,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,k,I;
3671 register double flux,ds;
3672 /*mwf debug*/
3673/* register double signDebug = -1.0; */
3674 /*first loop through and get element residual sums*/
3675 for (eN = 0; eN < nElements_global; eN++)
3676 {
3677 for (nN = 0; nN < nNodes_element; nN++)
3678 {
3679 conservationResidual[eN] += elementResidual[eN*nNodes_element + nN];
3680 }
3681 }
3682 /*now loop through element boundaries and update element sums*/
3683 /*interior*/
3684 for (ebNI = 0; ebNI < nInteriorElementBoundaries_global; ebNI++)
3685 {
3686 ebN = interiorElementBoundaries[ebNI];
3687 left_eN = elementBoundaryElements[ebN*2+0];
3688 right_eN = elementBoundaryElements[ebN*2+1];
3689 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3690 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
3691
3692 flux = 0.0;
3693 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3694 {
3695 ds = dS[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
3696 left_ebN_element*nQuadraturePoints_elementBoundary+
3697 k];
3698 for (I = 0; I < nSpace; I++)
3699 {
3700 flux+= velocity[ebN*nQuadraturePoints_elementBoundary*nSpace+
3701 k*nSpace+
3702 I]
3703 *
3704 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
3705 k*nSpace+
3706 I]
3707 *
3708 ds;
3709 }
3710 }/*k*/
3711 conservationResidual[left_eN] += flux;
3712 conservationResidual[right_eN]-= flux;
3713
3714 }/*ebNI*/
3715 /*exterior*/
3716 for (ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3717 {
3718 if (!exteriorElementBoundariesToSkip[ebNE])
3719 {
3720 ebN = exteriorElementBoundaries[ebNE];
3721 eN = elementBoundaryElements[ebN*2+0];
3722 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3723 flux = 0.0;
3724 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3725 {
3726 ds = dS[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
3727 ebN_element*nQuadraturePoints_elementBoundary+
3728 k];
3729 for (I = 0; I < nSpace; I++)
3730 {
3731 flux+= velocity[ebN*nQuadraturePoints_elementBoundary*nSpace+
3732 k*nSpace+
3733 I]
3734 *
3735 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
3736 k*nSpace+
3737 I]
3738 *
3739 ds;
3740 }
3741 }/*k*/
3742 conservationResidual[eN] += flux;
3743 }/*not skipping*/
3744 }/*ebNE*/
3745}
3746
3747void sunWheelerGSsweep(int nElements_global,
3748 int nElementBoundaries_global,
3749 int nInteriorElementBoundaries_global,
3750 int nExteriorElementBoundaries_global,
3751 int nElementBoundaries_element,
3752 int nQuadraturePoints_elementBoundary,
3753 int nSpace,
3754 int* interiorElementBoundaries,
3755 int* exteriorElementBoundaries,
3756 int* elementBoundaryElements,
3757 int* elementBoundaryLocalElementBoundaries,
3758 double* dS,
3759 double* normal,
3760 double* sqrt_det_g,
3761 double* alpha,
3762 double* fluxCorrection,
3763 double* conservationResidual)
3764{
3765 int ebNI,ebN,left_eN,right_eN,left_ebN_element,right_ebN_element;
3766 register double area,areaFact,F_ebN;
3767 areaFact = 1.0;
3768 if (nSpace == 3)
3769 areaFact = 0.5;
3770 /*interior faces*/
3771 for (ebNI = 0; ebNI < nInteriorElementBoundaries_global; ebNI++)
3772 {
3773 ebN = interiorElementBoundaries[ebNI];
3774 left_eN = elementBoundaryElements[ebN*2+0];
3775 right_eN = elementBoundaryElements[ebN*2+1];
3776 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3777 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
3778 /*assumed affine*/
3779 area = areaFact*sqrt_det_g[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary +
3780 left_ebN_element*nQuadraturePoints_elementBoundary + 0];
3781 F_ebN = alpha[ebN*2+0]*conservationResidual[left_eN]
3782 + alpha[ebN*2+1]*conservationResidual[right_eN];
3783
3784 fluxCorrection[ebN] += F_ebN;
3785
3786 /*Our residual is (R,1)_e-<Fn_f,n_e>_e which is opposite of our paper formulation*/
3787 conservationResidual[left_eN] -= area*F_ebN;
3788 conservationResidual[right_eN]+= area*F_ebN;
3789 }
3790 /*need to figure out what to do about exterior faces*/
3791 /*what about reversing?*/
3792
3793
3794}
3795
3796void fluxCorrectionVelocityUpdate(int nElements_global,
3797 int nElementBoundaries_global,
3798 int nInteriorElementBoundaries_global,
3799 int nExteriorElementBoundaries_global,
3800 int nElementBoundaries_element,
3801 int nQuadraturePoints_elementBoundary,
3802 int nSpace,
3803 int* interiorElementBoundaries,
3804 int* exteriorElementBoundaries,
3805 int* elementBoundaryElements,
3806 int* elementBoundaryLocalElementBoundaries,
3807 double* dS,
3808 double* normal,
3809 double* fluxCorrection,
3810 double* vConservative,
3811 double* vConservative_element)
3812{
3813 int eN,ebN_element,ebNI,ebNE,ebN,left_eN,right_eN,left_ebN_element,right_ebN_element,k,I;
3814
3815 for (ebN = 0; ebN < nElementBoundaries_global; ebN++)
3816 {
3817 for (k=0; k < nQuadraturePoints_elementBoundary; k++)
3818 {
3819 for (I=0; I < nSpace; I++)
3820 {
3821 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace +
3822 k*nSpace + I]
3823 -= fluxCorrection[ebN]
3824 * normal[ebN*nQuadraturePoints_elementBoundary*nSpace +
3825 k*nSpace + I];
3826 }
3827 }/*k*/
3828 }/*ebN*/
3829 /*copy the global velocity onto the elements*/
3830 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
3831 {
3832 ebN = interiorElementBoundaries[ebNI];
3833 left_eN = elementBoundaryElements[ebN*2+0];
3834 right_eN = elementBoundaryElements[ebN*2+1];
3835 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3836 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
3837 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3838 for(I=0;I<nSpace;I++)
3839 {
3840 vConservative_element[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3841 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3842 k*nSpace+
3843 I]
3844 =
3845 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
3846 k*nSpace+
3847 I];
3848 vConservative_element[right_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3849 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3850 k*nSpace+
3851 I]
3852 =
3853 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
3854 k*nSpace+
3855 I];
3856 }/*I,k*/
3857 }/*ebNI*/
3858 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
3859 {
3860 ebN = exteriorElementBoundaries[ebNE];
3861 eN = elementBoundaryElements[ebN*2+0];
3862 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3863 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3864 for(I=0;I<nSpace;I++)
3865 {
3866 vConservative_element[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3867 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3868 k*nSpace+
3869 I]
3870 =
3871 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
3872 k*nSpace+
3873 I];
3874 }
3875 }/*ebNE*/
3876}
3877
3878void computeFluxCorrectionPWC(int nElementBoundaries_global,
3879 int nInteriorElementBoundaries_global,
3880 int nExteriorElementBoundaries_global,
3881 int* interiorElementBoundaries,
3882 int* exteriorElementBoundaries,
3883 int* elementBoundaryElements,
3884 double* pwcW,
3885 double* pwcV,
3886 double* fluxCorrection)
3887{
3888 /*
3889 generate flux correction from element V's
3890 correction on face f = \partial \Omega_l \cap \partial \Omega_{r}
3891 \Delta_f = (V_l - V_r)/|\gamma_f|
3892 */
3893
3894 int ebNI,ebNE,ebN,eN_left,eN_right;
3895 double V_left,V_right,w_left,w_right;
3896 for (ebNI = 0; ebNI < nInteriorElementBoundaries_global; ebNI++)
3897 {
3898 ebN = interiorElementBoundaries[ebNI];
3899 eN_left = elementBoundaryElements[ebN*2 + 0];
3900 eN_right= elementBoundaryElements[ebN*2 + 1];
3901 V_left = pwcV[eN_left];
3902 V_right = pwcV[eN_right];
3903 w_left = pwcW[ebN*2 + 0];
3904 w_right = pwcW[ebN*2 + 1];
3905 fluxCorrection[ebN] = V_left*w_left + V_right*w_right;
3906 }
3907 for (ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3908 {
3909 ebN = exteriorElementBoundaries[ebNE];
3910 eN_left = elementBoundaryElements[ebN*2 + 0];
3911 V_left = pwcV[eN_left];
3912 w_left = pwcW[ebN*2 + 0];
3913 fluxCorrection[ebN] = V_left*w_left;
3914 }
3915}
3916
3917/***********************************************************************
3918 node star solver data type stuff
3919 ***********************************************************************/
3920
3921
3922int nodeStar_init(int nElements_global,
3923 int nNodes_element,
3924 int nNodes_global,
3925 int* nElements_node,
3926 int* nodeStarElementsArray,
3927 int* nodeStarElementNeighborsArray,
3928 int* N_p,
3929 int** subdomain_dim_p,
3930 double *** subdomain_L_p,
3931 double *** subdomain_R_p,
3932 double *** subdomain_U_p,
3933 PROTEUS_LAPACK_INTEGER*** subdomain_pivots_p,
3934 PROTEUS_LAPACK_INTEGER*** subdomain_column_pivots_p)
3935{
3936 int I;
3937
3938 int N;
3939 int* subdomain_dim;
3940 double** subdomain_R;
3941 double** subdomain_U;
3942 double** subdomain_L;
3943 PROTEUS_LAPACK_INTEGER** subdomain_pivots;
3944 PROTEUS_LAPACK_INTEGER** subdomain_column_pivots;
3945
3946 N = nNodes_global;
3947
3948 *N_p = N;
3949 *subdomain_dim_p = (int*) malloc(N*sizeof(int));
3950 *subdomain_pivots_p = (PROTEUS_LAPACK_INTEGER**)malloc(N*sizeof(PROTEUS_LAPACK_INTEGER*));
3951 *subdomain_column_pivots_p = (PROTEUS_LAPACK_INTEGER**)malloc(N*sizeof(PROTEUS_LAPACK_INTEGER*));
3952 *subdomain_R_p = (double**)malloc(N*sizeof(double*));
3953 *subdomain_U_p = (double**)malloc(N*sizeof(double*));
3954 *subdomain_L_p = (double**)malloc(N*sizeof(double*));
3955
3956 if ( (*subdomain_dim_p == NULL) ||
3957 (*subdomain_R_p == NULL) ||
3958 (*subdomain_U_p == NULL) ||
3959 (*subdomain_L_p == NULL) ||
3960 (*subdomain_pivots_p == NULL) ||
3961 (*subdomain_column_pivots_p == NULL))
3962 {
3963 return 1;
3964 }
3965
3966 subdomain_dim = *subdomain_dim_p;
3967 subdomain_pivots = *subdomain_pivots_p;
3968 subdomain_column_pivots = *subdomain_column_pivots_p;
3969 subdomain_R = *subdomain_R_p;
3970 subdomain_U = *subdomain_U_p;
3971 subdomain_L = *subdomain_L_p;
3972
3973 /*setup local node star system sizes*/
3974 for (I = 0; I < N; I++)
3975 {
3976 subdomain_dim[I] = nElements_node[I];
3977 subdomain_pivots[I] = (PROTEUS_LAPACK_INTEGER*) malloc(subdomain_dim[I]*sizeof(PROTEUS_LAPACK_INTEGER));
3978 subdomain_column_pivots[I] = (PROTEUS_LAPACK_INTEGER*) malloc(subdomain_dim[I]*sizeof(PROTEUS_LAPACK_INTEGER));
3979 subdomain_R[I] = (double*) malloc(subdomain_dim[I]*sizeof(double));
3980 subdomain_U[I] = (double*) malloc(subdomain_dim[I]*sizeof(double));
3981 subdomain_L[I] = (double*) malloc(subdomain_dim[I]*subdomain_dim[I]*sizeof(double));
3982
3983 if ((subdomain_pivots[I] == NULL) ||
3984 (subdomain_column_pivots[I] == NULL) ||
3985 (subdomain_R[I] == NULL) ||
3986 (subdomain_U[I] == NULL) ||
3987 (subdomain_L[I] == NULL))
3988 {
3989 return 1;
3990 }
3991
3992 }/*I*/
3993
3994 return 0;
3995}
3996
3997int nodeStar_free(int N,
3998 int * subdomain_dim,
3999 double ** subdomain_L,
4000 double ** subdomain_R,
4001 double ** subdomain_U,
4002 PROTEUS_LAPACK_INTEGER** subdomain_pivots,
4003 PROTEUS_LAPACK_INTEGER** subdomain_column_pivots)
4004{
4005 int I;
4006 free(subdomain_dim);
4007 for (I = 0; I < N; I++)
4008 {
4009 free(subdomain_pivots[I]);
4010 free(subdomain_column_pivots[I]);
4011 free(subdomain_R[I]);
4012 free(subdomain_U[I]);
4013 free(subdomain_L[I]);
4014 }
4015 free(subdomain_pivots);
4016 free(subdomain_column_pivots);
4017 free(subdomain_R);
4018 free(subdomain_U);
4019 free(subdomain_L);
4020
4021 subdomain_pivots = 0;
4022 subdomain_column_pivots = 0;
4023 subdomain_R = 0;
4024 subdomain_U = 0;
4025 subdomain_L = 0;
4026 return 0;
4027}
4028
4029int nodeStar_setU(NodeStarFactorStruct* nodeStarFactor, double val)
4030{
4031 int I,i;
4032 assert(nodeStarFactor);
4033 for (I = 0; I < nodeStarFactor->N; I++)
4034 for(i=0; i < nodeStarFactor->subdomain_dim[I]; i++)
4035 nodeStarFactor->subdomain_U[I][i] = val;
4036 return 0;
4037}
4038
4039int nodeStar_copy(int other_N,
4040 int * other_subdomain_dim,
4041 double ** other_subdomain_L,
4042 double ** other_subdomain_R,
4043 double ** other_subdomain_U,
4044 PROTEUS_LAPACK_INTEGER** other_subdomain_pivots,
4045 PROTEUS_LAPACK_INTEGER** other_subdomain_column_pivots,
4046 int* N_p,
4047 int** subdomain_dim_p,
4048 double *** subdomain_L_p,
4049 double *** subdomain_R_p,
4050 double *** subdomain_U_p,
4051 PROTEUS_LAPACK_INTEGER*** subdomain_pivots_p,
4052 PROTEUS_LAPACK_INTEGER*** subdomain_column_pivots_p)
4053{
4054 /*assumes both sets of structure have been allocated*/
4055 int I,i,N;
4056 int* subdomain_dim;
4057 double** subdomain_R;
4058 double** subdomain_U;
4059 double** subdomain_L;
4060 PROTEUS_LAPACK_INTEGER** subdomain_pivots;
4061 PROTEUS_LAPACK_INTEGER** subdomain_column_pivots;
4062
4063 int failed = 0;
4064 int realloc = 0;
4065 realloc = other_N != *N_p;
4066 if (!realloc)
4067 {
4068 for (I=0; I < other_N; I++)
4069 {
4070 realloc = realloc || (other_subdomain_dim[I] != (*subdomain_dim_p)[I]);
4071 }
4072 }
4073 if (realloc)
4074 {
4075 /*mwf debug*/
4076 printf("nodeStar_copy needs to reaalloc self_N=%d other_N=%d \n",*N_p,other_N);
4077 failed = nodeStar_free(*N_p,
4078 *subdomain_dim_p,
4079 *subdomain_L_p,
4080 *subdomain_R_p,
4081 *subdomain_U_p,
4082 *subdomain_pivots_p,
4083 *subdomain_column_pivots_p);
4084
4085 if (failed)
4086 return 1;
4087 N = other_N;
4088 *N_p = N;
4089 *subdomain_dim_p = (int*)malloc(N*sizeof(int));
4090 *subdomain_pivots_p = (PROTEUS_LAPACK_INTEGER**)malloc(N*sizeof(PROTEUS_LAPACK_INTEGER*));
4091 *subdomain_column_pivots_p = (PROTEUS_LAPACK_INTEGER**)malloc(N*sizeof(PROTEUS_LAPACK_INTEGER*));
4092 *subdomain_R_p = (double**)malloc(N*sizeof(double*));
4093 *subdomain_U_p = (double**)malloc(N*sizeof(double*));
4094 *subdomain_L_p = (double**)malloc(N*sizeof(double*));
4095
4096 if ( (*subdomain_dim_p == NULL) ||
4097 (*subdomain_R_p == NULL) ||
4098 (*subdomain_U_p == NULL) ||
4099 (*subdomain_L_p == NULL) ||
4100 (*subdomain_pivots_p == NULL) ||
4101 (*subdomain_column_pivots_p == NULL) )
4102 {
4103 return 1;
4104 }
4105 subdomain_dim = *subdomain_dim_p;
4106 subdomain_pivots = *subdomain_pivots_p;
4107 subdomain_column_pivots = *subdomain_column_pivots_p;
4108 subdomain_R = *subdomain_R_p;
4109 subdomain_U = *subdomain_U_p;
4110 subdomain_L = *subdomain_L_p;
4111
4112 /*setup local node star system sizes*/
4113 for (I = 0; I < N; I++)
4114 {
4115 subdomain_dim[I] = other_subdomain_dim[I];
4116 subdomain_pivots[I] = (PROTEUS_LAPACK_INTEGER*)malloc(subdomain_dim[I]*sizeof(PROTEUS_LAPACK_INTEGER));
4117 subdomain_column_pivots[I] = (PROTEUS_LAPACK_INTEGER*)malloc(subdomain_dim[I]*sizeof(PROTEUS_LAPACK_INTEGER));
4118 subdomain_R[I] = (double*)malloc(subdomain_dim[I]*sizeof(double));
4119 subdomain_U[I] = (double*)malloc(subdomain_dim[I]*sizeof(double));
4120 subdomain_L[I] = (double*)malloc(subdomain_dim[I]*subdomain_dim[I]*sizeof(double));
4121
4122 if ((subdomain_pivots[I] == NULL) ||
4123 (subdomain_column_pivots[I] == NULL) ||
4124 (subdomain_R[I] == NULL) ||
4125 (subdomain_U[I] == NULL) ||
4126 (subdomain_L[I] == NULL))
4127 {
4128 return 1;
4129 }
4130
4131 }/*I*/
4132 }/*had to realloc*/
4133 N = *N_p;
4134 subdomain_dim = *subdomain_dim_p;
4135 subdomain_pivots = *subdomain_pivots_p;
4136 subdomain_column_pivots = *subdomain_column_pivots_p;
4137 subdomain_R = *subdomain_R_p;
4138 subdomain_U = *subdomain_U_p;
4139 subdomain_L = *subdomain_L_p;
4140
4141 /*now copy*/
4142 /*mwf debug
4143 printf("nodeStar_copy data now N=%d \n",N);
4144 */
4145 for (I = 0; I < N; I++)
4146 {
4147 assert(subdomain_dim[I] == other_subdomain_dim[I]);
4148 for (i = 0; i < subdomain_dim[I]; i++)
4149 {
4150 subdomain_pivots[I][i] = other_subdomain_pivots[I][i];
4151 subdomain_column_pivots[I][i] = other_subdomain_column_pivots[I][i];
4152 subdomain_R[I][i] = other_subdomain_R[I][i];
4153 subdomain_U[I][i] = other_subdomain_U[I][i];
4154 }
4155 for (i = 0; i < subdomain_dim[I]*subdomain_dim[I]; i++)
4156 subdomain_L[I][i] = other_subdomain_L[I][i];
4157
4158 }/*I*/
4159
4160 return 0;
4161}
4162/***********************************************************************
4163 end node star solver data type stuff
4164 ***********************************************************************/
4165
4166/***********************************************************************
4167 try to put in version of Swedish Postprocessing with Neumann boundaries
4168 enforced explicitly and NodeStarFactorStruct data type
4169
4170 This version,
4171 skips flux boundaries in element integrals: calculateConservationResidualPWL,
4172 calculateConservationJacobianPWL
4173
4174 sets row 0 of flux boundary nodes jacobian to Id and rhs to 0:
4175 calculateConservationJacobianPWL, calculateConservationFluxPWL
4176
4177 When using this version,
4178 do not remove boundary flux terms from element residual
4179 must load in boundary fluxes into ebq_global velocity and ebq velocity
4180 after calling
4181
4182 Note,
4183 fluxElementBoundaries holds global exterior element boundary ids for flux bcs
4184 fluxElementBoundaryNodes holds global node numbers for flux bcs
4185 ***********************************************************************/
4186void calculateConservationResidualPWL(int nElements_global,
4187 int nInteriorElementBoundaries_global,
4188 int nExteriorElementBoundaries_global,
4189 int nElementBoundaries_element,
4190 int nQuadraturePoints_elementBoundary,
4191 int nNodes_element,
4192 int nDOF_element,
4193 int nSpace,
4194 int* interiorElementBoundaries,
4195 int* exteriorElementBoundaries,
4196 int* elementBoundaryElements,
4197 int* elementBoundaryLocalElementBoundaries,
4198 int* elementNodes,
4199 int* dofMapl2g,
4200 int* nodeStarElements,
4201 int* nodeStarElementNeighbors,
4202 int* nElements_node,
4203 int* fluxElementBoundaries,
4204 double* elementResidual,
4205 double* vAverage,
4206 double* dX,
4207 double* w,
4208 double* normal,
4209 NodeStarFactorStruct* nodeStarFactor,
4210 double* conservationResidual,
4211 double* vConservative,
4212 double* vConservative_element)
4213{
4214 int ebNI,ebNE,ebN,eN,eN_star,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,left_eN_star,right_eN_star,nN,nN_global,k,I;
4215 register double flux,fluxAverage,fluxCorrection,dx=0.0;
4216 /*mwf now access node star system using NodeStarFactorStruct*/
4217 double ** starR, ** starU;
4218 /*mwf add for debugging
4219 double * fluxSumDebug;
4220 double resSum;
4221 fluxSumDebug = calloc(nElements_global,sizeof(double));
4222 memset(fluxSumDebug,0,sizeof(double)*nElements_global);
4223 */
4224 /*mwf end for debugging*/
4225 memset(conservationResidual,0,sizeof(double)*nElements_global);
4226 /*mwf debug
4227 printf("calcConsResidPWL nQuadraturePoints_elementBoundary=%d \n",
4228 nQuadraturePoints_elementBoundary);
4229 */
4230 assert(nodeStarFactor);
4231 starR = nodeStarFactor->subdomain_R;
4232 starU = nodeStarFactor->subdomain_U;
4233 /*initial residual with element residual*/
4234 for (eN=0;eN<nElements_global;eN++)
4235 {
4236 for (nN=0;nN<nNodes_element;nN++)
4237 {
4238 nN_global = dofMapl2g[eN*nDOF_element+
4239 nN];
4240 eN_star = nodeStarElements[eN*nNodes_element+
4241 nN];
4242 starR[nN_global][eN_star]
4243 =
4244 elementResidual[eN*nNodes_element+
4245 nN];
4246 conservationResidual[eN]
4247 +=
4248 elementResidual[eN*nNodes_element+
4249 nN];
4250 /*mwf debug
4251 printf("calcConsResPWL eN=%d nN=%d starR=%g\n",eN,nN,
4252 starR[nN_global][eN_star]);
4253 // */
4254 }
4255 // /*mwf debug
4256 // printf("calcConsResPWL eN=%d consRes=%g\n",eN,
4257 // conservationResidual[eN]);
4258 // */
4259 }
4260 /*calculate interior element boundary fluxes and update residual*/
4261 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
4262 {
4263 ebN = interiorElementBoundaries[ebNI];
4264 left_eN = elementBoundaryElements[ebN*2+0];
4265 right_eN = elementBoundaryElements[ebN*2+1];
4266 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4267 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
4268 /* printf("ebN = %d, left_eN = %d, right_eN = %d, left_ebN_element = %d, right_ebN_element = %d\n", */
4269 /* ebN, */
4270 /* left_eN, */
4271 /* right_eN, */
4272 /* left_ebN_element, */
4273 /* right_ebN_element); */
4274 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4275 {
4276 fluxAverage = 0.0;
4277 /* get the integration weight */
4278 dx = dX[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4279 left_ebN_element*nQuadraturePoints_elementBoundary+
4280 k];
4281 /*mwf debug
4282 printf("calcConsResPWL ebNI=%d k=%d dx=%g \n",ebNI,k,dx);
4283 */
4284 for (I=0;I<nSpace;I++)
4285 {
4286 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4287 k*nSpace+
4288 I]
4289 =
4290 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
4291 k*nSpace+
4292 I];
4293 /*mwf debug
4294 printf("pwl get flux vAverage int I=%d, val=%g \n",I,
4295 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
4296 k*nSpace+
4297 I]);
4298 */
4299 fluxAverage
4300 +=
4301 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
4302 k*nSpace+
4303 I]
4304 *
4305 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
4306 k*nSpace+
4307 I];
4308 }
4309 for (nN=0;nN<nNodes_element;nN++)
4310 {
4311 // ARB CHANGE
4312 nN_global = dofMapl2g[left_eN*nDOF_element+
4313 nN];
4314 left_eN_star = nodeStarElements[left_eN*nNodes_element+
4315 nN];
4316 /* check if node is opposite element boundary we're computing and ignore 0 contribution */
4317 /* this shouldn't be necessary */
4318 if (nN != left_ebN_element)
4319 {
4320 right_eN_star = nodeStarElementNeighbors[left_eN*nNodes_element*nElementBoundaries_element+
4321 nN*nElementBoundaries_element+
4322 left_ebN_element];
4323 fluxCorrection = (starU[nN_global][left_eN_star]
4324 -
4325 starU[nN_global][right_eN_star])
4326 *
4327 w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4328 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4329 k*nNodes_element+
4330 nN];
4331 for (I=0;I<nSpace;I++)
4332 {
4333 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4334 k*nSpace+
4335 I]
4336 +=
4337 fluxCorrection
4338 *
4339 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
4340 k*nSpace+
4341 I];
4342 }
4343 flux = (fluxAverage*w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4344 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4345 k*nNodes_element+
4346 nN]
4347 + fluxCorrection)*dx;
4348 starR[nN_global][left_eN_star]
4349 += flux;
4350 starR[nN_global][right_eN_star]
4351 -= flux;
4352 conservationResidual[left_eN] += flux;
4353 conservationResidual[right_eN] -= flux;
4354 /*mwf debug
4355 printf("ppwl nN_global=%d left_eN=%d right_eN=%d fluxAvg=%g fluxCorr=%g flux=%g \n",
4356 nN_global,left_eN,right_eN,fluxAverage,fluxCorrection,flux);
4357 */
4358 /*mwf debug
4359 fluxSumDebug[left_eN] += flux;
4360 fluxSumDebug[right_eN]-= flux;
4361 */
4362 }
4363 }
4364 }
4365 }
4366 /*calculate fluxes on exterior element boundaries and update residual*/
4367 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4368 {
4369 ebN = exteriorElementBoundaries[ebNE];
4370 eN = elementBoundaryElements[ebN*2+0];
4371 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4372
4373 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4374 {
4375 /* get the integration weight, I don't think this was here before */
4376 dx = dX[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4377 ebN_element*nQuadraturePoints_elementBoundary+
4378 k];
4379
4380 fluxAverage=0.0;
4381 for (I=0;I<nSpace;I++)
4382 {
4383 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4384 k*nSpace+
4385 I]
4386 =
4387 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
4388 k*nSpace+
4389 I];
4390 /*mwf debug
4391 printf("pwl get flux vAverage ext ebN=%d free=%d I=%d, val=%g \n",ebN,fluxElementBoundaries[ebNE],
4392 I,vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
4393 k*nSpace+
4394 I]);
4395 */
4396 fluxAverage +=
4397 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
4398 k*nSpace+
4399 I]
4400 *
4401 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
4402 k*nSpace+
4403 I];
4404 }
4405 for (nN=0;nN<nNodes_element;nN++)
4406 {
4407 nN_global = dofMapl2g[eN*nDOF_element+
4408 nN];
4409 eN_star = nodeStarElements[eN*nNodes_element+
4410 nN];
4411 /* check if this node lies opposite the element boundary whose contribution we're computing */
4412 /* in that case there is no flux contribution because the test function is zero*/
4413 /* however, we may be able to remove this conditional for speed */
4414 if (nN != ebN_element && !fluxElementBoundaries[ebNE])/*mwf add skip for flux*/
4415 {
4416 fluxCorrection = starU[nN_global][eN_star]
4417 *
4418 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4419 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4420 k*nNodes_element+
4421 nN];
4422 for (I=0;I<nSpace;I++)
4423 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4424 k*nSpace+
4425 I] +=
4426 fluxCorrection
4427 *
4428 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
4429 k*nSpace+
4430 I];
4431 flux = (fluxAverage
4432 *
4433 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4434 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4435 k*nNodes_element+
4436 nN]
4437 +
4438 fluxCorrection)*dx;
4439 starR[nN_global][eN_star]
4440 += flux;
4441 conservationResidual[eN] += flux;
4442 /*mwf debug
4443 fluxSumDebug[eN] += flux;
4444 */
4445 /*mwf debug
4446 printf("pwl get flux corr ext ebN=%d eN=%d fluxBC=%d fluxCorr=%g flux=%g \n",ebN,eN,
4447 fluxElementBoundaries[ebNE],
4448 fluxCorrection,flux);
4449 */
4450
4451
4452 }
4453 }
4454 }
4455 }
4456 /*copy the global velocity onto the elements*/
4457 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
4458 {
4459 ebN = interiorElementBoundaries[ebNI];
4460 left_eN = elementBoundaryElements[ebN*2+0];
4461 right_eN = elementBoundaryElements[ebN*2+1];
4462 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4463 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
4464 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4465 for(I=0;I<nSpace;I++)
4466 {
4467 vConservative_element[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4468 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4469 k*nSpace+
4470 I]
4471 =
4472 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4473 k*nSpace+
4474 I];
4475 vConservative_element[right_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4476 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4477 k*nSpace+
4478 I]
4479 =
4480 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4481 k*nSpace+
4482 I];
4483 /*mwf debug
4484 printf("pwl int copy ebN=%d eN=[%d,%d] ebN_element=[%d,%d] k=%d vl[%d]=%g, vr[%d]=%g \n",
4485 ebN,left_eN,right_eN,left_ebN_element,right_ebN_element,k,I,
4486 vConservative_element[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4487 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4488 k*nSpace+
4489 I],I,
4490 vConservative_element[right_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4491 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4492 k*nSpace+
4493 I]);
4494 */
4495
4496 }
4497 }
4498
4499 /*copy the global velocity onto the elements*/
4500 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4501 {
4502 ebN = exteriorElementBoundaries[ebNE];
4503 eN = elementBoundaryElements[ebN*2+0];
4504 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4505 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4506 for(I=0;I<nSpace;I++)
4507 {
4508 vConservative_element[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4509 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4510 k*nSpace+
4511 I]
4512 =
4513 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4514 k*nSpace+
4515 I];
4516 /*mwf debug
4517 printf("pwl ext copy ebN=%d eN=%d ebN_element=%d k=%d v[%d]=%g \n",
4518 ebN,eN,ebN_element,k,I,
4519 vConservative_element[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4520 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4521 k*nSpace+
4522 I]);
4523 */
4524
4525 }
4526 }
4527 /*mwf debug
4528 for (eN=0; eN < nElements_global; eN++)
4529 {
4530 printf("leaving calcConsResPWL eN=%d consResid=%g \n",eN,conservationResidual[eN]);
4531 }
4532 */
4533 /*mwf debug
4534 for (eN=0; eN < nElements_global; eN++)
4535 {
4536 resSum = 0;
4537 for (ebN = 0; ebN < nDOF_element; ebN++)
4538 resSum += elementResidual[eN*nDOF_element+ebN];
4539 printf("eN=%d fluxSum=%g elemResid=%g consResid=%g \n",eN,fluxSumDebug[eN],resSum,
4540 conservationResidual[eN]);
4541 }
4542 */
4543 /*mwf debug
4544 free(fluxSumDebug);
4545 */
4546}
4547
4548void calculateConservationJacobianPWL(int nDOF_global,
4549 int nNodes_internal,
4550 int nElements_global,
4551 int nInteriorElementBoundaries_global,
4552 int nExteriorElementBoundaries_global,
4553 int nElementBoundaries_element,
4554 int nQuadraturePoints_elementBoundary,
4555 int nNodes_element,
4556 int nDOF_element,
4557 int nSpace,
4558 int* interiorElementBoundaries,
4559 int* exteriorElementBoundaries,
4560 int* elementBoundaryElements,
4561 int* elementBoundaryLocalElementBoundaries,
4562 int* elementNodes,
4563 int* dofMapl2g,
4564 int* dofStarElements,
4565 int* dofStarElementNeighbors,
4566 int* nElements_node,
4567 int* internalNodes,
4568 int* fluxElementBoundaries,
4569 int* fluxBoundaryNodes,
4570 double* w,
4571 double* normal,
4572 NodeStarFactorStruct* nodeStarFactor)
4573
4574{
4575 int eN,ebNI,ebNE,ebN,eN_star,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,left_eN_star,right_eN_star,nN,nN_global,nNI,k;
4576 PROTEUS_LAPACK_INTEGER INFO=0;
4577 register double wflux;
4578 /*mwf add for boundaries*/
4579 int ii,jj;
4580 double ** starJacobian = nodeStarFactor->subdomain_L;
4581 int * subdomain_dim = nodeStarFactor->subdomain_dim;
4582 PROTEUS_LAPACK_INTEGER ** starPivots = nodeStarFactor->subdomain_pivots;
4583 /*for full pivoting, to avoid pathological bow-tie nodes*/
4584 PROTEUS_LAPACK_INTEGER ** starColPivots = nodeStarFactor->subdomain_column_pivots;
4585
4586 /*zero everything for safety*/
4587 assert(nodeStarFactor);
4588 for (nN = 0; nN < nDOF_global; nN++)
4589 {
4590 for(ii=0; ii < subdomain_dim[nN]; ii++)
4591 {
4592 starPivots[nN][ii]=0;
4593 starColPivots[nN][ii]=0;
4594 for (jj=0; jj < subdomain_dim[nN]; jj++)
4595 starJacobian[nN][ii + jj*subdomain_dim[nN]] = 0.0;
4596 }
4597 }
4598 // printf("break1\n");
4599 /*Load Jacobian entries arising from iterior element boundaries*/
4600 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
4601 {
4602 ebN = interiorElementBoundaries[ebNI];
4603 left_eN = elementBoundaryElements[ebN*2+0];
4604 right_eN = elementBoundaryElements[ebN*2+1];
4605 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4606 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
4607 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4608 {
4609 for (nN=0;nN<nNodes_element;nN++)
4610 {
4611 nN_global = dofMapl2g[left_eN*nDOF_element+
4612 nN];
4613 left_eN_star = dofStarElements[left_eN*nNodes_element+
4614 nN];
4615 /* check if node is opposite element boundary we're computing and ignore 0 contribution */
4616 /* this shouldn't be necessary */
4617 if (nN != left_ebN_element)
4618 {
4619 right_eN_star = dofStarElementNeighbors[left_eN*nNodes_element*nElementBoundaries_element+
4620 nN*nElementBoundaries_element+
4621 left_ebN_element];
4622 wflux = w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4623 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4624 k*nNodes_element+
4625 nN];
4626
4627 starJacobian[nN_global][left_eN_star+
4628 left_eN_star*subdomain_dim[nN_global]]
4629 += wflux;
4630 starJacobian[nN_global][left_eN_star+
4631 right_eN_star*subdomain_dim[nN_global]]
4632 -= wflux;
4633 starJacobian[nN_global][right_eN_star+
4634 left_eN_star*subdomain_dim[nN_global]]
4635 -= wflux;
4636 starJacobian[nN_global][right_eN_star+
4637 right_eN_star*subdomain_dim[nN_global]]
4638 += wflux;
4639 }
4640 }
4641 }
4642 }
4643 // printf("break2\n");
4644 /*Load Jacobian entries arising from exterior element boundaries*/
4645 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4646 {
4647 ebN = exteriorElementBoundaries[ebNE];
4648 eN = elementBoundaryElements[ebN*2+0];
4649 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4650 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4651 {
4652 for (nN=0;nN<nNodes_element;nN++)
4653 {
4654 nN_global = dofMapl2g[eN*nDOF_element+
4655 nN];
4656 eN_star = dofStarElements[eN*nNodes_element+
4657 nN];
4658 /* check if this node lies opposite the element boundary whose contribution we're computing */
4659 /* in that case there is no flux contribution because the test function is zero*/
4660 /* however, we may be able to remove this conditional for speed */
4661 /* cek modified to be: also exclude if the flux is specified on this boundary AND nN is a free node */
4662 /* mwf try to change definition of fluxElementBoundaries so that they are excluded if they contain a node
4663 that is Dirichlet but has no "Dirichlet" faces associated with it
4664 if (nN != ebN_element && !(fluxElementBoundaries[ebNE] && fluxBoundaryNodes[nN_global]) )*/
4665 if (nN != ebN_element && !fluxElementBoundaries[ebNE])
4666 {
4667 starJacobian[nN_global][eN_star+
4668 eN_star*nElements_node[nN_global]]
4669 +=
4670 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4671 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4672 k*nNodes_element+
4673 nN];
4674 }
4675 }
4676 }
4677 }
4678 /*set Dirichlet boundary conditions at one element on interior node stars*/
4679 for (nNI=0;nNI<nNodes_internal;nNI++)
4680 {
4681 nN = internalNodes[nNI];
4682 starJacobian[nN][0]=1.0;
4683 for (eN=1;eN<nElements_node[nN];eN++)
4684 starJacobian[nN][eN*subdomain_dim[nN]]
4685 =0.0;
4686 }
4687 /*repeat for Neumann boundary node stars*/
4688 for (nN=0; nN < nDOF_global; nN++)
4689 {
4690 if (fluxBoundaryNodes[nN]==1)
4691 {
4692 starJacobian[nN][0]=1.0;
4693 for (eN=1;eN<nElements_node[nN];eN++)
4694 starJacobian[nN][eN*subdomain_dim[nN]]
4695 =0.0;
4696
4697 /*mwf debug
4698 printf("jacobian fluxNode=%d nElements_node=%d \n",nN,nElements_node[nN]);
4699 */
4700 }
4701 }
4702
4703 /*factor with lapack*/
4704 for (nN=0;nN<nDOF_global;nN++)
4705 {
4706 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER)subdomain_dim[nN]);
4707/* dgetrf_(&nE_n, */
4708/* &nE_n, */
4709/* starJacobian[nN], */
4710/* &nE_n, */
4711/* starPivots[nN], */
4712/* &INFO); */
4713 /*could try using dgetc2 and dgesc2 to solve when have a zero row because of
4714 isolated nodes*/
4715 dgetc2_(&nE_n,
4716 starJacobian[nN],
4717 &nE_n,
4718 starPivots[nN],
4719 starColPivots[nN],
4720 &INFO);
4721
4722 /*mwf debug*/
4723 /* if (INFO > 0) */
4724 /* { */
4725 /* printf("velPP jac dgetrf INFO=%d nN=%d \n",(int)(INFO),nN); */
4726 /* for (ii=0;ii<nE_n;ii++) */
4727 /* { */
4728 /* for(jj=0;jj<nE_n;jj++) */
4729 /* { */
4730
4731 /* printf("%12.5e \t",starJacobian[nN][ii*nE_n + jj]); */
4732 /* } */
4733 /* printf("\n"); */
4734 /* } */
4735 /* } */
4736 /*assert(INFO == 0);*//*need to turn off if use dgetc2*/
4737 }
4738}
4739void calculateConservationFluxPWL(int nNodes_global,
4740 int nNodes_internal,
4741 int* nElements_node,
4742 int* internalNodes,
4743 int* fluxBoundaryNodes,
4744 NodeStarFactorStruct* nodeStarFactor)
4745{
4746 int nN,nNI,eN;
4747 PROTEUS_LAPACK_INTEGER NRHS=1,INFO=0;
4748 double ** starJ = nodeStarFactor->subdomain_L;
4749 double ** starR = nodeStarFactor->subdomain_R;
4750 double ** starU = nodeStarFactor->subdomain_U;
4751 PROTEUS_LAPACK_INTEGER ** starPivots = nodeStarFactor->subdomain_pivots;
4752 /*for full pivoting, to avoid pathological bow-tie nodes*/
4753 PROTEUS_LAPACK_INTEGER ** starColPivots = nodeStarFactor->subdomain_column_pivots;
4754 double scale;
4755 int * subdomain_dim = nodeStarFactor->subdomain_dim;
4756 char TRANS='N';
4757 /*load -R into U*/
4758 for (nN=0;nN<nNodes_global;nN++)
4759 for(eN=0;eN<subdomain_dim[nN];eN++)
4760 starU[nN][eN] = -starR[nN][eN];
4761 /*set Dirichlet boundary conditions on interior node stars*/
4762 for (nNI=0;nNI<nNodes_internal;nNI++)
4763 {
4764 nN = internalNodes[nNI];
4765 starU[nN][0]=0.0;
4766 }
4767 /*repeat for Neumann boundary node stars*/
4768 for (nN=0; nN < nNodes_global; nN++)
4769 {
4770 if (fluxBoundaryNodes[nN] == 1)
4771 {
4772 starU[nN][0]=0.0;
4773 }
4774 }
4775 /*solve with lapack*/
4776 for (nN=0;nN<nNodes_global;nN++)
4777 {
4778 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER)subdomain_dim[nN]);
4779 /*mwf orig, partial pivoting*/
4780/* dgetrs_(&TRANS, */
4781/* &nE_n, */
4782/* &NRHS, */
4783/* starJ[nN], */
4784/* &nE_n, */
4785/* starPivots[nN], */
4786/* starU[nN], */
4787/* &nE_n, */
4788/* &INFO); */
4789 /*could try using dgetc2 and dgesc2 to solve when have a zero row because of
4790 isolated nodes*/
4791
4792 dgesc2_(&nE_n,
4793 starJ[nN],
4794 &nE_n,
4795 starU[nN],
4796 starPivots[nN],
4797 starColPivots[nN],
4798 &scale);
4799 }
4800}
4801
4802void calculateConservationFluxPWL_noNeumannFix(int nNodes_global,
4803 int* nElements_node,
4804 NodeStarFactorStruct* nodeStarFactor)
4805{
4806 int nN,nNI,eN;
4807 PROTEUS_LAPACK_INTEGER NRHS=1,INFO=0;
4808 double ** starJ = nodeStarFactor->subdomain_L;
4809 double ** starR = nodeStarFactor->subdomain_R;
4810 double ** starU = nodeStarFactor->subdomain_U;
4811 PROTEUS_LAPACK_INTEGER ** starPivots = nodeStarFactor->subdomain_pivots;
4812 /*for full pivoting, to avoid pathological bow-tie nodes*/
4813 PROTEUS_LAPACK_INTEGER ** starColPivots = nodeStarFactor->subdomain_column_pivots;
4814 double scale;
4815 int * subdomain_dim = nodeStarFactor->subdomain_dim;
4816 char TRANS='N';
4817 /*load -R into U*/
4818 for (nN=0;nN<nNodes_global;nN++)
4819 for(eN=0;eN<subdomain_dim[nN];eN++)
4820 starU[nN][eN] = -starR[nN][eN];
4821 /*solve with lapack*/
4822 for (nN=0;nN<nNodes_global;nN++)
4823 {
4824 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER)subdomain_dim[nN]);
4825 dgesc2_(&nE_n,
4826 starJ[nN],
4827 &nE_n,
4828 starU[nN],
4829 starPivots[nN],
4830 starColPivots[nN],
4831 &scale);
4832 }
4833}
4834
4835void calculateConservationResidualPWL_opt(int nNodes_owned,
4836 int nElements_global,
4837 int nInteriorElementBoundaries_global,
4838 int nExteriorElementBoundaries_global,
4839 int nElementBoundaries_element,
4840 int nQuadraturePoints_elementBoundary,
4841 int nNodes_element,
4842 int nSpace,
4843 int* interiorElementBoundaries,
4844 int* exteriorElementBoundaries,
4845 int* elementBoundaryElements,
4846 int* elementBoundaryLocalElementBoundaries,
4847 int* elementNodes,
4848 int* nodeStarElements,
4849 int* nodeStarElementNeighbors,
4850 int* nElements_node,
4851 int* fluxElementBoundaries,
4852 double* elementResidual,
4853 double* vAverage,
4854 double* dX,
4855 double* w,
4856 double* normal,
4857 NodeStarFactorStruct* nodeStarFactor,
4858 double* conservationResidual,
4859 double* vConservative,
4860 double* vConservative_element)
4861{
4862 int foundNonzeroR;
4863 int ebNI,ebNE,ebN,eN,eN_star,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,left_eN_star,right_eN_star,nN,nN_global,k,I;
4864 register double flux,fluxAverage,fluxCorrection,dx=0.0;
4865 /*mwf now access node star system using NodeStarFactorStruct*/
4866 double ** starR, ** starU;
4867 /*mwf add for debugging
4868 double * fluxSumDebug;
4869 double resSum;
4870 fluxSumDebug = calloc(nElements_global,sizeof(double));
4871 memset(fluxSumDebug,0,sizeof(double)*nElements_global);
4872 */
4873 /*mwf end for debugging*/
4874 memset(conservationResidual,0,sizeof(double)*nElements_global);
4875 memset(vConservative,0,sizeof(double)*(nInteriorElementBoundaries_global+nExteriorElementBoundaries_global)*nQuadraturePoints_elementBoundary*nSpace);
4876 /*mwf debug
4877 printf("calcConsResidPWL nQuadraturePoints_elementBoundary=%d \n",
4878 nQuadraturePoints_elementBoundary);
4879 */
4880 assert(nodeStarFactor);
4881 starR = nodeStarFactor->subdomain_R;
4882 starU = nodeStarFactor->subdomain_U;
4883
4884 for (nN=nNodes_owned;nN<nodeStarFactor->N;nN++)
4885 for(eN=0; eN < nodeStarFactor->subdomain_dim[nN]; eN++)
4886 {
4887 starU[nN][eN]=0.0;
4888 }
4889 /*initial residual with element residual*/
4890 for (eN=0;eN<nElements_global;eN++)
4891 {
4892 for (nN=0;nN<nNodes_element;nN++)
4893 {
4894 nN_global = elementNodes[eN*nNodes_element+
4895 nN];
4896 eN_star = nodeStarElements[eN*nNodes_element+
4897 nN];
4898 starR[nN_global][eN_star]
4899 =
4900 elementResidual[eN*nNodes_element+
4901 nN];
4902 conservationResidual[eN]
4903 +=
4904 elementResidual[eN*nNodes_element+
4905 nN];
4906 /*mwf debug
4907 printf("calcConsResPWL eN=%d nN=%d starR=%g\n",eN,nN,
4908 starR[nN_global][eN_star]);
4909 */
4910 }
4911 /*mwf debug
4912 printf("calcConsResPWL eN=%d consRes=%g\n",eN,
4913 conservationResidual[eN]);
4914 */
4915 }
4916 /*calculate interior element boundary fluxes and update residual*/
4917 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
4918 {
4919 ebN = interiorElementBoundaries[ebNI];
4920 left_eN = elementBoundaryElements[ebN*2+0];
4921 right_eN = elementBoundaryElements[ebN*2+1];
4922 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4923 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
4924 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4925 {
4926 fluxAverage = 0.0;
4927 /* get the integration weight */
4928 dx = dX[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4929 left_ebN_element*nQuadraturePoints_elementBoundary+
4930 k];
4931 /*mwf debug
4932 printf("calcConsResPWL ebNI=%d k=%d dx=%g \n",ebNI,k,dx);
4933 */
4934 for (I=0;I<nSpace;I++)
4935 {
4936
4937 fluxAverage
4938 +=
4939 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
4940 k*nSpace+
4941 I]
4942 *
4943 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
4944 k*nSpace+
4945 I];
4946 }
4947 for (nN=0;nN<nNodes_element;nN++)
4948 {
4949 nN_global = elementNodes[left_eN*nNodes_element+
4950 nN];
4951 left_eN_star = nodeStarElements[left_eN*nNodes_element+
4952 nN];
4953 /* check if node is opposite element boundary we're computing and ignore 0 contribution */
4954 /* this shouldn't be necessary */
4955 if (nN != left_ebN_element)
4956 {
4957 right_eN_star = nodeStarElementNeighbors[left_eN*nNodes_element*nElementBoundaries_element+
4958 nN*nElementBoundaries_element+
4959 left_ebN_element];
4960 fluxCorrection = (starU[nN_global][left_eN_star]
4961 -
4962 starU[nN_global][right_eN_star])
4963 *
4964 w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4965 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4966 k*nNodes_element+
4967 nN];
4968 for (I=0;I<nSpace;I++)
4969 {
4970 if (nN_global < nNodes_owned)
4971 {
4972 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
4973 k*nSpace+
4974 I]
4975 +=
4976 fluxCorrection
4977 *
4978 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
4979 k*nSpace+
4980 I];
4981 }
4982 }
4983 flux = (fluxAverage*w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
4984 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
4985 k*nNodes_element+
4986 nN]
4987 + fluxCorrection)*dx;
4988 starR[nN_global][left_eN_star]
4989 += flux;
4990 starR[nN_global][right_eN_star]
4991 -= flux;
4992 conservationResidual[left_eN] += flux;
4993 conservationResidual[right_eN] -= flux;
4994 /*mwf debug
4995 printf("ppwl nN_global=%d left_eN=%d right_eN=%d fluxAvg=%g fluxCorr=%g flux=%g \n",
4996 nN_global,left_eN,right_eN,fluxAverage,fluxCorrection,flux);
4997 */
4998 /*mwf debug
4999 fluxSumDebug[left_eN] += flux;
5000 fluxSumDebug[right_eN]-= flux;
5001 */
5002 }
5003 }
5004 }
5005 }
5006 /*calculate fluxes on exterior element boundaries and update residual*/
5007 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5008 {
5009 ebN = exteriorElementBoundaries[ebNE];
5010 eN = elementBoundaryElements[ebN*2+0];
5011 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5012
5013 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5014 {
5015 /* get the integration weight, I don't think this was here before */
5016 dx = dX[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5017 ebN_element*nQuadraturePoints_elementBoundary+
5018 k];
5019
5020 fluxAverage=0.0;
5021 for (I=0;I<nSpace;I++)
5022 {
5023 fluxAverage +=
5024 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
5025 k*nSpace+
5026 I]
5027 *
5028 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5029 k*nSpace+
5030 I];
5031 }
5032 for (nN=0;nN<nNodes_element;nN++)
5033 {
5034 nN_global = elementNodes[eN*nNodes_element+
5035 nN];
5036 eN_star = nodeStarElements[eN*nNodes_element+
5037 nN];
5038 /* check if this node lies opposite the element boundary whose contribution we're computing */
5039 /* in that case there is no flux contribution because the test function is zero*/
5040 /* however, we may be able to remove this conditional for speed */
5041 if (nN != ebN_element && !fluxElementBoundaries[ebNE])/*mwf add skip for flux*/
5042 {
5043 fluxCorrection = starU[nN_global][eN_star]
5044 *
5045 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5046 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5047 k*nNodes_element+
5048 nN];
5049 if (nN_global < nNodes_owned)
5050 {
5051 for (I=0;I<nSpace;I++)
5052 {
5053 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5054 k*nSpace+
5055 I] +=
5056 fluxCorrection
5057 *
5058 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5059 k*nSpace+
5060 I];
5061 }
5062 }
5063 flux = (fluxAverage
5064 *
5065 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5066 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5067 k*nNodes_element+
5068 nN]
5069 +
5070 fluxCorrection)*dx;
5071 starR[nN_global][eN_star]
5072 += flux;
5073 conservationResidual[eN] += flux;
5074 }
5075 }
5076 }
5077 }
5078}
5079
5080void calculateConservationResidualPWL_primative(int nElements_global,
5081 int nInteriorElementBoundaries_global,
5082 int nExteriorElementBoundaries_global,
5083 int nElementBoundaries_element,
5084 int nQuadraturePoints_elementBoundary,
5085 int nNodes_element,
5086 int nSpace,
5087 int* interiorElementBoundaries,
5088 int* exteriorElementBoundaries,
5089 int* elementBoundaryElements,
5090 int* elementBoundaryLocalElementBoundaries,
5091 int* skipflag_elementBoundaries,
5092 double* elementResidual,
5093 double* dX,
5094 double* normal,
5095 double* conservationResidual,
5096 double* vConservative)
5097{
5098 int ebNI,ebNE,ebN,eN,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,nN,k,I;
5099 register double flux,divergence=0.0;
5100 memset(conservationResidual,0,sizeof(double)*nElements_global);
5101 /* /\*initialize residual with element residual and assume external flux terms are in the element residual*\/ */
5102 /* for (eN=0;eN<nElements_global;eN++) */
5103 /* { */
5104 /* for (nN=0;nN<nNodes_element;nN++) */
5105 /* { */
5106 /* conservationResidual[eN] */
5107 /* += */
5108 /* elementResidual[eN*nNodes_element+ */
5109 /* nN]; */
5110 /* } */
5111 /* } */
5112 /*calculate interior element boundary fluxes and update residual*/
5113 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
5114 {
5115 ebN = interiorElementBoundaries[ebNI];
5116 left_eN = elementBoundaryElements[ebN*2+0];
5117 right_eN = elementBoundaryElements[ebN*2+1];
5118 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5119 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
5120 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5121 {
5122 flux = 0.0;
5123 for (I=0;I<nSpace;I++)
5124 {
5125 flux
5126 +=
5127 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5128 k*nSpace+
5129 I]
5130 *
5131 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5132 k*nSpace+
5133 I];
5134 }
5135 flux*=dX[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5136 left_ebN_element*nQuadraturePoints_elementBoundary+
5137 k];
5138 conservationResidual[left_eN] += flux;
5139 conservationResidual[right_eN] -= flux;
5140 }
5141 }
5142 /*calculate fluxes on exterior element boundaries and update residual*/
5143 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5144 {
5145 ebN = exteriorElementBoundaries[ebNE];
5146 eN = elementBoundaryElements[ebN*2+0];
5147 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5148 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5149 {
5150 flux=0.0;
5151 for (I=0;I<nSpace;I++)
5152 {
5153 flux +=
5154 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5155 k*nSpace+
5156 I]
5157 *
5158 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5159 k*nSpace+
5160 I];
5161 }
5162 flux*=dX[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5163 ebN_element*nQuadraturePoints_elementBoundary+
5164 k];
5165 conservationResidual[eN] += flux;
5166 divergence += flux;
5167 }
5168 }
5169}
5170
5171void calculateConservationJacobianPWL_opt(int nNodes_owned,
5172 int nNodes_global,
5173 int nNodes_internal,
5174 int nElements_global,
5175 int nInteriorElementBoundaries_global,
5176 int nExteriorElementBoundaries_global,
5177 int nElementBoundaries_element,
5178 int nQuadraturePoints_elementBoundary,
5179 int nNodes_element,
5180 int nSpace,
5181 int* interiorElementBoundaries,
5182 int* exteriorElementBoundaries,
5183 int* elementBoundaryElements,
5184 int* elementBoundaryLocalElementBoundaries,
5185 int* elementNodes,
5186 int* nodeStarElements,
5187 int* nodeStarElementNeighbors,
5188 int* nElements_node,
5189 int* internalNodes,
5190 int* fluxElementBoundaries,
5191 int* fluxBoundaryNodes,
5192 double* w,
5193 double* normal,
5194 NodeStarFactorStruct* nodeStarFactor)
5195
5196{
5197 int eN,ebNI,ebNE,ebN,eN_star,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,left_eN_star,right_eN_star,nN,nN_global,nNI,k;
5198 PROTEUS_LAPACK_INTEGER INFO=0;
5199 register double wflux;
5200 /*mwf add for boundaries*/
5201 int ii,jj;
5202 double ** starJacobian = nodeStarFactor->subdomain_L;
5203 int * subdomain_dim = nodeStarFactor->subdomain_dim;
5204 PROTEUS_LAPACK_INTEGER ** starPivots = nodeStarFactor->subdomain_pivots;
5205 /*for full pivoting, to avoid pathological bow-tie nodes*/
5206 PROTEUS_LAPACK_INTEGER ** starColPivots = nodeStarFactor->subdomain_column_pivots;
5207
5208 /*zero everything for safety*/
5209 assert(nodeStarFactor);
5210 for (nN = 0; nN < nNodes_global; nN++)
5211 {
5212 for(ii=0; ii < subdomain_dim[nN]; ii++)
5213 {
5214 starPivots[nN][ii]=0;
5215 starColPivots[nN][ii]=0;
5216 for (jj=0; jj < subdomain_dim[nN]; jj++)
5217 starJacobian[nN][ii + jj*subdomain_dim[nN]] = 0.0;
5218 }
5219 }
5220 /*Load Jacobian entries arising from iterior element boundaries*/
5221 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
5222 {
5223 ebN = interiorElementBoundaries[ebNI];
5224 left_eN = elementBoundaryElements[ebN*2+0];
5225 right_eN = elementBoundaryElements[ebN*2+1];
5226 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5227 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
5228 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5229 {
5230 for (nN=0;nN<nNodes_element;nN++)
5231 {
5232 nN_global = elementNodes[left_eN*nNodes_element+
5233 nN];
5234 left_eN_star = nodeStarElements[left_eN*nNodes_element+
5235 nN];
5236 /* check if node is opposite element boundary we're computing and ignore 0 contribution */
5237 /* this shouldn't be necessary */
5238 if (nN != left_ebN_element && nN_global < nNodes_owned)
5239 {
5240 right_eN_star = nodeStarElementNeighbors[left_eN*nNodes_element*nElementBoundaries_element+
5241 nN*nElementBoundaries_element+
5242 left_ebN_element];
5243 wflux = w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5244 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5245 k*nNodes_element+
5246 nN];
5247
5248 starJacobian[nN_global][left_eN_star+
5249 left_eN_star*subdomain_dim[nN_global]]
5250 += wflux;
5251 starJacobian[nN_global][left_eN_star+
5252 right_eN_star*subdomain_dim[nN_global]]
5253 -= wflux;
5254 starJacobian[nN_global][right_eN_star+
5255 left_eN_star*subdomain_dim[nN_global]]
5256 -= wflux;
5257 starJacobian[nN_global][right_eN_star+
5258 right_eN_star*subdomain_dim[nN_global]]
5259 += wflux;
5260 }
5261 }
5262 }
5263 }
5264 /*Load Jacobian entries arising from exterior element boundaries*/
5265 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5266 {
5267 ebN = exteriorElementBoundaries[ebNE];
5268 eN = elementBoundaryElements[ebN*2+0];
5269 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5270 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5271 {
5272 for (nN=0;nN<nNodes_element;nN++)
5273 {
5274 nN_global = elementNodes[eN*nNodes_element+
5275 nN];
5276 eN_star = nodeStarElements[eN*nNodes_element+
5277 nN];
5278 /* check if this node lies opposite the element boundary whose contribution we're computing */
5279 /* in that case there is no flux contribution because the test function is zero*/
5280 /* however, we may be able to remove this conditional for speed */
5281 /* cek modified to be: also exclude if the flux is specified on this boundary AND nN is a free node */
5282 /* mwf try to change definition of fluxElementBoundaries so that they are excluded if they contain a node
5283 that is Dirichlet but has no "Dirichlet" faces associated with it
5284 if (nN != ebN_element && !(fluxElementBoundaries[ebNE] && fluxBoundaryNodes[nN_global]) )*/
5285 if (nN != ebN_element && !fluxElementBoundaries[ebNE] && nN_global < nNodes_owned)
5286 {
5287 starJacobian[nN_global][eN_star+
5288 eN_star*nElements_node[nN_global]]
5289 +=
5290 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5291 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5292 k*nNodes_element+
5293 nN];
5294 }
5295 }
5296 }
5297 }
5298 /*set Dirichlet boundary conditions at one element on interior node stars*/
5299 for (nNI=0;nNI<nNodes_internal;nNI++)
5300 {
5301 nN = internalNodes[nNI];
5302 if (nN < nNodes_owned)
5303 {
5304 starJacobian[nN][0]=1.0;
5305 for (eN=1;eN<nElements_node[nN];eN++)
5306 starJacobian[nN][eN*subdomain_dim[nN]]
5307 =0.0;
5308 }
5309 }
5310 /*repeat for Neumann boundary node stars*/
5311 for (nN=0; nN < nNodes_owned; nN++)
5312 {
5313 if (fluxBoundaryNodes[nN]==1)
5314 {
5315 starJacobian[nN][0]=1.0;
5316 for (eN=1;eN<nElements_node[nN];eN++)
5317 starJacobian[nN][eN*subdomain_dim[nN]]
5318 =0.0;
5319
5320 /*mwf debug
5321 printf("jacobian fluxNode=%d nElements_node=%d \n",nN,nElements_node[nN]);
5322 */
5323 }
5324 }
5325
5326 /*factor with lapack*/
5327 for (nN=0;nN<nNodes_owned;nN++)
5328 {
5329 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER)subdomain_dim[nN]);
5330/* dgetrf_(&nE_n, */
5331/* &nE_n, */
5332/* starJacobian[nN], */
5333/* &nE_n, */
5334/* starPivots[nN], */
5335/* &INFO); */
5336 /*could try using dgetc2 and dgesc2 to solve when have a zero row because of
5337 isolated nodes*/
5338 dgetc2_(&nE_n,
5339 starJacobian[nN],
5340 &nE_n,
5341 starPivots[nN],
5342 starColPivots[nN],
5343 &INFO);
5344
5345 /*mwf debug*/
5346 if (INFO > 0)
5347 {
5348 printf("velPP jac dgetrf INFO=%d nN=%d \n",(int)(INFO),nN);
5349 for (ii=0;ii<nE_n;ii++)
5350 {
5351 for(jj=0;jj<nE_n;jj++)
5352 {
5353
5354 printf("%12.5e \t",starJacobian[nN][ii*nE_n + jj]);
5355 }
5356 printf("\n");
5357 }
5358 }
5359 /*assert(INFO == 0);*//*need to turn off if use dgetc2*/
5360 }
5361}
5362void calculateConservationFluxPWL_opt(int nNodes_owned,
5363 int nNodes_global,
5364 int nNodes_internal,
5365 int* nElements_node,
5366 int* internalNodes,
5367 int* fluxBoundaryNodes,
5368 NodeStarFactorStruct* nodeStarFactor)
5369{
5370 int nN,nNI,eN;
5371 PROTEUS_LAPACK_INTEGER NRHS=1,INFO=0;
5372 double ** starJ = nodeStarFactor->subdomain_L;
5373 double ** starR = nodeStarFactor->subdomain_R;
5374 double ** starU = nodeStarFactor->subdomain_U;
5375 PROTEUS_LAPACK_INTEGER ** starPivots = nodeStarFactor->subdomain_pivots;
5376 /*for full pivoting, to avoid pathological bow-tie nodes*/
5377 PROTEUS_LAPACK_INTEGER ** starColPivots = nodeStarFactor->subdomain_column_pivots;
5378 double scale;
5379 int * subdomain_dim = nodeStarFactor->subdomain_dim;
5380 char TRANS='N';
5381 /*load -R into U*/
5382 for (nN=0;nN<nNodes_owned;nN++)
5383 for(eN=0;eN<subdomain_dim[nN];eN++)
5384 starU[nN][eN] = -starR[nN][eN];
5385 /*set Dirichlet boundary conditions on interior node stars*/
5386 for (nNI=0;nNI<nNodes_internal;nNI++)
5387 {
5388 nN = internalNodes[nNI];
5389 if (nN < nNodes_owned)
5390 starU[nN][0]=0.0;
5391 }
5392 /*repeat for Neumann boundary node stars*/
5393 for (nN=0; nN < nNodes_owned; nN++)
5394 {
5395 if (fluxBoundaryNodes[nN] == 1)
5396 {
5397 starU[nN][0]=0.0;
5398 }
5399 }
5400 /*solve with lapack*/
5401 for (nN=0;nN<nNodes_owned;nN++)
5402 {
5403 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER)subdomain_dim[nN]);
5404 /*mwf orig, partial pivoting*/
5405/* dgetrs_(&TRANS, */
5406/* &nE_n, */
5407/* &NRHS, */
5408/* starJ[nN], */
5409/* &nE_n, */
5410/* starPivots[nN], */
5411/* starU[nN], */
5412/* &nE_n, */
5413/* &INFO); */
5414 /*could try using dgetc2 and dgesc2 to solve when have a zero row because of
5415 isolated nodes*/
5416
5417 dgesc2_(&nE_n,
5418 starJ[nN],
5419 &nE_n,
5420 starU[nN],
5421 starPivots[nN],
5422 starColPivots[nN],
5423 &scale);
5424 }
5425}
5426void subdomain_U_copy_global2local(int max_nN_owned,
5427 int nElements_global,
5428 int nNodes_element,
5429 int* elementNodes,
5430 int* nodeStarElements,
5431 NodeStarFactorStruct* nodeStarFactor,
5432 double* subdomain_U)
5433{
5434 int eN,nN,eN_star,nN_global;
5435 double ** starU = nodeStarFactor->subdomain_U;
5436 for (eN=0;eN<nElements_global;eN++)
5437 for (nN=0;nN<nNodes_element;nN++)
5438 {
5439 nN_global = elementNodes[eN*nNodes_element+
5440 nN];
5441 eN_star = nodeStarElements[eN*nNodes_element+
5442 nN];
5443 starU[nN_global][eN_star]
5444 =
5445 subdomain_U[eN*nNodes_element+
5446 nN];
5447 }
5448}
5449
5450void subdomain_U_copy_local2global(int max_nN_owned,
5451 int nElements_global,
5452 int nNodes_element,
5453 int* elementNodes,
5454 int* nodeStarElements,
5455 NodeStarFactorStruct* nodeStarFactor,
5456 double* subdomain_U)
5457{
5458 int eN,nN,eN_star,nN_global;
5459 double ** starU = nodeStarFactor->subdomain_U;
5460 for (eN=0;eN<nElements_global;eN++)
5461 for (nN=0;nN<nNodes_element;nN++)
5462 {
5463
5464 nN_global = elementNodes[eN*nNodes_element+
5465 nN];
5466 eN_star = nodeStarElements[eN*nNodes_element+
5467 nN];
5468 if (nN_global < max_nN_owned)
5469 subdomain_U[eN*nNodes_element+
5470 nN] =
5471 starU[nN_global][eN_star];
5472 else//throw out correction if this node star isn't owned
5473 {
5474 subdomain_U[eN*nNodes_element+
5475 nN] = 0.0;
5476 }
5477 }
5478}
5479
5483void updateSelectedExteriorElementBoundaryFlux(int nExteriorElementBoundaries_global,
5484 int nElementBoundaries_element,
5485 int nQuadraturePoints_elementBoundary,
5486 int nDOF_test_element,
5487 int* exteriorElementBoundaries,
5488 int* elementBoundaryElements,
5489 int* elementBoundaryLocalElementBoundaries,
5490 int* skipflag_elementBoundaries,
5491 double* flux,
5492 double* w_dS,
5493 double* residual)
5494{
5495 int ebNE,ebN,eN_global,ebN_element,i,k;
5496 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5497 {
5498 ebN = exteriorElementBoundaries[ebNE];
5499 eN_global = elementBoundaryElements[ebN*2+0];
5500 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5501 if (skipflag_elementBoundaries[ebNE] == 0)
5502 {
5503 for(i=0;i<nDOF_test_element;i++)
5504 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5505 {
5506 /*mwf debug
5507 if (fabs(flux[ebN*nQuadraturePoints_elementBoundary+
5508 k]) > 0.0)
5509 {
5510 printf("postproc updateSelectedExtBnd ebN=%d k=%d eN = %d flux=%g\n",ebN,k,eN_global,
5511 flux[ebN*nQuadraturePoints_elementBoundary+
5512 k]);
5513 }
5514 mwf end debug*/
5515 residual[eN_global*nDOF_test_element+
5516 i]
5517 +=
5518 flux[ebN*nQuadraturePoints_elementBoundary+
5519 k]
5520 *
5521 w_dS[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_test_element+
5522 ebN_element*nQuadraturePoints_elementBoundary*nDOF_test_element+
5523 k*nDOF_test_element+
5524 i];
5525 }
5526 }
5527 }/*ebNE*/
5528}
5529
5530void postprocessAdvectiveVelocityPointEval(int nPoints,
5531 int nSpace,
5532 double updateCoef,
5533 const double* f,
5534 double * velocity)
5535{
5536 /*advective portion of pointwise evaluation of velocity
5537 v = -\ten{a}_h\grad phi_h + \vec f_h
5538
5539 q_a and or q_f may be null if not defined for system
5540 (only reason why calculateFlowVelocity wouldnt work)
5541 */
5542 int k,I;
5543 for (k = 0; k < nPoints; k++)
5544 {
5545 for (I = 0; I < nSpace; I++)
5546 {
5547 velocity[k*nSpace + I] *= updateCoef;
5548 velocity[k*nSpace + I] += f[k*nSpace + I];
5549 }/*I*/
5550 }/*k*/
5551
5552}
5553void postprocessDiffusiveVelocityPointEval(int nPoints,
5554 int nSpace,
5555 double updateCoef,
5556 const double* a,
5557 const double* grad_phi,
5558 double * velocity)
5559{
5560 /*diffusive velocity part of pointwise evaluation of velocity
5561 v = -\ten{a}_h\grad phi_h + \vec f_h
5562
5563 q_a and or q_f may be null if not defined for system
5564 (only reason why calculateFlowVelocity wouldnt work)
5565 */
5566 int eN,k,I,J;
5567 const int nSpace2 = nSpace*nSpace;
5568 for (k = 0; k < nPoints; k++)
5569 {
5570 for (I = 0; I < nSpace; I++)
5571 {
5572 velocity[k*nSpace + I] *= updateCoef;
5573 for (J=0; J < nSpace; J++)
5574 {
5575 velocity[k*nSpace + I] -=
5576 a[k*nSpace2 + I*nSpace + J]
5577 *
5578 grad_phi[k*nSpace + J];
5579 }/*J*/
5580 }/*I*/
5581 }/*k*/
5582}
5583
5584void postprocessDiffusiveVelocityPointEval_sd(int nPoints,
5585 int nSpace,
5586 double updateCoef,
5587 int* rowptr,
5588 int* colind,
5589 const double* a,
5590 const double* grad_phi,
5591 double * velocity)
5592{
5593 /*diffusive velocity part of pointwise evaluation of velocity
5594 v = -\ten{a}_h\grad phi_h + \vec f_h
5595
5596 q_a and or q_f may be null if not defined for system
5597 (only reason why calculateFlowVelocity wouldnt work)
5598 */
5599 int eN,k,I,m;
5600 const int nnz=rowptr[nSpace];
5601 for (k = 0; k < nPoints; k++)
5602 {
5603 for (I = 0; I < nSpace; I++)
5604 {
5605 velocity[k*nSpace + I] *= updateCoef;
5606 for(m=rowptr[I];m<rowptr[I+1];m++)
5607 {
5608 velocity[k*nSpace + I] -=
5609 a[k*nnz+m]
5610 *
5611 grad_phi[k*nSpace + colind[m]];
5612 }/*J*/
5613 }/*I*/
5614 }/*k*/
5615}
5616
5617void calculateElementResidualPWL(int nElements, int nDOF_element_res, int nDOF_element_resPWL,double* alpha, double* elementResidual, double* elementResidualPWL)
5618{
5619 int eN,i,j;
5620 memset(elementResidualPWL,0,sizeof(double)*nElements*nDOF_element_resPWL);
5621 for(eN=0;eN<nElements;eN++)
5622 for(i=0;i<nDOF_element_resPWL;i++)
5623 for(j=0;j<nDOF_element_res;j++)
5624 elementResidualPWL[eN*nDOF_element_resPWL+i] += alpha[i*nDOF_element_res+j]*elementResidual[eN*nDOF_element_res+j];
5625}
5626
5627
5628/***********************************************************************
5629 version of Swedish Postprocessing with internal Neumann boundaries
5630 enforced explicitly
5631
5632 This version,
5633 skips flux boundaries in element integrals: calculateConservationResidualPWL,
5634 calculateConservationJacobianPWL
5635
5636
5637 also will use this with routines that do NOT modify jacobians for
5638 node-stars since it uses dgetc2 Lapack routines that have column
5639 pivoting and can handle the simple singularity
5640
5641 When using this version,
5642 do not remove boundary flux terms from element residual
5643 must load in boundary fluxes into ebq_global velocity and ebq velocity
5644 after calling
5645
5646 Note,
5647 fluxElementBoundaries holds global element boundary ids for flux bcs
5648
5649 ***********************************************************************/
5650void calculateConservationResidualPWL_interiorBoundaries(int nElements_global,
5651 int nInteriorElementBoundaries_global,
5652 int nExteriorElementBoundaries_global,
5653 int nElementBoundaries_element,
5654 int nQuadraturePoints_elementBoundary,
5655 int nNodes_element,
5656 int nSpace,
5657 int* interiorElementBoundaries,
5658 int* exteriorElementBoundaries,
5659 int* elementBoundaryElements,
5660 int* elementBoundaryLocalElementBoundaries,
5661 int* elementNodes,
5662 int* nodeStarElements,
5663 int* nodeStarElementNeighbors,
5664 int* nElements_node,
5665 int* fluxElementBoundaries,
5666 double* elementResidual,
5667 double* vAverage,
5668 double* dX,
5669 double* w,
5670 double* normal,
5671 NodeStarFactorStruct* nodeStarFactor,
5672 double* conservationResidual,
5673 double* vConservative,
5674 double* vConservative_element)
5675{
5676 int ebNI,ebNE,ebN,eN,eN_star,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,left_eN_star,right_eN_star,nN,nN_global,k,I;
5677 register double flux,fluxAverage,fluxCorrection,dx=0.0;
5678 /*mwf now access node star system using NodeStarFactorStruct*/
5679 double ** starR, ** starU;
5680 /*mwf add for debugging
5681 double * fluxSumDebug;
5682 double resSum;
5683 fluxSumDebug = calloc(nElements_global,sizeof(double));
5684 memset(fluxSumDebug,0,sizeof(double)*nElements_global);
5685 */
5686 /*mwf end for debugging*/
5687 memset(conservationResidual,0,sizeof(double)*nElements_global);
5688 /*mwf debug
5689 printf("calcConsResidPWL nQuadraturePoints_elementBoundary=%d \n",
5690 nQuadraturePoints_elementBoundary);
5691 */
5692 assert(nodeStarFactor);
5693 starR = nodeStarFactor->subdomain_R;
5694 starU = nodeStarFactor->subdomain_U;
5695
5696 /*initial residual with element residual*/
5697 for (eN=0;eN<nElements_global;eN++)
5698 {
5699 for (nN=0;nN<nNodes_element;nN++)
5700 {
5701 nN_global = elementNodes[eN*nNodes_element+
5702 nN];
5703 eN_star = nodeStarElements[eN*nNodes_element+
5704 nN];
5705 starR[nN_global][eN_star]
5706 =
5707 elementResidual[eN*nNodes_element+
5708 nN];
5709 conservationResidual[eN]
5710 +=
5711 elementResidual[eN*nNodes_element+
5712 nN];
5713 /*mwf debug
5714 printf("calcConsResPWL eN=%d nN=%d starR=%g\n",eN,nN,
5715 starR[nN_global][eN_star]);
5716 */
5717 }
5718 /*mwf debug
5719 printf("calcConsResPWL eN=%d consRes=%g\n",eN,
5720 conservationResidual[eN]);
5721 */
5722 }
5723 /*calculate interior element boundary fluxes and update residual*/
5724 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
5725 {
5726 ebN = interiorElementBoundaries[ebNI];
5727 left_eN = elementBoundaryElements[ebN*2+0];
5728 right_eN = elementBoundaryElements[ebN*2+1];
5729 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5730 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
5731 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5732 {
5733 fluxAverage = 0.0;
5734 /* get the integration weight */
5735 dx = dX[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5736 left_ebN_element*nQuadraturePoints_elementBoundary+
5737 k];
5738 /*mwf debug
5739 printf("calcConsResPWL ebNI=%d k=%d dx=%g \n",ebNI,k,dx);
5740 */
5741 for (I=0;I<nSpace;I++)
5742 {
5743 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5744 k*nSpace+
5745 I]
5746 =
5747 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
5748 k*nSpace+
5749 I];
5750 /*mwf debug
5751 printf("pwl get flux vAverage int I=%d, val=%g \n",I,
5752 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
5753 k*nSpace+
5754 I]);
5755 */
5756 fluxAverage
5757 +=
5758 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
5759 k*nSpace+
5760 I]
5761 *
5762 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5763 k*nSpace+
5764 I];
5765 }
5766 for (nN=0;nN<nNodes_element;nN++)
5767 {
5768 nN_global = elementNodes[left_eN*nNodes_element+
5769 nN];
5770 left_eN_star = nodeStarElements[left_eN*nNodes_element+
5771 nN];
5772 /* check if node is opposite element boundary we're computing and ignore 0 contribution */
5773 /* also skip if the face is an internal boundary*/
5774 if (nN != left_ebN_element && !fluxElementBoundaries[ebN])
5775 {
5776 right_eN_star = nodeStarElementNeighbors[left_eN*nNodes_element*nElementBoundaries_element+
5777 nN*nElementBoundaries_element+
5778 left_ebN_element];
5779 fluxCorrection = (starU[nN_global][left_eN_star]
5780 -
5781 starU[nN_global][right_eN_star])
5782 *
5783 w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5784 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5785 k*nNodes_element+
5786 nN];
5787 for (I=0;I<nSpace;I++)
5788 {
5789 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5790 k*nSpace+
5791 I]
5792 +=
5793 fluxCorrection
5794 *
5795 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5796 k*nSpace+
5797 I];
5798 }
5799 flux = (fluxAverage*w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5800 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5801 k*nNodes_element+
5802 nN]
5803 + fluxCorrection)*dx;
5804 starR[nN_global][left_eN_star]
5805 += flux;
5806 starR[nN_global][right_eN_star]
5807 -= flux;
5808 conservationResidual[left_eN] += flux;
5809 conservationResidual[right_eN] -= flux;
5810 /*mwf debug
5811 printf("ppwl nN_global=%d left_eN=%d right_eN=%d fluxAvg=%g fluxCorr=%g flux=%g \n",
5812 nN_global,left_eN,right_eN,fluxAverage,fluxCorrection,flux);
5813 */
5814 /*mwf debug
5815 fluxSumDebug[left_eN] += flux;
5816 fluxSumDebug[right_eN]-= flux;
5817 */
5818 }
5819 }
5820 }
5821 }
5822 /*calculate fluxes on exterior element boundaries and update residual*/
5823 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5824 {
5825 ebN = exteriorElementBoundaries[ebNE];
5826 eN = elementBoundaryElements[ebN*2+0];
5827 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5828
5829 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5830 {
5831 /* get the integration weight, I don't think this was here before */
5832 dx = dX[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5833 ebN_element*nQuadraturePoints_elementBoundary+
5834 k];
5835
5836 fluxAverage=0.0;
5837 for (I=0;I<nSpace;I++)
5838 {
5839 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5840 k*nSpace+
5841 I]
5842 =
5843 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
5844 k*nSpace+
5845 I];
5846 /*mwf debug
5847 printf("pwl get flux vAverage ext ebN=%d free=%d I=%d, val=%g \n",ebN,fluxElementBoundaries[ebNE],
5848 I,vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
5849 k*nSpace+
5850 I]);
5851 */
5852 fluxAverage +=
5853 vAverage[ebN*nQuadraturePoints_elementBoundary*nSpace+
5854 k*nSpace+
5855 I]
5856 *
5857 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5858 k*nSpace+
5859 I];
5860 }
5861 for (nN=0;nN<nNodes_element;nN++)
5862 {
5863 nN_global = elementNodes[eN*nNodes_element+
5864 nN];
5865 eN_star = nodeStarElements[eN*nNodes_element+
5866 nN];
5867 /* check if this node lies opposite the element boundary whose contribution we're computing */
5868 /* in that case there is no flux contribution because the test function is zero*/
5869 /* however, we may be able to remove this conditional for speed */
5870 if (nN != ebN_element && !fluxElementBoundaries[ebN])/*skip for flux note ebN not ebNE*/
5871 {
5872 fluxCorrection = starU[nN_global][eN_star]
5873 *
5874 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5875 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5876 k*nNodes_element+
5877 nN];
5878 for (I=0;I<nSpace;I++)
5879 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5880 k*nSpace+
5881 I] +=
5882 fluxCorrection
5883 *
5884 normal[ebN*nQuadraturePoints_elementBoundary*nSpace+
5885 k*nSpace+
5886 I];
5887 flux = (fluxAverage
5888 *
5889 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
5890 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
5891 k*nNodes_element+
5892 nN]
5893 +
5894 fluxCorrection)*dx;
5895 starR[nN_global][eN_star]
5896 += flux;
5897 conservationResidual[eN] += flux;
5898 /*mwf debug
5899 fluxSumDebug[eN] += flux;
5900 */
5901 /*mwf debug
5902 printf("pwl get flux corr ext ebN=%d eN=%d fluxBC=%d fluxCorr=%g flux=%g \n",ebN,eN,
5903 fluxElementBoundaries[ebNE],
5904 fluxCorrection,flux);
5905 */
5906
5907
5908 }
5909 }
5910 }
5911 }
5912 /*copy the global velocity onto the elements*/
5913 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
5914 {
5915 ebN = interiorElementBoundaries[ebNI];
5916 left_eN = elementBoundaryElements[ebN*2+0];
5917 right_eN = elementBoundaryElements[ebN*2+1];
5918 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5919 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
5920 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5921 for(I=0;I<nSpace;I++)
5922 {
5923 vConservative_element[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5924 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5925 k*nSpace+
5926 I]
5927 =
5928 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5929 k*nSpace+
5930 I];
5931 vConservative_element[right_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5932 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5933 k*nSpace+
5934 I]
5935 =
5936 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5937 k*nSpace+
5938 I];
5939 /*mwf debug
5940 printf("pwl int copy ebN=%d eN=[%d,%d] ebN_element=[%d,%d] k=%d vl[%d]=%g, vr[%d]=%g \n",
5941 ebN,left_eN,right_eN,left_ebN_element,right_ebN_element,k,I,
5942 vConservative_element[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5943 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5944 k*nSpace+
5945 I],I,
5946 vConservative_element[right_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5947 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5948 k*nSpace+
5949 I]);
5950 */
5951
5952 }
5953 }
5954 /*copy the global velocity onto the elements*/
5955 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5956 {
5957 ebN = exteriorElementBoundaries[ebNE];
5958 eN = elementBoundaryElements[ebN*2+0];
5959 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5960 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5961 for(I=0;I<nSpace;I++)
5962 {
5963 vConservative_element[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5964 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5965 k*nSpace+
5966 I]
5967 =
5968 vConservative[ebN*nQuadraturePoints_elementBoundary*nSpace+
5969 k*nSpace+
5970 I];
5971 /*mwf debug
5972 printf("pwl ext copy ebN=%d eN=%d ebN_element=%d k=%d v[%d]=%g \n",
5973 ebN,eN,ebN_element,k,I,
5974 vConservative_element[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5975 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5976 k*nSpace+
5977 I]);
5978 */
5979
5980 }
5981 }
5982 /*mwf debug
5983 for (eN=0; eN < nElements_global; eN++)
5984 {
5985 printf("leaving calcConsResPWL eN=%d consResid=%g \n",eN,conservationResidual[eN]);
5986 }
5987 */
5988 /*mwf debug
5989 for (eN=0; eN < nElements_global; eN++)
5990 {
5991 resSum = 0;
5992 for (ebN = 0; ebN < nNodes_element; ebN++)
5993 resSum += elementResidual[eN*nNodes_element+ebN];
5994 printf("eN=%d fluxSum=%g elemResid=%g consResid=%g \n",eN,fluxSumDebug[eN],resSum,
5995 conservationResidual[eN]);
5996 }
5997 */
5998 /*mwf debug
5999 free(fluxSumDebug);
6000 */
6001}
6002/***************************************************
6003 has internal boundaries and does not "fix"
6004 node star systems for pure-neumann degeneracy
6005 **************************************************/
6006void calculateConservationJacobianPWL_interiorBoundaries(int nNodes_global,
6007 int nElements_global,
6008 int nInteriorElementBoundaries_global,
6009 int nExteriorElementBoundaries_global,
6010 int nElementBoundaries_element,
6011 int nQuadraturePoints_elementBoundary,
6012 int nNodes_element,
6013 int nSpace,
6014 int* interiorElementBoundaries,
6015 int* exteriorElementBoundaries,
6016 int* elementBoundaryElements,
6017 int* elementBoundaryLocalElementBoundaries,
6018 int* elementNodes,
6019 int* nodeStarElements,
6020 int* nodeStarElementNeighbors,
6021 int* nElements_node,
6022 int* fluxElementBoundaries,
6023 double* w,
6024 double* normal,
6025 NodeStarFactorStruct* nodeStarFactor)
6026
6027{
6028 int eN,ebNI,ebNE,ebN,eN_star,left_eN,right_eN,ebN_element,left_ebN_element,right_ebN_element,left_eN_star,right_eN_star,nN,nN_global,nNI,k;
6029 PROTEUS_LAPACK_INTEGER INFO=0;
6030 register double wflux;
6031 /*mwf add for boundaries*/
6032 int ii,jj;
6033 double ** starJacobian = nodeStarFactor->subdomain_L;
6034 int * subdomain_dim = nodeStarFactor->subdomain_dim;
6035 PROTEUS_LAPACK_INTEGER ** starPivots = nodeStarFactor->subdomain_pivots;
6036 /*for full pivoting, to avoid pathological bow-tie nodes*/
6037 PROTEUS_LAPACK_INTEGER ** starColPivots = nodeStarFactor->subdomain_column_pivots;
6038
6039 /*zero everything for safety*/
6040 assert(nodeStarFactor);
6041 for (nN = 0; nN < nNodes_global; nN++)
6042 {
6043 for(ii=0; ii < subdomain_dim[nN]; ii++)
6044 {
6045 starPivots[nN][ii]=0;
6046 starColPivots[nN][ii]=0;
6047 for (jj=0; jj < subdomain_dim[nN]; jj++)
6048 starJacobian[nN][ii + jj*subdomain_dim[nN]] = 0.0;
6049 }
6050 }
6051 /*Load Jacobian entries arising from iterior element boundaries*/
6052 for (ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
6053 {
6054 ebN = interiorElementBoundaries[ebNI];
6055 left_eN = elementBoundaryElements[ebN*2+0];
6056 right_eN = elementBoundaryElements[ebN*2+1];
6057 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
6058 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
6059 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
6060 {
6061 for (nN=0;nN<nNodes_element;nN++)
6062 {
6063 nN_global = elementNodes[left_eN*nNodes_element+
6064 nN];
6065 left_eN_star = nodeStarElements[left_eN*nNodes_element+
6066 nN];
6067 /* check if node is opposite element boundary we're computing and ignore 0 contribution */
6068 /* this shouldn't be necessary */
6069 /* skip if flux boundary (now interior or exterior */
6070 if (nN != left_ebN_element && !fluxElementBoundaries[ebN])
6071 {
6072 right_eN_star = nodeStarElementNeighbors[left_eN*nNodes_element*nElementBoundaries_element+
6073 nN*nElementBoundaries_element+
6074 left_ebN_element];
6075 wflux = w[left_eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
6076 left_ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
6077 k*nNodes_element+
6078 nN];
6079
6080 starJacobian[nN_global][left_eN_star+
6081 left_eN_star*subdomain_dim[nN_global]]
6082 += wflux;
6083 starJacobian[nN_global][left_eN_star+
6084 right_eN_star*subdomain_dim[nN_global]]
6085 -= wflux;
6086 starJacobian[nN_global][right_eN_star+
6087 left_eN_star*subdomain_dim[nN_global]]
6088 -= wflux;
6089 starJacobian[nN_global][right_eN_star+
6090 right_eN_star*subdomain_dim[nN_global]]
6091 += wflux;
6092 }
6093 }
6094 }
6095 }
6096 /*Load Jacobian entries arising from exterior element boundaries*/
6097 for (ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
6098 {
6099 ebN = exteriorElementBoundaries[ebNE];
6100 eN = elementBoundaryElements[ebN*2+0];
6101 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
6102 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
6103 {
6104 for (nN=0;nN<nNodes_element;nN++)
6105 {
6106 nN_global = elementNodes[eN*nNodes_element+
6107 nN];
6108 eN_star = nodeStarElements[eN*nNodes_element+
6109 nN];
6110 /* check if this node lies opposite the element boundary whose contribution we're computing */
6111 /* in that case there is no flux contribution because the test function is zero*/
6112 /* however, we may be able to remove this conditional for speed */
6113 /* skip if flux boundary (now interior or exterior */
6114 if (nN != ebN_element && !fluxElementBoundaries[ebN])
6115 {
6116 starJacobian[nN_global][eN_star+
6117 eN_star*nElements_node[nN_global]]
6118 +=
6119 w[eN*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nNodes_element+
6120 ebN_element*nQuadraturePoints_elementBoundary*nNodes_element+
6121 k*nNodes_element+
6122 nN];
6123 }
6124 }
6125 }
6126 }
6127
6128 /*factor with lapack*/
6129 for (nN=0;nN<nNodes_global;nN++)
6130 {
6131 PROTEUS_LAPACK_INTEGER nE_n = ((PROTEUS_LAPACK_INTEGER)subdomain_dim[nN]);
6132 dgetc2_(&nE_n,
6133 starJacobian[nN],
6134 &nE_n,
6135 starPivots[nN],
6136 starColPivots[nN],
6137 &INFO);
6138
6139 /*mwf debug
6140 if (INFO > 0)
6141 {
6142 printf("velPP jac dgetrf INFO=%d nN=%d \n",(int)(INFO),nN);
6143 for (ii=0;ii<nE_n;ii++)
6144 {
6145 for(jj=0;jj<nE_n;jj++)
6146 {
6147
6148 printf("%12.5e \t",starJacobian[nN][ii*nE_n + jj]);
6149 }
6150 printf("\n");
6151 }
6152 }
6153 */
6154
6155 }
6156}
6157
n
Int n
Definition Headers.h:28
r
Double r
Definition Headers.h:83
f
Double f
Definition Headers.h:64
s
Double s
Definition Headers.h:84
u
Double u
Definition Headers.h:89
phi
Double phi
Definition Headers.h:76
calculateConservationResidualPWL
void calculateConservationResidualPWL(int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *elementNodes, int *nodeStarElements, int *nodeStarElementNeighbors, int *nodeStarOffsets, int *nElements_node, double *elementResidual, double *vAverage, double *starU, double *dX, double *w, double *normal, double *conservationResidual, double *starR, double *vConservative, double *vConservative_element)
calculateConservationFluxPWL
void calculateConservationFluxPWL(int nNodes_global, int nNodes_internal, int *nElements_node, int *nodeStarOffsets, int *nodeStarJacobianOffsets, int *internalNodes, double *starR, double *starJ, double *starU)
calculateConservationResidualGlobalBoundaries
void calculateConservationResidualGlobalBoundaries(int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, double *dS, double *normal, double *elementResidual, double *velocity, double *conservationResidual)
calculateConservationJacobianPWL
void calculateConservationJacobianPWL(int nNodes_global, int nNodes_internal, int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *elementNodes, int *nodeStarElements, int *nodeStarElementNeighbors, int *nodeStarOffsets, int *nodeStarJacobianOffsets, int *nElements_node, int *internalNodes, double *w, double *normal, double *starJacobian)
calculateConservationResidualPWL_opt
void calculateConservationResidualPWL_opt(int nNodes_owned, int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *elementNodes, int *nodeStarElements, int *nodeStarElementNeighbors, int *nElements_node, int *fluxElementBoundaries, double *elementResidual, double *vAverage, double *dX, double *w, double *normal, NodeStarFactorStruct *nodeStarFactor, double *conservationResidual, double *vConservative, double *vConservative_element)
Definition postprocessing.c:4835
calculateConservationResidualPWL_primative
void calculateConservationResidualPWL_primative(int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *skipflag_elementBoundaries, double *elementResidual, double *dX, double *normal, double *conservationResidual, double *vConservative)
Definition postprocessing.c:5080
updateSelectedExteriorElementBoundaryFlux
void updateSelectedExteriorElementBoundaryFlux(int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nDOF_test_element, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *skipflag_elementBoundaries, double *flux, double *w_dS, double *residual)
Update the element boundary flux on exterior element boundaries.
Definition postprocessing.c:5483
calculateConservationJacobianPWL_interiorBoundaries
void calculateConservationJacobianPWL_interiorBoundaries(int nNodes_global, int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *elementNodes, int *nodeStarElements, int *nodeStarElementNeighbors, int *nElements_node, int *fluxElementBoundaries, double *w, double *normal, NodeStarFactorStruct *nodeStarFactor)
Definition postprocessing.c:6006
getElementBoundaryBDM1velocityValuesLagrangeRep
void getElementBoundaryBDM1velocityValuesLagrangeRep(int nElements_global, int nBoundaries_Element, int nQuadraturePoints_elementBoundary, int nSpace, int nDOF_trial_element, int nVDOF_element, int *elementBoundaryElementsArray, int *exteriorElementBoundariesArray, double *ebq_v, double *p1_velocity_dofs, double *ebq_velocity)
Definition postprocessing.c:3598
calculateConservationFluxPWL_noNeumannFix
void calculateConservationFluxPWL_noNeumannFix(int nNodes_global, int *nElements_node, NodeStarFactorStruct *nodeStarFactor)
Definition postprocessing.c:4802
getElementBoundaryRT0velocityValues
void getElementBoundaryRT0velocityValues(int nElements_global, int nElementBoundaries_element, int nPoints_elementBoundary, int nSpace, double *x_elementBoundary, double *rt0vdofs_element, double *v_elementBoundary)
Definition postprocessing.c:1544
getGlobalExteriorElementBoundaryBDM1velocityValuesLagrangeRep
void getGlobalExteriorElementBoundaryBDM1velocityValuesLagrangeRep(int nExteriorElementBoundaries_global, int nQuadraturePoints_elementBoundary, int nSpace, int nDOF_trial_element, int nVDOF_element, int *elementBoundaryElementsArray, int *exteriorElementBoundariesArray, double *ebqe_v, double *p1_velocity_dofs, double *ebqe_velocity)
Definition postprocessing.c:3507
updateRT0velocityWithAveragedPotentialP1nc
void updateRT0velocityWithAveragedPotentialP1nc(int nElements_global, int nQuadraturePoints_element, int nSpace, double *detJ, double *quad_a, double *phi, double *gradphi, double *a, double *rt0vdofs)
Definition postprocessing.c:1316
nodeStar_init
int nodeStar_init(int nElements_global, int nNodes_element, int nNodes_global, int *nElements_node, int *nodeStarElementsArray, int *nodeStarElementNeighborsArray, int *N_p, int **subdomain_dim_p, double ***subdomain_L_p, double ***subdomain_R_p, double ***subdomain_U_p, PROTEUS_LAPACK_INTEGER ***subdomain_pivots_p, PROTEUS_LAPACK_INTEGER ***subdomain_column_pivots_p)
Definition postprocessing.c:3922
solveLocalBDM2projection
void solveLocalBDM2projection(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int nDOFs_test_element, int nVDOFs_element, double *BDMprojectionMatFact_element, int *BDMprojectionMatPivots_element, double *w_dS_f, double *ebq_n, double *w_interior_gradients, double *q_velocity, double *ebq_velocity, double *p1_velocity_dofs)
Definition postprocessing.c:3214
postProcessRT0velocityFromP1nc
void postProcessRT0velocityFromP1nc(int nElements_global, int nQuadraturePoints_element, int nDOF_test_element, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int *nFreeDOF_element, int *freeLocal_element, double *detJ, double *sqrt_det_g, double *n, double *elementBarycenters, double *quad_a, double *quad_f, double *w_dV_r, double *w_dV_m, double *u, double *gradu, double *a, double *f, double *r, double *mt, double *rt0vdofs)
Definition postprocessing.c:52
postProcessRT0velocityFromP1ncNoMass_sd
void postProcessRT0velocityFromP1ncNoMass_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_test_element, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int *rowptr, int *colind, int *nFreeDOF_element, int *freeLocal_element, double *detJ, double *sqrt_det_g, double *n, double *elementBarycenters, double *quad_a, double *quad_f, double *w_dV_r, double *u, double *gradu, double *a, double *f, double *r, double *rt0vdofs)
Definition postprocessing.c:1047
nodeStar_copy
int nodeStar_copy(int other_N, int *other_subdomain_dim, double **other_subdomain_L, double **other_subdomain_R, double **other_subdomain_U, PROTEUS_LAPACK_INTEGER **other_subdomain_pivots, PROTEUS_LAPACK_INTEGER **other_subdomain_column_pivots, int *N_p, int **subdomain_dim_p, double ***subdomain_L_p, double ***subdomain_R_p, double ***subdomain_U_p, PROTEUS_LAPACK_INTEGER ***subdomain_pivots_p, PROTEUS_LAPACK_INTEGER ***subdomain_column_pivots_p)
Definition postprocessing.c:4039
solveLocalBDM1projectionFromFlux
void solveLocalBDM1projectionFromFlux(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nDOFs_test_element, int nVDOFs_element, double *BDMprojectionMatFact_element, int *BDMprojectionMatPivots_element, int *elementBoundaryElementsArray, int *elementBoundariesArray, double *w_dS_f, double *ebq_global_flux, double *p1_velocity_dofs)
Definition postprocessing.c:3262
computeFluxCorrectionPWC
void computeFluxCorrectionPWC(int nElementBoundaries_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, double *pwcW, double *pwcV, double *fluxCorrection)
Definition postprocessing.c:3878
fluxCorrectionVelocityUpdate
void fluxCorrectionVelocityUpdate(int nElements_global, int nElementBoundaries_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, double *dS, double *normal, double *fluxCorrection, double *vConservative, double *vConservative_element)
Definition postprocessing.c:3796
postprocessDiffusiveVelocityPointEval
void postprocessDiffusiveVelocityPointEval(int nPoints, int nSpace, double updateCoef, const double *a, const double *grad_phi, double *velocity)
Definition postprocessing.c:5553
getGlobalExteriorElementBoundaryRT0velocityValues
void getGlobalExteriorElementBoundaryRT0velocityValues(int nExteriorElementBoundaries_global, int nPoints_elementBoundary, int nSpace, int *elementBoundaryElementsArray, int *exteriorElementBoundariesArray, double *x_elementBoundary_global, double *rt0vdofs_element, double *v_elementBoundary_global)
Definition postprocessing.c:1637
getGlobalExteriorElementBoundaryRT0velocityValuesFluxRep
void getGlobalExteriorElementBoundaryRT0velocityValuesFluxRep(int nExteriorElementBoundaries_global, int nPoints_elementBoundary_global, int nSpace, int nDetVals_element, double *nodeArray, int *elementNodesArray, int *elementBoundaryElementsArray, int *exteriorElementBoundariesArray, double *abs_det_J, double *x_ebqe, double *rt0vdofs_element, double *v_ebqe)
Definition postprocessing.c:2341
getElementRT0velocityValues
void getElementRT0velocityValues(int nElements_global, int nPoints_element, int nSpace, double *x_element, double *rt0vdofs_element, double *v_element)
Definition postprocessing.c:1505
nodeStar_setU
int nodeStar_setU(NodeStarFactorStruct *nodeStarFactor, double val)
Definition postprocessing.c:4029
calculateElementResidualPWL
void calculateElementResidualPWL(int nElements, int nDOF_element_res, int nDOF_element_resPWL, double *alpha, double *elementResidual, double *elementResidualPWL)
Definition postprocessing.c:5617
postProcessRT0velocityFromP1ncNoMass
void postProcessRT0velocityFromP1ncNoMass(int nElements_global, int nQuadraturePoints_element, int nDOF_test_element, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int *nFreeDOF_element, int *freeLocal_element, double *detJ, double *sqrt_det_g, double *n, double *elementBarycenters, double *quad_a, double *quad_f, double *w_dV_r, double *u, double *gradu, double *a, double *f, double *r, double *rt0vdofs)
Definition postprocessing.c:779
subdomain_U_copy_local2global
void subdomain_U_copy_local2global(int max_nN_owned, int nElements_global, int nNodes_element, int *elementNodes, int *nodeStarElements, NodeStarFactorStruct *nodeStarFactor, double *subdomain_U)
Definition postprocessing.c:5450
postprocessAdvectiveVelocityPointEval
void postprocessAdvectiveVelocityPointEval(int nPoints, int nSpace, double updateCoef, const double *f, double *velocity)
Definition postprocessing.c:5530
getElementBDM2velocityValuesLagrangeRep
void getElementBDM2velocityValuesLagrangeRep(int nElements_global, int nQuadraturePoints_element, int nSpace, int nDOF_trial_element, int nVDOF_element, double *q_v, double *p1_velocity_dofs, double *q_velocity)
Definition postprocessing.c:3428
getElementBoundaryRT0velocityValuesFluxRep
void getElementBoundaryRT0velocityValuesFluxRep(int nElements_global, int nElementBoundaries_element, int nPoints_elementBoundary, int nSpace, int nDetVals_element, double *nodeArray, int *elementNodesArray, double *abs_det_J, double *x_elementBoundary, double *rt0vdofs_element, double *v_elementBoundary)
Definition postprocessing.c:2210
postProcessRT0velocityFromP1nc_sd
void postProcessRT0velocityFromP1nc_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_test_element, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int *rowptr, int *colind, int *nFreeDOF_element, int *freeLocal_element, double *detJ, double *sqrt_det_g, double *n, double *elementBarycenters, double *quad_a, double *quad_f, double *w_dV_r, double *w_dV_m, double *u, double *gradu, double *a, double *f, double *r, double *mt, double *rt0vdofs)
Definition postprocessing.c:415
getElementBDM1velocityValuesLagrangeRep
void getElementBDM1velocityValuesLagrangeRep(int nElements_global, int nQuadraturePoints_element, int nSpace, int nDOF_trial_element, int nVDOF_element, double *q_v, double *p1_velocity_dofs, double *q_velocity)
Definition postprocessing.c:3387
calculateConservationResidualPWL_interiorBoundaries
void calculateConservationResidualPWL_interiorBoundaries(int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *elementNodes, int *nodeStarElements, int *nodeStarElementNeighbors, int *nElements_node, int *fluxElementBoundaries, double *elementResidual, double *vAverage, double *dX, double *w, double *normal, NodeStarFactorStruct *nodeStarFactor, double *conservationResidual, double *vConservative, double *vConservative_element)
Definition postprocessing.c:5650
getGlobalElementBoundaryRT0velocityValuesFluxRep
void getGlobalElementBoundaryRT0velocityValuesFluxRep(int nElementBoundaries_global, int nPoints_elementBoundary_global, int nSpace, int nDetVals_element, double *nodeArray, int *elementNodesArray, int *elementBoundaryElementsArray, double *abs_det_J, double *x_elementBoundary_global, double *rt0vdofs_element, double *v_elementBoundary_global)
Definition postprocessing.c:2278
sunWheelerGSsweep
void sunWheelerGSsweep(int nElements_global, int nElementBoundaries_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, double *dS, double *normal, double *sqrt_det_g, double *alpha, double *fluxCorrection, double *conservationResidual)
Definition postprocessing.c:3747
buildBDM2rhs
void buildBDM2rhs(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nQuadraturePoints_elementInterior, int nSpace, int nDOFs_test_element, int nVDOFs_element, int nDOFs_trial_interior_element, double *BDMprojectionMatFact_element, int *BDMprojectionMatPivots_element, int *edgeFlags, double *w_dS_f, double *ebq_n, double *w_interior_grads, double *w_interior_divfree, double *ebq_velocity, double *q_velocity, double *p1_velocity_dofs)
Definition postprocessing.c:3055
postProcessRT0potentialFromP1nc_sd
void postProcessRT0potentialFromP1nc_sd(int nElements_global, int nQuadraturePoints_element, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int *rowptr, int *colind, double *uQuadratureWeights_element, double *elementBarycenters, double *aElementQuadratureWeights, double *detJ, double *uQuadratureWeights_elementBoundary, double *x, double *u, double *gradu, double *x_elementBoundary, double *u_elementBoundary, double *n, double *a, double *f, double *r, double *rt0vdofs, double *rt0potential)
Definition postprocessing.c:1843
buildLocalBDM2projectionMatrices
void buildLocalBDM2projectionMatrices(int degree, int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nQuadraturePoints_elementInterior, int nSpace, int nDOFs_test_element, int nDOFs_trial_boundary_element, int nDOFs_trial_interior_element, int nVDOFs_element, int *edgeFlags, double *w_dS_f, double *ebq_n, double *ebq_v, double *BDMprojectionMat_element, double *q_basis_vals, double *w_int_test_grads, double *w_int_div_free, double *piola_trial_fun)
Definition postprocessing.c:2710
subdomain_U_copy_global2local
void subdomain_U_copy_global2local(int max_nN_owned, int nElements_global, int nNodes_element, int *elementNodes, int *nodeStarElements, NodeStarFactorStruct *nodeStarFactor, double *subdomain_U)
Definition postprocessing.c:5426
getRT0velocityValuesFluxRep_arbitraryElementMembership
void getRT0velocityValuesFluxRep_arbitraryElementMembership(int nElements_global, int nElementBoundaries_element, int nPoints, int nSpace, int nDetVals_element, const double *nodeArray, const int *elementNodesArray, const double *abs_det_J, const double *x, const int *element_locations, const double *rt0vdofs_element, double *v_element)
Definition postprocessing.c:2406
updateRT0velocityWithAveragedPotentialP1nc_sd
void updateRT0velocityWithAveragedPotentialP1nc_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, int *rowptr, int *colind, double *detJ, double *quad_a, double *phi, double *gradphi, double *a, double *rt0vdofs)
Definition postprocessing.c:1410
projectElementBoundaryFluxToRT0fluxRep
void projectElementBoundaryFluxToRT0fluxRep(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nDOF_RT0V_element, int *elementBoundaryElementsArray, int *elementBoundariesArray, double *elementBoundaryQuadratureWeights, double *flux_elementBoundary, double *rt0vdofs_element)
Definition postprocessing.c:2085
solveLocalBDM1projection
void solveLocalBDM1projection(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int nDOFs_test_element, int nVDOFs_element, double *BDMprojectionMatFact_element, int *BDMprojectionMatPivots_element, double *w_dS_f, double *ebq_n, double *ebq_velocity, double *p1_velocity_dofs)
Definition postprocessing.c:2965
factorLocalBDM1projectionMatrices
void factorLocalBDM1projectionMatrices(int nElements_global, int nVDOFs_element, double *BDMprojectionMat_element, int *BDMprojectionMatPivots_element)
Definition postprocessing.c:2907
buildLocalBDM1projectionMatrices
void buildLocalBDM1projectionMatrices(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int nDOFs_test_element, int nDOFs_trial_element, int nVDOFs_element, double *w_dS_f, double *ebq_n, double *ebq_v, double *BDMprojectionMat_element)
Definition postprocessing.c:2572
getElementLDGvelocityValuesLagrangeRep
void getElementLDGvelocityValuesLagrangeRep(int nElements_global, int nQuadraturePoints_element, int nSpace, int nDOF_trial_element, int nVDOF_element, double *q_v, double *velocity_dofs, double *q_velocity)
Definition postprocessing.c:3469
getElementRT0velocityValuesFluxRep
void getElementRT0velocityValuesFluxRep(int nElements_global, int nElementBoundaries_element, int nPoints_element, int nSpace, int nDetVals_element, double *nodeArray, int *elementNodesArray, double *abs_det_J, double *x_element, double *rt0vdofs_element, double *v_element)
Definition postprocessing.c:2144
postprocessDiffusiveVelocityPointEval_sd
void postprocessDiffusiveVelocityPointEval_sd(int nPoints, int nSpace, double updateCoef, int *rowptr, int *colind, const double *a, const double *grad_phi, double *velocity)
Definition postprocessing.c:5584
factorLocalBDM2projectionMatrices
void factorLocalBDM2projectionMatrices(int nElements_global, int nVDOFs_element, double *BDMprojectionMat_element, int *BDMprojectionMatPivots_element)
Definition postprocessing.c:2936
projectElementBoundaryVelocityToRT0fluxRep
void projectElementBoundaryVelocityToRT0fluxRep(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, double *elementBoundaryQuadratureWeights, double *n, double *v_elementBoundary, double *rt0vdofs_element)
Definition postprocessing.c:2004
calculateConservationJacobianPWL_opt
void calculateConservationJacobianPWL_opt(int nNodes_owned, int nNodes_global, int nNodes_internal, int nElements_global, int nInteriorElementBoundaries_global, int nExteriorElementBoundaries_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nNodes_element, int nSpace, int *interiorElementBoundaries, int *exteriorElementBoundaries, int *elementBoundaryElements, int *elementBoundaryLocalElementBoundaries, int *elementNodes, int *nodeStarElements, int *nodeStarElementNeighbors, int *nElements_node, int *internalNodes, int *fluxElementBoundaries, int *fluxBoundaryNodes, double *w, double *normal, NodeStarFactorStruct *nodeStarFactor)
Definition postprocessing.c:5171
nodeStar_free
int nodeStar_free(int N, int *subdomain_dim, double **subdomain_L, double **subdomain_R, double **subdomain_U, PROTEUS_LAPACK_INTEGER **subdomain_pivots, PROTEUS_LAPACK_INTEGER **subdomain_column_pivots)
Definition postprocessing.c:3997
getGlobalElementBoundaryRT0velocityValues
void getGlobalElementBoundaryRT0velocityValues(int nElementBoundaries_global, int nPoints_elementBoundary, int nSpace, int *elementBoundaryElementsArray, double *x_elementBoundary_global, double *rt0vdofs_element, double *v_elementBoundary_global)
Definition postprocessing.c:1594
postProcessRT0potentialFromP1nc
void postProcessRT0potentialFromP1nc(int nElements_global, int nQuadraturePoints_element, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, double *uQuadratureWeights_element, double *elementBarycenters, double *aElementQuadratureWeights, double *detJ, double *uQuadratureWeights_elementBoundary, double *x, double *u, double *gradu, double *x_elementBoundary, double *u_elementBoundary, double *n, double *a, double *f, double *r, double *rt0vdofs, double *rt0potential)
Definition postprocessing.c:1683
calculateConservationFluxPWL_opt
void calculateConservationFluxPWL_opt(int nNodes_owned, int nNodes_global, int nNodes_internal, int *nElements_node, int *internalNodes, int *fluxBoundaryNodes, NodeStarFactorStruct *nodeStarFactor)
Definition postprocessing.c:5362
sign
#define sign(x, y)
Definition jf.h:44
w
#define w(x)
Definition jf.h:22
nnz
#define nnz
Definition m_comp_co2.h:19
getGlobalElementBoundaryBDM1velocityValuesLagrangeRep
void getGlobalElementBoundaryBDM1velocityValuesLagrangeRep(int nExteriorElementBoundaries_global, int nQuadraturePoints_elementBoundary, int nSpace, int nDOF_trial_element, int nVDOF_element, int *elementBoundaryElementsArray, int *exteriorElementBoundariesArray, double *ebqe_v, double *p1_velocity_dofs, double *ebq_global_velocity)
Definition postprocessing.c:3552
invertLocal
void invertLocal(int nSpace, double A[3][3], double AI[3][3])
Definition postprocessing.c:13
getElementBDM1velocityValuesLagrangeRep_orig
void getElementBDM1velocityValuesLagrangeRep_orig(int nElements_global, int nQuadraturePoints_element, int nSpace, int nDOF_trial_element, int nVDOF_element, double *q_v, double *p1_velocity_dofs, double *q_velocity)
Definition postprocessing.c:3346
buildLocalBDM1projectionMatrices_orig
void buildLocalBDM1projectionMatrices_orig(int nElements_global, int nElementBoundaries_element, int nQuadraturePoints_elementBoundary, int nSpace, int nDOFs_test_element, int nDOFs_trial_element, int nVDOFs_element, double *w_dS_f, double *ebq_n, double *ebq_v, double *BDMprojectionMat_element)
Definition postprocessing.c:2470
postprocessing.h
Python interface to velocity postprocessing library.
dgetrf_
int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info)
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
dgetrs_
int dgetrs_(char *trans, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info)
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
NodeStarFactorStruct::subdomain_column_pivots
PROTEUS_LAPACK_INTEGER ** subdomain_column_pivots
Definition postprocessing.h:27
NodeStarFactorStruct::subdomain_pivots
PROTEUS_LAPACK_INTEGER ** subdomain_pivots
Definition postprocessing.h:26
Generated on for proteus by doxygen 1.18.0