working remeshing for a fixed domain diffusion problem

remeshing_diffusion_problem/191223-101432_parameters.json

+{
+    "resolution": 10,
+    "L": 1.0,
+    "D": 1,
+    "timestep": 0.1,
+    "maxtime": 0.5,
+    "timestamp": "191223-101432"
+}
\ No newline at end of file

remeshing_diffusion_problem/compare.py

+import dolfin as df
+import numpy as np
+np.set_printoptions(precision=16)
+
+c1File = df.XDMFFile('171223-155711_concentration.xdmf')
+c2File = df.XDMFFile('171223-155743_concentration.xdmf')
+
+mesh = df.IntervalMesh(128, -0.5, 0.5)
+times = np.arange(0.0, 1.1, 0.1)
+SFS = df.FunctionSpace(mesh, 'P', 1)
+
+for steps in range(len(times)):
+    c1 = df.Function(SFS)
+    c2 = df.Function(SFS)
+    c1File.read_checkpoint(c1, 'concentration', steps)
+    c2File.read_checkpoint(c2, 'concentration', steps)
+    print(df.assemble((c2 - c1)**2 * df.dx(mesh))**0.5)
+
+c1File.close()
+c2File.close()

remeshing_diffusion_problem/parameters.json

+{
+    "resolution": 10,
+    "L": 1.0, 
+    "D": 1,
+    "timestep": 0.1,
+    "maxtime": 0.5
+
+}   
\ No newline at end of file

remeshing_diffusion_problem/remeshing_diffusion.py

+import numpy as np
+import dolfin as df
+import progressbar
+
+np.set_printoptions(precision=16)
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+class diffusion(object):
+    def __init__(self, parameters):
+        # read in parameters
+        for key in parameters:
+            setattr(self, key, parameters[key])
+
+        self.mesh = df.IntervalMesh(self.resolution, 0.0, self.L)
+
+    def setup(self):
+        self.function_space = df.FunctionSpace(self.mesh,'P', 1)
+        self.c = df.Function(self.function_space)
+        self.tc = df.TestFunction(self.function_space)
+        self.c0 = df.Function(self.function_space)
+        
+    def setup_initial_conditions(self):
+        self.c0.interpolate(df.Expression('1 + 0.2 * cos(pi*x[0]/L)', pi=np.pi, L=self.L, degree=1))
+        print(self.c0.compute_vertex_values(self.mesh))
+
+    def setup_weak_forms(self):
+        self.form = (df.inner((self.c - self.c0)/self.timestep, self.tc) 
+            + self.D * df.inner(df.nabla_grad(self.c), df.nabla_grad(self.tc))
+            ) * df.dx
+        
+    def refine_mesh(self, mesh):
+        ## if selective refinement
+        # cell_markers = df.MeshFunction('bool', mesh, mesh.topology().dim())
+        # cell_markers.set_all(False)
+        # for cell in df.cells(mesh):
+        #     if cell.h() > 0.0:
+        #         cell_markers[cell] = True
+        # mesh = df.refine(mesh, cell_markers)
+
+        return df.refine(mesh)
+
+    def solve(self):
+
+        # function space and functions
+        self.setup()
+
+        # create file if it doesn't exist, open if it exists
+        self.cFile = df.XDMFFile('%s_concentration.xdmf' % self.timestamp)
+        
+        # time-variables
+        self.time = 0.0
+        maxsteps = round(self.maxtime/self.timestep)
+
+        # initial conditions
+        self.setup_initial_conditions()
+        self.cFile.write_checkpoint(self.c0, 'concentration', self.time) # save initial condition
+
+        for steps in range(1, maxsteps + 1):
+            # save initial condition in a new function to be used later
+            oldc = df.Function(self.function_space)
+            oldc.assign(self.c0)
+           
+           # refine mesh
+            self.mesh = self.refine_mesh(self.mesh)
+            # redefine function space and functions
+            self.setup()
+            # project old iniial condition onto new function space
+            self.c0 = df.project(oldc, self.function_space)
+            # define weak forms on new mesh
+            self.setup_weak_forms()
+
+            # solve
+            df.solve(self.form == 0, self.c)
+            # save
+            self.cFile.write_checkpoint(self.c, 'concentration', self.time, append=True)
+            # update
+            self.c0.assign(self.c)
+             
+        self.cFile.close()
+
+if __name__ ==  '__main__':
+    import json, datetime
+    with open('parameters.json') as jsonfile:
+        parameters = json.load(jsonfile)
+    
+    timestamp = datetime.datetime.now().strftime('%d%m%y-%H%M%S')
+    parameters['timestamp'] = timestamp
+
+    diff = diffusion(parameters)
+    diff.solve()
+
+    with open(parameters['timestamp']+'_parameters.json', 'w') as fp:
+        json.dump(parameters, fp, indent=4)
+  
+

remeshing_diffusion_problem/viz_remeshing_diffusion.py

+import dolfin as df
+import numpy as np
+import os
+from matplotlib.widgets import Slider
+import progressbar
+import matplotlib.pyplot as plt
+import h5py
+
+def visualize(params, DIR=''):
+    savesteps = round(params['maxtime']/params['timestep'])
+    times = np.arange(savesteps+1) * params['timestep']
+    
+    # read mesh geometry and topology from h5 file
+    h5 = h5py.File(os.path.join(DIR, '%s_concentration.h5' % params['timestamp']), 'r')
+
+    #NOTE: geometry and topology are lists of length len(times)
+    #NOTE geometry[i][:, 0] is an array of dimension (num_vertices,) for time step i
+    #NOTE topology[i] is an array of dimension (num_cells, 2) for time step i 
+      
+    geometry = []
+    topology = []
+    for i in progressbar.progressbar(range(len(times))):
+        geometry.append(np.array(h5['concentration/concentration_%d/mesh/geometry' %i])[:, 0])
+        topology.append(np.array(h5['concentration/concentration_%d/mesh/topology' %i]))
+
+    # create mesh at every timestep with this geometry and topology
+    meshes = []
+    for k in range(len(times)):
+        # Initialize mesh and mesh editor
+        mesh = df.Mesh()
+        editor = df.MeshEditor()
+        editor.open(mesh, "interval", 1, 1)
+
+        # Initialize vertices and cells
+        editor.init_vertices(len(geometry[k]))
+        editor.init_cells(len(topology[k]))
+
+        # Add vertices
+        for j in range(len(geometry[k])):
+            editor.add_vertex(j, [geometry[k][j]])
+
+        # Add cells
+        for j in range(len(topology[k])):
+            editor.add_cell(j, topology[k][j])
+
+        # Close editor
+        editor.close()
+
+        # append to list
+        meshes.append(mesh)
+
+    # read field values
+    c = [] 
+    cFile   = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp']))
+
+    # Reading data
+    print('Reading data...')
+    for steps in progressbar.progressbar(range(savesteps+1)):
+        SFS = df.FunctionSpace(meshes[steps], 'P', 1)
+        ci = df.Function(SFS)    
+        cFile.read_checkpoint(ci, 'concentration', steps)        
+        c.append(ci.compute_vertex_values(meshes[steps]))
+    cFile.close()
+
+    # plot
+    plt.rcParams.update({'font.size': 15})
+
+    fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    axes.set_xlabel(r'$x$')
+    axes.set_ylabel(r'$c$')
+
+    # spatial coordinate
+    x = []
+    for i in range(len(times)):
+        x.append(meshes[i].coordinates()[:,0])
+
+    cplot, = axes.plot(x[0], c[0], ms=5, marker = '.', linestyle = 'None')
+    
+    def update(value):
+        ti = np.abs(times-value).argmin()
+        cplot.set_xdata(x[ti])
+        cplot.set_ydata(c[ti])
+        plt.draw()
+        
+    sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+    slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+            valinit=min(times), valfmt='%3.1f',
+            fc='#999999')
+    slider.drawon = False
+    slider.on_changed(update)
+    plt.show()
+
+if __name__ == '__main__':  
+    import argparse, json
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
+    args = parser.parse_args()
+
+    with open(args.jsonfile) as jsonFile:
+        params = json.load(jsonFile)
+
+    visualize(params, DIR=os.path.dirname(args.jsonfile))