Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Orphan Nodes #183

Merged
merged 3 commits into from
Apr 11, 2019
Merged

Orphan Nodes #183

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4ec9736
Merge branch 'master' of github.com:/stoiver/anuga_core
stoiver Apr 9, 2019
1c46050
Merge branch 'master' of github.com:/GeoscienceAustralia/anuga_core
stoiver Apr 9, 2019
5f7c6b9
vertex value bug when mesh has orphan nodes
stoiver Apr 11, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
50 changes: 29 additions & 21 deletions anuga/abstract_2d_finite_volumes/general_mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,15 +82,17 @@ def __init__(self,

"""

self.verbose = verbose

if verbose: log.critical('General_mesh: Building basic mesh structure')

self.use_inscribed_circle = use_inscribed_circle

self.triangles = num.array(triangles, num.int)

if verbose:
if verbose:
log.timingInfo("numTriangles, " + str(self.triangles.shape[0]))

self.nodes = num.array(nodes, num.float)

# Register number of elements and nodes
Expand Down Expand Up @@ -163,11 +165,11 @@ def __init__(self,
y2 = V2[:,1]

self.areas[:] = -((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2.0

#areas = -((x0-x1)*(y2-y1) - (y0-y1)*(x2-x1))/2.0

#assert num.allclose(self.areas, areas)

ind = num.where(self.areas <= 0.0)
msg = 'Degenerate Triangle(s) '+str(ind[0])
assert num.all(self.areas > 0.0), msg
Expand Down Expand Up @@ -225,7 +227,7 @@ def __init__(self,

self.normals[:,4] = yn2
self.normals[:,5] = -xn2

self.edgelengths[:,0] = l0
self.edgelengths[:,1] = l1
self.edgelengths[:,2] = l2
Expand Down Expand Up @@ -371,9 +373,9 @@ def __init__(self,
# Build structure listing which triangles belong to which node.
if verbose: log.critical('General Mesh: Building inverted triangle structure')
self.build_inverted_triangle_structure()

if verbose: log.timingInfo("aoi, '%s'" % self.get_area())


def __len__(self):

Expand All @@ -399,18 +401,18 @@ def get_normal(self, i, j):
"""

return self.normals[i, 2*j:2*j+2]

def get_edgelength(self, i, j):
"""Return length of j'th edge of the i'th triangle.
Return value is the numeric array slice [x, y]
"""
return self.edgelengths[i, j]


def get_number_of_triangles(self):
return self.number_of_triangles


def get_number_of_nodes(self):
return self.number_of_nodes

Expand Down Expand Up @@ -444,7 +446,7 @@ def get_node(self, i, absolute=False):
Default is False as many parts of ANUGA expects relative coordinates.
(To see which, switch to default absolute=True and run tests).

Note: This method returns a modified _copy_ of the nodes slice if
Note: This method returns a modified _copy_ of the nodes slice if
absolute is True. If absolute is False, just return the slice.
This is related to the ensure_numeric() returning a copy problem.
"""
Expand Down Expand Up @@ -490,8 +492,8 @@ def get_vertex_coordinates(self, triangle_id=None, absolute=False):
if absolute is True and not self.geo_reference.is_absolute():
offset=num.array([self.geo_reference.get_xllcorner(),
self.geo_reference.get_yllcorner()], num.float)
return V[i3:i3+3,:] + offset

return V[i3:i3+3,:] + offset
else:
return V[i3:i3+3,:]

Expand Down Expand Up @@ -557,7 +559,7 @@ def get_edge_midpoint_coordinates(self, triangle_id=None, absolute=False):
"""

E = self.edge_midpoint_coordinates

if triangle_id is None:
if absolute is True:
if not self.geo_reference.is_absolute():
Expand All @@ -574,7 +576,7 @@ def get_edge_midpoint_coordinates(self, triangle_id=None, absolute=False):
offset=num.array([self.geo_reference.get_xllcorner(),
self.geo_reference.get_yllcorner()], num.float)

return E[i3:i3+3,:] + offset
return E[i3:i3+3,:] + offset
else:
return E[i3:i3+3,:]

Expand All @@ -590,7 +592,7 @@ def get_edge_midpoint_coordinate(self, i, j, absolute=False):
E = self.get_edge_midpoint_coordinates(triangle_id=i, absolute=absolute)
return E[j,:] # Return (x, y) for edge mid point


def compute_edge_midpoint_coordinates(self):
"""Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
where the jth edge midpoint of the ith triangle is located in row 3*i+j.
Expand All @@ -610,7 +612,7 @@ def compute_edge_midpoint_coordinates(self):
V1 = V[1:3*M:3, :]
V2 = V[2:3*M:3, :]


#print V.shape, V0.shape, V1.shape, V2.shape

#print E.shape, E[0:3*M:3, :].shape, E[1:3*M:3, :].shape, E[2:3*M:3, :].shape
Expand Down Expand Up @@ -791,7 +793,14 @@ def build_inverted_triangle_structure(self):
number_of_triangles_per_node = num.bincount(self.triangles.flat).astype(num.int)
number_of_lone_nodes = self.number_of_nodes - len(number_of_triangles_per_node)

#print number_of_lone_nodes

orphan_nodes = num.argwhere(number_of_triangles_per_node==0)
number_of_orphan_nodes = len(orphan_nodes)

if number_of_orphan_nodes > 0 and self.verbose:
msg = 'Node(s) %d not associated to a triangle.' % orphan_nodes[0]
print msg

if number_of_lone_nodes > 0:
number_of_triangles_per_node = \
num.append(number_of_triangles_per_node,num.zeros(number_of_lone_nodes,num.int))
Expand All @@ -802,7 +811,7 @@ def build_inverted_triangle_structure(self):
number_of_entries = num.sum(number_of_triangles_per_node)

assert number_of_entries == 3*self.number_of_triangles

#vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#

# Array of vertex_indices (3*vol_id+vertex_id) sorted into contiguous
Expand Down Expand Up @@ -863,4 +872,3 @@ def set_georeference(self, g):

def get_georeference(self):
return self.geo_reference

12 changes: 8 additions & 4 deletions anuga/abstract_2d_finite_volumes/mesh_factory.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,9 @@ def rectangular_cross(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
len1 = float(len1)
len2 = float(len2)

m = int(m)
n = int(n)

params = []
params.append(m)
params.append(n)
Expand All @@ -158,7 +161,8 @@ def rectangular_cross(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):

arrParams = num.array(params)
arrOrigin = num.array(origin)
points = num.empty([(m+1)*(n+1)+m*n,2])

points = num.empty([(m+1)*(n+1)+m*n,2], dtype=num.float)
elements = num.empty([4*m*n,3], dtype=num.int)

from mesh_factory_ext import rectangular_cross_construct
Expand Down Expand Up @@ -326,7 +330,7 @@ def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)
processor = 0
numproc = 1


n = n_g
m_low = -1
m_high = m_g +1
Expand Down Expand Up @@ -414,7 +418,7 @@ def __call__(self, i,j):
boundary[nt, 2] = 'right'
else:
boundary[nt, 2] = 'ghost'

if j == 0:
boundary[nt, 1] = 'bottom'
elements[nt,:] = [i4,i3,i2]
Expand Down Expand Up @@ -447,7 +451,7 @@ def __call__(self, i,j):

Idfl = num.array(Idfl, num.int)
Idgr = num.array(Idgr, num.int)

full_send_dict[processor] = [Idfl, Idfl]
ghost_recv_dict[processor] = [Idgr, Idgr]

Expand Down