Re: [Getfem-users] periodic boundary conditions example and source term question

[Yves Renard]

On lundi 5 avril 2010, Kajetan Sikorski wrote:
> Hello,
>
> I'm new to getfem, but really like it so far. Great Job!
>
> I was browsing the old threads and could not find much information on
> how to setup periodic boundary conditions. I figured it out eventually
> so I'm attaching a cooked up matlab example which might hopefully help
> someone in the future.

Thank you.

>
> I also have a question. In the stokes problem implemented in the code
> below I would like to have an additional div(S) forcing term where S
> is a second order extra stress tensor. In the weak formulation this
> should show up as \int_{\omega} S: grad(v) dx where : is the double
> contraction product of two tensors. How could I do this? Thanks.


There is no specific brick for the addition of such a source term.
You can make the assembly of this term using a generic assembly which should 
be something like

V = gf_asm('volumic' , 'S=data(qdim(#1),#1);V(#1)+=comp(vGrad(#1))
(i,j,i,j).S(i,j)"", mim, mf, S);

and then use 
gf_model_set(model M, 'add explicit rhs', string varname, vec L)
to explicitely add your term to the model.

Yves.

>
> Kai
>
> function stokes
> %compute solution of stokes problem on semi-periodic domain
>
> M = 30;                                                 %grid resolution
> N = 30;
> left = 0;                                               %grid boundaries
> right = 1;
> top = 1;
> bottom = 0;
> gf_workspace('clear all');
>
> m = gf_mesh('cartesian',linspace(left,right,N),linspace(bottom,top,M));
> femu = gf_mesh_fem(m,2);
> femp = gf_mesh_fem(m,1);
> femd = gf_mesh_fem(m,1);
>
> gf_mesh_fem_set(femu,'fem',gf_fem('FEM_QK(2,2)'));          %assign FEs
> gf_mesh_fem_set(femp,'fem',gf_fem('FEM_QK(2,1)'));
> gf_mesh_fem_set(femd,'fem',gf_fem('FEM_QK(2,2)'));
>                                  %assign integration type
> mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,10)'));
>
> md = gf_model('real');          %this holdes the model
>
> gf_model_set(md, 'add fem variable', 'u',femu);
> gf_model_set(md, 'add fem variable', 'p',femp);
> gf_model_set(md, 'add linear incompressibility brick', mim, 'u', 'p');
> gf_model_set(md, 'add Laplacian brick', mim, 'u');
> gf_model_set(md, 'add variable', 'mult_spec', 1);
> gf_model_set(md, 'add constraint with multipliers', 'p',
> 'mult_spec', ...
>                 sparse(ones(1, gf_mesh_fem_get(femp, 'nbdof'))), 0);
> P=gf_mesh_get(m,'pts');
>                                                          %find
> boundaries
> pidtop=find(abs(P(2,:)-top)<1e-6);
> pidbottom = find(abs(P(2,:)-bottom)<1e-6);
> pidleft = find(abs(P(1,:)-left)<1e-6);
> pidright = find(abs(P(1,:)-right)<1e-6);
> ftop=gf_mesh_get(m,'faces from pid',pidtop);
> fbot=gf_mesh_get(m,'faces from pid',pidbottom);
> fleft=gf_mesh_get(m,'faces from pid',pidleft);
> fright=gf_mesh_get(m,'faces from pid',pidright);
> gf_mesh_set(m,'boundary',1,[ftop,fbot]);
> gf_mesh_set(m,'boundary',2,fleft);
> gf_mesh_set(m,'boundary',3,fright);
>
>                                  %apply periodic conditions on left
> and right
> leftDof = gf_mesh_fem_get(femu, 'basic dof on region', 2);
> rightDof = gf_mesh_fem_get(femu, 'basic dof on region',3);
> ConstraintMatrix = sparse(zeros(1,gf_mesh_fem_get(femu, 'nbdof')));
> for i=3:length(leftDof)-2,
>      ConstraintMatrix(1,leftDof(i))=1;
>      ConstraintMatrix(1,rightDof(i))=-1;
>      gf_model_set(md, 'add variable', strcat('mult_spec',num2str(i)),
> 1);
>      gf_model_set(md, 'add constraint with multipliers', 'u',
> strcat('mult_spec',num2str(i)), ...
>                      ConstraintMatrix(1,:), 0);
> end
>
>                                  %apply no-slip at top and bottom and
>
> F = gf_mesh_fem_get(femd, 'eval', {'(-1).*sin (2.*pi.*x).*sin
> (4.*pi.*y)';'cos (2.*pi.*x).*sin (2.*pi.*y).^2'});
> F2 = gf_mesh_fem_get(femd, 'eval',
> {'(-20
> ).*pi
> .^
> 2
> .*sin
> (2
> .*pi
> .*x
> ).*sin
> (4
> .*pi
> .*y
> )';'(-8
> ).*pi
> .^
> 2
> .*cos
> (2
> .*pi
> .*x).*cos(2.*pi.*y).^2+12.*pi.^2.*cos(2.*pi.*x).*sin(2.*pi.*y).^2'});
> gf_model_set(md, 'add initialized fem data', 'Dirdata', femd, ...
>      gf_mesh_fem_get(femd,'eval',{'0';'0'}));
> gf_model_set(md, 'add initialized fem data', 'Dirdata2',femd,...
>      F);
> gf_model_set(md, 'add initialized fem data', 'VolumicData',femd,...
>      F2);
>
> gf_model_set(md, 'add Dirichlet condition with multipliers', ...
>      mim, 'u', femu, 1, 'Dirdata2');
>
>                                  %add forcing
>
> gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
>
> gf_model_get(md, 'solve');
> P = gf_model_get(md, 'variable', 'p');
> U = gf_model_get(md, 'variable', 'u');
> gf_plot(femp,P,'mesh','on'); hold on
> gf_plot(femu,U,'mesh','on');
> colorbar; title('computed solution'); hold off



-- 

  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

---------