Changeset 22299

Timestamp:

12/22/17 08:13:13 (7 years ago)

Author:

bdef

Message:

NEW: implemeting domain orientation in python bamg (holes and subdomain too)

File:

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

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

@@ -22 +22 @@
 
            - 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)
@@ -74 +78 @@
         bamg_mesh=bamgmesh()
 
-        subdomain_ref = 1
-        hole_ref = 1
+        subdomain_ref = 1
+        hole_ref = 1
 
         # Bamg Geometry parameters {{{
         if options.exist('domain'):
-
                 #Check that file exists
                 domainfile=options.getfieldvalue('domain')
@@ -92 +95 @@
                 count=0
                 for i,domaini in enumerate(domain):
-
                         #Check that the domain is closed
                         if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]):
                                 raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")
 
-                        #Checks that all holes are INSIDE the principle domain outline
+                        #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
                         if i:
                                 flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0]
@@ -103 +105 @@
                                         raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
 
+                        #Check orientation
+                        nods=domaini['nods']-1#the domain are closed 1=end;
+
+                        test = np.sum((domaini['x'][1:nods+1] - domaini['x'][0:nods])*(domaini['y'][1:nods+1] + domaini['y'][0:nods]))
+                        if (i==0 and test>0) or (i>0 and test<0):
+                                print('At least one contour was not correctly oriented and has been re-oriented')
+                                domaini['x'] = np.flipud(domaini['x'])
+                                domaini['y'] = np.flipud(domaini['y'])
+
+
                         #Add all points to bamg_geometry
                         nods=domaini['nods']-1    #the domain are closed 0=end
                         bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T))
                         bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
-                        if i:
-                            bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,-subdomain_ref]))
-                            subdomain_ref = subdomain_ref+1;
-                        else:
-                            bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,0]))
-
-                        # bamg_geometry.Vertices=np.hstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods].reshape(-1),domaini['y'][0:nods].reshape(-1),np.ones((nods))))))
-                        # bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.hstack((np.arange(count+1,count+nods+1).reshape(-1,),np.hstack((np.arange(count+2,count+nods+1),count+1)).reshape(-1,),1.*np.ones((nods,))))))
-                        # if i:
-                        #       bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,1]))
-
+                        if i:
+                                bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,-subdomain_ref]))
+                                subdomain_ref = subdomain_ref+1;
+                        else:
+                                bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,0]))
+                        
                         #update counter
                         count+=nods
 
+                #Deal with domain holes
+                if options.exist('holes'):
+                        holesfile=options.getfieldvalue('holes')
+                        if type(holesfile) == str:
+                                if not os.path.exists(holesfile):
+                                        raise IOError("bamg error message: file '%s' not found" % holesfile)
+                                holes=expread(holesfile)
+                        else:
+                                holes=holesfile
+
+                        #Build geometry 
+                        for i,holei in enumerate(holes):
+                                #Check that the hole is closed
+                                if (holei['x'][0]!=holei['x'][-1] or holei['y'][0]!=holei['y'][-1]):
+                                        raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed")
+
+                                #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
+                                flags=ContourToNodes(holei['x'],holei['y'],domainfile,0)[0]
+                                if np.any(np.logical_not(flags)):
+                                        raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
+
+                                #Check orientation
+                                nods=holei['nods']-1#the hole are closed 1=end;
+                                test = np.sum((holei['x'][1:nods+1] - holei['x'][0:nods])*(holei['y'][1:nods+1] + holei['y'][0:nods]))
+                                if test<0:
+                                        print('At least one hole was not correctly oriented and has been re-oriented')
+                                        holei['x'] = np.flipud(holei['x'])
+                                        holei['y'] = np.flipud(holei['y'])
+
+                                #Add all points to bamg_geometry
+                                nods=holei['nods']-1    #the hole are closed 0=end
+                                bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((holei['x'][0:nods],holei['y'][0:nods],np.ones((nods)))).T))
+                                bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
+                                #update counter
+                                count+=nods
+
+                #And subdomains
+                if options.exist('subdomains'):
+                        subdomainfile=options.getfieldvalue('subdomains')
+                        if type(subdomainfile) == str:
+                                if not os.path.exists(subdomainfile):
+                                        raise IOError("bamg error message: file '%s' not found" % subdomainfile)
+                                subdomains=expread(subdomainfile)
+                        else:
+                                subdomains=subdomainfile
+
+                        #Build geometry 
+                        for i,subdomaini in enumerate(subdomains):
+                                #Check that the subdomain is closed
+                                if (subdomaini['x'][0]!=subdomaini['x'][-1] or subdomaini['y'][0]!=subdomaini['y'][-1]):
+                                        raise RuntimeError("bamg error message: all contours provided in ''subdomain'' should be closed")
+
+                                #Checks that all subdomains are INSIDE the principle subdomain outline princial domain should be index 0
+                                if i:
+                                        flags=ContourToNodes(subdomaini['x'],subdomaini['y'],domainfile,0)[0]
+                                        if np.any(np.logical_not(flags)):
+                                                raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal subdomain")
+
+                                #Check orientation
+                                nods=subdomaini['nods']-1#the subdomain are closed 1=end;
+
+                                test = np.sum((subdomaini['x'][1:nods+1] - subdomaini['x'][0:nods])*(subdomaini['y'][1:nods+1] + subdomaini['y'][0:nods]))
+                                if test>0:
+                                        print('At least one subcontour was not correctly oriented and has been re-oriented')
+                                        subdomaini['x'] = np.flipud(subdomaini['x'])
+                                        subdomaini['y'] = np.flipud(subdomaini['y'])
+
+                                #Add all points to bamg_geometry
+                                nods=subdomaini['nods']-1    #the subdomain are closed 0=end
+                                bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((subdomaini['x'][0:nods],subdomaini['y'][0:nods],np.ones((nods)))).T))
+                                bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
+                                #update counter
+                                count+=nods
+                        
                 if options.getfieldvalue('vertical',0):
                         if np.size(options.getfieldvalue('Markers',[]))!=np.size(bamg_geometry.Edges,0):