Python ufl.inner函数代码示例

 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