Code covered by the BSD License  

Highlights from
Double porosity model

image thumbnail

Double porosity model

by

 

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

[coordinates,NeueElemente,dirichlet,neumann,fatherelement]=...
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
 

Contact us