Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17651)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17652)
@@ -7,5 +7,5 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 2
+//#define FSANALYTICAL 21
 
 /*Model processing*/
@@ -3039,97 +3039,4 @@
 	return Ke;
 }/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
-
-	if(element->IsFloating() || !element->IsOnBase()) return NULL;
-
-	/*If on water or not FS, skip stiffness: */
-	int approximation;
-	element->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
-
-	/*Intermediaries*/
-	bool        mainlyfloating;
-	int         i,meshtype,dim,epssize;
-	int         migration_style,point1;
-	IssmDouble  alpha2,Jdet,fraction1,fraction2;
-	IssmDouble  gllevelset,phi=1.;
-	IssmDouble *xyz_list_base = NULL;
-	Gauss*      gauss         = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&meshtype,MeshTypeEnum);
-	switch(meshtype){
-		case Mesh2DverticalEnum: dim = 2;break;
-		case Mesh3DEnum:         dim = 3;break;
-		case Mesh3DtetrasEnum:   dim = 3;break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumnodes = element->NumberofNodesPressure();
-	int numdof    = vnumnodes*dim + pnumnodes;
-
-	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
-	IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
-	IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	element->FindParam(&migration_style,GroundinglineMigrationEnum);
-	Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
-	Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
-	Input* vz_input         = NULL;
-	if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
-	Input* gllevelset_input = NULL;
-
-	/*build friction object, used later on: */
-	Friction* friction=new Friction(element,dim==3?3:1);
-
-	/*Recover portion of element that is grounded*/
-	if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
-	if(migration_style==SubelementMigration2Enum){
-		if(meshtype==Mesh2DverticalEnum) _error_("Subelement Migration 2 not implemented yet for Flowline");
-		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
-		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
-		//gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
-		gauss=element->NewGaussBase(3);
-	}
-	else{
-		gauss=element->NewGaussBase(3);
-	}
-
-	/* Start  looping on the number of gaussian points: */
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
-		if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
-		if(migration_style==SubelementMigration2Enum){
-			gllevelset_input->GetInputValue(&gllevelset, gauss);
-			if(gllevelset<0.) alpha2=0.;
-		}
-
-		this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
-		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		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)
-
-		TripleMultiply(B,dim-1,numdof,1,
-					D,dim-1,dim-1,0,
-					B,dim-1,numdof,0,
-					&Ke->values[0],1);
-	}
-
-	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
-
-	/*Clean up and return*/
-	delete gauss;
-	delete friction;
-	xDelete<IssmDouble>(xyz_list_base);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(D);
-	return Ke;
-}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/
 
@@ -3207,4 +3114,77 @@
 }/*}}}*/
 #ifdef FSANALYTICAL
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
+
+	if(element->IsFloating() || !element->IsOnBase()) return NULL;
+
+	/*If on water or not FS, skip stiffness: */
+	int approximation;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+
+	/*Intermediaries*/
+	int         i,meshtype,dim;
+	IssmDouble  alpha2,Jdet;
+	IssmDouble  x_coord,y_coord,z_coord;
+	IssmDouble *xyz_list_base = NULL;
+	Gauss*      gauss         = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2;break;
+		case Mesh3DEnum:         dim = 3;break;
+		case Mesh3DtetrasEnum:   dim = 3;break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+	int numdof    = vnumnodes*dim + pnumnodes;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
+	IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
+	Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
+	Input* vz_input         = NULL;
+	if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=element->NewGaussBase(10);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		x_coord=element->GetXcoord(gauss);
+		y_coord=element->GetYcoord(gauss);
+		if(dim==3) z_coord=element->GetZcoord(gauss);
+		else z_coord=0.;
+
+		alpha2=alpha(x_coord,y_coord,z_coord,FSANALYTICAL);
+
+		this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		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)
+
+		TripleMultiply(B,dim-1,numdof,1,
+					D,dim-1,dim-1,0,
+					B,dim-1,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(D);
+	return Ke;
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
@@ -3282,4 +3262,97 @@
 }/*}}}*/
 #else
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
+
+	if(element->IsFloating() || !element->IsOnBase()) return NULL;
+
+	/*If on water or not FS, skip stiffness: */
+	int approximation;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+
+	/*Intermediaries*/
+	bool        mainlyfloating;
+	int         i,meshtype,dim,epssize;
+	int         migration_style,point1;
+	IssmDouble  alpha2,Jdet,fraction1,fraction2;
+	IssmDouble  gllevelset,phi=1.;
+	IssmDouble *xyz_list_base = NULL;
+	Gauss*      gauss         = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2;break;
+		case Mesh3DEnum:         dim = 3;break;
+		case Mesh3DtetrasEnum:   dim = 3;break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+	int numdof    = vnumnodes*dim + pnumnodes;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
+	IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&migration_style,GroundinglineMigrationEnum);
+	Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
+	Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
+	Input* vz_input         = NULL;
+	if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
+	Input* gllevelset_input = NULL;
+
+	/*build friction object, used later on: */
+	Friction* friction=new Friction(element,dim==3?3:1);
+
+	/*Recover portion of element that is grounded*/
+	if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
+	if(migration_style==SubelementMigration2Enum){
+		if(meshtype==Mesh2DverticalEnum) _error_("Subelement Migration 2 not implemented yet for Flowline");
+		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+		//gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+		gauss=element->NewGaussBase(3);
+	}
+	else{
+		gauss=element->NewGaussBase(3);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+		if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
+		if(migration_style==SubelementMigration2Enum){
+			gllevelset_input->GetInputValue(&gllevelset, gauss);
+			if(gllevelset<0.) alpha2=0.;
+		}
+
+		this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		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)
+
+		TripleMultiply(B,dim-1,numdof,1,
+					D,dim-1,dim-1,0,
+					B,dim-1,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(D);
+	return Ke;
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp	(revision 17651)
+++ /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp	(revision 17652)
@@ -25,5 +25,5 @@
 			z=y;
 			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;
-		case 2: 
+		case 2: case 11:
 			z=y;
 			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);
@@ -39,5 +39,5 @@
 			z=y;
 			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;
-		case 8: 
+		case 8: case 18:
 			return 1.0;
 		case 9:
@@ -45,4 +45,9 @@
 		case 10:
 			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);
+		case 21:
+			z=y;
+			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);
+		case 24:
+			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);
 
 		case 101: 
@@ -62,5 +67,5 @@
 			z=y;
 			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;
-		case 2: 
+		case 2: case 11:
 			z=y;
 			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);
@@ -76,5 +81,5 @@
 			z=y;
 			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);
-		case 8: 
+		case 8: case 18:
 			return 1.0;
 		case 9: 
@@ -82,4 +87,9 @@
 		case 10:
 			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);
+		case 21:
+			z=y;
+			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);
+		case 24:
+			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); 
 
 		case 101: 
@@ -96,5 +106,5 @@
 
 	switch(testid){
-		case 1: case 2: case 7: case 8:
+		case 1: case 2: case 7: case 8: case 11: case 18: case 21:
 			return 0.;
 		case 3: 
@@ -110,4 +120,6 @@
 		case 10:
 			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);
+		case 24:
+			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);
 		default:
 			_error_("FS analytical solution"<<testid<<" not implemented yet");
@@ -115,2 +127,20 @@
 }
 /*}}}*/
+IssmDouble alpha(IssmDouble x,IssmDouble y,IssmDouble z,int testid){ /*{{{*/ 
+	IssmDouble p = 2.0;  
+	IssmDouble q = 2.0;   
+
+	switch(testid){
+		case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: case 9: case 10: case 21:
+			return 0.;
+		case 11:
+			z=y;
+			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));
+		case 18:
+			return (120.0L/109.0L)*sqrt(109)/(10*x - 3*z);
+		default:
+			_error_("FS analytical solution"<<testid<<" not implemented yet");
+	}
+}
+/*}}}*/
+
Index: /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h	(revision 17651)
+++ /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.h	(revision 17652)
@@ -11,4 +11,5 @@
 IssmDouble fy(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord, int testid);
 IssmDouble fz(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord, int testid);
+IssmDouble alpha(IssmDouble x_coord, IssmDouble y_coord, IssmDouble z_coord, int testid);
 
 #endif //ifndef _SHARED_ANALYTICALS_H_