Add Mesh.to_dict and generalize *_satisfying. Fixes #194, #158, #51.

docs/examples/ex02.py

@@ -16,18 +16,20 @@ def bilinf(u, du, ddu, v, dv, ddv, w):
 
     def C(T):
         trT = T[0,0] + T[1,1]
-        return np.array([[E/(1.0+nu)*(T[0, 0]+nu/(1.0-nu)*trT), E/(1.0+nu)*T[0, 1]],
-                         [E/(1.0+nu)*T[1, 0], E/(1.0+nu)*(T[1, 1]+nu/(1.0-nu)*trT)]])
+        return np.array([[E/(1.0+nu)*(T[0, 0]+nu/(1.0-nu)*trT),
+                          E/(1.0+nu)*T[0, 1]],
+                         [E/(1.0+nu)*T[1, 0],
+                          E/(1.0+nu)*(T[1, 1]+nu/(1.0-nu)*trT)]])
 
     def Eps(ddw):
         return np.array([[ddw[0][0], ddw[0][1]],
                          [ddw[1][0], ddw[1][1]]])
 
     def ddot(T1, T2):
-        return T1[0, 0]*T2[0, 0] +\
-               T1[0, 1]*T2[0, 1] +\
-               T1[1, 0]*T2[1, 0] +\
-               T1[1, 1]*T2[1, 1]
+        return (T1[0, 0]*T2[0, 0] +
+                T1[0, 1]*T2[0, 1] +
+                T1[1, 0]*T2[1, 0] +
+                T1[1, 1]*T2[1, 1])
 
     return d**3/12.0*ddot(C(Eps(ddu)), Eps(ddv))
 
@@ -39,9 +41,9 @@ def linf(v, dv, ddv, w):
 f = asm(linf, ib)
 
 boundary = {
-        'left':  m.facets_satisfying(lambda x, y: x==0),
-        'right': m.facets_satisfying(lambda x, y: x==1),
-        'top':   m.facets_satisfying(lambda x, y: y==1),
+        'left':  m.facets_satisfying(lambda x: x[0]==0),
+        'right': m.facets_satisfying(lambda x: x[0]==1),
+        'top':   m.facets_satisfying(lambda x: x[1]==1),
         }
 
 dofs = ib.get_dofs(boundary)

docs/examples/ex03.py

@@ -3,8 +3,8 @@
 from scipy.sparse.linalg import eigsh
 import numpy as np
 
-m1 = MeshLine(np.linspace(0,5,50))
-m2 = MeshLine(np.linspace(0,1,10))
+m1 = MeshLine(np.linspace(0, 5, 50))
+m2 = MeshLine(np.linspace(0, 1, 10))
 m = m1*m2
 
 e1 = ElementQuad1()
@@ -23,7 +23,7 @@ def mass(u, du, v, dv, w):
 
 M = asm(mass, gb)
 
-D = gb.get_dofs(lambda x, y: x==0.0).all()
+D = gb.get_dofs(lambda x: x[0]==0.0).all()
 y = np.zeros(gb.N)
 
 I = gb.complement_dofs(D)
@@ -33,6 +33,6 @@ def mass(u, du, v, dv, w):
 y[I] = x[:, 4]
 
 if __name__ == "__main__":
-    MeshQuad(np.array([m.p[0, :] + y[gb.nodal_dofs[0, :]],\
+    MeshQuad(np.array([m.p[0, :] + y[gb.nodal_dofs[0, :]],
                        m.p[1, :] + y[gb.nodal_dofs[1, :]]]), m.t).draw()
     m.show()

docs/examples/ex05.py

@@ -30,7 +30,7 @@ def facetlinf(v, dv, w):
 
 b = asm(facetlinf, fb)
 
-D = ib.get_dofs(lambda x, y: (y == 0.0)).all()
+D = ib.get_dofs(lambda x: (x[1] == 0.0)).all()
 I = ib.complement_dofs(D)
 
 import scipy.sparse

docs/examples/ex11.py

@@ -13,8 +13,8 @@
 K = asm(linear_elasticity(*lame_parameters(1e3, 0.3)), ib)
 
 dofs = {
-    'left' : ib.get_dofs(lambda x,y,z: x==0.0),
-    'right': ib.get_dofs(lambda x,y,z: x==1.0),
+    'left' : ib.get_dofs(lambda x: x[0]==0.0),
+    'right': ib.get_dofs(lambda x: x[0]==1.0),
     }
 
 u = np.zeros(K.shape[0])

docs/examples/ex25.py

@@ -23,8 +23,8 @@ def advection(u, du, v, dv, w):
     return v * velocity_0 * du[0]
 
 
-dofs = {'inlet': basis.get_dofs(lambda x, y: x == 0.),
-        'floor': basis.get_dofs(lambda x, y: y == 0.)}
+dofs = {'inlet': basis.get_dofs(lambda x: x[0] == 0.),
+        'floor': basis.get_dofs(lambda x: x[1] == 0.)}
 D = np.concatenate([d.all() for d in dofs.values()])
 interior = basis.complement_dofs(D)
 

skfem/mesh/mesh.py

@@ -62,20 +62,24 @@ class Mesh():
     p: ndarray = np.array([])
     t: ndarray = np.array([])
 
-    subdomains: Dict[str, ndarray] = {}
-    boundaries: Dict[str, ndarray] = {}
+    subdomains: Optional[Dict[str, ndarray]] = None
+    boundaries: Optional[Dict[str, ndarray]] = None
     external: Any = None
 
     def __init__(self):
         """Check that p and t are C_CONTIGUOUS as this leads
         to better performance."""
         if self.p is not None:
+            if not isinstance(self.p, ndarray):
+                self.p = np.array(self.p, dtype=np.float_)
             if self.p.flags['F_CONTIGUOUS']:
                 if self.p.shape[1] > 1000:
                     warnings.warn("Mesh.__init__(): Transforming "
                                   "over 100 vertices to C_CONTIGUOUS.")
                 self.p = np.ascontiguousarray(self.p)
         if self.t is not None:
+            if not isinstance(self.t, ndarray):
+                self.t = np.array(self.t, dtype=np.intp)
             if self.t.flags['F_CONTIGUOUS']:
                 if self.t.shape[1] > 1000:
                     warnings.warn("Mesh.__init__(): Transforming "
@@ -320,6 +324,64 @@ def element_finder(self) -> Callable[[ndarray], ndarray]:
         raise NotImplementedError("element_finder not implemented "
                                   "for the given Mesh type.")
 
+    def nodes_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
+        """Return nodes that satisfy some condition.
+
+        Parameters
+        ----------
+        test
+            A function which returns True for the set of nodes that are to be
+            included in the return set.
+
+        """
+        return np.nonzero(test(self.p))[0]
+
+    def facets_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
+        """Return facets whose midpoints satisfy some condition.
+
+        Parameters
+        ----------
+        test
+            A function which returns True for the facet midpoints that are to
+            be included in the return set.
+
+        """
+        midp = [np.sum(self.p[itr, self.facets], axis=0) / self.facets.shape[0]
+                for itr in range(self.p.shape[0])]
+        return np.nonzero(test(np.array(midp)))[0]
+
+    def elements_satisfying(self,
+                            test: Callable[[ndarray], bool]) -> ndarray:
+        """Return elements whose midpoints satisfy some condition.
+
+        Parameters
+        ----------
+        test
+            A function which returns True for the element midpoints that are to
+            be included in the return set.
+
+        """
+        midp = [np.sum(self.p[itr, self.t], axis=0) / self.t.shape[0]
+                for itr in range(self.p.shape[0])]
+        return np.nonzero(test(np.array(midp)))[0]
+
+    def to_dict(self) -> Dict[str, ndarray]:
+        """Return json serializable dictionary."""
+        if self.boundaries is not None:
+            boundaries = {k: v.tolist() for k, v in self.boundaries.items()}
+        else:
+            boundaries = self.boundaries
+        if self.subdomains is not None:
+            subdomains = {k: v.tolist() for k, v in self.subdomains.items()}
+        else:
+            subdomains = self.subdomains
+        return {
+            'p': self.p.tolist(),
+            't': self.t.tolist(),
+            'boundaries': boundaries,
+            'subdomains': subdomains,
+        }
+
     @staticmethod
     def strip_extra_coordinates(p: ndarray) -> ndarray:
         return p

skfem/mesh/mesh2d/mesh2d.py

@@ -23,47 +23,6 @@ class Mesh2D(Mesh):
     f2t: ndarray = np.array([])
     t2f: ndarray = np.array([])
 
-    def nodes_satisfying(self, test: Callable[[float, float], bool]) -> ndarray:
-        """Return nodes that satisfy some condition.
-
-        Parameters
-        ----------
-        test
-            A function which returns True for the set of nodes that are to be
-            included in the return set.
-
-        """
-        return np.nonzero(test(self.p[0, :], self.p[1, :]))[0]
-
-    def facets_satisfying(self, test: Callable[[float, float], bool]) -> ndarray:
-        """Return facets whose midpoints satisfy some condition.
-
-        Parameters
-        ----------
-        test
-            A function which returns True for the facet midpoints that are to
-            be included in the return set.
-
-        """
-        mx = 0.5*(self.p[0, self.facets[0, :]] + self.p[0, self.facets[1, :]])
-        my = 0.5*(self.p[1, self.facets[0, :]] + self.p[1, self.facets[1, :]])
-        return np.nonzero(test(mx, my))[0]
-
-    def elements_satisfying(self,
-                            test: Callable[[float, float], bool]) -> ndarray:
-        """Return elements whose midpoints satisfy some condition.
-
-        Parameters
-        ----------
-        test
-            A function which returns True for the element midpoints that are to
-            be included in the return set.
-
-        """
-        mx = np.sum(self.p[0, self.t], axis=0)/self.t.shape[0]
-        my = np.sum(self.p[1, self.t], axis=0)/self.t.shape[0]
-        return np.nonzero(test(mx, my))[0]
-
     def draw(self,
              ax: Optional[Axes] = None,
              node_numbering: Optional[bool] = False,

skfem/mesh/mesh2d/mesh_quad.py

@@ -75,10 +75,10 @@ def __init__(self,
         self.t = t
         self.boundaries = boundaries
         self.subdomains = subdomains
+        super(MeshQuad, self).__init__()
         if validate:
             self._validate()
         self._build_mappings()
-        super(MeshQuad, self).__init__()
 
     @classmethod
     def init_tensor(cls: Type[MeshType],

skfem/mesh/mesh2d/mesh_tri.py

@@ -121,10 +121,10 @@ def __init__(self,
         self.t = t
         self.boundaries = boundaries
         self.subdomains = subdomains
+        super(MeshTri, self).__init__()
         if validate:
             self._validate()
         self._build_mappings(sort_t=sort_t)
-        super(MeshTri, self).__init__()
 
     @classmethod
     def init_tensor(cls: Type[MeshType],

skfem/mesh/mesh3d/mesh3d.py

@@ -7,9 +7,6 @@
 from numpy import ndarray
 
 
-BoolFun = Callable[[float, float, float], bool]
-
-
 class Mesh3D(Mesh):
     """Three dimensional meshes, common methods.
 
@@ -20,33 +17,7 @@ class Mesh3D(Mesh):
 
     """
 
-    def nodes_satisfying(self, test: BoolFun) -> ndarray:
-        """Return nodes that satisfy some condition.
-
-        Parameters
-        ----------
-        test
-            Evaluates to 1 or True for nodes belonging to the output set.
-
-        """
-        return np.nonzero(test(self.p[0, :], self.p[1, :], self.p[2, :]))[0]
-
-    def facets_satisfying(self, test: BoolFun) -> ndarray:
-        """Return facets whose midpoints satisfy some condition.
-        
-        Parameters
-        ----------
-        test
-            Evaluates to 1 or True for facet midpoints of the facets belonging
-            to the output set.
-
-        """
-        mx = np.sum(self.p[0, self.facets], axis=0) / self.facets.shape[0]
-        my = np.sum(self.p[1, self.facets], axis=0) / self.facets.shape[0]
-        mz = np.sum(self.p[2, self.facets], axis=0) / self.facets.shape[0]
-        return np.nonzero(test(mx, my, mz))[0]
-
-    def edges_satisfying(self, test: BoolFun) -> ndarray:
+    def edges_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
         """Return edges whose midpoints satisfy some condition.
 
         Parameters
@@ -61,21 +32,6 @@ def edges_satisfying(self, test: BoolFun) -> ndarray:
         mz = 0.5 * (self.p[2, self.edges[0, :]] + self.p[2, self.edges[1, :]])
         return np.nonzero(test(mx, my, mz))[0]
 
-    def elements_satisfying(self, test: BoolFun) -> ndarray:
-        """Return elements whose midpoints satisfy some condition.
-
-        Parameters
-        ----------
-        test
-            Evaluates to 1 or True for element midpoints of the elements
-            belonging to the output set.
-
-        """
-        mx = np.sum(self.p[0, self.t], axis=0) / self.t.shape[0]
-        my = np.sum(self.p[1, self.t], axis=0) / self.t.shape[0]
-        mz = np.sum(self.p[2, self.t], axis=0) / self.t.shape[0]
-        return np.nonzero(test(mx, my, mz))[0]
-
     def boundary_edges(self) -> ndarray:
         """Return an array of boundary edge indices."""
         return np.nonzero(np.isin(self.edges,

skfem/mesh/mesh3d/mesh_hex.py

@@ -77,10 +77,10 @@ def __init__(self,
         self.t = t
         self.boundaries = boundaries
         self.subdomains = subdomains
+        super(MeshHex, self).__init__()
         if validate:
             self._validate()
         self._build_mappings()
-        super(MeshHex, self).__init__()
 
     @classmethod
     def init_tensor(cls: Type[MeshType],

skfem/mesh/mesh3d/mesh_tet.py

@@ -70,11 +70,11 @@ def __init__(self,
         self.t = t
         self.boundaries = boundaries
         self.subdomains = subdomains
+        super(MeshTet, self).__init__()
         if validate:
             self._validate()
         self.enable_facets = True
         self._build_mappings()
-        super(MeshTet, self).__init__()
 
     @classmethod
     def init_tensor(cls: Type[MeshType],

skfem/mesh/mesh_line.py

@@ -36,11 +36,11 @@ def __init__(self,
         self.facets = np.arange(self.p.shape[1])[None, :]
         self.t = np.vstack([self.facets[0, :-1],
                             self.facets[0, 1:]]) if t is None else t
+        super(MeshLine, self).__init__()
         self._build_mappings()
 
         if validate:
             self._validate()
-        super(MeshLine, self).__init__()
 
     @classmethod
     def init_refdom(cls: Type[MeshType]) -> MeshType:

tests/test_assembly.py

@@ -70,7 +70,7 @@ def runTest(self):
         mtype, etype = self.case
         m = mtype()
         m.refine(3)
-        bnd = m.facets_satisfying(lambda x,y,z: x==1.0)
+        bnd = m.facets_satisfying(lambda x: x[0]==1.0)
         fb = FacetBasis(m, etype(), facets=bnd)
 
         @bilinear_form
@@ -109,7 +109,8 @@ def runTest(self):
         x = self.initOnes(ib)
         f = ib.interpolator(x)
 
-        self.assertTrue(np.sum(f(np.array([np.sin(m.p[0, :]), np.sin(3.0*m.p[1, :])]))-1.0) < 1e-10)
+        self.assertTrue(np.sum(f(np.array([np.sin(m.p[0, :]),
+                                           np.sin(3.0*m.p[1, :])]))-1.0) < 1e-10)
 
 class BasisInterpolatorMorley(BasisInterpolator):
     case = (MeshTri, ElementTriMorley)