diff --git a/requirements-dev.txt b/requirements-dev.txt
index 385dbdab..8dac3c99 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,8 +1,8 @@
-black == 23.*
+black == 24.*
deepdiff
flake8
isort
-mypy >= 0.981
+mypy
pytest >= 4.6
pytest-cov
pytest-runner
diff --git a/src/pygeon/__init__.py b/src/pygeon/__init__.py
index caca007a..8347e6cc 100644
--- a/src/pygeon/__init__.py
+++ b/src/pygeon/__init__.py
@@ -18,9 +18,14 @@
from pygeon.discretizations.vec_discretization import VecDiscretization
from pygeon.discretizations.fem.hcurl import Nedelec0, Nedelec1
-from pygeon.discretizations.fem.hdiv import RT0, BDM1, VecBDM1
+from pygeon.discretizations.fem.hdiv import RT0, BDM1, VecBDM1, VecRT0
from pygeon.discretizations.fem.h1 import Lagrange1, VecLagrange1
-from pygeon.discretizations.fem.l2 import PwConstants, PwLinears, VecPwConstants
+from pygeon.discretizations.fem.l2 import (
+ PwConstants,
+ PwLinears,
+ VecPwConstants,
+ VecPwLinears,
+)
from pygeon.discretizations.vem.hdiv import VRT0, VBDM1
from pygeon.discretizations.vem.h1 import VLagrange1, VecVLagrange1
@@ -44,6 +49,7 @@
SpanningTree,
SpanningWeightedTrees,
SpanningTreeElasticity,
+ SpanningTreeCosserat,
)
import pygeon.utils.bmat as bmat
diff --git a/src/pygeon/discretizations/discretization.py b/src/pygeon/discretizations/discretization.py
index 1fb87d1e..2da02f22 100644
--- a/src/pygeon/discretizations/discretization.py
+++ b/src/pygeon/discretizations/discretization.py
@@ -265,6 +265,7 @@ def error_l2(
ana_sol: Callable[[np.ndarray], np.ndarray],
relative: Optional[bool] = True,
etype: Optional[str] = "standard",
+ data: dict = None,
) -> float:
"""
Returns the l2 error computed against an analytical solution given as a function.
@@ -281,7 +282,7 @@ def error_l2(
float: The computed error.
"""
int_sol = self.interpolate(sd, ana_sol)
- mass = self.assemble_mass_matrix(sd)
+ mass = self.assemble_mass_matrix(sd, data)
norm = (int_sol @ mass @ int_sol.T) if relative else 1
diff --git a/src/pygeon/discretizations/fem/hdiv.py b/src/pygeon/discretizations/fem/hdiv.py
index fe6ca7d7..7d5c2bf9 100644
--- a/src/pygeon/discretizations/fem/hdiv.py
+++ b/src/pygeon/discretizations/fem/hdiv.py
@@ -41,7 +41,7 @@ class RT0(pg.Discretization, pp.RT0):
assemble_nat_bc(sd: pg.Grid, func: Callable[[np.ndarray], np.ndarray],
b_faces: np.ndarray) -> np.ndarray:
- Assembles the natural boundary condition term (n dot q, func)_\Gamma
+ Assembles the natural boundary condition term (n dot q, func)_Gamma
get_range_discr_class(dim: int) -> pg.Discretization:
Returns the range discretization class for the given dimension.
@@ -121,8 +121,9 @@ def assemble_mass_matrix(
Returns:
sps.csc_matrix: The mass matrix.
"""
-
+ # create dummy data, unitary permeability, in case not present
data = self.create_dummy_data(sd, data)
+ # perform the rt0 discretization
pp.RT0.discretize(self, sd, data)
return data[pp.DISCRETIZATION_MATRICES][self.keyword][
self.mass_matrix_key
@@ -211,7 +212,7 @@ def assemble_nat_bc(
) -> np.ndarray:
"""
Assembles the natural boundary condition term
- (n dot q, func)_\Gamma
+ (n dot q, func)_Gamma
Args:
sd (pg.Grid): The grid object representing the computational domain.
@@ -584,7 +585,7 @@ def assemble_nat_bc(
) -> np.ndarray:
"""
Assembles the natural boundary condition term
- (n dot q, func)_\Gamma
+ (n dot q, func)_Gamma
Args:
sd (pg.Grid): The grid object representing the computational domain.
@@ -606,7 +607,7 @@ def assemble_nat_bc(
sign = np.sum(sd.cell_faces.tocsr()[face, :])
loc_vals = np.array(
[func(sd.nodes[:, node]) for node in sd.face_nodes[:, face].indices]
- )
+ ).ravel()
vals[face + np.arange(sd.dim) * sd.num_faces] = sign * local_mass @ loc_vals
@@ -812,7 +813,7 @@ def assemble_mass_matrix(self, sd: pg.Grid, data: dict) -> sps.csc_matrix:
# Assemble the trace part
B = self.assemble_trace_matrix(sd)
- # Assemble the piecewise linear mass mastrix, to assemble the term
+ # Assemble the piecewise linear mass matrix, to assemble the term
# (Trace(sigma), Trace(tau))
discr = pg.PwLinears(self.keyword)
M = discr.assemble_mass_matrix(sd, data_for_PwL)
@@ -820,6 +821,46 @@ def assemble_mass_matrix(self, sd: pg.Grid, data: dict) -> sps.csc_matrix:
# Compose all the parts and return them
return D - B.T @ M @ B
+ def assemble_mass_matrix_cosserat(self, sd: pg.Grid, data: dict) -> sps.csc_matrix:
+ """
+ Assembles and returns the mass matrix for vector BDM1 discretizing the Cosserat inner
+ product, which is given by (A sigma, tau) where
+ A sigma = (sym(sigma) - coeff * Trace(sigma) * I) / (2 mu) + skw(sigma) / (2 mu_c)
+ with mu and lambda the Lamé constants, coeff = lambda / (2*mu + dim*lambda), and
+ mu_c the coupling Lamé modulus.
+
+ Args:
+ sd (pg.Grid): The grid.
+ data (dict): Data for the assembly.
+
+ Returns:
+ sps.csc_matrix: The mass matrix obtained from the discretization.
+ """
+ M = self.assemble_mass_matrix(sd, data)
+
+ # Extract the data
+ mu = data[pp.PARAMETERS][self.keyword]["mu"]
+ mu_c = data[pp.PARAMETERS][self.keyword]["mu_c"]
+
+ coeff = 0.25 * (1 / mu_c - 1 / mu)
+
+ # If coeff is a scalar, replace it by a vector so that it can be accessed per cell
+ if isinstance(coeff, np.ScalarType):
+ coeff = np.full(sd.num_cells, coeff)
+
+ data_for_R = pp.initialize_data(sd, {}, self.keyword, {"weight": coeff})
+
+ if sd.dim == 2:
+ R_space = pg.PwLinears(self.keyword)
+ elif sd.dim == 3:
+ R_space = pg.VecPwLinears(self.keyword)
+
+ R_mass = R_space.assemble_mass_matrix(sd, data_for_R)
+
+ asym = self.assemble_asym_matrix(sd, False)
+
+ return M + asym.T @ R_mass @ asym
+
def assemble_trace_matrix(self, sd: pg.Grid) -> sps.csc_matrix:
"""
Assembles and returns the trace matrix for the vector BDM1.
@@ -876,69 +917,200 @@ def assemble_trace_matrix(self, sd: pg.Grid) -> sps.csc_matrix:
# Construct the global matrices
return sps.csc_matrix((data_IJ[:idx], (rows_I[:idx], cols_J[:idx])))
- def assemble_asym_matrix(self, sd: pg.Grid) -> sps.csc_matrix:
+ def assemble_asym_matrix(self, sd: pg.Grid, as_pwconstant=True) -> sps.csc_matrix:
"""
Assembles and returns the asymmetric matrix for the vector BDM1.
The asymmetric operator `as' for a tensor is a scalar and it is defined in 2d as
- as(tau) = tau_xy - tau_yx
+ as(tau) = tau_yx - tau_xy
while for a tensor in 3d it is a vector and given by
as(tau) = [tau_zy - tau_yz, tau_xz - tau_zx, tau_yx - tau_xy]^T
- Note: We assume that the as(tau) is a cell variable.
+ Note: We assume that the as(tau) is a piecewise linear.
Args:
sd (pg.Grid): The grid.
+ as_pwconstant (bool): Compute the operator with the range on the piece-wise
+ constant (default), otherwise the mapping is on the piece-wise linears.
Returns:
sps.csc_matrix: The asymmetric matrix obtained from the discretization.
"""
- # We need to map the BDM1 to the cell center, thus we construct the proper operator
- P = self.eval_at_cell_centers(sd)
- nc, cv = sd.num_cells, sd.cell_volumes
- # We consider a different approach if sd is 2d or 3d
- if sd.dim == 2:
- # Every cell has two entries, one for the t_xy and one for the t_yx
- rows_I = np.tile(np.arange(nc), 2)
- # Given a 2d tensor, represented as a vector, the t_xy component is the second
- # and the t_yx is the third
- cols_J = np.hstack(
- (
- np.arange(nc, 2 * nc),
- np.arange(3 * nc, 4 * nc),
- )
+
+ # overestimate the size
+ size = np.square((sd.dim + 1) * sd.dim) * sd.num_cells
+ rows_I = np.empty(size, dtype=int)
+ cols_J = np.empty(size, dtype=int)
+ data_IJ = np.empty(size)
+ idx = 0
+
+ # Helper functions for inside the loop
+ negate_col = [2, 0, 1]
+ zeroed_col = [0, 1, 2]
+
+ if sd.dim == 3:
+ ind_list = np.arange(3)
+ shift = ind_list
+ rot_space = pg.VecPwLinears(self.keyword)
+ scaling = sps.diags(np.tile(sd.cell_volumes, 3))
+ elif sd.dim == 2:
+ ind_list = [2]
+ shift = [0, 0, 0]
+ rot_space = pg.PwLinears(self.keyword)
+ scaling = sps.diags(sd.cell_volumes)
+ else:
+ raise ValueError("The grid should be either two or three-dimensional")
+
+ cell_nodes = sd.cell_nodes()
+ for c in np.arange(sd.num_cells):
+ # For the current cell retrieve its faces and
+ # determine the location of the dof
+ loc = slice(sd.cell_faces.indptr[c], sd.cell_faces.indptr[c + 1])
+ faces_loc = sd.cell_faces.indices[loc]
+ dof_loc = np.reshape(
+ sd.face_nodes[:, faces_loc].indices, (sd.dim, -1), order="F"
)
- # t_xy gets a +1 and t_yx a -1, represented as p0 dofs so with the cell volume
- data_IJ = np.hstack((cv, -cv))
- # Assemble the matrix with a given shape, t_xx and t_yy are not considered
- T = sps.csc_matrix((data_IJ, (rows_I, cols_J)), shape=(nc, P.shape[0]))
- return T @ P
- elif sd.dim == 3:
- # Define an help function
- enum = lambda i: np.arange(i * nc, (i + 1) * nc)
- # Since as(tau) is a piecewise vector we need to arrange the entries accordingly
- # to the second component of the tensor
- # as(tau)_x = tau_zy - tau_yz -> z and y with -1 and +1
- # as(tau)_y = tau_xz - tau_zx -> z and x with +1 and -1
- # as(tau)_z = tau_yx - tau_xy -> x and y with -1 and +1
- rows_I = np.hstack((enum(2), enum(1), enum(2), enum(0), enum(1), enum(0)))
- # Given a 3d tensor, represented as a vector, the involved components are
- # (t_xy, t_xz, t_yx, t_yz, t_zx, t_zy) which are in positions (1, 2, 3) and
- # (5, 6, 7)
- cols_J = np.hstack(
- (
- np.arange(nc, 4 * nc),
- np.arange(5 * nc, 8 * nc),
- )
+
+ cell_nodes_loc = cell_nodes[:, c].toarray()
+ Psi = self.scalar_discr.eval_basis_at_node(
+ sd, dof_loc, cell_nodes_loc, faces_loc
)
- # Assording to the ordering of rows and cols given abose we have the following
- # signs (-1, +1, +1, -1, -1, +1), weighted by the cell volumes since it is
- # represented as a p0 vector element
- data_IJ = np.hstack((-cv, cv, cv, -cv, -cv, cv))
- T = sps.csc_matrix((data_IJ, (rows_I, cols_J)), shape=(3 * nc, P.shape[0]))
- return T @ P
+
+ # Get all the components of the basis at node
+ Psi_i, Psi_j, Psi_v = sps.find(Psi)
+
+ for ind in ind_list:
+ Psi_v_copy = Psi_v.copy()
+ Psi_v_copy[np.mod(Psi_j, 3) == negate_col[ind]] *= -1
+ Psi_v_copy[np.mod(Psi_j, 3) == zeroed_col[ind]] *= 0
+
+ loc_ind = np.hstack([faces_loc] * sd.dim)
+ loc_ind += np.repeat(np.arange(sd.dim), sd.dim + 1) * sd.num_faces
+
+ cols = np.tile(loc_ind, (3, 1))
+ cols[0, :] += np.mod(-ind, 3) * self.scalar_discr.ndof(sd)
+ cols[1, :] += np.mod(-ind - 1, 3) * self.scalar_discr.ndof(sd)
+ cols[2, :] += np.mod(-ind - 2, 3) * self.scalar_discr.ndof(sd)
+
+ cols = np.tile(cols, (sd.dim + 1, 1)).T
+
+ cols = cols[Psi_i, Psi_j]
+
+ nodes_loc = np.arange((sd.dim + 1) * c, (sd.dim + 1) * (c + 1))
+ rows = np.repeat(nodes_loc, 3)[Psi_j]
+
+ # Save values of the local matrix in the global structure
+ loc_idx = slice(idx, idx + cols.size)
+ rows_I[loc_idx] = rows + shift[ind] * (sd.dim + 1) * sd.num_cells
+ cols_J[loc_idx] = cols
+ data_IJ[loc_idx] = Psi_v_copy
+ idx += cols.size
+
+ # Construct the global matrices
+ asym = sps.csc_matrix((data_IJ[:idx], (rows_I[:idx], cols_J[:idx])))
+
+ # Return the operator that maps to the piece-wise constant
+ if as_pwconstant:
+ return scaling @ rot_space.eval_at_cell_centers(sd) @ asym
else:
- raise ValueError("The grid should be either bi or three-dimensional")
+ return asym
+
+ def assemble_lumped_matrix(
+ self, sd: pg.Grid, data: Optional[dict] = None
+ ) -> sps.csc_matrix:
+ """
+ Assembles the lumped matrix for the given grid.
+
+ Args:
+ sd (pg.Grid): The grid object.
+ data (Optional[dict]): Optional data dictionary.
+
+ Returns:
+ sps.csc_matrix: The assembled lumped matrix.
+ """
+ # Assemble the block diagonal mass matrix for the base discretization class
+ D = super().assemble_lumped_matrix(sd)
+
+ # Assemble the trace part
+ B = self.assemble_trace_matrix(sd)
+
+ # Assemble the piecewise linear mass matrix, to assemble the term
+ # (Trace(sigma), Trace(tau))
+ discr = pg.PwLinears(self.keyword)
+ M = discr.assemble_lumped_matrix(sd)
+
+ # Extract the data and compute the coefficient for the trace part
+ mu = data[pp.PARAMETERS][self.keyword]["mu"]
+ lambda_ = data[pp.PARAMETERS][self.keyword]["lambda"]
+ coeff = lambda_ / (2 * mu + sd.dim * lambda_)
+
+ # Compose all the parts and return them
+ return (D - coeff * B.T @ M @ B) / (2 * mu)
+
+ def assemble_lumped_matrix_cosserat(
+ self, sd: pg.Grid, data: dict
+ ) -> sps.csc_matrix:
+ """
+ Assembles the lumped matrix with cosserat terms for the given grid.
+
+ Args:
+ sd (pg.Grid): The grid object.
+ data (Optional[dict]): Optional data dictionary.
+
+ Returns:
+ sps.csc_matrix: The assembled lumped matrix.
+ """
+
+ M = self.assemble_lumped_matrix(sd, data)
+
+ # Extract the data
+ mu = data[pp.PARAMETERS][self.keyword]["mu"]
+ mu_c = data[pp.PARAMETERS][self.keyword]["mu_c"]
+
+ coeff = 0.25 * (1 / mu_c - 1 / mu)
+
+ # If coeff is a scalar, replace it by a vector so that it can be accessed per cell
+ if isinstance(coeff, np.ScalarType):
+ coeff = np.full(sd.num_cells, coeff)
+
+ data_for_R = pp.initialize_data(sd, {}, self.keyword, {"weight": coeff})
+
+ if sd.dim == 2:
+ R_space = pg.PwLinears(self.keyword)
+ elif sd.dim == 3:
+ R_space = pg.VecPwLinears(self.keyword)
+
+ R_mass = R_space.assemble_lumped_matrix(sd, data_for_R)
+
+ asym = self.assemble_asym_matrix(sd, False)
+
+ return M + asym.T @ R_mass @ asym
+
+ def proj_to_RT0(self, sd: pg.Grid) -> sps.csc_matrix:
+ """
+ Project the function space to the lowest order Raviart-Thomas (RT0) space.
+
+ Args:
+ sd (pg.Grid): The grid object representing the computational domain.
+
+ Returns:
+ sps.csc_matrix: The projection matrix to the RT0 space.
+ """
+ proj = self.scalar_discr.proj_to_RT0(sd)
+ return sps.block_diag([proj] * sd.dim, format="csc")
+
+ def proj_from_RT0(self, sd: pg.Grid) -> sps.csc_matrix:
+ """
+ Project the RT0 finite element space onto the faces of the given grid.
+
+ Args:
+ sd (pg.Grid): The grid on which the projection is performed.
+
+ Returns:
+ sps.csc_matrix: The projection matrix.
+ """
+ proj = self.scalar_discr.proj_from_RT0(sd)
+ return sps.block_diag([proj] * sd.dim, format="csc")
def get_range_discr_class(self, dim: int) -> object:
"""
@@ -952,3 +1124,163 @@ def get_range_discr_class(self, dim: int) -> object:
differential
"""
return pg.VecPwConstants
+
+
+class VecRT0(pg.VecDiscretization):
+ def __init__(self, keyword: str) -> None:
+ """
+ Initialize the vector RT0 discretization class.
+ The scalar discretization class is pg.RT0.
+
+ We are considering the following structure of the stress tensor in 2d
+
+ sigma = [[sigma_xx, sigma_xy],
+ [sigma_yx, sigma_yy]]
+
+ which is represented in the code unrolled row-wise as a vector of length 4
+
+ sigma = [sigma_xx, sigma_xy,
+ sigma_yx, sigma_yy]
+
+ While in 3d the stress tensor can be written as
+
+ sigma = [[sigma_xx, sigma_xy, sigma_xz],
+ [sigma_yx, sigma_yy, sigma_yz],
+ [sigma_zx, sigma_zy, sigma_zz]]
+
+ where its vectorized structure of length 9 is given by
+
+ sigma = [sigma_xx, sigma_xy, sigma_xz,
+ sigma_yx, sigma_yy, sigma_yz,
+ sigma_zx, sigma_zy, sigma_zz]
+
+ Args:
+ keyword (str): The keyword for the vector discretization class.
+
+ Returns:
+ None
+ """
+ super().__init__(keyword, pg.RT0)
+
+ def assemble_mass_matrix(self, sd: pg.Grid, data: dict) -> sps.csc_matrix:
+ """
+ Assembles and returns the mass matrix for vector RT0, which is given by
+ (A sigma, tau) where A sigma = (sigma - coeff * Trace(sigma) * I) / (2 mu)
+ with mu and lambda the Lamé constants and coeff = lambda / (2*mu + dim*lambda)
+
+ Args:
+ sd (pg.Grid): The grid.
+ data (dict): Data for the assembly.
+
+ Returns:
+ sps.csc_matrix: The mass matrix obtained from the discretization.
+ """
+ # Assemble the block diagonal mass matrix for the base discretization class,
+ # with unitary data. The data are handled afterwards
+ D = super().assemble_mass_matrix(sd)
+
+ # Assemble the trace part
+ B = self.assemble_trace_matrix(sd)
+
+ # Assemble the piecewise linear mass matrix, to assemble the term
+ # (Trace(sigma), Trace(tau))
+ discr = pg.PwLinears(self.keyword)
+ M = discr.assemble_mass_matrix(sd)
+
+ # Extract the data and compute the coefficient for the trace part
+ mu = data[pp.PARAMETERS][self.keyword]["mu"]
+ lambda_ = data[pp.PARAMETERS][self.keyword]["lambda"]
+ coeff = lambda_ / (2 * mu + sd.dim * lambda_)
+
+ # Compose all the parts and return them
+ return (D - coeff * B.T @ M @ B) / (2 * mu)
+
+ def assemble_mass_matrix_cosserat(self, sd: pg.Grid, data: dict) -> sps.csc_matrix:
+ """
+ Assembles and returns the mass matrix for vector BDM1 discretizing the Cosserat inner
+ product, which is given by (A sigma, tau) where
+ A sigma = (sym(sigma) - coeff * Trace(sigma) * I) / (2 mu) + skw(sigma) / (2 mu_c)
+ with mu and lambda the Lamé constants, coeff = lambda / (2*mu + dim*lambda), and
+ mu_c the coupling Lamé modulus.
+
+ Args:
+ sd (pg.Grid): The grid.
+ data (dict): Data for the assembly.
+
+ Returns:
+ sps.csc_matrix: The mass matrix obtained from the discretization.
+
+ TODO: Consider using inheritance from VecBDM1.assemble_mass_matrix_cosserat
+ """
+
+ M = self.assemble_mass_matrix(sd, data)
+
+ # Extract the data
+ mu = data[pp.PARAMETERS][self.keyword]["mu"]
+ mu_c = data[pp.PARAMETERS][self.keyword]["mu_c"]
+
+ coeff = 0.25 * (1 / mu_c - 1 / mu)
+
+ # If coeff is a scalar, replace it by a vector so that it can be accessed per cell
+ if isinstance(coeff, np.ScalarType):
+ coeff = np.full(sd.num_cells, coeff)
+
+ data_for_R = pp.initialize_data(sd, {}, self.keyword, {"weight": coeff})
+
+ if sd.dim == 2:
+ R_space = pg.PwLinears(self.keyword)
+ elif sd.dim == 3:
+ R_space = pg.VecPwLinears(self.keyword)
+
+ R_mass = R_space.assemble_mass_matrix(sd, data_for_R)
+
+ asym = self.assemble_asym_matrix(sd, False)
+
+ return M + asym.T @ R_mass @ asym
+
+ def assemble_trace_matrix(self, sd: pg.Grid) -> sps.csc_matrix:
+ """
+ Assembles and returns the trace matrix for the vector BDM1.
+
+ Args:
+ sd (pg.Grid): The grid.
+
+ Returns:
+ sps.csc_matrix: The trace matrix obtained from the discretization.
+ """
+ vec_bdm1 = VecBDM1(self.keyword)
+ proj = vec_bdm1.proj_from_RT0(sd)
+ return vec_bdm1.assemble_trace_matrix(sd) @ proj
+
+ def assemble_asym_matrix(self, sd: pg.Grid, as_pwconstant=True) -> sps.csc_matrix:
+ """
+ Assembles and returns the asymmetric matrix for the vector RT0.
+
+ The asymmetric operator `as' for a tensor is a scalar and it is defined in 2d as
+ as(tau) = tau_xy - tau_yx
+ while for a tensor in 3d it is a vector and given by
+ as(tau) = [tau_zy - tau_yz, tau_xz - tau_zx, tau_yx - tau_xy]^T
+
+ Note: We assume that the as(tau) is a cell variable.
+
+ Args:
+ sd (pg.Grid): The grid.
+
+ Returns:
+ sps.csc_matrix: The asymmetric matrix obtained from the discretization.
+ """
+ vec_bdm1 = VecBDM1(self.keyword)
+ proj = vec_bdm1.proj_from_RT0(sd)
+ return vec_bdm1.assemble_asym_matrix(sd, as_pwconstant) @ proj
+
+ def get_range_discr_class(self, dim: int) -> object:
+ """
+ Returns the range discretization class for the given dimension.
+
+ Args:
+ dim (int): The dimension of the range space.
+
+ Returns:
+ pg.Discretization: The range discretization class.
+ """
+ return pg.VecPwConstants
diff --git a/src/pygeon/discretizations/fem/l2.py b/src/pygeon/discretizations/fem/l2.py
index bf851c9f..ced9a6f2 100644
--- a/src/pygeon/discretizations/fem/l2.py
+++ b/src/pygeon/discretizations/fem/l2.py
@@ -341,6 +341,27 @@ def assemble_mass_matrix(
# Construct the global matrix
return sps.csc_matrix((data_IJ, (rows_I, cols_J)))
+ def assemble_lumped_matrix(
+ self, sd: pg.Grid, data: Optional[dict] = None
+ ) -> sps.csc_matrix:
+ """
+ Assembles the lumped matrix for the given grid.
+
+ Args:
+ sd (pg.Grid): The grid object.
+ data (Optional[dict]): Optional data dictionary.
+
+ Returns:
+ sps.csc_matrix: The assembled lumped matrix.
+ """
+ try:
+ weight = data[pp.PARAMETERS][self.keyword]["weight"]
+ except Exception:
+ weight = np.ones(sd.num_cells)
+
+ diag = np.repeat(weight * sd.cell_volumes, sd.dim + 1) / (sd.dim + 1)
+ return sps.diags(diag, format="csc")
+
def assemble_diff_matrix(self, sd: pg.Grid) -> sps.csc_matrix:
"""
Assembles the matrix corresponding to the differential operator.
@@ -393,7 +414,12 @@ def eval_at_cell_centers(self, sd: pg.Grid) -> np.ndarray:
Returns:
np.ndarray: The evaluation matrix.
"""
- raise NotImplementedError
+
+ rows = np.repeat(np.arange(sd.num_cells), sd.dim + 1)
+ cols = np.arange(self.ndof(sd))
+ data = np.ones(self.ndof(sd)) / (sd.dim + 1)
+
+ return sps.csc_matrix((data, (rows, cols)))
def interpolate(
self, sd: pg.Grid, func: Callable[[np.ndarray], np.ndarray]
@@ -468,3 +494,95 @@ def assemble_nat_bc(
np.ndarray: The assembled natural boundary condition vector.
"""
return np.zeros(self.ndof(sd))
+
+ def error_l2(
+ self,
+ sd: pg.Grid,
+ num_sol: np.ndarray,
+ ana_sol: Callable[[np.ndarray], np.ndarray],
+ relative: Optional[bool] = True,
+ etype: Optional[str] = "specific",
+ ) -> float:
+ """
+ Returns the l2 error computed against an analytical solution given as a function.
+
+ Args:
+ sd (pg.Grid): Grid, or a subclass.
+ num_sol (np.ndarray): Vector of the numerical solution.
+ ana_sol (Callable[[np.ndarray], np.ndarray]): Function that represents the
+ analytical solution.
+ relative (Optional[bool], optional): Compute the relative error or not.
+ Defaults to True.
+ etype (Optional[str], optional): Type of error computed. Defaults to "specific".
+
+ Returns:
+ float: The computed error.
+ """
+
+ err2 = 0
+ num_sol = num_sol.reshape((sd.dim, -1))
+ for d in np.arange(sd.dim):
+ ana_sol_dim = lambda x: ana_sol(x)[d]
+ num_sol_dim = num_sol[d]
+
+ err2_dim = self.scalar_discr.error_l2(
+ sd, num_sol_dim, ana_sol_dim, relative, etype
+ )
+ err2 += err2_dim**2
+ return np.sqrt(err2)
+
+
+class VecPwLinears(pg.VecDiscretization):
+ """
+ A class representing the discretization using vector piecewise linear functions.
+
+ Attributes:
+ keyword (str): The keyword for the vector discretization class.
+
+ Methods:
+ get_range_discr_class(self, dim: int) -> pg.Discretization:
+ Returns the discretization class for the range of the differential.
+
+ """
+
+ def __init__(self, keyword: str) -> None:
+ """
+ Initialize the vector discretization class.
+ The scalar discretization class is pg.PwLinears.
+
+ Args:
+ keyword (str): The keyword for the vector discretization class.
+
+ Returns:
+ None
+ """
+ super().__init__(keyword, pg.PwLinears)
+
+ def get_range_discr_class(self, dim: int) -> pg.Discretization:
+ """
+ Returns the discretization class for the range of the differential.
+
+ Args:
+ dim (int): The dimension of the range space.
+
+ Returns:
+ pg.Discretization: The discretization class for the range of the differential.
+ """
+ return self.scalar_discr.get_range_discr_class(dim)
+
+ def assemble_nat_bc(
+ self, sd: pg.Grid, func: Callable[[np.ndarray], np.ndarray], b_faces: np.ndarray
+ ) -> np.ndarray:
+ """
+ Assembles the natural boundary condition vector, equal to zero.
+
+ Args:
+ sd (pg.Grid): The grid object.
+ func (Callable[[np.ndarray], np.ndarray]): The function defining the
+ natural boundary condition.
+ b_faces (np.ndarray): The array of boundary faces.
+
+ Returns:
+ np.ndarray: The assembled natural boundary condition vector.
+ """
+ return np.zeros(self.ndof(sd))
diff --git a/src/pygeon/discretizations/vem/h1.py b/src/pygeon/discretizations/vem/h1.py
index a5561b9a..955f8679 100644
--- a/src/pygeon/discretizations/vem/h1.py
+++ b/src/pygeon/discretizations/vem/h1.py
@@ -150,7 +150,7 @@ def assemble_loc_proj_to_mon(
) -> np.ndarray:
"""
Computes the local projection onto the monomials
- Returns the coefficients {a_i} in a_0 + [a_1, a_2] \dot (x - c) / d
+ Returns the coefficients {a_i} in a_0 + [a_1, a_2] dot (x - c) / d
for each VL1 basis function.
Args:
diff --git a/src/pygeon/discretizations/vem/hdiv.py b/src/pygeon/discretizations/vem/hdiv.py
index 645e3dc6..846ebb9c 100644
--- a/src/pygeon/discretizations/vem/hdiv.py
+++ b/src/pygeon/discretizations/vem/hdiv.py
@@ -167,7 +167,7 @@ def assemble_nat_bc(
) -> np.ndarray:
"""
Assembles the natural boundary condition term
- (n dot q, func)_\Gamma.
+ (n dot q, func)_Gamma.
Args:
sd (pg.Grid): The grid object representing the computational domain.
@@ -408,7 +408,7 @@ def assemble_nat_bc(
) -> np.ndarray:
"""
Assembles the natural boundary condition term
- (n dot q, func)_\Gamma
+ (n dot q, func)_Gamma
Args:
sd (pg.Grid): The grid object representing the computational domain.
diff --git a/src/pygeon/grids/create_grid.py b/src/pygeon/grids/create_grid.py
index 5d9e66b2..f75a5254 100644
--- a/src/pygeon/grids/create_grid.py
+++ b/src/pygeon/grids/create_grid.py
@@ -25,12 +25,14 @@ def grid_from_domain(
Returns:
Either a pg.MixedDimensionalGrid or a pg.Grid, depending on the value of as_mdg.
"""
+ as_mdg = kwargs.get("as_mdg", True)
mesh_size_min = kwargs.get("mesh_size_min", mesh_size / 10)
+
mesh_kwargs = {"mesh_size_frac": mesh_size, "mesh_size_min": mesh_size_min}
- mdg = pp.create_fracture_network(domain=domain).mesh(mesh_kwargs)
+ mdg = pp.create_fracture_network(domain=domain).mesh(mesh_kwargs, **kwargs)
pg.convert_from_pp(mdg)
- if kwargs.get("as_mdg", True):
+ if as_mdg:
return mdg
else:
return mdg.subdomains(dim=mdg.dim_max())[0]
diff --git a/src/pygeon/grids/graph.py b/src/pygeon/grids/graph.py
index e61e6d70..b9c0c122 100644
--- a/src/pygeon/grids/graph.py
+++ b/src/pygeon/grids/graph.py
@@ -12,7 +12,7 @@ class Graph(pp.Grid):
def __init__(self, graph, dim=2):
"""Construct a pp.Grid like representation of a graph.
- The following informations are stored
+ The following information are stored
dim: the dimension of the graph, always
nodes: the coordinate of the nodes
cell_faces: the map from the cells (graph nodes) to the faces (graph edges)
@@ -89,7 +89,7 @@ def compute_ridges(self):
This function create mapping between the faces and the ridges in the graph,
the following data are created
- num_ridges: the number of rideges (graph cicles)
+ num_ridges: the number of ridges (graph cycles)
num_peaks: the number of peaks (always 0)
face_ridges: the map from the faces to the ridges
ridge_peaks: an empty matrix that maps the ridges to peaks
@@ -182,7 +182,7 @@ def edges_of_nodes(self, nodes):
def collapse(self, dim):
"""Collapse the graph removing the nodes with a given dimension.
- From a bipartite graph it is possible to create a dim = 2 graph (only fracures)
+ From a bipartite graph it is possible to create a dim = 2 graph (only fractures)
or a dim = 1 (only fracture intersections) graph.
"""
@@ -224,7 +224,7 @@ def shortest_paths(self, start, end):
return np.array(list(sp), dtype=np.object)
def not_shortest_paths(self, start, end, sp=None, cutoff=None):
- """Compyte all the not shortest paths between two graph nodes.
+ """Compute all the not shortest paths between two graph nodes.
It is possible to set a cutoff to make the computation faster leaving out long paths.
"""
@@ -291,8 +291,8 @@ def to_file(self, file_name):
"""
# make sure that an edge is sorted by dimension
- sort = (
- lambda e: e
+ sort = lambda e: (
+ e
if self.graph.nodes[e[0]]["dim"] > self.graph.nodes[e[1]]["dim"]
else np.flip(e)
)
diff --git a/src/pygeon/grids/md_grid.py b/src/pygeon/grids/md_grid.py
index 565af939..91bba4e4 100644
--- a/src/pygeon/grids/md_grid.py
+++ b/src/pygeon/grids/md_grid.py
@@ -89,26 +89,27 @@ def initialize_data(self) -> None:
Returns:
None
"""
+ unit_keyword = "unit"
for sd, data in self.subdomains(return_data=True):
perm = pp.SecondOrderTensor(np.ones(sd.num_cells))
data[pp.PARAMETERS] = pp.Parameters(
sd,
- ["unit"],
+ [unit_keyword],
[
{"second_order_tensor": perm},
],
)
- data[pp.DISCRETIZATION_MATRICES] = {}
+ data[pp.DISCRETIZATION_MATRICES] = {unit_keyword: {}}
for _, data in self.interfaces(return_data=True):
data[pp.PARAMETERS] = pp.Parameters(
sd,
- ["unit"],
+ [unit_keyword],
[
{"normal_diffusivity": 1.0},
],
)
- data[pp.DISCRETIZATION_MATRICES] = {}
+ data[pp.DISCRETIZATION_MATRICES] = {unit_keyword: {}}
def num_subdomain_faces(
self, cond: Optional[Callable[[pg.Grid], bool]] = None
diff --git a/src/pygeon/numerics/projections.py b/src/pygeon/numerics/projections.py
index e739725f..fd4e4247 100644
--- a/src/pygeon/numerics/projections.py
+++ b/src/pygeon/numerics/projections.py
@@ -72,7 +72,7 @@ def proj_faces_to_cells(
)
if discr is None:
- discr = pg.RT0("flow")
+ discr = pg.RT0("unit")
# Local mass matrices
for nn_sd, (sd, d_sd) in enumerate(mdg.subdomains(return_data=True)):
diff --git a/src/pygeon/numerics/spanningtree.py b/src/pygeon/numerics/spanningtree.py
index 7899e4b5..99dcb364 100644
--- a/src/pygeon/numerics/spanningtree.py
+++ b/src/pygeon/numerics/spanningtree.py
@@ -497,6 +497,70 @@ def compute_system(self, sd: pg.Grid) -> sps.csc_matrix:
return sps.vstack((-div, -asym)) @ self.expand
+class SpanningTreeCosserat(SpanningTreeElasticity):
+ """
+ Represents a class for computing the spanning tree for the Cosserat problem.
+
+ Attributes:
+ expand (sps.csc_matrix): The expanded matrix for spanning tree computation.
+ system (sps.csc_matrix): The computed system matrix.
+
+ Methods:
+ setup_system(self, mdg: pg.MixedDimensionalGrid, flagged_faces: np.ndarray) -> None:
+ Set up the system for the spanning tree algorithm.
+
+ compute_expand(self, sd: pg.Grid, flagged_faces: np.ndarray) -> sps.csc_matrix:
+ Compute the expanded matrix for spanning tree computation.
+
+ compute_system(self, sd: pg.Grid) -> sps.csc_matrix:
+ Computes the system matrix for the given grid.
+ """
+
+ def compute_expand(self, sd: pg.Grid, flagged_faces: np.ndarray) -> sps.csc_matrix:
+ """
+ Compute the expanded matrix for spanning tree computation.
+
+ Args:
+ sd (pg.Grid): The grid object.
+ flagged_faces (np.ndarray): Array of flagged faces.
+
+ Returns:
+ sps.csc_matrix: The expanded matrix for spanning tree computation.
+ """
+ dim_sig_omega = sd.dim * (sd.dim + 1) // 2
+
+ expand = pg.numerics.linear_system.create_restriction(flagged_faces).T.tocsc()
+
+ return sps.block_diag([expand] * dim_sig_omega, format="csc")
+
+ def compute_system(self, sd: pg.Grid) -> sps.csc_matrix:
+ """
+ Computes the system matrix for the given grid.
+
+ Args:
+ sd (pg.Grid): The grid object representing the domain.
+
+ Returns:
+ sps.csc_matrix: The computed system matrix.
+ """
+ assert sd.dim == 3, "Only implemented the 3D Raviart-Thomas version."
+
+ # first we assemble the B matrix
+ key = "tree"
+ vec_rt0 = pg.VecRT0(key)
+ vec_p0 = pg.VecPwConstants(key)
+
+ M = vec_p0.assemble_mass_matrix(sd)
+
+ div = M @ vec_rt0.assemble_diff_matrix(sd)
+ asym = M @ vec_rt0.assemble_asym_matrix(sd)
+
+ B = sps.bmat([[-div, None], [-asym, -div]], format="csc")
+
+ # create the solution operator
+ return B @ self.expand
+
+
class SpanningWeightedTrees:
"""
Class that can perform a spanning weighted trees solve, based on one of
diff --git a/tests/integration/test_elasticity.py b/tests/integration/test_elasticity.py
index 8c8dfc26..bf8d135d 100644
--- a/tests/integration/test_elasticity.py
+++ b/tests/integration/test_elasticity.py
@@ -180,9 +180,9 @@ def run_elasticity_2d(self, u_boundary, N):
asym = Mr @ vec_bdm1.assemble_asym_matrix(sd)
# fmt: off
- spp = sps.bmat([[ Ms, div.T, asym.T],
- [ -div, None, None],
- [-asym, None, None]], format = "csc")
+ spp = sps.bmat([[ Ms, div.T, -asym.T],
+ [-div, None, None],
+ [asym, None, None]], format = "csc")
# fmt: on
b_faces = sd.tags["domain_boundary_faces"]
diff --git a/tests/unit/test_projections.py b/tests/unit/test_projections.py
index 5763f5bd..83cce460 100644
--- a/tests/unit/test_projections.py
+++ b/tests/unit/test_projections.py
@@ -294,5 +294,4 @@ def test3(self):
if __name__ == "__main__":
- # ProjectionsUnitTest().test3()
unittest.main()
diff --git a/tests/unit/test_pwlinear.py b/tests/unit/test_pwlinear.py
index db457936..f63056b5 100644
--- a/tests/unit/test_pwlinear.py
+++ b/tests/unit/test_pwlinear.py
@@ -137,7 +137,11 @@ def test_eval_at_cell_centers(self):
discr = pg.PwLinears("P1b")
- self.assertRaises(NotImplementedError, discr.eval_at_cell_centers, sd)
+ known_func = np.ones(discr.ndof(sd))
+
+ P = discr.eval_at_cell_centers(sd)
+
+ self.assertTrue(np.allclose(P @ known_func, np.ones(sd.num_cells)))
def test_assemble_nat_bc(self):
dim = 2
diff --git a/tests/unit/test_spanningtree.py b/tests/unit/test_spanningtree.py
index 4208c8bd..a52fe34f 100644
--- a/tests/unit/test_spanningtree.py
+++ b/tests/unit/test_spanningtree.py
@@ -230,5 +230,57 @@ def test_for_errors(self):
self.assertRaises(NotImplementedError, pg.SpanningTreeElasticity, mdg)
+class SpanningTreeCosseratTest(unittest.TestCase):
+ def check(self, sd):
+ mdg = pg.as_mdg(sd)
+ pg.convert_from_pp(mdg)
+ mdg.compute_geometry()
+
+ B = self.assemble_B(mdg)
+ f = np.random.rand(B.shape[0])
+
+ sptr = pg.SpanningTreeCosserat(mdg)
+
+ s_f = sptr.solve(f)
+ self.assertTrue(np.allclose(B @ s_f, f))
+
+ def assemble_B(self, mdg):
+ sd = mdg.subdomains(dim=mdg.dim_max())[0]
+
+ key = "tree"
+ vec_rt0 = pg.VecRT0(key)
+ vec_p0 = pg.VecPwConstants(key)
+
+ M = vec_p0.assemble_mass_matrix(sd)
+
+ div = M @ vec_rt0.assemble_diff_matrix(sd)
+ asym = M @ vec_rt0.assemble_asym_matrix(sd)
+
+ return sps.bmat([[-div, None], [-asym, -div]])
+
+ def test_elasticity_struct_tet_grid(self):
+ sd = pp.StructuredTetrahedralGrid([1] * 3)
+ self.check(sd)
+
+ def test_elasticity_unstruct_tet_grid(self):
+ sd = pg.unit_grid(3, 1.0)
+ self.check(sd)
+
+ def test_assemble_SI(self):
+ N, dim = 3, 3
+ sd = pp.StructuredTetrahedralGrid([N] * dim, [1] * dim)
+ mdg = pg.as_mdg(sd)
+ pg.convert_from_pp(mdg)
+ mdg.compute_geometry()
+
+ sptr = pg.SpanningTreeCosserat(mdg)
+
+ SI = sptr.assemble_SI()
+ B = self.assemble_B(mdg)
+ check = sps.eye_array(B.shape[0]) - B @ SI
+
+ self.assertTrue(np.allclose(check.data, 0))
+
+
if __name__ == "__main__":
unittest.main()
diff --git a/tests/unit/test_vecbdm1.py b/tests/unit/test_vecbdm1.py
index 07e36705..c20b8802 100644
--- a/tests/unit/test_vecbdm1.py
+++ b/tests/unit/test_vecbdm1.py
@@ -6,6 +6,15 @@
class VecBDM1Test(unittest.TestCase):
+
+ def test_1d(self):
+ sd = pp.CartGrid([1])
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ vec_bdm1 = pg.VecBDM1("vecbdm1")
+ self.assertRaises(ValueError, vec_bdm1.assemble_asym_matrix, sd)
+
def test_trace_2d(self):
sd = pp.StructuredTriangleGrid([1] * 2, [1] * 2)
pg.convert_from_pp(sd)
@@ -52,6 +61,61 @@ def test_assemble_mass_matrix_2d(self):
M = vec_bdm1.assemble_mass_matrix(sd, data)
self.assertAlmostEqual(u.T @ M @ u, 30)
+ def test_assemble_mass_matrix_cosserat_2d(self):
+ N = 10
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5, "mu_c": 0.25}}}
+ M = vec_bdm1.assemble_mass_matrix_cosserat(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0]])
+ u = vec_bdm1.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 28)
+
+ def test_assemble_lumped_matrix_2d(self):
+ N = 10
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5}}}
+ M = vec_bdm1.assemble_lumped_matrix(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0]])
+ u = vec_bdm1.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 26)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0}}}
+ M = vec_bdm1.assemble_lumped_matrix(sd, data)
+ self.assertAlmostEqual(u.T @ M @ u, 30)
+
+ def test_assemble_lumped_matrix_cosserat_2d(self):
+ N = 10
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5, "mu_c": 0.25}}}
+ M = vec_bdm1.assemble_lumped_matrix_cosserat(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0]])
+ u = vec_bdm1.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 28)
+
def test_eval_at_cell_centers_2d(self):
N = 1
sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
@@ -72,6 +136,24 @@ def linear(x):
self.assertAlmostEqual(np.linalg.norm(eval - known), 0)
+ def test_proj_to_and_from_rt0_2d(self):
+ N = 1
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ def linear(x):
+ return np.array([x, 2 * x])
+
+ interp = vec_bdm1.interpolate(sd, linear)
+ interp_to_rt0 = vec_bdm1.proj_to_RT0(sd) @ interp
+ interp_from_rt0 = vec_bdm1.proj_from_RT0(sd) @ interp_to_rt0
+
+ self.assertAlmostEqual(np.linalg.norm(interp - interp_from_rt0), 0)
+
def test_range(self):
key = "vecbdm1"
vec_bdm1 = pg.VecBDM1(key)
@@ -88,19 +170,19 @@ def test_assemble_asym_matrix_2d(self):
fun = lambda _: np.array([[1, 2, 0], [4, 3, 0]])
u = vec_bdm1.interpolate(sd, fun)
- asym = vec_bdm1.assemble_asym_matrix(sd)
+ asym = vec_bdm1.assemble_asym_matrix(sd, False)
- p0 = pg.PwConstants("p0")
- cell_asym_u = p0.eval_at_cell_centers(sd) @ (asym @ u)
+ p1 = pg.PwLinears("p1")
+ cell_asym_u = p1.eval_at_cell_centers(sd) @ (asym @ u)
- self.assertTrue(np.allclose(cell_asym_u, -2))
+ self.assertTrue(np.allclose(cell_asym_u, 2))
def test_trace_3d(self):
sd = pp.StructuredTetrahedralGrid([1] * 3, [1] * 3)
pg.convert_from_pp(sd)
sd.compute_geometry()
- vec_bdm1 = pg.VecBDM1("vec_bdm1")
+ vec_bdm1 = pg.VecBDM1("vecbdm1")
B = vec_bdm1.assemble_trace_matrix(sd)
@@ -118,7 +200,7 @@ def test_ndof_3d(self):
pg.convert_from_pp(sd)
sd.compute_geometry()
- vec_bdm1 = pg.VecBDM1("vec_bdm1")
+ vec_bdm1 = pg.VecBDM1("vecbdm1")
self.assertEqual(vec_bdm1.ndof(sd), 162)
@@ -128,7 +210,7 @@ def test_assemble_mass_matrix_3d(self):
pg.convert_from_pp(sd)
sd.compute_geometry()
- key = "vec_bdm1"
+ key = "vecbdm1"
vec_bdm1 = pg.VecBDM1(key)
data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5}}}
@@ -143,13 +225,72 @@ def test_assemble_mass_matrix_3d(self):
M = vec_bdm1.assemble_mass_matrix(sd, data)
self.assertAlmostEqual(u.T @ M @ u, 32)
+ def test_assemble_mass_matrix_cosserat_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5, "mu_c": 0.25}}}
+ M = vec_bdm1.assemble_mass_matrix_cosserat(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0], [0, 1, 1]])
+ u = vec_bdm1.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 29.5)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0, "mu_c": 0.25}}}
+ M = vec_bdm1.assemble_mass_matrix_cosserat(sd, data)
+ self.assertAlmostEqual(u.T @ M @ u, 34.5)
+
+ def test_assemble_lumped_matrix_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5}}}
+ M = vec_bdm1.assemble_lumped_matrix(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0], [0, 1, 1]])
+ u = vec_bdm1.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 27)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0}}}
+ M = vec_bdm1.assemble_lumped_matrix(sd, data)
+ self.assertAlmostEqual(u.T @ M @ u, 32)
+
+ def test_assemble_lumped_matrix_cosserat_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5, "mu_c": 0.25}}}
+ M = vec_bdm1.assemble_lumped_matrix_cosserat(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0], [0, 1, 1]])
+ u = vec_bdm1.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 29.5)
+
def test_eval_at_cell_centers_3d(self):
N = 1
sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
pg.convert_from_pp(sd)
sd.compute_geometry()
- key = "vec_bdm1"
+ key = "vecbdm1"
vec_bdm1 = pg.VecBDM1(key)
def linear(x):
@@ -168,21 +309,39 @@ def test_assemble_asym_matrix_3d(self):
pg.convert_from_pp(sd)
sd.compute_geometry()
- key = "vec_bdm1"
+ key = "vecbdm1"
vec_bdm1 = pg.VecBDM1(key)
fun = lambda _: np.array([[1, 2, -1], [4, 3, 2], [1, 1, 1]])
u = vec_bdm1.interpolate(sd, fun)
- asym = vec_bdm1.assemble_asym_matrix(sd)
+ asym = vec_bdm1.assemble_asym_matrix(sd, False)
- p0 = pg.VecPwConstants("p0")
- cell_asym_u = p0.eval_at_cell_centers(sd) @ (asym @ u)
+ p1 = pg.VecPwLinears("p1")
+ cell_asym_u = p1.eval_at_cell_centers(sd) @ (asym @ u)
cell_asym_u = cell_asym_u.reshape((3, -1))
self.assertTrue(np.allclose(cell_asym_u[0], -1))
self.assertTrue(np.allclose(cell_asym_u[1], -2))
self.assertTrue(np.allclose(cell_asym_u[2], 2))
+ def test_proj_to_and_from_rt0_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecbdm1"
+ vec_bdm1 = pg.VecBDM1(key)
+
+ def linear(x):
+ return np.array([x, 2 * x, 3 * x])
+
+ interp = vec_bdm1.interpolate(sd, linear)
+ interp_to_rt0 = vec_bdm1.proj_to_RT0(sd) @ interp
+ interp_from_rt0 = vec_bdm1.proj_from_RT0(sd) @ interp_to_rt0
+
+ self.assertAlmostEqual(np.linalg.norm(interp - interp_from_rt0), 0)
+
if __name__ == "__main__":
unittest.main()
diff --git a/tests/unit/test_vecpwconstant.py b/tests/unit/test_vecpwconstant.py
index ee4b0f73..58a74714 100644
--- a/tests/unit/test_vecpwconstant.py
+++ b/tests/unit/test_vecpwconstant.py
@@ -13,7 +13,7 @@ def test_ndof(self):
sd.compute_geometry()
discr = pg.VecPwConstants("P0")
- assert discr.ndof(sd) == sd.num_cells * dim
+ self.assertTrue(discr.ndof(sd) == sd.num_cells * dim)
def test_assemble_mass_matrix(self):
dim = 2
@@ -179,8 +179,13 @@ def test_error_l2(self):
# fmt: on
ana_sol = lambda x: np.sin(x)
+ err = discr.error_l2(sd, num_sol, ana_sol, etype="standard")
+ err_known = 7.951651025069553e-08
+
+ self.assertTrue(np.allclose(err, err_known))
+
err = discr.error_l2(sd, num_sol, ana_sol)
- err_known = 5.622666361455642e-08
+ err_known = 0.5254595621486393
self.assertTrue(np.allclose(err, err_known))
diff --git a/tests/unit/test_vecpwlinear.py b/tests/unit/test_vecpwlinear.py
new file mode 100644
index 00000000..90cdbc1d
--- /dev/null
+++ b/tests/unit/test_vecpwlinear.py
@@ -0,0 +1,167 @@
+import unittest
+import numpy as np
+import scipy.sparse as sps
+
+import pygeon as pg
+import porepy as pp
+
+
+class VecPwLinearsTest(unittest.TestCase):
+ def test_ndof(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([2] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+ self.assertTrue(discr.ndof(sd) == sd.num_cells * dim * (dim + 1))
+
+ def test_assemble_mass_matrix(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+ M = discr.assemble_mass_matrix(sd)
+
+ # fmt: off
+ M_known_data = np.array(
+ [0.08333333, 0.04166667, 0.04166667, 0.04166667, 0.08333333,
+ 0.04166667, 0.04166667, 0.04166667, 0.08333333, 0.08333333,
+ 0.04166667, 0.04166667, 0.04166667, 0.08333333, 0.04166667,
+ 0.04166667, 0.04166667, 0.08333333, 0.08333333, 0.04166667,
+ 0.04166667, 0.04166667, 0.08333333, 0.04166667, 0.04166667,
+ 0.04166667, 0.08333333, 0.08333333, 0.04166667, 0.04166667,
+ 0.04166667, 0.08333333, 0.04166667, 0.04166667, 0.04166667,
+ 0.08333333]
+ )
+
+ M_known_indices = np.array(
+ [ 0, 1, 2, 0, 1, 2, 0, 1, 2, 3, 4, 5, 3, 4, 5, 3, 4,
+ 5, 6, 7, 8, 6, 7, 8, 6, 7, 8, 9, 10, 11, 9, 10, 11, 9,
+ 10, 11]
+ )
+
+ M_known_indptr = np.array(
+ [ 0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36]
+ )
+ # fmt: on
+
+ self.assertTrue(np.allclose(M.data, M_known_data))
+ self.assertTrue(np.allclose(M.indptr, M_known_indptr))
+ self.assertTrue(np.allclose(M.indices, M_known_indices))
+
+ def test_assemble_lumped_matrix(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+
+ M = discr.assemble_lumped_matrix(sd)
+
+ # fmt: off
+ M_known_data = np.array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]) / 6
+
+ M_known_indices = np.array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
+
+ M_known_indptr = np.array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])
+ # fmt: on
+
+ self.assertTrue(np.allclose(M.data, M_known_data))
+ self.assertTrue(np.allclose(M.indptr, M_known_indptr))
+ self.assertTrue(np.allclose(M.indices, M_known_indices))
+
+ def test_assemble_diff_matrix(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+
+ self.assertRaises(NotImplementedError, discr.assemble_diff_matrix, sd)
+
+ def test_assemble_stiff_matrix(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+
+ self.assertRaises(NotImplementedError, discr.assemble_stiff_matrix, sd)
+
+ def test_interpolate(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+
+ func = lambda x: np.sin(x) # Example function
+ self.assertRaises(NotImplementedError, discr.interpolate, sd, func)
+
+ def test_eval_at_cell_centers(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+
+ P = discr.eval_at_cell_centers(sd)
+
+ # fmt: off
+ P_known_data = np.ones(discr.ndof(sd)) / (sd.dim + 1)
+
+ P_known_indices = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3])
+
+ P_known_indptr = np.arange(discr.ndof(sd)+1)
+ # fmt: on
+
+ self.assertTrue(np.allclose(P.data, P_known_data))
+ self.assertTrue(np.allclose(P.indptr, P_known_indptr))
+ self.assertTrue(np.allclose(P.indices, P_known_indices))
+
+ def test_assemble_nat_bc(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+
+ func = lambda x: np.sin(x[0]) # Example function
+
+ b_faces = np.array([0, 1, 3]) # Example boundary faces
+ vals = discr.assemble_nat_bc(sd, func, b_faces)
+
+ vals_known = np.zeros(discr.ndof(sd))
+
+ self.assertTrue(np.allclose(vals, vals_known))
+
+ def test_get_range_discr_class(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+
+ self.assertRaises(
+ NotImplementedError,
+ discr.get_range_discr_class,
+ dim,
+ )
+
+ def test_error_l2(self):
+ dim = 2
+ sd = pp.StructuredTriangleGrid([1] * dim, [1] * dim)
+ sd.compute_geometry()
+
+ discr = pg.VecPwLinears("P1")
+ # fmt: off
+ num_sol = np.empty(0)
+ # fmt: on
+ ana_sol = lambda x: np.sin(x)
+
+ self.assertRaises(NotImplementedError, discr.error_l2, sd, num_sol, ana_sol)
+
+
+if __name__ == "__main__":
+ unittest.main()
diff --git a/tests/unit/test_vecrt0.py b/tests/unit/test_vecrt0.py
new file mode 100644
index 00000000..9b8bbabe
--- /dev/null
+++ b/tests/unit/test_vecrt0.py
@@ -0,0 +1,222 @@
+import unittest
+import numpy as np
+
+import porepy as pp
+import pygeon as pg
+
+
+class VecRT0Test(unittest.TestCase):
+ def test_trace_2d(self):
+ sd = pp.StructuredTriangleGrid([1] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ vec_rt0 = pg.VecRT0("vecrt0")
+
+ B = vec_rt0.assemble_trace_matrix(sd)
+
+ fun = lambda _: np.array([[1, 2, 0], [3, 4, 0]])
+ u = vec_rt0.interpolate(sd, fun)
+
+ trace = B @ u
+
+ self.assertTrue(np.allclose(trace, 5))
+
+ def test_ndof_2d(self):
+ sd = pp.StructuredTriangleGrid([1] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ vec_rt0 = pg.VecRT0("vecrt0")
+
+ self.assertEqual(vec_rt0.ndof(sd), 10)
+
+ def test_assemble_mass_matrix_2d(self):
+ N = 10
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5}}}
+ M = vec_rt0.assemble_mass_matrix(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0]])
+ u = vec_rt0.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 26)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0}}}
+ M = vec_rt0.assemble_mass_matrix(sd, data)
+ self.assertAlmostEqual(u.T @ M @ u, 30)
+
+ def test_assemble_mass_matrix_cosserat_2d(self):
+ N = 10
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5, "mu_c": 0.25}}}
+ M = vec_rt0.assemble_mass_matrix_cosserat(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0]])
+ u = vec_rt0.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 28)
+
+ def test_eval_at_cell_centers_2d(self):
+ N = 1
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ def linear(x):
+ return np.array([x, 2 * x])
+
+ interp = vec_rt0.interpolate(sd, linear)
+ eval = vec_rt0.eval_at_cell_centers(sd) @ interp
+ eval = np.reshape(eval, (3, -1), order="F")
+ eval = np.vstack(np.split(eval, sd.dim, axis=1))
+
+ known = np.array([linear(x).ravel() for x in sd.cell_centers.T]).T
+
+ self.assertAlmostEqual(np.linalg.norm(eval - known), 0)
+
+ def test_range(self):
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+ self.assertTrue(vec_rt0.get_range_discr_class(2) is pg.VecPwConstants)
+
+ def test_assemble_asym_matrix_2d(self):
+ N = 1
+ sd = pp.StructuredTriangleGrid([N] * 2, [1] * 2)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0]])
+ u = vec_rt0.interpolate(sd, fun)
+ asym = vec_rt0.assemble_asym_matrix(sd)
+
+ p0 = pg.PwConstants("p0")
+ cell_asym_u = p0.eval_at_cell_centers(sd) @ (asym @ u)
+
+ self.assertTrue(np.allclose(cell_asym_u, 2))
+
+ def test_trace_3d(self):
+ sd = pp.StructuredTetrahedralGrid([1] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ vec_rt0 = pg.VecRT0("vecrt0")
+
+ B = vec_rt0.assemble_trace_matrix(sd)
+
+ fun = lambda _: np.array([[1, 2, 3], [4, 5, 6], [0, 0, 7]])
+ u = vec_rt0.interpolate(sd, fun)
+
+ trace = B @ u
+
+ self.assertTrue(np.allclose(trace, 13))
+
+ def test_ndof_3d(self):
+ sd = pp.StructuredTetrahedralGrid([1] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ vec_rt0 = pg.VecRT0("vecrt0")
+
+ self.assertEqual(vec_rt0.ndof(sd), 54)
+
+ def test_assemble_mass_matrix_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5}}}
+ M = vec_rt0.assemble_mass_matrix(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0], [0, 1, 1]])
+ u = vec_rt0.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 27)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0}}}
+ M = vec_rt0.assemble_mass_matrix(sd, data)
+ self.assertAlmostEqual(u.T @ M @ u, 32)
+
+ def test_assemble_mass_matrix_cosserat_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ data = {pp.PARAMETERS: {key: {"mu": 0.5, "lambda": 0.5, "mu_c": 0.25}}}
+ M = vec_rt0.assemble_mass_matrix_cosserat(sd, data)
+
+ fun = lambda _: np.array([[1, 2, 0], [4, 3, 0], [0, 1, 1]])
+ u = vec_rt0.interpolate(sd, fun)
+
+ self.assertAlmostEqual(u.T @ M @ u, 29.5)
+
+ def test_eval_at_cell_centers_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ def linear(x):
+ return np.array([x, 2 * x, -x])
+
+ interp = vec_rt0.interpolate(sd, linear)
+ eval = vec_rt0.eval_at_cell_centers(sd) @ interp
+ eval = np.reshape(eval, (3, -1), order="F")
+ eval = np.vstack(np.split(eval, sd.dim, axis=1))
+
+ known = np.array([linear(x).ravel() for x in sd.cell_centers.T]).T
+ self.assertAlmostEqual(np.linalg.norm(eval - known), 0)
+
+ def test_assemble_asym_matrix_3d(self):
+ N = 1
+ sd = pp.StructuredTetrahedralGrid([N] * 3, [1] * 3)
+ pg.convert_from_pp(sd)
+ sd.compute_geometry()
+
+ key = "vecrt0"
+ vec_rt0 = pg.VecRT0(key)
+
+ fun = lambda _: np.array([[1, 2, -1], [4, 3, 2], [1, 1, 1]])
+ u = vec_rt0.interpolate(sd, fun)
+ asym = vec_rt0.assemble_asym_matrix(sd)
+
+ p0 = pg.VecPwConstants("p0")
+ cell_asym_u = p0.eval_at_cell_centers(sd) @ (asym @ u)
+ cell_asym_u = cell_asym_u.reshape((3, -1))
+
+ self.assertTrue(np.allclose(cell_asym_u[0], -1))
+ self.assertTrue(np.allclose(cell_asym_u[1], -2))
+ self.assertTrue(np.allclose(cell_asym_u[2], 2))
+
+
+if __name__ == "__main__":
+ unittest.main()
diff --git a/tutorials/cosserat.ipynb b/tutorials/cosserat.ipynb
new file mode 100644
index 00000000..dff6017f
--- /dev/null
+++ b/tutorials/cosserat.ipynb
@@ -0,0 +1,310 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "essential-american",
+ "metadata": {},
+ "source": [
+ "# Cosserat equation in mixed form\n",
+ "\n",
+ "In this tutorial we present how to solve the Cosserat equation with [PyGeoN](https://github.com/compgeo-mox/pygeon) in mixed form. The unknowns are the stress $\\sigma$, the micro-stress $\\omega$, the displacement $u$, and the rotation $r$.\n",
+ "\n",
+ "Let $\\Omega=(0,1)^2$ with boundary $\\partial \\Omega$ and outward unit normal ${\\nu}$. Given the macroscopic \n",
+ "$\\lambda_\\sigma$ Lamé constant and $\\mu_\\sigma$ the Kirchhoff modulus and their microscopic version $\\lambda_\\omega$ and $\\mu_\\omega$, \n",
+ "we want to solve the following problem: find $(\\sigma, \\omega, u, r)$ such that\n",
+ "\\begin{align*}\n",
+ "&A_\\sigma \\sigma - \\epsilon(u) + {\\rm asym}^* r = 0\\\\\n",
+ "&A_\\omega \\omega - \\nabla r = 0\\\\\n",
+ "&\\nabla \\cdot \\sigma = - b\\\\\n",
+ "&{\\rm asym}\\, \\sigma - \\nabla \\cdot \\omega = 0\n",
+ "\\end{align*}\n",
+ "with $\\epsilon$ the symmetric gradient and $b$ a body force, which is set to $0$ in this example.\n",
+ "The operator $A_\\sigma$ is defined as\n",
+ "$$\n",
+ "A_\\sigma \\sigma = \\frac1{2\\mu_\\sigma} \\left( \\mathrm{sym} \\sigma - \\frac{\\lambda_\\sigma}{2\\mu_\\sigma + d \\lambda_\\sigma} \\mathrm{Tr} (\\sigma) I\\right) + \\frac1{2\\mu_\\sigma^c} \\mathrm{skw} \\sigma\n",
+ "$$\n",
+ "with $d$ the dimension of the domain, here $d=2$. While the operator $A_\\omega$ depends on the dimension \n",
+ "and is defined as \n",
+ "$$\n",
+ "A_\\omega \\omega = \\frac1{2\\mu_\\omega} \\left( \\mathrm{sym} \\omega - \\frac{\\lambda_\\omega}{2\\mu_\\omega + d \\lambda_\\omega} \\mathrm{Tr} (\\omega) I\\right) + \\frac1{2\\mu_\\omega^c} \\mathrm{skw} \\omega\n",
+ " \\qquad d = 3\\\\\n",
+ "A_\\omega \\omega = \\frac{1}{2\\mu_\\omega}\\omega \\qquad d = 2\n",
+ "$$\n",
+ "Finally, ${\\rm asym}$ is the asymmetric operator, which is defined as followed\n",
+ "$$\n",
+ " {\\rm asym }\\, \\sigma = \\begin{bmatrix}\n",
+ " \\sigma_{32} - \\sigma_{23} \\\\\n",
+ " \\sigma_{13} - \\sigma_{31} \\\\\n",
+ " \\sigma_{21} - \\sigma_{12}\n",
+ " \\end{bmatrix} \n",
+ " \\qquad\n",
+ " {\\rm asym}^* r =\n",
+ " \\begin{bmatrix}\n",
+ " 0 & -r_3 & r_2 \\\\\n",
+ " r_3 & 0 & -r_1 \\\\\n",
+ " -r_2 & r_1 & 0\n",
+ " \\end{bmatrix}\n",
+ " \\qquad\n",
+ " d = 3\\\\\n",
+ " {\\rm asym} \\, \\sigma = \\sigma_{21} - \\sigma_{12} \\qquad\n",
+ " {\\rm asym}^* r =\n",
+ " \\begin{bmatrix}\n",
+ " 0 & -r \\\\\n",
+ " r & 0\n",
+ " \\end{bmatrix} \\qquad\n",
+ " d = 2.\n",
+ "$$\n",
+ "For this test case we set the following boundary conditions related to the so-called footstep problem:\n",
+ "$$ u = 0 \\text{ on } \\partial_{bottom} \\Omega \\qquad \\nu \\cdot \\sigma = [0, 0]^\\top \\text{ on } \\partial_{left} \\Omega \\cup \\partial_{right} \\Omega \\qquad \\nu \\cdot \\sigma = [0, -1]^\\top \\text{ on } \\partial_{top} \\Omega$$\n",
+ "\n",
+ "We present *step-by-step* how to create the grid, declare the problem data, and finally solve the problem."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "planned-danger",
+ "metadata": {},
+ "source": [
+ "First we import some of the standard modules, like `numpy` and `scipy.sparse`. Since PyGeoN is based on [PorePy](https://github.com/pmgbergen/porepy) we import both modules."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "dietary-perth",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import scipy.sparse as sps\n",
+ "\n",
+ "import porepy as pp\n",
+ "import pygeon as pg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "roman-glossary",
+ "metadata": {},
+ "source": [
+ "We create now the grid, since we use vector RT0 for ${\\sigma}$, RT0 for $\\omega$, vector piece-wise constant for $u$ and piece-wise constant for $r$, we are restricted to simplices. In this example we consider a 2-dimensional structured grid, but the presented code will work also in 3d with the only attention that in 3d the micro-stress is a vector RT0 and the rotations is a vector piece-wise constant."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "spectacular-saturn",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "mesh_size = 0.05\n",
+ "dim = 2\n",
+ "\n",
+ "sd = pg.unit_grid(dim, mesh_size, as_mdg=False)\n",
+ "sd.compute_geometry()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d625fca2",
+ "metadata": {},
+ "source": [
+ "We define now the finite element spaces for all the variables."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "52400a88",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "key = \"cosserat\"\n",
+ "\n",
+ "# finite element spaces\n",
+ "vec_rt0 = pg.VecRT0(key)\n",
+ "rt0 = pg.RT0(key)\n",
+ "vec_p0 = pg.VecPwConstants(key)\n",
+ "p0 = pg.PwConstants(key)\n",
+ "\n",
+ "# build the degrees of freedom\n",
+ "dofs = np.array([vec_rt0.ndof(sd), rt0.ndof(sd), vec_p0.ndof(sd), p0.ndof(sd)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "precious-belle",
+ "metadata": {},
+ "source": [
+ "With the following code we set the data, in particular the Lamé and the Kirchhoff modulus, and the boundary conditions. Since we need to identify each side of $\\partial \\Omega$ we need few steps."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "spare-person",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# data for sigma\n",
+ "data_sigma = {pp.PARAMETERS: {key: {\"lambda\": 1, \"mu\": 0.5, \"mu_c\": 0.25}}}\n",
+ "\n",
+ "# data for omega\n",
+ "sot = pp.SecondOrderTensor(np.ones(sd.num_cells))\n",
+ "data_omega = {pp.PARAMETERS: {key: {\"second_order_tensor\": sot}}}\n",
+ "\n",
+ "# select the portions of the boundary where we want to apply the boundary conditions\n",
+ "bottom = np.isclose(sd.face_centers[1, :], 0)\n",
+ "top = np.isclose(sd.face_centers[1, :], 1)\n",
+ "left = np.isclose(sd.face_centers[0, :], 0)\n",
+ "right = np.isclose(sd.face_centers[0, :], 1)\n",
+ "\n",
+ "# select the part of the boundary where we apply the essential boundary conditions\n",
+ "b_faces = np.logical_or.reduce((top, left, right))\n",
+ "ess_dof = np.tile(b_faces, sd.dim)\n",
+ "\n",
+ "# function for the essential boundary conditions\n",
+ "val = np.array([[0, 0, 0], [0, -1, 0]])\n",
+ "fct = lambda pt: val if np.isclose(pt[1], 1) else 0 * val\n",
+ "\n",
+ "# interpolate the essential boundary conditions\n",
+ "ess_val = vec_rt0.interpolate(sd, fct)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "secure-flesh",
+ "metadata": {},
+ "source": [
+ "Once the data are assigned to the grid, we construct the matrices. Once the latter is created, we also construct the right-hand side containing the boundary conditions."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "romance-findings",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# assemble mass matrices\n",
+ "Ms = vec_rt0.assemble_mass_matrix_cosserat(sd, data_sigma)\n",
+ "Mw = rt0.assemble_mass_matrix(sd, data_omega)\n",
+ "Mu = vec_p0.assemble_mass_matrix(sd)\n",
+ "Mr = p0.assemble_mass_matrix(sd)\n",
+ "\n",
+ "# assemble matrices associated with the divergence and asymmetry operators\n",
+ "div_s = Mu @ vec_rt0.assemble_diff_matrix(sd)\n",
+ "asym = Mr @ vec_rt0.assemble_asym_matrix(sd)\n",
+ "div_w = Mr @ rt0.assemble_diff_matrix(sd)\n",
+ "\n",
+ "# assemble the saddle point problem matrix\n",
+ "# fmt: off\n",
+ "A = sps.block_diag([Ms, Mw])\n",
+ "B = sps.bmat([[-div_s, None],\n",
+ " [ asym, -div_w]])\n",
+ "spp = sps.bmat([[A, -B.T],\n",
+ " [B, None]], format=\"csc\")\n",
+ "# fmt: on"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "mobile-nirvana",
+ "metadata": {},
+ "source": [
+ "We need to solve the linear system, PyGeoN provides a framework for that. The actual imposition of essential boundary conditions (stress boundary conditions) might change the symmetry of the global system, the class `pg.LinearSystem` preserves this structure by internally eliminating these degrees of freedom."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "subtle-wonder",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# solve the linear system\n",
+ "ls = pg.LinearSystem(spp, np.zeros(spp.shape[0]))\n",
+ "ls.flag_ess_bc(np.where(ess_dof)[0], ess_val)\n",
+ "x = ls.solve()\n",
+ "\n",
+ "# split the solution into the components\n",
+ "idx = np.cumsum(dofs[:-1])\n",
+ "sigma, omega, u, r = np.split(x, idx)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "pacific-alpha",
+ "metadata": {},
+ "source": [
+ "We finally export the solution to be visualized by [ParaView](https://www.paraview.org/)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "satisfactory-jerusalem",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# post process of the displacement\n",
+ "proj_u = vec_p0.eval_at_cell_centers(sd)\n",
+ "cell_u = (proj_u @ u).reshape((sd.dim, -1), order=\"C\")\n",
+ "# since we are in 2d we need to add the z component for the exporting\n",
+ "cell_u = np.vstack((cell_u, np.zeros(cell_u.shape[1])))\n",
+ "\n",
+ "# post process of the rotation\n",
+ "proj_r = p0.eval_at_cell_centers(sd)\n",
+ "cell_r = proj_r @ r\n",
+ "\n",
+ "save = pp.Exporter(sd, \"sol\", folder_name=key)\n",
+ "save.write_vtu([(\"cell_u\", cell_u), (\"cell_r\", cell_r)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "developing-mobile",
+ "metadata": {},
+ "source": [
+ "A representation of the computed solution is given below.
\n",
+ ""
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.12.4"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "e4cc1db98167c7fd7d55a1da8057731abc6cd6fe154328a2ae319df8aab4e24d"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tutorials/elasticity_mixed.ipynb b/tutorials/elasticity_mixed.ipynb
new file mode 100644
index 00000000..40ec2f36
--- /dev/null
+++ b/tutorials/elasticity_mixed.ipynb
@@ -0,0 +1,294 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "essential-american",
+ "metadata": {},
+ "source": [
+ "# Elasticity equation in mixed form\n",
+ "\n",
+ "In this tutorial we present how to solve the elasticity equation with [PyGeoN](https://github.com/compgeo-mox/pygeon) in mixed form. The unknowns are the stress $\\sigma$, the displacement $u$, and the rotation $r$.\n",
+ "\n",
+ "Let $\\Omega=(0,1)^2$ with boundary $\\partial \\Omega$ and outward unit normal ${\\nu}$. Given \n",
+ "$\\lambda$ Lamé constant and $\\mu$ the Kirchhoff modulus, we want to solve the following problem: find $(\\sigma, u, r)$ such that\n",
+ "\\begin{align*}\n",
+ "&A \\sigma - \\epsilon(u) - {\\rm asym}^* r = 0\\\\\n",
+ "&\\nabla \\cdot \\sigma = - b\\\\\n",
+ "&{\\rm asym}\\, \\sigma = 0\n",
+ "\\end{align*}\n",
+ "with $\\epsilon$ the symmetric gradient and $b$ a body force, which is set to $0$ in this example.\n",
+ "The operator $A$ is defined as\n",
+ "$$\n",
+ "A \\sigma = \\frac{1}{2\\mu}\\left( \\sigma -\\frac{\\lambda}{2\\mu + d \\lambda} {\\rm Tr}(\\sigma) I\\right)\n",
+ "$$\n",
+ "with $d$ the dimension of the domain, here $d=2$.\n",
+ "Finally, ${\\rm asym}$ is the asymmetric operator, which is defined as followed\n",
+ "$$\n",
+ " {\\rm asym }\\, \\sigma = \\begin{bmatrix}\n",
+ " \\sigma_{32} - \\sigma_{23} \\\\\n",
+ " \\sigma_{13} - \\sigma_{31} \\\\\n",
+ " \\sigma_{21} - \\sigma_{12}\n",
+ " \\end{bmatrix} \n",
+ " \\qquad\n",
+ " {\\rm asym}^* r =\n",
+ " \\begin{bmatrix}\n",
+ " 0 & -r_3 & r_2 \\\\\n",
+ " r_3 & 0 & -r_1 \\\\\n",
+ " -r_2 & r_1 & 0\n",
+ " \\end{bmatrix}\n",
+ " \\qquad\n",
+ " d = 3\\\\\n",
+ " {\\rm asym} \\, \\sigma = \\sigma_{21} - \\sigma_{12} \\qquad\n",
+ " {\\rm asym}^* r =\n",
+ " \\begin{bmatrix}\n",
+ " 0 & -r \\\\\n",
+ " r & 0\n",
+ " \\end{bmatrix} \\qquad\n",
+ " d = 2.\n",
+ "$$\n",
+ "For this test case we set the following boundary conditions related to the so-called footstep problem:\n",
+ "$$ u = 0 \\text{ on } \\partial_{bottom} \\Omega \\qquad \\nu \\cdot \\sigma = [0, 0]^\\top \\text{ on } \\partial_{left} \\Omega \\cup \\partial_{right} \\Omega \\qquad \\nu \\cdot \\sigma = [0, -1]^\\top \\text{ on } \\partial_{top} \\Omega$$\n",
+ "\n",
+ "We present *step-by-step* how to create the grid, declare the problem data, and finally solve the problem."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "planned-danger",
+ "metadata": {},
+ "source": [
+ "First we import some of the standard modules, like `numpy` and `scipy.sparse`. Since PyGeoN is based on [PorePy](https://github.com/pmgbergen/porepy) we import both modules."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "dietary-perth",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import scipy.sparse as sps\n",
+ "\n",
+ "import porepy as pp\n",
+ "import pygeon as pg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "roman-glossary",
+ "metadata": {},
+ "source": [
+ "We create now the grid, since we use vector BDM1 for ${\\sigma}$, vector piece-wise constant for $u$ and piece-wise constant for $r$, we are restricted to simplices. In this example we consider a 2-dimensional structured grid, but the presented code will work also in 3d with the only attention that in 3d the rotations should be discretized with vector piece-wise constant finite elements."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "spectacular-saturn",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "mesh_size = 0.05\n",
+ "dim = 2\n",
+ "\n",
+ "sd = pg.unit_grid(dim, mesh_size, as_mdg=False)\n",
+ "sd.compute_geometry()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d625fca2",
+ "metadata": {},
+ "source": [
+ "We define now the finite element spaces for all the variables."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "52400a88",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "key = \"elasticity\"\n",
+ "\n",
+ "# finite element spaces\n",
+ "vec_bdm1 = pg.VecBDM1(key)\n",
+ "vec_p0 = pg.VecPwConstants(key)\n",
+ "p0 = pg.PwConstants(key)\n",
+ "\n",
+ "# build the degrees of freedom\n",
+ "dofs = np.array([vec_bdm1.ndof(sd), vec_p0.ndof(sd), p0.ndof(sd)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "precious-belle",
+ "metadata": {},
+ "source": [
+ "With the following code we set the data, in particular the Lamé and the Kirchhoff modulus, and the boundary conditions. Since we need to identify each side of $\\partial \\Omega$ we need few steps."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "spare-person",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "data = {pp.PARAMETERS: {key: {\"lambda\": 1, \"mu\": 0.5}}}\n",
+ "\n",
+ "# select the portions of the boundary where we want to apply the boundary conditions\n",
+ "bottom = np.isclose(sd.face_centers[1, :], 0)\n",
+ "top = np.isclose(sd.face_centers[1, :], 1)\n",
+ "left = np.isclose(sd.face_centers[0, :], 0)\n",
+ "right = np.isclose(sd.face_centers[0, :], 1)\n",
+ "\n",
+ "# select the part of the boundary where we apply the essential boundary conditions\n",
+ "b_faces = np.logical_or.reduce((top, left, right))\n",
+ "ess_dof = np.tile(b_faces, sd.dim**2)\n",
+ "\n",
+ "# function for the essential boundary conditions\n",
+ "val = np.array([[0, 0, 0], [0, -1, 0]])\n",
+ "fct = lambda pt: val if np.isclose(pt[1], 1) else 0 * val\n",
+ "\n",
+ "# interpolate the essential boundary conditions\n",
+ "ess_val = vec_bdm1.interpolate(sd, fct)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "secure-flesh",
+ "metadata": {},
+ "source": [
+ "Once the data are assigned to the grid, we construct the matrices. Once the latter is created, we also construct the right-hand side containing the boundary conditions."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "romance-findings",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# assemble mass matrices\n",
+ "Ms = vec_bdm1.assemble_mass_matrix(sd, data)\n",
+ "Mu = vec_p0.assemble_mass_matrix(sd)\n",
+ "Mr = p0.assemble_mass_matrix(sd)\n",
+ "\n",
+ "# assemble matrices associated with the divergence and asymmetry operators\n",
+ "div = Mu @ vec_bdm1.assemble_diff_matrix(sd)\n",
+ "asym = Mr @ vec_bdm1.assemble_asym_matrix(sd)\n",
+ "\n",
+ "# assemble the saddle point problem matrix\n",
+ "# fmt: off\n",
+ "spp = sps.bmat([[ Ms, div.T, -asym.T],\n",
+ " [-div, None, None],\n",
+ " [asym, None, None]], format = \"csc\")\n",
+ "# fmt: on"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "mobile-nirvana",
+ "metadata": {},
+ "source": [
+ "We need to solve the linear system, PyGeoN provides a framework for that. The actual imposition of essential boundary conditions (stress boundary conditions) might change the symmetry of the global system, the class `pg.LinearSystem` preserves this structure by internally eliminating these degrees of freedom."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "subtle-wonder",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# solve the linear system\n",
+ "ls = pg.LinearSystem(spp, np.zeros(spp.shape[0]))\n",
+ "ls.flag_ess_bc(np.where(ess_dof)[0], ess_val)\n",
+ "x = ls.solve()\n",
+ "\n",
+ "# split the solution into the components\n",
+ "idx = np.cumsum(dofs[:-1])\n",
+ "sigma, u, r = np.split(x, idx)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "pacific-alpha",
+ "metadata": {},
+ "source": [
+ "We finally export the solution to be visualized by [ParaView](https://www.paraview.org/)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "satisfactory-jerusalem",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# post process of the displacement\n",
+ "proj_u = vec_p0.eval_at_cell_centers(sd)\n",
+ "cell_u = (proj_u @ u).reshape((sd.dim, -1), order=\"C\")\n",
+ "# since we are in 2d we need to add the z component for the exporting\n",
+ "cell_u = np.vstack((cell_u, np.zeros(cell_u.shape[1])))\n",
+ "\n",
+ "# post process of the rotation\n",
+ "proj_r = p0.eval_at_cell_centers(sd)\n",
+ "cell_r = proj_r @ r\n",
+ "\n",
+ "save = pp.Exporter(sd, \"sol\")\n",
+ "save.write_vtu([(\"cell_u\", cell_u), (\"cell_r\", cell_r)])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "developing-mobile",
+ "metadata": {},
+ "source": [
+ "A representation of the computed solution is given below.
\n",
+ "\n",
+ "
\n",
+ ""
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.11.1"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "e4cc1db98167c7fd7d55a1da8057731abc6cd6fe154328a2ae319df8aab4e24d"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tutorials/fig/cosserat.png b/tutorials/fig/cosserat.png
new file mode 100644
index 00000000..1d764e3a
Binary files /dev/null and b/tutorials/fig/cosserat.png differ