Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 20944) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 20945) @@ -3,6 +3,6 @@ */ -#ifndef MATEARL_H_ -#define MATEARL_H_ +#ifndef MATESTAR_H_ +#define MATESTAR_H_ /*Headers:*/ @@ -85,3 +85,3 @@ }; -#endif /* _MATEARL_H_ */ +#endif /* _MATESTAR_H_ */ Index: /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 20944) +++ /issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 20945) @@ -6,5 +6,5 @@ /*Intermediaries*/ - IssmDouble omega[3],omega_norm; + IssmDouble omega[3]; IssmDouble nrsp[3],nrsp_norm; IssmDouble eps[3][3],epso; @@ -81,4 +81,5 @@ /*Intermediaries*/ IssmDouble omega_norm; + IssmDouble omega_rigid[3]; /*Create vorticity vector*/ @@ -87,13 +88,18 @@ /*Create vorticity vector, corrected for rigid body rotation - * \overline{\omega} =\omega - \Omega - * =-\nabla\times{\bf v} - 2*\nabla U\times{\bf v}/U; + * \overline{\omega} =\omega - \omega_rigid + * =\nabla\times{\bf v} - 2*U*\kappa*\hat{b}; + * =\nabla\times{\bf v} - (2*{\bf v}\times(({\bf v}\cdot\nabla)*{\bf v}))/U^2 * check the magnitude of the second term -- if it is small, then the two * vorticities (omega and first term in omega) are approx. equal * * */ - omega[0] = -(dvz[1] - dvy[2]) + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag; - omega[1] = -(dvx[2] - dvz[0]) + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag; - omega[2] = -(dvy[0] - dvx[1]) + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag; + omega_rigid[0] = 2*(vy*(vx*dvz[0]+vy*dvz[1]+vz*dvz[2]) - vz*(vx*dvy[0]+vy*dvy[1]+vz*dvy[2]))/(vmag*vmag); + omega_rigid[1] = 2*(vz*(vx*dvx[0]+vy*dvx[1]+vz*dvx[2]) - vx*(vx*dvz[0]+vy*dvz[1]+vz*dvz[2]))/(vmag*vmag); + omega_rigid[2] = 2*(vx*(vx*dvy[0]+vy*dvy[1]+vz*dvy[2]) - vy*(vx*dvx[0]+vy*dvx[1]+vz*dvx[2]))/(vmag*vmag); + + omega[0] = (dvz[1] - dvy[2]) - omega_rigid[0]; + omega[1] = (dvx[2] - dvz[0]) - omega_rigid[1]; + omega[2] = (dvy[0] - dvx[1]) - omega_rigid[2]; /*Take out vorticity component aligned with the velocity*/ Index: /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 20944) +++ /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 20945) @@ -9,6 +9,6 @@ * case 6: 3d test with sinusoidal functions, non homogeneous Dirichlet conditions * - * case 201: 3d test quadratic functions FS, EARL flow law - * case 202: 3d test quadratic functions HO, EARL flow law + * case 201: 3d test quadratic functions FS, ESTAR flow law + * case 202: 3d test quadratic functions HO, ESTAR flow law */ @@ -83,4 +83,5 @@ case 212: return (16.0L/9.0L)*x*(4*x + 1)/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 1.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 4.0L/3.0L)) + (64.0L/3.0L)*x*(-(4*x + 1)*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/pow(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11, 2) + 2*(pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (y + 3*z)/fabs(y + 3*z))*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)) + (4.0L/9.0L)*(4*x + 1)/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 1.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 4.0L/3.0L)) - 4/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 1.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)) + (16.0L/3.0L)*(-(4*x + 1)*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/pow(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11, 2) + 2*(pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (y + 3*z)/fabs(y + 3*z))*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)); + default: _error_("FS analytical solution"<