Code covered by the BSD License

# Double porosity model

### Jan Valdman (view profile)

07 Aug 2009 (Updated )

Finite element solution of the double porosity model and its a posteriori error estimate

Verfeinerung_Uniform(coordinates,elements3,dirichlet,neumann,nodes2edge)
```function [coordinates,NeueElemente,dirichlet,neumann,fatherelement]=...
Verfeinerung_Uniform(coordinates,elements3,dirichlet,neumann,nodes2edge)
edges_number=full(max(max(nodes2edge)));
nodes_number=size(coordinates,1);
elements3_number=size(elements3,1);

%New array created, it returns two nodes belonging to a given edge
NodenzuKante=zeros(edges_number,2);
[a,b,c]=find(nodes2edge);
NodenzuKante(c,:)=[a b];

% Koordinaten anpassen
coordinates(nodes_number+1:nodes_number+edges_number,:)=(coordinates(NodenzuKante(:,1),:)+coordinates(NodenzuKante(:,2),:))/2;

% Elemente anpassen
NeueElemente=zeros(4*size(elements3,1),3);

fatherelement=zeros(size(NeueElemente,1),1);

ell=0;
for j=1:elements3_number
I=[nodes2edge(elements3(j,1),elements3(j,2)); nodes2edge(elements3(j,2),elements3(j,3)); nodes2edge(elements3(j,3),elements3(j,1))];
NeuKnoten=(I+nodes_number)';
NeueElemente(ell+1:ell+4,:)=[NeuKnoten([2 3 1]);                        %elem1
elements3(j,1), NeuKnoten(1), NeuKnoten(3); %elem2
NeuKnoten(1), elements3(j,2), NeuKnoten(2); %elem3
NeuKnoten(3), NeuKnoten(2), elements3(j,3)];%elem4
fatherelement(ell+1:ell+4)=j;
ell=ell+4;
end

% Dirichlet und Neumann anpassen
tmp=0;
dirichletnr=diag(nodes2edge(dirichlet(:,1),dirichlet(:,2)));
for k=1:size(dirichlet,1)
%if VK(dirichletnr(k))
dirichlet=[dirichlet(1:k-1+tmp,:);dirichlet(k+tmp,1),dirichletnr(k)+nodes_number;
dirichletnr(k)+nodes_number,dirichlet(k+tmp,2);
dirichlet(k+1+tmp:size(dirichlet,1),:)];
tmp=tmp+1;
%end
end
tmp=0;
if ~isempty(neumann)
neumannnr=diag(nodes2edge(neumann(:,1),neumann(:,2)));
end
for k=1:size(neumann,1)
%if VK(Neumannnr(k))
neumann = [neumann(1:k-1+tmp,:);neumann(k+tmp,1),neumannnr(k)+nodes_number;
neumannnr(k)+nodes_number,neumann(k+tmp,2);
neumann(k+1+tmp:size(neumann,1),:)];
tmp=tmp+1;
%end
end

```