Modeling gravity field with finite elements

Question

Well, my question is rather from applied maths area, not pure mathematics, so I am not sure that it's a place on this board for one.

I want to solve a direct gravity gradiometry problem on 3D triangular mesh in the $[-1,1]^3$ cube.

Let the area $z>0$ be overground, $z<0$ --- underground. Now push a solid object into the ground and see, how changed the gravity field above the ground.

According to the divergence theorem, we have the field equation: $\nabla \bar{G} = \rho,$ where $\rho$ is the density of mass.

Next, we can build mesh. I used gmsh, my .geo file:

Point(1) = {-1, -1, 0, 1.0};
Point(2) = {-1, 1, 0, 1.0};
Line(1) = {2, 1};
Extrude {2, 0, 0} {
  Line{1};
}
Extrude {0, 0, 1} {
  Surface{5};
}
Extrude {0, 0, -1} {
  Surface{5};
}
Point(25) = {-0.2, -0.2, -0.4, 1.0};
Point(26) = {-0.2, 0.2, -0.4, 1.0};
Line(50) = {25, 26};
Extrude {0.4, 0, 0} {
  Line{50};
}
Extrude {0, 0, -0.4} {
  Surface{54};
}

Physical Volume(84) = {1};  // overground
Physical Volume(85) = {2};  // underground
Physical Volume(86) = {3};  // brick

So, after building mesh and choosing the test functions class, we can build a weak formulation of problem:

$$\int_\Omega \nabla \bar{G} \cdot \bar{v} = \int_\Omega \rho \cdot \bar{v}$$

Now we can run a fem solver. I used sfepy, my .py file:

import numpy as nm

filename_mesh = 'my.mesh'

regions = {
    'Omega'       : ('all', {}),
    'Overground'  : ('nodes of group 1', {}),
    'Underground' : ('nodes of group 2', {}),
    'Brick'       : ('nodes of group 3', {}),
}

field_1 = {
    'name' : 'gravity',
    'dtype' : nm.float64,
    'shape' : (3,),
    'region' : 'Omega',
    'approx_order' : 1,
}

variables = {
    'G'       : ('unknown field',   'gravity',  0 ),
    'g'       : ('test field',      'gravity', 'G'),
}

ebcs = {
}

materials = {
    'm' : ({'rho': {
              'Overground':  1.0e-7,
              'Underground': 1.0e+0,
              'Brick':       1.0e+5
           }},
    ),
    'n' : ({'G' : 1.0 }, )
}

equations = {
    'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_div.1.Omega( m.rho, g )"""
}

solvers = {
    'ls' : ('ls.scipy_direct', {}),
    'newton' : ('nls.newton', {
        'i_max'      : 1,
        'eps_a'      : 1e-10,
    }),
}

The problem is that the results a little bit unadequate. I expected the field, uniformly directed to the brick, but got that:

What I have to do to correct the model?

Thanks in advance.