How to find the value at the cell-centered

General
1

I have this code that im working on. I only know how to apply ‘c.vector().get_local()’ to get the solution. But i want to know if there is a way to find the value at the cell-centered


def solver(f, c0, t_0, t_N, dt, radius, height, xmin, ymin, fs_order, u_in, u_x, D_xx, D_xy, D_yy, solution_dict):
    diameter = radius * 2
    xmax = diameter
    ymax = height

    domain = Rectangle(Point(xmin, ymin), Point(xmax, ymax))
    mesh = generate_mesh(domain, 4)

    element = FiniteElement(family='CG', cell=mesh.ufl_cell(),
                            degree=fs_order, quad_scheme='quadrature')
    V = FunctionSpace(mesh, element)

    c = TrialFunction(V)
    w = TestFunction(V)

    u_y = u_in
    D_yx = D_xy
    u = as_vector([u_x, u_y])
    D = as_matrix(((D_xx, D_xy), (D_yx, D_yy)))

    c_in = c0 / 100  # [ppm]
    g_in = u_in * c_in
    c_0 = Constant(c0)

    if t_0 in solution_dict:
        c_init = solution_dict[t_0]
    else:
        c_init = Function(V)
        c_init.assign(c_0)

    domains = MeshFunction("size_t", mesh, mesh.topology().dim())
    boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

    class BottomBoundary(SubDomain):
        def inside(self, x, on_boundary):
            tol = 1E-14
            return on_boundary and abs(x[1] - ymin) < tol

    Gamma_1 = BottomBoundary()
    Gamma_1.mark(boundary_parts, 1)
    c_B = Constant(c_in)
    bc1 = DirichletBC(V, c_B, boundary_parts, 1)

    class RightBoundary(SubDomain):
        def inside(self, x, on_boundary):
            tol = 1E-14
            return on_boundary and abs(x[0] - xmax) < tol

    Gamma_2 = RightBoundary()
    Gamma_2.mark(boundary_parts, 2)
    c_R = Constant(0.0)
    bc2 = DirichletBC(V, c_R, boundary_parts, 2)

    class TopBoundary(SubDomain):
        def inside(self, x, on_boundary):
            tol = 1E-14
            return on_boundary and abs(x[1] - ymax) < tol

    Gamma_3 = TopBoundary()
    Gamma_3.mark(boundary_parts, 3)
    c_T = Constant(0.0)
    bc3 = DirichletBC(V, c_T, boundary_parts, 3)

    class LeftBoundary(SubDomain):
        def inside(self, x, on_boundary):
            tol = 1E-14
            return on_boundary and abs(x[0] - xmin) < tol

    Gamma_4 = LeftBoundary()
    Gamma_4.mark(boundary_parts, 4)
    c_L = Constant(0.0)
    bc4 = DirichletBC(V, c_L, boundary_parts, 4)

    bcs = [bc1]

    dx = Measure('dx', domain=mesh, subdomain_data=domains)
    ds = Measure('ds', domain=mesh, subdomain_data=boundary_parts)

    a = c * w * dx + dt * c * w * dx + dt * inner(grad(w), dot(D, grad(c)) - u * c) * dx + dt * c * w * u_y * ds(3)
    L = c_init * w * dx + dt * f * w * dx - dt * w * g_in * ds(1)

    c = Function(V)
    c.rename('concentración', 'c')
    c_file = File('Transporte_Prueba_new/concentracion_new.pvd')

    c.assign(c_init)

    A = assemble(a)
    b = assemble(L)

    for bc in bcs:
        bc.apply(A, b)

    solve(A, c.vector(), b, 'lu')

    print(f"Time step: {t_0}, Concentration: {c.vector().get_local()[0]}")

    solution_dict[t_0 + dt] = c.copy(deepcopy=True)

    c_file << (c, t_0)

    return c

def run_solver():
    radius = 0.01  # m
    height = 0.05  # m
    xmin = 0.
    ymin = 0.
    fs_order = 1
    u_in = 2.5e-07
    u_x = 0
    D_xx = 0
    D_xy = 0
    D_yy = 5e-10
    f = 0
    c0 = 1000
    t_0 = 0.
    t_N = 90
    dt = 10

    solution_dict = {}

    while True:
        c = solver(f, c0, t_0, t_N, dt, radius, height, xmin, ymin, fs_order, u_in, u_x, D_xx, D_xy, D_yy, solution_dict)
        t_0 += dt

        user_input = input("Enter 'q' to quit or press Enter to continue with the next run: ")
        if user_input.lower() == 'q':
            break

if __name__ == '__main__':
    run_solver()
2

The easiest way to get the value at the cell center, you could interpolate your solution into a DG 0 space, and then get the values throguh vector().get_local(), and then the midpoints with dolfin.FunctionSpace.tabulate_dof_coordinates()