-
Notifications
You must be signed in to change notification settings - Fork 0
/
hifi_model.m
88 lines (68 loc) · 2.37 KB
/
hifi_model.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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
function [Mesh,Pf,NumNodes,k,mask] = hifi_model(h)
templatefile = 'MeshTemplate';
%% GMSH OPTIONS
options.gmsh = 'gmsh'; %GMSH
options.delete = false;
options.meshopts = '-1 -2 -order 1 -format msh2';
options.meshfile = 'DarcyBlockfine2.msh';
%% Create and Import Mesh
Mesh = CreateGMSHfromTemplate(templatefile,options,'#h#',h);
%% Coordinates of the Inlet Boundary
Coord = GetCoord(Mesh,'Inlet');
CoordY=Coord(:,2);
[Y,Ind]=sort(CoordY);
%% Connectivity of the domain
CONECT{1}=GetConnectivity(Mesh,'Omega1');
CONECT{2}=GetConnectivity(Mesh,'Omega2');
CONECT{3}=GetConnectivity(Mesh,'Omega3');
CONECT{4}=GetConnectivity(Mesh,'Omega4');
CONECT{5}=GetConnectivity(Mesh,'Omega5');
CONECT{6}=GetConnectivity(Mesh,'Omega6');
CONECT{7}=GetConnectivity(Mesh,'Omega7');
CONECT{8}=GetConnectivity(Mesh,'Omega8');
CONECT{9}=GetConnectivity(Mesh,'Omega9');
%%% Connectivity of the boundaries
CInlet = GetConnectivity(Mesh,'Inlet');
CWall = GetConnectivity(Mesh,'Wall');
COutlet = GetConnectivity(Mesh,'Outlet');
%% Boundary Nodes
NumNodes= size(Mesh.X,1);
dofsDir = GetDofs(Mesh,'Outlet');
dofsNB1= GetDofs(Mesh,'Inlet');
dofsNB2= GetDofs(Mesh,'Wall');
%for k=1:size(Permeability,1)
ValNB1=ones(length(dofsNB1),1);% Neumann Boundary Value
%end
%% Permeability
k1=[CONECT{1}, ones(size(CONECT{1},1),1)];
k2=[CONECT{2},ones(size(CONECT{2},1),1)];
k3=[CONECT{3}, ones(size(CONECT{3},1),1)];
k4=[CONECT{4}, ones(size(CONECT{4},1),1)];
k5=[CONECT{5}, 0.0005*ones(size(CONECT{5},1),1)];
k6=[CONECT{6}, ones(size(CONECT{6},1),1)];
k7= [CONECT{7}, ones(size(CONECT{7},1),1)];
k8=[CONECT{8}, ones(size(CONECT{8},1),1)];
k9= [CONECT{9}, ones(size(CONECT{9},1),1)];
k=[k1;k2;k3;k4;k5;k6;k7;k8;k9];
%% Assembly of stiffness and mass matrices and R.H.S flux vector
mask=true(NumNodes,1);
mask(dofsDir)=false; %%% Eliminate the nodes containing the Dirichlet Boundary Conditions
%%%% Assemble matrices for each sub domain
for i=1:9
Kx{i} = Matrix2Dv5(Mesh.X(:,1:2),CONECT{i},1,0,1,0);
Ky{i} = Matrix2Dv5(Mesh.X(:,1:2),CONECT{i},0,1,0,1);
K{i}=Kx{i}+Ky{i};
Pf.K{i}=K{i}(mask,mask);
%M{i}= Matrix2Dv5(Mesh.X(:,1:2),CONECT{i},0,0,0,0);
%P.M{i}=M{i}(mask,mask);
end
%% Mass Matrix on the inlet boundary
M = FEM_mat_1D(Y,0,0);
%% Flux vector
f=zeros(NumNodes,1);
aux=M*ValNB1;
aux2=zeros(size(aux));
aux2(Ind)= aux;
f(dofsNB1)=aux2; %% Apply Neumann BC
Pf.rhs = f(mask) ;
end