From 659e7466a74d926adc85e98b926d8c717d7b74f2 Mon Sep 17 00:00:00 2001
From: Gonzalo Revilla <gonzalo.revilla-mut@inria.fr>
Date: Thu, 30 Jan 2025 19:14:46 +0100
Subject: [PATCH] Refractor the entire domain submodule
---
doc/library_description.rst | 4 +-
doc/user_guide.md | 4 +-
src/bvpy/domains/__init__.py | 7 +-
src/bvpy/domains/abstract.py | 854 ++++--------------------
src/bvpy/domains/custom_domain.py | 134 ----
src/bvpy/domains/custom_gmsh.py | 935 +++++++++++++++++++++++++--
src/bvpy/domains/custom_polygonal.py | 246 -------
src/bvpy/domains/fixed_domain.py | 82 +++
src/bvpy/domains/fixed_gmsh.py | 58 ++
src/bvpy/domains/geometry.py | 25 +-
src/bvpy/domains/primitives.py | 894 -------------------------
src/bvpy/utils/examples.py | 10 +-
test/data/mesh_example.msh | 754 ++++++++-------------
test/test_bvp.py | 5 +-
test/test_domains.py | 297 +++++----
test/test_elasticity.py | 6 +-
test/test_gbvp.py | 5 +-
test/test_geometry.py | 43 +-
test/test_gmshio.py | 2 +
test/test_helmholtz.py | 7 +-
test/test_ibvp.py | 5 +-
test/test_io.py | 29 +-
test/test_periodic.py | 5 +-
test/test_plasticity.py | 7 +-
test/test_poisson.py | 17 +-
test/test_post_process.py | 4 +-
test/test_transport.py | 26 +-
test/test_visu.py | 7 +-
test/test_visu_pyvista.py | 5 +-
29 files changed, 1710 insertions(+), 2767 deletions(-)
delete mode 100644 src/bvpy/domains/custom_domain.py
delete mode 100644 src/bvpy/domains/custom_polygonal.py
create mode 100644 src/bvpy/domains/fixed_domain.py
create mode 100644 src/bvpy/domains/fixed_gmsh.py
delete mode 100644 src/bvpy/domains/primitives.py
diff --git a/doc/library_description.rst b/doc/library_description.rst
index 2d7014e..ab4fa89 100644
--- a/doc/library_description.rst
+++ b/doc/library_description.rst
@@ -60,8 +60,8 @@ The notion of *geometry* and the related ``geometry()`` method are of paramount
* The classes contained within the ``bvpy.domains.primitives`` sub-module are built around ``geometry()`` methods that construct basic geometric elements (such as rectangles, disk, spheres...) exclusively with **Gmsh** tools.
.. note:: Basic operators such as ``__add__()``, ``__sub__()`` and ``__and__()`` have been surcharged so primitives can be combined together to form more complex domains.
-* The ``CustomPolygonalDomain`` class (contained within the ``bvpy.domains.custom_polygonal`` sub-module) is built around a ``geometry()`` method that construct piecewise polygonal structures from lists of oriented points. This is also done thanks to **Gmsh** tools. The purpose of this class is to enable the use of cellularized structures as integration domains. This is an important feature in the context of modeling biological tissues and other multicellular continua.
-* Finally, the ``CustomDomain`` class (within the ``bvpy.domains.custom_domain`` sub-module) features an empty ``geometry()`` method for its purpose is to define a domain from an already existing mesh, sorted as a ``.ply``-formatted file.
+* The ``ConformalPolygon`` class (contained within the ``bvpy.domains.conformal_polygon`` sub-module) is built around a ``geometry()`` method that construct piecewise polygonal structures from lists of oriented points. This is also done thanks to **Gmsh** tools. The purpose of this class is to enable the use of cellularized structures as integration domains. This is an important feature in the context of modeling biological tissues and other multicellular continua.
+* Finally, the ``FixedDomain`` class (within the ``bvpy.domains.fixed_domain`` sub-module) its purpose is to define a domain from an already existing mesh, sorted as a ``.ply``-formatted file.
Besides a *geometry* and a *mesh*, domain classes also encompasses ``fenics.Measures`` elements defined within the domain (``self.dx``) and onto its border (``self.ds``) and some *setter* methods to parametrize
diff --git a/doc/user_guide.md b/doc/user_guide.md
index df1cd2e..98dfa65 100644
--- a/doc/user_guide.md
+++ b/doc/user_guide.md
@@ -98,8 +98,8 @@ The `bvpy.domains` module contains the classes that can be used to implement the
The notion of *geometry* and the related `geometry()` method are of paramount importance, for each domain class corresponds to a specific implementation of the `geometry()` method:
* The classes contained within the `bvpy.domains.primitives` sub-module are built around `geometry()` methods that construct basic geometric elements (such as rectangles, disk, spheres...) exclusively with **Gmsh** tools.
-* The `CustomPolygonalDomain` class (contained within the `bvpy.domains.custom_polygonal` sub-module) is built around a `geometry()` method that construct piecewise polygonal structures from lists of oriented points. This is also done thanks to **Gmsh** tools. The purpose of this class is to enable the use of cellularized structures as integration domains. This is an important feature in the context of modeling biological tissues and other multicellular continua.
-* Finally, the `CustomDomain` class (within the `bvpy.domains.custom_domain` sub-module) features an empty `geometry()` method for its purpose is to define a domain from an already existing mesh, sorted as a `.ply`-formatted file.
+* The `ConformalPolygon` class (contained within the `bvpy.domains.conformal_polygon` sub-module) is built around a `geometry()` method that construct piecewise polygonal structures from lists of oriented points. This is also done thanks to **Gmsh** tools. The purpose of this class is to enable the use of cellularized structures as integration domains. This is an important feature in the context of modeling biological tissues and other multicellular continua.
+* Finally, the `FixedDomain` class (within the `bvpy.domains.fixed_domain` sub-module) its purpose is to define a domain from an already existing mesh, sorted as a `.ply`-formatted file.
> Besides a *geometry* and a *mesh*, domain classes also encompasses `fenics.Measures` elements defined within the domain (`self.dx`) and onto its border (`self.ds`) and some *setter* methods to parametrize
diff --git a/src/bvpy/domains/__init__.py b/src/bvpy/domains/__init__.py
index 922b336..e4cea4a 100644
--- a/src/bvpy/domains/__init__.py
+++ b/src/bvpy/domains/__init__.py
@@ -1,4 +1,3 @@
-from .custom_domain import *
-from .custom_polygonal import *
-from .custom_gmsh import *
-from .primitives import *
+from .fixed_domain import FixedDomain
+from .fixed_gmsh import FixedGMSH
+from .custom_gmsh import CustomGMSH, OccModel, BuiltInModel
diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py
index f060b98..dac3e82 100644
--- a/src/bvpy/domains/abstract.py
+++ b/src/bvpy/domains/abstract.py
@@ -9,6 +9,7 @@
# File contributor(s):
# Florian Gacon <florian.gacon@inria.fr>
# Olivier Ali <olivier.ali@inria.fr>
+# Gonzalo Revilla <gonzalo.revilla-mut@inria.fr>
#
# File maintainer(s):
# Olivier Ali <olivier.ali@inria.fr>
@@ -19,179 +20,78 @@
# https://www.gnu.org/licenses/lgpl-3.0.en.html
#
# -----------------------------------------------------------------------
-import os
-import gmsh
-import sys
-
-from bvpy import logger as log
-import meshio as mio
-from abc import ABC, abstractmethod
-import tempfile
-
-import fenics as fe
import numpy as np
+import fenics as fe
-# Available meshing algorithm from gmsh (experimental algorithms are removed)
-# See https://gmsh.info/doc/texinfo/gmsh.html#Choosing-the-right-unstructured-algorithm
-# to have more information about what algo fits your need. Briefy :
-# - complex curved surfaces (2D) --> MeshAdapt
-# - high element quality (2D/3D) --> Frontal-Delaunay
-# - very large meshes --> Delaunay (2D) / HXT (3D)
-AVAILABLE_MESHING_ALGO_2D = {'MeshAdapt': 1, 'Automatic': 2, 'Delaunay': 5, 'Frontal-Delaunay': 6}
-AVAILABLE_MESHING_ALGO_3D = {'Delaunay': 1, 'Frontal-Delaunay': 4, 'HXT': 10}
+from bvpy import logger
# Available cell type for the mesh (hybrid mesh are not allowed)
-AVAILABLE_CELL_TYPE = {'tetra': 3, 'triangle': 2, 'line': 1, 'vertex': 0} # cell_type : topological_dimension
-ELEMENT_GMSH_TYPES = {'line': 1, 'triangle': 2, 'tetra': 4} # Gmsh code for a given element type
+AVAILABLE_CELL_TYPE = {'tetrahedron': 3, 'tetra': 3, 'triangle': 2, 'interval': 1, 'line': 1, 'vertex': 0}
-def _generate_default_boundary_data(mesh: fe.Mesh) -> fe.MeshFunction:
- """Sets the boundary of a mesh.
-
- Parameters
- ----------
- mesh : :class:`Mesh<Dolfin.cpp.mesh.Mesh>`
- The mesh on which we want to define boundaries.
-
- Returns
- -------
- :class:`MeshFunction<Dolfin.cpp.mesh.MeshFunction>`
- A fenics MeshFunction that marks (with the index 0)
- the boundary of the considered mesh.
- """
- mf = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
- mf.set_all(sys.maxsize)
- subdomain = fe.CompiledSubDomain('on_boundary')
- subdomain.mark(mf, 0)
- return mf
-
-
-def _generate_default_cell_data(mesh: fe.Mesh) -> fe.MeshFunction:
- """Sets the body (i.e. inner part) of a mesh.
-
- Parameters
- ----------
- mesh : :class:`Mesh<Dolfin.cpp.mesh.Mesh>`
- The mesh on which we want to mark the body.
-
- Returns
- -------
- :class:`MeshFunction<Dolfin.cpp.mesh.MeshFunction>`
- A fenics MeshFunction that marks (with the index 0 ?)
- the boundary of the considered mesh.
- """
- mf = fe.MeshFunction("size_t", mesh, mesh.topology().dim())
- mf.set_all(0)
- return mf
-
-
-class AbstractDomain(ABC):
+class AbstractDomain:
"""The most general class that defines an integration domain.
Attributes
----------
+ tdim : int
+ Topological dimension of the domain.
+ gdim : int
+ Dimension of the embedded space, also called geometrical or physical dimension
+ cell_type : str
+ The type of cell used for the mesh.
+ mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
+ A mesh representation of the considered domain.
cdata : :class:`MeshFunctionSizet<dolfin.cpp.mesh.MeshFunctionSizet>`
- Description of attribute `cdata`.
+ Function that assigns to each cell a label
bdata : :class:`MeshFunctionSizet<dolfin.cpp.mesh.MeshFunctionSizet>`
- Description of attribute `bdata`.
+ Function that assigns to each facet a label
+ sub_domain_names : dict
+ A dictionary that maps the cdata labels to a name
+ sub_boundary_names : dict
+ A dictionary that maps the bdata labels to a name
dx : :class:`Measure<ufl.measure.Measure>`
An integral measure within the domain of interest.
ds : :class:`Measure<ufl.measure.Measure>`
An integral measure on the border of the domain of interest.
- mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
- A mesh representation (topology and geometry) of the considered domain.
- model : :class:`Model<gmsh.model>`
- A purely geometrical description of the domain.
- domain_type : str
- Name of the geometrical model of the domain.
- details : str
- A printable description of the domain main characteristics.
+
+ Methods
+ -------
+ info()
+ Prints the description of the domain.
+ set_mesh_functions()
+ Sets cdata, bdata, sub_domain_names, sub_boundary_names, ds and dx attributes.
+ move(displacement, return_cdata=True)
+ Create a FeniCS mesh copy and move it using a displacement field.
+ plot()
+ Plots the mesh.
"""
- _model = gmsh.model
- _domain_instances = 0
def __init__(self,
- cell_type: str = 'triangle',
- cell_size=1,
- verbosity: bool = False,
- algorithm: str = 'Delaunay',
- dim: int = 3,
- clear: bool = False):
- """Generates an instance of the AbstractDomain class.
-
+ cell_type: str,
+ gdim: int):
+ """
Parameters
----------
cell_type : str
Type of element used to tile the domain.
- (the default is 'triangle').
- cell_size : float, [float, float]
- Characteristic size of the tiling elements (the default is 1). You can relax this constraint by defining
- cell_size as [min_cell_size, max_cell_size].
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements.
- (the default is 'Delaunay').
-
- Note
- ----
- The following cell types can be used to tile the considered domain:
- * 'line' (1D element)
- * 'triangle' (2D)
- * 'tetra' (3D).
-
- Accepted meshing algorithms depends on the cell_type dimension:
- * 2D ('triangle') : 'Delaunay', 'Frontal-Delaunay', 'MeshAdapt', 'Automatic'
- * 3D ('tetra') : 'Delaunay', 'Frontal-Delaunay', 'HXT'
+ It can be 'tetrahedron', 'tetra', 'triangle', 'interval' or 'line'.
+ gdim : int
+ Dimension of the embedded space, also called geometrical or physical dimension.
"""
- super().__init__(cell_type=cell_type, cell_size=cell_size,
- verbosity=verbosity, algorithm=algorithm)
-
- AbstractDomain._domain_instances += 1
- if AbstractDomain._domain_instances == 1:
- gmsh.initialize()
-
- if clear:
- gmsh.clear()
-
- self._characteristic_size = None # TODO : do we need to init these values here ?
- self._cell_type = cell_type
- self._cell_size = cell_size
- self._algorithm = algorithm
-
- self._mesh_dir = tempfile.TemporaryDirectory()
- self._dir_name = self._mesh_dir.name
-
- self.set_cell_type(cell_type)
- self.set_verbosity(verbosity)
- self.set_cell_size(cell_size)
- self.set_meshing_algorithm(algorithm)
- self.set_dim(dim)
+ self.__cell_type = self.__check_cell_type(cell_type)
+ self.__tdim = self.__set_tdim()
+ self.__gdim = self.__check_gdim(gdim)
- self.dim = dim
self.mesh = None
self.cdata = None
self.bdata = None
- self.sub_domain_names = None # specify explicit names for labels in cdata (MeshFunctionSizet)
- self.sub_boundary_names = None # specify explicit names for labels in bdata (MeshFunctionSizet)
-
-
- self.dx = fe.dx
- self.ds = fe.ds
-
- self.surfaces = dict()
- self.volumes = dict()
-
- def __del__(self):
- """Clears the temporary file containing the domain.
-
- """
- self._mesh_dir.cleanup()
- AbstractDomain._domain_instances -= 1
- if AbstractDomain._domain_instances == 0:
- gmsh.finalize()
+ self.sub_domain_names = {} # cdata labels names
+ self.sub_boundary_names = {} # bdata labels names
+ self.dx = None
+ self.ds = None
def __repr__(self):
"""Returns a short description of the considered domain.
@@ -215,48 +115,37 @@ class AbstractDomain(ABC):
"""
return self._details
+ def info(self):
+ """Prints the description of the domain."""
+ print(self._details)
+
@property
def domain_type(self):
- """str: Name of the geometrical model of the domain.
-
- """
+ """Class name of the domain."""
return self.__class__.__name__
- @property
- def model(self):
- """:class:`Model<gmsh.model>`: The (gmsh) geometrical model used
- to define the domain.
-
- """
- return AbstractDomain._model
-
@property
def _details(self):
- """str: A detailed printable description of the domain.
-
- """
+ """A detailed printable description of the domain."""
description = f'''Domain
------
- * Shape: {self.domain_type}
- * Dimension: {self.dim}\n'''
-
- if self._characteristic_size is not None:
- for name, val in self._characteristic_size.items():
- description += f' * {name}: {val:.2e}\n'
+ * Type: {self.domain_type}
+ * Geometric dimension: {self.gdim}
+ * Topological dimension: {self.tdim}\n'''
if self.mesh is not None:
- description += '\n' + f'''Meshing
--------
- * Algorithm: {self._algorithm}
- * Cell type: {self._cell_type}
- * Number: {len(self.mesh.cells())}
- * Resolution (i.e. prescribed element size): {self._cell_size:.2e}
- * Actual size (min, max): ({self.mesh.hmin():.2e},\
+ description += '\n' + f'''Mesh
+----
+ * Cell type: {self.cell_type}
+ * Number of vertices: {self.mesh.topology().size(0)}
+ * Number of cells: {self.mesh.topology().size(self.tdim)}
+ * Number of facets: {self.mesh.topology().size(self.tdim - 1)}
+ * Number of facets on the boundary: {fe.BoundaryMesh(self.mesh, "exterior").topology().size(self.tdim - 1)}
+ * Mesh size (min, max): ({self.mesh.hmin():.2e},\
{self.mesh.hmax():.2e})
'''
- cells_quality = np.array([c.radius_ratio()
- for c in fe.cells(self.mesh)])
+ cells_quality = np.array([c.radius_ratio() for c in fe.cells(self.mesh)])
mn_c_qual = cells_quality.mean()
std_c_qual = cells_quality.std()
@@ -264,178 +153,104 @@ class AbstractDomain(ABC):
description += f" {mn_c_qual:.2e} +/- {std_c_qual:.2e} \n"
else:
- description = "Mesh not generate yet."
- description += " Use self.discretize() to to so."
+ description += "\nMesh not generated.\n"
return description
- def synchronize(self):
- """synchronizes the CAD model with the gmsh model.
-
- """
- self.factory.synchronize()
-
- @abstractmethod
- def geometry(self, *args, **kwargs):
- """Computes the geometry of the domain.
-
- Parameters
- ----------
- None
-
- Returns
- -------
- None
-
- Note
- ----
- This abstract method must be defined within daughter class
- depending on the specific need.
-
- """
- pass
-
- def run_gui(self):
- """Generates a quick gmsh gui with the current model in it.
-
- """
- gmsh.fltk.run()
-
- def set_dim(self, dim: int):
- """Sets the dimension of the embedded space
-
- Parameters
- ----------
- dim : int
- Dimension of the embedded space ?
-
- """
- if (dim < AVAILABLE_CELL_TYPE[self._cell_type]) or (type(dim) is not int):
- raise ValueError(f'Dimension of the embedded space need to be equal or larger than the dimension'
- f' of the finite element. Here, you are using "{self._cell_type}" meaning'
- f' that dim >= {AVAILABLE_CELL_TYPE[self._cell_type]}.')
-
- self.dim = dim
-
- def set_cell_size(self, cell_size):
- """Sets the characteristic size of the domain elements.
-
- Parameters
- ----------
- cell_size : float, [float, float]
- The cell_size value we want to use.
+ @property
+ def cell_type(self):
+ """The type of cell used for the mesh.
+ Notes
+ -----
+ Currently bvpy does not accept hybrid meshes. The cell type should be one of the following:
+ - 'tetrahedron' or 'tetra' for 3D meshes
+ - 'triangle' for 2D meshes
+ - 'interval' or 'line' for 1D meshes
"""
- if isinstance(cell_size, (float, int)):
- gmsh.option.setNumber("Mesh.CharacteristicLengthMin", cell_size)
- gmsh.option.setNumber("Mesh.CharacteristicLengthMax", cell_size)
- elif (isinstance(cell_size, list) and (len(cell_size) == 2)
- and all(isinstance(ele, (int, float)) for ele in cell_size)):
- gmsh.option.setNumber("Mesh.CharacteristicLengthMin", cell_size[0])
- gmsh.option.setNumber("Mesh.CharacteristicLengthMax", cell_size[1])
- else:
- raise ValueError('cell_size can be a number or a list of two numbers')
+ return self.__cell_type
- self._cell_size = cell_size
- def set_cell_type(self, cell_type: str):
- """Sets the type of finite element (triangle, tetra) you want to use for your mesh
+ @staticmethod
+ def __check_cell_type(cell_type: str):
+ """Checks that the cell_type is valid
Parameters
----------
- cell_type: str
- The cell_type you want to use
+ cell_type : int
"""
if cell_type not in AVAILABLE_CELL_TYPE:
raise ValueError(f'Invalid cell_type, should be one of {AVAILABLE_CELL_TYPE}')
+ return cell_type
- self._cell_type = cell_type
-
- def set_cell_size_from_curvature(self, n_elements_per_circle: int):
-
- if type(n_elements_per_circle) is not int:
- raise ValueError('n_elements should be an integer')
+ @property
+ def gdim(self):
+ """Dimension of the embedded space, also called geometrical or physical dimension."""
+ return self.__gdim
- gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", n_elements_per_circle)
+ @property
+ def dim(self):
+ """Old geometry dimension"""
+ logger.warning("DEPRECATION WARNING: dim will be deprecated in future bvpy versions, use gdim instead")
+ return self.__gdim
- def set_meshing_algorithm(self, algorithm: str):
- """Sets the algorithm to use to discretize the domain elements.
+ def __check_gdim(self, gdim: int):
+ """Checks that the dimension of the embedded space is consistent with the topological dimension
Parameters
----------
- algorithm : str
- Name of the tiling algorithm to use.
+ gdim : int
+ Dimension of the embedded space
"""
- if AVAILABLE_CELL_TYPE[self._cell_type] == 2: # topological_dimension = 2
- if algorithm not in AVAILABLE_MESHING_ALGO_2D.keys():
- raise ValueError(f'Invalid meshing algorithm, can be {list(AVAILABLE_MESHING_ALGO_2D)}')
-
- gmsh.option.setNumber("Mesh.Algorithm", AVAILABLE_MESHING_ALGO_2D[algorithm])
-
- elif AVAILABLE_CELL_TYPE[self._cell_type] == 3: # topological_dimension = 3
- if algorithm not in AVAILABLE_MESHING_ALGO_3D.keys():
- raise ValueError(f'Invalid meshing algorithm, can be {list(AVAILABLE_MESHING_ALGO_3D)}')
-
- gmsh.option.setNumber("Mesh.Algorithm", AVAILABLE_MESHING_ALGO_3D[algorithm])
-
- self._algorithm = algorithm
+ if (gdim < self.tdim) or (type(gdim) is not int) or gdim > 3:
+ raise ValueError(f'Dimension of the embedded space need to be equal or larger than the dimension'
+ f' of the topological space, and less than 4. Here, you are using '
+ f'"{self.tdim}" as tdim, and {gdim} as gdim.')
+ return gdim
+ @property
+ def tdim(self):
+ """Topological dimension of the domain."""
+ return self.__tdim
- def set_verbosity(self, verbosity: bool):
- """Sets the verbosity of the considered domain.
+ def __set_tdim(self):
+ return AVAILABLE_CELL_TYPE[self.__cell_type]
- Parameters
- ----------
- verbosity : bool
-
- """
- if not isinstance(verbosity, bool):
- raise ValueError('verbosity should be a boolean {True, False}')
+ def set_mesh_functions(self):
+ """Sets cdata, bdata, sub_domain_names, sub_boundary_names, ds and dx attributes."""
+ assert self.mesh is not None, "The mesh has not been generated yet."
- if verbosity:
- gmsh.option.setNumber("General.Verbosity", 5)
+ # - Assert that the entities are allowed
+ cell_type = self.mesh.ufl_cell().cellname()
+ if cell_type not in AVAILABLE_CELL_TYPE.keys():
+ raise RuntimeError(f'Input mesh contains an entity that is not handled : {cell_type}')
+
+ if self.cdata is None:
+ self.cdata = fe.MeshFunction("size_t", self.mesh, self.tdim)
+ self.cdata.set_all(1)
+ self.sub_domain_names = {1: 'domain'}
+
+ boundary = fe.CompiledSubDomain('on_boundary')
+ if self.bdata is None:
+ self.bdata = fe.MeshFunction("size_t", self.mesh, self.tdim - 1)
+ self.bdata.set_all(0)
+ boundary.mark(self.bdata, 1)
+ self.sub_boundary_names = {1: 'boundary'}
else:
- gmsh.option.setNumber("General.Verbosity", 0)
-
- def discretize(self):
- """Generates a fenics mesh representing the domain.
-
- Yields
- ------
- mesh : :class:`Mesh <dolfin.cpp.mesh.Mesh>`
- A fenics mesh corresponding to the domain.
- bdata : :class:`MeshFunctionSizet <dolfin.cpp.mesh.MeshFunction>`
- A MeshFunction that labels the borders (i.e., outermost facets)
- of the mesh of the domain.
- cdata : :class:`MeshFunctionSizet <dolfin.cpp.mesh.MeshFunction>`
- A MeshFunction that labels the inner elements (i.e., cells)
- of the mesh of the domain.
- ds : :class:`Measure <ufl.measure.Measure>`
- An integral measure on the border of the domain of interest.
- dx : :class:`Measure <ufl.measure.Measure>`
- An integral measure within the domain of interest.
-
- Notes
- -----
- - The method also generates a temporary gmsh mesh
- and stores it temporarily on disk.
- - MeshFunctionSizet: `Sizet` means that the output of
- the function is a positive integer.
- """
- if self.mesh is not None:
- self.model.mesh.clear()
-
- self.generate()
- self.convert_gmsh_to_fenics()
+ # Check that the boundary and non-boundary tags are distinct
+ boundary_mask = fe.MeshFunction("bool", self.mesh, self.tdim - 1)
+ boundary.mark(boundary_mask, True) # True on the boundary, False elsewhere
+ boundary_values = set(self.bdata.array()[boundary_mask.array()])
+ non_boundary_values = set(self.bdata.array()[~boundary_mask.array()])
+ assert not (boundary_values & non_boundary_values), (
+ f"Some boundary facets and non-boundary facets have the same tag. Please correct it."
+ )
self.ds = fe.Measure('ds', domain=self.mesh, subdomain_data=self.bdata)
self.dx = fe.Measure('dx', domain=self.mesh, subdomain_data=self.cdata)
- def rotate(self, point, axe, angle):
- dimTags = [(2, val) for val in self.surfaces.values()] + [(3, val) for val in self.volumes.values()]
- self.factory.rotate(dimTags, *point, *axe, angle)
- self.synchronize()
+ return
def move(self, displacement, return_cdata=True):
""" Create a FeniCS mesh copy and move it using a displacement field
@@ -470,409 +285,12 @@ class AbstractDomain(ABC):
else:
return new_mesh
- def generate(self):
- """ Generate the Gmsh mesh from the current model and write it in a temporary file.
+ def plot(self, **kwargs):
"""
- # - Get the finite element topological dimension
- ele_dim = AVAILABLE_CELL_TYPE[self._cell_type]
- log.debug(f"Convert a mesh with topological element of dimension {ele_dim}")
-
- self.model.mesh.generate(ele_dim)
- gmsh.write(self._dir_name + '/mesh.msh')
-
- def convert_gmsh_to_fenics(self, fname=None):
- """ Converts a temporary gmsh mesh into a fenics (dolphin) one.
-
- Other Parameters
- ----------------
- fname : str, optional
- The path to the gmsh mesh. Default is the automatic path generated by `generate()` method.
-
- Returns
- -------
- mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
- Fenics mesh corresponding to the domain.
- bdata : :class:`MeshFunctionSizet<dolfin.cpp.mesh.MeshFunction>`
- MeshFunction that labels the borders (i.e. outermost facets)
- of the mesh of the domain.
- cdata : :class:`MeshFunctionSizet<dolfin.cpp.mesh.MeshFunction>`
- MeshFunction that labels the inner elements (i.e. cells)
- of the mesh of the domain.
-
+ Plots the mesh.
+ For a list of possible arguments, see the documentation of the function `plot` in `bvpy.utils.visu`.
"""
- if fname is None:
- fname = self._dir_name + '/mesh.msh'
- cell_type = self._cell_type # we know the cell_type in advance
- elif os.path.isfile(fname):
- cell_type = None
- if not fname.endswith('.msh'):
- raise ValueError('The input filename should be a gmsh file with extension ".msh"')
- else:
- raise ValueError('The input filename should be a gmsh file with extension ".msh"')
-
- # - Read the gmsh file with meshio
- mesh_io = mio.read(fname, file_format='gmsh')
-
- # - Assert that the entities are allowed
- unknown_entites = set(mesh_io.cells_dict.keys()) - set(AVAILABLE_CELL_TYPE.keys())
-
- if unknown_entites:
- raise RuntimeError(f'Input mesh contains entities that are not handled : {unknown_entites}')
-
- # - Guess the dimension of the mesh if no model was provided
- if cell_type is None:
- if 'tetra' in mesh_io.cells_dict.keys():
- log.debug('Input mesh is a 3D tetra mesh')
- cell_type = 'tetra'
- elif 'triangle' in mesh_io.cells_dict.keys():
- log.debug('Input mesh is a 2D triangular mesh')
- cell_type = 'triangle'
- elif 'line' in mesh_io.cells_dict.keys():
- log.debug('Input mesh is a 1D line mesh')
- cell_type = 'line'
- else:
- raise RuntimeError(f'No entities were detected after the meshio conversion in the current input mesh!')
-
- ele_dim = AVAILABLE_CELL_TYPE[cell_type]
- else:
- ele_dim = AVAILABLE_CELL_TYPE[cell_type]
- element_type, element_tags, _ = self.model.mesh.getElements(ele_dim)
- n_entities = len(element_tags[np.where(element_type == ELEMENT_GMSH_TYPES[cell_type])[0][0]])
-
- if cell_type not in mesh_io.cells_dict:
- raise RuntimeError(f'None elements of the required topological dimension (element_dim = {ele_dim} for '
- f'{cell_type}) have been found during the conversion from Gmsh to MeshIo.'
- f'\nKnown issue 1: no Entities of the required topological dimension ({ele_dim}) '
- f'are defined in your Gmsh mesh. Use YourAbstractObject.factory.getEntities(dim={ele_dim}) '
- f'to check if it is the case.\nKnown issue 2: some PhysicalGroups have been set for '
- f'Entities of topological dimension < {ele_dim}. If PhysicalGroups are defined, be sure '
- f'to define them for all the Entities at topological dimension = {ele_dim}. If not, '
- f'all the unmarked Entities will be ignored.\nKnown issue 3: ???')
- elif len(mesh_io.cells_dict[cell_type]) != n_entities:
- raise Warning(f'Some elements have been lost during the conversion of the gmsh mesh, one possible '
- f'reason is the partial labelling of entities at topological dimension {ele_dim}. Make '
- f'sure that you label all the entities ! Unlabelled entities at dim {ele_dim} will be'
- f' ignored...')
-
- mesh_io.points = mesh_io.points[:, :self.dim]
-
- mesh = fe.Mesh()
-
- # - Extract physical properties (and possible name) for the domain
- if ('gmsh:physical' in mesh_io.cell_data_dict.keys()
- and cell_type in mesh_io.cell_data_dict["gmsh:physical"].keys()):
- log.debug(f"Found labels at the dimension {ele_dim}")
-
- # Notes:
- # 'gmsh: physical' lists the tags given to group of entities of same dimension. Here, the dimension will
- # be defined by ele_dim (from the _cell_type).
- # There tags can be defined by addPhysicalGroup(ele_dim, ele_tags, tagOfGroup) in gmsh.
- # It is possible to set a name for a given PhysicalGroup using setPhysicalName in gmsh.
- # Gmsh does not allowed intersection of PhysicalGroup, the first tag for a given entity
- # will be the one you will find here.
- # Ex: if elements [0, 1] are tagged with 4,
- # then 1 is tagged with 2,
- # you will see ele. 1 with tag 4 here.
- # Important warning : In the case where some entities have a PhysicialGroup, if it exists entities without
- # any PhysicialGroup, they will be ignored !
- # TO BE CLEAR : you will loose part of your mesh !
- # Conclusion : using PhysicalGroup, the resulting mesh will always contains PhysicialGroups that form a
- # partition of the mesh (at the given topological dimension)
- # However, it is possible that some of the PhysicialGroup do not have a PhysicalName, we will
- # deal with this situation by setting a default name.
-
- labels = mesh_io.cell_data_dict["gmsh:physical"][cell_type]
-
- # - Catch PhysicalName associated with PhysicalGroup
- # - reformat field_data as {tag : name} and filter PhysicalName (`ele_dim` topological dimension)
- if mesh_io.field_data:
- field_data = {lab_dim[0]: lab_name
- for lab_name, lab_dim in mesh_io.field_data.items()
- if lab_dim[1] == (ele_dim)} # tag: name
- log.debug(f"domains : field data => {field_data} found")
- else:
- field_data = {i: f'sub_domain_{i}' for i in np.unique(labels)}
- log.debug(f"domains : default field data => {field_data}")
-
- if set(labels) != set(field_data):
- log.warning('Domain : the number of PhysicalGroup does not correspond to the number of PhysicalName !')
- log.warning('Domain : some of your PhysicalGroup will be automatically named...')
-
- # - set automatic names for the other labels
- field_data.update({lab: f'sub_domain_{i}' for i, lab in enumerate(set(labels) - set(field_data))})
-
- mio.write(self._dir_name + '/' + 'mesh.xdmf',
- mio.Mesh(points=mesh_io.points,
- cells={cell_type: mesh_io.cells_dict[cell_type]},
- cell_data={f"{cell_type}_labels": [labels]}))
-
- with fe.XDMFFile(self._dir_name + '/' + 'mesh.xdmf') as f:
- f.read(mesh)
-
- mvc = fe.MeshValueCollection("size_t", mesh, ele_dim)
- with fe.XDMFFile(self._dir_name + '/' + 'mesh.xdmf') as infile:
- infile.read(mvc, f'{cell_type}_labels')
- cdata = fe.cpp.mesh.MeshFunctionSizet(mesh, mvc)
-
- self.mesh = mesh
- self.cdata = cdata
- self.sub_domain_names = field_data
-
- else:
- log.debug(f"No labels found at the dimension {ele_dim}, use default.")
-
- mio.write(self._dir_name + '/' + 'mesh.xdmf',
- mio.Mesh(points=mesh_io.points,
- cells={cell_type: mesh_io.cells_dict[cell_type]}))
-
- with fe.XDMFFile(self._dir_name + '/' + 'mesh.xdmf') as f:
- f.read(mesh)
-
- self.mesh = mesh
- self.cdata = _generate_default_cell_data(mesh)
- self.sub_domain_names = {0: 'domain'} # default naming for the domain
-
- # - Extract physical properties (and possible names) for the boundaries (ele_dim - 1)
- bnd_cell_type = [cell_type for cell_type, d in AVAILABLE_CELL_TYPE.items() if d == (ele_dim - 1)][0]
-
- if ('gmsh:physical' in mesh_io.cell_data_dict.keys()
- and bnd_cell_type in mesh_io.cell_data_dict["gmsh:physical"].keys()):
-
- labels = mesh_io.cell_data_dict["gmsh:physical"][bnd_cell_type]
-
- # - Catch PhysicalName associated with PhysicalGroup
- # - reformat field_data as {tag : name} and filter PhysicalName (`ele_dim` - 1 topological dimension)
- if mesh_io.field_data:
- field_data = {lab_dim[0]: lab_name
- for lab_name, lab_dim in mesh_io.field_data.items()
- if lab_dim[1] == (ele_dim - 1)} # tag: name
- log.debug(f"boundary : field data => {field_data} found")
- else:
- field_data = {i: f'sub_boundary_{i}' for i in np.unique(labels)}
- log.debug(f"boundary : default field data => {field_data}")
-
- if set(labels) != set(field_data):
- log.warning('Boundary: the number of PhysicalGroup does not correspond to the number of PhysicalName !')
- log.warning('Boundary: some of your PhysicalGroup will be automatically named...')
-
- # - set automatic names for the other labels
- field_data.update({lab: f'sub_boundary_{i}' for i, lab in enumerate(set(labels)-set(field_data))})
-
- bnd_mesh = mio.Mesh(points=mesh_io.points,
- cells={bnd_cell_type:
- mesh_io.cells_dict[bnd_cell_type]},
- cell_data={f"boundary": [labels]})
- mio.xdmf.write(self._dir_name + '/mesh_boundary.xdmf', bnd_mesh)
-
- mvc = fe.MeshValueCollection("size_t", mesh, ele_dim - 1)
- with fe.XDMFFile(self._dir_name + '/mesh_boundary.xdmf') as infile:
- infile.read(mvc, 'boundary')
- bdata = fe.cpp.mesh.MeshFunctionSizet(mesh, mvc)
-
- self.bdata = bdata
- self.sub_boundary_names = field_data
-
- else:
- self.bdata = _generate_default_boundary_data(mesh)
- self.sub_boundary_names = {0: 'boundary'}
-
- def info(self):
- """Prints the descrition of the domain.
-
- Returns
- -------
- None
-
- """
- print(self._details)
-
-
-class BuiltInModel(object):
- """A model based on the Gmsh `geo` factory.
-
- Attributes
- ----------
- factory : :class:`occ<gmsh.model.geo>`
- A gmsh model `geo` factory to discretize mesh from geometry.
-
- """
-
- def __init__(self, *args, **kwargs):
- """Sets the attribute factory to a specific type of gmsh factory.
-
- Returns
- -------
- factory : :class:`occ<gmsh.model.geo>`
- A gmsh model `geo` factory to discretize mesh from geometry.
-
- """
- super().__init__()
- self.factory = gmsh.model.geo
-
-
-class OccModel(object):
- """A model based on the Gmsh `Occ` factory.
-
- Attributes
- ----------
- factory : :class:`occ<gmsh.model.occ>`
- A gmsh model `occ` factory to discretize mesh from geometry.
-
- """
-
- def __init__(self, *args, **kwargs):
- """Sets the attribute factory to a specific type of gmsh factory.
-
- Returns
- -------
- factory : :class:`occ<gmsh.model.occ>`
- A gmsh model `occ` factory to discretize mesh from geometry.
-
- """
- super().__init__()
- self.factory = gmsh.model.occ
-
- def __sub__(self, other):
- """Cutting the left hand side Domain with the right hand side domain.
-
- Parameters
- ----------
- other : instance of subclass(AbstractDomain)
-
- Returns
- -------
- CSGDomain
- Returns th result of the cut operation
-
- """
-
- csg = CSGDomain(cell_type=self._cell_type,
- cell_size=self._cell_size,
- algorithm=self._algorithm,
- dim=self.dim)
-
- self.compute_new_size(other, csg)
-
- if not self.volumes:
- dim = 2
- objectdimtags = [(dim, tag) for tag in self.surfaces.values()]
- tooldimtags = [(dim, tag) for tag in other.surfaces.values()]
- outdimtag, _ = self.factory.cut(objectdimtags, tooldimtags)
- csg.surfaces[0] = outdimtag[0][1]
- else:
- dim = 3
- objectdimtags = [(dim, tag) for tag in self.volumes.values()]
- tooldimtags = [(dim, tag) for tag in other.volumes.values()]
- outdimtag, _ = self.factory.cut(objectdimtags, tooldimtags)
- csg.volumes[0] = outdimtag[0][1]
-
- self.factory.synchronize()
-
- return csg
-
- def __add__(self, other):
- """Adding the left hand side Domain with the right hand side domain.
-
- Parameters
- ----------
- other : instance of subclass(AbstractDomain)
-
- Returns
- -------
- CSGDomain
- Returns th result of the addition operation
-
- """
-
- csg = CSGDomain(cell_type=self._cell_type,
- cell_size=self._cell_size,
- algorithm=self._algorithm,
- dim=self.dim)
-
- self.compute_new_size(other, csg)
-
- if not self.volumes:
- dim = 2
- objectdimtags = [(dim, tag) for tag in self.surfaces.values()]
- tooldimtags = [(dim, tag) for tag in other.surfaces.values()]
- outdimtag, _ = self.factory.fuse(objectdimtags, tooldimtags)
- csg.surfaces[0] = outdimtag[0][1]
- else:
- dim = 3
- objectdimtags = [(dim, tag) for tag in self.volumes.values()]
- tooldimtags = [(dim, tag) for tag in other.volumes.values()]
- outdimtag, _ = self.factory.fuse(objectdimtags, tooldimtags)
- csg.volumes[0] = outdimtag[0][1]
-
- self.factory.synchronize()
-
- return csg
-
- def __and__(self, other):
- """Compute the intersection of the two domains.
-
- Parameters
- ----------
- other : instance of subclass(AbstractDomain)
-
- Returns
- -------
- CSGDomain
- Return the result of the intersection.
-
- """
-
- csg = CSGDomain(cell_type=self._cell_type,
- cell_size=self._cell_size,
- algorithm=self._algorithm,
- dim=self.dim)
-
- self.compute_new_size(other, csg)
-
- if not self.volumes:
- dim = 2
- objectdimtags = [(dim, tag) for tag in self.surfaces.values()]
- tooldimtags = [(dim, tag) for tag in other.surfaces.values()]
- outdimtag, _ = self.factory.intersect(objectdimtags, tooldimtags)
- csg.surfaces[0] = outdimtag[0][1]
- else:
- dim = 3
- objectdimtags = [(dim, tag) for tag in self.volumes.values()]
- tooldimtags = [(dim, tag) for tag in other.volumes.values()]
- outdimtag, _ = self.factory.intersect(objectdimtags, tooldimtags)
- csg.volumes[0] = outdimtag[0][1]
-
- self.factory.synchronize()
-
- return csg
-
- def compute_new_size(self, other, csg):
- """Computes a characteristic size for the new structure.
-
- Parameters
- ----------
- other : instance of subclass(AbstractDomain)
- The domain to combine with the considered one.
- csg : `CSGDomain<bvpy.domains.abstract.CSGDomain>`
- The combined domain that must be updated.
-
- Returns
- -------
- None
-
- """
-
- all_sizes = list(self._characteristic_size.values())\
- + list(other._characteristic_size.values())
- csg._characteristic_size = {'width': max(all_sizes)}
-
-
-class CSGDomain(AbstractDomain, OccModel):
- def __init__(self, *args, **kwargs):
- """The resulting class of an add/sub or intersection of Domains.
-
- """
- super(CSGDomain, self).__init__(**kwargs)
-
- def geometry(self, *args, **kwargs):
- pass
+ if not self.mesh:
+ raise ValueError("The mesh has not been generated yet.")
+ from bvpy.utils.visu import plot
+ plot(self, **kwargs)
diff --git a/src/bvpy/domains/custom_domain.py b/src/bvpy/domains/custom_domain.py
deleted file mode 100644
index 86eff0c..0000000
--- a/src/bvpy/domains/custom_domain.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-#
-# bvpy.templates.custom_polygonal
-#
-# File author(s):
-# Florian Gacon <florian.gacon@inria.fr>
-#
-# File contributor(s):
-# Florian Gacon <florian.gacon@inria.fr>
-#
-# File maintainer(s):
-# Olivier Ali <olivier.ali@inria.fr>
-#
-# Copyright © by Inria
-# Distributed under the LGPL License..
-# See accompanying file LICENSE.txt or copy at
-# https://www.gnu.org/licenses/lgpl-3.0.en.html
-#
-# -----------------------------------------------------------------------
-import gmsh
-import fenics as fe
-import meshio as mio
-import numpy as np
-
-from bvpy.domains.abstract import (AbstractDomain, BuiltInModel,
- _generate_default_cell_data,
- _generate_default_boundary_data)
-
-__all__ = ['CustomDomain']
-
-
-class CustomDomain(AbstractDomain, BuiltInModel):
- """A Domain subclass to describe custom surfaces.
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- """
-
- def __init__(self, points=[], cells=[], **kwargs):
- """Generates a cellularized domain structure.
-
- Parameters
- ----------
- path : str
- Location of the file to read to get the mesh.
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- """
-
- super(CustomDomain, self).__init__(**kwargs)
-
- self._characteristic_size = {}
-
- points = np.array(points)
- cells = np.array(cells)
- self.mesh = fe.Mesh()
-
- editor = fe.MeshEditor()
- if self._cell_type == 'line':
- top = 1
- elif self._cell_type == 'triangle':
- top = 2
- elif self._cell_type == 'tetra':
- top = 3
- editor.open(self.mesh, self._cell_type, top, self.dim)
- editor.init_vertices(len(points))
- editor.init_cells(len(cells))
-
- [editor.add_vertex(ind, point[:self.dim]) for ind, point in enumerate(points)]
- [editor.add_cell(ind, cell) for ind, cell in enumerate(cells)]
-
- editor.close()
-
- self.cdata = _generate_default_cell_data(self.mesh)
- self.bdata = _generate_default_boundary_data(self.mesh)
-
- self.ds = fe.Measure('ds', domain=self.mesh, subdomain_data=self.bdata)
- self.dx = fe.Measure('dx', domain=self.mesh, subdomain_data=self.cdata)
-
- def geometry(self):
- """Adds a custom structure to the Gmsh factory.
-
- WARNING: For custom domains generated from existing triangular mesh,
- The Gmsh step is not required and therefore the corresponding domains
- will not contain any `geometry`, the corresponding method is therefore
- useless. It remains in the code at a token of heritage.
-
- """
- pass
-
- @staticmethod
- def read(path):
- """Convert a mesh recorded on disk into a fenics one.
-
- Parameters
- ----------
- path : str
- The location of the mesh file to convert.
-
- Notes
- -----
- Gmsh is not involved in the domain production workflow.
-
- """
- mesh_io = mio.read(path)
- _, dim = mesh_io.points.shape
-
- points = mesh_io.points[:, :dim]
- cells = mesh_io.cells_dict["triangle"]
-
- return CustomDomain(points, cells)
diff --git a/src/bvpy/domains/custom_gmsh.py b/src/bvpy/domains/custom_gmsh.py
index 52808a7..ed88ef1 100644
--- a/src/bvpy/domains/custom_gmsh.py
+++ b/src/bvpy/domains/custom_gmsh.py
@@ -1,66 +1,917 @@
-# -*- python -*-
-# -*- coding: utf-8 -*-
-#
-# bvpy.templates.custom_gmsh
-#
-# file is an adaptation of customDomain class
-# File contributor(s)
-# Adrien Heymans <adrien.heymans@slu.se>
-# -----------------------------------------------------------------------
+import sys
-# -----------------------------------------------------------------------
-import fenics as fe
+from typing import Union, Tuple
+import numpy as np
+import numpy.linalg as lng
-from bvpy.domains.abstract import AbstractDomain, BuiltInModel
+import gmsh
+from bvpy.domains.abstract import AbstractDomain
+from bvpy import logger as log
+from bvpy import MPI_COMM_WORLD
+from bvpy.utils.gmshio import model_to_mesh
+from bvpy.utils.io import read_polygonal_data
-__all__ = ['CustomDomainGmsh']
+__all__ = ['CustomGMSH', 'OccModel', 'BuiltInModel']
+# Available meshing algorithm from gmsh (experimental algorithms are removed)
+# See https://gmsh.info/doc/texinfo/gmsh.html#Choosing-the-right-unstructured-algorithm
+# to have more information about what algo fits your need. Briefy :
+# - complex curved surfaces (2D) --> MeshAdapt
+# - high element quality (2D/3D) --> Frontal-Delaunay
+# - very large meshes --> Delaunay (2D) / HXT (3D)
+AVAILABLE_MESHING_ALGO_2D = {'MeshAdapt': 1, 'Automatic': 2, 'Delaunay': 5, 'Frontal-Delaunay': 6}
+AVAILABLE_MESHING_ALGO_3D = {'Delaunay': 1, 'Frontal-Delaunay': 4, 'HXT': 10}
+
+
+def _set_current_domain(func):
+ """Decorator to set the current domain in GMSH before calling any method."""
+ def wrapper(self, *args, **kwargs):
+ if not sys.is_finalizing() and gmsh.isInitialized() and (gmsh.model.getCurrent() != self.model_id):
+ gmsh.model.setCurrent(self.model_id)
+ return func(self, *args, **kwargs)
+
+ return wrapper
+
+
+class CustomGMSH(AbstractDomain):
+ """A Domain subclass to create customized meshes from GMSH geometries
+
+ Attributes
+ ----------
+ model_id : str
+ GMSH model name
+ cell_size : (float, float)
+ The cell_size min and max.
+ cell_size_from_curvature : int
+ GMSH Mesh.CharacteristicLengthFromCurvature option
+ algorithm : str
+ Name of the meshing algorithm to use.
+ verbosity : bool
+ GMSH verbosity
+ entities : list
+ GMSH Entities (tags) of tdim
+
+ Methods
+ -------
+ set_cell_size(cell_size)
+ Sets the characteristic size of the domain elements.
+ set_cell_size_from_curvature(n_elements_per_circle)
+ Sets the GMSH Mesh.CharacteristicLengthFromCurvature option
+ set_meshing_algorithm(algorithm)
+ Sets the algorithm to use to discretize the domain elements.
+ set_verbosity(verbosity)
+ Sets the GMSH verbosity.
+ synchronize()
+ synchronizes the CAD model with the gmsh model.
+ discretize()
+ Generates a fenics mesh from a GMSH geometry.
+ run_gui()
+ Opens GMSH GUI with the current model
+ rotate(point, axe, angle)
+ Rotates the GMSH model around a given point, axe and angle.
-class CustomDomainGmsh(AbstractDomain, BuiltInModel):
- """A Domain subclass to create bvpy meshes from existing file meshes
"""
+ def __init__(self,
+ cell_type: str,
+ gdim: int,
+ cell_size: Union[float, Tuple[float, float]],
+ verbosity: bool,
+ algorithm: str):
+ """
+ Parameters
+ ----------
+ cell_type : str
+ Type of element used to tile the domain.
+ It can be 'tetrahedron', 'tetra', 'triangle', 'interval' or 'line'.
+ gdim : int
+ Dimension of the embedded space, also called geometrical or physical dimension
+ cell_size : float, [float, float]
+ The cell_size value we want to use.
+ verbosity : bool
+ GMSH verbosity
+ algorithm : str
+ Name of the meshing algorithm to use.
+ Accepted meshing algorithms depends on the cell_type dimension:
+ * 2D ('triangle') : 'Delaunay', 'Frontal-Delaunay', 'MeshAdapt', 'Automatic'
+ * 3D ('tetra') : 'Delaunay', 'Frontal-Delaunay', 'HXT'
+ """
+
+ if not gmsh.isInitialized():
+ gmsh.initialize()
+
+ self.model_id = str(hex(id(self))) # Unique instance id based on the memory address
+ gmsh.model.add(self.model_id)
+
+ super().__init__(cell_type=cell_type, gdim=gdim)
+
+ # GMSH specs
+ self.__verbosity = None
+ self.set_verbosity(verbosity)
+ self.__cell_size = None
+ self.set_cell_size(cell_size)
+ self.__algorithm = None
+ self.set_meshing_algorithm(algorithm)
+ self.__cell_size_from_curvature = None
+
+ # Occ or BuiltIn factory
+ self.factory = None
+
+ # GMSH Entities (tags) of tdim
+ self.entities = []
+
+ @_set_current_domain
+ def __del__(self):
+ if not sys.is_finalizing() and gmsh.isInitialized():
+ gmsh.model.remove()
+
+ @property
+ def _details(self):
+ """A detailed printable description of the domain."""
+
+ description = super()._details
+
+ if self.mesh is None:
+ description += "Use self.discretize() to generate it.\n"
+
+# description += '''\nGMSH geometry specs
+# -------------------\n'''
+# if self._characteristic_size is not None:
+# for name, val in self._characteristic_size.items():
+# description += f' * {name}: {val:.2e}\n'
+
+ description += '\n' + f'''GMSH meshing specs
+------------------
+ * Algorithm: {self.algorithm}
+ * Resolution (cell size): min:{self.cell_size[0]:.2e}, max:{self.cell_size[1]:.2e}
+ * Verbosity: {self.verbosity}
+ '''
+
+ return description
+
+ @property
+ def cell_size(self):
+ """Cell size given to GMSH to compute the mesh.
+ It can be a float or a list of two floats (min/max)."""
+ return self.__cell_size
+
+ def set_cell_size(self, cell_size):
+ """Sets the characteristic size of the domain elements.
+
+ Parameters
+ ----------
+ cell_size : float, [float, float]
+ The cell_size value we want to use.
- def __init__(self, fname=None, **kwargs):
- """Initialize the custom domain.
+ """
+ if isinstance(cell_size, (float, int)):
+ self.__cell_size = (cell_size, cell_size)
+ elif (isinstance(cell_size, (list, tuple)) and (len(cell_size) == 2)
+ and all(isinstance(ele, (int, float)) for ele in cell_size)):
+ self.__cell_size = cell_size
+ else:
+ raise ValueError('cell_size must be a number or a list/tuple of two numbers')
+
+ @property
+ def cell_size_from_curvature(self):
+ """GMSH Mesh.CharacteristicLengthFromCurvature option"""
+ return self.__cell_size_from_curvature
+
+ def set_cell_size_from_curvature(self, n_elements_per_circle: int):
+ """Sets the GMSH Mesh.CharacteristicLengthFromCurvature option
+ """
+
+ if type(n_elements_per_circle) is not int:
+ raise ValueError('n_elements should be an integer')
+
+ gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", n_elements_per_circle)
+
+ self.__cell_size_from_curvature = n_elements_per_circle
+
+ @property
+ def algorithm(self):
+ """GMSH algorithm used to compute the mesh"""
+ return self.__algorithm
+
+ def set_meshing_algorithm(self, algorithm: str):
+ """Sets the algorithm to use to discretize the domain elements.
Parameters
----------
- fname : str
- Location of the file to read to get the mesh.
+ algorithm : str
+ Name of the meshing algorithm to use.
+ Accepted meshing algorithms depends on the cell_type dimension:
+ * 2D ('triangle') : 'Delaunay', 'Frontal-Delaunay', 'MeshAdapt', 'Automatic'
+ * 3D ('tetra') : 'Delaunay', 'Frontal-Delaunay', 'HXT'
- Other parameters
- ----------------
- dim : float
- Dimension of the embedding space to consider.
"""
- super(CustomDomainGmsh, self).__init__(**kwargs)
+ if self.tdim == 2 and algorithm not in AVAILABLE_MESHING_ALGO_2D.keys():
+ raise ValueError(f'Invalid meshing algorithm, can be {list(AVAILABLE_MESHING_ALGO_2D)}')
+
+ elif self.tdim == 3 and algorithm not in AVAILABLE_MESHING_ALGO_3D.keys():
+ raise ValueError(f'Invalid meshing algorithm, can be {list(AVAILABLE_MESHING_ALGO_3D)}')
+
+ elif self.tdim not in [2, 3]:
+ algorithm = None
+
+ self.__algorithm = algorithm
- if fname is not None:
- self.read_gmsh(fname)
+ @property
+ def verbosity(self):
+ """GMSH verbosity"""
+ return self.__verbosity
- def read_gmsh(self, fname):
- """Read a gmsh mesh recorded on disk.
+ def set_verbosity(self, verbosity: bool):
+ """Sets the GMSH verbosity.
Parameters
----------
- fname : str
- The location of the mesh file to convert.
+ verbosity : bool
+
"""
- if self.mesh is not None:
- self.model.mesh.clear()
+ if not isinstance(verbosity, bool):
+ raise ValueError('verbosity should be a boolean {True, False}')
+
+ self.__verbosity = verbosity
+
+ @_set_current_domain
+ def __set_gmsh_specs(self):
+ """Sets the GMSH specs before generating the mesh."""
+ # cell_size
+ gmsh.option.setNumber("Mesh.CharacteristicLengthMin", self.cell_size[0])
+ gmsh.option.setNumber("Mesh.CharacteristicLengthMax", self.cell_size[1])
- self.convert_gmsh_to_fenics(fname)
+ # cell_size_from_curvature
+ if self.__cell_size_from_curvature is not None:
+ gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", self.cell_size_from_curvature)
- self.ds = fe.Measure('ds', domain=self.mesh, subdomain_data=self.bdata)
- self.dx = fe.Measure('dx', domain=self.mesh, subdomain_data=self.cdata)
+ # algorithm
+ if self.tdim == 2:
+ gmsh.option.setNumber("Mesh.Algorithm", AVAILABLE_MESHING_ALGO_2D[self.algorithm])
+ elif self.tdim == 3:
+ gmsh.option.setNumber("Mesh.Algorithm", AVAILABLE_MESHING_ALGO_3D[self.algorithm])
- def geometry(self):
- """Adds a custom structure to the Gmsh factory.
+ # verbosity
+ if self.verbosity:
+ gmsh.option.setNumber("General.Verbosity", 5)
+ else:
+ gmsh.option.setNumber("General.Verbosity", 0)
- WARNING: For custom domains generated from existing triangular mesh,
- The Gmsh step is not required and therefore the corresponding domains
- will not contain any `geometry`, the corresponding method is therefore
- useless. It remains in the code at a token of heritage.
+ @_set_current_domain
+ def synchronize(self):
+ """synchronizes the CAD model with the gmsh model.
"""
- pass
\ No newline at end of file
+ self.factory.synchronize()
+
+ @_set_current_domain
+ def discretize(self):
+ """Generates a fenics mesh from a GMSH geometry.
+ Fills the mesh, cdata, bdata, sub_domain_names, sub_boundary_names, ds and dx attributes.
+
+ """
+ log.debug(f"Computing the mesh with topological element of dimension {self.tdim}"
+ f"and physical dimension of {self.gdim}")
+
+ # Adding missing PhysicalGroups for dimension tdim
+ c_ent = gmsh.model.getEntities(self.tdim) # total entities in the model
+ c_ent_pg = list({entity for group in gmsh.model.getPhysicalGroups(self.tdim) # entities already in a physical group
+ for entity in gmsh.model.getEntitiesForPhysicalGroup(group[0], group[1])})
+ if not c_ent:
+ raise ValueError(f'No entities have been defined for the topological dimension {self.tdim}.\n'
+ f'Please define them before generating the mesh.')
+ elif len(c_ent_pg) < len(c_ent):
+ log.info(f'PhysicalGroups are not defined for all the entities of topological dimension {self.tdim}.\n'
+ f'Defining the remaining entities PhysicalGroups as domain_automatically_generated.')
+ tags = [entity[1] for entity in c_ent if entity[1] not in c_ent_pg]
+ gmsh.model.addPhysicalGroup(self.tdim, tags, name='domain_automatically_generated')
+ # # Adding missing PhysicalGroups for dimension tdim - 1
+ # b_ent = gmsh.model.getEntities(self.tdim - 1) # total entities in the model
+ # b_ent_pg = list({entity for group in gmsh.model.getPhysicalGroups(self.tdim - 1) # entities already in a physical group
+ # for entity in gmsh.model.getEntitiesForPhysicalGroup(group[0], group[1])})
+ # if len(b_ent_pg) < len(b_ent) and len(b_ent_pg) > 0:
+ # log.info(f'PhysicalGroups are not defined for all the entities at topological dimension {self.tdim - 1}.\n'
+ # f'Defining the remaining entities PhysicalGroups as boundary_automatically_generated.')
+ # tags = [entity[1] for entity in b_ent if entity[1] not in b_ent_pg]
+ # gmsh.model.addPhysicalGroup(self.tdim - 1, tags, name='boundary_automatically_generated')
+
+ # Generate the mesh
+ self.__set_gmsh_specs()
+ gmsh.model.mesh.generate(self.tdim)
+ self.mesh, self.cdata, self.bdata, _, _, _ = model_to_mesh(gmsh.model, MPI_COMM_WORLD, MPI_COMM_WORLD.rank, gdim=self.gdim)
+
+ self.sub_domain_names = {
+ tag: gmsh.model.getPhysicalName(dim, tag) if gmsh.model.getPhysicalName(dim, tag) != ""
+ else f"sub_domain_{tag}"
+ for dim, tag in gmsh.model.getPhysicalGroups(self.tdim)
+ }
+ if self.bdata:
+ self.sub_boundary_names = {
+ tag: gmsh.model.getPhysicalName(dim, tag) if gmsh.model.getPhysicalName(dim, tag) != ""
+ else f"sub_boundary_{tag}"
+ for dim, tag in gmsh.model.getPhysicalGroups(self.tdim - 1)
+ }
+
+ self.set_mesh_functions()
+ gmsh.model.mesh.clear() # Clear the gmsh mesh to save memory
+
+ @staticmethod
+ @_set_current_domain
+ def run_gui():
+ """Opens GMSH GUI with the current model
+
+ Warnings
+ --------
+ The GMSH GUI should be previously installed on your machine.
+
+ """
+ gmsh.fltk.run()
+
+ @_set_current_domain
+ def rotate(self, point, axe, angle):
+ """Rotates the GMSH model around a given point, axe and angle.
+ """
+ dimTags = ([(self.tdim, val) for val in self.entities])
+ self.factory.rotate(dimTags, *point, *axe, angle)
+ self.synchronize()
+
+
+class OccModel(CustomGMSH):
+ """A subclass to create customized meshes from GMSH OpenCASCADE geometries
+
+ Attributes
+ ----------
+ factory : gmsh.model.occ
+ Factory to use to create the mesh
+
+ Methods
+ -------
+ add_rectangle(x=0, y=0, length=1, width=1)
+ Adds a rectangular domain to the existing GMSH model.
+ add_disk(center=(0,0), radius=1)
+ Adds a disk domain to the existing GMSH model.
+ add_cube(x=0, y=0, z=0, dx=1, dy=1, dz=1)
+ Adds a cube domain to the existing GMSH model.
+ add_cylinder(x=0, y=0, z=0, dx=0, dy=0, dz=1, r=0.5)
+ Adds a cylinder domain to the existing GMSH model.
+ add_sphere(center=(0,0,0), radius=1)
+ Adds a sphere domain to the existing GMSH model.
+ add_hemisphere(center=(0,0,0), radius=1)
+ Adds a hemisphere domain to the existing GMSH model.
+ add_ellipsoid(center=(0,0,0), radius_xy=1, radius_z=2)
+ Adds an ellipsoid domain to the existing GMSH model.
+ add_squeezed_ellipsoid(center=(0,0,0), radius_x=1, radius_y=1.5, radius_z=2)
+ Adds a squeezed ellipsoid domain to the existing GMSH model.
+ add_torus(center=(0,0,0), main_radius=2, tube_radius=1)
+ Adds a torus domain to the existing GMSH model.
+ fuse(tag1, tag2)
+ Fuses two GMSH entities
+ intersect(tag1, tag2)
+ Intersects two GMSH entities
+ cut(tag1, tag2)
+ Cuts two GMSH entities: Removes the part of the first entity that is inside the second entity
+
+ """
+
+ def __init__(self,
+ cell_type: str,
+ gdim: int,
+ cell_size: Union[float, Tuple[float, float]] = 1.0,
+ verbosity: bool = False,
+ algorithm: str = 'Delaunay'):
+ """
+ Parameters
+ ----------
+ cell_type : str
+ Type of element used to tile the domain.
+ It can be 'tetrahedron', 'tetra', 'triangle', 'interval' or 'line'.
+ gdim : int
+ Dimension of the embedded space, also called geometrical or physical dimension
+ cell_size : float, [float, float]
+ The cell_size value we want to use.
+ verbosity : bool
+ GMSH verbosity
+ algorithm : str
+ Name of the meshing algorithm to use.
+ Accepted meshing algorithms depends on the cell_type dimension:
+ * 2D ('triangle') : 'Delaunay', 'Frontal-Delaunay', 'MeshAdapt', 'Automatic'
+ * 3D ('tetra') : 'Delaunay', 'Frontal-Delaunay', 'HXT'
+ """
+
+ super().__init__(cell_type=cell_type,
+ gdim=gdim,
+ cell_size=cell_size,
+ verbosity=verbosity,
+ algorithm=algorithm)
+
+ self.factory = gmsh.model.occ
+
+ @_set_current_domain
+ def add_rectangle(self, x=0, y=0, length=1, width=1):
+ """
+ Adds a rectangular domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ x : float
+ Horizontal position of the lower left corner.
+ y : float
+ Vertical position of the lower left corner.
+ length : float
+ Horizontal extension of the rectangle.
+ width : float
+ Vertical extension of the rectangle.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added rectangle
+
+ """
+ assert self.gdim == 2, 'The rectangle can only be added when gdim is 2'
+
+ result_entity = gmsh.model.occ.addRectangle(x, y, 0, length, width)
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_disk(self, center=(0,0), radius=1):
+ """
+ Adds a disk domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ center : 2-tuple
+ The center of the disk.
+ radius : float
+ The radius of the disk.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added disk
+
+ """
+ assert self.gdim == 2, 'The disk can only be added when gdim is 2'
+
+ result_entity = gmsh.model.occ.addDisk(center[0], center[1], 0, radius, radius)
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_cube(self, x=0, y=0, z=0, dx=1, dy=1, dz=1):
+ """
+ Adds a cube domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ x : float
+ Horizontal position of the lower left corner.
+ y : float
+ Vertical position of the lower left corner.
+ z : float
+ Depth position of the lower left corner.
+ dx : float
+ Horizontal extension of the cube.
+ dy : float
+ Vertical extension of the cube.
+ dz : float
+ Depth extension of the cube.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added cube
+
+ """
+ assert self.gdim == 3, 'The cube can only be added when gdim is 3'
+
+ result_entity = gmsh.model.occ.addBox(x, y, z, dx, dy, dz)
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_cylinder(self, x=0, y=0, z=0, dx=0, dy=0, dz=1, r=0.5):
+ """
+ Adds a cylinder domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ x : float
+ Horizontal position of the lower left corner.
+ y : float
+ Vertical position of the lower left corner.
+ z : float
+ Depth position of the lower left corner.
+ dx : float
+ Horizontal extension of the cylinder.
+ dy : float
+ Vertical extension of the cylinder.
+ dz : float
+ Depth extension of the cylinder.
+ r : float
+ Radius of the cylinder.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added cylinder
+
+ Notes
+ -----
+ If the topological dimension is 2, the top and bottom faces are removed.
+
+ """
+ assert self.gdim == 3, 'The cylinder can only be added when gdim is 3'
+
+ result_entity = gmsh.model.occ.addCylinder(x, y, z, dx, dy, dz, r)
+ self.synchronize()
+ if self.tdim == 2: # If the cylinder is 2D, we remove the top and bottom faces
+ bnd = gmsh.model.getBoundary([(3, result_entity)], oriented=False)
+ gmsh.model.removeEntities([(3, result_entity)])
+ gmsh.model.removeEntities([bnd[1], bnd[2]])
+ result_entity = bnd[0][1]
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_sphere(self, center=(0,0,0), radius=1):
+ """
+ Adds a sphere domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ center : tuple
+ The position vector of the center of the sphere.
+ radius : float
+ The radius of the sphere.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added sphere
+
+ """
+ assert self.gdim == 3, 'The sphere can only be added when gdim is 3'
+
+ result_entity = gmsh.model.occ.addSphere(center[0], center[1], center[2], radius)
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_hemisphere(self, center=(0,0,0), radius=1):
+ """
+ Adds a hemisphere domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ center : tuple
+ The position vector of the center of the complete sphere.
+ radius : float
+ The radius of the hemisphere.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added hemisphere
+
+ Notes
+ -----
+ If the topological dimension is 2, the bottom face is removed.
+
+ """
+ assert self.gdim == 3, 'The hemisphere can only be added when gdim is 3'
+
+ result_entity = gmsh.model.occ.addSphere(center[0], center[1], center[2], radius, angle1=0)
+ self.synchronize()
+ if self.tdim == 2: # If the hemisphere is 2D, we remove the disk
+ surfaces = gmsh.model.getBoundary([(3, result_entity)], oriented=False)
+ gmsh.model.occ.remove([(3, result_entity), surfaces[1]], recursive=False)
+ self.synchronize()
+ result_entity = surfaces[0][1]
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_ellipsoid(self, center=(0,0,0), radius_xy=1, radius_z=2):
+ """
+ Adds an ellipsoid domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ center : tuple
+ Position vector of the center of the sphere.
+ radius_xy : float
+ Radius of the ellipsoid within the Oxy plane.
+ radius_z : float
+ Radius of the ellipsoid along the Oz axis.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added ellipsoid
+
+ """
+ assert self.gdim == 3, 'The ellipsoid can only be added when gdim is 3'
+
+ ratio = radius_xy/radius_z
+ result_entity = gmsh.model.occ.addSphere(center[0], center[1], center[2], radius_z)
+ gmsh.model.occ.dilate([(3, result_entity)], *center, ratio, ratio, 1)
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_squeezed_ellipsoid(self, center=(0,0,0), radius_x=1, radius_y=1.5, radius_z=2):
+ """
+ Adds a squeezed ellipsoid domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ center : tuple
+ Position vector of the center of the sphere.
+ radius_x : float
+ Radius of the ellipsoid along the Ox.
+ radius_y : float
+ Radius of the ellipsoid along the Oy.
+ radius_z : float
+ Radius of the ellipsoid along the Oz axis.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added ellipsoid
+
+ """
+ assert self.gdim == 3, 'The ellipsoid can only be added when gdim is 3'
+
+ ratio_xy = radius_y/radius_x
+ ratio_xz = radius_z/radius_x
+ result_entity = gmsh.model.occ.addSphere(center[0], center[1], center[2], radius_x)
+ gmsh.model.occ.dilate([(3, result_entity)], *center, 1, ratio_xy, ratio_xz)
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_torus(self, center=(0,0,0), main_radius=2, tube_radius=1):
+ """
+ Adds a torus domain to the existing GMSH model.
+
+ Parameters
+ ----------
+ center : tuple
+ The position vector of the center of the torus, i.e. the center of the main circle.
+ main_radius : float
+ The radius of circle around which the tube is wrapped.
+ tube_radius : float
+ The tube radius.
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the added torus
+
+ """
+ assert self.gdim == 3, 'The torus can only be added when gdim is 3'
+
+ result_entity = gmsh.model.occ.addTorus(center[0], center[1], center[2], main_radius, tube_radius)
+ self.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def add_conformal_polygon(self, points=None, cells=None, path=None, markers=None, isoriented=True):
+ if path and points:
+ raise ValueError("You cannot provide both a file path and points.")
+ elif not path and not points:
+ raise ValueError("You must provide either a file path or (points and cells).")
+ elif path:
+ points, cells, labels = read_polygonal_data(path)
+
+ assert self.gdim == len(points[0]), 'Incorrect dimension of the points'
+ assert self.cell_type == 'triangle', 'The cell type must be triangle' # TODO: Adapt this to other cell types
+
+ # TODO: This method should be optimized and better written
+
+ cell_widths = [_max_width([pts for vid, pts in enumerate(points)
+ if vid in cell]) for cell in cells]
+ avg_cell_width = np.mean(cell_widths)
+
+ if not isoriented:
+ cells = _orient_cell_contours(points, cells)
+
+ # Add the geometry to the GMSH model
+ pts = np.asarray(points)
+ cells = np.asarray(cells, dtype=object)
+
+ points = {id: gmsh.model.occ.addPoint(*point, tag=id)
+ for id, point in enumerate(pts)}
+
+ [gmsh.model.occ.addSurfaceFilling(
+ gmsh.model.occ.addCurveLoop([gmsh.model.occ.addLine(points[c[ind_p]],
+ points[c[(ind_p + 1) % len(c)]])
+ for ind_p, p in enumerate(c)]), tag=ind_c+1)
+ for ind_c, c in enumerate(cells)]
+
+ gmsh.model.occ.removeAllDuplicates()
+
+ self.synchronize()
+
+ if markers is None:
+ [gmsh.model.addPhysicalGroup(2, [ind+1], ind)
+ for ind in range(len(cells))]
+ else:
+ markers = np.asarray(markers)
+ indexes = np.unique(markers)
+ for i in indexes:
+ tags = np.argwhere(markers==i).flatten()+1
+ gmsh.model.addPhysicalGroup(2, tags, i)
+
+ @_set_current_domain
+ def fuse(self, tag1, tag2):
+ """
+ Fuses two GMSH entities
+
+ Parameters
+ ----------
+ tag1 : int
+ A GMSH entity of topological dimension tdim
+ tag2 : int
+ A GMSH entity of topological dimension tdim
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the result of the fusion
+
+ """
+ result_entity, _ = gmsh.model.occ.fuse([(self.tdim, tag1)], [(self.tdim, tag2)])
+ gmsh.model.occ.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def intersect(self, tag1, tag2):
+ """
+ Intersects two GMSH entities
+
+ Parameters
+ ----------
+ tag1 : int
+ A GMSH entity of topological dimension tdim
+ tag2 : int
+ A GMSH entity of topological dimension tdim
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the result of the intersection
+
+ """
+ result_entity, _ = gmsh.model.occ.intersect([(self.tdim, tag1)], [(self.tdim, tag2)])
+ gmsh.model.occ.synchronize()
+
+ return result_entity
+
+ @_set_current_domain
+ def cut(self, tag1, tag2):
+ """
+ Cuts two GMSH entities: Removes the part of the first entity that is inside the second entity
+
+ Parameters
+ ----------
+ tag1 : int
+ A GMSH entity of topological dimension tdim
+ tag2 : int
+ A GMSH entity of topological dimension tdim
+
+ Returns
+ -------
+ int
+ The GMSH entity tag of the result of the cut
+
+ """
+ result_entity, _ = gmsh.model.occ.cut([(self.tdim, tag1)], [(self.tdim, tag2)])
+ gmsh.model.occ.synchronize()
+
+ return result_entity
+
+
+class BuiltInModel(CustomGMSH):
+ """A subclass to create customized meshes from GMSH built-in geometries
+
+ Attributes
+ ----------
+ factory : gmsh.model.geo
+ Factory to use to create the mesh
+
+ """
+
+ def __init__(self,
+ cell_type: str,
+ gdim: int,
+ cell_size: Union[float, Tuple[float, float]] = 1.0,
+ verbosity: bool = False,
+ algorithm: str = 'Delaunay'):
+ """
+ Parameters
+ ----------
+ cell_type : str
+ Type of element used to tile the domain.
+ It can be 'tetrahedron', 'tetra', 'triangle', 'interval' or 'line'.
+ gdim : int
+ Dimension of the embedded space, also called geometrical or physical dimension
+ cell_size : float, [float, float]
+ The cell_size value we want to use.
+ verbosity : bool
+ GMSH verbosity
+ algorithm : str
+ Name of the meshing algorithm to use.
+ Accepted meshing algorithms depends on the cell_type dimension:
+ * 2D ('triangle') : 'Delaunay', 'Frontal-Delaunay', 'MeshAdapt', 'Automatic'
+ * 3D ('tetra') : 'Delaunay', 'Frontal-Delaunay', 'HXT'
+ """
+ super().__init__(cell_type=cell_type,
+ gdim=gdim,
+ cell_size=cell_size,
+ verbosity=verbosity,
+ algorithm=algorithm)
+
+ self.factory = gmsh.model.geo
+
+
+def _orient_cell_contours(points, cells):
+ """Arranges the vertex labels surrounding cells in a coherent order.
+ To be used in the add_conformal_polygon method
+
+ Parameters
+ ----------
+ points : list of list of float
+ The list of the position vectors of the vertices.
+ cells : list of list of int
+ The lists of vertex indices surrounding each cell.
+
+ Returns
+ -------
+ oriented_cells : list of list of int
+ An oriented version of the `cells` argument: Within each list
+ the vertices labels are arranged in a continuous manner.
+
+ """
+ # TODO: Adapt properly this algo !!!
+ oriented_cells = []
+ for i, cell in enumerate(cells):
+ vrtcs = [points[vid] for vid in cell]
+ vrtcs = [np.array(v) - np.array(vrtcs).mean(axis=0) for v in vrtcs]
+
+ o_cell, o_vrtcs = [cell.pop(0)], [vrtcs.pop(0)]
+ while len(cell) > 0:
+ if len(o_cell) == 1:
+ orient = np.zeros(len(vrtcs))
+ else:
+ last_nrml = np.cross(o_vrtcs[-2], o_vrtcs[-1])
+ nrmls = [np.cross(o_vrtcs[-1], v) for v in vrtcs]
+ orient = [np.sign(np.dot(n, last_nrml)) for n in nrmls]
+
+ good_vrtcs = [(k, v) for k, (v, o)
+ in enumerate(zip(vrtcs, orient)) if o >= 0]
+
+ dist = [np.dot(v / lng.norm(v),
+ o_vrtcs[-1] / lng.norm(o_vrtcs[-1]))
+ for _, v in good_vrtcs]
+
+ vid_2_keep = good_vrtcs[dist.index(max(dist))][0]
+
+ o_cell.append(cell.pop(vid_2_keep))
+ o_vrtcs.append(vrtcs.pop(vid_2_keep))
+
+ oriented_cells.append(o_cell)
+
+ return oriented_cells
+
+
+def _max_width(points):
+ """Computes the maximal distance to the center in a point cloud.
+ To be used in the add_conformal_polygon method
+
+ Parameters
+ ----------
+ points : list of float
+ List of the position vectors of the vertices.
+
+ Returns
+ -------
+ float
+ Twice the maximal distance between any point within the considered
+ list of points and the center of this list.
+ """
+ center = np.array(points).mean(axis=0)
+ max_width = 2 * max([np.linalg.norm(pos - center)
+ for pos in np.array(points)])
+ return max_width
diff --git a/src/bvpy/domains/custom_polygonal.py b/src/bvpy/domains/custom_polygonal.py
deleted file mode 100644
index ea52993..0000000
--- a/src/bvpy/domains/custom_polygonal.py
+++ /dev/null
@@ -1,246 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-#
-# bvpy.templates.custom_polygonal
-#
-# File author(s):
-# Florian Gacon <florian.gacon@inria.fr>
-#
-# File contributor(s):
-# Florian Gacon <florian.gacon@inria.fr>
-#
-# File maintainer(s):
-# Olivier Ali <olivier.ali@inria.fr>
-#
-# Copyright © by Inria
-# Distributed under the LGPL License..
-# See accompanying file LICENSE.txt or copy at
-# https://www.gnu.org/licenses/lgpl-3.0.en.html
-#
-# -----------------------------------------------------------------------
-import numpy as np
-import numpy.linalg as lng
-
-from bvpy.domains.abstract import AbstractDomain, OccModel
-from bvpy.utils.io import read_polygonal_data
-
-__all__ = ['CustomPolygonalDomain']
-
-
-class CustomPolygonalDomain(AbstractDomain, OccModel):
- """A Domain subclass to describe cellularized strucutres.
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- Notes
- -----
- We choose to describe the characteristic size of the structure here with
- three numbers: its maximal width, the number of polygons it contains and
- the average width of these polygons.
- """
-
- def __init__(self, points=[], cells=[],
- markers=None, isoriented=True, **kwargs):
- """Generates a cellularized domain structure.
-
- Parameters
- ----------
- points : list of list of float
- List of the position vectors of the vertices (the default is []).
- cells : list of list of int
- Oriented lists of vertex indices
- surrounding each cell (the default is []).
- markers : list of int
- label values given to each cell (the default is None).
- isoriented : bool
- if True this means that the lists of vertex indices describing the
- cells are all sorted in a consistent order (either clockwise
- or anti-clockwize). If False, a routine in applied
- to orient them properly.
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- ------
- :class:`CustomPolygonalDomain<bvpy.templates.domains.custom_polygonal.CustomPolygonalDomain>`
- A custom piecewise-polygonal domain instance.
-
- """
- super(CustomPolygonalDomain, self).__init__(**kwargs)
-
- cell_widths = [self.max_width([pts for vid, pts in enumerate(points)
- if vid in cell]) for cell in cells]
- avg_cell_width = np.mean(cell_widths)
-
- self._characteristic_size = {'Maximal width': self.max_width(points),
- 'Number of polygons': len(cells),
- 'Averaged polygon width': avg_cell_width}
-
- if not isoriented:
- cells = self.orient_cell_contours(points, cells)
-
- self.geometry(points, cells, markers)
-
- def max_width(self, points):
- """Computes the maximal distance to the center in a point cloud.
-
- Parameters
- ----------
- points : list of float
- List of the position vectors of the vertices.
-
- Returns
- -------
- float
- Twice the maximal distance between any point within the considered
- list of points and the center of this list.
- """
- center = np.array(points).mean(axis=0)
- max_width = 2 * max([np.linalg.norm(pos - center)
- for pos in np.array(points)])
- return max_width
-
- def geometry(self, points, cells, markers):
- """Adds a custom piecewise-polygonal structure to the Gmsh factory.
-
- Parameters
- ----------
- points : list of list of float
- List of the position vectors of the vertices.
- cells : list of list of int
- Oriented lists of vertex indices
- surrounding each cell.
- markers : list of int
- label values given to each cell.
-
- Returns
- -------
- None
-
- """
- pts = np.asarray(points)
- cells = np.asarray(cells, dtype=object)
-
- self.points = {id: self.factory.addPoint(*point, tag=id)
- for id, point in enumerate(pts)}
-
- [self.factory.addSurfaceFilling(
- self.factory.addCurveLoop([self.factory.addLine(self.points[c[ind_p]],
- self.points[c[(ind_p + 1) % len(c)]])
- for ind_p, p in enumerate(c)]), tag=ind_c+1)
- for ind_c, c in enumerate(cells)]
-
- self.factory.removeAllDuplicates()
-
- self.synchronize()
-
- if markers is None:
- [self.model.addPhysicalGroup(2, [ind+1], ind)
- for ind in range(len(cells))]
- else:
- markers = np.asarray(markers)
- indexes = np.unique(markers)
- for i in indexes:
- tags = np.argwhere(markers==i).flatten()+1
- self.model.addPhysicalGroup(2, tags, i)
-
-
- def orient_cell_contours(self, points, cells):
- """Arranges the vertex labels surrounding cells in a coherent order.
-
- Parameters
- ----------
- points : list of list of float
- The list of the position vectors of the vertices.
- cells : list of list of int
- The lists of vertex indices surrounding each cell.
-
- Returns
- -------
- oriented_cells : list of list of int
- An oriented version of the `cells` argument: Within each list
- the vertices labels are arranged in a continuous manner.
-
- """
- # TODO: Adapt properly this algo !!!
- oriented_cells = []
- for i, cell in enumerate(cells):
- vrtcs = [points[vid] for vid in cell]
- vrtcs = [np.array(v) - np.array(vrtcs).mean(axis=0) for v in vrtcs]
-
- o_cell, o_vrtcs = [cell.pop(0)], [vrtcs.pop(0)]
- while len(cell) > 0:
- if len(o_cell) == 1:
- orient = np.zeros(len(vrtcs))
- else:
- last_nrml = np.cross(o_vrtcs[-2], o_vrtcs[-1])
- nrmls = [np.cross(o_vrtcs[-1], v) for v in vrtcs]
- orient = [np.sign(np.dot(n, last_nrml)) for n in nrmls]
-
- good_vrtcs = [(k, v) for k, (v, o)
- in enumerate(zip(vrtcs, orient)) if o >= 0]
-
- dist = [np.dot(v / lng.norm(v),
- o_vrtcs[-1] / lng.norm(o_vrtcs[-1]))
- for _, v in good_vrtcs]
-
- vid_2_keep = good_vrtcs[dist.index(max(dist))][0]
-
- o_cell.append(cell.pop(vid_2_keep))
- o_vrtcs.append(vrtcs.pop(vid_2_keep))
-
- oriented_cells.append(o_cell)
-
- return oriented_cells
-
- @staticmethod
- def read(path, cell_size=.1, isoriented=True, clear=True, dim=3):
- """Instantiates a cellularized domain from a mesh file.
-
- Parameters
- ----------
- path : str
- location of the file to read.
- isoriented : bool optional
- if false, the vertex orientation around each cell is recomputed
- (the default is True).
- clear : bool
- Optional. If truen, clear all the previously defined Gmsh
- geometries. The default is False.
- dim : int
- Optional. Dimension of the problem. The default is 3.
-
- Returns
- -------
- :class: `CustomPolygonalDomain<bvpy.templates.domains.custom_polygonal.CustomPolygonalDomain>`
- A cellularized do5main.
-
- """
- points, cells, labels = read_polygonal_data(path)
- return CustomPolygonalDomain(points, cells,
- cell_size=cell_size,
- markers=labels,
- isoriented=isoriented,
- clear=clear,
- dim=dim)
diff --git a/src/bvpy/domains/fixed_domain.py b/src/bvpy/domains/fixed_domain.py
new file mode 100644
index 0000000..e3d78b3
--- /dev/null
+++ b/src/bvpy/domains/fixed_domain.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# bvpy.templates.custom_polygonal
+#
+# File author(s):
+# Florian Gacon <florian.gacon@inria.fr>
+#
+# File contributor(s):
+# Florian Gacon <florian.gacon@inria.fr>
+# Gonzalo Revilla <gonzalo.revilla-mut@inria.fr>
+#
+# File maintainer(s):
+# Olivier Ali <olivier.ali@inria.fr>
+#
+# Copyright © by Inria
+# Distributed under the LGPL License..
+# See accompanying file LICENSE.txt or copy at
+# https://www.gnu.org/licenses/lgpl-3.0.en.html
+#
+# -----------------------------------------------------------------------
+import fenics as fe
+import meshio as mio
+import numpy as np
+
+from bvpy.domains.abstract import AbstractDomain
+
+__all__ = ['FixedDomain']
+
+
+class FixedDomain(AbstractDomain):
+ """
+ Subclass to describe domains from meshes provided by the user.
+ They can be either provided as a file path or as points and cells.
+
+ """
+
+ def __init__(self, points=[], cells=[], path=None):
+ """
+ Parameters
+ ----------
+ points : list of list of float
+ List of the position vectors of the vertices.
+ cells : list of list of int
+ Lists of vertex indices surrounding each cell.
+ path : str
+ location of the file to read.
+
+ Warnings
+ --------
+ - Currently, the only cell type supported is the triangle.
+ - If path is given, the mesh must be in a format readable by meshio.
+
+ """
+ # TODO : Make this function generic for all mesh types. Currently only admits triangular meshes
+
+ if path and points:
+ raise ValueError("You cannot provide both a file path and points.")
+ elif not path and len(points)==0:
+ raise ValueError("You must provide either a file path or (points and cells).")
+ elif path:
+ # Read the mesh
+ mesh_io = mio.read(path)
+ points = mesh_io.points
+ cells = mesh_io.cells_dict["triangle"]
+
+ gdim = len(points[0])
+ points = np.array(points)
+ cells = np.array(cells)
+
+ super().__init__(cell_type="triangle", gdim=gdim)
+
+ self.mesh = fe.Mesh()
+ editor = fe.MeshEditor()
+ editor.open(self.mesh, self.cell_type, 2, self.gdim)
+ editor.init_vertices(len(points))
+ editor.init_cells(len(cells))
+ [editor.add_vertex(ind, point[:self.gdim]) for ind, point in enumerate(points)]
+ [editor.add_cell(ind, cell) for ind, cell in enumerate(cells)]
+ editor.close()
+
+ self.set_mesh_functions()
diff --git a/src/bvpy/domains/fixed_gmsh.py b/src/bvpy/domains/fixed_gmsh.py
new file mode 100644
index 0000000..0a5ff3d
--- /dev/null
+++ b/src/bvpy/domains/fixed_gmsh.py
@@ -0,0 +1,58 @@
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+# bvpy.templates.custom_gmsh
+#
+# file is an adaptation of customDomain class
+# File contributor(s)
+# Adrien Heymans <adrien.heymans@slu.se>
+# Gonzalo Revilla <gonzalo.revilla-mut@inria.fr>
+# -----------------------------------------------------------------------
+
+# -----------------------------------------------------------------------
+import gmsh
+from bvpy import MPI_COMM_WORLD
+from bvpy.domains.abstract import AbstractDomain
+from bvpy.utils.gmshio import model_to_mesh
+
+__all__ = ['FixedGMSH']
+
+
+class FixedGMSH(AbstractDomain):
+ """A Domain subclass to instanciate a bvpy domain loading a GMSH mesh from a file.
+ """
+
+ def __init__(self, fname, gdim):
+ """
+ Parameters
+ ----------
+ fname : str
+ Location of the file to read to get the mesh with .msh extension.
+ gdim : int
+ Dimension of the embedded space, also called geometrical or physical dimension.
+
+ Warnings
+ --------
+ - The mesh must have at least one physical group defined on all the cells.
+ """
+ if not fname.endswith('.msh'):
+ raise ValueError('The input filename should be a gmsh file with extension ".msh"')
+
+ # TODO: Mesh should have physicalGroup
+ if not gmsh.isInitialized():
+ gmsh.initialize()
+ model_id = str(hex(id(self)))
+ gmsh.model.add(model_id)
+ gmsh.merge(str(fname))
+ mesh, cdata, bdata, _, _, _ = model_to_mesh(gmsh.model, MPI_COMM_WORLD, MPI_COMM_WORLD.rank, gdim=gdim)
+ gmsh.model.remove()
+
+ cell_type = mesh.ufl_cell().cellname()
+ gdim = mesh.geometry().dim()
+ super().__init__(cell_type=cell_type, gdim=gdim)
+
+ self.mesh = mesh
+ self.cdata = cdata
+ self.bdata = bdata
+
+ self.set_mesh_functions()
diff --git a/src/bvpy/domains/geometry.py b/src/bvpy/domains/geometry.py
index df9c99c..e5cef89 100644
--- a/src/bvpy/domains/geometry.py
+++ b/src/bvpy/domains/geometry.py
@@ -257,34 +257,23 @@ def boundary_normal(domain, scale=1, name='boundary_normal'):
"""
if isinstance(domain, fe.Mesh):
mesh = domain
- dim = len(domain.coordinates()[0])
elif issubclass(type(domain), AbstractDomain):
mesh = domain.mesh
- dim = domain.dim
else:
raise ValueError("domain must be a mesh or a subclass of AbstractDomain")
+ gdim = mesh.geometry().dim()
+ tdim = mesh.topology().dim()
mesh.init()
bnd_facets = get_boundary_facets(mesh)
facets_index = [f.index() for f in bnd_facets]
- facet_normals = np.zeros((mesh.num_entities(1), 3))
- normals = np.zeros((mesh.num_entities(0), 3))
-
- #for fct in bnd_facets:
- # cell = [c for c in fe.cells(fct)][0]
- # vertex = [v for v in fe.vertices(fct)][0]
- # normal_cell = cell.cell_normal().array()
- # midpoint_facet = fct.midpoint().array()
- # midpoint_vertex = vertex.midpoint().array()
- # vec = midpoint_vertex-midpoint_facet
- # vec = vec/np.linalg.norm(vec)
- # n = np.cross(normal_cell, vec)
- # facet_normals[fct.index()] = n
+ facet_normals = np.zeros((mesh.num_entities(tdim - 1), gdim))
+ normals = np.zeros((mesh.num_entities(0), gdim))
for fct in bnd_facets:
- if fct.dim() == dim - 1: # Manifold without curvature
- facet_normals[fct.index()] = fct.normal().array()
- elif fct.dim() == 1 and dim == 3: # 2-manifold embedded in a 3-dimension space
+ if gdim == tdim: # Manifold without curvature
+ facet_normals[fct.index()] = fct.normal().array()[:gdim]
+ elif tdim == 2 and gdim == 3:
cell = [c for c in fe.cells(fct)][0]
vertex = [v for v in fe.vertices(fct)][0]
normal_cell = cell.cell_normal().array()
diff --git a/src/bvpy/domains/primitives.py b/src/bvpy/domains/primitives.py
deleted file mode 100644
index a2ae31d..0000000
--- a/src/bvpy/domains/primitives.py
+++ /dev/null
@@ -1,894 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-#
-# bvpy.domains.primitives
-#
-# File author(s):
-# Florian Gacon <florian.gacon@inria.fr>
-#
-# File contributor(s):
-# Florian Gacon <florian.gacon@inria.fr>
-# Olivier Ali <olivier.ali@inria.fr>
-#
-# File maintainer(s):
-# Olivier Ali <olivier.ali@inria.fr>
-#
-# Copyright © by Inria
-# Distributed under the LGPL License..
-# See accompanying file LICENSE.txt or copy at
-# https://www.gnu.org/licenses/lgpl-3.0.en.html
-#
-# -----------------------------------------------------------------------
-
-from bvpy.domains.abstract import AbstractDomain, BuiltInModel, OccModel
-
-import numpy as np
-
-
-class Rectangle(AbstractDomain, OccModel):
- """Instanciates a rectangular domain.
-
- Parameters
- ----------
- x : float
- Horizontal position of the lower left corder (the default is 0).
- y : float
- Vertical position of the lower left corder (the default is 0).
- length : float
- Horizontal extention of the rectangle (the default is 1).
- width : float
- Vertical extention of the rectangle (the default is 1).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- Methods
- -------
- geometry :
- Adds a rectangle element to the Gmsh factory
- in charge of computing the domain.
- discretize :
- Generates a fenics mesh representing the domain.
- See the doc of the mother class for more info:
- :class:`AbstractDomain<bvpy.domains.abstract.AbstractDomain>`
-
- Notes
- -----
- We choose to describe the characteristic size of the structure here with
- two numbers: its length and width parameters.
- """
-
- def __init__(self, x=0, y=0, length=1, width=1, **kwargs):
- """Generates an intance of a rectangular domain.
-
- Parameters
- ----------
- x : float
- Horizontal position of the lower left corder (the default is 0).
- y : float
- Vertical position of the lower left corder (the default is 0).
- length : float
- Horizontal extention of the rectangle (the default is 1).
- width : float
- Vertical extention of the rectangle (the default is 1).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- -------
- :class:`Rectangle<bvpy.templates.domains.primitives.Rectangle>`
- A rectangular domain instance.
-
- """
- super(Rectangle, self).__init__(**kwargs)
- self._characteristic_size = {'Length': length,
- 'Width': width}
- self.geometry(x, y, length, width)
-
- def geometry(self, x, y, length, width):
- """Adds a rectangle to the Gmsh factory computing the domain.
-
- Parameters
- ----------
- x : float
- Horizontal position of the lower left corder (the default is 0).
- y : float
- Vertical position of the lower left corder (the default is 0).
- length : float
- Horizontal extention of the rectangle (the default is 1).
- width : float
- Vertical extention of the rectangle (the default is 1).
-
- Returns
- -------
- None
-
- """
-
- self.surfaces[0] = self.factory.addRectangle(x, y, 0, length, width)
- self.synchronize()
-
-
-class Disk(AbstractDomain, OccModel):
- """Instantiates a circular domain.
-
- Parameters
- ----------
- radius : float
- The radius of the disk (the default is 1).
- angle_x : float
- Rotation angle around the Ox axis (the default is 0).
- angle_y : float
- Rotation angle around the Oy axis (the default is 0).
- angle_z : float
- Rotation angle around the Oz axis (the default is 0).
-
- Other Parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain (the default is 'triangle').
- Accepted values are 'line' (1D element), 'triangle' (2D), or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider (the default is 3).
-
- Attributes
- ----------
- _characteristic_size : dict
- A dictionary describing the characteristic sizes of the structure.
-
- - **keys** (str): The names of characteristics (e.g., "length", "width", "radius").
- - **values** (float): The corresponding size value for each characteristic.
-
- Methods
- -------
- geometry
- Adds a rectangle element to the Gmsh factory in charge of computing the domain.
- discretize
- Generates a fenics mesh representing the domain. See the documentation
- of the parent class :class:`AbstractDomain <bvpy.domains.abstract.AbstractDomain>`.
-
- Notes
- -----
- The characteristic size of the structure is described here with a single parameter, its radius.
- """
-
- def __init__(self, center=[0, 0, 0],
- radius=1, angle_x=0, angle_y=0, angle_z=0, **kwargs):
- """Generates a disk instance.
-
- Parameters
- ----------
- radius : float
- The radius of the disk (the default is 1).
- angle_x : float
- Rotation angle around the Ox axis (the default is 0).
- angle_y : float
- Rotation angle around the Oy axis (the default is 0).
- angle_z : float
- Rotation angle around the Oz axis (the default is 0).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- ------
- :class:`Rectangle<bvpy.templates.domains.primitives.Disk>`
- A disk domain instance.
-
- """
-
- super(Disk, self).__init__(**kwargs)
- self._characteristic_size = {'Radius': radius}
- self.geometry(center, radius)
-
- def geometry(self, center, radius):
- """Adds a disk to the Gmsh factory computing the domain.
-
- Parameters
- ----------
- radius : float
- The radius of the disk (the default is 1).
- angle_x : float
- Rotation angle around the Ox axis (the default is 0).
- angle_y : float
- Rotation angle around the Oy axis (the default is 0).
- angle_z : float
- Rotation angle around the Oz axis (the default is 0).
-
- Returns
- -------
- None
-
- """
-
- circle = self.factory.addCircle(*center, radius)
- curve_loop = self.factory.addCurveLoop([circle])
- self.surfaces[0] = self.factory.addPlaneSurface([curve_loop])
- self.synchronize()
-
-
-class HemiSphere(AbstractDomain, OccModel):
- """Instanciates an hemispherical domain.
-
- Parameters
- ----------
- radius : float
- The radius of the hemisphere (the default is 1).
- center : list of floats
- The position vector of the center of the sphere,
- the hemisphere is extracted from (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- Methods
- -------
- geometry :
- Adds a rectangle element to the Gmsh factory
- in charge of computing the domain.
- discretize :
- Generates a fenics mesh representing the domain.
- See the doc of the mother class for more info:
- :class:`AbstractDomain<bvpy.domains.abstract.AbstractDomain>`
-
- Notes
- -----
- We choose to describe the characteristic size of the structure here with
- one number: its radius parameter.
-
- """
-
- def __init__(self, radius=1, center=[0, 0, 0], **kwargs):
- """Generates an instance of hemisphere domain.
-
- Parameters
- ----------
- radius : float
- The radius of the hemisphere (the default is 1).
- center : list of floats
- The position vector of the center of the sphere,
- the hemisphere is extracted from (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- ------
- :class:`HemiSphere<bvpy.templates.domains.primitives.HemiSphere>`
- An hemispherical domain instance.
-
- """
- super(HemiSphere, self).__init__(**kwargs)
- self._characteristic_size = {'Radius': radius}
- self.geometry(radius, center)
-
- def geometry(self, radius, center):
- """Adds an hemisphere to the Gmsh factory computing the domain.
-
- Parameters
- ----------
- radius : float
- The radius of the hemisphere (the default is 1).
- center : list of floats
- The position vector of the center of the sphere,
- the hemisphere is extracted from (the default is [0, 0, 0]).
-
- Returns
- -------
- None
-
- """
- self.volumes[0] = self.factory.addSphere(*center, radius, angle1=0)
- self.synchronize()
-
- if self._cell_type == 'triangle': # Case where tdim=2 and gdim=3
- # In this case we just want the exterior boundary
- surfaces = self.model.getBoundary([(3, self.volumes[0])], oriented=False)
- self.factory.remove([(3, self.volumes[0]), surfaces[1]], recursive=False)
- self.synchronize()
-
-
-class Cube(AbstractDomain, OccModel):
- def __init__(self, x=0, y=0, z=0, dx=1, dy=1, dz=1, **kwargs):
- super(Cube, self).__init__(**kwargs)
- self._characteristic_size = {'Length': max(dx, dy, dz),
- 'Width': min(dx, dy, dz)}
- self.geometry(x, y, z, dx, dy, dz)
-
- def geometry(self, x, y, z, dx, dy, dz):
- self.volumes[0] = self.factory.addBox(x, y, z, dx, dy, dz)
- self.factory.synchronize()
-
-
-class Cylinder(AbstractDomain, OccModel):
- def __init__(self, x=0, y=0, z=0,
- dx=0, dy=0, dz=1, r=0.5, open=False, **kwargs):
- super(Cylinder, self).__init__(**kwargs)
- self._characteristic_size = {'Length': max(dx, dy, dz),
- 'Width': min(dx, dy, dz)}
- self.geometry(x, y, z, dx, dy, dz, r, open)
-
- def geometry(self, x, y, z, dx, dy, dz, r, open):
- self.volumes[0] = self.factory.addCylinder(x, y, z, dx, dy, dz, r)
- self.factory.synchronize()
-
- if open:
- self._cell_type = "triangle"
- bnd = self.model.getBoundary([(3, self.volumes[0])])
- self.model.removeEntities([(3, self.volumes[0])])
- self.model.removeEntities([bnd[1], bnd[2]])
-
-
-class Sphere(AbstractDomain, OccModel):
- """Instanciates an spherical domain.
-
- Parameters
- ----------
- radius : float
- The radius of the sphere (the default is 1).
- center : list of floats
- The position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- Methods
- -------
- geometry :
- Adds a rectangle element to the Gmsh factory
- in charge of computing the domain.
- discretize :
- Generates a fenics mesh representing the domain.
- See the doc of the mother class for more info:
- :class:`AbstractDomain<bvpy.domains.abstract.AbstractDomain>`
-
- Notes
- -----
- We choose to describe the characteristic size of the structure here with
- one number: its radius parameter.
-
- """
-
- def __init__(self, radius=1, center=[0, 0, 0], **kwargs):
- """Generates an instance of spherical domain.
-
- Parameters
- ----------
- radius : float
- The radius of the sphere (the default is 1).
- center : list of floats
- The position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- ------
- :class:`Sphere<bvpy.templates.domains.primitives.Sphere>`
- A spherical domain instance.
-
- """
- super(Sphere, self).__init__(**kwargs)
- self._characteristic_size = {'Radius': radius}
- self.geometry(radius, center)
-
- def geometry(self, radius, center):
- """Adds a sphere to the Gmsh factory computing the domain.
-
- Parameters
- ----------
- radius : float
- The radius of the sphere (the default is 1).
- center : list of floats
- The position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Returns
- -------
- None
-
- """
- self.volumes[0] = self.factory.addSphere(*center, radius)
- self.synchronize()
-
-
-class Ellipsoid(AbstractDomain, OccModel):
- """Instanciates an Ellispoid domain.
-
- Parameters
- ----------
- radius_xy : float
- The radius of the ellipsoid within the Oxy plane (the default is 1).
- radius_z : float
- The radius of the ellipsoid along the Oz axis (the default is 2).
- center : list of floats
- The position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- Methods
- -------
- geometry :
- Adds an ellipsoid element to the Gmsh factory
- in charge of computing the domain.
- discretize :
- Generates a fenics mesh representing the domain.
- See the doc of the mother class for more info:
- :class:`AbstractDomain<bvpy.domains.abstract.AbstractDomain>`
-
- Notes
- -----
- We choose to describe the characteristic size of the structure here with
- two numbers: its two radius parameters.
-
- """
-
- def __init__(self, radius_xy=1, radius_z=2, center=[0, 0, 0], **kwargs):
- """Generates an instance of ellipsoid domain.
-
- Parameters
- ----------
- radius_xy : float
- Radius of the ellipsoid within the Oxy plane (the default is 1).
- radius_z : float
- Radius of the ellipsoid along the Oz axis (the default is 2).
- center : list of floats
- Position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- ------
- :class:`Ellipsoid<bvpy.templates.domains.primitives.Ellipsoid>`
- An ellipsoid domain instance.
-
- """
- super(Ellipsoid, self).__init__(**kwargs)
- self._characteristic_size = {'Oxy radius': radius_xy,
- 'Oz radius': radius_z}
- self.geometry(radius_xy, radius_z, center)
-
- def geometry(self, radius_xy, radius_z, center=[0, 0, 0]):
- """Adds an ellipsoid to the Gmsh factory computing the domain.
-
- Parameters
- ----------
- radius_xy : float
- Radius of the ellipsoid within the Oxy plane (the default is 1).
- radius_z : float
- Radius of the ellipsoid along the Oz axis (the default is 2).
- center : list of floats
- Position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Returns
- -------
- None
-
- """
- assert len(center) == 3
- ratio = radius_xy/radius_z
- self.volumes[0] = self.factory.addSphere(*center, radius_z)
- self.factory.dilate([(3, self.volumes[0])], *center, ratio, ratio, 1)
- self.synchronize()
-
-class Squeezed_Ellipsoid(AbstractDomain, OccModel):
- """Instanciates an Squeezed Ellispoid domain.
-
- Parameters
- ----------
- radius_xy : float
- The radius of the ellipsoid within the Oxy plane (the default is 1).
- radius_z : float
- The radius of the ellipsoid along the Oz axis (the default is 2).
- center : list of floats
- The position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- Methods
- -------
- geometry :
- Adds an ellipsoid element to the Gmsh factory
- in charge of computing the domain.
- discretize :
- Generates a fenics mesh representing the domain.
- See the doc of the mother class for more info:
- :class:`AbstractDomain<bvpy.domains.abstract.AbstractDomain>`
-
- Notes
- -----
- We choose to describe the characteristic size of the structure here with
- two numbers: its two radius parameters.
-
- """
-
- def __init__(self, radius_x=1, radius_y=1.5, radius_z=2,
- center=[0, 0, 0], **kwargs):
- """Generates an instance of ellipsoid domain.
-
- Parameters
- ----------
- radius_x : float
- Radius of the ellipsoid along the Ox (the default is 1).
- radius_y : float
- Radius of the ellipsoid along the Oy (the default is 1.5).
- radius_z : float
- Radius of the ellipsoid along the Oz axis (the default is 2).
- center : list of floats
- Position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- ------
- :class:`Ellipsoid<bvpy.templates.domains.primitives.Ellipsoid>`
- An ellipsoid domain instance.
-
- """
- super(Squeezed_Ellipsoid, self).__init__(**kwargs)
- self._characteristic_size = {'Ox radius': radius_x,
- 'Oy radius': radius_y,
- 'Oz radius': radius_z}
- self.geometry(radius_x, radius_y, radius_z, center)
-
- def geometry(self, radius_x, radius_y, radius_z, center=[0, 0, 0]):
- """Adds an ellipsoid to the Gmsh factory computing the domain.
-
- Parameters
- ----------
- radius_x : float
- Radius of the ellipsoid within along the Ox axis (default is 1).
- radius_y : float
- Radius of the ellipsoid within along the Oy axis (default is 1.5).
- radius_z : float
- Radius of the ellipsoid along the Oz axis (default is 2).
- center : list of floats
- Position vector of the center of the sphere
- (the default is [0, 0, 0]).
-
- Returns
- -------
- None
-
- """
- assert len(center) == 3
- ratio_xy = radius_y/radius_x
- ratio_xz = radius_z/radius_x
- self.volumes[0] = self.factory.addSphere(*center, radius_x)
- self.factory.dilate([(3, self.volumes[0])], *center, 1,
- ratio_xy, ratio_xz)
- self.synchronize()
-
-
-class Torus(AbstractDomain, OccModel):
- """Instanciates a toroid domain.
-
- Parameters
- ----------
- main_radius : float
- The radius of circle around which the tube is wrapped
- (the default is 2).
- tube_radius : float
- The tube radius (the default is 1).
- center : list of floats
- The position vector of the center of the torus,
- i.e. the center of the main circle (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Attributes
- ----------
- _characteristic_size : dict
- - keys : str
- The name of the characteristics (i.e. length, width, radius...)
- - values : float
- The corresponding value.
-
- Methods
- -------
- geometry :
- Adds an ellipsoid element to the Gmsh factory
- in charge of computing the domain.
- discretize :
- Generates a fenics mesh representing the domain.
- See the doc of the mother class for more info:
- :class:`AbstractDomain<bvpy.domains.abstract.AbstractDomain>`
-
- Notes
- -----
- We choose to describe the characteristic size of the structure here with
- two numbers: the main_radius and the tube_radius parameters.
-
- """
-
- def __init__(self, main_radius=2, tube_radius=1, center=[0, 0, 0],
- **kwargs):
- """Generates an instance of toroidal domain.
-
- Parameters
- ----------
- main_radius : float
- The radius of circle around which the tube is wrapped
- (the default is 2).
- tube_radius : float
- The tube radius (the default is 1).
- center : list of floats
- The position vector of the center of the torus,
- i.e. the center of the main circle (the default is [0, 0, 0]).
-
- Other parameters
- ----------------
- cell_type : str
- Type of element used to tile the domain
- (the default is 'triangle'). The accepted values can be:
- 'line' (1D element), 'triangle' (2D) or 'tetra' (3D).
- cell_size : float
- Characteristic size of the tiling elements (the default is 1).
- verbosity : bool
- If `True`, log is printed (the default is False).
- algorithm : str
- Name of the algorithm used to compute the elements
- (the default is 'Delaunay').
- dimension : float
- Dimension of the embedding space to consider.
- (the default is 3)
-
- Yields
- ------
- :class:`Torus<bvpy.templates.domains.primitives.Torus>`
- A toroidal domain instance.
-
- """
- super(Torus, self).__init__(**kwargs)
- self._characteristic_size = {'Main radius': main_radius,
- 'Tube radius': tube_radius}
- self.geometry(center, main_radius, tube_radius)
-
- def geometry(self, center, main_radius, tube_radius):
- """Adds a torus to the Gmsh factory computing the domain.
-
- Parameters
- ----------
- center : list of floats
- The position vector of the center of the torus,
- i.e. the center of the main circle (the default is [0, 0, 0]).
- main_radius : float
- The radius of circle around which the tube is wrapped
- (the default is 2).
- tube_radius : float
- The tube radius (the default is 1).
-
- Returns
- -------
- None
-
- """
- self.factory.addTorus(*center, main_radius, tube_radius)
- self.synchronize()
diff --git a/src/bvpy/utils/examples.py b/src/bvpy/utils/examples.py
index 1bb1b3a..e4352d0 100644
--- a/src/bvpy/utils/examples.py
+++ b/src/bvpy/utils/examples.py
@@ -19,7 +19,7 @@
#
# -----------------------------------------------------------------------
-from bvpy.domains import Cylinder
+from bvpy.domains import OccModel
from bvpy.boundary_conditions import dirichlet
from bvpy.vforms.elasticity import LinearElasticForm
from bvpy.utils.visu import plot
@@ -84,9 +84,9 @@ def cantilever_beam(visualizer: str='pyvista'):
assert visualizer in ['pyvista', 'plotly'], "Acceptable values are 'pyvista' (default) or 'plotly'."
# - Create a cylinder and set up some parameters for the mesh construction
- cylinder = Cylinder(x=0, dx=2, dz=0, r=0.2, cell_type='tetra',
- cell_size=.1, algorithm='Frontal-Delaunay', clear=True)
- cylinder.discretize() # discretize (from Gmsh -> FeniCS)
+ domain = OccModel(cell_type='tetra', gdim=3, cell_size=.1, algorithm='Frontal-Delaunay', clear=True)
+ domain.add_cylinder(x=0, dx=2, dz=0, r=0.2)
+ domain.discretize() # discretize (from Gmsh -> FeniCS)
# - Create boundary conditions : clamp the beam at x=0
bc = [dirichlet(val=[0, 0, 0], boundary="near(x, 0, 1e-2)")]
@@ -99,7 +99,7 @@ def cantilever_beam(visualizer: str='pyvista'):
source=[0, 0, -rho * g])
# - Create the boundary-value problem and solve it
- bvp = BVP(domain=cylinder, vform=vform, bc=bc)
+ bvp = BVP(domain=domain, vform=vform, bc=bc)
bvp.solve()
# - Visualize the resulting displacement field
diff --git a/test/data/mesh_example.msh b/test/data/mesh_example.msh
index 076883c..ecdd631 100644
--- a/test/data/mesh_example.msh
+++ b/test/data/mesh_example.msh
@@ -11,10 +11,10 @@ $Entities
2 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 2 2 -3
3 -9.999999994736442e-08 0.9999999000000001 -1e-07 1.0000001 1.0000001 1e-07 0 2 3 -4
4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 1.0000001 1e-07 0 2 4 -1
-1 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 4 1 2 3 4
+1 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 1 1 4 1 2 3 4
$EndEntities
$Nodes
-9 143 1 143
+9 98 1 98
0 1 0 1
1
0 0 0
@@ -27,7 +27,7 @@ $Nodes
0 4 0 1
4
0 1 0
-1 1 0 9
+1 1 0 7
5
6
7
@@ -35,57 +35,59 @@ $Nodes
9
10
11
+0.1249999999999998 0 0
+0.2499999999999998 0 0
+0.3749999999999994 0 0
+0.4999999999999988 0 0
+0.6249999999999988 0 0
+0.7499999999999991 0 0
+0.8749999999999994 0 0
+1 2 0 7
12
13
-0.1 0 0
-0.2 0 0
-0.3 0 0
-0.4 0 0
-0.5 0 0
-0.6 0 0
-0.7 0 0
-0.8 0 0
-0.9 0 0
-1 2 0 9
14
15
16
17
18
+1 0.1249999999999998 0
+1 0.2499999999999998 0
+1 0.3749999999999994 0
+1 0.4999999999999988 0
+1 0.6249999999999988 0
+1 0.7499999999999991 0
+1 0.8749999999999994 0
+1 3 0 7
19
20
21
22
-1 0.1 0
-1 0.2 0
-1 0.3 0
-1 0.4 0
-1 0.5 0
-1 0.6 0
-1 0.7 0
-1 0.8 0
-1 0.9 0
-1 3 0 9
23
24
25
+0.8750000000000001 1 0
+0.7500000000000002 1 0
+0.6250000000000007 1 0
+0.5000000000000011 1 0
+0.3750000000000012 1 0
+0.2500000000000009 1 0
+0.1250000000000006 1 0
+1 4 0 7
26
27
28
29
30
31
-0.9 1 0
-0.8 1 0
-0.7 1 0
-0.6 1 0
-0.5 1 0
-0.4 1 0
-0.3 1 0
-0.2 1 0
-0.09999999999999998 1 0
-1 4 0 9
32
+0 0.8750000000000001 0
+0 0.7500000000000002 0
+0 0.6250000000000007 0
+0 0.5000000000000011 0
+0 0.3750000000000012 0
+0 0.2500000000000009 0
+0 0.1250000000000006 0
+2 1 0 66
33
34
35
@@ -94,16 +96,6 @@ $Nodes
38
39
40
-0 0.9 0
-0 0.8 0
-0 0.7 0
-0 0.6 0
-0 0.5 0
-0 0.4 0
-0 0.3 0
-0 0.2 0
-0 0.09999999999999998 0
-2 1 0 103
41
42
43
@@ -162,452 +154,236 @@ $Nodes
96
97
98
-99
-100
-101
-102
-103
-104
-105
-106
-107
-108
-109
-110
-111
-112
-113
-114
-115
-116
-117
-118
-119
-120
-121
-122
-123
-124
-125
-126
-127
-128
-129
-130
-131
-132
-133
-134
-135
-136
-137
-138
-139
-140
-141
-142
-143
-0.5 0.5 0
-0.7071428571428572 0.2928571428571428 0
-0.7071428571428572 0.7071428571428572 0
-0.2813605605204175 0.701969500996966 0
-0.2928571428571428 0.2928571428571428 0
-0.7758579074798024 0.5147174832368676 0
-0.493149962973225 0.7759084652631476 0
-0.506850037026775 0.2240915347368524 0
-0.2240915347368524 0.493149962973225 0
-0.8294973544973545 0.1705026455026455 0
-0.8294973544973545 0.8294973544973545 0
-0.1705026455026455 0.1705026455026454 0
-0.1705026455026454 0.8294973544973545 0
-0.3592898873606312 0.1511339309159103 0
-0.6407947253692821 0.1518248097738525 0
-0.8488660690840897 0.6407101126393689 0
-0.3607435781083035 0.8484499785903453 0
-0.1511339309159103 0.6407101126393688 0
-0.1518248097738525 0.359205274630718 0
-0.8488660690840897 0.3592898873606312 0
-0.6407101126393688 0.8488660690840897 0
-0.6401567944250871 0.4330139372822299 0
-0.5669860627177701 0.6401567944250871 0
-0.3598432055749129 0.5669860627177701 0
-0.4330139372822299 0.3598432055749129 0
-0.6733856581720996 0.5770645560214576 0
-0.4489409534615607 0.6663521796462133 0
-0.333826635273307 0.4413012504750118 0
-0.5586987495249882 0.333826635273307 0
-0.5335396291359238 0.109843205574913 0
-0.109843205574913 0.4664603708640763 0
-0.4664603708640763 0.890156794425087 0
-0.890156794425087 0.5335396291359238 0
-0.7403047640904488 0.09944092595747056 0
-0.9023348642131718 0.2609263968561355 0
-0.9006094296444203 0.7402985969818545 0
-0.7390736031438645 0.9023348642131718 0
-0.09939057035557974 0.7402985969818544 0
-0.09944092595747056 0.2596952359095513 0
-0.2609263968561355 0.09766513578682823 0
-0.2610747326467143 0.9022924059995243 0
-0.745282971718187 0.3958793750760887 0
-0.6041206249239113 0.745282971718187 0
-0.3958793750760887 0.2547170282818129 0
-0.2547170282818129 0.6041206249239113 0
-0.384704891613151 0.7505210099667106 0
-0.615295108386849 0.2494789900332894 0
-0.2494789900332894 0.384704891613151 0
-0.5978296212305664 0.5245716233160158 0
-0.9059804230793693 0.09401957692063068 0
-0.09401957692063068 0.9059804230793693 0
-0.09688906079890419 0.09688906079890419 0
-0.9031109392010959 0.9031109392010959 0
-0.4299212442790752 0.08324867532408453 0
-0.9167513246759155 0.4299212442790752 0
-0.08324867532408453 0.5700787557209248 0
-0.5700787557209248 0.9167513246759155 0
-0.2003834335631949 0.7370236452034065 0
-0.7953061572217176 0.725001607673958 0
-0.7223833510181306 0.2060517816310846 0
-0.2060517816310846 0.2776166489818695 0
-0.2722061756621068 0.8012652966219919 0
-0.8012652966219918 0.2722061756621067 0
-0.7277938243378933 0.8012652966219918 0
-0.2722061756621067 0.1987347033780082 0
-0.7523637681935527 0.6266755117910972 0
-0.542335390056664 0.421352693305455 0
-0.4501658790294238 0.5758506912721821 0
-0.4239643555425228 0.4574794944307918 0
-0.3575006807940414 0.6519524273377049 0
-0.6921239612473346 0.4945204197407769 0
-0.6424993192059586 0.3480475726622952 0
-0.3480475726622952 0.3575006807940414 0
-0.3500275056991684 0.06640954840536456 0
-0.9335904515946354 0.3500275056991684 0
-0.07222178826124714 0.3570721762808691 0
-0.6429278237191309 0.07222178826124714 0
-0.3576557363238189 0.9281798358029913 0
-0.06675463531911485 0.6502174930684296 0
-0.9279264586307194 0.6429096677514294 0
-0.6499724943008316 0.9335904515946354 0
-0.07997140572242911 0.8201017201097869 0
-0.8238071187089842 0.92184034909451 0
-0.1761928812910159 0.07815965090548999 0
-0.9200285942775709 0.1798982798902132 0
-0.9200285942775709 0.8201017201097869 0
-0.1798982798902131 0.9200285942775709 0
-0.8201017201097869 0.07997140572242911 0
-0.07997140572242911 0.1798982798902131 0
-0.519231323232856 0.7031172999754856 0
-0.4736105247275205 0.2931196009667213 0
-0.6516259541526819 0.6560383883487352 0
-0.2931196009667213 0.5263894752724796 0
-0.1646068749667147 0.5549039654243012 0
-0.5549039654243012 0.8353931250332853 0
-0.4450960345756987 0.1646068749667147 0
-0.8313685985821415 0.4403319049504258 0
-0.1828291668772978 0.4276562509996817 0
-0.4279309480492506 0.8172199024734518 0
-0.5723437490003183 0.1828291668772978 0
-0.8176142401408163 0.5757657555868232 0
-0.5224615528514942 0.5758647740932955 0
-0.2194965262590453 0.6679809534620259 0
+0.4375000000000012 0.8917468245269454 0
+0.5624999999999989 0.1082531754730547 0
+0.8787408590474761 0.5688788296333867 0
+0.1107249391757802 0.4359051368937709 0
+0.3098918535636471 0.1070150426096389 0
+0.8917468245269455 0.3124999999999996 0
+0.6871891340151978 0.8955068139135626 0
+0.1182626851017918 0.6692259739230109 0
+0.1094851921887065 0.2033420655976694 0
+0.8119703288398898 0.1058465056319792 0
+0.1875000000000007 0.8917468245269449 0
+0.879230463236312 0.8130906423764597 0
+0.687411721473314 0.1078520638328755 0
+0.6222382678779044 0.2148535068271494 0
+0.4995397113129832 0.2162308769262824 0
+0.5592159774975001 0.3228521478837691 0
+0.4350399505186866 0.3232496269229517 0
+0.4995115585047765 0.4327221170255846 0
+0.6130571419252585 0.4327008212914639 0
+0.5569896182521401 0.5431506844806554 0
+0.4365001961261519 0.5415315810733969 0
+0.4989149690630473 0.6498774713089112 0
+0.6204622660663317 0.6518781417959358 0
+0.3748733700916653 0.6495569611339232 0
+0.3123122610363024 0.541316479365892 0
+0.5624999999999974 0.757772228311381 0
+0.249999999999999 0.6495190528383259 0
+0.3133176815799774 0.7523710733550599 0
+0.685881360403847 0.7503445567222475 0
+0.7558138451570435 0.6422638761822286 0
+0.1875710532781787 0.5449975168824133 0
+0.2499999999999996 0.4330127018922171 0
+0.1778732193721211 0.3232140480521861 0
+0.3014840143016392 0.3178829819971146 0
+0.2375681941678718 0.2131232466900194 0
+0.7646672502010368 0.229066006106814 0
+0.891746824526945 0.6874999999999991 0
+0.7774226411965255 0.5265527835888375 0
+0.8628769461960453 0.4382017528180521 0
+0.7819074654780986 0.3451910717748128 0
+0.3125000000000009 0.891746824526945 0
+0.4369885941461045 0.1080009076591809 0
+0.5625000000000008 0.8917468245269453 0
+0.371489561667968 0.2148949910804583 0
+0.3724746634145927 0.4316192480461928 0
+0.4374999999999981 0.7577722283113807 0
+0.8959004857786097 0.1879321779033213 0
+0.8081462788186125 0.9070743513265296 0
+0.1034157830572141 0.8066379931605242 0
+0.1911333561284868 0.1028269135594174 0
+0.7923273839860845 0.74371081557924 0
+0.2010889394410992 0.7499266504670594 0
+0.6739266620600561 0.3223988573944038 0
+0.08453570123751114 0.5625000000000009 0
+0.671840423222428 0.5422756641907901 0
+0.07733755452579338 0.3125000000000011 0
+0.09150635094610993 0.09150635094611009 0
+0.90849364905389 0.09150635094610993 0
+0.09150635094611009 0.9084936490538901 0
+0.9084936490538901 0.90849364905389 0
+0.7301718800130687 0.43455349184306 0
+0.3750000000000013 0.8247595264191645 0
+0.5000000000000009 0.8247595264191647 0
+0.6250000000000007 0.8247595264191652 0
+0.2601608324359393 0.8246932292153027 0
+0.7489644266062752 0.8253550026152656 0
$EndNodes
$Elements
-9 288 1 288
-0 1 15 1
-1 1
-0 2 15 1
-2 2
-0 3 15 1
-3 3
-0 4 15 1
-4 4
-1 1 1 10
-5 1 5
-6 5 6
-7 6 7
-8 7 8
-9 8 9
-10 9 10
-11 10 11
-12 11 12
-13 12 13
-14 13 2
-1 2 1 10
-15 2 14
-16 14 15
-17 15 16
-18 16 17
-19 17 18
-20 18 19
-21 19 20
-22 20 21
-23 21 22
-24 22 3
-1 3 1 10
-25 3 23
-26 23 24
-27 24 25
-28 25 26
-29 26 27
-30 27 28
-31 28 29
-32 29 30
-33 30 31
-34 31 4
-1 4 1 10
-35 4 32
-36 32 33
-37 33 34
-38 34 35
-39 35 36
-40 36 37
-41 37 38
-42 38 39
-43 39 40
-44 40 1
-2 1 2 244
-45 3 93 22
-46 2 90 13
-47 4 91 31
-48 5 92 1
-49 1 92 40
-50 14 90 2
-51 32 91 4
-52 23 93 3
-53 64 109 108
-54 107 109 65
-55 68 109 64
-56 69 107 65
-57 46 106 66
-58 46 141 106
-59 83 104 61
-60 82 103 60
-61 84 105 54
-62 42 103 82
-63 45 105 84
-64 43 104 83
-65 63 89 66
-66 47 86 67
-67 48 87 69
-68 49 88 68
-69 55 100 87
-70 59 101 88
-71 67 142 63
-72 28 72 27
-73 37 71 36
-74 19 73 18
-75 10 70 9
-76 19 120 73
-77 28 118 72
-78 10 117 70
-79 37 116 71
-80 67 130 47
-81 69 131 48
-82 68 133 49
-83 66 132 63
-84 63 130 67
-85 64 133 68
-86 65 131 69
-87 87 100 42
-88 88 101 45
-89 66 111 46
-90 108 142 67
-91 41 109 107
-92 108 109 41
-93 70 117 55
-94 73 120 56
-95 72 118 57
-96 71 116 59
-97 7 80 6
-98 16 75 15
-99 34 78 33
-100 12 74 11
-101 21 76 20
-102 25 77 24
-103 30 81 29
-104 39 79 38
-105 57 102 86
-106 86 102 44
-107 70 94 9
-108 73 95 18
-109 72 97 27
-110 71 96 36
-111 58 98 78
-112 76 99 56
-113 79 101 59
-114 74 100 55
-115 46 137 73
-116 47 135 72
-117 49 134 71
-118 48 136 70
-119 99 106 56
-120 64 110 85
-121 65 113 84
-122 62 112 82
-123 82 111 62
-124 77 104 51
-125 53 102 81
-126 75 103 50
-127 80 105 52
-128 85 110 44
-129 84 113 45
-130 82 112 42
-131 46 111 82
-132 107 112 62
-133 67 110 108
-134 68 113 109
-135 108 110 64
-136 109 113 65
-137 69 112 107
-138 63 142 89
-139 36 96 35
-140 27 97 26
-141 9 94 8
-142 18 95 17
-143 48 140 87
-144 47 139 86
-145 49 138 88
-146 77 123 24
-147 78 122 33
-148 80 124 6
-149 30 127 81
-150 75 125 15
-151 21 126 76
-152 12 128 74
-153 39 129 79
-154 60 103 75
-155 81 102 57
-156 54 105 80
-157 61 104 77
-158 78 98 53
-159 52 101 79
-160 50 100 74
-161 51 99 76
-162 43 106 99
-163 41 107 89
-164 89 107 62
-165 24 123 23
-166 33 122 32
-167 31 127 30
-168 6 124 5
-169 40 129 39
-170 22 126 21
-171 15 125 14
-172 13 128 12
-173 70 136 94
-174 73 137 95
-175 71 134 96
-176 72 135 97
-177 90 128 13
-178 93 126 22
-179 5 124 92
-180 91 127 31
-181 92 129 40
-182 14 125 90
-183 32 122 91
-184 23 123 93
-185 98 102 53
-186 50 103 100
-187 51 104 99
-188 52 105 101
-189 85 133 64
-190 84 131 65
-191 83 130 63
-192 63 132 83
-193 56 141 73
-194 57 139 72
-195 55 140 70
-196 59 138 71
-197 73 141 46
-198 70 140 48
-199 72 139 47
-200 71 138 49
-201 78 119 58
-202 74 117 11
-203 55 117 74
-204 59 116 79
-205 56 120 76
-206 57 118 81
-207 75 115 60
-208 34 119 78
-209 79 116 38
-210 80 114 54
-211 7 114 80
-212 76 120 20
-213 81 118 29
-214 77 121 61
-215 25 121 77
-216 16 115 75
-217 43 132 106
-218 106 132 66
-219 83 132 43
-220 48 131 84
-221 47 130 83
-222 49 133 85
-223 79 129 52
-224 81 127 53
-225 74 128 50
-226 50 125 75
-227 52 124 80
-228 53 122 78
-229 51 123 77
-230 76 126 51
-231 83 135 47
-232 60 137 82
-233 85 134 49
-234 84 136 48
-235 54 136 84
-236 82 137 46
-237 58 134 85
-238 61 135 83
-239 87 112 69
-240 44 110 86
-241 45 113 88
-242 86 110 67
-243 88 113 68
-244 42 112 87
-245 99 104 43
-246 100 103 42
-247 101 105 45
-248 44 102 98
-249 89 111 66
-250 62 111 89
-251 85 143 58
-252 44 143 85
-253 58 143 98
-254 93 123 51
-255 51 126 93
-256 53 127 91
-257 90 125 50
-258 50 128 90
-259 92 124 52
-260 52 129 92
-261 91 122 53
-262 61 121 97
-263 54 114 94
-264 60 115 95
-265 58 119 96
-266 94 114 8
-267 95 115 17
-268 96 119 35
-269 97 121 26
-270 17 115 16
-271 20 120 19
-272 11 117 10
-273 8 114 7
-274 29 118 28
-275 38 116 37
-276 26 121 25
-277 35 119 34
-278 86 139 57
-279 88 138 59
-280 87 140 55
-281 89 142 41
-282 95 137 60
-283 94 136 54
-284 97 135 61
-285 96 134 58
-286 98 143 44
-287 106 141 56
-288 41 142 108
+1 162 37 198
+2 1 2 162
+37 68 79 38
+38 38 72 68
+39 35 69 62
+40 40 63 59
+41 62 70 35
+42 40 86 63
+43 15 71 14
+44 41 67 65
+45 37 76 67
+46 65 88 41
+47 35 71 15
+48 14 71 38
+49 36 65 64
+50 36 64 63
+51 49 77 66
+52 67 76 66
+53 42 79 68
+54 70 71 35
+55 66 77 64
+56 66 76 49
+57 70 93 71
+58 71 93 72
+59 69 83 62
+60 59 84 40
+61 65 66 64
+62 81 84 43
+63 84 97 43
+64 72 85 68
+65 71 72 38
+66 65 67 66
+67 67 82 37
+68 16 69 35
+69 23 73 33
+70 45 68 46
+71 9 45 34
+72 45 46 34
+73 10 45 9
+74 17 69 16
+75 42 45 10
+76 24 73 23
+77 7 74 37
+78 9 34 8
+79 44 69 17
+80 42 68 45
+81 16 35 15
+82 33 75 22
+83 14 38 13
+84 25 43 24
+85 43 73 24
+86 47 76 74
+87 34 74 8
+88 22 75 21
+89 21 75 39
+90 30 36 29
+91 74 76 37
+92 49 76 47
+93 8 74 7
+94 32 41 31
+95 7 37 6
+96 23 33 22
+97 46 47 34
+98 46 48 47
+99 51 52 50
+100 18 44 17
+101 11 42 10
+102 47 74 34
+103 48 51 50
+104 52 53 50
+105 54 56 53
+106 55 58 54
+107 21 39 20
+108 48 49 47
+109 53 77 50
+110 52 55 54
+111 56 57 53
+112 56 59 57
+113 28 40 27
+114 50 77 49
+115 52 54 53
+116 59 63 57
+117 48 50 49
+118 57 77 53
+119 58 78 54
+120 63 64 57
+121 55 61 58
+122 64 77 57
+123 55 62 61
+124 56 78 60
+125 56 60 59
+126 54 78 56
+127 62 87 70
+128 80 98 44
+129 44 98 83
+130 20 80 19
+131 27 81 26
+132 13 79 12
+133 63 86 36
+134 6 82 5
+135 37 82 6
+136 36 88 65
+137 38 79 13
+138 40 81 27
+139 39 80 20
+140 51 93 87
+141 41 82 67
+142 68 85 46
+143 48 85 51
+144 36 86 29
+145 28 86 40
+146 51 87 52
+147 55 87 62
+148 41 88 31
+149 30 88 36
+150 1 89 32
+151 4 91 25
+152 3 92 18
+153 5 89 1
+154 2 90 11
+155 12 90 2
+156 26 91 4
+157 19 92 3
+158 85 93 51
+159 44 83 69
+160 62 83 61
+161 60 84 59
+162 87 93 70
+163 46 85 48
+164 40 84 81
+165 29 86 28
+166 52 87 55
+167 72 93 85
+168 79 90 12
+169 80 92 19
+170 81 91 26
+171 31 88 30
+172 25 91 43
+173 11 90 42
+174 18 92 44
+175 32 89 41
+176 82 89 5
+177 94 97 60
+178 58 96 95
+179 73 94 33
+180 78 95 94
+181 58 95 78
+182 78 94 60
+183 95 96 75
+184 61 96 58
+185 33 95 75
+186 75 96 39
+187 94 95 33
+188 73 97 94
+189 43 97 73
+190 42 90 79
+191 44 92 80
+192 43 91 81
+193 61 98 96
+194 96 98 39
+195 41 89 82
+196 39 98 80
+197 60 97 84
+198 83 98 61
$EndElements
diff --git a/test/test_bvp.py b/test/test_bvp.py
index 8b00ee1..fff76c1 100644
--- a/test/test_bvp.py
+++ b/test/test_bvp.py
@@ -1,7 +1,7 @@
import unittest
from bvpy import BVP, logger
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
from bvpy.vforms import PoissonForm
from bvpy.boundary_conditions import dirichlet, neumann
from bvpy.utils.pre_processing import create_expression
@@ -12,7 +12,8 @@ class TestPoisson(unittest.TestCase):
"""Initiates the test.
"""
logger.setLevel(50)
- self.domain = Rectangle(x=-2, y=-2, length=4, width=4, cell_size=0.05)
+ self.domain = OccModel(cell_type="triangle", gdim=2, cell_size=0.05)
+ self.domain.add_rectangle(x=-2, y=-2, length=4, width=4)
self.vform = PoissonForm(source=0)
self.bc = dirichlet(0, boundary="all")
diff --git a/test/test_domains.py b/test/test_domains.py
index 0c9d0f0..245f6c6 100644
--- a/test/test_domains.py
+++ b/test/test_domains.py
@@ -2,8 +2,8 @@ import unittest
# -- test class
-from bvpy.domains import (Rectangle, Disk, Sphere, HemiSphere, Ellipsoid,
- Torus, CustomPolygonalDomain, Cylinder, Cube)
+import gmsh
+from bvpy.domains import CustomGMSH, OccModel
import fenics as fe
import numpy as np
@@ -20,97 +20,106 @@ class TestPrimitive(unittest.TestCase):
"""
def test_input(self):
- domain = Rectangle(cell_size=[0.05, 0.1], clear=True)
+ domain = OccModel(cell_type="triangle", gdim=2, cell_size=(0.05, 0.1))
+ domain.add_rectangle()
domain.set_cell_size_from_curvature(2)
def test_AbstractDomain_throws_exception(self):
# - try to put dummy input in an AbstractDomain class
- self.assertRaises(ValueError, Rectangle,
- cell_type='tri') # invalid cell_type, can be {'line', 'triangle', 'tetra'}
- self.assertRaises(ValueError, Rectangle,
- cell_size='1') # invalid cell_size, should be a list of two numbers or one number
- self.assertRaises(ValueError, Rectangle, cell_size=[1]) # invalid cell_size
- self.assertRaises(ValueError, Rectangle, cell_size=[1, 2, 4]) # invalid cell_size
- self.assertRaises(ValueError, Rectangle, verbosity=1) # invalid verbosity, can be {True, False}
- self.assertRaises(ValueError, Rectangle, dim=2,
- cell_type='tetra') # invalid embedded dimension : should be equal or larger to the dimension of finite element
- self.assertRaises(ValueError, Rectangle, cell_type='triangle', algorithm='MyAlgo') # wrong algo for 2D
- self.assertRaises(ValueError, Rectangle, cell_type='tetra', algorithm='MyAlgo') # wrong algo for 3D
-
- domain = Rectangle()
+ self.assertRaises(ValueError, OccModel, gdim=2, cell_type='tri') # invalid cell_type
+ self.assertRaises(ValueError, OccModel, gdim=2, cell_type='triangle', cell_size='1') # invalid cell_size
+ self.assertRaises(ValueError, OccModel, gdim=2, cell_type='triangle', cell_size=[1]) # invalid cell_size
+ self.assertRaises(ValueError, OccModel, gdim=4, cell_type='triangle') # invalid gdim
+ self.assertRaises(ValueError, OccModel, gdim=2, cell_type='triangle', cell_size=[1, 2, 4]) # invalid cell_size
+ self.assertRaises(ValueError, OccModel, gdim=2, cell_type='triangle', verbosity=1) # invalid verbosity
+ self.assertRaises(ValueError, OccModel, gdim=2, cell_type='triangle', algorithm='MyAlgo') # wrong algo for 2D
+ self.assertRaises(ValueError, OccModel, gdim=2, cell_type='tetra', algorithm='MyAlgo') # wrong algo for 3D
+
+ domain = OccModel(cell_type="triangle", gdim=2)
self.assertRaises(ValueError, domain.set_cell_size_from_curvature, 'a')
def test_rectangle(self):
- domain = Rectangle(cell_size=0.1, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=2, cell_size=0.1)
+ domain.add_rectangle()
domain.discretize()
domain.discretize()
self.assertEqual(domain.mesh.coordinates()[:, 0].max(), 1)
- domain._details
+ domain.info()
def test_hemisphere(self):
- domain = HemiSphere(cell_size=0.1, clear=True)
+ domain = OccModel(cell_type="tetra", gdim=3, cell_size=0.1)
+ domain.add_hemisphere()
domain.discretize()
self.assertEqual(domain.mesh.coordinates()[:, 2].min(), 0)
def test_disk(self):
- domain = Disk(cell_size=0.05, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=2, cell_size=0.05)
+ domain.add_disk()
domain.discretize()
self.assertIsNotNone(domain.mesh)
self.assertAlmostEqual(domain.mesh.coordinates()[:, 1].min(), -1, delta=1e-3)
def test_sphere(self):
- domain = Sphere(cell_size=0.05, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=3, cell_size=0.05)
+ domain.add_sphere()
domain.discretize()
self.assertIsNotNone(domain.mesh)
self.assertAlmostEqual(domain.mesh.coordinates()[:, 1].max(), 1, delta=1e-3)
def test_ellipsoid(self):
- domain = Ellipsoid(cell_size=0.05, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=3, cell_size=0.05)
+ domain.add_ellipsoid()
domain.discretize()
self.assertIsNotNone(domain.mesh)
self.assertAlmostEqual(domain.mesh.coordinates()[:, 0].min(), -1, delta=1e-3)
def test_torus(self):
- domain = Torus(cell_size=0.1, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=3, cell_size=0.1)
+ domain.add_torus()
domain.discretize()
self.assertIsNotNone(domain.mesh)
self.assertAlmostEqual(domain.mesh.coordinates()[:, 0].min(), -3, delta=1e-3)
def test_cylinder(self):
- domain = Cylinder(open=True, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=3)
+ domain.add_cylinder()
domain.discretize()
self.assertIsNotNone(domain.mesh)
def test_cube(self):
- domain = Cube(clear=True)
+ domain = OccModel(cell_type="tetra", gdim=3)
+ domain.add_cube()
domain.discretize()
self.assertIsNotNone(domain.mesh)
domain.rotate((0, 0, 0), (1, 0, 0), 0.1) # try rotation
def test_csg_sub(self):
- d1 = Rectangle(x=0, y=0, length=1, width=1, clear=True)
- d2 = Rectangle(x=0.5, y=0, length=0.5, width=1)
- d = d1 - d2
- d.discretize()
- self.assertEqual(d.mesh.coordinates()[:, 0].max(), 0.5)
+ domain = OccModel(cell_type="triangle", gdim=2)
+ rect1 = domain.add_rectangle(x=0, y=0, length=1, width=1)
+ rect2 = domain.add_rectangle(x=0.5, y=0, length=0.5, width=1)
+ domain.cut(rect1, rect2)
+ domain.discretize()
+ self.assertEqual(domain.mesh.coordinates()[:, 0].max(), 0.5)
def test_csg_add(self):
- d1 = Rectangle(x=0, y=0, length=0.5, width=1, clear=True)
- d2 = Rectangle(x=0.5, y=0, length=0.5, width=1)
- d = d1 + d2
- d.discretize()
- self.assertEqual(d.mesh.coordinates()[:, 0].max(), 1)
+ domain = OccModel(cell_type="triangle", gdim=2)
+ rect1 = domain.add_rectangle(x=0, y=0, length=0.5, width=1)
+ rect2 = domain.add_rectangle(x=0.5, y=0, length=0.5, width=1)
+ domain.fuse(rect1, rect2)
+ domain.discretize()
+ self.assertEqual(domain.mesh.coordinates()[:, 0].max(), 1)
def test_csg_intersect(self):
- d1 = Rectangle(x=0, y=0, length=1.5, width=1, clear=True)
- d2 = Rectangle(x=1, y=0, length=1, width=1)
- d = d1 & d2
- d.discretize()
- self.assertEqual(d.mesh.coordinates()[:, 0].min(), 1)
+ domain = OccModel(cell_type="triangle", gdim=2)
+ rect1 = domain.add_rectangle(x=0, y=0, length=1.5, width=1)
+ rect2 = domain.add_rectangle(x=1, y=0, length=1, width=1)
+ domain.intersect(rect1, rect2)
+ domain.discretize()
+ self.assertEqual(domain.mesh.coordinates()[:, 0].min(), 1)
- def test_primitive_custompolygonal(self):
+ def test_primitive_conformalpolygon(self):
p = [[0, 0, 0],
[1, 0, 0],
[1, 1, 0],
@@ -119,7 +128,8 @@ class TestPrimitive(unittest.TestCase):
c = [[0, 1, 2, 3],
[1, 4, 2]]
- domain = CustomPolygonalDomain(p, c, cell_size=0.05, isoriented=True, clear=True)
+ domain = OccModel(gdim=3, cell_type="triangle", cell_size=0.05)
+ domain.add_conformal_polygon(p, c, isoriented=True)
domain.discretize()
self.assertEqual(domain.mesh.coordinates()[:, 0].max(), 2)
@@ -132,7 +142,8 @@ class TestPrimitive(unittest.TestCase):
c = [[0, 3, 2, 1],
[1, 2, 4]]
- domain = CustomPolygonalDomain(p, c, cell_size=0.05, isoriented=False, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=3)
+ domain.add_conformal_polygon(p, c, isoriented=False)
domain.discretize()
self.assertEqual(domain.mesh.coordinates()[:, 0].max(), 2)
@@ -145,7 +156,8 @@ class TestPrimitive(unittest.TestCase):
c = [[0, 3, 2, 1],
[1, 2, 4]]
- domain = CustomPolygonalDomain(p, c, markers=[0, 1], cell_size=0.05, isoriented=False, clear=True)
+ domain = OccModel(cell_type="triangle", gdim=3, cell_size=0.05)
+ domain.add_conformal_polygon(p, c, markers=[0, 1], isoriented=False)
domain.discretize()
ind = 10
label = domain.cdata.array()[ind]
@@ -156,33 +168,36 @@ class TestPrimitive(unittest.TestCase):
self.assertEqual(label, 1)
def test_domain_repr_str(self):
- dom = Rectangle(verbosity=True)
- self.assertEqual(dom.__repr__()[1:10], 'Rectangle')
- self.assertEqual(dom.__str__()[:22], 'Mesh not generate yet.')
+ domain = OccModel(cell_type="triangle", gdim=2, verbosity=True)
+ self.assertEqual(domain.__repr__()[1:9], 'OccModel')
+ self.assertTrue("Mesh not generated" in str(domain))
def test_discretize_square_cube(self):
- sqr = Rectangle(cell_type='line')
- sqr.discretize()
- sqr.info()
-
- def test_sphere_sub(self):
- big = Sphere(radius=10, cell_type='tetra')
- sml = Sphere(radius=9, cell_type='tetra')
- shl = big - sml
- shl.geometry()
- shl.discretize()
-
- def test_sphere_add(self):
- sph = Sphere(radius=2, cell_type='tetra')
- sph_2 = Sphere(radius=2, center=[0, 0, 2], cell_type='tetra')
- intsec = sph + sph_2
- intsec.discretize()
-
- def test_sphere_intersec(self):
- sph = Sphere(radius=2, cell_type='tetra')
- sph_2 = Sphere(radius=2, center=[0, 0, 2], cell_type='tetra')
- intsec = sph & sph_2
- intsec.discretize()
+ domain = OccModel(cell_type="line", gdim=2)
+ domain.add_rectangle()
+ domain.discretize()
+ domain.info()
+
+ def test_sphere_cut(self):
+ domain = OccModel(cell_type="tetra", gdim=3)
+ big = domain.add_sphere(radius=10)
+ sml = domain.add_sphere(radius=6)
+ domain.cut(big, sml)
+ domain.discretize()
+
+ def test_sphere_fuse(self):
+ domain = OccModel(cell_type="tetra", gdim=3)
+ sph1 = domain.add_sphere(radius=2)
+ sph2 = domain.add_sphere(center=(0,0,2), radius=2)
+ domain.fuse(sph1, sph2)
+ domain.discretize()
+
+ def test_sphere_intersect(self):
+ domain = OccModel(cell_type="tetra", gdim=3)
+ sph1 = domain.add_sphere(radius=2)
+ sph2 = domain.add_sphere(center=(0,0,2), radius=2)
+ domain.intersect(sph1, sph2)
+ domain.discretize()
def test_mesh_conversion(self):
# - Assert that we can convert labelled meshes (entities in the gmsh mesh has an associated label + name)
@@ -193,135 +208,137 @@ class TestPrimitive(unittest.TestCase):
### - 1D element dimension : line - ###
###########################################
# - Create a domain as cylinder with two sub-domains
- geometry = Disk(radius=L, cell_type='line', cell_size=0.4, clear=True)
- geometry.factory.remove([(2, 1)]) # no need surface
- circle = geometry.factory.getEntities(dim=1)
+ domain = OccModel(cell_type="line", gdim=2, cell_size=0.4)
+ domain.add_disk(radius=L)
+ domain.factory.remove([(2, 1)]) # no need surface
+ circle = domain.factory.getEntities(dim=1)
# - Cut the circle by adding first two points
- geometry.factory.addPoint(x=L, y=0, z=0, tag=2)
- geometry.factory.addPoint(x=-L, y=0, z=0, tag=3)
+ domain.factory.addPoint(x=L, y=0, z=0, tag=2)
+ domain.factory.addPoint(x=-L, y=0, z=0, tag=3)
# - Divide the circle according to the points : create sub-domains
- geometry.factory.fragment(objectDimTags=circle, toolDimTags=[(0, 2), (0, 3)])
+ domain.factory.fragment(objectDimTags=circle, toolDimTags=[(0, 2), (0, 3)])
- geometry.factory.synchronize() # synchronize prior sub-domains labellisation
+ domain.factory.synchronize() # synchronize prior sub-domains labellisation
# - Label the surface entities to create two sub-domains (0 <= x <= L/2 and L/2 <= x <= L)
- geometry.model.addPhysicalGroup(1, [1], 1) # add label 1 to line 1
- geometry.model.setPhysicalName(1, 1, 'left_side') # add label title
- geometry.model.addPhysicalGroup(1, [2], 2) # add label 2 ti line 2
- geometry.model.setPhysicalName(1, 2, 'right_side') # add label title
+ gmsh.model.addPhysicalGroup(1, [1], 1) # add label 1 to line 1
+ gmsh.model.setPhysicalName(1, 1, 'left_side') # add label title
+ gmsh.model.addPhysicalGroup(1, [2], 2) # add label 2 ti line 2
+ gmsh.model.setPhysicalName(1, 2, 'right_side') # add label title
# - Discretize
- geometry.discretize()
+ domain.discretize()
# - Assert that the label and label_names are stored in the AbstractModel
- self.assertTrue(set(geometry.cdata.array()) == {1, 2}) # 1D : label are stored
+ self.assertTrue(set(domain.cdata.array()) == {1, 2}) # 1D : label are stored
self.assertTrue(
- set(list(geometry.sub_domain_names.keys())) == {1, 2}) # 1D: label are also stored in sub_domain_names
+ set(list(domain.sub_domain_names.keys())) == {1, 2}) # 1D: label are also stored in sub_domain_names
self.assertTrue(
- set(list(geometry.sub_domain_names.values())) == {'left_side', 'right_side'}) # 1D: name of the sub-domains
+ set(list(domain.sub_domain_names.values())) == {'left_side', 'right_side'}) # 1D: name of the sub-domains
### - 2D element dimension : triangle - ###
###########################################
# - Create a domain as cylinder with two sub-domains
- geometry = Cylinder(x=0, dx=L, dz=0, r=r, cell_type='triangle', cell_size=0.2, clear=True)
- geometry.factory.remove([(3, 1)]) # no need volume
- surfaces = geometry.factory.getEntities(dim=2)
+ domain = OccModel(cell_type="triangle", gdim=3, cell_size=0.2)
+ domain.add_cylinder(x=0, dx=L, dz=0, r=r)
+ domain.factory.remove([(3, 1)]) # no need volume
+ surfaces = domain.factory.getEntities(dim=2)
# - Create a circle at the middle of the cylinder
- circle = geometry.factory.addCircle(x=L / 2, y=0, z=0, r=r)
- geometry.factory.rotate([(1, circle)], x=L / 2, y=0, z=0, ax=0, ay=1, az=0,
+ circle = domain.factory.addCircle(x=L / 2, y=0, z=0, r=r)
+ domain.factory.rotate([(1, circle)], x=L / 2, y=0, z=0, ax=0, ay=1, az=0,
angle=np.pi / 2) # rotate the circle
# - Divide the cylinder according to the circle : create sub-domains
- geometry.factory.fragment(objectDimTags=surfaces, toolDimTags=[(1, circle)])
+ domain.factory.fragment(objectDimTags=surfaces, toolDimTags=[(1, circle)])
- geometry.factory.synchronize() # synchronize prior sub-domains labellisation
+ domain.factory.synchronize() # synchronize prior sub-domains labellisation
# - Label the surface entities to create two sub-domains (0 <= x <= L/2 and L/2 <= x <= L)
- geometry.model.addPhysicalGroup(2, [2, 4], 1) # add label 1 to surfaces 2 and 4
- geometry.model.setPhysicalName(2, 1, 'left_side') # add label title
- geometry.model.addPhysicalGroup(2, [3, 5], 2) # add label 1 to surfaces 3 and 5
- geometry.model.setPhysicalName(2, 2, 'right_side') # add label title
+ gmsh.model.addPhysicalGroup(2, [2, 4], 1) # add label 1 to surfaces 2 and 4
+ gmsh.model.setPhysicalName(2, 1, 'left_side') # add label title
+ gmsh.model.addPhysicalGroup(2, [3, 5], 2) # add label 1 to surfaces 3 and 5
+ gmsh.model.setPhysicalName(2, 2, 'right_side') # add label title
# - Discretize
- geometry.discretize()
+ domain.discretize()
# - Assert that the label and label_names are stored in the AbstractModel
- self.assertTrue(set(geometry.cdata.array()) == {1, 2}) # 2D: label are stored
+ self.assertTrue(set(domain.cdata.array()) == {1, 2}) # 2D: label are stored
self.assertTrue(
- set(list(geometry.sub_domain_names.keys())) == {1, 2}) # 2D: label are also stored in sub_domain_names
+ set(list(domain.sub_domain_names.keys())) == {1, 2}) # 2D: label are also stored in sub_domain_names
self.assertTrue(
- set(list(geometry.sub_domain_names.values())) == {'left_side', 'right_side'}) # 2D: name of the sub-domains
+ set(list(domain.sub_domain_names.values())) == {'left_side', 'right_side'}) # 2D: name of the sub-domains
### - 3D element dimension : tetra - ###
###########################################
# - Create a domain as cylinder with two sub-domains
- geometry = Cylinder(x=0, dx=L, dz=0, r=r, cell_type='tetra', cell_size=0.2, clear=True)
- volumes = geometry.factory.getEntities(dim=3)
+ domain = OccModel(cell_type="tetra", gdim=3, cell_size=0.2)
+ domain.add_cylinder(x=0, dx=L, dz=0, r=r)
+ volumes = domain.factory.getEntities(dim=3)
# - Create a disk at the middle of the cylinder
- circle = geometry.factory.addCircle(x=L / 2, y=0, z=0, r=r)
- curve_loop = geometry.factory.addCurveLoop([circle])
- disk = geometry.factory.addPlaneSurface([curve_loop])
+ circle = domain.factory.addCircle(x=L / 2, y=0, z=0, r=r)
+ curve_loop = domain.factory.addCurveLoop([circle])
+ disk = domain.factory.addPlaneSurface([curve_loop])
- geometry.factory.rotate([(2, disk)], x=L / 2, y=0, z=0, ax=0, ay=1, az=0,
+ domain.factory.rotate([(2, disk)], x=L / 2, y=0, z=0, ax=0, ay=1, az=0,
angle=np.pi / 2) # rotate the circle
# - Divide the cylinder according to the circle : create sub-domains
- geometry.factory.fragment(objectDimTags=volumes, toolDimTags=[(2, disk)])
+ domain.factory.fragment(objectDimTags=volumes, toolDimTags=[(2, disk)])
- geometry.factory.synchronize() # synchronize prior sub-domains labellisation
+ domain.factory.synchronize() # synchronize prior sub-domains labellisation
# - Label the volume entities to create two sub-domains (0 <= x <= L/2 and L/2 <= x <= L)
- geometry.model.addPhysicalGroup(3, [1], 1) # add label 1 to volume 1
- geometry.model.setPhysicalName(3, 1, 'left_side') # add label title
- geometry.model.addPhysicalGroup(3, [2], 2) # add label 2 to volume 2
- geometry.model.setPhysicalName(3, 2, 'right_side') # add label title
+ gmsh.model.addPhysicalGroup(3, [1], 1) # add label 1 to volume 1
+ gmsh.model.setPhysicalName(3, 1, 'left_side') # add label title
+ gmsh.model.addPhysicalGroup(3, [2], 2) # add label 2 to volume 2
+ gmsh.model.setPhysicalName(3, 2, 'right_side') # add label title
# - Label also the surface entities to create two sub-boundaries (0 <= x <= L/2 and L/2 <= x <= L)
- geometry.model.addPhysicalGroup(2, [8], 1) # add label 1 to surfaces 8
- geometry.model.setPhysicalName(2, 1, 'left_boundary') # add label title
- geometry.model.addPhysicalGroup(2, [6], 2) # add label 2 to surfaces 6
- geometry.model.setPhysicalName(2, 2, 'right_boundary') # add label title
- geometry.model.addPhysicalGroup(2, [5, 7], 3) # add label 3 to surfaces [5, 7]
- geometry.model.setPhysicalName(2, 3, 'other_boundary') # add label title
+ gmsh.model.addPhysicalGroup(2, [8], 1) # add label 1 to surfaces 8
+ gmsh.model.setPhysicalName(2, 1, 'left_boundary') # add label title
+ gmsh.model.addPhysicalGroup(2, [6], 2) # add label 2 to surfaces 6
+ gmsh.model.setPhysicalName(2, 2, 'right_boundary') # add label title
+ gmsh.model.addPhysicalGroup(2, [5, 7], 3) # add label 3 to surfaces [5, 7]
+ gmsh.model.setPhysicalName(2, 3, 'other_boundary') # add label title
# - Discretize
- geometry.discretize()
+ domain.discretize()
# - Assert that the label and label_names are stored in the AbstractModel
- self.assertTrue(set(geometry.cdata.array()) == {1, 2}) # 3D: label are stored
+ self.assertTrue(set(domain.cdata.array()) == {1, 2}) # 3D: label are stored
self.assertTrue(
- set(list(geometry.sub_domain_names.keys())) == {1, 2}) # 3D: label are also stored in sub_domain_names
+ set(list(domain.sub_domain_names.keys())) == {1, 2}) # 3D: label are also stored in sub_domain_names
self.assertTrue(
- set(list(geometry.sub_domain_names.values())) == {'left_side', 'right_side'}) # 3D: name of the sub-domains
- self.assertTrue(len({1, 2, 3} - set(geometry.bdata.array())) == 0) # 3D: all boundary labels are stored
+ set(list(domain.sub_domain_names.values())) == {'left_side', 'right_side'}) # 3D: name of the sub-domains
+ self.assertTrue(len({1, 2, 3} - set(domain.bdata.array())) == 0) # 3D: all boundary labels are stored
self.assertTrue(
- set(list(geometry.sub_boundary_names.keys())) == {1, 2, 3}) # 3D: label are also stored in sub_domain_names
- self.assertTrue(set(list(geometry.sub_boundary_names.values())) == {'left_boundary', 'right_boundary',
+ set(list(domain.sub_boundary_names.keys())) == {1, 2, 3}) # 3D: label are also stored in sub_domain_names
+ self.assertTrue(set(list(domain.sub_boundary_names.values())) == {'left_boundary', 'right_boundary',
'other_boundary'}) # 3D: name of the sub-domains
### - RunTimeError - ###
###########################################
- # Let ele_dim be the topological dimension of the finite element of the gmsh mesh (ele_dim can be defined from
- # _cell_type in the AbstractDomain class).
- # If no PhysicalGroup is defined for entities at dim. ele_dim whereas for some other entities
- # at dim. < ele_dim PhysicialGroup have been set => error with RunTimeError
- def error_test():
- geometry = Cylinder(x=0, dx=L, dz=0, r=r, cell_type='tetra', clear=True)
+ # If not PhysicalGroups are defined, a default one is computed
+ # We check if the default one is computed correctly
+ domain = OccModel(cell_type="tetra", gdim=3)
+ domain.add_cylinder(x=0, dx=L, dz=0, r=r)
- geometry.factory.synchronize() # synchronize prior sub-domains labellisation
+ domain.factory.synchronize() # synchronize prior sub-domains labellisation
- # - Label surfaces entities WITHOUT labelling volume entities
- #geometry.model.addPhysicalGroup(3, [1], 1) # add label 1 to volume 1
- #geometry.model.setPhysicalName(3, 1, 'my_domain') # add label title
+ # - Label surfaces entities WITHOUT labelling volume entities
+ gmsh.model.addPhysicalGroup(2, [1, 2, 3], 1) # add label 1 to surfaces 1, 2 and 3
+ gmsh.model.setPhysicalName(2, 1, 'my_boundary') # add label title
- geometry.model.addPhysicalGroup(2, [1, 2, 3], 1) # add label 1 to surfaces 1, 2 and 3
- geometry.model.setPhysicalName(2, 1, 'my_boundary') # add label title
-
- # - Discretize
- geometry.discretize()
+ # Compute mesh
+ domain.discretize()
- self.assertRaises(RuntimeError, error_test)
+ # Assertions
+ self.assertTrue(set(list(domain.sub_domain_names.keys())) == {2})
+ self.assertTrue(set(list(domain.sub_domain_names.values())) == {'domain_automatically_generated'})
+ self.assertTrue(set(list(domain.sub_boundary_names.keys())) == {1})
+ self.assertTrue(set(list(domain.sub_boundary_names.values())) == {'my_boundary'})
diff --git a/test/test_elasticity.py b/test/test_elasticity.py
index 3c8bbbd..5d857c9 100644
--- a/test/test_elasticity.py
+++ b/test/test_elasticity.py
@@ -3,7 +3,7 @@ import unittest
import numpy as np
from bvpy import BVP, logger
-from bvpy.domains import CustomDomain
+from bvpy.domains import FixedDomain
from bvpy.vforms.elasticity import LinearElasticForm, HyperElasticForm
from bvpy.vforms.elasticity import NeoHookeanPotential, StVenantKirchoffPotential, QuadraticGreenLagrangePotential
from bvpy.boundary_conditions import dirichlet
@@ -21,7 +21,7 @@ class TestLinearElasticityProblem(unittest.TestCase):
mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(self.L, self.W),
100, 10, "crossed")
- self.domain = CustomDomain(mesh.coordinates(), mesh.cells(), dim=2)
+ self.domain = FixedDomain(mesh.coordinates(), mesh.cells())
dummy_mesh = fe.UnitSquareMesh(1, 1)
self.dummy_V = fe.VectorFunctionSpace(dummy_mesh, "P", 1)
@@ -94,7 +94,7 @@ class TestLinearElasticityProblem(unittest.TestCase):
# - create deformed mesh
mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(self.L, self.W), 100, 10, "crossed")
- domain = CustomDomain(mesh.coordinates(), mesh.cells(), dim=2)
+ domain = FixedDomain(mesh.coordinates(), mesh.cells())
domain.move(bvp.solution)
def test_linear_hyper(self):
diff --git a/test/test_gbvp.py b/test/test_gbvp.py
index 3361dad..5ee521f 100644
--- a/test/test_gbvp.py
+++ b/test/test_gbvp.py
@@ -3,7 +3,7 @@ import tempfile
import numpy as np
from bvpy.gbvp import GBVP
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
from bvpy.vforms.elasticity import StVenantKirchoffPotential
from bvpy.vforms.plasticity import MorphoElasticForm
from bvpy.vforms.plasticity import ConstantGrowth
@@ -14,7 +14,8 @@ class TestGBVP(unittest.TestCase):
def setUp(self):
"""Initiates the test.
"""
- self.domain = Rectangle(x=-2, y=-2, length=4, width=4, cell_size=0.05, dim=2)
+ self.domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.05)
+ self.domain.add_rectangle(x=-2, y=-2, length=4, width=4)
self.bc = dirichlet(val=[0, 0], boundary="near(x, -2)")
potential_energy = StVenantKirchoffPotential()
diff --git a/test/test_geometry.py b/test/test_geometry.py
index 75ee0ee..38d580c 100644
--- a/test/test_geometry.py
+++ b/test/test_geometry.py
@@ -3,8 +3,7 @@ import unittest
import numpy as np
import fenics as fe
-from bvpy.domains import Disk, Sphere
-
+from bvpy.domains import OccModel
# -- module to test
from bvpy.domains.geometry import *
@@ -24,30 +23,30 @@ class TestGeometry(unittest.TestCase):
def test_boundary_normal(self):
# 3D mesh
- disk = Sphere(cell_size=0.1, dim=3, cell_type="tetra", clear=True)
- disk.discretize()
- bnd_nmrl = boundary_normal(disk)
- np.testing.assert_allclose(bnd_nmrl(1, 0, 0),
- np.array([1, 0, 0]), atol=1e-2)
+ domain = OccModel(gdim=3, cell_type="tetra", cell_size=0.3)
+ domain.add_sphere()
+ domain.discretize()
+ bnd_nmrl = boundary_normal(domain)
+ np.testing.assert_allclose(bnd_nmrl(0, 0, -1), np.array([0, 0, -1]), atol=1e-1)
# 2D mesh
- disk = Disk(cell_size=0.05, dim=2, clear=True)
- disk.discretize()
- bnd_nmrl = boundary_normal(disk)
- np.testing.assert_allclose(bnd_nmrl(1, 0),
- np.array([1, 0]), atol=1e-15)
+ domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.05)
+ domain.add_disk()
+ domain.discretize()
+ bnd_nmrl = boundary_normal(domain)
+ np.testing.assert_allclose(bnd_nmrl(1, 0), np.array([1, 0]), atol=1e-15)
# 2D mesh embedded in 3D
- disk = Disk(cell_size=0.05, dim=3, clear=True)
- disk.discretize()
- bnd_nmrl = boundary_normal(disk)
- np.testing.assert_allclose(bnd_nmrl(1, 0, 0),
- np.array([1, 0, 0]), atol=1e-15)
+ domain = OccModel(gdim=3, cell_type="triangle", cell_size=0.05)
+ domain.add_hemisphere()
+ domain.discretize()
+ bnd_nmrl = boundary_normal(domain)
+ np.testing.assert_allclose(bnd_nmrl(1, 0, 0), np.array([0, 0, -1]), atol=1e-1)
def test_normal(self):
- disk = Disk(cell_size=0.05, clear=True)
- disk.discretize()
- n = normal(disk, interior=[0, 0, -1])
- np.testing.assert_allclose(n(1, 0, 0),
- np.array([0, 0, 1]), atol=1e-15)
+ domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.05)
+ domain.add_disk()
+ domain.discretize()
+ n = normal(domain, interior=[0, 0, -1])
+ np.testing.assert_allclose(n(1, 0), np.array([0, 0, 1]), atol=1e-15)
def test_init_orientation_3D(self):
mesh = fe.UnitSquareMesh(1, 1)
diff --git a/test/test_gmshio.py b/test/test_gmshio.py
index 700be98..1cd678f 100644
--- a/test/test_gmshio.py
+++ b/test/test_gmshio.py
@@ -24,6 +24,8 @@ class TestGMSHio(unittest.TestCase):
"""Runs the test. Creates a gmsh mesh, writes it to a file, reads it with read_from_msh, and asserts the values.
"""
# Create and write the gmsh mesh
+ if gmsh.isInitialized(): # The cell size is keep even when gmsh is cleared
+ gmsh.finalize()
gmsh.initialize()
create_gmsh_mesh(self.tdim, self.gdim)
gmsh.write(str(self.file_path))
diff --git a/test/test_helmholtz.py b/test/test_helmholtz.py
index 66e2992..86721d8 100644
--- a/test/test_helmholtz.py
+++ b/test/test_helmholtz.py
@@ -3,7 +3,7 @@ import unittest
# -- test class
from bvpy.vforms import HelmholtzForm
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
from bvpy.utils.pre_processing import create_expression
from bvpy.boundary_conditions import dirichlet, neumann
from bvpy import BVP
@@ -23,9 +23,10 @@ class TestHelmholtz(unittest.TestCase):
pass
def test_bvp_creator(self):
- rect = Rectangle(cell_size=0.1)
+ domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.1)
+ domain.add_rectangle()
hlmz = HelmholtzForm(linear_coef=50**2)
hlmz.add_point_source([.5, .5, 0.], 1.)
diri = dirichlet(0, boundary='all')
- prblm = BVP(domain=rect, vform=hlmz, bc=diri)
+ prblm = BVP(domain=domain, vform=hlmz, bc=diri)
prblm.solve()
diff --git a/test/test_ibvp.py b/test/test_ibvp.py
index 18d2a70..a06a320 100644
--- a/test/test_ibvp.py
+++ b/test/test_ibvp.py
@@ -1,7 +1,7 @@
import unittest
from bvpy import IBVP, logger
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
from bvpy.vforms import PoissonForm
from bvpy.solvers.time_integration import ExplicitEuler
from bvpy.boundary_conditions import dirichlet
@@ -15,7 +15,8 @@ class TestPoisson(unittest.TestCase):
"""Initiates the test.
"""
logger.setLevel(50)
- self.domain = Rectangle(x=-2, y=-2, length=4, width=4, cell_size=0.05)
+ self.domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.05)
+ self.domain.add_rectangle(x=-2, y=-2, length=4, width=4)
self.vform = PoissonForm(source=0)
self.bc = dirichlet(0, boundary="all")
self.dir = tempfile.TemporaryDirectory()
diff --git a/test/test_io.py b/test/test_io.py
index d293ceb..ed65d4f 100644
--- a/test/test_io.py
+++ b/test/test_io.py
@@ -4,10 +4,10 @@ from pathlib import Path
from fenics import UnitSquareMesh
from bvpy import BVP
-from bvpy.domains.primitives import Rectangle
+from bvpy.domains import OccModel
from bvpy.vforms.poisson import PoissonForm
from bvpy.boundary_conditions.dirichlet import dirichlet
-from bvpy.domains.custom_gmsh import CustomDomainGmsh
+from bvpy.domains.fixed_gmsh import FixedGMSH
# -- test class
from bvpy.utils.io import *
@@ -30,7 +30,7 @@ class TestIO(unittest.TestCase):
def test_read_gmsh_data(self):
path = str(self.test_data_dir / "mesh_example.msh")
- dom = CustomDomainGmsh(path)
+ dom = FixedGMSH(path, gdim=2)
def test_read_polygonal_data_ply(self):
path = str(self.test_data_dir / "ply_example.ply")
@@ -77,13 +77,15 @@ class TestIO(unittest.TestCase):
write_ply(mesh, path)
def test_write_ply_domain(self):
- domain = Rectangle()
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
domain.discretize()
path = self.tmp.name+'/domain.ply'
write_ply(domain, path)
def test_write_ply_wrong_type(self):
- domain = Rectangle()
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
domain.discretize()
path = self.tmp.name+'/domain.ply'
write_ply(domain.cdata, path)
@@ -94,31 +96,38 @@ class TestIO(unittest.TestCase):
self.assertEqual(numerize('three'), 'three')
def test_save_domain(self):
- domain = Rectangle()
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
domain.discretize()
path = self.tmp.name+'/save_test.xdmf'
save(domain, path)
def test_save_cdata(self):
- domain = Rectangle()
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
domain.discretize()
path = self.tmp.name+'/save_test.xdmf'
save(domain.cdata, path)
def test_save_mesh(self):
- domain = Rectangle()
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
path = self.tmp.name+'/save_test.xdmf'
save(domain.mesh, path)
def test_save_problem(self):
- prblm = BVP(domain=Rectangle(),
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
+ prblm = BVP(domain=domain,
vform=PoissonForm(),
bc=dirichlet(0, 'all'))
path = self.tmp.name+'/save_test'
save(prblm, path)
def test_save_solution(self):
- prblm = BVP(domain=Rectangle(),
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
+ prblm = BVP(domain=domain,
vform=PoissonForm(),
bc=dirichlet(0, 'all'))
path = self.tmp.name+'/save_test.xdmf'
diff --git a/test/test_periodic.py b/test/test_periodic.py
index ac54403..a60fb2c 100644
--- a/test/test_periodic.py
+++ b/test/test_periodic.py
@@ -1,7 +1,7 @@
import unittest
import fenics as fe
import numpy as np
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
from bvpy.boundary_conditions.periodic import SquarePeriodicBoundary, LinearPeriodicBoundary
from bvpy.vforms.elasticity import LinearElasticForm
from bvpy.bvp import BVP
@@ -9,7 +9,8 @@ from bvpy.bvp import BVP
class TestPeriodic(unittest.TestCase):
def setUp(self):
- self.dom = Rectangle(length=1, width=1, dim=2, cell_size=0.5)
+ self.dom = OccModel(gdim=2, cell_type="triangle", cell_size=0.5)
+ self.dom.add_rectangle()
self.dom.discretize()
self.vform = LinearElasticForm(source=[0, 0])
diff --git a/test/test_plasticity.py b/test/test_plasticity.py
index fd398c2..4e9d465 100644
--- a/test/test_plasticity.py
+++ b/test/test_plasticity.py
@@ -1,12 +1,11 @@
import unittest
from bvpy import BVP, logger
-from bvpy.domains import Rectangle, CustomDomain
+from bvpy.domains import FixedDomain
from bvpy.vforms.plasticity import MorphoElasticForm
from bvpy.vforms.elasticity import StVenantKirchoffPotential
from bvpy.boundary_conditions import dirichlet
-from bvpy.utils.pre_processing import create_expression
-import tempfile
+
import fenics as fe
import numpy as np
@@ -22,7 +21,7 @@ class TestLinearElasticityProblem(unittest.TestCase):
mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(self.L, self.W),
100, 10, "crossed")
- self.domain = CustomDomain(mesh.coordinates(), mesh.cells(), dim=2)
+ self.domain = FixedDomain(mesh.coordinates(), mesh.cells())
self.rho_g = 1e-2
young = 1e5
diff --git a/test/test_poisson.py b/test/test_poisson.py
index 04e8354..01e0e15 100644
--- a/test/test_poisson.py
+++ b/test/test_poisson.py
@@ -3,7 +3,7 @@ import unittest
# -- test class
from bvpy.vforms import PoissonForm
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
from bvpy.utils.pre_processing import create_expression
from bvpy.boundary_conditions import dirichlet, neumann
from bvpy import BVP
@@ -35,22 +35,25 @@ class TestPoisson(unittest.TestCase):
self.assertEqual(poisson._elements[0].family, "DG")
def test_bvp_creator(self):
- rectangle = Rectangle(cell_size=0.1, clear=True)
+ domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.1)
+ domain.add_rectangle()
poisson = PoissonForm(source=-6)
diri = dirichlet('1 + x*x + 2*y*y', boundary='all')
- bvp = BVP(domain=rectangle, vform=poisson, bc=diri)
+ bvp = BVP(domain=domain, vform=poisson, bc=diri)
def test_bvp_solve(self):
- rectangle = Rectangle(cell_size=0.1, clear=True)
+ domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.1)
+ domain.add_rectangle()
source = create_expression('x*x', degree=2)
poisson = PoissonForm(source=source)
poisson.add_point_source([0.5, 0.5, 0], 1)
diri = dirichlet('1 + x*x + 2*y*y', boundary='all')
- BVP(domain=rectangle, vform=poisson, bc=diri).solve()
+ BVP(domain=domain, vform=poisson, bc=diri).solve()
def test_neumann(self):
- rectangle = Rectangle(cell_size=0.1, clear=True)
+ domain = OccModel(gdim=2, cell_type="triangle", cell_size=0.1)
+ domain.add_rectangle()
poisson = PoissonForm(source=-6)
neum = neumann(100, 'near(y,0)')
diri = dirichlet('1 + x*x + 2*y*y', boundary='near(x,0)')
- BVP(domain=rectangle, vform=poisson, bc=[neum, diri]).solve()
+ BVP(domain=domain, vform=poisson, bc=[neum, diri]).solve()
diff --git a/test/test_post_process.py b/test/test_post_process.py
index 018f6b3..3fec4df 100644
--- a/test/test_post_process.py
+++ b/test/test_post_process.py
@@ -6,7 +6,7 @@ import ufl
from bvpy.utils.post_processing import *
from bvpy import BVP
from bvpy.vforms import PoissonForm
-from bvpy.domains import CustomDomain
+from bvpy.domains import FixedDomain
from fenics import (UnitSquareMesh, FunctionSpace, VectorFunctionSpace, project, Constant)
@@ -102,7 +102,7 @@ class TestPreProcess(unittest.TestCase):
def test_estimate_precision(self):
sol2 = project(Constant(1.1), self.V)
- domain = CustomDomain(self.mesh.coordinates(), self.mesh.cells(), dim=2)
+ domain = FixedDomain(self.mesh.coordinates(), self.mesh.cells())
vform = PoissonForm()
bvp = BVP(domain, vform)
bvp.solution = self.sol
diff --git a/test/test_transport.py b/test/test_transport.py
index c932ed6..3eb4bc3 100644
--- a/test/test_transport.py
+++ b/test/test_transport.py
@@ -6,7 +6,7 @@ import fenics as fe
from bvpy.vforms import TransportForm, CoupledTransportForm
from bvpy.utils.pre_processing import create_expression
from bvpy.solvers.time_integration import ImplicitEuler
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
from bvpy.boundary_conditions import dirichlet
from bvpy.boundary_conditions.periodic import SquarePeriodicBoundary
from bvpy import BVP, IBVP
@@ -37,10 +37,12 @@ class TestTransport(unittest.TestCase):
self.assertEqual(vform._parameters['reaction'].values()[0], 1)
def test_TransportFrom_init_heterogeneous_diffusion(self):
+ domain = OccModel(gdim=2, cell_type="triangle")
+ domain.add_rectangle()
vform = TransportForm(diffusion_coef=create_expression('x', degree=1),
source=create_expression('x', degree=1),
- velocity=[1, 0, 0])
- prblm = BVP(domain=Rectangle(),
+ velocity=[1, 0])
+ prblm = BVP(domain=domain,
vform=vform,
bc=dirichlet(0, boundary='all'))
prblm.solve()
@@ -48,7 +50,9 @@ class TestTransport(unittest.TestCase):
self.assertIs(type(vform._parameters['source']), fe.Expression)
def test_TransportForm_PoissonLike(self):
- prblm = BVP(domain=Rectangle(),
+ domain = OccModel(gdim=2, cell_type="triangle")
+ domain.add_rectangle()
+ prblm = BVP(domain=domain,
vform=TransportForm(source=-6),
bc=dirichlet('1 + x*x + 2*y*y', boundary='all'))
prblm.solve()
@@ -62,14 +66,18 @@ class TestTransport(unittest.TestCase):
vform = TransportForm(reaction=reaction)
vform.set_element_order(4)
- prblm = BVP(domain=Rectangle(cell_size=.05),
+ domain = OccModel(gdim=2, cell_type="triangle", cell_size=.05)
+ domain.add_rectangle()
+
+ prblm = BVP(domain=domain,
vform=vform,
bc=dirichlet('sin(pi*x)*sin(pi*y)', boundary='all'))
prblm.solve()
def test_CoupledTransportForm_init(self):
- dom = Rectangle(length=10, width=10, dim=2)
- dom.discretize()
+ domain = OccModel(gdim=2, cell_type="triangle")
+ domain.add_rectangle(length=10, width=10)
+ domain.discretize()
k1 = 9
k2 = 11
@@ -84,7 +92,7 @@ class TestTransport(unittest.TestCase):
random = np.random.uniform(low=-1,
high=1,
- size=dom.mesh.num_entities(2))
+ size=domain.mesh.num_entities(2))
state = [[1+0.04*k1**2+0.1*r, 0.2*k2+0.1*r] for r in random]
init = create_expression(state)
@@ -95,7 +103,7 @@ class TestTransport(unittest.TestCase):
periodic_bc.map(x, z)
periodic_bc.map(y, z)
- ibvp = IBVP(domain=dom,
+ ibvp = IBVP(domain=domain,
vform=vform,
bc=periodic_bc,
scheme=ImplicitEuler,
diff --git a/test/test_visu.py b/test/test_visu.py
index 38686aa..cf1e09e 100644
--- a/test/test_visu.py
+++ b/test/test_visu.py
@@ -8,7 +8,7 @@ from bvpy.utils.visu import (_cone_plot, _wireframe_plot_mesh,
_surface_plot_function, _surface_plot_mesh,
_plot_dofs, _plot_dirichlet_bc,
_surface_plot_meshfunc, plot)
-from bvpy.domains import Rectangle
+from bvpy.domains import OccModel
import fenics as fe
class TestVisu(unittest.TestCase):
@@ -83,10 +83,11 @@ class TestVisu(unittest.TestCase):
self.assertTrue(os.path.exists(file_path))
def test_plot_domain(self):
- rec = Rectangle()
+ domain = OccModel(gdim=2, cell_type="triangle")
+ domain.add_rectangle()
with tempfile.TemporaryDirectory() as tmpdir:
file_path = os.path.join(tmpdir, "temp_figure.html")
- plot(rec, save_path=file_path)
+ plot(domain, save_path=file_path)
self.assertTrue(os.path.exists(file_path))
def test_plot_meshfunc(self):
diff --git a/test/test_visu_pyvista.py b/test/test_visu_pyvista.py
index 3492a68..6081d96 100644
--- a/test/test_visu_pyvista.py
+++ b/test/test_visu_pyvista.py
@@ -6,7 +6,7 @@ from bvpy.utils.visu_pyvista import (_check_is3D, _check_scalar_bar_title, _visu
_visu_feFunctionSizet, _visu_feMesh, visualize)
from bvpy.bvp import BVP
from bvpy.vforms.elasticity import LinearElasticForm
-from bvpy.domains.primitives import Rectangle
+from bvpy.domains import OccModel
from bvpy.boundary_conditions import dirichlet
import pyvista as pv
import fenics as fe
@@ -22,7 +22,8 @@ class TestVisu(unittest.TestCase):
self.mf = fe.MeshFunction("size_t", self.mesh, 2)
# - Create dummy BVP file with vector field solution
- domain = Rectangle(dim=2)
+ domain = OccModel(cell_type="triangle", gdim=2)
+ domain.add_rectangle()
diri = dirichlet([0, 0], boundary='near(x, 0)')
vform = LinearElasticForm(source=[0, -1e-2],
young=1e5,
--
GitLab