1

I have a (discontinous) function f from R^n to {0, 1}, provided by the user. E.g. something like this "circle inquality":

def charfunc(x):
    # assume x is a numpy array
    return x[0]**2 + x[1]**2  < 1

I want to find the border of the region(s) where f == 1 holds by sampling function values on a regular grid and refine that grid (iteratively) in regions where the value changes. I think that should to a be not so exotic problem and I bet that it has already been solved in the python universe. However I only found complicated PDE/FEM-Packages and NDTAMR. The former seem too much overhead. NDTAMR looks very promising but delivers strange results (at least with the default parameters from the examples):

Original Circle Circle after refinement

There are cells refined which are not close to the border, and some cells which are close, are not refined. Code for this example lives in this repo.

Question: How to use NDTAMR or some other Python package to get a mesh which is refined everywhere near the border but nowhere else?

cknoll
  • 2,130
  • 4
  • 18
  • 34

2 Answers2

1

I finally implemented it by myself. For simplicity the code lives in some other package, I already maintain, but should be usable independently with little or no effort:

https://github.com/TUD-RST/symbtools/blob/master/symbtools/meshtools.py

Some unittests might serve as substitute for yet missing docs.

enter image description here

cknoll
  • 2,130
  • 4
  • 18
  • 34
-1

This should be easy enough to do by yourself.

  • For all the corners of your quads, check if they are inside or outside the circle.
  • Refine those quads which have nodes both inside and outside, and compute the status for the new points.

I might add this kind of border-finding is not typical in numerical analysis. Usually your domain is specified by a continuous function that is negative inside and positive outside of the domain, i.e.,

def f(x):
    return x[0] ** 2 + x[1] ** 2 - 1

along with the gradient

def grad(x):
    return 2 * x

(The gradient can also be computed automatically.) With this info, starting from any point, you can then iterate towards the boundary.

Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
  • That general formulation is indeed very simple. However, implementing the details for n dimensions and general function seems to be nontrivial. Thanks for the hint, that my problem is not typical. The question now explicitly mentions the discontinuous character of the given function. Using the gradient is not an option. – cknoll Dec 20 '19 at 14:23