source: issm/oecreview/Archive/20545-21336/ISSM-20648-20649.diff@ 21337

Last change on this file since 21337 was 21337, checked in by Mathieu Morlighem, 8 years ago

CHG: added Archive/20545-21336

File size: 61.3 KB
RevLine 
[21337]1Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
2===================================================================
3--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20648)
4+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20649)
5@@ -232,8 +232,11 @@
6 MaterialsRheologyLawEnum,
7 MaterialsRheologyNEnum,
8 MaterialsRheologyKoEnum,
9+ MaterialsRheologyKobarEnum,
10 MaterialsRheologyEcEnum,
11+ MaterialsRheologyEcbarEnum,
12 MaterialsRheologyEsEnum,
13+ MaterialsRheologyEsbarEnum,
14 DamageIsdamageEnum,
15 DamageDEnum,
16 DamageFEnum,
17Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
18===================================================================
19--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20648)
20+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20649)
21@@ -236,8 +236,11 @@
22 case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw";
23 case MaterialsRheologyNEnum : return "MaterialsRheologyN";
24 case MaterialsRheologyKoEnum : return "MaterialsRheologyKo";
25+ case MaterialsRheologyKobarEnum : return "MaterialsRheologyKobar";
26 case MaterialsRheologyEcEnum : return "MaterialsRheologyEc";
27+ case MaterialsRheologyEcbarEnum : return "MaterialsRheologyEcbar";
28 case MaterialsRheologyEsEnum : return "MaterialsRheologyEs";
29+ case MaterialsRheologyEsbarEnum : return "MaterialsRheologyEsbar";
30 case DamageIsdamageEnum : return "DamageIsdamage";
31 case DamageDEnum : return "DamageD";
32 case DamageFEnum : return "DamageF";
33Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
34===================================================================
35--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20648)
36+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20649)
37@@ -239,8 +239,11 @@
38 else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
39 else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
40 else if (strcmp(name,"MaterialsRheologyKo")==0) return MaterialsRheologyKoEnum;
41+ else if (strcmp(name,"MaterialsRheologyKobar")==0) return MaterialsRheologyKobarEnum;
42 else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum;
43+ else if (strcmp(name,"MaterialsRheologyEcbar")==0) return MaterialsRheologyEcbarEnum;
44 else if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
45+ else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
46 else if (strcmp(name,"DamageIsdamage")==0) return DamageIsdamageEnum;
47 else if (strcmp(name,"DamageD")==0) return DamageDEnum;
48 else if (strcmp(name,"DamageF")==0) return DamageFEnum;
49@@ -256,13 +259,13 @@
50 else if (strcmp(name,"DamageKappa")==0) return DamageKappaEnum;
51 else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum;
52 else if (strcmp(name,"DamageMaxiter")==0) return DamageMaxiterEnum;
53- else if (strcmp(name,"DamageSpcdamage")==0) return DamageSpcdamageEnum;
54- else if (strcmp(name,"DamageMaxDamage")==0) return DamageMaxDamageEnum;
55- else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
56 else stage=3;
57 }
58 if(stage==3){
59- if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum;
60+ if (strcmp(name,"DamageSpcdamage")==0) return DamageSpcdamageEnum;
61+ else if (strcmp(name,"DamageMaxDamage")==0) return DamageMaxDamageEnum;
62+ else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
63+ else if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum;
64 else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
65 else if (strcmp(name,"Damage")==0) return DamageEnum;
66 else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
67@@ -379,13 +382,13 @@
68 else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
69 else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
70 else if (strcmp(name,"TimesteppingInterpForcings")==0) return TimesteppingInterpForcingsEnum;
71- else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum;
72- else if (strcmp(name,"TransientIscoupler")==0) return TransientIscouplerEnum;
73- else if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum;
74 else stage=4;
75 }
76 if(stage==4){
77- if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
78+ if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum;
79+ else if (strcmp(name,"TransientIscoupler")==0) return TransientIscouplerEnum;
80+ else if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum;
81+ else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
82 else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
83 else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
84 else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
85@@ -502,13 +505,13 @@
86 else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
87 else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
88 else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
89- else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
90- else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
91- else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
92 else stage=5;
93 }
94 if(stage==5){
95- if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
96+ if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
97+ else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
98+ else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
99+ else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
100 else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
101 else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
102 else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
103@@ -625,13 +628,13 @@
104 else if (strcmp(name,"Masscon")==0) return MassconEnum;
105 else if (strcmp(name,"MassconName")==0) return MassconNameEnum;
106 else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum;
107- else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
108- else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
109- else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
110 else stage=6;
111 }
112 if(stage==6){
113- if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
114+ if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
115+ else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
116+ else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
117+ else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
118 else if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum;
119 else if (strcmp(name,"MassconaxpbyNamey")==0) return MassconaxpbyNameyEnum;
120 else if (strcmp(name,"MassconaxpbyAlpha")==0) return MassconaxpbyAlphaEnum;
121@@ -748,13 +751,13 @@
122 else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
123 else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
124 else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
125- else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
126- else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
127- else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
128 else stage=7;
129 }
130 if(stage==7){
131- if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
132+ if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
133+ else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
134+ else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
135+ else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
136 else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
137 else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
138 else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
139@@ -871,13 +874,13 @@
140 else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
141 else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
142 else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
143- else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
144- else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
145- else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
146 else stage=8;
147 }
148 if(stage==8){
149- if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
150+ if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
151+ else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
152+ else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
153+ else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
154 else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
155 else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
156 else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
157@@ -994,13 +997,13 @@
158 else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
159 else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
160 else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
161- else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
162- else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
163- else if (strcmp(name,"Regular")==0) return RegularEnum;
164 else stage=9;
165 }
166 if(stage==9){
167- if (strcmp(name,"Scaled")==0) return ScaledEnum;
168+ if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
169+ else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
170+ else if (strcmp(name,"Regular")==0) return RegularEnum;
171+ else if (strcmp(name,"Scaled")==0) return ScaledEnum;
172 else if (strcmp(name,"Separate")==0) return SeparateEnum;
173 else if (strcmp(name,"Sset")==0) return SsetEnum;
174 else if (strcmp(name,"Verbose")==0) return VerboseEnum;
175Index: ../trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp
176===================================================================
177--- ../trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 20648)
178+++ ../trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp (revision 20649)
179@@ -7,6 +7,8 @@
180 * case 4: 3d test with sinusoidal functions
181 * case 5: 3d test with sinusoidal functions, no variations with z, non homogeneous Dirichlet conditions
182 * case 6: 3d test with sinusoidal functions, non homogeneous Dirichlet conditions
183+ *
184+ * case 201: 3d test quadratic functions, BAFL flow law
185 */
186
187 #include <math.h>
188@@ -70,6 +72,9 @@
189
190 case 101:
191 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);
192+ case 201:
193+ 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);
194+
195 default:
196 _error_("FS analytical solution"<<testid<<" not implemented yet");
197 }
198@@ -130,6 +135,8 @@
199
200 case 101:
201 return 12*pow(x, 2)*y*z*(x - 1)*(y - 1)*(z - 1) + 6*pow(x, 2)*y*z*(x - 1)*(2*y - 1)*(z - 1) + 6*pow(x, 2)*z*(x - 1)*(y - 1)*(2*y - 1)*(z - 1) - 8*x*pow(y, 2)*z*(x - 1)*(2*x - 1)*(z - 1) - 4*x*pow(y, 2)*z*pow(y - 1, 2)*(z - 1) - 2*x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2) + 12*x*y*z*pow(x - 1, 2)*(y - 1)*(z - 1) + 6*x*y*z*pow(x - 1, 2)*(2*y - 1)*(z - 1) - 32*x*y*z*(x - 1)*(2*x - 1)*(y - 1)*(z - 1) + 6*x*z*pow(x - 1, 2)*(y - 1)*(2*y - 1)*(z - 1) - 8*x*z*(x - 1)*(2*x - 1)*pow(y - 1, 2)*(z - 1) - 4*pow(y, 2)*z*(x - 1)*pow(y - 1, 2)*(z - 1) - 2*pow(y, 2)*z*(2*x - 1)*pow(y - 1, 2)*(z - 1);
202+ case 201:
203+ return -x*pow(y, 2)*(x - 1)*(2*x - 1)*pow(y - 1, 2)/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) + x + (1.0L/18.0L)*(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*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)*(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)*(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)))/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*(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))/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*((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)))*(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*(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), 1.0L/3.0L);
204 default:
205 _error_("FS analytical solution"<<testid<<" not implemented yet");
206 }
207@@ -157,6 +164,9 @@
208 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);
209 case 24:
210 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);
211+ case 201:
212+ return (1.0L/18.0L)*(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)*(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)))/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/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))*((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) - (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((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) + (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))/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);
213+
214 default:
215 _error_("FS analytical solution"<<testid<<" not implemented yet");
216 }
217Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
218===================================================================
219--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 20648)
220+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 20649)
221@@ -83,6 +83,17 @@
222 iomodel->FetchDataToInput(elements,MaterialsRheologyEcEnum);
223 iomodel->FetchDataToInput(elements,MaterialsRheologyEsEnum);
224 for(i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matbafl(i+1,i,iomodel));
225+ switch(iomodel->domaindim){
226+ case 2:
227+ elements->InputDuplicate(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum);
228+ elements->InputDuplicate(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum);
229+ elements->InputDuplicate(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum);
230+ break;
231+ case 3:
232+ break;
233+ default:
234+ _error_("Mesh not supported yet");
235+ }
236 break;
237 default:
238 _error_("Materials "<<EnumToStringx(materials_type)<<" not supported");
239Index: ../trunk-jpl/src/c/classes/Materials/Matbafl.h
240===================================================================
241--- ../trunk-jpl/src/c/classes/Materials/Matbafl.h (revision 20648)
242+++ ../trunk-jpl/src/c/classes/Materials/Matbafl.h (revision 20649)
243@@ -81,6 +81,7 @@
244 void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
245 void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
246 /*}}}*/
247+ IssmDouble GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
248 };
249
250 #endif /* _MATBAFL_H_ */
251Index: ../trunk-jpl/src/c/classes/Materials/Matbafl.cpp
252===================================================================
253--- ../trunk-jpl/src/c/classes/Materials/Matbafl.cpp (revision 20648)
254+++ ../trunk-jpl/src/c/classes/Materials/Matbafl.cpp (revision 20649)
255@@ -274,27 +274,143 @@
256 /*}}}*/
257 void Matbafl::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
258
259- _assert_(dim==3);
260+ /*Intermediaries*/
261+ IssmDouble vx,vy,vz;
262+ IssmDouble dvx[3],dvy[3],dvz[3];
263+ IssmDouble ko,Ec,Es;
264
265+ /*Get velocity derivatives in all directions*/
266+ _assert_(dim>1);
267+ _assert_(vx_input);
268+ vx_input->GetInputValue(&vx,gauss);
269+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
270+ _assert_(vy_input);
271+ vy_input->GetInputValue(&vy,gauss);
272+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
273+ if(dim==3){
274+ _assert_(vz_input);
275+ vz_input->GetInputValue(&vz,gauss);
276+ vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
277+ }
278+ else{
279+ vz = 0.;
280+ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
281+ }
282+
283+ /*Get enhancement factors and ko*/
284+ Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
285+ Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
286+ Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
287+ ec_input->GetInputValue(&Ec,gauss);
288+ es_input->GetInputValue(&Es,gauss);
289+ ko_input->GetInputValue(&ko,gauss);
290+
291+ /*Compute viscosity*/
292+ *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
293+}
294+/*}}}*/
295+void Matbafl::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
296+ this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
297+}/*}}}*/
298+void Matbafl::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
299+
300 /*Intermediaries*/
301- IssmDouble viscosity;
302 IssmDouble vx,vy,vz;
303 IssmDouble dvx[3],dvy[3],dvz[3];
304+ IssmDouble ko,Ec,Es;
305+
306+ /*Get velocity derivatives in all directions*/
307+ _assert_(dim==2 || dim==3);
308+ _assert_(vx_input);
309+ vx_input->GetInputValue(&vx,gauss);
310+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
311+ if(dim==3){
312+ _assert_(vy_input);
313+ vy_input->GetInputValue(&vy,gauss);
314+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
315+ }
316+ else{
317+ vy = 0.;
318+ dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
319+ }
320+ vz = 0.;
321+ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
322+
323+ /*Get enhancement factors and ko*/
324+ Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
325+ Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
326+ Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
327+ ec_input->GetInputValue(&Ec,gauss);
328+ es_input->GetInputValue(&Es,gauss);
329+ ko_input->GetInputValue(&ko,gauss);
330+
331+ /*Compute viscosity*/
332+ *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
333+}/*}}}*/
334+void Matbafl::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
335+ _error_("not implemented yet");
336+}/*}}}*/
337+void Matbafl::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
338+ _error_("not implemented yet");
339+}/*}}}*/
340+void Matbafl::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
341+ /*Intermediaries*/
342+ IssmDouble vx,vy,vz;
343+ IssmDouble dvx[3],dvy[3],dvz[3];
344+ IssmDouble ko,Ec,Es;
345+
346+ /*Get velocity derivatives in all directions*/
347+ _assert_(dim==1 || dim==2);
348+ _assert_(vx_input);
349+ vx_input->GetInputValue(&vx,gauss);
350+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
351+ if(dim==2){
352+ _assert_(vy_input);
353+ vy_input->GetInputValue(&vy,gauss);
354+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
355+ }
356+ else{
357+ vy = 0.;
358+ dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
359+ }
360+ vz = 0.;
361+ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
362+
363+ /*Get enhancement factors and ko*/
364+ Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input);
365+ Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input);
366+ Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
367+ ec_input->GetInputValue(&Ec,gauss);
368+ es_input->GetInputValue(&Es,gauss);
369+ ko_input->GetInputValue(&ko,gauss);
370+
371+ /*Compute viscosity*/
372+ *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
373+}/*}}}*/
374+void Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
375+ _error_("not implemented yet");
376+}/*}}}*/
377+void Matbafl::ResetHooks(){/*{{{*/
378+
379+ this->element=NULL;
380+
381+ /*Get Element type*/
382+ this->helement->reset();
383+
384+}
385+/*}}}*/
386+IssmDouble Matbafl::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
387+
388+ /*Intermediaries*/
389+ IssmDouble viscosity;
390 IssmDouble vorticity[3],vorticity_norm;
391 IssmDouble nrsp[3],nrsp_norm;
392 IssmDouble eps[3][3],epso;
393 IssmDouble epsprime[3],epsprime_norm;
394- IssmDouble ko,E,Ec,Es,lambdas;
395+ IssmDouble E,lambdas;
396
397- /*Get velocity derivatives in all directions*/
398- vx_input->GetInputValue(&vx,gauss);
399- vy_input->GetInputValue(&vy,gauss);
400- vz_input->GetInputValue(&vz,gauss);
401- vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
402- vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
403- vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
404-
405 /*Create vorticity vector*/
406+ _assert_(dvx && dvy && dvz);
407 vorticity[0] = dvz[1] - dvy[2];
408 vorticity[1] = dvx[2] - dvz[0];
409 vorticity[2] = dvy[0] - dvx[1];
410@@ -376,16 +492,8 @@
411 lambdas = sqrt(2*epsprime_norm/(3*epso));
412
413 /*Get total enhancement factor E(lambdas)*/
414- Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
415- Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
416- ec_input->GetInputValue(&Ec,gauss);
417- es_input->GetInputValue(&Es,gauss);
418 E = Ec + (Es-Ec)*lambdas*lambdas;
419
420- /*Get flow law parameter*/
421- Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
422- ko_input->GetInputValue(&ko,gauss);
423-
424 /*Compute viscosity*/
425 _assert_(E>0.);
426 _assert_(ko>0.);
427@@ -393,34 +501,6 @@
428 viscosity = 1./(2*pow(ko*E*epso*epso,1./3.));
429
430 /*Assign output pointer*/
431- *pviscosity=viscosity;
432+ return viscosity;
433 }
434 /*}}}*/
435-void Matbafl::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
436- this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
437-}/*}}}*/
438-void Matbafl::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
439-
440- _error_("not implemented yet");
441-}/*}}}*/
442-void Matbafl::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
443- _error_("not implemented yet");
444-}/*}}}*/
445-void Matbafl::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
446- _error_("not implemented yet");
447-}/*}}}*/
448-void Matbafl::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
449- _error_("not implemented yet");
450-}/*}}}*/
451-void Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
452- _error_("not implemented yet");
453-}/*}}}*/
454-void Matbafl::ResetHooks(){/*{{{*/
455-
456- this->element=NULL;
457-
458- /*Get Element type*/
459- this->helement->reset();
460-
461-}
462-/*}}}*/
463Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
464===================================================================
465--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20648)
466+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20649)
467@@ -1410,7 +1410,7 @@
468 void Element::InputDuplicate(int original_enum,int new_enum){/*{{{*/
469
470 /*Call inputs method*/
471- if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
472+ if(IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
473
474 }
475 /*}}}*/
476@@ -1530,6 +1530,12 @@
477 name==MaterialsRheologyBEnum ||
478 name==MaterialsRheologyBbarEnum ||
479 name==MaterialsRheologyNEnum ||
480+ name==MaterialsRheologyKoEnum ||
481+ name==MaterialsRheologyKobarEnum ||
482+ name==MaterialsRheologyEcEnum ||
483+ name==MaterialsRheologyEcbarEnum ||
484+ name==MaterialsRheologyEsEnum ||
485+ name==MaterialsRheologyEsbarEnum ||
486 name==SealevelEnum ||
487 name==SealevelEustaticEnum ||
488 name==SealevelriseDeltathicknessEnum ||
489Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
490===================================================================
491--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20648)
492+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20649)
493@@ -2250,6 +2250,9 @@
494 if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
495 break;
496 case MatbaflEnum:
497+ this->InputDepthAverageAtBase(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum);
498+ this->InputDepthAverageAtBase(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum);
499+ this->InputDepthAverageAtBase(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum);
500 break;
501 default:
502 _error_("not supported yet");
Note: See TracBrowser for help on using the repository browser.