diff --git a/.zenodo.json b/.zenodo.json index c8d766fe..51a9db6f 100755 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,8 +1,8 @@ { "description": "

cashocs is a computational, adjoint-based shape optimization and optimal control software.

", "license": "GPL-3.0+", - "title": "cashocs v1.8.9", - "version": "v1.8.9", + "title": "cashocs v1.8.10", + "version": "v1.8.10", "upload_type": "software", "creators": [ { diff --git a/cashocs/__init__.py b/cashocs/__init__.py index 65540fc0..879db1eb 100755 --- a/cashocs/__init__.py +++ b/cashocs/__init__.py @@ -59,7 +59,7 @@ from cashocs.nonlinear_solvers import newton_solve from cashocs.nonlinear_solvers import picard_iteration -__version__ = "1.8.9" +__version__ = "1.8.10" __citation__ = """ @Article{Blauth2021cashocs, diff --git a/cashocs/geometry/mesh_quality.py b/cashocs/geometry/mesh_quality.py index 1b39610f..ada7679f 100644 --- a/cashocs/geometry/mesh_quality.py +++ b/cashocs/geometry/mesh_quality.py @@ -274,6 +274,27 @@ def __init__(self) -> None: """Initializes self.""" pass + @classmethod + def _skewness(cls, mesh: fenics.Mesh) -> np.ndarray: + r"""Computes the skewness of the mesh. + + Args: + mesh: The mesh whose quality shall be computed. + + Returns: + The element wise skewness of the mesh on process 0. + + """ + comm = fenics.MPI.comm_world + skewness_array = cls._quality_object.skewness(mesh).array() + skewness_list: np.ndarray = comm.gather(skewness_array, root=0) + if comm.rank == 0: + skewness_list = np.concatenate(skewness_list, axis=None) + else: + skewness_list = np.zeros(1) + + return skewness_list + @classmethod def min_skewness(cls, mesh: fenics.Mesh) -> float: r"""Computes the minimal skewness of the mesh. @@ -303,7 +324,16 @@ def min_skewness(cls, mesh: fenics.Mesh) -> float: The minimum skewness of the mesh. """ - quality: float = np.min(cls._quality_object.skewness(mesh).array()) + skewness_list = cls._skewness(mesh) + comm = fenics.MPI.comm_world + + if comm.rank == 0: + qual = float(np.min(skewness_list)) + else: + qual = None + + quality: float = comm.bcast(qual, root=0) + return quality @classmethod @@ -335,9 +365,39 @@ def avg_skewness(cls, mesh: fenics.Mesh) -> float: The average skewness of the mesh. """ - quality: float = np.average(cls._quality_object.skewness(mesh).array()) + skewness_list = cls._skewness(mesh) + comm = fenics.MPI.comm_world + + if comm.rank == 0: + qual = float(np.average(skewness_list)) + else: + qual = None + + quality: float = comm.bcast(qual, root=0) + return quality + @classmethod + def _maximum_angle(cls, mesh: fenics.Mesh) -> np.ndarray: + r"""Computes the largest angle of each element. + + Args: + mesh: The mesh, whose quality shall be computed. + + Returns: + The maximum angle quality measure for each element on process 0. + + """ + comm = fenics.MPI.comm_world + maximum_angle_array = cls._quality_object.maximum_angle(mesh).array() + maximum_angle_list: np.ndarray = comm.gather(maximum_angle_array, root=0) + if comm.rank == 0: + maximum_angle_list = np.concatenate(maximum_angle_list, axis=None) + else: + maximum_angle_list = np.zeros(1) + + return maximum_angle_list + @classmethod def min_maximum_angle(cls, mesh: fenics.Mesh) -> float: r"""Computes the minimal quality measure based on the largest angle. @@ -361,7 +421,16 @@ def min_maximum_angle(cls, mesh: fenics.Mesh) -> float: The minimum value of the maximum angle quality measure. """ - quality: float = np.min(cls._quality_object.maximum_angle(mesh).array()) + maximum_angle_list = cls._maximum_angle(mesh) + comm = fenics.MPI.comm_world + + if comm.rank == 0: + qual = float(np.min(maximum_angle_list)) + else: + qual = None + + quality: float = comm.bcast(qual, root=0) + return quality @classmethod @@ -387,9 +456,39 @@ def avg_maximum_angle(cls, mesh: fenics.Mesh) -> float: The average quality, based on the maximum angle measure. """ - quality: float = np.average(cls._quality_object.maximum_angle(mesh).array()) + maximum_angle_list = cls._maximum_angle(mesh) + comm = fenics.MPI.comm_world + + if comm.rank == 0: + qual = float(np.average(maximum_angle_list)) + else: + qual = None + + quality: float = comm.bcast(qual, root=0) + return quality + @classmethod + def _radius_ratios(cls, mesh: fenics.Mesh) -> np.ndarray: + r"""Computes the radius ratios of the mesh elements. + + Args: + mesh: The mesh, whose quality shall be computed. + + Returns: + The radius ratios of the mesh elements on process 0. + + """ + comm = fenics.MPI.comm_world + radius_ratios_array = fenics.MeshQuality.radius_ratios(mesh).array() + radius_ratios_list: np.ndarray = comm.gather(radius_ratios_array, root=0) + if comm.rank == 0: + radius_ratios_list = np.concatenate(radius_ratios_list, axis=None) + else: + radius_ratios_list = np.zeros(1) + + return radius_ratios_list + @classmethod def min_radius_ratios(cls, mesh: fenics.Mesh) -> float: r"""Computes the minimal radius ratio of the mesh. @@ -410,7 +509,16 @@ def min_radius_ratios(cls, mesh: fenics.Mesh) -> float: The minimal radius ratio of the mesh. """ - quality: float = np.min(fenics.MeshQuality.radius_ratios(mesh).array()) + radius_ratios_list = cls._radius_ratios(mesh) + comm = fenics.MPI.comm_world + + if comm.rank == 0: + qual = float(np.min(radius_ratios_list)) + else: + qual = None + + quality: float = comm.bcast(qual, root=0) + return quality @classmethod @@ -433,22 +541,31 @@ def avg_radius_ratios(cls, mesh: fenics.Mesh) -> float: The average radius ratio of the mesh. """ - quality: float = np.average(fenics.MeshQuality.radius_ratios(mesh).array()) + radius_ratios_list = cls._radius_ratios(mesh) + comm = fenics.MPI.comm_world + + if comm.rank == 0: + qual = float(np.average(radius_ratios_list)) + else: + qual = None + + quality: float = comm.bcast(qual, root=0) + return quality @classmethod - def min_condition_number(cls, mesh: fenics.Mesh) -> float: - r"""Computes quality based on the condition number of the reference mapping. + def _cell_condition_number(cls, mesh: fenics.Mesh) -> fenics.Function: + r"""Computes the condition number quality for each cell. This quality criterion uses the condition number (in the Frobenius norm) of the (linear) mapping from the elements of the mesh to the reference element. - Computes the minimum of the condition number over all elements. Args: mesh: The mesh, whose quality shall be computed. Returns: - The minimal condition number quality measure. + A fenics.Function of a piecewise constant function space which holds the + cell's condition number quality measure. """ function_space_dg0 = fenics.FunctionSpace(mesh, "DG", 0) @@ -488,61 +605,44 @@ def min_condition_number(cls, mesh: fenics.Mesh) -> float: cond.vector().vec().scale(np.sqrt(mesh.geometric_dimension())) cond.vector().apply("") - quality: float = cond.vector().vec().min()[1] - return quality + return cond @classmethod - def avg_condition_number(cls, mesh: fenics.Mesh) -> float: - """Computes quality based on the condition number of the reference mapping. + def min_condition_number(cls, mesh: fenics.Mesh) -> float: + r"""Computes quality based on the condition number of the reference mapping. This quality criterion uses the condition number (in the Frobenius norm) of the (linear) mapping from the elements of the mesh to the reference element. - Computes the average of the condition number over all elements. + Computes the minimum of the condition number over all elements. Args: mesh: The mesh, whose quality shall be computed. Returns: - The average mesh quality based on the condition number. + The minimal condition number quality measure. """ - function_space_dg0 = fenics.FunctionSpace(mesh, "DG", 0) - jac = ufl.Jacobian(mesh) - inv = ufl.JacobianInverse(mesh) + cond = cls._cell_condition_number(mesh) + quality: float = cond.vector().vec().min()[1] - options: List[List[Union[str, int, float]]] = [ - ["ksp_type", "preonly"], - ["pc_type", "jacobi"], - ["pc_jacobi_type", "diagonal"], - ["ksp_rtol", 1e-16], - ["ksp_atol", 1e-20], - ["ksp_max_it", 1000], - ] + return quality - dx = measure.NamedMeasure("dx", mesh) - lhs = ( - fenics.TrialFunction(function_space_dg0) - * fenics.TestFunction(function_space_dg0) - * dx - ) - rhs = ( - fenics.sqrt(fenics.inner(jac, jac)) - * fenics.sqrt(fenics.inner(inv, inv)) - * fenics.TestFunction(function_space_dg0) - * dx - ) + @classmethod + def avg_condition_number(cls, mesh: fenics.Mesh) -> float: + """Computes quality based on the condition number of the reference mapping. - cond = fenics.Function(function_space_dg0) + This quality criterion uses the condition number (in the Frobenius norm) of the + (linear) mapping from the elements of the mesh to the reference element. + Computes the average of the condition number over all elements. - _utils.assemble_and_solve_linear( - lhs, rhs, x=cond.vector().vec(), ksp_options=options - ) - cond.vector().apply("") + Args: + mesh: The mesh, whose quality shall be computed. - cond.vector().vec().reciprocal() - cond.vector().apply("") - cond.vector().vec().scale(np.sqrt(mesh.geometric_dimension())) - cond.vector().apply("") + Returns: + The average mesh quality based on the condition number. + """ + cond = cls._cell_condition_number(mesh) quality: float = cond.vector().vec().sum() / cond.vector().vec().getSize() + return quality diff --git a/docs/source/conf.py b/docs/source/conf.py index fa1ffd5e..c6f6c4a9 100755 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -23,7 +23,7 @@ author = "Sebastian Blauth" # The full version, including alpha/beta/rc tags -release = "1.8.9" +release = "1.8.10" # -- General configuration --------------------------------------------------- diff --git a/setup.cfg b/setup.cfg index f3ae453e..71b38049 100755 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = cashocs -version = 1.8.9 +version = 1.8.10 author = Sebastian Blauth author_email = sebastian.blauth@itwm.fraunhofer.de description = Computational Adjoint-Based Shape Optimization and Optimal Control Software diff --git a/tests/sm_mesh/fine.h5 b/tests/sm_mesh/fine.h5 index 13c5d42a..2cba8ad3 100644 Binary files a/tests/sm_mesh/fine.h5 and b/tests/sm_mesh/fine.h5 differ diff --git a/tests/sm_mesh/fine_boundaries.h5 b/tests/sm_mesh/fine_boundaries.h5 index 68496484..9682eb1e 100644 Binary files a/tests/sm_mesh/fine_boundaries.h5 and b/tests/sm_mesh/fine_boundaries.h5 differ diff --git a/tests/sm_mesh/fine_subdomains.h5 b/tests/sm_mesh/fine_subdomains.h5 index 881ef4f1..c04b3e0a 100644 Binary files a/tests/sm_mesh/fine_subdomains.h5 and b/tests/sm_mesh/fine_subdomains.h5 differ