More fixes to make Oasis work with FEniCS 2018.1.0

oasis/common/io.py

@@ -91,16 +91,14 @@ def save_tstep_solution_h5(tstep, q_, u_, newfolder, tstepfiles, constrained_dom
         # project or store velocity to vector function space
         for comp, tstepfile in six.iteritems(tstepfiles):
             if comp == "u":
-                V = q_['u0'].function_space()
-                # First time around create vector function and assigners
-                if not hasattr(tstepfile, 'uv'):
-                    tstepfile.uv = AssignedVectorFunction(u_)
+                # Create vector function and assigners
+                uv = AssignedVectorFunction(u_)
 
                 # Assign solution to vector
-                tstepfile.uv()
+                uv()
 
                 # Store solution vector
-                tstepfile.write(tstepfile.uv, float(tstep))
+                tstepfile.write(uv, float(tstep))
 
             elif comp in q_:
                 tstepfile.write(q_[comp], float(tstep))

oasis/problems/NSCoupled/Nozzle2D.py

@@ -26,7 +26,7 @@ def create_bcs(VQ, mesh, sys_comp, re_high, **NS_namespce):
     # Analytical, could be more exact numerical, different r_0
     u_maks = Q / (4. * r_0 * r_0 * (1. - 2. / pi))
     #inn = Expression(("u_maks * cos(sqrt(pow(x[1],2))/r_0/2.*pi)", "0"), u_maks=u_maks, r_0=r_0)
-    inn = Expression(("u_maks * (1-x[1]*x[1]/r_0/r_0)", "0"), u_maks=u_maks, r_0=r_0)
+    inn = Expression(("u_maks * (1-x[1]*x[1]/r_0/r_0)", "0"), u_maks=u_maks, r_0=r_0, degree=2)
 
     bc0 = DirichletBC(VQ.sub(0),    inn,  inlet)
     bc1 = DirichletBC(VQ.sub(0), (0, 0),  walls)
@@ -42,7 +42,7 @@ def pre_solve_hook(mesh, V, **NS_namespace):
     Outlet = AutoSubDomain(outlet)
     Walls = AutoSubDomain(walls)
     Centerline = AutoSubDomain(centerline)
-    facets = FacetFunction('size_t', mesh, 0)
+    facets = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
     Inlet.mark(facets, 1)
     Outlet.mark(facets, 2)
     Walls.mark(facets, 3)

oasis/problems/NSCoupled/SkewedFlow.py

@@ -20,7 +20,7 @@ def problem_parameters(NS_parameters, **NS_namespace):
 def create_bcs(V, VQ, mesh, **NS_namespace):
     # Create inlet profile by solving Poisson equation on boundary
     bmesh = BoundaryMesh(mesh, 'exterior')
-    cc = CellFunction('size_t', bmesh, 0)
+    cc = MeshFunction('size_t', bmesh, bmesh.topology().dim(), 0)
     ii = AutoSubDomain(inlet)
     ii.mark(cc, 1)
     smesh = SubMesh(bmesh, cc, 1)
@@ -32,7 +32,7 @@ def create_bcs(V, VQ, mesh, **NS_namespace):
           bcs=[DirichletBC(Vu, Constant(0), DomainBoundary())])
 
     # Wrap the boundary function in an Expression to avoid the need to interpolate it back to V
-    class MyExp(Expression):
+    class MyExp(UserExpression):
         def eval(self, values, x):
             try:
                 values[0] = su(x)

oasis/problems/NSfracStep/Channel.py

@@ -116,9 +116,9 @@ def pre_solve_hook(V, u_, mesh, AssignedVectorFunction, newfolder, MPI,
     stats = ChannelGrid(V, [Nx, Ny + 1, Nz], [tol, -Ly / 2., -Lz / 2. + tol], [[1., 0., 0.], [
                         0., 1., 0.], [0., 0., 1.]], [Lx - Lx / Nx, Ly, Lz - Lz / Nz], statistics=True)
 
-    # Create FacetFunction to compute flux
+    # Create MeshFunction to compute flux
     Inlet = AutoSubDomain(inlet)
-    facets = FacetFunction('size_t', mesh, 0)
+    facets = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
     Inlet.mark(facets, 1)
     normal = FacetNormal(mesh)
 
@@ -137,7 +137,7 @@ def walls(x, on_bnd):
     return bcs
 
 
-class RandomStreamVector(Expression):
+class RandomStreamVector(UserExpression):
     random.seed(2 + MPI.rank(MPI.comm_world))
 
     def eval(self, values, x):

oasis/problems/NSfracStep/FlowPastSphere3D.py

@@ -49,7 +49,7 @@ def create_bcs(V, Q, mesh, **NS_namespace):
     outlet = "x[0] > 6-1e-8 && on_boundary"
 
     bmesh = BoundaryMesh(mesh, 'exterior')
-    cc = CellFunction('size_t', bmesh, 0)
+    cc = MeshFunction('size_t', bmesh, bmesh.topology().dim(), 0)
     ii = AutoSubDomain(lambda x, on_bnd: near(x[0], -3))
     ii.mark(cc, 1)
     smesh = SubMesh(bmesh, cc, 1)

oasis/problems/NSfracStep/LaminarChannel.py

@@ -6,7 +6,7 @@
 
 from ..NSfracStep import *
 from numpy import pi, arctan, array
-set_log_active(False)
+#set_log_active(False)
 
 
 # Override some problem specific parameters

oasis/problems/NSfracStep/SkewedFlow.py

@@ -8,7 +8,7 @@
 from ..SkewedFlow import *
 from numpy import cos, pi, cosh
 
-warning("""
+print("""
 This problem does not work well with IPCS since the outflow
 boundary condition
 
@@ -37,7 +37,7 @@ def problem_parameters(NS_parameters, **NS_namespace):
 def create_bcs(V, Q, mesh, **NS_namespace):
     # Create inlet profile by solving Poisson equation on boundary
     bmesh = BoundaryMesh(mesh, 'exterior')
-    cc = CellFunction('size_t', bmesh, 0)
+    cc = MeshFunction('size_t', bmesh, bmesh.topology().dim(), 0)
     ii = AutoSubDomain(inlet)
     ii.mark(cc, 1)
     smesh = SubMesh(bmesh, cc, 1)
@@ -49,7 +49,7 @@ def create_bcs(V, Q, mesh, **NS_namespace):
           bcs=[DirichletBC(Vu, Constant(0), DomainBoundary())])
 
     # Wrap the boundary function in an Expression to avoid the need to interpolate it back to V
-    class MyExp(Expression):
+    class MyExp(UserExpression):
         def eval(self, values, x):
             try:
                 values[0] = su(x)

oasis/solvers/NSfracStep/IPCS_ABCN.py

@@ -138,12 +138,10 @@ def get_solvers(use_krylov_solvers, krylov_solvers, bcs,
     else:
         ## tentative velocity solver ##
         u_sol = LUSolver()
-        u_sol.parameters['same_nonzero_pattern'] = True
+        #u_sol.parameters['same_nonzero_pattern'] = True
         ## pressure solver ##
         p_sol = LUSolver()
-        p_sol.parameters['reuse_factorization'] = True
-        if bcs['p'] == []:
-            p_sol.normalize = True
+        #p_sol.parameters['reuse_factorization'] = True
         sols = [u_sol, p_sol]
         ## scalar solver ##
         if len(scalar_components) > 0:
@@ -260,7 +258,7 @@ def pressure_solve(dp_, x_, Ap, b, p_sol, bcs, **NS_namespace):
     p_sol.solve(Ap, x_['p'], b['p'])
     t1.stop()
     # LUSolver use normalize directly for normalization of pressure
-    if hasattr(p_sol, 'normalize'):
+    if bcs['p'] == []:
         normalize(x_['p'])
 
     dpv = dp_.vector()

oasis/solvers/NSfracStep/IPCS_ABE.py

@@ -138,10 +138,10 @@ def get_solvers(use_krylov_solvers, krylov_solvers, bcs,
     else:
         ## tentative velocity solver ##
         u_sol = LUSolver('mumps')
-        u_sol.parameters['same_nonzero_pattern'] = True
+        #u_sol.parameters['same_nonzero_pattern'] = True
         ## pressure solver ##
         p_sol = LUSolver('mumps')
-        p_sol.parameters['reuse_factorization'] = True
+        #p_sol.parameters['reuse_factorization'] = True
         if bcs['p'] == []:
             p_sol.normalize = True
         sols = [u_sol, p_sol]

oasis/solvers/NSfracStep/LES/DynamicLagrangian.py

@@ -5,7 +5,7 @@
 
 from dolfin import (Function, FunctionSpace, TestFunction, sym, grad, dx, inner,
     sqrt, TrialFunction, project, CellVolume, as_vector, solve, Constant,
-    LagrangeInterpolator, assemble, FacetFunction, DirichletBC)
+    LagrangeInterpolator, assemble, MeshFunction, DirichletBC)
 from .DynamicModules import (tophatfilter, lagrange_average, compute_Lij,
                             compute_Mij)
 import numpy as np
@@ -36,7 +36,7 @@ def les_setup(u_, mesh, assemble_matrix, CG1Function, nut_krylov_solver, bcs, **
     Cs = Function(CG1)
     nut_form = Cs**2 * delta**2 * magS
     # Create nut_ BCs
-    ff = FacetFunction("size_t", mesh, 0)
+    ff = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
     bcs_nut = []
     for i, bc in enumerate(bcs['u0']):
         bc.apply(u_[0].vector())  # Need to initialize bc

oasis/solvers/NSfracStep/LES/ScaleDepDynamicLagrangian.py

@@ -4,7 +4,7 @@
 __license__ = 'GNU Lesser GPL version 3 or any later version'
 
 from dolfin import (Function, assemble, TestFunction, dx, solve, Constant,
-    FacetFunction, DirichletBC)
+    MeshFunction, DirichletBC)
 from .DynamicModules import (tophatfilter, lagrange_average, compute_Lij,
     compute_Mij, compute_Qij, compute_Nij)
 from . import DynamicLagrangian

oasis/solvers/NSfracStep/LES/Smagorinsky.py

@@ -4,7 +4,7 @@
 __license__ = 'GNU Lesser GPL version 3 or any later version'
 
 from dolfin import (Function, FunctionSpace, assemble, TestFunction, sym, grad, dx, inner, sqrt,
-    FacetFunction, DirichletBC, Constant)
+    MeshFunction, DirichletBC, Constant)
 
 from .common import derived_bcs
 

oasis/solvers/NSfracStep/LES/Wale.py

@@ -4,7 +4,7 @@
 __license__ = 'GNU Lesser GPL version 3 or any later version'
 
 from dolfin import (Function, FunctionSpace, assemble, TestFunction, sym, grad, tr,
-    Identity, dx, inner, Max, FacetFunction, DirichletBC, Constant)
+    Identity, dx, inner, Max, MeshFunction, DirichletBC, Constant)
 
 from .common import derived_bcs
 
@@ -29,7 +29,7 @@ def les_setup(u_, mesh, Wale, bcs, CG1Function, nut_krylov_solver, **NS_namespac
     Sd = sym(Gij * Gij) - 1. / 3. * Identity(mesh.geometry().dim()) * Skk * Skk
     nut_form = (Wale['Cw']**2 * pow(delta, 2. / dim) * pow(inner(Sd, Sd), 1.5)
                 / (Max(pow(inner(Sij, Sij), 2.5) + pow(inner(Sd, Sd), 1.25), 1e-6)))
-    ff = FacetFunction("size_t", mesh, 0)
+    ff = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
     bcs_nut = derived_bcs(CG1, bcs['u0'], u_)
     nut_ = CG1Function(nut_form, mesh, method=nut_krylov_solver,
                        bcs=bcs_nut, name='nut', bounded=True)

oasis/solvers/NSfracStep/LES/common.py

@@ -4,7 +4,7 @@
 __license__ = 'GNU Lesser GPL version 3 or any later version'
 
 import warnings
-from dolfin import FacetFunction, DirichletBC, Constant
+from dolfin import MeshFunction, DirichletBC, Constant
 
 
 def derived_bcs(V, original_bcs, u_):
@@ -15,7 +15,7 @@ def derived_bcs(V, original_bcs, u_):
     subdomain = original_bcs[0].user_sub_domain()
     if subdomain is None:
         mesh = V.mesh()
-        ff = FacetFunction("size_t", mesh, 0)
+        ff = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
         for i, bc in enumerate(original_bcs):
             bc.apply(u_[0].vector())  # Need to initialize bc
             m = bc.markers()  # Get facet indices of boundary