Changeset 21958

Timestamp:

08/15/17 14:17:24 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: enable plot for amr

Location:

issm/trunk-jpl/src/m/plot

Files:

• processdata.m (modified) (10 diffs)
• processmesh.m (modified) (2 diffs)

issm/trunk-jpl/src/m/plot/processdata.m

@@ -47 +47 @@
         numberofelements2d=NaN;
 end
+numberofvertices = md.mesh.numberofvertices;
+numberofelements = md.mesh.numberofelements;
+if exist(options,'amr'),
+        step = getfieldvalue(options,'amr');
+        numberofvertices = numel(md.results.TransientSolution(step).MeshX);
+        numberofelements = numel(md.results.TransientSolution(step).MeshElements);
+end
 
 %initialize datatype
@@ -67 +74 @@
 
 %check length
-if datasize(1)~=md.mesh.numberofvertices & datasize(1)~=md.mesh.numberofelements & datasize(1)~=md.mesh.numberofvertices*6 & (strcmp(md.mesh.domaintype(),'3D') & ~(datasize(1)==numberofelements2d | datasize(1)==numberofvertices2d))
+if datasize(1)~=numberofvertices & datasize(1)~=numberofelements & datasize(1)~=numberofvertices*6 & (strcmp(md.mesh.domaintype(),'3D') & ~(datasize(1)==numberofelements2d | datasize(1)==numberofvertices2d))
         error('plotmodel error message: data not supported yet');
 end
@@ -86 +93 @@
 
 %treat the case datasize(1)=6*nodes
-if datasize(1)==6*md.mesh.numberofvertices
+if datasize(1)==6*numberofvertices
         %keep the only norm of data
-        data1=data(1:6:md.mesh.numberofvertices*6,:);
-        data2=data(2:6:md.mesh.numberofvertices*6,:);
+        data1=data(1:6:numberofvertices*6,:);
+        data2=data(2:6:numberofvertices*6,:);
         data=sqrt(data1.^2+data2.^2);
-        datasize(1)=md.mesh.numberofvertices;
+        datasize(1)=numberofvertices;
         %---> go to node data
 end
@@ -98 +105 @@
 if (dimension(md.mesh)==3 & datasize(1)==numberofvertices2d),
         data=project3d(md,'vector',data,'type','node');
-        datasize(1)=md.mesh.numberofvertices;
+        datasize(1)=numberofvertices;
         %---> go to node data
 end
@@ -105 +112 @@
 if (dimension(md.mesh)==3 & datasize(1)==numberofelements2d),
         data=project3d(md,'vector',data,'type','element');
-        datasize(1)=md.mesh.numberofelements;
+        datasize(1)=numberofelements;
         %---> go to node data
 end
@@ -112 +119 @@
 if exist(options,'smooth')
         data=averaging(md,data,getfieldvalue(options,'smooth'));
-        datasize(1)=md.mesh.numberofvertices;
+        datasize(1)=numberofvertices;
         %---> go to node data
 end
 
 %element data
-if (datasize(1)==md.mesh.numberofelements & datasize(2)==1),
+if (datasize(1)==numberofelements & datasize(2)==1),
 
         %Initialize datatype if non patch
@@ -129 +136 @@
                 maskvalue=getfieldvalue(options,'maskvalue',NaN);
                 pos=find(~flags);
-                if length(flags)==md.mesh.numberofvertices,
+                if length(flags)==numberofvertices,
                         [pos2 dummy]=find(ismember(md.mesh.elements,pos));
                         data(pos2,:)=maskvalue;
-                elseif length(flags)==md.mesh.numberofelements
+                elseif length(flags)==numberofelements
                         data(pos,:)=maskvalue;
                 else
-                        disp('plotmodel warning: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements');
+                        disp('plotmodel warning: mask length not supported yet (supported length are numberofvertices and numberofelements');
                 end
         end
@@ -152 +159 @@
 
 %node data
-if (datasize(1)==md.mesh.numberofvertices & datasize(2)==1),
+if (datasize(1)==numberofvertices & datasize(2)==1),
         datatype=2;
 
@@ -160 +167 @@
                 maskvalue=getfieldvalue(options,'maskvalue',NaN);
                 pos=find(~flags);
-                if length(flags)==md.mesh.numberofvertices,
+                if length(flags)==numberofvertices,
                         data(pos,:)=maskvalue;
-                elseif length(flags)==md.mesh.numberofelements
+                elseif length(flags)==numberofelements
                         data(md.mesh.elements(pos,:),:)=maskvalue;
                 else
-                        disp('plotmodel warning: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements');
+                        disp('plotmodel warning: mask length not supported yet (supported length are numberofvertices and numberofelements');
                 end
         end
@@ -199 +206 @@
                 flags=getfieldvalue(options,'mask');
                 pos=find(~flags);
-                if length(flags)==md.mesh.numberofvertices,
+                if length(flags)==numberofvertices,
                         data(pos,:)=NaN;
-                elseif length(flags)==md.mesh.numberofelements
+                elseif length(flags)==numberofelements
                         data(md.mesh.elements(pos,:),:)=NaN;
                 else
-                        disp('plotmodel warning: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements');
+                        disp('plotmodel warning: mask length not supported yet (supported length are numberofvertices and numberofelements');
                 end
         end

issm/trunk-jpl/src/m/plot/processmesh.m

@@ -27 +27 @@
 end
 
-if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon'),
-        x=md.mesh.x;
-        if isprop(md.mesh,'x2d'), x2d=md.mesh.x2d; end
-        y=md.mesh.y;
-        if isprop(md.mesh,'y2d'), y2d=md.mesh.y2d; end
+if isprop(md.mesh,'elements2d'), elements2d=md.mesh.elements2d; end
+
+if exist(options,'amr'),
+        step = getfieldvalue(options,'amr');
+        x = md.results.TransientSolution(step).MeshX;
+        y = md.results.TransientSolution(step).MeshY;
+        elements = md.results.TransientSolution(step).MeshElements;
 else
-        x=md.mesh.long;
-        y=md.mesh.lat;
+        elements=md.mesh.elements;
+        if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon') ,
+                x=md.mesh.x;
+                if isprop(md.mesh,'x2d'), x2d=md.mesh.x2d; end
+                y=md.mesh.y;
+                if isprop(md.mesh,'y2d'), y2d=md.mesh.y2d; end
+        else
+                x=md.mesh.long;
+                y=md.mesh.lat;
+        end
 end
 
@@ -46 +56 @@
         z=md.(z);
 end
-
-if isprop(md.mesh,'elements2d'), elements2d=md.mesh.elements2d; end
-elements=md.mesh.elements;
 
 %is it a 2d plot?