Changeset 21542

Timestamp:

02/10/17 10:32:55 (8 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added Higher order finite elements for enthalpy

Location:

issm/trunk-jpl/src

Files:

• c/analyses/EnthalpyAnalysis.cpp (modified) (3 diffs)
• c/analyses/ThermalAnalysis.cpp (modified) (3 diffs)
• c/classes/IoModel.cpp (modified) (2 diffs)
• m/classes/thermal.m (modified) (5 diffs)
• m/classes/thermal.py (modified) (3 diffs)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -98 +98 @@
 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");
 }/*}}}*/
@@ -121 +124 @@
         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;
@@ -126 +133 @@
                 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++;
                 }

issm/trunk-jpl/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");
@@ -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++){

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -656 +656 @@
                                                 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){
@@ -797 +799 @@
                                                 } 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));

issm/trunk-jpl/src/m/classes/thermal.m

@@ -15 +15 @@
                 isenthalpy        = 0;
                 isdynamicbasalspc = 0;
+                fe                = 'P1';
                 requested_outputs = {};
         end
@@ -66 +67 @@
                         self.isdynamicbasalspc=0;
 
+                        %Linear elements by default
+                        self.fe='P1';
+
                         %default output
                         self.requested_outputs={'default'};
@@ -76 +80 @@
                         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),
 
@@ -110 +115 @@
                         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');
 
@@ -122 +128 @@
                         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');
 

issm/trunk-jpl/src/m/classes/thermal.py

@@ -23 +23 @@
                 self.isenthalpy        = 0
                 self.isdynamicbasalspc = 0;
+                self.fe                = 'P1';
                 self.requested_outputs = []
 
@@ -79 +80 @@
 
                 #will basal boundary conditions be set dynamically
-                self.isdynamicbasalspc=0;
+                self.isdynamicbasalspc=0
+
+                #Finite element interpolation
+                self.fe='P1'
 
                 #default output
@@ -124 +128 @@
                 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');