Changeset 22212

Timestamp:

11/02/17 14:18:56 (7 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added holes and subdomains options, provided as extra .exp files

File:

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

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

@@ -5 +5 @@
 %
 %   - domain :            followed by an ARGUS file that prescribes the domain outline
+%   - holes :             followed by an ARGUS file that prescribes the holes
+%   - subdomains :        followed by an ARGUS file that prescribes the list of
+%                         subdomains (that need to be inside domain)
+%
 %   - hmin :              minimum edge length (default is 10^-100)
 %   - hmax :              maximum edge length (default is 10^100)
@@ -59 +63 @@
 
 subdomain_ref = 1;
+hole_ref = 1;
 % Bamg Geometry parameters {{{
 if exist(options,'domain'),
@@ -82 +87 @@
         end
 
+        holes = [];
+        if exist(options,'holes'),
+                holesfile=getfieldvalue(options,'holes');
+                if ischar(holesfile),
+                        if ~exist(holesfile,'file') error(['bamg error message: file ' holesfile ' not found']); end
+
+                        %read holes according to its extension: 
+                        [path,name,ext]=fileparts(holesfile);
+                        if strcmp(ext,'.exp'),
+                                holes=expread(holesfile);
+                        elseif strcmp(ext,'.shp'),
+                                holes=shpread(holesfile);
+                        else
+                                error(['bamg error message: file ' holesfile ' format not supported (.shp or .exp)']);
+                        end
+                elseif isstruct(holesfile),
+                        holes = holesfile;
+                else
+                        error('''holes'' type not supported yet');
+                end
+        end
+        subdomains = [];
+        if exist(options,'subdomains'),
+                subdomainsfile=getfieldvalue(options,'subdomains');
+                if ischar(subdomainsfile),
+                        if ~exist(subdomainsfile,'file') error(['bamg error message: file ' subdomainsfile ' not found']); end
+
+                        %read subdomains according to its extension: 
+                        [path,name,ext]=fileparts(subdomainsfile);
+                        if strcmp(ext,'.exp'),
+                                subdomains=expread(subdomainsfile);
+                        elseif strcmp(ext,'.shp'),
+                                subdomains=shpread(subdomainsfile);
+                        else
+                                error(['bamg error message: file ' subdomainsfile ' format not supported (.shp or .exp)']);
+                        end
+                elseif isstruct(subdomainsfile),
+                        subdomains = subdomainsfile;
+                else
+                        error('''subdomains'' type not supported yet');
+                end
+        end
+
         %Build geometry 
         count=0;
@@ -99 +147 @@
                 end
 
+                %Check orientation
+                nods=domain(i).nods-1; %the domain are closed 1=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),
+                        disp('At least one contour was not correctly oriented and has been re-oriented');
+                        domain(i).x = flipud(domain(i).x); domain(i).y = flipud(domain(i).y);
+                end
+
                 %Add all points to bamg_geometry
-                nods=domain(i).nods-1; %the domain are closed 1=end;
                 bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
+                bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1*ones(nods,1)]];
+                if i>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];
+                end
+
+                %update counter
+                count=count+nods;
+        end
+        for i=1:length(holes),
+
+                %Check that the subdomains 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
+                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;
+                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');
+                        holes(i).x = flipud(holes(i).x); holes(i).y = flipud(holes(i).y);
+                end
+
+                %Add all points to bamg_geometry
+                bamg_geometry.Vertices=[bamg_geometry.Vertices; [holes(i).x(1:nods) holes(i).y(1:nods) ones(nods,1)]];
                 bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1.*ones(nods,1)]];
-                if i>1, bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 subdomain_ref]; subdomain_ref = subdomain_ref+1; end
+                bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 -hole_ref]; hole_ref = hole_ref+1;
+
+                %update counter
+                count=count+nods;
+        end
+        for i=1:length(subdomains),
+
+                %Check that the subdomains 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
+                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;
+                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);
+                end
+
+                %Add all points to bamg_geometry
+                bamg_geometry.Vertices=[bamg_geometry.Vertices; [subdomains(i).x(1:nods) subdomains(i).y(1:nods) ones(nods,1)]];
+                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;
 
                 %update counter