Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 5441) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 5442) @@ -84,5 +84,7 @@ fields={'edges'}; checksize(md,fields,[NaN 4]); - checknan(md,fields); + fields={'edges(:,1:3)'}; + checknan(md,fields); + checkgreaterstrict(md,fields,0); end %}}} @@ -429,5 +431,6 @@ %CHECKSIZE - check length of a field for i=1:length(fields), - if length(eval(['md.' fields{i}]))~=fieldlength, + eval(['field=md.' fields{i} ';']); + if length(field)~=fieldlength, error(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]); end @@ -439,30 +442,15 @@ %CHECKSIZE - check size of a field for i=1:length(fields), + eval(['field=md.' fields{i} ';']); if isnan(fieldsize(1)), - if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)), - %LOG warnings: - if strcmpi(fields{i},'pressureload'), - disp(' '); - disp(' If pressureload does not have the right size (one column missing), this is due to a change in md.pressure.'); - disp(' To update your model, use the following line'); - disp(' md.pressureload=[md.pressureload WaterEnum*md.elementoniceshelf(md.pressureload(:,end))+AirEnum*md.elementonicesheet(md.pressureload(:,end))];'); - disp(' '); - end + if (size(field,2)~=fieldsize(2)), error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']); end elseif isnan(fieldsize(2)), - if (size(eval(['md.' fields{i}]),1)~=fieldsize(1)), + if (size(field,1)~=fieldsize(1)), error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']); end else - if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) | (size(eval(['md.' fields{i}]),2)~=fieldsize(2))) - %LOG warnings: - if strcmpi(fields{i},'weights'), - disp(' '); - disp(' ''weights'' is a new field that should be of length numberofgrids.'); - disp(' To update your model, use the following line'); - disp(' md.weights=ones(md.numberofgrids,1)'); - disp(' '); - end + if ((size(field)~=fieldsize(1)) | (size(field,2)~=fieldsize(2))) error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]); end @@ -475,6 +463,7 @@ %CHECKNAN - check nan values of a field for i=1:length(fields), - if any(isnan(eval(['md.' fields{i}]))), - error(['model not consistent: NaN values in field ' fields{i}]); + eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); + if any(isnan(field)), + error(['model not consistent: NaN values found in field ' fields{i}]); end end @@ -485,6 +474,7 @@ %CHECKREAL - check real values of a field for i=1:length(fields), - if any(eval(['~isreal(md.' fields{i} ')'])), - error(['model not consistent: complex values in field ' fields{i}]); + eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); + if any(~isreal(field)), + error(['model not consistent: complex values found in field ' fields{i}]); end end @@ -495,5 +485,6 @@ %CHECKGREATERSTRICT - check values of a field for i=1:length(fields), - if any(eval(['md.' fields{i} '<=' num2str(lowerbound) ])), + eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); + if any(field<=lowerbound), error(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]); end @@ -505,5 +496,6 @@ %CHECKGREATER - check values of a field for i=1:length(fields), - if any(eval(['md.' fields{i} '<' num2str(lowerbound) ])), + eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); + if any(field=' num2str(upperbound) ])), + eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); + if any(field>=upperbound), error(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]); end @@ -525,5 +518,6 @@ %CHECKLESS - check values of a field for i=1:length(fields), - if any(eval(['md.' fields{i} '>' num2str(upperbound) ])), + eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); + if any(field>upperbound), error(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]); end @@ -535,5 +529,6 @@ %CHECKVALUE - check that a field has specified values for i=1:length(fields), - if eval(['any(~ismember( md.' fields{i} ',values))']), + eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); + if any(~ismember(field,values)), error(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]); end Index: /issm/trunk/src/m/classes/public/modelextract.m =================================================================== --- /issm/trunk/src/m/classes/public/modelextract.m (revision 5441) +++ /issm/trunk/src/m/classes/public/modelextract.m (revision 5442) @@ -148,16 +148,24 @@ if size(md2.edges,2)>1, %do not use ~isnan because there are some NaNs... %renumber first two columns - md2.edges(:,1)=Pgrid(md2.edges(:,1)); - md2.edges(:,2)=Pgrid(md2.edges(:,2)); - md2.edges(:,3)=Pelem(md2.edges(:,3)); - %remove column where the first element is not in the domain - md2.edges=md2.edges(find(md2.edges(:,1) & md2.edges(:,2) & md2.edges(:,3)),:); - %now we must renumber the second column and put NaN if the element is not in the domain - A=md2.edges(:,4); - A(find(~isnan(A)))=Pelem(A(find(~isnan(A)))); - A(find(A==0))=NaN; - md2.edges(:,4)=A; clear A; - end - + pos=find(~isnan(md2.edges(:,4))); + md2.edges(: ,1)=Pgrid(md2.edges(:,1)); + md2.edges(: ,2)=Pgrid(md2.edges(:,2)); + md2.edges(: ,3)=Pelem(md2.edges(:,3)); + md2.edges(pos,4)=Pelem(md2.edges(pos,4)); + %remove edges when the 2 vertices are not in the domain. + md2.edges=md2.edges(find(md2.edges(:,1) & md2.edges(:,2)),:); + %Replace all zeros by NaN in the last two columns; + pos=find(md2.edges(:,3)==0); + md2.edges(pos,3)=NaN; + pos=find(md2.edges(:,4)==0); + md2.edges(pos,4)=NaN; + %Invert NaN of the third column with last column + pos=find(isnan(md2.edges(:,3))); + md2.edges(pos,3)=md2.edges(pos,4); + md2.edges(pos,4)=NaN; + %Finally remove edges that do not belong to any element + pos=find(isnan(md2.edges(:,3)) & isnan(md2.edges(:,4))); + md2.edges(pos,:)=[]; + end %Penalties