def compute_err(self, is_tent, velocity, t):
     if self.doErrControl:
         er_list_L2 = self.listDict['u2L2' if is_tent else 'u_L2']['list']
         er_list_H1 = self.listDict['u2H1' if is_tent else 'u_H1']['list']
         self.tc.start('errorV')
         # assemble is faster than errornorm
         errorL2_sq = assemble(inner(velocity - self.solution, velocity - self.solution) * dx)
         errorH1seminorm_sq = assemble(inner(grad(velocity - self.solution), grad(velocity - self.solution)) * dx)
         info('  H1 seminorm error: %f' % sqrt(errorH1seminorm_sq))
         errorL2 = sqrt(errorL2_sq)
         errorH1 = sqrt(errorL2_sq + errorH1seminorm_sq)
         info("  Relative L2 error in velocity = %f" % (errorL2 / self.analytic_v_norm_L2))
         self.last_error = errorH1 / self.analytic_v_norm_H1
         self.last_status_functional = self.last_error
         info("  Relative H1 error in velocity = %f" % self.last_error)
         er_list_L2.append(errorL2)
         er_list_H1.append(errorH1)
         self.tc.end('errorV')
         if self.testErrControl:
             er_list_test_H1 = self.listDict['u2H1test' if is_tent else 'u_H1test']['list']
             er_list_test_L2 = self.listDict['u2L2test' if is_tent else 'u_L2test']['list']
             self.tc.start('errorVtest')
             er_list_test_L2.append(errornorm(velocity, self.solution, norm_type='L2', degree_rise=0))
             er_list_test_H1.append(errornorm(velocity, self.solution, norm_type='H1', degree_rise=0))
             self.tc.end('errorVtest')
         # stopping criteria for detecting diverging solution
         if self.last_error > self.divergence_treshold:
             raise RuntimeError('STOPPED: Failed divergence test!')

  def __init__(self, coarse_mesh, nref, p_coarse, p_fine):
    super(LaplaceEigenvalueProblem, self).__init__(coarse_mesh, nref, p_coarse, p_fine)

    print0("Assembling fine-mesh problem")

    self.dirichlet_bdry = lambda x,on_boundary: on_boundary

    bc = DirichletBC(self.V_fine, 0.0, self.dirichlet_bdry)
    u = TrialFunction(self.V_fine)
    v = TestFunction(self.V_fine)
    a = inner(grad(u), grad(v))*dx
    m = u*v*dx

    # Assemble the stiffness matrix and the mass matrix.
    b = v*dx # just need this to feed an argument to assemble_system
    assemble_system(a, b, bc, A_tensor=self.A_fine)
    assemble_system(m, b, bc, A_tensor=self.B_fine)
    # set the diagonal elements of M corresponding to boundary nodes to zero to
    # remove spurious eigenvalues.
    bc.zero(self.B_fine)

    print0("Assembling coarse-mesh problem")

    self.bc_coarse = DirichletBC(self.V_coarse, 0.0, self.dirichlet_bdry)
    u = TrialFunction(self.V_coarse)
    v = TestFunction(self.V_coarse)
    a = inner(grad(u), grad(v))*dx
    m = u*v*dx

    # Assemble the stiffness matrix and the mass matrix, without Dirichlet BCs. Dirichlet DOFs will be removed later.
    assemble(a, tensor=self.A_coarse)
    assemble(m, tensor=self.B_coarse)

    def update_time(self, actual_time, step_number):
        super(Problem, self).update_time(actual_time, step_number)
        if self.actual_time > 0.5 and int(round(self.actual_time * 1000)) % 1000 == 0:
            self.isWholeSecond = True
            seconds = int(round(self.actual_time))
            self.second_list.append(seconds)
            self.N1 = seconds*self.stepsInSecond
            self.N0 = (seconds-1)*self.stepsInSecond
        else:
            self.isWholeSecond = False

        self.solution = self.assemble_solution(self.actual_time)

        # Update boundary condition
        self.tc.start('updateBC')
        self.v_in.assign(self.solution)
        self.tc.end('updateBC')

        # construct analytic pressure (used for computing pressure and force errors)
        self.tc.start('analyticP')
        analytic_pressure = womersleyBC.analytic_pressure(self.factor, self.actual_time)
        self.sol_p = interpolate(analytic_pressure, self.pSpace)
        self.tc.end('analyticP')

        self.tc.start('analyticVnorms')
        self.analytic_v_norm_L2 = norm(self.solution, norm_type='L2')
        self.analytic_v_norm_H1 = norm(self.solution, norm_type='H1')
        self.analytic_v_norm_H1w = sqrt(assemble((inner(grad(self.solution), grad(self.solution)) +
                                                  inner(self.solution, self.solution)) * self.dsWall))
        self.listDict['av_norm_L2']['list'].append(self.analytic_v_norm_L2)
        self.listDict['av_norm_H1']['list'].append(self.analytic_v_norm_H1)
        self.listDict['av_norm_H1w']['list'].append(self.analytic_v_norm_H1w)
        self.tc.end('analyticVnorms')

    def __init__(self, v, v_out, solver_parameters=None):

        if isinstance(v, expression.Expression) or \
           not isinstance(v, (ufl.core.expr.Expr, function.Function)):
            raise ValueError("Can only project UFL expression or Functions not '%s'" % type(v))

        self._same_fspace = (isinstance(v, function.Function) and v.function_space() ==
                             v_out.function_space())
        self.v = v
        self.v_out = v_out

        if not self._same_fspace:
            V = v_out.function_space()

            p = ufl_expr.TestFunction(V)
            q = ufl_expr.TrialFunction(V)

            a = ufl.inner(p, q)*ufl.dx
            L = ufl.inner(p, v)*ufl.dx

            problem = vs.LinearVariationalProblem(a, L, v_out)

            if solver_parameters is None:
                solver_parameters = {}

            solver_parameters.setdefault("ksp_type", "cg")

            self.solver = vs.LinearVariationalSolver(problem,
                                                     solver_parameters=solver_parameters)

 def compute_err(self, is_tent, velocity, t):
     super(Problem, self).compute_err(is_tent, velocity, t)
     er_list_H1w = self.listDict['u2H1w' if is_tent else 'u_H1w']['list']
     errorH1wall = sqrt(assemble((inner(grad(velocity - self.solution), grad(velocity - self.solution)) +
                                  inner(velocity - self.solution, velocity - self.solution)) * self.dsWall))
     er_list_H1w.append(errorH1wall)
     print('  Relative H1wall error:', errorH1wall / self.analytic_v_norm_H1w)

    def compute_functionals(self, velocity, pressure, t):
        if self.args.wss:
            info('Computing stress tensor')
            I = Identity(velocity.geometric_dimension())
            T = TensorFunctionSpace(self.mesh, 'Lagrange', 1)
            stress = project(-pressure*I + 2*sym(grad(velocity)), T)
            info('Generating boundary mesh')
            wall_mesh = BoundaryMesh(self.mesh, 'exterior')
            # wall_mesh = SubMesh(self.mesh, self.facet_function, 1)   # QQ why does not work?
            # plot(wall_mesh, interactive=True)
            info('  Boundary mesh geometric dim: %d' % wall_mesh.geometry().dim())
            info('  Boundary mesh topologic dim: %d' % wall_mesh.topology().dim())
            info('Projecting stress to boundary mesh')
            Tb = TensorFunctionSpace(wall_mesh, 'Lagrange', 1)
            stress_b = interpolate(stress, Tb)
            self.fileDict['wss']['file'] << stress_b


            if False:  # does not work
                info('Computing WSS')
                n = FacetNormal(wall_mesh)
                info(stress_b, True)
                # wss = stress_b*n - inner(stress_b*n, n)*n
                wss = dot(stress_b, n) - inner(dot(stress_b, n), n)*n   # equivalent
                Vb = VectorFunctionSpace(wall_mesh, 'Lagrange', 1)
                Sb = FunctionSpace(wall_mesh, 'Lagrange', 1)
                # wss_func = project(wss, Vb)
                wss_norm = project(sqrt(inner(wss, wss)), Sb)
                plot(wss_norm, interactive=True)

    def initialize(self, V, Q, PS, D):
        super(Problem, self).initialize(V, Q, PS, D)

        print("IC type: " + self.ic)
        print("Velocity scale factor = %4.2f" % self.factor)
        reynolds = 728.761 * self.factor
        print("Computing with Re = %f" % reynolds)

        # set constants for
        self.area = assemble(interpolate(Expression("1.0"), Q) * self.dsIn)  # inflow area

        self.solution = interpolate(Expression(("0.0", "0.0", "factor*(1081.48-43.2592*(x[0]*x[0]+x[1]*x[1]))"),
                                               factor=self.factor), self.vSpace)
        analytic_pressure = womersleyBC.average_analytic_pressure_expr(self.factor)
        self.sol_p = interpolate(analytic_pressure, self.pSpace)
        self.analytic_gradient = womersleyBC.average_analytic_pressure_grad(self.factor)
        self.analytic_pressure_norm = norm(self.sol_p, norm_type='L2')
        self.analytic_v_norm_L2 = norm(self.solution, norm_type='L2')
        self.analytic_v_norm_H1 = norm(self.solution, norm_type='H1')
        self.analytic_v_norm_H1w = sqrt(assemble((inner(grad(self.solution), grad(self.solution)) +
                                                  inner(self.solution, self.solution)) * self.dsWall))
        print("Prepared analytic solution.")

        self.pg_normalization_factor.append(womersleyBC.average_analytic_pressure_grad(self.factor))
        self.p_normalization_factor.append(self.analytic_pressure_norm)
        self.vel_normalization_factor.append(norm(self.solution, norm_type='L2'))

        print('Normalisation factors (vel, p, pg):', self.vel_normalization_factor[0], self.p_normalization_factor[0],
              self.pg_normalization_factor[0])

        one = (interpolate(Expression('1.0'), Q))
        self.outflow_area = assemble(one*self.dsOut)
        print('Outflow area:', self.outflow_area)

 def compute_force(self, velocity, pressure, t):
     self.tc.start('errorForce')
     I = Identity(3)  # Identity tensor
     def T(p, v):
         return -p * I + 2.0 * self.nu * sym(grad(v))
     error_force = sqrt(
             assemble(inner((T(pressure, velocity) - T(self.sol_p, self.solution)) * self.normal,
                            (T(pressure, velocity) - T(self.sol_p, self.solution)) * self.normal) * self.dsWall))
     an_force = sqrt(assemble(inner(T(self.sol_p, self.solution) * self.normal,
                                         T(self.sol_p, self.solution) * self.normal) * self.dsWall))
     an_f_normal = sqrt(assemble(inner(inner(T(self.sol_p, self.solution) * self.normal, self.normal),
                                            inner(T(self.sol_p, self.solution) * self.normal, self.normal)) * self.dsWall))
     error_f_normal = sqrt(
             assemble(inner(inner((T(self.sol_p, self.solution) - T(pressure, velocity)) * self.normal, self.normal),
                            inner((T(self.sol_p, self.solution) - T(pressure, velocity)) * self.normal, self.normal)) * self.dsWall))
     an_f_shear = sqrt(
             assemble(inner((I - outer(self.normal, self.normal)) * T(self.sol_p, self.solution) * self.normal,
                            (I - outer(self.normal, self.normal)) * T(self.sol_p, self.solution) * self.normal) * self.dsWall))
     error_f_shear = sqrt(
             assemble(inner((I - outer(self.normal, self.normal)) *
                            (T(self.sol_p, self.solution) - T(pressure, velocity)) * self.normal,
                            (I - outer(self.normal, self.normal)) *
                            (T(self.sol_p, self.solution) - T(pressure, velocity)) * self.normal) * self.dsWall))
     self.listDict['a_force_wall']['list'].append(an_force)
     self.listDict['a_force_wall_normal']['list'].append(an_f_normal)
     self.listDict['a_force_wall_shear']['list'].append(an_f_shear)
     self.listDict['force_wall']['list'].append(error_force)
     self.listDict['force_wall_normal']['list'].append(error_f_normal)
     self.listDict['force_wall_shear']['list'].append(error_f_shear)
     if self.isWholeSecond:
         self.listDict['force_wall']['slist'].append(
             sqrt(sum([i*i for i in self.listDict['force_wall']['list'][self.N0:self.N1]])/self.stepsInSecond))
     print('  Relative force error:', error_force/an_force)
     self.tc.end('errorForce')

    def update_time(self, actual_time, step_number):
        super(Problem, self).update_time(actual_time, step_number)
        if self.actual_time > 0.5 and abs(math.modf(actual_time)[0]) < 0.5*self.metadata['dt']:
            self.second_list.append(int(round(self.actual_time)))

        self.solution = self.assemble_solution(self.actual_time)

        # Update boundary condition
        self.tc.start('updateBC')
        self.v_in.assign(self.onset_factor * self.solution)
        self.tc.end('updateBC')

        # construct analytic pressure (used for computing pressure and force errors)
        self.tc.start('analyticP')
        analytic_pressure = womersleyBC.analytic_pressure(self.factor, self.actual_time)
        self.sol_p = interpolate(analytic_pressure, self.pSpace)
        self.tc.end('analyticP')

        self.tc.start('analyticVnorms')
        self.analytic_v_norm_L2 = norm(self.solution, norm_type='L2')
        self.analytic_v_norm_H1 = norm(self.solution, norm_type='H1')
        self.analytic_v_norm_H1w = sqrt(assemble((inner(grad(self.solution), grad(self.solution)) +
                                                  inner(self.solution, self.solution)) * self.dsWall))
        self.listDict['av_norm_L2']['list'].append(self.analytic_v_norm_L2)
        self.listDict['av_norm_H1']['list'].append(self.analytic_v_norm_H1)
        self.listDict['av_norm_H1w']['list'].append(self.analytic_v_norm_H1w)
        self.tc.end('analyticVnorms')

 def nonlinearity(function):
     if self.use_ema:
        return 2*inner(dot(sym(grad(function)), u_ext), v1) * dx + inner(div(function)*u_ext, v1) * dx
         # return 2*inner(dot(sym(grad(function)), u_ext), v) * dx + inner(div(u_ext)*function, v) * dx
         # QQ implement this way?
     else:
         return inner(dot(grad(function), u_ext), v1) * dx

    def update_time(self, actual_time, step_number):
        super(Problem, self).update_time(actual_time, step_number)
        if self.actual_time > 0.5 and int(round(self.actual_time * 1000)) % 1000 == 0:
            self.isWholeSecond = True
            seconds = int(round(self.actual_time))
            self.second_list.append(seconds)
            self.N1 = seconds*self.stepsInSecond
            self.N0 = (seconds-1)*self.stepsInSecond
        else:
            self.isWholeSecond = False

        # Update boundary condition
        self.tc.start('updateBC')
        if not self.ic == 'correct':
            self.v_in.t = self.actual_time
        self.tc.end('updateBC')

        self.tc.start('analyticVnorms')
        self.analytic_v_norm_L2 = norm(self.solution, norm_type='L2')
        self.analytic_v_norm_H1 = norm(self.solution, norm_type='H1')
        self.analytic_v_norm_H1w = sqrt(assemble((inner(grad(self.solution), grad(self.solution)) +
                                                  inner(self.solution, self.solution)) * self.dsWall))
        self.listDict['av_norm_L2']['list'].append(self.analytic_v_norm_L2)
        self.listDict['av_norm_H1']['list'].append(self.analytic_v_norm_H1)
        self.listDict['av_norm_H1w']['list'].append(self.analytic_v_norm_H1w)
        self.tc.end('analyticVnorms')

 def compute_err(self, is_tent, velocity, t):
     super(Problem, self).compute_err(is_tent, velocity, t)
     er_list_H1w = self.listDict['u2H1w' if is_tent else 'u_H1w']['list']
     errorH1wall = sqrt(assemble((inner(grad(velocity - self.solution), grad(velocity - self.solution)) +
                                  inner(velocity - self.solution, velocity - self.solution)) * self.dsWall))
     er_list_H1w.append(errorH1wall)
     print('  Relative H1wall error:', errorH1wall / self.analytic_v_norm_H1w)
     if self.isWholeSecond:
         self.listDict['u2H1w' if is_tent else 'u_H1w']['slist'].append(
             sqrt(sum([i*i for i in er_list_H1w[self.N0:self.N1]])/self.stepsInSecond))

 def diffusion(fce):
     if self.useLaplace:
         return nu * inner(grad(fce), grad(v1)) * dx
     else:
         form = inner(nu * 2 * sym(grad(fce)), sym(grad(v1))) * dx
         if self.bcv == 'CDN':
             return form
         if self.bcv == 'LAP':
             return form - inner(nu * dot(grad(fce).T, n), v1) * problem.get_outflow_measure_form()
         if self.bcv == 'DDN':
             return form  # additional term must be added to non-constant part

 def diffusion(fce):
     if self.useLaplace:
         return nu*inner(grad(fce), grad(v1)) * dx
     else:
         form = inner(nu * 2 * sym(grad(fce)), sym(grad(v1))) * dx
         if self.bcv == 'CDN':
             # IMP will work only if p=0 on output, or we must add term
             # inner(p0*n, v)*problem.get_outflow_measure_form() to avoid boundary layer
             return form
         if self.bcv == 'LAP':
             return form - inner(nu*dot(grad(fce).T, n), v1)  * problem.get_outflow_measure_form()
         if self.bcv == 'DDN':
             # IMP will work only if p=0 on output, or we must add term
             # inner(p0*n, v)*problem.get_outflow_measure_form() to avoid boundary layer
             return form  # additional term must be added to non-constant part

def form(V, itype, request):
    if request.param == "functional":
        u = ufl.Coefficient(V)
        v = ufl.Coefficient(V)
    elif request.param == "1-form":
        u = ufl.Coefficient(V)
        v = ufl.TestFunction(V)
    elif request.param == "2-form":
        u = ufl.TrialFunction(V)
        v = ufl.TestFunction(V)

    if itype == "cell":
        return ufl.inner(u, v)*ufl.dx
    elif itype == "ext_facet":
        return ufl.inner(u, v)*ufl.ds
    elif itype == "int_facet":
        return ufl.inner(u('+'), v('-'))*ufl.dS

def forms(arguments, coefficients):
    v, u = arguments
    c, f = coefficients
    n = FacetNormal(triangle)
    a = u * v * dx
    L = f * v * dx
    b = u * v * dx(0) + inner(c * grad(u), grad(v)) * \
        dx(1) + dot(n, grad(u)) * v * ds + f * v * dx
    return (a, L, b)

    def error_indicators(self):
        """
        Generate and return linear form defining error indicators
        """
        # Extract these to increase readability
        R_T = self._R_T
        R_dT = self._R_dT
        z = self._Ez_h
        z_h = self._z_h

        # Define linear form for computing error indicators
        v = self.module.TestFunction(self._DG0)
        eta_T = (v * inner(R_T, z - z_h) * dx(self.domain) +
                 avg(v)*(inner(R_dT('+'), (z - z_h)('+')) +
                         inner(R_dT('-'), (z - z_h)('-'))) * dS(self.domain) +
                 v * inner(R_dT, z - z_h) * ds(self.domain))

        return eta_T

    def compute_functionals(self, velocity, pressure, t, step):
        if self.args.wss == 'all' or \
                (step >= self.stepsInCycle and self.args.wss == 'peak' and
                 (self.distance_from_chosen_steps < 0.5 * self.metadata['dt'])):
            # TODO check if choosing time steps works properly
            # QQ might skip step? change 0.5 to 0.51?
            self.tc.start('WSS')
            begin('WSS (%dth step)' % step)
            if self.args.wss_method == 'expression':
                stress = project(self.nu*2*sym(grad(velocity)), self.T)
                # pressure is not used as it contributes only to the normal component
                stress.set_allow_extrapolation(True)   # need because of some inaccuracies in BoundaryMesh coordinates
                stress_b = interpolate(stress, self.Tb)    # restrict stress to boundary mesh
                # self.fileDict['stress']['file'].write(stress_b, self.actual_time)
                # info('Saved stress tensor')
                info('Computing WSS')
                wss = dot(stress_b, self.nb) - inner(dot(stress_b, self.nb), self.nb)*self.nb
                wss_func = project(wss, self.Vb)
                wss_norm = project(sqrt_ufl(inner(wss, wss)), self.Sb)
                info('Saving WSS')
                self.fileDict['wss']['file'].write(wss_func, self.actual_time)
                self.fileDict['wss_norm']['file'].write(wss_norm, self.actual_time)
            if self.args.wss_method == 'integral':
                wss_norm = Function(self.SDG)
                mS = TestFunction(self.SDG)
                scaling = 1/FacetArea(self.mesh)
                stress = self.nu*2*sym(grad(velocity))
                wss = dot(stress, self.normal) - inner(dot(stress, self.normal), self.normal)*self.normal
                wss_norm_form = scaling*mS*sqrt_ufl(inner(wss, wss))*ds   # ds is integral over exterior facets only
                assemble(wss_norm_form, tensor=wss_norm.vector())
                self.fileDict['wss_norm']['file'].write(wss_norm, self.actual_time)

                # to get vector WSS values:
                # NT this works, but in ParaView for (DG,1)-vector space glyphs are displayed in cell centers
                # wss_vector = []
                # for i in range(3):
                #     wss_component = Function(self.SDG)
                #     wss_vector_form = scaling*wss[i]*mS*ds
                #     assemble(wss_vector_form, tensor=wss_component.vector())
                #     wss_vector.append(wss_component)
                # wss_func = project(as_vector(wss_vector), self.VDG)
                # self.fileDict['wss']['file'].write(wss_func, self.actual_time)
            self.tc.end('WSS')
            end()

 def save_vel(self, is_tent, field, t):
     self.vFunction.assign(field)
     self.fileDict['u2' if is_tent else 'u']['file'] << self.vFunction
     if self.doSaveDiff:
         self.vFunction.assign((1.0 / self.vel_normalization_factor[0]) * (field - self.solution))
         self.fileDict['u2D' if is_tent else 'uD']['file'] << self.vFunction
     if self.args.ldsg:
         # info(div(2.*sym(grad(field))-grad(field)).ufl_shape)
         form = div(2.*sym(grad(field))-grad(field))
         self.pFunction.assign(project(sqrt_ufl(inner(form, form)), self.pSpace))
         self.fileDict['ldsg2' if is_tent else 'ldsg']['file'] << self.pFunction

def test_gradient_error(cell, degree):
    """Test that tabulating gradient evaluations of the trace
    element triggers `gem.Failure` to raise the TraceError
    exception.
    """
    trace_element = FiniteElement("HDiv Trace", cell, degree)
    lambdar = TrialFunction(trace_element)
    gammar = TestFunction(trace_element)

    with pytest.raises(TraceError):
        compile_form(inner(grad(lambdar('+')), grad(gammar('+'))) * dS)