How to obtain the basis functions in DG space

[Chen-03]

Hi everyone,

I am solving a complex problem that requires computations on each basis in the DG space. Is there any way to accomplish this?

I have read the content of the link Accessing the basis function of a Finite Element, but it is about querying CG.

I tried using the following approach:

from dolfin import *
from petsc4py import PETSc
import numpy as np
import sympy as sp
from sympy import pi , sin ,cos
import matplotlib.pyplot as plt
import math

RectangleMesh(Point(0, 0), Point(1, 1), size, size) #size=2, 4, 8

P1 = VectorElement("DG", triangle, pdeg) #pdeg=1, 2, 3
P2 = VectorElement("DG", triangle, pdeg)
P3 = FiniteElement("DG", triangle, pdeg)
P1P2P3 = MixedElement([P1, P2, P3])
W = FunctionSpace(mesh, P1P2P3)

V1 = FunctionSpace(mesh, P1)
V2 = FunctionSpace(mesh, P2)
V3 = FunctionSpace(mesh, P3)

(q,sigma,u) = TrialFunctions(W)
(zeta,tau,v) = TestFunctions(W)

n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

for j in range(V3.dim()):
    phi_j = Function(V3)
    phi_j.vector()[j] = 1

However, this seems to only work in certain situations.

[dokken]

@Chen-03

I have read the content of the link Accessing the basis function of a Finite Element , but it is about querying CG.

What happens if you try the same approach with DG?

@Chen-03

However, this seems to only work in certain situations.

What do you mean by this? Do you have an example where this doesn’t work? If it doesn’t work, what output do you get?

[Chen-03]

I am working on an LDG/minimization problem, and I want to solve the following equation for the trial function u and the test function v in the space V_3. Here, \Sigma_h is V_1, D_{DG} belongs to V_1 and is the function that needs to be solved. However, I cannot directly substitute u and v into the equation, so I thought of using the looping method mentioned above to solve it.

\begin{aligned} \forall \boldsymbol{\zeta}_h \in \Sigma_h, & \left(D_{D G}(v ; g), \boldsymbol{\zeta}_h\right)_{\Omega} \\ = & -\left(v, \nabla_h \cdot \boldsymbol{\zeta}_h\right)_{\Omega}+\left\langle\{\{v\}\}+\boldsymbol{C}_{12} \cdot [[ v ]], [[ \boldsymbol{\zeta}_h ]]\right\rangle_{\Gamma^o}\\&+\left\langle v, \boldsymbol{\zeta}_h \cdot \boldsymbol{n}\right\rangle_{\Gamma^N}+\left\langle g, \boldsymbol{\zeta}_h \cdot \boldsymbol{n}\right\rangle_{\Gamma^D} .\end{aligned}
where
\boldsymbol{C}_{12} \cdot \boldsymbol{n} = \frac{1}{2} \operatorname{sign}\left((1,0)^T \cdot \boldsymbol{n}\right)

I tested this approach in SIP/DG by replacing (\nabla u, \nabla v)_{\Omega} with (\nabla \phi_i, \nabla \phi_j)_{\Omega}, then assembling them to solve the Poisson equation, and it seems to work fine.

[Chen-03]

Hi dokken.

Here are two ways I wrote dot(grad(u), grad(v)) * dx. I’m curious why these two results are the same. Can I consider phi as a basis function?

from dolfin import *
from petsc4py import PETSc
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

size = 2
pdeg = 3
mesh = UnitSquareMesh(size, size)
V = FunctionSpace(mesh, "DG", pdeg)
u = TrialFunction(V)
v = TestFunction(V)
alpha = Constant(10.0)
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

#test1
a1 = assemble(dot( grad(u) , grad (v))*dx)
a1 = as_backend_type(a1).mat()
a1 = a1.getValues(range(0, a1.getSize()[0]), range(0, a1.getSize()[1]))
print("A1:",a1)

#test2
Amat = np.zeros((V.dim(), V.dim()))
dofmap = V.dofmap()
global_dofs = dofmap.dofs()

# dot( grad (u), grad (v) )*dx 
for i in range(V.dim()):
    phi_i = Function(V)
    phi_i.vector()[i] = 1
    for j in range(V.dim()):
        phi_j = Function(V)
        phi_j.vector()[j] = 1
        phi_ij = assemble(inner(grad(phi_i), grad(phi_j)) * dx)
        Amat[global_dofs[i], global_dofs[j]] = phi_ij

print("A2:",Amat)

[dokken]

Well, they will of course give the same result, as by setting the ith and jth degree of freedom to one, you are making a global basis function.

This only works for specific spaces (Lagrange, DG), and not necessarily for Nedelec/RT, as covered in: Assemble with dolfin VS dolfinx - #2 by dokken

Note that the second approach is doing alot more work than the first approach.
It will integrate over the whole domain M**2 times (M=V.dim()), while the first approach does the integration once.

I still do not understand the need for the tabulated values of a basis function here?

You have also not explained what is wrong with the approach in:

Accessing the basis function of a Finite Element

Hello to everybody, I am working with a problem where I need to directly access the basis functions of a Finite Element. A minimal example would be: >>> from fenics import * >>> mesh = UnitSquareMesh(8,8) >>> V = FunctionSpace(mesh, 'P', 1) >>> V.dim() displays 81. Then I’d like to have a Python list of all these 81 functions, and being able to evaluate/use them. The precise reason beyond this is very long to describe, since it involves connections with other mathematical objects I am studyin…

[Chen-03]

Sorry for the unclear expression. I meant that I substituted \phi into this equation.

\begin{aligned} \forall \boldsymbol{\zeta}_h \in \Sigma_h, & \left(D_{D G}(v ; g), \boldsymbol{\zeta}_h\right)_{\Omega} \\ = & -\left(v, \nabla_h \cdot \boldsymbol{\zeta}_h\right)_{\Omega}+\left\langle\{\{v\}\}+\boldsymbol{C}_{12} \cdot[[v]],\left[\left[\boldsymbol{\zeta}_h\right]\right]\right\rangle_{\Gamma^o} \\ & \quad+\left\langle g, \boldsymbol{\zeta}_h \cdot \boldsymbol{n}\right\rangle_{\Gamma^D} . \end{aligned}

That is, substituting v with \phi , with g as the Dirichlet boundary condition, and solving D_{DG} through the Riesz representation.

The reason for doing this is that I am unable to directly substitute the trial function and test function to solve.

Additionally, I’m also curious about how to convert the Lambda functions obtained from ’ Accessing the basis function of a Finite Element ’ into the regular FEniCS functions we usually use.

[dokken]

@Chen-03

The reason for doing this is that I am unable to directly substitute the trial function and test function to solve.

Wouldn’t it be better to assemble the matrix where v is a TrialFunction once, and then apply a matrix-vector product with the appropriate unit vector (i.e. any of the vectors you created above).
Then the assembly is done once, and a matrix vector product isn’t that expensive.

[Chen-03]

Sorry, I don’t quite understand what you mean. Here’s the code I have written so far, but I’m not entirely sure if it’s correct.

from dolfin import *
import numpy as np

pdeg = 1
size = 2

mesh = RectangleMesh(Point(0, 0), Point(1, 1), size, size)

P1 = VectorElement("DG", triangle, pdeg)
P2 = VectorElement("DG", triangle, pdeg)
P3 = FiniteElement("DG", triangle, pdeg)
P1P2P3 = MixedElement([P1, P2, P3])
W = FunctionSpace(mesh, P1P2P3)

V1 = FunctionSpace(mesh, P1)
V2 = FunctionSpace(mesh, P2)
V3 = FunctionSpace(mesh, P3)

# (q,sigma,u) = TrialFunctions(W)
# (zeta,tau,v) = TestFunctions(W)

n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

g = Constant(0.0)
u = TrialFunction(V3)

phi_j = Function(V3)
phi_j.vector()[0] = 1

print(" phi_j :",phi_j.vector().get_local())
def Discrete_weak_gradient(v,g):
  (q,sigma,uu) = TrialFunctions(W)
  (tau,zeta,vv) = TestFunctions(W)

  v = project(v,V3)
  g = project(g,V3)

  C12 = Constant((1.0,0))
  dot_product = inner(C12,n('+'))
  sign = conditional(dot_product > 0, 1.0, -1.0)

  # bc = DirichletBC(W.sub(2), g, "on_boundary")

  L = inner(q,zeta)*dx

  R = -v*div(zeta)*dx + avg(v)*inner(jump(zeta),n('+'))*dS\
    + inner(0.5*sign,jump(v))*inner(jump(zeta),n('+'))*dS + g*inner(zeta,n)*ds

  F = L - R
  L = lhs(F)
  R = rhs(F)

  w = Function(W)
  solve(L==R,w)
  (q,sigma,u) = w.split()

  return q


q = Discrete_weak_gradient(phi_j,g)
print(" Dwg phi_j :",q.vector().get_local())

The problem I’m encountering is that I don’t know how to handle the case with

v = TrialFunction}(V3)

and

v = TestFunction(V3)

.

[dokken]

Consider the minimal example that illustrates an efficient apporach for your assembly:

from dolfin import *
import time
import numpy as np

mesh = UnitSquareMesh(50, 50)

V = FunctionSpace(mesh, "CG", 1)


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

a = inner(jump(grad(u)), avg(grad(v))) * dS

uh = Function(V)
L = inner(jump(grad(uh)), avg(grad(v))) * dS

matrix_time = 0
start = time.perf_counter()
A = assemble(a)
end = time.perf_counter()
matrix_time += end - start

vector_time = 0

L1 = Vector()
for i in range(V.dim()):
    uh.vector()[:] = 0.0
    uh.vector()[i] = 1.0

    vector_start = time.perf_counter()
    L0 = assemble(L)
    vector_end = time.perf_counter()
    vector_time += vector_end - vector_start

    matrix_start = time.perf_counter()
    A.mult(uh.vector(), L1)
    matrix_end = time.perf_counter()
    matrix_time += matrix_end - matrix_start
    np.testing.assert_allclose(L0.get_local(), L1.get_local(), atol=1e-10)

print("Matrix assembly time: ", matrix_time)
print("Vector assembly time: ", vector_time)

which yields

Matrix assembly time:  0.14058589899786966
Vector assembly time:  7.7443401540010655

i.e. assembling the matrix with a trial function, then using the matrix vector product to get the appropriate rhs vector is way faster than assembling each vector.

[Chen-03]

Hi Dokken,

Sorry for the delayed response. I roughly understand your point, but I’m not sure how to apply it to my situation

I want to solve the following equation, but I’m not sure how to implement it using the method you suggested.

\text { find } w_h \in V_h \text { such that } \forall v_h \in V_h \text { : }

-J_h^{\prime}\left(u_h\right)\left(v_h\right)=\left\{\begin{array}{l} \left(D_{D G}\left(w_h ; 0\right), D_{D G}\left(v_h ; 0\right)\right)_{\Omega} \\ +\left\langle\eta h^{-1} [[ w_h\right],\left[[v_h ]]\right\rangle_{\Gamma^o}+\left\langle\eta h^{-1} w_h, v_h\right\rangle_{\Gamma^D}, \end{array}\right.

where

\begin{aligned} \forall \boldsymbol{\zeta}_h \in \Sigma_h, & \left(D_{D G}(v ; g), \boldsymbol{\zeta}_h\right)_{\Omega} \\ = & -\left(v, \nabla_h \cdot \boldsymbol{\zeta}_h\right)_{\Omega}+\left\langle\{\{v\}\}+\boldsymbol{C}_{12} \cdot [[ v ]], [[\boldsymbol{\zeta}_h ]]\right\rangle_{\Gamma^{\circ}}+\left\langle g, \boldsymbol{\zeta}_h \cdot \boldsymbol{n}\right\rangle_{\Gamma^D} . \end{aligned}

from dolfin import *
from petsc4py import PETSc
import numpy as np
import matplotlib.pyplot as plt
import math


x, y = sp.symbols('x[0] x[1]')
p = 2.0
test = 2
pdeg = 1
size = 2

P1 = VectorElement("DG", triangle, pdeg)
P2 = VectorElement("DG", triangle, pdeg)
P3 = FiniteElement("DG", triangle, pdeg)
P1P2P3 = MixedElement([P1, P2, P3])
W = FunctionSpace(mesh, P1P2P3)

V1 = FunctionSpace(mesh, P1)
V2 = FunctionSpace(mesh, P2)
V3 = FunctionSpace(mesh, P3)

(q,sigma,u) = TrialFunctions(W)
(zeta,tau,v) = TestFunctions(W)

n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

def Discrete_weak_gradient(v,g):
  (q,sigma,uu) = TrialFunctions(W)
  (zeta,tau,vv) = TestFunctions(W)
  # v = interpolate(v,V3)
  g = project(g,V3) 

  C12 = Constant((1.0,0))
  dot_product = inner(C12,n('+'))
  sign = conditional(dot_product > 0, 1.0, -1.0)
  # sign = conditional(dot_product > 0, 1.0,conditional(dot_product < 0, -1.0, 0.0) )

  # bc = DirichletBC(W.sub(2), g, "on_boundary")

  L = inner(q,zeta)*dx

  R = -v*div(zeta)*dx + avg(v)*inner(jump(zeta),n('+'))*dS\
    + inner(0.5*sign,jump(v))*inner(jump(zeta),n('+'))*dS + g*inner(zeta,n)*ds

  F = L - R
  L = lhs(F)
  R = rhs(F)

  w = Function(W)
  solve(L==R,w)
  (q,sigma,u) = w.split()

  return q