Thanks for the informative links. I read the discussion fully, but I doubt if my problem was addressed there.
Dr. Mosher mentioned that primarily the Gain matrix is computed using an infinity (non-physical) reference, but what I concern here is that in computing FEM head model (either by Simbio or DUNEuro) we first consider the first channel as reference and set its share in all sources as zero no matter where it is located and how is connected to all sources in the brain anatomy, then after completing the computation of Gain matrix we apply the average-reference operator to give a bit of a value to the first channel thanks to the non-zero average value of every column (other channels' share in sources). In contrary, when computing a head model like 3 layers spheres we do not need to trivially consider a channel as reference and as a result no such problem will arise. So that's where I say computing the shares is nonlinear and I am skeptical if we can ignore one channel first and then try to rebuild it at last.
I know my question sounds complicated and I myself have no clue how to either mathematically or spatially explain what happens in details while computing FEM head model like this, but only thought maybe someone had any resources for this.
Thanks again.