From 0a6f2fba45d667ee889960abe2fa7624fdb11d6f Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sat, 21 Jun 2014 18:37:19 -0500
Subject: [PATCH] Make boundary restriction work, add unit test
---
meshmode/discretization/connection.py | 103 ++++++++++++++++++++----
meshmode/discretization/poly_element.py | 16 +++-
meshmode/mesh/__init__.py | 4 +-
test/blob-2d.step | 95 ++++++++++++++++++++++
test/test_meshmode.py | 91 +++++++++++++++++++++
5 files changed, 287 insertions(+), 22 deletions(-)
create mode 100644 test/blob-2d.step
create mode 100644 test/test_meshmode.py
diff --git a/meshmode/discretization/connection.py b/meshmode/discretization/connection.py
index bbdfd21..974a18c 100644
--- a/meshmode/discretization/connection.py
+++ b/meshmode/discretization/connection.py
@@ -26,7 +26,7 @@ import numpy as np
import modepy as mp
import pyopencl as cl
import pyopencl.array # noqa
-from pytools import memoize_method, memoize_method_nested
+from pytools import memoize_method, memoize_method_nested, Record
__doc__ = """
@@ -34,6 +34,8 @@ __doc__ = """
.. autofunction:: make_same_mesh_connection
+.. autofunction:: make_boundary_restriction
+
Implementation details
^^^^^^^^^^^^^^^^^^^^^^
@@ -54,13 +56,13 @@ class InterpolationBatch(object):
.. attribute:: source_element_indices
- A :class:`numpy.ndarray` of length ``nelements``, containing the
+ 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:`numpy.ndarray` of length ``nelements``, containing the
+ A :class:`pyopencl.array.Arrays` of length ``nelements``, containing the
element index to which this "*to*" element's data will be
interpolated.
@@ -72,6 +74,7 @@ class InterpolationBatch(object):
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
@@ -136,7 +139,18 @@ class DiscretizationConnection(object):
0<=k<nelements and
0<=i<n_to_nodes and
0<=j<n_from_nodes}""",
- "result[k,i] = sum(j, resample_mat[i, j] * vec[k, j])",
+ "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")
@@ -153,9 +167,12 @@ class DiscretizationConnection(object):
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):
- knl()(queue,
- resample_mat=self._resample_matrix(i_grp, i_batch),
- result=sgrp.view(result), vec=tgrp.view(vec))
+ 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
@@ -193,7 +210,47 @@ def make_same_mesh_connection(queue, to_discr, from_discr):
# {{{ boundary restriction constructor
-def make_boundary_extractor(queue, discr, group_factory):
+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 xrange(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)``
"""
@@ -205,6 +262,7 @@ def make_boundary_extractor(queue, discr, group_factory):
for igrp, mgrp in enumerate(discr.mesh.groups):
grp_face_vertex_indices = mgrp.face_vertex_indices()
+
for iel_grp in xrange(mgrp.nelements):
for fid, loc_face_vertices in enumerate(grp_face_vertex_indices):
face_vertices = frozenset(
@@ -237,6 +295,8 @@ def make_boundary_extractor(queue, discr, group_factory):
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 = [
@@ -274,7 +334,7 @@ def make_boundary_extractor(queue, discr, group_factory):
batch_base = 0
for face_id in xrange(len(grp_face_vertex_indices)):
- batch_boundary_el_numbers_in_vol = np.array(
+ batch_boundary_el_numbers_in_grp = np.array(
[
ibface_el
for ibface_el, ibface_face in group_boundary_faces
@@ -283,7 +343,7 @@ def make_boundary_extractor(queue, discr, group_factory):
new_el_numbers = np.arange(
batch_base,
- batch_base + len(batch_boundary_el_numbers_in_vol))
+ batch_base + len(batch_boundary_el_numbers_in_grp))
# {{{ no per-element axes in these computations
@@ -309,9 +369,11 @@ def make_boundary_extractor(queue, discr, group_factory):
# }}}
+ # {{{ build information for mesh element group
+
# Find vertex_indices
glob_face_vertices = mgrp.vertex_indices[
- batch_boundary_el_numbers_in_vol][:, loc_face_vertices]
+ batch_boundary_el_numbers_in_grp][:, loc_face_vertices]
vertex_indices[new_el_numbers] = \
vol_to_bdry_vertices[glob_face_vertices]
@@ -319,9 +381,18 @@ def make_boundary_extractor(queue, discr, group_factory):
nodes[:, new_el_numbers, :] = np.einsum(
"ij,dej->dei",
resampling_mat,
- mgrp.nodes[:, batch_boundary_el_numbers_in_vol, :])
+ 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_vol)
+ batch_base += len(batch_boundary_el_numbers_in_grp)
bdry_mesh_group = SimplexElementGroup(
mgrp.order, vertex_indices, nodes, unit_nodes=bdry_unit_nodes)
@@ -333,10 +404,8 @@ def make_boundary_extractor(queue, discr, group_factory):
bdry_discr = Discretization(
discr.cl_context, bdry_mesh, group_factory)
- # FIXME
- connection = None
-
- return bdry_mesh, bdry_discr, connection
+ return bdry_mesh, bdry_discr, _build_boundary_connection(
+ queue, discr, bdry_discr, connection_data)
# }}}
diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py
index 43445c6..c15a1b1 100644
--- a/meshmode/discretization/poly_element.py
+++ b/meshmode/discretization/poly_element.py
@@ -96,7 +96,13 @@ class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
@property
@memoize_method
def unit_nodes(self):
- return self._quadrature_rule().nodes
+ result = self._quadrature_rule().nodes
+ if len(result.shape) == 1:
+ result = np.array([result])
+
+ dim2, nunit_nodes = result.shape
+ assert dim2 == self.mesh_el_group.dim
+ return result
@property
@memoize_method
@@ -108,8 +114,12 @@ class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase):
@property
@memoize_method
def unit_nodes(self):
- dims = self.mesh_el_group.dim
- return mp.warp_and_blend_nodes(dims, self.order)
+ dim = self.mesh_el_group.dim
+ result = mp.warp_and_blend_nodes(dim, self.order)
+
+ dim2, nunit_nodes = result.shape
+ assert dim2 == dim
+ return result
@property
@memoize_method
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 91f3a7b..ef2c104 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -80,7 +80,7 @@ class MeshElementGroup(Record):
vertex_indices=vertex_indices,
nodes=nodes,
unit_nodes=unit_nodes,
- element_nr_base=None, node_nr_base=None)
+ element_nr_base=element_nr_base, node_nr_base=node_nr_base)
@property
def dim(self):
@@ -217,7 +217,7 @@ class Mesh(Record):
el_nr += ng.nelements
node_nr += ng.nnodes
- Record.__init__(self, vertices=vertices, groups=groups)
+ Record.__init__(self, vertices=vertices, groups=new_groups)
@property
def ambient_dim(self):
diff --git a/test/blob-2d.step b/test/blob-2d.step
new file mode 100644
index 0000000..f0b8bf0
--- /dev/null
+++ b/test/blob-2d.step
@@ -0,0 +1,95 @@
+ISO-10303-21;
+HEADER;
+FILE_DESCRIPTION(('FreeCAD Model'),'2;1');
+FILE_NAME('/home/andreas/own/graphics/freecad/blob.step',
+ '2014-06-16T23:10:06',('FreeCAD'),('FreeCAD'),
+ 'Open CASCADE STEP processor 6.5','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('BSpline','BSpline','',(#8));
+#8 = MECHANICAL_CONTEXT('',#2,'mechanical');
+#9 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design');
+#10 = MANIFOLD_SURFACE_SHAPE_REPRESENTATION('',(#11,#15),#53);
+#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),#36,.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('',(-4.72543090582E-04,1.64651736617E-04,0.));
+#24 = SURFACE_CURVE('',#25,(#35),.PCURVE_S1.);
+#25 = B_SPLINE_CURVE_WITH_KNOTS('',3,(#26,#27,#28,#29,#30,#31,#32,#33,
+ #34),.UNSPECIFIED.,.T.,.F.,(1,1,2,1,1,1,1,1,2,1,1),(-0.574540704647,
+ -0.371453182518,0.,0.263932510809,0.520877807931,0.808151393473,
+ 1.190527579206,1.393615101335,1.765068283853,2.029000794662,
+ 2.285946091784),.UNSPECIFIED.);
+#26 = CARTESIAN_POINT('',(-4.451395471938E-04,3.893567516391E-04,0.));
+#27 = CARTESIAN_POINT('',(-4.920144165846E-04,4.989734519847E-06,0.));
+#28 = CARTESIAN_POINT('',(-4.111934608223E-04,-1.366538732583E-04,0.));
+#29 = CARTESIAN_POINT('',(-1.187450244609E-04,-2.333067459084E-05,0.));
+#30 = CARTESIAN_POINT('',(2.486504062773E-04,-3.199830111205E-04,0.));
+#31 = CARTESIAN_POINT('',(1.570695404062E-04,2.640870359319E-04,0.));
+#32 = CARTESIAN_POINT('',(-1.562373623195E-04,2.334791360325E-04,0.));
+#33 = CARTESIAN_POINT('',(-4.451395471938E-04,3.893567516391E-04,0.));
+#34 = CARTESIAN_POINT('',(-4.920144165846E-04,4.989734519847E-06,0.));
+#35 = PCURVE('',#36,#41);
+#36 = PLANE('',#37);
+#37 = AXIS2_PLACEMENT_3D('',#38,#39,#40);
+#38 = CARTESIAN_POINT('',(-4.451395471938E-04,3.893567516391E-04,0.));
+#39 = DIRECTION('',(0.,0.,-1.));
+#40 = DIRECTION('',(-1.,0.,0.));
+#41 = DEFINITIONAL_REPRESENTATION('',(#42),#52);
+#42 = B_SPLINE_CURVE_WITH_KNOTS('',3,(#43,#44,#45,#46,#47,#48,#49,#50,
+ #51),.UNSPECIFIED.,.F.,.F.,(1,1,2,1,1,1,1,1,2,1,1),(-0.574540704647,
+ -0.371453182518,0.,0.263932510809,0.520877807931,0.808151393473,
+ 1.190527579206,1.393615101335,1.765068283853,2.029000794662,
+ 2.285946091784),.UNSPECIFIED.);
+#43 = CARTESIAN_POINT('',(0.,0.));
+#44 = CARTESIAN_POINT('',(4.687486939076E-05,-3.843670171192E-04));
+#45 = CARTESIAN_POINT('',(-3.394608637149E-05,-5.260106248974E-04));
+#46 = CARTESIAN_POINT('',(-3.263945227329E-04,-4.126874262299E-04));
+#47 = CARTESIAN_POINT('',(-6.937899534711E-04,-7.093397627596E-04));
+#48 = CARTESIAN_POINT('',(-6.022090876E-04,-1.252697157072E-04));
+#49 = CARTESIAN_POINT('',(-2.889021848743E-04,-1.558776156066E-04));
+#50 = CARTESIAN_POINT('',(0.,0.));
+#51 = CARTESIAN_POINT('',(4.687486939076E-05,-3.843670171192E-04));
+#52 = ( GEOMETRIC_REPRESENTATION_CONTEXT(2)
+PARAMETRIC_REPRESENTATION_CONTEXT() REPRESENTATION_CONTEXT('2D SPACE',''
+ ) );
+#53 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3)
+GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#57)) GLOBAL_UNIT_ASSIGNED_CONTEXT(
+(#54,#55,#56)) REPRESENTATION_CONTEXT('Context #1',
+ '3D Context with UNIT and UNCERTAINTY') );
+#54 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT($,.METRE.) );
+#55 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) );
+#56 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() );
+#57 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-10),#54,
+ 'distance_accuracy_value','confusion accuracy');
+#58 = PRODUCT_TYPE('part',$,(#7));
+#59 = MECHANICAL_DESIGN_GEOMETRIC_PRESENTATION_REPRESENTATION('',(#60),
+ #53);
+#60 = STYLED_ITEM('color',(#61),#17);
+#61 = PRESENTATION_STYLE_ASSIGNMENT((#62));
+#62 = SURFACE_STYLE_USAGE(.BOTH.,#63);
+#63 = SURFACE_SIDE_STYLE('',(#64));
+#64 = SURFACE_STYLE_FILL_AREA(#65);
+#65 = FILL_AREA_STYLE('',(#66));
+#66 = FILL_AREA_STYLE_COLOUR('',#67);
+#67 = 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
new file mode 100644
index 0000000..a1e0d4a
--- /dev/null
+++ b/test/test_meshmode.py
@@ -0,0 +1,91 @@
+from __future__ import division
+
+__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 numpy.linalg as la
+import pyopencl as cl
+import pyopencl.array # noqa
+import pyopencl.clmath # noqa
+
+from pyopencl.tools import ( # noqa
+ pytest_generate_tests_for_pyopencl
+ as pytest_generate_tests)
+
+
+def test_boundary_interpolation(ctx_getter):
+ 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 \
+ QuadratureSimplexGroupFactory
+ from pytools.convergence import EOCRecorder
+ from meshmode.discretization.connection import make_boundary_restriction
+
+ 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_dimension=2,
+ other_options=[
+ "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
+ )
+ print "END GEN"
+
+ vol_discr = Discretization(cl_ctx, mesh,
+ QuadratureSimplexGroupFactory(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)
+
+ bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
+ queue, vol_discr, QuadratureSimplexGroupFactory(order))
+
+ 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)
+
+ 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
+
+
+if __name__ == "__main__":
+ import sys
+ if len(sys.argv) > 1:
+ exec(sys.argv[1])
+ else:
+ from py.test.cmdline import main
+ main([__file__])
+
+# vim: fdm=marker
--
GitLab