From 82ab9adac95f7326d4fad00b55949e0465610dca Mon Sep 17 00:00:00 2001
From: karamoko <karamoko.samassa@inria.fr>
Date: Tue, 18 Oct 2022 09:06:21 +0200
Subject: [PATCH] toward upgrading all primitive domain 2
---
src/bvpy/domains/abstract.py | 13 ++-
src/bvpy/domains/primitives.py | 9 ++-
src/bvpy/utils/post_processing.py | 3 +-
src/bvpy/utils/visu.py | 25 +++++-
tutorials/bvpy_tutorial_1_hello_world.ipynb | 22 +++--
tutorials/bvpy_tutorial_2_domains.ipynb | 89 +++++++++++++++++----
6 files changed, 127 insertions(+), 34 deletions(-)
diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py
index cb7fcd3..8babc7f 100644
--- a/src/bvpy/domains/abstract.py
+++ b/src/bvpy/domains/abstract.py
@@ -304,7 +304,7 @@ class AbstractDomain(ABC):
"""
gmsh.fltk.run()
- def set_cell_size(self, cell_size):
+ def set_cell_size(self, cell_size, cell_size_changed=False):
"""Sets the characteristic size of the domain elements.
Parameters
@@ -316,6 +316,8 @@ class AbstractDomain(ABC):
self._cell_size = cell_size
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", cell_size)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", cell_size)
+ if cell_size_changed:
+ self.synchronize()
def set_meshing_algorithm(self, algorithm):
"""Sets the algorithm to use to discretize the domain elements.
@@ -379,8 +381,11 @@ class AbstractDomain(ABC):
if self.mesh is not None:
self.model.mesh.clear()
- self.model.mesh.generate(3)
+ self.model.mesh.generate(dim)
self.mesh, cell_markers, facet_markers = gmshio.model_to_mesh(self.model, MPI.COMM_SELF, 0)
+ # self.bdata = _generate_default_boundary_data(self.mesh)
+ self.ds = ufl.Measure('ds', domain=self.mesh)
+ self.dx = ufl.Measure('dx', domain=self.mesh)
# self._mesh_conversion(dim=dim)
@@ -396,8 +401,8 @@ class AbstractDomain(ABC):
# with dolfinx.io.XDMFFile(MPI.COMM_WORLD, self._dir_name+'/'+'mesh.xdmf', "w") as xdmf:
# xdmf.write_mesh(self.mesh)
- self.ds = ufl.ds
- self.dx = ufl.dx
+ # self.ds = ufl.ds
+ # self.dx = ufl.dx
def rotate(self, point, axe, angle):
dimTags = [(2, val) for val in self.surfaces.values()] + [(3, val) for val in self.volumes.values()]
diff --git a/src/bvpy/domains/primitives.py b/src/bvpy/domains/primitives.py
index 3505d49..ee57dc3 100644
--- a/src/bvpy/domains/primitives.py
+++ b/src/bvpy/domains/primitives.py
@@ -618,8 +618,15 @@ class Sphere(AbstractDomain, OccModel):
None
"""
- self.volumes[0] = self.factory.addSphere(*center, radius)
+ # self.volumes[0] = self.factory.addSphere(*center, radius)
+ # self.synchronize()
+
+ _dim = 2
+ _tag = [1]
+ self.surfaces[0] = self.factory.addSphere(*center, radius)
self.synchronize()
+ self.model.addPhysicalGroup(_dim, _tag)
+ self.model.mesh.generate(_dim)
class Ellipsoid(AbstractDomain, OccModel):
diff --git a/src/bvpy/utils/post_processing.py b/src/bvpy/utils/post_processing.py
index 33ec16d..c7bf8e6 100644
--- a/src/bvpy/utils/post_processing.py
+++ b/src/bvpy/utils/post_processing.py
@@ -321,8 +321,9 @@ def estimate_precision(prblm, prediction, output='array', **kwargs):
# error_max = np.max(np.abs(expected_sol.x.array-computed_sol.x.array))
- # error_array = fem.Function(func_space)
+ error_array = fem.Function(func_space)
+ error_array.interpolate(expected_sol - computed_sol)
# error_array.x.array = np.abs(expected_sol.x.array-computed_sol.x.array)
error_array = expected_sol - computed_sol
diff --git a/src/bvpy/utils/visu.py b/src/bvpy/utils/visu.py
index 0274f34..5233da0 100644
--- a/src/bvpy/utils/visu.py
+++ b/src/bvpy/utils/visu.py
@@ -20,12 +20,13 @@
#
# -----------------------------------------------------------------------
# import fenics as fe
+import logging
from tokenize import cookie_re
import numpy as np
import plotly.io as pio
import plotly.graph_objects as go
-from bvpy.domains.abstract import AbstractDomain
+from bvpy.domains.abstract import AbstractDomain, OccModel
from bvpy.utils.post_processing import SolutionExplorer
from dolfinx import fem, mesh
@@ -252,10 +253,26 @@ def plot(obj):
_plotter.show()
"""
- _topology, _cell_types, _geometry = dolfinx.plot.create_vtk_mesh(obj.function_space)
+ has_u_data = False
+
+ if issubclass(type(obj), AbstractDomain) :
+ if obj.mesh is None:
+ obj.discretize()
+ _topology, _cell_types, _geometry = dolfinx.plot.create_vtk_mesh(obj.mesh)
+ elif isinstance(obj, mesh.Mesh):
+ _topology, _cell_types, _geometry = dolfinx.plot.create_vtk_mesh(obj)
+ elif isinstance(obj, fem.Function):
+ _topology, _cell_types, _geometry = dolfinx.plot.create_vtk_mesh(obj.function_space)
+ has_u_data = True
+ else:
+ logging.warning("This type of data can't be plotted !")
+
_grid = pyvista.UnstructuredGrid(_topology, _cell_types, _geometry)
- _grid.point_data["u"] = obj.x.array.real
- _grid.set_active_scalars("u")
+
+ if has_u_data:
+ _grid.point_data["u"] = obj.x.array.real
+ _grid.set_active_scalars("u")
+
_plotter = pyvista.Plotter()
_plotter.add_mesh(_grid, show_edges=True)
_plotter.view_xy()
diff --git a/tutorials/bvpy_tutorial_1_hello_world.ipynb b/tutorials/bvpy_tutorial_1_hello_world.ipynb
index d957e81..67b279f 100644
--- a/tutorials/bvpy_tutorial_1_hello_world.ipynb
+++ b/tutorials/bvpy_tutorial_1_hello_world.ipynb
@@ -172,12 +172,12 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "f\n"
+ "<class 'dolfinx.fem.function.Function'>\n"
]
}
],
"source": [
- "print(prblm_1.solution)"
+ "print(type(prblm_1.solution))"
]
},
{
@@ -191,13 +191,13 @@
},
{
"cell_type": "code",
- "execution_count": 6,
+ "execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
- "model_id": "98645b5db9ed42ef87e1ca33d627382b",
+ "model_id": "68c9a98ec27d4ff98adea3a6524aac74",
"version_major": 2,
"version_minor": 0
},
@@ -227,7 +227,7 @@
},
{
"cell_type": "code",
- "execution_count": 8,
+ "execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
@@ -244,7 +244,7 @@
},
{
"cell_type": "code",
- "execution_count": 7,
+ "execution_count": 10,
"metadata": {},
"outputs": [
{
@@ -254,7 +254,7 @@
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)",
- "Cell \u001b[0;32mIn [7], line 5\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[1;32m 3\u001b[0m error \u001b[38;5;241m=\u001b[39m estimate_precision(prblm_1, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m1 + x*x + 2*y*y\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m----> 5\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mMaximal discrepency between expectation and computed solution: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mnp\u001b[38;5;241m.\u001b[39mmax(error\u001b[38;5;241m.\u001b[39mx\u001b[38;5;241m.\u001b[39marray)\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.2e\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n",
+ "Cell \u001b[0;32mIn [10], line 5\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[1;32m 3\u001b[0m error \u001b[38;5;241m=\u001b[39m estimate_precision(prblm_1, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m1 + x*x + 2*y*y\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m----> 5\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mMaximal discrepency between expectation and computed solution: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mnp\u001b[38;5;241m.\u001b[39mmax(error\u001b[38;5;241m.\u001b[39mx\u001b[38;5;241m.\u001b[39marray)\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.2e\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n",
"\u001b[0;31mAttributeError\u001b[0m: 'Sum' object has no attribute 'x'"
]
}
@@ -297,6 +297,14 @@
"File \u001b[0;32m~/miniconda3/envs/bvpyx-fenicsx-dev/lib/python3.10/site-packages/bvpy/utils/visu.py:255\u001b[0m, in \u001b[0;36mplot\u001b[0;34m(obj)\u001b[0m\n\u001b[1;32m 241\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mplot\u001b[39m(obj):\n\u001b[1;32m 242\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 243\u001b[0m \u001b[39m # Error plotting\u001b[39;00m\n\u001b[1;32m 244\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 253\u001b[0m \n\u001b[1;32m 254\u001b[0m \u001b[39m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 255\u001b[0m _topology, _cell_types, _geometry \u001b[39m=\u001b[39m dolfinx\u001b[39m.\u001b[39mplot\u001b[39m.\u001b[39mcreate_vtk_mesh(obj\u001b[39m.\u001b[39;49mfunction_space)\n\u001b[1;32m 256\u001b[0m _grid \u001b[39m=\u001b[39m pyvista\u001b[39m.\u001b[39mUnstructuredGrid(_topology, _cell_types, _geometry)\n\u001b[1;32m 257\u001b[0m _grid\u001b[39m.\u001b[39mpoint_data[\u001b[39m\"\u001b[39m\u001b[39mu\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m obj\u001b[39m.\u001b[39mx\u001b[39m.\u001b[39marray\u001b[39m.\u001b[39mreal\n",
"\u001b[0;31mAttributeError\u001b[0m: 'Sum' object has no attribute 'function_space'"
]
+ },
+ {
+ "ename": "",
+ "evalue": "",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[1;31mLe Kernel s’est bloqué lors de l’exécution du code dans la cellule active ou une cellule précédente. Veuillez vérifier le code dans la ou les cellules pour identifier une cause possible de l’échec. Cliquez <a href='https://aka.ms/vscodeJupyterKernelCrash'>ici</a> pour plus d’informations. Pour plus d’informations, consultez Jupyter <a href='command:jupyter.viewOutput'>log</a>."
+ ]
}
],
"source": [
diff --git a/tutorials/bvpy_tutorial_2_domains.ipynb b/tutorials/bvpy_tutorial_2_domains.ipynb
index a1a56fd..0b79446 100644
--- a/tutorials/bvpy_tutorial_2_domains.ipynb
+++ b/tutorials/bvpy_tutorial_2_domains.ipynb
@@ -142,16 +142,18 @@
"metadata": {},
"outputs": [
{
- "ename": "AttributeError",
- "evalue": "'Rectangle' object has no attribute 'function_space'",
- "output_type": "error",
- "traceback": [
- "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
- "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)",
- "Cell \u001b[0;32mIn [4], line 4\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mbvpy\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mutils\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mvisu\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m plot, set_renderer\n\u001b[1;32m 2\u001b[0m set_renderer(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mnotebook\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m----> 4\u001b[0m plot(domain)\n",
- "File \u001b[0;32m~/miniconda3/envs/bvpyx-fenicsx-dev/lib/python3.10/site-packages/bvpy/utils/visu.py:255\u001b[0m, in \u001b[0;36mplot\u001b[0;34m(obj)\u001b[0m\n\u001b[1;32m 241\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mplot\u001b[39m(obj):\n\u001b[1;32m 242\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 243\u001b[0m \u001b[39m # Error plotting\u001b[39;00m\n\u001b[1;32m 244\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 253\u001b[0m \n\u001b[1;32m 254\u001b[0m \u001b[39m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 255\u001b[0m _topology, _cell_types, _geometry \u001b[39m=\u001b[39m dolfinx\u001b[39m.\u001b[39mplot\u001b[39m.\u001b[39mcreate_vtk_mesh(obj\u001b[39m.\u001b[39;49mfunction_space)\n\u001b[1;32m 256\u001b[0m _grid \u001b[39m=\u001b[39m pyvista\u001b[39m.\u001b[39mUnstructuredGrid(_topology, _cell_types, _geometry)\n\u001b[1;32m 257\u001b[0m _grid\u001b[39m.\u001b[39mpoint_data[\u001b[39m\"\u001b[39m\u001b[39mu\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m obj\u001b[39m.\u001b[39mx\u001b[39m.\u001b[39marray\u001b[39m.\u001b[39mreal\n",
- "\u001b[0;31mAttributeError\u001b[0m: 'Rectangle' object has no attribute 'function_space'"
- ]
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "6bc813e37d6b4ad8b2064d2a454c1157",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
}
],
"source": [
@@ -182,11 +184,31 @@
"cell_type": "code",
"execution_count": null,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "ename": "",
+ "evalue": "",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[1;31mCanceled future for execute_request message before replies were done"
+ ]
+ },
+ {
+ "ename": "",
+ "evalue": "",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[1;31mLe Kernel s’est bloqué lors de l’exécution du code dans la cellule active ou une cellule précédente. Veuillez vérifier le code dans la ou les cellules pour identifier une cause possible de l’échec. Cliquez <a href='https://aka.ms/vscodeJupyterKernelCrash'>ici</a> pour plus d’informations. Pour plus d’informations, consultez Jupyter <a href='command:jupyter.viewOutput'>log</a>."
+ ]
+ }
+ ],
"source": [
"d1 = Rectangle(length=2, width=3, clear=True)\n",
"d1.set_cell_size(.1)\n",
- "plot(d1)"
+ "# d1.model.geo.mesh.set_size_from_boundary(2,1,1)\n",
+ "\n",
+ "plot(d1)\n",
+ "# help(d1.model.geo.mesh.set_size_from_boundary)\n"
]
},
{
@@ -206,9 +228,24 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 5,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "89d512655f71468e8b6bfd28c0495970",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
"source": [
"rec = Rectangle(cell_size=0.1, clear=True)\n",
"plot(rec)"
@@ -259,12 +296,30 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 1,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "8a76846c7d38485c96ef440ee387bc63",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
"source": [
"from bvpy.domains import Sphere\n",
- "sph = Sphere(cell_size=0.1, clear=True)\n",
+ "from bvpy.utils.visu import plot, set_renderer\n",
+ "# sph = Sphere(cell_size=0.1, clear=True)\n",
+ "sph = Sphere()\n",
+ "sph.set_cell_size(0.1)\n",
"plot(sph)"
]
},
--
GitLab