Re: [Getfem-users] number of Newton-Cotes points

[Edouard Oudet]

OK, Thanks !

Le 18/01/2019 à 09:39, Yves Renard a écrit :
> Dear Edouard,
>
> Yes, in Newton-Cotes integration methods, some  quadature points of
> the regular grid correspond to vanishing weights and these points are
> indeed eliminated. It would not be difficult to have an integration
> method with the complete grid (this is an option of the add_point
> fonction of getfem_integration.cc, in fact)  but this is not available
> for the moment.
>
> Best regards,
>
> Yves
>
>
>
>
>
>
> Le 16/01/2019 à 20:19, Edouard Oudet a écrit :
>> Dear all,
>>
>> I studied a simple situation : a triangular mesh and defined a MeshIm
>> class based on Newton-Cotes methods of degree "degi". I expected for
>> every triangle to have (degi + 1) * (degi + 2) / 2 quadrature points.
>>
>> Surprisingly for me, the following illustration below shows that for
>> a degree 4 for instance, I obtain 12 quadrature points which is
>> different from 15.. It seems that some vertices of the triangle are
>> not taken into account.
>>
>> Is it because the weights of these points is 0.? Or perhaps I didn't
>> call the class in the right way?
>>
>> Any comments are welcome,
>>
>> best,
>>
>> Edouard.
>>
>> /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
>>
>>
>> import getfem
>> import numpy as np
>>
>> def nb_imdata(degi):
>>     # creation mesh / fem / im
>>     NX = 6
>>     tmesh = getfem.Mesh('regular_simplices', np.arange(0, 1 + 1. /
>> NX, 1. / NX),
>>                                              np.arange(0, 1 + 1. /
>> NX, 1. / NX))
>>     tmeshfem = getfem.MeshFem(tmesh, 1)
>>     tmeshfem.set_fem(getfem.Fem('FEM_PK(2, 2)'))
>>     tmeshim = getfem.MeshIm(tmesh, getfem.Integ('IM_NC(2,' +
>> str(degi) + ')'))
>>     #print tmeshim.im_nodes()
>>     print 'degree ', degi, '|', \
>>           'nb of quad points ', tmeshim.im_nodes(2)[0].shape[0], '|', \
>>           'expected          ', ((degi + 1) * (degi + 2) / 2)
>>
>>
>> for  k in range(3, 10):
>>     nb_imdata(k)
>>
>

-- 
Edouard Oudet : http://www-ljk.imag.fr/membres/Edouard.Oudet/
IMAG - Bureau 164
700 avenue Centrale
38400 Saint Martin d'Hères
+33 (0)4 57 42 17 71  (office LJK)
+33 (0)4 79 68 82 06  (home)