%main
clear
% Définir les constantes
l = 20; % Remplacez par la valeur appropriée
h = 20; % Remplacez par la valeur appropriée
Qg=1;
k=1;
H=0.1;
Tex=18;
Tin=0;
nddlpN=1; % nombre de ddl par noeud
nddlpEl=4; % nombre de ddl par element
%paramettres pour mesh
meshParams.pointdecontrole=[0 l;0 h];
meshParams.NumElemsXYZ=[2 1];
NumElemsXDir = meshParams.NumElemsXYZ(1,1);
NumElemsYDir = meshParams.NumElemsXYZ(1,2);
nNXDir = meshParams.NumElemsXYZ(1,1) + 1; %nombre de noeuds suivant la
direction x
nNYDir = meshParams.NumElemsXYZ(1,2) + 1; %nombre de noeuds suivant la
direction y
nEl=NumElemsXDir*NumElemsYDir; %nombre total d'elements
nN=nNXDir*nNYDir;
%mesh
coor=[0 0
10 0
20 0
0 20
10 20
20 20];
conn=genererconn2D(meshParams);
% conn=[1 2 3 4
% 3 4 5 6
% 4 5 7 8];
% Définir les variables symboliques
syms e z;
%initialiser les matrices
Q=zeros(nN,1);
KT=zeros(nN);
A=[1;1;1;1];
% Définir les matrices
R = 1/4 * [-1+e, z-1; 1-e, -1-z; 1+e, 1+z; -1-e, 1-z];
S = [2/l, 0; 0, 2/h];
P = S; % P est la même que S dans votre description
Y = 1/4 * [-1+e, 1-e, 1+e, -1-e; z-1, -1-z, 1+z, 1-z];
detJ = (l*h)/4;
Qe = ((l*h)/4)*(Qg+k*H*(Tex-Tin))*A;
% Calculer le produit des quatre matrices
B = R * P * S * Y * detJ;
% Afficher le résultat
%disp('Le produit des quatre matrices est :');
%disp(B);
% Initialiser la matrice pour stocker les résultats des intégrations finales
K = zeros(size(B));
% Définir les points et les poids de la méthode de Gauss à 4 points
gauss_points = [-1/sqrt(3), -1/sqrt(3); 1/sqrt(3), -1/sqrt(3); -1/sqrt(3),
1/sqrt(3); 1/sqrt(3), 1/sqrt(3)];
gauss_weights = [1, 1, 1, 1];
% Calculer l'intégrale pour chaque élément de B
for i = 1:size(B, 1)
for j = 1:size(B, 2)
% Définir la fonction à intégrer pour l'élément B(i, j)
integrand = matlabFunction(B(i, j));
% Initialiser l'intégrale
integral_B = 0;
% Appliquer la méthode de Gauss à 4 points
for k = 1:length(gauss_weights)
e_val = gauss_points(k, 1);
z_val = gauss_points(k, 2);
integral_B = integral_B + gauss_weights(k) * integrand(e_val,
z_val);
end
% Multiplier par l'aire du quadrilatère
%integral_B = integral_B * (1 / sqrt(3))^2;
% Stocker le résultat dans la matrice finale K
K(i, j) = integral_B;
end
end
for i=1:nEl
Q(conn(i,:)) = Q(conn(i,:)) + Qe;
KT(conn(i,:), conn(i,:)) = KT(conn(i,:), conn(i,:)) + K;
end
% Assemblage
% for iel=1:nEl
% nods=conn(iel,:);
% index = globalIndex(nods);
% [Q, KT] = assVecMat(Q,KT,Qe,K,index);
% end
% Afficher le résultat final
disp('La matrice résultante K après intégration par méthode de Gauss à 4
points est :');
disp(K);
disp(Q);
disp(KT);
% % Used Functions
% function index = globalIndex(nods)
% n = length(nods);
% index = zeros(4,1);
% for nod = 1:n
% index(2*nod-1:2*nod) = [2*nods(nod)-1, 2*nods(nod)]; % Ajustement
des indices
% end
% end
% function [Q, kT] = assVecMat(Q,kT,Qe,K,index)
%
%
% % Ajouter Qe aux éléments de Q correspondants
% Q(index) = Q(index) + Qe;
%
% % Ajouter K aux éléments de kT correspondants
% kT(index, index) = kT(index, index) + K;
% end
Fonctions :
genererCoords2D(MeshParams)
function Coord = genererCoords2D(MeshParams)
% generate 2D coordinates for the structured domain
%
%---------------------inputs:
% MeshParams: geometry properties
% MeshParams.ControlPoints: Domain control points,
% format: 2D: [Xmin Xmax;Ymin Ymax]
% MeshParams.NumElemsXYZ: Number of elements in each direction
% in the domain, format: 2D: [XDir YDir]
%
%--------------------outputs:
% FeCoord: FE nodal coordinates [size: (NumFeNodes,2)]
%
%
%% step-1: get the geometry information
ControlPoints = MeshParams.ControlPoints;
NumNodesXDir = MeshParams.NumElemsXYZ(1,1) + 1;
NumNodesYDir = MeshParams.NumElemsXYZ(1,2) + 1;
%% step-2: generate equally distanced points using linspace
x_vals = linspace(ControlPoints(1,1),ControlPoints(1,2),NumNodesXDir);
y_vals = linspace(ControlPoints(2,1),ControlPoints(2,2),NumNodesYDir);
%% step-3: generate grid using meshgrid
[x_mesh,y_mesh] = meshgrid(x_vals,y_vals);
%% step-4: reshape and get the finalized coordinates
Coord = [(reshape(x_mesh',[],1)) (reshape(y_mesh',[],1))];
end
genererconn2D(MeshParams)
function conn = genererconn2D(MeshParams)
% generate 2D connectivity matrix for the structured domain
%
%---------------------inputs:
% MeshParams: geometry properties
% MeshParams.ControlPoints: Domain control points,
% format: 2D: [Xmin Xmax;Ymin Ymax]
% MeshParams.NumElemsXYZ: Number of elements in each direction
% in the domain, format: 2D: [XDir YDir]
%
%--------------------outputs:
% FeTopo: FE elemental connectivity [size: (NumFeElems,4)]
%
%
%% step-1: get the geometry information
NumElemsXDir = MeshParams.NumElemsXYZ(1,1);
NumElemsYDir = MeshParams.NumElemsXYZ(1,2);
NumNodesXDir = NumElemsXDir + 1;
%% step-2: generate nodal local pattern
incr_xdir = 1;
incr_ydir = NumNodesXDir;
node_pattern = [1 2 NumNodesXDir+2 NumNodesXDir+1];
meshParams.pointdecontrole=[0 20;0 20];
meshParams.NumElemsXYZ=[2 2];
NumElemsXDir = meshParams.NumElemsXYZ(1,1);
NumElemsYDir = meshParams.NumElemsXYZ(1,2);
nNXDir = meshParams.NumElemsXYZ(1,1) + 1; %nombre de noeuds suivant la
direction x
nNYDir = meshParams.NumElemsXYZ(1,2) + 1; %nombre de noeuds suivant la
direction y
nEl=NumElemsXDir*NumElemsYDir; %nombre total d'elements
NumNoeudEl=[1 2 3 4];
Coord = genererCoords2D(meshParams);
conn=genererconn2D(meshParams);
NumNE=NumNoeudEl(:,1);
for k=1:nEl
cordxk=Coord(conn(k,NumNE),1);
cordyk=Coord(conn(k,NumNE),2);
CoordXY=[cordxk cordyk];
end
