Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17302) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17303) @@ -178,7 +178,5 @@ ./shared/LatLong/Xy2llx.cpp\ ./shared/FSanalyticals/fsanalyticals.h\ - ./shared/FSanalyticals/FSanalytical1.cpp\ - ./shared/FSanalyticals/FSanalytical2.cpp\ - ./shared/FSanalyticals/FSanalytical3.cpp\ + ./shared/FSanalyticals/fsanalyticals.cpp\ ./shared/Enum/Enum.h\ ./shared/Enum/EnumDefinitions.h\ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17302) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17303) @@ -6,5 +6,5 @@ #include "../solutionsequences/solutionsequences.h" -//#define FSANALYTICAL 5 +//#define FSANALYTICAL 1 /*Model processing*/ @@ -2978,27 +2978,9 @@ y_coord=element->GetYcoord(gauss); if(dim==3) z_coord=element->GetZcoord(gauss); - - #if FSANALYTICAL == 1 - forcex=fx1(x_coord,y_coord,z_coord); - forcey=fy1(x_coord,y_coord,z_coord); - forcez=fz1(x_coord,y_coord,z_coord); - #elif FSANALYTICAL == 2 - forcex=fx2(x_coord,y_coord,z_coord); - forcey=fy2(x_coord,y_coord,z_coord); - forcez=fz2(x_coord,y_coord,z_coord); - #elif FSANALYTICAL == 3 - forcex=fx3(x_coord,y_coord); - forcey=fy3(x_coord,y_coord); - #elif FSANALYTICAL == 4 - forcex=fx4(x_coord,y_coord,z_coord); - forcey=fy4(x_coord,y_coord,z_coord); - forcez=fz4(x_coord,y_coord,z_coord); - #elif FSANALYTICAL == 5 - forcex=fx5(x_coord,y_coord,z_coord); - forcey=fy5(x_coord,y_coord,z_coord); - forcez=fz5(x_coord,y_coord,z_coord); - #else - _error_("FS analytical not implemented yet"); - #endif + else z_coord=0.; + + forcex=fx(x_coord,y_coord,z_coord,FSANALYTICAL); + forcey=fy(x_coord,y_coord,z_coord,FSANALYTICAL); + forcez=fz(x_coord,y_coord,z_coord,FSANALYTICAL); for(i=0;i -#include "../Numerics/types.h" -#include "../Numerics/constants.h" - -IssmDouble fx1(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble p = 2.0; - IssmDouble q = 2.0; - IssmDouble functionx; - - functionx = 4. * PI*PI * p*p * (cos(p*PI*x)-1) * sin(p*PI*y) * sin(p*PI*z)+ 2. * PI*PI*p*p* sin(p*PI*y) * sin(p*PI*z) * cos(p*PI*x) + q * PI * cos(q*PI*x) * sin(q*PI*y) * sin(q*PI*z); - - return functionx; -} -IssmDouble fy1(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble p = 2.0; - IssmDouble q = 2.0; - IssmDouble functiony; - - functiony = - 4. * PI*PI * p*p * (cos(p*PI*y)-1) * sin(p*PI*x) * sin(p*PI*z)- 2. * PI*PI * p*p * sin(p*PI*x) * cos(p*PI*y) * sin(p*PI*z)+ q * PI * sin(q*PI*x) * cos(q*PI*y) * sin(q*PI*z); - - return functiony; -} -IssmDouble fz1(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble p = 2.0; - IssmDouble q = 2.0; - IssmDouble functionz; - - functionz = - 2. * PI*PI * p*p * (cos(p*PI*z)-1) * sin(p*PI*x) * sin(p*PI*y)-PI*PI * p*p * sin(p*PI*x) * sin(p*PI*y) * cos(p*PI*z)+ q * PI * sin(q*PI*x) * sin(q*PI*y) * cos(q*PI*z); - - return functionz; -} Index: sm/trunk-jpl/src/c/shared/FSanalyticals/FSanalytical2.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/FSanalyticals/FSanalytical2.cpp (revision 17302) +++ (revision ) @@ -1,36 +1,0 @@ -/*!\file fx.cpp - * \brief: analytical fonction for FS - * test case for a cube with Dirichlet conditions on all faces and no variation in z - */ - -#include -#include "../Numerics/types.h" -#include "../Numerics/constants.h" - -IssmDouble fx2(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble p = 2.0; - IssmDouble q = 2.0; - IssmDouble functionx; - - functionx = 2.*(cos(PI*p*x) - 1)*PI*PI*p*p*sin(PI*p*y) + 3*PI*PI*p*p*sin(PI*p*y)*cos(PI*p*x) + PI*q*sin(PI*q*y)*cos(PI*q*x); - - return functionx; -} -IssmDouble fy2(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble p = 2.0; - IssmDouble q = 2.0; - IssmDouble functiony; - - functiony = -(cos(PI*p*y) - 1)*PI*PI*p*p*sin(PI*p*x) + PI*q*sin(PI*q*x)*cos(PI*q*y); - - return functiony; -} -IssmDouble fz2(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble p = 2.0; - IssmDouble q = 2.0; - IssmDouble functionz; - - functionz = 2*PI*PI*p*p*sin(PI*p*x)*sin(PI*p*y); - - return functionz; -} Index: sm/trunk-jpl/src/c/shared/FSanalyticals/FSanalytical3.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/FSanalyticals/FSanalytical3.cpp (revision 17302) +++ (revision ) @@ -1,72 +1,0 @@ -/*!\file fx.cpp - * \brief: analytical fonction for FS - * test case for a cube with Dirichlet conditions on all faces and no variation in z - */ - -#include -#include "../Numerics/types.h" -#include "../Numerics/constants.h" - -IssmDouble fx3(IssmDouble x,IssmDouble y){ - IssmDouble functionx; - IssmDouble z=y; - - functionx = 4*pow(x, 2)*z*pow(x - 1, 2) + 4*pow(x, 2)*z*(z - 1)*(2*z - 1) + 4*pow(x, 2)*pow(x - 1, 2)*(z - 1) + 2*pow(x, 2)*pow(x - 1, 2)*(2*z - 1) + 16*x*z*(x - 1)*(z - 1)*(2*z - 1) - 4*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) + 6*pow(z, 2)*(6*pow(x, 2) - 6*x + 1) + 4*z*pow(x - 1, 2)*(z - 1)*(2*z - 1) - 2*z*(6*pow(x, 2) - 6*x + 1) + z - 1.0L/2.0L; - return functionx; -} -IssmDouble fy3(IssmDouble x,IssmDouble y){ - IssmDouble functiony; - IssmDouble z=y; - - functiony = -8*pow(x, 3) + 4*pow(x, 2)*z*(x - 1)*(z - 1) + 2*pow(x, 2)*z*(x - 1)*(2*z - 1) + 2*pow(x, 2)*(x - 1)*(z - 1)*(2*z - 1) + 12*pow(x, 2) + 4*x*z*pow(x - 1, 2)*(z - 1) + 2*x*z*pow(x - 1, 2)*(2*z - 1) + 2*x*pow(x - 1, 2)*(z - 1)*(2*z - 1) - 3*x - 6*pow(z, 4)*(2*x - 1) + 12*pow(z, 3)*(2*x - 1) - 6*pow(z, 2)*(2*x - 1) - 24*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + 24*z*(2*pow(x, 3) - 3*pow(x, 2) + x) - 1.0L/2.0L; - - return functiony; -} -IssmDouble fx4(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble functionx; - - functionx = 4*pow(x, 2)*y*z*pow(x - 1, 2)*(z - 1) + 4*pow(x, 2)*y*z*(y - 1)*(2*y - 1)*(z - 1) + 2*pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1) + 4*pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 2*pow(x, 2)*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 4*x*pow(y, 2)*z*(x - 1)*(y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(2*x - 1)*(y - 1)*(z - 1) - 4*x*y*z*(x - 1)*pow(y - 1, 2)*(z - 1) + 16*x*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*y*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 2*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 4*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1) + y - 1.0L/2.0L; - - return functionx; -} -IssmDouble fy4(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble functiony; - - functiony = 4*pow(x, 2)*y*z*(x - 1)*(y - 1)*(z - 1) + 2*pow(x, 2)*y*z*(x - 1)*(2*y - 1)*(z - 1) + 2*pow(x, 2)*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 4*x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(z - 1) - 4*x*pow(y, 2)*z*pow(y - 1, 2)*(z - 1) - 2*x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2) + 4*x*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 2*x*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 16*x*y*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 2*x*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 4*x*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1) + x - 4*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - 2*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 1.0L/2.0L; - - return functiony; -} -IssmDouble fz4(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble functionz; - - functionz = 2*pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1) + 2*pow(x, 2)*y*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1) - 2*x*pow(y, 2)*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 2*x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) - 2*x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + 2*x*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*y*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1); - - return functionz; -} -IssmDouble fx5(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble functionx; - IssmDouble p=2; - IssmDouble q=2; - - functionx = 4*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) + 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x) + PI*q*sin(PI*q*y)*sin(PI*q*z)*cos(PI*q*x); - - return functionx; -} -IssmDouble fy5(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble functiony; - IssmDouble p=2; - IssmDouble q=2; - - functiony = -2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*sin(PI*p*x)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*z)*cos(PI*p*y) + PI*q*sin(PI*q*x)*sin(PI*q*z)*cos(PI*q*y); - - return functiony; -} -IssmDouble fz5(IssmDouble x,IssmDouble y,IssmDouble z){ - IssmDouble functionz; - IssmDouble p=2; - IssmDouble q=2; - - functionz = -2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*z) - 1)*sin(PI*p*x)*sin(PI*p*y) - pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*y)*cos(PI*p*z) + PI*q*sin(PI*q*x)*sin(PI*q*y)*cos(PI*q*z); - - return functionz; -} Index: /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 17303) +++ /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 17303) @@ -0,0 +1,229 @@ +/*!\file fx.cpp + * \brief: analytical fonction for FS + * test case for a cube with Dirichlet conditions on all faces and no variation in z + * case 1: 2d test with quadratic functions + * case 2: 2d test with sinusoidal functions + * case 3: 3d test with quadratic functions + * case 4: 3d test with sinusoidal functions + * case 5: 3d test with sinusoidal functions, no variations with z, non homogeneous Dirichlet conditions + * case 6: 3d test with sinusoidal functions, non homogeneous Dirichlet conditions + */ + +#include +#include "../Numerics/types.h" +#include "../Numerics/constants.h" +#include "./fsanalyticals.h" +#include "../shared.h" + +IssmDouble fx(IssmDouble x,IssmDouble y,IssmDouble z,int FSANALYTICAL){ /*{{{*/ + IssmDouble fx; + + switch(FSANALYTICAL){ + case 1: + return fx=fx1(x,y); + case 2: + return fx=fx2(x,y); + case 3: + return fx=fx3(x,y,z); + case 4: + return fx=fx4(x,y,z); + case 5: + return fx=fx5(x,y,z); + case 6: + return fx=fx6(x,y,z); + default: + _error_("FS analytical solution"<