Hi Richard,
Yes you are right, there was a bug in the face number bounds for
gf_mesh_im_get(mim,gf_eltm('normal'),cv,face)
I have fixed it in our svn repository , however I cannot recommend that
you use it right now because it looks like we have a (unrelated) bug in
the node searching code, which is quite annoying.
Best regards,
Julien
Richard George wrote:
> Dear Julien
>
> Thanks for your helpful answer to my gf_asm question, the solution you
> provided works with 2D meshes.
> for some reason I get a bad answer using the same command in a 3D case -
>
> I think I have found a bug with the 'normal' function, or I
> misunderstand it. Is normal() defined for 3D meshes?
>
> Here is an example that makes one convex in a 3D mesh, and tries to
> print the normals on the faces.
> I find two gf_ functions that behave as if the same convex has different
> numbers of faces?
>
> gf_mesh_get(m,'normal of face',cv,face) accepts cv=1, face =1,2,3,4
> and gives the unit normals on the four faces correctly.
>
> I tried:
>
> gf_mesh_im_get(mim,gf_eltm('normal'),cv,face)
>
> to provide the integral of the normal on the face, ie. the normal is a
> constant vector
> and the faces have a known area, so I'd expect
>
> Area_of_face * gf_mesh_get(m,'normal of face',cv,face) =
> gf_mesh_im_get(mim,gf_eltm('normal'),cv,face)
>
> This holds for face=1,2,3, but for the 3D mesh, gf_mesh_im_get does not
> accept face=4, is it meant to?
> please see attached examples.
>
>
> Thanks for any help you can offer
>
>
> Yours
>
> Richard George
>
>
>
>
> julien pommier wrote:
>
>> Hi Richard,
>>
>> Your expression for I1 is the correct one. The first index of a
>> Base(), vBase(), Grad etc is always the one refering to the
>> corresponding degree of freedom
>>
>> a=data(#1) and a=data$1(#1) are equivalent. The '$' is only useful
>> when you have more than one data argument (for example b=data$2(#1))
>> The #1 says that the corresponding dimension is the number of degrees
>> of freedom of the mesh_fem number 1
>>
>> if U lives on a mesh_fem mf, and a lives on a mesh_fem mfd (a scalar
>> mesh_fem whose qdim is 1), then the expression for I2 is:
>>
>> I2 = gf_asm('boundary',1,'u=data$1(#1);a=data$2(#2);
>> V()+=a(i).u(j).comp(vBase(#1).Base(#2).Normal())(j,k,i,k);',mim,mf,mfd,U,A);
>>
>>
>> In order to understand that, just write
>> a(x) = sum(a_i * Phi_i(x)) with Phi_i(x) being the scalar base
>> functions of mesh_fem mf1
>> U(x) = sum(u_j * Psi_j(x)) with Psi_j(x) being the vector base
>> functions of mf (so Psi_j(x)[k] is one of its components).
>>
>> Everything that is inside the 'comp' is in the integral, so you have
>>
>> sum_{i,j,k} a_i * U_j * integral(Phi_i(x) * Psi_j(x)[k] * Normal[k] dS)
>>
>> I hope it is more clear now !
>>
>> Best regards,
>> Julien
>>
>> Richard George wrote:
>>
>>> Hello
>>>
>>> I'd like to evaluate the integral of the normal component of a vector
>>> valued mesh_fem on a boundary,
>>>
>>> I have 'U' as a vector valued function represented by a FEM_PK(3,1)
>>> object, and 'a' being a scalar valued function
>>> that takes a constant value on each convex, represented by a
>>> FEM_DISCONTINUOUS_PK(3,0) object.
>>>
>>> term
>>>
>>> term 2
>>>
>>> The code for evaluating the first integral via gf_asm is possible I
>>> think by making a contraction of a vBase() with a Normal()
>>>
>>> I1 =
>>> gf_asm('boundary',1,'u=data(#1);V()+=u(i).comp(vBase(#1).Normal())(i,j,j);',mim,mf,U);
>>>
>>>
>>> This appears to give the right results in some simple tests - I'm
>>> assuming that i sums over nodes, while j,j sums over vector
>>> components and provides
>>> a dot product - but I don't really understand how are the indexes on
>>> the comp() function are determined ? When do the indexes all refer to
>>> cartesian
>>> vector components, and when are they local node numbers? am I using
>>> the Normal() option correctly?
>>>
>>> I think it should be possible to specify the second integral as a
>>> contraction of 'a', 'U' and a tensor but
>>> I don't think I grasp the syntax of the gf_asm command properly.
>>>
>>> Could you explain how to specify integral I2 via gf_asm, and when
>>> it's appropriate to use
>>>
>>> a=data(#1)
>>> a=data$1(#1)
>>> a=data(#1,qdim(#1))
>>>
>>> I2 = gf_asm('boundary',1,'u=data(#1);a=data(#2)
>>> ;V()+=a(i).u(j,k).comp(---??---.vBase(#1).Normal())(i,j,k,l,l);',mim,mf1,mf0,U,A);
>>>
>>>
>>> Thanks for any help you can provide
>>>
>>> Yours
>>>
>>> Richard George
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> ------------------------------------------------------------------------
>>>
>>> _______________________________________________
>>> Getfem-users mailing list
>>> [email protected]
>>> https://mail.gna.org/listinfo/getfem-users
>>>
>>>
>>
>
>
> ------------------------------------------------------------------------
>
> trace on;
> %% Make a mesh with a single convex
> m2=gf_mesh('empty',3);
> %% Add four coordinates
> pt1=gf_mesh_set(m2,'add point',[1;0;0]);
> pt2=gf_mesh_set(m2,'add point',[0;1;0]);
> pt3=gf_mesh_set(m2,'add point',[0;0;1]);
> pt4=gf_mesh_set(m2,'add point',[0;0;0]);
> %% Add a single tetrahedron
> gt=gf_geotrans('GT_PK(3,1)');
> pts=gf_mesh_get(m2,'pts',[pt1,pt2,pt3,pt4]);
> cv1=gf_mesh_set(m2,'add convex',gt,pts);
> gf_plot_mesh(m2);
> hold on;
> view(45,20);
> axis equal;
> %% Make a mesh fem
> mf2=gf_mesh_fem(m2,1);
> gf_mesh_fem_set(mf2,'fem',gf_fem('FEM_PK_DISCONTINUOUS(3,0)'));
> %% Provide an integration method
> mim=gf_mesh_im(m2,gf_integ('IM_NC(3,3)'));
> %% Normals via gf_mesh_get
> %
> disp('Normal of face 1 is:');
> n1=gf_mesh_get(m2,'normal of face',1,1)
> disp('Normal of face 2 is:');
> n2=gf_mesh_get(m2,'normal of face',1,2)
> disp('Normal of face 3 is:');
> n3=gf_mesh_get(m2,'normal of face',1,3)
> disp('Normal of face 4 is:');
> n4=gf_mesh_get(m2,'normal of face',1,4)
>
> %% Normals via gf_mesh_im_get
> %
> % This should also list integrals which are proportional to the normals,
> % (area of each face is 1/2)
> %
> gf=gf_eltm('normal');
>
> disp('normal 1');
> i1=gf_mesh_im_get(mim,'eltm',gf,1,1)'
> disp('normal 2');
> i2=gf_mesh_im_get(mim,'eltm',gf,1,2)'
> disp('normal 3');
> i3=gf_mesh_im_get(mim,'eltm',gf,1,3)'
>
> %% Try the fourth face
> disp('normal 4');
> i4=gf_mesh_im_get(mim,'eltm',gf,1,4)'
>
--
Julien Pommier, bureau 111
GMM, INSA Toulouse, tél:05 61 55 93 42
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
