3D Mesh generated from imported surface volume exhibits irregularities

mesh
1

I have extracted a brain surface from MRI data using a MatLab program (MRIcroS) into an stl file. I convert the resultant stl file to an ascII stl file format within python to load into mshr using Surface3D(). Upon meshing the volume within FEniCS, the mesh produced has 7 dense areas of extremely small cells on the surface of the mesh. (The interior mesh looks beautiful.). There is no such behavior in either stl file. I have attempted to use the smooth_boundary() attribute of the Mesh class to smooth the surface vertices of the mesh, but get the following error: TypeError: smooth_boundary(): incompatible function arguments. The following argument types are supported:
1. (self: dolfin.cpp.mesh.Mesh, arg0: int, arg1: bool) -> None
This error occurs even for a very simple mesh, such as:
import dolfin
mesh = dolfin.UnitSquareMesh(4,4)
mesh.smooth_boundary()
Am I using this attribute incorrectly?
I have successfully used the smooth() attribute on the brain mesh, but it has a limited effect on the surface mesh. Would appreciate any insight into using smooth_boundary() or handling the import of the surface file for a better mesh outcome. The relevant part of my program is below.
Thank you, Lori

from dolfin import *
from mshr import *

#—Build 3D mesh from surface extracted from MRI Data—

Load surface from file

geometry = Surface3D(‘Mouse_Mesh_03_sm12_asc.stl’)
domain = CSGCGALDomain3D(geometry)

Generate mesh preserving surface

mesh = generate_mesh(domain, 160)
mesh.smooth_boundary(20)
print(“max mesh size =”, mesh.hmax())
print(“min mesh size =”, mesh.hmin())
print(“number of cells =”, mesh.num_cells())
print(“number of vertices =”, mesh.num_vertices())

save mesh to file

mfile = File(‘mesh_sm12_160_mshrsmb1.pvd’)
mesh_file = File(‘mesh_sm12_160_mshrsmb1.xml.gz’)
mfile << mesh
mesh_file << mesh

2

Following the documentation of this function in C++:

  /// Smooth boundary vertices of mesh by local averaging.
    ///
    /// @param  num_iterations (std::size_t)
    ///         Number of iterations to perform smoothing,
    ///         default value is 1.
    ///
    /// @param  harmonic_smoothing (bool)
    ///         Flag to turn on harmonics smoothing, default
    ///         value is true.
    void smooth_boundary(std::size_t num_iterations=1,
                         bool harmonic_smoothing=true);
3

Dokken,
Thank you for sharing the code. From looking at it, I’m interpreting there are two parameters, which each have default values, and the ‘void’ indicates it does not return anything. I assume it performs the the smoothing operation on the mesh, so that the mesh is changed after the code is executed, but it does not return anything. I have tried:

  1. mesh.smooth_boundary(), invoking the default values
  2. mesh.smooth-boundary(20), specifying a value for the first parameter
  3. mesh.smooth-boundary(num_iterations=20), specifying a value for the first parameter
  4. mesh.smooth-boundary(20, True), specifying a value for each parameter
    They all give the same error reported below. Based on the code below, I do not understand why I would get this error. Do you have any ideas?
    Thank you,
    Lori
4

A minimal example turning on the log level:

set_log_level(LogLevel.PROGRESS)
mesh = UnitSquareMesh(10,10)

mesh.smooth_boundary(10,True)

Gives the output

Smoothing boundary of mesh: <Mesh of topological dimension 2 (triangles) with 121 vertices and 200 cells, ordered>
Smoothing mesh: <Mesh of topological dimension 1 (intervals) with 40 vertices and 40 cells, ordered>
Mesh smoothing repeated 10 times.
PETSc Krylov solver starting to solve 121 x 121 system.
PETSc Krylov solver starting to solve 121 x 121 system.

to reproduce this, you can use docker, this is the docker command (and image) I used

docker run -ti -v $(pwd):/home/fenics/shared -w /home/fenics/shared/ --rm quay.io/fenicsproject/dev:latest
5

Thank you Dokken,
The minimal example you give below worked great using your Docker container. If I change the unit square (2D mesh) to a unit cube (3D mesh), then I get the error below. smooth_boundary() does not work in 3-dimensions? Is there another way to smooth the 2D surface mesh of my 3D mesh?
Best,
Lori

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

6

If you had looked at the whole error message, it is clear that smooth_boundary only works in 2D.

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to find normal.
*** Reason:  Normal vector is not defined in dimension 3 (only defined when the triangle is in R^2.
*** Where:   This error was encountered inside TriangleCell.cpp.
*** Process: 0
***

In general, I do not recommend using mshr to generate meshes, as this module is no longer maintained. I would rather recommend using GMSH.

7

Thank you for the suggestion.