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/b820c2e0aaf525034498f9740f19cb >bd7b172056/ 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.
>>>
>>
>>
>