Mentions légales du service

Skip to content
Snippets Groups Projects
Commit b2d13c18 authored by ALI Olivier's avatar ALI Olivier :monkey_face:
Browse files

Fixed bug in CustomDomain to produce 3D meshes.

parent fdd7ebc7
No related branches found
No related tags found
1 merge request!10Draft: Feature/upgrade to fenicsx
......@@ -86,14 +86,19 @@ class CustomDomain(AbstractDomain, BuiltInModel):
# self.mesh = fe.Mesh()
self.mesh = None
self.dim = points.shape[1]
print(f'self.dim = {self.dim}')
# editor = fe.MeshEditor()
if self._cell_type == 'line':
top = 1
topo_dim = 1 # OA: changed variable name top -> topo_dim
# OA: what is the use of this variable ?
elif self._cell_type == 'triangle':
top = 2
topo_dim = 2
elif self._cell_type == 'tetra':
top = 3
# editor.open(self.mesh, self._cell_type, top, self.dim)
topo_dim = 3
# editor.open(self.mesh, self._cell_type, topo_dim, self.dim)
# editor.init_vertices(len(points))
# editor.init_cells(len(cells))
......@@ -101,13 +106,15 @@ class CustomDomain(AbstractDomain, BuiltInModel):
# [editor.add_cell(ind, cell) for ind, cell in enumerate(cells)]
# editor.close()
_cell = ufl.Cell(self._cell_type, geometric_dimension=top)
print(f'pts used to create mesh = {points[:, :self.dim]}')
_cell = ufl.Cell(self._cell_type, geometric_dimension=self.dim)
_domain = ufl.Mesh(ufl.VectorElement("CG", _cell, 1))
self.mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, points[:, :top], _domain)
self.mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells,
points[:, :self.dim], _domain)
# self.cdata = _generate_default_cell_data(self.mesh, )
# self.bdata = _generate_default_boundary_data(self.mesh)
......@@ -143,9 +150,7 @@ class CustomDomain(AbstractDomain, BuiltInModel):
"""
mesh_io = mio.read(path)
_, dim = mesh_io.points.shape
points = mesh_io.points[:, :dim]
points = mesh_io.points
cells = mesh_io.cells_dict["triangle"]
return CustomDomain(points, cells)
%% Cell type:markdown id: tags:
# Example of Domain's Primitives
In this tutorial, we will explore in particular the `bvpy.domains` module. The goal is to get familiar with the concept of integration domain and in particular with its implementation in the **Bvpy** library.
**Covered topics:**
- Basic operations on domains: instanciation, `.info()`, `.discretize()` and `set_cell_size()` methods.
- Quick domain visualization: `plot()` function from the `bvpy.utils.visu` sub-module.
- Basic geometrical domains overview: `Rectangle()`, `Disk()` and other primitives.
- Constructive Solid Geometry (*CSG*) manipulation of geometrical primitives.
- Generation of domains from existing files: `CustomDomain` and `CustomPolygonalDomain` classes.
Download the notebook: [bvpy_tutorial_2](https://gitlab.inria.fr/mosaic/publications/bvpy/-/raw/develop/tutorials/bvpy_tutorial_2_domains.ipynb?inline=false).
%% Cell type:markdown id: tags:
## Basic operations with domains
%% Cell type:markdown id: tags:
### Geometrical domain instanciation
Let's first see how to implement a simple domain and perform basic task with it. Various domain classes are gathered in the `bvpy.domains` modules. Within this module they are sorted in various submodules depending on their characteristics. The `bvpy.domains.primitives` submodule gathers purely geometrical domains built with the **Gmsh** library.
**Remark:** To import a specific domain class you do not need to call the corresponding submodule, calling the `bvpy.domains` module is enough.
%% Cell type:code id: tags:
``` python
from bvpy.domains import Rectangle
domain = Rectangle()
```
%% Cell type:markdown id: tags:
### Domain meshing
With the previous line, we have justed instanciated the default rectangular domain. As mentioned in the documentation [click](file:///Users/oali/Documents/Work/Research/Devlp/bvpy/doc/build/html/library_description.html#domains), in order to mesh the domain with triangles, one needs to call the `discretize()` method:
%% Cell type:code id: tags:
``` python
try:
# domain.mesh.num_cells()
# print(f"Number of triangles within the domain BEFORE discretization: {domain.mesh.num_cells()}")
number_triangles = domain.mesh.topology.index_map(2).size_global
print(f"Number of triangles within the domain BEFORE discretization: {number_triangles}")
except AttributeError:
domain.discretize()
# print(f"Number of triangles within the domain AFTER discretization: {domain.mesh.num_cells()}")
number_triangles = domain.mesh.topology.index_map(2).size_global
print(f"Number of triangles within the domain AFTER discretization: {number_triangles}")
```
%% Cell type:markdown id: tags:
### Get information about the considered domain
Useful informations about the domain are available through the `info()` method:
%% Cell type:code id: tags:
``` python
domain.info()
```
%% Cell type:markdown id: tags:
### Quick visualization of the meshed domains
A quick visualization can be performed thanks to the `plot()` function from the `bvpy.utils.visu` sub-module:
%% Cell type:code id: tags:
``` python
from bvpy.utils.visu import plot, set_renderer
set_renderer('notebook')
plot(domain)
```
%% Cell type:markdown id: tags:
**Remarks:**
* The `plot()`function requires the existence of the mesh within the domain but can be called upon a domain that has not been discretized... In this case the plot function calls itself the `domain.discretize()` method.
* The `set_renderer('notebook')` is a required instruction to enable the encapsulation of plotly objects in notebooks.
%% Cell type:markdown id: tags:
### Resizing domains
Size and resolution can be adjusted at will, either at instanciation or afterwards:
%% Cell type:code id: tags:
``` python
from bvpy.utils.visu import plot, set_renderer
from bvpy.domains import Rectangle
set_renderer('notebook')
d1 = Rectangle(length=2, width=3, clear=True)
d1.set_cell_size(.05)
plot(d1)
```
%% Cell type:markdown id: tags:
## Available geometrical primitives
Let's now have a look at the geometrical primitives that we can instantiate thanls to the `bvpy.domains.primitives` sub-module.
%% Cell type:markdown id: tags:
### Rectangle
%% Cell type:code id: tags:
``` python
rec = Rectangle(length=2, width=3, cell_size=0.1, clear=True)
plot(rec)
```
%% Cell type:markdown id: tags:
### Disk
%% Cell type:code id: tags:
``` python
from bvpy.domains import Disk
dsk = Disk(cell_size=0.1, clear=True)
plot(dsk)
```
%% Cell type:markdown id: tags:
### HemiSphere
%% Cell type:code id: tags:
``` python
from bvpy.domains import HemiSphere
from bvpy.utils.visu import plot, set_renderer
hsph = HemiSphere(cell_size=0.1, clear=True)
plot(hsph)
```
%% Cell type:markdown id: tags:
### Sphere
%% Cell type:code id: tags:
``` python
from bvpy.domains import Sphere
from bvpy.utils.visu import plot, set_renderer
sph = Sphere(cell_size=0.1, clear=True)
plot(sph)
```
%% Cell type:markdown id: tags:
### Ellipsoid
%% Cell type:code id: tags:
``` python
from bvpy.domains import Ellipsoid
elp = Ellipsoid(cell_size=0.1, clear=True)
plot(elp)
```
%% Cell type:markdown id: tags:
### Torus
%% Cell type:code id: tags:
``` python
from bvpy.domains import Torus
tor = Torus(cell_size=0.2, clear=True)
plot(tor)
```
%% Cell type:markdown id: tags:
## Combining geometrical primitives
We included Constructive Solid Geometry (CSG) tools from **Gmsh**, in order to easily combine geometrical primitives into more complex forms.
Three operators from **Gmsh** have been implemented within the `bvpy.domains` module:
* Addition: `__add__()`
* Substraction: `__sub__()`
* Intersection: `__and__()`
**Remark:** The classic operators `+`, `-` and `&` have been surcharged to ease even more the combinaison of shapes.
%% Cell type:markdown id: tags:
### Combining flat surfaces
%% Cell type:code id: tags:
``` python
from bvpy.utils.visu import plot, set_renderer
from bvpy.domains import Rectangle
from bvpy.domains import Disk
set_renderer('notebook')
dsk1 = Disk(center=[0, -.25, 0], clear=True, cell_size=.1)
dsk2 = Disk(center=[0, .25, 0])
dsk3 = Disk(radius=.3)
rec1 = Rectangle(y=-.05, x=-.5, length=1, width=.1)
cmplx_dom = dsk1 & dsk2 - (dsk3 + rec1)
plot(cmplx_dom)
```
%% Cell type:markdown id: tags:
> **Note:** You can see in the first line above we added the argument `clear=True` at the first domain instanciation. When multiple domains are instanciated, they share the same **Gmsh** factory and are, so to speak, “super-imposed" on one another. This can create conflicts latter on when domains are processed. For instance, when the `plot` function probes the factory for a specific domain, It might not return the desired one but all the previously instanced ones contained within the factory. to avoid this, we added a call to the `gmsh.clear()` (triggered with the argument `clear` is set to `True`) method within our `AbstractDomain` class. This **Gmsh** method resets the domain factory and avoids visualization conflits.
%% Cell type:markdown id: tags:
### Combining curved surfaces:
%% Cell type:code id: tags:
``` python
from bvpy.domains import Sphere
from bvpy.utils.visu import plot, set_renderer
sph1 = Sphere(clear=True)
sph2 = Sphere(center=[0,1,0], radius=.5)
pin = sph1 + sph2
pin.set_cell_size(.1)
plot(pin)
```
%% Cell type:markdown id: tags:
### Combining volumes
This CSG feature can also be used to combine volumes together and carve complex 3D domains.
%% Cell type:code id: tags:
``` python
big = Sphere(radius=10, cell_type='tetra', clear=True)
sml = Sphere(radius=9, cell_type='tetra')
shell = big - sml
plot(shell)
```
%% Cell type:markdown id: tags:
## Importing domains from elsewhere
One can be able to import within **bvpy** his/her own custom domain of interest. To do so, we enabled the generation of domains from existing meshes, recorded in the widespread `.ply` format.
The `bvpy.domains.custom_domain` sub-module contains the `CustomDomain` class that takes as input such an existing mesh structure and generate a proper *bvpy domain* from it.
**Remarks:**
* Once again, the subdomain `custom_domain` does not have to be mentioned within the import command.
* The `read` static method of the `CustomDomain` class relies on the **meshio** library. The `.ply` file used as imput should not contain too much properties. And their header should look like:
```
ply
format ascii 1.0
element vertex 7359
property float x
property float y
property float z
element face 14304
property list int int vertex_indices
end_header
```
%% Cell type:code id: tags:
``` python
from bvpy.domains import CustomDomain
from bvpy.utils.visu import plot
path = '../data/example_meristem.ply'
cd = CustomDomain.read(path)
cd.info()
plot(cd)
```
%% Output
Domain
------
* Shape: CustomDomain
* Dimension: 3
Meshing
-------
* Algorithm: Delaunay
* Cell type: triangle
* Number: 3
* Resolution (i.e. prescribed element size): 1.00e+00
* Actual size (min, max): (2.24e+00, 2.24e+00)
* Cells quality: 4.50e+00 +/- 2.87e+00
%% Cell type:markdown id: tags:
## Importing piecewise polygonal surfaces
In a biological context it is important to be able to handle data structure that account for cellularized tissues.
This can be done with the `CustomPolygonalDomain` class from the `bvpy.domains.custom_polygonal` sub-module.
%% Cell type:code id: tags:
``` python
from bvpy.utils.visu import plot
from bvpy.domains import CustomPolygonalDomain
path ='../data/tissue_example_curved_big.txt'
tissue = CustomPolygonalDomain.read(path)
tissue.discretize()
tissue.info()
plot(tissue)
```
%% Output
WARNING:root:Adding physical group of Hemisphere...
Domain
------
* Shape: CustomPolygonalDomain
* Dimension: 3
* Maximal width: 9.19e+00
* Number of polygons: 6.10e+01
* Averaged polygon width: 1.58e+00
Meshing
-------
* Algorithm: Delaunay
* Cell type: triangle
* Number: 3
* Resolution (i.e. prescribed element size): 1.00e-01
* Actual size (min, max): (9.40e-02, 9.40e-02)
* Cells quality: 4.50e+00 +/- 2.87e+00
Le 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>.
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment