-
Notifications
You must be signed in to change notification settings - Fork 0
/
MFEM_preProcess_Macro.m
73 lines (60 loc) · 2.56 KB
/
MFEM_preProcess_Macro.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
%*********************************************************************
% Preprocess matrices for computing the macro-scale solution
%
% This code is based on:
% Bahriawati, C., & Carstensen, C. (2005).
% Three MATLAB implementations of the lowest-order Raviart-Thomas
% MFEM with a posteriori error control.
% Computational Methods in Applied Mathematics, 5(4), 333-359
%*********************************************************************
%
%***------------------------------------
% Manuela Bastidas - 2020
% Hasselt University, Belgium
function [C,D,geometry] = MFEM_preProcess_Macro(geometry,N_macro)
%%
%*********************************************************************
%* *
%* MFEM - MACRO SOLUTION *
%* *
%*********************************************************************
% Assemble matrices B and C
% B = sparse(geometry.noedges, geometry.noedges);
C = sparse(geometry.noedges,geometry.nElement);
D = sparse(geometry.nElement,geometry.nElement);
geometry.size = 0;
geometry.bari = zeros(geometry.nElement,2);
% Imposition of the source as dirichlet conditions
% This part can be improved
N_real = [1 1];
int = 1/N_macro;
P = [ 0 0;int 0;0 int
N_real(1) N_real(2); N_real(1)-int N_real(2); N_real(1) N_real(2)-int];
TR = triangulation([1 2 3; 4 5 6],P);
ele1 = []; ele2=[];
for j = 1:geometry.nElement
coord = geometry.coordinate(geometry.element(j,:),:)';
I = full(diag(geometry.nodes2edge(geometry.element(j,[2 3 1]),geometry.element(j,[3 1 2]))));
signum = ones(1,3);
signum((j==geometry.edge2element(I,4)))= -1;
geometry.area(j)=polyarea(coord(1,:),coord(2,:));
n = coord(:,[3,1,2])-coord(:,[2,3,1]);
C(I,j) = diag(signum)*[norm(n(:,1)) norm(n(:,2)) norm(n(:,3))]';
D(j,j) = geometry.area(j);
geometry.size = max(geometry.size,max([norm(n(:,1)) norm(n(:,2)) norm(n(:,3))]));
% B(I,I)= B(I,I)+ diag(signum)*...
% stimaB(coord,eye(2))*diag(signum);
geometry.bari(j,:) = sum(coord,2)/3;
%% Source
% THIS PART SHOULD BE IMPROVED!
stbary= cartesianToBarycentric(TR,1,geometry.bari(j,:));
if stbary(1)>0 && stbary(2)>0
ele1 = [ele1;j];
end
stbary= cartesianToBarycentric(TR,2,geometry.bari(j,:));
if stbary(1)>0 && stbary(2)>0
ele2 = [ele2;j];
end
end
geometry.source1 = ele1;
geometry.source2 = ele2;