On_boundary Clarification

[tdiorio]

Hey everyone,

I’ve been running into issues defining boundaries on a complex geometry (2D Brain slice) using the on_boundary function. To see the root of my issue, I’m running a test case using a 2D hole in square plate with linear elasticity (code below).
rect = Rectangle(Point(0,0), Point(1,1))
circ=Circle(Point(0.5,0.5), 0.25, 100)
domain = rect-circ
mesh1 = generate_mesh(domain, 20)

When defining the boundaries, I am using a combination of x[0] and x[1] inequalities and calling the on_boundary function:
class Inner(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0]>0.2 and x[0] < 0.8 and x[1] < 0.8 and x[1]> 0.2

However, this seems to be returning True for more than just the boundary elements; it grabs everything within the specified inequality box.

Is anyone able to provide some clarification on how the above boundary definition could be improved to just return True for the boundary elements that lie within the specified domain? Or even just some clarity on how the on_boundary function is working?

Though I’m using the simple example of a hole-in-plate, this is ultimately to help me identify the boundaries of a complex geometry (2D Brain) shown below:

Thank you!
Tyler

[bleyerj]

Hi, here is what I get with

V = FunctionSpace(mesh1, "CG", 1)
u = Function(V)
bc = DirichletBC(V, Constant(1), Inner())
bc.apply(u.vector())
plot(u)

so on_boundary behaves as expected, it gives only facets belonging to the boundary

[bleyerj]

As regards your ultimate goal, I imagine that you will most likely get your brain mesh from an external mesher and not mshr which has only limited capabilities. External meshers provide ways to mark different parts of the boundaries which can then be read directly in dolfin without the need to use mathematical definition of the boundary.

[tdiorio]

Thank you for your help, that quick code was great for me to check my boundaries using simple geometry.

In regards to the brain mesh, for the time being I’m unable to use my external mesher to mark boundary parts unfortunately. So I’ve utilized a similar method to above where I define boxed regions of interest then call the on_boundary function:

class Leftv(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and x[0]>-0.0075-1e-7 and x[0]<0.0018+1e-7 and x[1]>-0.019-1e-7 and x[1 <0.029+1e-7) or (on_boundary and x[0]>-0.02-1e-7 and x[0]<-0.0075+1e-7 and x[1]>-0.027-1e-7 and x[1]<0.036+1e-7)

When I go to check if the correct regions are being identified by applying the same code you provided:

V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
bc = DirichletBC(V, Constant(1), Leftv())
bc.apply(u.vector())
plot(u)

I get a graphic highlighting the boxed region of interest (left) rather than just the boundary elements that I am hoping to see (right: drawn in pink outline). Below are graphics of the plot(u) output as well as a diagram of my intended regions:

I apologize for the odd boundary definition, but is there something I am missing here? I’m not entirely understanding why the output of plot(u) is grabbing more than just the boundary elements.

[dokken]

Could you check mesh.geometry().dim()?

I have a hunch you have a three-dimensional mesh as input (aka the points in the mesh are of the form (x,y,z), and not just (x,y)). This will confuse dolfin.

Additionally, Im curious, which mesh generator are you using?

[tdiorio]

I actually had the same thought a bit ago, but in fact print(mesh.geometry().dim()) yields 2.

Also in the .xdmf file, I am seeing correctly that the .xdmf file is 2D as shown by the .xdmf script describing it below:

<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XY"><DataItem 
DataType="Float" Dimensions="18919 2" Format="HDF" 
Precision="8">SY_QuadBL_noZDone2D.h5:/data0</DataItem></Geometry><Topology 
NodesPerElement="3" NumberOfElements="19374" TopologyType="Triangle"><DataItem 
DataType="Int" Dimensions="19374 3" Format="HDF" 
Precision="8">SY_QuadBL_noZDone2D.h5:/data1</DataItem></Topology></Grid></Domain> 
</Xdmf>

Importantly though, the process of obtaining this 2D XDMF mesh was odd. I exported as STL from Hypermesh (3D coordinates), converted to XML using Meshio (still 3D coordinates), manually edited to 2D by removing all z-coordinates and formatting accordingly, then converted to XDMF using Meshio (2D coordinates now). It certainly sounds possible that something in these steps is causing troubles for dolfin when handling the mesh. I’m not entirely sure how to proceed and would appreciate any tips.

For mesh generation, I’m using Altair Hypermesh as our lab has a license for it and it features a nice GUI. Due to connection issues I’m temporarily locked out of remote access from our computer with Hypermesh while we have building restrictions (COVID-19). Otherwise I’d implement the boundary tagging in there as mentioned above. I’m hoping to try gmsh soon as I’ve been reading/hearing a lot about it from FEniCS forums.

[dokken]

One thing that can clearly simplify this is to directly convert the mesh from stl to xdmf and not go through VTK.

Another thing you can test is create a copy of the vtk file, remove almost all the cells (keep all the points in the geometry), and look at a subset for debugging. (The same can be done if you write your xdmffile with ascii encoding from meshio.)

[tdiorio]

Great, thank you for your help that was a very simple fix for the boundary issue. By going directly from stl to xdmf I was able to see that each boundary class was correctly identifying my regions of interest:

 V = FunctionSpace(mesh2, "CG", 1)
 u = Function(V)
 bc = DirichletBC(V, Constant(1), Rightv())
 #bcl = DirichletBC(V, Constant(1), Leftv())
 #bcf = DirichletBC(V, Constant(1), Fixed())
 bc.apply(u.vector())
 %matplotlib notebook
 plot(u)

However, I am now working in the 3D domain and would like to simplify to 2D if possible.

print(mesh2.geometry().dim())
3

Any further ideas on how I can get this planar mesh back to being truly 2D?

[dokken]

When you are converting your mesh with meshio, put remove the third column of the points in the mesh.

[tdiorio]

Thanks, I’ll give this a go. Really appreciate your help.

[tdiorio]

Update: A github contributor (keurfonluu) to Meshio provided a very simple method of removing the z-dimension during the conversion using the following line of code in Ubuntu terminal:

$ meshio-convert filename.stl filename.xdmf --input-format stl -z

This now allows me to correctly identifies my boundary facets (shown below in very faint coloring):

Thank you very much for your help, really appreciate it.

[nataliazon]

@bleyerj

External meshers provide ways to mark different parts of the boundaries

I’m looking for a tool like that, one which allows for marking different boundary regions on the mesh boundary. Specifically what I have is a liver segmentation with a vessel going through it, I want the outside of the liver to have a boundary condition, and the vessel “wall” a different boundary condition (to account for heat sink due to blood flow). Can you name/recommend any tool that would help with that? I tried SegmentMesher extension to Slicer3D but couldn’t figure out how to mark specific sub-boundaries

[bleyerj]

Hi,
sorry but I am not familiar with defining geometry based on image segmentation. I personally use Gmsh which provides a mechanism to mark boundaries but I am not sure how it can be used with image-based geometry.

[tdiorio]

Does Gmsh have the capability of marking boundaries on an imported geometry?

If so, a good pipeline might be to perform image segmentation, convert to a surface/geometry of interest, then mark boundaries and mesh with Gmsh before finally importing into the FEniCS environment.

As an example of image segmentation to geometry, I have used ITK Snap to get a surface mesh (2D brain slice in my case), then I use Geomagic DesignX to perform autosurfacing (mesh to geometry object). This usually provides a fairly clean surface but can get messy when dealing particularly complex shapes. Nevertheless, it provides a surface that you can clean up using software like Hypermesh or possibly Gmsh (unfamiliar with this).

[bleyerj]

Yes you can mark boundaries from imported geometries using Gmsh.

[Rafael_Ferro]

Esteemed,

I would like to take advantage of the topic and also the code made @bleyerj in post 2, where I use the code made by him also for an elastic analysis of a cantilever, however I insert 2 holes in the middle of the beam. See the error and the code below:

Error:

  File "/Users/.../2D_elasticity/2D_elasticity_withhole.py", line 77, in <module>
    solve(a == l, u, BCs)

  File "/Users/.../opt/anaconda3/envs/fenics2/lib/python3.8/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)

  File "/Users/.../opt/anaconda3/envs/fenics2/lib/python3.8/site-packages/dolfin/fem/solving.py", line 241, in _solve_varproblem
    problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,

  File "/Users/.../opt/anaconda3/envs/fenics2/lib/python3.8/site-packages/dolfin/fem/problem.py", line 59, in __init__
    cpp.fem.LinearVariationalProblem.__init__(self, a, L, u._cpp_object, bcs)

RuntimeError: 

*** -------------------------------------------------------------------------
*** 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 define linear variational problem a(u, v) = L(v) for all v.
*** Reason:  Expecting the boundary conditions to to live on (a subspace of) the trial space.
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  
*** -------------------------------------------------------------------------

Code:

#2D linear elasticity

## Introduction

from fenics import *
from mshr import *

L1 = 2.5
H = 1.
#Define geometry for background
x1 = 1.25
y1 = 0.25
x2 = 1.25
y2 = 0.75
r = 0.1

domain = Rectangle(Point(0, 0), Point(L1, H)) - Circle(Point(x1, y1), r) - Circle(Point(x2, y2), r)
         
mesh = Mesh(generate_mesh(domain, 200))

## Constitutive relation

def eps(v):
      return sym(grad(v))

E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):
       return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

## Variational formulation

rho_g = 1e-3
f = Constant((0, -rho_g))

V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx

## Resolution

def left(x, on_boundary):
      return near(x[0], 0.)

class Circle1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]>1.1 and x[0] < 1.4 and x[1] < 0.1 and x[1]> 0.4
    
class Circle2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]>1.1 and x[0] < 1.4 and x[1] < 0.6 and x[1]> 0.9

Vc = FunctionSpace(mesh, "CG", 1)
uc = Function(Vc)
bc1 = DirichletBC(Vc, Constant(1), Circle1())
bc1.apply(uc.vector())
bc2 = DirichletBC(Vc, Constant(1), Circle2())
bc2.apply(uc.vector())

bc_left = DirichletBC(V, Constant((0.,0.)), left)

BCs = [bc_left, bc1, bc2]

u = Function(V, name="Displacement")
solve(a == l, u, BCs)
plot(1e3*u, mode="displacement")

## Validation and post-processing

print("Maximal deflection:", -u(L,H/2.)[1])
print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3))
Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(sigma(u), Vsig))
print("Stress at (0,H):", sig(0, H))
file_results = XDMFFile("2D_elasticity_results_with_2holes.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)

Can anybody help me?

[bleyerj]

You should not use a different function space for the bc1 and bc2 boundary conditions as stated by the error message. Which component of the displacement would you like to set to 1 on the inner circle boundaries ? Use V.sub(0) or V.sub(1) depending on your choice, although this may seem strange bcs to me…

[Rafael_Ferro]

Thank you very much. V.sub (0) worked!

[NataliaKor]

Hello. For my work, I need mark subdomains in my sample, for this purpose I found the example from FEniCS tutorial (4.3 Defining subdomains for different materials). For my mesh I’m using Salome and get the mesh with two subdomains (I call it even and odd); magnum.msh convert my file to XML and shows two Domains name_ID = 1 for even and name_ID = 2 for odd.

From XML file I can see that I have 755 vertex index and number of nodes in my mesh is also 755.
I wrote a code for defining boundary:

k_0 = 1.1
k_1 = 2.1
tol = 1e-14
class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return (0 <= x[0] <= 40 + tol) or (60 - tol <= x[0] <= 100 + tol)

class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return (40 - tol <= x[0] <= 60 + tol) or (100 - tol <= x[0] <= 120 + tol)  
# Define conductivity for different subdomains
subdomains = CellFunction('size_t', mesh)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(subdomains, 0)
subdomain_1.mark(subdomains, 1)
# use the values of the mesh function materials to define the variable coefficient kappa 
class K(Expression):
    def __init__(self, subdomains, k_0, k_1, **kwargs):
        self.subdomains = subdomains
        self.k_0 = k_0
        self.k_1 = k_1
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 0:
            values[0] = self.k_0
        else:
            values[0] = self.k_1

kappa = K(subdomains, k_0, k_1, degree=1)

V = FunctionSpace(mesh, 'P', 1)

kappa_function = interpolate(kappa, V)

After i get the file with different kappa I can see that it marks nodes not in the range which I mention. For example node 485 has k = 1.1, but node 486 (which located in the same domain) has k = 2.1.
Lines from XML file:

<dof index="485" value="1.1000000000000001" cell_index="36" cell_dof_index="1"/>
<dof index="486" value="2.1000000000000001" cell_index="34" cell_dof_index="3"/>

This picture from Salome shows coordinates:

Could you please clarify what I did wrong in specifying boundary? I want that it marks all verices in even area with k = 2.1. Thanks in advance.

[dokken]

@NataliaKor

V = FunctionSpace(mesh, 'P', 1)

kappa_function = interpolate(kappa, V)

You are interpolating a function defined for each cell into a space whose degrees of freedom is at vertices.
You even instruct the expression to be defined at vertices:

@NataliaKor

kappa = K(subdomains, k_0, k_1, degree=1)

which is not well defined.
What should happen if a vertex is shared between subdomain 0 and 1?
For visualization, I would use

kappa = K(subdomains, k_0, k_1, degree=0)
V = FunctionSpace(mesh, 'DG', 0)
kappa_function = interpolate(kappa, V)