Make compatible with DOLFINx v0.10.0

.gitignore

@@ -80,3 +80,4 @@ docs/*/_autosummary
 # FEniCS failures
 jitfailure*
 
+_build

README.md

@@ -8,11 +8,11 @@
 
 ## Prerequisites
 **dolfinx_optim** requires: 
-* **FEniCSx** (v.0.8), see [installation instructions here](https://fenicsproject.org/download/).
+* **FEniCSx** (v.0.10), see [installation instructions here](https://fenicsproject.org/download/).
 * **MOSEK** (>= version 10 with its Python Fusion interface), see [installation instructions here](https://www.mosek.com/downloads/). The Python interface can be simply installed via `pip`:
 
 ```
-pip install -f https://download.mosek.com/stable/wheel/index.html Mosek
+pip install -f mosek
 ```
 
 Mosek is a commercial software so users need a valid Mosek license. Free unlimited licenses are available for education and research purposes, see the [Academic License section](https://www.mosek.com/products/academic-licenses/).

docs/_config.yml → _config.yml

@@ -20,7 +20,7 @@ latex:
 
 # Add a bibtex file so that we can create citations
 bibtex_bibfiles:
-  - references.bib
+  - ./docs/references.bib
 
 # Information about where the book exists on the web
 # repository:
@@ -32,7 +32,8 @@ exclude_patterns: [
     ...,
     docs/index.md,
     README.md,
-    "**/*.ipynb"]
+    "**/*.ipynb",
+    "*.pytest_cache"]
 
 # Sphinx configuration for custom theme
 sphinx:
@@ -58,7 +59,12 @@ sphinx:
   - 'sphinx.ext.viewcode'
   - 'sphinx.ext.autosummary'
   - 'docs.remove_docstring' # remove module docstring
-
+  - 'sphinx.ext.intersphinx'
+  - "sphinx_codeautolink"
+
+
+
+
 parse:
   myst_enable_extensions:
     - "amsmath"

docs/_toc.yml → _toc.yml

@@ -2,7 +2,7 @@
 # Learn more at https://jupyterbook.org/customize/toc.html
 
 format: jb-book
-root: index
+root: index.md
 title: dolfinx_optim Documentation
 
 
@@ -15,4 +15,4 @@ parts:
       - file: demos/cartoon_texture_decomposition/cartoon_texture_decomposition.md
   - caption: API
     chapters:
-      - file: api/api.rst
+      - file: docs/api/api.rst

demos/3D_soil_slope_limit_analysis/3D_soil_slope_limit_analysis.md

@@ -113,7 +113,7 @@ def border(x):
 
 
 gdim = 3
-V = fem.functionspace(domain, ("CG", 2, (gdim,)))
+V = fem.functionspace(domain, ("Lagrange", 2, (gdim,)))
 bc_dofs = fem.locate_dofs_geometrical(V, border)
 bcs = [fem.dirichletbc(np.zeros((gdim,)), bc_dofs, V)]
 ```

demos/cartoon_texture_decomposition/cartoon_texture_decomposition.md

@@ -83,7 +83,7 @@ colliding_cells = geometry.compute_colliding_cells(
 # represent image as DG0 on quad mesh
 V00 = fem.functionspace(domain0, ("DG", 0))
 y0 = fem.Function(V00)
-cells = [c for i in range(len(colliding_cells)) for c in colliding_cells.links(i)]
+cells = colliding_cells.array
 y0.x.array[cells] = image.ravel()
 ```
 
@@ -102,11 +102,11 @@ y = fem.Function(V0)
 fine_mesh_cell_map = domain.topology.index_map(domain.topology.dim)
 num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
 cells = np.arange(num_cells_on_proc, dtype=np.int32)
-interpolation_data = fem.create_nonmatching_meshes_interpolation_data(
-    V0.mesh.geometry, V0.element, V00.mesh, cells, padding=1e-14
+interpolation_data = fem.create_interpolation_data(
+    V0, V00, cells, padding=1e-14
 )
 # interpolate on non-matching mesh
-y.interpolate(y0, nmm_interpolation_data=interpolation_data)
+y.interpolate_nonmatching(y0, cells=cells, interpolation_data=interpolation_data)
 ```
 
 We now define the variational problem by creating the two optimization variables $u$ and $\boldsymbol{g}$.

demos/elastoplasticity/elastoplasticity_exponential_hardening.md

@@ -137,7 +137,7 @@ We create a rectangular plate with two circular notches on its sides using `gmsh
 
 def generate_notched_plate(W, H, R, mesh_size):
     import gmsh
-    from dolfinx.io.gmshio import model_to_mesh
+    from dolfinx.io.gmsh import model_to_mesh
 
     gmsh.initialize()
     gmsh.option.setNumber("General.Terminal", 0)  # to disable meshing info
@@ -189,15 +189,15 @@ def generate_notched_plate(W, H, R, mesh_size):
 
         gmsh.model.mesh.generate(gdim)
 
-        domain, markers, facets = model_to_mesh(
+        mesh_data = model_to_mesh(
             gmsh.model,
             mesh_comm,
             model_rank,
             gdim=gdim,
         )
 
     gmsh.finalize()
-    return domain, markers, facets
+    return mesh_data.mesh, mesh_data.cell_tags, mesh_data.facet_tags
 ```
 
 We generate the mesh and define the different physical constants for the problem.
@@ -368,7 +368,7 @@ We now define the relevant function spaces associated with the discretization of
 
 ```{code-cell} ipython3
 # P2 interpolation for velocity
-V = fem.functionspace(domain, ("CG", 2, (2,)))
+V = fem.functionspace(domain, ("Lagrange", 2, (2,)))
 Vuy, _ = V.sub(1).collapse()
 Vepsp = fem.functionspace(domain, ("DG", 1, (4,)))
 Vp,_ = Vepsp.sub(0).collapse()
@@ -432,11 +432,11 @@ for i, t in enumerate(t_list[1:]):
     prob.parameters["log_level"] = 0
     prob.optimize()
 
-    p.vector.copy(p_old.vector)
-    epsp.vector.copy(epsp_old.vector)
+    p.x.petsc_vec.copy(p_old.x.petsc_vec)
+    epsp.x.petsc_vec.copy(epsp_old.x.petsc_vec)
 
     sigma = sig(eps_el)
-    sig_exp = fem.Expression(sigma[1, 1], Vp.element.interpolation_points())
+    sig_exp = fem.Expression(sigma[1, 1], Vp.element.interpolation_points)
     sigyy.interpolate(sig_exp)
 
 

docs/references.bib

@@ -19,16 +19,7 @@ @article{halphen1975materiaux
   pages={39--63},
   year={1975}
 }
-@article{ortiz1999variational,
-  title={The variational formulation of viscoplastic constitutive updates},
-  author={Ortiz, Michael and Stainier, Laurent},
-  journal={Computer methods in applied mechanics and engineering},
-  volume={171},
-  number={3-4},
-  pages={419--444},
-  year={1999},
-  publisher={Elsevier}
-}
+
 
 
 % viscoplastic fluids
@@ -50,24 +41,7 @@ @article{coussot2017bingham
   year={2017},
   publisher={Springer}
 }
-@article{bleyer2018advances,
-  title={Advances in the simulation of viscoplastic fluid flows using interior-point methods},
-  author={Bleyer, Jeremy},
-  journal={Computer Methods in Applied Mechanics and Engineering},
-  volume={330},
-  pages={368--394},
-  year={2018},
-  publisher={Elsevier}
-}
-@article{bleyer2015efficient,
-  title={Efficient numerical computations of yield stress fluid flows using second-order cone programming},
-  author={Bleyer, Jeremy and Maillard, Mathilde and De Buhan, Patrick and Coussot, Philippe},
-  journal={Computer Methods in Applied Mechanics and Engineering},
-  volume={283},
-  pages={599--614},
-  year={2015},
-  publisher={Elsevier}
-}
+
 
 
 % limit analysis
@@ -244,14 +218,6 @@ @Article{hewitt2016obstructed
   publisher = {Cambridge University Press},
 }
 
-@Article{halphen1975materiaux,
-  author  = {Halphen, Bernard and Nguyen, Quoc Son},
-  title   = {Sur les mat{\'e}riaux standard g{\'e}n{\'e}ralis{\'e}s},
-  journal = {Journal de m{\'e}canique},
-  year    = {1975},
-  volume  = {14},
-  pages   = {39--63},
-}
 
 @Article{el2020elastoplastic,
   author    = {El Boustani, Chadi and Bleyer, Jeremy and Arquier, Mathieu and Ferradi, Mohammed-Khalil and Sab, Karam},
@@ -597,16 +563,6 @@ @article{makrodimopoulos2006lower
   year={2006},
   publisher={Wiley Online Library}
 }
-@article{makrodimopoulos2007upper,
-  title={Upper bound limit analysis using simplex strain elements and second-order cone programming},
-  author={Makrodimopoulos, A and Martin, CM1196},
-  journal={International journal for numerical and analytical methods in geomechanics},
-  volume={31},
-  number={6},
-  pages={835--865},
-  year={2007},
-  publisher={Wiley Online Library}
-}
 
 @article{nesterov1992conic,
   title={Conic formulation of a convex programming problem and duality},

dolfinx_optim/mosek_io.py

@@ -64,7 +64,7 @@ def _create_mass_matrix(V2, dx):
     # mass matrix on V2
     one = fem.Function(V2)
     v = ufl.TestFunction(V2)
-    one.vector.set(1.0)
+    one.x.petsc_vec.set(1.0)
     m_ufl = ufl.inner(one, v) * dx
     return fem.assemble_vector(fem.form(m_ufl)).array
 
@@ -151,10 +151,10 @@ def _default_parameters(self):
     def _create_variable_vector(self, var, name, ux, lx, cone):
         if cone is not None:
             n = cone.dim
-            m = len(var.vector.array) // n
+            m = len(var.x.petsc_vec.array) // n
             domain = mosek_cone_domain(cone)
             if cone.type == "sdp":
-                m = len(var.vector.array) // n**2
+                m = len(var.x.petsc_vec.array) // n**2
                 domain = mf.Domain.inPSDCone(n, m)
                 if name is not None:
                     vec = self.M.variable(name, domain)
@@ -193,11 +193,11 @@ def _create_variable_vector(self, var, name, ux, lx, cone):
     def _add_boundary_conditions(self, variable, vector, bcs):
         V = variable.function_space
         u_bc = fem.Function(V)
-        u_bc.vector.set(np.inf)
+        u_bc.x.petsc_vec.set(np.inf)
 
-        fem.set_bc(u_bc.vector, to_list(bcs))
-        dof_indices = np.where(u_bc.vector.array < np.inf)[0].astype(np.int32)
-        bc_values = u_bc.vector.array[dof_indices]
+        [bc.set(u_bc.x.array) for bc in to_list(bcs)]
+        dof_indices = np.where(u_bc.x.petsc_vec.array < np.inf)[0].astype(np.int32)
+        bc_values = u_bc.x.petsc_vec.array[dof_indices]
 
         self.M.constraint(vector.pick(dof_indices), mf.Domain.equalsTo(bc_values))
 
@@ -436,7 +436,7 @@ def add_convex_term(self, conv_fun: ConvexTerm):
             cone = conv_fun.cones[i]
 
             vector = self._create_variable_vector(var, name, ux, lx, cone)
-            assert len(var.vector.array) == vector.getSize()
+            assert len(var.x.petsc_vec.array) == vector.getSize()
             self.variables.append(var)
             self.vectors.append(vector)
 
@@ -516,11 +516,11 @@ def _apply_linear_constraints(self, conv_fun):
             if cons["bu"] is not None:
                 if isinstance(cons["bu"], float):
                     bu = fem.Function(cons["V"])
-                    xbu = bu.vector.array
+                    xbu = bu.x.petsc_vec.array
                     xbu[:] = cons["bu"]
                 elif cons["bu"] == 0:
                     bu = fem.Function(cons["V"])
-                    xbu = bu.vector.array
+                    xbu = bu.x.petsc_vec.array
                 else:
                     bu = fem.assemble_vector(
                         fem.form(ufl.inner(lamb_, cons["bu"]) * conv_fun.dx)
@@ -531,11 +531,11 @@ def _apply_linear_constraints(self, conv_fun):
             if cons["bl"] is not None:
                 if isinstance(cons["bl"], float):
                     bl = fem.Function(cons["V"])
-                    xbl = bl.vector.array
+                    xbl = bl.x.petsc_vec.array
                     xbl[:] = cons["bl"]
                 elif cons["bl"] == 0:
                     bl = fem.Function(cons["V"])
-                    xbl = bl.vector.array
+                    xbl = bl.x.petsc_vec.array
                 else:
                     bl = fem.assemble_vector(
                         fem.form(ufl.inner(lamb_, cons["bl"]) * conv_fun.dx)
@@ -715,7 +715,7 @@ def optimize(self, sense="min", dump=False):
         # Retrieve solution and save to file
         for var, vec in zip(self.variables, self.vectors):
             if isinstance(var, fem.Function):
-                var.vector.array[:] = vec.level()
+                var.x.petsc_vec.array[:] = vec.level()
             else:
                 var.value = vec.level()
 
@@ -735,5 +735,5 @@ def get_lagrange_multiplier(self, name):
         """Retrieves Lagrange multiplier function associated with constraint `name`."""
         constraint, V_cons = self.constraints[name]
         lag = fem.Function(V_cons, name=name)
-        lag.vector.array[:] = constraint.dual()
+        lag.x.petsc_vec.array[:] = constraint.dual()
         return lag

docs/index.md → index.md

@@ -6,11 +6,11 @@ This project is a rewrite of the [`fenics_optim` package](https://gitlab.enpc.fr
 
 ## Introduction
 
-```{include} ../README.md
+```{include} README.md
 :start-line: 1
 ```
 
-```{image} images/banner_tutelles.png
+```{image} docs/images/banner_tutelles.png
 :class: bg-primary mb-1
 :width: 600px
 :align: center

pyproject.toml

@@ -0,0 +1,25 @@
+[build-system]
+requires = ["setuptools>=61.2"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "dolfinx_optim"
+version = "0.2.0"
+authors = [{ name = "Jeremy Bleyer", email = "jeremy.bleyer@enpc.fr" }]
+description = "Automated Convex Optimization in FEniCSx"
+readme = "README.md"
+keywords = ["optimization", "PDF"]
+requires-python = ">=3.10"
+dependencies = ["mosek", "fenics-dolfinx>=0.10,<0.11"]
+
+[project.optional-dependencies]
+dev = ["pdbpp", "ipython", "mypy", "ruff"]
+docs = ["jupyter-book<2.0", "sphinx-codeautolink"]
+
+
+[tool.setuptools]
+zip-safe = false
+include-package-data = true
+
+[tool.setuptools.packages]
+find = { namespaces = false }

tests/test_interpolation_operator.py

@@ -24,7 +24,7 @@
 )
 
 deg=1
-V = fem.functionspace(domain, ("CG", deg, ()))
+V = fem.functionspace(domain, ("Lagrange", deg, ()))
 u = fem.Function(V)
 u.interpolate(lambda x: x[0]+2*x[1])