Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20648) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20649) @@ -1411,5 +1411,5 @@ /*Call inputs method*/ - if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); + if(IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); } @@ -1531,4 +1531,10 @@ name==MaterialsRheologyBbarEnum || name==MaterialsRheologyNEnum || + name==MaterialsRheologyKoEnum || + name==MaterialsRheologyKobarEnum || + name==MaterialsRheologyEcEnum || + name==MaterialsRheologyEcbarEnum || + name==MaterialsRheologyEsEnum || + name==MaterialsRheologyEsbarEnum || name==SealevelEnum || name==SealevelEustaticEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20648) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20649) @@ -2251,4 +2251,7 @@ break; case MatbaflEnum: + this->InputDepthAverageAtBase(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum); + this->InputDepthAverageAtBase(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum); + this->InputDepthAverageAtBase(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum); break; default: Index: /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp (revision 20648) +++ /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp (revision 20649) @@ -275,25 +275,141 @@ void Matbafl::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ - _assert_(dim==3); + /*Intermediaries*/ + IssmDouble vx,vy,vz; + IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble ko,Ec,Es; + + /*Get velocity derivatives in all directions*/ + _assert_(dim>1); + _assert_(vx_input); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + _assert_(vy_input); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + if(dim==3){ + _assert_(vz_input); + vz_input->GetInputValue(&vz,gauss); + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + } + else{ + vz = 0.; + dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.; + } + + /*Get enhancement factors and ko*/ + Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input); + Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input); + Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input); + ec_input->GetInputValue(&Ec,gauss); + es_input->GetInputValue(&Es,gauss); + ko_input->GetInputValue(&ko,gauss); + + /*Compute viscosity*/ + *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); +} +/*}}}*/ +void Matbafl::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ +void Matbafl::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble vx,vy,vz; + IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble ko,Ec,Es; + + /*Get velocity derivatives in all directions*/ + _assert_(dim==2 || dim==3); + _assert_(vx_input); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + if(dim==3){ + _assert_(vy_input); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + } + else{ + vy = 0.; + dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.; + } + vz = 0.; + dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.; + + /*Get enhancement factors and ko*/ + Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input); + Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input); + Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input); + ec_input->GetInputValue(&Ec,gauss); + es_input->GetInputValue(&Es,gauss); + ko_input->GetInputValue(&ko,gauss); + + /*Compute viscosity*/ + *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); +}/*}}}*/ +void Matbafl::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void Matbafl::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void Matbafl::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + /*Intermediaries*/ + IssmDouble vx,vy,vz; + IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble ko,Ec,Es; + + /*Get velocity derivatives in all directions*/ + _assert_(dim==1 || dim==2); + _assert_(vx_input); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + if(dim==2){ + _assert_(vy_input); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + } + else{ + vy = 0.; + dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.; + } + vz = 0.; + dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.; + + /*Get enhancement factors and ko*/ + Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input); + Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input); + Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input); + ec_input->GetInputValue(&Ec,gauss); + es_input->GetInputValue(&Es,gauss); + ko_input->GetInputValue(&ko,gauss); + + /*Compute viscosity*/ + *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); +}/*}}}*/ +void Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void Matbafl::ResetHooks(){/*{{{*/ + + this->element=NULL; + + /*Get Element type*/ + this->helement->reset(); + +} +/*}}}*/ +IssmDouble Matbafl::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/ /*Intermediaries*/ IssmDouble viscosity; - IssmDouble vx,vy,vz; - IssmDouble dvx[3],dvy[3],dvz[3]; IssmDouble vorticity[3],vorticity_norm; IssmDouble nrsp[3],nrsp_norm; IssmDouble eps[3][3],epso; IssmDouble epsprime[3],epsprime_norm; - IssmDouble ko,E,Ec,Es,lambdas; - - /*Get velocity derivatives in all directions*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + IssmDouble E,lambdas; /*Create vorticity vector*/ + _assert_(dvx && dvy && dvz); vorticity[0] = dvz[1] - dvy[2]; vorticity[1] = dvx[2] - dvz[0]; @@ -377,13 +493,5 @@ /*Get total enhancement factor E(lambdas)*/ - Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input); - Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input); - ec_input->GetInputValue(&Ec,gauss); - es_input->GetInputValue(&Es,gauss); E = Ec + (Es-Ec)*lambdas*lambdas; - - /*Get flow law parameter*/ - Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input); - ko_input->GetInputValue(&ko,gauss); /*Compute viscosity*/ @@ -394,33 +502,5 @@ /*Assign output pointer*/ - *pviscosity=viscosity; -} -/*}}}*/ -void Matbafl::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); -}/*}}}*/ -void Matbafl::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - _error_("not implemented yet"); -}/*}}}*/ -void Matbafl::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void Matbafl::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void Matbafl::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void Matbafl::ResetHooks(){/*{{{*/ - - this->element=NULL; - - /*Get Element type*/ - this->helement->reset(); - -} -/*}}}*/ + return viscosity; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matbafl.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matbafl.h (revision 20648) +++ /issm/trunk-jpl/src/c/classes/Materials/Matbafl.h (revision 20649) @@ -82,4 +82,5 @@ void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); /*}}}*/ + IssmDouble GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz); }; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 20648) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 20649) @@ -84,4 +84,15 @@ iomodel->FetchDataToInput(elements,MaterialsRheologyEsEnum); for(i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matbafl(i+1,i,iomodel)); + switch(iomodel->domaindim){ + case 2: + elements->InputDuplicate(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum); + elements->InputDuplicate(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum); + elements->InputDuplicate(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum); + break; + case 3: + break; + default: + _error_("Mesh not supported yet"); + } break; default: Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20648) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20649) @@ -233,6 +233,9 @@ MaterialsRheologyNEnum, MaterialsRheologyKoEnum, + MaterialsRheologyKobarEnum, MaterialsRheologyEcEnum, + MaterialsRheologyEcbarEnum, MaterialsRheologyEsEnum, + MaterialsRheologyEsbarEnum, DamageIsdamageEnum, DamageDEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20648) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20649) @@ -237,6 +237,9 @@ case MaterialsRheologyNEnum : return "MaterialsRheologyN"; case MaterialsRheologyKoEnum : return "MaterialsRheologyKo"; + case MaterialsRheologyKobarEnum : return "MaterialsRheologyKobar"; case MaterialsRheologyEcEnum : return "MaterialsRheologyEc"; + case MaterialsRheologyEcbarEnum : return "MaterialsRheologyEcbar"; case MaterialsRheologyEsEnum : return "MaterialsRheologyEs"; + case MaterialsRheologyEsbarEnum : return "MaterialsRheologyEsbar"; case DamageIsdamageEnum : return "DamageIsdamage"; case DamageDEnum : return "DamageD"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20648) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20649) @@ -240,6 +240,9 @@ else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; else if (strcmp(name,"MaterialsRheologyKo")==0) return MaterialsRheologyKoEnum; + else if (strcmp(name,"MaterialsRheologyKobar")==0) return MaterialsRheologyKobarEnum; else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum; + else if (strcmp(name,"MaterialsRheologyEcbar")==0) return MaterialsRheologyEcbarEnum; else if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum; + else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum; else if (strcmp(name,"DamageIsdamage")==0) return DamageIsdamageEnum; else if (strcmp(name,"DamageD")==0) return DamageDEnum; @@ -257,11 +260,11 @@ else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum; else if (strcmp(name,"DamageMaxiter")==0) return DamageMaxiterEnum; - else if (strcmp(name,"DamageSpcdamage")==0) return DamageSpcdamageEnum; - else if (strcmp(name,"DamageMaxDamage")==0) return DamageMaxDamageEnum; - else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum; + if (strcmp(name,"DamageSpcdamage")==0) return DamageSpcdamageEnum; + else if (strcmp(name,"DamageMaxDamage")==0) return DamageMaxDamageEnum; + else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum; + else if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum; else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum; else if (strcmp(name,"Damage")==0) return DamageEnum; @@ -380,11 +383,11 @@ else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum; else if (strcmp(name,"TimesteppingInterpForcings")==0) return TimesteppingInterpForcingsEnum; - else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum; - else if (strcmp(name,"TransientIscoupler")==0) return TransientIscouplerEnum; - else if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; + if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum; + else if (strcmp(name,"TransientIscoupler")==0) return TransientIscouplerEnum; + else if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum; + else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum; else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; - else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; - else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; - else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; + if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; + else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; + else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; + else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum; else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; @@ -626,11 +629,11 @@ else if (strcmp(name,"MassconName")==0) return MassconNameEnum; else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum; - else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; - else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; - else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum; + if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; + else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; + else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum; + else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum; else if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum; else if (strcmp(name,"MassconaxpbyNamey")==0) return MassconaxpbyNameyEnum; @@ -749,11 +752,11 @@ else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; - else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; - else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; - else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"StrainRate")==0) return StrainRateEnum; + if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; + else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; + else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; + else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; - else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; - else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum; - else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum; + if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; + else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum; + else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum; + else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum; else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum; else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; @@ -995,11 +998,11 @@ else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum; else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; - else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; - else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum; - else if (strcmp(name,"Regular")==0) return RegularEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"Scaled")==0) return ScaledEnum; + if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; + else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum; + else if (strcmp(name,"Regular")==0) return RegularEnum; + else if (strcmp(name,"Scaled")==0) return ScaledEnum; else if (strcmp(name,"Separate")==0) return SeparateEnum; else if (strcmp(name,"Sset")==0) return SsetEnum; Index: /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 20648) +++ /issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 20649) @@ -8,4 +8,6 @@ * case 5: 3d test with sinusoidal functions, no variations with z, non homogeneous Dirichlet conditions * case 6: 3d test with sinusoidal functions, non homogeneous Dirichlet conditions + * + * case 201: 3d test quadratic functions, BAFL flow law */ @@ -71,4 +73,7 @@ case 101: return 4*pow(x, 2)*y*z*pow(x - 1, 2)*(z - 1) + 8*pow(x, 2)*y*z*(y - 1)*(2*y - 1)*(z - 1) + 2*pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1) + 4*pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 2*pow(x, 2)*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 12*x*pow(y, 2)*z*(x - 1)*(y - 1)*(z - 1) - 6*x*pow(y, 2)*z*(2*x - 1)*(y - 1)*(z - 1) - 12*x*y*z*(x - 1)*pow(y - 1, 2)*(z - 1) + 32*x*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 6*x*y*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 6*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 8*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 6*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1); + case 201: + return pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)/pow((1.0L/6.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (1.0L/6.0L)*pow(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (4.0L/3.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (4.0L/3.0L)*pow(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (1.0L/6.0L)*pow(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2), 1.0L/3.0L) + y - 1.0L/18.0L*(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))*(2*pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1)) + 2*x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)) + 8*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1) + x*pow(y, 2)*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*y*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)) + 8*(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))*(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1) + pow(x, 2)*y*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + x*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1)) + (2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1) + 2*pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2) - 2*x*pow(y, 2)*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) - pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)))/pow((1.0L/6.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (1.0L/6.0L)*pow(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (4.0L/3.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (4.0L/3.0L)*pow(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (1.0L/6.0L)*pow(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2), 4.0L/3.0L) - 2.0L/9.0L*(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))*((x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2) + x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2) + 2*x*pow(y, 2)*(x - 1)*pow(y - 1, 2)*(z - 1) + x*pow(y, 2)*(2*x - 1)*pow(y - 1, 2)*(z - 1) + pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)) + 2*(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))*(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1) + pow(x, 2)*y*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + x*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1)) + 8*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(2*x*pow(y, 2)*z*(x - 1)*(y - 1)*(z - 1) + x*pow(y, 2)*z*(2*x - 1)*(y - 1)*(z - 1) + 2*x*y*z*(x - 1)*pow(y - 1, 2)*(z - 1) + x*y*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) + pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)) + 8*(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))*(pow(x, 2)*y*z*(y - 1)*(2*y - 1)*(z - 1) + 4*x*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1)) + 2*(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(2*pow(x, 2)*y*z*(x - 1)*(y - 1)*(z - 1) + pow(x, 2)*y*z*(x - 1)*(2*y - 1)*(z - 1) + pow(x, 2)*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*pow(y - 1, 2)*(z - 1) + 2*x*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + x*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1)))/pow((1.0L/6.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (1.0L/6.0L)*pow(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (4.0L/3.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (4.0L/3.0L)*pow(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (1.0L/6.0L)*pow(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2), 4.0L/3.0L) + 2*(pow(x, 2)*y*z*(y - 1)*(2*y - 1)*(z - 1) + 4*x*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))/pow((1.0L/6.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (1.0L/6.0L)*pow(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (4.0L/3.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (4.0L/3.0L)*pow(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (1.0L/6.0L)*pow(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2), 1.0L/3.0L) - 1.0L/18.0L*(2*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1) + x*pow(y, 2)*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*y*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)) + (pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))*(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1) + 2*pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1)) + 8*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(z - 1) + 4*x*y*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)) + 8*(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1))*(2*pow(x, 2)*y*z*(x - 1)*(y - 1)*(z - 1) + pow(x, 2)*y*z*(x - 1)*(2*y - 1)*(z - 1) + pow(x, 2)*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + 2*x*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + x*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1)) + 2*(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))*(2*pow(x, 2)*y*z*pow(x - 1, 2)*(z - 1) + 2*pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*(y - 1)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*(y - 1)*(z - 1) - 2*x*y*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*y*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) - y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1)))*(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))/pow((1.0L/6.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (1.0L/6.0L)*pow(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (4.0L/3.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (4.0L/3.0L)*pow(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (1.0L/6.0L)*pow(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2), 4.0L/3.0L) - 1.0L/2.0L + (2*pow(x, 2)*y*z*pow(x - 1, 2)*(z - 1) + 2*pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*(y - 1)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*(y - 1)*(z - 1) - 2*x*y*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*y*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) - y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1))/pow((1.0L/6.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2) + x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (1.0L/6.0L)*pow(pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1) + pow(x, 2)*y*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (4.0L/3.0L)*pow(x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + x*y*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2) + (4.0L/3.0L)*pow(pow(x, 2)*y*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) + x*y*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1), 2) + (1.0L/6.0L)*pow(2*pow(x, 2)*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + pow(x, 2)*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) + pow(x, 2)*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 2*x*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - x*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1) - pow(y, 2)*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1), 2), 1.0L/3.0L); + default: _error_("FS analytical solution"<