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

[Kajetan Sikorski]

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.

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.

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