Re: [Getfem-users] High level generic assembly procedures

[Yves Renard]

Dear Marco,

May be it is too much specific to add a new function to the langage (you
may if you think it is better). However, you also may use
(1 -cos(t))/t^2 = (sinc(t))^2 / (1+cos(t))


Yves.





Le 14/04/2014 11:23, Marco Pischedda a écrit :
> Dear Yves,
>
> thank you for your answer.
> Is there also something for the function (1-cos(phi))/(phi*phi) when
> phi tends to zero?
>
> Thanks in advance
>
> Marco
>
>
>
> 2014-04-13 10:44 GMT+02:00 Yves Renard <address@hidden>:
>> Dear Marco,
>>
>> The message "WARNING: detected wrong equivalent nodes" is a debugging 
>> message which can be ignored, it serves to control the efficiency of the 
>> hash codes for the semantic trees produced by the generic assembly.
>>
>> For the moment, their is no conditional expressions in the language. 
>> However, you can use pos_part and neg_part.
>> Concerning the cardinal sine function (sin(t)/t) I added it to the language. 
>> It is called "sinc" now in the generic assembly.
>>
>> Yves.
>>
>>
>>
>>
>>
>> ----- Original Message -----
>> From: "Marco Pischedda" <address@hidden>
>> To: "Yves Renard" <address@hidden>
>> Cc: address@hidden
>> Sent: Friday, April 11, 2014 9:28:22 AM
>> Subject: Re: [Getfem-users] High level generic assembly procedures
>>
>> Dear Yves,
>>
>> thank you for your answer.  I understand how to impose the boundary 
>> condition.
>> I have another question: in my formulation I'm assembling a term
>> A(rho)*f.Test_f, where A(rho) is a symmetric matrix defined as
>> follows:
>>
>>  " [-sqr(rho(2))-sqr(rho(3)), rho(1)*rho(2), rho(1)*rho(3);"
>>  " rho(1)*rho(2), -sqr(rho(1))-sqr(rho(3)), rho(2)*rho(3);"
>>  " rho(1)*rho(3), rho(2)*rho(3) , -sqr(rho(1))-sqr(rho(2))]";
>>
>> where rho and f are two FemVariables. When I'm assembling this term I
>> receive this warning message:  "WARNING: detected wrong equivalent
>> nodes".
>> What does it mean?
>>
>> Then I have an expression "sin(phi)/phi", how can I insert an "if
>> condition" (i.e if phi<1e-8) in the expression when "phi" tends to
>> zero?
>>
>> Thanks in advance
>>
>> Marco
>>
>>
>>
>> 2014-04-08 15:29 GMT+02:00 Yves Renard <address@hidden>:
>>>
>>> Dear Marco,
>>>
>>> No, "Grad_u.Test_p +p.Grad_Test_u+p.Test_u + p(0)u(0)-p(L)u(L)" is not
>>> valid. The expression of an assembly string is to be evaluated at each
>>> Gauss point of the concerned region. An expression such as "p(0)" will
>>> be understand to be the component 0 of p which is not valid because the
>>> first component is 1.
>>> If you want to add some expression on a boundary (an extremity of the
>>> intervall [0, L] here) then it is possible but not in the same assembly
>>> string. In your example you have to add the assembly string "p.u" at the
>>> boundary x=0 and "-p.u" at the boundary x=L.
>>>
>>> However, if you need to prescribe a Dirichlet boundary condition, the
>>> best would be to use a corresponding model brick (if you use the model
>>> system).
>>>
>>>
>>> Yves.
>>>
>>>
>>>
>>>
>>>
>>> Le 07/04/2014 15:43, Marco Pischedda a écrit :
>>>> Hi,
>>>>
>>>> how can I impose the boundary condition with high level generic
>>>> assembly? For example I want to assemble "Grad_u.Test_p +
>>>> p.Grad_Test_u+p.Test_u" and then imposing the value of p on the
>>>> boundaries of the 1d domain, i.e x=0 and x=L. It is possible
>>>> to write "Grad_u.Test_p +p.Grad_Test_u+p.Test_u + p(0)u(0)-p(L)u(L)",
>>>> i.e can I insert the boundaries conditions in the assembling string?
>>>>
>>>> Thanks in advance
>>>>
>>>> Marco
>>>>
>>>>
>>>> 2014-04-03 9:52 GMT+02:00 Marco Pischedda <address@hidden>:
>>>>> Dear Yves,
>>>>>
>>>>> thank you for bug correction, now it works.
>>>>>
>>>>>> Concerning your question on Grad_u for a vector field, it is usually an
>>>>>> order two  tensor (a matrix) except in 1D.
>>>>>   Ok that's good.
>>>>>
>>>>> I will let you know if there are other problems.
>>>>>
>>>>> Thank you
>>>>>
>>>>> Marco
>>>>>
>>>>> 2014-04-03 9:25 GMT+02:00 Yves Renard <address@hidden>:
>>>>>> Dear Marco,
>>>>>>
>>>>>> Concerning your question on Grad_u for a vector field, it is usually an
>>>>>> order two  tensor (a matrix) except in 1D.
>>>>>> I tried to to the best to make the most operations "dimension
>>>>>> indepedent" so there is some permitivity in the langage (components of
>>>>>> size 1 are sometimes ignored).
>>>>>> Remember also that you can have the expression of any term with the
>>>>>> commant Print. For instance "Print( Grad_u).Test_p" will print the
>>>>>> gradient of u on each Gauss point. It does a lot of print, but at least,
>>>>>> you can see the format of the term.
>>>>>>
>>>>>> Yves.
>>>>>>
>>>>>>
>>>>>> Le 02/04/2014 17:21, Marco Pischedda a écrit :
>>>>>>> Hi,
>>>>>>>
>>>>>>> I have other questions:
>>>>>>>
>>>>>>> - I tried to assemble separately the following terms:
>>>>>>>   "Grad_u.Test_p" , "p.Grad_Test_u", "p.Test_u" and the code works 
>>>>>>> correctly.
>>>>>>>   Then I want to assemble the sum of this terms, i.e: "Grad_u.Test_p +
>>>>>>> p.Grad_Test_u+p.Test_u"
>>>>>>>   but I receive the following error:
>>>>>>>
>>>>>>>    Addition or substraction of incompatible expressions or of different 
>>>>>>> sizes
>>>>>>>    terminate called after throwing an instance of 'gmm::gmm_error'
>>>>>>>    what():  Error in getfem_generic_assembly.cc, line 3591 :
>>>>>>>    Error in assembly string
>>>>>>>    Aborted
>>>>>>>
>>>>>>> - I'm working in 1d problem with vectors as unknowns. Grad_u is
>>>>>>> therefore a vector or is a        tensor? When I do Grad_u.Test_p the
>>>>>>> result is a scalar or a vector? Test_p is a vector or a scalar?
>>>>>>>
>>>>>>> Thanks in Advance
>>>>>>>
>>>>>>> Marco
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> 2014-04-01 17:51 GMT+02:00 Marco Pischedda <address@hidden>:
>>>>>>>> Ok thank you,
>>>>>>>>
>>>>>>>> I tried it and it works. I will let you know if there are other 
>>>>>>>> problems.
>>>>>>>>
>>>>>>>> Thanks
>>>>>>>>
>>>>>>>> Marco
>>>>>>>>
>>>>>>>> 2014-04-01 17:16 GMT+02:00 Yves Renard <address@hidden>:
>>>>>>>>> Dear Marco,
>>>>>>>>>
>>>>>>>>> All seems to me correct in your implementation. This is probably just
>>>>>>>>> the test line 2074 of getfem_generic_assembly.cc
>>>>>>>>> which is not correct for one-dimensionnal problems.
>>>>>>>>>
>>>>>>>>> I think the test should be
>>>>>>>>>
>>>>>>>>>  GA_DEBUG_ASSERT((qdim == 1 && t.sizes()[0] == N) ||
>>>>>>>>>                       (t.sizes()[1] == N && t.sizes()[0] == qdim) ||
>>>>>>>>>                       (N == 1 && t.sizes()[0] == qdim),
>>>>>>>>>                       "dimensions mismatch");
>>>>>>>>>
>>>>>>>>> May be you can try this. I will validate it if it works.
>>>>>>>>>
>>>>>>>>> Regards,
>>>>>>>>>
>>>>>>>>> Yves.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Le 01/04/2014 16:56, Marco Pischedda a écrit :
>>>>>>>>>> Dear Yves,
>>>>>>>>>>
>>>>>>>>>> thank you for your fast answer.
>>>>>>>>>>
>>>>>>>>>> I have another question:
>>>>>>>>>>
>>>>>>>>>> - I have a monodimensional problem but the unknowns are vectorial. I
>>>>>>>>>> set the vectorial
>>>>>>>>>> dimension of the unknows with mf_u.set_qdim(3) and with 
>>>>>>>>>> mf_p.set_qdim(3).
>>>>>>>>>> When I define the expression " Grad_u.Test_p "  i receive the 
>>>>>>>>>> following error:
>>>>>>>>>>
>>>>>>>>>> terminate called after throwing an instance of 'gmm::gmm_error'
>>>>>>>>>>   what():  Error in
>>>>>>>>>> ../../getfem-svn/getfem/src/getfem_generic_assembly.cc, line 2074 :
>>>>>>>>>> dimensions mismatch
>>>>>>>>>>
>>>>>>>>>> Grad_u should be a 2nd order tensor, while Test_p should be a vector.
>>>>>>>>>> I suppose the
>>>>>>>>>> problem is that Grad_u is interpreted as a vector while Test_p is
>>>>>>>>>> interpreted as a scalar.
>>>>>>>>>> How can I use the high level generic assembly procedures on 
>>>>>>>>>> vectorial problems
>>>>>>>>>> defined on a monodimensional computational domain?
>>>>>>>>>>
>>>>>>>>>> Thanks in advance
>>>>>>>>>>
>>>>>>>>>> Marco
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> 2014-04-01 8:42 GMT+02:00 Yves Renard <address@hidden>:
>>>>>>>>>>> Dear Marco,
>>>>>>>>>>>
>>>>>>>>>>> Yes, your code is correct and it should work correctly.
>>>>>>>>>>> Note that you can also use the high level generic assembly with the
>>>>>>>>>>> model object using a generic assembly brick.
>>>>>>>>>>> Note also that high level assembly is only available in the svn
>>>>>>>>>>> repository version of Getfem
>>>>>>>>>>> (and a very recent version, at least r4570 because a bug have been
>>>>>>>>>>> corrected on coupled problems recently).
>>>>>>>>>>>
>>>>>>>>>>> Yves.
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> Le 31/03/2014 14:00, Marco Pischedda a écrit :
>>>>>>>>>>>> Dear all,
>>>>>>>>>>>>
>>>>>>>>>>>> I want to use the high level generic assembly procedures in order 
>>>>>>>>>>>> to
>>>>>>>>>>>> use the automatic differentation for calculating the Jacobian of a
>>>>>>>>>>>> non-linear problem.
>>>>>>>>>>>>
>>>>>>>>>>>> For example I want to calculate the Jacobian of two vectorial
>>>>>>>>>>>> equations in the unknowns "u" and "p". The first vectorial 
>>>>>>>>>>>> equation is
>>>>>>>>>>>> (u+p).Test_u while the second is (u+p).Test_p.
>>>>>>>>>>>>
>>>>>>>>>>>>  It is correct the following code?
>>>>>>>>>>>>
>>>>>>>>>>>>  gmm::sub_interval Iu(0, ndofs_u);
>>>>>>>>>>>>  gmm::sub_interval Ip(ndofs_u,ndofs_p);
>>>>>>>>>>>>  workspace.add_fem_variable("u", mf_u, Iu, U);
>>>>>>>>>>>>  workspace.add_fem_variable("p", mf_p, Ip, P);
>>>>>>>>>>>>  workspace.add_expression("(u+p).Test_u + (u+p).Test_p",mim);
>>>>>>>>>>>>  getfem::model_real_sparse_matrix Jac(ndofs_u+ndofs_p, 
>>>>>>>>>>>> ndofs_u+ndofs_p);
>>>>>>>>>>>>  workspace.set_assembled_matrix(Jac);
>>>>>>>>>>>>  workspace.assembly(2);
>>>>>>>>>>>>
>>>>>>>>>>>> In the matrix Jac I have the Jacobian of the system?
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> Thanks in advance
>>>>>>>>>>>>
>>>>>>>>>>>> Marco
>>>>>>>>>>>>
>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>> Getfem-users mailing list
>>>>>>>>>>>> address@hidden
>>>>>>>>>>>> https://mail.gna.org/listinfo/getfem-users
>>>>>>>>>>> --
>>>>>>>>>>>
>>>>>>>>>>>   Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
>>>>>>>>>>>   Pole de Mathematiques, INSA-Lyon             fax : (33) 
>>>>>>>>>>> 04.72.43.85.29
>>>>>>>>>>>   20, rue Albert Einstein
>>>>>>>>>>>   69621 Villeurbanne Cedex, FRANCE
>>>>>>>>>>>   http://math.univ-lyon1.fr/~renard
>>>>>>>>>>>
>>>>>>>>>>> ---------
>>>>>>>>>>>
>>>>>>>>> --
>>>>>>>>>
>>>>>>>>>   Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
>>>>>>>>>   Pole de Mathematiques, INSA-Lyon             fax : (33) 
>>>>>>>>> 04.72.43.85.29
>>>>>>>>>   20, rue Albert Einstein
>>>>>>>>>   69621 Villeurbanne Cedex, FRANCE
>>>>>>>>>   http://math.univ-lyon1.fr/~renard
>>>>>>>>>
>>>>>>>>> ---------
>>>>>>>>>
>>>>>> --
>>>>>>
>>>>>>   Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
>>>>>>   Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
>>>>>>   20, rue Albert Einstein
>>>>>>   69621 Villeurbanne Cedex, FRANCE
>>>>>>   http://math.univ-lyon1.fr/~renard
>>>>>>
>>>>>> ---------
>>>>>>
>>>
>>> --
>>>
>>>   Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
>>>   Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
>>>   20, rue Albert Einstein
>>>   69621 Villeurbanne Cedex, FRANCE
>>>   http://math.univ-lyon1.fr/~renard
>>>
>>> ---------
>>>


-- 

  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------