source: issm/oecreview/Archive/21337-21723/ISSM-21541-21542.diff@ 21726

Index: ../trunk-jpl/src/m/classes/thermal.py
===================================================================
--- ../trunk-jpl/src/m/classes/thermal.py       (revision 21541)
+++ ../trunk-jpl/src/m/classes/thermal.py       (revision 21542)
@@ -22,6 +22,7 @@
                self.penalty_factor    = 0
                self.isenthalpy        = 0
                self.isdynamicbasalspc = 0;
+               self.fe                = 'P1';
                self.requested_outputs = []
 
                #set defaults
@@ -78,8 +79,11 @@
                self.isenthalpy=0
 
                #will basal boundary conditions be set dynamically
-               self.isdynamicbasalspc=0;
+               self.isdynamicbasalspc=0
 
+               #Finite element interpolation
+               self.fe='P1'
+
                #default output
                self.requested_outputs=['default']
                return self
@@ -123,6 +127,7 @@
                WriteData(fid,prefix,'object',self,'fieldname','penalty_lock','format','Integer')
                WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double')
                WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean')
+               WriteData(fid,prefix,'object',self,'fieldname','fe','format','String');
                WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean');
 
                #process requested outputs
Index: ../trunk-jpl/src/m/classes/thermal.m
===================================================================
--- ../trunk-jpl/src/m/classes/thermal.m        (revision 21541)
+++ ../trunk-jpl/src/m/classes/thermal.m        (revision 21542)
@@ -14,6 +14,7 @@
                penalty_factor    = 0;
                isenthalpy        = 0;
                isdynamicbasalspc = 0;
+               fe                = 'P1';
                requested_outputs = {};
        end
        methods
@@ -65,6 +66,9 @@
                        %will basal boundary conditions be set dynamically
                        self.isdynamicbasalspc=0;
 
+                       %Linear elements by default
+                       self.fe='P1';
+
                        %default output
                        self.requested_outputs={'default'};
                end % }}}
@@ -75,6 +79,7 @@
 
                        md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2]);
                        md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1,'>=',0);
+                       md = checkfield(md,'fieldname','thermal.fe','values',{'P1','P1xP2'});
                        if (ismember('EnthalpyAnalysis',analyses) & md.thermal.isenthalpy & dimension(md.mesh)==3),
 
                                %Make sure the spc are less than melting point
@@ -109,6 +114,7 @@
                        fielddisplay(self,'penalty_factor','scaling exponent (default is 3)');
                        fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
                        fielddisplay(self,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);
+                       fielddisplay(self,'fe','Finite Element type: ''P1'' (default), ''P1xP2''');
                        fielddisplay(self,'requested_outputs','additional outputs requested');
 
                end % }}}
@@ -121,6 +127,7 @@
                        WriteData(fid,prefix,'object',self,'fieldname','penalty_lock','format','Integer');
                        WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double');
                        WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean');
+                       WriteData(fid,prefix,'object',self,'fieldname','fe','format','String');
                        WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean');
 
                        %process requested outputs
Index: ../trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/IoModel.cpp      (revision 21541)
+++ ../trunk-jpl/src/c/classes/IoModel.cpp      (revision 21542)
@@ -656,6 +656,8 @@
                                                this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
                                        } else if(strcmp(record_name,"md.flowequation.fe_FS")==0){
                                                this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
+                                       } else if(strcmp(record_name,"md.thermal.fe")==0){
+                                               this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
                                        } else if(strcmp(record_name,"md.groundingline.migration")==0){
                                                this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
                                        } else if(strcmp(record_name,"md.masstransport.hydrostatic_adjustment")==0){
@@ -796,6 +798,8 @@
                                                        this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
                                                } else if(strcmp(record_name,"md.flowequation.fe_FS")==0){
                                                        this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
+                                               } else if(strcmp(record_name,"md.thermal.fe")==0){
+                                                       this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
                                                } else if(strcmp(record_name,"md.groundingline.migration")==0){
                                                        this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
                                                } else if(strcmp(record_name,"md.masstransport.hydrostatic_adjustment")==0){
Index: ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp    (revision 21541)
+++ ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp    (revision 21542)
@@ -97,8 +97,11 @@
 }/*}}}*/
 void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
+       int finiteelement;
+       iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+
        if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
-       ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
+       ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,finiteelement);
        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
 }/*}}}*/
 int  EnthalpyAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
@@ -120,12 +123,16 @@
        /*Fetch data needed: */
        iomodel->FetchData(3,"md.initialization.temperature","md.initialization.waterfraction","md.initialization.pressure");
 
+       /*Finite element type*/
+       int finiteelement;
+       iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+
        /*Update elements: */
        int counter=0;
        for(int i=0;i<iomodel->numberofelements;i++){
                if(iomodel->my_elements[i]){
                        Element* element=(Element*)elements->GetObjectByOffset(counter);
-                       element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
+                       element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
                        counter++;
                }
        }
Index: ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp     (revision 21541)
+++ ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp     (revision 21542)
@@ -8,7 +8,9 @@
 void ThermalAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
        /*Intermediary*/
-       int finiteelement = P1Enum;
+       int finiteelement;
+       iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+       _assert_(finiteelement==P1Enum); 
 
        /*Only 3d mesh supported*/
        if(iomodel->domaintype!=Domain3DEnum) _error_("not supported yet");
@@ -61,7 +63,9 @@
 }/*}}}*/
 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");
        ::CreateNodes(nodes,iomodel,ThermalAnalysisEnum,finiteelement);
@@ -78,7 +82,9 @@
        if(iomodel->domaintype==Domain2DhorizontalEnum)return;
 
        /*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++){
                if(iomodel->my_elements[i]){