Changeset 9075

Timestamp:

07/19/11 16:28:47 (14 years ago)

Author:

Mathieu Morlighem

Message:

Simplified spcthickness (use new NaN scheme)

Location:

issm/trunk/src

Files:

• c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (modified) (4 diffs)
• c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (modified) (2 diffs)
• c/objects/Elements/Tria.cpp (modified) (1 diff)
• m/classes/model.m (modified) (2 diffs)
• m/model/display/displaybc.m (modified) (1 diff)
• m/model/display/displayprognostic.m (modified) (1 diff)
• m/model/ismodelselfconsistent.m (modified) (4 diffs)
• m/model/outflow.m (modified) (1 diff)
• m/solvers/solver_linear.m (modified) (1 diff)
• m/utils/BC/SetIceSheetBC.m (modified) (1 diff)
• m/utils/BC/SetIceShelfBC.m (modified) (1 diff)
• m/utils/BC/SetMarineIceSheetBC.m (modified) (1 diff)
• m/utils/Interp/PatchToVec.m (modified) (1 diff)

issm/trunk/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

@@ -42 +42 @@
 
         /*Transient or static?:*/
-        if(M==iomodel->numberofvertices){ //{{{1
+        if(M==iomodel->numberofvertices){
                 /*static: just create Spc objects*/
-
-        
-                /*Initialize counter*/
                 count=0;
         
@@ -60 +57 @@
                                         count++;
                                 }
-
-                        } //if((my_nodes[i]==1))
+                        }
                 }
-        } //}}}
-        else{ //{{{1
+        }
+        else{
                 /*transient: create transient Spct objects. Same logic, except we need to retrieve 
                  * various times and values to initialize an Spct object: */
-
-                /*Initialize counter*/
                 count=0;
 
@@ -78 +72 @@
                 /*unit conversion: */
                 UnitConversion(times,N,ExtToIuEnum,TimeEnum);
-
-
 
                 /*Create spcs from x,y,z, as well as the spc values on those spcs: */
@@ -100 +92 @@
                                 }
                                 xfree((void**)&values);
-                        } //if((my_nodes[i]==1))
+                        }
                 }
-        } //}}}
+        }
 
         /*Free ressources:*/

issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp

@@ -1 @@
-/*
- * CreateConstraintsPrognostic.c:
- */
-
 #include "../../../Container/Container.h"
+#include "../../../modules/modules.h"
 #include "../../../toolkits/toolkits.h"
 #include "../../../EnumDefinitions/EnumDefinitions.h"
@@ -27 +24 @@
         if(!constraints) constraints = new Constraints(ConstraintsEnum);
 
-        /*Do not add constraints in DG*/
+        /*Do not add constraints in DG, they are weakly imposed*/
         if(!iomodel->prognostic_DG){
-
-                /*Fetch data: */
-                IoModelFetchData(&iomodel->spcthickness,NULL,NULL,iomodel_handle,SpcthicknessEnum);
-
-                /*Initialize counter*/
-                count=0;
-
-                /*Create spcs from x,y,z, as well as the spc values on those spcs: */
-                for (i=0;i<iomodel->numberofvertices;i++){
-                        /*keep only this partition's nodes:*/
-                        if((iomodel->my_vertices[i])){
-
-                                if ((int)iomodel->spcthickness[2*i]){
-
-                                        constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,*(iomodel->spcthickness+2*i+1),PrognosticAnalysisEnum));
-                                        count++;
-                                }
-                        } //if((my_vertices[i]))
-                }
-
-                /*Free data: */
-                xfree((void**)&iomodel->spcthickness);
+                IoModelToConstraintsx(constraints,iomodel,iomodel_handle,SpcthicknessEnum,PrognosticAnalysisEnum);
         }
 

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -775 +775 @@
         /*Constants*/
         const int  numdof=NDOF2*NUMVERTICES;
+        _error_("STOP");
 
         /*Intermediaries*/

issm/trunk/src/m/classes/model.m

@@ -15 +15 @@
                  name            = {'',true,'String'};
                  runtimename     = {'',false}; %name used when running a parallel job
-                 counter         = {0,false}; %1 mesh, 2 geography, 3 parameterize
+                 counter         = {0,false};
                  domainoutline   = {NaN,false};
                  riftoutline     = {NaN,false};
@@ -346 +346 @@
          end
          methods (Static)
-                 function md = loadobj(md) % {{{
+                 function md = lllloadobj(md) % {{{
                          % This function is directly called by matlab when a model object is
                          % loaded. If the input is a struct it is an old version of model and

issm/trunk/src/m/model/display/displaybc.m

@@ -25 +25 @@
 
 disp(sprintf('\n      prognostic:'));
-fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)');
+fielddisplay(md,'spcthickness','thickness constraints (NaN means no constraint)');
 
 disp(sprintf('\n      thermal:'));
-fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)');
+fielddisplay(md,'spctemperature','temperature constraints (NaN means no constraint)');
 fielddisplay(md,'basal_melting_rate','basal melting rate [m/a]');

issm/trunk/src/m/model/display/displayprognostic.m

@@ -18 +18 @@
 
 disp(sprintf('\n      boundary conditions:'));
-fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)');
+fielddisplay(md,'spcthickness','thickness constraints (NaN means no constraint)');

issm/trunk/src/m/model/ismodelselfconsistent.m

@@ -224 +224 @@
         if any(md.thickness<=0),
                 pos=find(md.thickness<=0);
-                if any(find(md.spcthickness(pos,1)==0)),
+                if any(isnan(md.spcthickness(pos))),
                         message(['model not consistent: model ' md.name ' has some nodes with 0 thickness']);
                 end
@@ -369 +369 @@
                         if any(md.thickness<=0),
                                 pos=find(md.thickness<=0);
-                                if any(find(md.spcthickness(pos,1)==0)),
+                                if any(isnan(md.spcthickness(pos))),
                                         message(['model not consistent: model ' md.name ' has some nodes with 0 thickness']);
                                 end
@@ -458 +458 @@
                         checksize(md,fields,[md.numberofnodes 1]);
                         checknan(md,fields);
-
-                        %CHECK THAT WE ARE NOT FULLY CONSTRAINED
-                        if (md.dim==2),
-                                if isnan(find(~md.spcthickness(:,1))),
-                                        message(['model not consistent: model ' md.name ' is totally constrained for prognostic, no need to solve!']);
-                                end
-                        end
                         %}}}
                 case HydrologyAnalysisEnum,
@@ -560 +553 @@
                         checksize(md,fields,[md.numberofnodes 1]);
                         checknan(md,fields);
-
-                        %SPC                             
-                        %       if ~md.prognostic_DG,
-                        %               if any(md.spcthickness(find(md.nodeonboundary))~=1),             
-                        %                       message(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);                    
-                        %               end 
-                        %       end
 
                         %Triangle with zero velocity

issm/trunk/src/m/model/outflow.m

@@ -1 +1 @@
 function flag=outflow(md);
+%OUTFLOW - flag nodes on outflux boundary
+%
+%   Usage:
+%      flag=outflow(md);
 
 A=md.segments(:,1);

issm/trunk/src/m/solvers/solver_linear.m

@@ -9 +9 @@
 
         [K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+        full(K_ff)
+        condest(K_ff)
         ys  =CreateNodalConstraints(femmodel.nodes,configuration_type);
         p_f =Reduceload( p_f, K_fs, ys);
+        full(p_f)
 
         issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
         u_f=Solver(K_ff,p_f,[],df,femmodel.parameters);
+        u_f
+        error('stop')
         u_g= Mergesolutionfromftog( u_f, ys, femmodel.nodes,femmodel.parameters); 
         [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,u_g);

issm/trunk/src/m/utils/BC/SetIceSheetBC.m

@@ -58 +58 @@
 end
 
-md.spcthickness=zeros(md.numberofnodes,2);
+md.spcthickness=NaN*ones(md.numberofnodes,1);
 
 if (length(md.observed_temperature)==md.numberofnodes),

issm/trunk/src/m/utils/BC/SetIceShelfBC.m

@@ -90 +90 @@
 end
 
-md.spcthickness=zeros(md.numberofnodes,2);
+md.spcthickness=NaN*ones(md.numberofnodes,1);
 
 if (length(md.observed_temperature)==md.numberofnodes),

issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m

@@ -101 +101 @@
 end
 
-md.spcthickness=zeros(md.numberofnodes,2);
+md.spcthickness=NaN*ones(md.numberofnodes,1);
 
 if (length(md.observed_temperature)==md.numberofnodes),

issm/trunk/src/m/utils/Interp/PatchToVec.m

@@ -8 +8 @@
 switch(Patch.interpolation(1)),
         case P0Enum,
-                vec(Patch.element)=Patch.value;
+                %vec(Patch.element)=Patch.value;
+                connectivity=sparse(Patch.index(:),1,1);
+                value       =sparse(Patch.index(:),1,Patch.value(:));
+                vec=full(value./connectivity);
         case P1Enum,
                 connectivity=sparse(Patch.index(:),1,1);