Integral Matrix

General
1

Suppose “mesh” and “fes” are defined.

When I want to intergrate the multiple of two GridFunctions, it works with the following

B = GridFunction(fes)
H.Set(sin(x))
H = GridFunction(fes)
B.Set(x)
BH = Integrate(B*H*dx, mesh)

When I want to obtain the integral matrix with one specific function multiplied by an unknown funtion, what comes in my mind is the following, but it is really slow.
Anyone has better idea of doing it?

B = GridFunction(fes)
B.Set(sin(x))
S = zeros((len(B.vec)),dtype=float64)
for n in range(len(B.vec)):
	H = GridFunction(fes)
	H.vec[n] = 1
	S[n] = Integrate(B*H*dx, mesh)
Integral = np.linalg.dot(S, Unknown_Vecotr)
2

If I understand you correctly you want to do:

f = LinearForm(sin(x) * v * dx).Assemble()
Integral = InnerProduct(f.vec, Unknown_Vector)
1 Like
3

Thank you, Christpher!
That was exactly what I wanted.

By the way, I have defined two spaces and tried to optain the interploation matrix.
The following script was ok. I could obtain the consistent B3 with B2.

fes1 = H1(mesh, order=3)
u1 = fes1.TrialFunction()
v1 = fes1.TestFunction()

print(f'H1 space DoF = {fes1.ndof}')

fes2 = L2(mesh, order=0)
u2 = fes2.TrialFunction()
v2 = fes2.TestFunction()
print(f'L2 space DoF = {fes2.ndof}')

f22 = BilinearForm( v2*u2 *dx).Assemble()
print(shape(f22.mat))

f21 = BilinearForm( v2*u1 *dx).Assemble()
print(shape(f21.mat))

B1 = GridFunction(fes1)
B1.Set((2*x/d))

B2 = GridFunction(fes2)
B2.Set((2*x/d))

B3 = GridFunction(fes2)
B3.vec.data = f22.mat.Inverse() * (f21.mat *  B1.vec)

However, when I tried to obtain the matrix, it did not work out.

A = f22.mat.Inverse() * f21.mat
A = InnerProduct(f22.mat.Inverse() , f21.mat)

What are alternative commands?

4

To compose two linear operators you write

A = f22.mat.Inverse() @ f21.mat

In the second example,

A = InnerProduct(f22.mat.Inverse() , f21.mat)

do you mean the Frobenius inner product of two large matrices ? This is not supported, you have to reduce it to matrix-vector operations, best with a MulitVector.

However, in your case the matrix f22.mat is a diagonal matrix. If you give diagonal=True to your bilinear-form, it is stored in a vector, rather than in a CSR-matrix. That simplifies the .Inverse. Was this only a special case, or is that the general picture ?

Joachim

1 Like
5

Thank you, Joachim.

However, in your case the matrix f22.mat is a diagonal matrix. If you give diagonal=True to your bilinear-form, it is stored in a vector, rather than in a CSR-matrix. That simplifies the .Inverse . Was this only a special case, or is that the general picture ?

I am trying to get used to NGSolve which means I have to understand both cases.
“diagonal=True” is something I need to know. I can not yet see the whole picture or structure of NGSolve.

6

it took me 20 years to see the whole picture of NGSolve :slight_smile:

For most people it should be fine to see a certain fraction - I think you are already good

7

I see slight differences betweeen Integrate and InnderProduct, but that is something I say it is negligible. Thank you.

2.0000000000000001e-01 :: Analytical Integral
1.9678481603659206e-01 :: Integral by Integrate
1.9678481603659276e-01 :: Integral by LinearForm and InnerProduct
1.9678481603659279e-01 :: Integral by LinearFrom and numpy.dot

from numpy import *
from ngsolve import *
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

fes = L2(mesh, order=0)
u = fes.TrialFunction()
v = fes.TestFunction()

Az = GridFunction(fes)
Az.Set(x**4)

S_0 = 1/5
print(f'Integral = {S_0: 20.16e}')

S_1 = Integrate(Az*dx, mesh)
print(f'Integral = {S_1: 20.16e}')

Mat = LinearForm(v * dx).Assemble()
S_2 = dot(Mat.vec, Az.vec)
print(f'Integral = {S_2: 20.16e} ')