landau de gennes NI transition on flat 2D

isotropic_nematic_transition/iso_nem.py

+import matplotlib.pyplot as plt
+import dolfin as df
+import numpy as np
+import progressbar
+import sys
+np.set_printoptions(threshold=sys.maxsize)
+from ufl.constantvalue import PermutationSymbol as epsilon
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+
+class ActiveNematicManifold(object):
+    def __init__(self, parameters):
+        for key in parameters:
+            setattr(self, key, parameters[key])
+            
+        mesh = df.RectangleMesh(df.Point(-self.Lx/2, -self.Ly/2), 
+                                     df.Point(self.Lx/2, self.Ly/2),
+                                     self.Nx, self.Nx)
+        self.mesh = mesh
+            
+        # scalar and tensor finite elements
+        scalar_element = df.FiniteElement('P', cell=mesh.ufl_cell(), degree=self.degree)
+        tensor_element = df.TensorElement('P', cell=mesh.ufl_cell(), degree=self.degree)
+        
+        # Q, Lagrange mutiplier for trace of Q, Lagrange multiplier for antisymmetric part of Q
+        element = df.MixedElement([tensor_element, scalar_element, scalar_element])
+
+        # Construct the mixed function space
+        self.function_space = df.FunctionSpace(mesh, element)
+
+        # functions for values at t and t+1 respectively
+        self.function0 = df.Function(self.function_space)
+        self.function  = df.Function(self.function_space)
+
+    def setup_initial_conditions(self):
+        # Q0 = df.Constant(((0.0, 0.0), (0.0, 0.0)))
+        Q0 = df.Expression((('sin(x[0])* 0.01', 'cos(x[1])* 0.01'), ('cos(x[1])* 0.01', '- sin(x[0])* 0.01')), degree=1)
+        Q0 = df.interpolate(Q0, self.function_space.sub(0).collapse())
+        
+        l0 = df.interpolate(df.Constant(0.0), self.function_space.sub(1).collapse())
+        L0 = df.interpolate(df.Constant(0.0), self.function_space.sub(2).collapse())
+        
+        self.function = df.Function(self.function_space)
+        df.assign(self.function0, [Q0, l0, L0])
+
+    def setup_weak_forms(self, Qi=None):
+
+        Q0, _, _ = df.split(self.function0)
+        Q,  l, L = df.split(self.function)
+        tQ, tl, tL = df.TestFunctions(self.function_space)
+    
+        
+        Qform = (df.inner(tQ, (Q - Q0)/self.timestep)
+                + self.mobility * df.inner((2 * self.A * Q + 3 * self.B * df.dot(Q, Q) + 4 * self.C * Q * df.inner(Q, Q)), tQ)
+                + df.inner(l * df.Identity(2), tQ)
+                + df.inner(df.dot(epsilon(2), L * df.Identity(2)), tQ)
+                )
+        
+        lform = df.inner(tl, df.tr(Q))
+        Lform = df.inner(df.dot(epsilon(2), tL * df.Identity(2)), df.skew(Q))
+
+        self.form = (Qform + lform + Lform) * df.dx
+
+        
+    def solve(self):
+        QFile = df.XDMFFile('%s_Q.xdmf' % self.timestamp)
+
+        self.setup_initial_conditions()
+
+        time = 0.0
+        Q0, l0, L0 = self.function0.split(deepcopy=True) 
+        QFile.write_checkpoint(Q0, 'nematic', time)
+
+        self.setup_weak_forms()
+
+        maxsteps = round(self.maxtime/self.timestep)
+
+        for steps in progressbar.ProgressBar()(range(1, maxsteps+1)):
+            time += self.timestep
+
+            df.solve(self.form == 0, self.function)
+            Q, l, L = self.function.split(deepcopy=True) 
+            QFile.write_checkpoint(Q, 'nematic', time, append=True)
+
+            df.assign(self.function0, self.function)
+            
+        QFile.close()
+
+if __name__ == '__main__':
+    import json, datetime
+    import os
+    assert os.path.isfile('parameters.json'), 'parameters.json file not found'  # assert that parameters.json is a valid file, otherwise 
+                                                                                # give an error message parameters.json file not found
+    
+    # load the parameters
+    with open('parameters.json') as jsonFile:                    
+        params = json.load(jsonFile)
+        
+    
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+    params['timestamp'] = timestamp
+        
+    anm = ActiveNematicManifold(params)
+    anm.solve()
+    
+    with open(params['timestamp'] + '_parameters.json', "w") as fp:
+        json.dump(params, fp, indent=4)    

isotropic_nematic_transition/parameters.json

+{
+    "degree": 1,
+    "Lx": 1.0,
+    "Ly": 1.0,
+    "Nx": 32,
+    "Ny": 32,
+    "K":  1.0,
+    "A":  -0.5,
+    "B":  1,
+    "C": 1.0,
+    "maxtime": 10.0,
+    "timestep": 0.1,
+    "mobility": 1.0
+}
\ No newline at end of file

isotropic_nematic_transition/viz_iso_nematic_transition.py

+import vedo as vd
+import dolfin as df
+import numpy as np
+import h5py
+import progressbar
+np.set_printoptions(threshold=np.inf)
+from tempfile import TemporaryDirectory
+import os
+
+def get_mesh_from_xdmf(h5File):
+    # Read mesh geometry from h5 file
+    h5 = h5py.File(h5File, "r")
+    geometry = np.array(h5['nematic/nematic_0/mesh/geometry'])
+    topology = np.array(h5['nematic/nematic_0/mesh/topology'])
+    h5.close()
+
+    # create a mesh
+    mesh = df.Mesh()
+    editor = df.MeshEditor()
+    editor.open(mesh, type='triangle', tdim=2, gdim=2)
+    editor.init_vertices(len(geometry))
+    editor.init_cells(len(topology))
+    for i in range(len(geometry)):
+        editor.add_vertex(i, geometry[i])
+    for i in range(len(topology)):
+        editor.add_cell(i, topology[i])
+    editor.close()
+
+    return mesh
+
+import argparse, json
+
+# load the parameters
+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)
+
+gdim = 2
+xdmf = '%s_Q.xdmf' % params['timestamp']
+mesh = get_mesh_from_xdmf(xdmf.replace('xdmf', 'h5'))
+maxtime = params['maxtime']
+timestep = params['timestep']
+maxsteps = round(maxtime/timestep)
+
+n_cmap_vals = 8
+scalar_cmap = 'viridis'
+times = timestep * np.arange(maxsteps+1)
+
+print('Reading polarity....')
+nematic = np.zeros((len(times), mesh.num_vertices(), gdim, gdim))
+TFS = df.TensorFunctionSpace(mesh, 'P', 1)
+Q = df.Function(TFS)
+
+with df.XDMFFile(xdmf) as nFile:
+    for steps in progressbar.ProgressBar()(range(maxsteps + 1)):
+        nFile.read_checkpoint(Q, 'nematic', steps)
+        Qvals = Q.compute_vertex_values(mesh)
+        Qvals = Qvals.reshape(gdim, gdim, int(Qvals.shape[0]/gdim**2)).T
+        nematic[steps] = Qvals.transpose((0, 2, 1))  # exchange last two dimensions
+
+print(np.trace(nematic[5, 10]), np.trace(nematic[2, 1]))
+print( (nematic[10, 100] - np.transpose(nematic[10, 100]))/2)
+
+raise SystemExit
+Qmag = np.sqrt(np.einsum('tijk,tijk->ti', nematic, nematic)) # for each time and mesh point this is sqrt(Qjk Qjk)
+
+Qmin, Qmax = np.min(Qmag), np.max(Qmag)
+plotter = vd.plotter.Plotter(axes=4)
+poly = vd.utils.buildPolyData(mesh.coordinates(), mesh.cells())
+scalar_actor = vd.mesh.Mesh(poly)
+# scalar_actor.compute_normals(points=True, cells=True)  
+scalar_actor.pointdata['nematic'] = Qmag[0]
+scalar_actor.cmap('viridis', Qmag[0], vmin=Qmin, vmax=Qmax)
+scalar_actor.add_scalarbar(title="Qmag")
+plotter += scalar_actor
+
+director = np.zeros((len(times), mesh.num_vertices(), gdim))
+
+for steps in range(len(times)):
+    for i in range(mesh.num_vertices()):
+        l, v = np.linalg.eig(np.eye(gdim) + nematic[steps, i])  # why added identity?
+        idx = np.argmax(l)
+        director[steps, i] = np.real(v[idx]) * Qmag[steps, i]
+
+points = mesh.coordinates()
+
+scale = 0.1
+startPoints = points - scale * director[0]/2
+endPoints = points + scale * director[0]/2
+lines = vd.shapes.Lines(startPoints, endPoints, lw=3, c='k')
+plotter += lines
+
+def update(idx):
+    global plotter, lines
+    scalar_actor.pointdata['nematic'] = Qmag[idx]
+    scalar_actor.cmap('viridis', Qmag[idx], vmin=Qmin, vmax=Qmax)
+    plotter.remove(lines)
+    startPoints = points - scale * director[idx]/2
+    endPoints = points + scale * director[idx]/2
+    lines = vd.shapes.Lines(startPoints, endPoints, lw=3, c='k')
+    plotter += lines
+
+def slider_update(widget, event):
+    value = widget.GetRepresentation().GetValue()
+    idx = (abs(times-value)).argmin()
+    update(idx)
+    
+slider = plotter.add_slider3d(slider_update,
+                        pos1=(0.1, 0.9, 0.1), pos2=(0.9, 0.9, 0.1),
+                        xmin=times[0], xmax=times.max(), 
+                        value=times[0], title=r'$t$')
+slider.GetRepresentation().SetLabelFormat('%0.1f')
+
+
+
+vd.show(scalar_actor, interactive=False)
+
+FPS = 5
+movFile = '%s_nematic.mov' % params['timestamp']
+fps = float(FPS)
+command = "ffmpeg -y -r"
+options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
+tmp_dir = TemporaryDirectory()
+get_filename = lambda x: os.path.join(tmp_dir.name, x)
+for tt in range(len(times)):
+    idx = (abs(times-times[tt])).argmin()
+    update(idx)
+    slider.GetRepresentation().SetValue(times[tt])
+    fr = get_filename("%03d.png" % tt)
+    vd.screenshot(fr)
+os.system(command + " " + str(fps)
+            + " -i " + tmp_dir.name + os.sep
+            + "%03d.png " + options + " " + movFile)
+tmp_dir.cleanup()
+
+