How to estimate the gradient along the vector between the two centers of adjacent cells?

General
1

Hello, I am following this tutorial deal.II step-9 to estimate the errors of solution and refine the mesh. In this tutorial, the solution error is estimated by |\nabla_h u(K)|. It is the sum of gradients along the vectors which connect the centers of adjacent cells:

\underbrace{ \left(\sum_{K'} \frac{{\mathbf y}_{KK'} {\mathbf y}_{KK'}^T} {|{\mathbf y}_{KK'}|^2}\right)}_{=:Y} \nabla u \approx \sum_{K'} \frac{{\mathbf y}_{KK'}}{|{\mathbf y}_{KK'}|} \frac{u(K') - u(K)}{|{\mathbf y}_{KK'}|}.
\nabla u \approx Y^{-1} \left( \sum_{K'} \frac{{\mathbf y}_{KK'}}{|{\mathbf y}_{KK'}|} \frac{u(K') - u(K)}{|{\mathbf y}_{KK'}|} \right).

{\mathbf y}_{KK'} is the vector between two adjacent cell centers.

My question is:

  1. How to get the center of the cell, is there any convenient function to do this ?
  2. How to know which cells are adjacent to a given cell ?

Is it possible to do 1,2 efficiently by UFL ?
Thanks.

2

Here is an example of how to get the midpoints, and the values of u at the midpoints:

from mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.functionspace(mesh, ("CG", 2))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0] + x[1] ** 2)

Q = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.geometry.dim,)))
q = dolfinx.fem.Function(Q)
q.interpolate(lambda x: x[: mesh.geometry.dim])


num_cells_local = (
    mesh.topology.index_map(mesh.topology.dim).size_local
    + mesh.topology.index_map(mesh.topology.dim).num_ghosts
)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
exterior_facet_indices = dolfinx.mesh.exterior_facet_indices(mesh.topology)
num_facets_local = mesh.topology.index_map(mesh.topology.dim - 1).size_local
interior_facet_indices = np.arange(num_facets_local)
interior_facet_indices = np.delete(interior_facet_indices, exterior_facet_indices)

import basix

ref_cell_geometry = basix.cell.geometry(mesh.basix_cell())
midpoint = np.sum(ref_cell_geometry, axis=0) / ref_cell_geometry.shape[0]

u_at_midpoint = dolfinx.fem.Expression(u, midpoint.reshape(-1, len(midpoint)))
midpoint_values = u_at_midpoint.eval(mesh, np.arange(num_cells_local, dtype=np.int32))


f_to_c = mesh.topology.connectivity(mesh.topology.dim - 1, mesh.topology.dim)
q_r = q.x.array.reshape(-1, mesh.geometry.dim)
for interior_facet in interior_facet_indices:
    cells = f_to_c.links(interior_facet)

    u_values = midpoint_values[cells]
    q_values = q_r[cells]

    assert np.isclose(u_values[0][0], q_values[0][0] + q_values[0][1] ** 2)
    assert np.isclose(u_values[1][0], q_values[1][0] + q_values[1][1] ** 2)

This can now easily be extended to loop through each cell, then finding the facets of that cell, finding the neighbour, getting the midpoints to the adjacent cells and evaluate u in one cell minus u in another.

2 Likes
3

Hello, dokken, thanks for your help. I have tried to implement the code but there is a problem. I expect the facets always connect to two cells, but after the refinement, if I run the code with more than one process, some facets will connect to only one cell. In that case I cannot get the vector between the two cell centers. In the code below, this line will cause error

edge_index = np.vstack(edge_index)

Because some items of edge_index only has one cell index, while others have two cell indexes.

#!/mnt/d/software_install/fenics-0.9/bin/python3
from mpi4py import MPI
import numpy as np
from dolfinx.fem.petsc import LinearProblem
import dolfinx, ufl, basix
import pyvista
from dolfinx import plot
from dolfinx.cpp.mesh import DiagonalType

N = 8
domain = dolfinx.mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=[(-1.0, -1.0), (1.0, 1.0)],
    n = [N, N])

N_refine = 9
for i_iter in range(N_refine):
    print(f"++++++++{i_iter}++++++++")
    V = dolfinx.fem.functionspace(domain, ("Lagrange", 1))
    x = ufl.SpatialCoordinate(domain)
    
    beta = ufl.as_vector([2.0, 
                          1 + 4.0/5.0*ufl.sin(8.0*ufl.pi*x[0])])
    n = ufl.FacetNormal(domain)
    beta_dot_n = ufl.dot(beta, n)
    beta_dot_n_inflow = ufl.conditional(ufl.lt(beta_dot_n, 0), beta_dot_n, 0.0)

    x0 = [-3.0/4, -3.0/4]
    delta = 0.1 * ufl.Circumradius(domain)
    s = 0.1
    f = ufl.conditional(
        ufl.lt((x[0]-x0[0])**2 + (x[1]-x0[1])**2, s*s), 
        1.0/(10*s*s), 
        0)

    norm_x2 = x[0]**2 + x[1]**2
    g = ufl.exp(5*(1-norm_x2)) * ufl.sin(16.0*ufl.pi*norm_x2)

    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    a = ufl.dot(beta, ufl.grad(u)) * (v + delta*ufl.dot(beta, ufl.grad(v))) * ufl.dx - u*beta_dot_n_inflow*v*ufl.ds

    L = f*(v + delta*ufl.dot(beta, ufl.grad(v)))*ufl.dx - g*beta_dot_n_inflow*v*ufl.ds

    problem = LinearProblem(a, L, bcs=[], 
                            petsc_options={"ksp_type": "gmres", 
                                "pc_type":  "lu", 
                                "ksp_converged_reason":"",
                                "ksp_monitor":""})
    uh = problem.solve()

    filename = f"solution_iter_{i_iter}.xdmf"
    with dolfinx.io.XDMFFile(domain.comm, filename, "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(uh)
    
    if i_iter == 8:
        tdim = domain.topology.dim
        domain.topology.create_connectivity(tdim, tdim)
        topology, cell_types, geometry = plot.vtk_mesh(domain, tdim)
        grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
        plotter = pyvista.Plotter()
        plotter.add_mesh(grid, show_edges=True)
        plotter.view_xy()
        plotter.show()

    c_dim = domain.geometry.dim
    ###### Get the center of each cell ######
    Q = dolfinx.fem.functionspace(domain, ("DG", 0, (c_dim,)))
    q = dolfinx.fem.Function(Q)
    q.interpolate(lambda x: x[:c_dim])
    # q_r is the coordinate of the cell center
    q_r = q.x.array.reshape(-1, c_dim)

    ###### Get the index of all edges that are not on the boundary ###### 
    num_cells_local = domain.topology.index_map(c_dim).size_local \
                + domain.topology.index_map(c_dim).num_ghosts

    domain.topology.create_connectivity(c_dim - 1, c_dim)
    domain.topology.create_connectivity(c_dim, c_dim - 1)
    exterior_facet_indices = dolfinx.mesh.exterior_facet_indices(domain.topology)
    num_facets_local = domain.topology.index_map(domain.topology.dim - 1).size_local
    interior_facet_indices = np.arange(num_facets_local)
    interior_facet_indices = np.delete(interior_facet_indices, exterior_facet_indices)

    ###### Obtain the connection relationship between cells and facets ###### 
    f_to_c = domain.topology.connectivity(c_dim - 1, c_dim)

    ###### Retrieve all edges connecting the centers of adjacent cells, where the index is actually the cell ###### 
    edge_index = [f_to_c.links(interior_facet) for interior_facet in interior_facet_indices]
    #print(edge_index)
    edge_index = np.vstack(edge_index)
    edge_index = np.vstack((edge_index, edge_index[:, ::-1])).T

    y_kk = q_r[edge_index[0]] - q_r[edge_index[1]]
    L_y_kk = np.linalg.norm(y_kk, axis=1)[:, None]
    y_kk = y_kk / L_y_kk #vector Y_{KK'}, normalized
    
    # get the vetex coordinates of reference cell
    ref_cell_geometry = basix.cell.geometry(domain.basix_cell())
    # get the mid points
    midpoint = np.sum(ref_cell_geometry, axis=0) / ref_cell_geometry.shape[0]

    uh_at_midpoint = dolfinx.fem.Expression(uh, midpoint.reshape(-1, len(midpoint)))
    # (num_cells_local, 1)
    midpoint_values = uh_at_midpoint.eval(domain, np.arange(num_cells_local, dtype=np.int32))
    
    h_local = dolfinx.fem.Expression(ufl.Circumradius(domain), midpoint.reshape(-1, len(midpoint)))
    # (num_cells_local, 1)
    h_values = h_local.eval(domain, np.arange(num_cells_local, dtype=np.int32))

    d_uh = midpoint_values[edge_index[0]] - midpoint_values[edge_index[1]]
    grad_uh = y_kk * d_uh / L_y_kk

    grad_uh_cell = []
    for i in range(num_cells_local):
        index = np.where(edge_index[0] == i)[0]
        g = grad_uh[index].mean(axis=0)
        y_kk_sub = y_kk[index]
        Y = np.sum(np.einsum('bi, bj->bij', y_kk_sub, y_kk_sub), axis=0)
        if np.linalg.det(Y) < 1.0e-6:
            grad_uh_cell.append(100.0)
        else:
            grad_cell = np.matmul(np.linalg.inv(Y), g) 
            grad_cell = np.linalg.norm(grad_cell) * (h_values[i]**2)
            grad_uh_cell.append(grad_cell[0])
    
    grad_uh_cell = np.array(grad_uh_cell)

    sort_index = np.argsort(grad_uh_cell)

    cell_index = sort_index[-int(0.333*num_cells_local):]

    edge_index = dolfinx.mesh.compute_incident_entities(domain.topology, cell_index.astype(np.int32), 2, 1)

    new_domain, _, _ = dolfinx.mesh.refine(domain, edge_index)
    domain = new_domain
4

I solve this problem by adding a partitioner argument to refine function

partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode)

new_domain, _, _ = dolfinx.mesh.refine(domain, edge_index, partitioner=partitioner)
1 Like