Re: [Getfem-users] gf_asm : How to use?

[Richard George]

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
>> address@hidden
>> https://mail.gna.org/listinfo/getfem-users
>>   
>
>

mesh_fem_test.pdf
Description: Adobe PDF document

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)'