# Block Matrices and varf
**URL:** https://community.freefem.org/t/block-matrices-and-varf/3980
**Category:** General Discussion
**Created:** [June 22, 2025, 4:31pm UTC](https://community.freefem.org/t/block-matrices-and-varf/3980 "2025-06-22T16:31:43Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@kelvindsou91](https://community.freefem.org/u/kelvindsou91)
#### Post date: [June 22, 2025, 4:31pm UTC](https://community.freefem.org/t/block-matrices-and-varf/3980/1 "2025-06-22T16:31:44Z")
Hi,
I’m trying to solve an eigenvalue problem for a matrix that involves second-order derivatives, such as:
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.
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 24, 2025, 10:50am UTC](https://community.freefem.org/t/block-matrices-and-varf/3980/2 "2025-06-24T10:50:57Z")
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.
---
### Author:  [@kelvindsou91](https://community.freefem.org/u/kelvindsou91)
#### Post date: [June 24, 2025, 4:37pm UTC](https://community.freefem.org/t/block-matrices-and-varf/3980/3 "2025-06-24T16:37:26Z")
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
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) +2_dz(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) +2_dz(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) +2_dz(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) +2_dz(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)( u1_v1 + u2_v2 + u3_v3 + u4_v4 + u5_v5 + u6_v6 );
// Assemble
matrix A = aLK(Vh, Vh);
matrix M = m(Vh, Vh);
Thank you
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 25, 2025, 7:40am UTC](https://community.freefem.org/t/block-matrices-and-varf/3980/4 "2025-06-25T07:40:02Z")
I have not checked the details 😅, but in the principle it looks good.