diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..609ef1425d0a0df37ace94b784f8ba346efe7b15
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,36 @@
+Python 2.7 AMD CPU:
+ script:
+ - export PY_EXE=python2.7
+ - export PYOPENCL_TEST=amd:pu
+ - export EXTRA_INSTALL="numpy mako"
+ - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh
+ - ". ./build-and-test-py-project.sh"
+ tags:
+ - python2.7
+ - amd-cl-cpu
+ except:
+ - tags
+Python 2.7 POCL:
+ script:
+ - export PY_EXE=python2.7
+ - export PYOPENCL_TEST=portable
+ - export EXTRA_INSTALL="numpy mako"
+ - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh
+ - ". ./build-and-test-py-project.sh"
+ tags:
+ - python2.7
+ - pocl
+ except:
+ - tags
+Python 3.4 POCL:
+ script:
+ - export PY_EXE=python3.4
+ - export PYOPENCL_TEST=portable
+ - export EXTRA_INSTALL="numpy mako"
+ - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh
+ - ". ./build-and-test-py-project.sh"
+ tags:
+ - python3.4
+ - pocl
+ except:
+ - tags
diff --git a/README.rst b/README.rst
index 929ac20b35e2b6103772839755b69b76f44f1e3f..3ca2d6c34f44e78a074468ac94b3547b3c603742 100644
--- a/README.rst
+++ b/README.rst
@@ -1,4 +1,7 @@
meshmode
========
+* `Source code on Github <https://github.com/inducer/meshmode>`_
+* `Documentation <https://documen.tician.de/meshmode>`_
+
.. TODO
diff --git a/doc/_static/akdoc.css b/doc/_static/akdoc.css
deleted file mode 100644
index d8b61e3ff7a358e5d5c0f132b5040ec7c3f43e42..0000000000000000000000000000000000000000
--- a/doc/_static/akdoc.css
+++ /dev/null
@@ -1,59 +0,0 @@
-pre {
- line-height: 110%;
-}
-
-.footer {
- background-color: #eee;
-}
-
-body > div.container {
- margin-top:10px;
-}
-
-dd {
- margin-left: 40px;
-}
-
-tt.descname {
- font-size: 100%;
-}
-
-code {
- color: rgb(51,51,51);
-}
-
-h1 {
- padding-bottom:7px;
- border-bottom: 1px solid #ccc;
-}
-
-h2 {
- padding-bottom:5px;
- border-bottom: 1px solid #ccc;
-}
-
-h3 {
- padding-bottom:5px;
- border-bottom: 1px solid #ccc;
-}
-
-.rubric {
- font-size: 120%;
- padding-bottom:1px;
- border-bottom: 1px solid #ccc;
-}
-
-.headerlink {
- padding-left: 1ex;
- padding-right: 1ex;
-}
-
-a.headerlink:hover {
- text-decoration: none;
-}
-
-blockquote p {
- font-size: 100%;
- font-weight: normal;
- line-height: normal;
-};
diff --git a/doc/_templates/layout.html b/doc/_templates/layout.html
deleted file mode 100644
index 400e7ec1d49677537aff6bf744e2803ef5c01e9e..0000000000000000000000000000000000000000
--- a/doc/_templates/layout.html
+++ /dev/null
@@ -1,2 +0,0 @@
-{% extends "!layout.html" %}
-{% set bootswatch_css_custom = ['_static/akdoc.css']%}
diff --git a/doc/conf.py b/doc/conf.py
index f43c6d3596d715838132b0dfc6bfbe6d29eaa386..1e5a0227f1671b4097dfd3b92417b3c7729ac195 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -103,26 +103,24 @@ pygments_style = 'sphinx'
# -- Options for HTML output ----------------------------------------------
-try:
- import sphinx_bootstrap_theme
-except:
- from warnings import warn
- warn("I would like to use the sphinx bootstrap theme, but can't find it.\n"
- "'pip install sphinx_bootstrap_theme' to fix.")
-else:
- # Activate the theme.
- html_theme = 'bootstrap'
- html_theme_path = sphinx_bootstrap_theme.get_html_theme_path()
-
- # Theme options are theme-specific and customize the look and feel of a theme
- # further. For a list of options available for each theme, see the
- # documentation.
- html_theme_options = {
- "navbar_fixed_top": "true",
- "navbar_site_name": "Contents",
- 'bootstrap_version': '3',
- 'source_link_position': 'footer',
+html_theme = "alabaster"
+
+html_theme_options = {
+ "extra_nav_links": {
+ "🚀 Github": "https://github.com/inducer/meshmode",
+ "💾 Download Releases": "https://pypi.python.org/pypi/meshmode",
}
+ }
+
+html_sidebars = {
+ '**': [
+ 'about.html',
+ 'navigation.html',
+ 'relations.html',
+ 'searchbox.html',
+ ]
+}
+
# Add any paths that contain custom themes here, relative to this directory.
#html_theme_path = []
diff --git a/examples/multiple-meshes.py b/examples/multiple-meshes.py
new file mode 100644
index 0000000000000000000000000000000000000000..84918e1779a34932ed6269a8884e30e63bb846d0
--- /dev/null
+++ b/examples/multiple-meshes.py
@@ -0,0 +1,26 @@
+from __future__ import division
+
+import numpy as np # noqa
+
+
+order = 4
+
+
+def main():
+ from meshmode.mesh.generation import ( # noqa
+ make_curve_mesh, starfish)
+ mesh1 = make_curve_mesh(starfish, np.linspace(0, 1, 20), 4)
+
+ from meshmode.mesh.processing import affine_map, merge_disjoint_meshes
+ mesh2 = affine_map(mesh1, b=np.array([2, 3]))
+
+ mesh = merge_disjoint_meshes((mesh1, mesh2))
+
+ from meshmode.mesh.visualization import draw_2d_mesh
+ draw_2d_mesh(mesh, set_bounding_box=True)
+
+ import matplotlib.pyplot as pt
+ pt.show()
+
+if __name__ == "__main__":
+ main()
diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py
index 09d795418a05f16c033f3b969862a1eda3ab07e0..068983ca81672295af65cdf3c82112478edac1a5 100644
--- a/examples/plot-connectivity.py
+++ b/examples/plot-connectivity.py
@@ -2,47 +2,12 @@ from __future__ import division
import numpy as np # noqa
import pyopencl as cl
-import random
-import os
-order = 4
-
-def main():
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
- from meshmode.mesh.generation import ( # noqa
- generate_icosphere, generate_icosahedron,
- generate_torus)
- #mesh = generate_icosphere(1, order=order)
- #mesh = generate_icosahedron(1, order=order)
- mesh = generate_torus(3, 1, order=order)
- from meshmode.mesh.refinement import Refiner
- r = Refiner(mesh)
- #mesh = r.refine(0)
- from meshmode.discretization import Discretization
- from meshmode.discretization.poly_element import \
- PolynomialWarpAndBlendGroupFactory
-
- discr = Discretization(
- cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
-
- from meshmode.discretization.visualization import make_visualizer
- vis = make_visualizer(queue, discr, order=1)
- os.remove("geometry.vtu")
- os.remove("connectivity.vtu")
- vis.write_vtk_file("geometry.vtu", [
- ("f", discr.nodes()[0]),
- ])
-
- from meshmode.discretization.visualization import \
- write_mesh_connectivity_vtk_file
-
- write_mesh_connectivity_vtk_file("connectivity.vtu",
- mesh)
+order = 4
-def main2():
+def main():
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
@@ -50,17 +15,8 @@ def main2():
generate_icosphere, generate_icosahedron,
generate_torus)
#mesh = generate_icosphere(1, order=order)
- #mesh = generate_icosahedron(1, order=order)
- mesh = generate_torus(3, 1, order=order)
- from meshmode.mesh.refinement import Refiner
- r = Refiner(mesh)
-
- times = random.randint(1, 1)
- for time in xrange(times):
- flags = np.zeros(len(mesh.groups[0].vertex_indices))
- for i in xrange(0, len(flags)):
- flags[i] = random.randint(0, 1)
- mesh = r.refine(flags)
+ mesh = generate_icosahedron(1, order=order)
+ #mesh = generate_torus(3, 1, order=order)
from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import \
@@ -71,18 +27,17 @@ def main2():
from meshmode.discretization.visualization import make_visualizer
vis = make_visualizer(queue, discr, order)
- os.remove("geometry2.vtu")
- os.remove("connectivity2.vtu")
- vis.write_vtk_file("geometry2.vtu", [
+
+ vis.write_vtk_file("geometry.vtu", [
("f", discr.nodes()[0]),
])
from meshmode.discretization.visualization import \
- write_mesh_connectivity_vtk_file
+ write_nodal_adjacency_vtk_file
- write_mesh_connectivity_vtk_file("connectivity2.vtu",
+ write_nodal_adjacency_vtk_file("adjacency.vtu",
mesh)
+
if __name__ == "__main__":
main()
- main2()
diff --git a/examples/refinement-playground.py b/examples/refinement-playground.py
new file mode 100644
index 0000000000000000000000000000000000000000..9927d6b53a647c043ee706e948d9b3de44f77141
--- /dev/null
+++ b/examples/refinement-playground.py
@@ -0,0 +1,497 @@
+from __future__ import division, print_function
+
+import numpy as np # noqa
+import pyopencl as cl
+import random
+import os
+order = 1
+
+
+def main():
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron,
+ generate_torus, generate_box_mesh)
+ #mesh = generate_icosphere(1, order=order)
+ #mesh = generate_icosahedron(1, order=order)
+ #mesh = generate_torus(3, 1, order=order)
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),))
+ from meshmode.mesh.refinement import Refiner
+ r = Refiner(mesh)
+ #mesh = r.refine(0)
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory
+
+ discr = Discretization(cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, order)
+ os.remove("geometry.vtu")
+ os.remove("connectivity.vtu")
+ vis.write_vtk_file("geometry.vtu", [
+ ("f", discr.nodes()[0]),
+ ])
+
+ from meshmode.discretization.visualization import \
+ write_mesh_connectivity_vtk_file
+
+ write_mesh_connectivity_vtk_file("connectivity.vtu",
+ mesh)
+
+
+#construct vertex vertex_index
+def get_vertex(mesh, vertex_index):
+ vertex = np.empty([len(mesh.vertices)])
+ for i_cur_dim, cur_dim_coords in enumerate(mesh.vertices):
+ vertex[i_cur_dim] = cur_dim_coords[vertex_index]
+ return vertex
+
+
+def test_refiner_connectivity(mesh):
+ def group_and_iel_to_global_iel(igrp, iel):
+ return mesh.groups[igrp].element_nr_base + iel
+
+ from meshmode.mesh.tools import make_element_lookup_tree
+ tree = make_element_lookup_tree(mesh)
+
+ from meshmode.mesh.processing import find_bounding_box
+ bbox_min, bbox_max = find_bounding_box(mesh)
+
+ extent = bbox_max-bbox_min
+ cnx = mesh.element_connectivity
+ nvertices_per_element = len(mesh.groups[0].vertex_indices[0])
+ #stores connectivity obtained from geometry using barycentric coordinates
+ connected_to_element_geometry = np.zeros(
+ (mesh.nelements, mesh.nelements), dtype=bool)
+ #store connectivity obtained from mesh's connectivity
+ connected_to_element_connectivity = np.zeros(
+ (mesh.nelements, mesh.nelements), dtype=bool)
+
+ import six
+ for igrp, grp in enumerate(mesh.groups):
+ iel_base = grp.element_nr_base
+ for iel_grp in six.moves.range(grp.nelements):
+
+ iel_g = group_and_iel_to_global_iel(igrp, iel_grp)
+ nb_starts = cnx.neighbors_starts
+ for nb_iel_g in cnx.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]:
+ connected_to_element_connectivity[nb_iel_g, iel_g] = True
+ connected_to_element_connectivity[iel_g, nb_iel_g] = True
+
+ for vertex_index in grp.vertex_indices[iel_grp]:
+ vertex = get_vertex(mesh, vertex_index)
+ #check which elements touch this vertex
+ for bounding_igrp, bounding_iel in tree.generate_matches(vertex):
+ if bounding_igrp == igrp and bounding_iel == iel_grp:
+ continue
+ bounding_grp = mesh.groups[bounding_igrp]
+
+ last_bounding_vertex = get_vertex(mesh,
+ bounding_grp.vertex_indices[bounding_iel][nvertices_per_element-1])
+ transformation = np.empty([len(mesh.vertices), nvertices_per_element-1])
+ vertex_transformed = vertex - last_bounding_vertex
+ for ibounding_vertex_index, bounding_vertex_index in enumerate(
+ bounding_grp.vertex_indices[bounding_iel][:nvertices_per_element-1]):
+ bounding_vertex = get_vertex(mesh, bounding_vertex_index)
+ transformation[:,ibounding_vertex_index] = bounding_vertex - last_bounding_vertex
+ barycentric_coordinates = np.linalg.solve(transformation, vertex_transformed)
+ is_connected = True
+ sum_of_coords = 0.0
+ for coord in barycentric_coordinates:
+ if coord < 0:
+ is_connected = False
+ sum_of_coords += coord
+ if sum_of_coords > 1:
+ is_connected = False
+ if is_connected:
+ connected_to_element_geometry[group_and_iel_to_global_iel(bounding_igrp, bounding_iel),
+ group_and_iel_to_global_iel(igrp, iel_grp)] = True
+ connected_to_element_geometry[group_and_iel_to_global_iel(igrp, iel_grp),
+ group_and_iel_to_global_iel(bounding_igrp, bounding_iel)] = True
+ print("GEOMETRY: ")
+ print(connected_to_element_geometry)
+ print("CONNECTIVITY: ")
+ print(connected_to_element_connectivity)
+ cmpmatrix = (connected_to_element_geometry == connected_to_element_connectivity)
+ print(cmpmatrix)
+ print(np.where(cmpmatrix == False)) # noqa
+ if not (connected_to_element_geometry == connected_to_element_connectivity).all():
+ cmpmatrix = connected_to_element_connectivity == connected_to_element_geometry
+ for ii, i in enumerate(cmpmatrix):
+ for ij, j in enumerate(i):
+ if not j:
+ print(
+ ii, ij,
+ "GEOMETRY: ",
+ connected_to_element_geometry[ii][ij],
+ "CONNECTIVITY: ",
+ connected_to_element_connectivity[ii][ij])
+ assert (connected_to_element_geometry == connected_to_element_connectivity).all()
+
+
+def test_refiner_connectivity_efficient(mesh):
+ def group_and_iel_to_global_iel(igrp, iel):
+ return mesh.groups[igrp].element_nr_base + iel
+
+ from meshmode.mesh.tools import make_element_lookup_tree
+ tree = make_element_lookup_tree(mesh)
+
+ from meshmode.mesh.processing import find_bounding_box
+ bbox_min, bbox_max = find_bounding_box(mesh)
+
+ extent = bbox_max-bbox_min
+ cnx = mesh.element_connectivity
+ nvertices_per_element = len(mesh.groups[0].vertex_indices[0])
+ #stores connectivity obtained from geometry using barycentric coordinates
+ #connected_to_element_geometry = np.zeros([mesh.nelements, mesh.nelements], dtype=bool)
+ #store connectivity obtained from mesh's connectivity
+ #connected_to_element_connectivity = np.zeros([mesh.nelements, mesh.nelements], dtype=bool)
+ connected_to_element_geometry = []
+ connected_to_element_connectivity = []
+ import six
+ for i in six.moves.range(mesh.nelements):
+ connected_to_element_geometry.append([])
+ connected_to_element_connectivity.append([])
+
+ for igrp, grp in enumerate(mesh.groups):
+ iel_base = grp.element_nr_base
+ for iel_grp in six.moves.range(grp.nelements):
+
+ iel_g = group_and_iel_to_global_iel(igrp, iel_grp)
+ nb_starts = cnx.neighbors_starts
+ for nb_iel_g in cnx.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]:
+ if iel_g not in connected_to_element_connectivity[nb_iel_g]:
+ connected_to_element_connectivity[nb_iel_g].append(iel_g)
+ if nb_iel_g not in connected_to_element_connectivity[iel_g]:
+ connected_to_element_connectivity[iel_g].append(nb_iel_g)
+ #connected_to_element_connectivity[nb_iel_g, iel_g] = True
+ #connected_to_element_connectivity[iel_g, nb_iel_g] = True
+
+ for vertex_index in grp.vertex_indices[iel_grp]:
+ vertex = get_vertex(mesh, vertex_index)
+ #check which elements touch this vertex
+ for bounding_igrp, bounding_iel in tree.generate_matches(vertex):
+ if bounding_igrp == igrp and bounding_iel == iel_grp:
+ continue
+ bounding_grp = mesh.groups[bounding_igrp]
+
+ last_bounding_vertex = get_vertex(mesh, \
+ bounding_grp.vertex_indices[bounding_iel][nvertices_per_element-1])
+ transformation = np.empty([len(mesh.vertices), nvertices_per_element-1])
+ vertex_transformed = vertex - last_bounding_vertex
+ for ibounding_vertex_index, bounding_vertex_index in \
+ enumerate(bounding_grp.vertex_indices[bounding_iel][:nvertices_per_element-1]):
+ bounding_vertex = get_vertex(mesh, bounding_vertex_index)
+ transformation[:,ibounding_vertex_index] = bounding_vertex - last_bounding_vertex
+ barycentric_coordinates = np.linalg.solve(transformation, vertex_transformed)
+ is_connected = True
+ sum_of_coords = 0.0
+ for coord in barycentric_coordinates:
+ if coord < 0:
+ is_connected = False
+ sum_of_coords += coord
+ if sum_of_coords > 1:
+ is_connected = False
+ if is_connected:
+ el1 = group_and_iel_to_global_iel(bounding_igrp, bounding_iel)
+ el2 = group_and_iel_to_global_iel(igrp, iel_grp)
+ if el2 not in connected_to_element_geometry[el1]:
+ connected_to_element_geometry[el1].append(el2)
+ if el1 not in connected_to_element_geometry[el2]:
+ connected_to_element_geometry[el2].append(el1)
+
+ from collections import Counter
+ for i in six.moves.range(mesh.nelements):
+ assert Counter(connected_to_element_connectivity[i]) \
+ == Counter(connected_to_element_geometry[i])
+
+
+def linear_func(vert):
+ csum = 0
+ for i in vert:
+ csum += i
+ #print csum
+ return csum
+
+
+def sine_func(vert):
+ #print vert
+ import math
+ res = 1
+ for i in vert:
+ res *= math.sin(i*2.0*math.pi)
+ #print res
+ return abs(res) * 0.2 + 0.2
+
+
+def quadratic_func(vert):
+ csum = 0
+ for i in vert:
+ csum += i * i
+ return csum * 1.5
+
+
+def get_function_flags(mesh, function):
+ import six
+ from math import sqrt
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ for grp in mesh.groups:
+ iel_base = grp.element_nr_base
+ for iel_grp in six.moves.range(grp.nelements):
+ vertex_indices = grp.vertex_indices[iel_grp]
+ max_edge_len = 0
+ for i in six.moves.range(len(vertex_indices)):
+ for j in six.moves.range(i+1, len(vertex_indices)):
+ edge_len = 0
+ for k in six.moves.range(len(mesh.vertices)):
+ edge_len += (
+ (mesh.vertices[k, vertex_indices[i]]
+ - mesh.vertices[k, vertex_indices[j]])
+ * (mesh.vertices[k, vertex_indices[i]]
+ - mesh.vertices[k, vertex_indices[j]]))
+ edge_len = sqrt(edge_len)
+ max_edge_len = max(max_edge_len, edge_len)
+ #print(edge_lens[0], mesh.vertices[0, vertex_indices[i]], mesh.vertices[1, vertex_indices[i]], mesh.vertices[2, vertex_indices[i]]) # noqa
+ centroid = [0] * len(mesh.vertices)
+ for j in six.moves.range(len(mesh.vertices)):
+ centroid[j] += mesh.vertices[j, vertex_indices[i]]
+ for i in six.moves.range(len(mesh.vertices)):
+ centroid[i] /= len(vertex_indices)
+ yes = False
+ val = function(centroid)
+ if max_edge_len > val:
+ flags[iel_grp] = True
+ return flags
+
+
+def get_corner_flags(mesh):
+ import six
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ for grp in mesh.groups:
+ iel_base = grp.element_nr_base
+ for iel_grp in six.moves.range(grp.nelements):
+ is_corner_el = False
+ vertex_indices = grp.vertex_indices[iel_grp]
+ for i in six.moves.range(len(vertex_indices)):
+ cur_vertex_corner = True
+ for j in six.moves.range(len(mesh.vertices)):
+ print(iel_grp, i, mesh.vertices[j, vertex_indices[i]])
+ if mesh.vertices[j, vertex_indices[i]] != 0.0:
+ cur_vertex_corner = False
+ if cur_vertex_corner:
+ is_corner_el = True
+ break
+ if is_corner_el:
+ print(iel_grp)
+ flags[iel_grp] = True
+ return flags
+
+
+def get_random_flags(mesh):
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ for i in xrange(0, len(flags)):
+ flags[i] = random.randint(0, 1)
+ return flags
+
+
+def refine_and_generate_chart_function(mesh, filename, function):
+ import six
+ from time import clock
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+ print("NELEMENTS: ", mesh.nelements)
+ #print mesh
+ for i in six.moves.range(len(mesh.groups[0].vertex_indices[0])):
+ for k in six.moves.range(len(mesh.vertices)):
+ print(mesh.vertices[k, i])
+
+ #test_refiner_connectivity(mesh);
+ from meshmode.mesh.refinement import Refiner
+ r = Refiner(mesh)
+ #random.seed(0)
+ poss_flags = np.ones(len(mesh.groups[0].vertex_indices))
+ #times = 3
+ num_elements = []
+ time_t = []
+ #nelements = mesh.nelements
+ while True:
+ print("NELS:", mesh.nelements)
+ #flags = get_corner_flags(mesh)
+ flags = get_function_flags(mesh, function)
+ nels = 0
+ for i in flags:
+ if i:
+ nels += 1
+ if nels == 0:
+ break
+ print("LKJASLFKJALKASF:", nels)
+ num_elements.append(nels)
+ #flags = get_corner_flags(mesh)
+ beg = clock()
+ mesh = r.refine(flags)
+ end = clock()
+ time_taken = end - beg
+ time_t.append(time_taken)
+ #if nelements == mesh.nelements:
+ #break
+ #nelements = mesh.nelements
+ #from meshmode.mesh.visualization import draw_2d_mesh
+ #draw_2d_mesh(mesh, True, True, True, fill=None)
+ #import matplotlib.pyplot as pt
+ #pt.show()
+
+ #poss_flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ #for i in xrange(0, len(flags)):
+ # poss_flags[i] = flags[i]
+ #for i in xrange(len(flags), len(poss_flags)):
+ # poss_flags[i] = 1
+
+ import matplotlib.pyplot as plt
+ import matplotlib.pyplot as pt
+ pt.xlabel('Number of elements being refined')
+ pt.ylabel('Time taken')
+ pt.plot(num_elements, time_t, "o")
+ pt.savefig(filename, format='pdf')
+ pt.clf()
+ print('DONE REFINING')
+ '''
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ flags[0] = 1
+ flags[1] = 1
+ mesh = r.refine(flags)
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ flags[0] = 1
+ flags[1] = 1
+ flags[2] = 1
+ mesh = r.refine(flags)
+ '''
+ #test_refiner_connectivity_efficient(mesh)
+ #r.print_rays(70)
+ #r.print_rays(117)
+ #r.print_hanging_elements(10)
+ #r.print_hanging_elements(117)
+ #r.print_hanging_elements(757)
+ #from meshmode.mesh.visualization import draw_2d_mesh
+ #draw_2d_mesh(mesh, False, False, False, fill=None)
+ #import matplotlib.pyplot as pt
+ #pt.show()
+
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory
+ discr = Discretization(
+ cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, order)
+ os.remove("connectivity2.vtu")
+ os.remove("geometry2.vtu")
+ vis.write_vtk_file("geometry2.vtu", [
+ ("f", discr.nodes()[0]),
+ ])
+
+ from meshmode.discretization.visualization import \
+ write_mesh_connectivity_vtk_file
+
+ write_mesh_connectivity_vtk_file("connectivity2.vtu",
+ mesh)
+
+
+def main2():
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron,
+ generate_torus, generate_regular_rect_mesh,
+ generate_box_mesh)
+# mesh = generate_icosphere(1, order=order)
+ #mesh = generate_icosahedron(1, order=order)
+ #mesh = generate_torus(3, 1, order=order)
+ #mesh = generate_regular_rect_mesh()
+ #mesh = generate_box_mesh(3*(np.linspace(0, 3, 5),))
+ #mesh = generate_box_mesh(3*(np.linspace(0, 1, 3),))
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),))
+ refine_and_generate_chart_function(mesh, "plot.pdf", sine_func)
+
+
+def all_refine(num_mesh, depth, fname):
+ import six
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron,
+ generate_torus, generate_regular_rect_mesh,
+ generate_box_mesh)
+ import timeit
+ nelements = []
+ runtimes = []
+ for el_fact in six.moves.range(2, num_mesh+2):
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
+ from meshmode.mesh.refinement import Refiner
+ r = Refiner(mesh)
+ for time in xrange(depth):
+ flags = np.ones(len(mesh.groups[0].vertex_indices))
+ if time < depth-1:
+ mesh = r.refine(flags)
+ else:
+ start = timeit.default_timer()
+ mesh = r.refine(flags)
+ stop = timeit.default_timer()
+ nelements.append(mesh.nelements)
+ runtimes.append(stop-start)
+ #test_refiner_connectivity_efficient(mesh)
+ import matplotlib.pyplot as pt
+ pt.plot(nelements, runtimes, "o")
+ pt.savefig(fname)
+ pt.clf()
+ #pt.show()
+
+
+def uniform_refine(num_mesh, fract, depth, fname):
+ import six
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron,
+ generate_torus, generate_regular_rect_mesh,
+ generate_box_mesh)
+ import timeit
+ nelements = []
+ runtimes = []
+ for el_fact in six.moves.range(2, num_mesh+2):
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
+ from meshmode.mesh.refinement import Refiner
+ r = Refiner(mesh)
+ all_els = range(mesh.nelements)
+ for time in xrange(depth):
+ flags = np.zeros(mesh.nelements)
+ from random import shuffle
+ shuffle(all_els)
+ nels_this_round = 0
+ for i in six.moves.range(len(all_els)):
+ if i / len(flags) > fract:
+ break
+ flags[all_els[i]] = 1
+ nels_this_round += 1
+
+ if time < depth-1:
+ mesh = r.refine(flags)
+ else:
+ start = timeit.default_timer()
+ mesh = r.refine(flags)
+ stop = timeit.default_timer()
+ nelements.append(mesh.nelements)
+ runtimes.append(stop-start)
+ all_els = []
+ for i in six.moves.range(len(flags)):
+ if flags[i]:
+ all_els.append(i)
+ for i in six.moves.range(len(flags), mesh.nelements):
+ all_els.append(i)
+ test_refiner_connectivity(mesh)
+ import matplotlib.pyplot as pt
+ pt.plot(nelements, runtimes, "o")
+ pt.savefig(fname)
+ pt.clf()
+
+
+if __name__ == "__main__":
+ print("HEREERERE")
+ #all_refine(3, 2, 'all_a.pdf')
+ #all_refine(3, 3, 'all_b.pdf')
+ #uniform_refine(3, 0.2, 3, 'uniform_a.pdf')
+ main2()
diff --git a/meshmode/__init__.py b/meshmode/__init__.py
index e1a4c5e6aedbffd76ceed3f426ee347380d13cce..1c38d87894841834007489ac2a4eea5175ba1428 100644
--- a/meshmode/__init__.py
+++ b/meshmode/__init__.py
@@ -25,7 +25,7 @@ THE SOFTWARE.
__doc__ = """
.. exception:: Error
-.. exception:: ConnectivityUnavailable
+.. exception:: DataUnavailable
"""
@@ -33,5 +33,5 @@ class Error(RuntimeError):
pass
-class ConnectivityUnavailable(Error):
+class DataUnavailable(Error):
pass
diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index 6088753222fcdfde70e51b0cca237f2af2dd9f78..e1ed150f132a516d409757567d4b5219db669235 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -23,9 +23,10 @@ THE SOFTWARE.
"""
import numpy as np
-from pytools import memoize_method, memoize_method_nested
+from pytools import memoize_method, memoize_in
import loopy as lp
import pyopencl as cl
+import pyopencl.array # noqa
__doc__ = """
.. autoclass:: ElementGroupBase
@@ -37,11 +38,50 @@ __doc__ = """
# {{{ element group base
class ElementGroupBase(object):
- """Container for the :class:`QBXGrid` data corresponding to
+ """Container for the :class:`Discretization` data corresponding to
one :class:`meshmode.mesh.MeshElementGroup`.
.. attribute :: mesh_el_group
+ .. attribute :: order
.. attribute :: node_nr_base
+
+ .. autoattribute:: nelements
+ .. autoattribute:: nunit_nodes
+ .. autoattribute:: nnodes
+ .. autoattribute:: dim
+ .. automethod:: view
+
+ .. method:: unit_nodes()
+
+ Returns a :class:`numpy.ndarray` of shape ``(dim, nunit_nodes)``
+ of reference coordinates of interpolation nodes.
+
+ .. method:: basis()
+
+ Returns a :class:`list` of basis functions that take arrays
+ of shape ``(dim, n)`` and return an array of shape (n,)``
+ (which performs evaluation of the basis function).
+
+ .. method:: grad_basis()
+
+ :returns: a :class:`tuple` of functions, each of which
+ accepts arrays of shape *(dims, npts)*
+ and returns a :class:`tuple` of length *dims* containing
+ the derivatives along each axis as an array of size *npts*.
+ 'Scalar' evaluation, by passing just one vector of length *dims*,
+ is also supported.
+
+ .. method:: diff_matrices()
+
+ Return a :attr:`dim`-long :class:`tuple` of matrices of
+ shape ``(nunit_nodes, nunit_nodes)``, each of which,
+ when applied to an array of nodal values, take derivatives
+ in the reference (r,s,t) directions.
+
+ .. method:: weights()
+
+ Returns an array of length :attr:`nunit_nodes` containing
+ quadrature weights.
"""
def __init__(self, mesh_el_group, order, node_nr_base):
@@ -79,6 +119,11 @@ class ElementGroupBase(object):
(-1, -1))
def view(self, global_array):
+ """Return a view of *global_array* of shape ``(..., nelements,
+ nunit_nodes)`` where *global_array* is of shape ``(..., nnodes)``,
+ where *nnodes* is the global (per-discretization) node count.
+ """
+
return global_array[
..., self.node_nr_base:self.node_nr_base + self.nnodes] \
.reshape(
@@ -127,7 +172,9 @@ class Discretization(object):
.. attribute :: groups
- .. method:: empty(dtype, queue=None, extra_dims=None)
+ .. automethod:: empty
+
+ .. automethod:: zeros
.. method:: nodes()
@@ -141,11 +188,6 @@ class Discretization(object):
"""
def __init__(self, cl_ctx, mesh, group_factory, real_dtype=np.float64):
- """
- :arg order: A polynomial-order-like parameter passed unmodified to
- :attr:`group_class`. See subclasses for more precise definition.
- """
-
self.cl_context = cl_ctx
self.mesh = mesh
@@ -170,7 +212,20 @@ class Discretization(object):
def ambient_dim(self):
return self.mesh.ambient_dim
- def empty(self, dtype, queue=None, extra_dims=None):
+ def empty(self, queue=None, dtype=None, extra_dims=None, allocator=None):
+ """Return an empty DOF vector.
+
+ :arg dtype: type special value 'c' will result in a
+ vector of dtype :attr:`self.complex_dtype`. If
+ *None* (the default), a real vector will be returned.
+ """
+ if dtype is None:
+ dtype = self.real_dtype
+ elif dtype == "c":
+ dtype = self.complex_dtype
+ else:
+ dtype = np.dtype(dtype)
+
if queue is None:
first_arg = self.cl_context
else:
@@ -180,11 +235,14 @@ class Discretization(object):
if extra_dims is not None:
shape = extra_dims + shape
- return cl.array.empty(first_arg, shape, dtype=dtype)
+ return cl.array.empty(first_arg, shape, dtype=dtype, allocator=allocator)
+
+ def zeros(self, queue, dtype=None, extra_dims=None, allocator=None):
+ return self.empty(queue, dtype=dtype, extra_dims=extra_dims,
+ allocator=allocator).fill(0)
- def num_reference_derivative(
- self, queue, ref_axes, vec):
- @memoize_method_nested
+ def num_reference_derivative(self, queue, ref_axes, vec):
+ @memoize_in(self, "reference_derivative_knl")
def knl():
knl = lp.make_kernel(
"""{[k,i,j]:
@@ -196,9 +254,12 @@ class Discretization(object):
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
- result = self.empty(vec.dtype)
+ result = self.empty(dtype=vec.dtype)
for grp in self.groups:
+ if grp.nelements == 0:
+ continue
+
mat = None
for ref_axis in ref_axes:
next_mat = grp.diff_matrices()[ref_axis]
@@ -212,24 +273,28 @@ class Discretization(object):
return result
def quad_weights(self, queue):
- @memoize_method_nested
+ @memoize_in(self, "quad_weights_knl")
def knl():
knl = lp.make_kernel(
"{[k,i]: 0<=k<nelements and 0<=i<ndiscr_nodes}",
"result[k,i] = weights[i]",
- name="quad_weights")
+ name="quad_weights",
+ default_offset=lp.auto)
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
- result = self.empty(self.real_dtype)
+ result = self.empty(dtype=self.real_dtype)
for grp in self.groups:
+ if grp.nelements == 0:
+ continue
+
knl()(queue, result=grp.view(result), weights=grp.weights)
return result
@memoize_method
def nodes(self):
- @memoize_method_nested
+ @memoize_in(self, "nodes_knl")
def knl():
knl = lp.make_kernel(
"""{[d,k,i,j]:
@@ -250,13 +315,16 @@ class Discretization(object):
"stride:auto,stride:auto,stride:auto")
return knl
- result = self.empty(self.real_dtype, extra_dims=(self.ambient_dim,))
+ result = self.empty(dtype=self.real_dtype, extra_dims=(self.ambient_dim,))
with cl.CommandQueue(self.cl_context) as queue:
for grp in self.groups:
+ if grp.nelements == 0:
+ continue
+
meg = grp.mesh_el_group
knl()(queue,
- resampling_mat=grp.resampling_matrix(),
+ resampling_mat=grp.from_mesh_interp_matrix(),
result=grp.view(result), nodes=meg.nodes)
return result
diff --git a/meshmode/discretization/connection.py b/meshmode/discretization/connection.py
deleted file mode 100644
index 3f1614170e1deddb8e2aeca626c590114cae34c5..0000000000000000000000000000000000000000
--- a/meshmode/discretization/connection.py
+++ /dev/null
@@ -1,435 +0,0 @@
-from __future__ import division
-from __future__ import absolute_import
-import six
-from six.moves import range
-from six.moves import zip
-
-__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
-
-__license__ = """
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in
-all copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-THE SOFTWARE.
-"""
-
-import numpy as np
-import modepy as mp
-import pyopencl as cl
-import pyopencl.array # noqa
-from pytools import memoize_method, memoize_method_nested, Record
-
-import logging
-logger = logging.getLogger(__name__)
-
-
-__doc__ = """
-.. autoclass:: DiscretizationConnection
-
-.. autofunction:: make_same_mesh_connection
-
-.. autofunction:: make_boundary_restriction
-
-Implementation details
-^^^^^^^^^^^^^^^^^^^^^^
-
-.. autoclass:: InterpolationBatch
-
-.. autoclass:: DiscretizationConnectionElementGroup
-"""
-
-
-class InterpolationBatch(object):
- """One interpolation batch captures how a batch of elements *within* an
- element group should be an interpolated. Note that while it's possible that
- an interpolation batch takes care of interpolating an entire element group
- from source to target, that's not *necessarily* the case. Consider the case
- of extracting boundary values of a discretization. For, say, a triangle, at
- least three different interpolation batches are needed to cover boundary
- edges that fall onto each of the three edges of the unit triangle.
-
- .. attribute:: source_element_indices
-
- A :class:`pyopencl.array.Arrays` of length ``nelements``, containing the
- element index from which this "*to*" element's data will be
- interpolated.
-
- .. attribute:: target_element_indices
-
- A :class:`pyopencl.array.Arrays` of length ``nelements``, containing the
- element index to which this "*to*" element's data will be
- interpolated.
-
- .. attribute:: result_unit_nodes
-
- A :class:`numpy.ndarray` of shape
- ``(from_group.dim,to_group.nelements,to_group.nunit_nodes)``
- storing the coordinates of the nodes (in unit coordinates
- of the *from* reference element) from which the node
- locations of this element should be interpolated.
- """
-
- def __init__(self, source_element_indices,
- target_element_indices, result_unit_nodes):
- self.source_element_indices = source_element_indices
- self.target_element_indices = target_element_indices
- self.result_unit_nodes = result_unit_nodes
-
- @property
- def nelements(self):
- return len(self.source_element_indices)
-
-
-class DiscretizationConnectionElementGroup(object):
- """
- .. attribute:: batches
-
- A list of :class:`InterpolationBatch` instances.
- """
- def __init__(self, batches):
- self.batches = batches
-
-
-class DiscretizationConnection(object):
- """
- .. attribute:: from_discr
-
- .. attribute:: to_discr
-
- .. attribute:: groups
-
- a list of :class:`MeshConnectionGroup` instances, with
- a one-to-one correspondence to the groups in
- :attr:`from_discr` and :attr:`to_discr`.
- """
-
- def __init__(self, from_discr, to_discr, groups):
- if from_discr.cl_context != to_discr.cl_context:
- raise ValueError("from_discr and to_discr must live in the "
- "same OpenCL context")
-
- self.cl_context = from_discr.cl_context
-
- self.from_discr = from_discr
- self.to_discr = to_discr
- self.groups = groups
-
- @memoize_method
- def _resample_matrix(self, elgroup_index, ibatch_index):
- import modepy as mp
- ibatch = self.groups[elgroup_index].batches[ibatch_index]
- from_grp = self.from_discr.groups[elgroup_index]
-
- return mp.resampling_matrix(
- mp.simplex_onb(self.from_discr.dim, from_grp.order),
- ibatch.result_unit_nodes, from_grp.unit_nodes)
-
- def __call__(self, queue, vec):
- @memoize_method_nested
- def knl():
- import loopy as lp
- knl = lp.make_kernel(
- """{[k,i,j]:
- 0<=k<nelements and
- 0<=i<n_to_nodes and
- 0<=j<n_from_nodes}""",
- "result[target_element_indices[k], i] \
- = sum(j, resample_mat[i, j] \
- * vec[source_element_indices[k], j])",
- [
- lp.GlobalArg("result", None,
- shape="nelements_result, n_to_nodes"),
- lp.GlobalArg("vec", None,
- shape="nelements_vec, n_from_nodes"),
- lp.ValueArg("nelements_result", np.int32),
- lp.ValueArg("nelements_vec", np.int32),
- "...",
- ],
- name="oversample")
-
- knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
- return lp.tag_inames(knl, dict(k="g.0"))
-
- if not isinstance(vec, cl.array.Array):
- return vec
-
- result = self.to_discr.empty(vec.dtype)
-
- if vec.shape != (self.from_discr.nnodes,):
- raise ValueError("invalid shape of incoming resampling data")
-
- for i_grp, (sgrp, tgrp, cgrp) in enumerate(
- zip(self.to_discr.groups, self.from_discr.groups, self.groups)):
- for i_batch, batch in enumerate(cgrp.batches):
- if len(batch.source_element_indices):
- knl()(queue,
- resample_mat=self._resample_matrix(i_grp, i_batch),
- result=sgrp.view(result), vec=tgrp.view(vec),
- source_element_indices=batch.source_element_indices,
- target_element_indices=batch.target_element_indices)
-
- return result
-
- # }}}
-
-
-# {{{ same-mesh constructor
-
-def make_same_mesh_connection(queue, to_discr, from_discr):
- if from_discr.mesh is not to_discr.mesh:
- raise ValueError("from_discr and to_discr must be based on "
- "the same mesh")
-
- assert queue.context == from_discr.cl_context
- assert queue.context == to_discr.cl_context
-
- groups = []
- for fgrp, tgrp in zip(from_discr.groups, to_discr.groups):
- all_elements = cl.array.arange(queue,
- fgrp.nelements,
- dtype=np.intp).with_queue(None)
- ibatch = InterpolationBatch(
- source_element_indices=all_elements,
- target_element_indices=all_elements,
- result_unit_nodes=tgrp.unit_nodes)
-
- groups.append(
- DiscretizationConnectionElementGroup([ibatch]))
-
- return DiscretizationConnection(
- from_discr, to_discr, groups)
-
-# }}}
-
-
-# {{{ boundary restriction constructor
-
-class _ConnectionBatchData(Record):
- pass
-
-
-def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data):
- connection_groups = []
- for igrp, (vol_grp, bdry_grp) in enumerate(
- zip(vol_discr.groups, bdry_discr.groups)):
- connection_batches = []
- mgrp = vol_grp.mesh_el_group
-
- for face_id in range(len(mgrp.face_vertex_indices())):
- data = connection_data[igrp, face_id]
-
- bdry_unit_nodes_01 = (bdry_grp.unit_nodes + 1)*0.5
- result_unit_nodes = (np.dot(data.A, bdry_unit_nodes_01).T + data.b).T
-
- connection_batches.append(
- InterpolationBatch(
- source_element_indices=cl.array.to_device(
- queue,
- vol_grp.mesh_el_group.element_nr_base
- + data.group_source_element_indices)
- .with_queue(None),
- target_element_indices=cl.array.to_device(
- queue,
- bdry_grp.mesh_el_group.element_nr_base
- + data.group_target_element_indices)
- .with_queue(None),
- result_unit_nodes=result_unit_nodes,
- ))
-
- connection_groups.append(
- DiscretizationConnectionElementGroup(
- connection_batches))
-
- return DiscretizationConnection(
- vol_discr, bdry_discr, connection_groups)
-
-
-def make_boundary_restriction(queue, discr, group_factory):
- """
- :return: a tuple ``(bdry_mesh, bdry_discr, connection)``
- """
-
- logger.info("building boundary connection: start")
-
- # {{{ build face_map
-
- # maps (igrp, el_grp, face_id) to a frozenset of vertex IDs
- face_map = {}
-
- for igrp, mgrp in enumerate(discr.mesh.groups):
- grp_face_vertex_indices = mgrp.face_vertex_indices()
-
- for iel_grp in range(mgrp.nelements):
- for fid, loc_face_vertices in enumerate(grp_face_vertex_indices):
- face_vertices = frozenset(
- mgrp.vertex_indices[iel_grp, fvi]
- for fvi in loc_face_vertices
- )
- face_map.setdefault(face_vertices, []).append(
- (igrp, iel_grp, fid))
-
- del face_vertices
-
- # }}}
-
- boundary_faces = [
- face_ids[0]
- for face_vertices, face_ids in six.iteritems(face_map)
- if len(face_ids) == 1]
-
- from pytools import flatten
- bdry_vertex_vol_nrs = sorted(set(flatten(six.iterkeys(face_map))))
-
- vol_to_bdry_vertices = np.empty(
- discr.mesh.vertices.shape[-1],
- discr.mesh.vertices.dtype)
- vol_to_bdry_vertices.fill(-1)
- vol_to_bdry_vertices[bdry_vertex_vol_nrs] = np.arange(
- len(bdry_vertex_vol_nrs))
-
- bdry_vertices = discr.mesh.vertices[:, bdry_vertex_vol_nrs]
-
- from meshmode.mesh import Mesh, SimplexElementGroup
- bdry_mesh_groups = []
- connection_data = {}
-
- for igrp, grp in enumerate(discr.groups):
- mgrp = grp.mesh_el_group
- group_boundary_faces = [
- (ibface_el, ibface_face)
- for ibface_group, ibface_el, ibface_face in boundary_faces
- if ibface_group == igrp]
-
- if not isinstance(mgrp, SimplexElementGroup):
- raise NotImplementedError("can only take boundary of "
- "SimplexElementGroup-based meshes")
-
- # {{{ Preallocate arrays for mesh group
-
- ngroup_bdry_elements = len(group_boundary_faces)
- vertex_indices = np.empty(
- (ngroup_bdry_elements, mgrp.dim+1-1),
- mgrp.vertex_indices.dtype)
-
- bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order)
- bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5
-
- vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order)
- nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1]
- nodes = np.empty(
- (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes),
- dtype=np.float64)
-
- # }}}
-
- grp_face_vertex_indices = mgrp.face_vertex_indices()
- grp_vertex_unit_coordinates = mgrp.vertex_unit_coordinates()
-
- # batch by face_id
-
- batch_base = 0
-
- for face_id in range(len(grp_face_vertex_indices)):
- batch_boundary_el_numbers_in_grp = np.array(
- [
- ibface_el
- for ibface_el, ibface_face in group_boundary_faces
- if ibface_face == face_id],
- dtype=np.intp)
-
- new_el_numbers = np.arange(
- batch_base,
- batch_base + len(batch_boundary_el_numbers_in_grp))
-
- # {{{ no per-element axes in these computations
-
- # Find boundary vertex indices
- loc_face_vertices = list(grp_face_vertex_indices[face_id])
-
- # Find unit nodes for boundary element
- face_vertex_unit_coordinates = \
- grp_vertex_unit_coordinates[loc_face_vertices]
-
- # Find A, b such that A [e_1 e_2] + b = [r_1 r_2]
- # (Notation assumes that the volume is 3D and the face is 2D.
- # Code does not.)
-
- b = face_vertex_unit_coordinates[0]
- A = (
- face_vertex_unit_coordinates[1:]
- - face_vertex_unit_coordinates[0]).T
-
- face_unit_nodes = (np.dot(A, bdry_unit_nodes_01).T + b).T
-
- resampling_mat = mp.resampling_matrix(
- vol_basis,
- face_unit_nodes, mgrp.unit_nodes)
-
- # }}}
-
- # {{{ build information for mesh element group
-
- # Find vertex_indices
- glob_face_vertices = mgrp.vertex_indices[
- batch_boundary_el_numbers_in_grp][:, loc_face_vertices]
- vertex_indices[new_el_numbers] = \
- vol_to_bdry_vertices[glob_face_vertices]
-
- # Find nodes
- nodes[:, new_el_numbers, :] = np.einsum(
- "ij,dej->dei",
- resampling_mat,
- mgrp.nodes[:, batch_boundary_el_numbers_in_grp, :])
-
- # }}}
-
- connection_data[igrp, face_id] = _ConnectionBatchData(
- group_source_element_indices=batch_boundary_el_numbers_in_grp,
- group_target_element_indices=new_el_numbers,
- A=A,
- b=b,
- )
-
- batch_base += len(batch_boundary_el_numbers_in_grp)
-
- bdry_mesh_group = SimplexElementGroup(
- mgrp.order, vertex_indices, nodes, unit_nodes=bdry_unit_nodes)
- bdry_mesh_groups.append(bdry_mesh_group)
-
- bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups)
-
- from meshmode.discretization import Discretization
- bdry_discr = Discretization(
- discr.cl_context, bdry_mesh, group_factory)
-
- connection = _build_boundary_connection(
- queue, discr, bdry_discr, connection_data)
-
- logger.info("building boundary connection: done")
-
- return bdry_mesh, bdry_discr, connection
-
-# }}}
-
-
-# {{{ refinement connection
-
-def make_refinement_connection(refiner, coarse_discr):
- pass
-
-# }}}
-
-# vim: foldmethod=marker
diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..269ed5d14c61e3b33ef3d7e6d4354b105e9b5f50
--- /dev/null
+++ b/meshmode/discretization/connection/__init__.py
@@ -0,0 +1,444 @@
+from __future__ import division, print_function, absolute_import
+
+__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+from six.moves import range, zip
+
+import numpy as np
+import pyopencl as cl
+import pyopencl.array # noqa
+from pytools import memoize_method, memoize_in
+
+from meshmode.discretization.connection.same_mesh import \
+ make_same_mesh_connection
+from meshmode.discretization.connection.face import (
+ FRESTR_INTERIOR_FACES, FRESTR_ALL_FACES,
+ make_face_restriction, make_face_to_all_faces_embedding)
+from meshmode.discretization.connection.opposite_face import \
+ make_opposite_face_connection
+
+import logging
+logger = logging.getLogger(__name__)
+
+
+__all__ = [
+ "DiscretizationConnection",
+ "make_same_mesh_connection",
+ "FRESTR_INTERIOR_FACES", "FRESTR_ALL_FACES",
+ "make_face_restriction",
+ "make_face_to_all_faces_embedding",
+ "make_opposite_face_connection"
+ ]
+
+__doc__ = """
+.. autoclass:: DiscretizationConnection
+
+.. autofunction:: make_same_mesh_connection
+
+.. autofunction:: FRESTR_INTERIOR_FACES
+.. autofunction:: FRESTR_ALL_FACES
+.. autofunction:: make_face_restriction
+.. autofunction:: make_face_to_all_faces_embedding
+
+.. autofunction:: make_opposite_face_connection
+
+Implementation details
+^^^^^^^^^^^^^^^^^^^^^^
+
+.. autoclass:: InterpolationBatch
+
+.. autoclass:: DiscretizationConnectionElementGroup
+"""
+
+
+# {{{ interpolation batch
+
+class InterpolationBatch(object):
+ """One interpolation batch captures how a batch of elements *within* an
+ element group should be an interpolated. Note that while it's possible that
+ an interpolation batch takes care of interpolating an entire element group
+ from source to target, that's not *necessarily* the case. Consider the case
+ of extracting boundary values of a discretization. For, say, a triangle, at
+ least three different interpolation batches are needed to cover boundary
+ edges that fall onto each of the three edges of the unit triangle.
+
+ .. attribute:: from_group_index
+
+ An integer indicating from which element group in the *from* discretization
+ the data should be interpolated.
+
+ .. attribute:: from_element_indices
+
+ ``element_id_t [nelements]``. (a :class:`pyopencl.array.Array`)
+ This contains the (group-local) element index (relative to
+ :attr:`from_group_index` from which this "*to*" element's data will be
+ interpolated.
+
+ .. attribute:: to_element_indices
+
+ ``element_id_t [nelements]``. (a :class:`pyopencl.array.Array`)
+ This contains the (group-local) element index to which this "*to*"
+ element's data will be interpolated.
+
+ .. attribute:: result_unit_nodes
+
+ A :class:`numpy.ndarray` of shape
+ ``(from_group.dim,to_group.nelements,to_group.nunit_nodes)``
+ storing the coordinates of the nodes (in unit coordinates
+ of the *from* reference element) from which the node
+ locations of this element should be interpolated.
+
+ .. autoattribute:: nelements
+
+ .. attribute:: to_element_face
+
+ *int* or *None*. (a :class:`pyopencl.array.Array` if existent) If this
+ interpolation batch targets interpolation *to* a face, then this number
+ captures the face number (on all elements referenced by
+ :attr:`from_element_indices` to which this batch interpolates. (Since
+ there is a fixed set of "from" unit nodes per batch, one batch will
+ always go to a single face index.)
+ """
+
+ def __init__(self, from_group_index, from_element_indices,
+ to_element_indices, result_unit_nodes, to_element_face):
+ self.from_group_index = from_group_index
+ self.from_element_indices = from_element_indices
+ self.to_element_indices = to_element_indices
+ self.result_unit_nodes = result_unit_nodes
+ self.to_element_face = to_element_face
+
+ @property
+ def nelements(self):
+ return len(self.from_element_indices)
+
+# }}}
+
+
+# {{{ connection element group
+
+class DiscretizationConnectionElementGroup(object):
+ """
+ .. attribute:: batches
+
+ A list of :class:`InterpolationBatch` instances.
+ """
+ def __init__(self, batches):
+ self.batches = batches
+
+# }}}
+
+
+# {{{ connection class
+
+class DiscretizationConnection(object):
+ """Data supporting an interpolation-like operation that takes in data on
+ one discretization and returns it on another. Implemented applications
+ include:
+
+ * upsampling/downsampling on the same mesh
+ * restricition to the boundary
+ * interpolation to a refined/coarsened mesh
+ * interpolation onto opposing faces
+
+ .. attribute:: from_discr
+
+ .. attribute:: to_discr
+
+ .. attribute:: groups
+
+ a list of :class:`DiscretizationConnectionElementGroup`
+ instances, with a one-to-one correspondence to the groups in
+ :attr:`to_discr`.
+
+ .. attribute:: is_surjective
+
+ A :class:`bool` indicating whether every output degree
+ of freedom is set by the connection.
+
+ .. automethod:: __call__
+
+ .. automethod:: full_resample_matrix
+
+ """
+
+ def __init__(self, from_discr, to_discr, groups, is_surjective):
+ if from_discr.cl_context != to_discr.cl_context:
+ raise ValueError("from_discr and to_discr must live in the "
+ "same OpenCL context")
+
+ self.cl_context = from_discr.cl_context
+
+ if from_discr.mesh.vertex_id_dtype != to_discr.mesh.vertex_id_dtype:
+ raise ValueError("from_discr and to_discr must agree on the "
+ "vertex_id_dtype")
+
+ if from_discr.mesh.element_id_dtype != to_discr.mesh.element_id_dtype:
+ raise ValueError("from_discr and to_discr must agree on the "
+ "element_id_dtype")
+
+ self.cl_context = from_discr.cl_context
+
+ self.from_discr = from_discr
+ self.to_discr = to_discr
+ self.groups = groups
+
+ self.is_surjective = is_surjective
+
+ @memoize_method
+ def _resample_matrix(self, to_group_index, ibatch_index):
+ import modepy as mp
+ ibatch = self.groups[to_group_index].batches[ibatch_index]
+ from_grp = self.from_discr.groups[ibatch.from_group_index]
+
+ result = mp.resampling_matrix(
+ from_grp.basis(),
+ ibatch.result_unit_nodes, from_grp.unit_nodes)
+
+ with cl.CommandQueue(self.cl_context) as queue:
+ return cl.array.to_device(queue, result).with_queue(None)
+
+ @memoize_method
+ def _resample_point_pick_indices(self, to_group_index, ibatch_index,
+ tol_multiplier=None):
+ """If :meth:`_resample_matrix` *R* is a row subset of a permutation matrix *P*,
+ return the index subset I so that, loosely, ``x[I] == R @ x``.
+
+ Will return *None* if no such index array exists, or a
+ :class:`pyopencl.array.Array` containing the index subset.
+ """
+
+ with cl.CommandQueue(self.cl_context) as queue:
+ mat = self._resample_matrix(to_group_index, ibatch_index).get(
+ queue=queue)
+
+ nrows, ncols = mat.shape
+ result = np.zeros(nrows, dtype=self.to_discr.mesh.element_id_dtype)
+
+ if tol_multiplier is None:
+ tol_multiplier = 50
+
+ tol = np.finfo(mat.dtype).eps * tol_multiplier
+
+ for irow in range(nrows):
+ one_indices, = np.where(np.abs(mat[irow] - 1) < tol)
+ zero_indices, = np.where(np.abs(mat[irow]) < tol)
+
+ if len(one_indices) != 1:
+ return None
+ if len(zero_indices) != ncols - 1:
+ return None
+
+ one_index, = one_indices
+ result[irow] = one_index
+
+ with cl.CommandQueue(self.cl_context) as queue:
+ return cl.array.to_device(queue, result).with_queue(None)
+
+ def full_resample_matrix(self, queue):
+ """Build a dense matrix representing this discretization connection.
+
+ .. warning::
+
+ On average, this will be exceedingly expensive (:math:`O(N^2)` in
+ the number *N* of discretization points) in terms of memory usage
+ and thus not what you'd typically want.
+ """
+
+ @memoize_in(self, "oversample_mat_knl")
+ def knl():
+ import loopy as lp
+ knl = lp.make_kernel(
+ """{[k,i,j]:
+ 0<=k<nelements and
+ 0<=i<n_to_nodes and
+ 0<=j<n_from_nodes}""",
+ "result[itgt_base + to_element_indices[k]*n_to_nodes + i, \
+ isrc_base + from_element_indices[k]*n_from_nodes + j] \
+ = resample_mat[i, j]",
+ [
+ lp.GlobalArg("result", None,
+ shape="nnodes_tgt, nnodes_src",
+ offset=lp.auto),
+ lp.ValueArg("itgt_base,isrc_base", np.int32),
+ lp.ValueArg("nnodes_tgt,nnodes_src", np.int32),
+ "...",
+ ],
+ name="oversample_mat")
+
+ knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+ return lp.tag_inames(knl, dict(k="g.0"))
+
+ result = cl.array.zeros(
+ queue,
+ (self.to_discr.nnodes, self.from_discr.nnodes),
+ dtype=self.to_discr.real_dtype)
+
+ for i_tgrp, (tgrp, cgrp) in enumerate(
+ zip(self.to_discr.groups, self.groups)):
+ for i_batch, batch in enumerate(cgrp.batches):
+ if len(batch.from_element_indices):
+ if not len(batch.from_element_indices):
+ continue
+
+ sgrp = self.from_discr.groups[batch.from_group_index]
+
+ knl()(queue,
+ resample_mat=self._resample_matrix(i_tgrp, i_batch),
+ result=result,
+ itgt_base=tgrp.node_nr_base,
+ isrc_base=sgrp.node_nr_base,
+ from_element_indices=batch.from_element_indices,
+ to_element_indices=batch.to_element_indices)
+
+ return result
+
+ def __call__(self, queue, vec):
+ @memoize_in(self, "resample_by_mat_knl")
+ def mat_knl():
+ import loopy as lp
+ knl = lp.make_kernel(
+ """{[k,i,j]:
+ 0<=k<nelements and
+ 0<=i<n_to_nodes and
+ 0<=j<n_from_nodes}""",
+ "result[to_element_indices[k], i] \
+ = sum(j, resample_mat[i, j] \
+ * vec[from_element_indices[k], j])",
+ [
+ lp.GlobalArg("result", None,
+ shape="nelements_result, n_to_nodes",
+ offset=lp.auto),
+ lp.GlobalArg("vec", None,
+ shape="nelements_vec, n_from_nodes",
+ offset=lp.auto),
+ lp.ValueArg("nelements_result", np.int32),
+ lp.ValueArg("nelements_vec", np.int32),
+ "...",
+ ],
+ name="resample_by_mat")
+
+ knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+ return lp.tag_inames(knl, dict(k="g.0"))
+
+ @memoize_in(self, "resample_by_picking_knl")
+ def pick_knl():
+ import loopy as lp
+ knl = lp.make_kernel(
+ """{[k,i,j]:
+ 0<=k<nelements and
+ 0<=i<n_to_nodes}""",
+ "result[to_element_indices[k], i] \
+ = vec[from_element_indices[k], pick_list[i]]",
+ [
+ lp.GlobalArg("result", None,
+ shape="nelements_result, n_to_nodes",
+ offset=lp.auto),
+ lp.GlobalArg("vec", None,
+ shape="nelements_vec, n_from_nodes",
+ offset=lp.auto),
+ lp.ValueArg("nelements_result", np.int32),
+ lp.ValueArg("nelements_vec", np.int32),
+ lp.ValueArg("n_from_nodes", np.int32),
+ "...",
+ ],
+ name="resample_by_picking")
+
+ knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+ return lp.tag_inames(knl, dict(k="g.0"))
+
+ if not isinstance(vec, cl.array.Array):
+ return vec
+
+ if self.is_surjective:
+ result = self.to_discr.empty(dtype=vec.dtype)
+ else:
+ result = self.to_discr.zeros(queue, dtype=vec.dtype)
+
+ if vec.shape != (self.from_discr.nnodes,):
+ raise ValueError("invalid shape of incoming resampling data")
+
+ for i_tgrp, (tgrp, cgrp) in enumerate(
+ zip(self.to_discr.groups, self.groups)):
+ for i_batch, batch in enumerate(cgrp.batches):
+ sgrp = self.from_discr.groups[batch.from_group_index]
+
+ if not len(batch.from_element_indices):
+ continue
+
+ point_pick_indices = self._resample_point_pick_indices(
+ i_tgrp, i_batch)
+
+ if point_pick_indices is None:
+ mat_knl()(queue,
+ resample_mat=self._resample_matrix(i_tgrp, i_batch),
+ result=tgrp.view(result), vec=sgrp.view(vec),
+ from_element_indices=batch.from_element_indices,
+ to_element_indices=batch.to_element_indices)
+
+ else:
+ pick_knl()(queue,
+ pick_list=point_pick_indices,
+ result=tgrp.view(result), vec=sgrp.view(vec),
+ from_element_indices=batch.from_element_indices,
+ to_element_indices=batch.to_element_indices)
+
+ return result
+
+# }}}
+
+
+# {{{ check connection
+
+def check_connection(connection):
+ from_discr = connection.from_discr
+ to_discr = connection.to_discr
+
+ assert len(connection.groups) == len(to_discr.groups)
+
+ with cl.CommandQueue(to_discr.cl_context) as queue:
+ for cgrp, tgrp in zip(connection.groups, to_discr.groups):
+ for batch in cgrp.batches:
+ fgrp = from_discr.groups[batch.from_group_index]
+
+ from_element_indices = batch.from_element_indices.get(queue)
+ to_element_indices = batch.to_element_indices.get(queue)
+
+ assert (0 <= from_element_indices).all()
+ assert (0 <= to_element_indices).all()
+ assert (from_element_indices < fgrp.nelements).all()
+ assert (to_element_indices < tgrp.nelements).all()
+ if batch.to_element_face is not None:
+ assert 0 <= batch.to_element_face < fgrp.mesh_el_group.nfaces
+
+# }}}
+
+
+# {{{ refinement connection
+
+def make_refinement_connection(refiner, coarse_discr):
+ pass
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py
new file mode 100644
index 0000000000000000000000000000000000000000..52d7575d6dd83df9ed9868d063c6fa3adbe985cf
--- /dev/null
+++ b/meshmode/discretization/connection/face.py
@@ -0,0 +1,488 @@
+from __future__ import division, print_function, absolute_import
+
+__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import six
+from six.moves import range, zip
+
+from pytools import Record
+
+import numpy as np
+import pyopencl as cl
+import pyopencl.array # noqa
+import modepy as mp
+
+import logging
+logger = logging.getLogger(__name__)
+
+
+class FRESTR_INTERIOR_FACES:
+ """A special value to pass to
+ :func:`meshmode.discretization.connection.make_face_restriction`
+ to produce a discretization consisting of all interior faces
+ in a discretization.
+ """
+
+
+class FRESTR_ALL_FACES:
+ """A special value to pass to
+ :func:`meshmode.discretization.connection.make_face_restriction`
+ to produce a discretization consisting of all faces (interior and boundary)
+ faces in a discretization.
+ """
+
+
+# {{{ boundary connection
+
+class _ConnectionBatchData(Record):
+ pass
+
+
+def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data,
+ per_face_groups):
+ from meshmode.discretization.connection import (
+ InterpolationBatch, DiscretizationConnectionElementGroup,
+ DiscretizationConnection)
+
+ ibdry_grp = 0
+ batches = []
+
+ connection_groups = []
+ for igrp, vol_grp in enumerate(vol_discr.groups):
+ mgrp = vol_grp.mesh_el_group
+
+ for face_id in range(mgrp.nfaces):
+ bdry_grp = bdry_discr.groups[ibdry_grp]
+ data = connection_data[igrp, face_id]
+
+ bdry_unit_nodes_01 = (bdry_grp.unit_nodes + 1)*0.5
+ result_unit_nodes = (np.dot(data.A, bdry_unit_nodes_01).T + data.b).T
+
+ batches.append(
+ InterpolationBatch(
+ from_group_index=igrp,
+ from_element_indices=cl.array.to_device(
+ queue,
+ vol_grp.mesh_el_group.element_nr_base
+ + data.group_source_element_indices)
+ .with_queue(None),
+ to_element_indices=cl.array.to_device(
+ queue,
+ data.group_target_element_indices)
+ .with_queue(None),
+ result_unit_nodes=result_unit_nodes,
+ to_element_face=face_id
+ ))
+
+ is_last_face = face_id + 1 == mgrp.nfaces
+
+ if per_face_groups or is_last_face:
+ connection_groups.append(
+ DiscretizationConnectionElementGroup(batches))
+ batches = []
+
+ ibdry_grp += 1
+
+ assert ibdry_grp == len(bdry_discr.groups)
+
+ return DiscretizationConnection(
+ vol_discr, bdry_discr, connection_groups,
+ is_surjective=True)
+
+# }}}
+
+
+# {{{ pull together boundary vertices
+
+def _get_face_vertices(mesh, boundary_tag):
+ # a set of volume vertex numbers
+ bdry_vertex_vol_nrs = set()
+
+ if boundary_tag not in [FRESTR_INTERIOR_FACES, FRESTR_ALL_FACES]:
+ # {{{ boundary faces
+
+ btag_bit = mesh.boundary_tag_bit(boundary_tag)
+
+ for fagrp_map in mesh.facial_adjacency_groups:
+ bdry_grp = fagrp_map.get(None)
+ if bdry_grp is None:
+ continue
+
+ assert (bdry_grp.neighbors < 0).all()
+
+ grp = mesh.groups[bdry_grp.igroup]
+
+ nb_el_bits = -bdry_grp.neighbors
+ face_relevant_flags = (nb_el_bits & btag_bit) != 0
+
+ for iface, fvi in enumerate(grp.face_vertex_indices()):
+ bdry_vertex_vol_nrs.update(
+ grp.vertex_indices
+ [bdry_grp.elements[face_relevant_flags]]
+ [:, np.array(fvi, dtype=np.intp)]
+ .flat)
+
+ return np.array(sorted(bdry_vertex_vol_nrs), dtype=np.intp)
+
+ # }}}
+ else:
+ # For INTERIOR_FACES, this is likely every vertex in the book.
+ # Don't ever bother trying to cut the list down.
+ # For ALL_FACES, it literally is every single vertex.
+
+ return np.arange(mesh.nvertices, dtype=np.intp)
+
+# }}}
+
+
+def make_face_restriction(discr, group_factory, boundary_tag,
+ per_face_groups=False):
+ """Create a mesh, a discretization and a connection to restrict
+ a function on *discr* to its values on the edges of element faces
+ denoted by *boundary_tag*.
+
+ :arg boundary_tag: The boundary tag for which to create a face
+ restriction. May be
+ :class:`meshmode.discretization.connection.INTERIOR_FACES`
+ to indicate interior faces, or
+ :class:`meshmode.discretization.connection.ALL_FACES`
+ to make a discretization consisting of all (interior and
+ boundary) faces.
+
+ :arg per_face_groups: If *True*, the resulting discretization is
+ guaranteed to have groups organized as::
+
+ (grp0, face0), (grp0, face1), ... (grp0, faceN),
+ (grp1, face0), (grp1, face1), ... (grp1, faceN), ...
+
+ each with the elements in the same order as the originating
+ group. If *False*, volume and boundary groups correspond with
+ each other one-to-one, and an interpolation batch is created
+ per face.
+
+ :return: a :class:`meshmode.discretization.connection.DiscretizationConnection`
+ representing the new connection. The new boundary discretization can be
+ obtained from the
+ :attr:`meshmode.discretization.connection.DiscretizationConnection.to_discr`
+ attribute of the return value, and the corresponding new boundary mesh
+ from that.
+
+ """
+
+ if boundary_tag is None:
+ boundary_tag = FRESTR_INTERIOR_FACES
+ from warnings import warn
+ warn("passing *None* for boundary_tag is deprecated--pass "
+ "INTERIOR_FACES instead",
+ DeprecationWarning, stacklevel=2)
+
+ logger.info("building face restriction: start")
+
+ # {{{ gather boundary vertices
+
+ bdry_vertex_vol_nrs = _get_face_vertices(discr.mesh, boundary_tag)
+
+ vol_to_bdry_vertices = np.empty(
+ discr.mesh.vertices.shape[-1],
+ discr.mesh.vertices.dtype)
+ vol_to_bdry_vertices.fill(-1)
+ vol_to_bdry_vertices[bdry_vertex_vol_nrs] = np.arange(
+ len(bdry_vertex_vol_nrs), dtype=np.intp)
+
+ bdry_vertices = discr.mesh.vertices[:, bdry_vertex_vol_nrs]
+
+ # }}}
+
+ from meshmode.mesh import Mesh, SimplexElementGroup
+ bdry_mesh_groups = []
+ connection_data = {}
+
+ btag_bit = discr.mesh.boundary_tag_bit(boundary_tag)
+
+ for igrp, (grp, fagrp_map) in enumerate(
+ zip(discr.groups, discr.mesh.facial_adjacency_groups)):
+
+ mgrp = grp.mesh_el_group
+
+ if not isinstance(mgrp, SimplexElementGroup):
+ raise NotImplementedError("can only take boundary of "
+ "SimplexElementGroup-based meshes")
+
+ # {{{ pull together per-group face lists
+
+ group_boundary_faces = []
+
+ if boundary_tag is FRESTR_INTERIOR_FACES:
+ for fagrp in six.itervalues(fagrp_map):
+ if fagrp.ineighbor_group is None:
+ # boundary faces -> not looking for those
+ continue
+
+ group_boundary_faces.extend(
+ zip(fagrp.elements, fagrp.element_faces))
+
+ elif boundary_tag is FRESTR_ALL_FACES:
+ group_boundary_faces.extend(
+ (iel, iface)
+ for iface in range(grp.mesh_el_group.nfaces)
+ for iel in range(grp.nelements)
+ )
+
+ else:
+ bdry_grp = fagrp_map.get(None)
+ if bdry_grp is not None:
+ nb_el_bits = -bdry_grp.neighbors
+ face_relevant_flags = (nb_el_bits & btag_bit) != 0
+
+ group_boundary_faces.extend(
+ zip(
+ bdry_grp.elements[face_relevant_flags],
+ bdry_grp.element_faces[face_relevant_flags]))
+
+ # }}}
+
+ grp_face_vertex_indices = mgrp.face_vertex_indices()
+ grp_vertex_unit_coordinates = mgrp.vertex_unit_coordinates()
+
+ batch_base = 0
+
+ # group by face_id
+
+ for face_id in range(mgrp.nfaces):
+ batch_boundary_el_numbers_in_grp = np.array(
+ [
+ ibface_el
+ for ibface_el, ibface_face in group_boundary_faces
+ if ibface_face == face_id],
+ dtype=np.intp)
+
+ # {{{ Preallocate arrays for mesh group
+
+ nbatch_elements = len(batch_boundary_el_numbers_in_grp)
+
+ if per_face_groups or face_id == 0:
+ if per_face_groups:
+ ngroup_bdry_elements = nbatch_elements
+ else:
+ ngroup_bdry_elements = len(group_boundary_faces)
+
+ vertex_indices = np.empty(
+ (ngroup_bdry_elements, mgrp.dim+1-1),
+ mgrp.vertex_indices.dtype)
+
+ bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order)
+ bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5
+
+ vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order)
+ nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1]
+ nodes = np.empty(
+ (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes),
+ dtype=np.float64)
+
+ # }}}
+
+ new_el_numbers = batch_base + np.arange(nbatch_elements)
+ if not per_face_groups:
+ batch_base += nbatch_elements
+
+ # {{{ no per-element axes in these computations
+
+ # Find boundary vertex indices
+ loc_face_vertices = list(grp_face_vertex_indices[face_id])
+
+ # Find unit nodes for boundary element
+ face_vertex_unit_coordinates = \
+ grp_vertex_unit_coordinates[loc_face_vertices]
+
+ # Find A, b such that A [e_1 e_2] + b = [r_1 r_2]
+ # (Notation assumes that the volume is 3D and the face is 2D.
+ # Code does not.)
+
+ b = face_vertex_unit_coordinates[0]
+ A = ( # noqa
+ face_vertex_unit_coordinates[1:]
+ - face_vertex_unit_coordinates[0]).T
+
+ face_unit_nodes = (np.dot(A, bdry_unit_nodes_01).T + b).T
+
+ resampling_mat = mp.resampling_matrix(
+ vol_basis,
+ face_unit_nodes, mgrp.unit_nodes)
+
+ # }}}
+
+ # {{{ build information for mesh element group
+
+ # Find vertex_indices
+ glob_face_vertices = mgrp.vertex_indices[
+ batch_boundary_el_numbers_in_grp][:, loc_face_vertices]
+ vertex_indices[new_el_numbers] = \
+ vol_to_bdry_vertices[glob_face_vertices]
+
+ # Find nodes
+ nodes[:, new_el_numbers, :] = np.einsum(
+ "ij,dej->dei",
+ resampling_mat,
+ mgrp.nodes[:, batch_boundary_el_numbers_in_grp, :])
+
+ # }}}
+
+ connection_data[igrp, face_id] = _ConnectionBatchData(
+ group_source_element_indices=batch_boundary_el_numbers_in_grp,
+ group_target_element_indices=new_el_numbers,
+ A=A,
+ b=b,
+ )
+
+ is_last_face = face_id + 1 == mgrp.nfaces
+
+ if per_face_groups or is_last_face:
+ bdry_mesh_group = SimplexElementGroup(
+ mgrp.order, vertex_indices, nodes,
+ unit_nodes=bdry_unit_nodes)
+ bdry_mesh_groups.append(bdry_mesh_group)
+
+ bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups)
+
+ from meshmode.discretization import Discretization
+ bdry_discr = Discretization(
+ discr.cl_context, bdry_mesh, group_factory)
+
+ with cl.CommandQueue(discr.cl_context) as queue:
+ connection = _build_boundary_connection(
+ queue, discr, bdry_discr, connection_data,
+ per_face_groups)
+
+ logger.info("building face restriction: done")
+
+ return connection
+
+# }}}
+
+
+# {{{ face -> all_faces connection
+
+def make_face_to_all_faces_embedding(faces_connection, all_faces_discr):
+ """Return a
+ :class:`meshmode.discretization.connection.DiscretizationConnection`
+ connecting a discretization containing some faces of a discretization
+ to one containing all faces.
+
+ :arg faces_connection: must be the (connection) result of calling
+ :func:`meshmode.discretization.connection.make_face_restriction`
+ with
+ :class:`meshmode.discretization.connection.FRESTR_INTERIOR_FACES`
+ or a boundary tag.
+ :arg all_faces_discr: must be the (discretization) result of calling
+ :func:`meshmode.discretization.connection.make_face_restriction`
+ with
+ :class:`meshmode.discretization.connection.FRESTR_ALL_FACES`
+ for the same volume discretization as the one from which
+ *faces_discr* was obtained.
+ """
+
+ vol_discr = faces_connection.from_discr
+ faces_discr = faces_connection.to_discr
+
+ per_face_groups = (
+ len(vol_discr.groups) != len(faces_discr.groups))
+
+ if len(faces_discr.groups) != len(all_faces_discr.groups):
+ raise ValueError("faces_discr and all_faces_discr must have the "
+ "same number of groups")
+ if len(faces_connection.groups) != len(all_faces_discr.groups):
+ raise ValueError("faces_connection and all_faces_discr must have the "
+ "same number of groups")
+
+ from meshmode.discretization.connection import (
+ DiscretizationConnection,
+ DiscretizationConnectionElementGroup,
+ InterpolationBatch)
+
+ i_faces_grp = 0
+
+ with cl.CommandQueue(vol_discr.cl_context) as queue:
+ groups = []
+ for ivol_grp, vol_grp in enumerate(vol_discr.groups):
+ batches = []
+
+ nfaces = vol_grp.mesh_el_group.nfaces
+ for iface in range(nfaces):
+ faces_grp = faces_discr.groups[i_faces_grp]
+ all_faces_grp = all_faces_discr.groups[i_faces_grp]
+
+ if per_face_groups:
+ assert len(faces_connection.groups[i_faces_grp].batches) == 1
+ else:
+ assert (len(faces_connection.groups[i_faces_grp].batches)
+ == nfaces)
+
+ assert np.array_equal(
+ faces_grp.unit_nodes, all_faces_grp.unit_nodes)
+
+ # {{{ find src_batch
+
+ src_batches = faces_connection.groups[i_faces_grp].batches
+ if per_face_groups:
+ src_batch, = src_batches
+ else:
+ src_batch = src_batches[iface]
+ del src_batches
+
+ # }}}
+
+ if per_face_groups:
+ to_element_indices = src_batch.from_element_indices
+ else:
+ assert all_faces_grp.nelements == nfaces * vol_grp.nelements
+
+ to_element_indices = (
+ vol_grp.nelements*iface
+ + src_batch.from_element_indices.with_queue(queue)
+ ).with_queue(None)
+
+ batches.append(
+ InterpolationBatch(
+ from_group_index=i_faces_grp,
+ from_element_indices=src_batch.to_element_indices,
+ to_element_indices=to_element_indices,
+ result_unit_nodes=all_faces_grp.unit_nodes,
+ to_element_face=None))
+
+ is_last_face = iface + 1 == nfaces
+ if per_face_groups or is_last_face:
+ groups.append(
+ DiscretizationConnectionElementGroup(batches=batches))
+ batches = []
+
+ i_faces_grp += 1
+
+ return DiscretizationConnection(
+ faces_connection.to_discr,
+ all_faces_discr,
+ groups,
+ is_surjective=False)
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py
new file mode 100644
index 0000000000000000000000000000000000000000..70728522b2c06f119a5b20598616e2c464eea249
--- /dev/null
+++ b/meshmode/discretization/connection/opposite_face.py
@@ -0,0 +1,395 @@
+from __future__ import division, print_function, absolute_import
+
+__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+from six.moves import range
+
+import numpy as np
+import numpy.linalg as la
+import pyopencl as cl
+import pyopencl.array # noqa
+
+import logging
+logger = logging.getLogger(__name__)
+
+
+# {{{ opposite-face connection
+
+def _make_cross_face_batches(
+ queue, vol_discr, bdry_discr,
+ i_tgt_grp, i_src_grp,
+ i_face_tgt,
+ adj_grp,
+ vbc_tgt_grp_face_batch, src_grp_el_lookup):
+
+ # {{{ index wrangling
+
+ # Assert that the adjacency group and the restriction
+ # interpolation batch and the adjacency group have the same
+ # element ordering.
+
+ adj_grp_tgt_flags = adj_grp.element_faces == i_face_tgt
+
+ assert (
+ np.array_equal(
+ adj_grp.elements[adj_grp_tgt_flags],
+ vbc_tgt_grp_face_batch.from_element_indices
+ .get(queue=queue)))
+
+ # find to_element_indices
+
+ to_bdry_element_indices = (
+ vbc_tgt_grp_face_batch.to_element_indices
+ .get(queue=queue))
+
+ # find from_element_indices
+
+ from_vol_element_indices = adj_grp.neighbors[adj_grp_tgt_flags]
+ from_element_faces = adj_grp.neighbor_faces[adj_grp_tgt_flags]
+
+ from_bdry_element_indices = src_grp_el_lookup[
+ from_vol_element_indices, from_element_faces]
+
+ # }}}
+
+ # {{{ visualization (for debugging)
+
+ if 0:
+ print("TVE", adj_grp.elements[adj_grp_tgt_flags])
+ print("TBE", to_bdry_element_indices)
+ print("FVE", from_vol_element_indices)
+ from meshmode.mesh.visualization import draw_2d_mesh
+ import matplotlib.pyplot as pt
+ draw_2d_mesh(vol_discr.mesh, draw_element_numbers=True,
+ set_bounding_box=True,
+ draw_vertex_numbers=False,
+ draw_face_numbers=True,
+ fill=None)
+ pt.figure()
+
+ draw_2d_mesh(bdry_discr.mesh, draw_element_numbers=True,
+ set_bounding_box=True,
+ draw_vertex_numbers=False,
+ draw_face_numbers=True,
+ fill=None)
+
+ pt.show()
+ # }}}
+
+ # {{{ invert face map (using Gauss-Newton)
+
+ to_bdry_nodes = (
+ # FIXME: This should view-then-transfer (but PyOpenCL doesn't do
+ # non-contiguous transfers for now).
+ bdry_discr.groups[i_tgt_grp].view(
+ bdry_discr.nodes().get(queue=queue))
+ [:, to_bdry_element_indices])
+
+ tol = 1e4 * np.finfo(to_bdry_nodes.dtype).eps
+
+ from_mesh_grp = bdry_discr.mesh.groups[i_src_grp]
+ from_grp = bdry_discr.groups[i_src_grp]
+
+ dim = from_grp.dim
+ ambient_dim, nelements, nto_unit_nodes = to_bdry_nodes.shape
+
+ initial_guess = np.mean(from_mesh_grp.vertex_unit_coordinates(), axis=0)
+ from_unit_nodes = np.empty((dim, nelements, nto_unit_nodes))
+ from_unit_nodes[:] = initial_guess.reshape(-1, 1, 1)
+
+ import modepy as mp
+ from_vdm = mp.vandermonde(from_grp.basis(), from_grp.unit_nodes)
+ from_inv_t_vdm = la.inv(from_vdm.T)
+ from_nfuncs = len(from_grp.basis())
+
+ # (ambient_dim, nelements, nfrom_unit_nodes)
+ from_bdry_nodes = (
+ # FIXME: This should view-then-transfer (but PyOpenCL doesn't do
+ # non-contiguous transfers for now).
+ bdry_discr.groups[i_src_grp].view(
+ bdry_discr.nodes().get(queue=queue))
+ [:, from_bdry_element_indices])
+
+ def apply_map(unit_nodes):
+ # unit_nodes: (dim, nelements, nto_unit_nodes)
+
+ # basis_at_unit_nodes
+ basis_at_unit_nodes = np.empty((from_nfuncs, nelements, nto_unit_nodes))
+
+ for i, f in enumerate(from_grp.basis()):
+ basis_at_unit_nodes[i] = (
+ f(unit_nodes.reshape(dim, -1))
+ .reshape(nelements, nto_unit_nodes))
+
+ intp_coeffs = np.einsum("fj,jet->fet", from_inv_t_vdm, basis_at_unit_nodes)
+
+ # If we're interpolating 1, we had better get 1 back.
+ one_deviation = np.abs(np.sum(intp_coeffs, axis=0) - 1)
+ assert (one_deviation < tol).all(), np.max(one_deviation)
+
+ return np.einsum("fet,aef->aet", intp_coeffs, from_bdry_nodes)
+
+ def get_map_jacobian(unit_nodes):
+ # unit_nodes: (dim, nelements, nto_unit_nodes)
+
+ # basis_at_unit_nodes
+ dbasis_at_unit_nodes = np.empty(
+ (dim, from_nfuncs, nelements, nto_unit_nodes))
+
+ for i, df in enumerate(from_grp.grad_basis()):
+ df_result = df(unit_nodes.reshape(dim, -1))
+
+ for rst_axis, df_r in enumerate(df_result):
+ dbasis_at_unit_nodes[rst_axis, i] = (
+ df_r.reshape(nelements, nto_unit_nodes))
+
+ dintp_coeffs = np.einsum(
+ "fj,rjet->rfet", from_inv_t_vdm, dbasis_at_unit_nodes)
+
+ return np.einsum("rfet,aef->raet", dintp_coeffs, from_bdry_nodes)
+
+ # {{{ test map applier and jacobian
+
+ if 0:
+ u = from_unit_nodes
+ f = apply_map(u)
+ for h in [1e-1, 1e-2]:
+ du = h*np.random.randn(*u.shape)
+
+ f_2 = apply_map(u+du)
+
+ jf = get_map_jacobian(u)
+
+ f2_2 = f + np.einsum("raet,ret->aet", jf, du)
+
+ print(h, la.norm((f_2-f2_2).ravel()))
+
+ # }}}
+
+ # {{{ visualize initial guess
+
+ if 0:
+ import matplotlib.pyplot as pt
+ guess = apply_map(from_unit_nodes)
+ goals = to_bdry_nodes
+
+ from meshmode.discretization.visualization import draw_curve
+ draw_curve(bdry_discr)
+
+ pt.plot(guess[0].reshape(-1), guess[1].reshape(-1), "or")
+ pt.plot(goals[0].reshape(-1), goals[1].reshape(-1), "og")
+ pt.plot(from_bdry_nodes[0].reshape(-1), from_bdry_nodes[1].reshape(-1), "o",
+ color="purple")
+ pt.show()
+
+ # }}}
+
+ logger.info("make_opposite_face_connection: begin gauss-newton")
+
+ niter = 0
+ while True:
+ resid = apply_map(from_unit_nodes) - to_bdry_nodes
+
+ df = get_map_jacobian(from_unit_nodes)
+ df_inv_resid = np.empty_like(from_unit_nodes)
+
+ # For the 1D/2D accelerated versions, we'll use the normal
+ # equations and Cramer's rule. If you're looking for high-end
+ # numerics, look no further than meshmode.
+
+ if dim == 1:
+ # A is df.T
+ ata = np.einsum("iket,jket->ijet", df, df)
+ atb = np.einsum("iket,ket->iet", df, resid)
+
+ df_inv_resid = atb / ata[0, 0]
+
+ elif dim == 2:
+ # A is df.T
+ ata = np.einsum("iket,jket->ijet", df, df)
+ atb = np.einsum("iket,ket->iet", df, resid)
+
+ det = ata[0, 0]*ata[1, 1] - ata[0, 1]*ata[1, 0]
+
+ df_inv_resid = np.empty_like(from_unit_nodes)
+ df_inv_resid[0] = 1/det * (ata[1, 1] * atb[0] - ata[1, 0]*atb[1])
+ df_inv_resid[1] = 1/det * (-ata[0, 1] * atb[0] + ata[0, 0]*atb[1])
+
+ else:
+ # The boundary of a 3D mesh is 2D, so that's the
+ # highest-dimensional case we genuinely care about.
+ #
+ # This stinks, performance-wise, because it's not vectorized.
+ # But we'll only hit it for boundaries of 4+D meshes, in which
+ # case... good luck. :)
+ for e in range(nelements):
+ for t in range(nto_unit_nodes):
+ df_inv_resid[:, e, t], _, _, _ = \
+ la.lstsq(df[:, :, e, t].T, resid[:, e, t])
+
+ from_unit_nodes = from_unit_nodes - df_inv_resid
+
+ max_resid = np.max(np.abs(resid))
+ logger.debug("gauss-newton residual: %g" % max_resid)
+
+ if max_resid < tol:
+ logger.info("make_opposite_face_connection: gauss-newton: done, "
+ "final residual: %g" % max_resid)
+ break
+
+ niter += 1
+ if niter > 10:
+ raise RuntimeError("Gauss-Newton (for finding opposite-face reference "
+ "coordinates) did not converge")
+
+ # }}}
+
+ # {{{ find groups of from_unit_nodes
+
+ def to_dev(ary):
+ return cl.array.to_device(queue, ary, array_queue=None)
+
+ done_elements = np.zeros(nelements, dtype=np.bool)
+ while True:
+ todo_elements, = np.where(~done_elements)
+ if not len(todo_elements):
+ return
+
+ template_unit_nodes = from_unit_nodes[:, todo_elements[0], :]
+
+ unit_node_dist = np.max(np.max(np.abs(
+ from_unit_nodes[:, todo_elements, :]
+ -
+ template_unit_nodes.reshape(dim, 1, -1)),
+ axis=2), axis=0)
+
+ close_els = todo_elements[unit_node_dist < tol]
+ done_elements[close_els] = True
+
+ unit_node_dist = np.max(np.max(np.abs(
+ from_unit_nodes[:, todo_elements, :]
+ -
+ template_unit_nodes.reshape(dim, 1, -1)),
+ axis=2), axis=0)
+
+ from meshmode.discretization.connection import InterpolationBatch
+ yield InterpolationBatch(
+ from_group_index=i_src_grp,
+ from_element_indices=to_dev(from_bdry_element_indices[close_els]),
+ to_element_indices=to_dev(to_bdry_element_indices[close_els]),
+ result_unit_nodes=template_unit_nodes,
+ to_element_face=None)
+
+ # }}}
+
+
+def _find_ibatch_for_face(vbc_tgt_grp_batches, iface):
+ vbc_tgt_grp_face_batches = [
+ batch
+ for batch in vbc_tgt_grp_batches
+ if batch.to_element_face == iface]
+
+ assert len(vbc_tgt_grp_face_batches) == 1
+
+ vbc_tgt_grp_face_batch, = vbc_tgt_grp_face_batches
+
+ return vbc_tgt_grp_face_batch
+
+
+def _make_el_lookup_table(queue, connection, igrp):
+ from_nelements = connection.from_discr.groups[igrp].nelements
+ from_nfaces = connection.from_discr.mesh.groups[igrp].nfaces
+
+ iel_lookup = np.empty((from_nelements, from_nfaces),
+ dtype=connection.from_discr.mesh.element_id_dtype)
+ iel_lookup.fill(-1)
+
+ for ibatch, batch in enumerate(connection.groups[igrp].batches):
+ from_element_indices = batch.from_element_indices.get(queue=queue)
+ iel_lookup[from_element_indices, batch.to_element_face] = \
+ batch.to_element_indices.get(queue=queue)
+
+ return iel_lookup
+
+
+def make_opposite_face_connection(volume_to_bdry_conn):
+ """Given a boundary restriction connection *volume_to_bdry_conn*,
+ return a :class:`DiscretizationConnection` that performs data
+ exchange across opposite faces.
+ """
+
+ vol_discr = volume_to_bdry_conn.from_discr
+ vol_mesh = vol_discr.mesh
+ bdry_discr = volume_to_bdry_conn.to_discr
+
+ # make sure we were handed a volume-to-boundary connection
+ for i_tgrp, conn_grp in enumerate(volume_to_bdry_conn.groups):
+ for batch in conn_grp.batches:
+ assert batch.from_group_index == i_tgrp
+ assert batch.to_element_face is not None
+
+ ngrps = len(volume_to_bdry_conn.groups)
+ assert ngrps == len(vol_discr.groups)
+ assert ngrps == len(bdry_discr.groups)
+
+ # One interpolation batch in this connection corresponds
+ # to a key (i_tgt_grp,) (i_src_grp, i_face_tgt,)
+
+ with cl.CommandQueue(vol_discr.cl_context) as queue:
+ # a list of batches for each group
+ groups = [[] for i_tgt_grp in range(ngrps)]
+
+ for i_src_grp in range(ngrps):
+ src_grp_el_lookup = _make_el_lookup_table(
+ queue, volume_to_bdry_conn, i_src_grp)
+
+ for i_tgt_grp in range(ngrps):
+ vbc_tgt_grp_batches = volume_to_bdry_conn.groups[i_tgt_grp].batches
+
+ adj_grp = vol_mesh.facial_adjacency_groups[i_tgt_grp][i_src_grp]
+
+ for i_face_tgt in range(vol_mesh.groups[i_tgt_grp].nfaces):
+ vbc_tgt_grp_face_batch = _find_ibatch_for_face(
+ vbc_tgt_grp_batches, i_face_tgt)
+
+ groups[i_tgt_grp].extend(
+ _make_cross_face_batches(
+ queue, vol_discr, bdry_discr,
+ i_tgt_grp, i_src_grp,
+ i_face_tgt,
+ adj_grp,
+ vbc_tgt_grp_face_batch, src_grp_el_lookup))
+
+ from meshmode.discretization.connection import (
+ DiscretizationConnection, DiscretizationConnectionElementGroup)
+ return DiscretizationConnection(
+ from_discr=bdry_discr,
+ to_discr=bdry_discr,
+ groups=[
+ DiscretizationConnectionElementGroup(batches=batches)
+ for batches in groups],
+ is_surjective=True)
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/meshmode/discretization/connection/same_mesh.py b/meshmode/discretization/connection/same_mesh.py
new file mode 100644
index 0000000000000000000000000000000000000000..3b73e34fefd5115e01b97728ab414cc8c77fd178
--- /dev/null
+++ b/meshmode/discretization/connection/same_mesh.py
@@ -0,0 +1,65 @@
+from __future__ import division, print_function, absolute_import
+
+__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import numpy as np
+import pyopencl as cl
+import pyopencl.array # noqa
+
+
+# {{{ same-mesh constructor
+
+def make_same_mesh_connection(to_discr, from_discr):
+ from meshmode.discretization.connection import (
+ InterpolationBatch, DiscretizationConnectionElementGroup,
+ DiscretizationConnection)
+
+ if from_discr.mesh is not to_discr.mesh:
+ raise ValueError("from_discr and to_discr must be based on "
+ "the same mesh")
+
+ assert to_discr.cl_context == from_discr.cl_context
+
+ with cl.CommandQueue(to_discr.cl_context) as queue:
+ groups = []
+ for igrp, (fgrp, tgrp) in enumerate(zip(from_discr.groups, to_discr.groups)):
+ all_elements = cl.array.arange(queue,
+ fgrp.nelements,
+ dtype=np.intp).with_queue(None)
+ ibatch = InterpolationBatch(
+ from_group_index=igrp,
+ from_element_indices=all_elements,
+ to_element_indices=all_elements,
+ result_unit_nodes=tgrp.unit_nodes,
+ to_element_face=None)
+
+ groups.append(
+ DiscretizationConnectionElementGroup([ibatch]))
+
+ return DiscretizationConnection(
+ from_discr, to_discr, groups,
+ is_surjective=True)
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py
index d59346691c93f7e967a1174b0a4622fdfff05426..a5319bc84f47ab6689f502629ff6009893033832 100644
--- a/meshmode/discretization/poly_element.py
+++ b/meshmode/discretization/poly_element.py
@@ -38,6 +38,7 @@ Group types
.. autoclass:: InterpolatoryQuadratureSimplexElementGroup
.. autoclass:: QuadratureSimplexElementGroup
.. autoclass:: PolynomialWarpAndBlendElementGroup
+.. autoclass:: PolynomialEquidistantElementGroup
Group factories
^^^^^^^^^^^^^^^
@@ -45,6 +46,7 @@ Group factories
.. autoclass:: InterpolatoryQuadratureSimplexGroupFactory
.. autoclass:: QuadratureSimplexGroupFactory
.. autoclass:: PolynomialWarpAndBlendGroupFactory
+.. autoclass:: PolynomialEquidistantGroupFactory
"""
# FIXME Most of the loopy kernels will break as soon as we start using multiple
@@ -63,13 +65,16 @@ from meshmode.discretization import ElementGroupBase
class PolynomialSimplexElementGroupBase(ElementGroupBase):
def basis(self):
- return mp.simplex_onb(self.dim, self.order)
+ return mp.simplex_best_available_basis(self.dim, self.order)
+
+ def grad_basis(self):
+ return mp.grad_simplex_best_available_basis(self.dim, self.order)
@memoize_method
def from_mesh_interp_matrix(self):
meg = self.mesh_el_group
return mp.resampling_matrix(
- mp.simplex_onb(meg.dim, meg.order),
+ mp.simplex_best_available_basis(meg.dim, meg.order),
self.unit_nodes,
meg.unit_nodes)
@@ -77,7 +82,7 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase):
def diff_matrices(self):
result = mp.differentiation_matrices(
self.basis(),
- mp.grad_simplex_onb(self.dim, self.order),
+ self.grad_basis(),
self.unit_nodes)
if not isinstance(result, tuple):
@@ -85,13 +90,6 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase):
else:
return result
- @memoize_method
- def resampling_matrix(self):
- meg = self.mesh_el_group
- return mp.resampling_matrix(
- mp.simplex_onb(self.dim, meg.order),
- self.unit_nodes, meg.unit_nodes)
-
class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
"""Elemental discretization supplying a high-order quadrature rule
@@ -154,11 +152,20 @@ class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
return self._quadrature_rule().weights
-class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase):
+class _MassMatrixQuadratureElementGroup(PolynomialSimplexElementGroupBase):
+ @property
+ @memoize_method
+ def weights(self):
+ return np.dot(
+ mp.mass_matrix(self.basis(), self.unit_nodes),
+ np.ones(len(self.basis())))
+
+
+class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup):
"""Elemental discretization with a number of nodes matching the number of
polynomials in :math:`P^k`, hence usable for differentiation and
- interpolation. Interpolation nodes are present on the boundary of the
- simplex.
+ interpolation. Interpolation nodes edge-clustered for avoidance of Runge
+ phenomena. Nodes are present on the boundary of the simplex.
"""
@property
@memoize_method
@@ -170,12 +177,24 @@ class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase):
assert dim2 == dim
return result
+
+class PolynomialEquidistantElementGroup(_MassMatrixQuadratureElementGroup):
+ """Elemental discretization with a number of nodes matching the number of
+ polynomials in :math:`P^k`, hence usable for differentiation and
+ interpolation. Interpolation nodes are present on the boundary of the
+ simplex.
+
+ .. versionadded:: 2016.1
+ """
@property
@memoize_method
- def weights(self):
- return np.dot(
- mp.mass_matrix(self.basis(), self.unit_nodes),
- np.ones(len(self.basis())))
+ def unit_nodes(self):
+ dim = self.mesh_el_group.dim
+ result = mp.equidistant_nodes(dim, self.order)
+
+ dim2, nunit_nodes = result.shape
+ assert dim2 == dim
+ return result
# }}}
@@ -212,6 +231,15 @@ class PolynomialWarpAndBlendGroupFactory(OrderBasedGroupFactory):
mesh_group_class = _MeshSimplexElementGroup
group_class = PolynomialWarpAndBlendElementGroup
+
+class PolynomialEquidistantGroupFactory(OrderBasedGroupFactory):
+ """
+ .. versionadded:: 2016.1
+ """
+
+ mesh_group_class = _MeshSimplexElementGroup
+ group_class = PolynomialEquidistantElementGroup
+
# }}}
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index 3c2565fdb5a7ac206d76b0c2fb5ed44853939e54..d68ec0c58078e68417721c8205fc3f24021d0bad 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -34,7 +34,7 @@ __doc__ = """
.. autoclass:: Visualizer
-.. autofunction:: write_mesh_connectivity_vtk_file
+.. autofunction:: write_nodal_adjacency_vtk_file
"""
@@ -77,10 +77,10 @@ class Visualizer(object):
.. automethod:: write_vtk_file
"""
- def __init__(self, discr, vis_discr, connection):
- self.discr = discr
- self.vis_discr = vis_discr
+ def __init__(self, connection):
self.connection = connection
+ self.discr = connection.from_discr
+ self.vis_discr = connection.to_discr
def _resample_and_get(self, queue, vec):
from pytools.obj_array import with_object_array_or_scalar
@@ -240,16 +240,36 @@ def make_visualizer(queue, discr, vis_order):
real_dtype=discr.real_dtype)
from meshmode.discretization.connection import \
make_same_mesh_connection
- cnx = make_same_mesh_connection(queue, vis_discr, discr)
- return Visualizer(discr, vis_discr, cnx)
+ return Visualizer(make_same_mesh_connection(vis_discr, discr))
# }}}
-# {{{ connectivity
+# {{{ draw_curve
-def write_mesh_connectivity_vtk_file(file_name, mesh, compressor=None,):
+def draw_curve(discr):
+ mesh = discr.mesh
+
+ import matplotlib.pyplot as pt
+ pt.plot(mesh.vertices[0], mesh.vertices[1], "o")
+
+ color = pt.cm.rainbow(np.linspace(0, 1, len(discr.groups)))
+ with cl.CommandQueue(discr.cl_context) as queue:
+ for i, group in enumerate(discr.groups):
+ group_nodes = group.view(discr.nodes()).get(queue=queue)
+ pt.plot(
+ group_nodes[0].T,
+ group_nodes[1].T, "-x",
+ label="Group %d" % i,
+ color=color[i])
+
+# }}}
+
+
+# {{{ adjacency
+
+def write_nodal_adjacency_vtk_file(file_name, mesh, compressor=None,):
from pyvisfile.vtk import (
UnstructuredGrid, DataArray,
AppendedDataXMLGenerator,
@@ -266,16 +286,16 @@ def write_mesh_connectivity_vtk_file(file_name, mesh, compressor=None,):
np.sum(mesh.vertices[:, grp.vertex_indices], axis=-1)
/ grp.vertex_indices.shape[-1])
- cnx = mesh.element_connectivity
+ adj = mesh.nodal_adjacency
- nconnections = len(cnx.neighbors)
+ nconnections = len(adj.neighbors)
connections = np.empty((nconnections, 2), dtype=np.int32)
- nb_starts = cnx.neighbors_starts
+ nb_starts = adj.neighbors_starts
for iel in range(mesh.nelements):
connections[nb_starts[iel]:nb_starts[iel+1], 0] = iel
- connections[:, 1] = cnx.neighbors
+ connections[:, 1] = adj.neighbors
grid = UnstructuredGrid(
(mesh.nelements,
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 6bd4b92492de942de1ca0c388207d63d44257973..5c2b5c1cb6ec7be972ed1ca15e0eb2db1c4b852c 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -24,6 +24,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+import six
import numpy as np
import modepy as mp
import numpy.linalg as la
@@ -39,11 +40,59 @@ __doc__ = """
:members:
:undoc-members:
-.. autoclass:: ElementConnectivity
+.. autoclass:: NodalAdjacency
+.. autoclass:: FacialAdjacencyGroup
+.. autofunction:: as_python
+.. autofunction:: check_bc_coverage
+.. autofunction:: is_boundary_tag_empty
+
+Predefined Boundary tags
+------------------------
+
+.. autoclass:: BTAG_NONE
+.. autoclass:: BTAG_ALL
+.. autoclass:: BTAG_REALLY_ALL
+.. autoclass:: BTAG_NO_BOUNDARY
"""
+# {{{ element tags
+
+class BTAG_NONE(object):
+ """A boundary tag representing an empty boundary or volume."""
+ pass
+
+
+class BTAG_ALL(object):
+ """A boundary tag representing the entire boundary or volume.
+
+ In the case of the boundary, TAG_ALL does not include rank boundaries,
+ or, more generally, anything tagged with TAG_NO_BOUNDARY."""
+ pass
+
+
+class BTAG_REALLY_ALL(object):
+ """A boundary tag representing the entire boundary.
+
+ Unlike :class:`TAG_ALL`, this includes rank boundaries,
+ or, more generally, everything tagged with :class:`TAG_NO_BOUNDARY`."""
+ pass
+
+
+class BTAG_NO_BOUNDARY(object):
+ """A boundary tag indicating that this edge should not fall under
+ :class:`TAG_ALL`. Among other things, this is used to keep rank boundaries
+ out of :class:`BTAG_ALL`.
+ """
+ pass
+
+
+SYSTEM_TAGS = set([BTAG_NONE, BTAG_ALL, BTAG_REALLY_ALL, BTAG_NO_BOUNDARY])
+
+# }}}
+
+
# {{{ element group
class MeshElementGroup(Record):
@@ -78,6 +127,14 @@ class MeshElementGroup(Record):
The number of dimensions spanned by the element.
*Not* the ambient dimension, see :attr:`Mesh.ambient_dim`
for that.
+
+ .. automethod:: face_vertex_indices
+ .. automethod:: vertex_unit_coordinates
+
+ .. attribute:: nfaces
+
+ .. automethod:: __eq__
+ .. automethod:: __ne__
"""
def __init__(self, order, vertex_indices, nodes,
@@ -134,6 +191,36 @@ class MeshElementGroup(Record):
def nunit_nodes(self):
return self.unit_nodes.shape[-1]
+ def face_vertex_indices(self):
+ """Return a tuple of tuples indicating which vertices
+ (in mathematically positive ordering) make up each face
+ of an element in this group.
+ """
+ raise NotImplementedError()
+
+ def vertex_unit_coordinates(self):
+ """Return an array of shape ``(nfaces, dim)`` with the unit
+ coordinates of each vertex.
+ """
+ raise NotImplementedError()
+
+ @property
+ def nfaces(self):
+ return len(self.face_vertex_indices())
+
+ def __eq__(self, other):
+ return (
+ type(self) == type(other)
+ and self.order == other.order
+ and np.array_equal(self.vertex_indices, other.vertex_indices)
+ and np.array_equal(self.nodes, other.nodes)
+ and np.array_equal(self.unit_nodes, other.unit_nodes)
+ and self.element_nr_base == other.element_nr_base
+ and self.node_nr_base == other.node_nr_base)
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
+
class SimplexElementGroup(MeshElementGroup):
def __init__(self, order, vertex_indices, nodes,
@@ -164,7 +251,10 @@ class SimplexElementGroup(MeshElementGroup):
raise TypeError("'dim' must be passed "
"if 'unit_nodes' is not passed")
- unit_nodes = mp.warp_and_blend_nodes(dim, order)
+ if dim <= 3:
+ unit_nodes = mp.warp_and_blend_nodes(dim, order)
+ else:
+ unit_nodes = mp.equidistant_nodes(dim, order)
dims = unit_nodes.shape[0]
@@ -173,12 +263,15 @@ class SimplexElementGroup(MeshElementGroup):
"element. expected: %d, got: %d" % (dims+1,
vertex_indices.shape[-1]))
- MeshElementGroup.__init__(self, order, vertex_indices, nodes,
+ super(SimplexElementGroup, self).__init__(order, vertex_indices, nodes,
element_nr_base, node_nr_base, unit_nodes, dim)
def face_vertex_indices(self):
if self.dim == 1:
- return ((0, 1),)
+ return (
+ (0,),
+ (1,),
+ )
elif self.dim == 2:
return (
(0, 1),
@@ -196,36 +289,21 @@ class SimplexElementGroup(MeshElementGroup):
raise NotImplementedError("dim=%d" % self.dim)
def vertex_unit_coordinates(self):
- if self.dim == 1:
- return np.array([
- [-1], [1]
- ], dtype=np.float64)
- elif self.dim == 2:
- return np.array([
- [-1, -1],
- [1, -1],
- [-1, 1],
- ], dtype=np.float64)
- elif self.dim == 3:
- return np.array([
- [-1, -1, -1],
- [1, -1, -1],
- [-1, 1, -1],
- [-1, -1, 1],
- ], dtype=np.float64)
- else:
- raise NotImplementedError("dim=%d" % self.dim)
+ from modepy.tools import unit_vertices
+ return unit_vertices(self.dim)
# }}}
-# {{{ mesh
+# {{{ nodal adjacency
+
+class NodalAdjacency(Record):
+ """Describes nodal element adjacency information, i.e. information about
+ elements that touch in at least one point.
-class ElementConnectivity(Record):
- """
.. attribute:: neighbors_starts
- ``element_id_t [nelments+1]``
+ ``element_id_t [nelements+1]``
Use together with :attr:`neighbors`. ``neighbors_starts[iel]`` and
``neighbors_starts[iel+1]`` together indicate a ranges of element indices
@@ -236,8 +314,92 @@ class ElementConnectivity(Record):
``element_id_t []``
See :attr:`neighbors_starts`.
+
+ .. automethod:: __eq__
+ .. automethod:: __ne__
"""
+ def __eq__(self, other):
+ return (
+ type(self) == type(other)
+ and np.array_equal(self.neighbors_starts,
+ other.neighbors_starts)
+ and np.array_equal(self.neighbors, other.neighbors))
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
+
+# }}}
+
+
+# {{{ facial adjacency
+
+class FacialAdjacencyGroup(Record):
+ """Describes facial element adjacency information for one
+ :class:`MeshElementGroup`, i.e. information about elements that share (part
+ of) a face.
+
+ .. attribute:: igroup
+
+ .. attribute:: ineighbor_group
+
+ ID of neighboring group, or *None* for boundary faces. If identical
+ to :attr:`igroup`, then this contains the self-connectivity in this
+ group.
+
+ .. attribute:: elements
+
+ ``element_id_t [nfagrp_elements]``. ``elements[i]``
+ Group-local element numbers.
+
+ .. attribute:: element_faces
+
+ ``face_id_t [nfagrp_elements]``. ``element_faces[i]``
+ indicate what face index of the opposite element indicated in
+ ``neighbors[iel_grp][iface]`` touches face number *iface* of element
+ number *iel_grp* in this element group.
+
+ .. attribute:: neighbors
+
+ ``element_id_t [nfagrp_elements]``. ``neighbors[i]``
+ gives the element number within :attr:`ineighbor_group` of the element
+ opposite ``elements[i]``.
+
+ If this number is negative, then this indicates that this is a
+ boundary face, and the bits set in ``-neighbors[i]``
+ should be interpreted according to :attr:`Mesh.boundary_tags`.
+
+ .. attribute:: neighbor_faces
+
+ ``face_id_t [nfagrp_elements]``. ``neighbor_faces[i]`` indicate what
+ face index of the opposite element indicated in ``neighbors[i]`` has
+ facial contact with face number ``element_faces[i]`` of element number
+ ``elements[i]`` in element group *igroup*.
+
+ Zero if ``neighbors[i]`` is negative.
+
+ .. automethod:: __eq__
+ .. automethod:: __ne__
+ """
+
+ def __eq__(self, other):
+ return (
+ type(self) == type(other)
+ and self.igroup == other.igroup
+ and self.ineighbor_group == other.ineighbor_group
+ and np.array_equal(self.elements, other.elements)
+ and np.array_equal(self.element_faces, other.element_faces)
+ and np.array_equal(self.neighbors, other.neighbors)
+ and np.array_equal(self.neighbor_faces, other.neighbor_faces)
+ )
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
+
+# }}}
+
+
+# {{{ mesh
class Mesh(Record):
"""
@@ -250,20 +412,51 @@ class Mesh(Record):
A list of :class:`MeshElementGroup` instances.
- .. attribute:: element_connectivity
+ .. attribute:: nodal_adjacency
- An instance of :class:`ElementConnectivity`.
+ An instance of :class:`NodalAdjacency`.
Referencing this attribute may raise
- :exc:`meshmode.ConnectivityUnavailable`.
+ :exc:`meshmode.DataUnavailable`.
+
+ .. attribute:: facial_adjacency_groups
+
+ A list of mappings from neighbor group IDs to instances of
+ :class:`FacialAdjacencyGroup`.
+
+ ``facial_adjacency_groups[igrp][ineighbor_group]`` gives
+ the set of facial adjacency relations between group *igrp*
+ and *ineighbor_group*. *ineighbor_group* and *igrp* may be
+ identical, or *ineighbor_group* may be *None*, in which case
+ a group containing boundary faces is returned.
+
+ Referencing this attribute may raise
+ :exc:`meshmode.DataUnavailable`.
+
+ .. attribute:: boundary_tags
+
+ A tuple of boundary tag identifiers. :class:`BTAG_ALL` and
+ :class:`BTAG_REALLY_ALL` are guranateed to exist.
+
+ .. attribute:: btag_to_index
+
+ A mapping that maps boundary tag identifiers to their
+ corresponding index.
.. attribute:: vertex_id_dtype
.. attribute:: element_id_dtype
+
+ .. automethod:: __eq__
+ .. automethod:: __ne__
"""
+ face_id_dtype = np.int8
+
def __init__(self, vertices, groups, skip_tests=False,
- element_connectivity=False,
+ nodal_adjacency=False,
+ facial_adjacency_groups=False,
+ boundary_tags=None,
vertex_id_dtype=np.int32,
element_id_dtype=np.int32):
"""
@@ -271,7 +464,7 @@ class Mesh(Record):
:arg skip_tests: Skip mesh tests, in case you want to load a broken
mesh anyhow and then fix it inside of this data structure.
- :arg element_connectivity: One of three options:
+ :arg nodal_adjacency: One of three options:
*None*, in which case this information
will be deduced from vertex adjacency. *False*, in which case
this information will be marked unavailable (such as if there are
@@ -279,8 +472,16 @@ class Mesh(Record):
the full picture), and references to
:attr:`element_neighbors_starts` and :attr:`element_neighbors`
will result in exceptions. Lastly, a tuple
- *(element_neighbors_starts, element_neighbors)*, representing the
- correspondingly-named attributes.
+ :class:`NodalAdjacency` object.
+ :arg facial_adjacency_groups: One of three options:
+ *None*, in which case this information
+ will be deduced from vertex adjacency. *False*, in which case
+ this information will be marked unavailable (such as if there are
+ hanging nodes in the geometry, so that vertex adjacency does not convey
+ the full picture), and references to
+ :attr:`element_neighbors_starts` and :attr:`element_neighbors`
+ will result in exceptions. Lastly, a data structure as described in
+ :attr:`facial_adjacency_groups` may be passed.
"""
el_nr = 0
node_nr = 0
@@ -292,18 +493,48 @@ class Mesh(Record):
el_nr += ng.nelements
node_nr += ng.nnodes
- if element_connectivity is not False and element_connectivity is not None:
- nb_starts, nbs = element_connectivity
- element_connectivity = ElementConnectivity(
- neighbors_starts=nb_starts,
- neighbors=nbs)
+ # {{{ boundary tags
+
+ if boundary_tags is None:
+ boundary_tags = []
+ else:
+ boundary_tags = boundary_tags[:]
+
+ if BTAG_NONE in boundary_tags:
+ raise ValueError("BTAG_NONE is not allowed to be part of "
+ "boundary_tags")
+ if BTAG_ALL not in boundary_tags:
+ boundary_tags.append(BTAG_ALL)
+ if BTAG_REALLY_ALL not in boundary_tags:
+ boundary_tags.append(BTAG_REALLY_ALL)
+
+ max_boundary_tag_count = int(
+ np.log(np.iinfo(element_id_dtype).max)/np.log(2))
+ if len(boundary_tags) > max_boundary_tag_count:
+ raise ValueError("too few bits in element_id_dtype to represent all "
+ "boundary tags")
+
+ btag_to_index = dict(
+ (btag, i) for i, btag in enumerate(boundary_tags))
+
+ # }}}
+
+ if nodal_adjacency is not False and nodal_adjacency is not None:
+ if not isinstance(nodal_adjacency, NodalAdjacency):
+ nb_starts, nbs = nodal_adjacency
+ nodal_adjacency = NodalAdjacency(
+ neighbors_starts=nb_starts,
+ neighbors=nbs)
- del nb_starts
- del nbs
+ del nb_starts
+ del nbs
Record.__init__(
self, vertices=vertices, groups=new_groups,
- _element_connectivity=element_connectivity,
+ _nodal_adjacency=nodal_adjacency,
+ _facial_adjacency_groups=facial_adjacency_groups,
+ boundary_tags=boundary_tags,
+ btag_to_index=btag_to_index,
vertex_id_dtype=np.dtype(vertex_id_dtype),
element_id_dtype=np.dtype(element_id_dtype),
)
@@ -313,6 +544,32 @@ class Mesh(Record):
for g in self.groups:
assert g.vertex_indices.dtype == self.vertex_id_dtype
+ if nodal_adjacency:
+ assert nodal_adjacency.neighbor_starts.shape == (self.nelements+1,)
+ assert len(nodal_adjacency.neighbors.shape) == 1
+
+ assert nodal_adjacency.neighbor_starts.dtype == self.element_id_dtype
+ assert nodal_adjacency.neighbors.dtype == self.element_id_dtype
+
+ if facial_adjacency_groups:
+ assert len(facial_adjacency_groups) == len(self.groups)
+ for fagrp_map in facial_adjacency_groups:
+ for fagrp in six.itervalues(fagrp_map):
+ grp = self.groups[fagrp.igroup]
+
+ fvi = grp.face_vertex_indices()
+ assert fagrp.neighbors.dtype == self.element_id_dtype
+ assert fagrp.neighbors.shape == (
+ grp.nelements, len(fvi))
+ assert fagrp.neighbors.dtype == self.face_id_dtype
+ assert fagrp.neighbor_faces.shape == (
+ grp.nelements, len(fvi))
+
+ is_bdry = fagrp.neighbors < 0
+ assert ((1 << btag_to_index[BTAG_REALLY_ALL])
+ & fagrp.neighbors[is_bdry]).all(), \
+ "boundary faces without BTAG_REALLY_ALL found"
+
from meshmode.mesh.processing import \
test_volume_mesh_element_orientations
@@ -330,19 +587,63 @@ class Mesh(Record):
from pytools import single_valued
return single_valued(grp.dim for grp in self.groups)
+ @property
+ def nvertices(self):
+ return self.vertices.shape[-1]
+
@property
def nelements(self):
return sum(grp.nelements for grp in self.groups)
@property
- def element_connectivity(self):
- if self._element_connectivity is False:
- from meshmode import ConnectivityUnavailable
- raise ConnectivityUnavailable()
- elif self._element_connectivity is None:
- self._element_connectivity = _compute_connectivity_from_vertices(self)
+ def nodal_adjacency(self):
+ if self._nodal_adjacency is False:
+ from meshmode import DataUnavailable
+ raise DataUnavailable("nodal_adjacency")
+ elif self._nodal_adjacency is None:
+ self._nodal_adjacency = _compute_nodal_adjacency_from_vertices(self)
+
+ return self._nodal_adjacency
+
+ def nodal_adjacency_init_arg(self):
+ """Returns an 'nodal_adjacency' argument that can be
+ passed to a Mesh constructor.
+ """
- return self._element_connectivity
+ return self._nodal_adjacency
+
+ @property
+ def facial_adjacency_groups(self):
+ if self._facial_adjacency_groups is False:
+ from meshmode import DataUnavailable
+ raise DataUnavailable("facial_adjacency_groups")
+ elif self._facial_adjacency_groups is None:
+ self._facial_adjacency_groups = \
+ _compute_facial_adjacency_from_vertices(self)
+
+ return self._facial_adjacency_groups
+
+ def boundary_tag_bit(self, boundary_tag):
+ try:
+ return 1 << self.btag_to_index[boundary_tag]
+ except KeyError:
+ return 0
+
+ def __eq__(self, other):
+ return (
+ type(self) == type(other)
+ and np.array_equal(self.vertices, other.vertices)
+ and self.groups == other.groups
+ and self.vertex_id_dtype == other.vertex_id_dtype
+ and self.element_id_dtype == other.element_id_dtype
+ and (self._nodal_adjacency
+ == other._nodal_adjacency)
+ and (self._facial_adjacency_groups
+ == other._facial_adjacency_groups)
+ and self.boundary_tags == other.boundary_tags)
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
# Design experience: Try not to add too many global data structures to the
# mesh. Let the element groups be responsible for that at the mesh level.
@@ -355,10 +656,13 @@ class Mesh(Record):
# {{{ node-vertex consistency test
def _test_node_vertex_consistency_simplex(mesh, mgrp):
- from modepy.tools import UNIT_VERTICES
+ if mgrp.nelements == 0:
+ return True
+
resampling_mat = mp.resampling_matrix(
- mp.simplex_onb(mgrp.dim, mgrp.order),
- UNIT_VERTICES[mgrp.dim].T.copy(), mgrp.unit_nodes)
+ mp.simplex_best_available_basis(mgrp.dim, mgrp.order),
+ mgrp.vertex_unit_coordinates().T,
+ mgrp.unit_nodes)
# dim, nelments, nnvertices
map_vertices = np.einsum(
@@ -401,9 +705,9 @@ def _test_node_vertex_consistency(mesh):
# }}}
-# {{{ vertex-based connectivity
+# {{{ vertex-based nodal adjacency
-def _compute_connectivity_from_vertices(mesh):
+def _compute_nodal_adjacency_from_vertices(mesh):
# FIXME Native code would make this faster
_, nvertices = mesh.vertices.shape
@@ -436,10 +740,296 @@ def _compute_connectivity_from_vertices(mesh):
assert neighbors_starts[-1] == len(neighbors)
- return ElementConnectivity(
+ return NodalAdjacency(
neighbors_starts=neighbors_starts,
neighbors=neighbors)
# }}}
+
+# {{{ vertex-based facial adjacency
+
+def _compute_facial_adjacency_from_vertices(mesh):
+ # FIXME Native code would make this faster
+
+ # create face_map, which is a mapping of
+ # (vertices on a face) ->
+ # [(igrp, iel_grp, face_idx) for elements bordering that face]
+ face_map = {}
+ for igrp, grp in enumerate(mesh.groups):
+ for fid, face_vertex_indices in enumerate(grp.face_vertex_indices()):
+ all_fvi = grp.vertex_indices[:, face_vertex_indices]
+
+ for iel_grp, fvi in enumerate(all_fvi):
+ face_map.setdefault(
+ frozenset(fvi), []).append((igrp, iel_grp, fid))
+
+ # maps tuples (igrp, ineighbor_group) to number of elements
+ group_count = {}
+ for face_tuples in six.itervalues(face_map):
+ if len(face_tuples) == 2:
+ (igrp, _, _), (inb_grp, _, _) = face_tuples
+ group_count[igrp, inb_grp] = group_count.get((igrp, inb_grp), 0) + 1
+ group_count[inb_grp, igrp] = group_count.get((inb_grp, igrp), 0) + 1
+ elif len(face_tuples) == 1:
+ (igrp, _, _), = face_tuples
+ group_count[igrp, None] = group_count.get((igrp, None), 0) + 1
+ else:
+ raise RuntimeError("unexpected number of adjacent faces")
+
+ # {{{ build facial_adjacency_groups data structure, still empty
+
+ facial_adjacency_groups = []
+ for igroup in range(len(mesh.groups)):
+ grp_map = {}
+ facial_adjacency_groups.append(grp_map)
+
+ bdry_count = group_count.get((igroup, None))
+ if bdry_count is not None:
+ elements = np.empty(bdry_count, dtype=mesh.element_id_dtype)
+ element_faces = np.empty(bdry_count, dtype=mesh.face_id_dtype)
+ neighbors = np.empty(bdry_count, dtype=mesh.element_id_dtype)
+ neighbor_faces = np.zeros(bdry_count, dtype=mesh.face_id_dtype)
+
+ neighbors.fill(-(
+ mesh.boundary_tag_bit(BTAG_ALL)
+ | mesh.boundary_tag_bit(BTAG_REALLY_ALL)))
+
+ grp_map[None] = FacialAdjacencyGroup(
+ igroup=igroup,
+ ineighbor_group=None,
+ elements=elements,
+ element_faces=element_faces,
+ neighbors=neighbors,
+ neighbor_faces=neighbor_faces)
+
+ for ineighbor_group in range(len(mesh.groups)):
+ nb_count = group_count.get((igroup, ineighbor_group))
+ if nb_count is not None:
+ elements = np.empty(nb_count, dtype=mesh.element_id_dtype)
+ element_faces = np.empty(nb_count, dtype=mesh.face_id_dtype)
+ neighbors = np.empty(nb_count, dtype=mesh.element_id_dtype)
+ neighbor_faces = np.empty(nb_count, dtype=mesh.face_id_dtype)
+
+ grp_map[ineighbor_group] = FacialAdjacencyGroup(
+ igroup=igroup,
+ ineighbor_group=ineighbor_group,
+ elements=elements,
+ element_faces=element_faces,
+ neighbors=neighbors,
+ neighbor_faces=neighbor_faces)
+
+ # }}}
+
+ # maps tuples (igrp, ineighbor_group) to number of elements filled in group
+ fill_count = {}
+ for face_tuples in six.itervalues(face_map):
+ if len(face_tuples) == 2:
+ for (igroup, iel, iface), (inb_group, inb_el, inb_face) in [
+ (face_tuples[0], face_tuples[1]),
+ (face_tuples[1], face_tuples[0]),
+ ]:
+ idx = fill_count.get((igrp, inb_grp), 0)
+ fill_count[igrp, inb_grp] = idx + 1
+
+ fagrp = facial_adjacency_groups[igroup][inb_grp]
+ fagrp.elements[idx] = iel
+ fagrp.element_faces[idx] = iface
+ fagrp.neighbors[idx] = inb_el
+ fagrp.neighbor_faces[idx] = inb_face
+
+ elif len(face_tuples) == 1:
+ (igroup, iel, iface), = face_tuples
+
+ idx = fill_count.get((igrp, None), 0)
+ fill_count[igrp, None] = idx + 1
+
+ fagrp = facial_adjacency_groups[igroup][None]
+ fagrp.elements[idx] = iel
+ fagrp.element_faces[idx] = iface
+
+ else:
+ raise RuntimeError("unexpected number of adjacent faces")
+
+ return facial_adjacency_groups
+
+# }}}
+
+
+# {{{ as_python
+
+def _numpy_array_as_python(array):
+ return "np.array(%s, dtype=np.%s)" % (
+ repr(array.tolist()),
+ array.dtype.name)
+
+
+def as_python(mesh, function_name="make_mesh"):
+ """Return a snippet of Python code (as a string) that will
+ recreate the mesh given as an input parameter.
+ """
+
+ from pytools.py_codegen import PythonCodeGenerator, Indentation
+ cg = PythonCodeGenerator()
+ cg("""
+ # generated by meshmode.mesh.as_python
+
+ import numpy as np
+ from meshmode.mesh import (
+ Mesh, MeshElementGroup, FacialAdjacencyGroup,
+ BTAG_ALL, BTAG_REALLY_ALL)
+
+ """)
+
+ cg("def %s():" % function_name)
+ with Indentation(cg):
+ cg("vertices = " + _numpy_array_as_python(mesh.vertices))
+ cg("")
+ cg("groups = []")
+ cg("")
+ for group in mesh.groups:
+ cg("import %s" % type(group).__module__)
+ cg("groups.append(%s.%s(" % (
+ type(group).__module__,
+ type(group).__name__))
+ cg(" order=%s," % group.order)
+ cg(" vertex_indices=%s,"
+ % _numpy_array_as_python(group.vertex_indices))
+ cg(" nodes=%s,"
+ % _numpy_array_as_python(group.nodes))
+ cg(" unit_nodes=%s))"
+ % _numpy_array_as_python(group.unit_nodes))
+
+ # {{{ facial adjacency groups
+
+ def fagrp_params_str(fagrp):
+ params = {
+ "igroup": fagrp.igroup,
+ "ineighbor_group": repr(fagrp.ineighbor_group),
+ "elements": _numpy_array_as_python(fagrp.elements),
+ "element_faces": _numpy_array_as_python(fagrp.element_faces),
+ "neighbors": _numpy_array_as_python(fagrp.neighbors),
+ "neighbor_faces": _numpy_array_as_python(fagrp.neighbor_faces),
+ }
+ return ",\n ".join("%s=%s" % (k, v) for k, v in six.iteritems(params))
+
+ if mesh._facial_adjacency_groups:
+ cg("facial_adjacency_groups = []")
+
+ for igrp, fagrp_map in enumerate(mesh.facial_adjacency_groups):
+ cg("facial_adjacency_groups.append({%s})" % ",\n ".join(
+ "%r: FacialAdjacencyGroup(%s)" % (
+ inb_grp, fagrp_params_str(fagrp))
+ for inb_grp, fagrp in six.iteritems(fagrp_map)))
+
+ else:
+ cg("facial_adjacency_groups = %r" % mesh._facial_adjacency_groups)
+
+ # }}}
+
+ # {{{ boundary tags
+
+ def strify_boundary_tag(btag):
+ if isinstance(btag, type):
+ return btag.__name__
+ else:
+ return repr(btag)
+
+ btags_str = ", ".join(
+ strify_boundary_tag(btag) for btag in mesh.boundary_tags)
+
+ # }}}
+
+ cg("return Mesh(vertices, groups, skip_tests=True,")
+ cg(" vertex_id_dtype=np.%s," % mesh.vertex_id_dtype.name)
+ cg(" element_id_dtype=np.%s," % mesh.element_id_dtype.name)
+
+ if isinstance(mesh._nodal_adjacency, NodalAdjacency):
+ el_con_str = "(%s, %s)" % (
+ _numpy_array_as_python(
+ mesh._nodal_adjacency.neighbors_starts),
+ _numpy_array_as_python(
+ mesh._nodal_adjacency.neighbors),
+ )
+ else:
+ el_con_str = repr(mesh._nodal_adjacency)
+
+ cg(" nodal_adjacency=%s," % el_con_str)
+ cg(" facial_adjacency_groups=facial_adjacency_groups,")
+ cg(" boundary_tags=[%s])" % btags_str)
+
+ # FIXME: Handle facial adjacency, boundary tags
+
+ return cg.get()
+
+# }}}
+
+
+# {{{ check_bc_coverage
+
+def check_bc_coverage(mesh, boundary_tags, incomplete_ok=False):
+ """Verify boundary condition coverage.
+
+ Given a list of boundary tags as *boundary_tags*, this function verifies
+ that
+
+ 1. the union of all these boundaries gives the complete boundary,
+ 1. all these boundaries are disjoint.
+
+ :arg incomplete_ok: Do not report an error if some faces are not covered
+ by the boundary conditions.
+ """
+
+ for igrp, fagrp_map in enumerate(mesh.facial_adjacency_groups):
+ bdry_grp = fagrp_map.get(None)
+ if bdry_grp is None:
+ continue
+
+ nb_elements = bdry_grp.neighbors
+ assert (nb_elements < 0).all()
+
+ nb_el_bits = -nb_elements
+
+ seen = np.zeros_like(nb_el_bits, dtype=np.bool)
+
+ for btag in boundary_tags:
+ tag_bit = mesh.boundary_tag_bit(btag)
+ tag_set = (nb_el_bits & tag_bit) != 0
+
+ if (seen & tag_set).any():
+ raise RuntimeError("faces with multiple boundary conditions found")
+
+ seen = seen | tag_set
+
+ if not incomplete_ok and not seen.all():
+ raise RuntimeError("found faces without boundary conditions")
+
+# }}}
+
+
+# {{{ is_boundary_tag_empty
+
+def is_boundary_tag_empty(mesh, boundary_tag):
+ """Return *True* if the corresponding boundary tag does not occur as part of
+ *mesh*.
+ """
+
+ btag_bit = mesh.boundary_tag_bit(boundary_tag)
+ if not btag_bit:
+ return True
+
+ for igrp in range(len(mesh.groups)):
+ bdry_fagrp = mesh.facial_adjacency_groups[igrp].get(None, None)
+ if bdry_fagrp is None:
+ continue
+
+ neg = bdry_fagrp.neighbors < 0
+ if (-bdry_fagrp.neighbors[neg] & btag_bit).any():
+ return False
+
+ return True
+
+# }}}
+
+
# vim: foldmethod=marker
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 5b925dcb62dfacc73b90b1180800a737b1ba0a4d..65761586e01c8b0aa60f82d6f8b4bc00941eeed7 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -1,6 +1,4 @@
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
+from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
@@ -24,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+from six.moves import range
import numpy as np
import numpy.linalg as la
@@ -53,7 +52,13 @@ Surfaces
.. autofunction:: generate_icosahedron
.. autofunction:: generate_icosphere
.. autofunction:: generate_torus
+
+Volumes
+-------
+
+.. autofunction:: generate_box_mesh
.. autofunction:: generate_regular_rect_mesh
+.. autofunction:: generate_warped_rect_mesh
"""
@@ -207,8 +212,10 @@ def make_curve_mesh(curve_f, element_boundaries, order):
nodes=nodes,
unit_nodes=unodes)
- return Mesh(vertices=vertices, groups=[egroup],
- element_connectivity=None)
+ return Mesh(
+ vertices=vertices, groups=[egroup],
+ nodal_adjacency=None,
+ facial_adjacency_groups=None)
# }}}
@@ -227,7 +234,11 @@ def make_group_from_vertices(vertices, vertex_indices, order):
dim = nspan_vectors
# dim, nunit_nodes
- unit_nodes = mp.warp_and_blend_nodes(dim, order)
+ if dim <= 3:
+ unit_nodes = mp.warp_and_blend_nodes(dim, order)
+ else:
+ unit_nodes = mp.equidistant_nodes(dim, order)
+
unit_nodes_01 = 0.5 + 0.5*unit_nodes
nodes = np.einsum(
@@ -280,8 +291,10 @@ def generate_icosahedron(r, order):
grp = make_group_from_vertices(vertices, vertex_indices, order)
from meshmode.mesh import Mesh
- return Mesh(vertices, [grp],
- element_connectivity=None)
+ return Mesh(
+ vertices, [grp],
+ nodal_adjacency=None,
+ facial_adjacency_groups=None)
# }}}
@@ -297,8 +310,10 @@ def generate_icosphere(r, order):
nodes=grp.nodes * r / np.sqrt(np.sum(grp.nodes**2, axis=0)))
from meshmode.mesh import Mesh
- return Mesh(mesh.vertices, [grp],
- element_connectivity=None)
+ return Mesh(
+ mesh.vertices, [grp],
+ nodal_adjacency=None,
+ facial_adjacency_groups=None)
# }}}
@@ -360,7 +375,11 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner,
nodes[2] = b*np.sin(minor_theta)
from meshmode.mesh import Mesh
- return (Mesh(vertices, [grp.copy(nodes=nodes)], element_connectivity=None),
+ return (
+ Mesh(
+ vertices, [grp.copy(nodes=nodes)],
+ nodal_adjacency=None,
+ facial_adjacency_groups=None),
[idx(i, 0) for i in range(n_outer)],
[idx(0, j) for j in range(n_inner)])
@@ -373,6 +392,104 @@ def generate_torus(r_outer, r_inner, n_outer=20, n_inner=10, order=1):
return mesh
+# {{{ generate_box_mesh
+
+def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64):
+ """Create a semi-structured mesh.
+
+ :param axis_coords: a tuple with a number of entries corresponding
+ to the number of dimensions, with each entry a numpy array
+ specifying the coordinates to be used along that axis.
+ """
+
+ for iaxis, axc in enumerate(axis_coords):
+ if len(axc) < 2:
+ raise ValueError("need at least two points along axis %d"
+ % (iaxis+1))
+
+ dim = len(axis_coords)
+
+ shape = tuple(len(axc) for axc in axis_coords)
+
+ from pytools import product
+ nvertices = product(shape)
+
+ vertex_indices = np.arange(nvertices).reshape(*shape, order="F")
+
+ vertices = np.empty((dim,)+shape, dtype=coord_dtype)
+ for idim in range(dim):
+ vshape = (shape[idim],) + (1,)*idim
+ vertices[idim] = axis_coords[idim].reshape(*vshape)
+
+ vertices = vertices.reshape(dim, -1)
+
+ el_vertices = []
+
+ if dim == 1:
+ for i in range(shape[0]-1):
+ # a--b
+
+ a = vertex_indices[i]
+ b = vertex_indices[i+1]
+
+ el_vertices.append((a, b,))
+
+ elif dim == 2:
+ for i in range(shape[0]-1):
+ for j in range(shape[1]-1):
+
+ # c--d
+ # | |
+ # a--b
+
+ a = vertex_indices[i, j]
+ b = vertex_indices[i+1, j]
+ c = vertex_indices[i, j+1]
+ d = vertex_indices[i+1, j+1]
+
+ el_vertices.append((a, b, c))
+ el_vertices.append((d, c, b))
+
+ elif dim == 3:
+ for i in range(shape[0]-1):
+ for j in range(shape[1]-1):
+ for k in range(shape[2]-1):
+
+ a000 = vertex_indices[i, j, k]
+ a001 = vertex_indices[i, j, k+1]
+ a010 = vertex_indices[i, j+1, k]
+ a011 = vertex_indices[i, j+1, k+1]
+
+ a100 = vertex_indices[i+1, j, k]
+ a101 = vertex_indices[i+1, j, k+1]
+ a110 = vertex_indices[i+1, j+1, k]
+ a111 = vertex_indices[i+1, j+1, k+1]
+
+ el_vertices.append((a000, a100, a010, a001))
+ el_vertices.append((a101, a100, a001, a010))
+ el_vertices.append((a101, a011, a010, a001))
+
+ el_vertices.append((a100, a010, a101, a110))
+ el_vertices.append((a011, a010, a110, a101))
+ el_vertices.append((a011, a111, a101, a110))
+
+ else:
+ raise NotImplementedError("box meshes of dimension %d"
+ % dim)
+
+ el_vertices = np.array(el_vertices, dtype=np.int32)
+
+ grp = make_group_from_vertices(
+ vertices.reshape(dim, -1), el_vertices, order)
+
+ from meshmode.mesh import Mesh
+ return Mesh(vertices, [grp],
+ nodal_adjacency=None,
+ facial_adjacency_groups=None)
+
+# }}}
+
+
# {{{ generate_regular_rect_mesh
def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
@@ -386,42 +503,40 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
if min(n) < 2:
raise ValueError("need at least two points in each direction")
- node_dict = {}
- vertices = []
- points_1d = [np.linspace(a_i, b_i, n_i)
+ axis_coords = [np.linspace(a_i, b_i, n_i)
for a_i, b_i, n_i in zip(a, b, n)]
- for j in range(n[1]):
- for i in range(n[0]):
- node_dict[i, j] = len(vertices)
- vertices.append(np.array([points_1d[0][i], points_1d[1][j]]))
-
- vertices = np.array(vertices).T.copy()
-
- vertex_indices = []
-
- for i in range(n[0]-1):
- for j in range(n[1]-1):
-
- # c--d
- # | |
- # a--b
+ return generate_box_mesh(axis_coords, order=order)
- a = node_dict[i, j]
- b = node_dict[i+1, j]
- c = node_dict[i, j+1]
- d = node_dict[i+1, j+1]
+# }}}
- vertex_indices.append((a, b, c))
- vertex_indices.append((d, c, b))
- vertex_indices = np.array(vertex_indices, dtype=np.int32)
+# {{{ generate_warped_rect_mesh
- grp = make_group_from_vertices(vertices, vertex_indices, order)
+def generate_warped_rect_mesh(dim, order, n):
+ """Generate a mesh of a warped line/square/cube. Mainly useful for testing
+ functionality with curvilinear meshes.
+ """
- from meshmode.mesh import Mesh
- return Mesh(vertices, [grp],
- element_connectivity=None)
+ assert dim in [1, 2, 3]
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*dim, b=(0.5,)*dim,
+ n=(n,)*dim, order=order)
+
+ def m(x):
+ result = np.empty_like(x)
+ result[0] = (
+ 1.5*x[0] + np.cos(x[0])
+ + 0.1*np.sin(10*x[1]))
+ result[1] = (
+ 0.05*np.cos(10*x[0])
+ + 1.3*x[1] + np.sin(x[1]))
+ if len(x) == 3:
+ result[2] = x[2] + np.sin(x[0])
+ return result
+
+ from meshmode.mesh.processing import map_mesh
+ return map_mesh(mesh, m)
# }}}
diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py
index 292ee55077a178c70694bf8fd719b2fe42eef143..e27f8de352016e6f70c86be905fd4e247213d72b 100644
--- a/meshmode/mesh/io.py
+++ b/meshmode/mesh/io.py
@@ -1,8 +1,4 @@
-from __future__ import division
-from __future__ import absolute_import
-import six
-from six.moves import range
-from six.moves import zip
+from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
@@ -26,6 +22,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+import six
+from six.moves import range, zip
import numpy as np
from meshpy.gmsh_reader import ( # noqa
@@ -39,6 +37,8 @@ __doc__ = """
.. autofunction:: read_gmsh
.. autofunction:: generate_gmsh
+.. autofunction:: from_meshpy
+.. autofunction:: from_vertices_and_simplices
"""
@@ -171,11 +171,16 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
groups.append(group)
- return Mesh(vertices, groups, element_connectivity=None)
+ return Mesh(
+ vertices, groups,
+ nodal_adjacency=None,
+ facial_adjacency_groups=None)
# }}}
+# {{{ gmsh
+
def read_gmsh(filename, force_ambient_dim=None):
"""Read a gmsh mesh file from *filename* and return a
:class:`meshmode.mesh.Mesh`.
@@ -213,7 +218,7 @@ def generate_gmsh(source, dimensions, order=None, other_options=[],
mesh = recv.get_mesh()
if force_ambient_dim is None:
- AXIS_NAMES = "xyz"
+ AXIS_NAMES = "xyz" # noqa
dim = mesh.vertices.shape[0]
for idim in range(dim):
@@ -226,3 +231,66 @@ def generate_gmsh(source, dimensions, order=None, other_options=[],
break
return mesh
+
+# }}}
+
+
+# {{{ meshpy
+
+def from_meshpy(mesh_info, order=1):
+ """Imports a mesh from a :mod:`meshpy` *mesh_info* data structure,
+ which may be generated by either :mod:`meshpy.triangle` or
+ :mod:`meshpy_tet`.
+ """
+ from meshmode.mesh import Mesh
+ from meshmode.mesh.generation import make_group_from_vertices
+
+ vertices = np.array(mesh_info.points).T
+ elements = np.array(mesh_info.elements, np.int32)
+
+ grp = make_group_from_vertices(vertices, elements, order)
+
+ # FIXME: Should transfer boundary/volume markers
+
+ return Mesh(
+ vertices=vertices, groups=[grp],
+ nodal_adjacency=None,
+ facial_adjacency_groups=None)
+
+# }}}
+
+
+# {{{ from_vertices_and_simplices
+
+def from_vertices_and_simplices(vertices, simplices, order=1, fix_orientation=False):
+ """Imports a mesh from a numpy array of vertices and an array
+ of simplices.
+
+ :arg vertices:
+ An array of vertex coordinates with shape
+ *(ambient_dim, nvertices)*
+ :arg simplices:
+ An array *(nelements, nvertices)* of (mesh-wide)
+ vertex indices.
+ """
+ from meshmode.mesh import Mesh
+ from meshmode.mesh.generation import make_group_from_vertices
+
+ grp = make_group_from_vertices(vertices, simplices, order)
+
+ if fix_orientation:
+ from meshmode.mesh.processing import (
+ find_volume_mesh_element_group_orientation,
+ flip_simplex_element_group)
+ orient = find_volume_mesh_element_group_orientation(vertices, grp)
+ grp = flip_simplex_element_group(vertices, grp, orient < 0)
+
+ return Mesh(
+ vertices=vertices, groups=[grp],
+ nodal_adjacency=None,
+ facial_adjacency_groups=None)
+
+# }}}
+
+
+# vim: foldmethod=marker
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 94cd8b0951f5eb937c556bc239a2ebd7110bfb04..06631f9dedd6fdb1bcd3112a8c8157d4c72cb7f7 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1,7 +1,4 @@
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
-from functools import reduce
+from __future__ import division, absolute_import, print_function
__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
@@ -25,6 +22,9 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+from six.moves import range
+from functools import reduce
+
import numpy as np
import numpy.linalg as la
import modepy as mp
@@ -34,14 +34,15 @@ __doc__ = """
.. autofunction:: find_volume_mesh_element_orientations
.. autofunction:: perform_flips
.. autofunction:: find_bounding_box
-.. autofunction:: merge_dijsoint_meshes
+.. autofunction:: merge_disjoint_meshes
+.. autofunction:: map_mesh
.. autofunction:: affine_map
"""
# {{{ orientations
-def find_volume_mesh_element_group_orientation(mesh, grp):
+def find_volume_mesh_element_group_orientation(vertices, grp):
"""Return a positive floating point number for each positively
oriented element, and a negative floating point number for
each negatively oriented element.
@@ -56,11 +57,11 @@ def find_volume_mesh_element_group_orientation(mesh, grp):
"exclusively SimplexElementGroup-based meshes")
# (ambient_dim, nelements, nvertices)
- vertices = mesh.vertices[:, grp.vertex_indices]
+ my_vertices = vertices[:, grp.vertex_indices]
# (ambient_dim, nelements, nspan_vectors)
spanning_vectors = (
- vertices[:, :, 1:] - vertices[:, :, 0][:, :, np.newaxis])
+ my_vertices[:, :, 1:] - my_vertices[:, :, 0][:, :, np.newaxis])
ambient_dim = spanning_vectors.shape[0]
nspan_vectors = spanning_vectors.shape[-1]
@@ -81,6 +82,10 @@ def find_volume_mesh_element_group_orientation(mesh, grp):
from operator import xor
outer_prod = -reduce(xor, mvs)
+ if grp.dim == 1:
+ # FIXME: This is a little weird.
+ outer_prod = -outer_prod
+
return (outer_prod.I | outer_prod).as_scalar()
@@ -102,7 +107,8 @@ def find_volume_mesh_element_orientations(mesh, tolerate_unimplemented_checks=Fa
if tolerate_unimplemented_checks:
try:
signed_area_elements = \
- find_volume_mesh_element_group_orientation(mesh, grp)
+ find_volume_mesh_element_group_orientation(
+ mesh.vertices, grp)
except NotImplementedError:
result_grp_view[:] = float("nan")
else:
@@ -110,7 +116,8 @@ def find_volume_mesh_element_orientations(mesh, tolerate_unimplemented_checks=Fa
result_grp_view[:] = signed_area_elements
else:
signed_area_elements = \
- find_volume_mesh_element_group_orientation(mesh, grp)
+ find_volume_mesh_element_group_orientation(
+ mesh.vertices, grp)
assert not np.isnan(signed_area_elements).any()
result_grp_view[:] = signed_area_elements
@@ -158,7 +165,7 @@ def flip_simplex_element_group(vertices, grp, grp_flip_flags):
flipped_unit_nodes = barycentric_to_unit(flipped_bary_unit_nodes)
flip_matrix = mp.resampling_matrix(
- mp.simplex_onb(grp.dim, grp.order),
+ mp.simplex_best_available_basis(grp.dim, grp.order),
flipped_unit_nodes, grp.unit_nodes)
flip_matrix[np.abs(flip_matrix) < 1e-15] = 0
@@ -220,7 +227,7 @@ def find_bounding_box(mesh):
# {{{ merging
-def merge_dijsoint_meshes(meshes, skip_tests=False):
+def merge_disjoint_meshes(meshes, skip_tests=False):
if not meshes:
raise ValueError("must pass at least one mesh")
@@ -271,27 +278,48 @@ def merge_dijsoint_meshes(meshes, skip_tests=False):
# }}}
-# {{{ affine map
+# {{{ map
-def affine_map(mesh, A=None, b=None):
- """Apply the affine map *f(x)=Ax+b* to the geometry of *mesh*."""
+def map_mesh(mesh, f): # noqa
+ """Apply the map *f* to the mesh. *f* needs to accept and return arrays of
+ shape ``(ambient_dim, npoints)``."""
- vertices = np.einsum("ij,jv->iv", A, mesh.vertices) + b[:, np.newaxis]
+ vertices = f(mesh.vertices)
# {{{ assemble new groups list
new_groups = []
for group in mesh.groups:
- nodes = (
- np.einsum("ij,jen->ien", A, group.nodes)
- + b[:, np.newaxis, np.newaxis])
- new_groups.append(group.copy(nodes=nodes))
+ mapped_nodes = f(group.nodes.reshape(mesh.ambient_dim, -1))
+ new_groups.append(group.copy(
+ nodes=mapped_nodes.reshape(*group.nodes.shape)))
# }}}
from meshmode.mesh import Mesh
- return Mesh(vertices, new_groups, skip_tests=True)
+ return Mesh(vertices, new_groups, skip_tests=True,
+ nodal_adjacency=mesh.nodal_adjacency_init_arg(),
+ facial_adjacency_groups=mesh._facial_adjacency_groups)
+
+# }}}
+
+
+# {{{ affine map
+
+def affine_map(mesh, A=None, b=None): # noqa
+ """Apply the affine map *f(x)=Ax+b* to the geometry of *mesh*."""
+
+ if A is None:
+ A = np.eye(mesh.ambient_dim) # noqa
+
+ if b is None:
+ b = np.zeros(A.shape[0])
+
+ def f(x):
+ return np.einsum("ij,jv->iv", A, x) + b[:, np.newaxis]
+
+ return map_mesh(mesh, f)
# }}}
diff --git a/meshmode/mesh/refinement.py b/meshmode/mesh/refinement.py
index cd2ce945c7277274d0138ecff8bfedc8d12736a2..0a6518f0ec040c3b1a0adc8f18900a65ff197caf 100644
--- a/meshmode/mesh/refinement.py
+++ b/meshmode/mesh/refinement.py
@@ -1,4 +1,4 @@
-from __future__ import division
+from __future__ import division, print_function
__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
@@ -9,10 +9,8 @@ in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
-
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
-
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -22,153 +20,128 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-import numpy as np
-
-class VertexRay:
- def __init__(self, ray, pos):
- self.ray = ray
- self.pos = pos
-
-class _SplitFaceRecord(object):
- """
- .. attribute:: neighboring_elements
- .. attribute:: new_vertex_nr
-
- integer or None
- """
-class Adj:
- """One vertex-associated entry of a ray.
-
- .. attribute:: elements
-
- A list of numbers of elements adjacent to
- the edge following :attr:`vertex`
- along the ray.
-
- .. attribute:: velements
-
- A list of numbers of elements adjacent to
- :attr:`vertex`.
-
- """
- def __init__(self, vertex=None, elements=[None, None], velements
- =[None, None]):
- self.vertex = vertex
- self.elements = elements
- self.velements = velements
- def __str__(self):
- return 'vertex: ' + str(self.vertex) + ' ' + 'elements: ' +\
- str(self.elements) + ' velements: ' + str(self.velements)
-#map pair of vertices to ray and midpoint
-class PairMap:
- """Describes a segment of a ray between two vertices.
-
- .. attribute:: ray
-
- Index of the ray in the *rays* list.
-
- .. attribute:: d
-
- A :class:`bool` denoting direction in the ray,
- with *True* representing "positive" and *False*
- representing "negative".
+import numpy as np
+import itertools
+from six.moves import range
+
+
+class TreeRayNode(object):
+ """Describes a ray as a tree, this class represents each node in this tree
+ .. attribute:: left
+ Left child.
+ *None* if ray segment hasn't been split yet.
+ .. attribute:: right
+ Right child.
+ *None* if ray segment hasn't been split yet.
.. attribute:: midpoint
-
- Vertex index of the midpoint of this segment.
-
+ Vertex index of midpoint of this ray segment.
*None* if no midpoint has been assigned yet.
+ .. attribute:: adjacent_elements
+ List containing elements indices of elements adjacent
+ to this ray segment.
"""
+ def __init__(self, left_vertex, right_vertex, adjacent_elements=[]):
+ import copy
+ self.left = None
+ self.right = None
+ self.parent = None
+ self.left_vertex = left_vertex
+ self.right_vertex = right_vertex
+ self.midpoint = None
+ self.adjacent_elements = copy.deepcopy(adjacent_elements)
+ self.adjacent_add_diff = []
- def __init__(self, ray=None, d = True, midpoint=None):
- self.ray = ray
- #direction in ray, True means that second node (bigger index
- #) is after first in ray
- self.d = d
- self.midpoint = midpoint
- def __str__(self):
- return 'ray: ' + str(self.ray) + ' d: ' + str(self.d) + \
- ' midpoint: ' + str(self.midpoint)
class Refiner(object):
def __init__(self, mesh):
- from llist import dllist, dllistnode
- from meshmode.mesh.tesselate import tesselatetri
- self.tri_node_tuples, self.tri_result = tesselatetri()
- print self.tri_node_tuples
- print self.tri_result
+ from meshmode.mesh.tesselate import tesselatetet, tesselatetri
+ self.lazy = False
+ self.seen_tuple = {}
+ tri_node_tuples, tri_result = tesselatetri()
+ tet_node_tuples, tet_result = tesselatetet()
+ print(tri_node_tuples, tri_result)
+ #quadrilateral_node_tuples = [
+ #print tri_result, tet_result
+ self.simplex_node_tuples = [None, None, tri_node_tuples, tet_node_tuples]
+ self.simplex_result = [None, None, tri_result, tet_result]
+ #print tri_node_tuples, tri_result
+ #self.simplex_node_tuples, self.simplex_result = tesselatetet()
self.last_mesh = mesh
- # Indices in last_mesh that were split in the last round of
- # refinement. Only these elements may be split this time
- # around.
- self.last_split_elements = None
-
# {{{ initialization
- # a list of dllist instances containing Adj objects
- self.rays = []
-
- # mapping: (vertex_1, vertex_2) -> PairMap
+ # mapping: (vertex_1, vertex_2) -> TreeRayNode
# where vertex_i represents a vertex number
#
# Assumption: vertex_1 < vertex_2
self.pair_map = {}
nvertices = len(mesh.vertices[0])
-
- # list of dictionaries, with each entry corresponding to
- # one vertex.
- #
- # Each dictionary maps
- # ray number -> dllist node containing a :class:`Adj`,
- # (part of *rays*)
- self.vertex_to_ray = []
- import six
- for i in six.moves.range(nvertices):
- self.vertex_to_ray.append({})
+ #array containing element whose edge lies on corresponding vertex
+ self.hanging_vertex_element = []
+ for i in range(nvertices):
+ self.hanging_vertex_element.append([])
+
for grp in mesh.groups:
iel_base = grp.element_nr_base
- for iel_grp in six.moves.range(grp.nelements):
- #use six.moves.range
-
+ for iel_grp in range(grp.nelements):
vert_indices = grp.vertex_indices[iel_grp]
+ for i in range(len(vert_indices)):
+ for j in range(i+1, len(vert_indices)):
- for i in six.moves.range(len(vert_indices)):
- for j in six.moves.range(i+1, len(vert_indices)):
-
- #minimum and maximum of the two indices for storing
+ #minimum and maximum of the two indices for storing
#data in vertex_pair
- mn_idx = min(vert_indices[i], vert_indices[j])
- mx_idx = max(vert_indices[i], vert_indices[j])
+ min_index = min(vert_indices[i], vert_indices[j])
+ max_index = max(vert_indices[i], vert_indices[j])
- vertex_pair = (mn_idx, mx_idx)
+ vertex_pair = (min_index, max_index)
+ #print(vertex_pair)
if vertex_pair not in self.pair_map:
- els = []
- els.append(iel_base+iel_grp)
- fr = Adj(mn_idx, els)
- to = Adj(mx_idx, [None, None])
- self.rays.append(dllist([fr, to]))
- self.pair_map[vertex_pair] = PairMap(len(self.rays) - 1)
- self.vertex_to_ray[mn_idx][len(self.rays)-1]\
- = self.rays[len(self.rays)-1].nodeat(0)
- self.vertex_to_ray[mx_idx][len(self.rays)-1]\
- = self.rays[len(self.rays)-1].nodeat(1)
- else:
- self.rays[self.pair_map[vertex_pair].ray].nodeat(0).value.elements.append(iel_base+iel_grp)
+ self.pair_map[vertex_pair] = TreeRayNode(min_index, max_index)
+ self.pair_map[vertex_pair].adjacent_elements.append(iel_base+iel_grp)
+ elif (iel_base+iel_grp) not in self.pair_map[vertex_pair].adjacent_elements:
+ (self.pair_map[vertex_pair].
+ adjacent_elements.append(iel_base+iel_grp))
# }}}
- #print mesh.groups[0].vertex_indices
+ #print(vert_indices)
+ #generate reference tuples
+ self.index_to_node_tuple = []
+ self.index_to_midpoint_tuple = []
+ for d in range(len(vert_indices)):
+ dim = d + 1
+ cur_index_to_node_tuple = [()] * dim
+ for i in range(0, dim-1):
+ cur_index_to_node_tuple[0] = cur_index_to_node_tuple[0] + (0,)
+ for i in range(1, dim):
+ for j in range(1, dim):
+ if i == j:
+ cur_index_to_node_tuple[i] = cur_index_to_node_tuple[i] + (2,)
+ else:
+ cur_index_to_node_tuple[i] = cur_index_to_node_tuple[i] + (0,)
+ cur_index_to_midpoint_tuple = [()] * (int((dim * (dim - 1)) / 2))
+ curind = 0
+ for ind1 in range(0, len(cur_index_to_node_tuple)):
+ for ind2 in range(ind1+1, len(cur_index_to_node_tuple)):
+ i = cur_index_to_node_tuple[ind1]
+ j = cur_index_to_node_tuple[ind2]
+ #print(i, j)
+ for k in range(0, dim-1):
+ cur = int((i[k] + j[k]) / 2)
+ cur_index_to_midpoint_tuple[curind] = cur_index_to_midpoint_tuple[curind] + (cur,)
+ curind += 1
+ self.index_to_node_tuple.append(cur_index_to_node_tuple)
+ self.index_to_midpoint_tuple.append(cur_index_to_midpoint_tuple)
'''
- self.ray_vertices = np.empty([len(mesh.groups[0].vertex_indices) *
- len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1) / 2, 2],
+ self.ray_vertices = np.empty([len(mesh.groups[0].vertex_indices) *
+ len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1) / 2, 2],
dtype=np.int32)
- self.ray_elements = np.zeros([len(mesh.groups[0].vertex_indices) *
+ self.ray_elements = np.zeros([len(mesh.groups[0].vertex_indices) *
len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1)
/ 2, 1, 2], dtype=np.int32)
self.vertex_to_ray = []
-
for i in mesh.vertices[0]:
self.vertex_to_ray.append([]);
count = 0
@@ -184,7 +157,7 @@ class Refiner(object):
self.vertex_to_ray[grp.vertex_indices[i][k]].append(temp2)
count += 1
ind = 0
- #print self.ray_vertices
+ #print(self.ray_vertices)
for i in self.ray_vertices:
which = 0
for grp in mesh.groups:
@@ -214,16 +187,16 @@ class Refiner(object):
np.bool)
def get_current_mesh(self):
-
+
from meshmode.mesh import Mesh
#return Mesh(vertices, [grp], element_connectivity=self.generate_connectivity(len(self.last_mesh.groups[group].vertex_indices) \
# + count*3))
- groups = []
+ groups = []
grpn = 0
for grp in self.last_mesh.groups:
groups.append(np.empty([len(grp.vertex_indices),
len(self.last_mesh.groups[grpn].vertex_indices[0])], dtype=np.int32))
- for iel_grp in xrange(grp.nelements):
+ for iel_grp in range(grp.nelements):
for i in range(0, len(grp.vertex_indices[iel_grp])):
groups[grpn][iel_grp][i] = grp.vertex_indices[iel_grp][i]
grpn += 1
@@ -235,459 +208,552 @@ class Refiner(object):
self.last_mesh = Mesh(self.last_mesh.vertices, grp,\
element_connectivity=self.generate_connectivity(len(self.last_mesh.groups[0].vertex_indices),\
len(self.last_mesh.vertices[0])))
-
+
return self.last_mesh
+
+ def get_leaves(self, cur_node):
+ queue = [cur_node]
+ res = []
+ while queue:
+ vertex = queue.pop(0)
+ #if leaf node
+ if vertex.left is None and vertex.right is None:
+ res.append(vertex)
+ else:
+ queue.append(vertex.left)
+ queue.append(vertex.right)
+ return res
+
+ def get_subtree(self, cur_node):
+ queue = [cur_node]
+ res = []
+ while queue:
+ vertex = queue.pop(0)
+ res.append(vertex)
+ if not (vertex.left is None and vertex.right is None):
+ queue.append(vertex.left)
+ queue.append(vertex.right)
+ return res
+
+ def apply_diff(self, cur_node, new_hanging_vertex_elements):
+ for el in cur_node.adjacent_add_diff:
+ if el not in cur_node.adjacent_elements:
+ cur_node.adjacent_elements.append(el)
+ if el not in new_hanging_vertex_elements[cur_node.left_vertex]:
+ new_hanging_vertex_elements[cur_node.left_vertex].append(el)
+ if el not in new_hanging_vertex_elements[cur_node.right_vertex]:
+ new_hanging_vertex_elements[cur_node.right_vertex].append(el)
+ if cur_node.left is not None and cur_node.right is not None:
+ if el not in cur_node.left.adjacent_add_diff:
+ cur_node.left.adjacent_add_diff.append(el)
+ if el not in cur_node.right.adjacent_add_diff:
+ cur_node.right.adjacent_add_diff.append(el)
+ cur_node.adjacent_add_diff = []
+
+# def propagate(self, cur_node, new_hanging_vertex_elements):
+# if cur_node.parent is not None:
+# parent_node = cur_node.parent
+# self.propagate(parent_node, new_hanging_vertex_elements)
+# self.apply_diff(parent_node, new_hanging_vertex_elements)
+
+ def get_root(self, cur_node):
+ while(cur_node.parent is not None):
+ cur_node = cur_node.parent
+ return cur_node
+
+ def propagate_tree(self, cur_node, new_hanging_vertex_elements, element_to_element):
+ vertex_tuple = (min(cur_node.left_vertex, cur_node.right_vertex), max(cur_node.left_vertex, cur_node.right_vertex))
+ self.seen_tuple[vertex_tuple] = True
+ self.apply_diff(cur_node, new_hanging_vertex_elements)
+ if cur_node.left is not None and cur_node.right is not None:
+ self.propagate_tree(cur_node.left, new_hanging_vertex_elements, element_to_element)
+ self.propagate_tree(cur_node.right, new_hanging_vertex_elements, element_to_element)
+ else:
+ for el in cur_node.adjacent_elements:
+ element_to_element[el].update(cur_node.adjacent_elements)
+ for el in new_hanging_vertex_elements[cur_node.left_vertex]:
+ element_to_element[el].update(new_hanging_vertex_elements[cur_node.left_vertex])
+ for el in new_hanging_vertex_elements[cur_node.right_vertex]:
+ element_to_element[el].update(new_hanging_vertex_elements[cur_node.right_vertex])
+
+# def remove_from_subtree(self, cur_node, new_hanging_vertex_elements, to_remove):
+# if self.lazy:
+# self.propagate(cur_node, new_hanging_vertex_elements)
+# if to_remove in cur_node.adjacent_add_diff:
+# cur_node.adjacent_add_diff.remove(to_remove)
+# if to_remove not in cur_node.adjacent_remove_diff:
+# cur_node.adjacent_remove_diff.append(to_remove)
+# else:
+# subtree = self.get_subtree(cur_node)
+# for node in subtree:
+# if to_remove in node.adjacent_elements:
+# node.adjacent_elements.remove(to_remove)
+# if to_remove in new_hanging_vertex_elements[node.left_vertex]:
+# new_hanging_vertex_elements[node.left_vertex].remove(to_remove)
+# if to_remove in new_hanging_vertex_elements[node.right_vertex]:
+# new_hanging_vertex_elements[node.right_vertex].remove(to_remove)
+
+ def add_to_subtree(self, cur_node, new_hanging_vertex_elements, to_add):
+ if self.lazy:
+ if to_add not in cur_node.adjacent_add_diff:
+ cur_node.adjacent_add_diff.append(to_add)
+ else:
+ subtree = self.get_subtree(cur_node)
+ for node in subtree:
+ if to_add not in node.adjacent_elements:
+ node.adjacent_elements.append(to_add)
+ if to_add not in new_hanging_vertex_elements[node.left_vertex]:
+ new_hanging_vertex_elements[node.left_vertex].append(to_add)
+ if to_add not in new_hanging_vertex_elements[node.right_vertex]:
+ new_hanging_vertex_elements[node.right_vertex].append(to_add)
+
#refine_flag tells you which elements to split as a numpy array of bools
def refine(self, refine_flags):
- """
- :return: a refined mesh
- """
-
- # multi-level mapping:
- # {
- # dimension of intersecting entity (0=vertex,1=edge,2=face,...)
- # :
- # { frozenset of end vertices : _SplitFaceRecord }
- # }
- '''
- import numpy as np
- count = 0
- for i in refine_flags:
- if i:
- count += 1
-
- le = len(self.last_mesh.vertices[0])
-
- vertices = np.empty([len(self.last_mesh.vertices), len(self.last_mesh.groups[0].vertex_indices[0])
- * count + le])
- vertex_indices = np.empty([len(self.last_mesh.groups[0].vertex_indices)
- + count*3,
- len(self.last_mesh.groups[0].vertex_indices[0])], dtype=np.int32)
- indices_it = len(self.last_mesh.groups[0].vertex_indices)
- for i in range(0, len(self.last_mesh.vertices)):
- for j in range(0, len(self.last_mesh.vertices[i])):
- vertices[i][j] = self.last_mesh.vertices[i][j]
- for i in range(0, len(self.last_mesh.groups[0].vertex_indices)):
- for j in range(0, len(self.last_mesh.groups[0].vertex_indices[i])):
- vertex_indices[i][j] = self.last_mesh.groups[0].vertex_indices[i][j]
-
-
- import itertools
- for i in range(0, len(refine_flags)):
- if refine_flags[i]:
- for subset in itertools.combinations(self.last_mesh.groups[0].vertex_indices[i],
- len(self.last_mesh.groups[0].vertex_indices[i]) - 1):
- for j in range(0, len(self.last_mesh.vertices)):
- avg = 0
- for k in subset:
- avg += self.last_mesh.vertices[j][k]
- avg /= len(self.last_mesh.vertices)
- self.last_mesh.vertices[j][le] = avg
- le++
- vertex_indices[indices_it][0] = self.last_mesh.groups[0].vertex_indices[i][0]
- '''
- '''
import numpy as np
- count = 0
- for i in refine_flags:
- if i:
- count += 1
- #print count
- #print self.last_mesh.vertices
- #print vertices
- if(len(self.last_mesh.groups[0].vertex_indices[0]) == 3):
- for group in range(0, len(self.last_mesh.groups)):
- le = len(self.last_mesh.vertices[0])
- vertices = np.empty([len(self.last_mesh.vertices),
- len(self.last_mesh.groups[group].vertex_indices[0])
- * count + le])
- vertex_indices = np.empty([len(self.last_mesh.groups[group].vertex_indices)
- + count*3,
- len(self.last_mesh.groups[group].vertex_indices[0])], dtype=np.int32)
- indices_i = 0
- for i in range(0, len(self.last_mesh.vertices)):
- for j in range(0, len(self.last_mesh.vertices[i])):
- vertices[i][j] = self.last_mesh.vertices[i][j]
- #for i in range(0, len(self.last_mesh.groups[group].vertex_indices)):
- #for j in range(0, len(self.last_mesh.groups[group].vertex_indices[i])):
- #vertex_indices[i][j] = self.last_mesh.groups[group].vertex_indices[i][j]
- for fl in range(0, len(refine_flags)):
- if(refine_flags[fl]):
- i = self.last_mesh.groups[group].vertex_indices[fl]
- for j in range(0, len(i)):
- for k in range(j + 1, len(i)):
- for l in range(0, 3):
- #print self.last_mesh.vertices[l][i[j]], ' ', self.last_mesh.vertices[l][i[k]], '=', ((self.last_mesh.vertices[l][i[j]] + self.last_mesh.vertices[l][i[k]]) / 2)
- vertices[l][le]=((self.last_mesh.vertices[l][i[j]] + self.last_mesh.vertices[l][i[k]]) / 2)
- le += 1
- vertex_indices[indices_i][0] = i[0]
- vertex_indices[indices_i][1] = le-3
- vertex_indices[indices_i][2] = le-2
- indices_i += 1
- vertex_indices[indices_i][0] = i[1]
- vertex_indices[indices_i][1] = le-1
- vertex_indices[indices_i][2] = le-3
- indices_i += 1
- vertex_indices[indices_i][0] = i[2]
- vertex_indices[indices_i][1] = le-2
- vertex_indices[indices_i][2] = le-1
- indices_i += 1
- vertex_indices[indices_i][0] = le-3
- vertex_indices[indices_i][1] = le-2
- vertex_indices[indices_i][2] = le-1
- indices_i += 1
- else:
- for j in range(0, len(self.last_mesh.groups[group].vertex_indices[fl])):
- vertex_indices[indices_i][j] = self.last_mesh.groups[group].vertex_indices[fl][j]
- indices_i += 1
- '''
- import numpy as np
- import six
- # if (isinstance(self.last_mesh.groups[0], SimplexElementGroup) and
- # self.last_mesh.groups[0].dim == 2):
- print 'refining'
- if(len(self.last_mesh.groups[0].vertex_indices[0]) == 3):
- groups = []
- midpoint_already = {}
+ #vertices and groups for next generation
+ nvertices = len(self.last_mesh.vertices[0])
+
+ groups = []
+
+ midpoint_already = set()
+ grpn = 0
+ totalnelements = 0
+
+ for grp in self.last_mesh.groups:
+ iel_base = grp.element_nr_base
nelements = 0
- nvertices = len(self.last_mesh.vertices[0])
- grpn = 0
- totalnelements=0
-
- # {{{ create new vertices array and each group's vertex_indices
-
- for grp in self.last_mesh.groups:
- iel_base = grp.element_nr_base
- nelements = 0
- #print grp.nelements
- for iel_grp in six.moves.range(grp.nelements):
- nelements += 1
- vert_indices = grp.vertex_indices[iel_grp]
- if refine_flags[iel_base+iel_grp]:
- nelements += 3
- for i in six.moves.range(0, len(vert_indices)):
- for j in six.moves.range(i+1, len(vert_indices)):
- mn_idx = min(vert_indices[i], vert_indices[j])
- mx_idx = max(vert_indices[i], vert_indices[j])
- vertex_pair = (mn_idx, mx_idx)
- if vertex_pair not in midpoint_already and \
- self.pair_map[vertex_pair].midpoint is None:
+ for iel_grp in range(grp.nelements):
+ nelements += 1
+ vertex_indices = grp.vertex_indices[iel_grp]
+ if refine_flags[iel_base+iel_grp]:
+ cur_dim = len(grp.vertex_indices[iel_grp])-1
+ nelements += len(self.simplex_result[cur_dim]) - 1
+ for i in range(len(vertex_indices)):
+ for j in range(i+1, len(vertex_indices)):
+ i_index = vertex_indices[i]
+ j_index = vertex_indices[j]
+ index_tuple = (i_index, j_index) if i_index < j_index else (j_index, i_index)
+ if index_tuple not in midpoint_already and \
+ self.pair_map[index_tuple].midpoint is None:
nvertices += 1
- midpoint_already[vertex_pair] = True
- groups.append(np.empty([nelements,
- len(grp.vertex_indices[0])], dtype=np.int32))
- grpn += 1
- totalnelements += nelements
-
- vertices = np.empty([len(self.last_mesh.vertices), nvertices])
-
- # }}}
-
- # {{{ bring over vertex_indices and vertices info from previous generation
-
- for i in six.moves.range(0, len(self.last_mesh.vertices)):
- for j in six.moves.range(0, len(self.last_mesh.vertices[i])):
- # always use v[i,j]
- vertices[i,j] = self.last_mesh.vertices[i,j]
- grpn = 0
- for grp in self.last_mesh.groups:
- for iel_grp in six.moves.range(grp.nelements):
- for i in six.moves.range(0, len(grp.vertex_indices[iel_grp])):
- groups[grpn][iel_grp][i] = grp.vertex_indices[iel_grp][i]
- grpn += 1
- grpn = 0
-
- # }}}
-
- vertices_idx = len(self.last_mesh.vertices[0])
- for grp in self.last_mesh.groups:
- iel_base = grp.element_nr_base
- nelements_in_grp = len(grp.vertex_indices)
-
- # np.where
- for iel_grp in six.moves.range(grp.nelements):
- if refine_flags[iel_base+iel_grp]:
-
- # {{{ split element
-
- # {{{ go through vertex pairs in element
-
- # stores indices of all midpoints for this element
- # (in order of vertex pairs in elements)
- verts = []
-
- vert_indices = grp.vertex_indices[iel_grp]
- for i in six.moves.range(0, len(vert_indices)):
- for j in six.moves.range(i+1, len(vert_indices)):
- mn_idx = min(vert_indices[i], vert_indices[j])
- mx_idx = max(vert_indices[i], vert_indices[j])
- vertex_pair = (mn_idx, mx_idx)
-
- # {{{ create midpoint if it doesn't exist already
-
- cur_pmap = self.pair_map[vertex_pair]
- if cur_pmap.midpoint is None:
- self.pair_map[vertex_pair].midpoint = vertices_idx
- vertex_pair1 = (mn_idx, vertices_idx)
- vertex_pair2 = (mx_idx, vertices_idx)
- self.pair_map[vertex_pair1] =\
- PairMap(cur_pmap.ray, cur_pmap.d, None)
- self.pair_map[vertex_pair2] =\
- PairMap(cur_pmap.ray, not cur_pmap.d, None)
- self.vertex_to_ray.append({})
-
- # FIXME: Check where the new Adj.elements
- # gets populated.
-
- # try and collapse the two branches by setting up variables
- # ahead of time
+ midpoint_already.add(index_tuple)
+ groups.append(np.empty([nelements, len(grp.vertex_indices[0])], dtype=np.int32))
+ grpn += 1
+ totalnelements += nelements
+
+ vertices = np.empty([len(self.last_mesh.vertices), nvertices])
+
+ new_hanging_vertex_element = [
+ [] for i in range(nvertices)]
+
+# def remove_element_from_connectivity(vertices, new_hanging_vertex_elements, to_remove):
+# #print(vertices)
+# import itertools
+# if len(vertices) == 2:
+# min_vertex = min(vertices[0], vertices[1])
+# max_vertex = max(vertices[0], vertices[1])
+# ray = self.pair_map[(min_vertex, max_vertex)]
+# self.remove_from_subtree(ray, new_hanging_vertex_elements, to_remove)
+# return
+#
+# cur_dim = len(vertices)-1
+# element_rays = []
+# midpoints = []
+# split_possible = True
+# for i in range(len(vertices)):
+# for j in range(i+1, len(vertices)):
+# min_vertex = min(vertices[i], vertices[j])
+# max_vertex = max(vertices[i], vertices[j])
+# element_rays.append(self.pair_map[(min_vertex, max_vertex)])
+# if element_rays[len(element_rays)-1].midpoint is not None:
+# midpoints.append(element_rays[len(element_rays)-1].midpoint)
+# else:
+# split_possible = False
+
+ #for node in element_rays:
+ #self.remove_from_subtree(node, new_hanging_vertex_elements, to_remove)
+ #if split_possible:
+# if split_possible:
+# node_tuple_to_coord = {}
+# for node_index, node_tuple in enumerate(self.index_to_node_tuple[cur_dim]):
+# node_tuple_to_coord[node_tuple] = vertices[node_index]
+# for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]):
+# node_tuple_to_coord[midpoint_tuple] = midpoints[midpoint_index]
+# for i in range(len(self.simplex_result[cur_dim])):
+# next_vertices = []
+# for j in range(len(self.simplex_result[cur_dim][i])):
+# next_vertices.append(node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]])
+# all_rays_present = True
+# for v1 in range(len(next_vertices)):
+# for v2 in range(v1+1, len(next_vertices)):
+# vertex_tuple = (min(next_vertices[v1], next_vertices[v2]), max(next_vertices[v1], next_vertices[v2]))
+# if vertex_tuple not in self.pair_map:
+# all_rays_present = False
+# if all_rays_present:
+# remove_element_from_connectivity(next_vertices, new_hanging_vertex_elements, to_remove)
+# else:
+# split_possible = False
+# if not split_possible:
+# next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1))
+# for next_vertices in next_vertices_list:
+# remove_element_from_connectivity(next_vertices, new_hanging_vertex_elements, to_remove)
+
+ def add_element_to_connectivity(vertices, new_hanging_vertex_elements, to_add):
+ if len(vertices) == 2:
+ min_vertex = min(vertices[0], vertices[1])
+ max_vertex = max(vertices[0], vertices[1])
+ ray = self.pair_map[(min_vertex, max_vertex)]
+ self.add_to_subtree(ray, new_hanging_vertex_elements, to_add)
+ return
+
+ cur_dim = len(vertices)-1
+ element_rays = []
+ midpoints = []
+ split_possible = True
+ for i in range(len(vertices)):
+ for j in range(i+1, len(vertices)):
+ min_vertex = min(vertices[i], vertices[j])
+ max_vertex = max(vertices[i], vertices[j])
+ element_rays.append(self.pair_map[(min_vertex, max_vertex)])
+ if element_rays[len(element_rays)-1].midpoint is not None:
+ midpoints.append(element_rays[len(element_rays)-1].midpoint)
+ else:
+ split_possible = False
+ #for node in element_rays:
+ #self.add_to_subtree(node, new_hanging_vertex_elements, to_add)
+ if split_possible:
+ node_tuple_to_coord = {}
+ for node_index, node_tuple in enumerate(self.index_to_node_tuple[cur_dim]):
+ node_tuple_to_coord[node_tuple] = vertices[node_index]
+ for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]):
+ node_tuple_to_coord[midpoint_tuple] = midpoints[midpoint_index]
+ for i in range(len(self.simplex_result[cur_dim])):
+ next_vertices = []
+ for j in range(len(self.simplex_result[cur_dim][i])):
+ next_vertices.append(node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]])
+ all_rays_present = True
+ for v1 in range(len(next_vertices)):
+ for v2 in range(v1+1, len(next_vertices)):
+ vertex_tuple = (min(next_vertices[v1], next_vertices[v2]), max(next_vertices[v1], next_vertices[v2]))
+ if vertex_tuple not in self.pair_map:
+ all_rays_present = False
+ if all_rays_present:
+ add_element_to_connectivity(next_vertices, new_hanging_vertex_elements, to_add)
+ else:
+ split_possible = False
+ if not split_possible:
+ next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1))
+ for next_vertices in next_vertices_list:
+ add_element_to_connectivity(next_vertices, new_hanging_vertex_elements, to_add)
+# for node in element_rays:
+# self.add_element_to_connectivity(node, new_hanging_vertex_elements, to_add)
+ # leaves = self.get_subtree(node)
+ # for leaf in leaves:
+ # if to_add not in leaf.adjacent_elements:
+ # leaf.adjacent_elements.append(to_add)
+ # if to_add not in new_hanging_vertex_elements[leaf.left_vertex]:
+ # new_hanging_vertex_elements[leaf.left_vertex].append(to_add)
+ # if to_add not in new_hanging_vertex_elements[leaf.right_vertex]:
+ # new_hanging_vertex_elements[leaf.right_vertex].append(to_add)
+
+# next_element_rays = []
+# for i in range(len(element_rays)):
+# for j in range(i+1, len(element_rays)):
+# if element_rays[i].midpoint is not None and element_rays[j].midpoint is not None:
+# min_midpoint = min(element_rays[i].midpoint, element_rays[j].midpoint)
+# max_midpoint = max(element_rays[i].midpoint, element_rays[j].midpoint)
+# vertex_pair = (min_midpoint, max_midpoint)
+# if vertex_pair in self.pair_map:
+# next_element_rays.append(self.pair_map[vertex_pair])
+# cur_next_rays = []
+# if element_rays[i].left_vertex == element_rays[j].left_vertex:
+# cur_next_rays = [element_rays[i].left, element_rays[j].left, self.pair_map[vertex_pair]]
+# if element_rays[i].right_vertex == element_rays[j].right_vertex:
+# cur_next_rays = [element_rays[i].right, element_rays[j].right, self.pair_map[vertex_pair]]
+# if element_rays[i].left_vertex == element_rays[j].right_vertex:
+# cur_next_rays = [element_rays[i].left, element_rays[j].right, self.pair_map[vertex_pair]]
+# if element_rays[i].right_vertex == element_rays[j].left_vertex:
+# cur_next_rays = [element_rays[i].right, element_rays[j].left, self.pair_map[vertex_pair]]
+# assert (cur_next_rays != [])
+# #print cur_next_rays
+# add_element_to_connectivity(cur_next_rays, new_hanging_vertex_elements, to_add)
+# else:
+# return
+# else:
+# return
+# add_element_to_connectivity(next_element_rays, new_hanging_vertex_elements, to_add)
+
+ def add_hanging_vertex_el(v_index, el):
+ assert not (v_index == 37 and el == 48)
+
+ new_hanging_vertex_element[v_index].append(el)
+
+# def remove_ray_el(ray, el):
+# ray.remove(el)
+
+ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_grp):
+ for grp in groups:
+ iel_base = 0
+ for iel_grp in range(nelements_in_grp):
+ vertex_indices = grp[iel_grp]
+ for i in range(len(vertex_indices)):
+ for j in range(i+1, len(vertex_indices)):
+ min_index = min(vertex_indices[i], vertex_indices[j])
+ max_index = max(vertex_indices[i], vertex_indices[j])
+ cur_node = self.pair_map[(min_index, max_index)]
+ #print iel_base+iel_grp, cur_node.left_vertex, cur_node.right_vertex
+ #if (iel_base + iel_grp) not in cur_node.adjacent_elements:
+ #print min_index, max_index
+ #print iel_base + iel_grp, cur_node.left_vertex, cur_node.right_vertex, cur_node.adjacent_elements
+ #assert (4 in new_hanging_vertex_elements[cur_node.left_vertex] or 4 in new_hanging_vertex_elements[cur_node.right_vertex])
+ assert ((iel_base+iel_grp) in cur_node.adjacent_elements)
+ assert((iel_base+iel_grp) in new_hanging_vertex_elements[cur_node.left_vertex])
+ assert((iel_base+iel_grp) in new_hanging_vertex_elements[cur_node.right_vertex])
+
+ for i in range(len(self.last_mesh.vertices)):
+ for j in range(len(self.last_mesh.vertices[i])):
+ vertices[i,j] = self.last_mesh.vertices[i,j]
+ import copy
+ if i == 0:
+ new_hanging_vertex_element[j] = copy.deepcopy(self.hanging_vertex_element[j])
+ grpn = 0
+ for grp in self.last_mesh.groups:
+ for iel_grp in range(grp.nelements):
+ for i in range(len(grp.vertex_indices[iel_grp])):
+ groups[grpn][iel_grp][i] = grp.vertex_indices[iel_grp][i]
+ grpn += 1
+
+ grpn = 0
+ vertices_index = len(self.last_mesh.vertices[0])
+ nelements_in_grp = grp.nelements
+ for grp in self.last_mesh.groups:
+ iel_base = grp.element_nr_base
+ for iel_grp in range(grp.nelements):
+ if refine_flags[iel_base+iel_grp]:
+ midpoint_vertices = []
+ midpoint_tuples = []
+ vertex_indices = grp.vertex_indices[iel_grp]
+ #if simplex
+ if len(grp.vertex_indices[iel_grp]) == len(self.last_mesh.vertices)+1:
+ for i in range(len(vertex_indices)):
+ for j in range(i+1, len(vertex_indices)):
+ min_index = min(vertex_indices[i], vertex_indices[j])
+ max_index = max(vertex_indices[i], vertex_indices[j])
+ cur_node = self.pair_map[(min_index, max_index)]
+ if cur_node.midpoint is None:
+ cur_node.midpoint = vertices_index
import copy
- if self.pair_map[vertex_pair].d:
- velements = self.vertex_to_ray[mn_idx][cur_pmap.ray].value.elements
- self.rays[cur_pmap.ray].insert(Adj(vertices_idx, copy.deepcopy(velements),
- copy.deepcopy(velements)), \
- self.vertex_to_ray[mx_idx][cur_pmap.ray])
- self.vertex_to_ray[vertices_idx][cur_pmap.ray] = \
- self.vertex_to_ray[mx_idx][cur_pmap.ray].prev
- else:
- velements = self.vertex_to_ray[mx_idx][cur_pmap.ray].value.elements
- self.rays[cur_pmap.ray].insert(Adj(vertices_idx, copy.deepcopy(velements),
- copy.deepcopy(velements)), \
- self.vertex_to_ray[mn_idx][cur_pmap.ray])
- self.vertex_to_ray[vertices_idx][cur_pmap.ray] = \
- self.vertex_to_ray[mn_idx][cur_pmap.ray].prev
-
- # compute location of midpoint
- for k in six.moves.range(len(self.last_mesh.vertices)):
- vertices[k,vertices_idx] =\
- (self.last_mesh.vertices[k,vert_indices[i]]\
- +self.last_mesh.vertices[k,vert_indices[j]]) / 2.0
-
- verts.append(vertices_idx)
- vertices_idx += 1
+ cur_node.left = TreeRayNode(min_index, vertices_index,
+ copy.deepcopy(cur_node.adjacent_elements))
+ cur_node.left.parent = cur_node
+ cur_node.right = TreeRayNode(max_index, vertices_index,
+ copy.deepcopy(cur_node.adjacent_elements))
+ cur_node.right.parent = cur_node
+ vertex_pair1 = (min_index, vertices_index)
+ vertex_pair2 = (max_index, vertices_index)
+ self.pair_map[vertex_pair1] = cur_node.left
+ self.pair_map[vertex_pair2] = cur_node.right
+ for k in range(len(self.last_mesh.vertices)):
+ vertices[k, vertices_index] = \
+ (self.last_mesh.vertices[k, vertex_indices[i]] +
+ self.last_mesh.vertices[k, vertex_indices[j]]) / 2.0
+ midpoint_vertices.append(vertices_index)
+ vertices_index += 1
else:
- verts.append(self.pair_map[vertex_pair].midpoint)
-
- # }}}
-
- # }}}
-
- # {{{ fix connectivity
-
- # new elements will be nels+0 .. nels+2
- # (While those don't exist yet, we generate connectivity for them
- # anyhow.)
-
- nnew_elements = 0
-
- unique_vertex_pairs = [
- (i,j) for i in range(3) for j in range(i+1, 3)]
- for i, j in unique_vertex_pairs:
- mn_idx = min(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j])
- mx_idx = max(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j])
- if mn_idx == grp.vertex_indices[iel_grp][i]:
- min_element_index = i
- max_element_index = j
- else:
- min_element_index = j
- max_element_index = i
- vertex_pair = (mn_idx, mx_idx)
- cur_pmap = self.pair_map[vertex_pair]
- if cur_pmap.d:
- start_node =\
- self.vertex_to_ray[mn_idx][cur_pmap.ray]
- end_node = self.vertex_to_ray[mx_idx][cur_pmap.ray]
- first_element_index = min_element_index
- second_element_index = max_element_index
- else:
- start_node =\
- self.vertex_to_ray[mx_idx][cur_pmap.ray]
- end_node = self.vertex_to_ray[mn_idx][cur_pmap.ray]
- first_element_index = max_element_index
- second_element_index = min_element_index
- midpoint_node=\
- self.vertex_to_ray[cur_pmap.midpoint][cur_pmap.ray]
- # hop along ray from start node to midpoint node
- while start_node != midpoint_node:
- # replace old (big) element with new
- # (refined) element
- node_els = start_node.value.elements
- #print node_els
- iold_el = node_els.index(iel_base+iel_grp)
- node_els[iold_el] = iel_base+nelements_in_grp+first_element_index
- start_node = start_node.next
- # hop along ray from midpoint node to end node
- while start_node != end_node:
- #replace old (big) element with new
- # (refined element
- node_els = start_node.value.elements
- iold_el = node_els.index(iel_base+iel_grp)
- node_els[iold_el] = iel_base+nelements_in_grp+second_element_index
- start_node = start_node.next
-
- # }}}
+ cur_midpoint = cur_node.midpoint
+ midpoint_vertices.append(cur_midpoint)
+
#generate new rays
- from llist import dllist, dllistnode
- ind = 0
-
- for i in six.moves.range(0, len(verts)):
- for j in six.moves.range(i+1, len(verts)):
- mn_vert = min(verts[i], verts[j])
- mx_vert = max(verts[i], verts[j])
- vertex_pair = (mn_vert, mx_vert)
- els = []
- els.append(iel_base+iel_grp)
- els.append(iel_base+nelements_in_grp+ind)
-
- fr = Adj(mn_vert, els)
-
- # We're creating a new ray, and this is the end node
- # of it.
- to = Adj(mx_vert, [None, None])
-
- self.rays.append(dllist([fr, to]))
- self.pair_map[vertex_pair] = PairMap(len(self.rays) - 1)
- self.vertex_to_ray[mn_vert][len(self.rays)-1]\
- = self.rays[len(self.rays)-1].nodeat(0)
- self.vertex_to_ray[mx_vert][len(self.rays)-1]\
- = self.rays[len(self.rays)-1].nodeat(1)
- ind += 1
- ind = 0
-
- #map vertex indices to coordinates
+ cur_dim = len(grp.vertex_indices[0])-1
+ for i in range(len(midpoint_vertices)):
+ for j in range(i+1, len(midpoint_vertices)):
+ min_index = min(midpoint_vertices[i], midpoint_vertices[j])
+ max_index = max(midpoint_vertices[i], midpoint_vertices[j])
+ vertex_pair = (min_index, max_index)
+ if vertex_pair in self.pair_map:
+ continue
+ self.pair_map[vertex_pair] = TreeRayNode(min_index, max_index, [])
node_tuple_to_coord = {}
- node_tuple_to_coord[(0, 0)] = grp.vertex_indices[iel_grp][0]
- node_tuple_to_coord[(2, 0)] = grp.vertex_indices[iel_grp][1]
- node_tuple_to_coord[(0, 2)] = grp.vertex_indices[iel_grp][2]
- vertex_pair = (min(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][1]), \
- max(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][1]))
- node_tuple_to_coord[(1, 0)] = self.pair_map[vertex_pair].midpoint
- vertex_pair = (min(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][2]), \
- max(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][2]))
- node_tuple_to_coord[(0, 1)] = self.pair_map[vertex_pair].midpoint
- vertex_pair = (min(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][2]), \
- max(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][2]))
- node_tuple_to_coord[(1, 1)] = self.pair_map[vertex_pair].midpoint
- #generate actual elements
- #middle element
- for i in six.moves.range(0, len(self.tri_result[1])):
- groups[grpn][iel_grp][i] = \
- node_tuple_to_coord[self.tri_node_tuples[self.tri_result[1][i]]]
- for i in six.moves.range(0, 4):
- if i == 1:
- continue
- for j in six.moves.range(0, len(self.tri_result[i])):
- groups[grpn][nelements_in_grp][j] = \
- node_tuple_to_coord[self.tri_node_tuples[self.tri_result[i][j]]]
- nelements_in_grp += 1
- '''
- #middle element
- vertex_pair = (min(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][1]), \
- max(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][1]))
- groups[grpn][iel_grp][0] = self.pair_map[vertex_pair].midpoint
- vertex_pair = (min(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][2]), \
- max(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][2]))
- groups[grpn][iel_grp][1] = self.pair_map[vertex_pair].midpoint
- vertex_pair = (min(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][2]), \
- max(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][2]))
- groups[grpn][iel_grp][2] = self.pair_map[vertex_pair].midpoint
- #element 0
- groups[grpn][nelements_in_grp][0] = grp.vertex_indices[iel_grp][0]
- vertex_pair = (min(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][1]), \
- max(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][1]))
- groups[grpn][nelements_in_grp][1] = self.pair_map[vertex_pair].midpoint
- vertex_pair = (min(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][2]), \
- max(grp.vertex_indices[iel_grp][0], grp.vertex_indices[iel_grp][2]))
- groups[grpn][nelements_in_grp][2] = self.pair_map[vertex_pair].midpoint
- nelements_in_grp += 1
- #element 1
- groups[grpn][nelements_in_grp][0] = grp.vertex_indices[iel_grp][1]
- vertex_pair = (min(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][0]), \
- max(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][0]))
- groups[grpn][nelements_in_grp][1] = self.pair_map[vertex_pair].midpoint
- vertex_pair = (min(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][2]), \
- max(grp.vertex_indices[iel_grp][1], grp.vertex_indices[iel_grp][2]))
- groups[grpn][nelements_in_grp][2] = self.pair_map[vertex_pair].midpoint
- nelements_in_grp += 1
- #element 2
- groups[grpn][nelements_in_grp][0] = grp.vertex_indices[iel_grp][2]
- vertex_pair = (min(grp.vertex_indices[iel_grp][2], grp.vertex_indices[iel_grp][0]), \
- max(grp.vertex_indices[iel_grp][2], grp.vertex_indices[iel_grp][0]))
- groups[grpn][nelements_in_grp][1] = self.pair_map[vertex_pair].midpoint
- vertex_pair = (min(grp.vertex_indices[iel_grp][2], grp.vertex_indices[iel_grp][1]), \
- max(grp.vertex_indices[iel_grp][2], grp.vertex_indices[iel_grp][1]))
- groups[grpn][nelements_in_grp][2] = self.pair_map[vertex_pair].midpoint
- nelements_in_grp += 1
- '''
- # }}}
+ for node_index, node_tuple in enumerate(self.index_to_node_tuple[cur_dim]):
+ node_tuple_to_coord[node_tuple] = grp.vertex_indices[iel_grp][node_index]
+ for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]):
+ node_tuple_to_coord[midpoint_tuple] = midpoint_vertices[midpoint_index]
+ for i in range(len(self.simplex_result[cur_dim])):
+ for j in range(len(self.simplex_result[cur_dim][i])):
+ if i == 0:
+ groups[grpn][iel_grp][j] = \
+ node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]]
+ else:
+ groups[grpn][nelements_in_grp+i-1][j] = \
+ node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]]
+ nelements_in_grp += len(self.simplex_result[cur_dim])-1
+ #assuming quad otherwise
+ #else:
+ #quadrilateral
+# node_tuple_to_coord = {}
+# for node_index, node_tuple in enumerate(self.index_to_node_tuple[cur_dim]):
+# node_tuple_to_coord[node_tuple] = grp.vertex_indices[iel_grp][node_index]
+# def generate_all_tuples(cur_list):
+# if len(cur_list[len(cur_list)-1])
+
+
+ #clear connectivity data
+ for grp in self.last_mesh.groups:
+ iel_base = grp.element_nr_base
+ for iel_grp in range(grp.nelements):
+ for i in range(len(grp.vertex_indices[iel_grp])):
+ for j in range(i+1, len(grp.vertex_indices[iel_grp])):
+ min_vert = min(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j])
+ max_vert = max(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j])
+ vertex_pair = (min_vert, max_vert)
+ root_ray = self.get_root(self.pair_map[vertex_pair])
+ if root_ray not in self.seen_tuple:
+ self.seen_tuple[root_ray] = True
+ cur_tree = self.get_subtree(root_ray)
+ for node in cur_tree:
+ node.adjacent_elements = []
+ new_hanging_vertex_element[node.left_vertex] = []
+ new_hanging_vertex_element[node.right_vertex] = []
+
+ self.seen_tuple.clear()
+
+ nelements_in_grp = grp.nelements
+ for grp in groups:
+ for iel_grp in range(len(grp)):
+ add_verts = []
+ for i in range(len(grp[iel_grp])):
+ add_verts.append(grp[iel_grp][i])
+ add_element_to_connectivity(add_verts, new_hanging_vertex_element, iel_base+iel_grp)
+ #assert ray connectivity
+ #check_adjacent_elements(groups, new_hanging_vertex_element, nelements_in_grp)
- grpn += 1
- #print vertices
- #print vertex_indices
+ self.hanging_vertex_element = new_hanging_vertex_element
from meshmode.mesh.generation import make_group_from_vertices
- #grp = make_group_from_vertices(vertices, vertex_indices, 4)
grp = []
- grpn = 0
for grpn in range(0, len(groups)):
grp.append(make_group_from_vertices(vertices, groups[grpn], 4))
from meshmode.mesh import Mesh
- #return Mesh(vertices, [grp], element_connectivity=self.generate_connectivity(len(self.last_mesh.groups[group].vertex_indices) \
- # + count*3))
-
- self.last_mesh = Mesh(vertices, grp, element_connectivity=self.generate_connectivity(totalnelements, nvertices, groups))
- return self.last_mesh
- split_faces = {}
-
- ibase = self.get_refine_base_index()
- affected_group_indices = set()
-
- for grp in self.last_mesh.groups:
- iel_base
+ self.last_mesh = Mesh(vertices, grp, element_connectivity=self.generate_connectivity(
+ totalnelements, nvertices, groups))
+ return self.last_mesh
+ def print_rays(self, ind):
+ for i in range(len(self.last_mesh.groups[0].vertex_indices[ind])):
+ for j in range(i+1, len(self.last_mesh.groups[0].vertex_indices[ind])):
+ mn = min(self.last_mesh.groups[0].vertex_indices[ind][i],
+ self.last_mesh.groups[0].vertex_indices[ind][j])
+ mx = max(self.last_mesh.groups[0].vertex_indices[ind][j],
+ self.last_mesh.groups[0].vertex_indices[ind][i])
+ vertex_pair = (mn, mx)
+ print('LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex)
+ print('RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex)
+ print('ADJACENT:')
+ rays = self.get_leaves(self.pair_map[vertex_pair])
+ for k in rays:
+ print(k.adjacent_elements)
+ '''
+ def print_rays(self, groups, ind):
+ for i in range(len(groups[0][ind])):
+ for j in range(i+1, len(groups[0][ind])):
+ mn = min(groups[0][ind][i], groups[0][ind][j])
+ mx = max(groups[0][ind][i], groups[0][ind][j])
+ vertex_pair = (mn, mx)
+ print('LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex)
+ print('RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex)
+ print('ADJACENT:')
+ rays = self.get_leaves(self.pair_map[vertex_pair])
+ for k in rays:
+ print(k.adjacent_elements)
+ '''
+
+ def print_hanging_elements(self, ind):
+ for i in self.last_mesh.groups[0].vertex_indices[ind]:
+ print("IND:", i, self.hanging_vertex_element[i])
def generate_connectivity(self, nelements, nvertices, groups):
# medium-term FIXME: make this an incremental update
# rather than build-from-scratch
- vertex_to_element = [[] for i in xrange(nvertices)]
-
+ vertex_to_element = [[] for i in range(nvertices)]
element_index = 0
for grp in groups:
- for iel_grp in xrange(len(grp)):
+ for iel_grp in range(len(grp)):
for ivertex in grp[iel_grp]:
vertex_to_element[ivertex].append(element_index)
element_index += 1
- element_to_element = [set() for i in xrange(nelements)]
+ element_to_element = [set() for i in range(nelements)]
element_index = 0
- for grp in groups:
- for iel_grp in xrange(len(grp)):
- for ivertex in grp[iel_grp]:
- element_to_element[element_index].update(
- vertex_to_element[ivertex])
- element_index += 1
- #print self.ray_elements
+ if self.lazy:
+ for grp in groups:
+ for iel_grp in range(len(grp)):
+ for i in range(len(grp[iel_grp])):
+ for j in range(i+1, len(grp[iel_grp])):
+ vertex_pair = (min(grp[iel_grp][i], grp[iel_grp][j]), max(grp[iel_grp][i], grp[iel_grp][j]))
+ #print 'iel:', iel_grp, 'pair:', vertex_pair
+ if vertex_pair not in self.seen_tuple:
+ self.propagate_tree(self.get_root(self.pair_map[vertex_pair]), self.hanging_vertex_element, element_to_element)
+ #print self.pair_map[vertex_pair].left_vertex, self.pair_map[vertex_pair].right_vertex, self.pair_map[vertex_pair].adjacent_elements, self.hanging_vertex_element[self.pair_map[vertex_pair].left_vertex], self.hanging_vertex_element[self.pair_map[vertex_pair].right_vertex]
+
+ else:
+ for grp in groups:
+ for iel_grp in range(len(grp)):
+ for ivertex in grp[iel_grp]:
+ element_to_element[element_index].update(
+ vertex_to_element[ivertex])
+ if self.hanging_vertex_element[ivertex]:
+ for hanging_element in self.hanging_vertex_element[ivertex]:
+ if element_index != hanging_element:
+ element_to_element[element_index].update([hanging_element])
+ element_to_element[hanging_element].update([element_index])
+ for i in range(len(grp[iel_grp])):
+ for j in range(i+1, len(grp[iel_grp])):
+ vertex_pair = (min(grp[iel_grp][i], grp[iel_grp][j]),
+ max(grp[iel_grp][i], grp[iel_grp][j]))
+ #element_to_element[element_index].update(
+ #self.pair_map[vertex_pair].adjacent_elements)
+ queue = [self.pair_map[vertex_pair]]
+ while queue:
+ vertex = queue.pop(0)
+ #if leaf node
+ if vertex.left is None and vertex.right is None:
+ assert(element_index in vertex.adjacent_elements)
+ element_to_element[element_index].update(
+ vertex.adjacent_elements)
+ else:
+ queue.append(vertex.left)
+ queue.append(vertex.right)
+ '''
+ if self.hanging_vertex_element[ivertex] and element_index != self.hanging_vertex_element[ivertex][0]:
+ element_to_element[element_index].update([self.hanging_vertex_element[ivertex][0]])
+ element_to_element[self.hanging_vertex_element[ivertex][0]].update([element_index])
+ '''
+ element_index += 1
+ print(len(element_to_element))
+ for iel, neighbors in enumerate(element_to_element):
+ if iel in neighbors:
+ neighbors.remove(iel)
+ #print(self.ray_elements)
+ '''
for ray in self.rays:
curnode = ray.first
while curnode is not None:
- '''
if len(curnode.value.elements) >= 2:
if curnode.value.elements[0] is not None:
element_to_element[curnode.value.elements[0]].update(curnode.value.elements)
if curnode.value.elements[1] is not None:
element_to_element[curnode.value.elements[1]].update(curnode.value.elements)
- '''
if len(curnode.value.velements) >= 2:
if curnode.value.velements[0] is not None:
element_to_element[curnode.value.velements[0]].update(curnode.value.velements)
if curnode.value.velements[1] is not None:
element_to_element[curnode.value.velements[1]].update(curnode.value.velements)
curnode = curnode.next
-
+ '''
'''
for i in self.ray_elements:
for j in i:
@@ -717,5 +783,3 @@ class Refiner(object):
# vim: foldmethod=marker
-# vim: shiftwidth=4
-# vim: softtabstop=4
diff --git a/meshmode/mesh/tesselate.py b/meshmode/mesh/tesselate.py
new file mode 100644
index 0000000000000000000000000000000000000000..4dbfefab948decd297ad00abc303d5cf878f2054
--- /dev/null
+++ b/meshmode/mesh/tesselate.py
@@ -0,0 +1,73 @@
+from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \
+ as gnitstam
+
+
+def add_tuples(a, b):
+ return tuple(ac+bc for ac, bc in zip(a, b))
+
+
+def tesselatetri():
+ result = []
+
+ node_tuples = list(gnitstam(2, 2))
+ node_dict = dict(
+ (ituple, idx)
+ for idx, ituple in enumerate(node_tuples))
+
+ def try_add_tri(current, d1, d2, d3):
+ try:
+ result.append((
+ node_dict[add_tuples(current, d1)],
+ node_dict[add_tuples(current, d2)],
+ node_dict[add_tuples(current, d3)],
+ ))
+ except KeyError:
+ pass
+
+ if len(result) > 0:
+ return [node_tuples, result]
+ for current in node_tuples:
+ # this is a tesselation of a cube into six tets.
+ # subtets that fall outside of the master tet are simply not added.
+
+ # positively oriented
+ try_add_tri(current, (0, 0), (1, 0), (0, 1))
+ try_add_tri(current, (1, 0), (1, 1), (0, 1))
+ return [node_tuples, result]
+
+
+def tesselatetet():
+ node_tuples = list(gnitstam(2, 3))
+
+ node_dict = dict(
+ (ituple, idx)
+ for idx, ituple in enumerate(node_tuples))
+
+ def try_add_tet(current, d1, d2, d3, d4):
+ try:
+ result.append((
+ node_dict[add_tuples(current, d1)],
+ node_dict[add_tuples(current, d2)],
+ node_dict[add_tuples(current, d3)],
+ node_dict[add_tuples(current, d4)],
+ ))
+ except KeyError:
+ pass
+
+ result = []
+
+ if len(result) > 0:
+ return [node_tuples, result]
+ for current in node_tuples:
+ # this is a tesselation of a cube into six tets.
+ # subtets that fall outside of the master tet are simply not added.
+
+ # positively oriented
+ try_add_tet(current, (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1))
+ try_add_tet(current, (1, 0, 1), (1, 0, 0), (0, 0, 1), (0, 1, 0))
+ try_add_tet(current, (1, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 1))
+
+ try_add_tet(current, (1, 0, 0), (0, 1, 0), (1, 0, 1), (1, 1, 0))
+ try_add_tet(current, (0, 1, 1), (0, 1, 0), (1, 1, 0), (1, 0, 1))
+ try_add_tet(current, (0, 1, 1), (1, 1, 1), (1, 0, 1), (1, 1, 0))
+ return [node_tuples, result]
diff --git a/meshmode/mesh/tools.py b/meshmode/mesh/tools.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b46bb83d6022710aed5505c3c21646929de8426
--- /dev/null
+++ b/meshmode/mesh/tools.py
@@ -0,0 +1,95 @@
+from __future__ import division
+from __future__ import absolute_import
+from six.moves import range
+
+__copyright__ = "Copyright (C) 2010,2012,2013 Andreas Kloeckner, Michael Tom"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import numpy as np
+from pytools.spatial_btree import SpatialBinaryTreeBucket
+
+
+# {{{ make_element_lookup_tree
+
+def make_element_lookup_tree(mesh, eps=1e-12):
+ from meshmode.mesh.processing import find_bounding_box
+ bbox_min, bbox_max = find_bounding_box(mesh)
+ bbox_min -= eps
+ bbox_max += eps
+
+ tree = SpatialBinaryTreeBucket(bbox_min, bbox_max)
+
+ for igrp, grp in enumerate(mesh.groups):
+ for iel_grp in range(grp.nelements):
+ el_vertices = mesh.vertices[:, grp.vertex_indices[iel_grp]]
+
+ el_bbox_min = np.min(el_vertices, axis=-1) - eps
+ el_bbox_max = np.max(el_vertices, axis=-1) + eps
+
+ tree.insert((igrp, iel_grp), (el_bbox_min, el_bbox_max))
+
+ return tree
+
+# }}}
+
+
+# {{{ nd_quad_submesh
+
+def nd_quad_submesh(node_tuples):
+ """Return a list of tuples of indices into the node list that
+ generate a tesselation of the reference element.
+
+ :arg node_tuples: A list of tuples *(i, j, ...)* of integers
+ indicating node positions inside the unit element. The
+ returned list references indices in this list.
+
+ :func:`pytools.generate_nonnegative_integer_tuples_below`
+ may be used to generate *node_tuples*.
+
+ See also :func:`modepy.tools.simplex_submesh`.
+ """
+
+ from pytools import single_valued, add_tuples
+ dims = single_valued(len(nt) for nt in node_tuples)
+
+ node_dict = dict(
+ (ituple, idx)
+ for idx, ituple in enumerate(node_tuples))
+
+ from pytools import generate_nonnegative_integer_tuples_below as gnitb
+
+ result = []
+ for current in node_tuples:
+ try:
+ result.append(tuple(
+ node_dict[add_tuples(current, offset)]
+ for offset in gnitb(2, dims)))
+
+ except KeyError:
+ pass
+
+ return result
+
+# }}}
+
+
+# vim: foldmethod=marker
diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py
index 8966861b168f5d740a099e5ec4ad80baccdba4b0..caa9c1beb3cdb7123a3c9798ce04a3c882b14f93 100644
--- a/meshmode/mesh/visualization.py
+++ b/meshmode/mesh/visualization.py
@@ -30,7 +30,8 @@ import numpy as np
# {{{ draw_2d_mesh
def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True,
- draw_connectivity=False, **kwargs):
+ draw_nodal_adjacency=False, draw_face_numbers=False,
+ set_bounding_box=False, **kwargs):
assert mesh.ambient_dim == 2
import matplotlib.pyplot as pt
@@ -74,7 +75,7 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True,
ha="center", va="center", color="blue",
bbox=dict(facecolor='white', alpha=0.5, lw=0))
- if draw_connectivity:
+ if draw_nodal_adjacency:
def global_iel_to_group_and_iel(global_iel):
for igrp, grp in enumerate(mesh.groups):
if global_iel < grp.nelements:
@@ -83,7 +84,7 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True,
raise ValueError("invalid element nr")
- cnx = mesh.element_connectivity
+ cnx = mesh.nodal_adjacency
nb_starts = cnx.neighbors_starts
for iel_g in range(mesh.nelements):
@@ -110,6 +111,42 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True,
length_includes_head=True,
color="green", head_width=1e-2, lw=1e-2)
+ if draw_face_numbers:
+ for igrp, grp in enumerate(mesh.groups):
+ for iel, el in enumerate(grp.vertex_indices):
+ elverts = mesh.vertices[:, el]
+ el_center = np.mean(elverts, axis=-1)
+
+ for iface, fvi in enumerate(grp.face_vertex_indices()):
+ face_center = (
+ 0.3*el_center
+ +
+ 0.7*np.mean(elverts[:, fvi], axis=-1))
+
+ pt.text(face_center[0], face_center[1], str(iface), fontsize=12,
+ ha="center", va="center", color="purple",
+ bbox=dict(facecolor='white', alpha=0.5, lw=0))
+
+ if set_bounding_box:
+ from meshmode.mesh.processing import find_bounding_box
+ lower, upper = find_bounding_box(mesh)
+ pt.xlim([lower[0], upper[0]])
+ pt.ylim([lower[1], upper[1]])
+
+# }}}
+
+
+# {{{ draw_curve
+
+def draw_curve(mesh):
+ import matplotlib.pyplot as pt
+ pt.plot(mesh.vertices[0], mesh.vertices[1], "o")
+
+ for i, group in enumerate(mesh.groups):
+ pt.plot(
+ group.nodes[0].ravel(),
+ group.nodes[1].ravel(), "-x", label="Group %d" % i)
+
# }}}
# vim: foldmethod=marker
diff --git a/meshmode/version.py b/meshmode/version.py
index 92c29262fa4d3e5d7e00f6e6bef07af819a6d872..1da70a1e09bc16f4b7e677505af6cf4dd4502b86 100644
--- a/meshmode/version.py
+++ b/meshmode/version.py
@@ -1,2 +1,2 @@
-VERSION = (2014, 1)
+VERSION = (2015, 1)
VERSION_TEXT = ".".join(str(i) for i in VERSION)
diff --git a/requirements.txt b/requirements.txt
index a6902dd2d88736406f611ddf5d76b324987f5dcb..8e99cd7fdd25e578a5dc11c56f4c43506d18974f 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,6 +1,7 @@
numpy
git+git://github.com/inducer/modepy
git+git://github.com/pyopencl/pyopencl
+git+git://github.com/inducer/islpy
git+git://github.com/inducer/loopy
# required by pytential, which is in turn needed for some tests
diff --git a/setup.cfg b/setup.cfg
index 6faef2e65abe138f9ab22c94452c43be0e52c075..512924240ea5f744029e9dd0395769ca66267ca1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,3 +1,3 @@
[flake8]
-ignore = E126,E127,E128,E123,E226,E241,E242,E265
+ignore = E126,E127,E128,E123,E226,E241,E242,E265,W503,E402
max-line-length=85
diff --git a/test/circle.step b/test/circle.step
new file mode 100644
index 0000000000000000000000000000000000000000..78f370557b0fd9f525f484a7d6cf91ba1ea78e29
--- /dev/null
+++ b/test/circle.step
@@ -0,0 +1,85 @@
+ISO-10303-21;
+HEADER;
+FILE_DESCRIPTION(('FreeCAD Model'),'2;1');
+FILE_NAME('/home/andreas/own/graphics/freecad/circle.step',
+ '2015-04-28T14:41:10',('FreeCAD'),('FreeCAD'),
+ 'Open CASCADE STEP processor 6.7','FreeCAD','Unknown');
+FILE_SCHEMA(('AUTOMOTIVE_DESIGN_CC2 { 1 2 10303 214 -1 1 5 4 }'));
+ENDSEC;
+DATA;
+#1 = APPLICATION_PROTOCOL_DEFINITION('committee draft',
+ 'automotive_design',1997,#2);
+#2 = APPLICATION_CONTEXT(
+ 'core data for automotive mechanical design processes');
+#3 = SHAPE_DEFINITION_REPRESENTATION(#4,#10);
+#4 = PRODUCT_DEFINITION_SHAPE('','',#5);
+#5 = PRODUCT_DEFINITION('design','',#6,#9);
+#6 = PRODUCT_DEFINITION_FORMATION('','',#7);
+#7 = PRODUCT('Circle001','Circle001','',(#8));
+#8 = MECHANICAL_CONTEXT('',#2,'mechanical');
+#9 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design');
+#10 = MANIFOLD_SURFACE_SHAPE_REPRESENTATION('',(#11,#15),#46);
+#11 = AXIS2_PLACEMENT_3D('',#12,#13,#14);
+#12 = CARTESIAN_POINT('',(0.,0.,0.));
+#13 = DIRECTION('',(0.,0.,1.));
+#14 = DIRECTION('',(1.,0.,-0.));
+#15 = SHELL_BASED_SURFACE_MODEL('',(#16));
+#16 = OPEN_SHELL('',(#17));
+#17 = ADVANCED_FACE('',(#18),#31,.T.);
+#18 = FACE_BOUND('',#19,.F.);
+#19 = EDGE_LOOP('',(#20));
+#20 = ORIENTED_EDGE('',*,*,#21,.T.);
+#21 = EDGE_CURVE('',#22,#22,#24,.T.);
+#22 = VERTEX_POINT('',#23);
+#23 = CARTESIAN_POINT('',(2.E-04,0.,0.));
+#24 = SURFACE_CURVE('',#25,(#30),.PCURVE_S1.);
+#25 = CIRCLE('',#26,2.E-04);
+#26 = AXIS2_PLACEMENT_3D('',#27,#28,#29);
+#27 = CARTESIAN_POINT('',(0.,0.,0.));
+#28 = DIRECTION('',(0.,0.,1.));
+#29 = DIRECTION('',(1.,0.,0.));
+#30 = PCURVE('',#31,#36);
+#31 = PLANE('',#32);
+#32 = AXIS2_PLACEMENT_3D('',#33,#34,#35);
+#33 = CARTESIAN_POINT('',(2.E-04,0.,0.));
+#34 = DIRECTION('',(0.,0.,-1.));
+#35 = DIRECTION('',(-1.,0.,0.));
+#36 = DEFINITIONAL_REPRESENTATION('',(#37),#45);
+#37 = ( BOUNDED_CURVE() B_SPLINE_CURVE(2,(#38,#39,#40,#41,#42,#43,#44),
+.UNSPECIFIED.,.F.,.F.) B_SPLINE_CURVE_WITH_KNOTS((1,2,2,2,2,1),(
+ -2.094395102393,0.,2.094395102393,4.188790204786,6.28318530718,
+8.377580409573),.UNSPECIFIED.) CURVE() GEOMETRIC_REPRESENTATION_ITEM()
+RATIONAL_B_SPLINE_CURVE((1.,0.5,1.,0.5,1.,0.5,1.)) REPRESENTATION_ITEM(
+ '') );
+#38 = CARTESIAN_POINT('',(0.,0.));
+#39 = CARTESIAN_POINT('',(0.,3.464101615138E-04));
+#40 = CARTESIAN_POINT('',(3.E-04,1.732050807569E-04));
+#41 = CARTESIAN_POINT('',(6.E-04,4.898587196589E-20));
+#42 = CARTESIAN_POINT('',(3.E-04,-1.732050807569E-04));
+#43 = CARTESIAN_POINT('',(2.981555974335E-19,-3.464101615138E-04));
+#44 = CARTESIAN_POINT('',(0.,0.));
+#45 = ( GEOMETRIC_REPRESENTATION_CONTEXT(2)
+PARAMETRIC_REPRESENTATION_CONTEXT() REPRESENTATION_CONTEXT('2D SPACE',''
+ ) );
+#46 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3)
+GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#50)) GLOBAL_UNIT_ASSIGNED_CONTEXT(
+(#47,#48,#49)) REPRESENTATION_CONTEXT('Context #1',
+ '3D Context with UNIT and UNCERTAINTY') );
+#47 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT($,.METRE.) );
+#48 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) );
+#49 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() );
+#50 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-10),#47,
+ 'distance_accuracy_value','confusion accuracy');
+#51 = PRODUCT_TYPE('part',$,(#7));
+#52 = MECHANICAL_DESIGN_GEOMETRIC_PRESENTATION_REPRESENTATION('',(#53),
+ #46);
+#53 = STYLED_ITEM('color',(#54),#17);
+#54 = PRESENTATION_STYLE_ASSIGNMENT((#55));
+#55 = SURFACE_STYLE_USAGE(.BOTH.,#56);
+#56 = SURFACE_SIDE_STYLE('',(#57));
+#57 = SURFACE_STYLE_FILL_AREA(#58);
+#58 = FILL_AREA_STYLE('',(#59));
+#59 = FILL_AREA_STYLE_COLOUR('',#60);
+#60 = COLOUR_RGB('',0.800000011921,0.800000011921,0.800000011921);
+ENDSEC;
+END-ISO-10303-21;
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 6047113c956d1f616c57fc6610b8a75393b90db0..c7c2bc488932feec285147eaf2fc7fa7358deab4 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -1,7 +1,4 @@
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-from six.moves import range
+from __future__ import division, absolute_import, print_function
__copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
@@ -25,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+from six.moves import range
import numpy as np
import numpy.linalg as la
import pyopencl as cl
@@ -35,58 +33,342 @@ from pyopencl.tools import ( # noqa
pytest_generate_tests_for_pyopencl
as pytest_generate_tests)
+from meshmode.discretization.poly_element import (
+ InterpolatoryQuadratureSimplexGroupFactory,
+ PolynomialWarpAndBlendGroupFactory,
+ PolynomialEquidistantGroupFactory,
+ )
+from meshmode.mesh import BTAG_ALL
+from meshmode.discretization.connection import \
+ FRESTR_ALL_FACES, FRESTR_INTERIOR_FACES
+
import pytest
import logging
logger = logging.getLogger(__name__)
-def test_boundary_interpolation(ctx_getter):
+# {{{ circle mesh
+
+def test_circle_mesh(do_plot=False):
+ from meshmode.mesh.io import generate_gmsh, FileSource
+ print("BEGIN GEN")
+ mesh = generate_gmsh(
+ FileSource("circle.step"), 2, order=2,
+ force_ambient_dim=2,
+ other_options=[
+ "-string", "Mesh.CharacteristicLengthMax = 0.05;"]
+ )
+ print("END GEN")
+ print(mesh.nelements)
+
+ from meshmode.mesh.processing import affine_map
+ mesh = affine_map(mesh, A=3*np.eye(2))
+
+ if do_plot:
+ from meshmode.mesh.visualization import draw_2d_mesh
+ draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True,
+ set_bounding_box=True)
+ import matplotlib.pyplot as pt
+ pt.show()
+
+# }}}
+
+
+# {{{ convergence of boundary interpolation
+
+@pytest.mark.parametrize("group_factory", [
+ InterpolatoryQuadratureSimplexGroupFactory,
+ PolynomialWarpAndBlendGroupFactory
+ ])
+@pytest.mark.parametrize("boundary_tag", [
+ BTAG_ALL,
+ FRESTR_ALL_FACES,
+ FRESTR_INTERIOR_FACES,
+ ])
+@pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
+ ("blob", 2, [1e-1, 8e-2, 5e-2]),
+ ("warp", 2, [10, 20, 30]),
+ ("warp", 3, [10, 20, 30]),
+ ])
+@pytest.mark.parametrize("per_face_groups", [False, True])
+def test_boundary_interpolation(ctx_getter, group_factory, boundary_tag,
+ mesh_name, dim, mesh_pars, per_face_groups):
cl_ctx = ctx_getter()
queue = cl.CommandQueue(cl_ctx)
- from meshmode.mesh.io import generate_gmsh, FileSource
from meshmode.discretization import Discretization
- from meshmode.discretization.poly_element import \
- InterpolatoryQuadratureSimplexGroupFactory
- from meshmode.discretization.connection import make_boundary_restriction
+ from meshmode.discretization.connection import (
+ make_face_restriction, check_connection)
from pytools.convergence import EOCRecorder
eoc_rec = EOCRecorder()
order = 4
- for h in [1e-1, 3e-2, 1e-2]:
- print("BEGIN GEN")
- mesh = generate_gmsh(
- FileSource("blob-2d.step"), 2, order=order,
- force_ambient_dim=2,
- other_options=[
- "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
- )
- print("END GEN")
+ def f(x):
+ return 0.1*cl.clmath.sin(30*x)
+
+ for mesh_par in mesh_pars:
+ # {{{ get mesh
+
+ if mesh_name == "blob":
+ assert dim == 2
+
+ h = mesh_par
+
+ from meshmode.mesh.io import generate_gmsh, FileSource
+ print("BEGIN GEN")
+ mesh = generate_gmsh(
+ FileSource("blob-2d.step"), 2, order=order,
+ force_ambient_dim=2,
+ other_options=[
+ "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
+ )
+ print("END GEN")
+ elif mesh_name == "warp":
+ from meshmode.mesh.generation import generate_warped_rect_mesh
+ mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par)
+
+ h = 1/mesh_par
+ else:
+ raise ValueError("mesh_name not recognized")
+
+ # }}}
vol_discr = Discretization(cl_ctx, mesh,
- InterpolatoryQuadratureSimplexGroupFactory(order))
+ group_factory(order))
print("h=%s -> %d elements" % (
h, sum(mgrp.nelements for mgrp in mesh.groups)))
x = vol_discr.nodes()[0].with_queue(queue)
- f = 0.1*cl.clmath.sin(30*x)
+ vol_f = f(x)
- bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
- queue, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(order))
+ bdry_connection = make_face_restriction(
+ vol_discr, group_factory(order),
+ boundary_tag, per_face_groups=per_face_groups)
+ check_connection(bdry_connection)
+ bdry_discr = bdry_connection.to_discr
bdry_x = bdry_discr.nodes()[0].with_queue(queue)
- bdry_f = 0.1*cl.clmath.sin(30*bdry_x)
- bdry_f_2 = bdry_connection(queue, f)
+ bdry_f = f(bdry_x)
+ bdry_f_2 = bdry_connection(queue, vol_f)
+
+ if mesh_name == "blob" and dim == 2:
+ mat = bdry_connection.full_resample_matrix(queue).get(queue)
+ bdry_f_2_by_mat = mat.dot(vol_f.get())
+
+ mat_error = la.norm(bdry_f_2.get(queue=queue) - bdry_f_2_by_mat)
+ assert mat_error < 1e-14, mat_error
err = la.norm((bdry_f-bdry_f_2).get(), np.inf)
eoc_rec.add_data_point(h, err)
print(eoc_rec)
- assert eoc_rec.order_estimate() >= order-0.5
+ assert (
+ eoc_rec.order_estimate() >= order-0.5
+ or eoc_rec.max_error() < 1e-14)
+
+# }}}
+
+
+# {{{ boundary-to-all-faces connecttion
+
+@pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
+ ("blob", 2, [1e-1, 8e-2, 5e-2]),
+ ("warp", 2, [10, 20, 30]),
+ ("warp", 3, [10, 20, 30]),
+ ])
+@pytest.mark.parametrize("per_face_groups", [False, True])
+def test_all_faces_interpolation(ctx_getter, mesh_name, dim, mesh_pars,
+ per_face_groups):
+ cl_ctx = ctx_getter()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.connection import (
+ make_face_restriction, make_face_to_all_faces_embedding,
+ check_connection)
+
+ from pytools.convergence import EOCRecorder
+ eoc_rec = EOCRecorder()
+
+ order = 4
+
+ def f(x):
+ return 0.1*cl.clmath.sin(30*x)
+
+ for mesh_par in mesh_pars:
+ # {{{ get mesh
+
+ if mesh_name == "blob":
+ assert dim == 2
+
+ h = mesh_par
+
+ from meshmode.mesh.io import generate_gmsh, FileSource
+ print("BEGIN GEN")
+ mesh = generate_gmsh(
+ FileSource("blob-2d.step"), 2, order=order,
+ force_ambient_dim=2,
+ other_options=[
+ "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
+ )
+ print("END GEN")
+ elif mesh_name == "warp":
+ from meshmode.mesh.generation import generate_warped_rect_mesh
+ mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par)
+
+ h = 1/mesh_par
+ else:
+ raise ValueError("mesh_name not recognized")
+
+ # }}}
+
+ vol_discr = Discretization(cl_ctx, mesh,
+ PolynomialWarpAndBlendGroupFactory(order))
+ print("h=%s -> %d elements" % (
+ h, sum(mgrp.nelements for mgrp in mesh.groups)))
+
+ all_face_bdry_connection = make_face_restriction(
+ vol_discr, PolynomialWarpAndBlendGroupFactory(order),
+ FRESTR_ALL_FACES, per_face_groups=per_face_groups)
+ all_face_bdry_discr = all_face_bdry_connection.to_discr
+
+ for ito_grp, ceg in enumerate(all_face_bdry_connection.groups):
+ for ibatch, batch in enumerate(ceg.batches):
+ assert np.array_equal(
+ batch.from_element_indices.get(queue),
+ np.arange(vol_discr.mesh.nelements))
+
+ if per_face_groups:
+ assert ito_grp == batch.to_element_face
+ else:
+ assert ibatch == batch.to_element_face
+
+ all_face_x = all_face_bdry_discr.nodes()[0].with_queue(queue)
+ all_face_f = f(all_face_x)
+
+ all_face_f_2 = all_face_bdry_discr.zeros(queue)
+
+ for boundary_tag in [
+ BTAG_ALL,
+ FRESTR_INTERIOR_FACES,
+ ]:
+ bdry_connection = make_face_restriction(
+ vol_discr, PolynomialWarpAndBlendGroupFactory(order),
+ boundary_tag, per_face_groups=per_face_groups)
+ bdry_discr = bdry_connection.to_discr
+
+ bdry_x = bdry_discr.nodes()[0].with_queue(queue)
+ bdry_f = f(bdry_x)
+
+ all_face_embedding = make_face_to_all_faces_embedding(
+ bdry_connection, all_face_bdry_discr)
+
+ check_connection(all_face_embedding)
+
+ all_face_f_2 += all_face_embedding(queue, bdry_f)
+
+ err = la.norm((all_face_f-all_face_f_2).get(), np.inf)
+ eoc_rec.add_data_point(h, err)
+
+ print(eoc_rec)
+ assert (
+ eoc_rec.order_estimate() >= order-0.5
+ or eoc_rec.max_error() < 1e-14)
+
+# }}}
+
+
+# {{{ convergence of opposite-face interpolation
+
+@pytest.mark.parametrize("group_factory", [
+ InterpolatoryQuadratureSimplexGroupFactory,
+ PolynomialWarpAndBlendGroupFactory
+ ])
+@pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
+ ("blob", 2, [1e-1, 8e-2, 5e-2]),
+ ("warp", 2, [3, 5, 7]),
+ ("warp", 3, [3, 5]),
+ ])
+def test_opposite_face_interpolation(ctx_getter, group_factory,
+ mesh_name, dim, mesh_pars):
+ logging.basicConfig(level=logging.INFO)
+
+ cl_ctx = ctx_getter()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.connection import (
+ make_face_restriction, make_opposite_face_connection,
+ check_connection)
+
+ from pytools.convergence import EOCRecorder
+ eoc_rec = EOCRecorder()
+
+ order = 5
+
+ def f(x):
+ return 0.1*cl.clmath.sin(30*x)
+
+ for mesh_par in mesh_pars:
+ # {{{ get mesh
+
+ if mesh_name == "blob":
+ assert dim == 2
+
+ h = mesh_par
+
+ from meshmode.mesh.io import generate_gmsh, FileSource
+ print("BEGIN GEN")
+ mesh = generate_gmsh(
+ FileSource("blob-2d.step"), 2, order=order,
+ force_ambient_dim=2,
+ other_options=[
+ "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
+ )
+ print("END GEN")
+ elif mesh_name == "warp":
+ from meshmode.mesh.generation import generate_warped_rect_mesh
+ mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par)
+ h = 1/mesh_par
+ else:
+ raise ValueError("mesh_name not recognized")
+
+ # }}}
+
+ vol_discr = Discretization(cl_ctx, mesh,
+ group_factory(order))
+ print("h=%s -> %d elements" % (
+ h, sum(mgrp.nelements for mgrp in mesh.groups)))
+
+ bdry_connection = make_face_restriction(
+ vol_discr, group_factory(order),
+ FRESTR_INTERIOR_FACES)
+ bdry_discr = bdry_connection.to_discr
+
+ opp_face = make_opposite_face_connection(bdry_connection)
+ check_connection(opp_face)
+
+ bdry_x = bdry_discr.nodes()[0].with_queue(queue)
+ bdry_f = f(bdry_x)
+
+ bdry_f_2 = opp_face(queue, bdry_f)
+
+ err = la.norm((bdry_f-bdry_f_2).get(), np.inf)
+ eoc_rec.add_data_point(h, err)
+
+ print(eoc_rec)
+ assert (
+ eoc_rec.order_estimate() >= order-0.5
+ or eoc_rec.max_error() < 1e-13)
+
+# }}}
+
+
+# {{{ element orientation
def test_element_orientation():
from meshmode.mesh.io import generate_gmsh, FileSource
@@ -116,6 +398,10 @@ def test_element_orientation():
assert ((mesh_orient < 0) == (flippy > 0)).all()
+# }}}
+
+
+# {{{ merge and map
def test_merge_and_map(ctx_getter, visualize=False):
from meshmode.mesh.io import generate_gmsh, FileSource
@@ -128,10 +414,10 @@ def test_merge_and_map(ctx_getter, visualize=False):
other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"]
)
- from meshmode.mesh.processing import merge_dijsoint_meshes, affine_map
+ from meshmode.mesh.processing import merge_disjoint_meshes, affine_map
mesh2 = affine_map(mesh, A=np.eye(2), b=np.array([5, 0]))
- mesh3 = merge_dijsoint_meshes((mesh2, mesh))
+ mesh3 = merge_disjoint_meshes((mesh2, mesh))
if visualize:
from meshmode.discretization import Discretization
@@ -147,6 +433,10 @@ def test_merge_and_map(ctx_getter, visualize=False):
vis = make_visualizer(queue, discr, 1)
vis.write_vtk_file("merged.vtu", [])
+# }}}
+
+
+# {{{ sanity checks: single element
@pytest.mark.parametrize("dim", [2, 3])
@pytest.mark.parametrize("order", [1, 3])
@@ -156,21 +446,21 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False):
cl_ctx = ctx_getter()
queue = cl.CommandQueue(cl_ctx)
- from modepy.tools import UNIT_VERTICES
- vertices = UNIT_VERTICES[dim].T.copy()
+ from modepy.tools import unit_vertices
+ vertices = unit_vertices(dim).T.copy()
center = np.empty(dim, np.float64)
center.fill(-0.5)
import modepy as mp
- from meshmode.mesh import SimplexElementGroup, Mesh
+ from meshmode.mesh import SimplexElementGroup, Mesh, BTAG_ALL
mg = SimplexElementGroup(
order=order,
vertex_indices=np.arange(dim+1, dtype=np.int32).reshape(1, -1),
nodes=mp.warp_and_blend_nodes(dim, order).reshape(dim, 1, -1),
dim=dim)
- mesh = Mesh(vertices, [mg])
+ mesh = Mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None)
from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import \
@@ -198,9 +488,11 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False):
# {{{ boundary discretization
- from meshmode.discretization.connection import make_boundary_restriction
- bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
- queue, vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3))
+ from meshmode.discretization.connection import make_face_restriction
+ bdry_connection = make_face_restriction(
+ vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3),
+ BTAG_ALL)
+ bdry_discr = bdry_connection.to_discr
# }}}
@@ -229,6 +521,60 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False):
assert normal_outward_check.get().all(), normal_outward_check.get()
+# }}}
+
+
+# {{{ sanity check: volume interpolation on scipy/qhull delaunay meshes in nD
+
+@pytest.mark.parametrize("dim", [2, 3, 4])
+@pytest.mark.parametrize("order", [3])
+def test_sanity_qhull_nd(ctx_getter, dim, order):
+ pytest.importorskip("scipy")
+
+ logging.basicConfig(level=logging.INFO)
+
+ ctx = ctx_getter()
+ queue = cl.CommandQueue(ctx)
+
+ from scipy.spatial import Delaunay
+ verts = np.random.rand(1000, dim)
+ dtri = Delaunay(verts)
+
+ from meshmode.mesh.io import from_vertices_and_simplices
+ mesh = from_vertices_and_simplices(dtri.points.T, dtri.simplices,
+ fix_orientation=True)
+
+ from meshmode.discretization import Discretization
+ low_discr = Discretization(ctx, mesh,
+ PolynomialEquidistantGroupFactory(order))
+ high_discr = Discretization(ctx, mesh,
+ PolynomialEquidistantGroupFactory(order+1))
+
+ from meshmode.discretization.connection import make_same_mesh_connection
+ cnx = make_same_mesh_connection(high_discr, low_discr)
+
+ def f(x):
+ return 0.1*cl.clmath.sin(x)
+
+ x_low = low_discr.nodes()[0].with_queue(queue)
+ f_low = f(x_low)
+
+ x_high = high_discr.nodes()[0].with_queue(queue)
+ f_high_ref = f(x_high)
+
+ f_high_num = cnx(queue, f_low)
+
+ err = (f_high_ref-f_high_num).get()
+
+ err = la.norm(err, np.inf)/la.norm(f_high_ref.get(), np.inf)
+
+ print(err)
+ assert err < 1e-2
+
+# }}}
+
+
+# {{{ sanity checks: ball meshes
# python test_meshmode.py 'test_sanity_balls(cl._csc, "disk-radius-1.step", 2, 2, visualize=True)' # noqa
@pytest.mark.parametrize(("src_file", "dim"), [
@@ -266,15 +612,15 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
# {{{ discretizations and connections
from meshmode.discretization import Discretization
- from meshmode.discretization.poly_element import \
- InterpolatoryQuadratureSimplexGroupFactory
vol_discr = Discretization(ctx, mesh,
InterpolatoryQuadratureSimplexGroupFactory(quad_order))
- from meshmode.discretization.connection import make_boundary_restriction
- bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
- queue, vol_discr,
- InterpolatoryQuadratureSimplexGroupFactory(quad_order))
+ from meshmode.discretization.connection import make_face_restriction
+ bdry_connection = make_face_restriction(
+ vol_discr,
+ InterpolatoryQuadratureSimplexGroupFactory(quad_order),
+ BTAG_ALL)
+ bdry_discr = bdry_connection.to_discr
# }}}
@@ -344,6 +690,10 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
print(surf_eoc_rec)
assert surf_eoc_rec.order_estimate() >= mesh_order
+# }}}
+
+
+# {{{ rect/box mesh generation
def test_rect_mesh(do_plot=False):
from meshmode.mesh.generation import generate_regular_rect_mesh
@@ -351,11 +701,250 @@ def test_rect_mesh(do_plot=False):
if do_plot:
from meshmode.mesh.visualization import draw_2d_mesh
- draw_2d_mesh(mesh, fill=None, draw_connectivity=True)
+ draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True)
import matplotlib.pyplot as pt
pt.show()
+def test_box_mesh(ctx_getter, visualize=False):
+ from meshmode.mesh.generation import generate_box_mesh
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),))
+
+ if visualize:
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory
+ cl_ctx = ctx_getter()
+ queue = cl.CommandQueue(cl_ctx)
+
+ discr = Discretization(cl_ctx, mesh,
+ PolynomialWarpAndBlendGroupFactory(1))
+
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, 1)
+ vis.write_vtk_file("box.vtu", [])
+
+# }}}
+
+
+# {{{ as_python stringification
+
+def test_as_python():
+ from meshmode.mesh.generation import generate_box_mesh
+ mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),))
+
+ # These implicitly compute these adjacency structures.
+ mesh.nodal_adjacency
+ mesh.facial_adjacency_groups
+
+ from meshmode.mesh import as_python
+ code = as_python(mesh)
+
+ print(code)
+ exec_dict = {}
+ exec(compile(code, "gen_code.py", "exec"), exec_dict)
+
+ mesh_2 = exec_dict["make_mesh"]()
+
+ assert mesh == mesh_2
+
+# }}}
+
+
+# {{{ test lookup tree for element finding
+
+def test_lookup_tree(do_plot=False):
+ from meshmode.mesh.generation import make_curve_mesh, cloverleaf
+ mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, 1000), order=3)
+
+ from meshmode.mesh.tools import make_element_lookup_tree
+ tree = make_element_lookup_tree(mesh)
+
+ from meshmode.mesh.processing import find_bounding_box
+ bbox_min, bbox_max = find_bounding_box(mesh)
+
+ extent = bbox_max-bbox_min
+
+ for i in range(20):
+ pt = bbox_min + np.random.rand(2) * extent
+ print(pt)
+ for igrp, iel in tree.generate_matches(pt):
+ print(igrp, iel)
+
+ if do_plot:
+ with open("tree.dat", "w") as outf:
+ tree.visualize(outf)
+
+# }}}
+
+
+# {{{ test refiner
+
+def _get_vertex(mesh, vertex_index):
+ vertex = np.empty([len(mesh.vertices)])
+ for i_cur_dim, cur_dim_coords in enumerate(mesh.vertices):
+ vertex[i_cur_dim] = cur_dim_coords[vertex_index]
+ return vertex
+
+
+def get_blobby_2d_mesh():
+ from meshmode.mesh.io import generate_gmsh, FileSource
+ return generate_gmsh(
+ FileSource("blob-2d.step"), 2, order=2,
+ force_ambient_dim=2,
+ other_options=[
+ "-string", "Mesh.CharacteristicLengthMax = 1e-1;"]
+ )
+
+
+def get_some_mesh():
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron,
+ generate_torus, generate_regular_rect_mesh)
+ import random
+ #mesh = generate_icosphere(1, order=4)
+ #mesh = generate_icosahedron(1, order=4)
+ mesh = generate_torus(3, 1, order=4)
+
+ #mesh = generate_regular_rect_mesh()
+
+ from meshmode.mesh.refinement import Refiner
+ r = Refiner(mesh)
+ times = random.randint(1, 2)
+ for time in range(times):
+ flags = np.zeros(len(mesh.groups[0].vertex_indices))
+ '''
+ if time == 0:
+ flags[0] = 1
+ if time == 1:
+ #flags[1] = 1
+ flags[33] = 1
+ if time == 2:
+ flags[3] = 1
+ pass
+ #flags[8] = 1
+ #flags[40] = 1
+ '''
+ for i in range(0, len(flags)):
+ flags[i] = random.randint(0, 1)
+ mesh = r.refine(flags)
+
+ return mesh
+
+
+@pytest.mark.parametrize("mesh_factory", [get_some_mesh])
+@pytest.mark.parametrize("num_rounds", [1, 2, 3])
+def test_refiner_connectivity(mesh_factory, num_rounds):
+ mesh = mesh_factory()
+
+ def group_and_iel_to_global_iel(igrp, iel):
+ return mesh.groups[igrp].element_nr_base + iel
+
+ from meshmode.mesh.tools import make_element_lookup_tree
+ tree = make_element_lookup_tree(mesh)
+
+ from meshmode.mesh.processing import find_bounding_box
+ bbox_min, bbox_max = find_bounding_box(mesh)
+
+ cnx = mesh.element_connectivity
+ nvertices_per_element = len(mesh.groups[0].vertex_indices[0])
+ #stores connectivity obtained from geometry using barycentric coordinates
+ connected_to_element_geometry = np.zeros(
+ (mesh.nelements, mesh.nelements), dtype=bool)
+ #store connectivity obtained from mesh's connectivity
+ connected_to_element_connectivity = np.zeros(
+ (mesh.nelements, mesh.nelements), dtype=bool)
+
+ import six
+ for igrp, grp in enumerate(mesh.groups):
+ for iel_grp in six.moves.range(grp.nelements):
+
+ iel_g = group_and_iel_to_global_iel(igrp, iel_grp)
+ nb_starts = cnx.neighbors_starts
+ for nb_iel_g in cnx.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]:
+ connected_to_element_connectivity[nb_iel_g, iel_g] = True
+ connected_to_element_connectivity[iel_g, nb_iel_g] = True
+
+ for vertex_index in grp.vertex_indices[iel_grp]:
+ vertex = _get_vertex(mesh, vertex_index)
+
+ # check which elements touch this vertex
+ for bounding_igrp, bounding_iel in tree.generate_matches(vertex):
+ if bounding_igrp == igrp and bounding_iel == iel_grp:
+ continue
+ bounding_grp = mesh.groups[bounding_igrp]
+
+ last_bounding_vertex = _get_vertex(mesh,
+ bounding_grp.vertex_indices[bounding_iel][nvertices_per_element-1])
+
+ transformation = np.empty([len(mesh.vertices), nvertices_per_element-1])
+ vertex_transformed = vertex - last_bounding_vertex
+ for ibounding_vertex_index, bounding_vertex_index in enumerate(
+ bounding_grp.vertex_indices[bounding_iel][:nvertices_per_element-1]):
+ bounding_vertex = _get_vertex(mesh, bounding_vertex_index)
+ transformation[:, ibounding_vertex_index] = \
+ bounding_vertex - last_bounding_vertex
+ barycentric_coordinates, resid, _, _ = np.linalg.lstsq(transformation, vertex_transformed)
+
+ bc = barycentric_coordinates
+ tol = 1e-12
+
+ if resid > tol:
+ continue
+
+ if ((bc > -tol).all() and np.sum(bc) < 1+tol):
+ connected_to_element_geometry[
+ group_and_iel_to_global_iel(bounding_igrp, bounding_iel),
+ group_and_iel_to_global_iel(igrp, iel_grp)] = True
+ connected_to_element_geometry[
+ group_and_iel_to_global_iel(igrp, iel_grp),
+ group_and_iel_to_global_iel(bounding_igrp, bounding_iel)] = True
+
+ '''
+ print ("GEOMETRY: ")
+ print (connected_to_element_geometry)
+ print ("CONNECTIVITY: ")
+ print (connected_to_element_connectivity)
+# cmpmatrix = (connected_to_element_geometry == connected_to_element_connectivity)
+ #print cmpmatrix
+# print np.where(cmpmatrix == False)
+ if not (connected_to_element_geometry == connected_to_element_connectivity).all():
+ cmpmatrix = connected_to_element_connectivity == connected_to_element_geometry
+ for ii, i in enumerate(cmpmatrix):
+ for ij, j in enumerate(i):
+ if j == False:
+ pass
+ print (ii, ij)
+ '''
+ assert (connected_to_element_geometry == connected_to_element_connectivity).all()
+
+# }}}
+
+
+# {{{ test_nd_quad_submesh
+
+@pytest.mark.parametrize("dims", [2, 3, 4])
+def test_nd_quad_submesh(dims):
+ from meshmode.mesh.tools import nd_quad_submesh
+ from pytools import generate_nonnegative_integer_tuples_below as gnitb
+
+ node_tuples = list(gnitb(3, dims))
+
+ for i, nt in enumerate(node_tuples):
+ print(i, nt)
+
+ assert len(node_tuples) == 3**dims
+
+ elements = nd_quad_submesh(node_tuples)
+
+ for e in elements:
+ print(e)
+
+ assert len(elements) == 2**dims
+
+# }}}
+
+
if __name__ == "__main__":
import sys
if len(sys.argv) > 1: