Add files via upload

BuildBCMaps2D.m

@@ -0,0 +1,23 @@
+function BuildBCMaps2D()
+
+% function BuildMaps2DBC
+% Purpose: Build specialized nodal maps for various types of
+%          boundary conditions, specified in BCType. 
+
+Globals2D;
+
+% create label of face nodes with boundary types from BCType
+bct    = BCType';
+bnodes = ones(Nfp, 1)*bct(:)';
+bnodes = bnodes(:);
+
+% find location of boundary nodes in face and volume node lists
+mapI = find(bnodes==In);           vmapI = vmapM(mapI);
+mapO = find(bnodes==Out);          vmapO = vmapM(mapO);
+mapW = find(bnodes==Wall);         vmapW = vmapM(mapW);
+mapF = find(bnodes==Far);          vmapF = vmapM(mapF);
+mapC = find(bnodes==Cyl);          vmapC = vmapM(mapC);
+mapD = find(bnodes==Dirichlet);    vmapD = vmapM(mapD);
+mapN = find(bnodes==Neuman);       vmapN = vmapM(mapN);
+mapS = find(bnodes==Slip);         vmapS = vmapM(mapS);
+return;

BuildCurvedOPS2D.m

@@ -0,0 +1,81 @@
+function [cinfo] = BuildCurvedOPS2D(intN)
+
+% function [cinfo] = BuildCurvedOPS2D(intN)
+% Purpose: build curved info for curvilinear elements
+
+Globals2D;
+
+% 1. Create cubature information
+
+% 1.1 Extract cubature nodes and weights
+[cR,cS,cW,Ncub] = Cubature2D(intN);
+
+% 1.1. Build interpolation matrix (nodes->cubature nodes)
+cV = InterpMatrix2D(cR, cS);
+
+% 1.2 Evaluate derivatives of Lagrange interpolants at cubature nodes
+[cDr,cDs] = Dmatrices2D(N, cR, cS, V);
+
+% 2. Create surface quadrature information
+
+% 2.1 Compute Gauss nodes and weights for 1D integrals
+[gz, gw] = JacobiGQ(0, 0, intN);
+
+% 2.2 Build Gauss nodes running counter-clockwise on element faces
+gR = [gz, -gz, -ones(size(gz))];
+gS = [-ones(size(gz)), gz, -gz];
+
+% 2.3 For each face 
+for f1=1:Nfaces
+  % 2.3.1 build nodes->Gauss quadrature interpolation and differentiation matrices
+  gV(:,:,f1) = InterpMatrix2D(gR(:,f1), gS(:,f1));
+  [gDr(:,:,f1),gDs(:,:,f1)] = Dmatrices2D(N, gR(:,f1), gS(:,f1), V);
+end
+
+% 3. For each curved element, evaluate custom operator matrices
+Ncurved = length(curved);
+
+% 3.1 Store custom information in array of Matlab structs
+cinfo = [];
+for c=1:Ncurved 
+  % find next curved element and the coordinates of its nodes
+  k1 = curved(c); x1 = x(:,k1); y1 = y(:,k1); cinfo(c).elmt = k1;
+
+  % compute geometric factors
+  [crx,csx,cry,csy,cJ] = GeometricFactors2D(x1,y1,cDr,cDs);
+
+  % build mass matrix
+  cMM = cV'*diag(cJ.*cW)*cV; cinfo(c).MM = cMM;
+
+  % build physical derivative matrices
+  cinfo(c).Dx = cMM\(cV'*diag(cW.*cJ)*(diag(crx)*cDr + diag(csx)*cDs));
+  cinfo(c).Dy = cMM\(cV'*diag(cW.*cJ)*(diag(cry)*cDr + diag(csy)*cDs));
+
+  % build individual lift matrices at each face
+  for f1=1:Nfaces
+    k2 = EToE(k1,f1); f2 = EToF(k1,f1);
+
+    % compute geometric factors
+    [grx,gsx,gry,gsy,gJ] = GeometricFactors2D(x1,y1,gDr(:,:,f1),gDs(:,:,f1));
+
+    % compute normals and surface Jacobian at Gauss points on face f1
+    if(f1==1) gnx = -gsx;     gny = -gsy;     end;
+    if(f1==2) gnx =  grx+gsx; gny =  gry+gsy; end;
+    if(f1==3) gnx = -grx;     gny = -gry;     end;
+
+    gsJ = sqrt(gnx.*gnx+gny.*gny);
+    gnx = gnx./gsJ;  gny = gny./gsJ;  gsJ = gsJ.*gJ;
+
+    % store normals and coordinates at Gauss nodes
+    cinfo(c).gnx(:,f1) = gnx;  cinfo(c).gx(:,f1)  = gV(:,:,f1)*x1;
+    cinfo(c).gny(:,f1) = gny;  cinfo(c).gy(:,f1)  = gV(:,:,f1)*y1;
+
+    % store Vandermondes for '-' and '+' traces
+    cinfo(c).gVM(:,:,f1) = gV(:,:,f1);
+    cinfo(c).gVP(:,:,f1) = gV(end:-1:1,:,f2);
+
+    % compute and store matrix to lift Gauss node data
+    cinfo(c).glift(:,:,f1) = cMM\(gV(:,:,f1)'*diag(gw.*gsJ));
+  end
+end
+return

BuildHNonCon2D.m

@@ -0,0 +1,96 @@
+function [neighbors] = BuildHNonCon2D(NGauss, tol)
+
+% function [neighbors] = BuildHNonCon2D(NGauss, tol)
+% Purpose: find element to element connections through non-conforming interfaces
+%          (** elements assumed straight sided **)
+
+Globals2D;
+
+% 1. Build Gauss nodes
+[gz, gw] = JacobiGQ(0, 0, NGauss-1);
+
+% 1.1 Find location of vertices of boundary faces
+vx1 = VX(EToV(:,[1,2,3])); vx2 = VX(EToV(:,[2,3,1]));
+vy1 = VY(EToV(:,[1,2,3])); vy2 = VY(EToV(:,[2,3,1]));
+
+idB = find(EToE==((1:K)'*ones(1,Nfaces)));
+x1 = vx1(idB)'; y1 = vy1(idB)';
+x2 = vx2(idB)'; y2 = vy2(idB)';
+
+% 1.2 Find those element-faces that are on boundary faces
+[elmtsB,facesB] = find(EToE==((1:K)'*ones(1,Nfaces)));
+Nbc = length(elmtsB);
+
+sk = 1;
+% 2.1 For each boundary face
+for b1=1:Nbc 
+  % 2.2 Find element and face of this boundary face
+  k1 = elmtsB(b1); f1 = facesB(b1);
+
+  % 2.3 Find end coordinates of b1'th boundary face
+  x11 = x1(b1);  y11 = y1(b1);  x12 = x2(b1);  y12 = y2(b1);
+
+  % 2.4 Compute areas, lengths and face coordinates used in intersection 
+  % tests comparing b1'th boundary face with all boundary faces
+  area1 = abs((x12-x11)*(y1-y11) - (y12-y11)*(x1-x11)); %scale
+  area2 = abs((x12-x11)*(y2-y11) - (y12-y11)*(x2-x11));
+  L   = (x12-x11)^2 + (y12-y11)^2 ; 
+  r21 = ((2*x1-x11-x12)*(x12-x11) + (2*y1-y11-y12)*(y12-y11))/L;
+  r22 = ((2*x2-x11-x12)*(x12-x11) + (2*y2-y11-y12)*(y12-y11))/L;
+
+  % 2.5 Find range of local face coordinate (bracketed between -1 and 1)
+  r1 = max(-1,min(r21,r22)); r2 = min(1,max(r21,r22));
+
+  % 2.6 Compute flag for overlap of b1 face with all other boundary faces
+  flag = area1+area2+(r1<= -1 & r2<= -1)+(r1>=1 & r2>=1)+(r2-r1<tol);
+
+  % 2.7 Find other faces with partial matches
+  matches = setdiff(find(flag < tol),b1); 
+  Nmatches = length(matches(:));
+
+  if(Nmatches>0)
+    % 3.1 Find matches
+    r1 = r1(matches); r2 = r2(matches); 
+
+    % 3.2 Find end points of boundary-boundary intersections
+    xy11 = 0.5*[x11;y11]*(1-r1) +  0.5*[x12;y12]*(1+r1);
+    xy12 = 0.5*[x11;y11]*(1-r2) +  0.5*[x12;y12]*(1+r2);
+
+    % 3.3 For each face-face match
+    for n=1:Nmatches
+
+      % 3.4 Store which elements intersect
+      k2 = elmtsB(matches(n)); f2 = facesB(matches(n));
+      neighbors{sk}.elmtM = k1; neighbors{sk}.faceM = f1;
+      neighbors{sk}.elmtP = k2; neighbors{sk}.faceP = f2;
+
+      % 3.5 Build physical Gauss nodes on face fragment
+      xg = 0.5*(1-gz)*xy11(1,n) + 0.5*(1+gz)*xy12(1,n);
+      yg = 0.5*(1-gz)*xy11(2,n) + 0.5*(1+gz)*xy12(2,n);
+
+      % 3.6 Find local coordinates of Gauss nodes
+      [rg1,sg1] = FindLocalCoords2D(k1, xg, yg);
+      [rg2,sg2] = FindLocalCoords2D(k2, xg, yg);
+
+      % 3.7 Build interpolation matrices for volume nodes ->Gauss nodes
+      gVM = InterpMatrix2D(rg1,sg1); neighbors{sk}.gVM  = gVM;
+      gVP = InterpMatrix2D(rg2,sg2); neighbors{sk}.gVP  = gVP;
+
+      % 3.8 Find face normal 
+      neighbors{sk}.nx = nx(1+(f1-1)*Nfp,k1);
+      neighbors{sk}.ny = ny(1+(f1-1)*Nfp,k1);
+
+      % 4.0 Build partial face data lift operator
+
+      % 4.1 Compute weights for lifting
+      partsJ = sqrt( (xy11(1,n)-xy12(1,n))^2 + (xy11(2,n)-xy12(2,n))^2 )/2;
+      dgw = gw*partsJ/J(1,k1); 
+
+      % 4.2 Build matrix to lift Gauss data to volume data
+      neighbors{sk}.lift = V*V'*(gVM')*diag(dgw);
+
+      sk = sk+1;
+    end
+  end
+end
+return

BuildMaps2D.m

@@ -0,0 +1,45 @@
+function [mapM, mapP, vmapM, vmapP, vmapB, mapB] = BuildMaps2D()
+
+% function [mapM, mapP, vmapM, vmapP, vmapB, mapB] = BuildMaps2D
+% Purpose: Connectivity and boundary tables in the K # of Np elements
+
+Globals2D;
+
+% number volume nodes consecutively
+nodeids = reshape(1:K*Np, Np, K);
+vmapM   = zeros(Nfp, Nfaces, K); vmapP   = zeros(Nfp, Nfaces, K); 
+mapM    = (1:K*Nfp*Nfaces)';     mapP    = reshape(mapM, Nfp, Nfaces, K);
+
+% find index of face nodes with respect to volume node ordering
+for k1=1:K
+  for f1=1:Nfaces
+    vmapM(:,f1,k1) = nodeids(Fmask(:,f1), k1);
+  end
+end
+
+one = ones(1, Nfp);
+for k1=1:K
+  for f1=1:Nfaces
+    % find neighbor
+    k2 = EToE(k1,f1); f2 = EToF(k1,f1);
+
+    % reference length of edge
+    v1 = EToV(k1,f1); v2 = EToV(k1, 1+mod(f1,Nfaces));
+    refd = sqrt( (VX(v1)-VX(v2))^2 + (VY(v1)-VY(v2))^2 );
+
+    % find find volume node numbers of left and right nodes 
+    vidM = vmapM(:,f1,k1); vidP = vmapM(:,f2,k2);    
+    x1 = x(vidM); y1 = y(vidM); x2 = x(vidP); y2 = y(vidP);
+    x1 = x1*one;  y1 = y1*one;  x2 = x2*one;  y2 = y2*one;
+
+    % Compute distance matrix
+    D = (x1 -x2').^2 + (y1-y2').^2;
+    [idM, idP] = find(sqrt(abs(D))<NODETOL*refd);
+    vmapP(idM,f1,k1) = vidP(idP); mapP(idM,f1,k1) = idP + (f2-1)*Nfp+(k2-1)*Nfaces*Nfp;
+  end
+end
+
+% reshape vmapM and vmapP to be vectors and create boundary node list
+vmapP = vmapP(:); vmapM = vmapM(:); mapP = mapP(:); 
+mapB = find(vmapP==vmapM); vmapB = vmapM(mapB);
+return

BuildNonCon2D.m

@@ -0,0 +1,96 @@
+function [neighbors] = BuildNonCon2D(NGauss, tol)
+
+% function [neighbors] = BuildNonCon2D(NGauss, tol)
+% Purpose: find element to element connections through non-conforming interfaces
+%          (** elements assumed straight sided **)
+
+Globals2D;
+
+% 1. Build Gauss nodes
+[gz, gw] = JacobiGQ(0, 0, NGauss-1);
+
+% 1.1 Find location of vertices of boundary faces
+vx1 = VX(EToV(:,[1,2,3])); vx2 = VX(EToV(:,[2,3,1]));
+vy1 = VY(EToV(:,[1,2,3])); vy2 = VY(EToV(:,[2,3,1]));
+
+idB = find(EToE==((1:K)'*ones(1,Nfaces)));
+x1 = vx1(idB)'; y1 = vy1(idB)';
+x2 = vx2(idB)'; y2 = vy2(idB)';
+
+% 1.2 Find those element-faces that are on boundary faces
+[elmtsB,facesB] = find(EToE==((1:K)'*ones(1,Nfaces)));
+Nbc = length(elmtsB);
+
+sk = 1;
+% 2.1 For each boundary face
+for b1=1:Nbc 
+  % 2.2 Find element and face of this boundary face
+  k1 = elmtsB(b1); f1 = facesB(b1);
+
+  % 2.3 Find end coordinates of b1'th boundary face
+  x11 = x1(b1);  y11 = y1(b1);  x12 = x2(b1);  y12 = y2(b1);
+
+  % 2.4 Compute areas, lengths and face coordinates used in intersection 
+  % tests comparing b1'th boundary face with all boundary faces
+  area1 = abs((x12-x11)*(y1-y11) - (y12-y11)*(x1-x11)); %scale
+  area2 = abs((x12-x11)*(y2-y11) - (y12-y11)*(x2-x11));
+  L   = (x12-x11)^2 + (y12-y11)^2 ; 
+  r21 = ((2*x1-x11-x12)*(x12-x11) + (2*y1-y11-y12)*(y12-y11))/L;
+  r22 = ((2*x2-x11-x12)*(x12-x11) + (2*y2-y11-y12)*(y12-y11))/L;
+
+  % 2.5 Find range of local face coordinate (bracketed between -1 and 1)
+  r1 = max(-1,min(r21,r22)); r2 = min(1,max(r21,r22));
+
+  % 2.6 Compute flag for overlap of b1 face with all other boundary faces
+  flag = area1+area2+(r1<= -1 & r2<= -1)+(r1>=1 & r2>=1)+(r2-r1<tol);
+
+  % 2.7 Find other faces with partial matches
+  matches = setdiff(find(flag < tol),b1); 
+  Nmatches = length(matches(:));
+
+  if(Nmatches>0)
+    % 3.1 Find matches
+    r1 = r1(matches); r2 = r2(matches); 
+
+    % 3.2 Find end points of boundary-boundary intersections
+    xy11 = 0.5*[x11;y11]*(1-r1) +  0.5*[x12;y12]*(1+r1);
+    xy12 = 0.5*[x11;y11]*(1-r2) +  0.5*[x12;y12]*(1+r2);
+
+    % 3.3 For each face-face match
+    for n=1:Nmatches
+
+      % 3.4 Store which elements intersect
+      k2 = elmtsB(matches(n)); f2 = facesB(matches(n));
+      neighbors{sk}.elmtM = k1; neighbors{sk}.faceM = f1;
+      neighbors{sk}.elmtP = k2; neighbors{sk}.faceP = f2;
+
+      % 3.5 Build physical Gauss nodes on face fragment
+      xg = 0.5*(1-gz)*xy11(1,n) + 0.5*(1+gz)*xy12(1,n);
+      yg = 0.5*(1-gz)*xy11(2,n) + 0.5*(1+gz)*xy12(2,n);
+
+      % 3.6 Find local coordinates of Gauss nodes
+      [rg1,sg1] = FindLocalCoords2D(k1, xg, yg);
+      [rg2,sg2] = FindLocalCoords2D(k2, xg, yg);
+
+      % 3.7 Build interpolation matrices for volume nodes ->Gauss nodes
+      gVM = InterpMatrix2D(rg1,sg1); neighbors{sk}.gVM  = gVM;
+      gVP = InterpMatrix2D(rg2,sg2); neighbors{sk}.gVP  = gVP;
+
+      % 3.8 Find face normal 
+      neighbors{sk}.nx = nx(1+(f1-1)*Nfp,k1);
+      neighbors{sk}.ny = ny(1+(f1-1)*Nfp,k1);
+
+      % 4.0 Build partial face data lift operator
+
+      % 4.1 Compute weights for lifting
+      partsJ = sqrt( (xy11(1,n)-xy12(1,n))^2 + (xy11(2,n)-xy12(2,n))^2 )/2;
+      dgw = gw*partsJ/J(1,k1); 
+
+      % 4.2 Build matrix to lift Gauss data to volume data
+      neighbors{sk}.lift = V*V'*(gVM')*diag(dgw);
+
+      sk = sk+1;
+    end
+  end
+end
+return