start.m

clear all
%format long

%fluid parameters
k1=1; k2=4; kappa=10;

%exact solution used for generation of right-hand sides h1 and h2 later
l1=2; %reference solution is: p1=sin(l1*pi*x)*sin(l1*pi*y)
l2=3; %reference solution is: p2=sin(l2*pi*x)*sin(l2*pi*y)

%maximal uniform refinement level
levels=5; %9; %maximum level for 2Gb memory

%Fridrichs' constant for H1_0 space on the unit square
%C_omega_exact=sqrt(1/(2*pi*pi)); %exact value for convergence comparison

%coarse mesh of the unit square
coordinates=[0 0; 1 0; 0 1; 1 1];
elements=[1 2 4; 1 4 3];  
dirichlet=[1 2; 2 4; 4 3; 3 1];
neumann=[];

%computation of double porosity problem on nested meshes
for level=1:levels
    fprintf('level=%d',level); 
    %uniform refinement   
    [element2edges, edge2nodes]=getEdges(elements);    
    nodes2edge=sparse(edge2nodes(:,1),edge2nodes(:,2),1:size(edge2nodes,1),size(coordinates,1),size(coordinates,1));
    nodes2edge=symetrizeMatrix(nodes2edge);
    %edge2elements=entryInWhichRows(element2edges); 
    [coordinates,elements,dirichlet,neumann]=refinement_uniform(coordinates,elements,dirichlet,neumann,nodes2edge);         
    number_elements=size(elements,1); 
    number_coordinates=size(coordinates,1);
    
    %stifness and mass matrices generation
    Xscalar=kron(ones(1,3),elements); Yscalar=kron(elements,ones(1,3)); Zscalar=zeros(size(Xscalar));
    areas=zeros(number_elements,1);   
    for j = 1:number_elements
        I=elements(j,:); 
        verticesTrans=coordinates(I,:)'; 
        area=det([1 1 1; verticesTrans])/2; areas(j)=area;
        PhiGrad = [1 1 1; verticesTrans]\[0 0; 1 0; 0 1];        
        Alocal = area * PhiGrad * PhiGrad'; 
        Zscalar(j,:)=reshape(Alocal,1,9); 
    end
    Zmassmatrix=kron(areas,reshape((ones(3)+eye(3))/12,1,9));
    STIF=sparse(Xscalar,Yscalar,Zscalar); 
    MASS=sparse(Xscalar,Yscalar,Zmassmatrix);
    
    all_nodes=(1:number_coordinates)';
    dirichlet_nodes=unique(dirichlet);
    free_nodes=setdiff(all_nodes,dirichlet_nodes); %internal + pure Neumann nodes 
    free_nodes2D=[free_nodes;free_nodes+number_coordinates];
    
    %system matrix 
    A=[k1*STIF+kappa*MASS -kappa*MASS; -kappa*MASS k2*STIF+kappa*MASS]; 
    
    %right hand side 
    p1_exact=sin(coordinates(:,1)*pi*l1).*sin(coordinates(:,2)*pi*l1);
    p2_exact=sin(coordinates(:,1)*pi*l2).*sin(coordinates(:,2)*pi*l2);
    h1=2*k1*pi^2*l1^2*p1_exact+kappa*(p1_exact-p2_exact);
    h2=2*k2*pi^2*l2^2*p2_exact+kappa*(p2_exact-p1_exact);
    b=[MASS*h1; MASS*h2];
    
    %solution of homegeneous Dirichlet problem
    p=b*0;
    p(free_nodes2D)=A(free_nodes2D,free_nodes2D)\b(free_nodes2D);
    p1=p(1:number_coordinates); 
    p2=p(number_coordinates+1:end);
    
    %approximation error in a problem related norm
    error=sqrt(k1*(p1-p1_exact)'*STIF*(p1-p1_exact)+k2*(p2-p2_exact)'*STIF*(p2-p2_exact)...
          +kappa*(p1-p1_exact)'*MASS*(p1-p1_exact)+kappa*(p2-p2_exact)'*MASS*(p2-p2_exact)...
          -2*kappa*(p1-p1_exact)'*MASS*(p2-p2_exact));
    fprintf(', approx. error=%f',error);  
    
    %functional majorant (not finished!)
    dual_F=kappa*(p1'*MASS*p1+p2'*MASS*p2-2*p1'*MASS*p2)/2 -h1'*MASS*p1-h2'*MASS*p2; %not complete 
    dual_G=k1*p1'*STIF*p1+k2*p2'*STIF*p2; %not complete 
    majorant=dual_F+dual_G;
    fprintf(', majorant (only part)=%f',majorant);  
    
    fprintf(', dof=%d \n',2*size(coordinates,1));  
    
    %visualization
    figure(1);
    subplot(2,1,1);
    show_piecewise_linear(p1(elements),elements,coordinates);title('pressure p_1');
    subplot(2,1,2);
    show_piecewise_linear(p2(elements),elements,coordinates);title('pressure p_2');
    
    %subplot(2,2,3);
    %show_piecewise_linear(h1(elements),elements,coordinates);title('right-hand side p1');
    %subplot(2,2,4);
    %show_piecewise_linear(h2(elements),elements,coordinates);title('right-hand side p2');
    
end