diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py index cb7fcd36f0347947868f5da8ba8d3059d7914f19..8babc7fafb9ef32f7e72f01f69fdeadcdbd8097f 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 3505d4900b44fba1c29e2e2cf4068b00873653a4..ee57dc38b41ae5ca7e7a30ef3f4fcae9c64e8149 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 33ec16d174588a3caf747a5625e456efac98a5b2..c7bf8e61f28ab9174b5161fecc22d34dd84a5831 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 0274f34bc851930113aa4f6f1b9d7160fbd84c7b..5233da0260796b07204e896595c47652758e4a1a 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 d957e81b9ce9d1be97db04b15c9fb63f2f423024..67b279f437168929a50d1ee722e997da3d850fb5 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 a1a56fd68bdccd972027c92dc3c81230053564da..0b79446d98c6ff0b0a0580ae7387241c5ec7dc1b 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)" ] },