Add files via upload

FEM/DirectFEM/FEM/ApplyDispPIV.m

@@ -0,0 +1,75 @@
+function [BC,nApplied]=ApplyDispPIV(dx,PIV,Set,X,X0)
+% Applies displacement of PIV data. PIV struct contains:
+% PIV.X(i,:): x annd y positions of PIV node i
+% PIV.U(i,:): x annd y displacement of PIV node i
+% X(i,:) : x coordinates of mesh node i. 3rd dimension is computed.
+nPIV=size(PIV.X,1);
+BCu=zeros(0,3);
+nnodes=1:size(X,1);
+nApplied=0;
+for i=1:nPIV
+    XP=PIV.X(i,:);
+    U=PIV.U(i,:);
+    nodes=nnodes(abs(XP(1)-X(:,1))<Set.PIVtol*dx(1) & abs(XP(2)-X(:,2))<Set.PIVtol*dx(2));
+    BCu=[BCu
+        nodes'   ones(length(nodes),1) ones(length(nodes),1)*U(1)];
+    if Set.ConstrainY
+        BCu=[BCu
+            nodes' 2*ones(length(nodes),1) ones(length(nodes),1)*U(2)];
+    end
+    nApplied=nApplied + length(nodes);
+end
+% Make sure that at least one dof is constrained along Y
+xmin=min(X0(:,1));
+xmax=max(X0(:,1));
+nodes=1:size(X0,1);
+nodeM=nodes(X0(:,1)==xmax);
+nodem=nodes(X0(:,1)==xmin);
+if ~Set.ConstrainY
+    BCu=[BCu
+        nodeM(1) 2 0 
+        nodem(1) 2 0];
+end
+% Make sure that at x=xmin, boundary remains flat (Sensitive to some
+% missing displacements). 
+% Apply mean displacements of measurements for points that x0 < xmin +dx(1)
+xmin=min(X0(:,1));
+nodes=nnodes(X0(:,1)-xmin<dx(1));
+u=0;
+k=0;
+iBC=zeros(size(BCu,1),1);
+for i=1:size(BCu,1)
+    if BCu(i,2)==1
+        if min(abs(BCu(i,1)-nodes))==0
+            u=u+BCu(i,3);
+            k=k+1;
+            iBC(k)=i;
+        end
+    end
+end
+iBC(k+1:end)=[];
+u=u/k;
+BCu(iBC,:)=[BCu(iBC,1) BCu(iBC,2) u*ones(size(iBC))];
+% Add previous displacements if accumulated (u accounts for all accumulated displacements)
+if Set.AccumulateU
+    uPrev=zeros(size(BCu,1),1);
+    for i=1:size(BCu,1)
+        uPrev(i)=X(BCu(i,1),BCu(i,2))-X0(BCu(i,1),BCu(i,2));
+    end
+    BC.u=[BCu(:,1) BCu(:,2) BCu(:,3)+uPrev];
+else
+    BC.u=BCu;
+end
+if nApplied ==0
+    warning('Could not apply displacements from PIV');
+    BC.u=[1 3 0];
+else
+    ymin=min(X0(:,2));
+    ymax=max(X0(:,2));
+    nodes=1:size(X0,1);
+    nodeM=nodes(X0(:,2)==ymax);
+    nodem=nodes(X0(:,2)==ymin);
+    BC.u=[BC.u
+        nodem(1) 3 0 
+        nodeM(1) 3 0];
+end

FEM/DirectFEM/FEM/Elem2Node.m

@@ -0,0 +1,16 @@
+function [Tx,Ty,Tz]=Elem2Node(dx,dy,Tx,Ty,Tz)
+%***************************************************
+% Transforms elemental tractions to nodal values
+% Syntax:
+%   Ae = Aeq4e(Xe,Ge)
+% Input:
+%   Tx,Ty   : elemental values of X and Y tractions
+% Output:
+%   Tx,Ty   : nodal forces in X and Y equivalent to Tx and Ty
+% Date:
+%   Version 1.0    04.05.14
+%***************************************************
+A=dx*dy;
+Tx=Tx*A;
+Ty=Ty*A;
+Tz=Tz*A;

FEM/DirectFEM/FEM/FEM.asv

@@ -0,0 +1,109 @@
+function [aux,En,Sn,Snv,u,uExp]=FEM(BC,C,En,Set,Sn,X,X0)
+%d,dim,E,h,Set,Tx,Ty,Tz,v)
+% Solves Forward problem using FEM.
+% Assumed (fully or partially) constrained problem.
+% INPUT:
+% E v  
+% En(e,gp,:) = Previous strains of element e at gauss  point gp. Exx,Eyy, Eyy, Exy,Exy,Exz,Eyz
+% Sn(e,gp,:) = Previous stresses of element e at gauss point. Sxx,Syy, Syy,Sxy,Sxy,Sxz,Syz
+% OUTPUT:
+% aux = dofs of constrained dofs (applied displacements of PIV)
+% uExp = experimental displacements
+% En(e,gp,:) = Strains of element e at gauss  point gp. Exx,Eyy, Eyy, Exy,Exy,Exz,Eyz
+% Sn(e,gp,:) = Stresses of element e at gauss point. Sxx,Syy, Szz,Sxy,Sxy,Sxz,Syz
+% Snv(e,gp,:) = Viscous Stresses of element e at gauss point. Sxx,Syy, Szz,Sxy,Sxy,Sxz,Syz
+if Set.AccumulateU % In viscoelasticity or accumulated elasticity
+    x=X0;
+else
+    x=X;
+end
+G=[Set.E Set.v];
+nodes=size(x,1);
+dim=3;
+nelem=size(C,1);
+nnod=size(C,2);
+%    A=sparse(nodes*dim,nelemb*dimT);
+%    Ae=Aeq4e(X(C(1,1:nnod),:),Set.Tz); % Matrix such that fe=Ae*[Tx ; Ty; Tz]
+dof=dim*nodes;
+f=zeros(dof,1);
+ig=zeros(nelem*nnod*dim*nnod*dim,1);
+jg=ig;
+vg=ig;
+k=0;
+for e=1:nelem
+    lnod=C(e,1:nnod);
+    if nnod==8
+        [Ke,qe]=keh8e(En(e,:,:),G,Sn(e,:,:),x(lnod,:),Set);
+    elseif nnod==4
+        [Ke,qe]=keq4e(En(e,:,:),G,Sn(e,:,:),x(lnod,:),Set);
+    else
+        error('elements with %i nodes not implemented',nnod)
+    end
+    dofe=kron(lnod-1,ones(1,dim))*dim+kron(ones(1,nnod),1:dim);
+    for i=1:nnod*dim
+        for j=1:nnod*dim
+            if abs(Ke(i,j))>eps
+                k=k+1;
+                ig(k)=dofe(i);
+                jg(k)=dofe(j);
+                vg(k)=Ke(i,j);
+            end
+        end
+    end
+    f(dofe)=f(dofe)-qe;
+end
+Kc=sparse(ig(1:k),jg(1:k),vg(1:k),dof,dof);
+% Apply constraints
+[aux,naux,uaux]=ApplyBC(BC,dim,dof);
+% FORWARD PROBLEM. Compute nodal displacements
+u=zeros(size(Kc,1),1);
+u(aux)=uaux; 
+uExp=ones(size(u))*(min(uaux)+max(uaux))/2; % All set to average imposed values
+uExp(aux)=uaux;
+uExp=reshape(uExp,size(x'))';
+f=f-Kc*u;
+f(aux)=[];
+Kc(aux,:)=[];
+Kc(:,aux)=[];
+us=Solve(Kc,f);
+u(naux)=us;
+u=reshape(u,size(x'))';
+Ep=En; % Previous strains
+Sp=Sn; % Previous stresses
+Snv=0*Sn;
+Een=zeros(size(En,2),size(En,3));
+% Compute Stresses and Strains
+for e=1:nelem
+    lnod=C(e,1:nnod);
+    Ue=u(lnod,:); % Total applied displacements at current time-step
+    Xe=x(lnod,:);
+    Een(:,:)=Ep(e,:,:); % elemental Previous strain
+    Sen(:,:)=Sp(e,:,:); % elemental Previous stresses
+    [~,Ee,Se,Sev] = qeh8e(Een,G,Ue,Sen,Set,Xe);
+    En(e,:,:)=Ee; 
+    Sn(e,:,:)=Se; 
+    Snv(e,:,:)=Sev; 
+end
+end
+%%
+function [aux,naux,uaux]=ApplyBC(BC,dim,dof)
+% OUTPUT:
+% aux = list of constratined dof
+% naux = listof free dof
+% uaxu = valus of displacement in constrainted dof
+ndof=size(BC.u,1);
+uaux=1:dof;
+aux=-uaux; % Constrained dof
+for i=1:ndof
+    node=BC.u(i,1);
+    dof=BC.u(i,2);
+    u=BC.u(i,3);
+    idof=(node-1)*dim+dof;
+    aux(idof)=idof;
+    uaux(idof)=u;
+end
+naux=-aux(aux<0);
+uaux(naux)=[];
+aux(naux)=[];
+end
+

FEM/DirectFEM/FEM/FEM.m

@@ -0,0 +1,113 @@
+function [aux,En,Sn,Snv,u,uExp]=FEM(BC,C,En,Set,Sn,X,X0)
+%d,dim,E,h,Set,Tx,Ty,Tz,v)
+% Solves Forward problem using FEM.
+% Assumed (fully or partially) constrained problem.
+% INPUT:
+% E v
+% En(e,gp,:) = Previous strains of element e at gauss  point gp. Exx,Eyy, Eyy, Exy,Exy,Exz,Eyz
+% Sn(e,gp,:) = Previous stresses of element e at gauss point. Sxx,Syy, Syy,Sxy,Sxy,Sxz,Syz
+% OUTPUT:
+% aux = dofs of constrained dofs (applied displacements of PIV)
+% uExp = experimental displacements
+% En(e,gp,:) = Strains of element e at gauss  point gp. Exx,Eyy, Eyy, Exy,Exy,Exz,Eyz
+% Sn(e,gp,:) = Stresses of element e at gauss point. Sxx,Syy, Szz,Sxy,Sxy,Sxz,Syz
+% Snv(e,gp,:) = Viscous Stresses of element e at gauss point. Sxx,Syy, Szz,Sxy,Sxy,Sxz,Syz
+if Set.AccumulateU % In viscoelasticity or accumulated elasticity
+    x=X0;
+else
+    x=X;
+end
+G=[Set.E Set.v];
+nodes=size(x,1);
+dim=size(X,2);
+nelem=size(C,1);
+nnod=size(C,2);
+%    A=sparse(nodes*dim,nelemb*dimT);
+%    Ae=Aeq4e(X(C(1,1:nnod),:),Set.Tz); % Matrix such that fe=Ae*[Tx ; Ty; Tz]
+dof=dim*nodes;
+f=zeros(dof,1);
+ig=zeros(nelem*nnod*dim*nnod*dim,1);
+jg=ig;
+vg=ig;
+k=0;
+for e=1:nelem
+    lnod=C(e,1:nnod);
+    if nnod==8
+        [Ke,qe]=keh8e(En(e,:,:),G,Sn(e,:,:),x(lnod,:),Set);
+    elseif nnod==4
+        [Ke,qe]=keq4e(En(e,:,:),G,Sn(e,:,:),x(lnod,:),Set);
+    else
+        error('elements with %i nodes not implemented',nnod)
+    end
+    dofe=kron(lnod-1,ones(1,dim))*dim+kron(ones(1,nnod),1:dim);
+    for i=1:nnod*dim
+        for j=1:nnod*dim
+            if abs(Ke(i,j))>eps
+                k=k+1;
+                ig(k)=dofe(i);
+                jg(k)=dofe(j);
+                vg(k)=Ke(i,j);
+            end
+        end
+    end
+    f(dofe)=f(dofe)-qe;
+end
+Kc=sparse(ig(1:k),jg(1:k),vg(1:k),dof,dof);
+% Apply constraints
+[aux,naux,uaux]=ApplyBC(BC,dim,dof);
+% FORWARD PROBLEM. Compute nodal displacements
+u=zeros(size(Kc,1),1);
+u(aux)=uaux;
+uExp=ones(size(u))*(min(uaux)+max(uaux))/2; % All set to average imposed values
+uExp(aux)=uaux;
+uExp=reshape(uExp,size(x'))';
+f=f-Kc*u;
+f(aux)=[];
+Kc(aux,:)=[];
+Kc(:,aux)=[];
+us=Solve(Kc,f);
+u(naux)=us;
+u=reshape(u,size(x'))';
+Ep=En; % Previous strains
+Sp=Sn; % Previous stresses
+Snv=0*Sn;
+Een=zeros(size(En,2),size(En,3));
+% Compute Stresses and Strains
+for e=1:nelem
+    lnod=C(e,1:nnod);
+    Ue=u(lnod,:); % Total applied displacements at current time-step
+    Xe=x(lnod,:);
+    Een(:,:)=Ep(e,:,:); % elemental Previous strain
+    Sen(:,:)=Sp(e,:,:); % elemental Previous stresses
+    if nnod==8
+        [~,Ee,Se,Sev] = qeh8e(Een,G,Ue,Sen,Set,Xe);
+    elseif nnod==4
+        [~,Ee,Se,Sev] = qeq4e(Een,G,Ue,Sen,Set,Xe);
+    end
+    En(e,:,:)=Ee;
+    Sn(e,:,:)=Se;
+    Snv(e,:,:)=Sev;
+end
+end
+%%
+function [aux,naux,uaux]=ApplyBC(BC,dim,dof)
+% OUTPUT:
+% aux = list of constratined dof
+% naux = listof free dof
+% uaxu = valus of displacement in constrainted dof
+ndof=size(BC.u,1);
+uaux=1:dof;
+aux=-uaux; % Constrained dof
+for i=1:ndof
+    node=BC.u(i,1);
+    dof=BC.u(i,2);
+    u=BC.u(i,3);
+    idof=(node-1)*dim+dof;
+    aux(idof)=idof;
+    uaux(idof)=u;
+end
+naux=-aux(aux<0);
+uaux(naux)=[];
+aux(naux)=[];
+end
+

FEM/DirectFEM/FEM/MainArtificalDispPIV.m

@@ -0,0 +1,62 @@
+function resultats=MainArtificalDispPIV(Case)
+% INPUT:
+% Case :
+%   =1: Initial linear (10 stesp) and the constant displacements (10
+%   steps). For creep test.
+%	=2: Constant displacements (20 steps). For stress relaxation test.
+% OUTPUT:
+% resultats.Xlag: x coordinates of applied x displacement
+% resultats.Ylag: y coordinates of applied y displacement
+% resultats.Ulag: x velocities
+% resultats.Vlag: y velocities
+% resultats.time: time points
+% SINTAX:
+% clearvars;resultats=MainArtificalDispPIV(1);save('results_creep.mat')
+% clearvars;resultats=MainArtificalDispPIV(2);save('results_relaxation.mat')
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Settings
+nt=20; % number time points
+Utot=10; % final applied displacement
+Xmin=90;
+Xmax=400;
+Ymin=-30;
+Ymax=30;
+dY=2;
+resultats.time=1:nt;
+Y=Ymin:dY:Ymax;
+nY=length(Y);
+Xlag=zeros(2*nY,nt);
+Ylag=zeros(2*nY,nt);
+Ulag=zeros(2*nY,nt);
+Vlag=zeros(2*nY,nt);
+X1=Xmin*ones(1,nY);
+X2=Xmax*ones(1,nY);
+Zeros=zeros(1,nY);
+Ones=ones(size(Zeros));
+for t=1:nt
+    if Case==1 % linear then constant, for creep testing.
+        Ylag(:,t) = [Y Y];
+        Vlag(:,t) = [Zeros Zeros];
+        if t<=nt/2
+            Xlag(:,t) = [X1 X2+Utot/(nt/2)*(t-1)];
+            Ulag(:,t) = [Zeros Utot/(nt/2)*Ones];
+        else
+            Xlag(:,t) = [X1 X2+Utot];
+            Ulag(:,t) = [Zeros Zeros];
+        end
+    elseif Case==2 % all constant, for stress relaxation.
+        Ylag(:,t) = [Y Y];
+        Vlag(:,t) = [Zeros Zeros];
+        if t==1
+            Xlag(:,t) = [X1 X2];
+            Ulag(:,t) = [Zeros Utot*Ones];
+        else
+            Xlag(:,t) = [X1 X2+Utot];
+            Ulag(:,t) = [Zeros Zeros];
+        end
+    end
+end
+resultats.Xlag=Xlag;
+resultats.Ylag=Ylag;
+resultats.Ulag=Ulag;
+resultats.Vlag=Vlag;