![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 13495 for issm/trunk-jpl/src/m/mesh/bamg.py
- Timestamp:
- 10/01/12 12:06:53 (12 years ago)
- File:
-
- 1 edited
-
issm/trunk-jpl/src/m/mesh/bamg.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/mesh/bamg.py
r13493 r13495 7 7 from expread import * 8 8 from expwrite import * 9 from SegIntersect import * 9 10 from MatlabFuncs import * 10 11 from BamgMesher import * … … 121 122 122 123 elif numpy.any(numpy.logical_not(flags)): 123 raise RuntimeError("bamg.py for rifts is not complete.")124 124 #We LOTS of work to do 125 125 print("Rift tip outside of or on the domain has been detected and is being processed...") 126 """ 127 128 %check that only one point is outside (for now) 129 if sum(~flags)~=1, 130 error('bamg error message: only one point outside of the domain is supported yet'); 131 end 132 133 %Move tip outside to the first position 134 if flags(1)==0, 135 %OK, first point is outside (do nothing), 136 elseif (flags(end)==0), 137 rift(i).x=flipud(rift(i).x); 138 rift(i).y=flipud(rift(i).y); 139 else 140 error('bamg error message: only a rift tip can be outside of the domain'); 141 end 142 143 %Get cordinate of intersection point 144 x1=rift(i).x(1); y1=rift(i).y(1); 145 x2=rift(i).x(2); y2=rift(i).y(2); 146 for j=1:length(domain[0]['x'])-1; 147 if SegIntersect([x1 y1; x2 y2],[domain[0]['x'](j) domain[0]['y'](j); domain[0]['x'](j+1) domain[0]['y'](j+1)]), 148 149 %Get position of the two nodes of the edge in domain 150 i1=j; 151 i2=mod(j+1,domain[0]['nods']); 152 153 %rift is crossing edge [i1 i2] of the domain 154 %Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) 155 x3=domain[0]['x'](i1); y3=domain[0]['y'](i1); 156 x4=domain[0]['x'](i2); y4=domain[0]['y'](i2); 157 x=det([det([x1 y1; x2 y2]) x1-x2;det([x3 y3; x4 y4]) x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); 158 y=det([det([x1 y1; x2 y2]) y1-y2;det([x3 y3; x4 y4]) y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); 159 160 segdis= sqrt((x4-x3)**2+(y4-y3)**2); 161 tipdis=[sqrt((x-x3)**2+(y-y3)**2) sqrt((x-x4)**2+(y-y4)**2)]; 162 163 if (min(tipdis)/segdis) < getfieldvalue(options,'toltip',0), 164 disp('moving tip-domain intersection point'); 165 166 %Get position of the closest point 167 if tipdis(1)>tipdis(2), 168 pos=i2; 169 else 170 pos=i1; 171 end 172 173 %This point is only in Vertices (number pos). 174 %OK, now we can add our own rift 175 nods=rift(i).nods-1; 176 bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]]; 177 bamg_geometry.Edges=[bamg_geometry.Edges;... 178 pos count+1 (1+i);... 179 [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) (1+i)*ones(nods-1,1)]]; 126 127 #check that only one point is outside (for now) 128 if numpy.sum(numpy.logical_not(flags).astype(int))!=1: 129 raise RuntimeError("bamg error message: only one point outside of the domain is supported yet") 130 131 #Move tip outside to the first position 132 if not flags[0]: 133 #OK, first point is outside (do nothing), 134 pass 135 elif not flags[-1]: 136 rifti['x']=numpy.flipud(rifti['x']) 137 rifti['y']=numpy.flipud(rifti['y']) 138 else: 139 raise RuntimeError("bamg error message: only a rift tip can be outside of the domain") 140 141 #Get cordinate of intersection point 142 x1=rifti['x'][0] 143 y1=rifti['y'][0] 144 x2=rifti['x'][1] 145 y2=rifti['y'][1] 146 for j in xrange(0,numpy.size(domain[0]['x'])-1): 147 if SegIntersect(numpy.array([[x1,y1],[x2,y2]]),numpy.array([[domain[0]['x'][j],domain[0]['y'][j]],[domain[0]['x'][j+1],domain[0]['y'][j+1]]])): 148 149 #Get position of the two nodes of the edge in domain 150 i1=j 151 i2=j+1 152 153 #rift is crossing edge [i1 i2] of the domain 154 #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) 155 x3=domain[0]['x'][i1] 156 y3=domain[0]['y'][i1] 157 x4=domain[0]['x'][i2] 158 y4=domain[0]['y'][i2] 159 # x=det([det([x1 y1; x2 y2]) x1-x2;det([x3 y3; x4 y4]) x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); 160 # y=det([det([x1 y1; x2 y2]) y1-y2;det([x3 y3; x4 y4]) y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); 161 x=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),x3-x4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])); 162 y=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),y3-y4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])); 163 164 segdis= sqrt((x4-x3)**2+(y4-y3)**2) 165 tipdis=numpy.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)]) 166 167 if numpy.min(tipdis)/segdis < options.getfieldvalue('toltip',0): 168 print("moving tip-domain intersection point") 169 170 #Get position of the closer point 171 if tipdis[0]>tipdis[1]: 172 pos=i2 173 else: 174 pos=i1 175 176 #This point is only in Vertices (number pos). 177 #OK, now we can add our own rift 178 nods=rifti['nods']-1 179 bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1)))))) 180 bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\ 181 numpy.array([[pos,count+1,(1+i)]]),\ 182 numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1)))))) 183 count+=nods 184 185 break 186 187 else: 188 #Add intersection point to Vertices 189 bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.array([[x,y,1]]))) 190 count+=1 191 192 #Decompose the crossing edge into 2 subedges 193 pos=numpy.nonzero(numpy.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0] 194 if not pos: 195 raise RuntimeError("bamg error message: a problem occurred...") 196 bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges[0:pos-1,:],\ 197 numpy.array([[bamg_geometry.Edges[pos,0],count ,bamg_geometry.Edges[pos,2]]]),\ 198 numpy.array([[count ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\ 199 bamg_geometry.Edges[pos+1:,:])) 200 201 #OK, now we can add our own rift 202 nods=rifti['nods']-1 203 bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1)))))) 204 bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\ 205 numpy.array([[count,count+1,2]]),\ 206 numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1)))))) 180 207 count=count+nods; 181 208 182 break; 183 184 else 185 %Add intersection point to Vertices 186 bamg_geometry.Vertices=[bamg_geometry.Vertices; x y 1]; 187 count=count+1; 188 189 %Decompose the crossing edge in 2 subedges 190 pos=find(bamg_geometry.Edges(:,1)==i1 & bamg_geometry.Edges(:,2)==i2); 191 if isempty(pos) error('bamg error message: a problem occured...'); end 192 bamg_geometry.Edges=[bamg_geometry.Edges(1:pos-1,:);... 193 bamg_geometry.Edges(pos,1) count bamg_geometry.Edges(pos,3);... 194 count bamg_geometry.Edges(pos,2) bamg_geometry.Edges(pos,3);... 195 bamg_geometry.Edges(pos+1:end,:)]; 196 197 %OK, now we can add our own rift 198 nods=rift(i).nods-1; 199 bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]]; 200 bamg_geometry.Edges=[bamg_geometry.Edges;... 201 count count+1 2 ;... 202 [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) (1+i)*ones(nods-1,1)]]; 203 count=count+nods; 204 205 break; 206 end 207 end 208 end 209 """ 209 break 210 210 211 else: 211 212 nods=rifti['nods']-1 … … 384 385 385 386 #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) 386 x= det(numpy.array([det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[det(numpy.array([[x3,y3],[x4,y4]])),x3-x4])/det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))387 y= det(numpy.array([det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[det(numpy.array([[x3,y3],[x4,y4]])),y3-y4])/det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))387 x=numpy.linalg.det(numpy.array([numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),x3-x4])/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))) 388 y=numpy.linalg.det(numpy.array([numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),y3-y4])/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))) 388 389 389 390 #Add vertex to the list of vertices
Note:
See TracChangeset
for help on using the changeset viewer.
