proteus  1.8.1
C/C++/Fortran libraries
transportCoefficients.c
Go to the documentation of this file.
1 #include "transportCoefficients.h"
2 
8 /*define relaxation function according to Jacobsen et al 2012, INJNMF*/
9 double relaxationFunction(double phi, double phiStart, double phiEnd)
10 {
11  double H;
12  double x;
13  double Length;
14 
15  if(phiStart < phiEnd)
16  {
17  Length = phiEnd - phiStart;
18  x = (phi - phiStart)/Length;
19  H = 1 - (exp(pow(x,3.5)) - 1.)/ (exp(1) - 1.);
20  }
21  else
22  {
23  Length = -(phiEnd - phiStart);
24  x = 1 - (phi - phiStart)/Length;
25  H = 1 - (exp(pow(x,3.5)) - 1.)/ (exp(1) - 1.);
26  }
27  return H;
28 
29 
30 
31 }
32 /*#define SCALAR_DIFFUSION*/
33 double smoothedHeaviside(double eps, double phi)
34 {
35  double H;
36  if (phi > eps)
37  H=1.0;
38  else if (phi < -eps)
39  H=0.0;
40  else if (phi==0.0)
41  H=0.5;
42  else
43  H = 0.5*(1.0 + phi/eps + sin(M_PI*phi/eps)/M_PI);
44  return H;
45 }
46 
47 double smoothedHeaviside_integral(double eps, double phi)
48 {
49  double HI;
50  if (phi > eps)
51  {
52  HI= phi - eps + 0.5*(eps + 0.5*eps*eps/eps - eps*cos(M_PI*eps/eps)/(M_PI*M_PI)) - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
53  }
54  else if (phi < -eps)
55  {
56  HI=0.0;
57  }
58  else
59  {
60  HI = 0.5*(phi + 0.5*phi*phi/eps - eps*cos(M_PI*phi/eps)/(M_PI*M_PI)) - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
61  }
62  return HI;
63 }
64 
65 double smoothedDirac(double eps, double phi)
66 {
67  double d;
68  if (phi > eps)
69  d=0.0;
70  else if (phi < -eps)
71  d=0.0;
72  else
73  d = 0.5*(1.0 + cos(M_PI*phi/eps))/eps;
74  return d;
75 }
76 
77 double linearHeaviside(double eps, double phi)
78 {
79  double H;
80  if (phi > eps)
81  H=1.0;
82  else if (phi < -eps)
83  H=0.0;
84  else
85  H = 0.5*((phi+eps)/eps);
86  return H;
87 }
88 
89 double linearDirac(double eps, double phi)
90 {
91  double d;
92  if (phi > eps)
93  d=0.0;
94  else if (phi < -eps)
95  d=0.0;
96  else
97  d = 0.5/eps;
98  return d;
99 }
100 
101 void linearADR_ConstantCoefficientsEvaluate(const int nPoints,
102  const int nSpace,
103  const double M,
104  const double *A,
105  const double *B,
106  const double C,
107  const double t,
108  const double *x,
109  const double *u,
110  double *m,
111  double *dm,
112  double *f,
113  double *df,
114  double *a,
115  double *r,
116  double *dr)
117 {
118  int k,I,J;
119  const int nSpace2=nSpace*nSpace;
120  for (k=0;k<nPoints;k++)
121  {
122  m[k]=M*u[k];
123  dm[k]=M;
124 
125  for (I=0;I<nSpace;I++)
126  {
127  f[k*nSpace+I]=B[I]*u[k];
128  df[k*nSpace+I]=B[I];
129  for (J=0;J<nSpace;J++)
130  {
131  a[k*nSpace2+I*nSpace+J]=A[I*nSpace+J];
132  }
133  }
134 
135  r[k]=C*u[k];
136  dr[k]=C;
137  }
138 }
139 
140 void groundwaterTransportCoefficientsEvaluate(const int nPoints,
141  const int nSpace,
142  const double omega,
143  const double d,
144  const double alpha_L,
145  const double alpha_T,
146  const double *v,
147  const double *u,
148  double *m,
149  double *dm,
150  double *f,
151  double *df,
152  double *a)
153 {
154  int k,I,J;
155  const int nSpace2=nSpace*nSpace;
156  double norm_v;
157 
158  for (k=0;k<nPoints;k++)
159  {
160  m[k]=omega*u[k];
161  dm[k]=omega;
162  norm_v = 0.0;
163  for (I=0;I<nSpace;I++)
164  {
165  f[k*nSpace+I]=v[k*nSpace+I]*u[k];
166  df[k*nSpace+I]=v[k*nSpace+I];
167  norm_v += v[k*nSpace+I]*v[k*nSpace+I];
168  }
169  norm_v = sqrt(norm_v);
170  if (norm_v > 0.0)
171  {
172  for (I=0;I<nSpace;I++)
173  {
174  a[k*nSpace2+I*nSpace+I]=omega*d + alpha_T*norm_v + (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
175  for (J=I+1;J<nSpace;J++)
176  {
177  a[k*nSpace2+I*nSpace+J]=(alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
178  a[k*nSpace2+J*nSpace+I]=a[k*nSpace2+I*nSpace+J];
179  }
180  }
181  }
182  else
183  for (I=0;I<nSpace;I++)
184  a[k*nSpace2+I*nSpace+I]=omega*d;
185  }
186 }
187 void groundwaterBiodegradation01EvaluateFC(const int nPoints,
188  const int nSpace,
189  const double omega,
190  const double d_c,
191  const double d_e,
192  const double alpha_L,
193  const double alpha_T,
194  const double Kox_max,
195  const double Kox_C,
196  const double Kox_E,
197  const double Kox_X,
198  const double Yield,
199  const double k_d,
200  const double *v,
201  const double *c_c,
202  const double *c_e,
203  const double *c_x,
204  double *m_c,
205  double *dm_c,
206  double *m_e,
207  double *dm_e,
208  double *m_x,
209  double *dm_x,
210  double *f_c,
211  double *df_c,
212  double *f_e,
213  double *df_e,
214  double *a_c,
215  double *a_e,
216  double *r_c,
217  double *dr_c_dc,
218  double *dr_c_de,
219  double *dr_c_dx,
220  double *r_e,
221  double *dr_e_dc,
222  double *dr_e_de,
223  double *dr_e_dx,
224  double *r_x,
225  double *dr_x_dc,
226  double *dr_x_de,
227  double *dr_x_dx)
228 {
229  int k,I,J;
230  const int nSpace2=nSpace*nSpace;
231  double norm_v;
232  double C,E,X,denomC,denomE,denomX,rox,drox_dC,drox_dE,drox_dX;
233  for (k=0;k<nPoints;k++)
234  {
235  C = c_c[k]; E = c_e[k]; X = c_x[k];
236  m_c[k]=omega*C;
237  dm_c[k]=omega;
238  m_e[k]=omega*E;
239  dm_e[k]=omega;
240  m_x[k]=omega*X;
241  dm_x[k]=omega;
242 
243  norm_v = 0.0;
244  for (I=0;I<nSpace;I++)
245  {
246  f_c[k*nSpace+I]=v[k*nSpace+I]*C;
247  df_c[k*nSpace+I]=v[k*nSpace+I];
248  f_e[k*nSpace+I]=v[k*nSpace+I]*E;
249  df_e[k*nSpace+I]=v[k*nSpace+I];
250 
251  norm_v += v[k*nSpace+I]*v[k*nSpace+I];
252  }
253  norm_v = sqrt(norm_v);
254  if (norm_v > 0.0)
255  {
256  for (I=0;I<nSpace;I++)
257  {
258  a_c[k*nSpace2+I*nSpace+I]=omega*d_c + alpha_T*norm_v + (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
259  for (J=I+1;J<nSpace;J++)
260  {
261  a_c[k*nSpace2+I*nSpace+J]= (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
262  a_c[k*nSpace2+J*nSpace+I]=a_c[k*nSpace2+I*nSpace+J];
263  }
264  a_e[k*nSpace2+I*nSpace+I]=omega*d_e + alpha_T*norm_v + (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
265  for (J=I+1;J<nSpace;J++)
266  {
267  a_e[k*nSpace2+I*nSpace+J]= (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
268  a_e[k*nSpace2+J*nSpace+I]=a_e[k*nSpace2+I*nSpace+J];
269  }
270  }
271  }
272  else
273  {
274  for (I=0;I<nSpace;I++)
275  {
276  a_c[k*nSpace2+I*nSpace+I]=omega*d_c;
277  a_e[k*nSpace2+I*nSpace+I]=omega*d_e;
278  }
279  }/*dispersion calc*/
280  /*reactions*/
281  denomC = C + Kox_C; denomE = E + Kox_E; denomX = X + Kox_X;
282  rox = Kox_max*X*(C/denomC)*(E/denomE)*(Kox_X/denomX);
283  drox_dC = Kox_max*X*(E/denomE)*(Kox_X/denomX)*(1.0/denomC -C/(denomC*denomC));
284  drox_dE = Kox_max*X*(C/denomC)*(Kox_X/denomX)*(1.0/denomE -E/(denomE*denomE));
285  drox_dX = Kox_max*(C/denomC)*(E/denomE)*(Kox_X/denomX -Kox_X*X/(denomX*denomX));
286 
287  r_c[k] = omega*rox;
288  r_e[k] = 3.0*omega*rox;
289  r_x[k] = -Yield*omega*rox + X*omega*k_d;
290 
291  dr_c_dc[k] = omega*drox_dC; dr_c_de[k] = omega*drox_dE; dr_c_dx[k] = omega*drox_dX;
292 
293  dr_e_dc[k] = 3.0*omega*drox_dC; dr_e_de[k] = 3.0*omega*drox_dE; dr_e_dx[k] = 3.0*omega*drox_dX;
294 
295  dr_x_dc[k] = -Yield*omega*drox_dC;
296  dr_x_de[k] = -Yield*omega*drox_dE;
297  dr_x_dx[k] = -Yield*omega*drox_dX + omega*k_d;
298 /* /\*mwf debug*\/ */
299 /* printf("bio01eval k=%d C=%g E=%g X=%g rox=%g \n",k,C,E,X,rox); */
300  }
301 }
302 
303 void groundwaterBryantDawsonIonExEvaluateFC(const int nPoints,
304  const int nSpace,
305  const double omega,
306  const double d_m,
307  const double d_h,
308  const double alpha_L,
309  const double alpha_T,
310  const double K_m,
311  const double K_h,
312  const double K_w,
313  const double Z_tot,
314  const double *v,
315  const double *c_m,
316  const double *c_h,
317  double *m_m,
318  double *dm_m_m,
319  double *dm_m_h,
320  double *m_h,
321  double *dm_h_m,
322  double *dm_h_h,
323  double *f_m,
324  double *df_m,
325  double *f_h,
326  double *df_h,
327  double *a_m,
328  double *a_h,
329  double *phi_h,
330  double *dphi_h,
331  double *r_m,
332  double *dr_m_dm,
333  double *dr_m_dh,
334  double *r_h,
335  double *dr_h_dm,
336  double *dr_h_dh)
337 
338 {
339  int k,I,J;
340  const int nSpace2=nSpace*nSpace;
341  double norm_v;
342  double C_m,C_h,C_oh,C_a,Z_m,Z_h,dC_a_dC_h,denomZ;
343  const double eps = 1.0e-12;
344  for (k=0;k<nPoints;k++)
345  {
346  /*metal, proton, hydroxyl*/
347  C_m = c_m[k]; C_h= c_h[k]; C_oh = K_w/(C_h+eps);
348  /*acidity*/
349  C_a = C_h - C_oh; dC_a_dC_h = 1.0 + K_w/(C_h*C_h+eps);
350  /*sorbed concentrations*/
351  denomZ = 1.0 + K_m*C_m + K_h*C_h;
352  Z_m = K_m*C_m*Z_tot/denomZ;
353  Z_h = K_h*C_h*Z_tot/denomZ;
354 
355  m_m[k] =omega*(C_m + Z_m);
356  dm_m_m[k]=omega*(1.0 + K_m*Z_tot/denomZ - K_m*C_m*Z_tot/(denomZ*denomZ)*K_m*C_m);
357  dm_m_h[k]=omega*( - K_m*C_m*Z_tot/(denomZ*denomZ)*K_h*C_h);
358  m_h[k] =omega*(C_a + Z_h);
359  dm_h_m[k]=omega*( - K_h*C_h*Z_tot/(denomZ*denomZ)*K_m*C_m);
360  dm_h_h[k]=omega*(dC_a_dC_h + K_h*Z_tot/denomZ - K_h*C_h*Z_tot/(denomZ*denomZ)*K_h*C_h);
361 
362  /*nonlinear potential for C_h*/
363  phi_h[k] = C_a;
364  dphi_h[k] = dC_a_dC_h;
365  norm_v = 0.0;
366  for (I=0;I<nSpace;I++)
367  {
368  f_m[k*nSpace+I]=v[k*nSpace+I]*C_m;
369  df_m[k*nSpace+I]=v[k*nSpace+I];
370  f_h[k*nSpace+I]=v[k*nSpace+I]*C_a;
371  df_h[k*nSpace+I]=v[k*nSpace+I]*dC_a_dC_h;
372 
373  norm_v += v[k*nSpace+I]*v[k*nSpace+I];
374  }
375  norm_v = sqrt(norm_v);
376  if (norm_v > 0.0)
377  {
378  for (I=0;I<nSpace;I++)
379  {
380  a_m[k*nSpace2+I*nSpace+I]=omega*d_m + alpha_T*norm_v + (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
381  for (J=I+1;J<nSpace;J++)
382  {
383  a_m[k*nSpace2+I*nSpace+J]= (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
384  a_m[k*nSpace2+J*nSpace+I]=a_m[k*nSpace2+I*nSpace+J];
385  }
386  a_h[k*nSpace2+I*nSpace+I]=omega*d_h + alpha_T*norm_v + (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
387  for (J=I+1;J<nSpace;J++)
388  {
389  a_h[k*nSpace2+I*nSpace+J]= (alpha_L - alpha_T)*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
390  a_h[k*nSpace2+J*nSpace+I]=a_h[k*nSpace2+I*nSpace+J];
391  }
392  }
393  }
394  else
395  {
396  for (I=0;I<nSpace;I++)
397  {
398  a_m[k*nSpace2+I*nSpace+I]=omega*d_m;
399  a_h[k*nSpace2+I*nSpace+I]=omega*d_h;
400  }
401  }/*dispersion calc*/
402  r_m[k] = 0.0;
403  r_h[k] = 0.0;
404  dr_m_dm[k] = 0.0;
405  dr_m_dh[k] = 0.0;
406  dr_h_dm[k] = 0.0;
407  dr_h_dh[k] = 0.0;
408 
409 /* /\*mwf debug*\/ */
410 /* printf("ionExeval k=%d C_m=%g C_h=%g \n",k,C_m,C_h); */
411  }
412 }
413 void groundwaterTransportCoefficientsEvaluate_hetMat(const int nSimplex,
414  const int nPointsPerSimplex,
415  const int nSpace,
416  const double d,
417  const int* materialTypes,
418  const double *omega_types,
419  const double *alpha_L_types,
420  const double *alpha_T_types,
421  const double *v,
422  const double *u,
423  double *m,
424  double *dm,
425  double *f,
426  double *df,
427  double *a)
428 {
429  int i,j,k,I,J,matID;
430  const int nSpace2=nSpace*nSpace;
431  double norm_v;
432 
433  for (i=0;i<nSimplex;i++)
434  {
435  matID = materialTypes[i];
436  for (j=0; j < nPointsPerSimplex; j++)
437  {
438  k = i*nPointsPerSimplex+j;
439 
440  m[k]=omega_types[matID]*u[k];
441  dm[k]=omega_types[matID];
442  norm_v = 0.0;
443  for (I=0;I<nSpace;I++)
444  {
445  f[k*nSpace+I]=v[k*nSpace+I]*u[k];
446  df[k*nSpace+I]=v[k*nSpace+I];
447  norm_v += v[k*nSpace+I]*v[k*nSpace+I];
448  }
449  norm_v = sqrt(norm_v);
450  if (norm_v > 0.0)
451  {
452  for (I=0;I<nSpace;I++)
453  {
454  a[k*nSpace2+I*nSpace+I]=omega_types[matID]*d + alpha_T_types[matID]*norm_v + (alpha_L_types[matID] - alpha_T_types[matID])*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
455  for (J=I+1;J<nSpace;J++)
456  {
457  a[k*nSpace2+I*nSpace+J]=(alpha_L_types[matID] - alpha_T_types[matID])*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
458  a[k*nSpace2+J*nSpace+I]=a[k*nSpace2+I*nSpace+J];
459  }
460  }
461  }
462  else
463  for (I=0;I<nSpace;I++)
464  a[k*nSpace2+I*nSpace+I]=omega_types[matID]*d;
465  }
466  }
467 }
468 void variablySaturatedGroundwaterTransportCoefficientsEvaluate_hetMat(const int nSimplex,
469  const int nPointsPerSimplex,
470  const int nSpace,
471  const double d,
472  const int* materialTypes,
473  const double *theta, /*phase volume fraction*/
474  const double *alpha_L_types,
475  const double *alpha_T_types,
476  const double *v,/*phase darcy velocity*/
477  const double *u,
478  double *m,
479  double *dm,
480  double *f,
481  double *df,
482  double *a)
483 {
484  int i,j,k,I,J,matID;
485  const int nSpace2=nSpace*nSpace;
486  double norm_v;
487 
488  for (i=0;i<nSimplex;i++)
489  {
490  matID = materialTypes[i];
491  for (j=0; j < nPointsPerSimplex; j++)
492  {
493  k = i*nPointsPerSimplex+j;
494 
495  m[k]=theta[k]*u[k];
496  dm[k]=theta[k];
497  norm_v = 0.0;
498  for (I=0;I<nSpace;I++)
499  {
500  f[k*nSpace+I]=v[k*nSpace+I]*u[k];
501  df[k*nSpace+I]=v[k*nSpace+I];
502  norm_v += v[k*nSpace+I]*v[k*nSpace+I];
503  }
504  norm_v = sqrt(norm_v);
505  if (norm_v > 0.0)
506  {
507  for (I=0;I<nSpace;I++)
508  {
509  a[k*nSpace2+I*nSpace+I]=theta[k]*d + alpha_T_types[matID]*norm_v + (alpha_L_types[matID] - alpha_T_types[matID])*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
510  for (J=I+1;J<nSpace;J++)
511  {
512  a[k*nSpace2+I*nSpace+J]=(alpha_L_types[matID] - alpha_T_types[matID])*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
513  a[k*nSpace2+J*nSpace+I]=a[k*nSpace2+I*nSpace+J];
514  }
515  }
516  }
517  else
518  for (I=0;I<nSpace;I++)
519  a[k*nSpace2+I*nSpace+I]=theta[k]*d;
520  }
521  }
522 }
523 void variablySaturatedGroundwaterEnergyTransportCoefficientsEvaluate_hetMat(const int nSimplex,
524  const int nPointsPerSimplex,
525  const int nSpace,
526  const double rho_w,
527  const double rho_n,
528  const double specificHeat_w,
529  const double specificHeat_n,
530  const int* materialTypes,
531  const double *theta, /*phase volume fraction*/
532  const double *thetaS_types,
533  const double *alpha_L_types,
534  const double *alpha_T_types,
535  const double *rho_s_types,
536  const double *specificHeat_s_types,
537  const double *lambda_sat_types,
538  const double *lambda_dry_types,
539  const double *lambda_aniso_types,
540  const double *v,/*phase darcy velocity*/
541  const double *u,
542  double *m,
543  double *dm,
544  double *f,
545  double *df,
546  double *a)
547 {
548  int i,j,k,I,J,matID;
549  const int nSpace2=nSpace*nSpace;
550  double norm_v,tmp,sw,Ke,lambda;
551 
552  for (i=0;i<nSimplex;i++)
553  {
554  matID = materialTypes[i];
555  for (j=0; j < nPointsPerSimplex; j++)
556  {
557  k = i*nPointsPerSimplex+j;
558  tmp = theta[k]*rho_w*specificHeat_w + (thetaS_types[matID]-theta[k])*rho_n*specificHeat_n +
559  (1.0-thetaS_types[matID])*rho_s_types[matID]*specificHeat_s_types[matID];
560  m[k]=tmp*u[k];
561  dm[k]=tmp;
562  norm_v = 0.0;
563  for (I=0;I<nSpace;I++)
564  {
565  f[k*nSpace+I]=v[k*nSpace+I]*rho_w*specificHeat_w*u[k];
566  df[k*nSpace+I]=v[k*nSpace+I]*rho_w*specificHeat_w;
567  norm_v += v[k*nSpace+I]*v[k*nSpace+I];
568  }
569  norm_v = sqrt(norm_v);
570  sw = theta[k]/thetaS_types[matID];
571  if (norm_v > 0.0)
572  {
573  for (I=0;I<nSpace;I++)
574  {
575  a[k*nSpace2+I*nSpace+I]= rho_w*specificHeat_w*alpha_T_types[matID]*norm_v + rho_w*specificHeat_w*(alpha_L_types[matID] - alpha_T_types[matID])*v[k*nSpace+I]*v[k*nSpace+I]/norm_v;
576  for (J=I+1;J<nSpace;J++)
577  {
578  a[k*nSpace2+I*nSpace+J]=rho_w*specificHeat_w*(alpha_L_types[matID] - alpha_T_types[matID])*v[k*nSpace+I]*v[k*nSpace+J]/norm_v;
579  a[k*nSpace2+J*nSpace+I]=a[k*nSpace2+I*nSpace+J];
580  }
581  }
582  }
583 
584  /*todo check with Stacy for right form*/
585  if (sw > 0.05)
586  Ke = 0.7*log(sw) + 1.0;
587  lambda = (lambda_sat_types[matID]-lambda_dry_types[matID])*Ke + lambda_dry_types[matID];
588  for (I=0;I<nSpace;I++)
589  {
590  a[k*nSpace2+I*nSpace+I] += lambda*lambda_aniso_types[matID*nSpace+I];
591  }
592 
593  }
594  }
595 }
596 
597 void nonlinearADR_pqrstEvaluate(const int nPoints,
598  const int nSpace,
599  const double M,
600  const double* A,
601  const double* B,
602  const double C,
603  const double p_pow,
604  const double q_pow,
605  const double r_pow,
606  const double s_pow,
607  const double t_pow,
608  const double t,
609  const double *x,
610  const double *u,
611  double *m,
612  double *dm,
613  double *f,
614  double *df,
615  double *a,
616  double *da,
617  double *phi,
618  double *dphi,
619  double *r,
620  double *dr)
621 {
622  int k,I,J;
623  const int nSpace2=nSpace*nSpace;
624  double uPlus,
625  tmp_f,
626  tmp_df,
627  tmp_a,
628  tmp_da;
629  const double pM1_pow=p_pow-1.0,
630  qM1_pow=q_pow-1.0,
631  rM1_pow=r_pow-1.0,
632  sM1_pow=s_pow-1.0,
633  tM1_pow=t_pow-1.0;
634  for (k=0;k<nPoints;k++)
635  {
636  uPlus = u[k] > 0.0 ? u[k] : 0.0;
637  m[k] = p_pow > 1.0 ? M*pow(uPlus,p_pow) : M*u[k];
638  dm[k] = p_pow > 1.0 ? p_pow*M*pow(uPlus,pM1_pow) : M;
639 
640  tmp_f = q_pow > 1.0 ? pow(uPlus,q_pow) : u[k];
641  tmp_df = q_pow > 1.0 ? q_pow*pow(uPlus,qM1_pow) : 1.0;
642 
643  tmp_a = t_pow > 0.0 ? pow(uPlus,t_pow) : 1.0;
644  tmp_da = t_pow > 0.0 ? t_pow*pow(uPlus,tM1_pow) : 0.0;
645 
646  for (I=0;I<nSpace;I++)
647  {
648  f[k*nSpace+I] = B[I]*tmp_f;
649  df[k*nSpace+I] = B[I]*tmp_df;
650  for (J=0;J<nSpace;J++)
651  {
652  a[k*nSpace2+I*nSpace+J] = A[I*nSpace + J]*tmp_a;
653  da[k*nSpace2+I*nSpace+J] = A[I*nSpace + J]*tmp_da;
654  }
655  }
656  phi[k] = r_pow > 1.0 ? pow(uPlus,r_pow) : u[k];
657  dphi[k] = r_pow > 1.0 ? r_pow*pow(uPlus,rM1_pow) : 1.0;
658 
659  r[k] = s_pow > 1.0 ? C*pow(uPlus,s_pow) : C*u[k];
660  dr[k] = s_pow > 1.0 ? s_pow*C*pow(uPlus,sM1_pow) : C;
661  }
662 }
663 
664 void nonlinearADR_pqrstDualEvaluate(const int nPoints,
665  const int nSpace,
666  const double M,
667  const double* A,
668  const double* B,
669  const double C,
670  const double p1,
671  const double q1,
672  const double r1,
673  const double s1,
674  const double t1,
675  const double p2,
676  const double q2,
677  const double r2,
678  const double s2,
679  const double t2,
680  const double t,
681  const double *x,
682  const double *u,
683  double *m,
684  double *dm,
685  double *f,
686  double *df,
687  double *a,
688  double *da,
689  double *phi,
690  double *dphi,
691  double *r,
692  double *dr)
693 {
694  int k,I,J;
695  const int nSpace2=nSpace*nSpace;
696  double max_1mu_0,atmp,datmp;
697  const double p2M1=p2-1.0,
698  q2M1=q2-1.0,
699  r2M1=r2-1.0,
700  s2M1=s2-1.0,
701  t2M1=t2-1.0;
702 
703  nonlinearADR_pqrstEvaluate(nPoints,
704  nSpace,
705  M,
706  A,
707  B,
708  C,
709  p1,q1,r1,s1,t1,
710  t,
711  x,
712  u,
713  m,dm,
714  f,df,
715  a,da,
716  phi,dphi,
717  r,
718  dr);
719 
720  for (k=0; k < nPoints; k++)
721  {
722  max_1mu_0 = 1.0-u[k] > 0.0 ? 1.0-u[k] : 0.0;
723 
724  if (p2 > 1.0)
725  {
726  m[k] *= pow(max_1mu_0,p2);
727  dm[k] *= pow(max_1mu_0,p2M1)*p2;
728  }
729 
730  if (q2 > 1.0)
731  {
732  for (I=0; I < nSpace; I++)
733  {
734  f[k*nSpace+I] *= pow(max_1mu_0,q2);
735  df[k*nSpace+I] *= pow(max_1mu_0,q2M1)*q2;
736  }
737  }
738 
739  if (t2 > 0.0)
740  {
741  atmp = pow(max_1mu_0,t2);
742  datmp = pow(max_1mu_0,t2M1);
743  for (I=0; I < nSpace; I++)
744  {
745  for (J=0; J < nSpace; J++)
746  {
747  a[k*nSpace2+I*nSpace+J] *= atmp;
748  da[k*nSpace2+I*nSpace+J] *= datmp*t2;
749  }
750  }
751  }
752 
753  if (r2 > 1.0)
754  {
755  phi[k] *= pow(max_1mu_0,r2);
756  dphi[k] *= pow(max_1mu_0,r2M1)*r2;
757  }
758 
759  if (s2 > 1.0)
760  {
761  r[k] *= pow(max_1mu_0,s2*r[k]);
762  dr[k] *= pow(max_1mu_0,s2M1)*s2;
763  }
764  }
765 }
766 
767 void unitSquareRotationEvaluate(const int nPoints,
768  const int nSpace,
769  const double *x,
770  const double *u,
771  double *m,
772  double *dm,
773  double *f,
774  double *df)
775 {
776  double vx, vy;
777  int k;
778  for (k=0; k < nPoints; k++)
779  {
780  m[k] = u[k];
781  dm[k] = 1.0;
782  vx = 2.0*M_PI*(x[k*3+1] - 0.5);
783  vy = 2.0*M_PI*(0.5 - x[k*3]);
784  f[k*nSpace] = vx*u[k];
785  f[k*nSpace+1] = vy*u[k];
786  df[k*nSpace] = vx;
787  df[k*nSpace+1] = vy;
788  }
789 }
790 
791 void unitCubeRotationEvaluate(const int nPoints,
792  const int nSpace,
793  const double *x,
794  const double *u,
795  double *m,
796  double *dm,
797  double *f,
798  double *df)
799 {
800  double vx, vy, vz;
801  int k;
802  for (k=0; k < nPoints; k++)
803  {
804  m[k] = u[k];
805  dm[k] = 1.0;
806  vx = 2.0*M_PI*(x[k*3+1] - 0.5);
807  vy = 2.0*M_PI*(0.5 - x[k*3]);
808  vz = 0.0;
809  f[k*nSpace] = vx*u[k];
810  f[k*nSpace+1] = vy*u[k];
811  f[k*nSpace+2] = vz*u[k];
812  df[k*nSpace] = vx;
813  df[k*nSpace+1] = vy;
814  df[k*nSpace+2] = vz;
815  }
816 }
817 
818 void rotatingPulseVelEvaluate(const int nPoints,
819  const int nSpace,
820  const double self_a,
821  const double *x,
822  const double *u,
823  double *m,
824  double *dm,
825  double *f,
826  double *df,
827  double *a,
828  double *da,
829  double *phi,
830  double *dphi)
831 {
832  /*mwf add variable declarations*/
833  double vx,vy;
834  int k,I;
835  const int nSpace2 = nSpace*nSpace;
836  memset(da, 0, nPoints * nSpace2 * sizeof(double));
837  memcpy(m, u, nPoints * sizeof(double));
838  memcpy(phi, u, nPoints * sizeof(double));
839 
840 
841  for (k=0; k < nPoints; k++)
842  {
843  dm[k] = dphi[k] = 1.0;
844 
845  vx = -4.0*(x[k*3+1] - 0.5);
846  vy = 4.0*(x[k*3] - 0.5);
847  f[k*nSpace] = vx*u[k];
848  f[k*nSpace+1] = vy*u[k];
849  df[k*nSpace] = vx;
850  df[k*nSpace+1] = vy;
851 
852  for (I=0; I < nSpace; I++)
853  {
854  a[k*nSpace2 + I*nSpace + I] = self_a;
855  }
856  }
857 }
858 
859 void disRotatingPulseVelEvaluate(const int nPoints,
860  const int nSpace,
861  const double self_a,
862  const double *x,
863  const double *u,
864  double *m,
865  double *dm,
866  double *f,
867  double *df,
868  double *a,
869  double *da,
870  double *phi,
871  double *dphi)
872 {
873  double X,Y,vx,vy;
874  int k,I;
875  const int nSpace2 = nSpace*nSpace;
876 
877  memcpy(m, u, nPoints * sizeof(double));
878  memcpy(phi, u, nPoints * sizeof(double));
879  memset(da, 0, nPoints * nSpace2 * sizeof(double));
880 
881  for (k=0; k < nPoints; k++)
882  {
883  dm[k] = dphi[k] = 1.0;
884 
885  X = x[k*3+1] - 0.5;
886  Y = x[k*3] - 0.5;
887  vx = -4.0*X;
888  vy = 4.0*Y;
889  f[k*nSpace] = vx*u[k];
890  f[k*nSpace+1] = vy*u[k];
891  df[k*nSpace] = vx;
892  df[k*nSpace+1] = vy;
893  if (sqrt(X*X+Y*Y) < 0.25)
894  {
895  f[k*nSpace] *= 0.001;
896  f[k*nSpace+1] *= 0.001;
897  df[k*nSpace] *= 0.001;
898  df[k*nSpace+1] *= 0.001;
899  }
900 
901  for (I=0; I < nSpace; I++)
902  {
903  a[k*nSpace2 + I*nSpace + I] = self_a;
904  }
905  }
906 }
907 
908 void disVelEvaluate(const int nPoints,
909  const int nSpace,
910  const double self_a,
911  const double *x,
912  const double *u,
913  double *m,
914  double *dm,
915  double *f,
916  double *df,
917  double *a,
918  double *da,
919  double *phi,
920  double *dphi)
921 {
922  int k,I,J;
923  const int nSpace2 = nSpace*nSpace;
924 
925  for (k=0; k < nPoints; k++)
926  {
927  m[k] = u[k];
928  phi[k] = 0.0;
929  dm[k] = dphi[k] = 1.0;
930 
931  f[k*nSpace] = u[k];
932  f[k*nSpace+1] = 0.0;
933  df[k*nSpace] = 1.0;
934  df[k*nSpace+1] = 0.0;
935  if (x[k*3+1] > 0.5)
936  {
937  f[k*nSpace] *= 0.25;
938  df[k*nSpace] *= 0.25;
939  }
940 
941  for (I=0; I < nSpace; I++)
942  {
943  a[k*nSpace2 + I*nSpace + I] = self_a;
944 
945  for (J=0; J < nSpace; J++)
946  {
947  da[k*nSpace2 + I*nSpace + J] = 0.0;
948  }
949  }
950  }
951 }
952 
953 void burgersDiagonalVelEvaluate(const int nPoints,
954  const int nSpace,
955  const double self_a,
956  const double *self_v,
957  const double *u,
958  double *m,
959  double *dm,
960  double *f,
961  double *df,
962  double *a,
963  double *phi,
964  double *dphi)
965 {
966  double u2;
967  int k,I;
968  const int nSpace2 = nSpace*nSpace;
969  /*mwf changed to remove maxu0*/
970  for (k=0; k < nPoints; k++)
971  {
972  m[k] = phi[k] = u[k];
973  dm[k] = dphi[k] = 1.0;
974  u2 = u[k]*u[k];
975  for (I=0; I < nSpace; I++)
976  {
977  a[k*nSpace2 + I*nSpace + I] = self_a;
978  f[k*nSpace+I] = self_v[I] * u[k] * u[k] * 0.5;
979  df[k*nSpace+I]= self_v[I] * u[k];
980 
981 /* for (J=0; J < nSpace; J++) */
982 /* { */
983 /* da[k*nSpace2 + I*nSpace + J] = 0.0; */
984 /* } */
985  }
986  }
987 }
988 
989 void burgersDiagonalVelHJEvaluate(const int nPoints,
990  const int nSpace,
991  const double self_a,
992  const double *self_v,
993  const double *u,
994  const double *grad_u,
995  double *m,
996  double *dm,
997  double *H,
998  double *dH,
999  double *a,
1000  double *phi,
1001  double *dphi)
1002 {
1003  double u2;
1004  int k,I;
1005  const int nSpace2 = nSpace*nSpace;
1006  /*mwf changed to remove maxu0*/
1007  for (k=0; k < nPoints; k++)
1008  {
1009  m[k] = phi[k] = u[k];
1010  dm[k] = dphi[k] = 1.0;
1011  u2 = u[k]*u[k];
1012  H[k] = 0.0;
1013  for (I=0; I < nSpace; I++)
1014  {
1015  a[k*nSpace2 + I*nSpace + I] = self_a;
1016  H[k] += self_v[I] * grad_u[k*nSpace+I] * u[k];
1017  dH[k*nSpace+I]= self_v[I] * u[k];
1018 
1019 /* for (J=0; J < nSpace; J++) */
1020 /* { */
1021 /* da[k*nSpace2 + I*nSpace + J] = 0.0; */
1022 /* } */
1023  }
1024  }
1025 }
1026 
1027 void twophasePotentialFlowEvaluate(int nPoints,
1028  int nSpace,
1029  double *M,
1030  double *A,
1031  double *B,
1032  double *Bcon,
1033  double *C,
1034  double t,
1035  double *x,
1036  double *u,
1037  double *m,
1038  double *dm,
1039  double *f,
1040  double *df,
1041  double *a,
1042  double *da,
1043  double *phi,
1044  double *dphi,
1045  double *r,
1046  double *dr)
1047 {
1048  int k,I,J;
1049  const int nSpace2=nSpace*nSpace;
1050  for (k=0;k<nPoints;k++)
1051  {
1052  m[k]=M[k]*u[k];
1053  dm[k]=M[k];
1054  for (I=0;I<nSpace;I++)
1055  {
1056  f[k*nSpace+I]=B[k*nSpace+I]*u[k] + Bcon[k*nSpace+I];
1057  df[k*nSpace+I]=B[k*nSpace+I];
1058  for (J=0;J<nSpace;J++)
1059  {
1060  a[k*nSpace2+I*nSpace+J]=A[k*nSpace2 + I*nSpace+J];
1061  da[k*nSpace2+I*nSpace+J]=0.0;
1062  }
1063  }
1064 
1065  phi[k]=u[k];
1066  dphi[k]=1.0;
1067 
1068  r[k]=C[k]*u[k];
1069  dr[k]=C[k];
1070  }
1071 }
1072 
1073 void twophasePotentialFlowUpdateFreeSurface(int nPoints,
1074  int nSpace,
1075  double eps,
1076  double* u_levelSet,
1077  double M1, double M2, double *M,
1078  double* A1, double* A2, double *A,
1079  double* B1, double* B2, double *B,
1080  double* Bcon1, double* Bcon2, double *Bcon,
1081  double C1, double C2, double *C)
1082 {
1083  int k,I,J;
1084  const int nSpace2=nSpace*nSpace;
1085  double H,oneMinusH;
1086  for (k=0;k<nPoints;k++)
1087  {
1088  if (u_levelSet[k] > eps)
1089  {
1090  M[k] = M1;
1091  for (I=0;I<nSpace;I++)
1092  {
1093  B[k*nSpace+I]=B1[I];
1094  Bcon[k*nSpace+I]=Bcon1[I];
1095  for (J=0;J<nSpace;J++)
1096  {
1097  A[k*nSpace2+I*nSpace+J]=A1[I*nSpace+J];
1098  }
1099  }
1100  C[k] = C1;
1101  }
1102  else if (u_levelSet[k] < -eps)
1103  {
1104  M[k] = M2;
1105  for (I=0;I<nSpace;I++)
1106  {
1107  B[k*nSpace+I]=B2[I];
1108  Bcon[k*nSpace+I]=Bcon2[I];
1109  for (J=0;J<nSpace;J++)
1110  {
1111  A[k*nSpace2+I*nSpace+J]=A2[I*nSpace+J];
1112  }
1113  }
1114  C[k] = C2;
1115  }
1116  else
1117  {
1118  H = 0.5*(1.0 + u_levelSet[k]/eps + sin((M_PI*u_levelSet[k])/eps)/M_PI);
1119  oneMinusH=1.0-H;
1120  M[k] = oneMinusH*M2 + H*M1;
1121  for (I=0;I<nSpace;I++)
1122  {
1123  B[k*nSpace+I]=oneMinusH*B2[I] + H*B1[I];
1124  Bcon[k*nSpace+I]=oneMinusH*Bcon2[I] + H*Bcon1[I];
1125  for (J=0;J<nSpace;J++)
1126  {
1127  A[k*nSpace2+I*nSpace+J]=oneMinusH*A2[I*nSpace+J] + H*A1[I*nSpace+J];
1128  }
1129  }
1130  C[k]=oneMinusH*C2+H*C1;
1131  }
1132  }
1133 }
1134 
1135 void twophaseLevelSetCoefficientsUpdateVelocity(int nPoints,
1136  int nSpace,
1137  double v_scale,
1138  double* vIn,
1139  double* vOut)
1140 {
1141  int i,I;
1142  for (i=0;i<nPoints;i++)
1143  for (I=0;I<nSpace;I++)
1144  vOut[i*nSpace+I]=v_scale*vIn[i*nSpace+I];
1145 }
1146 
1147 void twophaseLevelSetCoefficientsEvaluate(int nPoints,
1148  int nSpace,
1149  double* B,
1150  double t,
1151  double* x,
1152  double* u,
1153  double* grad_u,
1154  double* m, double* dm,
1155  double* h, double* dh,
1156  double* rh)
1157 {
1158  int i,I;
1159  for (i=0;i<nPoints;i++)
1160  {
1161  rh[i]=0.0;
1162  h[i]=0.0;
1163  m[i]=u[i];
1164  dm[i]=1.0;
1165  for (I=0;I<nSpace;I++)
1166  {
1167  h[i] += B[i*nSpace+I]*grad_u[i*nSpace+I];
1168  dh[i*nSpace+I]=B[i*nSpace+I];
1169  }
1170  }
1171 }
1172 
1173 void twophaseLevelSetCoefficientsEvaluateCI(int nPoints,
1174  int nSpace,
1175  double* B,
1176  double t,
1177  double* x,
1178  double* u,
1179  double* m, double* dm,
1180  double* f, double* df,
1181  double* a, double* da,
1182  double* phi, double* dphi,
1183  double* r, double* dr)
1184 {
1185  int i,I;
1186  for (i=0;i<nPoints;i++)
1187  {
1188  m[i]=u[i];
1189  dm[i]=1.0;
1190  for (I=0;I<nSpace;I++)
1191  {
1192  f[i*nSpace+I] = B[i*nSpace+I]*u[i];
1193  df[i*nSpace+I]=B[i*nSpace+I];
1194  }
1195  }
1196 }
1197 
1198 void ncLevelSetCoefficientsEvaluate(int nPoints,
1199  int nSpace,
1200  double* v,
1201  double* u,
1202  double* grad_u,
1203  double* m,
1204  double* dm,
1205  double* H,
1206  double* dH)
1207 {
1208  int i,I;
1209  for (i=0;i<nPoints;i++)
1210  {
1211  m[i]=u[i];
1212  dm[i]=1.0;
1213  H[i] = 0.0;
1214  for (I=0;I<nSpace;I++)
1215  {
1216  H[i] += v[i*nSpace+I]*grad_u[i*nSpace+I];
1217  dH[i*nSpace+I] = v[i*nSpace+I];
1218  }
1219  }
1220 }
1221 
1222 void cLevelSetCoefficientsEvaluate(int nPoints,
1223  int nSpace,
1224  double* v,
1225  double* u,
1226  double* m,
1227  double* dm,
1228  double* f,
1229  double* df)
1230 {
1231  int i,I;
1232  for (i=0;i<nPoints;i++)
1233  {
1234  m[i]=u[i];
1235  dm[i]=1.0;
1236  for (I=0;I<nSpace;I++)
1237  {
1238  f[i*nSpace+I] = v[i*nSpace+I]*u[i];
1239  df[i*nSpace+I] = v[i*nSpace+I];
1240  }
1241  }
1242 }
1243 
1244 void VOFCoefficientsEvaluate(int nPoints,
1245  int nSpace,
1246  double eps,
1247  double* v,
1248  double* phi,
1249  double* u,
1250  double* m,
1251  double* dm,
1252  double* f,
1253  double* df)
1254 {
1255  int i,I;
1256  for (i=0;i<nPoints;i++)
1257  {
1258  m[i]=u[i];
1259  dm[i]=1.0;
1260  for (I=0;I<nSpace;I++)
1261  {
1262 /* f[i*nSpace+I] = v[i*nSpace+I]*smoothedHeaviside(eps,phi[i]); */
1263 /* df[i*nSpace+I] = 0.0; */
1264  f[i*nSpace+I] = v[i*nSpace+I]*u[i];
1265  df[i*nSpace+I] = v[i*nSpace+I];
1266  }
1267  }
1268 }
1269 
1270 void levelSetCurvatureCoefficientsEvaluate(int nPoints,
1271  int nSpace,
1272  double *grad_phi,
1273  double *u,
1274  double *f,
1275  double *r,
1276  double *dr)
1277 {
1278  int i,I;
1279  double norm_grad_phi;
1280  for (i=0;i<nPoints;i++)
1281  {
1282  r[i] = u[i];
1283  dr[i] = 1.0;
1284  norm_grad_phi = 0.0;
1285  for (I=0;I<nSpace;I++)
1286  norm_grad_phi += grad_phi[i*nSpace+I]*grad_phi[i*nSpace+I];
1287  norm_grad_phi = sqrt(norm_grad_phi);
1288  if (norm_grad_phi > 0.0)
1289  {
1290  for (I=0;I<nSpace;I++)
1291  f[i*nSpace+I] = grad_phi[i*nSpace+I]/norm_grad_phi;
1292  }
1293  else
1294  f[i*nSpace+I] = 0.0;
1295  }
1296 }
1297 
1298 void twophaseSignedDistanceCoefficientsUpdateSignFunction(int nPoints,
1299  double eps,
1300  double* u_levelSet,
1301  double* S)
1302 {
1303  int k;
1304  double H;
1305  for (k=0;k<nPoints;k++)
1306  {
1307  if (u_levelSet[k] > eps)
1308  S[k]= 1.0;
1309  else if (u_levelSet[k] < -eps)
1310  S[k] = -1.0;
1311  else
1312  {
1313  H = 0.5*(1.0 + u_levelSet[k]/eps + sin((M_PI*u_levelSet[k])/eps)/M_PI);
1314  S[k]=2.0*H - 1.0;
1315  }
1316  }
1317 /* for (k=0;k<nPoints;k++) */
1318 /* S[k] = u_levelSet[k]/sqrt(u_levelSet[k]*u_levelSet[k] + eps*eps); */
1319 /* if(u_levelSet[k] > 0.0) */
1320 /* S[k] = 1.0; */
1321 /* else */
1322 /* S[k] = -1.0; */
1323 }
1324 
1325 void twophaseSignedDistanceCoefficientsEvaluate(int nPoints,
1326  int nSpace,
1327  double* S,
1328  double* u,
1329  double* grad_u,
1330  double* m, double* dm,
1331  double* h, double* dh,
1332  double* rh)
1333 {
1334  int i,I;
1335  for (i=0;i<nPoints;i++)
1336  {
1337  m[i]=u[i];
1338  dm[i]=1.0;
1339  rh[i]=-1.0;
1340  h[i]=0.0;
1341  for (I=0;I<nSpace;I++)
1342  h[i] += grad_u[i*nSpace+I]*grad_u[i*nSpace+I];
1343  h[i] = sqrt(h[i]);
1344  for (I=0;I<nSpace;I++)
1345  if(h[i]>= fabs(S[i]*grad_u[i*nSpace+I])*1.0e-8)
1346  dh[i*nSpace+I] = (S[i]*grad_u[i*nSpace+I])/h[i];
1347  else
1348  dh[i*nSpace+I] = (S[i]*grad_u[i*nSpace+I])/1.0e-8;
1349  h[i]+=rh[i];
1350  h[i]*=S[i];
1351  rh[i]*=S[i];
1352  }
1353 }
1354 void eikonalEquationEvaluate(int nPoints,
1355  int nSpace,
1356  double rhs,
1357  double* u,
1358  double* grad_u,
1359  double* m,
1360  double* dm,
1361  double* H,
1362  double* dH,
1363  double* r)
1364 {
1365  int i,I;
1366  double normGradU;
1367  for (i=0;i<nPoints;i++)
1368  {
1369  m[i]=u[i];
1370  dm[i]=1.0;
1371  H[i] = 0.0;
1372  r[i]=-rhs;
1373  normGradU=0.0;
1374  for (I=0;I<nSpace;I++)
1375  normGradU+= grad_u[i*nSpace+I]*grad_u[i*nSpace+I];
1376  normGradU = sqrt(normGradU);
1377  H[i] = normGradU;
1378  for (I=0;I<nSpace;I++)
1379  {
1380  dH[i*nSpace+I] = grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1381  }/*I*/
1382  }/*i*/
1383 }
1384 void redistanceLevelSetCoefficientsEvaluate(int nPoints,
1385  int nSpace,
1386  double eps,
1387  double* u_levelSet,
1388  double* u,
1389  double* grad_u,
1390  double* m,
1391  double* dm,
1392  double* H,
1393  double* dH,
1394  double* r)
1395 {
1396  int i,I;
1397  double Si,normGradU;/* He */
1398  /*mwf debug
1399  printf("redistanceLS nPoints= %d nSpace= %d eps= %g \n",nPoints,nSpace,eps);
1400  */
1401  for (i=0;i<nPoints;i++)
1402  {
1403  m[i]=u[i];
1404  dm[i]=1.0;
1405  H[i] = 0.0;
1406  Si = -1.0+2.0*smoothedHeaviside(eps,u_levelSet[i]);
1407 /* if (u_levelSet[i] > eps) */
1408 /* Si=1.0; */
1409 /* else if (u_levelSet[i] < -eps) */
1410 /* Si=-1.0; */
1411 /* else */
1412 /* { */
1413 /* He=0.5*(1.0 + u_levelSet[i]/eps + sin(M_PI*u_levelSet[i]/eps)/M_PI); */
1414 /* Si= 2.0*He-1.0; */
1415 /* } */
1416  /*mwf now try just straight sign with small eps*/
1417  /*Si = u_levelSet[i]/sqrt(u_levelSet[i]*u_levelSet[i]+1.0e-12)*/;
1418  /*
1419  r =-S
1420  H =S*|\grad d|
1421  dH=S*\grad d/|\grad d|
1422  */
1423  r[i]=-Si;
1424  normGradU=0.0;
1425  for (I=0;I<nSpace;I++)
1426  normGradU+= grad_u[i*nSpace+I]*grad_u[i*nSpace+I];
1427  normGradU = sqrt(normGradU);
1428  H[i] = Si*normGradU;
1429  /*
1430  mwf debug what about solving with r= 0 and H=S*(|\grad d|-1)?
1431  no longer homogeneous of order 1, gets Hamiltonian wrong in stabilization
1432  */
1433  /*
1434  r[i] = 0.0;
1435  H[i] = Si*(normGradU-1.0);
1436  */
1437  for (I=0;I<nSpace;I++)
1438  {
1439  dH[i*nSpace+I] = Si*grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1440  }/*I*/
1441  }/*i*/
1442 }
1443 void redistanceLevelSetCoefficientsWithWeakPenaltyEvaluate(int nPoints,
1444  int nSpace,
1445  double eps,
1446  double lambda_penalty,
1447  double* u_levelSet,
1448  double* u,
1449  double* grad_u,
1450  double* m,
1451  double* dm,
1452  double* H,
1453  double* dH,
1454  double* r,
1455  double* dr)
1456 {
1457  int i,I;
1458  double Si,normGradU;/* He */
1459  /*mwf debug
1460  printf("redistanceLS nPoints= %d nSpace= %d eps= %g \n",nPoints,nSpace,eps);
1461  */
1462  for (i=0;i<nPoints;i++)
1463  {
1464  m[i]=u[i];
1465  dm[i]=1.0;
1466  H[i] = 0.0;
1467  Si = -1.0+2.0*smoothedHeaviside(eps,u_levelSet[i]);
1468 /* if (u_levelSet[i] > eps) */
1469 /* Si=1.0; */
1470 /* else if (u_levelSet[i] < -eps) */
1471 /* Si=-1.0; */
1472 /* else */
1473 /* { */
1474 /* He=0.5*(1.0 + u_levelSet[i]/eps + sin(M_PI*u_levelSet[i]/eps)/M_PI); */
1475 /* Si= 2.0*He-1.0; */
1476 /* } */
1477  /*mwf now try just straight sign with small eps*/
1478  /*Si = u_levelSet[i]/sqrt(u_levelSet[i]*u_levelSet[i]+1.0e-12)*/;
1479  /*
1480  r =-S
1481  H =S*|\grad d|
1482  dH=S*\grad d/|\grad d|
1483  */
1484  r[i]=-Si;
1485  normGradU=0.0;
1486  for (I=0;I<nSpace;I++)
1487  normGradU+= grad_u[i*nSpace+I]*grad_u[i*nSpace+I];
1488  normGradU = sqrt(normGradU);
1489  H[i] = Si*normGradU;
1490  /*
1491  mwf debug what about solving with r= 0 and H=S*(|\grad d|-1)?
1492  no longer homogeneous of order 1, gets Hamiltonian wrong in stabilization
1493  */
1494  /*
1495  r[i] = 0.0;
1496  H[i] = Si*(normGradU-1.0);
1497  */
1498  for (I=0;I<nSpace;I++)
1499  {
1500  dH[i*nSpace+I] = Si*grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1501  }/*I*/
1502  /*add in weak penalty*/
1503  r[i] += (u[i]-u_levelSet[i])*lambda_penalty*smoothedDirac(eps,u_levelSet[i]);
1504  dr[i] = lambda_penalty*smoothedDirac(eps,u_levelSet[i]);
1505  }/*i*/
1506 }
1507 
1508 /*
1509  redistance as before, but include lagrange multiplier source term to improve volume
1510  conservation following Sussman and Fatemi 99
1511 */
1512 void redistanceLevelSetSandFCoefficientsEvaluate(int nSimplex,
1513  int nPointsPerSimplex,
1514  int nSpace,
1515  double eps,
1516  double* u_levelSet,
1517  double* dV,
1518  double* u,
1519  double* grad_u,
1520  double* m,
1521  double* dm,
1522  double* H,
1523  double* dH,
1524  double* r)
1525 {
1526  int ie,iq,i,I;
1527  double He,Si,normGradU,
1528  dHe,Li,lambda;
1529  /*mwf debug
1530  printf("redistanceLSSandF nSimplex=%d nPointsPerSimplex= %d nSpace= %d eps= %g \n",nSimplex,
1531  nPointsPerSimplex,nSpace,eps);
1532  */
1533  for (ie=0; ie < nSimplex; ie++)
1534  {
1535  /*first loop through simplex and compute normal coefficients,
1536  then loop back through and include lagrange multiplier term
1537  */
1538  for (iq=0;iq<nPointsPerSimplex;iq++)
1539  {
1540  i = ie*nPointsPerSimplex + iq;
1541  m[i]=u[i];
1542  dm[i]=1.0;
1543  H[i] = 0.0;
1544  dHe = 0.0;
1545  if (u_levelSet[i] > eps)
1546  {
1547  Si=1.0; He = 1.0;
1548  }
1549  else if (u_levelSet[i] < -eps)
1550  {
1551  Si=-1.0; He = 0.0;
1552  }
1553  else
1554  {
1555  He =0.5*(1.0 + u_levelSet[i]/eps + sin(M_PI*u_levelSet[i]/eps)/M_PI);
1556  dHe=0.5*(0.0 + 1.0/eps + cos(M_PI*u_levelSet[i]/eps)/eps);
1557  Si= 2.0*He-1.0;
1558  }
1559  normGradU=0.0;
1560  for (I=0;I<nSpace;I++)
1561  normGradU+= grad_u[i*nSpace+I]*grad_u[i*nSpace+I];
1562  normGradU = sqrt(normGradU);
1563  /*now loop through simplex and compute averages*/
1564  r[i]=-Si;
1565  H[i] = Si*normGradU;
1566  for (I=0;I<nSpace;I++)
1567  {
1568  dH[i*nSpace+I] = Si*grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1569  }/*I*/
1570  }/*iq loop 1*/
1571  /*the following is wasteful*/
1572  lambda = 0.0;
1573  for (iq = 0; iq < nPointsPerSimplex; iq++)
1574  {
1575  i = ie*nPointsPerSimplex + iq;
1576  Li= -r[i] - H[i];/* Si(1-|gradU|) */
1577  dHe = 0.0;
1578  if (u_levelSet[i] > eps)
1579  {
1580  Si=1.0; He = 1.0;
1581  }
1582  else if (u_levelSet[i] < -eps)
1583  {
1584  Si=-1.0; He = 0.0;
1585  }
1586  else
1587  {
1588  He =0.5*(1.0 + u_levelSet[i]/eps + sin(M_PI*u_levelSet[i]/eps)/M_PI);
1589  dHe=0.5*(0.0 + 1.0/eps + cos(M_PI*u_levelSet[i]/eps)/eps);
1590  Si= 2.0*He-1.0;
1591  }
1592  normGradU = H[i]/Si;
1593  lambda -= dHe*Li*dV[i]/(dHe*dHe*normGradU+1.0e-8);
1594  }/*iq*/
1595  /*go back through and update source term*/
1596  for (iq = 0; iq < nPointsPerSimplex; iq++)
1597  {
1598  i = ie*nPointsPerSimplex + iq;
1599  dHe = 0.0; He=1.0;
1600  if (u_levelSet[i] > eps)
1601  {Si=1.0; He = 1.0;}
1602  else if (u_levelSet[i] < -eps)
1603  {Si=-1.0; He = 0.0;}
1604  else
1605  {
1606  He =0.5*(1.0 + u_levelSet[i]/eps + sin(M_PI*u_levelSet[i]/eps)/M_PI);
1607  dHe=0.5*(0.0 + 1.0/eps + cos(M_PI*u_levelSet[i]/eps)/eps);
1608  Si= 2.0*He-1.0;
1609  }
1610  normGradU = H[i]/Si;
1611  r[i] -= lambda*dHe*normGradU;/*recall r is minus right hand side*/
1612  /*mwf debug
1613  printf("redistSandF ie=%d iq=%d u=%g He=%g dHe=%g Si=%g |gradu|=%g lam=%g r=%g\n",
1614  ie,iq,u_levelSet[i],He,dHe,Si,normGradU,lambda,r[i]);
1615  */
1616  }/*iq*/
1617  }/*ie*/
1618 }/*func*/
1619 
1620 void setWeakDirichletConditionsForLevelSet(int nElements_global,
1621  int nDOF_trial_element,
1622  double epsilon_freeze_factor,
1623  const double *elementDiameter,
1624  const int * u_l2g,
1625  const double *u_dof,
1626  int * freeze_nodes_tmp,
1627  int * weakDirichletConditionFlags)
1628 {
1629  int eN,j,jj,signU,j0,J0,J;
1630  double eps;
1631  for (eN = 0; eN < nElements_global; eN++)
1632  {
1633  for (j = 0; j < nDOF_trial_element; j++)
1634  freeze_nodes_tmp[j] = 0;
1635  eps = epsilon_freeze_factor*elementDiameter[eN];
1636  signU = 0; j0 = 0;
1637  while (signU == 0 && j0 < nDOF_trial_element)
1638  {
1639  J0 = u_l2g[eN*nDOF_trial_element + j0];
1640  if (u_dof[J0] < -eps)
1641  signU = -1;
1642  else if (u_dof[J0] > eps)
1643  signU = 1;
1644  else
1645  freeze_nodes_tmp[j0] = 1;
1646  j0++;
1647  }
1648  for (j = j0; j < nDOF_trial_element; j++)
1649  {
1650  J = u_l2g[eN*nDOF_trial_element + j];
1651  if ((u_dof[J] < -eps && signU == 1) ||
1652  (u_dof[J] > eps && signU == -1))
1653  {
1654  for (jj = 0; jj < nDOF_trial_element; jj++)
1655  freeze_nodes_tmp[jj] = 1;
1656  break;
1657  }
1658  else if (fabs(u_dof[J]) < eps)
1659  freeze_nodes_tmp[j] = 1;
1660  }
1661  for (j = 0; j < nDOF_trial_element; j++)
1662  {
1663  if (freeze_nodes_tmp[j] == 1)
1664  {
1665  J = u_l2g[eN*nDOF_trial_element + j];
1666  weakDirichletConditionFlags[J] = 1;
1667  }
1668  }
1669  }//eN
1670 }
1671 
1672 void setSimpleWeakDirichletConditionsForLevelSet(int nElements_global,
1673  int nDOF_trial_element,
1674  double epsilon_freeze_factor,
1675  const double *elementDiameter,
1676  const int * u_l2g,
1677  const double *u_dof,
1678  int * freeze_nodes_tmp,
1679  int * weakDirichletConditionFlags)
1680 {
1681  int eN,j,J;
1682  double eps;
1683  for (eN = 0; eN < nElements_global; eN++)
1684  {
1685  eps = epsilon_freeze_factor*elementDiameter[eN];
1686  for (j = 0; j < nDOF_trial_element; j++)
1687  {
1688  J = u_l2g[eN*nDOF_trial_element + j];
1689  if (fabs(u_dof[J]) < eps)
1690  weakDirichletConditionFlags[J] = 1;
1691  }
1692  }
1693 }
1694 
1695 /***********************************************************************
1696  begin two phase potential flow duplication efforts by mwf
1697  ***********************************************************************/
1698 void darcySharpInterfaceFlowEvaluate(int nPoints,
1699  int nSpace,
1700  double Km, double rhoM,
1701  double Kp, double rhoP,
1702  double eps,
1703  double * gravity_u,
1704  double * u,
1705  double * gradu,
1706  double * u_levelSet,
1707  double * phi_pot,
1708  double * a,
1709  double * f,
1710  double * r,
1711  double * m,
1712  double * dphi_pot,
1713  double * da,
1714  double * df,
1715  double * dr,
1716  double * dm)
1717 {
1718  int k,I,J,nSpace2;
1719  double He;
1720 
1721  nSpace2 = nSpace*nSpace;
1722  for (k = 0; k < nPoints; k++)
1723  {
1724  m[k] = 0.0; dm[k] = 1.0;
1725  r[k] = 0.0; dr[k] = 0.0;
1726  phi_pot[k] =u[k];
1727  dphi_pot[k]=1.0;
1728 
1729  He = smoothedHeaviside(eps,u_levelSet[k]);
1730  /*mwf debug
1731  printf("u_ls[%d]= %g He=%g \n",k,u_levelSet[k],He);
1732  */
1733  for (I = 0; I < nSpace; I++)
1734  {
1735  for (J = 0; J < nSpace; J++)
1736  {
1737  if (I==J)
1738  a[k*nSpace2 + I*nSpace+J] = Km + He*(Kp-Km);
1739  else
1740  a[k*nSpace2 + I*nSpace+J] = 0.0;
1741  da[k*nSpace2 + I*nSpace+J] = 0.0;
1742  }
1743  f[k*nSpace + I] = (Km*rhoM + He*(Kp*rhoP-Km*rhoM))*gravity_u[I];
1744  df[k*nSpace + I]= 0.0;
1745  }/*I*/
1746  }/*k*/
1747 
1748 }
1749 
1750 void darcySharpInterfaceFlowImEvaluate(int nPoints,
1751  int nSpace,
1752  double Km, double rhoM,
1753  double Kp, double rhoP,
1754  double eps,
1755  double * gravity_u,
1756  double * u,
1757  double * gradu,
1758  double * u_levelSet,
1759  double * phi_pot,
1760  double * a,
1761  double * f,
1762  double * r,
1763  double * m,
1764  double * dphi_pot,
1765  double * da,
1766  double * df,
1767  double * dr,
1768  double * dm)
1769 {
1770  int k,I,J,nSpace2;
1771  double He,rhoHe;
1772  double rhoEps = 0.0; /*do not smear density*/
1773  nSpace2 = nSpace*nSpace;
1774  for (k = 0; k < nPoints; k++)
1775  {
1776  m[k] = 0.0; dm[k] = 1.0;
1777  r[k] = 0.0; dr[k] = 0.0;
1778  phi_pot[k] =u[k];
1779  dphi_pot[k]=1.0;
1780 
1781  He = smoothedHeaviside(eps,u_levelSet[k]);
1782  rhoHe = smoothedHeaviside(rhoEps,u_levelSet[k]);
1783  /*mwf debug
1784  printf("u_ls[%d]= %g He=%g \n",k,u_levelSet[k],He);
1785  */
1786  for (I = 0; I < nSpace; I++)
1787  {
1788  for (J = 0; J < nSpace; J++)
1789  {
1790  if (I==J)
1791  a[k*nSpace2 + I*nSpace+J] = Km + He*(Kp-Km);
1792  else
1793  a[k*nSpace2 + I*nSpace+J] = 0.0;
1794  da[k*nSpace2 + I*nSpace+J] = 0.0;
1795  }
1796  f[k*nSpace + I] = (Km*rhoM + He*(Kp*rhoP-Km*rhoM))*gravity_u[I];
1797  df[k*nSpace + I]= 0.0;
1798  }/*I*/
1799  }/*k*/
1800 
1801 }
1802 
1803 void Laplace_Evaluate2D(const int nPoints,
1804  double *mom_p_diff_ten,
1805  double *mom_u_diff_ten,
1806  double *mom_v_diff_ten)
1807 {
1808  int k;
1809  for (k=0; k<nPoints; k++)
1810  {
1811  mom_p_diff_ten[k*2+0] = 1.0;
1812  mom_p_diff_ten[k*2+1] = 1.0;
1813 
1814  mom_u_diff_ten[k*2+0] = 1.0;
1815  mom_u_diff_ten[k*2+1] = 1.0;
1816 
1817  mom_v_diff_ten[k*2+0] = 1.0;
1818  mom_v_diff_ten[k*2+1] = 1.0;
1819  }
1820 }
1821 
1822 void Laplace_Evaluate3D(const int nPoints,
1823  double *mom_p_diff_ten,
1824  double *mom_u_diff_ten,
1825  double *mom_v_diff_ten,
1826  double *mom_w_diff_ten)
1827 {
1828  int k;
1829  for (k=0; k<nPoints; k++)
1830  {
1831  mom_p_diff_ten[k*3+0] = 1.0;
1832  mom_p_diff_ten[k*3+1] = 1.0;
1833  mom_p_diff_ten[k*3+2] = 1.0;
1834 
1835  mom_u_diff_ten[k*3+0] = 1.0;
1836  mom_u_diff_ten[k*3+1] = 1.0;
1837  mom_u_diff_ten[k*3+2] = 1.0;
1838 
1839  mom_v_diff_ten[k*3+0] = 1.0;
1840  mom_v_diff_ten[k*3+1] = 1.0;
1841  mom_v_diff_ten[k*3+2] = 1.0;
1842 
1843  mom_w_diff_ten[k*3+0] = 1.0;
1844  mom_w_diff_ten[k*3+1] = 1.0;
1845  mom_w_diff_ten[k*3+2] = 1.0;
1846  }
1847 }
1848 
1849 
1850 void NavierStokes_2D_Evaluate(const int nPoints,
1851  const double rho,
1852  const double nu,
1853  const double *g,
1854  const double *p,
1855  const double *grad_p,
1856  const double *u,
1857  const double *v,
1858  double *mom_u_acc,
1859  double *dmom_u_acc_u,
1860  double *mom_v_acc,
1861  double *dmom_v_acc_v,
1862  double *mass_adv,
1863  double *dmass_adv_u,
1864  double *dmass_adv_v,
1865  double *mom_u_adv,
1866  double *dmom_u_adv_u,
1867  double *dmom_u_adv_v,
1868  double *mom_v_adv,
1869  double *dmom_v_adv_u,
1870  double *dmom_v_adv_v,
1871  double *mom_u_diff_ten,
1872  double *mom_v_diff_ten,
1873  double *mom_u_source,
1874  double *mom_v_source,
1875  double *mom_u_ham,
1876  double *dmom_u_ham_grad_p,
1877  double *mom_v_ham,
1878  double *dmom_v_ham_grad_p)
1879 {
1880  int k;
1881  for (k=0;k<nPoints;k++)
1882  {
1884  //momentum accumulation
1885  mom_u_acc[k]=u[k];
1886  dmom_u_acc_u[k]=1.0;
1887 
1888  mom_v_acc[k]=v[k];
1889  dmom_v_acc_v[k]=1.0;
1890 
1891  //mass advective flux
1892  mass_adv[k*2+0]=u[k];
1893  mass_adv[k*2+1]=v[k];
1894 
1895  dmass_adv_u[k*2+0]=1.0;
1896  dmass_adv_v[k*2+1]=1.0;
1897 
1898  //u momentum advective flux
1899  mom_u_adv[k*2+0]=u[k]*u[k];
1900  mom_u_adv[k*2+1]=u[k]*v[k];
1901 
1902  dmom_u_adv_u[k*2+0]=2.0*u[k];
1903  dmom_u_adv_u[k*2+1]=v[k];
1904 
1905  dmom_u_adv_v[k*2+1]=u[k];
1906 
1907  //v momentum advective_flux
1908  mom_v_adv[k*2+0]=v[k]*u[k];
1909  mom_v_adv[k*2+1]=v[k]*v[k];
1910 
1911  dmom_v_adv_u[k*2+0]=v[k];
1912 
1913  dmom_v_adv_v[k*2+0]=u[k];
1914  dmom_v_adv_v[k*2+1]=2.0*v[k];
1915 
1916  //u momentum diffusion tensor
1917  mom_u_diff_ten[k*4+0] = nu;
1918  mom_u_diff_ten[k*4+3] = nu;
1919 
1920  //v momentum diffusion tensor
1921  mom_v_diff_ten[k*4+0] = nu;
1922  mom_v_diff_ten[k*4+3] = nu;
1923 
1924  //momentum sources
1925  mom_u_source[k] = -g[0];
1926  mom_v_source[k] = -g[1];
1927 
1928  //u momentum Hamiltonian (pressure)
1929  mom_u_ham[k] = grad_p[k*2+0]/rho;
1930  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
1931 
1932  //v momentum Hamiltonian (pressure)
1933  mom_v_ham[k] = grad_p[k*2+1]/rho;
1934  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
1935  }
1936 }
1937 
1938 void NavierStokes_3D_Evaluate(const int nPoints,
1939  const double rho,
1940  const double nu,
1941  const double *g,
1942  const double *p,
1943  const double *grad_p,
1944  const double *u,
1945  const double *v,
1946  const double *w,
1947  double *mom_u_acc,
1948  double *dmom_u_acc_u,
1949  double *mom_v_acc,
1950  double *dmom_v_acc_v,
1951  double *mom_w_acc,
1952  double *dmom_w_acc_w,
1953  double *mass_adv,
1954  double *dmass_adv_u,
1955  double *dmass_adv_v,
1956  double *dmass_adv_w,
1957  double *mom_u_adv,
1958  double *dmom_u_adv_u,
1959  double *dmom_u_adv_v,
1960  double *dmom_u_adv_w,
1961  double *mom_v_adv,
1962  double *dmom_v_adv_u,
1963  double *dmom_v_adv_v,
1964  double *dmom_v_adv_w,
1965  double *mom_w_adv,
1966  double *dmom_w_adv_u,
1967  double *dmom_w_adv_v,
1968  double *dmom_w_adv_w,
1969  double *mom_u_diff_ten,
1970  double *mom_v_diff_ten,
1971  double *mom_w_diff_ten,
1972  double *mom_u_source,
1973  double *mom_v_source,
1974  double *mom_w_source,
1975  double *mom_u_ham,
1976  double *dmom_u_ham_grad_p,
1977  double *mom_v_ham,
1978  double *dmom_v_ham_grad_p,
1979  double *mom_w_ham,
1980  double *dmom_w_ham_grad_p)
1981 {
1982  int k;
1983  for (k=0;k<nPoints;k++)
1984  {
1986  //momentum accumulation
1987 
1988  mom_u_acc[k]=u[k];
1989  dmom_u_acc_u[k]=1.0;
1990 
1991  mom_v_acc[k]=v[k];
1992  dmom_v_acc_v[k]=1.0;
1993 
1994  mom_w_acc[k]=w[k];
1995  dmom_w_acc_w[k]=1.0;
1996 
1997  //mass advective flux
1998  mass_adv[k*3+0]=u[k];
1999  mass_adv[k*3+1]=v[k];
2000  mass_adv[k*3+2]=w[k];
2001 
2002  dmass_adv_u[k*3+0]=1.0;
2003  dmass_adv_v[k*3+1]=1.0;
2004  dmass_adv_w[k*3+2]=1.0;
2005 
2006  //u momentum advective flux
2007  mom_u_adv[k*3+0]=u[k]*u[k];
2008  mom_u_adv[k*3+1]=u[k]*v[k];
2009  mom_u_adv[k*3+2]=u[k]*w[k];
2010 
2011  dmom_u_adv_u[k*3+0]=2.0*u[k];
2012  dmom_u_adv_u[k*3+1]=v[k];
2013  dmom_u_adv_u[k*3+2]=w[k];
2014 
2015  dmom_u_adv_v[k*3+1]=u[k];
2016 
2017  dmom_u_adv_w[k*3+2]=u[k];
2018 
2019  //v momentum advective_flux
2020  mom_v_adv[k*3+0]=v[k]*u[k];
2021  mom_v_adv[k*3+1]=v[k]*v[k];
2022  mom_v_adv[k*3+2]=v[k]*w[k];
2023 
2024  dmom_v_adv_u[k*3+0]=v[k];
2025 
2026  dmom_v_adv_v[k*3+0]=u[k];
2027  dmom_v_adv_v[k*3+1]=2.0*v[k];
2028  dmom_v_adv_v[k*3+2]=w[k];
2029 
2030  dmom_v_adv_w[k*3+2]=v[k];
2031 
2032  //w momentum advective_flux
2033  mom_w_adv[k*3+0]=w[k]*u[k];
2034  mom_w_adv[k*3+1]=w[k]*v[k];
2035  mom_w_adv[k*3+2]=w[k]*w[k];
2036 
2037  dmom_w_adv_u[k*3+0]=w[k];
2038 
2039  dmom_w_adv_v[k*3+0]=w[k];
2040 
2041  dmom_w_adv_w[k*3+0]=u[k];
2042  dmom_w_adv_w[k*3+1]=v[k];
2043  dmom_w_adv_w[k*3+2]=2.0*w[k];
2044 
2045  //u momentum diffusion tensor
2046  mom_u_diff_ten[k*9+0] = nu;
2047  mom_u_diff_ten[k*9+4] = nu;
2048  mom_u_diff_ten[k*9+8] = nu;
2049 
2050  //v momentum diffusion tensor
2051  mom_v_diff_ten[k*9+0] = nu;
2052  mom_v_diff_ten[k*9+4] = nu;
2053  mom_v_diff_ten[k*9+8] = nu;
2054 
2055  //w momentum diffusion tensor
2056  mom_w_diff_ten[k*9+0] = nu;
2057  mom_w_diff_ten[k*9+4] = nu;
2058  mom_w_diff_ten[k*9+8] = nu;
2059 
2060  //momentum sources
2061  mom_u_source[k] = -g[0];
2062  mom_v_source[k] = -g[1];
2063  mom_w_source[k] = -g[2];
2064 
2065  //u momentum Hamiltonian (pressure)
2066  mom_u_ham[k] = grad_p[k*3+0]/rho;
2067  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
2068 
2069  //v momentum Hamiltonian (pressure)
2070  mom_v_ham[k] = grad_p[k*3+1]/rho;
2071  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
2072 
2073  //w momentum Hamiltonian (pressure)
2074  mom_w_ham[k] = grad_p[k*3+2]/rho;
2075  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
2076  }
2077 }
2078 
2079 void Stokes_2D_Evaluate(const int nPoints,
2080  const double rho,
2081  const double nu,
2082  const double *g,
2083  const double *p,
2084  const double *grad_p,
2085  const double *u,
2086  const double *v,
2087  double *mom_u_acc,
2088  double *dmom_u_acc_u,
2089  double *mom_v_acc,
2090  double *dmom_v_acc_v,
2091  double *mass_adv,
2092  double *dmass_adv_u,
2093  double *dmass_adv_v,
2094  double *mom_u_diff_ten,
2095  double *mom_v_diff_ten,
2096  double *mom_u_source,
2097  double *mom_v_source,
2098  double *mom_u_ham,
2099  double *dmom_u_ham_grad_p,
2100  double *mom_v_ham,
2101  double *dmom_v_ham_grad_p)
2102 {
2103  int k;
2104  for (k=0;k<nPoints;k++)
2105  {
2107  //momentum accumulation
2108  mom_u_acc[k]=u[k];
2109  dmom_u_acc_u[k]=1.0;
2110 
2111  mom_v_acc[k]=v[k];
2112  dmom_v_acc_v[k]=1.0;
2113 
2114  //mass advective flux
2115  mass_adv[k*2+0]=u[k];
2116  mass_adv[k*2+1]=v[k];
2117 
2118  dmass_adv_u[k*2+0]=1.0;
2119  dmass_adv_v[k*2+1]=1.0;
2120 
2121  //u momentum diffusion tensor
2122  mom_u_diff_ten[k*4+0] = nu;
2123  mom_u_diff_ten[k*4+3] = nu;
2124 
2125  //v momentum diffusion tensor
2126  mom_v_diff_ten[k*4+0] = nu;
2127  mom_v_diff_ten[k*4+3] = nu;
2128 
2129  //momentum sources
2130  mom_u_source[k] = -g[0];
2131  mom_v_source[k] = -g[1];
2132 
2133  //u momentum Hamiltonian (pressure)
2134  mom_u_ham[k] = grad_p[k*2+0]/rho;
2135  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2136 
2137  //v momentum Hamiltonian (pressure)
2138  mom_v_ham[k] = grad_p[k*2+1]/rho;
2139  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2140  }
2141 }
2142 
2143 void StokesP_2D_Evaluate(const int nPoints,
2144  const double rho,
2145  const double nu,
2146  const double *g,
2147  const double *p,
2148  const double *u,
2149  const double *v,
2150  double *mom_u_acc,
2151  double *dmom_u_acc_u,
2152  double *mom_v_acc,
2153  double *dmom_v_acc_v,
2154  double *mass_adv,
2155  double *dmass_adv_u,
2156  double *dmass_adv_v,
2157  double *mom_u_adv,
2158  double *dmom_u_adv_p,
2159  double *mom_v_adv,
2160  double *dmom_v_adv_p,
2161  double *mom_u_diff_ten,
2162  double *mom_v_diff_ten,
2163  double *mom_u_source,
2164  double *mom_v_source)
2165 {
2166  int k;
2167  for (k=0;k<nPoints;k++)
2168  {
2170  //momentum accumulation
2171  mom_u_acc[k]=u[k];
2172  dmom_u_acc_u[k]=1.0;
2173 
2174  mom_v_acc[k]=v[k];
2175  dmom_v_acc_v[k]=1.0;
2176 
2177  //mass advective flux
2178  mass_adv[k*2+0]=u[k];
2179  mass_adv[k*2+1]=v[k];
2180 
2181  dmass_adv_u[k*2+0]=1.0;
2182  dmass_adv_v[k*2+1]=1.0;
2183 
2184  //u momentum advective flux
2185  mom_u_adv[k*2+0]=p[k]/rho;
2186  dmom_u_adv_p[k*2+0]=1.0/rho;
2187 
2188  //v momentum advective flux
2189  mom_v_adv[k*2+1]=p[k]/rho;
2190  dmom_v_adv_p[k*2+1]=1.0/rho;
2191 
2192  //u momentum diffusion tensor
2193  mom_u_diff_ten[k*4+0] = nu;
2194  mom_u_diff_ten[k*4+3] = nu;
2195 
2196  //v momentum diffusion tensor
2197  mom_v_diff_ten[k*4+0] = nu;
2198  mom_v_diff_ten[k*4+3] = nu;
2199 
2200  //momentum sources
2201  mom_u_source[k] = -g[0];
2202  mom_v_source[k] = -g[1];
2203  }
2204 }
2205 
2206 void Stokes_3D_Evaluate(const int nPoints,
2207  const double rho,
2208  const double nu,
2209  const double *g,
2210  const double *p,
2211  const double *grad_p,
2212  const double *u,
2213  const double *v,
2214  const double *w,
2215  double *mom_u_acc,
2216  double *dmom_u_acc_u,
2217  double *mom_v_acc,
2218  double *dmom_v_acc_v,
2219  double *mom_w_acc,
2220  double *dmom_w_acc_w,
2221  double *mass_adv,
2222  double *dmass_adv_u,
2223  double *dmass_adv_v,
2224  double *dmass_adv_w,
2225  double *mom_u_diff_ten,
2226  double *mom_v_diff_ten,
2227  double *mom_w_diff_ten,
2228  double *mom_u_source,
2229  double *mom_v_source,
2230  double *mom_w_source,
2231  double *mom_u_ham,
2232  double *dmom_u_ham_grad_p,
2233  double *mom_v_ham,
2234  double *dmom_v_ham_grad_p,
2235  double *mom_w_ham,
2236  double *dmom_w_ham_grad_p)
2237 {
2238  int k;
2239  for (k=0;k<nPoints;k++)
2240  {
2242  //momentum accumulation
2243  mom_u_acc[k]=u[k];
2244  dmom_u_acc_u[k]=1.0;
2245 
2246  mom_v_acc[k]=v[k];
2247  dmom_v_acc_v[k]=1.0;
2248 
2249  mom_w_acc[k]=w[k];
2250  dmom_w_acc_w[k]=1.0;
2251 
2252  //mass advective flux
2253  mass_adv[k*3+0]=u[k];
2254  mass_adv[k*3+1]=v[k];
2255  mass_adv[k*3+2]=w[k];
2256 
2257  dmass_adv_u[k*3+0]=1.0;
2258  dmass_adv_v[k*3+1]=1.0;
2259  dmass_adv_w[k*3+2]=1.0;
2260 
2261  //u momentum diffusion tensor
2262  mom_u_diff_ten[k*9+0] = nu;
2263  mom_u_diff_ten[k*9+4] = nu;
2264  mom_u_diff_ten[k*9+8] = nu;
2265 
2266  //v momentum diffusion tensor
2267  mom_v_diff_ten[k*9+0] = nu;
2268  mom_v_diff_ten[k*9+4] = nu;
2269  mom_v_diff_ten[k*9+8] = nu;
2270 
2271  //w momentum diffusion tensor
2272  mom_w_diff_ten[k*9+0] = nu;
2273  mom_w_diff_ten[k*9+4] = nu;
2274  mom_w_diff_ten[k*9+8] = nu;
2275 
2276  //momentum sources
2277  mom_u_source[k] = -g[0];
2278  mom_v_source[k] = -g[1];
2279  mom_w_source[k] = -g[2];
2280 
2281  //u momentum Hamiltonian (pressure)
2282  mom_u_ham[k] = grad_p[k*3+0]/rho;
2283  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
2284 
2285  //v momentum Hamiltonian (pressure)
2286  mom_v_ham[k] = grad_p[k*3+1]/rho;
2287  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
2288 
2289  //w momentum Hamiltonian (pressure)
2290  mom_w_ham[k] = grad_p[k*3+2]/rho;
2291  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
2292  }
2293 }
2294 
2295 void StokesP_3D_Evaluate(const int nPoints,
2296  const double rho,
2297  const double nu,
2298  const double *g,
2299  const double *p,
2300  const double *u,
2301  const double *v,
2302  const double *w,
2303  double *mom_u_acc,
2304  double *dmom_u_acc_u,
2305  double *mom_v_acc,
2306  double *dmom_v_acc_v,
2307  double *mom_w_acc,
2308  double *dmom_w_acc_w,
2309  double *mass_adv,
2310  double *dmass_adv_u,
2311  double *dmass_adv_v,
2312  double *dmass_adv_w,
2313  double *mom_u_adv,
2314  double *dmom_u_adv_p,
2315  double *mom_v_adv,
2316  double *dmom_v_adv_p,
2317  double *mom_w_adv,
2318  double *dmom_w_adv_p,
2319  double *mom_u_diff_ten,
2320  double *mom_v_diff_ten,
2321  double *mom_w_diff_ten,
2322  double *mom_u_source,
2323  double *mom_v_source,
2324  double *mom_w_source)
2325 {
2326  int k;
2327  for (k=0;k<nPoints;k++)
2328  {
2330  //momentum accumulation
2331  mom_u_acc[k]=u[k];
2332  dmom_u_acc_u[k]=1.0;
2333 
2334  mom_v_acc[k]=v[k];
2335  dmom_v_acc_v[k]=1.0;
2336 
2337  mom_w_acc[k]=w[k];
2338  dmom_w_acc_w[k]=1.0;
2339 
2340  //mass advective flux
2341  mass_adv[k*3+0]=u[k];
2342  mass_adv[k*3+1]=v[k];
2343  mass_adv[k*3+2]=w[k];
2344 
2345  dmass_adv_u[k*3+0]=1.0;
2346  dmass_adv_v[k*3+1]=1.0;
2347  dmass_adv_w[k*3+2]=1.0;
2348 
2349  //u momentum advective flux
2350  mom_u_adv[k*3+0]=p[k]/rho;
2351  dmom_u_adv_p[k*3+0]=1.0/rho;
2352 
2353  //v momentum advective flux
2354  mom_v_adv[k*3+1]=p[k]/rho;
2355  dmom_v_adv_p[k*3+1]=1.0/rho;
2356 
2357  //w momentum advective flux
2358  mom_w_adv[k*3+2]=p[k]/rho;
2359  dmom_w_adv_p[k*3+2]=1.0/rho;
2360 
2361  //u momentum diffusion tensor
2362  mom_u_diff_ten[k*9+0] = nu;
2363  mom_u_diff_ten[k*9+4] = nu;
2364  mom_u_diff_ten[k*9+8] = nu;
2365 
2366  //v momentum diffusion tensor
2367  mom_v_diff_ten[k*9+0] = nu;
2368  mom_v_diff_ten[k*9+4] = nu;
2369  mom_v_diff_ten[k*9+8] = nu;
2370 
2371  //w momentum diffusion tensor
2372  mom_w_diff_ten[k*9+0] = nu;
2373  mom_w_diff_ten[k*9+4] = nu;
2374  mom_w_diff_ten[k*9+8] = nu;
2375 
2376  //momentum sources
2377  mom_u_source[k] = -g[0];
2378  mom_v_source[k] = -g[1];
2379  mom_w_source[k] = -g[2];
2380  }
2381 }
2382 
2383 void TwophaseNavierStokes_LS_SO_2D_Evaluate(const int nPoints,
2384  const double eps,
2385  const double rho_0,
2386  const double nu_0,
2387  const double rho_1,
2388  const double nu_1,
2389  const double* g,
2390  const double* phi,
2391  const double *p,
2392  const double *grad_p,
2393  const double *u,
2394  const double *v,
2395  double *mom_u_acc,
2396  double *dmom_u_acc_u,
2397  double *mom_v_acc,
2398  double *dmom_v_acc_v,
2399  double *mass_adv,
2400  double *dmass_adv_u,
2401  double *dmass_adv_v,
2402  double *mom_u_adv,
2403  double *dmom_u_adv_u,
2404  double *dmom_u_adv_v,
2405  double *mom_v_adv,
2406  double *dmom_v_adv_u,
2407  double *dmom_v_adv_v,
2408  double *mom_u_diff_ten,
2409  double *mom_v_diff_ten,
2410  double *mom_u_source,
2411  double *mom_v_source,
2412  double *mom_u_ham,
2413  double *dmom_u_ham_grad_p,
2414  double *mom_v_ham,
2415  double *dmom_v_ham_grad_p)
2416 {
2417  int k;
2418  double rho,nu,H;
2419  for (k=0;k<nPoints;k++)
2420  {
2422  H = smoothedHeaviside(eps,phi[k]);
2423  rho = rho_0*(1.0-H)+rho_1*H;
2424  nu = nu_0*(1.0-H)+nu_1*H;
2425 
2426  //u momentum accumulation
2427  mom_u_acc[k]=u[k];
2428  dmom_u_acc_u[k]=1.0;
2429 
2430  //v momentum accumulation
2431  mom_v_acc[k]=v[k];
2432  dmom_v_acc_v[k]=1.0;
2433 
2434  //mass advective flux
2435  mass_adv[k*2+0]=u[k];
2436  mass_adv[k*2+1]=v[k];
2437 
2438  dmass_adv_u[k*2+0]=1.0;
2439  dmass_adv_v[k*2+1]=1.0;
2440 
2441  //u momentum advective flux
2442  mom_u_adv[k*2+0]=u[k]*u[k];
2443  mom_u_adv[k*2+1]=u[k]*v[k];
2444 
2445  dmom_u_adv_u[k*2+0]=2.0*u[k];
2446  dmom_u_adv_u[k*2+1]=v[k];
2447 
2448  dmom_u_adv_v[k*2+1]=u[k];
2449 
2450  //v momentum advective_flux
2451  mom_v_adv[k*2+0]=v[k]*u[k];
2452  mom_v_adv[k*2+1]=v[k]*v[k];
2453 
2454  dmom_v_adv_u[k*2+0]=v[k];
2455 
2456  dmom_v_adv_v[k*2+0]=u[k];
2457  dmom_v_adv_v[k*2+1]=2.0*v[k];
2458 
2459  //u momentum diffusion tensor
2460  mom_u_diff_ten[k*4+0] = nu;
2461  mom_u_diff_ten[k*4+3] = nu;
2462 
2463  //v momentum diffusion tensor
2464  mom_v_diff_ten[k*4+0] = nu;
2465  mom_v_diff_ten[k*4+3] = nu;
2466 
2467  //momentum sources
2468  mom_u_source[k] = -g[0];
2469  mom_v_source[k] = -g[1];
2470 
2471  //u momentum Hamiltonian (pressure)
2472  mom_u_ham[k] = grad_p[k*2+0]/rho;
2473  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2474 
2475  //v momentum Hamiltonian (pressure)
2476  mom_v_ham[k] = grad_p[k*2+1]/rho;
2477  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2478  }
2479 }
2480 
2481 void TwophaseNavierStokes_ST_LS_SO_2D_Evaluate(const int nPoints,
2482  const double eps_rho,
2483  const double eps_mu,
2484  const double sigma,
2485  const double rho_0,
2486  const double nu_0,
2487  const double rho_1,
2488  const double nu_1,
2489  const double* g,
2490  const double* phi,
2491  const double* n,
2492  const double* kappa,
2493  const double *p,
2494  const double *grad_p,
2495  const double *u,
2496  const double *v,
2497  double *mom_u_acc,
2498  double *dmom_u_acc_u,
2499  double *mom_v_acc,
2500  double *dmom_v_acc_v,
2501  double *mass_adv,
2502  double *dmass_adv_u,
2503  double *dmass_adv_v,
2504  double *mom_u_adv,
2505  double *dmom_u_adv_u,
2506  double *dmom_u_adv_v,
2507  double *mom_v_adv,
2508  double *dmom_v_adv_u,
2509  double *dmom_v_adv_v,
2510  double *mom_u_diff_ten,
2511  double *mom_v_diff_ten,
2512  double *mom_uv_diff_ten,
2513  double *mom_vu_diff_ten,
2514  double *mom_u_source,
2515  double *mom_v_source,
2516  double *mom_u_ham,
2517  double *dmom_u_ham_grad_p,
2518  double *mom_v_ham,
2519  double *dmom_v_ham_grad_p)
2520 {
2521  int k;
2522  double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
2523  for (k=0;k<nPoints;k++)
2524  {
2526  H_rho = smoothedHeaviside(eps_rho,phi[k]);
2527  d_rho = smoothedDirac(eps_rho,phi[k]);
2528  H_mu = smoothedHeaviside(eps_mu,phi[k]);
2529  d_mu = smoothedDirac(eps_mu,phi[k]);
2530 
2531  rho = rho_0*(1.0-H_rho)+rho_1*H_rho;
2532  nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
2533  mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
2534 
2535 /* //u momentum accumulation */
2536 /* mom_u_acc[k]=rho*u[k]; */
2537 /* dmom_u_acc_u[k]=rho; */
2538 
2539 /* //v momentum accumulation */
2540 /* mom_v_acc[k]=rho*v[k]; */
2541 /* dmom_v_acc_v[k]=rho; */
2542 
2543 /* //mass advective flux */
2544 /* mass_adv[k*2+0]=u[k]; */
2545 /* mass_adv[k*2+1]=v[k]; */
2546 
2547 /* dmass_adv_u[k*2+0]=1.0; */
2548 /* dmass_adv_v[k*2+1]=1.0; */
2549 
2550 /* //u momentum advective flux */
2551 /* mom_u_adv[k*2+0]=rho*u[k]*u[k]; */
2552 /* mom_u_adv[k*2+1]=rho*u[k]*v[k]; */
2553 
2554 /* dmom_u_adv_u[k*2+0]=2.0*rho*u[k]; */
2555 /* dmom_u_adv_u[k*2+1]=rho*v[k]; */
2556 
2557 /* dmom_u_adv_v[k*2+1]=rho*u[k]; */
2558 
2559 /* //v momentum advective_flux */
2560 /* mom_v_adv[k*2+0]=rho*v[k]*u[k]; */
2561 /* mom_v_adv[k*2+1]=rho*v[k]*v[k]; */
2562 
2563 /* dmom_v_adv_u[k*2+0]=rho*v[k]; */
2564 
2565 /* dmom_v_adv_v[k*2+0]=rho*u[k]; */
2566 /* dmom_v_adv_v[k*2+1]=2.0*rho*v[k]; */
2567 
2568 /* #ifdef SCALAR_DIFFUSION */
2569 /* //u momentum diffusion tensor */
2570 /* mom_u_diff_ten[k*4+0] = mu; */
2571 /* mom_u_diff_ten[k*4+3] = mu; */
2572 
2573 /* //v momentum diffusion tensor */
2574 /* mom_v_diff_ten[k*4+0] = mu; */
2575 /* mom_v_diff_ten[k*4+3] = mu; */
2576 /* #else */
2577 /* //u momentum diffusion tensor */
2578 /* mom_u_diff_ten[k*4+0] = 2.0*mu; */
2579 /* mom_u_diff_ten[k*4+3] = mu; */
2580 /* mom_uv_diff_ten[k*4+2]=mu; */
2581 
2582 /* //v momentum diffusion tensor */
2583 /* mom_v_diff_ten[k*4+0] = mu; */
2584 /* mom_v_diff_ten[k*4+3] = 2.0*mu; */
2585 /* mom_vu_diff_ten[k*4+1] = mu; */
2586 /* #endif */
2587 
2588 /* //momentum sources */
2589 /* norm_n = sqrt(n[k*2+0]*n[k*2+0]+n[k*2+1]*n[k*2+1]); */
2590 /* if (norm_n < 1.0e-8) */
2591 /* norm_n = 1.0e-8; */
2592 /* mom_u_source[k] = -rho*g[0] - d_mu*sigma*kappa[k]*n[k*2+0]/(norm_n); */
2593 /* mom_v_source[k] = -rho*g[1] - d_mu*sigma*kappa[k]*n[k*2+1]/(norm_n); */
2594 
2595 /* //u momentum Hamiltonian (pressure) */
2596 
2597 /* mom_u_ham[k] = grad_p[k*2+0]; */
2598 /* dmom_u_ham_grad_p[k*2+0]=1.0; */
2599 
2600 /* //v momentum Hamiltonian (pressure) */
2601 /* mom_v_ham[k] = grad_p[k*2+1]; */
2602 /* dmom_v_ham_grad_p[k*21]=1.0; */
2603 
2604  //cek incomp form
2605  //u momentum accumulation
2606  mom_u_acc[k]=u[k];
2607  dmom_u_acc_u[k]=1.0;
2608 
2609  //v momentum accumulation
2610  mom_v_acc[k]=v[k];
2611  dmom_v_acc_v[k]=1.0;
2612 
2613  //mass advective flux
2614  mass_adv[k*2+0]=u[k];
2615  mass_adv[k*2+1]=v[k];
2616 
2617  dmass_adv_u[k*2+0]=1.0;
2618  dmass_adv_v[k*2+1]=1.0;
2619 
2620  //u momentum advective flux
2621  mom_u_adv[k*2+0]=u[k]*u[k];
2622  mom_u_adv[k*2+1]=u[k]*v[k];
2623 
2624  dmom_u_adv_u[k*2+0]=2.0*u[k];
2625  dmom_u_adv_u[k*2+1]=v[k];
2626 
2627  dmom_u_adv_v[k*2+1]=u[k];
2628 
2629  //v momentum advective_flux
2630  mom_v_adv[k*2+0]=v[k]*u[k];
2631  mom_v_adv[k*2+1]=v[k]*v[k];
2632 
2633  dmom_v_adv_u[k*2+0]=v[k];
2634 
2635  dmom_v_adv_v[k*2+0]=u[k];
2636  dmom_v_adv_v[k*2+1]=2.0*v[k];
2637 
2638 #ifdef SCALAR_DIFFUSION
2639  //u momentum diffusion tensor
2640  mom_u_diff_ten[k*4+0] = nu;
2641  mom_u_diff_ten[k*4+3] = nu;
2642 
2643  //v momentum diffusion tensor
2644  mom_v_diff_ten[k*4+0] = nu;
2645  mom_v_diff_ten[k*4+3] = nu;
2646 #else
2647  //u momentum diffusion tensor
2648  mom_u_diff_ten[k*4+0] = 2.0*nu;
2649  mom_u_diff_ten[k*4+3] = nu;
2650  mom_uv_diff_ten[k*4+2]=nu;
2651 
2652  //v momentum diffusion tensor
2653  mom_v_diff_ten[k*4+0] = nu;
2654  mom_v_diff_ten[k*4+3] = 2.0*nu;
2655  mom_vu_diff_ten[k*4+1] = nu;
2656 #endif
2657 
2658  //momentum sources
2659 /* mom_u_source[k] = -g[0]; */
2660 /* mom_v_source[k] = -g[1]; */
2661  norm_n = sqrt(n[k*2+0]*n[k*2+0]+n[k*2+1]*n[k*2+1]);
2662  mom_u_source[k] = -g[0] - d_mu*sigma*kappa[k]*n[k*2+0]/(rho*(norm_n+1.0e-8));
2663  mom_v_source[k] = -g[1] - d_mu*sigma*kappa[k]*n[k*2+1]/(rho*(norm_n+1.0e-8));
2664 
2665  //u momentum Hamiltonian (pressure)
2666 
2667  mom_u_ham[k] = grad_p[k*2+0]/rho;
2668  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2669 
2670  //v momentum Hamiltonian (pressure)
2671  mom_v_ham[k] = grad_p[k*2+1]/rho;
2672  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2673  }
2674 }
2675 
2676 
2677 void TwophaseNavierStokes_ST_LS_SO_2D_Evaluate_sd(const int nPoints,
2678  const double eps_rho,
2679  const double eps_mu,
2680  const double sigma,
2681  const double rho_0,
2682  const double nu_0,
2683  const double rho_1,
2684  const double nu_1,
2685  const double* g,
2686  const double* phi,
2687  const double* n,
2688  const double* kappa,
2689  const double *p,
2690  const double *grad_p,
2691  const double *u,
2692  const double *v,
2693  double *mom_u_acc,
2694  double *dmom_u_acc_u,
2695  double *mom_v_acc,
2696  double *dmom_v_acc_v,
2697  double *mass_adv,
2698  double *dmass_adv_u,
2699  double *dmass_adv_v,
2700  double *mom_u_adv,
2701  double *dmom_u_adv_u,
2702  double *dmom_u_adv_v,
2703  double *mom_v_adv,
2704  double *dmom_v_adv_u,
2705  double *dmom_v_adv_v,
2706  double *mom_u_diff_ten,
2707  double *mom_v_diff_ten,
2708  double *mom_uv_diff_ten,
2709  double *mom_vu_diff_ten,
2710  double *mom_u_source,
2711  double *mom_v_source,
2712  double *mom_u_ham,
2713  double *dmom_u_ham_grad_p,
2714  double *mom_v_ham,
2715  double *dmom_v_ham_grad_p)
2716 {
2717  int k;
2718  double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
2719  for (k=0;k<nPoints;k++)
2720  {
2722  H_rho = smoothedHeaviside(eps_rho,phi[k]);
2723  d_rho = smoothedDirac(eps_rho,phi[k]);
2724  H_mu = smoothedHeaviside(eps_mu,phi[k]);
2725  d_mu = smoothedDirac(eps_mu,phi[k]);
2726 
2727  rho = rho_0*(1.0-H_rho)+rho_1*H_rho;
2728  nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
2729  mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
2730 
2731  //u momentum accumulation
2732  mom_u_acc[k]=u[k];
2733  dmom_u_acc_u[k]=1.0;
2734 
2735  //v momentum accumulation
2736  mom_v_acc[k]=v[k];
2737  dmom_v_acc_v[k]=1.0;
2738 
2739  //mass advective flux
2740  mass_adv[k*2+0]=u[k];
2741  mass_adv[k*2+1]=v[k];
2742 
2743  dmass_adv_u[k*2+0]=1.0;
2744  dmass_adv_v[k*2+1]=1.0;
2745 
2746  //u momentum advective flux
2747  mom_u_adv[k*2+0]=u[k]*u[k];
2748  mom_u_adv[k*2+1]=u[k]*v[k];
2749 
2750  dmom_u_adv_u[k*2+0]=2.0*u[k];
2751  dmom_u_adv_u[k*2+1]=v[k];
2752 
2753  dmom_u_adv_v[k*2+1]=u[k];
2754 
2755  //v momentum advective_flux
2756  mom_v_adv[k*2+0]=v[k]*u[k];
2757  mom_v_adv[k*2+1]=v[k]*v[k];
2758 
2759  dmom_v_adv_u[k*2+0]=v[k];
2760 
2761  dmom_v_adv_v[k*2+0]=u[k];
2762  dmom_v_adv_v[k*2+1]=2.0*v[k];
2763 
2764  //u momentum diffusion tensor
2765  mom_u_diff_ten[k*2+0] = 2.0*nu;
2766  mom_u_diff_ten[k*2+1] = nu;
2767  mom_uv_diff_ten[k]=nu;
2768 
2769  //v momentum diffusion tensor
2770  mom_v_diff_ten[k*2+0] = nu;
2771  mom_v_diff_ten[k*2+1] = 2.0*nu;
2772  mom_vu_diff_ten[k] = nu;
2773 
2774  //momentum sources
2775  /* printf("rho = %.2f\n",rho); */
2776  /* printf("mom_u_source = %f\n",mom_u_source[k]); */
2777  /* printf("mom_v_source = %f\n",mom_u_source[k]); */
2778  /* printf("d_mu = %f\n" , d_mu); */
2779  /* printf("sigma = %f\n " , sigma); */
2780 
2781  norm_n = sqrt(n[k*2+0]*n[k*2+0]+n[k*2+1]*n[k*2+1]);
2782  mom_u_source[k] = -g[0] - d_mu*sigma*kappa[k]*n[k*2+0]/(rho*(norm_n+1.0e-8));
2783  mom_v_source[k] = -g[1] - d_mu*sigma*kappa[k]*n[k*2+1]/(rho*(norm_n+1.0e-8));
2784  //u momentum Hamiltonian (pressure)
2785 
2786  mom_u_ham[k] = grad_p[k*2+0]/rho;
2787  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2788 
2789  //v momentum Hamiltonian (pressure)
2790  mom_v_ham[k] = grad_p[k*2+1]/rho;
2791  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2792 
2793  /* //compressible form */
2794  /* //u momentum accumulation */
2795  /* mom_u_acc[k]=rho*u[k]; */
2796  /* dmom_u_acc_u[k]=rho; */
2797 
2798  /* //v momentum accumulation */
2799  /* mom_v_acc[k]=rho*v[k]; */
2800  /* dmom_v_acc_v[k]=rho; */
2801 
2802  /* //mass advective flux */
2803  /* mass_adv[k*2+0]=u[k]; */
2804  /* mass_adv[k*2+1]=v[k]; */
2805 
2806  /* dmass_adv_u[k*2+0]=1.0; */
2807  /* dmass_adv_v[k*2+1]=1.0; */
2808 
2809  /* //u momentum advective flux */
2810  /* mom_u_adv[k*2+0]=rho*u[k]*u[k]; */
2811  /* mom_u_adv[k*2+1]=rho*u[k]*v[k]; */
2812 
2813  /* dmom_u_adv_u[k*2+0]=rho*2.0*u[k]; */
2814  /* dmom_u_adv_u[k*2+1]=rho*v[k]; */
2815 
2816  /* dmom_u_adv_v[k*2+1]=rho*u[k]; */
2817 
2818  /* //v momentum advective_flux */
2819  /* mom_v_adv[k*2+0]=rho*v[k]*u[k]; */
2820  /* mom_v_adv[k*2+1]=rho*v[k]*v[k]; */
2821 
2822  /* dmom_v_adv_u[k*2+0]=rho*v[k]; */
2823 
2824  /* dmom_v_adv_v[k*2+0]=rho*u[k]; */
2825  /* dmom_v_adv_v[k*2+1]=rho*2.0*v[k]; */
2826 
2827  /* //u momentum diffusion tensor */
2828  /* mom_u_diff_ten[k*2+0] = 2.0*mu; */
2829  /* mom_u_diff_ten[k*2+1] = mu; */
2830  /* mom_uv_diff_ten[k]=mu; */
2831 
2832  /* //v momentum diffusion tensor */
2833  /* mom_v_diff_ten[k*2+0] = mu; */
2834  /* mom_v_diff_ten[k*2+1] = 2.0*mu; */
2835  /* mom_vu_diff_ten[k] = mu; */
2836 
2837  /* //momentum sources */
2838  /* norm_n = sqrt(n[k*2+0]*n[k*2+0]+n[k*2+1]*n[k*2+1]); */
2839  /* mom_u_source[k] = -rho*g[0] - d_mu*sigma*kappa[k]*n[k*2+0]/(norm_n+1.0e-8); */
2840  /* mom_v_source[k] = -rho*g[1] - d_mu*sigma*kappa[k]*n[k*2+1]/(norm_n+1.0e-8); */
2841 
2842  /* //u momentum Hamiltonian (pressure) */
2843 
2844  /* mom_u_ham[k] = grad_p[k*2+0]; */
2845  /* dmom_u_ham_grad_p[k*2+0]=1.0; */
2846 
2847  /* //v momentum Hamiltonian (pressure) */
2848  /* mom_v_ham[k] = grad_p[k*2+1]; */
2849  /* dmom_v_ham_grad_p[k*2+1]=1.0; */
2850  }
2851 }
2852 
2853 void ThreephaseNavierStokes_ST_LS_SO_2D_EvaluateOrig(const int nPoints,
2854  const double eps_rho,
2855  const double eps_mu,
2856  const double sigma,
2857  const double rho_0,
2858  const double nu_0,
2859  const double rho_1,
2860  const double nu_1,
2861  const double rho_s,
2862  const double nu_s,
2863  const double* g,
2864  const double* phi,
2865  const double* n,
2866  const double* kappa,
2867  const double* phi_s,
2868  const double* n_s,
2869  const double *p,
2870  const double *grad_p,
2871  const double *u,
2872  const double *v,
2873  double *mom_u_acc,
2874  double *dmom_u_acc_u,
2875  double *mom_v_acc,
2876  double *dmom_v_acc_v,
2877  double *mass_adv,
2878  double *dmass_adv_u,
2879  double *dmass_adv_v,
2880  double *mom_u_adv,
2881  double *dmom_u_adv_u,
2882  double *dmom_u_adv_v,
2883  double *mom_v_adv,
2884  double *dmom_v_adv_u,
2885  double *dmom_v_adv_v,
2886  double *mom_u_diff_ten,
2887  double *mom_v_diff_ten,
2888  double *mom_uv_diff_ten,
2889  double *mom_vu_diff_ten,
2890  double *mom_u_source,
2891  double *dmom_u_source_u,
2892  double *dmom_u_source_v,
2893  double *mom_v_source,
2894  double *dmom_v_source_u,
2895  double *dmom_v_source_v,
2896  double *mom_u_ham,
2897  double *dmom_u_ham_grad_p,
2898  double *mom_v_ham,
2899  double *dmom_v_ham_grad_p)
2900 {
2901  int k;
2902  double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,
2903  H_rho_s,d_rho_s,H_mu_s,d_mu_s,norm_n,norm_n_s;
2904  for (k=0;k<nPoints;k++)
2905  {
2907  H_rho = smoothedHeaviside(eps_rho,phi[k]);
2908  d_rho = smoothedDirac(eps_rho,phi[k]);
2909  H_mu = smoothedHeaviside(eps_mu,phi[k]);
2910  d_mu = smoothedDirac(eps_mu,phi[k]);
2911 
2912  rho = rho_0*(1.0-H_rho)+rho_1*H_rho;
2913  nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
2914  mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
2915  rho = rho_0;
2916  nu = nu_0;
2917  mu = rho_0*nu_0;
2918 
2919 /* H_rho_s = smoothedHeaviside(eps_rho,phi_s[k]); */
2920 /* d_rho_s = smoothedDirac(eps_rho,phi_s[k]); */
2921 /* H_mu_s = smoothedHeaviside(eps_mu,phi_s[k]); */
2922 /* d_mu_s = smoothedDirac(eps_mu,phi_s[k]); */
2923 
2924 /* rho = rho_s*(1.0-H_rho_s)+rho*H_rho_s; */
2925 /* nu = nu_s*(1.0-H_mu_s)+nu*H_mu_s; */
2926 /* mu = rho_s*nu_s*(1.0-H_mu_s)+rho*nu*H_mu_s; */
2927 
2928  //u momentum accumulation
2929  mom_u_acc[k]=rho*u[k];
2930  dmom_u_acc_u[k]=rho;
2931 
2932  //v momentum accumulation
2933  mom_v_acc[k]=rho*v[k];
2934  dmom_v_acc_v[k]=rho;
2935 
2936  //mass advective flux
2937  mass_adv[k*2+0]=u[k];
2938  mass_adv[k*2+1]=v[k];
2939 
2940  dmass_adv_u[k*2+0]=1.0;
2941  dmass_adv_v[k*2+1]=1.0;
2942 
2943  //u momentum advective flux
2944  mom_u_adv[k*2+0]=rho*u[k]*u[k];
2945  mom_u_adv[k*2+1]=rho*u[k]*v[k];
2946 
2947  dmom_u_adv_u[k*2+0]=2.0*rho*u[k];
2948  dmom_u_adv_u[k*2+1]=rho*v[k];
2949 
2950  dmom_u_adv_v[k*2+1]=rho*u[k];
2951 
2952  //v momentum advective_flux
2953  mom_v_adv[k*2+0]=rho*v[k]*u[k];
2954  mom_v_adv[k*2+1]=rho*v[k]*v[k];
2955 
2956  dmom_v_adv_u[k*2+0]=rho*v[k];
2957 
2958  dmom_v_adv_v[k*2+0]=rho*u[k];
2959  dmom_v_adv_v[k*2+1]=2.0*rho*v[k];
2960 
2961  //u momentum diffusion tensor
2962  mom_u_diff_ten[k*4+0] = 2.0*mu;
2963  mom_u_diff_ten[k*4+3] = mu;
2964  mom_uv_diff_ten[k*4+2]=mu;
2965 
2966  //v momentum diffusion tensor
2967  mom_v_diff_ten[k*4+0] = mu;
2968  mom_v_diff_ten[k*4+3] = 2.0*mu;
2969  mom_vu_diff_ten[k*4+1] = mu;
2970 
2971  //momentum sources
2972  norm_n = sqrt(n[k*2+0]*n[k*2+0]+n[k*2+1]*n[k*2+1]);
2973  norm_n_s = sqrt(n_s[k*2+0]*n_s[k*2+0]+n_s[k*2+1]*n_s[k*2+1]);
2974 
2975 /* mom_u_source[k] = -rho*g[0] */
2976 /* - d_mu*sigma*kappa[k]*n[k*2+0]/norm_n */
2977 /* +rho*d_mu_s*(u[k]*u[k]*n_s[k*2+0] + u[k]*v[k]*n_s[k*2+1])/norm_n_s */
2978 /* +rho*(1.0-H_mu_s)*u[k]; */
2979 /* dmom_u_source_u[k] = rho*d_mu_s*(2.0*u[k]*n_s[k*2+0] */
2980 /* + v[k]*n_s[k*2+1])/norm_n_s */
2981 /* +rho*(1.0-H_mu_s); */
2982 /* dmom_u_source_v[k] = rho*d_mu_s*u[k]*n_s[k*2+1]/norm_n_s; */
2983 
2984 /* mom_v_source[k] = -rho*g[1] */
2985 /* - d_mu*sigma*kappa[k]*n[k*2+1]/norm_n */
2986 /* +rho*d_mu_s*(v[k]*u[k]*n_s[k*2+0] */
2987 /* + v[k]*v[k]*n_s[k*2+1])/norm_n_s */
2988 /* +rho*(1.0-H_mu_s)*v[k]; */
2989 /* dmom_v_source_u[k] = rho*d_mu_s*v[k]*n_s[k*2+0]/norm_n_s; */
2990 /* dmom_v_source_v[k] = rho*d_mu_s*(u[k]*n_s[k*2+0] */
2991 /* + 2.0*v[k]*n_s[k*2+1])/norm_n_s */
2992 /* +rho*(1.0-H_mu_s); */
2993 
2994  /* mom_u_source[k] = -rho*g[0] */
2995  /* - d_mu*sigma*kappa[k]*n[k*2+0]/norm_n */
2996  /* +rho*2.0*(1.0-H_mu_s)*u[k]; */
2997  /* dmom_u_source_u[k] = rho*2.0*(1.0-H_mu_s); */
2998  /* dmom_u_source_v[k] = 0.0; */
2999 
3000  /* mom_v_source[k] = -rho*g[1] */
3001  /* - d_mu*sigma*kappa[k]*n[k*2+1]/norm_n */
3002  /* +rho*2.0*(1.0-H_mu_s)*v[k]; */
3003  /* dmom_v_source_u[k] = 0.0; */
3004  /* dmom_v_source_v[k] = rho*2.0*(1.0-H_mu_s); */
3005 
3006  /* mom_u_source[k] = -rho*g[0] */
3007  /* - d_mu*sigma*kappa[k]*n[k*2+0]/norm_n */
3008  /* - rho*d_mu_s*(u[k]*u[k]*n_s[k*2+0] + u[k]*v[k]*n_s[k*2+1])/norm_n_s; */
3009 
3010  /* dmom_u_source_u[k] = -rho*d_mu_s*(2.0*u[k]*n_s[k*2+0] */
3011  /* + v[k]*n_s[k*2+1])/norm_n_s; */
3012 
3013  /* dmom_u_source_v[k] = -rho*d_mu_s*u[k]*n_s[k*2+1]/norm_n_s; */
3014 
3015  /* mom_v_source[k] = -rho*g[1] */
3016  /* - d_mu*sigma*kappa[k]*n[k*2+1]/norm_n */
3017  /* - rho*d_mu_s*(v[k]*u[k]*n_s[k*2+0] */
3018  /* + v[k]*v[k]*n_s[k*2+1])/norm_n_s; */
3019 
3020  /* dmom_v_source_u[k] = -rho*d_mu_s*v[k]*n_s[k*2+0]/norm_n_s; */
3021  /* dmom_v_source_v[k] = -rho*d_mu_s*(u[k]*n_s[k*2+0] */
3022  /* + 2.0*v[k]*n_s[k*2+1])/norm_n_s; */
3023 
3024  /* mom_u_source[k] = -rho*g[0] +(1.0-H_mu_s)*u[k]; */
3025  /* dmom_u_source_u[k] = (1.0-H_mu_s); */
3026  /* dmom_u_source_v[k] = 0.0; */
3027 
3028  /* mom_v_source[k] = -rho*g[1] +(1.0-H_mu_s)*v[k]; */
3029  /* dmom_v_source_u[k] = 0.0; */
3030  /* dmom_v_source_v[k] = (1.0-H_mu_s); */
3031  mom_u_source[k] = -rho*g[0];
3032  dmom_u_source_u[k] = 0.0;
3033  dmom_u_source_v[k] = 0.0;
3034 
3035  mom_v_source[k] = -rho*g[1];
3036  dmom_v_source_u[k] = 0.0;
3037  dmom_v_source_v[k] = 0.0;
3038 
3039  //u momentum Hamiltonian (pressure)
3040 
3041  mom_u_ham[k] = grad_p[k*2+0];
3042  dmom_u_ham_grad_p[k*2+0]=1.0;
3043 
3044  //v momentum Hamiltonian (pressure)
3045  mom_v_ham[k] = grad_p[k*2+1];
3046  dmom_v_ham_grad_p[k*2+1]=1.0;
3047  }
3048 }
3049 
3050 void ThreephaseNavierStokes_ST_LS_SO_2D_Evaluate(const int nPoints,
3051  const double boundaryPenaltyCoef,
3052  const double volumePenaltyCoef,
3053  const double eps_rho,
3054  const double eps_mu,
3055  const double sigma,
3056  const double rho_0,
3057  const double nu_0,
3058  const double rho_1,
3059  const double nu_1,
3060  const double rho_s,
3061  const double nu_s,
3062  const double* g,
3063  const double* phi,
3064  const double* n,
3065  const double* kappa,
3066  const double* phi_s,
3067  const double* n_s,
3068  const double *p,
3069  const double *grad_p,
3070  const double *u,
3071  const double *v,
3072  double *mom_u_acc,
3073  double *dmom_u_acc_u,
3074  double *mom_v_acc,
3075  double *dmom_v_acc_v,
3076  double *mass_adv,
3077  double *dmass_adv_u,
3078  double *dmass_adv_v,
3079  double *mom_u_adv,
3080  double *dmom_u_adv_u,
3081  double *dmom_u_adv_v,
3082  double *mom_v_adv,
3083  double *dmom_v_adv_u,
3084  double *dmom_v_adv_v,
3085  double *mom_u_diff_ten,
3086  double *mom_v_diff_ten,
3087  double *mom_uv_diff_ten,
3088  double *mom_vu_diff_ten,
3089  double *mom_u_source,
3090  double *dmom_u_source_u,
3091  double *dmom_u_source_v,
3092  double *mom_v_source,
3093  double *dmom_v_source_u,
3094  double *dmom_v_source_v,
3095  double *mom_u_ham,
3096  double *dmom_u_ham_grad_p,
3097  double *mom_v_ham,
3098  double *dmom_v_ham_grad_p)
3099 {
3100  int k;
3101  double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,
3102  H_rho_s,d_rho_s,H_mu_s,d_mu_s,norm_n,norm_n_s,volumeFlux,sp;
3103  int quadraticPenalty = 0,sipgPenalty=1.0;
3104  sp=(double)(sipgPenalty);
3105  //sp=0.0;
3106  for (k=0;k<nPoints;k++)
3107  {
3109  H_rho = smoothedHeaviside(eps_rho,phi[k]);
3110  d_rho = smoothedDirac(eps_rho,phi[k]);
3111  H_mu = smoothedHeaviside(eps_mu,phi[k]);
3112  d_mu = smoothedDirac(eps_mu,phi[k]);
3113 
3114  rho = rho_0*(1.0-H_rho)+rho_1*H_rho;
3115  nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
3116  mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
3117 
3118  H_rho_s = smoothedHeaviside(eps_rho,phi_s[k]);
3119  d_rho_s = smoothedDirac(eps_rho,phi_s[k]);
3120  H_mu_s = smoothedHeaviside(eps_mu,phi_s[k]);
3121  d_mu_s = smoothedDirac(eps_mu,phi_s[k]);
3122 
3123  norm_n = sqrt(n[k*2+0]*n[k*2+0]+n[k*2+1]*n[k*2+1]);
3124  norm_n_s = sqrt(n_s[k*2+0]*n_s[k*2+0]+n_s[k*2+1]*n_s[k*2+1]);
3125 
3126  //mwf start hacking here ...
3127  //make coefficients just be fluid values first?
3128  /* rho = rho_s*(1.0-H_rho_s)+rho*H_rho_s; */
3129  /* nu = nu_s*(1.0-H_mu_s)+nu*H_mu_s; */
3130  /* mu = rho_s*nu_s*(1.0-H_mu_s)+rho*nu*H_mu_s; */
3131 
3132  //u momentum accumulation
3133  mom_u_acc[k]=H_mu_s*rho*u[k];
3134  dmom_u_acc_u[k]=H_mu_s*rho;
3135 
3136  //v momentum accumulation
3137  mom_v_acc[k]=H_mu_s*rho*v[k];
3138  dmom_v_acc_v[k]=H_mu_s*rho;
3139 
3140  //mwf volume conservation holds in both solid and fluid
3141  //mass advective flux
3142  mass_adv[k*2+0]=u[k];
3143  mass_adv[k*2+1]=v[k];
3144 
3145  dmass_adv_u[k*2+0]=1.0;
3146  dmass_adv_v[k*2+1]=1.0;
3147 
3148  //mwf continue killing acceleration in solid phase
3149  //u momentum advective flux
3150  mom_u_adv[k*2+0]=H_mu_s*rho*u[k]*u[k] - sp*(u[k]-0.0)*d_mu_s*n_s[k*2+0]/norm_n_s;
3151  mom_u_adv[k*2+1]=H_mu_s*rho*u[k]*v[k] - sp*(u[k]-0.0)*d_mu_s*n_s[k*2+1]/norm_n_s;
3152 
3153  dmom_u_adv_u[k*2+0]=2.0*H_mu_s*rho*u[k] - sp*d_mu_s*n_s[k*2+0]/norm_n_s;
3154  dmom_u_adv_u[k*2+1]=H_mu_s*rho*v[k] - sp*d_mu_s*n_s[k*2+1]/norm_n_s;
3155 
3156  dmom_u_adv_v[k*2+1]=H_mu_s*rho*u[k];
3157 
3158  //v momentum advective_flux
3159  mom_v_adv[k*2+0]=H_mu_s*rho*v[k]*u[k] - sp*(v[k]-0.0)*d_mu_s*n_s[k*2+0]/norm_n_s;
3160  mom_v_adv[k*2+1]=H_mu_s*rho*v[k]*v[k] - sp*(v[k]-0.0)*d_mu_s*n_s[k*2+1]/norm_n_s;
3161 
3162  dmom_v_adv_u[k*2+0]=H_mu_s*rho*v[k];
3163 
3164  dmom_v_adv_v[k*2+0]=H_mu_s*rho*u[k] - sp*d_mu_s*n_s[k*2+0]/norm_n_s;
3165  dmom_v_adv_v[k*2+1]=2.0*H_mu_s*rho*v[k] - sp*d_mu_s*n_s[k*2+1]/norm_n_s;
3166 
3167  //mwf no diffusion of momentum in solid phase either,
3168  //mwf also need to switch to Laplace form temporarily because of bc's?
3169  //u momentum diffusion tensor
3170  mom_u_diff_ten[k*4+0] = 2.0*H_mu_s*mu;
3171  mom_u_diff_ten[k*4+3] = H_mu_s*mu;
3172  mom_uv_diff_ten[k*4+2]=H_mu_s*mu;
3173 
3174  //v momentum diffusion tensor
3175  mom_v_diff_ten[k*4+0] = H_mu_s*mu;
3176  mom_v_diff_ten[k*4+3] = 2.0*H_mu_s*mu;
3177  mom_vu_diff_ten[k*4+1] = H_mu_s*mu;
3178 /* //u momentum diffusion tensor */
3179 /* mom_u_diff_ten[k*4+0] = H_mu_s*mu; */
3180 /* mom_u_diff_ten[k*4+3] = H_mu_s*mu; */
3181 /* mom_uv_diff_ten[k*4+2]= 0.0; */
3182 
3183 /* //v momentum diffusion tensor */
3184 /* mom_v_diff_ten[k*4+0] = H_mu_s*mu; */
3185 /* mom_v_diff_ten[k*4+3] = H_mu_s*mu; */
3186 /* mom_vu_diff_ten[k*4+1] = 0.0; */
3187 
3188  //momentum sources
3189 
3190  //mwf momentum in solid phase is just \grad p = 0 (i.e. pressure is constant) with penalties to enforce
3191  //mwf velocity is equal to input (v^s = 0 for now) on boundary and in solid region
3192  if (quadraticPenalty)
3193  {
3194  mom_u_source[k] = -H_mu_s*rho*g[0]
3195  - H_mu_s*d_mu*sigma*kappa[k]*n[k*2+0]/norm_n
3196  +boundaryPenaltyCoef*rho*d_mu_s*(u[k] - 0.0)
3197  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(u[k] - 0.0)*(u[k] - 0.0);
3198 
3199  dmom_u_source_u[k] = boundaryPenaltyCoef*rho*d_mu_s
3200  +volumePenaltyCoef*rho*(1.0-H_mu_s)*2.0*(u[k]-0.0);
3201  dmom_u_source_v[k] = 0.0;
3202 
3203  mom_v_source[k] = -H_mu_s*rho*g[1]
3204  - H_mu_s*d_mu*sigma*kappa[k]*n[k*2+1]/norm_n
3205  +boundaryPenaltyCoef*rho*d_mu_s*(v[k] - 0.0)
3206  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(v[k] - 0.0)*(v[k] - 0.0);
3207 
3208  dmom_v_source_u[k] = 0.0;
3209  dmom_v_source_v[k] = boundaryPenaltyCoef*rho*d_mu_s
3210  +volumePenaltyCoef*rho*(1.0-H_mu_s)*2.0*(v[k] - 0.0);
3211  }
3212  else if (sipgPenalty)
3213  {
3214  mom_u_source[k] = -H_mu_s*rho*g[0]
3215  - H_mu_s*d_mu*sigma*kappa[k]*n[k*2+0]/norm_n
3216  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(u[k] - 0.0);
3217 
3218  dmom_u_source_u[k] = volumePenaltyCoef*rho*(1.0-H_mu_s);
3219  dmom_u_source_v[k] = 0.0;
3220 
3221  mom_v_source[k] = -H_mu_s*rho*g[1]
3222  - H_mu_s*d_mu*sigma*kappa[k]*n[k*2+1]/norm_n
3223  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(v[k] - 0.0);
3224 
3225  dmom_v_source_u[k] = 0.0;
3226  dmom_v_source_v[k] = volumePenaltyCoef*rho*(1.0-H_mu_s);
3227 
3228  volumeFlux = u[k]*n_s[k*2+0] + v[k]*n_s[k*2+1];
3229  if (volumeFlux < 0.0)
3230  {
3231  mom_u_source[k] += rho*d_mu_s*u[k]*(u[k]*n_s[k*2+0] + v[k]*n_s[k*2+1])/norm_n_s;
3232 
3233  dmom_u_source_u[k] += rho*d_mu_s*(2.0*u[k]*n_s[k*2+0]
3234  + v[k]*n_s[k*2+1])/norm_n_s;
3235 
3236  dmom_u_source_v[k] = rho*d_mu_s*u[k]*n_s[k*2+1]/norm_n_s;
3237 
3238  mom_v_source[k] += rho*d_mu_s*(v[k]*u[k]*n_s[k*2+0]
3239  + v[k]*v[k]*n_s[k*2+1])/norm_n_s;
3240 
3241  dmom_v_source_u[k] += rho*d_mu_s*v[k]*n_s[k*2+0]/norm_n_s;
3242 
3243  dmom_v_source_v[k] += rho*d_mu_s*(u[k]*n_s[k*2+0]
3244  + 2.0*v[k]*n_s[k*2+1])/norm_n_s;
3245  }
3246  }
3247  else
3248  {
3249  mom_u_source[k] = -H_mu_s*rho*g[0]
3250  - H_mu_s*d_mu*sigma*kappa[k]*n[k*2+0]/norm_n
3251  +boundaryPenaltyCoef*rho*d_mu_s*(u[k] - 0.0)
3252  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(u[k] - 0.0);
3253 
3254  dmom_u_source_u[k] = boundaryPenaltyCoef*rho*d_mu_s
3255  +volumePenaltyCoef*rho*(1.0-H_mu_s);
3256  dmom_u_source_v[k] = 0.0;
3257 
3258  mom_v_source[k] = -H_mu_s*rho*g[1]
3259  - H_mu_s*d_mu*sigma*kappa[k]*n[k*2+1]/norm_n
3260  +boundaryPenaltyCoef*rho*d_mu_s*(v[k] - 0.0)
3261  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(v[k] - 0.0);
3262 
3263  dmom_v_source_u[k] = 0.0;
3264  dmom_v_source_v[k] = boundaryPenaltyCoef*rho*d_mu_s
3265  +volumePenaltyCoef*rho*(1.0-H_mu_s);
3266 
3267  }
3268  //mwf include grad_p term in solid phase to enforce p = constant (or integral around boundary is zero?)
3269  //u momentum Hamiltonian (pressure)
3270  mom_u_ham[k] = grad_p[k*2+0];
3271  dmom_u_ham_grad_p[k*2+0]=1.0;
3272 
3273  //v momentum Hamiltonian (pressure)
3274  mom_v_ham[k] = grad_p[k*2+1];
3275  dmom_v_ham_grad_p[k*2+1]=1.0;
3276  }
3277 }
3278 
3279 void TwophaseNavierStokes_ST_LS_SO_3D_Evaluate(const int nPoints,
3280  const double eps_rho,
3281  const double eps_mu,
3282  const double sigma,
3283  const double rho_0,
3284  const double nu_0,
3285  const double rho_1,
3286  const double nu_1,
3287  const double* g,
3288  const double* phi,
3289  const double* n,
3290  const double* kappa,
3291  const double *p,
3292  const double *grad_p,
3293  const double *u,
3294  const double *v,
3295  const double *w,
3296  double *mom_u_acc,
3297  double *dmom_u_acc_u,
3298  double *mom_v_acc,
3299  double *dmom_v_acc_v,
3300  double *mom_w_acc,
3301  double *dmom_w_acc_w,
3302  double *mass_adv,
3303  double *dmass_adv_u,
3304  double *dmass_adv_v,
3305  double *dmass_adv_w,
3306  double *mom_u_adv,
3307  double *dmom_u_adv_u,
3308  double *dmom_u_adv_v,
3309  double *dmom_u_adv_w,
3310  double *mom_v_adv,
3311  double *dmom_v_adv_u,
3312  double *dmom_v_adv_v,
3313  double *dmom_v_adv_w,
3314  double *mom_w_adv,
3315  double *dmom_w_adv_u,
3316  double *dmom_w_adv_v,
3317  double *dmom_w_adv_w,
3318  double *mom_u_diff_ten,
3319  double *mom_v_diff_ten,
3320  double *mom_w_diff_ten,
3321  double *mom_uv_diff_ten,
3322  double *mom_uw_diff_ten,
3323  double *mom_vu_diff_ten,
3324  double *mom_vw_diff_ten,
3325  double *mom_wu_diff_ten,
3326  double *mom_wv_diff_ten,
3327  double *mom_u_source,
3328  double *mom_v_source,
3329  double *mom_w_source,
3330  double *mom_u_ham,
3331  double *dmom_u_ham_grad_p,
3332  double *mom_v_ham,
3333  double *dmom_v_ham_grad_p,
3334  double *mom_w_ham,
3335  double *dmom_w_ham_grad_p)
3336 {
3337  int k;
3338  double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
3339  for (k=0;k<nPoints;k++)
3340  {
3342  /*H = smoothedHeaviside(eps,phi[k]);*/
3343  H_rho = smoothedHeaviside(eps_rho,phi[k]);
3344  d_rho = smoothedDirac(eps_rho,phi[k]);
3345  H_mu = smoothedHeaviside(eps_mu,phi[k]);
3346  d_mu = smoothedDirac(eps_mu,phi[k]);
3347 
3348  rho = rho_0*(1.0-H_rho)+rho_1*H_rho;
3349  nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
3350  mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
3351 
3352 /* //u momentum accumulation */
3353 /* mom_u_acc[k]=rho*u[k]; */
3354 /* dmom_u_acc_u[k]=rho; */
3355 
3356 /* //v momentum accumulation */
3357 /* mom_v_acc[k]=rho*v[k]; */
3358 /* dmom_v_acc_v[k]=rho; */
3359 
3360 /* //w momentum accumulation */
3361 /* mom_w_acc[k]=rho*w[k]; */
3362 /* dmom_w_acc_w[k]=rho; */
3363 
3364 
3365 /* //mass advective flux */
3366 /* mass_adv[k*3+0]=u[k]; */
3367 /* mass_adv[k*3+1]=v[k]; */
3368 /* mass_adv[k*3+2]=w[k]; */
3369 
3370 /* dmass_adv_u[k*3+0]=1.0; */
3371 /* dmass_adv_v[k*3+1]=1.0; */
3372 /* dmass_adv_w[k*3+2]=1.0; */
3373 
3374 /* //u momentum advective flux */
3375 /* mom_u_adv[k*3+0]=rho*u[k]*u[k]; */
3376 /* mom_u_adv[k*3+1]=rho*u[k]*v[k]; */
3377 /* mom_u_adv[k*3+2]=rho*u[k]*w[k]; */
3378 
3379 /* dmom_u_adv_u[k*3+0]=2.0*rho*u[k]; */
3380 /* dmom_u_adv_u[k*3+1]=rho*v[k]; */
3381 /* dmom_u_adv_u[k*3+2]=rho*w[k]; */
3382 
3383 /* dmom_u_adv_v[k*3+1]=rho*u[k]; */
3384 
3385 /* dmom_u_adv_w[k*3+2]=rho*u[k]; */
3386 
3387 /* //v momentum advective_flux */
3388 /* mom_v_adv[k*3+0]=rho*v[k]*u[k]; */
3389 /* mom_v_adv[k*3+1]=rho*v[k]*v[k]; */
3390 /* mom_v_adv[k*3+2]=rho*v[k]*w[k]; */
3391 
3392 /* dmom_v_adv_u[k*3+0]=rho*v[k]; */
3393 
3394 /* dmom_v_adv_w[k*3+2]=rho*v[k]; */
3395 
3396 /* dmom_v_adv_v[k*3+0]=rho*u[k]; */
3397 /* dmom_v_adv_v[k*3+1]=2.0*rho*v[k]; */
3398 /* dmom_v_adv_v[k*3+2]=rho*w[k]; */
3399 
3400 /* //w momentum advective_flux */
3401 /* mom_w_adv[k*3+0]=rho*w[k]*u[k]; */
3402 /* mom_w_adv[k*3+1]=rho*w[k]*v[k]; */
3403 /* mom_w_adv[k*3+2]=rho*w[k]*w[k]; */
3404 
3405 /* dmom_w_adv_u[k*3+0]=rho*w[k]; */
3406 
3407 /* dmom_w_adv_v[k*3+1]=rho*w[k]; */
3408 
3409 /* dmom_w_adv_w[k*3+0]=rho*u[k]; */
3410 /* dmom_w_adv_w[k*3+1]=rho*v[k]; */
3411 /* dmom_w_adv_w[k*3+2]=2.0*rho*w[k]; */
3412 
3413 /* //u momentum diffusion tensor */
3414 /* mom_u_diff_ten[k*9+0] = 2.0*mu; */
3415 /* mom_u_diff_ten[k*9+4] = mu; */
3416 /* mom_u_diff_ten[k*9+8] = mu; */
3417 
3418 /* mom_uv_diff_ten[k*9+3]=mu; */
3419 
3420 /* mom_uw_diff_ten[k*9+6]=mu; */
3421 
3422 /* //v momentum diffusion tensor */
3423 /* mom_v_diff_ten[k*9+0] = mu; */
3424 /* mom_v_diff_ten[k*9+4] = 2.0*mu; */
3425 /* mom_v_diff_ten[k*9+8] = mu; */
3426 
3427 /* mom_vu_diff_ten[k*9+1]=mu; */
3428 
3429 /* mom_vw_diff_ten[k*9+7]=mu; */
3430 
3431 /* //w momentum diffusion tensor */
3432 /* mom_w_diff_ten[k*9+0] = mu; */
3433 /* mom_w_diff_ten[k*9+4] = mu; */
3434 /* mom_w_diff_ten[k*9+8] = 2.0*mu; */
3435 
3436 /* mom_wu_diff_ten[k*9+2]=mu; */
3437 
3438 /* mom_wv_diff_ten[k*9+5]=mu; */
3439 
3440 /* //momentum sources */
3441 /* norm_n = sqrt(n[k*3+0]*n[k*3+0]+n[k*3+1]*n[k*3+1]+n[k*3+2]*n[k*3+2]); */
3442 /* mom_u_source[k] = -rho*g[0] - d_mu*sigma*kappa[k]*n[k*3+0]/(norm_n); */
3443 /* mom_v_source[k] = -rho*g[1] - d_mu*sigma*kappa[k]*n[k*3+1]/(norm_n); */
3444 /* mom_w_source[k] = -rho*g[2] - d_mu*sigma*kappa[k]*n[k*3+2]/(norm_n); */
3445 
3446 
3447 /* //u momentum Hamiltonian (pressure) */
3448 /* mom_u_ham[k] = grad_p[k*3+0]; */
3449 /* dmom_u_ham_grad_p[k*3+0]=1.0; */
3450 
3451 /* //v momentum Hamiltonian (pressure) */
3452 /* mom_v_ham[k] = grad_p[k*3+1]; */
3453 /* dmom_v_ham_grad_p[k*3+1]=1.0; */
3454 
3455 /* //w momentum Hamiltonian (pressure) */
3456 /* mom_w_ham[k] = grad_p[k*3+2]; */
3457 /* dmom_w_ham_grad_p[k*3+2]=1.0; */
3458 
3459  //cek "incompressible" form
3460  //u momentum accumulation
3461  mom_u_acc[k]=u[k];
3462  dmom_u_acc_u[k]=1.0;
3463 
3464  //v momentum accumulation
3465  mom_v_acc[k]=v[k];
3466  dmom_v_acc_v[k]=1.0;
3467 
3468  //w momentum accumulation
3469  mom_w_acc[k]=w[k];
3470  dmom_w_acc_w[k]=1.0;
3471 
3472 
3473  //mass advective flux
3474  mass_adv[k*3+0]=u[k];
3475  mass_adv[k*3+1]=v[k];
3476  mass_adv[k*3+2]=w[k];
3477 
3478  dmass_adv_u[k*3+0]=1.0;
3479  dmass_adv_v[k*3+1]=1.0;
3480  dmass_adv_w[k*3+2]=1.0;
3481 
3482  //u momentum advective flux
3483  mom_u_adv[k*3+0]=u[k]*u[k];
3484  mom_u_adv[k*3+1]=u[k]*v[k];
3485  mom_u_adv[k*3+2]=u[k]*w[k];
3486 
3487  dmom_u_adv_u[k*3+0]=2.0*u[k];
3488  dmom_u_adv_u[k*3+1]=v[k];
3489  dmom_u_adv_u[k*3+2]=w[k];
3490 
3491  dmom_u_adv_v[k*3+1]=u[k];
3492 
3493  dmom_u_adv_w[k*3+2]=u[k];
3494 
3495  //v momentum advective_flux
3496  mom_v_adv[k*3+0]=v[k]*u[k];
3497  mom_v_adv[k*3+1]=v[k]*v[k];
3498  mom_v_adv[k*3+2]=v[k]*w[k];
3499 
3500  dmom_v_adv_u[k*3+0]=v[k];
3501 
3502  dmom_v_adv_w[k*3+2]=v[k];
3503 
3504  dmom_v_adv_v[k*3+0]=u[k];
3505  dmom_v_adv_v[k*3+1]=2.0*v[k];
3506  dmom_v_adv_v[k*3+2]=w[k];
3507 
3508  //w momentum advective_flux
3509  mom_w_adv[k*3+0]=w[k]*u[k];
3510  mom_w_adv[k*3+1]=w[k]*v[k];
3511  mom_w_adv[k*3+2]=w[k]*w[k];
3512 
3513  dmom_w_adv_u[k*3+0]=w[k];
3514 
3515  dmom_w_adv_v[k*3+1]=w[k];
3516 
3517  dmom_w_adv_w[k*3+0]=u[k];
3518  dmom_w_adv_w[k*3+1]=v[k];
3519  dmom_w_adv_w[k*3+2]=2.0*w[k];
3520 
3521  //u momentum diffusion tensor
3522  mom_u_diff_ten[k*9+0] = 2.0*nu;
3523  mom_u_diff_ten[k*9+4] = nu;
3524  mom_u_diff_ten[k*9+8] = nu;
3525 
3526  mom_uv_diff_ten[k*9+3]=nu;
3527 
3528  mom_uw_diff_ten[k*9+6]=nu;
3529 
3530  //v momentum diffusion tensor
3531  mom_v_diff_ten[k*9+0] = nu;
3532  mom_v_diff_ten[k*9+4] = 2.0*nu;
3533  mom_v_diff_ten[k*9+8] = nu;
3534 
3535  mom_vu_diff_ten[k*9+1]=nu;
3536 
3537  mom_vw_diff_ten[k*9+7]=nu;
3538 
3539  //w momentum diffusion tensor
3540  mom_w_diff_ten[k*9+0] = nu;
3541  mom_w_diff_ten[k*9+4] = nu;
3542  mom_w_diff_ten[k*9+8] = 2.0*nu;
3543 
3544  mom_wu_diff_ten[k*9+2]=nu;
3545 
3546  mom_wv_diff_ten[k*9+5]=nu;
3547 
3548  //momentum sources
3549  norm_n = sqrt(n[k*3+0]*n[k*3+0]+n[k*3+1]*n[k*3+1]+n[k*3+2]*n[k*3+2]);
3550  mom_u_source[k] = -g[0] - d_mu*sigma*kappa[k]*n[k*3+0]/(rho*(norm_n+1.0e-8));
3551  mom_v_source[k] = -g[1] - d_mu*sigma*kappa[k]*n[k*3+1]/(rho*(norm_n+1.0e-8));
3552  mom_w_source[k] = -g[2] - d_mu*sigma*kappa[k]*n[k*3+2]/(rho*(norm_n+1.0e-8));
3553 
3554 
3555  //u momentum Hamiltonian (pressure)
3556  mom_u_ham[k] = grad_p[k*3+0]/rho;
3557  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
3558 
3559  //v momentum Hamiltonian (pressure)
3560  mom_v_ham[k] = grad_p[k*3+1]/rho;
3561  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
3562 
3563  //w momentum Hamiltonian (pressure)
3564  mom_w_ham[k] = grad_p[k*3+2]/rho;
3565  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
3566  }
3567 }
3568 void TwophaseNavierStokes_ST_LS_SO_3D_Evaluate_sd(const int nPoints,
3569  const double eps_rho,
3570  const double eps_mu,
3571  const double sigma,
3572  const double rho_0,
3573  const double nu_0,
3574  const double rho_1,
3575  const double nu_1,
3576  const double* g,
3577  const double* phi,
3578  const double* n,
3579  const double* kappa,
3580  const double *p,
3581  const double *grad_p,
3582  const double *u,
3583  const double *v,
3584  const double *w,
3585  double *mom_u_acc,
3586  double *dmom_u_acc_u,
3587  double *mom_v_acc,
3588  double *dmom_v_acc_v,
3589  double *mom_w_acc,
3590  double *dmom_w_acc_w,
3591  double *mass_adv,
3592  double *dmass_adv_u,
3593  double *dmass_adv_v,
3594  double *dmass_adv_w,
3595  double *mom_u_adv,
3596  double *dmom_u_adv_u,
3597  double *dmom_u_adv_v,
3598  double *dmom_u_adv_w,
3599  double *mom_v_adv,
3600  double *dmom_v_adv_u,
3601  double *dmom_v_adv_v,
3602  double *dmom_v_adv_w,
3603  double *mom_w_adv,
3604  double *dmom_w_adv_u,
3605  double *dmom_w_adv_v,
3606  double *dmom_w_adv_w,
3607  double *mom_u_diff_ten,
3608  double *mom_v_diff_ten,
3609  double *mom_w_diff_ten,
3610  double *mom_uv_diff_ten,
3611  double *mom_uw_diff_ten,
3612  double *mom_vu_diff_ten,
3613  double *mom_vw_diff_ten,
3614  double *mom_wu_diff_ten,
3615  double *mom_wv_diff_ten,
3616  double *mom_u_source,
3617  double *mom_v_source,
3618  double *mom_w_source,
3619  double *mom_u_ham,
3620  double *dmom_u_ham_grad_p,
3621  double *mom_v_ham,
3622  double *dmom_v_ham_grad_p,
3623  double *mom_w_ham,
3624  double *dmom_w_ham_grad_p)
3625 {
3626  int k;
3627  double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
3628  for (k=0;k<nPoints;k++)
3629  {
3631  /*H = smoothedHeaviside(eps,phi[k]);*/
3632  H_rho = smoothedHeaviside(eps_rho,phi[k]);
3633  d_rho = smoothedDirac(eps_rho,phi[k]);
3634  H_mu = smoothedHeaviside(eps_mu,phi[k]);
3635  d_mu = smoothedDirac(eps_mu,phi[k]);
3636 
3637  rho = rho_0*(1.0-H_rho)+rho_1*H_rho;
3638  nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
3639  mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
3640 
3641  //u momentum accumulation
3642  mom_u_acc[k]=u[k];
3643  dmom_u_acc_u[k]=1.0;
3644 
3645  //v momentum accumulation
3646  mom_v_acc[k]=v[k];
3647  dmom_v_acc_v[k]=1.0;
3648 
3649  //w momentum accumulation
3650  mom_w_acc[k]=w[k];
3651  dmom_w_acc_w[k]=1.0;
3652 
3653 
3654  //mass advective flux
3655  mass_adv[k*3+0]=u[k];
3656  mass_adv[k*3+1]=v[k];
3657  mass_adv[k*3+2]=w[k];
3658 
3659  dmass_adv_u[k*3+0]=1.0;
3660  dmass_adv_v[k*3+1]=1.0;
3661  dmass_adv_w[k*3+2]=1.0;
3662 
3663  //u momentum advective flux
3664  mom_u_adv[k*3+0]=u[k]*u[k];
3665  mom_u_adv[k*3+1]=u[k]*v[k];
3666  mom_u_adv[k*3+2]=u[k]*w[k];
3667 
3668  dmom_u_adv_u[k*3+0]=2.0*u[k];
3669  dmom_u_adv_u[k*3+1]=v[k];
3670  dmom_u_adv_u[k*3+2]=w[k];
3671 
3672  dmom_u_adv_v[k*3+1]=u[k];
3673 
3674  dmom_u_adv_w[k*3+2]=u[k];
3675 
3676  //v momentum advective_flux
3677  mom_v_adv[k*3+0]=v[k]*u[k];
3678  mom_v_adv[k*3+1]=v[k]*v[k];
3679  mom_v_adv[k*3+2]=v[k]*w[k];
3680 
3681  dmom_v_adv_u[k*3+0]=v[k];
3682 
3683  dmom_v_adv_w[k*3+2]=v[k];
3684 
3685  dmom_v_adv_v[k*3+0]=u[k];
3686  dmom_v_adv_v[k*3+1]=2.0*v[k];
3687  dmom_v_adv_v[k*3+2]=w[k];
3688 
3689  //w momentum advective_flux
3690  mom_w_adv[k*3+0]=w[k]*u[k];
3691  mom_w_adv[k*3+1]=w[k]*v[k];
3692  mom_w_adv[k*3+2]=w[k]*w[k];
3693 
3694  dmom_w_adv_u[k*3+0]=w[k];
3695 
3696  dmom_w_adv_v[k*3+1]=w[k];
3697 
3698  dmom_w_adv_w[k*3+0]=u[k];
3699  dmom_w_adv_w[k*3+1]=v[k];
3700  dmom_w_adv_w[k*3+2]=2.0*w[k];
3701 
3702  //u momentum diffusion tensor
3703  mom_u_diff_ten[k*3+0] = 2.0*nu;
3704  mom_u_diff_ten[k*3+1] = nu;
3705  mom_u_diff_ten[k*3+2] = nu;
3706 
3707  mom_uv_diff_ten[k]=nu;
3708 
3709  mom_uw_diff_ten[k]=nu;
3710 
3711  //v momentum diffusion tensor
3712  mom_v_diff_ten[k*3+0] = nu;
3713  mom_v_diff_ten[k*3+1] = 2.0*nu;
3714  mom_v_diff_ten[k*3+2] = nu;
3715 
3716  mom_vu_diff_ten[k]=nu;
3717 
3718  mom_vw_diff_ten[k]=nu;
3719 
3720  //w momentum diffusion tensor
3721  mom_w_diff_ten[k*3+0] = nu;
3722  mom_w_diff_ten[k*3+1] = nu;
3723  mom_w_diff_ten[k*3+2] = 2.0*nu;
3724 
3725  mom_wu_diff_ten[k]=nu;
3726 
3727  mom_wv_diff_ten[k]=nu;
3728 
3729  //momentum sources
3730  norm_n = sqrt(n[k*3+0]*n[k*3+0]+n[k*3+1]*n[k*3+1]+n[k*3+2]*n[k*3+2]);
3731  mom_u_source[k] = -g[0] - d_mu*sigma*kappa[k]*n[k*3+0]/(rho*(norm_n+1.0e-8));
3732  mom_v_source[k] = -g[1] - d_mu*sigma*kappa[k]*n[k*3+1]/(rho*(norm_n+1.0e-8));
3733  mom_w_source[k] = -g[2] - d_mu*sigma*kappa[k]*n[k*3+2]/(rho*(norm_n+1.0e-8));
3734 
3735 
3736  //u momentum Hamiltonian (pressure)
3737  mom_u_ham[k] = grad_p[k*3+0]/rho;
3738  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
3739 
3740  //v momentum Hamiltonian (pressure)
3741  mom_v_ham[k] = grad_p[k*3+1]/rho;
3742  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
3743 
3744  //w momentum Hamiltonian (pressure)
3745  mom_w_ham[k] = grad_p[k*3+2]/rho;
3746  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
3747  }
3748 }
3749 void ThreephaseNavierStokes_ST_LS_SO_3D_Evaluate(const int nPoints,
3750  const double boundaryPenaltyCoef,
3751  const double volumePenaltyCoef,
3752  const double eps_rho,
3753  const double eps_mu,
3754  const double sigma,
3755  const double rho_0,
3756  const double nu_0,
3757  const double rho_1,
3758  const double nu_1,
3759  const double rho_s,
3760  const double nu_s,
3761  const double* g,
3762  const double* phi,
3763  const double* n,
3764  const double* kappa,
3765  const double* phi_s,
3766  const double* n_s,
3767  const double *p,
3768  const double *grad_p,
3769  const double *u,
3770  const double *v,
3771  const double *w,
3772  double *mom_u_acc,
3773  double *dmom_u_acc_u,
3774  double *mom_v_acc,
3775  double *dmom_v_acc_v,
3776  double *mom_w_acc,
3777  double *dmom_w_acc_w,
3778  double *mass_adv,
3779  double *dmass_adv_u,
3780  double *dmass_adv_v,
3781  double *dmass_adv_w,
3782  double *mom_u_adv,
3783  double *dmom_u_adv_u,
3784  double *dmom_u_adv_v,
3785  double *dmom_u_adv_w,
3786  double *mom_v_adv,
3787  double *dmom_v_adv_u,
3788  double *dmom_v_adv_v,
3789  double *dmom_v_adv_w,
3790  double *mom_w_adv,
3791  double *dmom_w_adv_u,
3792  double *dmom_w_adv_v,
3793  double *dmom_w_adv_w,
3794  double *mom_u_diff_ten,
3795  double *mom_v_diff_ten,
3796  double *mom_w_diff_ten,
3797  double *mom_uv_diff_ten,
3798  double *mom_uw_diff_ten,
3799  double *mom_vu_diff_ten,
3800  double *mom_vw_diff_ten,
3801  double *mom_wu_diff_ten,
3802  double *mom_wv_diff_ten,
3803  double *mom_u_source,
3804  double *dmom_u_source_u,
3805  double *dmom_u_source_v,
3806  double *dmom_u_source_w,
3807  double *mom_v_source,
3808  double *dmom_v_source_u,
3809  double *dmom_v_source_v,
3810  double *dmom_v_source_w,
3811  double *mom_w_source,
3812  double *dmom_w_source_u,
3813  double *dmom_w_source_v,
3814  double *dmom_w_source_w,
3815  double *mom_u_ham,
3816  double *dmom_u_ham_grad_p,
3817  double *mom_v_ham,
3818  double *dmom_v_ham_grad_p,
3819  double *mom_w_ham,
3820  double *dmom_w_ham_grad_p)
3821 {
3822  int k;
3823  double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,
3824  H_rho_s,d_rho_s,H_mu_s,d_mu_s,norm_n,norm_n_s;
3825 
3826  for (k=0;k<nPoints;k++)
3827  {
3829  /*H = smoothedHeaviside(eps,phi[k]);*/
3830  H_rho = smoothedHeaviside(eps_rho,phi[k]);
3831  d_rho = smoothedDirac(eps_rho,phi[k]);
3832  H_mu = smoothedHeaviside(eps_mu,phi[k]);
3833  d_mu = smoothedDirac(eps_mu,phi[k]);
3834 
3835  rho = rho_0*(1.0-H_rho)+rho_1*H_rho;
3836  nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
3837  mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
3838 
3839  H_rho_s = smoothedHeaviside(eps_rho,phi_s[k]);
3840  d_rho_s = smoothedDirac(eps_rho,phi_s[k]);
3841  H_mu_s = smoothedHeaviside(eps_mu,phi_s[k]);
3842  d_mu_s = smoothedDirac(eps_mu,phi_s[k]);
3843 
3844  rho = rho_s*(1.0-H_rho_s)+rho*H_rho_s;
3845  nu = nu_s*(1.0-H_mu_s)+nu*H_mu_s;
3846  mu = rho_s*nu_s*(1.0-H_mu_s)+rho*nu*H_mu_s;
3847 
3848  //u momentum accumulation
3849  mom_u_acc[k]=H_mu_s*rho*u[k];
3850  dmom_u_acc_u[k]=H_mu_s*rho;
3851 
3852  //v momentum accumulation
3853  mom_v_acc[k]=H_mu_s*rho*v[k];
3854  dmom_v_acc_v[k]=H_mu_s*rho;
3855 
3856  //w momentum accumulation
3857  mom_w_acc[k]=H_mu_s*rho*w[k];
3858  dmom_w_acc_w[k]=H_mu_s*rho;
3859 
3860 
3861  //mass advective flux
3862  mass_adv[k*3+0]=u[k];
3863  mass_adv[k*3+1]=v[k];
3864  mass_adv[k*3+2]=w[k];
3865 
3866  dmass_adv_u[k*3+0]=1.0;
3867  dmass_adv_v[k*3+1]=1.0;
3868  dmass_adv_w[k*3+2]=1.0;
3869 
3870  //u momentum advective flux
3871  mom_u_adv[k*3+0]=H_mu_s*rho*u[k]*u[k];
3872  mom_u_adv[k*3+1]=H_mu_s*rho*u[k]*v[k];
3873  mom_u_adv[k*3+2]=H_mu_s*rho*u[k]*w[k];
3874 
3875  dmom_u_adv_u[k*3+0]=2.0*H_mu_s*rho*u[k];
3876  dmom_u_adv_u[k*3+1]=H_mu_s*rho*v[k];
3877  dmom_u_adv_u[k*3+2]=H_mu_s*rho*w[k];
3878 
3879  dmom_u_adv_v[k*3+1]=H_mu_s*rho*u[k];
3880 
3881  dmom_u_adv_w[k*3+2]=H_mu_s*rho*u[k];
3882 
3883  //v momentum advective_flux
3884  mom_v_adv[k*3+0]=H_mu_s*rho*v[k]*u[k];
3885  mom_v_adv[k*3+1]=H_mu_s*rho*v[k]*v[k];
3886  mom_v_adv[k*3+2]=H_mu_s*rho*v[k]*w[k];
3887 
3888  dmom_v_adv_u[k*3+0]=H_mu_s*rho*v[k];
3889 
3890  dmom_v_adv_w[k*3+2]=H_mu_s*rho*v[k];
3891 
3892  dmom_v_adv_v[k*3+0]=H_mu_s*rho*u[k];
3893  dmom_v_adv_v[k*3+1]=2.0*H_mu_s*rho*v[k];
3894  dmom_v_adv_v[k*3+2]=H_mu_s*rho*w[k];
3895 
3896  //w momentum advective_flux
3897  mom_w_adv[k*3+0]=H_mu_s*rho*w[k]*u[k];
3898  mom_w_adv[k*3+1]=H_mu_s*rho*w[k]*v[k];
3899  mom_w_adv[k*3+2]=H_mu_s*rho*w[k]*w[k];
3900 
3901  dmom_w_adv_u[k*3+0]=H_mu_s*rho*w[k];
3902 
3903  dmom_w_adv_v[k*3+1]=H_mu_s*rho*w[k];
3904 
3905  dmom_w_adv_w[k*3+0]=H_mu_s*rho*u[k];
3906  dmom_w_adv_w[k*3+1]=H_mu_s*rho*v[k];
3907  dmom_w_adv_w[k*3+2]=2.0*H_mu_s*rho*w[k];
3908 
3909  //u momentum diffusion tensor
3910  mom_u_diff_ten[k*9+0] = 2.0*H_mu_s*mu;
3911  mom_u_diff_ten[k*9+4] = H_mu_s*mu;
3912  mom_u_diff_ten[k*9+8] = H_mu_s*mu;
3913 
3914  mom_uv_diff_ten[k*9+3]=H_mu_s*mu;
3915 
3916  mom_uw_diff_ten[k*9+6]=H_mu_s*mu;
3917 
3918  //v momentum diffusion tensor
3919  mom_v_diff_ten[k*9+0] = H_mu_s*mu;
3920  mom_v_diff_ten[k*9+4] = 2.0*H_mu_s*mu;
3921  mom_v_diff_ten[k*9+8] = H_mu_s*mu;
3922 
3923  mom_vu_diff_ten[k*9+1]=H_mu_s*mu;
3924 
3925  mom_vw_diff_ten[k*9+7]=H_mu_s*mu;
3926 
3927  //w momentum diffusion tensor
3928  mom_w_diff_ten[k*9+0] = H_mu_s*mu;
3929  mom_w_diff_ten[k*9+4] = H_mu_s*mu;
3930  mom_w_diff_ten[k*9+8] = 2.0*H_mu_s*mu;
3931 
3932  mom_wu_diff_ten[k*9+2]=H_mu_s*mu;
3933 
3934  mom_wv_diff_ten[k*9+5]=H_mu_s*mu;
3935 
3936  //momentum sources
3937  norm_n = sqrt(n[k*3+0]*n[k*3+0]+n[k*3+1]*n[k*3+1]+n[k*3+2]*n[k*3+2]);
3938  mom_u_source[k] = -H_mu_s*rho*g[0] - H_mu_s*d_mu*sigma*kappa[k]*n[k*3+0]/(norm_n)
3939  +boundaryPenaltyCoef*rho*d_mu_s*(u[k] - 0.0)
3940  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(u[k] - 0.0);
3941 
3942  dmom_u_source_u[k] = boundaryPenaltyCoef*rho*d_mu_s
3943  +volumePenaltyCoef*rho*(1.0-H_mu_s);
3944  dmom_u_source_v[k] = 0.0;
3945  dmom_u_source_w[k] = 0.0;
3946 
3947  mom_v_source[k] = -H_mu_s*rho*g[1] - H_mu_s*d_mu*sigma*kappa[k]*n[k*3+1]/(norm_n)
3948  +boundaryPenaltyCoef*rho*d_mu_s*(v[k] - 0.0)
3949  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(v[k] - 0.0);
3950 
3951  dmom_v_source_u[k] = 0.0;
3952  dmom_v_source_v[k] = boundaryPenaltyCoef*rho*d_mu_s
3953  +volumePenaltyCoef*rho*(1.0-H_mu_s);
3954  dmom_v_source_w[k] = 0.0;
3955 
3956 
3957  mom_w_source[k] = -H_mu_s*rho*g[2] - H_mu_s*d_mu*sigma*kappa[k]*n[k*3+2]/(norm_n)
3958  +boundaryPenaltyCoef*rho*d_mu_s*(w[k] - 0.0)
3959  +volumePenaltyCoef*rho*(1.0-H_mu_s)*(w[k] - 0.0);
3960  dmom_w_source_u[k] = 0.0;
3961  dmom_w_source_v[k] = 0.0;
3962  dmom_w_source_w[k] = boundaryPenaltyCoef*rho*d_mu_s
3963  +volumePenaltyCoef*rho*(1.0-H_mu_s);
3964 
3965 
3966 
3967  //u momentum Hamiltonian (pressure)
3968  mom_u_ham[k] = grad_p[k*3+0];
3969  dmom_u_ham_grad_p[k*3+0]=1.0;
3970 
3971  //v momentum Hamiltonian (pressure)
3972  mom_v_ham[k] = grad_p[k*3+1];
3973  dmom_v_ham_grad_p[k*3+1]=1.0;
3974 
3975  //w momentum Hamiltonian (pressure)
3976  mom_w_ham[k] = grad_p[k*3+2];
3977  dmom_w_ham_grad_p[k*3+2]=1.0;
3978  }
3979 }
3980 
3981 void TwophaseNavierStokes_LS_SO_3D_Evaluate(const int nPoints,
3982  const double eps,
3983  const double rho_0,
3984  const double nu_0,
3985  const double rho_1,
3986  const double nu_1,
3987  const double* g,
3988  const double* phi,
3989  const double *p,
3990  const double *grad_p,
3991  const double *u,
3992  const double *v,
3993  const double *w,
3994  double *mom_u_acc,
3995  double *dmom_u_acc_u,
3996  double *mom_v_acc,
3997  double *dmom_v_acc_v,
3998  double *mom_w_acc,
3999  double *dmom_w_acc_w,
4000  double *mass_adv,
4001  double *dmass_adv_u,
4002  double *dmass_adv_v,
4003  double *dmass_adv_w,
4004  double *mom_u_adv,
4005  double *dmom_u_adv_u,
4006  double *dmom_u_adv_v,
4007  double *dmom_u_adv_w,
4008  double *mom_v_adv,
4009  double *dmom_v_adv_u,
4010  double *dmom_v_adv_v,
4011  double *dmom_v_adv_w,
4012  double *mom_w_adv,
4013  double *dmom_w_adv_u,
4014  double *dmom_w_adv_v,
4015  double *dmom_w_adv_w,
4016  double *mom_u_diff_ten,
4017  double *mom_v_diff_ten,
4018  double *mom_w_diff_ten,
4019  double *mom_u_source,
4020  double *mom_v_source,
4021  double *mom_w_source,
4022  double *mom_u_ham,
4023  double *dmom_u_ham_grad_p,
4024  double *mom_v_ham,
4025  double *dmom_v_ham_grad_p,
4026  double *mom_w_ham,
4027  double *dmom_w_ham_grad_p)
4028 {
4029  int k;
4030  double rho,nu,H;
4031  for (k=0;k<nPoints;k++)
4032  {
4033  H = smoothedHeaviside(eps,phi[k]);
4034  rho = rho_0*(1.0-H)+rho_1*H;
4035  nu = nu_0*(1.0-H)+nu_1*H;
4036 
4037  //momentum accumulation
4038  mom_u_acc[k]=u[k];
4039  dmom_u_acc_u[k]=1.0;
4040 
4041  mom_v_acc[k]=v[k];
4042  dmom_v_acc_v[k]=1.0;
4043 
4044  mom_w_acc[k]=w[k];
4045  dmom_w_acc_w[k]=1.0;
4046 
4047  //mass advective flux
4048  mass_adv[k*3+0]=u[k];
4049  mass_adv[k*3+1]=v[k];
4050  mass_adv[k*3+2]=w[k];
4051 
4052  dmass_adv_u[k*3+0]=1.0;
4053  dmass_adv_v[k*3+1]=1.0;
4054  dmass_adv_w[k*3+2]=1.0;
4055 
4056  //u momentum advective flux
4057  mom_u_adv[k*3+0]=u[k]*u[k];
4058  mom_u_adv[k*3+1]=u[k]*v[k];
4059  mom_u_adv[k*3+2]=u[k]*w[k];
4060 
4061  dmom_u_adv_u[k*3+0]=2.0*u[k];
4062  dmom_u_adv_u[k*3+1]=v[k];
4063  dmom_u_adv_u[k*3+2]=w[k];
4064 
4065  dmom_u_adv_v[k*3+1]=u[k];
4066 
4067  dmom_u_adv_w[k*3+2]=u[k];
4068 
4069  //v momentum advective_flux
4070  mom_v_adv[k*3+0]=v[k]*u[k];
4071  mom_v_adv[k*3+1]=v[k]*v[k];
4072  mom_v_adv[k*3+2]=v[k]*w[k];
4073 
4074  dmom_v_adv_u[k*3+0]=v[k];
4075 
4076  dmom_v_adv_v[k*3+0]=u[k];
4077  dmom_v_adv_v[k*3+1]=2.0*v[k];
4078  dmom_v_adv_v[k*3+2]=w[k];
4079 
4080  dmom_v_adv_w[k*3+2]=v[k];
4081 
4082  //w momentum advective_flux
4083  mom_w_adv[k*3+0]=w[k]*u[k];
4084  mom_w_adv[k*3+1]=w[k]*v[k];
4085  mom_w_adv[k*3+2]=w[k]*w[k];
4086 
4087  dmom_w_adv_u[k*3+0]=w[k];
4088 
4089  dmom_w_adv_v[k*3+0]=w[k];
4090 
4091  dmom_w_adv_w[k*3+0]=u[k];
4092  dmom_w_adv_w[k*3+1]=v[k];
4093  dmom_w_adv_w[k*3+2]=2.0*w[k];
4094 
4095  //u momentum diffusion tensor
4096  mom_u_diff_ten[k*9+0] = nu;
4097  mom_u_diff_ten[k*9+4] = nu;
4098  mom_u_diff_ten[k*9+8] = nu;
4099 
4100  //v momentum diffusion tensor
4101  mom_v_diff_ten[k*9+0] = nu;
4102  mom_v_diff_ten[k*9+4] = nu;
4103  mom_v_diff_ten[k*9+8] = nu;
4104 
4105  //w momentum diffusion tensor
4106  mom_w_diff_ten[k*9+0] = nu;
4107  mom_w_diff_ten[k*9+4] = nu;
4108  mom_w_diff_ten[k*9+8] = nu;
4109 
4110  //momentum sources
4111  mom_u_source[k] = -g[0];
4112  mom_v_source[k] = -g[1];
4113  mom_w_source[k] = -g[2];
4114 
4115  //u momentum Hamiltonian (pressure)
4116  mom_u_ham[k] = grad_p[k*3+0]/rho;
4117  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4118 
4119  //v momentum Hamiltonian (pressure)
4120  mom_v_ham[k] = grad_p[k*3+1]/rho;
4121  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4122 
4123  //w momentum Hamiltonian (pressure)
4124  mom_w_ham[k] = grad_p[k*3+2]/rho;
4125  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4126  }
4127 }
4128 
4129 void TwophaseStokes_LS_SO_2D_Evaluate(const int nPoints,
4130  const double eps,
4131  const double rho_0,
4132  const double nu_0,
4133  const double rho_1,
4134  const double nu_1,
4135  const double* g,
4136  const double* phi,
4137  const double *p,
4138  const double *grad_p,
4139  const double *u,
4140  const double *v,
4141  double *mom_u_acc,
4142  double *dmom_u_acc_u,
4143  double *mom_v_acc,
4144  double *dmom_v_acc_v,
4145  double *mass_adv,
4146  double *dmass_adv_u,
4147  double *dmass_adv_v,
4148  double *mom_u_diff_ten,
4149  double *mom_v_diff_ten,
4150  double *mom_u_source,
4151  double *mom_v_source,
4152  double *mom_u_ham,
4153  double *dmom_u_ham_grad_p,
4154  double *mom_v_ham,
4155  double *dmom_v_ham_grad_p)
4156 {
4157  int k;
4158  double rho,nu,H;
4159  for (k=0;k<nPoints;k++)
4160  {
4162  H = smoothedHeaviside(eps,phi[k]);
4163  rho = rho_0*(1.0-H)+rho_1*H;
4164  nu = nu_0*(1.0-H)+nu_1*H;
4165 
4166  //u momentum accumulation
4167  mom_u_acc[k]=u[k];
4168  dmom_u_acc_u[k]=1.0;
4169 
4170  //v momentum accumulation
4171  mom_v_acc[k]=v[k];
4172  dmom_v_acc_v[k]=1.0;
4173 
4174  //mass advective flux
4175  mass_adv[k*2+0]=u[k];
4176  mass_adv[k*2+1]=v[k];
4177 
4178  dmass_adv_u[k*2+0]=1.0;
4179  dmass_adv_v[k*2+1]=1.0;
4180 
4181  //u momentum diffusion tensor
4182  mom_u_diff_ten[k*4+0] = nu;
4183  mom_u_diff_ten[k*4+3] = nu;
4184 
4185  //v momentum diffusion tensor
4186  mom_v_diff_ten[k*4+0] = nu;
4187  mom_v_diff_ten[k*4+3] = nu;
4188 
4189  //momentum sources
4190  mom_u_source[k] = -g[0];
4191  mom_v_source[k] = -g[1];
4192 
4193  //u momentum Hamiltonian (pressure)
4194  mom_u_ham[k] = grad_p[k*2+0]/rho;
4195  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
4196 
4197  //v momentum Hamiltonian (pressure)
4198  mom_v_ham[k] = grad_p[k*2+1]/rho;
4199  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
4200  }
4201 }
4202 
4203 void TwophaseStokes_LS_SO_3D_Evaluate(const int nPoints,
4204  const double eps,
4205  const double rho_0,
4206  const double nu_0,
4207  const double rho_1,
4208  const double nu_1,
4209  const double* g,
4210  const double* phi,
4211  const double *p,
4212  const double *grad_p,
4213  const double *u,
4214  const double *v,
4215  const double *w,
4216  double *mom_u_acc,
4217  double *dmom_u_acc_u,
4218  double *mom_v_acc,
4219  double *dmom_v_acc_v,
4220  double *mom_w_acc,
4221  double *dmom_w_acc_w,
4222  double *mass_adv,
4223  double *dmass_adv_u,
4224  double *dmass_adv_v,
4225  double *dmass_adv_w,
4226  double *mom_u_diff_ten,
4227  double *mom_v_diff_ten,
4228  double *mom_w_diff_ten,
4229  double *mom_u_source,
4230  double *mom_v_source,
4231  double *mom_w_source,
4232  double *mom_u_ham,
4233  double *dmom_u_ham_grad_p,
4234  double *mom_v_ham,
4235  double *dmom_v_ham_grad_p,
4236  double *mom_w_ham,
4237  double *dmom_w_ham_grad_p)
4238 {
4239  int k;
4240  double rho,nu,H;
4241  for (k=0;k<nPoints;k++)
4242  {
4243  H = smoothedHeaviside(eps,phi[k]);
4244  rho = rho_0*(1.0-H)+rho_1*H;
4245  nu = nu_0*(1.0-H)+nu_1*H;
4246 
4247  //momentum accumulation
4248  mom_u_acc[k]=u[k];
4249  dmom_u_acc_u[k]=1.0;
4250 
4251  mom_v_acc[k]=v[k];
4252  dmom_v_acc_v[k]=1.0;
4253 
4254  mom_w_acc[k]=w[k];
4255  dmom_w_acc_w[k]=1.0;
4256 
4257  //mass advective flux
4258  mass_adv[k*3+0]=u[k];
4259  mass_adv[k*3+1]=v[k];
4260  mass_adv[k*3+2]=w[k];
4261 
4262  dmass_adv_u[k*3+0]=1.0;
4263  dmass_adv_v[k*3+1]=1.0;
4264  dmass_adv_w[k*3+2]=1.0;
4265 
4266  //u momentum diffusion tensor
4267  mom_u_diff_ten[k*9+0] = nu;
4268  mom_u_diff_ten[k*9+4] = nu;
4269  mom_u_diff_ten[k*9+8] = nu;
4270 
4271  //v momentum diffusion tensor
4272  mom_v_diff_ten[k*9+0] = nu;
4273  mom_v_diff_ten[k*9+4] = nu;
4274  mom_v_diff_ten[k*9+8] = nu;
4275 
4276  //w momentum diffusion tensor
4277  mom_w_diff_ten[k*9+0] = nu;
4278  mom_w_diff_ten[k*9+4] = nu;
4279  mom_w_diff_ten[k*9+8] = nu;
4280 
4281  //momentum sources
4282  mom_u_source[k] = -g[0];
4283  mom_v_source[k] = -g[1];
4284  mom_w_source[k] = -g[2];
4285 
4286  //u momentum Hamiltonian (pressure)
4287  mom_u_ham[k] = grad_p[k*3+0]/rho;
4288  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4289 
4290  //v momentum Hamiltonian (pressure)
4291  mom_v_ham[k] = grad_p[k*3+1]/rho;
4292  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4293 
4294  //w momentum Hamiltonian (pressure)
4295  mom_w_ham[k] = grad_p[k*3+2]/rho;
4296  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4297  }
4298 }
4299 
4300 void TwophaseNavierStokes_VOF_SO_2D_Evaluate(const int nPoints,
4301  const double eps,
4302  const double rho_0,
4303  const double nu_0,
4304  const double rho_1,
4305  const double nu_1,
4306  const double* g,
4307  const double* vof,
4308  const double *p,
4309  const double *grad_p,
4310  const double *u,
4311  const double *v,
4312  double *mom_u_acc,
4313  double *dmom_u_acc_u,
4314  double *mom_v_acc,
4315  double *dmom_v_acc_v,
4316  double *mass_adv,
4317  double *dmass_adv_u,
4318  double *dmass_adv_v,
4319  double *mom_u_adv,
4320  double *dmom_u_adv_u,
4321  double *dmom_u_adv_v,
4322  double *mom_v_adv,
4323  double *dmom_v_adv_u,
4324  double *dmom_v_adv_v,
4325  double *mom_u_diff_ten,
4326  double *mom_v_diff_ten,
4327  double *mom_u_source,
4328  double *mom_v_source,
4329  double *mom_u_ham,
4330  double *dmom_u_ham_grad_p,
4331  double *mom_v_ham,
4332  double *dmom_v_ham_grad_p)
4333 {
4334  int k;
4335  double rho,nu,H;
4336  for (k=0;k<nPoints;k++)
4337  {
4339  H = fmax(0.0,fmin(1.0,vof[k]));
4340  rho = rho_0*(1.0-H)+rho_1*H;
4341  nu = nu_0*(1.0-H)+nu_1*H;
4342 
4343  //u momentum accumulation
4344  mom_u_acc[k]=u[k];
4345  dmom_u_acc_u[k]=1.0;
4346 
4347  //v momentum accumulation
4348  mom_v_acc[k]=v[k];
4349  dmom_v_acc_v[k]=1.0;
4350 
4351  //mass advective flux
4352  mass_adv[k*2+0]=u[k];
4353  mass_adv[k*2+1]=v[k];
4354 
4355  dmass_adv_u[k*2+0]=1.0;
4356  dmass_adv_v[k*2+1]=1.0;
4357 
4358  //u momentum advective flux
4359  mom_u_adv[k*2+0]=u[k]*u[k];
4360  mom_u_adv[k*2+1]=u[k]*v[k];
4361 
4362  dmom_u_adv_u[k*2+0]=2.0*u[k];
4363  dmom_u_adv_u[k*2+1]=v[k];
4364 
4365  dmom_u_adv_v[k*2+1]=u[k];
4366 
4367  //v momentum advective_flux
4368  mom_v_adv[k*2+0]=v[k]*u[k];
4369  mom_v_adv[k*2+1]=v[k]*v[k];
4370 
4371  dmom_v_adv_u[k*2+0]=v[k];
4372 
4373  dmom_v_adv_v[k*2+0]=u[k];
4374  dmom_v_adv_v[k*2+1]=2.0*v[k];
4375 
4376  //u momentum diffusion tensor
4377  mom_u_diff_ten[k*4+0] = nu;
4378  mom_u_diff_ten[k*4+3] = nu;
4379 
4380  //v momentum diffusion tensor
4381  mom_v_diff_ten[k*4+0] = nu;
4382  mom_v_diff_ten[k*4+3] = nu;
4383 
4384  //momentum sources
4385  mom_u_source[k] = -g[0];
4386  mom_v_source[k] = -g[1];
4387 
4388  //u momentum Hamiltonian (pressure)
4389  mom_u_ham[k] = grad_p[k*2+0]/rho;
4390  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
4391 
4392  //v momentum Hamiltonian (pressure)
4393  mom_v_ham[k] = grad_p[k*2+1]/rho;
4394  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
4395  }
4396 }
4397 
4398 void TwophaseNavierStokes_VOF_SO_3D_Evaluate(const int nPoints,
4399  const double eps,
4400  const double rho_0,
4401  const double nu_0,
4402  const double rho_1,
4403  const double nu_1,
4404  const double* g,
4405  const double* vof,
4406  const double *p,
4407  const double *grad_p,
4408  const double *u,
4409  const double *v,
4410  const double *w,
4411  double *mom_u_acc,
4412  double *dmom_u_acc_u,
4413  double *mom_v_acc,
4414  double *dmom_v_acc_v,
4415  double *mom_w_acc,
4416  double *dmom_w_acc_w,
4417  double *mass_adv,
4418  double *dmass_adv_u,
4419  double *dmass_adv_v,
4420  double *dmass_adv_w,
4421  double *mom_u_adv,
4422  double *dmom_u_adv_u,
4423  double *dmom_u_adv_v,
4424  double *dmom_u_adv_w,
4425  double *mom_v_adv,
4426  double *dmom_v_adv_u,
4427  double *dmom_v_adv_v,
4428  double *dmom_v_adv_w,
4429  double *mom_w_adv,
4430  double *dmom_w_adv_u,
4431  double *dmom_w_adv_v,
4432  double *dmom_w_adv_w,
4433  double *mom_u_diff_ten,
4434  double *mom_v_diff_ten,
4435  double *mom_w_diff_ten,
4436  double *mom_u_source,
4437  double *mom_v_source,
4438  double *mom_w_source,
4439  double *mom_u_ham,
4440  double *dmom_u_ham_grad_p,
4441  double *mom_v_ham,
4442  double *dmom_v_ham_grad_p,
4443  double *mom_w_ham,
4444  double *dmom_w_ham_grad_p)
4445 {
4446  int k;
4447  double rho,nu,H;
4448  for (k=0;k<nPoints;k++)
4449  {
4450  H = fmax(0.0,fmin(1.0,vof[k]));
4451  rho = rho_0*(1.0-H)+rho_1*H;
4452  nu = nu_0*(1.0-H)+nu_1*H;
4453 
4454  //momentum accumulation
4455  mom_u_acc[k]=u[k];
4456  dmom_u_acc_u[k]=1.0;
4457 
4458  mom_v_acc[k]=v[k];
4459  dmom_v_acc_v[k]=1.0;
4460 
4461  mom_w_acc[k]=w[k];
4462  dmom_w_acc_w[k]=1.0;
4463 
4464  //mass advective flux
4465  mass_adv[k*3+0]=u[k];
4466  mass_adv[k*3+1]=v[k];
4467  mass_adv[k*3+2]=w[k];
4468 
4469  dmass_adv_u[k*3+0]=1.0;
4470  dmass_adv_v[k*3+1]=1.0;
4471  dmass_adv_w[k*3+2]=1.0;
4472 
4473  //u momentum advective flux
4474  mom_u_adv[k*3+0]=u[k]*u[k];
4475  mom_u_adv[k*3+1]=u[k]*v[k];
4476  mom_u_adv[k*3+2]=u[k]*w[k];
4477 
4478  dmom_u_adv_u[k*3+0]=2.0*u[k];
4479  dmom_u_adv_u[k*3+1]=v[k];
4480  dmom_u_adv_u[k*3+2]=w[k];
4481 
4482  dmom_u_adv_v[k*3+1]=u[k];
4483 
4484  dmom_u_adv_w[k*3+2]=u[k];
4485 
4486  //v momentum advective_flux
4487  mom_v_adv[k*3+0]=v[k]*u[k];
4488  mom_v_adv[k*3+1]=v[k]*v[k];
4489  mom_v_adv[k*3+2]=v[k]*w[k];
4490 
4491  dmom_v_adv_u[k*3+0]=v[k];
4492 
4493  dmom_v_adv_v[k*3+0]=u[k];
4494  dmom_v_adv_v[k*3+1]=2.0*v[k];
4495  dmom_v_adv_v[k*3+2]=w[k];
4496 
4497  dmom_v_adv_w[k*3+2]=v[k];
4498 
4499  //w momentum advective_flux
4500  mom_w_adv[k*3+0]=w[k]*u[k];
4501  mom_w_adv[k*3+1]=w[k]*v[k];
4502  mom_w_adv[k*3+2]=w[k]*w[k];
4503 
4504  dmom_w_adv_u[k*3+0]=w[k];
4505 
4506  dmom_w_adv_v[k*3+0]=w[k];
4507 
4508  dmom_w_adv_w[k*3+0]=u[k];
4509  dmom_w_adv_w[k*3+1]=v[k];
4510  dmom_w_adv_w[k*3+2]=2.0*w[k];
4511 
4512  //u momentum diffusion tensor
4513  mom_u_diff_ten[k*9+0] = nu;
4514  mom_u_diff_ten[k*9+4] = nu;
4515  mom_u_diff_ten[k*9+8] = nu;
4516 
4517  //v momentum diffusion tensor
4518  mom_v_diff_ten[k*9+0] = nu;
4519  mom_v_diff_ten[k*9+4] = nu;
4520  mom_v_diff_ten[k*9+8] = nu;
4521 
4522  //w momentum diffusion tensor
4523  mom_w_diff_ten[k*9+0] = nu;
4524  mom_w_diff_ten[k*9+4] = nu;
4525  mom_w_diff_ten[k*9+8] = nu;
4526 
4527  //momentum sources
4528  mom_u_source[k] = -g[0];
4529  mom_v_source[k] = -g[1];
4530  mom_w_source[k] = -g[2];
4531 
4532  //u momentum Hamiltonian (pressure)
4533  mom_u_ham[k] = grad_p[k*3+0]/rho;
4534  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4535 
4536  //v momentum Hamiltonian (pressure)
4537  mom_v_ham[k] = grad_p[k*3+1]/rho;
4538  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4539 
4540  //w momentum Hamiltonian (pressure)
4541  mom_w_ham[k] = grad_p[k*3+2]/rho;
4542  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4543  }
4544 }
4545 
4546 void TwophaseStokes_VOF_SO_2D_Evaluate(const int nPoints,
4547  const double eps,
4548  const double rho_0,
4549  const double nu_0,
4550  const double rho_1,
4551  const double nu_1,
4552  const double* g,
4553  const double* vof,
4554  const double *p,
4555  const double *grad_p,
4556  const double *u,
4557  const double *v,
4558  double *mom_u_acc,
4559  double *dmom_u_acc_u,
4560  double *mom_v_acc,
4561  double *dmom_v_acc_v,
4562  double *mass_adv,
4563  double *dmass_adv_u,
4564  double *dmass_adv_v,
4565  double *mom_u_diff_ten,
4566  double *mom_v_diff_ten,
4567  double *mom_u_source,
4568  double *mom_v_source,
4569  double *mom_u_ham,
4570  double *dmom_u_ham_grad_p,
4571  double *mom_v_ham,
4572  double *dmom_v_ham_grad_p)
4573 {
4574  int k;
4575  double rho,nu,H;
4576  for (k=0;k<nPoints;k++)
4577  {
4579  H = fmax(0.0,fmin(1.0,vof[k]));
4580  rho = rho_0*(1.0-H)+rho_1*H;
4581  nu = nu_0*(1.0-H)+nu_1*H;
4582 
4583  //u momentum accumulation
4584  mom_u_acc[k]=u[k];
4585  dmom_u_acc_u[k]=1.0;
4586 
4587  //v momentum accumulation
4588  mom_v_acc[k]=v[k];
4589  dmom_v_acc_v[k]=1.0;
4590 
4591  //mass advective flux
4592  mass_adv[k*2+0]=u[k];
4593  mass_adv[k*2+1]=v[k];
4594 
4595  dmass_adv_u[k*2+0]=1.0;
4596  dmass_adv_v[k*2+1]=1.0;
4597 
4598  //u momentum diffusion tensor
4599  mom_u_diff_ten[k*4+0] = nu;
4600  mom_u_diff_ten[k*4+3] = nu;
4601 
4602  //v momentum diffusion tensor
4603  mom_v_diff_ten[k*4+0] = nu;
4604  mom_v_diff_ten[k*4+3] = nu;
4605 
4606  //momentum sources
4607  mom_u_source[k] = -g[0];
4608  mom_v_source[k] = -g[1];
4609 
4610  //u momentum Hamiltonian (pressure)
4611  mom_u_ham[k] = grad_p[k*2+0]/rho;
4612  dmom_u_ham_grad_p[k*2+0]=1.0/rho;
4613 
4614  //v momentum Hamiltonian (pressure)
4615  mom_v_ham[k] = grad_p[k*2+1]/rho;
4616  dmom_v_ham_grad_p[k*2+1]=1.0/rho;
4617  }
4618 }
4619 
4620 void TwophaseStokes_VOF_SO_3D_Evaluate(const int nPoints,
4621  const double eps,
4622  const double rho_0,
4623  const double nu_0,
4624  const double rho_1,
4625  const double nu_1,
4626  const double* g,
4627  const double* vof,
4628  const double *p,
4629  const double *grad_p,
4630  const double *u,
4631  const double *v,
4632  const double *w,
4633  double *mom_u_acc,
4634  double *dmom_u_acc_u,
4635  double *mom_v_acc,
4636  double *dmom_v_acc_v,
4637  double *mom_w_acc,
4638  double *dmom_w_acc_w,
4639  double *mass_adv,
4640  double *dmass_adv_u,
4641  double *dmass_adv_v,
4642  double *dmass_adv_w,
4643  double *mom_u_diff_ten,
4644  double *mom_v_diff_ten,
4645  double *mom_w_diff_ten,
4646  double *mom_u_source,
4647  double *mom_v_source,
4648  double *mom_w_source,
4649  double *mom_u_ham,
4650  double *dmom_u_ham_grad_p,
4651  double *mom_v_ham,
4652  double *dmom_v_ham_grad_p,
4653  double *mom_w_ham,
4654  double *dmom_w_ham_grad_p)
4655 {
4656  int k;
4657  double rho,nu,H;
4658  for (k=0;k<nPoints;k++)
4659  {
4660  H = fmax(0.0,fmin(1.0,vof[k]));
4661  rho = rho_0*(1.0-H)+rho_1*H;
4662  nu = nu_0*(1.0-H)+nu_1*H;
4663 
4664  //momentum accumulation
4665  mom_u_acc[k]=u[k];
4666  dmom_u_acc_u[k]=1.0;
4667 
4668  mom_v_acc[k]=v[k];
4669  dmom_v_acc_v[k]=1.0;
4670 
4671  mom_w_acc[k]=w[k];
4672  dmom_w_acc_w[k]=1.0;
4673 
4674  //mass advective flux
4675  mass_adv[k*3+0]=u[k];
4676  mass_adv[k*3+1]=v[k];
4677  mass_adv[k*3+2]=w[k];
4678 
4679  dmass_adv_u[k*3+0]=1.0;
4680  dmass_adv_v[k*3+1]=1.0;
4681  dmass_adv_w[k*3+2]=1.0;
4682 
4683  //u momentum diffusion tensor
4684  mom_u_diff_ten[k*9+0] = nu;
4685  mom_u_diff_ten[k*9+4] = nu;
4686  mom_u_diff_ten[k*9+8] = nu;
4687 
4688  //v momentum diffusion tensor
4689  mom_v_diff_ten[k*9+0] = nu;
4690  mom_v_diff_ten[k*9+4] = nu;
4691  mom_v_diff_ten[k*9+8] = nu;
4692 
4693  //w momentum diffusion tensor
4694  mom_w_diff_ten[k*9+0] = nu;
4695  mom_w_diff_ten[k*9+4] = nu;
4696  mom_w_diff_ten[k*9+8] = nu;
4697 
4698  //momentum sources
4699  mom_u_source[k] = -g[0];
4700  mom_v_source[k] = -g[1];
4701  mom_w_source[k] = -g[2];
4702 
4703  //u momentum Hamiltonian (pressure)
4704  mom_u_ham[k] = grad_p[k*3+0]/rho;
4705  dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4706 
4707  //v momentum Hamiltonian (pressure)
4708  mom_v_ham[k] = grad_p[k*3+1]/rho;
4709  dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4710 
4711  //w momentum Hamiltonian (pressure)
4712  mom_w_ham[k] = grad_p[k*3+2]/rho;
4713  dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4714  }
4715 }
4716 
4717 void unitSquareVortexEvaluate(const int nPoints,
4718  const int nSpace,
4719  double t,
4720  const double *x,
4721  const double *u,
4722  double *m,
4723  double *dm,
4724  double *f,
4725  double *df)
4726 {
4727  double vx, vy, xk, yk;
4728  int k;
4729  double one8 = 1.0/8.0;
4730  for (k=0; k < nPoints; k++)
4731  {
4732  m[k] = u[k];
4733  dm[k] = 1.0;
4734  xk = x[k*3]; yk = x[k*3+1];
4735  vx = cos(M_PI*one8*t)*sin(2.0*M_PI*yk)*sin(M_PI*xk)*sin(M_PI*xk);
4736  vy =-cos(M_PI*one8*t)*sin(2.0*M_PI*xk)*sin(M_PI*yk)*sin(M_PI*yk);
4737  f[k*nSpace] = vx*u[k];
4738  f[k*nSpace+1] = vy*u[k];
4739  df[k*nSpace] = vx;
4740  df[k*nSpace+1] = vy;
4741  }
4742 }
4743 
4744 /*for HJ testing*/
4745 void constantVelocityLevelSetEvaluate(const int nPoints,
4746  const int nSpace,
4747  const double *b,
4748  const double *x,
4749  const double *u,
4750  const double *gradu,
4751  double *m,
4752  double *dm,
4753  double *f,
4754  double *df,
4755  double *H,
4756  double *dH)
4757 {
4758  double bdotgrad;
4759  int k,id;
4760  for (k=0; k < nPoints; k++)
4761  {
4762  m[k] = u[k];
4763  dm[k] = 1.0;
4764  bdotgrad=0.0;
4765  for (id=0; id < nSpace; id++)
4766  {
4767  f[k*nSpace+id] = 0.0; /*just for now since helps get all the right quad terms*/
4768  df[k*nSpace+id]= 0.0; /*just for now since helps get all the right quad terms*/
4769 
4770  bdotgrad+= gradu[k*nSpace+id]*b[id];
4771  dH[k*nSpace+id]=b[id];
4772  }
4773  H[k] = bdotgrad;
4774  }
4775 }
4776 void constantNormalVelocityLevelSetEvaluate(const int nPoints,
4777  const int nSpace,
4778  double b,
4779  const double *x,
4780  const double *u,
4781  const double *gradu,
4782  double *m,
4783  double *dm,
4784  double *f,
4785  double *df,
4786  double *H,
4787  double *dH)
4788 {
4789  double normgradu;
4790  int k,id;
4791  for (k=0; k < nPoints; k++)
4792  {
4793  m[k] = u[k];
4794  dm[k] = 1.0;
4795  normgradu=0.0;
4796  for (id=0; id < nSpace; id++)
4797  {
4798  f[k*nSpace+id] = 0.0; /*just for now since helps get all the right quad terms*/
4799  df[k*nSpace+id]= 0.0; /*just for now since helps get all the right quad terms*/
4800 
4801  normgradu+= gradu[k*nSpace+id]*gradu[k*nSpace+id];
4802 
4803  }
4804  normgradu = sqrt(normgradu);
4805  H[k] = b*normgradu;
4806  for (id=0; id < nSpace; id++)
4807  {
4808  dH[k*nSpace+id]=gradu[k*nSpace+id]/(normgradu+1.0e-8);
4809  }
4810  }
4811 }
4812 void unitSquareVortexLevelSetEvaluate(const int nPoints,
4813  const int nSpace,
4814  double t,
4815  const double *x,
4816  const double *u,
4817  const double *gradu,
4818  double *m,
4819  double *dm,
4820  double *f,
4821  double *df,
4822  double *H,
4823  double *dH)
4824 {
4825  double v[3], xk, yk, vdotgrad, one8;
4826  int k,id;
4827  v[2] = 0.0;
4828  one8 = 1.0/8.0;
4829  for (k=0; k < nPoints; k++)
4830  {
4831  m[k] = u[k];
4832  dm[k] = 1.0;
4833  xk = x[k*3]; yk = x[k*3+1];
4834  v[0] = cos(M_PI*one8*t)*sin(2.0*M_PI*yk)*sin(M_PI*xk)*sin(M_PI*xk);
4835  v[1] =-cos(M_PI*one8*t)*sin(2.0*M_PI*xk)*sin(M_PI*yk)*sin(M_PI*yk);
4836  vdotgrad=0.0;
4837  for (id=0; id < nSpace; id++)
4838  {
4839  f[k*nSpace+id] = 0.0; /*just for now since helps get all the right quad terms*/
4840  df[k*nSpace+id]= 0.0; /*just for now since helps get all the right quad terms*/
4841 
4842  vdotgrad+= gradu[k*nSpace+id]*v[id];
4843  dH[k*nSpace+id]=v[id];
4844  }
4845  H[k] = vdotgrad;
4846  }
4847 }
4848 void unitSquareRotationLevelSetEvaluate(const int nPoints,
4849  const int nSpace,
4850  double t,
4851  const double *x,
4852  const double *u,
4853  const double *gradu,
4854  double *m,
4855  double *dm,
4856  double *f,
4857  double *df,
4858  double *H,
4859  double *dH)
4860 {
4861  double v[3], xk, yk, vdotgrad, one8;
4862  int k,id;
4863  v[2] = 0.0;
4864  one8 = 1.0/8.0;
4865  for (k=0; k < nPoints; k++)
4866  {
4867  m[k] = u[k];
4868  dm[k] = 1.0;
4869  xk = x[k*3]; yk = x[k*3+1];
4870  v[0] = 2.0*M_PI*(x[k*3+1] - 0.5);
4871  v[1] =2.0*M_PI*(0.5 - x[k*3]);
4872  vdotgrad=0.0;
4873  for (id=0; id < nSpace; id++)
4874  {
4875  f[k*nSpace+id] = 0.0; /*just for now since helps get all the right quad terms*/
4876  df[k*nSpace+id]= 0.0; /*just for now since helps get all the right quad terms*/
4877 
4878  vdotgrad+= gradu[k*nSpace+id]*v[id];
4879  dH[k*nSpace+id]=v[id];
4880  }
4881  H[k] = vdotgrad;
4882  }
4883 }
4884 
4885 void HJBurgersEvaluate(const int nPoints,
4886  const int nSpace,
4887  const double offset,
4888  const double *u,
4889  const double *gradu,
4890  double *m,
4891  double *dm,
4892  double *H,
4893  double *dH)
4894 {
4895  int k,id,jd;
4896  double tmp;
4897  for (k=0; k < nPoints; k++)
4898  {
4899  m[k] = u[k];
4900  dm[k] = 1.0;
4901  tmp = offset;
4902  for (id=0; id < nSpace; id++)
4903  {
4904  tmp += gradu[k*nSpace+id];
4905  dH[k*nSpace+id]= offset;
4906  for (jd=0; jd < nSpace; jd++)
4907  dH[k*nSpace+id] += gradu[k*nSpace+jd];
4908  }
4909  H[k] = tmp*tmp*0.5;
4910  }
4911 }
4912 
4913 /* /\** Coefficients for the mass conservative head-based Richards' equation using Mualem-Van Genuchten. */
4914 /* *\/ */
4915 /* void conservativeHeadRichardsMualemVanGenuchtenEvaluate(const int nPoints, */
4916 /* const int nSpace, */
4917 /* const double rho, */
4918 /* const double* gravity, */
4919 /* const double* alpha, */
4920 /* const double* n, */
4921 /* const double* m, */
4922 /* const double* thetaS, */
4923 /* const double* thetaR, */
4924 /* const double* thetaSR, */
4925 /* const double* KWs, */
4926 /* double *u, */
4927 /* double *mass, */
4928 /* double *dmass, */
4929 /* double *f, */
4930 /* double *df, */
4931 /* double *a, */
4932 /* double *da, */
4933 /* double *phi, */
4934 /* double *dphi) */
4935 /* { */
4936 /* int k,I,J; */
4937 /* const int nSpace2=nSpace*nSpace; */
4938 /* register double psiC,alphaPsiC,alphaPsiC_n,alphaPsiC_nM1,onePlus_alphaPsiC_n,sBar,sBarByOnePlus_alphaPsiC_n,sqrt_sBar,sqrt_1minusSbar,thetaW,DsBar_DpC,DthetaW_DpC,vBar,uBar,krW,krN, */
4939 /* alphaPsiC_nM2,sBarBy_onePlus_alphaPsiC_n_2,DDsBar_DDpC,DDthetaW_DDpC,DkrW_DpC,rho2=rho*rho; */
4940 /* double KW[9],KN[9],DKW_DpC[9]; */
4941 /* for (k=0;k<nPoints;k++) */
4942 /* { */
4943 /* psiC = -u[k]; */
4944 /* if (psiC > 0.0) */
4945 /* { */
4946 /* alphaPsiC = alpha[k]*psiC; */
4947 /* alphaPsiC_n = pow(alphaPsiC,n[k]); */
4948 /* alphaPsiC_nM1 = alphaPsiC_n/alphaPsiC; */
4949 /* alphaPsiC_nM2 = alphaPsiC_nM1/alphaPsiC; */
4950 /* onePlus_alphaPsiC_n = 1.0 + alphaPsiC_n; */
4951 /* sBar = pow(onePlus_alphaPsiC_n,-m[k]); */
4952 /* sBarByOnePlus_alphaPsiC_n = sBar/onePlus_alphaPsiC_n; */
4953 /* sqrt_sBar = sqrt(sBar); */
4954 /* sqrt_1minusSbar = sqrt(1.0 - sBar); */
4955 /* thetaW = thetaSR[k]*sBar + thetaR[k]; */
4956 /* DsBar_DpC = -alpha[k]*(n[k]-1.0)*alphaPsiC_nM1 */
4957 /* *sBarByOnePlus_alphaPsiC_n; */
4958 /* DthetaW_DpC = thetaSR[k] * DsBar_DpC; */
4959 /* vBar = 1.0-alphaPsiC_nM1*sBar; */
4960 /* uBar = alphaPsiC_nM1*sBar; */
4961 /* krW = sqrt_sBar*vBar*vBar; */
4962 /* DkrW_DpC = (0.5/sqrt_sBar)*DsBar_DpC*vBar*vBar */
4963 /* - */
4964 /* 2.0*sqrt_sBar*vBar* */
4965 /* (alpha[k]*(n[k]-1.0)*alphaPsiC_nM2*sBar */
4966 /* + alphaPsiC_nM1 * DsBar_DpC); */
4967 /* for (I=0;I<nSpace;I++) */
4968 /* for(J=0;J<nSpace;J++) */
4969 /* { */
4970 /* KW[I*nSpace + J] = KWs[k*nSpace2 + I*nSpace + J]*krW; */
4971 /* DKW_DpC[I*nSpace + J] = KWs[k*nSpace2 + I*nSpace + J]*DkrW_DpC; */
4972 /* } */
4973 /* } */
4974 /* else */
4975 /* { */
4976 /* sBar = 1.0; */
4977 /* thetaW = thetaS[k]; */
4978 /* DsBar_DpC = 0.0; */
4979 /* DthetaW_DpC = 0.0; */
4980 /* krW = 1.0; */
4981 /* for (I=0;I<nSpace;I++) */
4982 /* for(J=0;J<nSpace;J++) */
4983 /* { */
4984 /* KW[I*nSpace + J] = KWs[k*nSpace2 + I*nSpace + J]; */
4985 /* DKW_DpC[I*nSpace + J] = 0.0; */
4986 /* } */
4987 /* } */
4988 /* mass[k] = rho*thetaW; */
4989 /* dmass[k] = -rho*DthetaW_DpC; */
4990 /* for (I=0;I<nSpace;I++) */
4991 /* { */
4992 /* f[k*nSpace+I] = 0.0; */
4993 /* df[k*nSpace+I] = 0.0; */
4994 /* for (J=0;J<nSpace;J++) */
4995 /* { */
4996 /* f[k*nSpace+I] += rho2*KW[I*nSpace + J]*gravity[J]; */
4997 /* df[k*nSpace+I] -= rho2*DKW_DpC[I*nSpace + J]*gravity[J]; */
4998 /* a[k*nSpace2+I*nSpace+J] = rho*KW[I*nSpace + J]; */
4999 /* da[k*nSpace2+I*nSpace+J] = -rho*DKW_DpC[I*nSpace + J]; */
5000 /* } */
5001 /* } */
5002 /* } */
5003 /* } */
5004 
5007 void conservativeHeadRichardsMualemVanGenuchtenHomEvaluate(const int nPoints,
5008  const int nSpace,
5009  const double rho,
5010  const double beta,
5011  const double* gravity,
5012  const double* x,
5013  const double alpha,
5014  const double n,
5015  const double m,
5016  const double thetaR,
5017  const double thetaSR,
5018  const double KWs,
5019  double *u,
5020  double *mass,
5021  double *dmass,
5022  double *f,
5023  double *df,
5024  double *a,
5025  double *da,
5026  double *phi,
5027  double *dphi)
5028 {
5029  int k,I;
5030  const int nSpace2=nSpace*nSpace;
5031  register double psiC,
5032  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5033  onePlus_pcBar_n,
5034  sBar,sqrt_sBar,DsBar_DpsiC,
5035  thetaW,DthetaW_DpsiC,
5036  vBar,vBar2,DvBar_DpsiC,
5037  KW,DKW_DpsiC,
5038  rho2=rho*rho,
5039  thetaS=thetaR+thetaSR,
5040  rhom,drhom;
5041  for (k=0;k<nPoints;k++)
5042  {
5043  psiC = -u[k];
5044  if (psiC > 0.0)
5045  {
5046  pcBar = alpha*psiC;
5047  pcBar_nM2 = pow(pcBar,n-2);
5048  pcBar_nM1 = pcBar_nM2*pcBar;
5049  pcBar_n = pcBar_nM1*pcBar;
5050  onePlus_pcBar_n = 1.0 + pcBar_n;
5051 
5052  sBar = pow(onePlus_pcBar_n,-m);
5053  /* using -mn = 1-n */
5054  DsBar_DpsiC = alpha*(1.0-n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5055 
5056  vBar = 1.0-pcBar_nM1*sBar;
5057  vBar2 = vBar*vBar;
5058  DvBar_DpsiC = -alpha*(n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5059 
5060  thetaW = thetaSR*sBar + thetaR;
5061  DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5062 
5063  sqrt_sBar = sqrt(sBar);
5064  KW= KWs*sqrt_sBar*vBar2;
5065  DKW_DpsiC= KWs*
5066  ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5067  +
5068  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5069  }
5070  else
5071  {
5072  thetaW = thetaS;
5073  DthetaW_DpsiC = 0.0;
5074  KW = KWs;
5075  DKW_DpsiC = 0.0;
5076  }
5077  //slight compressibility
5078  rhom = rho*exp(beta*u[k]);
5079  drhom = beta*rhom;
5080 
5081  mass[k] = rhom*thetaW;
5082  dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
5083  for (I=0;I<nSpace;I++)
5084  {
5085  f[k*nSpace+I] = rho2*KW*gravity[I];
5086  df[k*nSpace+I] = -rho2*DKW_DpsiC*gravity[I];
5087 
5088  a[k*nSpace2+I*nSpace+I] = rho*KW;
5089  da[k*nSpace2+I*nSpace+I] = -rho*DKW_DpsiC;
5090  }
5091  }
5092 }
5093 /* mwf begin unnecessary RE additions */
5094 /* TODO: figure out how to handle dmass term to get right jacobian*/
5095 void conservativeHeadRichardsL2projMualemVanGenuchtenHomEvaluate(const int nSimplices,
5096  const int nPointsPerSimplex,
5097  const int nSpace,
5098  const double rho,
5099  const double* gravity,
5100  const double alpha,
5101  const double n,
5102  const double m,
5103  const double thetaR,
5104  const double thetaSR,
5105  const double KWs,
5106  double *dV,
5107  double *u,
5108  double *mass,
5109  double *dmass,
5110  double *f,
5111  double *df,
5112  double *a,
5113  double *da)
5114 {
5115  int k,I,J,eN;
5116  const int nSpace2=nSpace*nSpace;
5117  register double psiC,
5118  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5119  onePlus_pcBar_n,
5120  sBar,sqrt_sBar,DsBar_DpsiC,
5121  thetaW,DthetaW_DpsiC,
5122  vBar,vBar2,DvBar_DpsiC,
5123  KW,DKW_DpsiC,
5124  rho2=rho*rho,
5125  thetaS=thetaR+thetaSR;
5126 
5127  double mavg,dmavg,vol;
5128  double favg[3] = {0.0,0.0,0.0};
5129  double dfavg[3] = {0.0,0.0,0.0};
5130  double aavg[3][3]= {{0.0,0.0,0.0},
5131  {0.0,0.0,0.0},
5132  {0.0,0.0,0.0}};
5133  double daavg[3][3]= {{0.0,0.0,0.0},
5134  {0.0,0.0,0.0},
5135  {0.0,0.0,0.0}};
5136  /*mwf debug
5137  printf("revgL2proj nSimplices=%d nPerSimp=%d nSpace=%d\n",nSimplices,nPointsPerSimplex,
5138  nSpace);
5139  */
5140  for (eN = 0; eN < nSimplices; eN++)
5141  {
5142  mavg = 0.0; dmavg = 0.0; vol = 0.0;
5143  for (I=0; I < nSpace; I++)
5144  {
5145  favg[I] =0.0; dfavg[I] = 0.0;
5146  for (J=0; J < nSpace;J++)
5147  {
5148  aavg[I][J] = 0.0; daavg[I][J] = 0.0;
5149  }
5150  }
5151  for (k=0;k<nPointsPerSimplex;k++)
5152  {
5153  psiC = -u[eN*nPointsPerSimplex + k];
5154  if (psiC > 0.0)
5155  {
5156  pcBar = alpha*psiC;
5157  pcBar_nM2 = pow(pcBar,n-2);
5158  pcBar_nM1 = pcBar_nM2*pcBar;
5159  pcBar_n = pcBar_nM1*pcBar;
5160  onePlus_pcBar_n = 1.0 + pcBar_n;
5161 
5162  sBar = pow(onePlus_pcBar_n,-m);
5163  /* using -mn = 1-n */
5164  DsBar_DpsiC = alpha*(1.0-n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5165 
5166  vBar = 1.0-pcBar_nM1*sBar;
5167  vBar2 = vBar*vBar;
5168  DvBar_DpsiC = -alpha*(n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5169 
5170  thetaW = thetaSR*sBar + thetaR;
5171  DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5172 
5173  sqrt_sBar = sqrt(sBar);
5174  KW= KWs*sqrt_sBar*vBar2;
5175  DKW_DpsiC= KWs*
5176  ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5177  +
5178  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5179  }
5180  else
5181  {
5182  thetaW = thetaS;
5183  DthetaW_DpsiC = 0.0;
5184  KW = KWs;
5185  DKW_DpsiC = 0.0;
5186  }
5187  /*mwf debug
5188  printf("revgmL2 eN=%d k=%d psiC=%g thetaW=%g KW=%g rho=%g DKW_DpsiC=%g\n",
5189  eN,k,psiC,thetaW,KW,rho,DKW_DpsiC);
5190  */
5191  vol += dV[eN*nPointsPerSimplex + k];
5192  mavg += rho*thetaW*dV[eN*nPointsPerSimplex + k];
5193  dmavg +=-rho*DthetaW_DpsiC*dV[eN*nPointsPerSimplex + k];
5194  /*go ahead and assign point values for derivs since this is what's actually
5195  appropriate at least for stiffness and advection terms,
5196  but it's not correct for mass term*/
5197  dmass[eN*nPointsPerSimplex+ k] =-rho*DthetaW_DpsiC;
5198 
5199  for (I=0; I < nSpace; I++)
5200  {
5201  favg[I] += rho2*KW*gravity[I]*dV[eN*nPointsPerSimplex + k];
5202  dfavg[I]+=-rho2*DKW_DpsiC*gravity[I]*dV[eN*nPointsPerSimplex + k];
5203  J=I;/*assume isotropic for now*/
5204  aavg[I][J] += rho*KW*dV[eN*nPointsPerSimplex + k];
5205  daavg[I][J] +=-rho*DKW_DpsiC*dV[eN*nPointsPerSimplex + k];
5206 
5207  df[eN*nPointsPerSimplex*nSpace+ k*nSpace + I] = -rho2*DKW_DpsiC*gravity[I];
5208  da[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J] = -rho*DKW_DpsiC;
5209  }
5210  }/*end k 1*/
5211  assert(vol > 0.0);
5212  /*mwf debug
5213  printf("eN=%d vol=%g mavg=%g dmavg=%g favg=[%g,%g,%g] aavg=[%g,%g,%g] \n",
5214  eN,vol,mavg,dmavg,favg[0],favg[1],favg[2],aavg[0][0],aavg[1][1],aavg[2][2]);
5215  */
5216  for (k=0; k < nPointsPerSimplex; k++)
5217  {
5218  mass[eN*nPointsPerSimplex + k] = mavg/vol;
5219  /*dmass[eN*nPointsPerSimplex+ k] = dmavg/vol;*/
5220  for (I=0; I < nSpace; I++)
5221  {
5222  f[eN*nPointsPerSimplex*nSpace + k*nSpace + I] = favg[I]/vol;
5223  /*df[eN*nPointsPerSimplex*nSpace+ k*nSpace + I] = dfavg[I]/vol;*/
5224  for (J=0; J < nSpace; J++)
5225  {
5226  if (J == I)
5227  {
5228  a[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J] = aavg[I][J]/vol;
5229  /*da[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J]= daavg[I][J]/vol;*/
5230  }
5231  else
5232  {
5233  a[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J] = 0.0;
5234  /*da[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J]= 0.0;*/
5235  }
5236  }
5237 
5238  }/*I*/
5239  }/*k*/
5240  }
5241 }
5242 /* TODO: figure out how to handle dmass term to get right jacobian*/
5243 void conservativeHeadRichardsL2projBndMualemVanGenuchtenHomEvaluate(const int nElements,
5244  const int nElementBoundaries_element,
5245  const int nPointsPerElementBoundary,
5246  const int nSpace,
5247  const double rho,
5248  const double* gravity,
5249  const double alpha,
5250  const double n,
5251  const double m,
5252  const double thetaR,
5253  const double thetaSR,
5254  const double KWs,
5255  double *dV,
5256  double *u,
5257  double *mass,
5258  double *dmass,
5259  double *f,
5260  double *df,
5261  double *a,
5262  double *da)
5263 {
5264  int k,I,J,eN,ebN;
5265  const int nSpace2=nSpace*nSpace;
5266  register double psiC,
5267  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5268  onePlus_pcBar_n,
5269  sBar,sqrt_sBar,DsBar_DpsiC,
5270  thetaW,DthetaW_DpsiC,
5271  vBar,vBar2,DvBar_DpsiC,
5272  KW,DKW_DpsiC,
5273  rho2=rho*rho,
5274  thetaS=thetaR+thetaSR;
5275 
5276  double mavg,dmavg,vol;
5277  double favg[3] = {0.0,0.0,0.0};
5278  double dfavg[3] = {0.0,0.0,0.0};
5279  double aavg[3][3]= {{0.0,0.0,0.0},
5280  {0.0,0.0,0.0},
5281  {0.0,0.0,0.0}};
5282  double daavg[3][3]= {{0.0,0.0,0.0},
5283  {0.0,0.0,0.0},
5284  {0.0,0.0,0.0}};
5285  /*mwf debug
5286  printf("revgL2projBnd nElements=%d nElemBndp=%d nPtsPerBnd=%d nSpace=%d\n",nElements,nElementBoundaries_element,
5287  nPointsPerElementBoundary,
5288  nSpace);
5289  */
5290  for (eN = 0; eN < nElements; eN++)
5291  {
5292  for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
5293  {
5294  mavg = 0.0; dmavg = 0.0; vol = 0.0;
5295  for (I=0; I < nSpace; I++)
5296  {
5297  favg[I] =0.0; dfavg[I] = 0.0;
5298  for (J=0; J < nSpace;J++)
5299  {
5300  aavg[I][J] = 0.0; daavg[I][J] = 0.0;
5301  }
5302  }
5303  for (k=0;k<nPointsPerElementBoundary;k++)
5304  {
5305  psiC = -u[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k];
5306  if (psiC > 0.0)
5307  {
5308  pcBar = alpha*psiC;
5309  pcBar_nM2 = pow(pcBar,n-2);
5310  pcBar_nM1 = pcBar_nM2*pcBar;
5311  pcBar_n = pcBar_nM1*pcBar;
5312  onePlus_pcBar_n = 1.0 + pcBar_n;
5313 
5314  sBar = pow(onePlus_pcBar_n,-m);
5315  /* using -mn = 1-n */
5316  DsBar_DpsiC = alpha*(1.0-n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5317 
5318  vBar = 1.0-pcBar_nM1*sBar;
5319  vBar2 = vBar*vBar;
5320  DvBar_DpsiC = -alpha*(n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5321 
5322  thetaW = thetaSR*sBar + thetaR;
5323  DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5324 
5325  sqrt_sBar = sqrt(sBar);
5326  KW= KWs*sqrt_sBar*vBar2;
5327  DKW_DpsiC= KWs*
5328  ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5329  +
5330  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5331  }
5332  else
5333  {
5334  thetaW = thetaS;
5335  DthetaW_DpsiC = 0.0;
5336  KW = KWs;
5337  DKW_DpsiC = 0.0;
5338  }
5339  /*mwf debug
5340  printf("revgmL2 eN=%d ebN=%d k=%d psiC=%g thetaW=%g KW=%g rho=%g DKW_DpsiC=%g\n",
5341  eN,ebN,k,psiC,thetaW,KW,rho,DKW_DpsiC);
5342  */
5343  vol += dV[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k];
5344  mavg += rho*thetaW*dV[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k];
5345  dmavg +=-rho*DthetaW_DpsiC*dV[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k];
5346  /*mwf go ahead and assign point values for derivs since this is what's needed for stiffness terms but
5347  even though it's not the right thing to do for mass term*/
5348  dmass[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k] =
5349  -rho*DthetaW_DpsiC;
5350  for (I=0; I < nSpace; I++)
5351  {
5352  favg[I] += rho2*KW*gravity[I]*dV[eN*nElementBoundaries_element*nPointsPerElementBoundary +
5353  ebN*nPointsPerElementBoundary + k];
5354  dfavg[I]+=-rho2*DKW_DpsiC*gravity[I]*dV[eN*nElementBoundaries_element*nPointsPerElementBoundary +
5355  ebN*nPointsPerElementBoundary + k];
5356  J = I; /*assume isotropic*/
5357  aavg[I][J] += rho*KW*dV[eN*nElementBoundaries_element*nPointsPerElementBoundary +
5358  ebN*nPointsPerElementBoundary + k];
5359  daavg[I][J] +=-rho*DKW_DpsiC*dV[eN*nElementBoundaries_element*nPointsPerElementBoundary +
5360  ebN*nPointsPerElementBoundary + k];
5361 
5362  df[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace + ebN*nPointsPerElementBoundary*nSpace +
5363  k*nSpace + I] = -rho2*DKW_DpsiC*gravity[I];
5364  da[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace2+ ebN*nPointsPerElementBoundary*nSpace2 +
5365  k*nSpace2 + I*nSpace + J]= -rho*DKW_DpsiC;
5366 
5367  }
5368  }/*end k 1*/
5369  assert(vol > 0.0);
5370  /*mwf debug
5371  printf("eN=%d vol=%g mavg=%g dmavg=%g favg=[%g,%g,%g] aavg=[%g,%g,%g] \n",
5372  eN,vol,mavg,dmavg,favg[0],favg[1],favg[2],aavg[0][0],aavg[1][1],aavg[2][2]);
5373  */
5374  for (k=0; k < nPointsPerElementBoundary; k++)
5375  {
5376  mass[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k] = mavg/vol;
5377  /*dmass[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k] =
5378  dmavg/vol;*/
5379  for (I=0; I < nSpace; I++)
5380  {
5381  f[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace + ebN*nPointsPerElementBoundary*nSpace +
5382  k*nSpace + I] = favg[I]/vol;
5383  /*df[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace + ebN*nPointsPerElementBoundary*nSpace +
5384  k*nSpace + I] = dfavg[I]/vol;*/
5385  for (J=0; J < nSpace; J++)
5386  {
5387  /*assume diagonal for now*/
5388  if (J == I)
5389  {
5390  a[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace2 + ebN*nPointsPerElementBoundary*nSpace2 +
5391  k*nSpace2 + I*nSpace + J] = aavg[I][J]/vol;
5392  /*da[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace2+
5393  ebN*nPointsPerElementBoundary*nSpace2 +
5394  k*nSpace2 + I*nSpace + J]= daavg[I][J]/vol;*/
5395 
5396  }
5397  else
5398  {
5399  a[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace2 + ebN*nPointsPerElementBoundary*nSpace2 +
5400  k*nSpace2 + I*nSpace + J] = 0.0;
5401  /*da[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace2+
5402  ebN*nPointsPerElementBoundary*nSpace2 +
5403  k*nSpace2 + I*nSpace + J]= 0.0;*/
5404 
5405  }
5406  }
5407  }/*I*/
5408  }/*k*/
5409  }/*ebN*/
5410  }/*eN*/
5411 }
5412 void conservativeHeadRichardsL2projMualemVanGenuchtenHetEvaluate(const int nSimplices,
5413  const int nPointsPerSimplex,
5414  const int nSpace,
5415  const double rho,
5416  const double *gravity,
5417  const double *alpha,
5418  const double *n,
5419  const double *thetaR,
5420  const double *thetaSR,
5421  const double *KWs,
5422  double *dV,
5423  double *u,
5424  double *mass,
5425  double *dmass,
5426  double *f,
5427  double *df,
5428  double *a,
5429  double *da)
5430 {
5431  int k,I,J,eN;
5432  const int nSpace2=nSpace*nSpace;
5433  register double psiC,
5434  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5435  onePlus_pcBar_n,
5436  sBar,sqrt_sBar,DsBar_DpsiC,
5437  thetaW,DthetaW_DpsiC,
5438  vBar,vBar2,DvBar_DpsiC,
5439  KW,DKW_DpsiC,
5440  rho2=rho*rho;
5441  register double thetaS,m;
5442 
5443  double mavg,dmavg,vol;
5444  double favg[3] = {0.0,0.0,0.0};
5445  double dfavg[3] = {0.0,0.0,0.0};
5446  double aavg[3][3]= {{0.0,0.0,0.0},
5447  {0.0,0.0,0.0},
5448  {0.0,0.0,0.0}};
5449  double daavg[3][3]= {{0.0,0.0,0.0},
5450  {0.0,0.0,0.0},
5451  {0.0,0.0,0.0}};
5452  /*mwf debug
5453  printf("revgL2proj nSimplices=%d nPerSimp=%d nSpace=%d\n",nSimplices,nPointsPerSimplex,
5454  nSpace);
5455  */
5456  for (eN = 0; eN < nSimplices; eN++)
5457  {
5458  mavg = 0.0; dmavg = 0.0; vol = 0.0;
5459  for (I=0; I < nSpace; I++)
5460  {
5461  favg[I] =0.0; dfavg[I] = 0.0;
5462  for (J=0; J < nSpace;J++)
5463  {
5464  aavg[I][J] = 0.0; daavg[I][J] = 0.0;
5465  }
5466  }
5467  for (k=0;k<nPointsPerSimplex;k++)
5468  {
5469  psiC = -u[eN*nPointsPerSimplex + k];
5470  m = 1.0 - 1.0/n[eN*nPointsPerSimplex + k];
5471  thetaS = thetaR[eN*nPointsPerSimplex + k] + thetaSR[eN*nPointsPerSimplex + k];
5472  if (psiC > 0.0)
5473  {
5474  pcBar = alpha[eN*nPointsPerSimplex + k]*psiC;
5475  pcBar_nM2 = pow(pcBar,n[eN*nPointsPerSimplex + k]-2);
5476  pcBar_nM1 = pcBar_nM2*pcBar;
5477  pcBar_n = pcBar_nM1*pcBar;
5478  onePlus_pcBar_n = 1.0 + pcBar_n;
5479 
5480  sBar = pow(onePlus_pcBar_n,-m);
5481  /* using -mn = 1-n */
5482  DsBar_DpsiC = alpha[eN*nPointsPerSimplex + k]*(1.0-n[eN*nPointsPerSimplex + k])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5483 
5484  vBar = 1.0-pcBar_nM1*sBar;
5485  vBar2 = vBar*vBar;
5486  DvBar_DpsiC = -alpha[eN*nPointsPerSimplex + k]*(n[eN*nPointsPerSimplex + k]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5487 
5488  thetaW = thetaSR[eN*nPointsPerSimplex + k]*sBar + thetaR[eN*nPointsPerSimplex + k];
5489  DthetaW_DpsiC = thetaSR[eN*nPointsPerSimplex + k] * DsBar_DpsiC;
5490 
5491  sqrt_sBar = sqrt(sBar);
5492  KW= KWs[eN*nPointsPerSimplex + k]*sqrt_sBar*vBar2;
5493  DKW_DpsiC= KWs[eN*nPointsPerSimplex + k]*
5494  ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5495  +
5496  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5497  }
5498  else
5499  {
5500  thetaW = thetaS;
5501  DthetaW_DpsiC = 0.0;
5502  KW = KWs[eN*nPointsPerSimplex + k];
5503  DKW_DpsiC = 0.0;
5504  }
5505  /*mwf debug
5506  printf("revgmL2 eN=%d k=%d psiC=%g thetaW=%g KW=%g rho=%g DKW_DpsiC=%g\n",
5507  eN,k,psiC,thetaW,KW,rho,DKW_DpsiC);
5508  */
5509  vol += dV[eN*nPointsPerSimplex + k];
5510  mavg += rho*thetaW*dV[eN*nPointsPerSimplex + k];
5511  dmavg +=-rho*DthetaW_DpsiC*dV[eN*nPointsPerSimplex + k];
5512  /*go ahead and assign point values for deriv terms since this is right thing to do
5513  for stiffness terms even though it's not right for mass terms*/
5514  dmass[eN*nPointsPerSimplex+ k] = -rho*DthetaW_DpsiC;
5515  for (I=0; I < nSpace; I++)
5516  {
5517  favg[I] += rho2*KW*gravity[I]*dV[eN*nPointsPerSimplex + k];
5518  dfavg[I]+=-rho2*DKW_DpsiC*gravity[I]*dV[eN*nPointsPerSimplex + k];
5519 
5520  df[eN*nPointsPerSimplex*nSpace+ k*nSpace + I] = -rho2*DKW_DpsiC*gravity[I];
5521  J=I;
5522  aavg[I][J] += rho*KW*dV[eN*nPointsPerSimplex + k];
5523  daavg[I][J] +=-rho*DKW_DpsiC*dV[eN*nPointsPerSimplex + k];
5524  da[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + I] = -rho*DKW_DpsiC;
5525  }
5526  }/*end k 1*/
5527  assert(vol > 0.0);
5528  /*mwf debug
5529  printf("eN=%d vol=%g mavg=%g dmavg=%g favg=[%g,%g,%g] aavg=[%g,%g,%g] \n",
5530  eN,vol,mavg,dmavg,favg[0],favg[1],favg[2],aavg[0][0],aavg[1][1],aavg[2][2]);
5531  */
5532  for (k=0; k < nPointsPerSimplex; k++)
5533  {
5534  mass[eN*nPointsPerSimplex + k] = mavg/vol;
5535  /*dmass[eN*nPointsPerSimplex+ k] = dmavg/vol;*/
5536  for (I=0; I < nSpace; I++)
5537  {
5538  /*assume diagonal*/
5539  f[eN*nPointsPerSimplex*nSpace + k*nSpace + I] = favg[I]/vol;
5540  /*df[eN*nPointsPerSimplex*nSpace+ k*nSpace + I] = dfavg[I]/vol;*/
5541  a[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + I] = aavg[I][I]/vol;
5542  /*da[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + I]= daavg[I][I]/vol;*/
5543  }/*I*/
5544  }/*k*/
5545  }
5546 }
5547 
5550 void conservativeTotalHeadRichardsMualemVanGenuchtenHomEvaluate(const int nPoints,
5551  const int nSpace,
5552  const double rho,
5553  const double* gravity,
5554  const double* x,
5555  const double alpha,
5556  const double n,
5557  const double m,
5558  const double thetaR,
5559  const double thetaSR,
5560  const double KWs,
5561  double *u,
5562  double *mass,
5563  double *dmass,
5564  double *f,
5565  double *df,
5566  double *a,
5567  double *da,
5568  double *phi,
5569  double *dphi)
5570 {
5571  int k,I;
5572  const int nSpace2=nSpace*nSpace;
5573  register double psiC,
5574  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5575  onePlus_pcBar_n,
5576  sBar,sqrt_sBar,DsBar_DpsiC,
5577  thetaW,DthetaW_DpsiC,
5578  vBar,vBar2,DvBar_DpsiC,
5579  KW,DKW_DpsiC,
5580  thetaS=thetaR+thetaSR;
5581  /*mwf elevation */
5582  register double elev;
5583  for (k=0;k<nPoints;k++)
5584  {
5585  elev = 0.0;
5586  for (I=0; I < nSpace; I++)
5587  elev += gravity[I]*x[k*3+I];
5588  /*mwf if unknown is h
5589  psiC = -u[k]-elev;
5590  */
5591  /*mwf if unknown is psi*/
5592  psiC = -u[k];
5593 
5594  if (psiC > 0.0)
5595  {
5596  pcBar = alpha*psiC;
5597  pcBar_nM2 = pow(pcBar,n-2);
5598  pcBar_nM1 = pcBar_nM2*pcBar;
5599  pcBar_n = pcBar_nM1*pcBar;
5600  onePlus_pcBar_n = 1.0 + pcBar_n;
5601 
5602  sBar = pow(onePlus_pcBar_n,-m);
5603  /* using -mn = 1-n */
5604  DsBar_DpsiC = alpha*(1.0-n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5605 
5606  vBar = 1.0-pcBar_nM1*sBar;
5607  vBar2 = vBar*vBar;
5608  DvBar_DpsiC = -alpha*(n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5609 
5610  thetaW = thetaSR*sBar + thetaR;
5611  DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5612 
5613  sqrt_sBar = sqrt(sBar);
5614  KW= KWs*sqrt_sBar*vBar2;
5615  DKW_DpsiC= KWs*
5616  ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5617  +
5618  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5619  }
5620  else
5621  {
5622  thetaW = thetaS;
5623  DthetaW_DpsiC = 0.0;
5624  KW = KWs;
5625  DKW_DpsiC = 0.0;
5626  }
5627  mass[k] = rho*thetaW;
5628  dmass[k] = -rho*DthetaW_DpsiC;
5629  /* remember to turn on nonlinear potential if you want to use total head */
5630  /*mwf unknown is psi*/
5631  phi[k] = -psiC;
5632  dphi[k] = 1.0;
5633 
5634  /*mwf unknown is h
5635  phi[k] = u[k];
5636  dphi[k] = 1.0;
5637  */
5638  for (I=0;I<nSpace;I++)
5639  {
5640  f[k*nSpace+I] = 0.0;
5641  df[k*nSpace+I] = 0.0;
5642  /*mwf unknown is psi*/
5643  phi[k] -= rho*gravity[I]*x[k*3+I];
5644 
5645  a[k*nSpace2+I*nSpace+I] = rho*KW;
5646  da[k*nSpace2+I*nSpace+I] = -rho*DKW_DpsiC;
5647  }
5648  }
5649 }
5650 
5651 void l2projectScalar(const int nSimplices,
5652  const int nPointsPerSimplex,
5653  double * dV,
5654  double * r)
5655 {
5656  int eN,k;
5657  double ravg,vol;
5658 
5659  for (eN = 0; eN < nSimplices; eN++)
5660  {
5661  ravg = 0.0; vol = 0.0;
5662  for (k=0; k < nPointsPerSimplex; k++)
5663  {
5664  vol += dV[eN*nPointsPerSimplex+k];
5665  ravg += r[eN*nPointsPerSimplex + k]*dV[eN*nPointsPerSimplex+k];
5666  }
5667 
5668  assert(vol > 0.0);
5669  ravg = ravg / vol;
5670  for (k=0; k < nPointsPerSimplex; k++)
5671  {
5672  r[eN*nPointsPerSimplex + k] = ravg;
5673  }
5674  }
5675 }
5676 void l2projectVector(const int nSimplices,
5677  const int nPointsPerSimplex,
5678  const int nSpace,
5679  double * dV,
5680  double * r)
5681 {
5682  int eN,k,I;
5683  double vol;
5684  /*need different max size, take a chance on variable length array?*/
5685  double ravg[nSpace];
5686  for (eN = 0; eN < nSimplices; eN++)
5687  {
5688  vol = 0.0;
5689  for (I=0; I < nSpace; I++)
5690  ravg[I] = 0.0;
5691  for (k=0; k < nPointsPerSimplex; k++)
5692  {
5693  vol += dV[eN*nPointsPerSimplex+k];
5694  for (I=0; I < nSpace; I++)
5695  ravg[I] += r[eN*nPointsPerSimplex*nSpace + k*nSpace + I]*dV[eN*nPointsPerSimplex+k];
5696  }
5697  assert(vol > 0.0);
5698  for (k=0; k < nPointsPerSimplex; k++)
5699  {
5700  for (I=0; I < nSpace; I++)
5701  r[eN*nPointsPerSimplex*nSpace + k*nSpace + I] = ravg[I]/vol;
5702  }
5703  }
5704 }
5705 void l2project2Tensor(const int nSimplices,
5706  const int nPointsPerSimplex,
5707  const int nSpace,
5708  double * dV,
5709  double * r)
5710 {
5711  int eN,k,I,J,nSpace2;
5712  double vol;
5713  /*need different max size, take a chance on variable length array?*/
5714  double ravg[nSpace*nSpace];
5715  nSpace2 = nSpace*nSpace;
5716  for (eN = 0; eN < nSimplices; eN++)
5717  {
5718  vol = 0.0;
5719  for (I=0; I < nSpace2; I++)
5720  ravg[I] = 0.0;
5721  for (k=0; k < nPointsPerSimplex; k++)
5722  {
5723  vol += dV[eN*nPointsPerSimplex+k];
5724  for (I=0; I < nSpace; I++)
5725  for (J=0; J < nSpace; J++)
5726  ravg[I*nSpace+J] += r[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J]
5727  *
5728  dV[eN*nPointsPerSimplex+k];
5729  }
5730  assert(vol > 0.0);
5731  for (k=0; k < nPointsPerSimplex; k++)
5732  {
5733  for (I=0; I < nSpace; I++)
5734  for (J=0; J < nSpace; J++)
5735  r[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J] = ravg[I*nSpace+J]/vol;
5736  }
5737  }
5738 }
5739 
5740 void conservativeHeadRichardsMualemVanGenuchten_sd_het(const int nSimplex,
5741  const int nPointsPerSimplex,
5742  const int nSpace,
5743  double pc_eps,
5744  const int * rowptr,
5745  const int * colind,
5746  const int* materialTypes,
5747  const double rho,
5748  const double beta,
5749  const double* gravity,
5750  const double* alpha,
5751  const double* n,
5752  const double* thetaR,
5753  const double* thetaSR,
5754  const double* KWs,
5755  double *u,
5756  double *mass,
5757  double *dmass,
5758  double *f,
5759  double *df,
5760  double *a,
5761  double *da,
5762  double* vol_frac)
5763 {
5764  int i,j,k,I,matID,ii;
5765  const int nSpace2=nSpace*nSpace;
5766  register double psiC,pcBarStar,
5767  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5768  onePlus_pcBar_n,
5769  sBar,sqrt_sBar,DsBar_DpsiC,
5770  thetaW,DthetaW_DpsiC,
5771  vBar,vBar2,DvBar_DpsiC,
5772  KWr,DKWr_DpsiC,
5773  rho2=rho*rho,
5774  thetaS,
5775  rhom,drhom,m,
5776  betauStar;
5777  const int nnz = rowptr[nSpace];
5778  for (i=0; i < nSimplex; i++)
5779  {
5780  matID= materialTypes[i];
5781  for (j=0;j<nPointsPerSimplex;j++)
5782  {
5783  k = i*nPointsPerSimplex + j;
5784  psiC = -u[k];
5785  m = 1.0 - 1.0/n[matID];
5786  thetaS = thetaR[matID] + thetaSR[matID];
5787  if (psiC > 0.0)
5788  {
5789  pcBar = alpha[matID]*psiC;
5790  pcBarStar = fmax(pcBar,pc_eps);
5791 
5792  pcBar_nM2 = pow(pcBarStar,n[matID]-2);
5793  pcBar_nM1 = pcBar_nM2*pcBar;
5794  pcBar_n = pcBar_nM1*pcBar;
5795  onePlus_pcBar_n = 1.0 + pcBar_n;
5796 
5797  sBar = pow(onePlus_pcBar_n,-m);
5798  /* using -mn = 1-n */
5799  DsBar_DpsiC = alpha[matID]*(1.0-n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5800 
5801  vBar = 1.0-pcBar_nM1*sBar;
5802  vBar2 = vBar*vBar;
5803  DvBar_DpsiC = -alpha[matID]*(n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5804 
5805  thetaW = thetaSR[matID]*sBar + thetaR[matID];
5806  DthetaW_DpsiC = thetaSR[matID] * DsBar_DpsiC;
5807 
5808  sqrt_sBar = sqrt(sBar);
5809  KWr= sqrt_sBar*vBar2;
5810  DKWr_DpsiC= ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5811  +
5812  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5813  }
5814  else
5815  {
5816  thetaW = thetaS;
5817  DthetaW_DpsiC = 0.0;
5818  KWr = 1.0;
5819  DKWr_DpsiC = 0.0;
5820  }
5821  //slight compressibility
5822  betauStar = fmin(beta*u[k],1000.0);
5823  rhom = rho*exp(betauStar);
5824  drhom = beta*rhom;
5825  mass[k] = rhom*thetaW;
5826  dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
5827  vol_frac[k] = thetaW;
5828  //mass[k] = rho*thetaW;
5829  //dmass[k] = -rho*DthetaW_DpsiC;
5830  for (I=0;I<nSpace;I++)
5831  {
5832  f[k*nSpace+I] = 0.0;
5833  df[k*nSpace+I] = 0.0;
5834  for (ii=rowptr[I]; ii < rowptr[I+1]; ii++)
5835  {
5836  f[k*nSpace+I] += rho2*KWr*KWs[matID*nnz+ii]*gravity[colind[ii]];
5837  df[k*nSpace+I] += -rho2*DKWr_DpsiC*KWs[matID*nnz+ii]*gravity[colind[ii]];
5838  a[k*nnz+ii] = rho*KWr*KWs[matID*nnz+ii];
5839  da[k*nnz+ii] = -rho*DKWr_DpsiC*KWs[matID*nnz+ii];
5840  }/*m*/
5841  }/*I*/
5842  }/*k*/
5843  }/*j*/
5844 
5845 }
5846 
5847 void conservativeHeadRichardsMualemVanGenuchten_sd_het_linearized_at_saturation(const int nSimplex,
5848  const int nPointsPerSimplex,
5849  const int nSpace,
5850  double linear_break,
5851  const int * rowptr,
5852  const int * colind,
5853  const int* materialTypes,
5854  const double rho,
5855  const double beta,
5856  const double* gravity,
5857  const double* alpha,
5858  const double* n,
5859  const double* thetaR,
5860  const double* thetaSR,
5861  const double* KWs,
5862  double *u,
5863  double *mass,
5864  double *dmass,
5865  double *f,
5866  double *df,
5867  double *a,
5868  double *da,
5869  double* vol_frac)
5870 {
5871  int i,j,k,I,matID,ii;
5872  const int nSpace2=nSpace*nSpace;
5873  register double psiC,pcBarStar,sBarStar,KWrStar,
5874  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5875  onePlus_pcBar_n,
5876  sBar,sqrt_sBar,DsBar_DpsiC,
5877  thetaW,DthetaW_DpsiC,
5878  vBar,vBar2,DvBar_DpsiC,
5879  KWr,DKWr_DpsiC,
5880  rho2=rho*rho,
5881  thetaS,
5882  rhom,drhom,m,
5883  betauStar;
5884  const int nnz = rowptr[nSpace];
5885  pcBarStar = linear_break;
5886 
5887  for (i=0; i < nSimplex; i++)
5888  {
5889  matID= materialTypes[i];
5890  for (j=0;j<nPointsPerSimplex;j++)
5891  {
5892  k = i*nPointsPerSimplex + j;
5893  psiC = -u[k];
5894  m = 1.0 - 1.0/n[matID];
5895  thetaS = thetaR[matID] + thetaSR[matID];
5896  if (psiC > 0.0)
5897  {
5898  pcBar = alpha[matID]*psiC;
5899  if (psiC <= pcBarStar)
5900  {
5901  /*calculate regularization parameters again because n varies*/
5902  pcBar_nM2 = pow(pcBarStar,n[matID]-2.);
5903  pcBar_nM1 = pcBar_nM2*pcBar;
5904  pcBar_n = pcBar_nM1*pcBar;
5905  onePlus_pcBar_n = 1.0 + pcBar_n;
5906 
5907  sBar = pow(onePlus_pcBar_n,-m);
5908  /* using -mn = 1-n */
5909  DsBar_DpsiC = alpha[matID]*(1.0-n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5910 
5911  vBar = 1.0-pcBar_nM1*sBar;
5912  vBar2 = vBar*vBar;
5913  DvBar_DpsiC = -alpha[matID]*(n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5914 
5915  sBarStar = sBar;
5916  sqrt_sBar = sqrt(sBar);
5917  KWrStar= sqrt_sBar*vBar2;
5918 
5919  /*now compute linearized solution*/
5920  DsBar_DpsiC = (sBarStar-1.0)/(pcBarStar-0.0);
5921  sBar = DsBar_DpsiC*(psiC-0.0) + 1.0;
5922  thetaW = thetaSR[matID]*sBar + thetaR[matID];
5923  DthetaW_DpsiC = thetaSR[matID] * DsBar_DpsiC;
5924 
5925  DKWr_DpsiC= (KWrStar - 1.0)/(pcBarStar-0.0);
5926  KWr = DKWr_DpsiC*(psiC-0.0) + 1.0;
5927  }
5928  else
5929  {
5930  pcBar = alpha[matID]*psiC;
5931  pcBarStar = fmax(pcBar,1.0e-8);
5932 
5933  pcBar_nM2 = pow(pcBarStar,n[matID]-2);
5934  pcBar_nM1 = pcBar_nM2*pcBar;
5935  pcBar_n = pcBar_nM1*pcBar;
5936  onePlus_pcBar_n = 1.0 + pcBar_n;
5937 
5938  sBar = pow(onePlus_pcBar_n,-m);
5939  /* using -mn = 1-n */
5940  DsBar_DpsiC = alpha[matID]*(1.0-n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5941 
5942  vBar = 1.0-pcBar_nM1*sBar;
5943  vBar2 = vBar*vBar;
5944  DvBar_DpsiC = -alpha[matID]*(n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5945 
5946  thetaW = thetaSR[matID]*sBar + thetaR[matID];
5947  DthetaW_DpsiC = thetaSR[matID] * DsBar_DpsiC;
5948 
5949  sqrt_sBar = sqrt(sBar);
5950  KWr= sqrt_sBar*vBar2;
5951  DKWr_DpsiC= ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5952  +
5953  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5954 
5955  }
5956  }
5957  else
5958  {
5959  thetaW = thetaS;
5960  DthetaW_DpsiC = 0.0;
5961  KWr = 1.0;
5962  DKWr_DpsiC = 0.0;
5963  }
5964  //slight compressibility
5965  betauStar = fmin(beta*u[k],1000.0);
5966  rhom = rho*exp(betauStar);
5967  drhom = beta*rhom;
5968  mass[k] = rhom*thetaW;
5969  dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
5970  vol_frac[k] = thetaW;
5971  //mass[k] = rho*thetaW;
5972  //dmass[k] = -rho*DthetaW_DpsiC;
5973  for (I=0;I<nSpace;I++)
5974  {
5975  f[k*nSpace+I] = 0.0;
5976  df[k*nSpace+I] = 0.0;
5977  for (ii=rowptr[I]; ii < rowptr[I+1]; ii++)
5978  {
5979  f[k*nSpace+I] += rho2*KWr*KWs[matID*nnz+ii]*gravity[colind[ii]];
5980  df[k*nSpace+I] += -rho2*DKWr_DpsiC*KWs[matID*nnz+ii]*gravity[colind[ii]];
5981  a[k*nnz+ii] = rho*KWr*KWs[matID*nnz+ii];
5982  da[k*nnz+ii] = -rho*DKWr_DpsiC*KWs[matID*nnz+ii];
5983  }/*m*/
5984  }/*I*/
5985  }/*k*/
5986  }/*j*/
5987 
5988 }
5989 
5990 void conservativeHeadRichardsMualemVanGenuchtenHetEvaluateV2(const int nSimplex,
5991  const int nPointsPerSimplex,
5992  const int nSpace,
5993  const int* materialTypes,
5994  const double rho,
5995  const double beta,
5996  const double* gravity,
5997  const double* alpha,
5998  const double* n,
5999  const double* thetaR,
6000  const double* thetaSR,
6001  const double* KWs,
6002  double *u,
6003  double *mass,
6004  double *dmass,
6005  double *f,
6006  double *df,
6007  double *a,
6008  double *da)
6009 {
6010  int i,j,k,I,matID;
6011  const int nSpace2=nSpace*nSpace;
6012  register double psiC,
6013  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
6014  onePlus_pcBar_n,
6015  sBar,sqrt_sBar,DsBar_DpsiC,
6016  thetaW,DthetaW_DpsiC,
6017  vBar,vBar2,DvBar_DpsiC,
6018  KW,DKW_DpsiC,
6019  rho2=rho*rho,
6020  thetaS,
6021  rhom,drhom,m;
6022  for (i=0; i < nSimplex; i++)
6023  {
6024  matID= materialTypes[i];
6025  for (j=0;j<nPointsPerSimplex;j++)
6026  {
6027  k = i*nPointsPerSimplex + j;
6028  psiC = -u[k];
6029  m = 1.0 - 1.0/n[matID];
6030  thetaS = thetaR[matID] + thetaSR[matID];
6031  if (psiC > 0.0)
6032  {
6033  pcBar = alpha[matID]*psiC;
6034  pcBar_nM2 = pow(pcBar,n[matID]-2);
6035  pcBar_nM1 = pcBar_nM2*pcBar;
6036  pcBar_n = pcBar_nM1*pcBar;
6037  onePlus_pcBar_n = 1.0 + pcBar_n;
6038 
6039  sBar = pow(onePlus_pcBar_n,-m);
6040  /* using -mn = 1-n */
6041  DsBar_DpsiC = alpha[matID]*(1.0-n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
6042 
6043  vBar = 1.0-pcBar_nM1*sBar;
6044  vBar2 = vBar*vBar;
6045  DvBar_DpsiC = -alpha[matID]*(n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
6046 
6047  thetaW = thetaSR[matID]*sBar + thetaR[matID];
6048  DthetaW_DpsiC = thetaSR[matID] * DsBar_DpsiC;
6049 
6050  sqrt_sBar = sqrt(sBar);
6051  KW= KWs[matID]*sqrt_sBar*vBar2;
6052  DKW_DpsiC= KWs[matID]*
6053  ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
6054  +
6055  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
6056  }
6057  else
6058  {
6059  thetaW = thetaS;
6060  DthetaW_DpsiC = 0.0;
6061  KW = KWs[matID];
6062  DKW_DpsiC = 0.0;
6063  }
6064  //slight compressibility
6065  rhom = rho*exp(beta*u[k]);
6066  drhom = beta*rhom;
6067  mass[k] = rhom*thetaW;
6068  dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
6069  //mass[k] = rho*thetaW;
6070  //dmass[k] = -rho*DthetaW_DpsiC;
6071  for (I=0;I<nSpace;I++)
6072  {
6073  f[k*nSpace+I] = rho2*KW*gravity[I];
6074  df[k*nSpace+I] = -rho2*DKW_DpsiC*gravity[I];
6075  a[k*nSpace2+I*nSpace+I] = rho*KW;
6076  da[k*nSpace2+I*nSpace+I] = -rho*DKW_DpsiC;
6077  }/*I*/
6078  }/*k*/
6079  }/*j*/
6080 
6081 }
6082 
6083 void seepageBrezis(const int nSimplex,
6084  const int nPointsPerSimplex,
6085  const int nSpace,
6086  const int* materialTypes,
6087  const double epsFact,
6088  const double rho,
6089  const double beta,
6090  const double* elementDiameter,
6091  const double* gravity,
6092  const double* alpha,
6093  const double* n,
6094  const double* thetaR,
6095  const double* thetaSR,
6096  const double* KWs,
6097  double *u,
6098  double *mass,
6099  double *dmass,
6100  double *f,
6101  double *df,
6102  double *a,
6103  double *da)
6104 {
6105  int i,j,k,I,matID;
6106  const int nSpace2=nSpace*nSpace;
6107  register double psiC,
6108  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
6109  onePlus_pcBar_n,
6110  sBar,sqrt_sBar,DsBar_DpsiC,
6111  thetaW,DthetaW_DpsiC,
6112  vBar,vBar2,DvBar_DpsiC,
6113  KW,DKW_DpsiC,
6114  rho2=rho*rho,
6115  thetaS,
6116  rhom,drhom,m,eps,hStar,dhStar;
6117  for (i=0; i < nSimplex; i++)
6118  {
6119  matID= materialTypes[i];
6120  eps = epsFact*elementDiameter[i];
6121  for (j=0;j<nPointsPerSimplex;j++)
6122  {
6123  k = i*nPointsPerSimplex + j;
6124  thetaS = thetaR[matID] + thetaSR[matID];
6125  /* hStar = linearHeaviside(eps,u[k]+eps)*(1.0-1.0e-5)+1.0e-5; */
6126  /* dhStar = linearDirac(eps,u[k]+eps)*(1.0-1.0e-5); */
6127  hStar = smoothedHeaviside(eps,u[k]+eps)*(1.0-1.0e-5)+1.0e-5;
6128  dhStar = smoothedDirac(eps,u[k]+eps)*(1.0-1.0e-5);
6129  mass[k] = thetaS*hStar;
6130  dmass[k] = thetaS*dhStar;
6131  for (I=0;I<nSpace;I++)
6132  {
6133  f[k*nSpace+I] = KWs[matID]*gravity[I]*hStar;
6134  df[k*nSpace+I] = KWs[matID]*gravity[I]*dhStar;
6135  a[k*nSpace2+I*nSpace+I] = KWs[matID]*hStar;
6136  da[k*nSpace2+I*nSpace+I] = KWs[matID]*dhStar;
6137  }/*I*/
6138  }/*k*/
6139  }/*j*/
6140 
6141 }
6142 
6143 void conservativeHeadRichardsJLeverett(const int nSimplex,
6144  const int nPointsPerSimplex,
6145  const int nSpace,
6146  const int* materialTypes,
6147  const double rho,
6148  const double beta,
6149  const double* gravity,
6150  const double* phi,
6151  const double* psiD,
6152  const double* ns,
6153  const double* nk,
6154  const double* S_wirr,
6155  const double* S_nwr,
6156  const double* kr0,
6157  double *u,
6158  double *mass,
6159  double *dmass,
6160  double *f,
6161  double *df,
6162  double *a,
6163  double *da)
6164 {
6165  int i,j,k,I,matID;
6166  const int nSpace2=nSpace*nSpace;
6167  register double psiC,
6168  Se,Sw,dSw,
6169  kr,dkr,
6170  rho2=rho*rho,
6171  rhom,drhom;
6172  const double reg_diff=1.0e-2;
6173  for (i=0; i < nSimplex; i++)
6174  {
6175  matID= materialTypes[i];
6176  for (j=0;j<nPointsPerSimplex;j++)
6177  {
6178  k = i*nPointsPerSimplex + j;
6179  psiC = -u[k];
6180  if (psiC > 0.0 && psiC < psiD[matID])
6181  {
6182  Se = 1.0 - pow(psiC/psiD[matID],ns[matID]);
6183  Sw = Se*(1.0 - S_wirr[matID] - S_nwr[matID]) + S_wirr[matID];
6184  dSw = (1.0 - S_wirr[matID] - S_nwr[matID])*pow(psiC/psiD[matID],ns[matID]-1.0)*ns[matID]/psiD[matID];
6185  kr = kr0[matID]*pow(Se,nk[matID])+reg_diff;
6186  dkr = kr0[matID]*pow(Se,nk[matID]-1.0)*nk[matID]*pow(psiC/psiD[matID],ns[matID]-1.0)*ns[matID]/psiD[matID];
6187  }
6188  else if (psiC <= 0.0)
6189  {
6190  Se = 1.0;
6191  Sw = 1.0 - S_nwr[matID];
6192  dSw = 0.0;
6193  kr = kr0[matID]+reg_diff;
6194  dkr = 0.0;
6195  }
6196  else
6197  {
6198  Se = 0.0;
6199  Sw = S_wirr[matID];
6200  dSw = 0.0;
6201  kr = reg_diff;
6202  dkr = 0.0;
6203  }
6204  //slight compressibility
6205  rhom = rho*exp(beta*u[k]);
6206  drhom = beta*rhom;
6207 
6208  mass[k] = rhom*Sw*phi[matID];
6209  dmass[k] = rhom*dSw*phi[matID]+drhom*Sw*phi[matID];
6210  for (I=0;I<nSpace;I++)
6211  {
6212  f[k*nSpace+I] = rho2*kr*gravity[I];
6213  df[k*nSpace+I] = rho2*dkr*gravity[I];
6214  a[k*nSpace2+I*nSpace+I] = rho*kr;
6215  da[k*nSpace2+I*nSpace+I] = rho*dkr;
6216  }/*I*/
6217  }/*k*/
6218  }/*j*/
6219 }
6220 void conservativeHeadRichardsJLeverettAni(const int nSimplex,
6221  const int nPointsPerSimplex,
6222  const int nSpace,
6223  const int* materialTypes,
6224  const double rho,
6225  const double beta,
6226  const double* gravity,
6227  const double* phi,
6228  const double* psiD,
6229  const double* ns,
6230  const double* nk,
6231  const double* S_wirr,
6232  const double* S_nwr,
6233  const double* kr0x,
6234  const double* kr0y,
6235  const double* kr0z,
6236  double *u,
6237  double *mass,
6238  double *dmass,
6239  double *f,
6240  double *df,
6241  double *a,
6242  double *da)
6243 {
6244  int i,j,k,I,matID;
6245  const int nSpace2=nSpace*nSpace;
6246  register double psiC,pcBar,
6247  Se,dSe,Sw,dSw,
6248  kr,dkr,
6249  rho2=rho*rho,
6250  rhom,drhom;
6251  const double reg_diff=1.0e-2;
6252  for (i=0; i < nSimplex; i++)
6253  {
6254  matID= materialTypes[i];
6255  for (j=0;j<nPointsPerSimplex;j++)
6256  {
6257  k = i*nPointsPerSimplex + j;
6258  psiC = -u[k];
6259  if (psiC >= psiD[matID])
6260  {
6261  pcBar = psiC/psiD[matID];
6262  Se = pow(pcBar, -ns[matID]);
6263  Sw = Se*(1.0 - S_wirr[matID] - S_nwr[matID]) + S_wirr[matID];
6264  dSe = ns[matID]*Se/(pcBar*psiD[matID]);
6265  dSw = (1.0 - S_wirr[matID] - S_nwr[matID])*dSe;
6266  kr = pow(Se, (2.0+3.0*ns[matID])/ns[matID]);
6267  dkr = ((2.0+3.0*ns[matID])/ns[matID]*pow(Se, 2.0/ns[matID]*(1.0+ns[matID]))*dSe);
6268  kr = pow(Se,nk[matID]);
6269  dkr = nk[matID]*pow(Se,nk[matID]-1)*dSe;
6270  }
6271  else
6272  {
6273  Se = 1.0;
6274  Sw = 1.0 - S_nwr[matID];
6275  dSw = 0.0;
6276  kr = 1.0;
6277  dkr = 0.0;
6278  }
6279 /* if (psiC > 0.0 && psiC < psiD[matID]) */
6280 /* { */
6281 /* Se = 1.0 - pow(psiC/psiD[matID],ns[matID]); */
6282 /* Sw = Se*(1.0 - S_wirr[matID] - S_nwr[matID]) + S_wirr[matID]; */
6283 /* dSw = (1.0 - S_wirr[matID] - S_nwr[matID])*pow(psiC/psiD[matID],ns[matID]-1.0)*ns[matID]/psiD[matID]; */
6284 /* kr = pow(Se,nk[matID])+reg_diff; */
6285 /* dkr = pow(Se,nk[matID]-1.0)*nk[matID]*pow(psiC/psiD[matID],ns[matID]-1.0)*ns[matID]/psiD[matID]; */
6286 /* } */
6287 /* else if (psiC <= 0.0) */
6288 /* { */
6289 /* Se = 1.0; */
6290 /* Sw = 1.0 - S_nwr[matID]; */
6291 /* dSw = 0.0; */
6292 /* kr = 1.0+reg_diff; */
6293 /* dkr = 0.0; */
6294 /* } */
6295 /* else */
6296 /* { */
6297 /* Se = 0.0; */
6298 /* Sw = S_wirr[matID]; */
6299 /* dSw = 0.0; */
6300 /* kr = reg_diff; */
6301 /* dkr = 0.0; */
6302 /* } */
6303  //slight compressibility
6304  rhom = rho*exp(beta*u[k]);
6305  drhom = beta*rhom;
6306 
6307  mass[k] = rhom*Sw*phi[matID];
6308  dmass[k] = rhom*dSw*phi[matID]+drhom*Sw*phi[matID];
6309  I=0;
6310  f[k*nSpace+I] = rho2*kr0x[matID]*kr*gravity[I];
6311  df[k*nSpace+I] = rho2*kr0x[matID]*dkr*gravity[I];
6312  a[k*nSpace2+I*nSpace+I] = rho*kr0x[matID]*kr;
6313  da[k*nSpace2+I*nSpace+I] = rho*kr0x[matID]*dkr;
6314  if (nSpace > 1)
6315  {
6316  I=1;
6317  f[k*nSpace+I] = rho2*kr0y[matID]*kr*gravity[I];
6318  df[k*nSpace+I] = rho2*kr0y[matID]*dkr*gravity[I];
6319  a[k*nSpace2+I*nSpace+I] = rho*kr0y[matID]*kr;
6320  da[k*nSpace2+I*nSpace+I] = rho*kr0y[matID]*dkr;
6321  if (nSpace > 2)
6322  {
6323  I=2;
6324  f[k*nSpace+I] = rho2*kr0z[matID]*kr*gravity[I];
6325  df[k*nSpace+I] = rho2*kr0z[matID]*dkr*gravity[I];
6326  a[k*nSpace2+I*nSpace+I] = rho*kr0z[matID]*kr;
6327  da[k*nSpace2+I*nSpace+I] = rho*kr0z[matID]*dkr;
6328  }
6329  }
6330  }/*k*/
6331  }/*j*/
6332 }
6333 
6334 
6335 /* mwf end unnecessary additions */
6336 
6337 void conservativeHeadRichardsMualemVanGenuchtenHetEvaluate(const int nPoints,
6338  const int nSpace,
6339  const double rho,
6340  const double* gravity,
6341  const double* alpha,
6342  const double* n,
6343  const double* thetaR,
6344  const double* thetaSR,
6345  const double* KWs,
6346  double *u,
6347  double *mass,
6348  double *dmass,
6349  double *f,
6350  double *df,
6351  double *a,
6352  double *da)
6353 {
6354  int k,I;
6355  const int nSpace2=nSpace*nSpace;
6356  register double psiC,
6357  pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
6358  onePlus_pcBar_n,
6359  sBar,sqrt_sBar,DsBar_DpsiC,
6360  thetaW,DthetaW_DpsiC,
6361  vBar,vBar2,DvBar_DpsiC,
6362  KW,DKW_DpsiC,
6363  rho2=rho*rho,
6364  thetaS,
6365  m;
6366  for (k=0;k<nPoints;k++)
6367  {
6368  psiC = -u[k];
6369  m = 1.0 - 1.0/n[k];
6370  thetaS = thetaR[k] + thetaSR[k];
6371  if (psiC > 0.0)
6372  {
6373  pcBar = alpha[k]*psiC;
6374  pcBar_nM2 = pow(pcBar,n[k]-2);
6375  pcBar_nM1 = pcBar_nM2*pcBar;
6376  pcBar_n = pcBar_nM1*pcBar;
6377  onePlus_pcBar_n = 1.0 + pcBar_n;
6378 
6379  sBar = pow(onePlus_pcBar_n,-m);
6380  /* using -mn = 1-n */
6381  DsBar_DpsiC = alpha[k]*(1.0-n[k])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
6382 
6383  vBar = 1.0-pcBar_nM1*sBar;
6384  vBar2 = vBar*vBar;
6385  DvBar_DpsiC = -alpha[k]*(n[k]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
6386 
6387  thetaW = thetaSR[k]*sBar + thetaR[k];
6388  DthetaW_DpsiC = thetaSR[k] * DsBar_DpsiC;
6389 
6390  sqrt_sBar = sqrt(sBar);
6391  KW= KWs[k]*sqrt_sBar*vBar2;
6392  DKW_DpsiC= KWs[k]*
6393  ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
6394  +
6395  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
6396  }
6397  else
6398  {
6399  thetaW = thetaS;
6400  DthetaW_DpsiC = 0.0;
6401  KW = KWs[k];
6402  DKW_DpsiC = 0.0;
6403  }
6404  mass[k] = rho*thetaW;
6405  dmass[k] = -rho*DthetaW_DpsiC;
6406  for (I=0;I<nSpace;I++)
6407  {
6408  f[k*nSpace+I] = rho2*KW*gravity[I];
6409  df[k*nSpace+I] = -rho2*DKW_DpsiC*gravity[I];
6410&#