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

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);

```