Commit 2bae2bbb authored by MartinA's avatar MartinA
Browse files

some modifs

parent 55d35e73
function Ms = domRegularize2D(data)
%+========================================================================+
%| |
%| OPENDOM - LIBRARY FOR NUMERICAL INTEGRATION |
%| openDom is part of the GYPSILAB toolbox for Matlab |
%| |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018. |
%| PROPERTY : Centre de Mathematiques Appliquees, Ecole polytechnique, |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved. |
%| LICENCE : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use, |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or |
%| later, http://www.gnu.org/licenses). For private use, dual licencing |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT : matthieu.aussal@polytechnique.edu |
%| francois.alouges@polytechnique.edu |
%| WEBSITE : www.cmap.polytechnique.fr/~aussal/gypsilab      |
%| |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it. |
%|________________________________________________________________________|
%| '&` | |
%| # | FILE : domRegularize2D.m |
%| # | VERSION : 0.53 |
%| _#_ | AUTHOR(S) : Matthieu Aussal & François Alouges |
%| ( # ) | CREATION : 25.11.2018 |
%| / 0 \ | LAST MODIF : 14.03.2019 |
%| ( === ) | SYNOPSIS : Finite element regularization matrix for |
%| `---' | singularities with Laplace kernel |
%+========================================================================+
%%% INPUT ANALYSIS
if isa(data{2},'Wdom')
Ms = WdomRegularize2D(data);
return
end
if (length(data) == 4)
X = data{1};
Ydom = data{2};
green = data{3};
v = data{4};
else
Xdom = data{1};
Ydom = data{2};
u = data{3};
green = data{4};
v = data{5};
end
%%% INITIALIZATION
% Mesh data from Y
vtx = Ydom.msh.vtx;
elt = Ydom.msh.elt;
ctr = Ydom.msh.ctr;
nrm = Ydom.msh.nrm;
stp = Ydom.msh.stp;
lgt = Ydom.msh.ndv;
tau = Ydom.msh.tgt;
nu = cell2mat({tau,tau});
Nelt = size(elt,1);
% Quadrature data from Y
[Y,Wy,elt2qud] = Ydom.qud;
Yun = ones(1,size(elt2qud,2));
% Degrees of freedom from Y
[~,elt2dof] = v.dof;
Nbas = size(elt2dof,2);
% Quadrature data from X
if (length(data) == 5)
[X,Wx] = Xdom.qud;
end
Nx = size(X,1);
% Rangesearch with max(|edge|)_Y
[Ielt,Relt] = rangeSearch(X,ctr,1.1*stp(2)); %%% DEBUG %%%
Mx = cell(Nelt,1);
%%% RIGHT INTEGRATION WITH REGULARIZATION
for el = 1:Nelt
% Edge data for Y
Sel = vtx(elt(el,:),:);
Nel = nrm(el,:);
Tel = tau(el,:);
NUel = reshape(nu(el,:),3,2)';
% Local size
rMin = 1.1*lgt(el); %%% DEBUG %%%
% Quadratures points in interaction
Iy = elt2qud(el,:);
Ix = sort(Ielt{el}(Relt{el}<rMin))';
% % Graphical representation %%% DEBUG %%%
% figure(10)
% plot(Ydom.msh.sub(el))
% hold on
% plot3(X(Ix,1),X(Ix,2),X(Ix,3),'*')
% axis equal
% hold off
% pause()
% If interactions
if ~isempty(Ix)
%%% CORRECTION WITH SEMI-ANALYTIC INTEGRATION
% Analytical integration
[logR,rlogR,gradlo gR] = domSemiAnalyticInt2D(X(Ix,:),Sel,Nel,Tel);
% logR(:) = 0; rlogR(:) = 0; gradlogR(:) = 0; %%% DEBUG %%%
% Vector yg-x
Xun = ones(length(Ix),1);
XY1 = Xun * Y(Iy,1)' - X(Ix,1) * Yun;
XY2 = Xun * Y(Iy,2)' - X(Ix,2) * Yun;
% Distance r = |yg-x|
Rxy = sqrt(XY1.^2 + XY2.^2);
logRxy = log(Rxy);
logRxy(Rxy<1e-12) = log(1e-12);
% Int_el(log(|r|)) - Sum_g log(|yg-x|)
logR = logR - logRxy * Wy(Iy);
% norm(logR) %%% DEBUG %%%
% Int_el(r*log(|r|)) - Sum_g (yg-x)*log(|yg-x|)
rlogR(:,1) = rlogR(:,1) - (XY1 .* logRxy) * Wy(Iy);
rlogR(:,2) = rlogR(:,2) - (XY2 .* logRxy) * Wy(Iy);
rlogR(:,3) = 0;
% norm(rlogR) %%% DEBUG %%%
% Int_el(r/|r|^2) - Sum_g (yg-x)/|yg-x|^2
Rxym12 = 1./(Rxy.^2);
Rxym12(Rxy<1e-12) = 0;
gradlogR(:,1) = gradlogR(:,1) - (XY1 .* Rxym12) * Wy(Iy);
gradlogR(:,2) = gradlogR(:,2) - (XY2 .* Rxym12) * Wy(Iy);
gradlogR(:,3) = 0;
% norm(gradlogR) %%% DEBUG %%%
% Nullify V
V = [];
%%% FINITE ELEMENT P0
if strcmp(v.typ,'P0')
% Correction
if strcmp(green,'[log(r)]') && strcmp(v.opr,'[psi]')
V = logR;
elseif strcmp(green,'grady[log(r)]') && strcmp(v.opr,'n*[psi]')
V = gradlogR * Nel';
elseif strcmp(green,'[log(r)]') && strcmp(v.opr,'n*[psi]')
V{1} = logR .* Nel(1);
V{2} = logR .* Nel(2);
V{3} = logR .* Nel(3);
elseif strcmp(green,'[log(r)]') && strcmp(v.opr,'nxgrad[psi]')
V{1} = zeros(size(logR));
V{2} = zeros(size(logR));
V{3} = zeros(size(logR));
else
error('domRegularize2D.m : unavailable case')
end
%%% FINITE ELEMENT P1
elseif strcmp(v.typ,'P1')
% For each basis function
for j = 1:Nbas
% Next dof
jp1 = mod(j,2) + 1;
% Height from j
hj = ((Sel(j,:)-Sel(jp1,:)) * NUel(j,:)');
% Scalar product (x-yk).nuj/hj
tmp = ( (X(Ix,1)-Sel(jp1,1))*NUel(j,1) + ...
(X(Ix,2)-Sel(jp1,2))*NUel(j,2) ) ./ hj;
% Correction
if strcmp(green,'[log(r)]') && strcmp(v.opr,'[psi]')
V(:,j) = logR.*tmp + rlogR*NUel(j,:)'/hj;
elseif strcmp(green,'[log(r)]') && strcmp(v.opr,'n*[psi]')
Vx = logR.*tmp + rlogR*NUel(j,:)'/hj;
V{1}(:,j) = Vx .* Nel(1);
V{2}(:,j) = Vx .* Nel(2);
V{3}(:,j) = Vx .* Nel(3);
elseif strcmp(green,'[log(r)]') && strcmp(v.opr,'nxgrad[psi]')
NxNUj = cross(Nel,NUel(j,:));
V{1}(:,j) = NxNUj(1)/hj .* logR;
V{2}(:,j) = NxNUj(2)/hj .* logR;
V{3}(:,j) = NxNUj(3)/hj .* logR;
elseif strcmp(green,'grady[log(r)]') && strcmp(v.opr,'n*[psi]')
V(:,j) = tmp .* (gradlogR * Nel');
else
error('domRegularize2D.m : unavailable case')
end
end
else
error('domRegularize2D.m : unavailable case')
end
% Matrix-Vector product
I = repmat(Ix,1,Nbas);
J = repmat(elt2dof(el,:),length(Ix),1);
if iscell(V)
Mx{el} = [I(:) J(:) V{1}(:) V{2}(:) V{3}(:)];
else
Mx{el} = [I(:) J(:) V(:)];
end
end
end
%%% LEFT INTEGRATION AND RIGHT REDUCTION
% Left integration matrix
Ndof = size(v.dof,1);
if (length(data) == 4)
Mu = speye(Nx,Nx);
Mw = speye(Nx,Nx);
Ms = sparse(Nx,Ndof);
else
Mu = u.uqm(Xdom);
Mw = spdiags(Wx,0,length(Wx),length(Wx));
Ms = sparse(length(u),Ndof);
end
% Regularization matrix
Mx = double(cell2mat(Mx));
if isempty(Mx)
if iscell(Mu)
Mx = [1 1 zeros(1,length(Mu))];
else
Mx = [1 1 0];
end
end
% Left integration
if iscell(Mu)
for i = 1:length(Mu)
Ms = Ms + Mu{i}' * Mw * sparse(Mx(:,1),Mx(:,2),Mx(:,2+i),Nx,Ndof);
end
else
Ms = Mu' * Mw * sparse(Mx(:,1),Mx(:,2),Mx(:,3),Nx,Ndof);
end
% Right reduction
[~,Mv] = v.unk;
Ms = Ms * Mv;
end
......@@ -4,18 +4,29 @@ clear all;
%% Geometry
% Create a non radially symetric mesh mesh whose boundary has two connected components:
N = 50; rad1 = 1; rad2 = 2;
m_Space = mshAssymetricRing(N,rad1,rad2);
m_Gamma = union(mshSphere(100,1),translatation([2,2,2],mshSphere(100,1)));
m_Gamma = m_Space.bnd;
c = m_Gamma.ctr;
indGamma1 = c(:,1).^2 + c(:,2).^2 > (rad1 + rad2)^2/4;
indGamma2 = ~indGamma1;
m_Gamma1 = m_Gamma.sub(indGamma1);
m_Gamma2 = m_Gamma.sub(indGamma2);
% m_Omega = mshAssymetricRing(N,rad1,rad2);
% m_Gamma = m_Omega.bnd;
% c = m_Gamma.ctr;
% indGamma1 = c(:,1).^2 + c(:,2).^2 > (rad1 + rad2)^2/4;
% indGamma2 = ~indGamma1;
% m_Gamma1 = m_Gamma.sub(indGamma1);
% m_Gamma2 = m_Gamma.sub(indGamma2);
N = 800;
m_Omega1 = mshDisk(N,0.1);
m_Omega2 = translate(mshDisk(N,0.3),[-0.5,0,0]);
m_Omega = union(m_Omega1,m_Omega2);
m_Gamma = m_Omega.bnd;
[~,ind_Gamma1] = intersect(m_Gamma,bnd(m_Omega1));
[~,ind_Gamma2] = intersect(m_Gamma,bnd(m_Omega2));
m_Gamma1 = m_Gamma.sub(ind_Gamma1);
m_Gamma2 = m_Gamma.sub(ind_Gamma2);
m_Space = mshDisk(3*N,2);
% Domains of integration
gss = 3; % Number of Gauss points
Gamma = dom(m_Gamma,gssm_Gamma = m_Space.bnd;);
Gamma = dom(m_Gamma,gss);
Gamma2 = dom(m_Gamma2,gss);
Gxy = @(X,Y)femGreenKernel(X,Y,'[log(r)]',0); % 0 wave number
......@@ -26,27 +37,40 @@ S1_Gamma2 = fem(m_Gamma2,'P1');
V = -1/(2*pi)*integral(Gamma,Gamma,S1_Gamma,Gxy,S1_Gamma);
V = V + -1/(2*pi)*regularize(Gamma,Gamma,S1_Gamma,'[log(r)]',S1_Gamma);
B = integral(Gamma,S1_Gamma);
sys = [V B;B' 0];
B = integral(Gamma,S1_Gamma); % Array containing the mean value of each basis function
sys = [V B;B' 0]; % Solve with zero mean value constraing
% right hand side:
P = restriction(S1_Gamma,m_Gamma2);
F = integral(Gamma,S1_Gamma,S1_Gamma)*(P'*ones(size(P,1),1));
P = restriction(S1_Gamma,m_Gamma2); % Operator restrictinng dofs from S^1,0(Gamma) to S^{1,0}(Gamma2). P' extends by 0 from Gamma2 to Gamma.
g = P'*ones(size(P,1),1); % Function equal to 1 on Gamma2 and 0 on Gamma1
F = integral(Gamma,S1_Gamma,S1_Gamma)*g;
rhs = [F;0];
sol = sys\rhs;
lambda = sol(1:end-1);
mu = sol(end);
dnu = P*lambda;
lambda = sol(1:end-1); % Normal trace
mu = sol(end); % Lagrange multiplier.
fprintf('Lagrange Multiplier: %s\n',num2str(mu))
dnu = P*lambda; % Restriction to Gamma2
figure;
surf(S1_Gamma2,dnu);
plot(m_Omega,'w');
hold on
surf(S1_Gamma,lambda);
axis equal;
title('Charge distribution on the conductor Omega2')
SL = (-1/(2*pi))*integral(m_Gamma1.vtx,Gamma,Gxy,S1_Gamma);
SL = (-1/(2*pi))*integral(m_Space.vtx,Gamma,Gxy,S1_Gamma);
% Array of int_{Gamma}G(x - y)phi_h(y) for x = x1 ... xN
SL = SL + (-1/(2*pi))*regularize(m_Gamma1.vtx,Gamma,'[log(r)]',S1_Gamma);
SL = SL + (-1/(2*pi))*regularize(m_Space.vtx,Gamma,'[log(r)]',S1_Gamma);
u = SL*lambda+mu;
u = SL*lambda+mu;
S1_Omega = fem(m_Space,'P1');
figure
graph(S1_Omega,u);
title('Electrostatic potential')
I = integral(Gamma2,S1_Gamma2,ntimes(S1_Gamma2));
F = 1/2*dnu'*I{1}*dnu;
\ No newline at end of file
F = 1/2*dnu'*I{1}*dnu;
fprintf('Force (in the x direction) %s\n',num2str(F))
\ No newline at end of file
......@@ -4,13 +4,22 @@ clear all;
%% Geometry
% Create a non radially symetric mesh mesh whose boundary has two connected components:
N = 50; rad1 = 1; rad2 = 2;
m_Space = mshAssymetricRing(N,rad1,rad2);
m_Gamma = m_Space.bnd;
c = m_Gamma.ctr;
indGamma1 = c(:,1).^2 + c(:,2).^2 > (rad1 + rad2)^2/4;
indGamma2 = ~indGamma1;
m_Gamma1 = m_Gamma.sub(indGamma1);
m_Gamma2 = m_Gamma.sub(indGamma2);
% m_Omega = mshAssymetricRing(N,rad1,rad2);
% m_Gamma = m_Omega.bnd;
% c = m_Gamma.ctr;
% indGamma1 = c(:,1).^2 + c(:,2).^2 > (rad1 + rad2)^2/4;
% indGamma2 = ~indGamma1;
% m_Gamma1 = m_Gamma.sub(indGamma1);
% m_Gamma2 = m_Gamma.sub(indGamma2);
m_Omega1 = mshDisk(100,1);
m_Omega2 = translate(mshSquare(100,[1,1]),[-3,0,0]);
m_Omega = union(m_Omega1,m_Omega2);
m_Gamma = m_Omega.bnd;
[~,ind_Gamma1] = intersect(m_Gamma,bnd(m_Omega1));
[~,ind_Gamma2] = intersect(m_Gamma,bnd(m_Omega2));
m_Gamma1 = m_Gamma.sub(
% Domains of integration
gss = 3; % Number of Gauss points
......@@ -18,34 +27,47 @@ Gamma = dom(m_Gamma,gss);
Gamma2 = dom(m_Gamma2,gss);
Gxy = @(X,Y)femGreenKernel(X,Y,'[log(r)]',0); % 0 wave number
Gxy = @(X,Y)(X(:,1))
S1_Gamma = fem(m_Gamma,'P1');
S1_Gamma2 = fem(m_Gamma2,'P1');
% Single layer potential
V = -1/(2*pi)*integral(Gamma,Gamma,S1_Gamma,Gxy,S1_Gamma);
% V = V + -1/(2*pi)*regularize(Gamma,Gamma,S1_Gamma,'[log(r)]',S1_Gamma);
V = V + -1/(2*pi)*regularize(Gamma,Gamma,S1_Gamma,'[log(r)]',S1_Gamma);
B = integral(Gamma,S1_Gamma);
sys = [V B;B' 0];
B = integral(Gamma,S1_Gamma); % Array containing the mean value of each basis function
sys = [V B;B' 0]; % Solve with zero mean value constraing
% right hand side:
P = restriction(S1_Gamma,m_Gamma2);
F = integral(Gamma,S1_Gamma,S1_Gamma)*(P'*ones(size(P,1),1));
P = restriction(S1_Gamma,m_Gamma2); % Operator restrictinng dofs from S^1,0(Gamma) to S^{1,0}(Gamma2). P' extends by 0 from Gamma2 to Gamma.
g = P'*ones(size(P,1),1); % Function equal to 1 on Gamma2 and 0 on Gamma1
F = integral(Gamma,S1_Gamma,S1_Gamma)*g;
rhs = [F;0];
sol = sys\rhs;
lambda = sol(1:end-1);
mu = sol(end);
dnu = P*lambda;
lambda = sol(1:end-1); % Normal trace
mu = sol(end); % Lagrange multiplier.
fprintf('Lagrange Multiplier: %s\n',num2str(mu))
dnu = P*lambda; % Restriction to Gamma2
figure;
surf(S1_Gamma2,dnu);
plot(m_Omega,'w');
hold on
surf(S1_Gamma,lambda);
axis equal;
title('Charge distribution on the conductor Omega2')
SL = (-1/(2*pi))*integral(m_Omega.vtx,Gamma,Gxy,S1_Gamma);
% Array of int_{Gamma}G(x - y)phi_h(y) for x = x1 ... xN
SL = SL + (-1/(2*pi))*regularize(m_Omega.vtx,Gamma,'[log(r)]',S1_Gamma);
SL = (-1/(2*pi))*integral(m_Gamma1.vtx,Gamma,Gxy,S1_Gamma);
SL = SL + (-1/(2*pi))*regularize(m_Gamma1.vtx,Gamma,'[log(r)]',S1_Gamma);
u = SL*lambda+mu;
u = SL*lambda+mu;
S1_Omega = fem(m_Omega,'P1');
figure
graph(S1_Omega,u);
title('Electrostatic potential')
I = integral(Gamma2,S1_Gamma2,ntimes(S1_Gamma2));
F = 1/2*dnu'*I{1}*dnu;
\ No newline at end of file
F = 1/2*dnu'*I{1}*dnu;
fprintf('Force (in the x direction) %s\n',num2str(F))
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment