function [M1, M2, F, b, x, e_conn] = assemb_co(param) % The FEM equation for the temp. dist at time t_{n+1} satisfies % (M_1 + M_2) z^{n+1} - M_1 z^n + F + b = 0 % The steady state solution is the solution of % M_2 z^{ss} + F + b = 0 %% Initialization %--------------------------------------------------------------------- x_1 = 0.000; y_1 = 0.000; x_2 = param.L; y_2 = param.W; n_nodesx = param.nodesx; % number of nodes in the x direction n_nodesy = param.nodesy; % number of nodes in the y direction n_gauss = 7; % number of point used in Gaussian integration %% Geometry Module %----------------------------------------------------------------------- % Generate a mesh of the unit square with quadratic elements n_nodes = n_nodesx*n_nodesy; n_elements = 2*(n_nodesx-1)*(n_nodesy-1)/4; x_nodes = linspace(x_1, x_2, n_nodesx); y_nodes = linspace(y_1, y_2, n_nodesy); % % Generate node coordinates. % x = zeros(n_nodes, 2); for i=1:n_nodesx x_i = x_nodes(i); k = i - n_nodesx; for j=1:n_nodesy k = k + n_nodesx; x(k,:) = [ x_i y_nodes(j)]; % k = i+(j-1)*n_nodesx; % x(k,1) = x_nodes(i); % x(k,2) = y_nodes(j); end end % Generate element connectivity e_conn = zeros(n_elements, 6); ie = 0; for j=1:2:n_nodesy-1 for i=1:2:n_nodesx-1 ie = ie + 1; k = i + (j-1)*n_nodesx; e_conn(ie,:) = [k, k+2+2*n_nodesx, k+2*n_nodesx, ... k + 1 + n_nodesx, k+1+2*n_nodesx, k+n_nodesx ]; ie = ie + 1; e_conn(ie,:) = [k, k+2, k+2+2*n_nodesx, ... k+1, k+2+n_nodesx, k+1+n_nodesx ]; end end % % Set up arrays for elements with faces on left/right boundary. % Omega_left = 1:(n_nodesx-1):n_elements; Omega_right = (n_nodesx-1):(n_nodesx-1):n_elements; % % With Neumann/Robins bc all nodes are unknown. % n_equations = n_nodes; % % Comment this line for PMODE, uncomment if for SPMD. % spmd % % Set up codistributed structure % % Column pointers and such for codistributed arrays. % Vc = codistributed.colon(1, n_equations); lP = getLocalPart(Vc); lP_1 = lP(1); lP_end = lP(end); dPM = getCodistributor(Vc.Partition); col_shft = [0 cumsum(dPM(1:end-1))]; % % local sparse arrays % M1_part = sparse(n_equations, dPM(labindex)); M2_part = M1_part; b_part = sparse(dPM(labindex), 1); F_part = b_part; fid = fopen(['out' int2str(labindex)], 'w'); fprintf(fid, 'Lab %d : column range %d %d \n', labindex, lP_1, lP_end); [n_elements, nel_dof] = size(e_conn); [r, s, w] = twod_gauss(n_gauss); % % Build the finite element matrices - Begin loop over elements % n_el_lab = 0; for n_el=1:n_elements nodes_local = e_conn(n_el,:);% which nodes are in this element % subset of nodes/columns on this lab lab_nodes_local = my_extract( nodes_local, lP_1, lP_end); % fprintf(fid, '\n \n n_el = %d \n', n_el); % fprintf(fid, 'lab_nodes_local has %d rows \n', size(lab_nodes_local, 1)); if ~isempty(lab_nodes_local) % continue the calculation for this elmnt n_el_lab = n_el_lab + 1; % % coordinates, shape functions & gradients evaluated at the Gauss pts. % x_local = x(nodes_local,:); [x_g, w_g, phi, p_x, p_y] = twod_shape(x_local, r, s, w); %------------------------------------------------------------------- %% Element contributions to the weak form of the equations % First we evaluate the surface integral (2D) terms %------------------------------------------------------------------- M1_loc = twod_bilinear( 1 , phi, phi, w_g); M2_loc = twod_bilinear(param.k_x, p_x, p_x, w_g) + ... twod_bilinear(param.k_y, p_y, p_y, w_g); src_g = source(x_g, param); % distributed source term F_loc = twod_f_int( src_g , phi, w_g ); % pre-allocate boundary arrays b_loc = zeros(nel_dof,1); %--------------------------------------------------------------------- %% Add boundary terms %--------------------------------------------------------------------- % % Add boundary-integral terms for left and right bounary elements % Left boundary n_el = [1:2*(n_nodesx-1):n_elements % Right boundary n_el = [2*(n_nodesx-1):2*(n_nodesx-1):n_elements % if any(Omega_left == n_el); % % Here we impose Robin boundary term. The node points on the boundary % are nodes_local([3 6 1]) and are at x_local([3 6 1],:) % [phi_g1, ip, x_g1, w_g1] = boundary2(3, x_local, 5); % n_gauss=5 alpha_g1 = steps(x_g1(:,2), param.yh, param.ep, param.alpha ); phi_a = diag(alpha_g1)*phi_g1; % mult. each row by alpha(y_g) beta_g1 = steps(x_g1(:,2), param.yh, param.ep, param.beta ); % % 1D integrals are from y = w to y = 0 % for ii=1:3 for jj=1:3 M2_loc(ip(ii),ip(jj)) = M2_loc(ip(ii),ip(jj)) + ... oned_f_int( phi_g1(:, jj), phi_a(:,ii), w_g1 ); end b_loc(ip(ii),1) = b_loc(ip(ii),1) + ... oned_f_int( beta_g1, phi_a(:,ii), w_g1 ); end elseif any(Omega_right == n_el); % % Here we impose boundary term for a specified flux % The node points on the boundary are nodes_local([2 3 5]) and % are at x_local([2 3 5],:) % We need \int f(y) \phi(L, y) dy along the right boundary. % [phi_g1, ip, x_g1, w_g1] = boundary2(2, x_local, 5); % n_gauss=5; for ii=1:3 b_loc(ip(ii),1) = b_loc(ip(ii),1) + ... oned_f_int( param.flux, phi_g1(:,ii), w_g1 ); end end %----------------------------------------------------------------- %% Assemble contributions into the global system matrices %----------------------------------------------------------------- % for n_t = 1:nel_dof % local DOF - test fcn t_glb = nodes_local(n_t); % global DOF - test fcn for n_u = 1:size(lab_nodes_local, 1) n_locj = lab_nodes_local(n_u, 1); % local DOF in current n_el n_glbj = lab_nodes_local(n_u, 2) ... -col_shft(labindex); % global DOF M1_part(t_glb,n_glbj) = M1_part(t_glb,n_glbj) ... + M1_loc(n_t,n_locj); M2_part(t_glb,n_glbj) = M2_part(t_glb,n_glbj) ... + param.dt*M2_loc(n_t,n_locj); end % if t_glb >= lP_1 && t_glb <= lP_end % why is this needed ? t_loc = t_glb - col_shft(labindex); b_part(t_loc,1) = b_part(t_loc,1) - param.dt*b_loc(n_t,1); F_part(t_loc,1) = F_part(t_loc,1) - param.dt*F_loc(n_t,1); end end % for n_t % fprintf(fid, 'b_loc(1:3) = %8.3f %8.3f %8.3f \n', b_loc(1:3)); % fprintf(fid, 'b_loc(4:6) = %8.3f %8.3f %8.3f \n', b_loc(4:6)); end % if not empty end % n_el fprintf(fid,'Fraction of elements evaluated on this lab is %8.6f \n',... n_el_lab/n_elements); % % Assemble the local parts in a codistributed format. % M1 = codistributed(M1_part, codistributor('1d', 2)); M2 = codistributed(M2_part, codistributor('1d', 2)); b = codistributed(b_part , codistributor('1d', 1)); F = codistributed(F_part , codistributor('1d', 1)); fclose(fid); end % spmd : comment this line for pmode, uncomment for spmd return end