Commit f3dcba4a authored by MartinA's avatar MartinA
Browse files

Laplace beltrami on multi-screen

parent 2bae2bbb
......@@ -42,7 +42,7 @@ methods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONSTRUCTOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fe = fem(mesh,str)
if nargin == 0
% Needed for child class. Martin A 22/09/2020
% Needed for child class
return
end
fe.typ = str;
......
......@@ -68,6 +68,11 @@ elseif strcmp(fe.typ,'RWG')
[face,elt2dof] = fe.msh.fce;
X = face.ctr;
end
elseif strcmp(fe.typ(1:5),'multi')
[X,elt2dof] = dofMS(fe);
% Others
else
......
......@@ -53,6 +53,8 @@ elseif strcmp(fe.typ,'RWG')
M = femRaoWiltonGlisson(fe,domain);
elseif strcmp(fe.typ,'NED')
M = femNedelec(fe,domain);
elseif strcmp(fe.typ(1:5),'multi')
M = femMultiScreen(fe,domain);
else
error('fem.m : unavailable case')
end
......
......@@ -618,13 +618,15 @@ classdef msh
end
% Rotation
function mesh = rotate(mesh,U,phi)
% Rotate by angle phi around vector U in trigonometric
function mesh = rotate(mesh,center,U,phi)
% Rotate by angle phi around axis center + kU in trigonometric
% direction.
N = U./norm(U);
mesh.vtx = cos(phi) * mesh.vtx + ...
(1-cos(phi)) .* ((mesh.vtx*N')*N) + ...
sin(phi) .* cross(ones(size(mesh.vtx,1),1)*N,mesh.vtx,2);
mtemp = translate(mesh,-center);
mtemp.vtx = cos(phi) * mtemp.vtx + ...
(1-cos(phi)) .* ((mtemp.vtx*N')*N) + ...
sin(phi) .* cross(ones(size(mtemp.vtx,1),1)*N,mtemp.vtx,2);
mesh = translate(mtemp,center);
end
% Swap
......
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_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
Gamma = dom(m_Gamma,gss);
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); % 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); % 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); % Normal trace
mu = sol(end); % Lagrange multiplier.
fprintf('Lagrange Multiplier: %s\n',num2str(mu))
dnu = P*lambda; % Restriction to Gamma2
figure;
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);
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;
fprintf('Force (in the x direction) %s\n',num2str(F))
\ No newline at end of file
function [dof,elt2dof,restrictions] = dofMS(Vh)
typ = Vh.typ(7:end);
ms = Vh.msh;
pans = ms.panels;
restrictions = cell(ms.npanels,1);
elt2dof= cell(ms.npanels,1);
X = cell(ms.npanels,1);
if strcmp(typ,'P1')
dof = [];
for i = 1:ms.npanels
femi = fem(pans{i},typ);
[X{i},elt2dof{i}] = femi.dof;
dof = [dof; [X{i},i*ones(size(X{i},1),1)]];
end
bound = ms.bnd;
[~,jdx] = ismembertol(dof(:,1:3),bound.vtx,'ByRows',true);
dof = dof(~jdx,:);
dofBound = zeros(bound.nvtx,4);
dofBound(:,1:3) = bound.vtx; dofBound(:,4) = 0;
dof = [dofBound;dof];
for i = 1:ms.npanels
Xi = [X{i},ones(size(X{i},1),1)*i];
[~,jdx] = ismembertol(Xi,dof,'ByRows',true);
idx1 = find(jdx~=0);
jdx1 = jdx(jdx~=0);
Xi(:,4) = 0;
[~,jdx] = ismembertol(Xi,dof,'ByRows',true);
idx2 = find(jdx~=0);
jdx2 = jdx(jdx~=0);
idx = [idx1;idx2];
jdx = [jdx1;jdx2];
elt2dof{i}(idx,:) = jdx(elt2dof{i}(idx,:));
restrictions{i} = sparse(idx,jdx,1,size(X{i},1),size(dof,1));
end
else
error('unavailable case')
end
end
function [dof,elt2dof,restrictions] = dofMS(Vh)
typ = Vh.typ(7:end);
ms = Vh.msh;
pans = ms.panels;
restrictions = cell(ms.npanels,1);
elt2dof= cell(ms.npanels,1);
X = cell(ms.npanels,1);
if strcmp(typ,'P1')
dof = [];
for i = 1:ms.npanels
femi = fem(pans{i},typ);
[X{i},elt2dof{i}] = femi.dof;
dof = [dof; [X{i},i*ones(size(X{i},1),1)]];
end
bound = ms.bnd;
[~,jdx] = ismembertol(dof(:,1:3),bound.vtx,'ByRows',true);
dof = dof(~jdx,:);
dofBound = zeros(bound.nvtx,4);
dofBound(:,1:3) = bound.vtx; dofBound(:,4) = 0;
dof = [dofBound;dof];
for i = 1:ms.npanels
Xi = [X{i},ones(size(X{i},1),1)*i];
[~,jdx] = ismembertol(Xi,dof,'ByRows',true);
idx1 = find(jdx~=0);
jdx1 = jdx(jdx~=0);
Xi(:,4) = 0;
[~,jdx] = ismembertol(Xi,dof,'ByRows',true);
idx2 = find(jdx~
elt2dof{i} = jdx(elt2dof{i});
restrictions{i} = sparse(1:size(X{i},1),jdx,1,size(X{i},1),size(dof,1));
end
else
error('unavailable case')
end
end
function [P] = femMultiScreen(Xfem,Xdom)
% Dof to quadrature matrix
ms = Xdom.msh;
gs = Xdom.gss;
[~,W,I] = qud(Xdom);
[~,~,restrictions] = dofMS(Xfem);
pani = ms.panels{1};
domi = dom(pani,gs);
femi = fem(pani,Xfem.typ(7:end));
femi.opr = Xfem.opr;
uqmi = femi.uqm(domi);
idx = I{1};
jdx = 1:size(domi.qud,1);
Q = sparse(idx,jdx,1,length(W),length(jdx));
if iscell(uqmi)
P = cell(length(uqmi),1);
for d = 1:length(uqmi)
P{d} = Q*uqmi{d}*restrictions{1};
end
else
P = Q*uqmi*restrictions{1};
end
for i = 2:ms.npanels
pani = ms.panels{i};
domi = dom(pani,gs);
femi = fem(pani,Xfem.typ(7:end));
femi.opr = Xfem.opr;
uqmi = femi.uqm(domi);
idx = I{i};
jdx = 1:size(domi.qud,1);
Q = sparse(idx,jdx,1,length(W),length(jdx));
if iscell(uqmi)
for d = 1:length(uqmi)
P{d} = P{d} + Q*uqmi{d}*restrictions{i};
end
else
P = P + Q*uqmi*restrictions{i};
end
end
end
\ No newline at end of file
classdef msDom < dom
%UNTITLED2 Summary of this class goes here
% Detailed explanation goes here
properties
end
methods
% Class constructor
function this = msDom(varargin)
if nargin == 0
else
this.msh = varargin{1};
this.gss = varargin{2};
end
end
function [X,W,I] = qud(this)
X = []; W = [];
gs = this.gss;
ms = this.msh;
I = cell(ms.npanels,1);
start = 0;
for i = 1:ms.npanels
pani = ms.panels{i};
domi = dom(pani,gs);
[x,w] = qud(domi);
X = [X;x]; W = [W;w];
I{i} = (start+1):(start+length(w));
start = start+length(w);
end
end
end
end
classdef msFem < fem
%UNTITLED Summary of this class goes here
% Detailed explanation goes here
properties
end
methods
function[] = surf(this,val)
[~,~,restrictions] = dofMS(this);
typ = this.typ(7:end);
for i = 1:this.msh.npanels
pani = this.msh.panels{i};
Nrm = pani.nrm;
n = mean(Nrm,1);
pani = translate(pani,-1e-3*n);
femi = fem(pani,typ);
surf(femi,restrictions{i}*val);
hold on;
end
end
end
end
classdef multiScreen
properties
panels; % List of meshes of each Gamma_j
vtx;
I
end
methods
% Class Constructor
function[this] = multiScreen(pan)
if nargin == 0
this.panels = cell(0);
else
this.panels = pan;
vtx = [];
for i = 1:length(pan)
vtx = [vtx; pan{i}.vtx];
end
[this.vtx,this.I] = unique(vtx,'rows');
end
end
% Methods
function np = npanels(this)
np = length(this.panels);
end
function[] = plot(varargin)
this = varargin{1};
plot(this.panels{1},varargin{2:end});
hold on
for i = 2:this.npanels
plot(this.panels{i},varargin{2:end})
end
end
function[this] = rotate(this,center,U,phi)
for i = 1:this.npanels
this.panels{i} = rotate(this.panels{i},center,U,phi);
end
end
function[E] = elt(this)
E = this.singleMesh.elt;
end
function[sg] = singleMesh(this)
sg = this.panels{1};
for i = 2:this.npanels
sg = union(sg,this.panels{i});
end
sg = mshClean(sg);
end
function[b] = bnd(this)
b = bnd(this.singleMesh);
end
end
end
close all;
ms = nFoldMultiScreen(500,3);
figure
plot(ms);
hold on
plotNrm(ms.panels{1},'r');
plotNrm(ms.panels{2},'g');
plotNrm(ms.panels{3},'b');
plot(ms.bnd,'y');
phi = msFem(ms,'multi_P1');
Gamma = msDom(ms,3);
M = integral(Gamma,phi,phi);
K = integral(Gamma,grad(phi),grad(phi));
[P,D] = eig(full(M\K));
[d,ind] = sort(diag(D));
figure;
surf(phi,P(:,ind(4)));
axis equal
\ No newline at end of file
function [ms] = nFoldMultiScreen(N,n)
Q = round(N/n);
originalPanel = mshSquare(Q,[1,2]);
singlePanel1 = originalPanel;
ang = 2*pi/n;
panels = cell(n,1);
for i = 1:n-1
singlePanel2 = swap(rotate(singlePanel1,[-0.5,0,0],[0,1,0],ang));
panels{i} = mshClean(union(singlePanel1,singlePanel2));
singlePanel1 = swap(singlePanel2);
end
singlePanel2 = originalPanel;
if mod(n,2)==1
singlePanel2 = swap(singlePanel2);
end
panels{end} = mshClean(union(singlePanel1,singlePanel2));
ms = multiScreen(panels);
ms = rotate(ms,[0,0,0],[1,0,0],pi/2);
end
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