I got a 3d tetraedrized space $\Omega$ and a Affine surface over the elements $\psi(x,y,z)=0$
I want to redefine the cuted elements maintaining the hierarchy between old elements and new elements (this is only a 2D example)
The idea is add new points $\psi \cap edges(\Omega)$ and in the next step add new edges for build a new elements but without to destroy the "father" elements.
then:
from dolfin import *
def psi(x, y, z):
#some level function
return x**2+y**2+z**2-1
omega = BoxMesh(-1, -1, -1, 1, 1, 1, 10, 10, 10)
cell_markers = CellFunction("bool", omega)
cell_markers.set_all(False)
coordinates = omega.coordinates()
for cell in cells(omega):
# the psi function cut some edge?
for edge in edges(cell):
v0 = psi(*coordinates[edge.entities(0)[0]])
v1 = psi(*coordinates[edge.entities(0)[1]])
if(v0*v1<0):
#if level function change the sign between the nodes then has a zero un that element
cell_markers[cell]=True
new_omega = refine(omega, cell_markers)
omega_file = File("omega.pvd")
omega_file << omega
new_omega_file = File("new_omega.pvd")
new_omega_file << new_omega
but I have not seen the difference between $\Omega$ and $\Omega_{\text{new}}$