Changeset 21759 for issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/ThermalAnalysis.cpp

Timestamp:

06/07/17 10:50:54 (8 years ago)

Author:

Eric.Larour

Message:

CHG: merged branch back to trunk-jpl 21754.

File:

• issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/ThermalAnalysis.cpp (modified) (14 diffs)

issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/ThermalAnalysis.cpp

@@ -9 +9 @@
 
         /*Intermediary*/
-        int finiteelement = P1Enum;
+        int finiteelement;
+        iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+        _assert_(finiteelement==P1Enum); 
 
         /*Only 3d mesh supported*/
@@ -62 +64 @@
 void ThermalAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
-        int finiteelement = P1Enum;
+        int finiteelement;
+        iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+        _assert_(finiteelement==P1Enum); 
         
         if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
@@ -73 +77 @@
 void ThermalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
-        int frictionlaw,basalforcing_model;
+        int frictionlaw,basalforcing_model,materialstype;
         int FrictionCoupling;
         /*Now, is the model 3d? otherwise, do nothing: */
@@ -79 +83 @@
 
         /*Update elements: */
-        int finiteelement = P1Enum;
+        int finiteelement;
+        iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+        _assert_(finiteelement==P1Enum); 
         int counter=0;
         for(int i=0;i<iomodel->numberofelements;i++){
@@ -93 +99 @@
         iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
         iomodel->FindConstant(&frictionlaw,"md.friction.law");
+        iomodel->FindConstant(&materialstype,"md.materials.type");
 
         iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
@@ -106 +113 @@
         iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
         iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
-        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
-        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
@@ -116 +121 @@
         InputUpdateFromConstantx(elements,0.,VyMeshEnum);
         InputUpdateFromConstantx(elements,0.,VzMeshEnum);
+
+        /*Rheology type*/
+        iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
+        switch(materialstype){
+                case MatenhancediceEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
+                        break;
+                case MatdamageiceEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+                        break;
+                case MatestarEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
+                        break;
+                case MaticeEnum:
+                        iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+                        break;
+                default:
+                        _error_("not supported");
+        }
         if(ismovingfront){
                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); // required for updating active nodes
@@ -131 +157 @@
         switch(frictionlaw){
                 case 1:
+                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
                         iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
                         iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
                         iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
+                        if (FrictionCoupling==1){
+                          iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
+                        }
                         break;
                 case 2:
@@ -144 +174 @@
                         iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
                         iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
-                        if (FrictionCoupling==0){
-                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
+                        if (FrictionCoupling==1){
+                          iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
                         }
                         break;
@@ -194 +224 @@
         iomodel->FindConstant(&frictionlaw,"md.friction.law");
         if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
-        if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
+        if(frictionlaw==3 || frictionlaw==1) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
+
 }/*}}}*/
 
@@ -722 +753 @@
         int        *doflist   = NULL;
         IssmDouble *xyz_list  = NULL;
+        IssmDouble  n=3.0;
         bool        hack      = false;
 
@@ -756 +788 @@
 
         /*Get all inputs and parameters*/
+        if(element->material->ObjectEnum()!=MatestarEnum) n=element->GetMaterialParameter(MaterialsRheologyNEnum);
         element->GetInputValue(&converged,ConvergedEnum);
         if(converged){
@@ -767 +800 @@
                         case NoneEnum:
                                 /*Do nothing: B is not temperature dependent*/
+                                break;
+                        case BuddJackaEnum:
+                                for(i=0;i<numnodes;i++) B[i]=BuddJacka(values[i]);
+                                element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
                                 break;
                         case CuffeyEnum:
@@ -778 +815 @@
                         case ArrheniusEnum:{
                                 element->GetVerticesCoordinates(&xyz_list);
-                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+                                for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],n);
                                 element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
                                 break;