Changeset 25455 for issm/trunk-jpl/src/m/mesh/bamg.m

Timestamp:

08/25/20 00:32:13 (5 years ago)

Author:

jdquinn

Message:

CHG: Saving chnages so that Basile has access to potential fix to solidearthmodel class

File:

• issm/trunk-jpl/src/m/mesh/bamg.m (modified) (14 diffs)

issm/trunk-jpl/src/m/mesh/bamg.m

@@ -80 +80 @@
                         domain=shpread(domainfile);
                 else
-                        error(['bamg error message: file ' domainfile ' format not supported (.shp or .exp)']);
+                        error(['bamg error message: file ' domainfile ' format not supported (.exp or .shp)']);
                 end
         elseif isstruct(domainfile),
@@ -140 +140 @@
                 end
 
-                %Checks that all holes are INSIDE the principle domain outline
+                %Check that all holes are INSIDE the principle domain outline
                 if i>1,
                         flags=ContourToNodes(domain(i).x,domain(i).y,domain(1),0);
@@ -149 +149 @@
 
                 %Check orientation
-                nods=domain(i).nods-1; %the domain are closed 1=end;
+                nods=domain(i).nods-1; %the domain is closed (domain[1] = domain[end])
                 test = sum([(domain(i).x(2:nods+1) - domain(i).x(1:nods)).*(domain(i).y(2:nods+1) + domain(i).y(1:nods))]);
                 if (i==1 && test>0) || (i>1 && test<0),
@@ -163 +163 @@
                 bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1*ones(nods,1)]];
 
-                %Flag how many edges we have now, that way we know which edges belong to the domain, will 
-                %be used later for required edges when NoBoundaryRefinement is one: 
+                % Flag how many edges we have now, that way we know which edges belong 
+                % to the domain. Will be used later for required edges if 
+                % NoBoundaryRefinement equals 1.
                 new_edge_length=length(bamg_geometry.Edges);
                 edges_required=(edge_length+1):new_edge_length;
 
                 if i>1,
-                        bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref]; subdomain_ref = subdomain_ref+1;
+                        bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -subdomain_ref]; 
+                        subdomain_ref = subdomain_ref+1;
                 else
                         bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 0];
@@ -179 +181 @@
         for i=1:length(holes),
 
-                %Check that the subdomains is closed
+                %Check that the hole is closed
                 if (holes(i).x(1)~=holes(i).x(end) | holes(i).y(1)~=holes(i).y(end)),
                         error('bamg error message: all contours provided in ''holes'' should be closed');
                 end
 
-                %Checks that all holes are INSIDE the principle domain outline
+                %Check that all holes are INSIDE the principal domain (principal domain should be index 0)
                 flags=ContourToNodes(holes(i).x,holes(i).y,domain(1),0);
                 if any(~flags), error('bamg error message: All holes should be strictly inside the principal domain'); end
 
                 %Check that hole is correctly oriented
-                nods=holes(i).nods-1; %the holes are closed 1=end;
+                nods=holes(i).nods-1; %the hole is closed (hole[1] = hole[end])
                 if(sum([(holes(i).x(2:nods+1) - holes(i).x(1:nods)).*(holes(i).y(2:nods+1) + holes(i).y(1:nods))]))<0
                         disp('At least one hole was not correctly oriented and has been re-oriented');
@@ -205 +207 @@
         for i=1:length(subdomains),
 
-                %Check that the subdomains is closed
+                %Check that the subdomain is closed
                 if (subdomains(i).x(1)~=subdomains(i).x(end) | subdomains(i).y(1)~=subdomains(i).y(end)),
                         error('bamg error message: all contours provided in ''subdomains'' should be closed');
                 end
 
-                %Checks that all holes are INSIDE the principle domain outline
+                %Checks that all subdomains are INSIDE the principal domain (principal domain should be index 0)
                 flags=ContourToNodes(subdomains(i).x,subdomains(i).y,domain(1),0);
                 if any(~flags),
-                        error('bamg error message: All holes should be strictly inside the principal domain');
-                end
-
-                %Check that hole is correctly oriented
-                nods=subdomains(i).nods-1; %the subdomains are closed 1=end;
+                        error('bamg error message: All subdomains should be strictly inside the principal domain');
+                end
+
+                %Check that subdomain is correctly oriented
+                nods=subdomains(i).nods-1; % the subdomains are closed (subdomains[1] = subdomains[end])
                 if(sum([(subdomains(i).x(2:nods+1) - subdomains(i).x(1:nods)).*(subdomains(i).y(2:nods+1) + subdomains(i).y(1:nods))]))>0
                         disp('At least one subdomain was not correctly oriented and has been re-oriented');
-                        subdomains(i).x = flipud(subdomains(i).x); subdomains(i).y = flipud(subdomains(i).y);
+                        subdomains(i).x = flipud(subdomains(i).x); 
+                        subdomains(i).y = flipud(subdomains(i).y);
                 end
 
@@ -227 +230 @@
                 bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1.*ones(nods,1)]];
                 
-                bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref]; subdomain_ref = subdomain_ref+1;
+                bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref];
+                subdomain_ref = subdomain_ref+1;
                 
                 %update counter
                 count=count+nods;
         end
+
         if getfieldvalue(options,'vertical',0),
                 if numel(getfieldvalue(options,'Markers',[]))~=size(bamg_geometry.Edges,1),
@@ -243 +248 @@
         %take care of rifts
         if exist(options,'rifts'),
-
-                %Check that file exists
+                %read rift file according to its extension: 
                 riftfile=getfieldvalue(options,'rifts');
-                [pathr,namer,extr]=fileparts(riftfile);
-                if ~exist(riftfile,'file')
-                        error(['bamg error message: file ' riftfile ' not found ']);
-                elseif strcmp(extr,'.exp'),
-                        rift=expread(riftfile);
-                elseif strcmp(extr,'.shp'),
-                        rift=shpread(riftfile);
-                end
-                %read rift file according to its extension: 
                 [path,name,ext]=fileparts(riftfile);
                 if strcmp(ext,'.exp'),
@@ -261 +256 @@
                         rift=shpread(riftfile);
                 else
-                        error(['bamg error message: file ' riftfile ' format not supported (.shp or .exp)']);
+                        error(['bamg error message: file ' riftfile ' format not supported (.exp or .shp)']);
                 end
 
@@ -272 +267 @@
 
                         elseif any(~flags),
-                                %We LOTS of work to do
+                                %We have LOTS of work to do
                                 disp('Rift tip outside of or on the domain has been detected and is being processed...');
 
                                 %check that only one point is outside (for now)
                                 if sum(~flags)~=1,
-                                        error('bamg error message: only one point outside of the domain is supported yet');
+                                        error('bamg error message: only one point outside of the domain is supported at this time');
                                 end
 
@@ -290 +285 @@
                                 end
 
-                                %Get cordinate of intersection point
-                                x1=rift(i).x(1); y1=rift(i).y(1);
-                                x2=rift(i).x(2); y2=rift(i).y(2);
+                                %Get coordinate of intersection point
+                                x1=rift(i).x(1);
+                                y1=rift(i).y(1);
+                                x2=rift(i).x(2);
+                                y2=rift(i).y(2);
                                 for j=1:length(domain(1).x)-1;
                                         if SegIntersect([x1 y1; x2 y2],[domain(1).x(j) domain(1).y(j); domain(1).x(j+1) domain(1).y(j+1)]),
@@ -489 +486 @@
         md.mesh.numberofvertices=length(md.mesh.x);
         md.mesh.numberofedges=size(md.mesh.edges,1);
-        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-
-        %Now, build the connectivity tables for this mesh.
-        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 
 elseif getfieldvalue(options,'3dsurface',0),
@@ -509 +504 @@
         md.mesh.numberofvertices=length(md.mesh.x);
         md.mesh.numberofedges=size(md.mesh.edges,1);
-        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 
 else 
@@ -524 +520 @@
         md.mesh.numberofvertices=length(md.mesh.x);
         md.mesh.numberofedges=size(md.mesh.edges,1);
-        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 end