Changeset 16743

Timestamp:

11/13/13 15:44:10 (11 years ago)

Author:

seroussi

Message:

CHG: better B in enthalpy solution

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (6 diffs)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶

@@ -4851 +4851 @@
         IssmDouble  values[numdof];
         IssmDouble  B[numdof],surface[numdof];
-        IssmDouble  B_average;
         int        *doflist = NULL;
         bool        hack    = false;
@@ -4890 +4889 @@
                  * otherwise the rheology could be negative*/
                 this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+                GetInputListOnVertices(&surface[0],SurfaceEnum);
                 switch(rheology_law){
                         case NoneEnum:
@@ -4895 +4895 @@
                                 break;
                         case PatersonEnum:
-                                B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0);
-                                for(i=0;i<numdof;i++) B[i]=B_average;
+                                for(i=0;i<numdof;i++) B[i]=Paterson(values[i]);
                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
                                 break;
                         case ArrheniusEnum:
-                                GetInputListOnVertices(&surface[0],SurfaceEnum);
                                 for(i=0;i<numdof;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i][2],material->GetN());
                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
@@ -4926 +4924 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble values[numdof];
-        IssmDouble pressure[NUMVERTICES];
+        IssmDouble pressure[numdof];
+        IssmDouble surface[numdof];
         IssmDouble temperatures[numdof];
         IssmDouble waterfraction[numdof];
@@ -4965 +4964 @@
                  * otherwise the rheology could be negative*/
                 this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+                GetInputListOnVertices(&surface[0],PressureEnum);
                 switch(rheology_law){
                         case NoneEnum:
@@ -4970 +4970 @@
                                 break;
                         case PatersonEnum:
-                                B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0);
-                                for(i=0;i<numdof;i++) B[i]=B_average;
+                                for(i=0;i<numdof;i++) B[i]=Paterson(values[i]);
                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
                                 break;
                         case ArrheniusEnum:
-                                surface_input->GetInputAverage(&s_average);
-                                B_average=Arrhenius((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0,
-                                                        s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
-                                                        material->GetN());
-                                for(i=0;i<numdof;i++) B[i]=B_average;
+                                for(i=0;i<numdof;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i][2],material->GetN());
                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
                                 break;
                         case LliboutryDuvalEnum:
-                                B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,
-                                                                                                                                 (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0,
-                                                                                                                                 material->GetN(), 
-                                                                                                                                 matpar->GetBeta(), matpar->GetReferenceTemperature(), matpar->GetHeatCapacity(), matpar->GetLatentHeat());
-                                for(i=0;i<numdof;i++) B[i]=B_average;
+                                for(i=0;i<numdof;i++) B[i]=LliboutryDuval(values[i],pressure[i],material->GetN(),matpar->GetBeta(),matpar->GetReferenceTemperature(),matpar->GetHeatCapacity(),matpar->GetLatentHeat());
                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
                                 break;