# Mesh with specified inside geometry

**URL:** <https://fenicsproject.discourse.group/t/mesh-with-specified-inside-geometry/1404>\
**Category:** mesh\
**Created:** [August 26, 2019, 2:02pm UTC](https://fenicsproject.discourse.group/t/mesh-with-specified-inside-geometry/1404 "2019-08-26T14:02:15Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Holzklotz](https://avatars.discourse-cdn.com/v4/letter/h/fbc32d/32.png) [@Holzklotz](https://fenicsproject.discourse.group/u/Holzklotz)\
**Post date:** [August 26, 2019, 2:02pm UTC](https://fenicsproject.discourse.group/t/mesh-with-specified-inside-geometry/1404/1 "2019-08-26T14:02:15Z")

</div>

I want to define a mesh with a subdomain. This is currently done by the following MWE and results in the plot below.

In my application I want to have a finer mesh around the boundary of the inner subdomain. One row of finer cells is enough. My idea of defining two extra rectangles (helper1 and helper2 in the code) is ignored, if I use the simple approach of adding helper1/2 to the domain.

**Question:** How do I force fenics/mshr to use my specified geometry in the inside of the domain?

Further thoughts:

- helper1/2 are ignored because the + operator is defined as the union. So in a mathematically sense the result is correct. The inner rectangle can be ignored because they are completely surrounded by the outer rectangle. Is there another command to combine/add two geometries but not as a union?
- I tried to use the set\_subdomain method, but it is restricted in such a way that a subdomain number of 0 can not be assigned. So a workaround would to be use subdomain ids \>0.

**MWE**

```
import fenics as fe
import mshr as mshr
    
#define geometry
domain = mshr.Rectangle( fe.Point(-1,-1), fe.Point(1,1) )
core = mshr.Rectangle( fe.Point(-0.5,-0.5), fe.Point(0.5,0.5) )
# this does not work
helper1 = mshr.Rectangle( fe.Point(-0.6,-0.6), fe.Point(0.6,0.6) )
helper2 = mshr.Rectangle( fe.Point(-0.4,-0.4), fe.Point(0.4,0.4) )
domain += helper1 + helper2
# mark subdomains
domain.set_subdomain(1, core)
# create mesh and meshfunction
mesh = mshr.generate_mesh(domain, 1)    
marker_subdomain = fe.MeshFunction('size_t', mesh, 2, mesh.domains())   
#plot
fe.plot(marker_subdomain)
fe.plot(mesh)

```

**Plot for the desired outcome** the red rectangles are missing in the MWE  
 ![Figure_2](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/da70953ae8313904e953978fa6036caf37a4f68b.png)

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [August 26, 2019, 2:16pm UTC](https://fenicsproject.discourse.group/t/mesh-with-specified-inside-geometry/1404/2 "2019-08-26T14:16:27Z")

</div>

Hi,  
For this sort of advanced mesh generation it is best to use a dedicated mesher such as Gmsh for instance. `mshr` is not designed to give control to local mesh sizes.

---

<div class="post-metadata">

**Author:** ![Holzklotz](https://avatars.discourse-cdn.com/v4/letter/h/fbc32d/32.png) [@Holzklotz](https://fenicsproject.discourse.group/u/Holzklotz)\
**Post date:** [August 28, 2019, 12:25pm UTC](https://fenicsproject.discourse.group/t/mesh-with-specified-inside-geometry/1404/3 "2019-08-28T12:25:07Z")

</div>

I found a workaround.

First define a subdomain for each geometry wanted.

```
domain.set_subdomain( 1, core_outer)
domain.set_subdomain(100, core)
domain.set_subdomain(101, core_inner)

```

I use number ranges of 0-99 which will later be mapped to 0.  
100-199 -\> 1  
200-299 -\> 2 and so on.

Create MeshFunction as usual

```
marker_subdomain = fe.MeshFunction('size_t', self.mesh, 2, self.mesh.domains()) 

```

And now ‘manually’ change the values.

```
temp = marker_subdomain.array()
temp = np.floor( temp/100 ).astype( np.uint64 )
marker_subdomain.set_values(temp.tolist())
```