In [1]:
%matplotlib notebook
import numpy as np
bathy=np.genfromtxt("FRF_FRF_20150915_1116_NAVD88_LARC_GPS_UTC.csv",
delimiter=",",
skiprows=1)
print bathy.shape, bathy.dtype
In [2]:
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(bathy[:,7], bathy[:,8], bathy[:,9],c=bathy[:,10],linewidths=0)
Out[2]:
In [3]:
xmax=bathy[:,7].max()
xmin=bathy[:,7].min()
ymax=bathy[:,8].max()
ymin=bathy[:,8].min()
zmax=bathy[:,9].max()
zmin=bathy[:,9].min()
print xmin,xmax
print ymin,ymax
print zmin,zmax
In [4]:
xmin = 65.0
xmax = 450.0
ymin = 900.0
ymax = 900.1
In [5]:
from proteus import Profiling, Comm
comm = Comm.init()
Profiling.verbose=1
Profiling.procID=comm.rank()
Profiling.openLog("frfDomain",level=7)
from proteus.Domain import InterpolatedBathymetryDomain
bathy_points = np.vstack((bathy[:,7],bathy[:,8],bathy[:,9])).transpose()
print bathy_points.shape
domain = InterpolatedBathymetryDomain(vertices=[[xmin,ymin],[xmin,ymax],[xmax,ymax],[xmax,ymin]],
vertexFlags=[1,1,2,2],
segments=[[0,1],[1,2],[2,3],[3,0]],
segmentFlags=[1,2,3,4],
regions=[(0.5*(xmax+xmin),0.5*(ymax+ymin))],
regionFlags=[1],
name="frfDomain2D",
units='m',
bathy = bathy_points)
domain.writePoly(domain.name)
In [6]:
from proteus.MeshTools import InterpolatedBathymetryMesh
mesh = InterpolatedBathymetryMesh(domain,
triangleOptions="gVApq30Dena%8.8f" % ((1000.0**2)/2.0,),
maxElementDiameter=5.0,
atol=1.0e-2,
rtol=1.0e-2,
maxLevels=25,
maxNodes=100000,
bathyType="points",
bathyAssignmentScheme="interpolation",
errorNormType="Linfty")
In [7]:
fig = plt.figure()
nodes = mesh.meshList[-1].nodeArray
x,y,z = nodes[:,0],nodes[:,1],nodes[:,2]
elements = mesh.meshList[-1].elementNodesArray
plt.triplot(nodes[:,0],nodes[:,1],elements)
Out[7]:
In [8]:
fit = plt.figure()
plt.tricontourf(nodes[:,0],nodes[:,1],elements,nodes[:,2])
plt.colorbar()
Out[8]:
In [9]:
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(x, y, z, cmap=cm.jet, linewidth=0.0)
plt.show()
In [10]:
x.shape
Out[10]:
In [11]:
print mesh.meshList[-1].nodeMaterialTypes
In [12]:
fineMesh = mesh.meshList[-1]
newNodes = {}
verticalEdges = {}
nN_start = fineMesh.nodeArray.shape[0]
nN = nN_start
for nN_bottom, n,f in zip(range(nN_start), fineMesh.nodeArray, fineMesh.nodeMaterialTypes):
if f > 0:
newNodes[nN] = (n[0],n[1],1.5*zmax)
verticalEdges[nN_bottom] = nN
nN += 1
newFacets = []
newFacetFlags = []
for t in fineMesh.elementNodesArray:
newFacets.append([[t[0],t[1],t[2]]])
newFacetFlags.append(6)
topConnectivity = {}
for edge, edgeF in zip(fineMesh.elementBoundaryNodesArray, fineMesh.elementBoundaryMaterialTypes):
if edgeF > 0:
n00 = edge[0]
n10 = edge[1]
n11 = verticalEdges[edge[1]]
n01 = verticalEdges[edge[0]]
newFacets.append([[n00,
n10,
n11,
n01]])
newFacetFlags.append(edgeF)
if topConnectivity.has_key(n11):
topConnectivity[n11].append(n01)
else:
topConnectivity[n11] = [n01]
if topConnectivity.has_key(n01):
topConnectivity[n01].append(n11)
else:
topConnectivity[n01] = [n11]
topFacet = [topConnectivity[nN_start][0],nN_start,topConnectivity[nN_start][1]]
while len(topFacet) < len(newNodes):
nN = topFacet[-1]
if topFacet[-2] == topConnectivity[nN][0]:
topFacet.append(topConnectivity[nN][1])
else:
topFacet.append(topConnectivity[nN][0])
newFacets.append([topFacet])
newFacetFlags.append(5)
newVertices = []
newVertexFlags = []
for n,nF in zip(fineMesh.nodeArray, fineMesh.nodeMaterialTypes):
newVertices.append([n[0],n[1],n[2]])
if nF > 0:
newVertexFlags.append(nF)
else:
newVertexFlags.append(6)
for nN in range(len(newNodes.keys())):
newVertices.append(newNodes[nN+nN_start])
newVertexFlags.append(5)
In [13]:
fig = plt.figure()
newNodes = np.vstack((fineMesh.nodeArray,np.array(newNodes.values())))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(newNodes[:,0], newNodes[:,1], newNodes[:,2],c = newNodes[:,2],linewidths=0)
Out[13]:
In [14]:
from proteus.Domain import PiecewiseLinearComplexDomain
xmin_new = fineMesh.nodeArray[:,0].min()
xmax_new = fineMesh.nodeArray[:,0].max()
ymin_new = fineMesh.nodeArray[:,1].min()
ymax_new = fineMesh.nodeArray[:,1].max()
regions = [[0.5*(xmin+xmax),0.5*(ymin+ymax),zmax-1.0]]
regionFlags = [1.0]
domain3D = PiecewiseLinearComplexDomain(vertices=newVertices,
facets=newFacets,
facetHoles=None,
holes=None,
regions=regions,
vertexFlags=newVertexFlags,
facetFlags=newFacetFlags,
regionFlags=regionFlags,
regionConstraints=None,
name="frfDomain3D",
units="m")
In [15]:
domain3D.writePoly("frfDomain3D")
domain3D.writePLY("frfDomain3D")
In [18]:
!tetgen -KVApq1.33q12feen frfDomain3D.poly
In [17]:
fig = plt.figure()
polyNodes = np.array(domain3D.vertices)
ax = fig.add_subplot(111, projection='3d')
ax.scatter(newNodes[:,0], newNodes[:,1], newNodes[:,2],c = newNodes[:,2],linewidths=0)
ax.scatter(regions[0][0],regions[0][1],regions[0][2],linewidths=3)
Out[17]:
In [ ]: