Block Matrices and varf

General Discussion
1

Hi,

I’m trying to solve an eigenvalue problem for a matrix that involves second-order derivatives, such as:

Screenshot 2025-06-22 at 12.26.56 PM
Screenshot 2025-06-22 at 12.26.56 PM696×538 12.2 KB

My initial approach was to implement each component of the matrix individually using separate varf definitions and then assemble the full matrix from these components. However, I suspect that this method introduces inconsistencies since the matrices are being constructed independently and may not align correctly.

What would be the correct way to construct a single, unified (block) matrix for this problem?

Thank you.

2

I understand from your post that x and y represent unknowns in a fespace. Then I don’t see any problem that could come from writing a varf for each block and then assembling the matrix by blocks.

3

Sorry for the confusion. Let me elaborate on the problem.
I am trying to solve the Luttinger-Kohn Equations over the finite element space. The matrix at each FE point is given by

Screenshot 2025-06-24 at 12.27.39 PM
Screenshot 2025-06-24 at 12.27.39 PM1060×756 39.7 KB

I have attached the code I tried to implement using the varf formulation by defining 6x6 finite element space, but not sure if it is the correct way to implement this.

int n = 10;
mesh3 Th = cube(n, n, n);

fespace Vh(Th, [P1, P1, P1, P1, P1, P1]);
Vh [u1, u2, u3, u4, u5, u6], [v1, v2, v3, v4, v5, v6];
//fespace Vh(Th, P1, complex);

//Vh u1, u2, u3, u4, u5, u6;
//Vh v1, v2, v3, v4, v5, v6;

// Constants (scaled units)
real gamma1 = 6.85;
real gamma2 = 2.1;
real gamma3 = 2.9;
real prefactor = 1.0; // Set hbar^2/2m = 1 for scaled units
real spindelta = 0.1;

// Weak form of full 6x6 LK Hamiltonian in 3D
varf aLK([u1,u2,u3,u4,u5,u6],[v1,v2,v3,v4,v5,v6]) =
int3d(Th)(
// Diagonal terms: P ± Q
prefactor * (gamma1 *(dx(u1)*dx(v1) + dy(u1)*dy(v1) + dz(u1)dz(v1)) + gamma2(-dx(u1)*dx(v1) - dy(u1)dy(v1) +2dz(u1)*dz(v1))) +// H11
prefactor * (gamma1 *(dx(u2)*dx(v2) + dy(u2)*dy(v2) + dz(u2)dz(v2)) + gamma2(-dx(u2)*dx(v2) - dy(u2)dy(v2) +2dz(u2)*dz(v2))) + // H22
prefactor * (gamma1 *(dx(u3)*dx(v3) + dy(u3)*dy(v3) + dz(u3)dz(v3)) + gamma2(-dx(u3)*dx(v3) - dy(u3)dy(v3) +2dz(u3)*dz(v3))) + // H33
prefactor * (gamma1 *(dx(u4)*dx(v4) + dy(u4)*dy(v4) + dz(u4)dz(v4)) + gamma2(-dx(u4)*dx(v4) - dy(u4)dy(v4) +2dz(u4)*dz(v4)))+ // H44
prefactor * gamma1 *(dx(u5)*dx(v5) + dy(u5)*dy(v5) + dz(u5)*dz(v5)) + // H55
prefactor * gamma1 *(dx(u6)*dx(v6) + dy(u6)*dy(v6) + dz(u1)*dz(v6)) - // H66

// Off-diagonal S terms: (-Sc)
1.732*prefactor*gamma3*(dx(u1)*dz(v2) -1i * dx(u1)*dz(v2)) - // H12 S*
1.732*prefactor*gamma3*(dx(u2)*dz(v1) +1i * dx(u2)*dz(v1)) + // H21

1.732*prefactor*(gamma3*(dx(u3)*dx(v1) - dy(u3)*dy(v1)) + 1i * 2*gamma3*dx(u3)*dy(v1)) - //H31

0.707*1.732*prefactor*gamma3*(dx(u3)*dz(v1) +1i*dx(u3)*dz(v1)) + //H51

1.414*1.732*prefactor*(gamma3*(dx(u6)*dx(v1)-dy(u6)*dy(v1)) + 1i*2*gamma3*dx(u6)*dy(v1)) + //H61

1.732*prefactor*gamma3*(dx(u2)*dz(v1) +1i*dx(u2)*dz(v1)) -//H42

1.414*prefactor*gamma2*(2*dz(u5)*dz(v2)-dx(u5)*dx(v2)-dy(u5)*dy(v2)) +//H52 u5v2

1.22*1.732*prefactor*gamma3*(dx(u6)*dz(v2) +1i*dx(u6)*dz(v2)) + //H62 u6v2 sqrt(3/2)*S

1.732*prefactor*(gamma3*(dx(u1)*dx(v3) - dy(u1)*dy(v3)) - 1i * 2*gamma3*dx(u1)*dy(v3)) +    // H13 u1 v3 R*

1.732*prefactor*gamma3*(dx(u4)*dz(v3) +1i * dx(u4)*dz(v3)) + //H43 u(4)v(3)S

1.22*1.732*prefactor*gamma3*(dx(u5)*dz(v3) -1i * dx(u5)*dz(v3)) + //sqrt(3/2)S* u(5)v(3)

1.414*prefactor*gamma2*(2*dz(u6)*dz(v3)-dx(u5)*dx(v3)-dy(u6)*dy(v3))+ //H63 sqrt(2)*Q u(6)v(3)

1.732*prefactor*(gamma3*(dx(u1)*dx(v4) - dy(u1)*dy(v4)) - 1i * 2*gamma3*dx(u1)*dy(v4)) + //H24 u(1) v(4) R*

1.732*prefactor*gamma3*(dx(u3)*dz(v4) -1i * dx(u3)*dz(v4)) - //H34 u(3)v(4) S*

1.414*1.732*prefactor*(gamma3*(dx(u5)*dx(v4) - dy(u5)*dy(v4)) - 1i * 2*gamma3*dx(u5)*dy(v4)) - // H54 u(5)v(4) -sqrt(2)R*

0.707*1.732*prefactor*gamma3*(dx(u6)*dz(v4) -1i * dx(u6)*dz(v4))- // H64 u(6)v(4) -S*/sqrt(2)

0.707*1.732*prefactor*gamma3*(dx(u1)*dz(v5) -1i * dx(u1)*dz(v5)) -//H15 u(1)v(5) -S*/sqrt(2))

1.414*prefactor*gamma2*(2*dz(u2)*dz(v5)-dx(u2)*dx(v5)-dy(u2)*dy(v5))+  //H25 u(2)v(5)    -sqrt(2)Q*

1.22*1.732*prefactor*gamma3*(dx(u3)*dz(v5) +1i*dx(u3)*dz(v5)) - //H35 u(3)v(5) sqrt(3/2)*S

1.414*1.732*prefactor*(gamma3*(dx(u4)*dx(v5) - dy(u4)*dy(v5)) + 1i * 2*gamma3*dx(u4)*dy(v5))+ //H45 u(4)v(5) -sqrt(2)R

1.414*1.732*prefactor*(gamma3*(dx(u1)*dx(v6) - dy(u1)*dy(v6)) - 1i * 2*gamma3*dx(u1)*dy(v6))+ //H16 u1) u6 sqt(2)*R* 

1.22*1.732*prefactor*gamma3*(dx(u2)*dz(v6) - 1i*dx(u2)*dz(v6)) + //H26 u2v6 sqrt(3/2)*S*

1.414*prefactor*gamma2*(2*dz(u3)*dz(v6)-dx(u3)*dx(v6)-dy(u3)*dy(v6)) - //H36 u3v6 sqrt(2)Q*

0.707*1.732*prefactor*gamma3*(dx(u4)*dz(v6) +1i * dx(u4)*dz(v6)) // H46 u(4)v(6) -S/sqrt(2)

);

// Mass matrix
varf m([u1,u2,u3,u4,u5,u6],[v1,v2,v3,v4,v5,v6]) =
int3d(Th)( u1v1 + u2v2 + u3v3 + u4v4 + u5v5 + u6v6 );

// Assemble
matrix A = aLK(Vh, Vh);
matrix M = m(Vh, Vh);
Thank you

4

I have not checked the details :sweat_smile:, but in the principle it looks good.