Changeset 17652


Ignore:
Timestamp:
04/04/14 16:07:14 (11 years ago)
Author:
seroussi
Message:

NEW: added friction fonction to test analytical solutions

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r17650 r17652  
    77#include "../cores/cores.h"
    88
    9 //#define FSANALYTICAL 2
     9//#define FSANALYTICAL 21
    1010
    1111/*Model processing*/
    … …  
    30393039        return Ke;
    30403040}/*}}}*/
    3041 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
    3042 
    3043         if(element->IsFloating() || !element->IsOnBase()) return NULL;
    3044 
    3045         /*If on water or not FS, skip stiffness: */
    3046         int approximation;
    3047         element->GetInputValue(&approximation,ApproximationEnum);
    3048         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    3049 
    3050         /*Intermediaries*/
    3051         bool        mainlyfloating;
    3052         int         i,meshtype,dim,epssize;
    3053         int         migration_style,point1;
    3054         IssmDouble  alpha2,Jdet,fraction1,fraction2;
    3055         IssmDouble  gllevelset,phi=1.;
    3056         IssmDouble *xyz_list_base = NULL;
    3057         Gauss*      gauss         = NULL;
    3058 
    3059         /*Get problem dimension*/
    3060         element->FindParam(&meshtype,MeshTypeEnum);
    3061         switch(meshtype){
    3062                 case Mesh2DverticalEnum: dim = 2;break;
    3063                 case Mesh3DEnum:         dim = 3;break;
    3064                 case Mesh3DtetrasEnum:   dim = 3;break;
    3065                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    3066         }
    3067 
    3068         /*Fetch number of nodes and dof for this finite element*/
    3069         int vnumnodes = element->NumberofNodesVelocity();
    3070         int pnumnodes = element->NumberofNodesPressure();
    3071         int numdof    = vnumnodes*dim + pnumnodes;
    3072 
    3073         /*Initialize Element matrix and vectors*/
    3074         ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
    3075         IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
    3076         IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
    3077 
    3078         /*Retrieve all inputs and parameters*/
    3079         element->GetVerticesCoordinatesBase(&xyz_list_base);
    3080         element->FindParam(&migration_style,GroundinglineMigrationEnum);
    3081         Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
    3082         Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
    3083         Input* vz_input         = NULL;
    3084         if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
    3085         Input* gllevelset_input = NULL;
    3086 
    3087         /*build friction object, used later on: */
    3088         Friction* friction=new Friction(element,dim==3?3:1);
    3089 
    3090         /*Recover portion of element that is grounded*/
    3091         if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
    3092         if(migration_style==SubelementMigration2Enum){
    3093                 if(meshtype==Mesh2DverticalEnum) _error_("Subelement Migration 2 not implemented yet for Flowline");
    3094                 gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    3095                 element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    3096                 //gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
    3097                 gauss=element->NewGaussBase(3);
    3098         }
    3099         else{
    3100                 gauss=element->NewGaussBase(3);
    3101         }
    3102 
    3103         /* Start  looping on the number of gaussian points: */
    3104         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3105                 gauss->GaussPoint(ig);
    3106 
    3107                 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    3108                 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
    3109                 if(migration_style==SubelementMigration2Enum){
    3110                         gllevelset_input->GetInputValue(&gllevelset, gauss);
    3111                         if(gllevelset<0.) alpha2=0.;
    3112                 }
    3113 
    3114                 this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
    3115                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    3116                 for(int i=0;i<dim-1;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
    3117 
    3118                 TripleMultiply(B,dim-1,numdof,1,
    3119                                         D,dim-1,dim-1,0,
    3120                                         B,dim-1,numdof,0,
    3121                                         &Ke->values[0],1);
    3122         }
    3123 
    3124         /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
    3125 
    3126         /*Clean up and return*/
    3127         delete gauss;
    3128         delete friction;
    3129         xDelete<IssmDouble>(xyz_list_base);
    3130         xDelete<IssmDouble>(B);
    3131         xDelete<IssmDouble>(D);
    3132         return Ke;
    3133 }/*}}}*/
    31343041ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/
    31353042
    … …  
    32073114}/*}}}*/
    32083115#ifdef FSANALYTICAL
     3116ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
     3117
     3118        if(element->IsFloating() || !element->IsOnBase()) return NULL;
     3119
     3120        /*If on water or not FS, skip stiffness: */
     3121        int approximation;
     3122        element->GetInputValue(&approximation,ApproximationEnum);
     3123        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     3124
     3125        /*Intermediaries*/
     3126        int         i,meshtype,dim;
     3127        IssmDouble  alpha2,Jdet;
     3128        IssmDouble  x_coord,y_coord,z_coord;
     3129        IssmDouble *xyz_list_base = NULL;
     3130        Gauss*      gauss         = NULL;
     3131
     3132        /*Get problem dimension*/
     3133        element->FindParam(&meshtype,MeshTypeEnum);
     3134        switch(meshtype){
     3135                case Mesh2DverticalEnum: dim = 2;break;
     3136                case Mesh3DEnum:         dim = 3;break;
     3137                case Mesh3DtetrasEnum:   dim = 3;break;
     3138                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     3139        }
     3140
     3141        /*Fetch number of nodes and dof for this finite element*/
     3142        int vnumnodes = element->NumberofNodesVelocity();
     3143        int pnumnodes = element->NumberofNodesPressure();
     3144        int numdof    = vnumnodes*dim + pnumnodes;
     3145
     3146        /*Initialize Element matrix and vectors*/
     3147        ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
     3148        IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
     3149        IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
     3150
     3151        /*Retrieve all inputs and parameters*/
     3152        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3153        Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
     3154        Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
     3155        Input* vz_input         = NULL;
     3156        if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
     3157
     3158        /* Start  looping on the number of gaussian points: */
     3159        gauss=element->NewGaussBase(10);
     3160        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3161                gauss->GaussPoint(ig);
     3162
     3163                x_coord=element->GetXcoord(gauss);
     3164                y_coord=element->GetYcoord(gauss);
     3165                if(dim==3) z_coord=element->GetZcoord(gauss);
     3166                else z_coord=0.;
     3167
     3168                alpha2=alpha(x_coord,y_coord,z_coord,FSANALYTICAL);
     3169
     3170                this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
     3171                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3172                for(int i=0;i<dim-1;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
     3173
     3174                TripleMultiply(B,dim-1,numdof,1,
     3175                                        D,dim-1,dim-1,0,
     3176                                        B,dim-1,numdof,0,
     3177                                        &Ke->values[0],1);
     3178        }
     3179
     3180        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     3181
     3182        /*Clean up and return*/
     3183        delete gauss;
     3184        xDelete<IssmDouble>(xyz_list_base);
     3185        xDelete<IssmDouble>(B);
     3186        xDelete<IssmDouble>(D);
     3187        return Ke;
     3188}/*}}}*/
    32093189ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
    32103190
    … …  
    32823262}/*}}}*/
    32833263#else
     3264ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
     3265
     3266        if(element->IsFloating() || !element->IsOnBase()) return NULL;
     3267
     3268        /*If on water or not FS, skip stiffness: */
     3269        int approximation;
     3270        element->GetInputValue(&approximation,ApproximationEnum);
     3271        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     3272
     3273        /*Intermediaries*/
     3274        bool        mainlyfloating;
     3275        int         i,meshtype,dim,epssize;
     3276        int         migration_style,point1;
     3277        IssmDouble  alpha2,Jdet,fraction1,fraction2;
     3278        IssmDouble  gllevelset,phi=1.;
     3279        IssmDouble *xyz_list_base = NULL;
     3280        Gauss*      gauss         = NULL;
     3281
     3282        /*Get problem dimension*/
     3283        element->FindParam(&meshtype,MeshTypeEnum);
     3284        switch(meshtype){
     3285                case Mesh2DverticalEnum: dim = 2;break;
     3286                case Mesh3DEnum:         dim = 3;break;
     3287                case Mesh3DtetrasEnum:   dim = 3;break;
     3288                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     3289        }
     3290
     3291        /*Fetch number of nodes and dof for this finite element*/
     3292        int vnumnodes = element->NumberofNodesVelocity();
     3293        int pnumnodes = element->NumberofNodesPressure();
     3294        int numdof    = vnumnodes*dim + pnumnodes;
     3295
     3296        /*Initialize Element matrix and vectors*/
     3297        ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
     3298        IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
     3299        IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
     3300
     3301        /*Retrieve all inputs and parameters*/
     3302        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3303        element->FindParam(&migration_style,GroundinglineMigrationEnum);
     3304        Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
     3305        Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
     3306        Input* vz_input         = NULL;
     3307        if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
     3308        Input* gllevelset_input = NULL;
     3309
     3310        /*build friction object, used later on: */
     3311        Friction* friction=new Friction(element,dim==3?3:1);
     3312
     3313        /*Recover portion of element that is grounded*/
     3314        if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
     3315        if(migration_style==SubelementMigration2Enum){
     3316                if(meshtype==Mesh2DverticalEnum) _error_("Subelement Migration 2 not implemented yet for Flowline");
     3317                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
     3318                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     3319                //gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
     3320                gauss=element->NewGaussBase(3);
     3321        }
     3322        else{
     3323                gauss=element->NewGaussBase(3);
     3324        }
     3325
     3326        /* Start  looping on the number of gaussian points: */
     3327        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3328                gauss->GaussPoint(ig);
     3329
     3330                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
     3331                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
     3332                if(migration_style==SubelementMigration2Enum){
     3333                        gllevelset_input->GetInputValue(&gllevelset, gauss);
     3334                        if(gllevelset<0.) alpha2=0.;
     3335                }
     3336
     3337                this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
     3338                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3339                for(int i=0;i<dim-1;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
     3340
     3341                TripleMultiply(B,dim-1,numdof,1,
     3342                                        D,dim-1,dim-1,0,
     3343                                        B,dim-1,numdof,0,
     3344                                        &Ke->values[0],1);
     3345        }
     3346
     3347        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     3348
     3349        /*Clean up and return*/
     3350        delete gauss;
     3351        delete friction;
     3352        xDelete<IssmDouble>(xyz_list_base);
     3353        xDelete<IssmDouble>(B);
     3354        xDelete<IssmDouble>(D);
     3355        return Ke;
     3356}/*}}}*/
    32843357ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
    32853358

  • issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp

    r17404 r17652  
    2525                        z=y;
    2626                        return  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;
    27                 case 2:
     27                case 2: case 11:
    2828                        z=y;
    2929                        return 10*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 3)*pow(sin(PI*p*z), 2)*cos(PI*p*z) - 2*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 3)*pow(cos(PI*p*z), 3) - 6*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*pow(sin(PI*p*z), 2)*pow(cos(PI*p*x), 2)*cos(PI*p*z) + PI*q*sin(PI*q*z)*cos(PI*q*x);
    … …  
    3939                        z=y;
    4040                        return 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) + 1;
    41                 case 8:
     41                case 8: case 18:
    4242                        return 1.0;
    4343                case 9:
    … …  
    4545                case 10:
    4646                        return 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*pow(z, a)*sin(PI*q*y)*cos(PI*q*x);
     47                case 21:
     48                        z=y;
     49                        return z - 4.0L/3.0L*pow(2, 1.0L/3.0L)*(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1))*(4*(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1))*(pow(x, 2)*z*(z - 1)*(2*z - 1) + 4*x*z*(x - 1)*(z - 1)*(2*z - 1) + z*pow(x - 1, 2)*(z - 1)*(2*z - 1)) + 4*(2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - 3*pow(z, 2)*(6*pow(x, 2) - 6*x + 1) + z*(6*pow(x, 2) - 6*x + 1))*(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x)) + (2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1))*(2*pow(x, 2)*z*(x - 1)*(z - 1) + pow(x, 2)*z*(x - 1)*(2*z - 1) + pow(x, 2)*(x - 1)*(z - 1)*(2*z - 1) + 2*x*z*pow(x - 1, 2)*(z - 1) + x*z*pow(x - 1, 2)*(2*z - 1) + x*pow(x - 1, 2)*(z - 1)*(2*z - 1) - 3*pow(z, 4)*(2*x - 1) + 6*pow(z, 3)*(2*x - 1) - 3*pow(z, 2)*(2*x - 1)))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 4.0L/3.0L) - 1.0L/3.0L*pow(2, 1.0L/3.0L)*(4*(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1))*(2*pow(x, 2)*z*(x - 1)*(z - 1) + pow(x, 2)*z*(x - 1)*(2*z - 1) + pow(x, 2)*(x - 1)*(z - 1)*(2*z - 1) + 2*x*z*pow(x - 1, 2)*(z - 1) + x*z*pow(x - 1, 2)*(2*z - 1) + x*pow(x - 1, 2)*(z - 1)*(2*z - 1)) + 4*(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x))*(2*pow(x, 3) - 3*pow(x, 2) + x + 6*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 6*z*(2*pow(x, 3) - 3*pow(x, 2) + x)) + (2*pow(x, 2)*z*pow(x - 1, 2) + 2*pow(x, 2)*pow(x - 1, 2)*(z - 1) + pow(x, 2)*pow(x - 1, 2)*(2*z - 1) - 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) + 3*pow(z, 2)*(6*pow(x, 2) - 6*x + 1) - z*(6*pow(x, 2) - 6*x + 1))*(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1)))*(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 4.0L/3.0L) + 2*pow(2, 1.0L/3.0L)*(pow(x, 2)*z*(z - 1)*(2*z - 1) + 4*x*z*(x - 1)*(z - 1)*(2*z - 1) + z*pow(x - 1, 2)*(z - 1)*(2*z - 1))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 1.0L/3.0L) - 1.0L/2.0L + pow(2, 1.0L/3.0L)*(2*pow(x, 2)*z*pow(x - 1, 2) + 2*pow(x, 2)*pow(x - 1, 2)*(z - 1) + pow(x, 2)*pow(x - 1, 2)*(2*z - 1) - 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) + 3*pow(z, 2)*(6*pow(x, 2) - 6*x + 1) - z*(6*pow(x, 2) - 6*x + 1))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 1.0L/3.0L);
     50                case 24:
     51                        return 2*pow(2, 1.0L/3.0L)*pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x)/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L) - 2.0L/3.0L*pow(2, 1.0L/3.0L)*PI*p*(12*pow(PI, 3)*pow(p, 3)*sin(PI*p*x)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2)*cos(PI*p*x) - (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*sin(PI*p*x)*sin(PI*p*z) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*z)*cos(PI*p*y)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*cos(PI*p*x)*cos(PI*p*y) + pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*cos(PI*p*x)*cos(PI*p*z)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*sin(PI*p*y) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*y)*cos(PI*p*z)))*sin(PI*p*x)*sin(PI*p*y)*sin(PI*p*z)/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) + PI*q*sin(PI*q*y)*sin(PI*q*z)*cos(PI*q*x) - 1.0L/6.0L*pow(2, 1.0L/3.0L)*(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(12*pow(PI, 3)*pow(p, 3)*pow(sin(PI*p*x), 2)*sin(PI*p*y)*pow(sin(PI*p*z), 2)*cos(PI*p*y) + (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*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)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*cos(PI*p*x)*cos(PI*p*y) - 2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*cos(PI*p*y)*cos(PI*p*z)))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) - 1.0L/6.0L*pow(2, 1.0L/3.0L)*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(12*pow(PI, 3)*pow(p, 3)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*sin(PI*p*z)*cos(PI*p*z) + (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(-2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*cos(PI*p*y)*cos(PI*p*z) + pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*cos(PI*p*x)*cos(PI*p*z)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(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)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x)))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) + pow(2, 1.0L/3.0L)*(2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L);
    4752
    4853                case 101:
    … …  
    6267                        z=y;
    6368                        return -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;
    64                 case 2:
     69                case 2: case 11:
    6570                        z=y;
    6671                        return -10*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*z), 3)*cos(PI*p*x) + 6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*sin(PI*p*z)*cos(PI*p*x)*pow(cos(PI*p*z), 2) + 2*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*z), 3)*pow(cos(PI*p*x), 3) + PI*q*sin(PI*q*x)*cos(PI*q*z);
    … …  
    7681                        z=y;
    7782                        return -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) - 4*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);
    78                 case 8:
     83                case 8: case 18:
    7984                        return 1.0;
    8085                case 9:
    … …  
    8287                case 10:
    8388                        return -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*pow(z, a)*sin(PI*q*x)*cos(PI*q*y);
     89                case 21:
     90                        z=y;
     91                        return x + (4.0L/3.0L)*pow(2, 1.0L/3.0L)*(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x))*(4*(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1))*(2*pow(x, 2)*z*(x - 1)*(z - 1) + pow(x, 2)*z*(x - 1)*(2*z - 1) + pow(x, 2)*(x - 1)*(z - 1)*(2*z - 1) + 2*x*z*pow(x - 1, 2)*(z - 1) + x*z*pow(x - 1, 2)*(2*z - 1) + x*pow(x - 1, 2)*(z - 1)*(2*z - 1)) + 4*(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x))*(2*pow(x, 3) - 3*pow(x, 2) + x + 6*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 6*z*(2*pow(x, 3) - 3*pow(x, 2) + x)) + (2*pow(x, 2)*z*pow(x - 1, 2) + 2*pow(x, 2)*pow(x - 1, 2)*(z - 1) + pow(x, 2)*pow(x - 1, 2)*(2*z - 1) - 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) + 3*pow(z, 2)*(6*pow(x, 2) - 6*x + 1) - z*(6*pow(x, 2) - 6*x + 1))*(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1)))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 4.0L/3.0L) - 1.0L/3.0L*pow(2, 1.0L/3.0L)*(4*(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1))*(pow(x, 2)*z*(z - 1)*(2*z - 1) + 4*x*z*(x - 1)*(z - 1)*(2*z - 1) + z*pow(x - 1, 2)*(z - 1)*(2*z - 1)) + 4*(2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - 3*pow(z, 2)*(6*pow(x, 2) - 6*x + 1) + z*(6*pow(x, 2) - 6*x + 1))*(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x)) + (2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1))*(2*pow(x, 2)*z*(x - 1)*(z - 1) + pow(x, 2)*z*(x - 1)*(2*z - 1) + pow(x, 2)*(x - 1)*(z - 1)*(2*z - 1) + 2*x*z*pow(x - 1, 2)*(z - 1) + x*z*pow(x - 1, 2)*(2*z - 1) + x*pow(x - 1, 2)*(z - 1)*(2*z - 1) - 3*pow(z, 4)*(2*x - 1) + 6*pow(z, 3)*(2*x - 1) - 3*pow(z, 2)*(2*x - 1)))*(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 4.0L/3.0L) - 1.0L/2.0L - 2*pow(2, 1.0L/3.0L)*(2*pow(x, 3) - 3*pow(x, 2) + x + 6*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 6*z*(2*pow(x, 3) - 3*pow(x, 2) + x))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 1.0L/3.0L) + pow(2, 1.0L/3.0L)*(2*pow(x, 2)*z*(x - 1)*(z - 1) + pow(x, 2)*z*(x - 1)*(2*z - 1) + pow(x, 2)*(x - 1)*(z - 1)*(2*z - 1) + 2*x*z*pow(x - 1, 2)*(z - 1) + x*z*pow(x - 1, 2)*(2*z - 1) + x*pow(x - 1, 2)*(z - 1)*(2*z - 1) - 3*pow(z, 4)*(2*x - 1) + 6*pow(z, 3)*(2*x - 1) - 3*pow(z, 2)*(2*x - 1))/pow(4*pow(pow(x, 2)*z*(x - 1)*(z - 1)*(2*z - 1) + x*z*pow(x - 1, 2)*(z - 1)*(2*z - 1), 2) + 4*pow(2*pow(z, 3)*(2*pow(x, 3) - 3*pow(x, 2) + x) - 3*pow(z, 2)*(2*pow(x, 3) - 3*pow(x, 2) + x) + z*(2*pow(x, 3) - 3*pow(x, 2) + x), 2) + (1.0L/2.0L)*pow(2*pow(x, 2)*z*pow(x - 1, 2)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*z - 1) + pow(x, 2)*pow(x - 1, 2)*(z - 1)*(2*z - 1) - pow(z, 4)*(6*pow(x, 2) - 6*x + 1) + 2*pow(z, 3)*(6*pow(x, 2) - 6*x + 1) - pow(z, 2)*(6*pow(x, 2) - 6*x + 1), 2), 1.0L/3.0L);
     92                case 24:
     93                        return -pow(2, 1.0L/3.0L)*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*z)*cos(PI*p*y)/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L) + (1.0L/3.0L)*pow(2, 1.0L/3.0L)*PI*p*(12*pow(PI, 3)*pow(p, 3)*pow(sin(PI*p*x), 2)*sin(PI*p*y)*pow(sin(PI*p*z), 2)*cos(PI*p*y) + (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*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)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*cos(PI*p*x)*cos(PI*p*y) - 2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*cos(PI*p*y)*cos(PI*p*z)))*sin(PI*p*x)*sin(PI*p*y)*sin(PI*p*z)/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) + PI*q*sin(PI*q*x)*sin(PI*q*z)*cos(PI*q*y) - 1.0L/6.0L*pow(2, 1.0L/3.0L)*(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(12*pow(PI, 3)*pow(p, 3)*sin(PI*p*x)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2)*cos(PI*p*x) - (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*sin(PI*p*x)*sin(PI*p*z) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*z)*cos(PI*p*y)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*cos(PI*p*x)*cos(PI*p*y) + pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*cos(PI*p*x)*cos(PI*p*z)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*sin(PI*p*y) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*y)*cos(PI*p*z)))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) - 1.0L/6.0L*pow(2, 1.0L/3.0L)*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(12*pow(PI, 3)*pow(p, 3)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*sin(PI*p*z)*cos(PI*p*z) + (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(-2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*cos(PI*p*y)*cos(PI*p*z) + pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*cos(PI*p*x)*cos(PI*p*z)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(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)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x)))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) - 1.0L/2.0L*pow(2, 1.0L/3.0L)*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*sin(PI*p*x)*sin(PI*p*z) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*z)*cos(PI*p*y))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L) - 1.0L/2.0L*pow(2, 1.0L/3.0L)*(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))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L);
    8494
    8595                case 101:
    … …  
    96106
    97107        switch(testid){
    98                 case 1: case 2: case 7: case 8:
     108                case 1: case 2: case 7: case 8: case 11: case 18: case 21:
    99109                        return 0.;
    100110                case 3:
    … …  
    110120                case 10:
    111121                        return a*pow(z, a - 1)*(sin(PI*q*x)*sin(PI*q*y) + 1) - 2*pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*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);
     122                case 24:
     123                        return -pow(2, 1.0L/3.0L)*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*y)*cos(PI*p*z)/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L) + (1.0L/3.0L)*pow(2, 1.0L/3.0L)*PI*p*(12*pow(PI, 3)*pow(p, 3)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*sin(PI*p*z)*cos(PI*p*z) + (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(-2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*cos(PI*p*y)*cos(PI*p*z) + pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*cos(PI*p*x)*cos(PI*p*z)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(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)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x)))*sin(PI*p*x)*sin(PI*p*y)*sin(PI*p*z)/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) + PI*q*sin(PI*q*x)*sin(PI*q*y)*cos(PI*q*z) - 1.0L/6.0L*pow(2, 1.0L/3.0L)*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(12*pow(PI, 3)*pow(p, 3)*pow(sin(PI*p*x), 2)*sin(PI*p*y)*pow(sin(PI*p*z), 2)*cos(PI*p*y) + (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*sin(PI*p*y)*sin(PI*p*z) - pow(PI, 2)*pow(p, 2)*sin(PI*p*y)*sin(PI*p*z)*cos(PI*p*x)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*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)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*cos(PI*p*x)*cos(PI*p*y) - 2*pow(PI, 2)*pow(p, 2)*(cos(PI*p*x) - 1)*cos(PI*p*y)*cos(PI*p*z)))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) - 1.0L/6.0L*pow(2, 1.0L/3.0L)*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(12*pow(PI, 3)*pow(p, 3)*sin(PI*p*x)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2)*cos(PI*p*x) - (-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x))*(pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*sin(PI*p*x)*sin(PI*p*z) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*z)*cos(PI*p*y)) + (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*cos(PI*p*x)*cos(PI*p*y) + pow(PI, 2)*pow(p, 2)*(cos(PI*p*y) - 1)*cos(PI*p*x)*cos(PI*p*z)) - (pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z))*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*sin(PI*p*y) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*y)*cos(PI*p*z)))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 4.0L/3.0L) - 1.0L/2.0L*pow(2, 1.0L/3.0L)*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*sin(PI*p*y) - 2*pow(PI, 2)*pow(p, 2)*sin(PI*p*x)*sin(PI*p*y)*cos(PI*p*z))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L) - 1.0L/2.0L*pow(2, 1.0L/3.0L)*(pow(PI, 3)*pow(p, 3)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*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))/pow(6*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(sin(PI*p*y), 2)*pow(sin(PI*p*z), 2) + (1.0L/2.0L)*pow(-2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*z)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*z)*cos(PI*p*x), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*x)*cos(PI*p*y) + PI*p*(cos(PI*p*y) - 1)*sin(PI*p*x)*cos(PI*p*z), 2) + (1.0L/2.0L)*pow(pow(PI, 2)*pow(p, 2)*(cos(PI*p*z)/(PI*p) - 1/(PI*p))*sin(PI*p*y)*cos(PI*p*x) - 2*PI*p*(cos(PI*p*x) - 1)*sin(PI*p*y)*cos(PI*p*z), 2), 1.0L/3.0L);
    112124                default:
    113125                        _error_("FS analytical solution"<<testid<<" not implemented yet");
    … …  
    115127}
    116128/*}}}*/
     129IssmDouble alpha(IssmDouble x,IssmDouble y,IssmDouble z,int testid){ /*{{{*/
     130        IssmDouble p = 2.0; 
     131        IssmDouble q = 2.0;   
     132
     133        switch(testid){
     134                case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: case 9: case 10: case 21:
     135                        return 0.;
     136                case 11:
     137                        z=y;
     138                        return -10*PI*p*sqrt((9.0L/25.0L)*pow(PI, 2)*pow(p, 2)*pow(cos(2*PI*p*x), 2) + 1)*(9*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*x), 2)*pow(cos(2*PI*p*x), 2)*pow(cos(PI*p*z), 2) - 9*pow(PI, 2)*pow(p, 2)*pow(sin(PI*p*z), 2)*pow(cos(PI*p*x), 2)*pow(cos(2*PI*p*x), 2) + 90*PI*p*sin(PI*p*x)*sin(PI*p*z)*cos(PI*p*x)*cos(2*PI*p*x)*cos(PI*p*z) - 25*pow(sin(PI*p*x), 2)*pow(cos(PI*p*z), 2) + 25*pow(sin(PI*p*z), 2)*pow(cos(PI*p*x), 2))/((9*pow(PI, 2)*pow(p, 2)*pow(cos(2*PI*p*x), 2) + 25)*(3*PI*p*sin(PI*p*z)*cos(PI*p*x)*cos(2*PI*p*x) - 5*sin(PI*p*x)*cos(PI*p*z))*sin(PI*p*x)*sin(PI*p*z));
     139                case 18:
     140                        return (120.0L/109.0L)*sqrt(109)/(10*x - 3*z);
     141                default:
     142                        _error_("FS analytical solution"<<testid<<" not implemented yet");
     143        }
     144}
     145/*}}}*/
     146

  • issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h

    r17360 r17652  
    1111IssmDouble fy(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord, int testid);
    1212IssmDouble fz(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord, int testid);
     13IssmDouble alpha(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord, int testid);
    1314
    1415#endif //ifndef _SHARED_ANALYTICALS_H_
Note: See TracChangeset for help on using the changeset viewer.