Robert,

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

On 06/19/2013 12:56 PM, Alexander Mihailov wrote:
> Yes, Robert, you are right, regions were empty.
>
> After correcting regions and setting post_process_hook, I have got some
> solution, slowly but surely it becomes believable, but the matrix is still
> singular and residual is high...
>
> sfepy: nls: iter: 0, residual: 1.528953e-02 (rel: 1.000000e+00)
> warning: (almost) singular matrix! (estimated cond. number: 6.06e+15)
> sfepy:   rezidual:    0.01 [s]
> sfepy:      solve:    0.01 [s]
> sfepy:     matrix:    0.00 [s]
> sfepy: linear system not solved! (err = 1.782581e-02 < 1.000000e-10)
> sfepy: linesearch: iter 1, (3.34778e-02 < 1.52894e-02) (new ls:
> 1.000000e-01)
> sfepy: nls: iter: 1, residual: 1.435926e-02 (rel: 9.391562e-01)
> sfepy: equation "tmp":
> sfepy: ev_grad.1.Omega( u )
> sfepy: updating materials...
> sfepy: ...done in 0.00 s
>
> There are still some subtleties that need to be considered?

The scalar field needs to be fixed somewhere, at least in a single vertex, as:

regions = {
     'Node' : ('node 0', {}),
}

ebcs = {
     'fix' : ('Node', {'u.0' : 0.0}),
}

Does it help?

r.

> <http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6husEB3Ixmr7KG5acBRI1Vq6DdhVj.jpg>
>
>
> среда, 19 июня 2013 г., 16:58:44 UTC+7 пользователь Robert Cimrman написал:
>>
>> On 06/19/2013 11:21 AM, Alexander Mihailov wrote:
>>> Thank you, Robert.
>>>
>>> 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.
>>
>> But it is not the same thing:
>> 1. the problem is more difficult to solve (1 more equation)
>> 2. the resulting G = grad(u) only in the weak sense, which is probably not
>> what
>> you want.
>>
>>> 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.
>>
>> First, check your region definitions by using:
>>
>> ./simple.py my.py --save-regions-as-groups --solve-not
>>
>> and inspect the resulting my_regions.vtk - it seems that the regions used
>> in
>> the material 'm' definition are empty. Maybe you want to use 'elements of
>> group
>> i' selector, where i = 1, 2, 3?
>>
>> r.
>>
>>> <
>> http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7ZoV5W4BY3DqAKi5o.jpg>
>>
>>>
>>>
>>> среда, 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-elementsbut
>>>>> 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.
>>>>>
>>>>
>>>>
>>>
>>
>>
>