Changeset 1361

Timestamp:

07/23/09 10:03:59 (16 years ago)

Author:

Mathieu Morlighem

Message:

fixed modelextract2

Location:

issm/trunk/src/m/classes/public

Files:

• ismodelselfconsistent.m (modified) (1 diff)
• modelextract2.m (modified) (5 diffs)

issm/trunk/src/m/classes/public/ismodelselfconsistent.m

@@ -328 +328 @@
                 end
         end
+
+        %SINGULAR
+        if ~any(md.gridondirichlet_diag),
+                disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
+                bool=0;return;
+        end
+
+        %DIRICHLET VALUES
+        if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
+                disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
+                bool=0;return;
+        end
+
+        %DIRICHLET IF THICKNESS <= 0
+        if any(md.thickness<=0),
+                pos=find(md.thickness<=0);
+                if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
+                        disp(['model ' md.name ' has some grids with 0 thickness']);
+                        bool=0; return;
+                end
+        end
 end
 

issm/trunk/src/m/classes/public/modelextract2.m

@@ -70 +70 @@
 end
 
-%element and grid position
+%kick out all elements with 3 dirichlets
+spc_elem=find(~flag_elem);
+spc_grid=sort(unique(md1.elements(spc_elem,:)));
+flag=ones(md1.numberofgrids,1);
+flag(spc_grid)=0;
+pos=find(sum(flag(md1.elements),2)==0);
+flag_elem(pos)=0;
+
+%extracted elements and nodes lists
 pos_elem=find(flag_elem);
 pos_grid=sort(unique(md1.elements(pos_elem,:)));
@@ -103 +111 @@
         md2=md1;
 
+        %automatically modify fields
+
+        %loop over model fields
+        model_fields=fields(md1);
+        for i=1:length(model_fields),
+
+                %get field
+                field=md1.(model_fields(i));
+                fieldsize=size(field);
+
+                %size = number of grids * n
+                if fieldsize(1)==numberofgrids1
+                        md2.(model_fields(i))=field(pos_grid,:);
+
+                %size = number of elements * n
+                elseif fieldsize(1)==numberofelements1
+                        md2.(model_fields(i))=field(pos_elem,:);
+                end
+        end
+
         %modify some specific fields
 
@@ -144 +172 @@
         if ~isnan(md2.elements_type), md2.elements_type=md1.elements_type(pos_elem,:);end
 
-        %Dirichlets Diag
+        %Boundary conditions: Dirichlets on new boundary
         %Catch the elements that have not been extracted
         orphans_elem=find(~flag_elem);
@@ -151 +179 @@
         gridstoflag1=intersect(orphans_grid,pos_grid);
         gridstoflag2=Pgrid(gridstoflag1);
-        if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
-                md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); 
-                md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
-        else
-                md2.dirichletvalues_diag=zeros(numberofgrids2,2);
-                disp(' ')
-                disp('!! modelextract warning: dirichlet values should be checked !!')
-                disp(' ')
+        if ~isnan(md1.gridondirichlet_diag),
+                md2.gridondirichlet_diag(gridstoflag2)=1;
+                if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
+                        md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); 
+                        md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
+                else
+                        md2.dirichletvalues_diag=zeros(numberofgrids2,2);
+                        disp(' ')
+                        disp('!! modelextract warning: dirichlet values should be checked !!')
+                        disp(' ')
+                end
+        end
+        if ~isnan(md1.gridondirichlet_thermal),
+                md2.gridondirichlet_thermal(gridstoflag2)=1;
+                if ~isnan(md1.observed_temperature)
+                        md2.dirichletvalues_thermal(gridstoflag2,1)=md2.observed_temperature(gridstoflag2); 
+                else
+                        md2.dirichletvalues_thermal=zeros(numberofgrids2,2);
+                        disp(' ')
+                        disp('!! modelextract warning: dirichlet values should be checked !!')
+                        disp(' ')
+                end
         end
 
@@ -215 +257 @@
         md2.strainrate=NaN;
 
-        %other fields
-
-        %loop over model fields
-        model_fields=fields(md1);
-        for i=1:length(model_fields),
-
-                %get field
-                field=md1.(model_fields(i));
-
-                %size = number of grids * 1
-                if (size(field,1)==numberofgrids1 & size(field,2)==1),
-                        md2.(model_fields(i))=field(pos_grid);
-
-                %size = number of grids * 2
-                elseif (size(field,1)==numberofgrids1 & size(field,2)==2), 
-                        md2.(model_fields(i))=[field(pos_grid,1) field(pos_grid,2)];
-
-                %size = number of elements * 1
-                elseif (size(field,1)==numberofelements1 & size(field,2)==1),
-                        md2.(model_fields(i))=field(pos_elem);
-
-                end
-        end
-
 %Keep track of pos_grid and pos_elem
 md2.extractedgrids=pos_grid;