add poisson 2d

README.md

@@ -13,8 +13,15 @@ pip install git+https://github.com/tataratat/feigen.git@main
 ```
 
 ## Quick start
+### iganet
 Current version supports iganet's BSplineSurface.
 Assuming that you have a server running,
 ```
 python3 -c "import feigen; feigen.BSpline2D('ws://localhost:9001').start()"
 ```
+
+### IGA examples
+#### Poisson problem 2D
+```
+python3 -c "import feigen; feigen.Poisson2D().start()"
+```

examples/poisson2d.py

@@ -0,0 +1,7 @@
+import feigen
+
+if __name__ == "__main__":
+    # feigen.log.configure(debug=True)
+
+    b = feigen.Poisson2D()
+    b.start()

feigen/__init__.py

@@ -1,11 +1,14 @@
-from feigen import bspline, comm, log
+from feigen import bspline, comm, log, poisson2d
 from feigen._version import __version__
 from feigen.bspline import BSpline2D
+from feigen.poisson2d import Poisson2D
 
 __all__ = [
     "__version__",
     "bspline",
     "comm",
     "log",
+    "poisson2d",
     "BSpline2D",
+    "Poisson2D",
 ]

feigen/poisson2d.py

@@ -1,10 +1,8 @@
 import numpy as np
 import splinepy
 import vedo
-from gustaf.utils.arr import enforce_len
 from scipy.sparse import csr_array, linalg
 
-from feigen import comm
 from feigen._base import FeigenBase
 
 # skip bound checks
@@ -63,8 +61,8 @@ def setup(self, boundary_splines):
         self.end_ids = (
             int(len(boundary_splines[0].cps) - 1),
             int(len(boundary_splines[0].cps) - 1),
-            int(len(boundary_splines[1].cps) - 1),
-            int(len(boundary_splines[1].cps) - 1),
+            int(len(boundary_splines[2].cps) - 1),
+            int(len(boundary_splines[2].cps) - 1),
         )
 
         # this is a bit verbose and can be also done with mod but
@@ -162,7 +160,7 @@ def _process_spline_actors(plt):
         return None
 
     plt._state["spline_cp_actors"][cp_id].SetPosition(
-        plt._state["spline"].cps[cp_id]
+        np.append(plt._state["spline"].cps[cp_id], [0])
     )
 
 
@@ -262,14 +260,39 @@ def _process_parametric_view(plt):
     plt._state["parametric_view_actors"] = actors
 
 
+def _bid_to_dof(spl, bid):
+    """ """
+    ortho_dim, extrema = divmod(bid, spl.para_dim)
+
+    all_ = slice(None)
+    indices = [all_] * spl.para_dim
+    indices[ortho_dim] = int(extrema * -1)
+
+    dir_flip = -1.0 if extrema else 1.0
+
+    return spl.multi_index[tuple(indices)], ortho_dim, dir_flip
+
+
+def refine(spl):
+    """
+    A predefined refinement rule, so that you can apply to both field and
+    boundaries.
+    """
+    spl.elevate_degrees(list(range(spl.para_dim)) * 2)
+    n_new_kvs = 10
+    for i, (l_bound, u_bound) in enumerate(zip(*spl.parametric_bounds)):
+        spl.insert_knots(
+            i,
+            np.linspace(l_bound, u_bound, int(n_new_kvs + 2))[1:-1],
+        )
+    return spl
+
+
 class Poisson2D(vedo.Plotter, FeigenBase):
     """
     Self contained interactive plotter base.
     Can move control points on the left plot and the right plot will show
-    the response from iganet.
-
-    Upon initialization, it will create bspline at iganet server and
-    setup interactive plotter.
+    the solution of simple laplace problem
 
     Parameters
     ----------
@@ -327,8 +350,21 @@ def __init__(self, spline=None):  # noqa PLR0915
         self.add_callback("LeftButtonRelease", self._left_release)
         self.add_callback("RightButtonPress", self._right_click)
 
-        # create matching spline
-        self._state["spline"] = splinepy.helpme.create.surface_circle(1)
+        # register spline
+        if spline is None:
+            self._state["spline"] = splinepy.helpme.create.surface_circle(
+                1
+            ).nurbs
+            self._state["spline"].insert_knots(0, [0.5])
+            self._state["spline"].insert_knots(1, [0.5])
+        else:
+            support_dim = 2
+            if spline.para_dim != support_dim or spline.dim != support_dim:
+                raise ValueError(
+                    "This plotter is built for splines with "
+                    "para_dim=2 and dim=2"
+                )
+            self._state["spline"] = spline
 
         # create boundary splines - will be used to control BCs
         conform_param_view = self._state["spline"].create.parametric_view(
@@ -581,18 +617,18 @@ def _right_click(self, evt):  # noqa ARG002
         solution = self._state.get("solution_spline", None)
         if solution is None:
             dict_spline = geometry.current_core_properties()
-            dict_spline["control_points"] = np.empty((len(geometry.cps), 1), dtype="float64")
+            dict_spline["control_points"] = np.zeros(
+                (len(geometry.cps), 1), dtype="float64"
+            )
             solution = type(geometry)(**dict_spline)
             self._state["solution_spline"] = solution
 
-            # we refine the solution here, we don't have to
-            solution.elevate_degrees([0,0,1,1])
-            n_new_kvs = 5
-            for i, (l_bound, u_bound) in enumerate(solution.parametric_bounds):
-                solution.insert_knots(i, np.linspace(l_bound, u_bound, int(n_new_kvs + 2))[1:-1])
+            refine(solution)
 
-            self._logd(f"created solution spline - degrees: {solution.ds}"
-                "knot_vectors: {solution.kvs}")
+            self._logd(
+                f"created solution spline - degrees: {solution.ds}"
+                f"knot_vectors: {solution.kvs}"
+            )
 
             # get greville points and sparsity pattern
             queries = solution.greville_abscissae
@@ -606,62 +642,68 @@ def _right_click(self, evt):  # noqa ARG002
             self._state["solution_row_ids"] = row_ids
             self._state["solution_n_col"] = n_col
             self._state["solution_col_ids"] = col_ids
-
 
         queries = self._state["solution_queries"]
 
         # get mapper and evaluate basis laplacian
         mapper = solution.mapper(reference=geometry)
         b_laplacian, _ = mapper.basis_laplacian_and_support(queries)
 
-        # evaluate rhs
-        # we could add source function here
-        rhs = np.ones(len(queries))
-
         # prepare lhs sparse matrix
         sp_laplacian = csr_array(
-            (b_laplacian.ravel(), (self._state["solution_row_ids"], self._state["solution_col_ids"])), shape=[self._state["solution_n_row"]] * 2
+            (
+                b_laplacian.ravel(),
+                (
+                    self._state["solution_row_ids"],
+                    self._state["solution_col_ids"],
+                ),
+            ),
+            shape=[self._state["solution_n_row"]] * 2,
         )
 
+        # prepare rhs
+        # we could add source function here
+        mi = solution.multi_index
+        boundary_dofs = np.concatenate((mi[[0, -1], :], mi[1:-1, [0, -1]]))
+        rhs = np.ones(len(queries))
+        rhs[boundary_dofs] = 0.0
+
+        # set zeros and ones on lhs
+        sp_laplacian[boundary_dofs] *= 0
+        sp_laplacian[boundary_dofs, boundary_dofs] = 1
+
+        # apply boundary condition
+        # we will loop and extract boundary values
+        for bid, (bdr_spline, bdr_area) in enumerate(
+            zip(
+                self._state["boundary_splines"],
+                self._state["boundary_spline_areas"],
+            )
+        ):
+            # now, based on bdr_area, we extract values for boundary condition
+            cp_half = int(len(bdr_area.cps) / 2)
+
+            bdr = bdr_spline.copy()
+            bdr.cps[:] = bdr_area.cps[:cp_half] - bdr_area.cps[cp_half:]
+            # self._state["solution_refinement"](bdr)
+            refine(bdr)
+            b_dof, axis, factor = _bid_to_dof(solution, bid)
+            rhs[b_dof] = bdr.cps[:, axis] * factor
+
+        solution.cps[:] = linalg.spsolve(-sp_laplacian, rhs).reshape(-1, 1)
+
         # eval field - currently, th
-        field = self._config["form"].evaluate_model(
-            self._config["spline_id"],
-            "ValueFieldMagnitude",
-            enforce_len(
-                self._config["sample_resolutions"],
-                self._state["spline"].para_dim,
-            ).tolist(),
-        )  # this is ravel
+        field = solution.sample(self._config["sample_resolutions"])
 
         # update vedo object
         # we should copy the geometry (spline),
         # because we don't want two plots to have same color
         geo_actor = self._state["spline_actors"]["spline"].clone()
-        geo_actor.cmap("plasma", field).alpha(0.9)
+        geo_actor.cmap("plasma", field).alpha(0.9).add_scalarbar3d()
         self._state["server_plot_actors"] = {
             "spline": geo_actor,
             "knots": self._state["spline_actors"]["knots"],
         }
-        self.add(
-            *self._state["server_plot_actors"].values(),
-            at=self._config["server_plot"],
-        )
-
-        # eval current geometry
-        geo_eval = self._config["form"].evaluate_model(
-            self._config["spline_id"],
-            "ValueField",
-            enforce_len(
-                int(self._config["sample_resolutions"] / 5),
-                self._state["spline"].para_dim,
-            ).tolist(),
-        )
-        eval_points = vedo.Points(geo_eval, c="white")
-        eval_point_ids = eval_points.labels("id", on="points", font="VTK")
-        self._state["server_plot_actors"]["evaluated_points"] = eval_points
-        self._state["server_plot_actors"][
-            "evaluated_point_ids"
-        ] = eval_point_ids.c("grey")
 
         # for some reason add() won't work here
         self.show(
@@ -712,7 +754,7 @@ def start(self):
         )
 
         self.show(
-            "IgaNet server - right click to sync",
+            "Solution - right click to sync",
             at=self._config["server_plot"],
             interactive=False,
             mode=self._config["plotter_mode"],
@@ -721,7 +763,4 @@ def start(self):
         # let's start
         self.interactive()
 
-        # TODO - close websocket here?
-        self._config["iganet_ws"].websocket.close()
-
         return self