Adding a new variable is a great way!

Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit and evident than hook magic.

field_1 = {
   'name' : 'gravity',
    'dtype' : nm.float64,
    'shape' : (3,),
    'region' : 'Omega',
    'approx_order' : 1,
}
field_2 = {
    'name' : 'gravity_potential',
    'dtype' : nm.float64,
    'shape' : (1,),
    'region' : 'Omega',
    'approx_order' : 1,
}
variables = {
    'G'       : ('unknown field',   'gravity',  1),
    'g'       : ('test field',      'gravity',  'G'),
    'u'       : ('unknown field',   'gravity_potential', 0),
    'v'       : ('test field',      'gravity_potential', 'u'),
}
equations = {
    'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""",
    'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """,
}

But here is another issue...

When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.

I think it's the wrong specifying or material, or the boundary conditions issue.

среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman написал:

Hi Alexander

On 06/19/2013 04:42 AM, Alexander Mihailov wrote:
> Hi all.
>
> I have already asked my question at
> http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-finite-elements
> and
> http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-finite-elements but
> seems there are no people there, who knows sfepy and finite elements, in
> particular.
>
> I want to model the gravity field on 3D rectangular area. It can be
> described by the equation: div(G) = rho. Here G is vector unknown function,
> rho is scalar parameter, which is constant at the selected point and fully
> determined by material.
>
> The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test
> function.
> Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) =
> dw_volume_lvf.1.Omega( m.rho, g )"""}

The terms you use do not correspond to the equation. As div(G) is a scalar, it
cannot be used with a vector test function, but a scalar test function - the
term that does that is dw_stokes, but then the matrix would not be square.

Try using a different approach - use the scalar gravity potential u as your
function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation,
use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields.
Check the equation in [1]. Then G could be computed in a postprocess hook
function by evaluating ev_grad term - see [2] for an example of a similar hook
function (post_process()).

r.

[1] http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html
[2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html

> I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to
> this message.
>
> When I try to run sfepy, it says that matrix is singular, so the solving
> process does not converge.
>
> What I have to do to correct the model?
>
> Thanks in advance.
>