Commit 6209ecc9 authored by Turnhout, M.C. van's avatar Turnhout, M.C. van
Browse files

add flux cdu/dx to fem1dcd (issue #8)

parent a561c715
......@@ -30,8 +30,8 @@ f = mat(imat,3);
etop=top(ielem,:);
% position of the local degrees of freedom in the global solution array
[vpos,vshp] = elm1d_i(1,[],etop,mat,5);
[n,dndxi,intw,intcoord] = elm1d_s([],[],etop,mat,vshp(1,:));
[vpos,vshp] = elm1dcd_i(1,[],etop,mat,5);
[n,dndxi,intw,intcoord] = elm1dcd_s([],[],etop,mat,vshp(1,:));
nint=length(intw);
ndofu=length(nonzeros(vpos));
......@@ -46,13 +46,13 @@ unod=sol(pos(ielem,:));
% stresses at the integration points
sigmaint=zeros(2,1);
for int=1:nint
for int = 1:nint
% get shape function derivatives with respect to global coordinates
[dndxy,detj]=getderiv(dndxi,nodcoord,int);
[dndxy, detj] = getderiv(dndxi, nodcoord, int);
% compute the strain displacement matrix
b=dndxy(:,1)';
b = dndxy(:, 1)';
% stress (c du/dx) at this integration point
sigmaint(int) = c*b*unod;
......
......@@ -33,17 +33,18 @@ M = K; % global mass matrix
C = K; % global damping matrix
rhs = zeros(ndof, 1); % global right-hand-side column
% assemble q and rhs
for ielem = 1:nelem
% compute the element stiffness matrix and right-hand-side column
[Me, Ce, Ke, rhse] = elm1dcd(ielem, coord, top, mat.mat);
% assemble the global stiffness matrix and the right-hand-side colomn
ii = pos(ielem, :);
K(ii, ii) = K(ii, ii) + Ke;
M(ii, ii) = M(ii, ii) + Me;
C(ii, ii) = C(ii, ii) + Ce;
rhs(ii) = rhs(ii) + rhse;
end
......@@ -54,18 +55,41 @@ if istat == 1
% solve the system of equations, stationairy problem
q = K + C;
sol = solvestat(q, rhs, bndcon, [], dest);
% calculate flux (c du/dx)
sigma = zeros(size(top, 1), size(top, 2)-2);
for ielem = 1:size(top, 1)
sigma(ielem, :) = elm1dcd_d(ielem, coord, top, mat.mat, ...
pos, sol);
end
elseif istat == 2
% unsteady solution
q = M/dt + theta*C + theta*K;
soln = sol(:, 1);
% calculate flux (c du/dx)
sigma = zeros(size(top, 1), size(top, 2)-2, ntime + 1);
for ielem = 1:size(top, 1)
sigma(ielem, :, 1) = elm1dcd_d(ielem, coord, top, ...
mat.mat, pos, soln);
end
for itime = 1:ntime
fprintf('** Start of time step %i\n', itime);
rhsn = rhs + (M/dt - (1 - theta)*(C + K))*soln;
sol(:, itime + 1) = solvestat(q, rhsn, bndcon, [], dest);
soln = sol(:, itime + 1);
% calculate flux (c du/dx)
for ielem = 1:size(top, 1)
sigma(ielem, :, itime + 1) = elm1dcd_d(ielem, coord, top, ...
mat.mat, pos, soln);
end
end
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment