I get a wrong result when using assemble_vector

dolfinx
,
1

Hi everyone! I want to get the integral of a given variable over all internal facets. When I use assemble_vector the size of the array is different from the number of internal boundaries and I get the wrong integral value. The following is MWE

import dolfinx
from dolfinx import *

from dolfinx.fem import (Constant,  Function, functionspace, assemble_scalar, 
                         assemble_vector,
                         dirichletbc, form, locate_dofs_topological)

from petsc4py import PETSc
from dolfinx.mesh import create_unit_square, locate_entities, meshtags

from joblib import Parallel, delayed
import ufl
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, dS, grad, inner, lhs, rhs,jump, avg)

import numpy as np
import time
from time import sleep
from mpi4py import MPI
N = 2
L = 100
# Generating mesh
msh = mesh.create_rectangle(comm = MPI.COMM_WORLD, points = ((0, 0), (L, L)), n = (N, N),
                            cell_type = mesh.CellType.triangle)
DG0 = fem.FunctionSpace(msh, ("Discontinuous Lagrange", 0))
num_cells = msh.topology.index_map(msh.topology.dim).size_local+ msh.topology.index_map(msh.topology.dim).num_ghosts
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
num_edges = msh.topology.index_map(1).size_local + msh.topology.index_map(1).num_ghosts
cell_list = np.arange(num_cells, dtype=np.int32)
edge_list = np.arange(num_edges, dtype=np.int32)
tdim = msh.topology.dim # Topology dimension
fdim = tdim - 1            # Facet dimension 
msh.topology.create_connectivity(fdim, tdim)
e4c_list = msh.topology.connectivity(tdim, fdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(msh.topology)
edge_list1 = np.delete(edge_list, boundary_facets)
sort_index = np.argsort(edge_list1)
cell_tag = meshtags(msh, msh.topology.dim, cell_list, cell_list)
edge_tag = meshtags(msh, 1, edge_list1[sort_index], edge_list1[sort_index])
dS = Measure("dS", domain=msh, subdomain_data = edge_tag)
dx = Measure("dx", domain=msh, subdomain_data = cell_tag)
n = FacetNormal(msh)
h = ufl.CellDiameter(msh)
h_avg = avg(h)
q = ufl.TestFunction(DG0)
for i in edge_tag.values:
    a111 = assemble_scalar(form(inner(1, q("+")) * dS(i)))
    print(f"a111_{i}",a111) 
    a11 =  assemble_vector(form(inner(1, q("+")) * dS(i))) 
    print(f"a11_{i}", a11.array, np.shape(a11.array))

a111_1 70.71067811865476
a11_1 [70.71067812  0.          0.          0.          0.          0.
  0.          0.        ] (8,)
a111_2 50.0
a11_2 [ 0. 50.  0.  0.  0.  0.  0.  0.] (8,)
a111_5 50.0
a11_5 [ 0. 50.  0.  0.  0.  0.  0.  0.] (8,)
a111_7 70.71067811865476
a11_7 [ 0.          0.         70.71067812  0.          0.          0.
  0.          0.        ] (8,)
a111_8 70.71067811865476
a11_8 [ 0.          0.          0.         70.71067812  0.          0.
  0.          0.        ] (8,)
a111_9 50.0
a11_9 [ 0.  0.  0.  0. 50.  0.  0.  0.] (8,)
a111_10 50.0
a11_10 [ 0.  0.  0.  0.  0. 50.  0.  0.] (8,)
a111_13 70.71067811865476
a11_13 [ 0.          0.          0.          0.          0.          0.
 70.71067812  0.        ] (8,)

Obviously the values ​​of the second and fifth sides are placed at the same position.

2

Have you searched this forum e.g. with keywords “interior”/“internal”, and sorted by date (most recent first)?

This topic has been discussed several times, including as recently as the last couple of weeks. Please make an effort to look for earlier posts.

3

Sorry for the late reply, thank you, I will check the relevant content of the forum carefully。