Commit 17a03c36 authored by MartinA's avatar MartinA
Browse files

Commuting Diagrams and NetForce test case

parent ef5d2dd2
classdef Florian
%UNTITLED2 Summary of this class goes here
% Detailed explanation goes here
properties
adress
tel
end
methods
function obj = getAdress(inputArg1,inputArg2)
%UNTITLED2 Construct an instance of this class
% Detailed explanation goes here
obj.Property1 = inputArg1 + inputArg2;
end
function outputArg = getTel(obj,inputArg)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
outputArg = obj.Property1 + inputArg;
end
end
end
% S = singleLayer(m);
clear all;
close all;
ps = 2:10;
Feng = 0*ps;
Fvol = 0*ps;
FBEM = 0*ps;
errEng = 0*ps;
errVol = 0*ps;
errBEM = 0*ps;
Ns = 0*ps;
Nvol = 0*ps;
Nbound = 0*ps;
refVal = 0.348731673835553;
% refVal = computeNetForce2D(2^(ps(end) +2));
for i = 1:length(ps)
N = 2^ps(i);
[FBEM(i),Fvol(i),Feng(i),Nvol(i),Nbound(i),h(i)] = computeNetForce2D(N);
errEng(i) = abs(Feng(i) - refVal);
errVol(i) = abs(Fvol(i) - refVal);
errBEM(i) = abs(FBEM(i) - refVal);
end
loglog(h,errEng);
hold on
loglog(h,errVol);
loglog(h,errBEM);
loglog(h,h*errEng(1)/h(1),'b--');
loglog(h,h.^2*errVol(1)/h(1)^2,'r--');
legend({'Surface formula','Volume formula','BEM formula',...
'O(h)','O(h^2)'})
xlabel('Mesh size')
ylabel('|Approx - Ref|')
function [FBEM,Fvol,Feng,NvtxVol,NvtxBound,h] = computeNetForce2D(N)
%% Geometry
% Create a non radially symetric mesh mesh whose boundary has two connected components:
rad1 = 1; rad2 = 2;
m_Omega = mshAssymetricRing(N,rad1,rad2);
m_Gamma = m_Omega.bnd; % boundary of the mesh
c = m_Gamma.ctr;
indGamma1 = c(:,1).^2 + c(:,2).^2 > (rad1 + rad2)^2/4;
m_Gamma1 = m_Gamma.sub(indGamma1);
m_Gamma2 = setdiff(m_Gamma,m_Gamma1);
% Domains of integration
gss = 3; % Number of Gauss points
Omega = dom(m_Omega,gss);
%% Variational setting
% Finite element spaces
S1_Omega = fem(m_Omega,'P1'); % Piecewise linear functions on Omega.
S10_Omega= dirichlet(S1_Omega,m_Gamma); % Pw linear Dirichlet condition
% Construct a function G on Omega which is equal to 1 on Gamma1 and 0 on
% Gamma2.
P= restriction(S1_Omega,m_Gamma2); % operator restriction on Gamma2
Q = P'; % Adjoint operator: extension by 0.
Gh= Q*ones(m_Gamma2.nvtx,1);
%% Computation of the electrostatic potential
% Linear form
rhs = -integral(Omega,grad(S10_Omega),grad(S1_Omega))*Gh;
% Bilinear form:
a = integral(Omega,grad(S10_Omega),grad(S10_Omega));
% Solution:
wh = a\rhs;
P = elimination(S1_Omega,m_Gamma);
Uh = P*wh + Gh;
%% Computation of the Electric field
S0_Omega = fem(m_Omega,'P0');
M0 = integral(Omega,S0_Omega,S0_Omega);
Omega = dom(m_Omega,gss); % We need just one Gauss point since everything is constant
E0 = cell(3,1);
E0{1} = M0\(integral(Omega,S0_Omega,grad(S1_Omega,1))*Uh);
E0{2} = M0\(integral(Omega,S0_Omega,grad(S1_Omega,2))*Uh);
E0{3} = M0\(integral(Omega,S0_Omega,grad(S1_Omega,3))*Uh);
%% Computation of the normal trace
gss = 1;
S0_Omega = fem(m_Omega,'P0');
M0 = integral(Omega,S0_Omega,S0_Omega);
Omega = dom(m_Omega,gss); % We need just one Gauss point since everything is constant
S0_Gamma2 = fem(m_Gamma2,'P0');
% Compute the gradient of uh at the center of each element, and restrict to
% Gamma2
[tr_Gamma_P0,S0_Gamma] = trace_P0P1(S0_Omega);
P = restriction(S0_Gamma,m_Gamma2);
tr_Gamma2_P0 = P*tr_Gamma_P0;
gradUh1 = tr_Gamma2_P0*(M0\integral(Omega,S0_Omega,grad(S1_Omega,1))*Uh);
gradUh2 = tr_Gamma2_P0*(M0\integral(Omega,S0_Omega,grad(S1_Omega,2))*Uh);
gradUh3 = tr_Gamma2_P0*(M0\integral(Omega,S0_Omega,grad(S1_Omega,3))*Uh);
% Scalar product by normals:
Nrm = m_Gamma2.nrm;
dnU = Nrm(:,1).*gradUh1 + Nrm(:,2).*gradUh2 + Nrm(:,3).*gradUh3;
%% Net force
Gamma2 = dom(m_Gamma2,1);
I = integral(Gamma2,S0_Gamma2,ntimes(S0_Gamma2));
Feng = 1/2*dnU'*I{1}*dnU;
E2 = E0{1}.^2 + E0{2}.^2 + E0{3}.^2;
Fvol = 1/2*sum(m_Omega.ndv.*E2.*E0{1});
NvtxVol= m_Omega.nvtx;
%% BEM approach
%% Geometry
% Create a non radially symetric mesh mesh whose boundary has two connected components:
% Domains of integration
gss = 3; % Number of Gauss points
Gamma = dom(m_Gamma,gss);
Gxy = @(X,Y)femGreenKernel(X,Y,'[log(r)]',0); % 0 wave number
S1_Gamma = fem(m_Gamma,'P0');
S1_Gamma2 = fem(m_Gamma2,'P0');
% 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);
B = integral(Gamma,S1_Gamma);
sys = [V B;B' 0];
% right hand side:
P = restriction(S1_Gamma,m_Gamma2);
F = integral(Gamma,S1_Gamma,S1_Gamma)*(P'*ones(size(P,1),1));
rhs = [F;0];
sol = sys\rhs;
lambda = sol(1:end-1);
dnu = P*lambda;
I = integral(Gamma2,S1_Gamma2,ntimes(S1_Gamma2));
FBEM = 1/2*dnu'*I{1}*dnu;
NvtxBound = m_Gamma.nvtx;
h = m_Omega.stp;
h = h(3);
end
function mesh = mshAssymetricRing(N,rad1,rad2)
rad = rad2-rad1;
% Radial discretisation
dr = sqrt(pi*rad^2/N);
dr = rad/ceil(rad/dr);
r = rad1:dr:rad2;
% Angular uniform discretization
rho = cell(length(r),1); theta = rho;
for ir = 1:length(r)
alpha = (r(ir) -rad1/2)/(rad2-rad1/2);
dtheta = dr/r(ir);
dtheta = 2*pi/ceil(2*pi/dtheta);
theta{ir} = (0:dtheta:2*pi-dtheta)';
rho{ir} = r(ir).*(1 + alpha*(1 + 0.5*cos(theta{ir})));
end
% Carthesian coordinates
[x,y] = pol2cart(cell2mat(theta),cell2mat(rho));
X = [0,0; x y];
% Unicity test
tmp = unique(X,'rows','stable');
if (max(abs(X-tmp)) > 1e-12)
error('mshDisk : non unicity of vertices')
end
% Delaunay triangulation
DT = delaunayTriangulation(X(:,1),X(:,2));
% Final mesh
elt = DT.ConnectivityList;
vtx = [DT.Points,zeros(size(DT.Points,1),1)];
mesh = msh(vtx,elt);
I = find(mesh.vtx(:,1).^2 + mesh.vtx(:,2).^2 < rad1.^2/2);
ind = and(and(~ismember(mesh.elt(:,1),I),~ismember(mesh.elt(:,2),I)),...
~ismember(mesh.elt(:,3),I));
mesh = mesh.sub(ind);
end
function mesh = mshBndAssymetricRing(N,rad1,rad2)
rad = rad2-rad1;
% Radial discretisation
dr = sqrt(pi*rad^2/N);
dr = rad/ceil(rad/dr);
r = rad1:dr:rad2;
% Angular uniform discretization
rho = cell(length(r),1); theta = rho;
alpha = rad1/2/(rad2-rad1/2);
dtheta = dr/rad1;
dtheta = 2*pi/ceil(2*pi/dtheta);
theta{1} = (0:dtheta:2*pi-dtheta)';
rho{1} = rad1.*(1 + alpha*(1 + 0.5*cos(theta{1})));
alpha = rad2/2/(rad2-rad1/2);
dtheta = dr/rad2;
dtheta = 2*pi/ceil(2*pi/dtheta);
theta{2} = (0:dtheta:2*pi-dtheta)';
rho{2} = rad2.*(1 + alpha*(1 + 0.5*cos(theta{2})));
% Carthesian coordinates
[x,y] = pol2cart(cell2mat(theta),cell2mat(rho));
X = [0,0; x y];
% Unicity test
tmp = unique(X,'rows','stable');
if (max(abs(X-tmp)) > 1e-12)
error('mshDisk : non unicity of vertices')
end
% Delaunay triangulation
DT = delaunayTriangulation(X(:,1),X(:,2));
% Final mesh
elt = DT.ConnectivityList;
vtx = [DT.Points,zeros(size(DT.Points,1),1)];
mesh = msh(vtx,elt);
I = find(mesh.vtx(:,1).^2 + mesh.vtx(:,2).^2 < rad1.^2/2);
ind = and(and(~ismember(mesh.elt(:,1),I),~ismember(mesh.elt(:,2),I)),...
~ismember(mesh.elt(:,3),I));
mesh = mesh.sub(ind);
mesh = mesh.bnd;
end
close all;
clear all;
%% Geometry
% Create a non radially symetric mesh mesh whose boundary has two connected components:
N = 30; rad1 = 1; rad2 = 2;
m_Omega = mshAssymetricRing(N,rad1,rad2);
m_Gamma = m_Omega.bnd; % boundary of the mesh
c = m_Gamma.ctr;
ind = c(:,1).^2 + c(:,2).^2 > (rad1 + rad2)^2/4;
m_Gamma1 = m_Gamma.sub(ind);
m_Gamma2 = setdiff(m_Gamma,m_Gamma1);
figure;
plot(m_Omega);
hold on;
plot(m_Gamma1,'g');
plot(m_Gamma2,'r');
axis off;
title('Mesh and components of the boundary');
% Domains of integration
gss = 3; % Number of Gauss points
Omega = dom(m_Omega,gss);
%% Variational setting
% Finite element spaces
S1_Omega = fem(m_Omega,'P1'); % Piecewise linear functions on Omega.
S10_Omega= dirichlet(S1_Omega,m_Gamma); % Pw linear Dirichlet condition
% Construct a function G on Omega which is equal to 1 on Gamma1 and 0 on
% Gamma2.
P= restriction(S1_Omega,m_Gamma2); % operator restriction on Gamma2
Q = P'; % Adjoint operator: extension by 0.
Gh= Q*ones(m_Gamma2.nvtx,1);
figure
graph(S1_Omega,Gh)
view(45,10)
axis equal
axis off
title('Function Gh')
%% Computation of the electrostatic potential
% Linear form
rhs = -integral(Omega,grad(S10_Omega),grad(S1_Omega))*Gh;
% Bilinear form:
a = integral(Omega,grad(S10_Omega),grad(S10_Omega));
% Solution:
wh = a\rhs;
P = elimination(S1_Omega,m_Gamma);
Uh = P*wh + Gh;
figure
graph(S1_Omega,Uh);
view(45,10)
axis equal
axis off
title('Electrostatic potential')
axis tight;
%% Computation of the Electric field
S0_Omega = fem(m_Omega,'P0');
M0 = integral(Omega,S0_Omega,S0_Omega);
Omega = dom(m_Omega,gss); % We need just one Gauss point since everything is constant
E0 = cell(3,1);
E0{1} = M0\(integral(Omega,S0_Omega,grad(S1_Omega,1))*Uh);
% int_{Omega} psi_h \partial_x(phi_h)
E0{2} = M0\(integral(Omega,S0_Omega,grad(S1_Omega,2))*Uh);
E0{3} = M0\(integral(Omega,S0_Omega,grad(S1_Omega,3))*Uh);
figure
title('Electric field')
plot(m_Omega,'w');
hold on
c = m_Omega.ctr;
quiver(c(:,1),c(:,2),E0{1},E0{2},'r')
%% Computation of the normal trace
S0_Gamma2 = fem(m_Gamma2,'P0');
gss = 1;
% Compute the gradient of uh at the center of each element, and restrict to
% Gamma2
[tr_Gamma_P0,S0_Gamma] = trace_P0P1(S0_Omega);
P = restriction(S0_Gamma,m_Gamma2);
tr_Gamma2_P0 = P*tr_Gamma_P0;
gradUh1 = tr_Gamma2_P0*(E0{1});
gradUh2 = tr_Gamma2_P0*(E0{2});
gradUh3 = tr_Gamma2_P0*(E0{3});
% Scalar product by normals:
Nrm = m_Gamma2.nrm;
dnU = Nrm(:,1).*gradUh1 + Nrm(:,2).*gradUh2 + Nrm(:,3).*gradUh3;
figure
plot(m_Omega,'w');
hold on
surf(S0_Gamma2,dnU)
%% Net force by engineering formula:
Gamma2 = dom(m_Gamma2,1);
I = integral(Gamma2,S0_Gamma2,ntimes(S0_Gamma2));
Fengineers = cell(3,1);
for i = 1:3
Fengineers{i} = 1/2*dnU'*I{i}*dnU;
end
disp(Fengineers{1}); % By symmetry, F_y = 0.
%% Net force by volume formula:
E2 = E0{1}.^2 + E0{2}.^2 + E0{3}.^2;
E2xE = cell(3,1);
Fvol = cell(3,1);
for i = 1:3
E2xE{i} = E2.*E0{i};
Fvol{i} = 1/2*sum(m_Omega.ndv.*E2xE{i});
end
close all;
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);
% Domains of integration
gss = 3; % Number of Gauss points
Gamma = dom(m_Gamma,gssm_Gamma = m_Space.bnd;);
Gamma2 = dom(m_Gamma2,gss);
Gxy = @(X,Y)femGreenKernel(X,Y,'[log(r)]',0); % 0 wave number
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);
B = integral(Gamma,S1_Gamma);
sys = [V B;B' 0];
% right hand side:
P = restriction(S1_Gamma,m_Gamma2);
F = integral(Gamma,S1_Gamma,S1_Gamma)*(P'*ones(size(P,1),1));
rhs = [F;0];
sol = sys\rhs;
lambda = sol(1:end-1);
mu = sol(end);
dnu = P*lambda;
figure;
surf(S1_Gamma2,dnu);
axis equal;
SL = (-1/(2*pi))*integral(m_Gamma1.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);
u = SL*lambda+mu;
I = integral(Gamma2,S1_Gamma2,ntimes(S1_Gamma2));
F = 1/2*dnu'*I{1}*dnu;
\ No newline at end of file
close all;
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);
% Domains of integration
gss = 3; % Number of Gauss points
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);
B = integral(Gamma,S1_Gamma);
sys = [V B;B' 0];
% right hand side:
P = restriction(S1_Gamma,m_Gamma2);
F = integral(Gamma,S1_Gamma,S1_Gamma)*(P'*ones(size(P,1),1));
rhs = [F;0];
sol = sys\rhs;
lambda = sol(1:end-1);
mu = sol(end);
dnu = P*lambda;
figure;
surf(S1_Gamma2,dnu);
axis equal;
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;
I = integral(Gamma2,S1_Gamma2,ntimes(S1_Gamma2));
F = 1/2*dnu'*I{1}*dnu;
\ No newline at end of file
%% Commuting diagrams:
% The following diagram is commuting:
% S^{1,0}(Omega) -> Ned(Omega) -> RT(Omega) -> S^{0,-1}(Omega)
% grad curl div
% | | |
% |gamma 1 | piT 2 | pi_n
% | | |
% V V V
% S^{1,0}(Gamma) -> Ned(Gamma) -> S^{0,-1}(Gamma)
% | grad | div(nx) |
% |Id 3 | nx 4 | Id
% V V V
% S^{1,0}(Gamma) -> RT(Gamma) -> S^{0,-1}(Gamma)
% nxgrad div
% In the sense that A o B = C o D whenever A and C are two arrows starting
% from the same space and B and D are two arrows pointing to the same
% space.
% When the topology of Omega is trivial, the horizontal seuquences are
% exact. A sequence
% A -> B -> C
% f g
% is exact if, for b in B, g(b) = 0 <=> b = f(a) for some a in A.
% Diagrams 3 and 4 are obvious!
%% Diagram 1:
m = mshCube(30,[1 1 1]);
dm = m.bnd;
Vh = fem(m,'P1');
Wh = fem(dm,'P1');
Xh = fem(m,'NED');
gamma = trace_P0P1(Vh);
piT = tangentialTrace_NED(Xh);
grad3D = grad_P1(Vh);
grad2D = grad_P1(Wh);
Comm = piT*grad3D - grad2D*gamma;
disp(max(max(abs(Comm))))
%% Diagram 2:
Vh = fem(m,'NED');
Wh = fem(dm,'NED');