17 Length = phiEnd - phiStart;
18 x = (
phi - phiStart)/Length;
19 H = 1 - (exp(pow(x,3.5)) - 1.)/ (exp(1) - 1.);
23 Length = -(phiEnd - phiStart);
24 x = 1 - (
phi - phiStart)/Length;
25 H = 1 - (exp(pow(x,3.5)) - 1.)/ (exp(1) - 1.);
43 H = 0.5*(1.0 +
phi/eps + sin(M_PI*
phi/eps)/M_PI);
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));
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));
73 d = 0.5*(1.0 + cos(M_PI*
phi/eps))/eps;
85 H = 0.5*((
phi+eps)/eps);
119 const int nSpace2=nSpace*nSpace;
120 for (k=0;k<nPoints;k++)
125 for (I=0;I<nSpace;I++)
127 f[k*nSpace+I]=
B[I]*
u[k];
129 for (J=0;J<nSpace;J++)
131 a[k*nSpace2+I*nSpace+J]=A[I*nSpace+J];
144 const double alpha_L,
145 const double alpha_T,
155 const int nSpace2=nSpace*nSpace;
158 for (k=0;k<nPoints;k++)
163 for (I=0;I<nSpace;I++)
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];
169 norm_v = sqrt(norm_v);
172 for (I=0;I<nSpace;I++)
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++)
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];
183 for (I=0;I<nSpace;I++)
184 a[k*nSpace2+I*nSpace+I]=omega*d;
192 const double alpha_L,
193 const double alpha_T,
194 const double Kox_max,
230 const int nSpace2=nSpace*nSpace;
232 double C,E,X,denomC,denomE,denomX,rox,drox_dC,drox_dE,drox_dX;
233 for (k=0;k<nPoints;k++)
235 C = c_c[k]; E = c_e[k]; X = c_x[k];
244 for (I=0;I<nSpace;I++)
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];
251 norm_v +=
v[k*nSpace+I]*
v[k*nSpace+I];
253 norm_v = sqrt(norm_v);
256 for (I=0;I<nSpace;I++)
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++)
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];
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++)
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];
274 for (I=0;I<nSpace;I++)
276 a_c[k*nSpace2+I*nSpace+I]=omega*d_c;
277 a_e[k*nSpace2+I*nSpace+I]=omega*d_e;
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));
288 r_e[k] = 3.0*omega*rox;
289 r_x[k] = -Yield*omega*rox + X*omega*k_d;
291 dr_c_dc[k] = omega*drox_dC; dr_c_de[k] = omega*drox_dE; dr_c_dx[k] = omega*drox_dX;
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;
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;
308 const double alpha_L,
309 const double alpha_T,
340 const int nSpace2=nSpace*nSpace;
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++)
347 C_m = c_m[k]; C_h= c_h[k]; C_oh = K_w/(C_h+eps);
349 C_a = C_h - C_oh; dC_a_dC_h = 1.0 + K_w/(C_h*C_h+eps);
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;
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);
364 dphi_h[k] = dC_a_dC_h;
366 for (I=0;I<nSpace;I++)
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;
373 norm_v +=
v[k*nSpace+I]*
v[k*nSpace+I];
375 norm_v = sqrt(norm_v);
378 for (I=0;I<nSpace;I++)
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++)
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];
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++)
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];
396 for (I=0;I<nSpace;I++)
398 a_m[k*nSpace2+I*nSpace+I]=omega*d_m;
399 a_h[k*nSpace2+I*nSpace+I]=omega*d_h;
414 const int nPointsPerSimplex,
417 const int* materialTypes,
418 const double *omega_types,
419 const double *alpha_L_types,
420 const double *alpha_T_types,
430 const int nSpace2=nSpace*nSpace;
433 for (i=0;i<nSimplex;i++)
435 matID = materialTypes[i];
436 for (j=0; j < nPointsPerSimplex; j++)
438 k = i*nPointsPerSimplex+j;
440 m[k]=omega_types[matID]*
u[k];
441 dm[k]=omega_types[matID];
443 for (I=0;I<nSpace;I++)
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];
449 norm_v = sqrt(norm_v);
452 for (I=0;I<nSpace;I++)
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++)
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];
463 for (I=0;I<nSpace;I++)
464 a[k*nSpace2+I*nSpace+I]=omega_types[matID]*d;
469 const int nPointsPerSimplex,
472 const int* materialTypes,
474 const double *alpha_L_types,
475 const double *alpha_T_types,
485 const int nSpace2=nSpace*nSpace;
488 for (i=0;i<nSimplex;i++)
490 matID = materialTypes[i];
491 for (j=0; j < nPointsPerSimplex; j++)
493 k = i*nPointsPerSimplex+j;
498 for (I=0;I<nSpace;I++)
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];
504 norm_v = sqrt(norm_v);
507 for (I=0;I<nSpace;I++)
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++)
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];
518 for (I=0;I<nSpace;I++)
519 a[k*nSpace2+I*nSpace+I]=theta[k]*d;
524 const int nPointsPerSimplex,
528 const double specificHeat_w,
529 const double specificHeat_n,
530 const int* materialTypes,
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,
549 const int nSpace2=nSpace*nSpace;
550 double norm_v,tmp,sw,Ke,lambda;
552 for (i=0;i<nSimplex;i++)
554 matID = materialTypes[i];
555 for (j=0; j < nPointsPerSimplex; j++)
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];
563 for (I=0;I<nSpace;I++)
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];
569 norm_v = sqrt(norm_v);
570 sw = theta[k]/thetaS_types[matID];
573 for (I=0;I<nSpace;I++)
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++)
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];
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++)
590 a[k*nSpace2+I*nSpace+I] += lambda*lambda_aniso_types[matID*nSpace+I];
623 const int nSpace2=nSpace*nSpace;
629 const double pM1_pow=p_pow-1.0,
634 for (k=0;k<nPoints;k++)
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;
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;
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;
646 for (I=0;I<nSpace;I++)
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++)
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;
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;
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;
695 const int nSpace2=nSpace*nSpace;
696 double max_1mu_0,atmp,datmp;
697 const double p2M1=p2-1.0,
720 for (k=0; k < nPoints; k++)
722 max_1mu_0 = 1.0-
u[k] > 0.0 ? 1.0-
u[k] : 0.0;
726 m[k] *= pow(max_1mu_0,p2);
727 dm[k] *= pow(max_1mu_0,p2M1)*p2;
732 for (I=0; I < nSpace; I++)
734 f[k*nSpace+I] *= pow(max_1mu_0,q2);
735 df[k*nSpace+I] *= pow(max_1mu_0,q2M1)*q2;
741 atmp = pow(max_1mu_0,t2);
742 datmp = pow(max_1mu_0,t2M1);
743 for (I=0; I < nSpace; I++)
745 for (J=0; J < nSpace; J++)
747 a[k*nSpace2+I*nSpace+J] *= atmp;
748 da[k*nSpace2+I*nSpace+J] *= datmp*t2;
755 phi[k] *= pow(max_1mu_0,r2);
756 dphi[k] *= pow(max_1mu_0,r2M1)*r2;
761 r[k] *= pow(max_1mu_0,s2*
r[k]);
762 dr[k] *= pow(max_1mu_0,s2M1)*s2;
778 for (k=0; k < nPoints; k++)
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];
802 for (k=0; k < nPoints; k++)
806 vx = 2.0*M_PI*(x[k*3+1] - 0.5);
807 vy = 2.0*M_PI*(0.5 - x[k*3]);
809 f[k*nSpace] =
vx*
u[k];
810 f[k*nSpace+1] =
vy*
u[k];
811 f[k*nSpace+2] = vz*
u[k];
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));
841 for (k=0; k < nPoints; k++)
843 dm[k] = dphi[k] = 1.0;
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];
852 for (I=0; I < nSpace; I++)
854 a[k*nSpace2 + I*nSpace + I] = self_a;
875 const int nSpace2 = nSpace*nSpace;
877 memcpy(m,
u, nPoints *
sizeof(
double));
878 memcpy(
phi,
u, nPoints *
sizeof(
double));
879 memset(da, 0, nPoints * nSpace2 *
sizeof(
double));
881 for (k=0; k < nPoints; k++)
883 dm[k] = dphi[k] = 1.0;
889 f[k*nSpace] =
vx*
u[k];
890 f[k*nSpace+1] =
vy*
u[k];
893 if (sqrt(X*X+
Y*
Y) < 0.25)
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;
901 for (I=0; I < nSpace; I++)
903 a[k*nSpace2 + I*nSpace + I] = self_a;
923 const int nSpace2 = nSpace*nSpace;
925 for (k=0; k < nPoints; k++)
929 dm[k] = dphi[k] = 1.0;
934 df[k*nSpace+1] = 0.0;
938 df[k*nSpace] *= 0.25;
941 for (I=0; I < nSpace; I++)
943 a[k*nSpace2 + I*nSpace + I] = self_a;
945 for (J=0; J < nSpace; J++)
947 da[k*nSpace2 + I*nSpace + J] = 0.0;
956 const double *self_v,
968 const int nSpace2 = nSpace*nSpace;
970 for (k=0; k < nPoints; k++)
972 m[k] =
phi[k] =
u[k];
973 dm[k] = dphi[k] = 1.0;
975 for (I=0; I < nSpace; I++)
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];
992 const double *self_v,
994 const double *grad_u,
1005 const int nSpace2 = nSpace*nSpace;
1007 for (k=0; k < nPoints; k++)
1009 m[k] =
phi[k] =
u[k];
1010 dm[k] = dphi[k] = 1.0;
1013 for (I=0; I < nSpace; I++)
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];
1049 const int nSpace2=nSpace*nSpace;
1050 for (k=0;k<nPoints;k++)
1054 for (I=0;I<nSpace;I++)
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++)
1060 a[k*nSpace2+I*nSpace+J]=A[k*nSpace2 + I*nSpace+J];
1061 da[k*nSpace2+I*nSpace+J]=0.0;
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)
1084 const int nSpace2=nSpace*nSpace;
1086 for (k=0;k<nPoints;k++)
1088 if (u_levelSet[k] > eps)
1091 for (I=0;I<nSpace;I++)
1093 B[k*nSpace+I]=B1[I];
1094 Bcon[k*nSpace+I]=Bcon1[I];
1095 for (J=0;J<nSpace;J++)
1097 A[k*nSpace2+I*nSpace+J]=A1[I*nSpace+J];
1102 else if (u_levelSet[k] < -eps)
1105 for (I=0;I<nSpace;I++)
1107 B[k*nSpace+I]=B2[I];
1108 Bcon[k*nSpace+I]=Bcon2[I];
1109 for (J=0;J<nSpace;J++)
1111 A[k*nSpace2+I*nSpace+J]=A2[I*nSpace+J];
1118 H = 0.5*(1.0 + u_levelSet[k]/eps + sin((M_PI*u_levelSet[k])/eps)/M_PI);
1120 M[k] = oneMinusH*M2 +
H*M1;
1121 for (I=0;I<nSpace;I++)
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++)
1127 A[k*nSpace2+I*nSpace+J]=oneMinusH*A2[I*nSpace+J] +
H*A1[I*nSpace+J];
1130 C[k]=oneMinusH*C2+
H*C1;
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];
1154 double* m,
double* dm,
1155 double* h,
double* dh,
1159 for (i=0;i<nPoints;i++)
1165 for (I=0;I<nSpace;I++)
1167 h[i] +=
B[i*nSpace+I]*grad_u[i*nSpace+I];
1168 dh[i*nSpace+I]=
B[i*nSpace+I];
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)
1186 for (i=0;i<nPoints;i++)
1190 for (I=0;I<nSpace;I++)
1192 f[i*nSpace+I] =
B[i*nSpace+I]*
u[i];
1193 df[i*nSpace+I]=
B[i*nSpace+I];
1209 for (i=0;i<nPoints;i++)
1214 for (I=0;I<nSpace;I++)
1216 H[i] +=
v[i*nSpace+I]*grad_u[i*nSpace+I];
1217 dH[i*nSpace+I] =
v[i*nSpace+I];
1232 for (i=0;i<nPoints;i++)
1236 for (I=0;I<nSpace;I++)
1238 f[i*nSpace+I] =
v[i*nSpace+I]*
u[i];
1239 df[i*nSpace+I] =
v[i*nSpace+I];
1256 for (i=0;i<nPoints;i++)
1260 for (I=0;I<nSpace;I++)
1264 f[i*nSpace+I] =
v[i*nSpace+I]*
u[i];
1265 df[i*nSpace+I] =
v[i*nSpace+I];
1279 double norm_grad_phi;
1280 for (i=0;i<nPoints;i++)
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)
1290 for (I=0;I<nSpace;I++)
1291 f[i*nSpace+I] = grad_phi[i*nSpace+I]/norm_grad_phi;
1294 f[i*nSpace+I] = 0.0;
1305 for (k=0;k<nPoints;k++)
1307 if (u_levelSet[k] > eps)
1309 else if (u_levelSet[k] < -eps)
1313 H = 0.5*(1.0 + u_levelSet[k]/eps + sin((M_PI*u_levelSet[k])/eps)/M_PI);
1330 double* m,
double* dm,
1331 double* h,
double* dh,
1335 for (i=0;i<nPoints;i++)
1341 for (I=0;I<nSpace;I++)
1342 h[i] += grad_u[i*nSpace+I]*grad_u[i*nSpace+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];
1348 dh[i*nSpace+I] = (S[i]*grad_u[i*nSpace+I])/1.0e-8;
1367 for (i=0;i<nPoints;i++)
1374 for (I=0;I<nSpace;I++)
1375 normGradU+= grad_u[i*nSpace+I]*grad_u[i*nSpace+I];
1376 normGradU = sqrt(normGradU);
1378 for (I=0;I<nSpace;I++)
1380 dH[i*nSpace+I] = grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1397 double Si,normGradU;
1401 for (i=0;i<nPoints;i++)
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;
1437 for (I=0;I<nSpace;I++)
1439 dH[i*nSpace+I] = Si*grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1446 double lambda_penalty,
1458 double Si,normGradU;
1462 for (i=0;i<nPoints;i++)
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;
1498 for (I=0;I<nSpace;I++)
1500 dH[i*nSpace+I] = Si*grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1503 r[i] += (
u[i]-u_levelSet[i])*lambda_penalty*
smoothedDirac(eps,u_levelSet[i]);
1513 int nPointsPerSimplex,
1527 double He,Si,normGradU,
1533 for (ie=0; ie < nSimplex; ie++)
1538 for (iq=0;iq<nPointsPerSimplex;iq++)
1540 i = ie*nPointsPerSimplex + iq;
1545 if (u_levelSet[i] > eps)
1549 else if (u_levelSet[i] < -eps)
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);
1560 for (I=0;I<nSpace;I++)
1561 normGradU+= grad_u[i*nSpace+I]*grad_u[i*nSpace+I];
1562 normGradU = sqrt(normGradU);
1565 H[i] = Si*normGradU;
1566 for (I=0;I<nSpace;I++)
1568 dH[i*nSpace+I] = Si*grad_u[i*nSpace+I]/(normGradU+1.0e-12);
1573 for (iq = 0; iq < nPointsPerSimplex; iq++)
1575 i = ie*nPointsPerSimplex + iq;
1578 if (u_levelSet[i] > eps)
1582 else if (u_levelSet[i] < -eps)
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);
1592 normGradU =
H[i]/Si;
1593 lambda -= dHe*Li*dV[i]/(dHe*dHe*normGradU+1.0e-8);
1596 for (iq = 0; iq < nPointsPerSimplex; iq++)
1598 i = ie*nPointsPerSimplex + iq;
1600 if (u_levelSet[i] > eps)
1602 else if (u_levelSet[i] < -eps)
1603 {Si=-1.0; He = 0.0;}
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);
1610 normGradU =
H[i]/Si;
1611 r[i] -= lambda*dHe*normGradU;
1621 int nDOF_trial_element,
1622 double epsilon_freeze_factor,
1623 const double *elementDiameter,
1625 const double *u_dof,
1626 int * freeze_nodes_tmp,
1627 int * weakDirichletConditionFlags)
1629 int eN,j,jj,signU,j0,J0,J;
1631 for (eN = 0; eN < nElements_global; eN++)
1633 for (j = 0; j < nDOF_trial_element; j++)
1634 freeze_nodes_tmp[j] = 0;
1635 eps = epsilon_freeze_factor*elementDiameter[eN];
1637 while (signU == 0 && j0 < nDOF_trial_element)
1639 J0 = u_l2g[eN*nDOF_trial_element + j0];
1640 if (u_dof[J0] < -eps)
1642 else if (u_dof[J0] > eps)
1645 freeze_nodes_tmp[j0] = 1;
1648 for (j = j0; j < nDOF_trial_element; j++)
1650 J = u_l2g[eN*nDOF_trial_element + j];
1651 if ((u_dof[J] < -eps && signU == 1) ||
1652 (u_dof[J] > eps && signU == -1))
1654 for (jj = 0; jj < nDOF_trial_element; jj++)
1655 freeze_nodes_tmp[jj] = 1;
1658 else if (fabs(u_dof[J]) < eps)
1659 freeze_nodes_tmp[j] = 1;
1661 for (j = 0; j < nDOF_trial_element; j++)
1663 if (freeze_nodes_tmp[j] == 1)
1665 J = u_l2g[eN*nDOF_trial_element + j];
1666 weakDirichletConditionFlags[J] = 1;
1673 int nDOF_trial_element,
1674 double epsilon_freeze_factor,
1675 const double *elementDiameter,
1677 const double *u_dof,
1678 int * freeze_nodes_tmp,
1679 int * weakDirichletConditionFlags)
1683 for (eN = 0; eN < nElements_global; eN++)
1685 eps = epsilon_freeze_factor*elementDiameter[eN];
1686 for (j = 0; j < nDOF_trial_element; j++)
1688 J = u_l2g[eN*nDOF_trial_element + j];
1689 if (fabs(u_dof[J]) < eps)
1690 weakDirichletConditionFlags[J] = 1;
1700 double Km,
double rhoM,
1701 double Kp,
double rhoP,
1706 double * u_levelSet,
1721 nSpace2 = nSpace*nSpace;
1722 for (k = 0; k < nPoints; k++)
1724 m[k] = 0.0; dm[k] = 1.0;
1725 r[k] = 0.0; dr[k] = 0.0;
1733 for (I = 0; I < nSpace; I++)
1735 for (J = 0; J < nSpace; J++)
1738 a[k*nSpace2 + I*nSpace+J] = Km + He*(Kp-Km);
1740 a[k*nSpace2 + I*nSpace+J] = 0.0;
1741 da[k*nSpace2 + I*nSpace+J] = 0.0;
1743 f[k*nSpace + I] = (Km*rhoM + He*(Kp*rhoP-Km*rhoM))*gravity_u[I];
1744 df[k*nSpace + I]= 0.0;
1752 double Km,
double rhoM,
1753 double Kp,
double rhoP,
1758 double * u_levelSet,
1772 double rhoEps = 0.0;
1773 nSpace2 = nSpace*nSpace;
1774 for (k = 0; k < nPoints; k++)
1776 m[k] = 0.0; dm[k] = 1.0;
1777 r[k] = 0.0; dr[k] = 0.0;
1786 for (I = 0; I < nSpace; I++)
1788 for (J = 0; J < nSpace; J++)
1791 a[k*nSpace2 + I*nSpace+J] = Km + He*(Kp-Km);
1793 a[k*nSpace2 + I*nSpace+J] = 0.0;
1794 da[k*nSpace2 + I*nSpace+J] = 0.0;
1796 f[k*nSpace + I] = (Km*rhoM + He*(Kp*rhoP-Km*rhoM))*gravity_u[I];
1797 df[k*nSpace + I]= 0.0;
1804 double *mom_p_diff_ten,
1805 double *mom_u_diff_ten,
1806 double *mom_v_diff_ten)
1809 for (k=0; k<nPoints; k++)
1811 mom_p_diff_ten[k*2+0] = 1.0;
1812 mom_p_diff_ten[k*2+1] = 1.0;
1814 mom_u_diff_ten[k*2+0] = 1.0;
1815 mom_u_diff_ten[k*2+1] = 1.0;
1817 mom_v_diff_ten[k*2+0] = 1.0;
1818 mom_v_diff_ten[k*2+1] = 1.0;
1823 double *mom_p_diff_ten,
1824 double *mom_u_diff_ten,
1825 double *mom_v_diff_ten,
1826 double *mom_w_diff_ten)
1829 for (k=0; k<nPoints; k++)
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;
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;
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;
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;
1855 const double *grad_p,
1859 double *dmom_u_acc_u,
1861 double *dmom_v_acc_v,
1863 double *dmass_adv_u,
1864 double *dmass_adv_v,
1866 double *dmom_u_adv_u,
1867 double *dmom_u_adv_v,
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,
1876 double *dmom_u_ham_grad_p,
1878 double *dmom_v_ham_grad_p)
1881 for (k=0;k<nPoints;k++)
1886 dmom_u_acc_u[k]=1.0;
1889 dmom_v_acc_v[k]=1.0;
1892 mass_adv[k*2+0]=
u[k];
1893 mass_adv[k*2+1]=
v[k];
1895 dmass_adv_u[k*2+0]=1.0;
1896 dmass_adv_v[k*2+1]=1.0;
1899 mom_u_adv[k*2+0]=
u[k]*
u[k];
1900 mom_u_adv[k*2+1]=
u[k]*
v[k];
1902 dmom_u_adv_u[k*2+0]=2.0*
u[k];
1903 dmom_u_adv_u[k*2+1]=
v[k];
1905 dmom_u_adv_v[k*2+1]=
u[k];
1908 mom_v_adv[k*2+0]=
v[k]*
u[k];
1909 mom_v_adv[k*2+1]=
v[k]*
v[k];
1911 dmom_v_adv_u[k*2+0]=
v[k];
1913 dmom_v_adv_v[k*2+0]=
u[k];
1914 dmom_v_adv_v[k*2+1]=2.0*
v[k];
1917 mom_u_diff_ten[k*4+0] = nu;
1918 mom_u_diff_ten[k*4+3] = nu;
1921 mom_v_diff_ten[k*4+0] = nu;
1922 mom_v_diff_ten[k*4+3] = nu;
1925 mom_u_source[k] = -g[0];
1926 mom_v_source[k] = -g[1];
1929 mom_u_ham[k] = grad_p[k*2+0]/rho;
1930 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
1933 mom_v_ham[k] = grad_p[k*2+1]/rho;
1934 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
1943 const double *grad_p,
1948 double *dmom_u_acc_u,
1950 double *dmom_v_acc_v,
1952 double *dmom_w_acc_w,
1954 double *dmass_adv_u,
1955 double *dmass_adv_v,
1956 double *dmass_adv_w,
1958 double *dmom_u_adv_u,
1959 double *dmom_u_adv_v,
1960 double *dmom_u_adv_w,
1962 double *dmom_v_adv_u,
1963 double *dmom_v_adv_v,
1964 double *dmom_v_adv_w,
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,
1976 double *dmom_u_ham_grad_p,
1978 double *dmom_v_ham_grad_p,
1980 double *dmom_w_ham_grad_p)
1983 for (k=0;k<nPoints;k++)
1989 dmom_u_acc_u[k]=1.0;
1992 dmom_v_acc_v[k]=1.0;
1995 dmom_w_acc_w[k]=1.0;
1998 mass_adv[k*3+0]=
u[k];
1999 mass_adv[k*3+1]=
v[k];
2000 mass_adv[k*3+2]=
w[k];
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;
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];
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];
2015 dmom_u_adv_v[k*3+1]=
u[k];
2017 dmom_u_adv_w[k*3+2]=
u[k];
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];
2024 dmom_v_adv_u[k*3+0]=
v[k];
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];
2030 dmom_v_adv_w[k*3+2]=
v[k];
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];
2037 dmom_w_adv_u[k*3+0]=
w[k];
2039 dmom_w_adv_v[k*3+0]=
w[k];
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];
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;
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;
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;
2061 mom_u_source[k] = -g[0];
2062 mom_v_source[k] = -g[1];
2063 mom_w_source[k] = -g[2];
2066 mom_u_ham[k] = grad_p[k*3+0]/rho;
2067 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
2070 mom_v_ham[k] = grad_p[k*3+1]/rho;
2071 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
2074 mom_w_ham[k] = grad_p[k*3+2]/rho;
2075 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
2084 const double *grad_p,
2088 double *dmom_u_acc_u,
2090 double *dmom_v_acc_v,
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,
2099 double *dmom_u_ham_grad_p,
2101 double *dmom_v_ham_grad_p)
2104 for (k=0;k<nPoints;k++)
2109 dmom_u_acc_u[k]=1.0;
2112 dmom_v_acc_v[k]=1.0;
2115 mass_adv[k*2+0]=
u[k];
2116 mass_adv[k*2+1]=
v[k];
2118 dmass_adv_u[k*2+0]=1.0;
2119 dmass_adv_v[k*2+1]=1.0;
2122 mom_u_diff_ten[k*4+0] = nu;
2123 mom_u_diff_ten[k*4+3] = nu;
2126 mom_v_diff_ten[k*4+0] = nu;
2127 mom_v_diff_ten[k*4+3] = nu;
2130 mom_u_source[k] = -g[0];
2131 mom_v_source[k] = -g[1];
2134 mom_u_ham[k] = grad_p[k*2+0]/rho;
2135 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2138 mom_v_ham[k] = grad_p[k*2+1]/rho;
2139 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2151 double *dmom_u_acc_u,
2153 double *dmom_v_acc_v,
2155 double *dmass_adv_u,
2156 double *dmass_adv_v,
2158 double *dmom_u_adv_p,
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)
2167 for (k=0;k<nPoints;k++)
2172 dmom_u_acc_u[k]=1.0;
2175 dmom_v_acc_v[k]=1.0;
2178 mass_adv[k*2+0]=
u[k];
2179 mass_adv[k*2+1]=
v[k];
2181 dmass_adv_u[k*2+0]=1.0;
2182 dmass_adv_v[k*2+1]=1.0;
2185 mom_u_adv[k*2+0]=p[k]/rho;
2186 dmom_u_adv_p[k*2+0]=1.0/rho;
2189 mom_v_adv[k*2+1]=p[k]/rho;
2190 dmom_v_adv_p[k*2+1]=1.0/rho;
2193 mom_u_diff_ten[k*4+0] = nu;
2194 mom_u_diff_ten[k*4+3] = nu;
2197 mom_v_diff_ten[k*4+0] = nu;
2198 mom_v_diff_ten[k*4+3] = nu;
2201 mom_u_source[k] = -g[0];
2202 mom_v_source[k] = -g[1];
2211 const double *grad_p,
2216 double *dmom_u_acc_u,
2218 double *dmom_v_acc_v,
2220 double *dmom_w_acc_w,
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,
2232 double *dmom_u_ham_grad_p,
2234 double *dmom_v_ham_grad_p,
2236 double *dmom_w_ham_grad_p)
2239 for (k=0;k<nPoints;k++)
2244 dmom_u_acc_u[k]=1.0;
2247 dmom_v_acc_v[k]=1.0;
2250 dmom_w_acc_w[k]=1.0;
2253 mass_adv[k*3+0]=
u[k];
2254 mass_adv[k*3+1]=
v[k];
2255 mass_adv[k*3+2]=
w[k];
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;
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;
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;
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;
2277 mom_u_source[k] = -g[0];
2278 mom_v_source[k] = -g[1];
2279 mom_w_source[k] = -g[2];
2282 mom_u_ham[k] = grad_p[k*3+0]/rho;
2283 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
2286 mom_v_ham[k] = grad_p[k*3+1]/rho;
2287 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
2290 mom_w_ham[k] = grad_p[k*3+2]/rho;
2291 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
2304 double *dmom_u_acc_u,
2306 double *dmom_v_acc_v,
2308 double *dmom_w_acc_w,
2310 double *dmass_adv_u,
2311 double *dmass_adv_v,
2312 double *dmass_adv_w,
2314 double *dmom_u_adv_p,
2316 double *dmom_v_adv_p,
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)
2327 for (k=0;k<nPoints;k++)
2332 dmom_u_acc_u[k]=1.0;
2335 dmom_v_acc_v[k]=1.0;
2338 dmom_w_acc_w[k]=1.0;
2341 mass_adv[k*3+0]=
u[k];
2342 mass_adv[k*3+1]=
v[k];
2343 mass_adv[k*3+2]=
w[k];
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;
2350 mom_u_adv[k*3+0]=p[k]/rho;
2351 dmom_u_adv_p[k*3+0]=1.0/rho;
2354 mom_v_adv[k*3+1]=p[k]/rho;
2355 dmom_v_adv_p[k*3+1]=1.0/rho;
2358 mom_w_adv[k*3+2]=p[k]/rho;
2359 dmom_w_adv_p[k*3+2]=1.0/rho;
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;
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;
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;
2377 mom_u_source[k] = -g[0];
2378 mom_v_source[k] = -g[1];
2379 mom_w_source[k] = -g[2];
2392 const double *grad_p,
2396 double *dmom_u_acc_u,
2398 double *dmom_v_acc_v,
2400 double *dmass_adv_u,
2401 double *dmass_adv_v,
2403 double *dmom_u_adv_u,
2404 double *dmom_u_adv_v,
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,
2413 double *dmom_u_ham_grad_p,
2415 double *dmom_v_ham_grad_p)
2419 for (k=0;k<nPoints;k++)
2428 dmom_u_acc_u[k]=1.0;
2432 dmom_v_acc_v[k]=1.0;
2435 mass_adv[k*2+0]=
u[k];
2436 mass_adv[k*2+1]=
v[k];
2438 dmass_adv_u[k*2+0]=1.0;
2439 dmass_adv_v[k*2+1]=1.0;
2442 mom_u_adv[k*2+0]=
u[k]*
u[k];
2443 mom_u_adv[k*2+1]=
u[k]*
v[k];
2445 dmom_u_adv_u[k*2+0]=2.0*
u[k];
2446 dmom_u_adv_u[k*2+1]=
v[k];
2448 dmom_u_adv_v[k*2+1]=
u[k];
2451 mom_v_adv[k*2+0]=
v[k]*
u[k];
2452 mom_v_adv[k*2+1]=
v[k]*
v[k];
2454 dmom_v_adv_u[k*2+0]=
v[k];
2456 dmom_v_adv_v[k*2+0]=
u[k];
2457 dmom_v_adv_v[k*2+1]=2.0*
v[k];
2460 mom_u_diff_ten[k*4+0] = nu;
2461 mom_u_diff_ten[k*4+3] = nu;
2464 mom_v_diff_ten[k*4+0] = nu;
2465 mom_v_diff_ten[k*4+3] = nu;
2468 mom_u_source[k] = -g[0];
2469 mom_v_source[k] = -g[1];
2472 mom_u_ham[k] = grad_p[k*2+0]/rho;
2473 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2476 mom_v_ham[k] = grad_p[k*2+1]/rho;
2477 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2482 const double eps_rho,
2483 const double eps_mu,
2492 const double* kappa,
2494 const double *grad_p,
2498 double *dmom_u_acc_u,
2500 double *dmom_v_acc_v,
2502 double *dmass_adv_u,
2503 double *dmass_adv_v,
2505 double *dmom_u_adv_u,
2506 double *dmom_u_adv_v,
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,
2517 double *dmom_u_ham_grad_p,
2519 double *dmom_v_ham_grad_p)
2522 double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
2523 for (k=0;k<nPoints;k++)
2607 dmom_u_acc_u[k]=1.0;
2611 dmom_v_acc_v[k]=1.0;
2614 mass_adv[k*2+0]=
u[k];
2615 mass_adv[k*2+1]=
v[k];
2617 dmass_adv_u[k*2+0]=1.0;
2618 dmass_adv_v[k*2+1]=1.0;
2621 mom_u_adv[k*2+0]=
u[k]*
u[k];
2622 mom_u_adv[k*2+1]=
u[k]*
v[k];
2624 dmom_u_adv_u[k*2+0]=2.0*
u[k];
2625 dmom_u_adv_u[k*2+1]=
v[k];
2627 dmom_u_adv_v[k*2+1]=
u[k];
2630 mom_v_adv[k*2+0]=
v[k]*
u[k];
2631 mom_v_adv[k*2+1]=
v[k]*
v[k];
2633 dmom_v_adv_u[k*2+0]=
v[k];
2635 dmom_v_adv_v[k*2+0]=
u[k];
2636 dmom_v_adv_v[k*2+1]=2.0*
v[k];
2638 #ifdef SCALAR_DIFFUSION
2640 mom_u_diff_ten[k*4+0] = nu;
2641 mom_u_diff_ten[k*4+3] = nu;
2644 mom_v_diff_ten[k*4+0] = nu;
2645 mom_v_diff_ten[k*4+3] = nu;
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;
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;
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));
2667 mom_u_ham[k] = grad_p[k*2+0]/rho;
2668 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2671 mom_v_ham[k] = grad_p[k*2+1]/rho;
2672 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2678 const double eps_rho,
2679 const double eps_mu,
2688 const double* kappa,
2690 const double *grad_p,
2694 double *dmom_u_acc_u,
2696 double *dmom_v_acc_v,
2698 double *dmass_adv_u,
2699 double *dmass_adv_v,
2701 double *dmom_u_adv_u,
2702 double *dmom_u_adv_v,
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,
2713 double *dmom_u_ham_grad_p,
2715 double *dmom_v_ham_grad_p)
2718 double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
2719 for (k=0;k<nPoints;k++)
2733 dmom_u_acc_u[k]=1.0;
2737 dmom_v_acc_v[k]=1.0;
2740 mass_adv[k*2+0]=
u[k];
2741 mass_adv[k*2+1]=
v[k];
2743 dmass_adv_u[k*2+0]=1.0;
2744 dmass_adv_v[k*2+1]=1.0;
2747 mom_u_adv[k*2+0]=
u[k]*
u[k];
2748 mom_u_adv[k*2+1]=
u[k]*
v[k];
2750 dmom_u_adv_u[k*2+0]=2.0*
u[k];
2751 dmom_u_adv_u[k*2+1]=
v[k];
2753 dmom_u_adv_v[k*2+1]=
u[k];
2756 mom_v_adv[k*2+0]=
v[k]*
u[k];
2757 mom_v_adv[k*2+1]=
v[k]*
v[k];
2759 dmom_v_adv_u[k*2+0]=
v[k];
2761 dmom_v_adv_v[k*2+0]=
u[k];
2762 dmom_v_adv_v[k*2+1]=2.0*
v[k];
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;
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;
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));
2786 mom_u_ham[k] = grad_p[k*2+0]/rho;
2787 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
2790 mom_v_ham[k] = grad_p[k*2+1]/rho;
2791 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
2854 const double eps_rho,
2855 const double eps_mu,
2866 const double* kappa,
2867 const double* phi_s,
2870 const double *grad_p,
2874 double *dmom_u_acc_u,
2876 double *dmom_v_acc_v,
2878 double *dmass_adv_u,
2879 double *dmass_adv_v,
2881 double *dmom_u_adv_u,
2882 double *dmom_u_adv_v,
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,
2897 double *dmom_u_ham_grad_p,
2899 double *dmom_v_ham_grad_p)
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++)
2929 mom_u_acc[k]=rho*
u[k];
2930 dmom_u_acc_u[k]=rho;
2933 mom_v_acc[k]=rho*
v[k];
2934 dmom_v_acc_v[k]=rho;
2937 mass_adv[k*2+0]=
u[k];
2938 mass_adv[k*2+1]=
v[k];
2940 dmass_adv_u[k*2+0]=1.0;
2941 dmass_adv_v[k*2+1]=1.0;
2944 mom_u_adv[k*2+0]=rho*
u[k]*
u[k];
2945 mom_u_adv[k*2+1]=rho*
u[k]*
v[k];
2947 dmom_u_adv_u[k*2+0]=2.0*rho*
u[k];
2948 dmom_u_adv_u[k*2+1]=rho*
v[k];
2950 dmom_u_adv_v[k*2+1]=rho*
u[k];
2953 mom_v_adv[k*2+0]=rho*
v[k]*
u[k];
2954 mom_v_adv[k*2+1]=rho*
v[k]*
v[k];
2956 dmom_v_adv_u[k*2+0]=rho*
v[k];
2958 dmom_v_adv_v[k*2+0]=rho*
u[k];
2959 dmom_v_adv_v[k*2+1]=2.0*rho*
v[k];
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;
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;
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]);
3031 mom_u_source[k] = -rho*g[0];
3032 dmom_u_source_u[k] = 0.0;
3033 dmom_u_source_v[k] = 0.0;
3035 mom_v_source[k] = -rho*g[1];
3036 dmom_v_source_u[k] = 0.0;
3037 dmom_v_source_v[k] = 0.0;
3041 mom_u_ham[k] = grad_p[k*2+0];
3042 dmom_u_ham_grad_p[k*2+0]=1.0;
3045 mom_v_ham[k] = grad_p[k*2+1];
3046 dmom_v_ham_grad_p[k*2+1]=1.0;
3051 const double boundaryPenaltyCoef,
3052 const double volumePenaltyCoef,
3053 const double eps_rho,
3054 const double eps_mu,
3065 const double* kappa,
3066 const double* phi_s,
3069 const double *grad_p,
3073 double *dmom_u_acc_u,
3075 double *dmom_v_acc_v,
3077 double *dmass_adv_u,
3078 double *dmass_adv_v,
3080 double *dmom_u_adv_u,
3081 double *dmom_u_adv_v,
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,
3096 double *dmom_u_ham_grad_p,
3098 double *dmom_v_ham_grad_p)
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);
3106 for (k=0;k<nPoints;k++)
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]);
3133 mom_u_acc[k]=H_mu_s*rho*
u[k];
3134 dmom_u_acc_u[k]=H_mu_s*rho;
3137 mom_v_acc[k]=H_mu_s*rho*
v[k];
3138 dmom_v_acc_v[k]=H_mu_s*rho;
3142 mass_adv[k*2+0]=
u[k];
3143 mass_adv[k*2+1]=
v[k];
3145 dmass_adv_u[k*2+0]=1.0;
3146 dmass_adv_v[k*2+1]=1.0;
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;
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;
3156 dmom_u_adv_v[k*2+1]=H_mu_s*rho*
u[k];
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;
3162 dmom_v_adv_u[k*2+0]=H_mu_s*rho*
v[k];
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;
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;
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;
3192 if (quadraticPenalty)
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);
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;
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);
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);
3212 else if (sipgPenalty)
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);
3218 dmom_u_source_u[k] = volumePenaltyCoef*rho*(1.0-H_mu_s);
3219 dmom_u_source_v[k] = 0.0;
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);
3225 dmom_v_source_u[k] = 0.0;
3226 dmom_v_source_v[k] = volumePenaltyCoef*rho*(1.0-H_mu_s);
3228 volumeFlux =
u[k]*n_s[k*2+0] +
v[k]*n_s[k*2+1];
3229 if (volumeFlux < 0.0)
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;
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;
3236 dmom_u_source_v[k] = rho*d_mu_s*
u[k]*n_s[k*2+1]/norm_n_s;
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;
3241 dmom_v_source_u[k] += rho*d_mu_s*
v[k]*n_s[k*2+0]/norm_n_s;
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;
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);
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;
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);
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);
3270 mom_u_ham[k] = grad_p[k*2+0];
3271 dmom_u_ham_grad_p[k*2+0]=1.0;
3274 mom_v_ham[k] = grad_p[k*2+1];
3275 dmom_v_ham_grad_p[k*2+1]=1.0;
3280 const double eps_rho,
3281 const double eps_mu,
3290 const double* kappa,
3292 const double *grad_p,
3297 double *dmom_u_acc_u,
3299 double *dmom_v_acc_v,
3301 double *dmom_w_acc_w,
3303 double *dmass_adv_u,
3304 double *dmass_adv_v,
3305 double *dmass_adv_w,
3307 double *dmom_u_adv_u,
3308 double *dmom_u_adv_v,
3309 double *dmom_u_adv_w,
3311 double *dmom_v_adv_u,
3312 double *dmom_v_adv_v,
3313 double *dmom_v_adv_w,
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,
3331 double *dmom_u_ham_grad_p,
3333 double *dmom_v_ham_grad_p,
3335 double *dmom_w_ham_grad_p)
3338 double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
3339 for (k=0;k<nPoints;k++)
3462 dmom_u_acc_u[k]=1.0;
3466 dmom_v_acc_v[k]=1.0;
3470 dmom_w_acc_w[k]=1.0;
3474 mass_adv[k*3+0]=
u[k];
3475 mass_adv[k*3+1]=
v[k];
3476 mass_adv[k*3+2]=
w[k];
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;
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];
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];
3491 dmom_u_adv_v[k*3+1]=
u[k];
3493 dmom_u_adv_w[k*3+2]=
u[k];
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];
3500 dmom_v_adv_u[k*3+0]=
v[k];
3502 dmom_v_adv_w[k*3+2]=
v[k];
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];
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];
3513 dmom_w_adv_u[k*3+0]=
w[k];
3515 dmom_w_adv_v[k*3+1]=
w[k];
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];
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;
3526 mom_uv_diff_ten[k*9+3]=nu;
3528 mom_uw_diff_ten[k*9+6]=nu;
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;
3535 mom_vu_diff_ten[k*9+1]=nu;
3537 mom_vw_diff_ten[k*9+7]=nu;
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;
3544 mom_wu_diff_ten[k*9+2]=nu;
3546 mom_wv_diff_ten[k*9+5]=nu;
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));
3556 mom_u_ham[k] = grad_p[k*3+0]/rho;
3557 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
3560 mom_v_ham[k] = grad_p[k*3+1]/rho;
3561 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
3564 mom_w_ham[k] = grad_p[k*3+2]/rho;
3565 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
3569 const double eps_rho,
3570 const double eps_mu,
3579 const double* kappa,
3581 const double *grad_p,
3586 double *dmom_u_acc_u,
3588 double *dmom_v_acc_v,
3590 double *dmom_w_acc_w,
3592 double *dmass_adv_u,
3593 double *dmass_adv_v,
3594 double *dmass_adv_w,
3596 double *dmom_u_adv_u,
3597 double *dmom_u_adv_v,
3598 double *dmom_u_adv_w,
3600 double *dmom_v_adv_u,
3601 double *dmom_v_adv_v,
3602 double *dmom_v_adv_w,
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,
3620 double *dmom_u_ham_grad_p,
3622 double *dmom_v_ham_grad_p,
3624 double *dmom_w_ham_grad_p)
3627 double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
3628 for (k=0;k<nPoints;k++)
3643 dmom_u_acc_u[k]=1.0;
3647 dmom_v_acc_v[k]=1.0;
3651 dmom_w_acc_w[k]=1.0;
3655 mass_adv[k*3+0]=
u[k];
3656 mass_adv[k*3+1]=
v[k];
3657 mass_adv[k*3+2]=
w[k];
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;
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];
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];
3672 dmom_u_adv_v[k*3+1]=
u[k];
3674 dmom_u_adv_w[k*3+2]=
u[k];
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];
3681 dmom_v_adv_u[k*3+0]=
v[k];
3683 dmom_v_adv_w[k*3+2]=
v[k];
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];
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];
3694 dmom_w_adv_u[k*3+0]=
w[k];
3696 dmom_w_adv_v[k*3+1]=
w[k];
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];
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;
3707 mom_uv_diff_ten[k]=nu;
3709 mom_uw_diff_ten[k]=nu;
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;
3716 mom_vu_diff_ten[k]=nu;
3718 mom_vw_diff_ten[k]=nu;
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;
3725 mom_wu_diff_ten[k]=nu;
3727 mom_wv_diff_ten[k]=nu;
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));
3737 mom_u_ham[k] = grad_p[k*3+0]/rho;
3738 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
3741 mom_v_ham[k] = grad_p[k*3+1]/rho;
3742 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
3745 mom_w_ham[k] = grad_p[k*3+2]/rho;
3746 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
3750 const double boundaryPenaltyCoef,
3751 const double volumePenaltyCoef,
3752 const double eps_rho,
3753 const double eps_mu,
3764 const double* kappa,
3765 const double* phi_s,
3768 const double *grad_p,
3773 double *dmom_u_acc_u,
3775 double *dmom_v_acc_v,
3777 double *dmom_w_acc_w,
3779 double *dmass_adv_u,
3780 double *dmass_adv_v,
3781 double *dmass_adv_w,
3783 double *dmom_u_adv_u,
3784 double *dmom_u_adv_v,
3785 double *dmom_u_adv_w,
3787 double *dmom_v_adv_u,
3788 double *dmom_v_adv_v,
3789 double *dmom_v_adv_w,
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,
3816 double *dmom_u_ham_grad_p,
3818 double *dmom_v_ham_grad_p,
3820 double *dmom_w_ham_grad_p)
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;
3826 for (k=0;k<nPoints;k++)
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;
3849 mom_u_acc[k]=H_mu_s*rho*
u[k];
3850 dmom_u_acc_u[k]=H_mu_s*rho;
3853 mom_v_acc[k]=H_mu_s*rho*
v[k];
3854 dmom_v_acc_v[k]=H_mu_s*rho;
3857 mom_w_acc[k]=H_mu_s*rho*
w[k];
3858 dmom_w_acc_w[k]=H_mu_s*rho;
3862 mass_adv[k*3+0]=
u[k];
3863 mass_adv[k*3+1]=
v[k];
3864 mass_adv[k*3+2]=
w[k];
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;
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];
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];
3879 dmom_u_adv_v[k*3+1]=H_mu_s*rho*
u[k];
3881 dmom_u_adv_w[k*3+2]=H_mu_s*rho*
u[k];
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];
3888 dmom_v_adv_u[k*3+0]=H_mu_s*rho*
v[k];
3890 dmom_v_adv_w[k*3+2]=H_mu_s*rho*
v[k];
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];
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];
3901 dmom_w_adv_u[k*3+0]=H_mu_s*rho*
w[k];
3903 dmom_w_adv_v[k*3+1]=H_mu_s*rho*
w[k];
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];
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;
3914 mom_uv_diff_ten[k*9+3]=H_mu_s*mu;
3916 mom_uw_diff_ten[k*9+6]=H_mu_s*mu;
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;
3923 mom_vu_diff_ten[k*9+1]=H_mu_s*mu;
3925 mom_vw_diff_ten[k*9+7]=H_mu_s*mu;
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;
3932 mom_wu_diff_ten[k*9+2]=H_mu_s*mu;
3934 mom_wv_diff_ten[k*9+5]=H_mu_s*mu;
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);
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;
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);
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;
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);
3968 mom_u_ham[k] = grad_p[k*3+0];
3969 dmom_u_ham_grad_p[k*3+0]=1.0;
3972 mom_v_ham[k] = grad_p[k*3+1];
3973 dmom_v_ham_grad_p[k*3+1]=1.0;
3976 mom_w_ham[k] = grad_p[k*3+2];
3977 dmom_w_ham_grad_p[k*3+2]=1.0;
3990 const double *grad_p,
3995 double *dmom_u_acc_u,
3997 double *dmom_v_acc_v,
3999 double *dmom_w_acc_w,
4001 double *dmass_adv_u,
4002 double *dmass_adv_v,
4003 double *dmass_adv_w,
4005 double *dmom_u_adv_u,
4006 double *dmom_u_adv_v,
4007 double *dmom_u_adv_w,
4009 double *dmom_v_adv_u,
4010 double *dmom_v_adv_v,
4011 double *dmom_v_adv_w,
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,
4023 double *dmom_u_ham_grad_p,
4025 double *dmom_v_ham_grad_p,
4027 double *dmom_w_ham_grad_p)
4031 for (k=0;k<nPoints;k++)
4039 dmom_u_acc_u[k]=1.0;
4042 dmom_v_acc_v[k]=1.0;
4045 dmom_w_acc_w[k]=1.0;
4048 mass_adv[k*3+0]=
u[k];
4049 mass_adv[k*3+1]=
v[k];
4050 mass_adv[k*3+2]=
w[k];
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;
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];
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];
4065 dmom_u_adv_v[k*3+1]=
u[k];
4067 dmom_u_adv_w[k*3+2]=
u[k];
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];
4074 dmom_v_adv_u[k*3+0]=
v[k];
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];
4080 dmom_v_adv_w[k*3+2]=
v[k];
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];
4087 dmom_w_adv_u[k*3+0]=
w[k];
4089 dmom_w_adv_v[k*3+0]=
w[k];
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];
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;
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;
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;
4111 mom_u_source[k] = -g[0];
4112 mom_v_source[k] = -g[1];
4113 mom_w_source[k] = -g[2];
4116 mom_u_ham[k] = grad_p[k*3+0]/rho;
4117 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4120 mom_v_ham[k] = grad_p[k*3+1]/rho;
4121 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4124 mom_w_ham[k] = grad_p[k*3+2]/rho;
4125 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4138 const double *grad_p,
4142 double *dmom_u_acc_u,
4144 double *dmom_v_acc_v,
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,
4153 double *dmom_u_ham_grad_p,
4155 double *dmom_v_ham_grad_p)
4159 for (k=0;k<nPoints;k++)
4168 dmom_u_acc_u[k]=1.0;
4172 dmom_v_acc_v[k]=1.0;
4175 mass_adv[k*2+0]=
u[k];
4176 mass_adv[k*2+1]=
v[k];
4178 dmass_adv_u[k*2+0]=1.0;
4179 dmass_adv_v[k*2+1]=1.0;
4182 mom_u_diff_ten[k*4+0] = nu;
4183 mom_u_diff_ten[k*4+3] = nu;
4186 mom_v_diff_ten[k*4+0] = nu;
4187 mom_v_diff_ten[k*4+3] = nu;
4190 mom_u_source[k] = -g[0];
4191 mom_v_source[k] = -g[1];
4194 mom_u_ham[k] = grad_p[k*2+0]/rho;
4195 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
4198 mom_v_ham[k] = grad_p[k*2+1]/rho;
4199 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
4212 const double *grad_p,
4217 double *dmom_u_acc_u,
4219 double *dmom_v_acc_v,
4221 double *dmom_w_acc_w,
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,
4233 double *dmom_u_ham_grad_p,
4235 double *dmom_v_ham_grad_p,
4237 double *dmom_w_ham_grad_p)
4241 for (k=0;k<nPoints;k++)
4249 dmom_u_acc_u[k]=1.0;
4252 dmom_v_acc_v[k]=1.0;
4255 dmom_w_acc_w[k]=1.0;
4258 mass_adv[k*3+0]=
u[k];
4259 mass_adv[k*3+1]=
v[k];
4260 mass_adv[k*3+2]=
w[k];
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;
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;
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;
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;
4282 mom_u_source[k] = -g[0];
4283 mom_v_source[k] = -g[1];
4284 mom_w_source[k] = -g[2];
4287 mom_u_ham[k] = grad_p[k*3+0]/rho;
4288 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4291 mom_v_ham[k] = grad_p[k*3+1]/rho;
4292 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4295 mom_w_ham[k] = grad_p[k*3+2]/rho;
4296 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4309 const double *grad_p,
4313 double *dmom_u_acc_u,
4315 double *dmom_v_acc_v,
4317 double *dmass_adv_u,
4318 double *dmass_adv_v,
4320 double *dmom_u_adv_u,
4321 double *dmom_u_adv_v,
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,
4330 double *dmom_u_ham_grad_p,
4332 double *dmom_v_ham_grad_p)
4336 for (k=0;k<nPoints;k++)
4339 H = fmax(0.0,fmin(1.0,vof[k]));
4345 dmom_u_acc_u[k]=1.0;
4349 dmom_v_acc_v[k]=1.0;
4352 mass_adv[k*2+0]=
u[k];
4353 mass_adv[k*2+1]=
v[k];
4355 dmass_adv_u[k*2+0]=1.0;
4356 dmass_adv_v[k*2+1]=1.0;
4359 mom_u_adv[k*2+0]=
u[k]*
u[k];
4360 mom_u_adv[k*2+1]=
u[k]*
v[k];
4362 dmom_u_adv_u[k*2+0]=2.0*
u[k];
4363 dmom_u_adv_u[k*2+1]=
v[k];
4365 dmom_u_adv_v[k*2+1]=
u[k];
4368 mom_v_adv[k*2+0]=
v[k]*
u[k];
4369 mom_v_adv[k*2+1]=
v[k]*
v[k];
4371 dmom_v_adv_u[k*2+0]=
v[k];
4373 dmom_v_adv_v[k*2+0]=
u[k];
4374 dmom_v_adv_v[k*2+1]=2.0*
v[k];
4377 mom_u_diff_ten[k*4+0] = nu;
4378 mom_u_diff_ten[k*4+3] = nu;
4381 mom_v_diff_ten[k*4+0] = nu;
4382 mom_v_diff_ten[k*4+3] = nu;
4385 mom_u_source[k] = -g[0];
4386 mom_v_source[k] = -g[1];
4389 mom_u_ham[k] = grad_p[k*2+0]/rho;
4390 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
4393 mom_v_ham[k] = grad_p[k*2+1]/rho;
4394 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
4407 const double *grad_p,
4412 double *dmom_u_acc_u,
4414 double *dmom_v_acc_v,
4416 double *dmom_w_acc_w,
4418 double *dmass_adv_u,
4419 double *dmass_adv_v,
4420 double *dmass_adv_w,
4422 double *dmom_u_adv_u,
4423 double *dmom_u_adv_v,
4424 double *dmom_u_adv_w,
4426 double *dmom_v_adv_u,
4427 double *dmom_v_adv_v,
4428 double *dmom_v_adv_w,
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,
4440 double *dmom_u_ham_grad_p,
4442 double *dmom_v_ham_grad_p,
4444 double *dmom_w_ham_grad_p)
4448 for (k=0;k<nPoints;k++)
4450 H = fmax(0.0,fmin(1.0,vof[k]));
4456 dmom_u_acc_u[k]=1.0;
4459 dmom_v_acc_v[k]=1.0;
4462 dmom_w_acc_w[k]=1.0;
4465 mass_adv[k*3+0]=
u[k];
4466 mass_adv[k*3+1]=
v[k];
4467 mass_adv[k*3+2]=
w[k];
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;
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];
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];
4482 dmom_u_adv_v[k*3+1]=
u[k];
4484 dmom_u_adv_w[k*3+2]=
u[k];
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];
4491 dmom_v_adv_u[k*3+0]=
v[k];
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];
4497 dmom_v_adv_w[k*3+2]=
v[k];
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];
4504 dmom_w_adv_u[k*3+0]=
w[k];
4506 dmom_w_adv_v[k*3+0]=
w[k];
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];
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;
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;
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;
4528 mom_u_source[k] = -g[0];
4529 mom_v_source[k] = -g[1];
4530 mom_w_source[k] = -g[2];
4533 mom_u_ham[k] = grad_p[k*3+0]/rho;
4534 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4537 mom_v_ham[k] = grad_p[k*3+1]/rho;
4538 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4541 mom_w_ham[k] = grad_p[k*3+2]/rho;
4542 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4555 const double *grad_p,
4559 double *dmom_u_acc_u,
4561 double *dmom_v_acc_v,
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,
4570 double *dmom_u_ham_grad_p,
4572 double *dmom_v_ham_grad_p)
4576 for (k=0;k<nPoints;k++)
4579 H = fmax(0.0,fmin(1.0,vof[k]));
4585 dmom_u_acc_u[k]=1.0;
4589 dmom_v_acc_v[k]=1.0;
4592 mass_adv[k*2+0]=
u[k];
4593 mass_adv[k*2+1]=
v[k];
4595 dmass_adv_u[k*2+0]=1.0;
4596 dmass_adv_v[k*2+1]=1.0;
4599 mom_u_diff_ten[k*4+0] = nu;
4600 mom_u_diff_ten[k*4+3] = nu;
4603 mom_v_diff_ten[k*4+0] = nu;
4604 mom_v_diff_ten[k*4+3] = nu;
4607 mom_u_source[k] = -g[0];
4608 mom_v_source[k] = -g[1];
4611 mom_u_ham[k] = grad_p[k*2+0]/rho;
4612 dmom_u_ham_grad_p[k*2+0]=1.0/rho;
4615 mom_v_ham[k] = grad_p[k*2+1]/rho;
4616 dmom_v_ham_grad_p[k*2+1]=1.0/rho;
4629 const double *grad_p,
4634 double *dmom_u_acc_u,
4636 double *dmom_v_acc_v,
4638 double *dmom_w_acc_w,
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,
4650 double *dmom_u_ham_grad_p,
4652 double *dmom_v_ham_grad_p,
4654 double *dmom_w_ham_grad_p)
4658 for (k=0;k<nPoints;k++)
4660 H = fmax(0.0,fmin(1.0,vof[k]));
4666 dmom_u_acc_u[k]=1.0;
4669 dmom_v_acc_v[k]=1.0;
4672 dmom_w_acc_w[k]=1.0;
4675 mass_adv[k*3+0]=
u[k];
4676 mass_adv[k*3+1]=
v[k];
4677 mass_adv[k*3+2]=
w[k];
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;
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;
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;
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;
4699 mom_u_source[k] = -g[0];
4700 mom_v_source[k] = -g[1];
4701 mom_w_source[k] = -g[2];
4704 mom_u_ham[k] = grad_p[k*3+0]/rho;
4705 dmom_u_ham_grad_p[k*3+0]=1.0/rho;
4708 mom_v_ham[k] = grad_p[k*3+1]/rho;
4709 dmom_v_ham_grad_p[k*3+1]=1.0/rho;
4712 mom_w_ham[k] = grad_p[k*3+2]/rho;
4713 dmom_w_ham_grad_p[k*3+2]=1.0/rho;
4727 double vx,
vy, xk, yk;
4729 double one8 = 1.0/8.0;
4730 for (k=0; k < nPoints; k++)
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];
4740 df[k*nSpace+1] =
vy;
4750 const double *gradu,
4760 for (k=0; k < nPoints; k++)
4765 for (
id=0;
id < nSpace;
id++)
4767 f[k*nSpace+id] = 0.0;
4768 df[k*nSpace+id]= 0.0;
4770 bdotgrad+= gradu[k*nSpace+id]*b[id];
4771 dH[k*nSpace+id]=b[id];
4781 const double *gradu,
4791 for (k=0; k < nPoints; k++)
4796 for (
id=0;
id < nSpace;
id++)
4798 f[k*nSpace+id] = 0.0;
4799 df[k*nSpace+id]= 0.0;
4801 normgradu+= gradu[k*nSpace+id]*gradu[k*nSpace+id];
4804 normgradu = sqrt(normgradu);
4806 for (
id=0;
id < nSpace;
id++)
4808 dH[k*nSpace+id]=gradu[k*nSpace+id]/(normgradu+1.0e-8);
4817 const double *gradu,
4825 double v[3], xk, yk, vdotgrad, one8;
4829 for (k=0; k < nPoints; k++)
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);
4837 for (
id=0;
id < nSpace;
id++)
4839 f[k*nSpace+id] = 0.0;
4840 df[k*nSpace+id]= 0.0;
4842 vdotgrad+= gradu[k*nSpace+id]*
v[id];
4843 dH[k*nSpace+id]=
v[id];
4853 const double *gradu,
4861 double v[3], xk, yk, vdotgrad, one8;
4865 for (k=0; k < nPoints; k++)
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]);
4873 for (
id=0;
id < nSpace;
id++)
4875 f[k*nSpace+id] = 0.0;
4876 df[k*nSpace+id]= 0.0;
4878 vdotgrad+= gradu[k*nSpace+id]*
v[id];
4879 dH[k*nSpace+id]=
v[id];
4887 const double offset,
4889 const double *gradu,
4897 for (k=0; k < nPoints; k++)
4902 for (
id=0;
id < nSpace;
id++)
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];
5011 const double* gravity,
5016 const double thetaR,
5017 const double thetaSR,
5030 const int nSpace2=nSpace*nSpace;
5031 register double psiC,
5032 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5034 sBar,sqrt_sBar,DsBar_DpsiC,
5035 thetaW,DthetaW_DpsiC,
5036 vBar,vBar2,DvBar_DpsiC,
5039 thetaS=thetaR+thetaSR,
5041 for (k=0;k<nPoints;k++)
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;
5052 sBar = pow(onePlus_pcBar_n,-m);
5054 DsBar_DpsiC = alpha*(1.0-
n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5056 vBar = 1.0-pcBar_nM1*sBar;
5058 DvBar_DpsiC = -alpha*(
n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5060 thetaW = thetaSR*sBar + thetaR;
5061 DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5063 sqrt_sBar = sqrt(sBar);
5064 KW= KWs*sqrt_sBar*vBar2;
5066 ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5068 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5073 DthetaW_DpsiC = 0.0;
5078 rhom = rho*exp(beta*
u[k]);
5081 mass[k] = rhom*thetaW;
5082 dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
5083 for (I=0;I<nSpace;I++)
5085 f[k*nSpace+I] = rho2*KW*gravity[I];
5086 df[k*nSpace+I] = -rho2*DKW_DpsiC*gravity[I];
5088 a[k*nSpace2+I*nSpace+I] = rho*KW;
5089 da[k*nSpace2+I*nSpace+I] = -rho*DKW_DpsiC;
5096 const int nPointsPerSimplex,
5099 const double* gravity,
5103 const double thetaR,
5104 const double thetaSR,
5116 const int nSpace2=nSpace*nSpace;
5117 register double psiC,
5118 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5120 sBar,sqrt_sBar,DsBar_DpsiC,
5121 thetaW,DthetaW_DpsiC,
5122 vBar,vBar2,DvBar_DpsiC,
5125 thetaS=thetaR+thetaSR;
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},
5133 double daavg[3][3]= {{0.0,0.0,0.0},
5140 for (eN = 0; eN < nSimplices; eN++)
5142 mavg = 0.0; dmavg = 0.0; vol = 0.0;
5143 for (I=0; I < nSpace; I++)
5145 favg[I] =0.0; dfavg[I] = 0.0;
5146 for (J=0; J < nSpace;J++)
5148 aavg[I][J] = 0.0; daavg[I][J] = 0.0;
5151 for (k=0;k<nPointsPerSimplex;k++)
5153 psiC = -
u[eN*nPointsPerSimplex + k];
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;
5162 sBar = pow(onePlus_pcBar_n,-m);
5164 DsBar_DpsiC = alpha*(1.0-
n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5166 vBar = 1.0-pcBar_nM1*sBar;
5168 DvBar_DpsiC = -alpha*(
n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5170 thetaW = thetaSR*sBar + thetaR;
5171 DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5173 sqrt_sBar = sqrt(sBar);
5174 KW= KWs*sqrt_sBar*vBar2;
5176 ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5178 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5183 DthetaW_DpsiC = 0.0;
5191 vol += dV[eN*nPointsPerSimplex + k];
5192 mavg += rho*thetaW*dV[eN*nPointsPerSimplex + k];
5193 dmavg +=-rho*DthetaW_DpsiC*dV[eN*nPointsPerSimplex + k];
5197 dmass[eN*nPointsPerSimplex+ k] =-rho*DthetaW_DpsiC;
5199 for (I=0; I < nSpace; I++)
5201 favg[I] += rho2*KW*gravity[I]*dV[eN*nPointsPerSimplex + k];
5202 dfavg[I]+=-rho2*DKW_DpsiC*gravity[I]*dV[eN*nPointsPerSimplex + k];
5204 aavg[I][J] += rho*KW*dV[eN*nPointsPerSimplex + k];
5205 daavg[I][J] +=-rho*DKW_DpsiC*dV[eN*nPointsPerSimplex + k];
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;
5216 for (k=0; k < nPointsPerSimplex; k++)
5218 mass[eN*nPointsPerSimplex + k] = mavg/vol;
5220 for (I=0; I < nSpace; I++)
5222 f[eN*nPointsPerSimplex*nSpace + k*nSpace + I] = favg[I]/vol;
5224 for (J=0; J < nSpace; J++)
5228 a[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J] = aavg[I][J]/vol;
5233 a[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + J] = 0.0;
5244 const int nElementBoundaries_element,
5245 const int nPointsPerElementBoundary,
5248 const double* gravity,
5252 const double thetaR,
5253 const double thetaSR,
5265 const int nSpace2=nSpace*nSpace;
5266 register double psiC,
5267 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5269 sBar,sqrt_sBar,DsBar_DpsiC,
5270 thetaW,DthetaW_DpsiC,
5271 vBar,vBar2,DvBar_DpsiC,
5274 thetaS=thetaR+thetaSR;
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},
5282 double daavg[3][3]= {{0.0,0.0,0.0},
5290 for (eN = 0; eN < nElements; eN++)
5292 for (ebN = 0; ebN < nElementBoundaries_element; ebN++)
5294 mavg = 0.0; dmavg = 0.0; vol = 0.0;
5295 for (I=0; I < nSpace; I++)
5297 favg[I] =0.0; dfavg[I] = 0.0;
5298 for (J=0; J < nSpace;J++)
5300 aavg[I][J] = 0.0; daavg[I][J] = 0.0;
5303 for (k=0;k<nPointsPerElementBoundary;k++)
5305 psiC = -
u[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k];
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;
5314 sBar = pow(onePlus_pcBar_n,-m);
5316 DsBar_DpsiC = alpha*(1.0-
n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5318 vBar = 1.0-pcBar_nM1*sBar;
5320 DvBar_DpsiC = -alpha*(
n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5322 thetaW = thetaSR*sBar + thetaR;
5323 DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5325 sqrt_sBar = sqrt(sBar);
5326 KW= KWs*sqrt_sBar*vBar2;
5328 ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5330 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5335 DthetaW_DpsiC = 0.0;
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];
5348 dmass[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k] =
5350 for (I=0; I < nSpace; I++)
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];
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];
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;
5374 for (k=0; k < nPointsPerElementBoundary; k++)
5376 mass[eN*nElementBoundaries_element*nPointsPerElementBoundary + ebN*nPointsPerElementBoundary + k] = mavg/vol;
5379 for (I=0; I < nSpace; I++)
5381 f[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace + ebN*nPointsPerElementBoundary*nSpace +
5382 k*nSpace + I] = favg[I]/vol;
5385 for (J=0; J < nSpace; J++)
5390 a[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace2 + ebN*nPointsPerElementBoundary*nSpace2 +
5391 k*nSpace2 + I*nSpace + J] = aavg[I][J]/vol;
5399 a[eN*nElementBoundaries_element*nPointsPerElementBoundary*nSpace2 + ebN*nPointsPerElementBoundary*nSpace2 +
5400 k*nSpace2 + I*nSpace + J] = 0.0;
5413 const int nPointsPerSimplex,
5416 const double *gravity,
5417 const double *alpha,
5419 const double *thetaR,
5420 const double *thetaSR,
5432 const int nSpace2=nSpace*nSpace;
5433 register double psiC,
5434 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5436 sBar,sqrt_sBar,DsBar_DpsiC,
5437 thetaW,DthetaW_DpsiC,
5438 vBar,vBar2,DvBar_DpsiC,
5441 register double thetaS,m;
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},
5449 double daavg[3][3]= {{0.0,0.0,0.0},
5456 for (eN = 0; eN < nSimplices; eN++)
5458 mavg = 0.0; dmavg = 0.0; vol = 0.0;
5459 for (I=0; I < nSpace; I++)
5461 favg[I] =0.0; dfavg[I] = 0.0;
5462 for (J=0; J < nSpace;J++)
5464 aavg[I][J] = 0.0; daavg[I][J] = 0.0;
5467 for (k=0;k<nPointsPerSimplex;k++)
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];
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;
5480 sBar = pow(onePlus_pcBar_n,-m);
5482 DsBar_DpsiC = alpha[eN*nPointsPerSimplex + k]*(1.0-
n[eN*nPointsPerSimplex + k])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5484 vBar = 1.0-pcBar_nM1*sBar;
5486 DvBar_DpsiC = -alpha[eN*nPointsPerSimplex + k]*(
n[eN*nPointsPerSimplex + k]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5488 thetaW = thetaSR[eN*nPointsPerSimplex + k]*sBar + thetaR[eN*nPointsPerSimplex + k];
5489 DthetaW_DpsiC = thetaSR[eN*nPointsPerSimplex + k] * DsBar_DpsiC;
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
5496 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5501 DthetaW_DpsiC = 0.0;
5502 KW = KWs[eN*nPointsPerSimplex + k];
5509 vol += dV[eN*nPointsPerSimplex + k];
5510 mavg += rho*thetaW*dV[eN*nPointsPerSimplex + k];
5511 dmavg +=-rho*DthetaW_DpsiC*dV[eN*nPointsPerSimplex + k];
5514 dmass[eN*nPointsPerSimplex+ k] = -rho*DthetaW_DpsiC;
5515 for (I=0; I < nSpace; I++)
5517 favg[I] += rho2*KW*gravity[I]*dV[eN*nPointsPerSimplex + k];
5518 dfavg[I]+=-rho2*DKW_DpsiC*gravity[I]*dV[eN*nPointsPerSimplex + k];
5520 df[eN*nPointsPerSimplex*nSpace+ k*nSpace + I] = -rho2*DKW_DpsiC*gravity[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;
5532 for (k=0; k < nPointsPerSimplex; k++)
5534 mass[eN*nPointsPerSimplex + k] = mavg/vol;
5536 for (I=0; I < nSpace; I++)
5539 f[eN*nPointsPerSimplex*nSpace + k*nSpace + I] = favg[I]/vol;
5541 a[eN*nPointsPerSimplex*nSpace2 + k*nSpace2 + I*nSpace + I] = aavg[I][I]/vol;
5553 const double* gravity,
5558 const double thetaR,
5559 const double thetaSR,
5572 const int nSpace2=nSpace*nSpace;
5573 register double psiC,
5574 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
5576 sBar,sqrt_sBar,DsBar_DpsiC,
5577 thetaW,DthetaW_DpsiC,
5578 vBar,vBar2,DvBar_DpsiC,
5580 thetaS=thetaR+thetaSR;
5582 register double elev;
5583 for (k=0;k<nPoints;k++)
5586 for (I=0; I < nSpace; I++)
5587 elev += gravity[I]*x[k*3+I];
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;
5602 sBar = pow(onePlus_pcBar_n,-m);
5604 DsBar_DpsiC = alpha*(1.0-
n)*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5606 vBar = 1.0-pcBar_nM1*sBar;
5608 DvBar_DpsiC = -alpha*(
n-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5610 thetaW = thetaSR*sBar + thetaR;
5611 DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
5613 sqrt_sBar = sqrt(sBar);
5614 KW= KWs*sqrt_sBar*vBar2;
5616 ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5618 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5623 DthetaW_DpsiC = 0.0;
5627 mass[k] = rho*thetaW;
5628 dmass[k] = -rho*DthetaW_DpsiC;
5638 for (I=0;I<nSpace;I++)
5640 f[k*nSpace+I] = 0.0;
5641 df[k*nSpace+I] = 0.0;
5643 phi[k] -= rho*gravity[I]*x[k*3+I];
5645 a[k*nSpace2+I*nSpace+I] = rho*KW;
5646 da[k*nSpace2+I*nSpace+I] = -rho*DKW_DpsiC;
5652 const int nPointsPerSimplex,
5659 for (eN = 0; eN < nSimplices; eN++)
5661 ravg = 0.0; vol = 0.0;
5662 for (k=0; k < nPointsPerSimplex; k++)
5664 vol += dV[eN*nPointsPerSimplex+k];
5665 ravg +=
r[eN*nPointsPerSimplex + k]*dV[eN*nPointsPerSimplex+k];
5670 for (k=0; k < nPointsPerSimplex; k++)
5672 r[eN*nPointsPerSimplex + k] = ravg;
5677 const int nPointsPerSimplex,
5685 double ravg[nSpace];
5686 for (eN = 0; eN < nSimplices; eN++)
5689 for (I=0; I < nSpace; I++)
5691 for (k=0; k < nPointsPerSimplex; k++)
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];
5698 for (k=0; k < nPointsPerSimplex; k++)
5700 for (I=0; I < nSpace; I++)
5701 r[eN*nPointsPerSimplex*nSpace + k*nSpace + I] = ravg[I]/vol;
5706 const int nPointsPerSimplex,
5711 int eN,k,I,J,nSpace2;
5714 double ravg[nSpace*nSpace];
5715 nSpace2 = nSpace*nSpace;
5716 for (eN = 0; eN < nSimplices; eN++)
5719 for (I=0; I < nSpace2; I++)
5721 for (k=0; k < nPointsPerSimplex; k++)
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]
5728 dV[eN*nPointsPerSimplex+k];
5731 for (k=0; k < nPointsPerSimplex; k++)
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;
5741 const int nPointsPerSimplex,
5746 const int* materialTypes,
5749 const double* gravity,
5750 const double* alpha,
5752 const double* thetaR,
5753 const double* thetaSR,
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,
5769 sBar,sqrt_sBar,DsBar_DpsiC,
5770 thetaW,DthetaW_DpsiC,
5771 vBar,vBar2,DvBar_DpsiC,
5777 const int nnz = rowptr[nSpace];
5778 for (i=0; i < nSimplex; i++)
5780 matID= materialTypes[i];
5781 for (j=0;j<nPointsPerSimplex;j++)
5783 k = i*nPointsPerSimplex + j;
5785 m = 1.0 - 1.0/
n[matID];
5786 thetaS = thetaR[matID] + thetaSR[matID];
5789 pcBar = alpha[matID]*psiC;
5790 pcBarStar = fmax(pcBar,pc_eps);
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;
5797 sBar = pow(onePlus_pcBar_n,-m);
5799 DsBar_DpsiC = alpha[matID]*(1.0-
n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5801 vBar = 1.0-pcBar_nM1*sBar;
5803 DvBar_DpsiC = -alpha[matID]*(
n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5805 thetaW = thetaSR[matID]*sBar + thetaR[matID];
5806 DthetaW_DpsiC = thetaSR[matID] * DsBar_DpsiC;
5808 sqrt_sBar = sqrt(sBar);
5809 KWr= sqrt_sBar*vBar2;
5810 DKWr_DpsiC= ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5812 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5817 DthetaW_DpsiC = 0.0;
5822 betauStar = fmin(beta*
u[k],1000.0);
5823 rhom = rho*exp(betauStar);
5825 mass[k] = rhom*thetaW;
5826 dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
5827 vol_frac[k] = thetaW;
5830 for (I=0;I<nSpace;I++)
5832 f[k*nSpace+I] = 0.0;
5833 df[k*nSpace+I] = 0.0;
5834 for (ii=rowptr[I]; ii < rowptr[I+1]; ii++)
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];
5848 const int nPointsPerSimplex,
5850 double linear_break,
5853 const int* materialTypes,
5856 const double* gravity,
5857 const double* alpha,
5859 const double* thetaR,
5860 const double* thetaSR,
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,
5876 sBar,sqrt_sBar,DsBar_DpsiC,
5877 thetaW,DthetaW_DpsiC,
5878 vBar,vBar2,DvBar_DpsiC,
5884 const int nnz = rowptr[nSpace];
5885 pcBarStar = linear_break;
5887 for (i=0; i < nSimplex; i++)
5889 matID= materialTypes[i];
5890 for (j=0;j<nPointsPerSimplex;j++)
5892 k = i*nPointsPerSimplex + j;
5894 m = 1.0 - 1.0/
n[matID];
5895 thetaS = thetaR[matID] + thetaSR[matID];
5898 pcBar = alpha[matID]*psiC;
5899 if (psiC <= pcBarStar)
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;
5907 sBar = pow(onePlus_pcBar_n,-m);
5909 DsBar_DpsiC = alpha[matID]*(1.0-
n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5911 vBar = 1.0-pcBar_nM1*sBar;
5913 DvBar_DpsiC = -alpha[matID]*(
n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5916 sqrt_sBar = sqrt(sBar);
5917 KWrStar= sqrt_sBar*vBar2;
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;
5925 DKWr_DpsiC= (KWrStar - 1.0)/(pcBarStar-0.0);
5926 KWr = DKWr_DpsiC*(psiC-0.0) + 1.0;
5930 pcBar = alpha[matID]*psiC;
5931 pcBarStar = fmax(pcBar,1.0e-8);
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;
5938 sBar = pow(onePlus_pcBar_n,-m);
5940 DsBar_DpsiC = alpha[matID]*(1.0-
n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
5942 vBar = 1.0-pcBar_nM1*sBar;
5944 DvBar_DpsiC = -alpha[matID]*(
n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
5946 thetaW = thetaSR[matID]*sBar + thetaR[matID];
5947 DthetaW_DpsiC = thetaSR[matID] * DsBar_DpsiC;
5949 sqrt_sBar = sqrt(sBar);
5950 KWr= sqrt_sBar*vBar2;
5951 DKWr_DpsiC= ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
5953 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
5960 DthetaW_DpsiC = 0.0;
5965 betauStar = fmin(beta*
u[k],1000.0);
5966 rhom = rho*exp(betauStar);
5968 mass[k] = rhom*thetaW;
5969 dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
5970 vol_frac[k] = thetaW;
5973 for (I=0;I<nSpace;I++)
5975 f[k*nSpace+I] = 0.0;
5976 df[k*nSpace+I] = 0.0;
5977 for (ii=rowptr[I]; ii < rowptr[I+1]; ii++)
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];
5991 const int nPointsPerSimplex,
5993 const int* materialTypes,
5996 const double* gravity,
5997 const double* alpha,
5999 const double* thetaR,
6000 const double* thetaSR,
6011 const int nSpace2=nSpace*nSpace;
6012 register double psiC,
6013 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
6015 sBar,sqrt_sBar,DsBar_DpsiC,
6016 thetaW,DthetaW_DpsiC,
6017 vBar,vBar2,DvBar_DpsiC,
6022 for (i=0; i < nSimplex; i++)
6024 matID= materialTypes[i];
6025 for (j=0;j<nPointsPerSimplex;j++)
6027 k = i*nPointsPerSimplex + j;
6029 m = 1.0 - 1.0/
n[matID];
6030 thetaS = thetaR[matID] + thetaSR[matID];
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;
6039 sBar = pow(onePlus_pcBar_n,-m);
6041 DsBar_DpsiC = alpha[matID]*(1.0-
n[matID])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
6043 vBar = 1.0-pcBar_nM1*sBar;
6045 DvBar_DpsiC = -alpha[matID]*(
n[matID]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
6047 thetaW = thetaSR[matID]*sBar + thetaR[matID];
6048 DthetaW_DpsiC = thetaSR[matID] * DsBar_DpsiC;
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
6055 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
6060 DthetaW_DpsiC = 0.0;
6065 rhom = rho*exp(beta*
u[k]);
6067 mass[k] = rhom*thetaW;
6068 dmass[k] = -rhom*DthetaW_DpsiC+drhom*thetaW;
6071 for (I=0;I<nSpace;I++)
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;
6084 const int nPointsPerSimplex,
6086 const int* materialTypes,
6087 const double epsFact,
6090 const double* elementDiameter,
6091 const double* gravity,
6092 const double* alpha,
6094 const double* thetaR,
6095 const double* thetaSR,
6106 const int nSpace2=nSpace*nSpace;
6107 register double psiC,
6108 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
6110 sBar,sqrt_sBar,DsBar_DpsiC,
6111 thetaW,DthetaW_DpsiC,
6112 vBar,vBar2,DvBar_DpsiC,
6116 rhom,drhom,m,eps,hStar,dhStar;
6117 for (i=0; i < nSimplex; i++)
6119 matID= materialTypes[i];
6120 eps = epsFact*elementDiameter[i];
6121 for (j=0;j<nPointsPerSimplex;j++)
6123 k = i*nPointsPerSimplex + j;
6124 thetaS = thetaR[matID] + thetaSR[matID];
6129 mass[k] = thetaS*hStar;
6130 dmass[k] = thetaS*dhStar;
6131 for (I=0;I<nSpace;I++)
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;
6144 const int nPointsPerSimplex,
6146 const int* materialTypes,
6149 const double* gravity,
6154 const double* S_wirr,
6155 const double* S_nwr,
6166 const int nSpace2=nSpace*nSpace;
6167 register double psiC,
6172 const double reg_diff=1.0e-2;
6173 for (i=0; i < nSimplex; i++)
6175 matID= materialTypes[i];
6176 for (j=0;j<nPointsPerSimplex;j++)
6178 k = i*nPointsPerSimplex + j;
6180 if (psiC > 0.0 && psiC < psiD[matID])
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];
6188 else if (psiC <= 0.0)
6191 Sw = 1.0 - S_nwr[matID];
6193 kr = kr0[matID]+reg_diff;
6205 rhom = rho*exp(beta*
u[k]);
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++)
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;
6221 const int nPointsPerSimplex,
6223 const int* materialTypes,
6226 const double* gravity,
6231 const double* S_wirr,
6232 const double* S_nwr,
6245 const int nSpace2=nSpace*nSpace;
6246 register double psiC,pcBar,
6251 const double reg_diff=1.0e-2;
6252 for (i=0; i < nSimplex; i++)
6254 matID= materialTypes[i];
6255 for (j=0;j<nPointsPerSimplex;j++)
6257 k = i*nPointsPerSimplex + j;
6259 if (psiC >= psiD[matID])
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;
6274 Sw = 1.0 - S_nwr[matID];
6304 rhom = rho*exp(beta*
u[k]);
6307 mass[k] = rhom*Sw*
phi[matID];
6308 dmass[k] = rhom*dSw*
phi[matID]+drhom*Sw*
phi[matID];
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;
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;
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;
6340 const double* gravity,
6341 const double* alpha,
6343 const double* thetaR,
6344 const double* thetaSR,
6355 const int nSpace2=nSpace*nSpace;
6356 register double psiC,
6357 pcBar,pcBar_n,pcBar_nM1,pcBar_nM2,
6359 sBar,sqrt_sBar,DsBar_DpsiC,
6360 thetaW,DthetaW_DpsiC,
6361 vBar,vBar2,DvBar_DpsiC,
6366 for (k=0;k<nPoints;k++)
6370 thetaS = thetaR[k] + thetaSR[k];
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;
6379 sBar = pow(onePlus_pcBar_n,-m);
6381 DsBar_DpsiC = alpha[k]*(1.0-
n[k])*(sBar/onePlus_pcBar_n)*pcBar_nM1;
6383 vBar = 1.0-pcBar_nM1*sBar;
6385 DvBar_DpsiC = -alpha[k]*(
n[k]-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
6387 thetaW = thetaSR[k]*sBar + thetaR[k];
6388 DthetaW_DpsiC = thetaSR[k] * DsBar_DpsiC;
6390 sqrt_sBar = sqrt(sBar);
6391 KW= KWs[k]*sqrt_sBar*vBar2;
6393 ((0.5/sqrt_sBar)*DsBar_DpsiC*vBar2
6395 2.0*sqrt_sBar*vBar*DvBar_DpsiC);
6400 DthetaW_DpsiC = 0.0;
6404 mass[k] = rho*thetaW;
6405 dmass[k] = -rho*DthetaW_DpsiC;
6406 for (I=0;I<nSpace;I++)
6408 f[k*nSpace+I] = rho2*KW*gravity[I];
6409 df[k*nSpace+I] = -rho2*DKW_DpsiC*gravity[I];