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