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
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
