Changeset 21546

Timestamp:

02/10/17 11:28:06 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: use IoModelToConstrain to create spcs

File:

• issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (modified) (3 diffs)

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

@@ -13 +13 @@
         int        M,N;
         bool       spcpresent = false;
+        int        finiteelement;
         IssmDouble heatcapacity;
         IssmDouble referencetemperature;
@@ -24 +25 @@
         iomodel->FindConstant(&heatcapacity,"md.materials.heatcapacity");
         iomodel->FindConstant(&referencetemperature,"md.constants.referencetemperature");
+        iomodel->FindConstant(&finiteelement,"md.thermal.fe");
 
         /*return if 2d mesh*/
@@ -31 +33 @@
         iomodel->FetchData(&spcvector,&M,&N,"md.thermal.spctemperature");
 
-        //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS 
-        /*Transient or static?:*/
-        if(M==iomodel->numberofvertices){
-                /*static: just create Constraints objects*/
-                count=0;
-
-                for(int i=0;i<iomodel->numberofvertices;i++){
-                        /*keep only this partition's nodes:*/
-                        if((iomodel->my_vertices[i])){
-
-                                if (!xIsNan<IssmDouble>(spcvector[i])){
-
-                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
-                                        count++;
-
-                                }
-                        }
-                }
-        }
-        else if (M==(iomodel->numberofvertices+1)){
-                /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve 
-                 * various times and values to initialize an SpcTransient object: */
-                count=0;
-
-                /*figure out times: */
-                times=xNew<IssmDouble>(N);
-                for(int j=0;j<N;j++){
-                        times[j]=spcvector[(M-1)*N+j];
-                }
-
-                /*Create constraints from x,y,z: */
-                for(int i=0;i<iomodel->numberofvertices;i++){
-
-                        /*keep only this partition's nodes:*/
-                        if((iomodel->my_vertices[i])){
-
-                                /*figure out times and values: */
-                                values=xNew<IssmDouble>(N);
-                                spcpresent=false;
-                                for(int j=0;j<N;j++){
-                                        values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
-                                        if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
-                                }
-
-                                if(spcpresent){
-                                        constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,N,times,values,EnthalpyAnalysisEnum));
-                                        count++;
-                                }
-                                xDelete<IssmDouble>(values);
-                        }
-                }
-        }
-        else{
-                _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
-        }
+        /*Convert spcs from temperatures to enthalpy*/
+        _assert_(N>0); _assert_(M>=iomodel->numberofvertices);
+        for(int i=0;i<iomodel->numberofvertices;i++){
+                for(int j=0;i<N;j++){
+                        spcvector[i*N+j] = heatcapacity*(spcvector[i]-referencetemperature);
+                }
+        }
+        IoModelToConstraintsx(constraints,iomodel,spcvector,M,N,EnthalpyAnalysisEnum,finiteelement,2);
 
         /*Free ressources:*/