294
views
0
12 months ago by
Hello,

it appears that adaptation of a function to a refined mesh does not preserve the integral of the function.

The following code produces two different results for me:

print(str(assemble(u0*dx)))

refined_mesh = refine(mesh, cell_markers)
dx = Measure("dx",domain =refined_mesh)

print(str(assemble(u0*dx)))
​

Is this a typical behavior or am I doing something wrong?

Community: FEniCS Project
1
Not sure about the adapt method, but have you tried interpolating u0 onto the refined functionspace first? Also, providing a complete MWE would be appreciated.
written 12 months ago by pf4d
Yes, I tried this, but it also failed.

I put together a 'MWE' based on the Cahn-Hilliard Demo.
The first few Iterations, the integral value stays the same and afterwards it changes, but only if I use Mesh adaptation.
import random
from dolfin import *

set_log_level(WARNING)

# Class representing the intial conditions
class InitialConditions(Expression):
def __init__(self, **kwargs):
random.seed(2 + MPI.rank(mpi_comm_world()))
def eval(self, values, x):
values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)

# Model parameters
lmbda  = 1.0e-02  # surface parameter
dt     = 5.0e-06  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

# Create mesh and build function space
mesh = UnitSquareMesh(96, 96)
mh = MeshHierarchy(mesh)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mh.finest(), P1*P1)

# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)

# Compute the chemical potential df/dc
c = variable(c)
f    = 100*c**2*(1-c)**2
dfdc = diff(f, c)

# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

# Weak statement of the equations
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = NonlinearVariationalProblem(L, u, J=a)
solver =  NonlinearVariationalSolver(problem)

# Output file
file = File("output.pvd", "compressed")

# Step in time
t = 0.0
T = 50*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()

mh = MeshHierarchy(mesh)
#to access values
c0, mu0 = u0.split()

for i in range(2):
c1 = project(c0,FunctionSpace(mh.finest(),P1))
cell_markers_refine = CellFunction("bool", mh.finest())
cell_markers_refine.set_all(False)
for cell in  cells(mh.finest()):
point = cell.midpoint()
tmp = abs(c1([point.x(), point.y()]))
if (cell.h() > .01) and (tmp > 0.1 and tmp < 0.9):
cell_markers_refine[cell] = True

mh = mh.refine(cell_markers_refine)

dx = Measure("dx",domain = mh.finest())
ds = Measure("ds",domain = mh.finest())

u0        = u0.child()
u         = u.child()
problem   = problem.child()
solver    = NonlinearVariationalSolver(problem.leaf_node())

solver.solve()

c0, mu0 = u0.split()
print(str(assemble(c0*dx)))
# file << (u.split()[0], t)
​
written 12 months ago by Lukas O.