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

Hmm, it might be caused by the second equation. I would really just solve the
Poisson equation and used the post-processing. But the solver converged even if
it complained - did you look at the solution?

r.

On 06/19/2013 05:40 PM, Alexander Mihailov wrote:
> Robert,
>
> I have set boundary conditions, but the matrix is still singular.
>
> regions = {
> 'Top' : ('nodes in (z > 0.9999)', {}),
> 'Bottom' : ('nodes in (z < -0.9999)', {}),
> 'Right' : ('nodes in (y > 0.9999)', {}),
> 'Left' : ('nodes in (y < -0.9999)', {}),
> 'Far' : ('nodes in (x > 0.9999)', {}),
> 'Near' : ('nodes in (x < -0.9999)', {}),
> }
>
> ebcs = {
> 'fix0' : ('Top', {'u.0' : 0.0},),
> 'fix1' : ('Bottom', {'u.0' : 0.0},),
> 'fix2' : ('Left', {'u.0' : 0.0},),
> 'fix3' : ('Right', {'u.0' : 0.0},),
> 'fix4' : ('Far', {'u.0' : 0.0},),
> 'fix5' : ('Near', {'u.0' : 0.0},),
> }
>
> sfepy: nls: iter: 0, residual: 1.528953e-01 (rel: 1.000000e+00)
> warning: (almost) singular matrix! (estimated cond. number: 3.76e+15)
> sfepy: rezidual: 0.00 [s]
> sfepy: solve: 0.00 [s]
> sfepy: matrix: 0.00 [s]
> sfepy: nls: iter: 1, residual: 4.692475e-17 (rel: 3.069077e-16)
> sfepy: equation "tmp":
> sfepy: ev_grad.1.Omega( u )
> sfepy: updating materials...
> sfepy: ...done in 0.00 s
>
>
>
> среда, 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.
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>
>>>>
>>>
>>
>>
>