diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py
index 69df9d9f52beffd6109978cdd2bc8b7dca7535a0..a39d6c7318b577aa0ff2d7656dbd30e95050fb11 100644
--- a/src/bvpy/domains/abstract.py
+++ b/src/bvpy/domains/abstract.py
@@ -104,20 +104,20 @@ class AbstractDomain(ABC):
bdata : :class:`MeshFunctionSizet<dolfin.cpp.mesh.MeshFunctionSizet>`
Description of attribute `bdata`.
dx : :class:`Measure<ufl.measure.Measure>`
- Description of attribute `dx`.
- ds : type
- Description of attribute `ds`.
- mesh : type
- Description of attribute `mesh`.
+ 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>`
- The type of mesh model used ????
+ A purely geometrical description of the domain.
model_name : str
- Name of the model.
+ Name of the geometrical model of the domain.
Methods
-------
_mesh_conversion
- Converts a mesh into ????
+ Converts a temporary gmsh mesh into a fenics (dolphin) one.
_set_factory
Sets the factory used to generate the domain elements?
geometry
@@ -220,8 +220,10 @@ class AbstractDomain(ABC):
Yields
-------
- ???
- Description of returned object.
+ :class:`occ<gmsh.model.occ>`
+ or
+ :class:`occ<gmsh.model.geo>`
+ A gmsh model factory to generate mesh from geometry.
"""
if factory == 'occ':
@@ -247,7 +249,7 @@ class AbstractDomain(ABC):
@property
def model_name(self):
- """Names the model.
+ """Name of the geometrical model of the domain.
Returns
-------
@@ -259,18 +261,18 @@ class AbstractDomain(ABC):
@property
def model(self):
- """The gmsh model used within the domain.
+ """The (gmsh) geometrical model used to define the domain.
Returns
-------
:class:`Model<gmsh.model>`
- The gmsh model used to tile the domain.
+ A purely geometrical description of the domain.
"""
return AbstractDomain._model
def synchronize(self):
- """synchronizes the element factory ????
+ """synchronizes the CAD model with the gmsh model.
Returns
-------
@@ -286,10 +288,7 @@ class AbstractDomain(ABC):
Parameters
----------
- *args : type
- Description of parameter `*args`.
- **kwargs : type
- Description of parameter `**kwargs`.
+ None
Returns
-------
@@ -303,7 +302,7 @@ class AbstractDomain(ABC):
"""
def run_gui(self):
- """???
+ """Generates a quick gmsh gui with the current model in it.
Returns
-------
@@ -351,16 +350,29 @@ class AbstractDomain(ABC):
gmsh.option.setNumber("General.Verbosity", 0)
def generate(self):
- """Generates the data structure and writes it on disc.
+ """Generates a fenics mesh representing the domain.
Yields
------
- :class:`Mesh<gmsh.mesh>` ???
- Tile mesh corresponding to the tiled domain.
- :class:'ds<fenics.mesh.MeshFunction>' ??
- The boundary surface incremental element.
- :class:'dx<fenics.mesh.MeshFunction>' ??
- The volume incremental element. ???
+ 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 temporary on disc.
+ * MeshFunctionSizet : `Sizet` means that the output of
+ the function is a positive integer.
"""
if self._cell_type == 'line':
@@ -380,12 +392,18 @@ class AbstractDomain(ABC):
self.dx = fe.Measure('dx', domain=self.mesh, subdomain_data=self.cdata)
def _mesh_conversion(self):
- """Short summary.
+ """Converts a temporary gmsh mesh into a fenics (dolphin) one.
Returns
-------
- type
- Description of returned object.
+ 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.
"""
path = self._dir_name + '/mesh.msh'
@@ -443,7 +461,8 @@ class AbstractDomain(ABC):
line_data =\
mesh_io.cell_data_dict["gmsh:physical"][key]
else:
- line_data = np.vstack([line_data, mesh_io.cell_data_dict["gmsh:physical"][key]])
+ line_data = np.vstack([line_data,
+ mesh_io.cell_data_dict["gmsh:physical"][key]])
line_mesh = meshio.Mesh(points=mesh_io.points,
cells=[("line", line_cells)],
@@ -460,12 +479,62 @@ class AbstractDomain(ABC):
class BuiltInModel(object):
+ """A model based on the Gmsh `geo` factory.
+
+ Parameters
+ ----------
+ None
+
+ Attributes
+ ----------
+ factory : :class:`occ<gmsh.model.geo>`
+ A gmsh model `geo` factory to generate mesh from geometry.
+
+ """
+
def __init__(self, *args, **kwargs):
+ """Sets the attribute factory to a specific type of gmsh factory.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ factory : :class:`occ<gmsh.model.geo>`
+ A gmsh model `geo` factory to generate mesh from geometry.
+
+ """
super().__init__()
self.factory = gmsh.model.geo
class OccModel(object):
+ """A model based on the Gmsh `Occ` factory.
+
+ Parameters
+ ----------
+ None
+
+ Attributes
+ ----------
+ factory : :class:`occ<gmsh.model.occ>`
+ A gmsh model `occ` factory to generate mesh from geometry.
+
+ """
+
def __init__(self, *args, **kwargs):
+ """Sets the attribute factory to a specific type of gmsh factory.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ factory : :class:`occ<gmsh.model.occ>`
+ A gmsh model `occ` factory to generate mesh from geometry.
+
+ """
super().__init__()
self.factory = gmsh.model.occ
diff --git a/src/bvpy/domains/geometry.py b/src/bvpy/domains/geometry.py
index 174de31a60d9fe5e14bd8752c315e6fd558c0952..6fc04069176520d560275856d245e8a838a8ba0f 100644
--- a/src/bvpy/domains/geometry.py
+++ b/src/bvpy/domains/geometry.py
@@ -22,7 +22,6 @@
import fenics as fe
import numpy as np
-from math import pi
class Init_orientation_3D(fe.UserExpression):
@@ -30,20 +29,37 @@ class Init_orientation_3D(fe.UserExpression):
Attributes
----------
- mesh : :class:`dolfin.cpp.mesh.Mesh<dolfin.cpp.mesh.Mesh>`
+ mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
Mesh to work with.
- interior : list(float)
+ barycenter : numpy.array of floats
+ The (3D) position vector of the barycenter of all the mesh vertices.
+ interior : list of floats
A 3D vector (x, y, z) to orient properly the mesh.
Methods
-------
- eval(value, x)
- Compute the relative position with respect to the barycenter.
- value_shape()
- Return the geometrical dimesion of the structure (always 3).
+ eval
+ Computes the relative position with respect to the barycenter.
+ value_shape
+ Returns the geometrical dimesion of the structure (always 3).
"""
def __init__(self, mesh, interior, **kwargs):
+ """Short summary.
+
+ Parameters
+ ----------
+ mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
+ Mesh to work with.
+ interior : list of floats
+ A 3D vector (x, y, z) to orient properly the mesh.
+
+ Returns
+ -------
+ barycenter : numpy.array of floats
+ The (3D) position vector of the barycenter of the mesh vertices.
+
+ """
if interior is None:
self.barycenter = mesh.coordinates().mean(axis=0)
else:
@@ -51,11 +67,38 @@ class Init_orientation_3D(fe.UserExpression):
super().__init__(**kwargs)
def eval(self, value, x):
+ """Computes the relative position with respect to the barycenter.
+
+ Parameters
+ ----------
+ value : numpy.array of float
+ the position we want to update (???).
+ x : list of float
+ The position vector to re-center.
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ This method is used by fenics and is not explicitly used within
+ the bvpy code but is nonetheless needed.
+
+ """
value[0] = x[0] - self.barycenter[0]
value[1] = x[1] - self.barycenter[1]
value[2] = x[2] - self.barycenter[2]
def value_shape(self):
+ """Returns the geometrical dimesion of the structure (always 3).
+
+ Returns
+ -------
+ tuple of int
+ This always equals 3.
+
+ """
return (3,)
@@ -64,13 +107,17 @@ def normal(domain, scale=1, interior=None, degree=1, name='normal'):
Parameters
----------
- domain : Domain
+ domain : subclass of :class:`Domain<bvpy.domains.AbstracDomain>`
The meshed surface on which we want to compute the normal field.
- scale : type
- Optional (Default : 1). The displacement intensity to consider.
- interior : list(float)
- Optional (Default : None). A 3D vector (x, y, z)
- to orient properly the mesh.
+ scale : float, optional
+ The displacement intensity to consider (the default is 1).
+ interior : list(float), optional
+ A 3D vector to set the mesh orientation (the default is None).
+ degree : int, optional
+ The interpolation degree to use to compute the vector space
+ on the border (the default is 1).
+ name : str, optional
+ the label we want to attach to the normal field.
Returns
-------
@@ -92,7 +139,7 @@ def normal(domain, scale=1, interior=None, degree=1, name='normal'):
V_scal = fe.FunctionSpace(mesh, 'CG', 1)
vf = np.zeros([mesh.num_entities(0), 3])
- mesh.init(0,2)
+ mesh.init(0, 2)
for v in fe.vertices(mesh):
f_vertex = np.zeros(3)
@@ -121,29 +168,75 @@ def normal(domain, scale=1, interior=None, degree=1, name='normal'):
return u
def get_boundary_facets(mesh):
+ """Selects the boundary elements of a given mesh.
+
+ Parameters
+ ----------
+ mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
+ The mesh from which we want to get the boundary elements.
+
+ Returns
+ -------
+ facets : list of :class:`facets<dolfin.cpp.mesh.facets>`
+ the list of the outermost elements of the considered mesh.
+
+ Notes
+ -----
+ The returned elements corresponds to mesh elements of topological
+ dimension (n-1) for a mesh of topological dimension n.
+
+ """
mesh.init()
facets = list()
for fct in fe.facets(mesh):
- if len(list(fe.cells(fct)))==1:
+ if len(list(fe.cells(fct))) == 1:
facets.append((fct))
+
return facets
def get_boundary_vertex(mesh):
+ """Short summary.
+
+ Parameters
+ ----------
+ mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
+ The mesh from which we want to get the boundary vertices.
+
+
+ Returns
+ -------
+ list of :class:`Vertex<dofin.cpp.mesh.Vertex>`
+ The list of the outermost vertices of the considered mesh.
+
+ """
bnd_facets = get_boundary_facets(mesh)
bnd_vertex = list()
for fct in bnd_facets:
[bnd_vertex.append(v.index()) for v in fe.vertices(fct)]
- bnd_vertex = list(set(bnd_vertex)) #Unique
+ bnd_vertex = list(set(bnd_vertex)) # Unique
bnd_vertex = [fe.Vertex(mesh, i) for i in bnd_vertex]
return bnd_vertex
def boundary_normal(mesh):
+ """Computes the normal vector field on the boundary of a mesh.
+
+ Parameters
+ ----------
+ mesh : :class:`Mesh<dolfin.cpp.mesh.Mesh>`
+ The mesh on which the boundary normals should be computed.
+
+ Returns
+ -------
+ :class:`Function<dolfin.functions.function.Function>`
+ a list of 3D vector defined on the outermost vertices of a mesh.
+
+ """
mesh.init()
bnd_facets = get_boundary_facets(mesh)
facets_index = [f.index() for f in bnd_facets]
@@ -151,9 +244,11 @@ def boundary_normal(mesh):
for fct in bnd_facets:
for vtx in fe.vertices(fct):
- norm = np.array([f.normal().array() for f in fe.facets(vtx) if f.index() in facets_index] )
+ norm = np.array([f.normal().array()
+ for f in fe.facets(vtx)
+ if f.index() in facets_index])
norm = np.sum(norm, axis=0)
- norm = norm/np.linalg.norm(norm)
+ norm = norm / np.linalg.norm(norm)
normals[vtx.index()] = norm
V = fe.VectorFunctionSpace(mesh, 'P', 1, dim=3)
diff --git a/src/bvpy/ivp.py b/src/bvpy/ivp.py
index 10008da76b51e973d3eff6006c3186d1c26a5e3b..491d33fd42e0056430afc8144c81d6b10c921798 100644
--- a/src/bvpy/ivp.py
+++ b/src/bvpy/ivp.py
@@ -33,22 +33,44 @@ class Time():
class Scheme():
- def __init__(self, vform):
+ def __init__(self, vform=None, dt=None):
self.vform = vform
+ self.previous_solution = None
+ self.dt = dt
- def implicit_euler(lhs, rhs):
- return lhs, Zero()
+ def set_previous_solution(self, previous):
+ self.previous_solution = previous
+
+ def implicit_euler(self, u, v, sol):
+ lhs, rhs = self.vform.get_form(u, v, sol)
+
+ lhs = fe.Constant(self.dt)*lhs
+ rhs = fe.Constant(self.dt)*rhs
+
+ lhs += u*v*self.dx
+ rhs += self.previous_solution*v*self.dx
+
+ return lhs, rhs
def explicit_euler(lhs, rhs):
return lhs, rhs
+ @property
+ def dx(self):
+ return self.vform.dx
+
+ @property
+ def ds(self):
+ return self.vform.ds
+
class IVP(BVP):
- def __init__(self, domain=None, vform=None, bc=None):
+ def __init__(self, domain=None, vform=None, bc=None,):
super().__init__(domain=domain, vform=vform, bc=bc)
self.time = Time()
+ self.scheme = Scheme()
def set_time(self, time):
self.time.current = time
diff --git a/tutorials/bvpy_tutorial_3_poissonEquation_FEniCSlike.ipynb b/tutorials/bvpy_tutorial_3_poissonEquation_FEniCSlike.ipynb
index e49cb576008e8dd6c22c38c15d98be8c72b7c48d..25a27622b9fb10ffa9daff86f86d836aeff8c60a 100644
--- a/tutorials/bvpy_tutorial_3_poissonEquation_FEniCSlike.ipynb
+++ b/tutorials/bvpy_tutorial_3_poissonEquation_FEniCSlike.ipynb
@@ -223,4 +223,4 @@
},
"nbformat": 4,
"nbformat_minor": 4
-}
\ No newline at end of file
+}
diff --git a/tutorials/bvpy_tutorial_5_linearElasticity_2DFlatCellularizedTissue.ipynb b/tutorials/bvpy_tutorial_5_linearElasticity_2DFlatCellularizedTissue.ipynb
index fc332e6faeb98e7577f7165222b02f4f508ade02..d9f095d27dbccfd06632501acc86363d26262dae 100644
--- a/tutorials/bvpy_tutorial_5_linearElasticity_2DFlatCellularizedTissue.ipynb
+++ b/tutorials/bvpy_tutorial_5_linearElasticity_2DFlatCellularizedTissue.ipynb
@@ -62,6 +62,18 @@
{
"cell_type": "code",
"execution_count": null,
+<<<<<<< HEAD
+=======
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "type(drchlt)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+>>>>>>> develop
"metadata": {},
"outputs": [],
"source": [
@@ -75,6 +87,18 @@
{
"cell_type": "code",
"execution_count": null,
+<<<<<<< HEAD
+=======
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p_force.vector().get_local().reshape(-1,3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+>>>>>>> develop
"metadata": {},
"outputs": [],
"source": [