Commit 4d6c9ee0 authored by MartinA's avatar MartinA
Browse files

Repaired the 3d regularization that was affected by my changes on the msh.tgt definition

parent f3dcba4a
......@@ -53,7 +53,9 @@ ctr = Ydom.msh.ctr;
nrm = Ydom.msh.nrm;
stp = Ydom.msh.stp;
srf = Ydom.msh.ndv;
tau = cell2mat(Ydom.msh.tgt);
tau = Ydom.msh.tgt;
tau2{1} = tau{2}; tau2{2} = tau{3}; tau2{3} = tau{1};
tau = cell2mat(tau2);
nu = cell2mat(Ydom.msh.nrmEdg);
Nelt = size(elt,1);
......
......@@ -341,6 +341,9 @@ classdef msh
% Edge normals
function Nu = nrmEdg(mesh,I)
if ~exist('I','var')||isempty(I)
I = 1:mesh.nelt;
end
switch mesh.type
case 'point'
error('No meaningful edge normal on point meshes')
......
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
......@@ -56,7 +56,9 @@ classdef multiScreen
function[b] = bnd(this)
b = bnd(this.singleMesh);
end
function[M] = regularize(varargin)
end
end
......
close all;
ms = nFoldMultiScreen(500,3);
figure
ms = nFoldMultiScreen(600,3);
figure
plot(ms);
hold on
plotNrm(ms.panels{1},'r');
......@@ -20,4 +20,30 @@ K = integral(Gamma,grad(phi),grad(phi));
figure;
surf(phi,P(:,ind(4)));
axis equal
\ No newline at end of file
axis equal
Gxy = @(X,Y) femGreenKernel(X,Y,'[1/r]',0);
N = 0;
[~,~,restrictions] = dofMS(phi);
for i = 1:ms.npanels
Gammai = dom(ms.panels{i},3);
Ri = restrictions{i};
phii = fem(ms.panels{i},'P1');
for j= 1:ms.npanels
Rj = restrictions{j};
phij = fem(ms.panels{j},'P1');
Gammaj = dom(ms.panels{j},3);
N = N + 1/(4*pi)*Ri'*integral(Gammai,Gammaj,nxgrad(phii),Gxy,nxgrad(phij))*Rj;
N = N + 1/(4*pi)*Ri'*regularize(Gammai,Gammaj,nxgrad(phii),'[1/r]',nxgrad(phij))*Rj;
end
end
[P,D] = eig(full(M\N));
[d,ind] = sort(diag(D));
i = find(d> 4,1,'first');
figure;
surf(phi,P(:,ind(i)));
axis equal;
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