Hi,
I was wondering how dolfinx integrates the matrix elements. Suppose we are dealing with Lagrange, order one functions and a weak form like
\iint \beta\,u\,v\,dx+\int \gamma\,u\,v \,ds
on a simple triangle ( (0,0),(1,0),(0,1) ). The line integral goes along the first 2 points. How is the integration finally done? I guess with gauss quadrature, as the result is correct when interpolating \beta,\gamma in DG0 or CG1 functions. Is it possible to change integration settings?
Below an example, computing the weak form in dolfinx and the second part computes the matrix elements with sympy.
from mpi4py import MPI
import dolfinx as dox
import numpy as np
import ufl
# make single triangle from point and cell list
p=np.array([[0,0],[1,0],[0,1]])
t=np.array([[0,1,2]])
gdim = 2
shape = "triangle"
degree = 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
cells = np.array(t, dtype=np.int32)
tri = dox.mesh.create_mesh(MPI.COMM_WORLD, cells, p, domain)
# Defintion function space
V = dox.fem.FunctionSpace(tri, ("Lagrange", 1))
Q=dox.fem.FunctionSpace(tri, ("DG", 0))
# defintion functions
x=ufl.SpatialCoordinate(tri)
beta=1/(x[0]+1)
gamma=1/(x[0]+x[1]+1)
# boundary
tri.topology.create_connectivity(1,0)
edges = tri.topology.connectivity(1,0)
e1 = dox.mesh.locate_entities_boundary(tri, 1 , lambda X: np.isclose(X[1],0) )
all_edges=np.arange(len(edges),dtype=np.int32)
marker = np.zeros(len(edges),dtype=np.int32)
marker[e1]=1
edge_tags = dox.mesh.meshtags(tri, 1 , all_edges ,marker)
ds = ufl.Measure("ds", domain=tri, subdomain_data=edge_tags)
# weak Form
# funktion
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
K=beta*u*v*ufl.dx + gamma*u*v*ds(1)
#K=gamma*u*v*ds(1)
Kform = dox.fem.form(K)
A = dox.fem.petsc.assemble_matrix(Kform)
A.assemble()
C = A.convert("dense")
print("Matrix : \n",C.getDenseArray(),"\n")
# Computation using sympy
import sympy as sy
x,y,t = sy.symbols('x y t')
L=[1-x-y,x,y]
N=[1-t,t]
beta=1/(x+1)
#gamma=1/(x+y+1), x=t, y=0 along integration line
gamma=1/(t+1)
T1=np.zeros((3,3))
for i,Li in enumerate(L):
for j,Lj in enumerate(L):
T1[i,j]=sy.integrate(sy.integrate(beta*Li*Lj,(y,0,1-x)),(x,0,1))
print(T1)
T2=np.zeros((3,3))
for i,Ni in enumerate(N):
for j,Nj in enumerate(N):
T2[i,j]=sy.integrate(gamma*Ni*Nj,(t,0,1))
print(T2)
print(T1+T2)