# Weird artifacts on DAGMC geometry plots

**URL:** <https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278>\
**Category:** User Support\
**Created:** [August 7, 2023, 1:30pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278 "2023-08-07T13:30:30Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![nlinden](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/nlinden/32/2897_2.png) [@nlinden](https://openmc.discourse.group/u/nlinden)\
**Post date:** [August 7, 2023, 1:30pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/1 "2023-08-07T13:30:30Z")

</div>

Hi everyone,

I have been using the DAGMC Cubit plugin to produce some DAGMC geometries that I could load into OpenMC. It seems to work well and I was able to successfully perform some eigenvalue calculations which yield sensible results.

However when trying to plot these geometries to check that the entire process went well, I noticed some weird artifacts on the resulting images. Since I could not trace these artifacts back to any of the DAGMC/Cubit operations I performed, I tried plotting the ready-made DAGMC geometries available in the openmc-notebooks repo. In particular I plotted the teapot running the following script:

```python
import openmc
import math
from matplotlib import pyplot as plt
import urllib.request

teapot_url = 'https://tinyurl.com/y4mcmc3u' # 29 MB
def download(url, filename="dagmc.h5m"):
    """
    Helper function for retrieving dagmc models
    """
    u = urllib.request.urlopen(url)
    
    if u.status != 200:
        raise RuntimeError("Failed to download file.")
    
    # save file as dagmc.h5m
    with open(filename, 'wb') as f:
        f.write(u.read())

download(teapot_url, "teapot.h5m")
dagmc_univ = openmc.DAGMCUniverse(filename="teapot.h5m")
geometry = openmc.Geometry(root=dagmc_univ)

iron = openmc.Material(name="iron")
iron.add_nuclide("Fe54", 0.0564555822608)
iron.add_nuclide("Fe56", 0.919015287728)
iron.add_nuclide("Fe57", 0.0216036861685)
iron.add_nuclide("Fe58", 0.00292544384231)
iron.set_density("g/cm3", 7.874)

water = openmc.Material(name="water")
water.add_nuclide('H1', 2.0, 'ao')
water.add_nuclide('O16', 1.0, 'ao')
water.set_density('g/cc', 1.0)
water.add_s_alpha_beta('c_H_in_H2O')
water.id = 41

materials = openmc.Materials([iron, water])

settings = openmc.Settings()
settings.batches = 10
settings.particles = 5000
settings.run_mode = "fixed source"

space = openmc.stats.Box(*geometry.bounding_box)
my_source = openmc.Source(space=space, domains=[water])
my_source.energy = openmc.stats.Discrete(x=[12.0,], p=[1.0,])

settings.source = my_source

mesh = openmc.RegularMesh() # could use from_domain here to avid setting lower_left and upper_right
mesh.dimension = (120, 1, 40)
mesh.lower_left = (-20.0, 0.0, -10.0)
mesh.upper_right = (20.0, 1.0, 4.0)

mesh_filter = openmc.MeshFilter(mesh)

pot_filter = openmc.CellFilter([1])
pot_tally = openmc.Tally()
pot_tally.filters = [mesh_filter, pot_filter]
pot_tally.scores = ['flux']

water_filter = openmc.CellFilter([5])
water_tally = openmc.Tally()
water_tally.filters = [mesh_filter, water_filter]
water_tally.scores = ['flux']

tallies = openmc.Tallies([pot_tally, water_tally])

model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies)

p = openmc.Plot()
p.basis = 'xy'
p.origin = (0.0, 0.0, 0.0)
p.width = (50., 50.)
p.pixels = (1000, 971)
p.color_by = 'cell'
p.colors = {iron: 'gray', water: 'blue'}
model.plots = openmc.Plots(plots=[p])
model.export_to_model_xml()
openmc.plot_geometry()

```

The resulting plot looks like this:

 ![plot_1](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/b/bfb95d0a82aef1bfb221a00d62d62caf358b487f.png)

Here we can clearly see the teapot’s body, spout, and handle, as well as the DAGMC graveyard. However we also see that there is some green tea bleeding into the implicit complement filled with void, we can also see some undefined white space.

Has anyone run into these problems before? I am wondering if these artifacts come from the plotting operation, or if there are some underlying issues with the CAD geometry.

Thanks,  
Nicolas

---

<div class="post-metadata">

**Author:** ![paulromano](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/paulromano/32/486_2.png) [@paulromano](https://openmc.discourse.group/u/paulromano)\
**Post date:** [August 7, 2023, 6:40pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/2 "2023-08-07T18:40:05Z")

</div>

@nlinden I’ve also run into this same issue recently. I don’t have a good idea of what’s going on yet. One thing that is slightly suspicious is that the artifacts lie at a 45 degree angle in the x-y plane. Internally in OpenMC’s plotting routines, it creates a particle at a given (x,y,z) location to determine what cell/material it corresponds to. It’s also necessary to define the particle’s direction, which is arbitrary but nevertheless right now the direction is assigned as (0.707, 0.707, 0), which corresponds to that 45 degree angle. So, that makes me suspect that it may be a plotting artifact but I don’t have a good explanation yet. Will let you know if I find out any more!

Pinging the DAGMC experts in case there are any more thoughts @pshriwise @Shimwell

---

<div class="post-metadata">

**Author:** ![Shimwell](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/shimwell/32/3876_2.png) [@Shimwell](https://openmc.discourse.group/u/Shimwell)\
**Post date:** [August 7, 2023, 11:08pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/3 "2023-08-07T23:08:00Z")

</div>

I think this teapot was made with the Cubit + Cubit plugin method which is the most mature method we have, but perhaps Patrick can confirm.

All I can offer is another way to view it

If you have Moab installed then you should have a command line tool called `mbconvert` also installed

You can convert it to a vtk file and load it up with paraview

`mbconvert dagmc.h5m dagmc.vtk`

This nice thing about this method is you see the actual mesh triangles that make up the geometry

---

<div class="post-metadata">

**Author:** ![gonuke](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/gonuke/32/1026_2.png) [@gonuke](https://openmc.discourse.group/u/gonuke)\
**Post date:** [August 15, 2023, 12:28pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/4 "2023-08-15T12:28:20Z")

</div>

We’ve also noticed this in some of our geometries. @camoreno may be able recall what we learned as we investigated

---

<div class="post-metadata">

**Author:** ![camoreno](https://avatars.discourse-cdn.com/v4/letter/c/13edae/32.png) [@camoreno](https://openmc.discourse.group/u/camoreno)\
**Post date:** [August 15, 2023, 8:44pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/5 "2023-08-15T20:44:47Z")

</div>

Yes, I ran into this issue when plotting stellarator geometries. Looking at the VTK in Paraview didn’t show any obvious geometry issues. In our experience, the artifacts always lie at a 45 degree angle in the XY plane but in the XZ plane, show parts of the geometry that should be out of plane (see the included pictures, material overlap shown in red). The orientation of the artifacts makes sense given @paulromano 's explanation of how the plotter works.

We also tested for other geometries and we found the same issue. However, when meshing the geometry via Gmsh, rather than Cubit, these artifacts disappeared, though strange material overlaps were found. Here’s a nice little breakdown of what we observed:

 ![Screenshot 2023-08-15 at 3.40.03 PM](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/8/87de303a72ab3fbf5377c4b6bb3db1a1c1e40c23.jpeg)

XZ plot of stellarator geometry:

 ![Screen Shot 2023-04-14 at 5.18.23 PM](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/3/3e72e9c8a45045f1626f0dd9099844f8f6056946.png)

We never ended up figuring out any solutions, though, and just hoped it was an artifact of the plotter.

---

<div class="post-metadata">

**Author:** ![camoreno](https://avatars.discourse-cdn.com/v4/letter/c/13edae/32.png) [@camoreno](https://openmc.discourse.group/u/camoreno)\
**Post date:** [August 17, 2023, 1:47pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/6 "2023-08-17T13:47:33Z")

</div>

Something else worth mentioning is that these geometries were made in OpenMC by reusing the same DAGMC universe. The universe is used four times, copying and rotating it, to make up each of the four quadrants. I think to some degree that explains why the artifacts “stop” at the end of each quadrant, corresponding to the boundary of each universe. When a full geometry is meshed and plotted in OpenMC, the artifacts appear differently:

 ![Screen Shot 2023-04-06 at 4.37.22 PM](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/5/5f08c50b5734d21c7781307e6910fa837fa07a70.png)

---

<div class="post-metadata">

**Author:** ![akoen](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/akoen/32/2832_2.png) [@akoen](https://openmc.discourse.group/u/akoen)\
**Post date:** [August 17, 2023, 11:04pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/7 "2023-08-17T23:04:56Z")

</div>

This issue has also been affecting our stellarator workflows.

The following script reproduces the issue with a simple layered torus meshed using both CAD-to-OpenMC and CAD-to-DAGMC. We’ve also observed the issue with CAD-to-OpenMC’s cadquery backend, so it seems not to be mesh backend specific.

Can also confirm that it seems to be mostly on 45 degree angles to the X-Y plane.

Source: [https://file.io/tlmxMXYyl7EB](https://file.io/tlmxMXYyl7EB) (not allowed to upload attachments as new user).

# Script (screenshot expands)

 ![CleanShot 2023-08-17 at 15.55.19](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/0/0d5cd7bfd8dc165c109db754a6ca2a74bce5a770.jpeg)

---

<div class="post-metadata">

**Author:** ![soha](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/soha/32/3658_2.png) [@soha](https://openmc.discourse.group/u/soha)\
**Post date:** [August 18, 2023, 4:28pm UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/8 "2023-08-18T16:28:21Z")

</div>

Adding these visuals to reflect @akoen 's comment above.

![CleanShot 2023-08-18 at 08.39.53](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/1/1d2ba8ae7005806d8bac24b6e2e7497a19eb210d.png)  
 ![CleanShot 2023-08-18 at 08.40.06](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/3/3c456341f494984964473ac834f46eeb09dc93a1.png)  
 ![CleanShot 2023-08-18 at 08.40.14](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/8/8062bbbc1a2f64832910ff344304161fd025b9df.png)

---

<div class="post-metadata">

**Author:** ![pshriwise](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/pshriwise/32/3863_2.png) [@pshriwise](https://openmc.discourse.group/u/pshriwise)\
**Post date:** [August 21, 2023, 7:46am UTC](https://openmc.discourse.group/t/weird-artifacts-on-dagmc-geometry-plots/3278/9 "2023-08-21T07:46:14Z")

</div>

Hi all!

Thanks for spurring investigation of these issues, everyone! I’ve seen it a few times myself and the problem has always been nagging at me for a while now.

* * *

As far as the teapot model goes, @nlinden, I happen to know that the water surface starts at z=0, sot he plot is slicing right on one of the planar triangle surfaces of the model. This causes problems for CSG geometry as well – the sample points for the image can’t quite decide what side of the surface to be on numerically. Changing the z value of the plot origin even slightly, say to `-0.0001`, fixes this problem. That said, the direction we use for plotting is zero in the z direction, which isn’t ideal as it increases the likelihood of problems with intersections on z-axis aligned planes, something pretty common in our geometries. That’s something we can easily update.

* * *

As far as issues where plot slices don’t line up with geometry planes, I did some digging and, while there may be more than one cause for these artifacts, here’s one (long-winded) explanation in the presence of unmerged coincident surfaces in the DAGMC model. The fix I implemented did seem to fix the problem in that case.

> <https://github.com/openmc-dev/openmc/issues/2660>
>
> \*\*TL;DR:\*\* Artifacts sometimes appear in the plotter in regions with unmerged su…rfaces in DAGMC models. The order in which we check volumes/cells during the \`find\_cell\` operation matters in these cases and we should defer checking DAGMC's implicit complement until the end.
> 
> \---
> \## The problem
> 
> "Streaks" sometimes appear in plotting images of DAGMC models. As in our CSG geometry, this can occur if the slicing plane is on the boundary of a region, but these artifacts also appear throughout volumes in DAGMC universes. This is an explanation as to why this occurs and a proposal for a fix.
> 
> !\[image\](https://github.com/openmc-dev/openmc/assets/4563941/2cdc817b-e611-46e8-bb1e-2165ecc67533)
> 
> Image from the Copenhagen Atomics open \[MSRE model\](https://github.com/openmsr/msre)
> 
> 
> I can confirm that there are unmerged, coincident surfaces in this model in this region that are causing these problems. It's possible that other models have different causes for these artifacts, but here's my explanation of why it occurs for unmerged, coincident surfaces.
> 
> \## DAGMC model description
> 
> Let's say a model is created with the following two volumes, 1 and 2.
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/bed4820e-d7d3-48d4-b07a-7412e083400e"\>
> 
> 
> Ideally, there would be a single surface, surface 1, between the two volumes. It would have a "forward" sense with respect to volume 1 and a "reverse" sense with respect to volume 2. This sense is used to determine whether a particle crossing is entering or exiting one of the volumes.
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/0e540d7b-6562-4627-80c3-a5368d82ff23"\>
> 
> 
> This shared surface is generated as part of the "imprint and merge" process in the CAD engine. Part of this process is to combine, or merge, coincident surfaces into one -- resulting in surface 1 in this case. In the case that this process fails, the model has two coincident surfaces, surfaces 1 and 2, between volumes 1 and 2. In this case, we get the following:
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/c4bcd327-0e5e-407f-af21-3948b57b8134"\>
> 
> 
> The surfaces are shown with some separation to ease visualization, but they are numerically coincident with one another. Surface 1 now has a "forward" sense with respect to volume 1 and nothing on the other side. Surface 2 has a "forward" sense with respect to volume 2 as well and nothing on the other side. When building the \[implicit complement\](https://svalinn.github.io/DAGMC/contribute/devtheory.html#implicit-complement) in DAGMC, the programatically constructed space claims the "nothing" side of these surfaces. A key thing to remember is that \_both\_ of these surfaces are now part of the IPC volume generated by 
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/1c07b283-13e5-4bb5-b332-1a72f90704f1"\>
> 
> 
> At first glance, it might appear that the implicit complement occupies the space between the two unmerged surfaces, but this isn't quite how it works in practice. 
> 
> \### The implicit complement from the perspective of volume 1:
> 
> In this case the implicit complement is bounded by surface 1.
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/7ca5f588-9862-4e1d-bb12-a4909629fdf9"\>
> 
> \### Implicit complement from the perspective of volume 2:
> 
> In this case the implicit complement is bounded by surface 2.
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/54b30567-87be-4a70-9545-a96a4edb5089"\>
> 
> \## How this affects the plotting process:
> 
> In OpenMC we check each volume in a DAGMC universe for point containment to establish the position of a particle. To perform this check, we fire a ray in any given direction from the particle's location in an attempt to find an \_exiting\_ intersection with the triangles of the DAGMC volume.
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/0f2d9054-c7f1-49cd-ba3c-9fe1eb0075c2"\>
> 
> 
> As we get into the subtleties of the problem, we'll recall that these surfaces are coincident and represent them as such in the image:
> 
> \<img width="677" alt="image" src="https://github.com/openmc-dev/openmc/assets/4563941/4fb6c845-d9e5-47fe-96f4-44dd4adcba57"\>
> 
> While looping over volumes searching for a successful point containment check, it's possible that the implicit complement may be checked before volume 1 and 2. This is where the problem occurs. Because both surfaces 1 and 2 are part of the implicit complement and are coincident, a nearest intersection with either surface is valid and results in two outcomes:
> 
> 1. The particle ray intersects with surface 1: Based on the image of the IPC from the perspective of volume 1, this is considered to be an \_entering\_ intersection as if the particle were moving from volume 1 to the IPC volume and therefore the point is considered to be outside the IPC. The search loop continues and eventually identifies that the particle is in volume 1 (only the triangles of surface 2 are part of this search, so there's no ambiguity in that case)
> 2. The particle ray intersects with surface 2: Based on the image of the IPC from the perspective of volume 2, this is considered to be an exiting intersection as if the particle were moving from the IPC to volume 2. The particle is considered to be in the IPC and the search loop stops \*\*resulting in the artifacts we see in the plotter\*\*.
> 
> \### A proposed solution
> 
> Instead of including the IPC volume in the initial search for a containing cell, we can first check all of the intentionally created volumes in the CAD model and fall back to a check on the IPC at the end. In this way we still check all cells in the DAGMC universe, but we do not inadvertently exit the loop too soon and miss a particle placement. An implementation of this fix can be found \[here\](https://github.com/pshriwise/openmc/tree/dagmc-find-cell-updates) and does seem to remove these artifacts. Note there is some special behavior depending on whether or not the \`DAGUniverse\` is the root universe of the geometry or not.
> 
> !\[image\](https://github.com/openmc-dev/openmc/assets/4563941/7871beaa-b12a-4c07-8b2d-780f7b3fa55d)
> 
> 
> These changes requires the \`DAGUniverse::find\_cell\` method to take ownership of the search loop (before we could rely on the \`Universe::find\_cell\` routine). That was nice, but I don't see a problem with having a universe containing a fundamentally different geometry type have its own search loop.
> 
> \---
> 
> I thought I'd tackle a couple of questions that could come to mind reading the post above.
> \### Why does this occur in some locations of a volume but not others? Why is it direction-dependent?
> 
> The intersection with coincident surfaces of triangles is going to come down to traversal order of the bounding volume hierarchy and floating point arithmetic during the triangle intersection calculation -- slight changes in the angle and position will result in selecting triangles from one surface or the other of the coincident surfaces. From the perspective of the ray tracer, a nearest intersection with either surface is valid.
> 
> \### How might this adversely affect transport?
> 
> During particle transport, as long as we start in the correct volume, the unmerged surfaces shouldn't have a negative affect due to how intersections are filtered. During point containment checks, DAGMC searches for the nearest \_exiting or entering\_ intersection. During transport, DAGMC performs a query to find the nearest \_exiting\_ intersection. So in the example above as a particle moves out of volume 1 it will:
> 
> 1. Intersect with surface 1
> 2. Cross into the IPC volume
> 3. Intersect with surface 2 at a zero distance (intersections with surface 1 are filtered b/c they would be entering intersections)
> 4. Cross into volume 2
> 
> So while there are additional steps involved, the unmerged surfaces shouldn't cause a problem with transport.
> 
> \### Why didn't we see this in testing?
> 
> The test models were created as analogs of eigenvalue problems. As a result the geometries are simple and don't contain unmerged surfaces. Going forward, it's worth considering adding some models with unmerged surfaces to ensure we can handle these cases as expected -- though this is arguably something that should be part of DAGMC testing as opposed to a downstream application.

**Before**

 ![image](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/1/1d998ab295432e09dace8ab7ae41e9c08e15845e.png)  
**After**  
 ![image](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/a/abd11ef959c4c9322f82e1c986c0402cb57fd496.png)

I’d be curious if the fix I’ve implemented would help with some of the cases in this thread. Is anyone willing/able to send some of their models my way or try out the branch I’ve linked in the issue above?

Cheers,

Patrick

Branch link here for reference: [GitHub - pshriwise/openmc at dagmc-find-cell-updates](https://github.com/pshriwise/openmc/tree/dagmc-find-cell-updates)
