Changeset 22758 for issm/trunk/src/m/mesh/bamg.py

Timestamp:

05/10/18 10:24:27 (7 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 22757

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/m/mesh/bamg.py (modified) (10 diffs)

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

@@ -1 +1 @@
 import os.path
 import numpy as  np
-from mesh2d import mesh2d
+from mesh2d import *
+from mesh2dvertical import *
+from mesh3dsurface import *
 from collections import OrderedDict
 from pairoptions import pairoptions
@@ -20 +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)
@@ -37 +43 @@
                                  1 -> use Green formula
            - KeepVertices :      try to keep initial vertices when adaptation is done on an existing mesh (default 1)
-           - MaxCornerAngle :    maximum angle of corners in degree (default is 10)
            - maxnbv :            maximum number of vertices used to allocate memory (default is 10^6)
            - maxsubdiv :         maximum subdivision of exisiting elements (default is 10)
@@ -49 +54 @@
            - power :             power applied to the metric (default is 1)
            - splitcorners :      split triangles whuch have 3 vertices on the outline (default is 1)
-           - geometricalmetric : take the geometry into account to generate the metric (default is 0)
            - verbose :           level of verbosity (default is 1)
 
@@ -74 +78 @@
         bamg_mesh=bamgmesh()
 
+        subdomain_ref = 1
+        hole_ref = 1
+
         # Bamg Geometry parameters {{{
         if options.exist('domain'):
-
                 #Check that file exists
                 domainfile=options.getfieldvalue('domain')
-                if not os.path.exists(domainfile):
-                        raise IOError("bamg error message: file '%s' not found" % domainfile)
-                domain=expread(domainfile)
+                if type(domainfile) == str:
+                        if not os.path.exists(domainfile):
+                                raise IOError("bamg error message: file '%s' not found" % domainfile)
+                        domain=expread(domainfile)
+                else:
+                        domain=domainfile
 
                 #Build geometry 
                 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]
                                 if np.any(np.logical_not(flags)):
                                         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
@@ -102 +120 @@
                         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,1]))
-
-                        # 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]))
-
+                                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):
+                                raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges,0)))
+                        if np.size(options.getfieldvalue('Markers',[]))==np.size(bamg_geometry.Edges,0):
+                                bamg_geometry.Edges[:,2]=options.getfieldvalue('Markers')
 
                 #take care of rifts
@@ -303 +399 @@
         bamg_options['hVertices']=options.getfieldvalue('hVertices',np.empty((0,1)))
         bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1)
-        bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.)
         bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6)
         bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.)
@@ -313 +408 @@
         bamg_options['power']=options.getfieldvalue('power',1.)
         bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1)
-        bamg_options['geometricalmetric']=options.getfieldvalue('geometricalmetric',0)
-        bamg_options['random']=options.getfieldvalue('rand',True)
         bamg_options['verbose']=options.getfieldvalue('verbose',1)
         #}}}
@@ -322 +415 @@
 
         # plug results onto model
+        if options.getfieldvalue('vertical',0):
+                md.mesh=mesh2dvertical()
+                md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+                md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+                md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+                md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+                md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+                md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+                #Fill in rest of fields:
+                md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+                md.mesh.numberofvertices=np.size(md.mesh.x)
+                md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+                md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+                md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+                #Now, build the connectivity tables for this mesh. Doubled in matlab for some reason
+                md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,)
+                md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=1
+
+        elif options.getfieldvalue('3dsurface',0):
+                md.mesh=mesh3dsurface()
+                md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+                md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+                md.mesh.z=md.mesh.x
+                md.mesh.z[:]=0
+                md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+                md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+                md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+                md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+                #Fill in rest of fields:
+                md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+                md.mesh.numberofvertices=np.size(md.mesh.x)
+                md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+                md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+                md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+        else:
+                md.mesh=mesh2d()
+                md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+                md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+                md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+                md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+                md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+                md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+                #Fill in rest of fields:
+                md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+                md.mesh.numberofvertices=np.size(md.mesh.x)
+                md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+                md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+                md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+        #Bamg private fields
         md.private.bamg=OrderedDict()
         md.private.bamg['mesh']=bamgmesh(bamgmesh_out)
         md.private.bamg['geometry']=bamggeom(bamggeom_out)
-        md.mesh = mesh2d()
-        md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
-        md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-        md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
-        md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
-        md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
-        md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
-
-        #Fill in rest of fields:
-        md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
-        md.mesh.numberofvertices=np.size(md.mesh.x)
-        md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
-        md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
-        md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
         md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
         md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0
@@ -345 +479 @@
         #Check for orphan
         if np.any(np.logical_not(np.in1d(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))):
-                raise RuntimeError("Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)")
+                raise RuntimeError("Output mesh has orphans. Check your Domain and/or RequiredVertices")
 
         return md