...
 
Commits (2)
......@@ -2,6 +2,6 @@ function u=Heat2D(f,g,p)
% HEAT2D solves the stationary heat equation over a square
% u=Heat2D(f,g,p) solves the heat equation over a square of size one
% by one. The source function is given by f, the boundary condition is
% given by b, the number of grid points for each direction is a quarter
% given by g, the number of grid points for each direction is a quarter
% of the length(g) and one chooses the iterative solver with p,
% i.e. p=0 uses Jacobi, p=1 uses Gauss-Seidel.
\documentclass[12pt,a4paper]{article}
\usepackage{amsmath,amssymb,amsthm,epsfig,amscd,verbatim}
\usepackage[french]{babel}
\usepackage{fancyvrb}
%\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\addtolength{\textwidth}{3cm}
\addtolength{\textheight}{4cm}
\addtolength{\hoffset}{-2cm}
\addtolength{\voffset}{-2cm}
\usepackage{hyperref}
%% commandes utiles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\tr}{\mbox{tr}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\noindent UNIVERSIT\'E DE GEN\`EVE \hfill Section de Math\'ematiques \\
\noindent Facult\'e des Sciences \hfill \`A rendre le mercredi 18 mars \\[-3mm] \hrule
\bigskip
\begin{center}
\textbf{\Large{Méthodes Itératives}} \\[1mm]\bigskip
\textbf{Série 2}
\end{center}
\bigskip
\begin{enumerate}
\item (2 pts) Ecrire les méthodes itératives de {\it Jacobi} et {\it Gauss-Seidel}
sous les deux formes standard, de correction, de résidu et de différences.
\item (4 pts) Soit $\vec{B}\in\C^{n\times n}$ une matrice carrée dont le
rayon spectral $\rho(\vec{B})$ vérifie $\rho(\vec{B})<1$. Démontrer que
$\vec{I}-\vec{B}$ est inversible et que
\[
(\vec{I}-\vec{B})^{-1} = \sum_{j=0}^\infty\vec{B}^j.
\]
Cette série est un cas particulier des séries de Neumann (généralisation aux
matrices carrées des séries entières).
\item (3 pts) On parle de matrice à diagonale strictement dominante $A$ si ses
coefficients diagonaux satisfont
\[
|a_{i,i}|>\sum_{j\neq i}|a_{i,j}|,\quad\forall i \in \{1,\ldots,n\}.
\]
Démontrer, pour une telle matrice, que la méthode de Jacobi converge.
\item (2 pts) Consid\'erer la matrice d'it\'eration $G=\begin{bmatrix}
\frac{5}{2} & -1 \\ 1 & 0
\end{bmatrix}$ et le vecteur initial $v_0 = \begin{bmatrix}
1 \\ 2
\end{bmatrix}$.
Calculer les it\'erations $v_1=G v_0$ et $v_2 = Gv_1$ et leurs normes.
Qu'est-ce que vous obs\'ervez?
R\'ep\'eter l'exp\'erience avec $v_0 = \begin{bmatrix}
2 \\ 1
\end{bmatrix}$.
Est-ce que le comportement est le m\^eme?
Pourquoi ou pourquoi pas?
\textit{Indice:} trouver les valeurs et vecteurs propres de $G$ et rappeler le th\'eor\`eme sur la convergence des m\'ethodes stationnaires.
\item (2+4+2 pts) Le produit de Kronecker $\otimes: \mathbb{R}^{m\times n} \times \mathbb{R}^{p \times q} \to \mathbb{R}^{mp \times nq}$ est d\'efini par
\begin{equation*}
A \otimes B = \begin{bmatrix}
a_{11} B & \dots & a_{1n} B \\ \vdots & \ddots & \vdots \\ a_{m1} B & \dots & a_{mn} B
\end{bmatrix}.
\end{equation*}
En Matlab, la fonction \texttt{kron(A,B)} calcule le produit de Kronecker.
\begin{enumerate}
\item Construire le Laplacien discret en dimension deux en utilisant \texttt{kron}.
Montrer qu'il est le m\^eme du r\'esultat de la fonction \texttt{FDLaplacian(n,2)} d'exercice 5 de s\'erie 1.
\item Soit $\left \{ \lambda_i \right \}$ l'ensemble des valeurs propres du Laplacien discret en dimension une.
D\'emontrer que $\left \{ \lambda_i + \lambda_j \right \}$ est l'ensemble des valeurs propres du Laplacien discret en dimension deux.
Quelle est la relation entre les vecteurs propres?
\item Dessiner les premiers cinq vecteurs propres du Laplacien discret en dimension deux et les comparer avec les vecteurs propres exacts.
En utilisant la relation de l'exercice 5b, construire les m\^emes vecteurs propres avec ceux du Laplacien discret en dimension une.
\end{enumerate}
\item (6 pts) Implémenter une fonction Matlab qui résout l'équation de la
chaleur stationnaire
\begin{equation*}
\begin{split}
\Delta u &= f \, \text{ dans $\Omega = (0,1)\times(0,1)$}, \\
u &= g \; \text{ sur $\partial \Omega$}, \\
\end{split}
\end{equation*}
en dimension deux à l'aide des méthodes de {\it Jacobi} et {\it
Gauss-Seidel}. Les solveurs itératifs seront implémentés dans des fonctions
externes au programme résolvant l'équation de la chaleur. Faire des en-têtes
claires et précises dans la description suivant le format ci-dessous.
La fonction Heat2D qui résout l'équation prendra en argument la source $f$, la condition au bord $g$,
et le choix du solveur et résoudra la chaleur
sur un domaine carré de côtés de longueur un. Son en-tête doit ressembler à
la suivante.
\VerbatimInput[firstline=0,lastline=7]{Heat2D.m}
\end{enumerate}
{\bf \'Evaluation:}
\begin{itemize}
\item[$\bullet$] Les exercices
\item[$\bullet$] Un examen oral durant la session d'examens sur le cours.
\end{itemize}
La note finale est de : $30 \%$ exercices et $70\%$ examen oral.\\
Les exercices sont notés $15\%$ pour la présentation de séries et
$15\%$ pour la présentation au tableau. Il faut passer au tableau au
moins 2 fois dans le semestre.
\end{document}
......@@ -191,11 +191,11 @@ Script d'ex
\noindent\\
Figures obtenues:
\begin{center}
\includegraphics[width=0.85\textwidth]{serie3fig1.eps}
\includegraphics[width=0.85\textwidth]{serie3fig2.eps}
\end{center}
%\begin{center}
%\includegraphics[width=0.85\textwidth]{serie3fig1.eps}
%
%\includegraphics[width=0.85\textwidth]{serie3fig2.eps}
%\end{center}
\item \textit{Sur 5 points} Il suffit de construire une suite de matrice
$(\mathbf{A}_k)_{0\leq k\leq n-2}$ et une suite de matrices de
......
function [P,Q,n] = Geometry(X,Y)
% GEOMETRY figures out the intersection of the triangles X and Y. It does
% so by first putting all variables in terms of two of the lines of the
% triangle Y, decides if the points of X lie within or outside Y, then
% calculates the intersections of the lines of X with those of Y and again
% deciding if they are inside or outside Y.
function [P,Q,R,n] = Geometry(U,V)
% GEOMETRY intersection of two triangles
% [P,Q,R,n] = Geometry(U,V) gives the vertices of U inside V (P),
% intersections between the edges of U and those of V (Q), the vertices
% of V inside U (R), and the neighbours of U that also intersect V (n).
ind=[0,0,0];
v0=Y(:,2)-Y(:,1); % vector between y1 and y2
v1=Y(:,3)-Y(:,1); % vector between y1 and y3
d00=v0'*v0; % norm of v0
d01=v0'*v1;
d11=v1'*v1; % norm of v1
v1=V(:,2)-V(:,1); % vector between V1 and V2
v2=V(:,3)-V(:,1); % vector between V1 and V3
d00=v1'*v1; % norm of v1
d01=v1'*v2;
d11=v2'*v2; % norm of v2
id =d00*d11 - d01*d01; % this value is used to normalize quantities
%---Interior points---%
U = zeros(2,3);
%---Change of Coordinates (sec.2)---%
X = zeros(2,3);
for i=1:3
r=X(:,i)-Y(:,1); % vector between ith vertex of triangle X and y1
d02=v0'*r; % projection of v2 in v0 direction
d12=v1'*r; % projection of v2 in v1 direction
U(1,i)=(d11*d02-d01*d12)/id; % amount of v2 in v0 direction
U(2,i)=(d00*d12-d01*d02)/id; % amount of v2 in v1 direction
% v2 = u*v0 + v*v1
% if u and v are positive and their sum less than 1 then v2 points
% inside the simplex created by the vectors v0 and v1
if U(1,i)>=0 && U(2,i)>=0 && sum(U(:,i))<=1 % also include nodes on the boundary
ind(i) = 1;
end
r=U(:,i)-V(:,1); % vector between ith vertex of triangle U and V1
d02=v1'*r; % projection of r in v1 direction
d12=v2'*r; % projection of r in v2 direction
X(1,i)=(d11*d02-d01*d12)/id; % amount of r in v1 direction
X(2,i)=(d00*d12-d01*d02)/id; % amount of r in v2 direction
% end result: r = X(1,i)*v1 + X(2,i)*v2
end
P = X(:,ind==1);
%---Intersections---%
%---Intersections (sec.3)---%
if nargout>1
ind = zeros(1,9);
Q = zeros(2,9);
n = [0,0,0];
Q = zeros(2,9);
indQ = zeros(1,9);
n = [0,0,0];
indR = n; % 1-(0,0); 2-(1,0); 3-(0,1)
indP = ones(1,3);
for i=1:3
x1=U(1,i);
x2=U(2,i);
y1=U(1,mod(i,3)+1);
y2=U(2,mod(i,3)+1);
x = -1; y = -1;
[Qtemp,ntemp,indRtemp,sp] = EdgeIntersect(X,i);
ind = (1:3) + 3*(i-1);
Q(:,ind) = Qtemp;
indQ(ind)= ntemp;
n = max(n,ntemp);
indR = max(indR,indRtemp);
indP = indP.*sp;
end
P = U(:,indP==1);
Q = V(:,1) + [v1,v2]*Q(:,indQ==1);
R = V(:,indR==1);
end
end
% with v1
if sign(x1)~=sign(y1) % check if intersection occurs
y = y2 - y1*(x2-y2)/(x1-y1);
elseif x1==0 && max(x2,y2)>1 && min(x2,y2)<1
% if the lines are not parallel then the only intersections
% that can occur are coincidental lines - if the entire line is
% coincident then the line is covered by Interior points above;
% if not, the only intersection not covered is at the point
% y=1, x=0.
y = 1;
end
if y>=0 && y<=1
Q(:,i) = Y(:,1) + y*v1;
ind(i) = 1;
n(i) = 1;
end
function [Q,ind,indR,sp]=EdgeIntersect(X,edge)
% EDGEINTERSECT intersection of edges of X with reference line of Y
% [Q,ind,indR,sp] = EdgeIntersect(X,edge) calculates the intersections
% between all edges of X and the reference line of Y indicated by edge
% and stores the values in Q.
% Boolean indicators for neighbours of U intersecting V, vertices of V
% inside U, and vertices of U inside V are stored in ind, indR and sp,
% respectively.
% with v0
if sign(x2)~=sign(y2)
x = y1 - y2*(x1-y1)/(x2-y2);
elseif x2==0 && max(x1,y1)>1 && min(x1,y1)<1
x = 1;
% Particulars for each edge of reference triangle
if edge==1 % v1
Param1 = @(x,y) y; Param2 = @(x,y) x; % parametrization of edge
Retrieve = @(q) [q;0]; % inversion of parametrization
R1 = 1; R2 = 2; % vertices of V on edge
elseif edge==2 % v2
Param1 = @(x,y) x; Param2 = @(x,y) y;
Retrieve = @(q) [0;q];
R1 = 1; R2 = 3;
elseif edge==3 % line connecting v1 and v2
Param1 = @(x,y) 1-x-y; Param2 = @(x,y) 0.5*(1-x+y);
Retrieve = @(q) [1-q;q];
R1 = 2; R2 = 3;
end
Q = zeros(2,3); ind = zeros(1,3); indR = ind;
p = Param1(X(1,:),X(2,:)); q = NaN*ones(1,3);
sp= sign(p)>=0;
for i = 1:3
j = mod(i,3) + 1;
if sp(i) ~= sp(j) % test if vertices of X lie on opposite sides of edge
if isnan(q(i))% test if q-coord. of vertex already calculated
q(i) = Param2(X(1,i),X(2,i));
end
if x>=0 && x<=1
Q(:,3+i) = Y(:,1) + x*v0;
ind(3+i) = 1;
n(i) = 1;
if isnan(q(j))
q(j) = Param2(X(1,j),X(2,j));
end
% with the line connecting v0 and v1
if sign(x1+x2-1)~=sign(y1+y2-1)
x = ( (x1-y1) + y1*(x2-y2) - y2*(x1-y1) )/( (x1-y1) + (x2-y2) );
y = ( (x2-y2) + y2*(x1-y1) - y1*(x2-y2) )/( (x1-y1) + (x2-y2) );
if x>=0 && y>=0 && x<=1 && y<=1
Q(:,6+i) = Y(:,1) + x*v0 + y*v1;
ind(6+i) = 1;
n(i) = 1;
end
q0=(q(i)*p(j) - q(j)*p(i))/(p(j)-p(i)); % calculate intersection
if q0>=0 && q0<=1 % test if intersection lies on edge of Y
Q(:,i) = Retrieve(q0);
ind(i) = 1;
end
if exist('qold','var') % test if 1st intersect. already calculated
indR(R1) = indR(R1) | sign(qold)~=sign(q0);
indR(R2) = indR(R2) | sign(qold-1)~=sign(q0-1);
end
qold = q0;
end
Q = Q(:,ind==1);
end
end
\ No newline at end of file
......@@ -4,8 +4,8 @@ function M=InterfaceMatrix(Na,Ta,Nb,Tb)
% and Tb with associated nodal coordinates in Na and Nb and
% computes the interface projection matrix M
bl=[1]; % bl: list of triangles of Tb to treat
bil=[1]; % bil: list of triangles Ta to start with
bl =1; % bl: list of triangles of Tb to treat
bil=1; % bil: list of triangles Ta to start with
bd=zeros(size(Tb,1)+1,1); % bd: flag for triangles in Tb treated
bd(end)=1; % guard, to treat boundaries
bd(1)=1; % mark first triangle in b list.
......
......@@ -4,18 +4,13 @@ function M=InterfaceMatrixNew(Na,Ta,Nb,Tb)
% and Tb with associated nodal coordinates in Na and Nb and
% computes the interface projection matrix M
% bl and bil need to be changed so that the first pair of triangles has an
% intersection - this seems like a real crutch
% replacing the new meshes in example with those from Jeronimo requires bl
% to be 15 (or so?)
bl =1; % bl: list of triangles of Tb to treat
bil=1; % bil: list of triangles Ta to start with
bd=zeros(size(Tb,1)+1,1); % bd: flag for triangles in Tb treated
bd(end)=1; % guard, to treat boundaries
bd(1)=1; % mark first triangle in b list.
M=sparse(size(Nb,2),size(Na,2));
while isempty(bl)~=1
while ~isempty(bl)
bc=bl(1); bl=bl(2:end); % bc: current triangle of Tb
al=bil(1); bil=bil(2:end); % triangle of Ta to start with
ad=zeros(size(Ta,1)+1,1); % same as for bd
......
......@@ -4,18 +4,13 @@ function M=InterfaceMatrixRF(Na,Ta,Nb,Tb)
% and Tb with associated nodal coordinates in Na and Nb and
% computes the interface projection matrix M
% bl and bil need to be changed so that the first pair of triangles has an
% intersection - this seems like a real crutch
% replacing the new meshes in example with those from Jeronimo requires bl
% to be 15 (or so?)
bl =1; % bl: list of triangles of Tb to treat
bil=1; % bil: list of triangles Ta to start with
bd=zeros(size(Tb,1)+1,1); % bd: flag for triangles in Tb treated
bd(end)=1; % guard, to treat boundaries
bd(1)=1; % mark first triangle in b list.
M=sparse(size(Nb,2),size(Na,2));
while isempty(bl)~=1
while ~isempty(bl)
bc=bl(1); bl=bl(2:end); % bc: current triangle of Tb
al=bil(1); bil=bil(2:end); % triangle of Ta to start with
ad=zeros(size(Ta,1)+1,1); % same as for bd
......@@ -52,16 +47,16 @@ function [P,n,M]=Intersect(X,Y)
% element Y. The numerical challenges are handled by including
% points on the boundary and removing duplicates at the end.
[P,Q,R,n] = RefFree(X,Y);
if size(P,2)>1 % if two or more interior points
[P,Q,R,n] = RefFree(Y,X);
if size(P,2)>1 %|| size(R,2)>1 % if two or more interior points
n=[1,1,1]; % the triangle is candidate for all
end % neighbors
P=[P Q R];
P=SortAndRemoveDoubles(P); % sort counter clock wise
[Pn,Qn,Rn,nn] = Geometry(X,Y);
Pn = SortAndRemoveDoubles([Pn,Qn,Rn]);
% [Pn,Qn,Rn,nn] = Geometry(X,Y);
% Pn = SortAndRemoveDoubles([Pn,Qn,Rn]);
M=zeros(3,3);
if size(P,2)>0
for j=2:size(P,2)-1 % compute interface matrix
......@@ -74,15 +69,15 @@ if size(P,2)>0
pause(0)
end
if size(P,2)>0 && size(Pn,2)>0
polyP = [ P, P(:,1)];
polyPn= [Pn,Pn(:,1)];
plot(polyP(1,:),polyP(2,:),'r*',polyPn(1,:),polyPn(2,:),'bo','markersize',10)
% triX = [ X , X(:,1)];
% triY = [ Y , Y(:,1)];
% plot(triX(1,:),triX(2,:),'-',triY(1,:),triY(2,:),'-')
pause(0.5)
end
% if size(P,2)>0 && size(Pn,2)>0
% polyP = [ P, P(:,1)];
% polyPn= [Pn,Pn(:,1)];
% plot(polyP(1,:),polyP(2,:),'r*',polyPn(1,:),polyPn(2,:),'bo','markersize',10)
% % triX = [ X , X(:,1)];
% % triY = [ Y , Y(:,1)];
% % plot(triX(1,:),triX(2,:),'-',triY(1,:),triY(2,:),'-')
% pause(0.5)
% end
end
function [P,Q,R,n]=RefFree(U,V)
......@@ -103,7 +98,7 @@ for i=1:3
q0=(qp(1,k)*qp(2,l)-qp(1,l)*qp(2,k))/(qp(2,l)-qp(2,k));
if q0>=0 && q0<=1 % intersection lies on edge of Y?
Q(:,m)=V(:,i)+q0*w; indQ(m)=1;
n(k)=1;
n(i)=1;
end
if isnan(qold) % 1st intersection not calculated?
qold=q0;
......
function M=QuadComplex(Na,Ta,Nb,Tb)
[m,~] = size(Ta);
[n,~] = size(Tb);
M=sparse(size(Nb,2),size(Na,2));
for i=1:m
for j=1:n
[P,nc,Mc] = Intersect(Nb(:,Tb(j,1:3)),Na(:,Ta(i,1:3)));
if ~isempty(P)
M(Tb(j,1:3),Ta(i,1:3))=M(Tb(j,1:3),Ta(i,1:3))+Mc;
end
end
end
end
function [P,n,M]=Intersect(X,Y)
% INTERSECT intersection of two triangles and mortar contribution
% [P,n,M]=Intersect(X,Y); computes for the two given triangles X and
% Y (point coordinates are stored column-wise, in counter clock
% order) the points P where they intersect, in n the indices of
% which neighbors of X are also intersecting with Y, and the local
% mortar matrix M of contributions of the P1 elements X on the P1
% element Y. The numerical challenges are handled by including
% points on the boundary and removing duplicates at the end.
[P,n]=EdgeIntersections(X,Y);
P1=PointsOfXInY(X,Y);
if size(P1,2)>1 % if two or more interior points
n=[1 1 1]; % the triangle is candidate for all
end % neighbors
P=[P P1];
P=[P PointsOfXInY(Y,X)];
P=SortAndRemoveDoubles(P); % sort counter clock wise
M=zeros(3,3);
if size(P,2)>0
for j=2:size(P,2)-1 % compute interface matrix
M=M+MortarInt(P(:,[1 j j+1]),X,Y);
end
patch(P(1,:),P(2,:),'m') % draw intersection for illustration
% H=line([P(1,:) P(1,1)],[P(2,:),P(2,1)]);
% set(H,'LineWidth',3,'Color','m');
pause(0)
end
end
function P=PointsOfXInY(X,Y)
% POINTSOFXINY finds corners of one triangle within another one
% P=PointsOfXInY(X,Y); computes for the two given triangles X
% and Y (point coordinates are stored column-wise, in counter clock
% order) the corners P of X which lie in the interior of Y.
k=0;P=[];
v0=Y(:,2)-Y(:,1); v1=Y(:,3)-Y(:,1); % find interior points of X in Y
d00=v0'*v0; d01=v0'*v1; d11=v1'*v1; % using baricentric coordinates
id=1/(d00*d11-d01*d01);
for i=1:3
v2=X(:,i)-Y(:,1); d02=v0'*v2; d12=v1'*v2;
u=(d11*d02-d01*d12)*id; v=(d00*d12-d01*d02)*id;
if u>=0 && v>=0 && u+v<=1 % also include nodes on the boundary
k=k+1; P(:,k)=X(:,i);
end
end
end
function [P,n]=EdgeIntersections(X,Y)
% EDGEINTERSECTIONS computes edge intersections of two triangles
% [P,n]=EdgeIntersections(X,Y) computes for the two given triangles X
% and Y (point coordinates are stored column-wise, in counter clock
% order) the points P where their edges intersect. In addition,
% in n the indices of which neighbors of X are also intersecting
% with Y are given.
P=[]; k=0; n=[0 0 0];
for i=1:3 % find all intersections of edges
for j=1:3
b=Y(:,j)-X(:,i);
A=[X(:,mod(i,3)+1)-X(:,i) -Y(:,mod(j,3)+1)+Y(:,j)];
if rank(A)==2 % edges not parallel
r=A\b;
if r(1)>=0 && r(1)<=1 && r(2)>=0 && r(2)<=1 % intersection found
k=k+1; P(:,k)=X(:,i)+r(1)*(X(:,mod(i,3)+1)-X(:,i)); n(i)=1;
end
end
end
end
end
function P=SortAndRemoveDoubles(P)
% SORTANDREMOVEDOUBLES sort points and remove duplicates
% P=SortAndRemoveDoubles(P); orders polygon corners in P counter
% clock wise and removes duplicates
ep=10*eps; % tolerance for identical nodes
m=size(P,2);
if m>0
c=sum(P')'/m; % order polygon corners counter
for i=1:m % clockwise
d=P(:,i)-c; ao(i)=angle(d(1)+sqrt(-1)*d(2));
end
[tmp,id]=sort(ao);
P=P(:,id);
i=1;j=2; % remove duplicates
while j<=m
if norm(P(:,i)-P(:,j))>ep
i=i+1;P(:,i)=P(:,j);j=j+1;
else
j=j+1;
end
end
P=P(:,1:i);
end
end
function M=MortarInt(T,X,Y)
% MORTARINT computes mortar contributions
% M=MortarInt(T,X,Y); computes for triangles X and Y with nodal coordinates
% stored columnwise the integrals of all P1 element shape function combinations
% between triangles X and Y on triangle T.
% The result is stored in the 3 by 3 matrix M
Jd=-T(1,1)*T(2,3)-T(1,2)*T(2,1)+T(1,2)*T(2,3)+T(1,1)*T(2,2)+...
T(1,3)*T(2,1)-T(1,3)*T(2,2);
a=T(1,1);
d=T(2,1);
e=-T(2,1)+T(2,2);
b=-T(1,1)+T(1,2);
c=-T(1,1)+T(1,3);
f=-T(2,1)+T(2,3);
T1=1/2*(X(1,1)*(X(2,2)-X(2,3))-X(1,2)*(X(2,1)-X(2,3))+X(1,3)*(X(2,1)-X(2,2)));
T2=1/2*(Y(1,1)*(Y(2,2)-Y(2,3))-Y(1,2)*(Y(2,1)-Y(2,3))+Y(1,3)*(Y(2,1)-Y(2,2)));
A11=1/(2*abs(T1))*((X(2,2)-X(2,3))*b-(X(1,2)-X(1,3))*e);
B11=1/(2*abs(T1))*((X(2,2)-X(2,3))*c-(X(1,2)-X(1,3))*f);
C11=1/(2*abs(T1))*((X(2,2)-X(2,3))*a-(X(1,2)-X(1,3))*d+(X(1,2)*X(2,3)-X(2,2)*X(1,3)));
A12=1/(2*abs(T1))*((X(2,3)-X(2,1))*b+(X(1,1)-X(1,3))*e);
B12=1/(2*abs(T1))*((X(2,3)-X(2,1))*c+(X(1,1)-X(1,3))*f);
C12=1/(2*abs(T1))*((X(2,3)-X(2,1))*a+(X(1,1)-X(1,3))*d-(X(1,1)*X(2,3)-X(2,1)*X(1,3)));
A13=1/(2*abs(T1))*((X(2,1)-X(2,2))*b-(X(1,1)-X(1,2))*e);
B13=1/(2*abs(T1))*((X(2,1)-X(2,2))*c-(X(1,1)-X(1,2))*f);
C13=1/(2*abs(T1))*((X(2,1)-X(2,2))*a-(X(1,1)-X(1,2))*d+(X(1,1)*X(2,2)-X(2,1)*X(1,2)));
A21=1/(2*abs(T2))*((Y(2,2)-Y(2,3))*b-(Y(1,2)-Y(1,3))*e);
B21=1/(2*abs(T2))*((Y(2,2)-Y(2,3))*c-(Y(1,2)-Y(1,3))*f);
C21=1/(2*abs(T2))*((Y(2,2)-Y(2,3))*a-(Y(1,2)-Y(1,3))*d+(Y(1,2)*Y(2,3)-Y(2,2)*Y(1,3)));
A22=1/(2*abs(T2))*((Y(2,3)-Y(2,1))*b+(Y(1,1)-Y(1,3))*e);
B22=1/(2*abs(T2))*((Y(2,3)-Y(2,1))*c+(Y(1,1)-Y(1,3))*f);
C22=1/(2*abs(T2))*((Y(2,3)-Y(2,1))*a+(Y(1,1)-Y(1,3))*d-(Y(1,1)*Y(2,3)-Y(2,1)*Y(1,3)));
A23=1/(2*abs(T2))*((Y(2,1)-Y(2,2))*b-(Y(1,1)-Y(1,2))*e);
B23=1/(2*abs(T2))*((Y(2,1)-Y(2,2))*c-(Y(1,1)-Y(1,2))*f);
C23=1/(2*abs(T2))*((Y(2,1)-Y(2,2))*a-(Y(1,1)-Y(1,2))*d+(Y(1,1)*Y(2,2)-Y(2,1)*Y(1,2)));
M(1,1)=1/24*Jd*(2*A11*A21+B11*A21+A11*B21+2*B11*B21+4*C11*A21+4*A11*C21+4*C11*B21+4*B11*C21+12*C11*C21);
M(1,2)=1/24*Jd*(2*A11*A22+B11*A22+A11*B22+2*B11*B22+4*C11*A22+4*A11*C22+4*C11*B22+4*B11*C22+12*C11*C22);
M(1,3)=1/24*Jd*(2*A11*A23+B11*A23+A11*B23+2*B11*B23+4*C11*A23+4*A11*C23+4*C11*B23+4*B11*C23+12*C11*C23);
M(2,1)=1/24*Jd*(2*A12*A21+B12*A21+A12*B21+2*B12*B21+4*C12*A21+4*A12*C21+4*C12*B21+4*B12*C21+12*C12*C21);
M(2,2)=1/24*Jd*(2*A12*A22+B12*A22+A12*B22+2*B12*B22+4*C12*A22+4*A12*C22+4*C12*B22+4*B12*C22+12*C12*C22);
M(2,3)=1/24*Jd*(2*A12*A23+B12*A23+A12*B23+2*B12*B23+4*C12*A23+4*A12*C23+4*C12*B23+4*B12*C23+12*C12*C23);
M(3,1)=1/24*Jd*(2*A13*A21+B13*A21+A13*B21+2*B13*B21+4*C13*A21+4*A13*C21+4*C13*B21+4*B13*C21+12*C13*C21);
M(3,2)=1/24*Jd*(2*A13*A22+B13*A22+A13*B22+2*B13*B22+4*C13*A22+4*A13*C22+4*C13*B22+4*B13*C22+12*C13*C22);
M(3,3)=1/24*Jd*(2*A13*A23+B13*A23+A13*B23+2*B13*B23+4*C13*A23+4*A13*C23+4*C13*B23+4*B13*C23+12*C13*C23);
end
\ No newline at end of file
......@@ -16,7 +16,7 @@ for i=1:3
q0=(qp(1,k)*qp(2,l)-qp(1,l)*qp(2,k))/(qp(2,l)-qp(2,k));
if q0>=0 && q0<=1 % intersection lies on edge of Y?
Q(:,m)=V(:,i)+q0*w; indQ(m)=1;
n(k)=1;
n(i)=1; % nhbrs of V intersect U?
end
if isnan(qold) % 1st intersection not calculated?
qold=q0;
......@@ -28,4 +28,5 @@ for i=1:3
end
indP=indP.*sp;
end
P=U(:,indP==1); Q=Q(:,indQ==1); R=V(:,indR==1);
\ No newline at end of file
P=U(:,indP==1); Q=Q(:,indQ==1); R=V(:,indR==1);
end
\ No newline at end of file
[N,T,Np]=CircleMesh(20,1e-16);
[N,T,Np]=CircleMesh(10,eps);
%%
figure(1)
......@@ -6,6 +6,7 @@ clf
PlotMesh(N,T,'b')
PlotMesh(Np,T,'r')
M = InterfaceMatrix(N,T,Np,T);
% M = QuadComplex(N,T,Np,T);
pause(1)
% savefig('CirclePANG.fig')
......
......@@ -5,4 +5,6 @@ set(fig,'DoubleBuffer','on');
PlotMesh(N_P,T_P,'b');
PlotMesh(N_S,T_S,'r');
perturb = 0 * (rand(2,3)-0.5);
M=InterfaceMatrixNew(N_P,T_P,N_S + perturb,T_S);
\ No newline at end of file
% M=InterfaceMatrixNew(N_P,T_P,N_S + perturb,T_S);
% M=QuadComplex(N_P,T_P,N_S+perturb,T_S);
M=InterfaceMatrixRF(N_P,T_P,N_S+perturb,T_S);
\ No newline at end of file