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

[DIVDIVmatrixRT0 MASSmatrixRT0]=stima_divdivRT0(elements3,coordinates,areas,element2edges,edge2elements)
function [DIVDIVmatrixRT0 MASSmatrixRT0]=stima_divdivRT0(elements3,coordinates,areas,element2edges,edge2elements)

orientation=edges_orientation(element2edges,edge2elements);

Mhelp=2*diag(ones(6,1))+diag(ones(4,1),2)+diag(ones(4,1),-2)+diag(ones(2,1),4)+diag(ones(2,1),-4);
X=kron(ones(1,3),element2edges); Y=kron(element2edges,ones(1,3)); Zmass=zeros(size(X)); Zdivdiv=zeros(size(X));
for j = 1:size(element2edges,1)
    orientationlocal=orientation(j,:);    
    area=areas(j);
    
    %divdiv matrix 
    DIVDIVmatrixlocal=ones(3,3)/area;
    Zdivdiv(j,:)=reshape(diag(orientationlocal)*DIVDIVmatrixlocal*diag(orientationlocal),1,9); 
    
    %mass matrix
    I=elements3(j,:); 
    vertices=coordinates(I,:);
    Nhelp=[vertices,vertices,vertices]-[vertices(1,:),vertices(2,:),vertices(3,:);vertices(1,:),vertices(2,:),vertices(3,:);vertices(1,:),vertices(2,:),vertices(3,:)];  
    MASSmatrixlocal=Nhelp*Mhelp*Nhelp'/48/area;
    Zmass(j,:)=reshape(diag(orientationlocal)*MASSmatrixlocal*diag(orientationlocal),1,9);   
end

MASSmatrixRT0=sparse(X,Y,Zmass); 
DIVDIVmatrixRT0=sparse(X,Y,Zdivdiv); 

Contact us