Interior DoF for Cubic Quadrilateral Elements Fixed and currently implemented for Quadrant 1 Mesh

Add dimension checks in the Get and Set methods of IntegrationPoint

CHANGELOG

@@ -13,6 +13,13 @@ Version 3.4.1 (development)
 - Added support for reading linear and quadratic 2D quadrilateral and triangular
   Cubit meshes.
 
+- The tetrahedral mesh refinement algorithm in serial and in parallel now
+  follows precisely the paper:
+     D. Arnold, A. Mukherjee, and L. Pouly, "Locally Adapted Tetrahedral Meshes
+     Using Bisection", SIAM J. Sci. Comput., 22(2), 431–448.
+  This guarantees that the shape regularity of the elements will be preserved
+  under refinement.
+
 
 Version 3.4, released on May 29, 2018
 =====================================

data/mixed-ele-meshes/annulus-trans-24.mesh

@@ -0,0 +1,218 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+74
+2 3 0 1 2 3
+2 3 1 5 6 2
+2 3 5 8 9 6
+2 3 8 11 12 9
+2 3 11 14 15 12
+2 3 14 17 18 15
+2 3 17 20 21 18
+2 3 20 23 24 21
+2 3 23 26 27 24
+2 3 26 29 30 27
+2 3 29 32 33 30
+2 3 32 35 36 33
+2 3 35 38 39 36
+2 3 38 41 42 39
+2 3 41 44 45 42
+2 3 44 47 48 45
+2 3 47 50 51 48
+2 3 50 53 54 51
+2 3 53 56 57 54
+2 3 56 59 60 57
+2 3 59 62 63 60
+2 3 62 65 66 63
+2 3 65 68 69 66
+2 3 68 71 72 69
+2 3 71 74 75 72
+1 2 2 3 4
+1 2 6 2 7
+1 2 9 6 10
+1 2 12 9 13
+1 2 15 12 16
+1 2 18 15 19
+1 2 21 18 22
+1 2 24 21 25
+1 2 27 24 28
+1 2 30 27 31
+1 2 33 30 34
+1 2 36 33 37
+1 2 39 36 40
+1 2 42 39 43
+1 2 45 42 46
+1 2 48 45 49
+1 2 51 48 52
+1 2 54 51 55
+1 2 57 54 58
+1 2 60 57 61
+1 2 63 60 64
+1 2 66 63 67
+1 2 69 66 70
+1 2 72 69 73
+1 2 75 72 76
+1 2 2 4 7
+1 2 6 7 10
+1 2 9 10 13
+1 2 12 13 16
+1 2 15 16 19
+1 2 18 19 22
+1 2 21 22 25
+1 2 24 25 28
+1 2 27 28 31
+1 2 30 31 34
+1 2 33 34 37
+1 2 36 37 40
+1 2 39 40 43
+1 2 42 43 46
+1 2 45 46 49
+1 2 48 49 52
+1 2 51 52 55
+1 2 54 55 58
+1 2 57 58 61
+1 2 60 61 64
+1 2 63 64 67
+1 2 66 67 70
+1 2 69 70 73
+1 2 72 73 76
+
+boundary
+53
+1 1 0 1
+1 1 1 5
+1 1 5 8
+1 1 8 11
+1 1 11 14
+1 1 14 17
+1 1 17 20
+1 1 20 23
+1 1 23 26
+1 1 26 29
+1 1 29 32
+1 1 32 35
+1 1 35 38
+1 1 38 41
+1 1 41 44
+1 1 44 47
+1 1 47 50
+1 1 50 53
+1 1 53 56
+1 1 56 59
+1 1 59 62
+1 1 62 65
+1 1 65 68
+1 1 68 71
+1 1 71 74
+1 1 74 75
+1 1 75 76
+1 1 76 73
+1 1 73 70
+1 1 70 67
+1 1 67 64
+1 1 64 61
+1 1 61 58
+1 1 58 55
+1 1 55 52
+1 1 52 49
+1 1 49 46
+1 1 46 43
+1 1 43 40
+1 1 40 37
+1 1 37 34
+1 1 34 31
+1 1 31 28
+1 1 28 25
+1 1 25 22
+1 1 22 19
+1 1 19 16
+1 1 16 13
+1 1 13 10
+1 1 10 7
+1 1 7 4
+1 1 4 3
+1 1 3 0
+
+vertices
+77
+2
+3.9788735773 0.0
+3.84329674785 1.02980825986
+2.88247256089 0.772356194895
+2.98415518297 0.0
+1.97241688113 0.259673608685
+3.44580559639 1.98943678865
+2.58435419729 1.49207759149
+1.83799993026 0.761324498753
+2.81348848799 2.81348848799
+2.11011636599 2.11011636599
+1.57832632157 1.21109238238
+1.98943678865 3.44580559639
+1.49207759149 2.58435419729
+1.21109238238 1.57832632157
+1.02980825986 3.84329674785
+0.772356194895 2.88247256089
+0.761324498753 1.83799993026
+2.43635739532e-16 3.9788735773
+1.82726804649e-16 2.98415518297
+0.259673608685 1.97241688113
+-1.02980825986 3.84329674785
+-0.772356194895 2.88247256089
+-0.259673608685 1.97241688113
+-1.98943678865 3.44580559639
+-1.49207759149 2.58435419729
+-0.761324498753 1.83799993026
+-2.81348848799 2.81348848799
+-2.11011636599 2.11011636599
+-1.21109238238 1.57832632157
+-3.44580559639 1.98943678865
+-2.58435419729 1.49207759149
+-1.57832632157 1.21109238238
+-3.84329674785 1.02980825986
+-2.88247256089 0.772356194895
+-1.83799993026 0.761324498753
+-3.9788735773 4.87271479065e-16
+-2.98415518297 3.65453609299e-16
+-1.97241688113 0.259673608685
+-3.84329674785 -1.02980825986
+-2.88247256089 -0.772356194895
+-1.97241688113 -0.259673608685
+-3.44580559639 -1.98943678865
+-2.58435419729 -1.49207759149
+-1.83799993026 -0.761324498753
+-2.81348848799 -2.81348848799
+-2.11011636599 -2.11011636599
+-1.57832632157 -1.21109238238
+-1.98943678865 -3.44580559639
+-1.49207759149 -2.58435419729
+-1.21109238238 -1.57832632157
+-1.02980825986 -3.84329674785
+-0.772356194895 -2.88247256089
+-0.761324498753 -1.83799993026
+-7.30907218597e-16 -3.9788735773
+-5.48180413948e-16 -2.98415518297
+-0.259673608685 -1.97241688113
+1.02980825986 -3.84329674785
+0.772356194895 -2.88247256089
+0.259673608685 -1.97241688113
+1.98943678865 -3.44580559639
+1.49207759149 -2.58435419729
+0.761324498753 -1.83799993026
+2.81348848799 -2.81348848799
+2.11011636599 -2.11011636599
+1.21109238238 -1.57832632157
+3.44580559639 -1.98943678865
+2.58435419729 -1.49207759149
+1.57832632157 -1.21109238238
+3.84329674785 -1.02980825986
+2.88247256089 -0.772356194895
+1.83799993026 -0.761324498753
+3.9788735773 -9.7454295813e-16
+2.98415518297 -7.30907218597e-16
+1.97241688113 -0.259673608685
+3.84329674785 1.02980825986
+2.88247256089 0.772356194895
+1.97241688113 0.259673608685

data/mixed-ele-meshes/annulus-trans-6.mesh

@@ -0,0 +1,74 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+20
+2 3 0 1 2 3
+2 3 1 5 6 2
+2 3 5 8 9 6
+2 3 8 11 12 9
+2 3 11 14 15 12
+2 3 14 17 18 15
+2 3 17 20 21 18
+1 2 2 3 4
+1 2 6 2 7
+1 2 9 6 10
+1 2 12 9 13
+1 2 15 12 16
+1 2 18 15 19
+1 2 21 18 22
+1 2 2 4 7
+1 2 6 7 10
+1 2 9 10 13
+1 2 12 13 16
+1 2 15 16 19
+1 2 18 19 22
+
+boundary
+17
+1 1 0 1
+1 1 1 5
+1 1 5 8
+1 1 8 11
+1 1 11 14
+1 1 14 17
+1 1 17 20
+1 1 20 21
+1 1 21 22
+1 1 22 19
+1 1 19 16
+1 1 16 13
+1 1 13 10
+1 1 10 7
+1 1 7 4
+1 1 4 3
+1 1 3 0
+
+vertices
+23
+2
+1.11408460164 0.0
+0.557042300822 0.964825566988
+0.417781725616 0.723619175241
+0.835563451232 0.0
+0.482412783494 0.278521150411
+-0.557042300822 0.964825566988
+-0.417781725616 0.723619175241
+3.41090035345e-17 0.557042300822
+-1.11408460164 1.36436014138e-16
+-0.835563451232 1.02327010604e-16
+-0.482412783494 0.278521150411
+-0.557042300822 -0.964825566988
+-0.417781725616 -0.723619175241
+-0.482412783494 -0.278521150411
+0.557042300822 -0.964825566988
+0.417781725616 -0.723619175241
+-1.02327010604e-16 -0.557042300822
+1.11408460164 -2.72872028276e-16
+0.835563451232 -2.04654021207e-16
+0.482412783494 -0.278521150411
+0.557042300822 0.964825566988
+0.417781725616 0.723619175241
+0.482412783494 0.278521150411

data/mixed-ele-meshes/circInSquare.py

@@ -0,0 +1,924 @@
+#Title:circInSquare.py
+#Author:T. M. McManus
+#Date:10-7-18
+#Purpose: Fill a circular sector with triangles and a bounding region,
+#defined by 3 nodes, with quads.  Then reflect/preserve QuadI twice to
+#create a complete disc bounded in a square.
+
+import scipy as sp
+import argparse
+import sys
+import subprocess
+import time
+
+parser=argparse.ArgumentParser(description='Fill a circular sector with triangles and a bounding region,\
+defined by 3 nodes, with quads.  Then reflect/preserve QuadI twice to create a complete disc bounded in a square.'
+,epilog='Sample run: python circInSquare.py -r 1 -e 2 -n 8 -g ../../../glvis/glvis')
+
+parser.add_argument('-r','--circRad', nargs='?',const=1, default = 1.0, type=float, help='Radius of circle')
+parser.add_argument('-e','--edgeLength', nargs='?',const=1,default=2.0,type=float,help='Edge-length of bounding square')
+parser.add_argument('-n','--numEdges',nargs='?',const=1,default=6,type=int,help='n-gon approximation of internal circle')
+parser.add_argument('-o','--outputFile',nargs='?',const=1,default='circInSquare', help='Output file name.')
+parser.add_argument('-g','--glvis',nargs='?',const=1,default='',type=str,help='Abs. or rel. path of glvis binary.')
+args=parser.parse_args()
+
+r=args.circRad 
+edgeLength=args.edgeLength
+numEdges=args.numEdges
+outputName=args.outputFile
+glvis=args.glvis
+
+visMesh=False;
+
+if glvis!='':
+    visMesh=True
+
+if r >= edgeLength:
+    print("Circle radius must be less than bounding square edge length")
+    sys.exit(1)
+
+if sp.mod(numEdges,2) != 0:
+    print("Currently this mixed element generator only supports an even numbers of edges.")
+    sys.exit(1)
+
+    
+#The basic idea:
+#1. Construct topology for regions
+#2. Combine topologies
+#3. Construct boundary
+#4. Construct geometry for regions
+#5. Combine geometries
+#6. Output
+
+def eleMatCirc(numEdges):
+    
+    nNodesSeq=sp.zeros([numEdges])
+    nNodesSeq[0]=3
+    if numEdges != 1:
+        for n in range(1,numEdges):
+            nNodesSeq[n]=nNodesSeq[n-1]+(2+n)
+            
+    numCircNodesTot =int(((numEdges+1)*(numEdges+2))/2)
+    
+    b=range(numCircNodesTot)
+    row_size=1
+    A=sp.zeros([numEdges+1,numEdges+1])
+    start=0;stop=1;
+    for m in range(numEdges+1):
+        if m==0:
+            A[m,range(m+1)]=b[0:1]
+            start=0
+            stop=1
+        else:
+            start=stop
+            stop=stop+m+1
+            A[m,range(m+1)]=b[start:stop]
+
+    M=sp.ones([numEdges**2,5])
+    m_row=0
+    for m in range(numEdges):
+        if m==0:
+            M[0,:]=[1,2,0,1,2]
+            m_row+=1
+        else:
+           holder=sp.size(sp.nonzero(A[m,:]))
+           for n in range(holder):
+               if n!=holder-1:
+                   M[m_row,:]=[1,2,A[m,n],A[m,n+1],A[m+1,n+1]]
+                   m_row+=1
+                   M[m_row,:]=[1,2,A[m,n],A[m+1,n],A[m+1,n+1]]
+                   m_row+=1
+               else:
+                   M[m_row,:]=[1,2,A[m,n],A[m+1,n],A[m+1,n+1]]
+                   m_row+=1
+            
+    return M.astype(int),numCircNodesTot
+
+def eleMatQuad(numEdges):
+    S0=numEdges*(numEdges+1)/(2.0)
+    A=sp.linspace(S0,(S0+(numEdges+1)**2)-1,(numEdges+1)**2)
+    A=A.reshape([numEdges+1,numEdges+1])
+    quadNode=sp.delete(A,-1,1)
+    quadNode=sp.delete(quadNode,-1,0)
+    quadNode=quadNode.flatten()
+    M=sp.zeros([numEdges**2,6])
+    for n in range(numEdges**2):
+        M[n,:]=[2,3,quadNode[n],quadNode[n]+1,quadNode[n]+numEdges+2,quadNode[n]+numEdges+1]
+    return M.astype(int)
+
+def boundMatTot(numEdges):
+    triS1=sp.zeros(numEdges+1)
+    triS3=sp.zeros(numEdges+1)
+    quadS1=sp.zeros(numEdges)
+    quadS2=sp.zeros(numEdges-1)
+    quadS3=sp.zeros(numEdges)
+
+    triS1[0]=0;
+    triS3[0]=0;
+    for n in range(1,numEdges+1):
+        triS1[n]=triS1[n-1]+n
+        triS3[n]=triS1[n]+n
+    ref1=triS3
+    
+    triS3=sp.flipud(triS3)
+    quadS1[0]=triS1[-1]+numEdges+1
+    quadS3[0]=triS1[-1]+2*numEdges+1
+
+    for n in range(1,numEdges):
+        quadS1[n]=quadS1[n-1]+(numEdges+1)
+        quadS3[n]=quadS3[n-1]+(numEdges+1)    
+    ref2=quadS3
+    xAxisRootRef=sp.concatenate([triS1.copy(),quadS1],axis=0)
+    quadS3=sp.flipud(quadS3)
+    quadS2=range(int(quadS1[-1]+1),int(quadS3[0]),1)
+    STOT=sp.concatenate([triS1,quadS1,quadS2,quadS3,triS3],axis=0)
+
+    filler=sp.zeros(1)
+    filler[0]=quadS3[0]
+    fillerFirst=sp.zeros(1)
+    fillerFirst[0]=quadS1[-1]
+    sTotRef=sp.concatenate([triS1,quadS1,quadS2,filler],axis=0)
+    newsTotRef=sp.concatenate([fillerFirst,quadS2,filler],axis=0)
+    boundMat=sp.zeros([STOT.size-1,4])
+    boundMatRef=sp.zeros([sTotRef.size-1,4])
+    new_boundMat_ref=sp.zeros([newsTotRef.size-1,4])
+    
+    for n in range(STOT.size-1):
+        boundMat[n,:]=[1,1,STOT[n],STOT[n+1]]
+    for n in range(sTotRef.size-1):
+        boundMatRef[n,:]=[1,1,sTotRef[n],sTotRef[n+1]]
+    for n in range(newsTotRef.size-1):
+        new_boundMat_ref[n,:]=[1,1,newsTotRef[n],newsTotRef[n+1]]
+        
+    ref=sp.concatenate([ref1,ref2],axis=0).astype(int)
+    return boundMat.astype(int),ref,boundMatRef.astype(int),xAxisRootRef.astype(int),new_boundMat_ref.astype(int)
+
+def vertMatCirc(numEdges):
+    r_o=sp.linspace(0,r,numEdges+1)
+    counter=0
+    vertMat=sp.zeros([numCircNodesTot,2])
+    for m in range(numEdges+1):
+        theta=sp.linspace(0,sp.pi/2.0,m+1)
+        for n in range(sp.size(theta)):
+            vertMat[counter,:]=[r_o[m]*sp.cos(theta[n]),r_o[m]*sp.sin(theta[n])]
+            counter+=1
+    return vertMat
+
+def vertMatQuad(numEdges):
+
+    theta=sp.linspace(0,sp.pi/2.0,numEdges+1)
+    AX=sp.zeros([numEdges+1,numEdges+1])
+    AY=sp.zeros([numEdges+1,numEdges+1])
+    AX[0,:]=r*sp.cos(theta)
+    AY[0,:]=r*sp.sin(theta)
+
+    vertLinSpace=sp.linspace(0,edgeLength,(numEdges/2)+1)
+    horzLineSpace=sp.linspace(edgeLength,0,(numEdges/2)+1)
+
+    #Assigning node locations along the boundary
+    vertCount=0
+    horzCount=1
+    for n in range(numEdges+1):
+        if n < (numEdges/2):
+            AX[-1,n]=edgeLength
+            AY[-1,n]=vertLinSpace[vertCount]
+            vertCount+=1
+        elif n == int(numEdges/2):
+            AX[-1,n]=edgeLength
+            AY[-1,n]=edgeLength
+        else:
+            AX[-1,n]=horzLineSpace[horzCount]
+            AY[-1,n]=edgeLength
+            horzCount+=1
+
+    #Linearly spacing nodes between the inner/outer boundaries
+    #One could then smooth this via r-based adaptivity
+    for col in range(numEdges+1):
+        for row in range(1,numEdges):
+            AX[row,col]=sp.linspace(AX[0,col],AX[-1,col],numEdges+1)[row]
+            AY[row,col]=sp.linspace(AY[0,col],AY[-1,col],numEdges+1)[row]
+
+    AX=sp.delete(AX,0,0)
+    AY=sp.delete(AY,0,0)
+    AX=AX.flatten()
+    AY=AY.flatten()
+    AX_reshape = AX.flatten()
+    
+    numQuadNodesTot=numEdges*(numEdges+1)
+    vertMat=sp.zeros([numQuadNodesTot,2])
+    for n in range(numQuadNodesTot):
+        vertMat[n,:]=[AX[n],AY[n]]    
+    return vertMat
+
+def orient(A):
+    aOrient=sp.zeros([A.shape[0],A.shape[1]])
+    triCounter=0
+    quadCounter=0
+    #Determine the number of triangle and quad elments in the given element matrix
+    for n in range(A.shape[0]):
+        if A[n,1]==2:
+            triCounter+=1
+        else:
+            quadCounter+=1
+    edgeMatTotal=sp.zeros([3*triCounter+4*quadCounter,2])
+    counter=0
+    for n in range(A.shape[0]):
+        detected=0
+        if A[n,1]==2:
+            for m in range(edgeMatTotal.shape[0]):
+                if detected != 1:
+                    if edgeMatTotal[m,0]==A[n,2] and edgeMatTotal[m,1]==A[n,3]:
+                        aOrient[n,:]=[1,2,A[n,2],A[n,4],A[n,3],0]
+                        detected=1
+                        #print("reorder:[{} {} {}] to [{} {} {}]".format(A[n,2],A[n,3],A[n,4],int(aOrient[n,2]),int(aOrient[n,3]),int(aOrient[n,4])))
+                    elif edgeMatTotal[m,0]==A[n,4] and edgeMatTotal[m,1]==A[n,2]:
+                        aOrient[n,:]=[1,2,A[n,2],A[n,4],A[n,3],0]
+                        detected=1
+                    else:
+                        aOrient[n,:]=A[n,:]
+
+            edgeMatTotal[counter,:]=[aOrient[n,2],aOrient[n,3]]
+            counter+=1
+            edgeMatTotal[counter,:]=[aOrient[n,3],aOrient[n,4]]
+            counter+=1
+            edgeMatTotal[counter,:]=[aOrient[n,4],aOrient[n,2]]
+            counter+=1
+        else:
+            for m in range(edgeMatTotal.shape[0]):
+                if detected != 1: 
+                    if edgeMatTotal[m,0]==A[n,2] and edgeMatTotal[m,1]==A[n,3]:
+                        aOrient[n,:]=[2,3,A[n,2],A[n,5],A[n,4],A[n,3]]
+                        detected=1
+                        #print("reorder:[{} {} {} {}] to [{} {} {} {}]".format(A[n,2],A[n,3],A[n,4],A[n,5],int(aOrient[n,2]),int(aOrient[n,3]),int(aOrient[n,4]),int(aOrient[n,5])))
+                    elif edgeMatTotal[m,0]==A[n,5] and edgeMatTotal[m,1]==A[n,2]:
+                        aOrient[n,:]=[2,3,A[n,2],A[n,5],A[n,4],A[n,3]]
+                        detected=1
+                    else:
+                        aOrient[n,:]=A[n,:]
+            edgeMatTotal[counter,:]=[aOrient[n,2],aOrient[n,3]]
+            counter+=1
+            edgeMatTotal[counter,:]=[aOrient[n,3],aOrient[n,4]]
+            counter+=1
+            edgeMatTotal[counter,:]=[aOrient[n,4],aOrient[n,5]]
+            counter+=1
+            edgeMatTotal[counter,:]=[aOrient[n,5],aOrient[n,2]]
+            counter+=1
+    return aOrient.astype(int)
+
+def gVis(_glvis,_meshFile):
+
+    if(_glvis==''):
+        print("Failure: Set glvis location via -g switch")
+        sys.exit(1)
+
+    colFuncFileName=_meshFile.replace('.mesh','.gf')
+    glvsScriptFileName=_meshFile.replace('.mesh','.glvs')
+    imageFileName=_meshFile.replace('.mesh','.png')
+    
+    #Create Coloring Function for mesh
+    _colFuncCommand=_glvis+ ' -m '+ _meshFile +' -sc -k q'
+    args=_colFuncCommand.split()
+    p=subprocess.Popen(args)#Create 'GLVis_coloring.gf'
+
+    _renameCommand='mv GLVis_coloring.gf {}'.format(colFuncFileName)
+    args=_renameCommand.split()
+    p=subprocess.Popen(args)
+
+    #Glvis script template
+    f=open(glvsScriptFileName,'w')
+    f.write('window 0 0 800 800\n'+'\n')
+    f.write('solution {} {}\n'.format(_meshFile,colFuncFileName)+'\n')
+    f.write('{\n'+'perspective off\n'+'zoom 1.5\n'+'keys gAeeRM\n'+'solution {} {} screenshot {}\n'.format(_meshFile,colFuncFileName,imageFileName)+'keys q\n'+'}\n')
+    f.close()
+
+    _runGlvisCommand=_glvis+' -run {}'.format(glvsScriptFileName)
+    args=_runGlvisCommand.split()
+    p=subprocess.Popen(args)
+    p.wait()
+    
+    return 0
+def quadInterDof(_edge,_linEleMat,_linVertMatRound):
+    _state=False
+    for n in range(_linEleMat.shape[0]):
+        if _linEleMat[n,1]==3:
+            if sp.any(_edge[0]==_linEleMat[n,2:6]) and sp.any(_edge[1]==_linEleMat[n,2:6]):
+                print("{} is possibly in {}".format(_edge,_linEleMat[n,2:6]))
+                _n1Loc=sp.where(_edge[0]==_linEleMat[n,2:6])[0][0]
+                _n2Loc=sp.where(_edge[1]==_linEleMat[n,2:6])[0][0]
+                if _n1Loc==sp.mod(_n2Loc+1,4) or _n1Loc==sp.mod(_n2Loc-1,4):
+                    _state=True
+                    xcent=(_linVertMatRound[_linEleMat[n,2],0]+_linVertMatRound[_linEleMat[n,3],0]+_linVertMatRound[_linEleMat[n,4],0]+_linVertMatRound[_linEleMat[n,5],0])/4.0
+                    ycent=(_linVertMatRound[_linEleMat[n,2],1]+_linVertMatRound[_linEleMat[n,3],1]+_linVertMatRound[_linEleMat[n,4],1]+_linVertMatRound[_linEleMat[n,5],1])/4.0
+                    _interDof=sp.zeros(2)
+                    _interDof[0]=sp.round_((_linVertMatRound[_edge[0],0]+_linVertMatRound[_edge[1],0]+xcent)/3.0,5)
+                    _interDof[1]=sp.round_((_linVertMatRound[_edge[0],1]+_linVertMatRound[_edge[1],1]+ycent)/3.0,5)
+                    print("dof loc is {},{}".format(_interDof[0],_interDof[1]))
+                    return(_state,_interDof[0],_interDof[1])
+
+    return(_state,0,0)
+
+[eleMatTriHolder,numCircNodesTot]=eleMatCirc(numEdges) #Construct tri element matrix for the region inside circular sector
+eleMatQuadHolder=eleMatQuad(numEdges) #Construct quad element matrix for region outside the circular sector
+
+#Combining eleMatTriHolder and eleMatQuadHolder
+linEleMat=sp.zeros([eleMatTriHolder.shape[0]+eleMatQuadHolder.shape[0],6])
+counter=0
+for n in range(eleMatTriHolder.shape[0]):
+    linEleMat[n,[0,1,2,3,4]]=eleMatTriHolder[n,:]
+    counter+=1
+for n in range(eleMatQuadHolder.shape[0]):
+    linEleMat[counter+n,:]=eleMatQuadHolder[n,:]
+
+linEleMat=linEleMat.astype(int)
+linBoundMat=boundMatTot(numEdges)[0] #Construct the boundary
+vertMatCircHolder = vertMatCirc(numEdges) #Construct vertex matrix for triang region
+vertMatQuadHolder = vertMatQuad(numEdges) #Construct vertex matrix for the quad region
+
+
+#Combining the two vertex matrices in Quadrant I (q1)
+linVertMat=sp.zeros([vertMatCircHolder.shape[0]+vertMatQuadHolder.shape[0],2])
+counter=0
+for n in range(vertMatCircHolder.shape[0]):
+    linVertMat[n,:]=vertMatCircHolder[n,:]
+    counter+=1
+for n in range(vertMatQuadHolder.shape[0]):
+    linVertMat[counter+n,:]=vertMatQuadHolder[n,:]
+
+#Outputting P1/Q1 mesh to a .mesh file
+g=open(outputName+'Lin.mesh','w')
+g.write('MFEM mesh v1.0\n'+'\n')
+g.write('dimension\n'+'2\n'+'\n')
+g.write('elements\n'+'{}\n'.format(linEleMat.shape[0]))
+for n in range(linEleMat.shape[0]):
+    if linEleMat[n,1]==2:
+        g.write('{} {} {} {} {}\n'.format(linEleMat[n,0],linEleMat[n,1],linEleMat[n,2],linEleMat[n,3],linEleMat[n,4]))
+    else:
+        g.write('{} {} {} {} {} {}\n'.format(linEleMat[n,0],linEleMat[n,1],linEleMat[n,2],linEleMat[n,3],linEleMat[n,4],linEleMat[n,5]))
+g.write('\n'+'boundary\n'+'{}\n'.format(linBoundMat.shape[0]))
+for n in range(linBoundMat.shape[0]):
+    g.write('{} {} {} {}\n'.format(linBoundMat[n,0],linBoundMat[n,1],linBoundMat[n,2],linBoundMat[n,3]))
+g.write('\n'+'vertices\n'+'{}\n'.format(linVertMat.shape[0])+'2\n')
+for n in range(linVertMat.shape[0]):
+    g.write('{} {}\n'.format(linVertMat[n,0],linVertMat[n,1]))
+g.close()
+
+if(visMesh==True):
+    gVis(glvis,outputName+'Lin.mesh')
+
+#Quadratic (P2/Q2) Element Generation
+
+#1.)Create Edge list from previously generated linear elements
+edgeMat=sp.zeros([3*eleMatTriHolder.shape[0]+4*eleMatQuadHolder.shape[0],2])
+linEleMat=orient(linEleMat)#Make sure that element orientation is in agreement with MFEM requirements
+
+counter=0
+for n in range(linEleMat.shape[0]):
+    if linEleMat[n,1]==2:
+        edgeMat[counter,:]=[linEleMat[n,2],linEleMat[n,3]]
+        counter+=1
+        edgeMat[counter,:]=[linEleMat[n,3],linEleMat[n,4]]
+        counter+=1
+        edgeMat[counter,:]=[linEleMat[n,4],linEleMat[n,2]]
+        counter+=1
+    else:
+        edgeMat[counter,:]=[linEleMat[n,2],linEleMat[n,3]]
+        counter+=1
+        edgeMat[counter,:]=[linEleMat[n,3],linEleMat[n,4]]
+        counter+=1
+        edgeMat[counter,:]=[linEleMat[n,4],linEleMat[n,5]]
+        counter+=1
+        edgeMat[counter,:]=[linEleMat[n,5],linEleMat[n,2]]
+        counter+=1
+
+#Remove duplicates
+holder=[]
+for n in range(edgeMat.shape[0]):
+    counter=0
+    for m in range(edgeMat.shape[0]):
+        if edgeMat[n,0]==edgeMat[m,0] and edgeMat[n,1]==edgeMat[m,1] and m!=n:
+            holder.append([n,m])
+        elif edgeMat[n,1]==edgeMat[m,0] and edgeMat[n,0]==edgeMat[m,1] and m!=n:
+            holder.append([n,m])
+
+removeIndices=sp.zeros(len(holder))
+for n in range(len(holder)):
+    if holder[n][0]>holder[n][1]:
+        removeIndices[n]=holder[n][0]
+    else:
+        removeIndices[n]=holder[n][1]
+removeIndices=sp.unique(removeIndices).astype(int)
+edgeMat=sp.delete(edgeMat,removeIndices,0)
+edgeMat=edgeMat.astype(int)
+
+edgeDofMat=sp.zeros([edgeMat.shape[0],2])#These will be the new DoFs that appear after the Element Vertices within the .mesh file
+linVertMatRound=sp.round_(linVertMat,5) 
+ 
+counter=0
+
+for n in edgeMat:
+    if linVertMatRound[n[0],1] == linVertMatRound[n[1],1]:
+        xmid=(linVertMatRound[n[0],0]+linVertMatRound[n[1],0])/2.0
+        ymid=linVertMatRound[n[0],1]
+        edgeDofMat[counter,:]=[xmid,ymid]
+    elif linVertMatRound[n[0],0] == linVertMatRound[n[1],0]:
+        xmid=linVertMatRound[n[0],0]
+        ymid=(linVertMatRound[n[0],1]+linVertMatRound[n[1],1])/2.0
+        edgeDofMat[counter,:]=[xmid,ymid]
+    else:
+        r0=sp.sqrt(linVertMatRound[n[0],0]**2+linVertMatRound[n[0],1]**2)
+        r1=sp.sqrt(linVertMatRound[n[1],0]**2+linVertMatRound[n[1],1]**2)
+        rmid = (r0+r1)/2.0 #should not be needed
+        xmidOld=(linVertMatRound[n[0],0]+linVertMatRound[n[1],0])/2.0
+        ymidOld=(linVertMatRound[n[0],1]+linVertMatRound[n[1],1])/2.0
+        midtheta=sp.arctan(ymidOld/xmidOld)
+        xmid=rmid*sp.cos(midtheta)
+        ymid=rmid*sp.sin(midtheta)
+        edgeDofMat[counter,:]=[xmid,ymid]
+    counter+=1
+edgeDofMat = sp.round_(edgeDofMat,5)
+
+#Determine midpoints of all Q1 elements:
+quadCentroidLoc=sp.zeros([eleMatQuadHolder.shape[0],2])
+for n in range(eleMatQuadHolder.shape[0]):
+    quadCentroidLoc[n,0]=(linVertMatRound[eleMatQuadHolder[n,2],0]+linVertMatRound[eleMatQuadHolder[n,3],0]+linVertMatRound[eleMatQuadHolder[n,4],0]+linVertMatRound[eleMatQuadHolder[n,5],0])/4.0
+    quadCentroidLoc[n,1]=(linVertMatRound[eleMatQuadHolder[n,2],1]+linVertMatRound[eleMatQuadHolder[n,3],1]+linVertMatRound[eleMatQuadHolder[n,4],1]+linVertMatRound[eleMatQuadHolder[n,5],1])/4.0
+
+quadCentroidLoc = sp.round_(quadCentroidLoc,5)
+
+#3.)Populate nodes section
+g=open(outputName+'Quad.mesh','w')
+g.write('MFEM mesh v1.0\n'+'\n')
+g.write('dimension\n'+'2\n'+'\n')
+g.write('elements\n'+'{}\n'.format(linEleMat.shape[0]))
+for n in range(linEleMat.shape[0]):
+    if linEleMat[n,1]==2:
+        g.write('{} {} {} {} {}\n'.format(linEleMat[n,0],linEleMat[n,1],linEleMat[n,2],linEleMat[n,3],linEleMat[n,4]))
+    else:
+        g.write('{} {} {} {} {} {}\n'.format(linEleMat[n,0],linEleMat[n,1],linEleMat[n,2],linEleMat[n,3],linEleMat[n,4],linEleMat[n,5]))
+g.write('\n'+'boundary\n'+'{}\n'.format(linBoundMat.shape[0]))
+for n in range(linBoundMat.shape[0]):
+    g.write('{} {} {} {}\n'.format(linBoundMat[n,0],linBoundMat[n,1],linBoundMat[n,2],linBoundMat[n,3]))
+g.write('\n'+'vertices\n'+'{}\n'.format(linVertMat.shape[0]))
+g.write('\n'+'nodes'+'\n'+'FiniteElementSpace'+'\n'+'FiniteElementCollection: H1_2D_P2'+'\n'+'VDim: 2'+'\n'+'Ordering: 1' +'\n\n')
+for n in range(linVertMatRound.shape[0]):
+    g.write('{} {}\n'.format(linVertMatRound[n,0],linVertMatRound[n,1]))
+for n in range(edgeDofMat.shape[0]):
+    g.write('{} {}\n'.format(edgeDofMat[n,0],edgeDofMat[n,1]))
+for n in range(quadCentroidLoc.shape[0]):
+    g.write('{} {}\n'.format(quadCentroidLoc[n,0],quadCentroidLoc[n,1]))
+g.close()
+
+if(visMesh==True):
+    gVis(glvis,outputName+'Quad.mesh')
+
+#Cubic (P3/Q3) Element Generation
+
+cubeDofMat=sp.zeros([2*edgeMat.shape[0],2])#These will be the new DoFs that appear after the Element Vertices within the .mesh file
+
+counter=0
+for n in edgeMat: #Here DoF ordering matters. 
+    if linVertMatRound[n[0],1] == linVertMatRound[n[1],1]:
+        xmid=(linVertMatRound[n[0],0]+linVertMatRound[n[1],0])/2.0
+        ymid=linVertMatRound[n[0],1]
+        xmid1=(linVertMatRound[n[0],0]+xmid)/2.0
+        ymid1=linVertMatRound[n[0],1]
+        xmid2=(linVertMatRound[n[1],0]+xmid)/2.0
+        ymid2=linVertMatRound[n[0],1]
+        if n[0] > n[1]:
+            cubeDofMat[counter,:]=[xmid2,ymid2]
+            counter+=1
+            cubeDofMat[counter,:]=[xmid1,ymid1]
+            counter+=1
+        else:
+            cubeDofMat[counter,:]=[xmid1,ymid1]
+            counter+=1
+            cubeDofMat[counter,:]=[xmid2,ymid2]
+            counter+=1
+            
+    elif linVertMatRound[n[0],0] == linVertMatRound[n[1],0]:
+        xmid=linVertMatRound[n[0],0]
+        ymid=(linVertMatRound[n[0],1]+linVertMatRound[n[1],1])/2.0
+        xmid1=linVertMatRound[n[0],0]
+        ymid1=(linVertMatRound[n[0],1]+ymid)/2.0
+        xmid2=linVertMatRound[n[0],0]
+        ymid2=(linVertMatRound[n[1],1]+ymid)/2.0
+        if n[0] > n[1]:
+            cubeDofMat[counter,:]=[xmid2,ymid2]
+            counter+=1
+            cubeDofMat[counter,:]=[xmid1,ymid1]
+            counter+=1
+        else:
+            cubeDofMat[counter,:]=[xmid1,ymid1]
+            counter+=1
+            cubeDofMat[counter,:]=[xmid2,ymid2]
+            counter+=1
+    else:
+        r0=sp.sqrt(linVertMatRound[n[0],0]**2+linVertMatRound[n[0],1]**2)
+        r1=sp.sqrt(linVertMatRound[n[1],0]**2+linVertMatRound[n[1],1]**2)
+        rmid = (r0+r1)/2.0 #should not be needed
+        xmidOld=(linVertMatRound[n[0],0]+linVertMatRound[n[1],0])/2.0
+        ymidOld=(linVertMatRound[n[0],1]+linVertMatRound[n[1],1])/2.0
+        midtheta=sp.arctan(ymidOld/xmidOld)
+        xmid=rmid*sp.cos(midtheta)
+        ymid=rmid*sp.sin(midtheta)
+        xmid1=(linVertMatRound[n[0],0]+xmid)/2.0
+        ymid1=(linVertMatRound[n[0],1]+ymid)/2.0
+        xmid2=(linVertMatRound[n[1],0]+xmid)/2.0
+        ymid2=(linVertMatRound[n[1],1]+ymid)/2.0
+        if n[0] > n[1]:
+            cubeDofMat[counter,:]=[xmid2,ymid2]
+            counter+=1
+            cubeDofMat[counter,:]=[xmid1,ymid1]
+            counter+=1
+        else:
+            cubeDofMat[counter,:]=[xmid1,ymid1]
+            counter+=1
+            cubeDofMat[counter,:]=[xmid2,ymid2]
+            counter+=1
+            
+cubeDofMat = sp.round_(cubeDofMat,5)
+
+triCentroidLoc=sp.zeros([eleMatTriHolder.shape[0],2])
+
+for n in range(eleMatTriHolder.shape[0]):
+    triCentroidLoc[n,0]=(linVertMatRound[eleMatTriHolder[n,2],0]+linVertMatRound[eleMatTriHolder[n,3],0]+linVertMatRound[eleMatTriHolder[n,4],0])/3.0
+    triCentroidLoc[n,1]=(linVertMatRound[eleMatTriHolder[n,2],1]+linVertMatRound[eleMatTriHolder[n,3],1]+linVertMatRound[eleMatTriHolder[n,4],1])/3.0
+
+quadCentroidLocCubic=sp.zeros([4*eleMatQuadHolder.shape[0],2])
+
+counter=0
+for n in range(eleMatQuadHolder.shape[0]):
+    xcent=quadCentroidLoc[n,0];ycent=quadCentroidLoc[n,1]
+    a=eleMatQuadHolder[n,2:6]
+    aMinIndex=sp.where(a[:]==a.min())[0][0]
+    dof0=0.5*sp.array([xcent+linVertMatRound[a[aMinIndex],0],ycent+linVertMatRound[a[aMinIndex],1]])
+    quadCentroidLocCubic[counter,:]=dof0
+    counter+=1
+    if aMinIndex==0:
+        aLeft=-1
+        aRight=1
+        aLast=2
+    else:
+        aLeft=aMinIndex-1
+        aRight=aMinIndex+1
+        aLast=sp.delete(a,[aMinIndex,aLeft,aRight])[0]
+    edge1=[a[aMinIndex], a[aLeft]]
+    edge2=[a[aMinIndex], a[aRight]]
+    edge1Index=0
+    edge2Index=0
+    edgeCounter=0
+    for edge in edgeMat:
+        if(edge[0]==edge1[0] and edge[1]==edge1[1]) or (edge[1]==edge1[0] and edge[0]==edge1[1]):
+            edge1Index=edgeCounter
+        if(edge[0]==edge2[0] and edge[1]==edge2[1]) or (edge[1]==edge2[0] and edge[0]==edge2[1]):
+            edge2Index=edgeCounter
+        edgeCounter+=1
+
+    if (edge1Index > edge2Index):
+        dof1=0.5*sp.array([xcent+linVertMatRound[a[aLeft],0],ycent+linVertMatRound[a[aLeft],1]])
+        quadCentroidLocCubic[counter,:]=dof1
+        counter+=1
+        dof2=0.5*sp.array([xcent+linVertMatRound[a[aRight],0],ycent+linVertMatRound[a[aRight],1]])
+        quadCentroidLocCubic[counter,:]=dof2
+        counter+=1
+        dof3=0.5*sp.array([xcent+linVertMatRound[a[aLast],0],ycent+linVertMatRound[a[aLast],1]])
+        quadCentroidLocCubic[counter,:]=dof3
+        counter+=1
+    else:
+        dof1=0.5*sp.array([xcent+linVertMatRound[a[aRight],0],ycent+linVertMatRound[a[aRight],1]])
+        quadCentroidLocCubic[counter,:]=dof1
+        counter+=1
+        dof2=0.5*sp.array([xcent+linVertMatRound[a[aLeft],0],ycent+linVertMatRound[a[aLeft],1]])
+        quadCentroidLocCubic[counter,:]=dof2
+        counter+=1
+        dof3=0.5*sp.array([xcent+linVertMatRound[a[aLast],0],ycent+linVertMatRound[a[aLast],1]])
+        quadCentroidLocCubic[counter,:]=dof3
+        counter+=1
+        
+truCentroidLoc=sp.round_(triCentroidLoc,5)
+    
+#3.)Populate nodes section
+g=open(outputName+'Cub.mesh','w')
+g.write('MFEM mesh v1.0\n'+'\n')
+g.write('dimension\n'+'2\n'+'\n')
+g.write('elements\n'+'{}\n'.format(linEleMat.shape[0]))
+for n in range(linEleMat.shape[0]):
+    if linEleMat[n,1]==2:
+        g.write('{} {} {} {} {}\n'.format(linEleMat[n,0],linEleMat[n,1],linEleMat[n,2],linEleMat[n,3],linEleMat[n,4]))
+    else:
+        g.write('{} {} {} {} {} {}\n'.format(linEleMat[n,0],linEleMat[n,1],linEleMat[n,2],linEleMat[n,3],linEleMat[n,4],linEleMat[n,5]))
+g.write('\n'+'boundary\n'+'{}\n'.format(linBoundMat.shape[0]))
+for n in range(linBoundMat.shape[0]):
+    g.write('{} {} {} {}\n'.format(linBoundMat[n,0],linBoundMat[n,1],linBoundMat[n,2],linBoundMat[n,3]))
+g.write('\n'+'vertices\n'+'{}\n'.format(linVertMat.shape[0]))
+g.write('\n'+'nodes'+'\n'+'FiniteElementSpace'+'\n'+'FiniteElementCollection: H1_2D_P3'+'\n'+'VDim: 2'+'\n'+'Ordering: 1' +'\n\n')
+for n in range(linVertMatRound.shape[0]):
+    g.write('{} {}\n'.format(linVertMatRound[n,0],linVertMatRound[n,1]))
+for n in range(cubeDofMat.shape[0]):
+    g.write('{} {}\n'.format(cubeDofMat[n,0],cubeDofMat[n,1]))
+for n in range(triCentroidLoc.shape[0]):
+    g.write('{} {}\n'.format(triCentroidLoc[n,0],triCentroidLoc[n,1]))
+for n in range(quadCentroidLocCubic.shape[0]):
+    g.write('{} {}\n'.format(quadCentroidLocCubic[n,0],quadCentroidLocCubic[n,1]))
+g.close()
+
+if(visMesh==True):
+    gVis(glvis,outputName+'Cub.mesh')
+#raw_input()
+#'Reflecting' topology about one of its edges and append it to itself
+upperPlaneEleMat = sp.zeros([2*linEleMat.shape[0],6])
+for n in range(linEleMat.shape[0]):
+    upperPlaneEleMat[n,:]=linEleMat[n,:]
+    
+#Create ele_mat_holder.shape[0]x2 matrix for mapping
+refEdge=boundMatTot(numEdges)[1]
+q1NumNodes=linVertMat.shape[0]
+
+mapping = sp.zeros([q1NumNodes])
+counter=0
+for n in range(q1NumNodes):
+     if (sp.any(refEdge == n)):
+         mapping[n]=n
+     else:
+         mapping[n]=counter+q1NumNodes
+         counter+=1
+
+mapping=mapping.astype(int)
+#Implement mapping
+
+counter=0
+for n in range(linEleMat.shape[0],2*linEleMat.shape[0]):
+    upperPlaneEleMat[n,0]=linEleMat[counter,0]
+    upperPlaneEleMat[n,1]=linEleMat[counter,1]
+    upperPlaneEleMat[n,2]=mapping[linEleMat[counter,2]]
+    upperPlaneEleMat[n,3]=mapping[linEleMat[counter,3]]
+    upperPlaneEleMat[n,4]=mapping[linEleMat[counter,4]]
+    upperPlaneEleMat[n,5]=mapping[linEleMat[counter,5]]
+    counter+=1
+    
+upperPlaneEleMat = upperPlaneEleMat.astype(int)
+
+#Reflecting boundary matrix
+origBound=boundMatTot(numEdges)[2]
+upperPlaneBoundMat=sp.zeros([2*origBound.shape[0],4])
+for n in range(origBound.shape[0]):
+    upperPlaneBoundMat[n,:]=origBound[n,:]
+counter=0
+newOrigBound=origBound.copy()
+newOrigBound[:,2]=sp.flipud(origBound[:,3])
+newOrigBound[:,3]=sp.flipud(origBound[:,2])
+for n in range(newOrigBound.shape[0],upperPlaneBoundMat.shape[0]):
+    upperPlaneBoundMat[n,0]=newOrigBound[counter,0]
+    upperPlaneBoundMat[n,1]=newOrigBound[counter,1]
+    upperPlaneBoundMat[n,2]=mapping[newOrigBound[counter,2]]
+    upperPlaneBoundMat[n,3]=mapping[newOrigBound[counter,3]]
+    counter+=1
+upperPlaneBoundMat=upperPlaneBoundMat.astype(int)
+
+#Reflecting vertex matrix about the y-axis and appending it to itself
+upperPlaneNumNodes=q1NumNodes+(q1NumNodes-refEdge.shape[0])
+upperPlaneVertMat = sp.zeros([upperPlaneNumNodes,2])
+for n in range(linVertMat.shape[0]):
+    upperPlaneVertMat[n,:]=linVertMat[n,:]
+counter=0
+for n in range(linVertMat.shape[0],upperPlaneNumNodes):
+        upperPlaneVertMat[n,0]=-1.0*linVertMat[sp.where(mapping==n)[0][0],0]
+        upperPlaneVertMat[n,1]=linVertMat[sp.where(mapping==n)[0][0],1]
+        counter+=1
+
+upperPlaneEleMat=orient(upperPlaneEleMat)
+g=open(outputName+'UpperPlaneLin.mesh','w')
+g.write('MFEM mesh v1.0\n'+'\n')
+g.write('dimension\n'+'2\n'+'\n')
+g.write('elements\n'+'{}\n'.format(upperPlaneEleMat.shape[0]))
+for n in range(upperPlaneEleMat.shape[0]):
+    if upperPlaneEleMat[n,1]==2:
+        g.write('{} {} {} {} {}\n'.format(upperPlaneEleMat[n,0],upperPlaneEleMat[n,1],upperPlaneEleMat[n,2],upperPlaneEleMat[n,3],upperPlaneEleMat[n,4]))
+    else:
+        g.write('{} {} {} {} {} {}\n'.format(upperPlaneEleMat[n,0],upperPlaneEleMat[n,1],upperPlaneEleMat[n,2],upperPlaneEleMat[n,3],upperPlaneEleMat[n,4],upperPlaneEleMat[n,5]))
+g.write('\n'+'boundary\n'+'{}\n'.format(upperPlaneBoundMat.shape[0]))
+for n in range(upperPlaneBoundMat.shape[0]):
+    g.write('{} {} {} {}\n'.format(upperPlaneBoundMat[n,0],upperPlaneBoundMat[n,1],upperPlaneBoundMat[n,2],upperPlaneBoundMat[n,3]))
+g.write('\n'+'vertices\n'+'{}\n'.format(upperPlaneVertMat.shape[0])+'2\n')
+for n in range(upperPlaneVertMat.shape[0]):
+    g.write('{} {}\n'.format(upperPlaneVertMat[n,0],upperPlaneVertMat[n,1]))
+g.close()
+
+if(visMesh==True):
+    gVis(glvis,outputName+'UpperPlaneLin.mesh')
+
+#'Reflecting' topology about one of its edges and append it to itself
+wholePlaneEleMat = sp.zeros([2*upperPlaneEleMat.shape[0],6])
+for n in range(upperPlaneEleMat.shape[0]):
+    wholePlaneEleMat[n,:]=upperPlaneEleMat[n,:]
+
+quad1Edge=boundMatTot(numEdges)[3]
+newRefEdge=sp.zeros(2*quad1Edge.shape[0]-1)
+for n in range(quad1Edge.shape[0]):
+    newRefEdge[n]=quad1Edge[n]
+counter=0
+for n in range(quad1Edge.shape[0],newRefEdge.shape[0]):
+    newRefEdge[n]=mapping[quad1Edge[counter]]
+    counter+=1
+newRefEdge=sp.unique(newRefEdge)
+newRefEdge=newRefEdge.astype(int)
+newTotNumNodes=upperPlaneVertMat.shape[0]
+
+newMapping=sp.zeros([newTotNumNodes])
+counter=0
+for n in range(newTotNumNodes):
+    if (sp.any(newRefEdge == n)):
+        newMapping[n]=n
+    else:
+        newMapping[n]=counter+newTotNumNodes
+        counter+=1
+newMapping=newMapping.astype(int)
+
+counter=0
+for n in range(upperPlaneEleMat.shape[0],2*upperPlaneEleMat.shape[0]):
+    wholePlaneEleMat[n,0]=upperPlaneEleMat[counter,0]
+    wholePlaneEleMat[n,1]=upperPlaneEleMat[counter,1]
+    wholePlaneEleMat[n,2]=newMapping[upperPlaneEleMat[counter,2]]
+    wholePlaneEleMat[n,3]=newMapping[upperPlaneEleMat[counter,3]]
+    wholePlaneEleMat[n,4]=newMapping[upperPlaneEleMat[counter,4]]
+    wholePlaneEleMat[n,5]=newMapping[upperPlaneEleMat[counter,5]]
+    counter+=1
+
+wholePlaneEleMat=wholePlaneEleMat.astype(int)
+
+#Reflecting boundary matrix
+newOrigBoundQuad1=boundMatTot(numEdges)[4]
+newFirstBoundMatHolder=sp.zeros([2*newOrigBoundQuad1.shape[0],4])
+for n in range(newOrigBoundQuad1.shape[0]):
+    newFirstBoundMatHolder[n,:]=newOrigBoundQuad1[n,:]
+
+newNewOrigBoundQuad1=newOrigBoundQuad1.copy()
+newNewOrigBoundQuad1[:,2]=sp.flipud(newOrigBoundQuad1[:,3])
+newNewOrigBoundQuad1[:,3]=sp.flipud(newOrigBoundQuad1[:,2])
+counter=0
+for n in range(newOrigBoundQuad1.shape[0],newFirstBoundMatHolder.shape[0]):
+    newFirstBoundMatHolder[n,0]=newNewOrigBoundQuad1[counter,0]
+    newFirstBoundMatHolder[n,1]=newNewOrigBoundQuad1[counter,1]
+    newFirstBoundMatHolder[n,2]=mapping[newNewOrigBoundQuad1[counter,2]]
+    newFirstBoundMatHolder[n,3]=mapping[newNewOrigBoundQuad1[counter,3]]
+    counter+=1
+
+upperQuadMat=newFirstBoundMatHolder.copy()
+wholePlaneBoundMat=sp.zeros([2*upperQuadMat.shape[0],4])
+for n in range(upperQuadMat.shape[0]):
+    wholePlaneBoundMat[n,:]=upperQuadMat[n,:]
+
+counter=0
+newNewOrigBound=upperQuadMat.copy()
+newNewOrigBound[:,2]=sp.flipud(upperQuadMat[:,3])
+newNewOrigBound[:,3]=sp.flipud(upperQuadMat[:,2])
+newNewOrigBound=newNewOrigBound.astype(int)
+for n in range(newNewOrigBound.shape[0],wholePlaneBoundMat.shape[0]):
+    wholePlaneBoundMat[n,0]=newNewOrigBound[counter,0]
+    wholePlaneBoundMat[n,1]=newNewOrigBound[counter,1]
+    wholePlaneBoundMat[n,2]=newMapping[newNewOrigBound[counter,2]]
+    wholePlaneBoundMat[n,3]=newMapping[newNewOrigBound[counter,3]]
+    counter+=1
+wholePlaneBoundMat=wholePlaneBoundMat.astype(int)
+
+wholePlaneNumNodes=newTotNumNodes+(newTotNumNodes-newRefEdge.shape[0])
+wholePlaneVertMat = sp.zeros([wholePlaneNumNodes,2])
+for n in range(upperPlaneVertMat.shape[0]):
+    wholePlaneVertMat[n,:]=upperPlaneVertMat[n,:]
+counter=0
+for n in range(upperPlaneVertMat.shape[0],wholePlaneNumNodes):
+        wholePlaneVertMat[n,0]=upperPlaneVertMat[sp.where(newMapping==n)[0][0],0]
+        wholePlaneVertMat[n,1]=-1.0*upperPlaneVertMat[sp.where(newMapping==n)[0][0],1]
+        counter+=1
+
+g=open(outputName+'WholePlaneLin.mesh','w')
+g.write('MFEM mesh v1.0\n'+'\n')
+g.write('dimension\n'+'2\n'+'\n')
+g.write('elements\n'+'{}\n'.format(wholePlaneEleMat.shape[0]))
+for n in range(wholePlaneEleMat.shape[0]):
+    if wholePlaneEleMat[n,1]==2:
+        g.write('{} {} {} {} {}\n'.format(wholePlaneEleMat[n,0],wholePlaneEleMat[n,1],wholePlaneEleMat[n,2],wholePlaneEleMat[n,3],wholePlaneEleMat[n,4]))
+    else:
+        g.write('{} {} {} {} {} {}\n'.format(wholePlaneEleMat[n,0],wholePlaneEleMat[n,1],wholePlaneEleMat[n,2],wholePlaneEleMat[n,3],wholePlaneEleMat[n,4],wholePlaneEleMat[n,5]))
+g.write('\n'+'boundary\n'+'{}\n'.format(wholePlaneBoundMat.shape[0]))
+for n in range(wholePlaneBoundMat.shape[0]):
+    g.write('{} {} {} {}\n'.format(wholePlaneBoundMat[n,0],wholePlaneBoundMat[n,1],wholePlaneBoundMat[n,2],wholePlaneBoundMat[n,3]))
+g.write('\n'+'vertices\n'+'{}\n'.format(wholePlaneVertMat.shape[0])+'2\n')
+for n in range(wholePlaneVertMat.shape[0]):
+    g.write('{} {}\n'.format(wholePlaneVertMat[n,0],wholePlaneVertMat[n,1]))
+g.close()
+
+if(visMesh==True):
+    gVis(glvis,outputName+'WholePlaneLin.mesh')
+
+#1.)Create Edge list from elements
+wholePlaneEleMat=orient(wholePlaneEleMat)
+triCounter=0;quadCounter=0;
+for n in range(wholePlaneEleMat.shape[0]):
+    if wholePlaneEleMat[n,1]==2:
+        triCounter+=1
+    else:
+        quadCounter+=1
+edgeMat=sp.zeros([3*triCounter+4*quadCounter,2])
+counter=0
+for n in range(wholePlaneEleMat.shape[0]):
+    if wholePlaneEleMat[n,1]==2:
+        edgeMat[counter,:]=[wholePlaneEleMat[n,2],wholePlaneEleMat[n,3]]
+        counter+=1
+        edgeMat[counter,:]=[wholePlaneEleMat[n,3],wholePlaneEleMat[n,4]]
+        counter+=1
+        edgeMat[counter,:]=[wholePlaneEleMat[n,4],wholePlaneEleMat[n,2]]
+        counter+=1
+    else:
+        edgeMat[counter,:]=[wholePlaneEleMat[n,2],wholePlaneEleMat[n,3]]
+        counter+=1
+        edgeMat[counter,:]=[wholePlaneEleMat[n,3],wholePlaneEleMat[n,4]]
+        counter+=1
+        edgeMat[counter,:]=[wholePlaneEleMat[n,4],wholePlaneEleMat[n,5]]
+        counter+=1
+        edgeMat[counter,:]=[wholePlaneEleMat[n,5],wholePlaneEleMat[n,2]]
+        counter+=1
+
+#Remove duplicates
+holder=[]
+for n in range(edgeMat.shape[0]):
+    counter=0
+    for m in range(edgeMat.shape[0]):
+        if edgeMat[n,0]==edgeMat[m,0] and edgeMat[n,1]==edgeMat[m,1] and m!=n:
+            holder.append([n,m])
+        elif edgeMat[n,1]==edgeMat[m,0] and edgeMat[n,0]==edgeMat[m,1] and m!=n:
+            holder.append([n,m])
+removeIndices=sp.zeros(len(holder))
+for n in range(len(holder)):
+    if holder[n][0]>holder[n][1]:
+        removeIndices[n]=holder[n][0]
+    else:
+        removeIndices[n]=holder[n][1]
+removeIndices=sp.unique(removeIndices).astype(int)
+edgeMat=sp.delete(edgeMat,removeIndices,0)
+edgeMat=edgeMat.astype(int)
+
+edgeDofMat=sp.zeros([edgeMat.shape[0],2])
+
+wholePlaneVertMatRound=sp.round_(wholePlaneVertMat,5)
+
+counter=0
+for n in edgeMat:
+    if wholePlaneVertMatRound[n[0],1] == wholePlaneVertMatRound[n[1],1]:
+        xmid=(wholePlaneVertMatRound[n[0],0]+wholePlaneVertMatRound[n[1],0])/2.0
+        ymid=wholePlaneVertMatRound[n[0],1]
+        edgeDofMat[counter,:]=[xmid,ymid]
+    elif wholePlaneVertMatRound[n[0],0] == wholePlaneVertMatRound[n[1],0]:
+        xmid=wholePlaneVertMatRound[n[0],0]
+        ymid=(wholePlaneVertMatRound[n[0],1]+wholePlaneVertMatRound[n[1],1])/2.0
+        edgeDofMat[counter,:]=[xmid,ymid]
+    else:
+        r0=sp.sqrt(wholePlaneVertMatRound[n[0],0]**2+wholePlaneVertMatRound[n[0],1]**2)
+        r1=sp.sqrt(wholePlaneVertMatRound[n[1],0]**2+wholePlaneVertMatRound[n[1],1]**2)
+        rmid = (r0+r1)/2.0 #should not be needed
+        xmidOld=(wholePlaneVertMatRound[n[0],0]+wholePlaneVertMatRound[n[1],0])/2.0
+        ymidOld=(wholePlaneVertMatRound[n[0],1]+wholePlaneVertMatRound[n[1],1])/2.0
+        midtheta=sp.arctan2(ymidOld,xmidOld)
+        xmid=rmid*sp.cos(midtheta)
+        ymid=rmid*sp.sin(midtheta)
+        edgeDofMat[counter,:]=[xmid,ymid]
+    counter+=1
+edgeDofMat = sp.round_(edgeDofMat,5)
+
+#2.)Create correct dof locations
+#Determine midpoints of all quads:
+quadCentroidLoc=sp.zeros([quadCounter,2])
+counter=0
+for n in range(wholePlaneEleMat.shape[0]):
+    if wholePlaneEleMat[n,1]==3:
+        quadCentroidLoc[counter,0]=(wholePlaneVertMatRound[wholePlaneEleMat[n,2],0]+wholePlaneVertMatRound[wholePlaneEleMat[n,3],0]+wholePlaneVertMatRound[wholePlaneEleMat[n,4],0]+wholePlaneVertMatRound[wholePlaneEleMat[n,5],0])/4.0
+        quadCentroidLoc[counter,1]=(wholePlaneVertMatRound[wholePlaneEleMat[n,2],1]+wholePlaneVertMatRound[wholePlaneEleMat[n,3],1]+wholePlaneVertMatRound[wholePlaneEleMat[n,4],1]+wholePlaneVertMatRound[wholePlaneEleMat[n,5],1])/4.0
+        counter+=1
+
+quadCentroidLoc = sp.round_(quadCentroidLoc,5)
+
+#3.)Populate nodes section
+g=open(outputName+'WholePlaneQuad.mesh','w')
+g.write('MFEM mesh v1.0\n'+'\n')
+g.write('dimension\n'+'2\n'+'\n')
+g.write('elements\n'+'{}\n'.format(wholePlaneEleMat.shape[0]))
+for n in range(wholePlaneEleMat.shape[0]):
+    if wholePlaneEleMat[n,1]==2:
+        g.write('{} {} {} {} {}\n'.format(wholePlaneEleMat[n,0],wholePlaneEleMat[n,1],wholePlaneEleMat[n,2],wholePlaneEleMat[n,3],wholePlaneEleMat[n,4]))
+    else:
+        g.write('{} {} {} {} {} {}\n'.format(wholePlaneEleMat[n,0],wholePlaneEleMat[n,1],wholePlaneEleMat[n,2],wholePlaneEleMat[n,3],wholePlaneEleMat[n,4],wholePlaneEleMat[n,5]))
+g.write('\n'+'boundary\n'+'{}\n'.format(wholePlaneBoundMat.shape[0]))
+for n in range(wholePlaneBoundMat.shape[0]):
+    g.write('{} {} {} {}\n'.format(wholePlaneBoundMat[n,0],wholePlaneBoundMat[n,1],wholePlaneBoundMat[n,2],wholePlaneBoundMat[n,3]))
+g.write('\n'+'vertices\n'+'{}\n'.format(wholePlaneVertMat.shape[0]))
+g.write('\n'+'nodes'+'\n'+'FiniteElementSpace'+'\n'+'FiniteElementCollection: H1_2D_P2'+'\n'+'VDim: 2'+'\n'+'Ordering: 1' +'\n\n')
+for n in range(wholePlaneVertMatRound.shape[0]):
+    g.write('{} {}\n'.format(wholePlaneVertMatRound[n,0],wholePlaneVertMatRound[n,1]))
+for n in range(edgeDofMat.shape[0]):
+    g.write('{} {}\n'.format(edgeDofMat[n,0],edgeDofMat[n,1]))
+for n in range(quadCentroidLoc.shape[0]):
+    g.write('{} {}\n'.format(quadCentroidLoc[n,0],quadCentroidLoc[n,1]))
+g.close()
+
+if(visMesh==True):
+    gVis(glvis,outputName+'WholePlaneQuad.mesh')

data/mixed-ele-meshes/comp_tri_quad_4_edges.mesh

@@ -0,0 +1,264 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+128
+1 2 0 1 2
+1 2 1 2 4
+1 2 1 3 4
+1 2 2 4 5
+1 2 3 4 7
+1 2 3 6 7
+1 2 4 5 8
+1 2 4 7 8
+1 2 5 8 9
+1 2 6 7 11
+1 2 6 10 11
+1 2 7 8 12
+1 2 7 11 12
+1 2 8 9 13
+1 2 8 12 13
+1 2 9 13 14
+2 3 10 15 16 11
+2 3 11 16 17 12
+2 3 12 17 18 13
+2 3 13 18 19 14
+2 3 15 20 21 16
+2 3 16 21 22 17
+2 3 17 22 23 18
+2 3 18 23 24 19
+2 3 20 25 26 21
+2 3 21 26 27 22
+2 3 22 27 28 23
+2 3 23 28 29 24
+2 3 25 30 31 26
+2 3 26 31 32 27
+2 3 27 32 33 28
+2 3 28 33 34 29
+1 2 0 35 2
+1 2 35 2 37
+1 2 35 36 37
+1 2 2 37 5
+1 2 36 37 39
+1 2 36 38 39
+1 2 37 5 40
+1 2 37 39 40
+1 2 5 40 9
+1 2 38 39 42
+1 2 38 41 42
+1 2 39 40 43
+1 2 39 42 43
+1 2 40 9 44
+1 2 40 43 44
+1 2 9 44 14
+2 3 41 45 46 42
+2 3 42 46 47 43
+2 3 43 47 48 44
+2 3 44 48 19 14
+2 3 45 49 50 46
+2 3 46 50 51 47
+2 3 47 51 52 48
+2 3 48 52 24 19
+2 3 49 53 54 50
+2 3 50 54 55 51
+2 3 51 55 56 52
+2 3 52 56 29 24
+2 3 53 57 58 54
+2 3 54 58 59 55
+2 3 55 59 60 56
+2 3 56 60 34 29
+1 2 0 1 61
+1 2 1 61 62
+1 2 1 3 62
+1 2 61 62 63
+1 2 3 62 64
+1 2 3 6 64
+1 2 62 63 65
+1 2 62 64 65
+1 2 63 65 66
+1 2 6 64 67
+1 2 6 10 67
+1 2 64 65 68
+1 2 64 67 68
+1 2 65 66 69
+1 2 65 68 69
+1 2 66 69 70
+2 3 10 15 71 67
+2 3 67 71 72 68
+2 3 68 72 73 69
+2 3 69 73 74 70
+2 3 15 20 75 71
+2 3 71 75 76 72
+2 3 72 76 77 73
+2 3 73 77 78 74
+2 3 20 25 79 75
+2 3 75 79 80 76
+2 3 76 80 81 77
+2 3 77 81 82 78
+2 3 25 30 83 79
+2 3 79 83 84 80
+2 3 80 84 85 81
+2 3 81 85 86 82
+1 2 0 35 61
+1 2 35 61 87
+1 2 35 36 87
+1 2 61 87 63
+1 2 36 87 88
+1 2 36 38 88
+1 2 87 63 89
+1 2 87 88 89
+1 2 63 89 66
+1 2 38 88 90
+1 2 38 41 90
+1 2 88 89 91
+1 2 88 90 91
+1 2 89 66 92
+1 2 89 91 92
+1 2 66 92 70
+2 3 41 45 93 90
+2 3 90 93 94 91
+2 3 91 94 95 92
+2 3 92 95 74 70
+2 3 45 49 96 93
+2 3 93 96 97 94
+2 3 94 97 98 95
+2 3 95 98 78 74
+2 3 49 53 99 96
+2 3 96 99 100 97
+2 3 97 100 101 98
+2 3 98 101 82 78
+2 3 53 102 103 99
+2 3 99 103 104 100
+2 3 100 104 105 101
+2 3 101 105 86 82
+
+boundary
+16
+1 1 30 31
+1 1 31 32
+1 1 32 33
+1 1 33 34
+1 1 34 60
+1 1 60 59
+1 1 59 58
+1 1 58 57
+1 1 102 103
+1 1 103 104
+1 1 104 105
+1 1 105 86
+1 1 86 85
+1 1 85 84
+1 1 84 83
+1 1 83 30
+
+vertices
+106
+2
+0.0 0.0
+0.125 0.0
+7.65404249467e-18 0.125
+0.25 0.0
+0.176776695297 0.176776695297
+1.53080849893e-17 0.25
+0.375 0.0
+0.324759526419 0.1875
+0.1875 0.324759526419
+2.2962127484e-17 0.375
+0.5 0.0
+0.461939766256 0.191341716183
+0.353553390593 0.353553390593
+0.191341716183 0.461939766256
+3.06161699787e-17 0.5
+0.625 0.0
+0.596454824692 0.268506287137
+0.515165042945 0.515165042945
+0.268506287137 0.596454824692
+2.2962127484e-17 0.625
+0.75 0.0
+0.730969883128 0.345670858091
+0.676776695297 0.676776695297
+0.345670858091 0.730969883128
+1.53080849893e-17 0.75
+0.875 0.0
+0.865484941564 0.422835429046
+0.838388347648 0.838388347648
+0.422835429046 0.865484941564
+7.65404249467e-18 0.875
+1.0 0.0
+1.0 0.5
+1.0 1.0
+0.5 1.0
+0.0 1.0
+-0.125 0.0
+-0.25 0.0
+-0.176776695297 0.176776695297
+-0.375 0.0
+-0.324759526419 0.1875
+-0.1875 0.324759526419
+-0.5 0.0
+-0.461939766256 0.191341716183
+-0.353553390593 0.353553390593
+-0.191341716183 0.461939766256
+-0.625 0.0
+-0.596454824692 0.268506287137
+-0.515165042945 0.515165042945
+-0.268506287137 0.596454824692
+-0.75 0.0
+-0.730969883128 0.345670858091
+-0.676776695297 0.676776695297
+-0.345670858091 0.730969883128
+-0.875 0.0
+-0.865484941564 0.422835429046
+-0.838388347648 0.838388347648
+-0.422835429046 0.865484941564
+-1.0 0.0
+-1.0 0.5
+-1.0 1.0
+-0.5 1.0
+7.65404249467e-18 -0.125
+0.176776695297 -0.176776695297
+1.53080849893e-17 -0.25
+0.324759526419 -0.1875
+0.1875 -0.324759526419
+2.2962127484e-17 -0.375
+0.461939766256 -0.191341716183
+0.353553390593 -0.353553390593
+0.191341716183 -0.461939766256
+3.06161699787e-17 -0.5
+0.596454824692 -0.268506287137
+0.515165042945 -0.515165042945
+0.268506287137 -0.596454824692
+2.2962127484e-17 -0.625
+0.730969883128 -0.345670858091
+0.676776695297 -0.676776695297
+0.345670858091 -0.730969883128
+1.53080849893e-17 -0.75
+0.865484941564 -0.422835429046
+0.838388347648 -0.838388347648
+0.422835429046 -0.865484941564
+7.65404249467e-18 -0.875
+1.0 -0.5
+1.0 -1.0
+0.5 -1.0
+0.0 -1.0
+-0.176776695297 -0.176776695297
+-0.324759526419 -0.1875
+-0.1875 -0.324759526419
+-0.461939766256 -0.191341716183
+-0.353553390593 -0.353553390593
+-0.191341716183 -0.461939766256
+-0.596454824692 -0.268506287137
+-0.515165042945 -0.515165042945
+-0.268506287137 -0.596454824692
+-0.730969883128 -0.345670858091
+-0.676776695297 -0.676776695297
+-0.345670858091 -0.730969883128
+-0.865484941564 -0.422835429046
+-0.838388347648 -0.838388347648
+-0.422835429046 -0.865484941564
+-1.0 -0.0
+-1.0 -0.5
+-1.0 -1.0
+-0.5 -1.0

data/mixed-ele-meshes/spec-3-mixed-type.mesh

@@ -0,0 +1,72 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+17
+1 2 0 1 2
+1 2 0 2 3
+1 2 0 3 4
+2 3 0 4 5 6
+2 3 0 6 7 1
+1 2 7 8 1
+1 2 1 8 9
+2 3 1 9 10 2
+1 2 2 10 11
+2 3 2 11 12 3
+1 2 3 12 13
+2 3 3 13 14 4
+1 2 4 14 15
+1 2 4 15 5
+1 2 5 16 6 
+1 2 6 16 17
+1 2 6 17 7
+
+boundary
+12
+1 1 7 8
+1 1 8 9
+1 1 9 10
+1 1 10 11
+1 1 11 12
+1 1 12 13
+1 1 13 14
+1 1 14 15
+1 1 15 5
+1 1 5 16
+1 1 16 17
+1 1 17 7
+
+vertices
+18
+2
+0 0
+1 0
+0.5 0.866025
+-0.5 0.866025
+-1 0
+-1 -1
+0 -1
+1 -1
+1.866025 -0.5
+1.866025 0.5
+1.366025 1.366025
+0.5 1.866025
+-0.5 1.866025
+-1.366025 1.366025
+-1.866025 0.5
+-1.866025 -0.5
+-0.5 -1.866025
+0.5 -1.866025

data/mixed-ele-meshes/spec-3-trans-and-ref.mesh

@@ -0,0 +1,362 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+154
+2 3 0 1 2 3
+2 3 1 5 6 2
+2 3 5 8 9 6
+2 3 8 11 12 9
+2 3 11 14 15 12
+2 3 14 17 18 15
+2 3 17 20 21 18
+2 3 20 23 24 21
+2 3 23 26 27 24
+2 3 26 29 30 27
+2 3 29 32 33 30
+2 3 32 35 36 33
+2 3 35 38 39 36
+2 3 38 41 42 39
+2 3 41 44 45 42
+2 3 44 47 48 45
+2 3 47 50 51 48
+2 3 50 53 54 51
+2 3 53 56 57 54
+2 3 56 59 60 57
+2 3 59 62 63 60
+2 3 62 65 66 63
+2 3 65 68 69 66
+2 3 68 71 72 69
+2 3 71 74 75 72
+2 3 74 77 78 75
+1 2 2 3 4
+1 2 6 2 7
+1 2 9 6 10
+1 2 12 9 13
+1 2 15 12 16
+1 2 18 15 19
+1 2 21 18 22
+1 2 24 21 25
+1 2 27 24 28
+1 2 30 27 31
+1 2 33 30 34
+1 2 36 33 37
+1 2 39 36 40
+1 2 42 39 43
+1 2 45 42 46
+1 2 48 45 49
+1 2 51 48 52
+1 2 54 51 55
+1 2 57 54 58
+1 2 60 57 61
+1 2 63 60 64
+1 2 66 63 67
+1 2 69 66 70
+1 2 72 69 73
+1 2 75 72 76
+1 2 78 75 79
+1 2 2 4 7
+1 2 6 7 10
+1 2 9 10 13
+1 2 12 13 16
+1 2 15 16 19
+1 2 18 19 22
+1 2 21 22 25
+1 2 24 25 28
+1 2 27 28 31
+1 2 30 31 34
+1 2 33 34 37
+1 2 36 37 40
+1 2 39 40 43
+1 2 42 43 46
+1 2 45 46 49
+1 2 48 49 52
+1 2 51 52 55
+1 2 54 55 58
+1 2 57 58 61
+1 2 60 61 64
+1 2 63 64 67
+1 2 66 67 70
+1 2 69 70 73
+1 2 72 73 76
+1 2 75 76 79
+2 3 80 81 82 83
+2 3 81 84 85 82
+2 3 84 86 87 85
+2 3 86 88 89 87
+2 3 88 90 91 89
+2 3 90 92 93 91
+2 3 92 94 95 93
+2 3 94 96 97 95
+2 3 96 98 99 97
+2 3 98 100 101 99
+2 3 100 102 103 101
+2 3 102 104 105 103
+2 3 104 106 107 105
+2 3 106 108 109 107
+2 3 108 110 111 109
+2 3 110 112 113 111
+2 3 112 114 115 113
+2 3 114 116 117 115
+2 3 116 118 119 117
+2 3 118 120 121 119
+2 3 120 122 123 121
+2 3 122 124 125 123
+2 3 124 126 127 125
+2 3 126 128 129 127
+2 3 128 130 131 129
+2 3 130 132 133 131
+1 2 82 83 4
+1 2 85 82 7
+1 2 87 85 10
+1 2 89 87 13
+1 2 91 89 16
+1 2 93 91 19
+1 2 95 93 22
+1 2 97 95 25
+1 2 99 97 28
+1 2 101 99 31
+1 2 103 101 34
+1 2 105 103 37
+1 2 107 105 40
+1 2 109 107 43
+1 2 111 109 46
+1 2 113 111 49
+1 2 115 113 52
+1 2 117 115 55
+1 2 119 117 58
+1 2 121 119 61
+1 2 123 121 64
+1 2 125 123 67
+1 2 127 125 70
+1 2 129 127 73
+1 2 131 129 76
+1 2 133 131 79
+1 2 82 4 7
+1 2 85 7 10
+1 2 87 10 13
+1 2 89 13 16
+1 2 91 16 19
+1 2 93 19 22
+1 2 95 22 25
+1 2 97 25 28
+1 2 99 28 31
+1 2 101 31 34
+1 2 103 34 37
+1 2 105 37 40
+1 2 107 40 43
+1 2 109 43 46
+1 2 111 46 49
+1 2 113 49 52
+1 2 115 52 55
+1 2 117 55 58
+1 2 119 58 61
+1 2 121 61 64
+1 2 123 64 67
+1 2 125 67 70
+1 2 127 70 73
+1 2 129 73 76
+1 2 131 76 79
+
+boundary
+60
+1 1 0 1
+1 1 1 5
+1 1 5 8
+1 1 8 11
+1 1 11 14
+1 1 14 17
+1 1 17 20
+1 1 20 23
+1 1 23 26
+1 1 26 29
+1 1 29 32
+1 1 32 35
+1 1 35 38
+1 1 38 41
+1 1 41 44
+1 1 44 47
+1 1 47 50
+1 1 50 53
+1 1 53 56
+1 1 56 59
+1 1 59 62
+1 1 62 65
+1 1 65 68
+1 1 68 71
+1 1 71 74
+1 1 74 77
+1 1 77 78
+1 1 78 79
+1 1 79 133
+1 1 133 132
+1 1 132 130
+1 1 130 128
+1 1 128 126
+1 1 126 124
+1 1 124 122
+1 1 122 120
+1 1 120 118
+1 1 118 116
+1 1 116 114
+1 1 114 112
+1 1 112 110
+1 1 110 108
+1 1 108 106
+1 1 106 104
+1 1 104 102
+1 1 102 100
+1 1 100 98
+1 1 98 96
+1 1 96 94
+1 1 94 92
+1 1 92 90
+1 1 90 88
+1 1 88 86
+1 1 86 84
+1 1 84 81
+1 1 81 80
+1 1 80 83
+1 1 83 4
+1 1 4 3
+1 1 3 0
+
+vertices
+134
+2
+0.0 0.0
+1.0 0.0
+1.0 1.0
+0.0 1.0
+0.5 1.86602540378
+2.0 0.0
+2.0 1.0
+1.5 1.86602540378
+3.0 0.0
+3.0 1.0
+2.5 1.86602540378
+4.0 0.0
+4.0 1.0
+3.5 1.86602540378
+5.0 0.0
+5.0 1.0
+4.5 1.86602540378
+6.0 0.0
+6.0 1.0
+5.5 1.86602540378
+7.0 0.0
+7.0 1.0
+6.5 1.86602540378
+8.0 0.0
+8.0 1.0
+7.5 1.86602540378
+9.0 0.0
+9.0 1.0
+8.5 1.86602540378
+10.0 0.0
+10.0 1.0
+9.5 1.86602540378
+11.0 0.0
+11.0 1.0
+10.5 1.86602540378
+12.0 0.0
+12.0 1.0
+11.5 1.86602540378
+13.0 0.0
+13.0 1.0
+12.5 1.86602540378
+14.0 0.0
+14.0 1.0
+13.5 1.86602540378
+15.0 0.0
+15.0 1.0
+14.5 1.86602540378
+16.0 0.0
+16.0 1.0
+15.5 1.86602540378
+17.0 0.0
+17.0 1.0
+16.5 1.86602540378
+18.0 0.0
+18.0 1.0
+17.5 1.86602540378
+19.0 0.0
+19.0 1.0
+18.5 1.86602540378
+20.0 0.0
+20.0 1.0
+19.5 1.86602540378
+21.0 0.0
+21.0 1.0
+20.5 1.86602540378
+22.0 0.0
+22.0 1.0
+21.5 1.86602540378
+23.0 0.0
+23.0 1.0
+22.5 1.86602540378
+24.0 0.0
+24.0 1.0
+23.5 1.86602540378
+25.0 0.0
+25.0 1.0
+24.5 1.86602540378
+26.0 0.0
+26.0 1.0
+25.5 1.86602540378
+0.0 3.73205080757
+1.0 3.73205080757
+1.0 2.73205080757
+0.0 2.73205080757
+2.0 3.73205080757
+2.0 2.73205080757
+3.0 3.73205080757
+3.0 2.73205080757
+4.0 3.73205080757
+4.0 2.73205080757
+5.0 3.73205080757
+5.0 2.73205080757
+6.0 3.73205080757
+6.0 2.73205080757
+7.0 3.73205080757
+7.0 2.73205080757
+8.0 3.73205080757
+8.0 2.73205080757
+9.0 3.73205080757
+9.0 2.73205080757
+10.0 3.73205080757
+10.0 2.73205080757
+11.0 3.73205080757
+11.0 2.73205080757
+12.0 3.73205080757
+12.0 2.73205080757
+13.0 3.73205080757
+13.0 2.73205080757
+14.0 3.73205080757
+14.0 2.73205080757
+15.0 3.73205080757
+15.0 2.73205080757
+16.0 3.73205080757
+16.0 2.73205080757
+17.0 3.73205080757
+17.0 2.73205080757
+18.0 3.73205080757
+18.0 2.73205080757
+19.0 3.73205080757
+19.0 2.73205080757
+20.0 3.73205080757
+20.0 2.73205080757
+21.0 3.73205080757
+21.0 2.73205080757
+22.0 3.73205080757
+22.0 2.73205080757
+23.0 3.73205080757
+23.0 2.73205080757
+24.0 3.73205080757
+24.0 2.73205080757
+25.0 3.73205080757
+25.0 2.73205080757
+26.0 3.73205080757
+26.0 2.73205080757

data/mixed-ele-meshes/spec-3-type-1-trans-5.mesh

@@ -0,0 +1,66 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+17
+2 3 0 1 2 3
+2 3 1 5 6 2
+2 3 5 8 9 6
+2 3 8 11 12 9
+2 3 11 14 15 12
+2 3 14 17 18 15
+1 2 2 3 4
+1 2 6 2 7
+1 2 9 6 10
+1 2 12 9 13
+1 2 15 12 16
+1 2 18 15 19
+1 2 2 4 7
+1 2 6 7 10
+1 2 9 10 13
+1 2 12 13 16
+1 2 15 16 19
+
+boundary
+15
+1 1 0 1
+1 1 1 5
+1 1 5 8
+1 1 8 11
+1 1 11 14
+1 1 14 17
+1 1 17 18
+1 1 18 19
+1 1 19 16
+1 1 16 13
+1 1 13 10
+1 1 10 7
+1 1 7 4
+1 1 4 3
+1 1 3 0
+
+vertices
+20
+2
+0.0 0.0
+1.0 0.0
+1.0 1.0
+0.0 1.0
+0.5 1.86602540378
+2.0 0.0
+2.0 1.0
+1.5 1.86602540378
+3.0 0.0
+3.0 1.0
+2.5 1.86602540378
+4.0 0.0
+4.0 1.0
+3.5 1.86602540378
+5.0 0.0
+5.0 1.0
+4.5 1.86602540378
+6.0 0.0
+6.0 1.0
+5.5 1.86602540378

data/mixed-ele-meshes/spec-3-type-1.mesh

@@ -0,0 +1,45 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+5
+1 2 0 1 2
+1 2 0 2 3
+1 2 0 3 4
+2 3 0 4 5 6
+2 3 0 6 7 1
+
+boundary
+7
+1 1 1 2
+1 1 2 3
+1 1 3 4
+1 1 4 5
+1 1 5 6
+1 1 6 7
+1 1 7 1
+
+vertices
+8
+2
+0 0
+1 0
+0.5 0.866025
+-0.5 0.866025
+-1 0
+-1 -1
+0 -1
+1 -1

data/mixed-ele-meshes/spec-3-type-2.mesh

@@ -0,0 +1,45 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+5
+1 2 0 1 2
+1 2 0 2 3
+1 2 0 5 6
+2 3 0 3 4 5
+2 3 0 6 7 1
+
+boundary
+7
+1 1 1 2
+1 1 2 3
+1 1 3 4
+1 1 4 5
+1 1 5 6
+1 1 6 7
+1 1 7 1
+
+vertices
+8
+2
+0 0
+0.866025 0.5
+0 1
+-0.866025 0.5
+-1.366025 -0.36602500
+-0.5 -0.866025
+0.5 -0.866025
+1.366025 -0.36602500

data/mixed-ele-meshes/tri_circ_sect_4.mesh

@@ -0,0 +1,57 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+16
+1 2 0 1 2
+1 2 1 2 4
+1 2 1 3 4
+1 2 2 4 5
+1 2 3 4 7
+1 2 3 6 7
+1 2 4 5 8
+1 2 4 7 8
+1 2 5 8 9
+1 2 6 7 11
+1 2 6 10 11
+1 2 7 8 12
+1 2 7 11 12
+1 2 8 9 13
+1 2 8 12 13
+1 2 9 13 14
+
+boundary
+12
+1 1 0 1
+1 1 1 3
+1 1 3 6
+1 1 6 10
+1 1 10 11
+1 1 11 12
+1 1 12 13
+1 1 13 14
+1 1 14 9
+1 1 9 5
+1 1 5 2
+1 1 2 0
+
+vertices
+15
+2
+0.0 0.0
+0.25 0.0
+1.53080849893e-17 0.25
+0.5 0.0
+0.353553390593 0.353553390593
+3.06161699787e-17 0.5
+0.75 0.0
+0.649519052838 0.375
+0.375 0.649519052838
+4.5924254968e-17 0.75
+1.0 0.0
+0.923879532511 0.382683432365
+0.707106781187 0.707106781187
+0.382683432365 0.923879532511
+6.12323399574e-17 1.0

data/mixed-ele-meshes/tri_circ_sect_8.mesh

@@ -0,0 +1,147 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+64
+1 2 0 1 2
+1 2 1 2 4
+1 2 1 3 4
+1 2 2 4 5
+1 2 3 4 7
+1 2 3 6 7
+1 2 4 5 8
+1 2 4 7 8
+1 2 5 8 9
+1 2 6 7 11
+1 2 6 10 11
+1 2 7 8 12
+1 2 7 11 12
+1 2 8 9 13
+1 2 8 12 13
+1 2 9 13 14
+1 2 10 11 16
+1 2 10 15 16
+1 2 11 12 17
+1 2 11 16 17
+1 2 12 13 18
+1 2 12 17 18
+1 2 13 14 19
+1 2 13 18 19
+1 2 14 19 20
+1 2 15 16 22
+1 2 15 21 22
+1 2 16 17 23
+1 2 16 22 23
+1 2 17 18 24
+1 2 17 23 24
+1 2 18 19 25
+1 2 18 24 25
+1 2 19 20 26
+1 2 19 25 26
+1 2 20 26 27
+1 2 21 22 29
+1 2 21 28 29
+1 2 22 23 30
+1 2 22 29 30
+1 2 23 24 31
+1 2 23 30 31
+1 2 24 25 32
+1 2 24 31 32
+1 2 25 26 33
+1 2 25 32 33
+1 2 26 27 34
+1 2 26 33 34
+1 2 27 34 35
+1 2 28 29 37
+1 2 28 36 37
+1 2 29 30 38
+1 2 29 37 38
+1 2 30 31 39
+1 2 30 38 39
+1 2 31 32 40
+1 2 31 39 40
+1 2 32 33 41
+1 2 32 40 41
+1 2 33 34 42
+1 2 33 41 42
+1 2 34 35 43
+1 2 34 42 43
+1 2 35 43 44
+
+boundary
+24
+1 1 0 1
+1 1 1 3
+1 1 3 6
+1 1 6 10
+1 1 10 15
+1 1 15 21
+1 1 21 28
+1 1 28 36
+1 1 36 37
+1 1 37 38
+1 1 38 39
+1 1 39 40
+1 1 40 41
+1 1 41 42
+1 1 42 43
+1 1 43 44
+1 1 44 35
+1 1 35 27
+1 1 27 20
+1 1 20 14
+1 1 14 9
+1 1 9 5
+1 1 5 2
+1 1 2 0
+
+vertices
+45
+2
+0.0 0.0
+0.125 0.0
+7.65404249467e-18 0.125
+0.25 0.0
+0.176776695297 0.176776695297
+1.53080849893e-17 0.25
+0.375 0.0
+0.324759526419 0.1875
+0.1875 0.324759526419
+2.2962127484e-17 0.375
+0.5 0.0
+0.461939766256 0.191341716183
+0.353553390593 0.353553390593
+0.191341716183 0.461939766256
+3.06161699787e-17 0.5
+0.625 0.0
+0.594410322684 0.193135621484
+0.505635621484 0.367365782683
+0.367365782683 0.505635621484
+0.193135621484 0.594410322684
+3.82702124734e-17 0.625
+0.75 0.0
+0.724444369717 0.194114283827
+0.649519052838 0.375
+0.53033008589 0.53033008589
+0.375 0.649519052838
+0.194114283827 0.724444369717
+4.5924254968e-17 0.75
+0.875 0.0
+0.853061923159 0.194705817212
+0.788347759415 0.379648271728
+0.68410254716 0.545553576626
+0.545553576626 0.68410254716
+0.379648271728 0.788347759415
+0.194705817212 0.853061923159
+5.35782974627e-17 0.875
+1.0 0.0
+0.980785280403 0.195090322016
+0.923879532511 0.382683432365
+0.831469612303 0.55557023302
+0.707106781187 0.707106781187
+0.55557023302 0.831469612303
+0.382683432365 0.923879532511
+0.195090322016 0.980785280403
+6.12323399574e-17 1.0

data/mixed-ele-meshes/tri_quad_2nd_order.mesh

@@ -0,0 +1,73 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+8
+1 2 0 1 2
+1 2 1 4 2
+1 2 1 3 4
+1 2 4 5 2
+2 3 3 6 7 4
+2 3 8 5 4 7
+2 3 6 9 10 7
+2 3 11 8 7 10
+
+boundary
+10
+1 1 0 1
+1 1 1 3
+1 1 3 6
+1 1 6 9
+1 1 9 10
+1 1 10 11
+1 1 11 8
+1 1 8 5
+1 1 5 2
+1 1 2 0
+
+vertices
+12
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: Quadratic
+VDim: 2
+Ordering: 1
+
+0.0 0.0
+0.5 0.0
+0.0 0.5
+1.0 0.0
+0.70711 0.70711
+0.0 1.0
+1.5 0.0
+1.35355 1.35355
+0.0 1.5
+2.0 0.0
+2.0 2.0
+0.0 2.0
+0.25 0.0
+0.35355 0.35355
+0.0 0.25
+0.64714 0.37909
+0.37909 0.64714
+0.75 0.0
+0.92388 0.38269
+0.38269 0.92388
+0.0 0.75
+1.25 0.0
+1.54238 0.73161
+1.03033 1.03033
+0.0 1.25
+0.73161 1.54238
+1.75 0.0
+2.0 1.0
+1.67678 1.67678
+0.0 1.75
+1.0 2.0
+1.14017 0.51516
+0.51516 1.14017
+1.71339 0.83839
+0.83839 1.71339

data/mixed-ele-meshes/whole_plane_tri_quad_2nd_order.mesh

@@ -0,0 +1,567 @@
+MFEM mesh v1.0
+
+dimension
+2
+
+elements
+128
+1 2 0 1 2
+1 2 1 4 2
+1 2 1 3 4
+1 2 2 4 5
+1 2 3 7 4
+1 2 3 6 7
+1 2 4 8 5
+1 2 4 7 8
+1 2 5 8 9
+1 2 6 11 7
+1 2 6 10 11
+1 2 7 12 8
+1 2 7 11 12
+1 2 8 13 9
+1 2 8 12 13
+1 2 9 13 14
+2 3 10 15 16 11
+2 3 11 16 17 12
+2 3 12 17 18 13
+2 3 13 18 19 14
+2 3 15 20 21 16
+2 3 16 21 22 17
+2 3 17 22 23 18
+2 3 18 23 24 19
+2 3 20 25 26 21
+2 3 21 26 27 22
+2 3 22 27 28 23
+2 3 23 28 29 24
+2 3 25 30 31 26
+2 3 26 31 32 27
+2 3 27 32 33 28
+2 3 28 33 34 29
+1 2 0 2 35
+1 2 35 2 37
+1 2 35 37 36
+1 2 2 5 37
+1 2 36 37 39
+1 2 36 39 38
+1 2 37 5 40
+1 2 37 40 39
+1 2 5 9 40
+1 2 38 39 42
+1 2 38 42 41
+1 2 39 40 43
+1 2 39 43 42
+1 2 40 9 44
+1 2 40 44 43
+1 2 9 14 44
+2 3 41 42 46 45
+2 3 42 43 47 46
+2 3 43 44 48 47
+2 3 44 14 19 48
+2 3 45 46 50 49
+2 3 46 47 51 50
+2 3 47 48 52 51
+2 3 48 19 24 52
+2 3 49 50 54 53
+2 3 50 51 55 54
+2 3 51 52 56 55
+2 3 52 24 29 56
+2 3 53 54 58 57
+2 3 54 55 59 58
+2 3 55 56 60 59
+2 3 56 29 34 60
+1 2 0 61 1
+1 2 1 61 62
+1 2 1 62 3
+1 2 61 63 62
+1 2 3 62 64
+1 2 3 64 6
+1 2 62 63 65
+1 2 62 65 64
+1 2 63 66 65
+1 2 6 64 67
+1 2 6 67 10
+1 2 64 65 68
+1 2 64 68 67
+1 2 65 66 69
+1 2 65 69 68
+1 2 66 70 69
+2 3 10 67 71 15
+2 3 67 68 72 71
+2 3 68 69 73 72
+2 3 69 70 74 73
+2 3 15 71 75 20
+2 3 71 72 76 75
+2 3 72 73 77 76
+2 3 73 74 78 77
+2 3 20 75 79 25
+2 3 75 76 80 79
+2 3 76 77 81 80
+2 3 77 78 82 81
+2 3 25 79 83 30
+2 3 79 80 84 83
+2 3 80 81 85 84
+2 3 81 82 86 85
+1 2 0 35 61
+1 2 35 87 61
+1 2 35 36 87
+1 2 61 87 63
+1 2 36 88 87
+1 2 36 38 88
+1 2 87 89 63
+1 2 87 88 89
+1 2 63 89 66
+1 2 38 90 88
+1 2 38 41 90
+1 2 88 91 89
+1 2 88 90 91
+1 2 89 92 66
+1 2 89 91 92
+1 2 66 92 70
+2 3 41 45 93 90
+2 3 90 93 94 91
+2 3 91 94 95 92
+2 3 92 95 74 70
+2 3 45 49 96 93
+2 3 93 96 97 94
+2 3 94 97 98 95
+2 3 95 98 78 74
+2 3 49 53 99 96
+2 3 96 99 100 97
+2 3 97 100 101 98
+2 3 98 101 82 78
+2 3 53 102 103 99
+2 3 99 103 104 100
+2 3 100 104 105 101
+2 3 101 105 86 82
+
+boundary
+16
+1 1 30 31
+1 1 31 32
+1 1 32 33
+1 1 33 34
+1 1 34 60
+1 1 60 59
+1 1 59 58
+1 1 58 57
+1 1 102 103
+1 1 103 104
+1 1 104 105
+1 1 105 86
+1 1 86 85
+1 1 85 84
+1 1 84 83
+1 1 83 30
+
+vertices
+106
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: H1_2D_P2
+VDim: 2
+Ordering: 1
+
+0.0 0.0
+0.25 0.0
+0.0 0.25
+0.5 0.0
+0.35355 0.35355
+0.0 0.5
+0.75 0.0
+0.64952 0.375
+0.375 0.64952
+0.0 0.75
+1.0 0.0
+0.92388 0.38268
+0.70711 0.70711
+0.38268 0.92388
+0.0 1.0
+1.25 0.0
+1.19291 0.53701
+1.03033 1.03033
+0.53701 1.19291
+0.0 1.25
+1.5 0.0
+1.46194 0.69134
+1.35355 1.35355
+0.69134 1.46194
+0.0 1.5
+1.75 0.0
+1.73097 0.84567
+1.67678 1.67678
+0.84567 1.73097
+0.0 1.75
+2.0 0.0
+2.0 1.0
+2.0 2.0
+1.0 2.0
+0.0 2.0
+-0.25 0.0
+-0.5 0.0
+-0.35355 0.35355
+-0.75 0.0
+-0.64952 0.375
+-0.375 0.64952
+-1.0 0.0
+-0.92388 0.38268
+-0.70711 0.70711
+-0.38268 0.92388
+-1.25 0.0
+-1.19291 0.53701
+-1.03033 1.03033
+-0.53701 1.19291
+-1.5 0.0
+-1.46194 0.69134
+-1.35355 1.35355
+-0.69134 1.46194
+-1.75 0.0
+-1.73097 0.84567
+-1.67678 1.67678
+-0.84567 1.73097
+-2.0 0.0
+-2.0 1.0
+-2.0 2.0
+-1.0 2.0
+0.0 -0.25
+0.35355 -0.35355
+0.0 -0.5
+0.64952 -0.375
+0.375 -0.64952
+0.0 -0.75
+0.92388 -0.38268
+0.70711 -0.70711
+0.38268 -0.92388
+0.0 -1.0
+1.19291 -0.53701
+1.03033 -1.03033
+0.53701 -1.19291
+0.0 -1.25
+1.46194 -0.69134
+1.35355 -1.35355
+0.69134 -1.46194
+0.0 -1.5
+1.73097 -0.84567
+1.67678 -1.67678
+0.84567 -1.73097
+0.0 -1.75
+2.0 -1.0
+2.0 -2.0
+1.0 -2.0
+0.0 -2.0
+-0.35355 -0.35355
+-0.64952 -0.375
+-0.375 -0.64952
+-0.92388 -0.38268
+-0.70711 -0.70711
+-0.38268 -0.92388
+-1.19291 -0.53701
+-1.03033 -1.03033
+-0.53701 -1.19291
+-1.46194 -0.69134
+-1.35355 -1.35355
+-0.69134 -1.46194
+-1.73097 -0.84567
+-1.67678 -1.67678
+-0.84567 -1.73097
+-2.0 -0.0
+-2.0 -1.0
+-2.0 -2.0
+-1.0 -2.0
+0.125 0.0
+0.17678 0.17678
+0.0 0.125
+0.32357 0.18954
+0.18954 0.32357
+0.375 0.0
+0.46194 0.19134
+0.19134 0.46194
+0.0 0.375
+0.59418 0.19384
+0.50569 0.36729
+0.625 0.0
+0.72444 0.19411
+0.36729 0.50569
+0.19384 0.59418
+0.53033 0.53033
+0.19411 0.72444
+0.0 0.625
+0.85299 0.19501
+0.78835 0.37964
+0.875 0.0
+0.98079 0.19509
+0.68405 0.54563
+0.54563 0.68405
+0.83147 0.55557
+0.37964 0.78835
+0.19501 0.85299
+0.55557 0.83147
+0.19509 0.98079
+0.0 0.875
+1.125 0.0
+1.24928 0.27462
+1.05851 0.4599
+1.13007 0.79668
+0.86872 0.86872
+0.79668 1.13007
+0.4599 1.05851
+0.27462 1.24928
+0.0 1.125
+1.375 0.0
+1.51779 0.35426
+1.32748 0.6142
+1.42864 1.03762
+1.19194 1.19194
+1.03762 1.42864
+0.6142 1.32748
+0.35426 1.51779
+0.0 1.375
+1.625 0.0
+1.78629 0.43396
+1.59649 0.76852
+1.72721 1.2785
+1.51516 1.51516
+1.2785 1.72721
+0.76852 1.59649
+0.43396 1.78629
+0.0 1.625
+1.875 0.0
+2.0 0.5
+1.8655 0.92284
+2.0 1.5
+1.83839 1.83839
+1.5 2.0
+0.92284 1.8655
+0.5 2.0
+0.0 1.875
+-0.17678 0.17678
+-0.125 0.0
+-0.18954 0.32357
+-0.32357 0.18954
+-0.46194 0.19134
+-0.375 0.0
+-0.19134 0.46194
+-0.50569 0.36729
+-0.59418 0.19384
+-0.72444 0.19411
+-0.625 0.0
+-0.19384 0.59418
+-0.36729 0.50569
+-0.53033 0.53033
+-0.19411 0.72444
+-0.78835 0.37964
+-0.85299 0.19501
+-0.98079 0.19509
+-0.875 0.0
+-0.54563 0.68405
+-0.68405 0.54563
+-0.83147 0.55557
+-0.19501 0.85299
+-0.37964 0.78835
+-0.55557 0.83147
+-0.19509 0.98079
+-1.05851 0.4599
+-1.24928 0.27462
+-1.125 0.0
+-0.86872 0.86872
+-1.13007 0.79668
+-0.4599 1.05851
+-0.79668 1.13007
+-0.27462 1.24928
+-1.32748 0.6142
+-1.51779 0.35426
+-1.375 0.0
+-1.19194 1.19194
+-1.42864 1.03762
+-0.6142 1.32748
+-1.03762 1.42864
+-0.35426 1.51779
+-1.59649 0.76852
+-1.78629 0.43396
+-1.625 0.0
+-1.51516 1.51516
+-1.72721 1.2785
+-0.76852 1.59649
+-1.2785 1.72721
+-0.43396 1.78629
+-1.8655 0.92284
+-2.0 0.5
+-1.875 0.0
+-1.83839 1.83839
+-2.0 1.5
+-0.92284 1.8655
+-1.5 2.0
+-0.5 2.0
+0.0 -0.125
+0.17678 -0.17678
+0.18954 -0.32357
+0.32357 -0.18954
+0.46194 -0.19134
+0.0 -0.375
+0.19134 -0.46194
+0.50569 -0.36729
+0.59418 -0.19384
+0.72444 -0.19411
+0.19384 -0.59418
+0.36729 -0.50569
+0.53033 -0.53033
+0.0 -0.625
+0.19411 -0.72444
+0.78835 -0.37964
+0.85299 -0.19501
+0.98079 -0.19509
+0.54563 -0.68405
+0.68405 -0.54563
+0.83147 -0.55557
+0.19501 -0.85299
+0.37964 -0.78835
+0.55557 -0.83147
+0.0 -0.875
+0.19509 -0.98079
+1.05851 -0.4599
+1.24928 -0.27462
+0.86872 -0.86872
+1.13007 -0.79668
+0.4599 -1.05851
+0.79668 -1.13007
+0.0 -1.125
+0.27462 -1.24928
+1.32748 -0.6142
+1.51779 -0.35426
+1.19194 -1.19194
+1.42864 -1.03762
+0.6142 -1.32748
+1.03762 -1.42864
+0.0 -1.375
+0.35426 -1.51779
+1.59649 -0.76852
+1.78629 -0.43396
+1.51516 -1.51516
+1.72721 -1.2785
+0.76852 -1.59649
+1.2785 -1.72721
+0.0 -1.625
+0.43396 -1.78629
+1.8655 -0.92284
+2.0 -0.5
+1.83839 -1.83839
+2.0 -1.5
+0.92284 -1.8655
+1.5 -2.0
+0.0 -1.875
+0.5 -2.0
+-0.17678 -0.17678
+-0.32357 -0.18954
+-0.18954 -0.32357
+-0.46194 -0.19134
+-0.19134 -0.46194
+-0.59418 -0.19384
+-0.50569 -0.36729
+-0.72444 -0.19411
+-0.36729 -0.50569
+-0.19384 -0.59418
+-0.53033 -0.53033
+-0.19411 -0.72444
+-0.85299 -0.19501
+-0.78835 -0.37964
+-0.98079 -0.19509
+-0.68405 -0.54563
+-0.54563 -0.68405
+-0.83147 -0.55557
+-0.37964 -0.78835
+-0.19501 -0.85299
+-0.55557 -0.83147
+-0.19509 -0.98079
+-1.24928 -0.27462
+-1.05851 -0.4599
+-1.13007 -0.79668
+-0.86872 -0.86872
+-0.79668 -1.13007
+-0.4599 -1.05851
+-0.27462 -1.24928
+-1.51779 -0.35426
+-1.32748 -0.6142
+-1.42864 -1.03762
+-1.19194 -1.19194
+-1.03762 -1.42864
+-0.6142 -1.32748
+-0.35426 -1.51779
+-1.78629 -0.43396
+-1.59649 -0.76852
+-1.72721 -1.2785
+-1.51516 -1.51516
+-1.2785 -1.72721
+-0.76852 -1.59649
+-0.43396 -1.78629
+-1.875 0.0
+-2.0 -0.5
+-1.8655 -0.92284
+-2.0 -1.5
+-1.83839 -1.83839
+-1.5 -2.0
+-0.92284 -1.8655
+-0.5 -2.0
+1.0917 0.22992
+0.96356 0.66428
+0.66428 0.96356
+0.22992 1.0917
+1.35121 0.30709
+1.25968 0.90306
+0.90306 1.25968
+0.30709 1.35121
+1.61073 0.38425
+1.55581 1.14184
+1.14184 1.55581
+0.38425 1.61073
+1.87024 0.46142
+1.85194 1.38061
+1.38061 1.85194
+0.46142 1.87024
+-1.0917 0.22992
+-0.96356 0.66428
+-0.66428 0.96356
+-0.22992 1.0917
+-1.35121 0.30709
+-1.25968 0.90306
+-0.90306 1.25968
+-0.30709 1.35121
+-1.61073 0.38425
+-1.55581 1.14184
+-1.14184 1.55581
+-0.38425 1.61073
+-1.87024 0.46142
+-1.85194 1.38061
+-1.38061 1.85194
+-0.46142 1.87024
+1.0917 -0.22992
+0.96356 -0.66428
+0.66428 -0.96356
+0.22992 -1.0917
+1.35121 -0.30709
+1.25968 -0.90306
+0.90306 -1.25968
+0.30709 -1.35121
+1.61073 -0.38425
+1.55581 -1.14184
+1.14184 -1.55581
+0.38425 -1.61073
+1.87024 -0.46142
+1.85194 -1.38061
+1.38061 -1.85194
+0.46142 -1.87024
+-1.0917 -0.22992
+-0.96356 -0.66428
+-0.66428 -0.96356
+-0.22992 -1.0917
+-1.35121 -0.30709
+-1.25968 -0.90306
+-0.90306 -1.25968
+-0.30709 -1.35121
+-1.61073 -0.38425
+-1.55581 -1.14184
+-1.14184 -1.55581
+-0.38425 -1.61073
+-1.87024 -0.46142
+-1.85194 -1.38061
+-1.38061 -1.85194
+-0.46142 -1.87024

examples/ex1.cpp

@@ -80,8 +80,8 @@ int main(int argc, char *argv[])
    //    largest number that gives a final mesh with no more than 50,000
    //    elements.
    {
-      int ref_levels =
-         (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
+     int ref_levels = 0; 
+       //(int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh->UniformRefinement();

examples/ex13p.cpp

@@ -110,6 +110,7 @@ int main(int argc, char *argv[])
    {
       pmesh->UniformRefinement();
    }
+   pmesh->ReorientTetMesh();
 
    // 6. Define a parallel finite element space on the parallel mesh. Here we
    //    use the Nedelec finite elements of the specified order.

fem/intrules.hpp

@@ -31,6 +31,7 @@ public:
 
    void Set(const double *p, const int dim)
    {
+      MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
       x = p[0];
       if (dim > 1)
       {
@@ -44,6 +45,7 @@ public:
 
    void Get(double *p, const int dim) const
    {
+      MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
       p[0] = x;
       if (dim > 1)
       {

general/sort_pairs.hpp

@@ -58,6 +58,11 @@ public:
    A one;
    B two;
    C three;
+
+   Triple() { }
+
+   Triple(const A &one, const B &two, const C &three)
+      : one(one), two(two), three(three) { }
 };
 
 /// @brief Lexicographic comparison operator for class Triple.

linalg/sparsemat.hpp

@@ -406,9 +406,6 @@ public:
    /// Print various sparse matrix staticstics.
    void PrintInfo(std::ostream &out) const;
 
-   /// Walks the sparse matrix
-   int Walk(int &i, int &j, double &a);
-
    /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
    double IsSymmetric() const;
 

mesh/element.hpp

@@ -17,6 +17,7 @@
 #include "../general/table.hpp"
 #include "../linalg/densemat.hpp"
 #include "../fem/geom.hpp"
+#include "../general/hash.hpp"
 
 namespace mfem
 {
@@ -80,7 +81,7 @@ public:
    virtual void MarkEdge(const DSTable &v_to_v, const int *length) {}
 
    /// Return 1 if the element needs refinement in order to get conforming mesh.
-   virtual int NeedRefinement(DSTable &v_to_v, int *middle) const { return 0; }
+   virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const { return 0; }
 
    /// Set current coarse-fine transformation number.
    virtual void ResetTransform(int tr) {}

mesh/mesh.cpp

@@ -5905,56 +5905,47 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
    }
    else if (Dim == 3) // ---------------------------------------------------
    {
-      // 1. Get table of vertex to vertex connections.
-      DSTable v_to_v(NumOfVertices);
-      GetVertexToVertexTable(v_to_v);
+      // 1. Hash table of vertex to vertex connections corresponding to refined
+      //    edges.
+      HashTable<Hashed2> v_to_v;
 
-      // 2. Get edge to element connections in arrays edge1 and edge2
-      nedges = v_to_v.NumberOfEntries();
-      int *middle = new int[nedges];
-
-      for (i = 0; i < nedges; i++)
-      {
-         middle[i] = -1;
-      }
-
-      // 3. Do the red refinement.
+      // 2. Do the red refinement.
       int ii;
       switch (type)
       {
          case 1:
             for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
             }
             break;
          case 2:
             for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
 
-               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(NumOfElements - 1, v_to_v);
+               Bisection(marked_el[i], v_to_v);
             }
             break;
          case 3:
             for (i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
 
                ii = NumOfElements - 1;
-               Bisection(ii, v_to_v, NULL, NULL, middle);
-               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
-               Bisection(ii, v_to_v, NULL, NULL, middle);
+               Bisection(ii, v_to_v);
+               Bisection(NumOfElements - 1, v_to_v);
+               Bisection(ii, v_to_v);
 
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
-               Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
+               Bisection(NumOfElements-1, v_to_v);
+               Bisection(marked_el[i], v_to_v);
             }
             break;
       }
 
-      // 4. Do the green refinement (to get conforming mesh).
+      // 3. Do the green refinement (to get conforming mesh).
       int need_refinement;
       // int need_refinement, onoe, max_gen = 0;
       do
@@ -5968,10 +5959,10 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
             // ((Tetrahedron *)elements[i])->
             // ParseRefinementFlag(redges, type, flag);
             // if (flag > max_gen)  max_gen = flag;
-            if (elements[i]->NeedRefinement(v_to_v, middle))
+            if (elements[i]->NeedRefinement(v_to_v))
             {
                need_refinement = 1;
-               Bisection(i, v_to_v, NULL, NULL, middle);
+               Bisection(i, v_to_v);
             }
          }
       }
@@ -5979,38 +5970,23 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       // mfem::out << "Maximum generation: " << max_gen << endl;
 
-      // 5. Update the boundary elements.
+      // 4. Update the boundary elements.
       do
       {
          need_refinement = 0;
          for (i = 0; i < NumOfBdrElements; i++)
-            if (boundary[i]->NeedRefinement(v_to_v, middle))
+            if (boundary[i]->NeedRefinement(v_to_v))
             {
                need_refinement = 1;
-               Bisection(i, v_to_v, middle);
+               BdrBisection(i, v_to_v);
             }
       }
       while (need_refinement == 1);
 
-      // 6. Un-mark the Pf elements.
-      int refinement_edges[2], type, flag;
-      for (i = 0; i < NumOfElements; i++)
-      {
-         Tetrahedron* el = (Tetrahedron*) elements[i];
-         el->ParseRefinementFlag(refinement_edges, type, flag);
-
-         if (type == Tetrahedron::TYPE_PF)
-         {
-            el->CreateRefinementFlag(refinement_edges, Tetrahedron::TYPE_PU,
-                                     flag);
-         }
-      }
-
+      NumOfVertices = vertices.Size();
       NumOfBdrElements = boundary.Size();
 
-      // 7. Free the allocated memory.
-      delete [] middle;
-
+      // 5. Update element-to-edge and element-to-face relations.
       DeleteLazyTables();
       if (el_to_edge != NULL)
       {
@@ -6588,7 +6564,21 @@ void Mesh::Bisection(int i, const DSTable &v_to_v,
       }
       NumOfElements++;
    }
-   else if (t == Element::TETRAHEDRON)
+   else
+   {
+      MFEM_ABORT("Bisection for now works only for triangles.");
+   }
+}
+
+void Mesh::Bisection(int i, HashTable<Hashed2> &v_to_v)
+{
+   int *vert;
+   int v[2][4], v_new, bisect, t;
+   Element *el = elements[i];
+   Vertex V;
+
+   t = el->GetType();
+   if (t == Element::TETRAHEDRON)
    {
       int j, type, new_type, old_redges[2], new_redges[2][2], flag;
       Tetrahedron *tet = (Tetrahedron *) el;
@@ -6599,32 +6589,19 @@ void Mesh::Bisection(int i, const DSTable &v_to_v,
       vert = tet->GetVertices();
 
       // 1. Get the index for the new vertex in v_new.
-      bisect = v_to_v(vert[0], vert[1]);
+      bisect = v_to_v.FindId(vert[0], vert[1]);
       if (bisect == -1)
       {
-         tet->ParseRefinementFlag(old_redges, type, flag);
-         mfem::err << "Error in Bisection(...) of tetrahedron!" << endl
-                   << "   redge[0] = " << old_redges[0]
-                   << "   redge[1] = " << old_redges[1]
-                   << "   type = " << type
-                   << "   flag = " << flag << endl;
-         mfem_error();
-      }
-
-      if (middle[bisect] == -1)
-      {
-         v_new = NumOfVertices++;
+         v_new = NumOfVertices + v_to_v.GetId(vert[0],vert[1]);
          for (j = 0; j < 3; j++)
          {
             V(j) = 0.5 * (vertices[vert[0]](j) + vertices[vert[1]](j));
          }
          vertices.Append(V);
-
-         middle[bisect] = v_new;
       }
       else
       {
-         v_new = middle[bisect];
+         v_new = NumOfVertices + bisect;
       }
 
       // 2. Set the node indices for the new elements in v[2][4] so that
@@ -6708,11 +6685,11 @@ void Mesh::Bisection(int i, const DSTable &v_to_v,
    }
    else
    {
-      MFEM_ABORT("Bisection for now works only for triangles & tetrahedra.");
+      MFEM_ABORT("Bisection with HashTable for now works only for tetrahedra.");
    }
 }
 
-void Mesh::Bisection(int i, const DSTable &v_to_v, int *middle)
+void Mesh::BdrBisection(int i, const HashTable<Hashed2> &v_to_v)
 {
    int *vert;
    int v[2][3], v_new, bisect, t;
@@ -6726,9 +6703,9 @@ void Mesh::Bisection(int i, const DSTable &v_to_v, int *middle)
       vert = tri->GetVertices();
 
       // 1. Get the index for the new vertex in v_new.
-      bisect = v_to_v(vert[0], vert[1]);
+      bisect = v_to_v.FindId(vert[0], vert[1]);
       MFEM_ASSERT(bisect >= 0, "");
-      v_new = middle[bisect];
+      v_new = NumOfVertices + bisect;
       MFEM_ASSERT(v_new != -1, "");
 
       // 2. Set the node indices for the new elements in v[0] and v[1] so that
@@ -6744,7 +6721,8 @@ void Mesh::Bisection(int i, const DSTable &v_to_v, int *middle)
    }
    else
    {
-      MFEM_ABORT("Bisection of boundary elements works only for triangles!");
+      MFEM_ABORT("Bisection of boundary elements with HashTable works only for"
+                 " triangles!");
    }
 }
 

mesh/mesh.hpp

@@ -250,11 +250,14 @@ protected:
                         int *edge1, int *edge2, int *middle)
    { Bisection(i, v_to_v, edge1, edge2, middle); }
 
-   /** Bisection. Element with index i is bisected. */
+   /// Bisect a triangle: element with index @a i is bisected.
    void Bisection(int i, const DSTable &, int *, int *, int *);
 
-   /** Bisection. Boundary element with index i is bisected. */
-   void Bisection(int i, const DSTable &, int *);
+   /// Bisect a tetrahedron: element with index @a i is bisected.
+   void Bisection(int i, HashTable<Hashed2> &);
+
+   /// Bisect a boundary triangle: boundary element with index @a i is bisected.
+   void BdrBisection(int i, const HashTable<Hashed2> &);
 
    /** Uniform Refinement. Element with index i is refined uniformly. */
    void UniformRefinement(int i, const DSTable &, int *, int *, int *);

mesh/pmesh.cpp

@@ -1214,48 +1214,76 @@ int ParMesh::GetEdgeSplittings(Element *edge, const DSTable &v_to_v,
    }
 }
 
-// For a triangular face with (correctly ordered) vertices v[0], v[1], v[2]
-// return a number with the following meaning:
-// 0 - the face was not refined
-// 1 - the face was refined once by splitting v[0],v[1]
-// 2 - the face was refined twice by splitting v[0],v[1] and then v[1],v[2]
-// 3 - the face was refined twice by splitting v[0],v[1] and then v[0],v[2]
-// 4 - the face was refined three times (as in 2+3)
-int ParMesh::GetFaceSplittings(Element *face, const DSTable &v_to_v,
-                               int *middle)
+void ParMesh::GetFaceSplittings(Element *face, const HashTable<Hashed2> &v_to_v,
+                                Array<unsigned> &codes)
 {
-   int m, right = 0;
-   int number_of_splittings = 0;
-   int *v = face->GetVertices();
+   const int *v = face->GetVertices();
+   typedef Triple<int,int,int> face_t;
+   Array<face_t> face_stack;
 
-   if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
+   unsigned code = 0;
+   face_stack.Append(face_t(v[0], v[1], v[2]));
+   for (unsigned bit = 0; face_stack.Size() > 0; bit++)
    {
-      number_of_splittings++;
-      if ((m = v_to_v(v[1], v[2])) != -1 && middle[m] != -1)
+      if (bit == 8*sizeof(unsigned))
       {
-         right = 1;
-         number_of_splittings++;
-      }
-      if ((m = v_to_v(v[2], v[0])) != -1 && middle[m] != -1)
-      {
-         number_of_splittings++;
+         codes.Append(code);
+         code = bit = 0;
       }
 
-      switch (number_of_splittings)
+      const face_t &f = face_stack.Last();
+      int mid = v_to_v.FindId(f.one, f.two);
+      if (mid == -1)
       {
-         case 2:
-            if (right == 0)
-            {
-               number_of_splittings++;
-            }
-            break;
-         case 3:
-            number_of_splittings++;
-            break;
+         // leave a 0 at bit 'bit'
+         face_stack.DeleteLast();
+      }
+      else
+      {
+         code += (1 << bit); // set bit 'bit' to 1
+         mid += NumOfVertices;
+         face_stack.Append(face_t(f.three, f.one, mid));
+         face_t &r = face_stack[face_stack.Size()-2];
+         r = face_t(r.two, r.three, mid);
       }
    }
+   codes.Append(code);
+}
 
-   return number_of_splittings;
+bool ParMesh::DecodeFaceSplittings(HashTable<Hashed2> &v_to_v, const int *v,
+                                   const Array<unsigned> &codes, int &pos)
+{
+   typedef Triple<int,int,int> face_t;
+   Array<face_t> face_stack;
+
+   bool need_refinement = 0;
+   face_stack.Append(face_t(v[0], v[1], v[2]));
+   for (unsigned bit = 0, code = codes[pos++]; face_stack.Size() > 0; bit++)
+   {
+      if (bit == 8*sizeof(unsigned))
+      {
+         code = codes[pos++];
+         bit = 0;
+      }
+
+      if ((code & (1 << bit)) == 0) { face_stack.DeleteLast(); continue; }
+
+      const face_t &f = face_stack.Last();
+      int mid = v_to_v.FindId(f.one, f.two);
+      if (mid == -1)
+      {
+         mid = v_to_v.GetId(f.one, f.two);
+         int ind[2] = { f.one, f.two };
+         vertices.Append(Vertex());
+         AverageVertices(ind, 2, vertices.Size()-1);
+         need_refinement = 1;
+      }
+      mid += NumOfVertices;
+      face_stack.Append(face_t(f.three, f.one, mid));
+      face_t &r = face_stack[face_stack.Size()-2];
+      r = face_t(r.two, r.three, mid);
+   }
+   return need_refinement;
 }
 
 void ParMesh::GenerateOffsets(int N, HYPRE_Int loc_sizes[],
@@ -2113,8 +2141,6 @@ void ParMesh::ReorientTetMesh()
 
 void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
 {
-   int i, j;
-
    if (pncmesh)
    {
       MFEM_ABORT("Local and nonconforming refinements cannot be mixed.");
@@ -2133,76 +2159,63 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
          uniform_refinement = 1;
       }
 
-      // 1. Get table of vertex to vertex connections.
-      DSTable v_to_v(NumOfVertices);
-      GetVertexToVertexTable(v_to_v);
+      // 1. Hash table of vertex to vertex connections corresponding to refined
+      //    edges.
+      HashTable<Hashed2> v_to_v;
 
-      // 2. Create a marker array for all edges (vertex to vertex connections).
-      Array<int> middle(v_to_v.NumberOfEntries());
-      middle = -1;
-
-      // 3. Do the red refinement.
+      // 2. Do the red refinement.
       switch (type)
       {
          case 1:
-            for (i = 0; i < marked_el.Size(); i++)
+            for (int i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
             }
             break;
          case 2:
-            for (i = 0; i < marked_el.Size(); i++)
+            for (int i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
 
-               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(NumOfElements - 1, v_to_v);
+               Bisection(marked_el[i], v_to_v);
             }
             break;
          case 3:
-            for (i = 0; i < marked_el.Size(); i++)
+            for (int i = 0; i < marked_el.Size(); i++)
             {
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
 
-               j = NumOfElements - 1;
-               Bisection(j, v_to_v, NULL, NULL, middle);
-               Bisection(NumOfElements - 1, v_to_v, NULL, NULL, middle);
-               Bisection(j, v_to_v, NULL, NULL, middle);
+               int j = NumOfElements - 1;
+               Bisection(j, v_to_v);
+               Bisection(NumOfElements - 1, v_to_v);
+               Bisection(j, v_to_v);
 
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
-               Bisection(NumOfElements-1, v_to_v, NULL, NULL, middle);
-               Bisection(marked_el[i], v_to_v, NULL, NULL, middle);
+               Bisection(marked_el[i], v_to_v);
+               Bisection(NumOfElements-1, v_to_v);
+               Bisection(marked_el[i], v_to_v);
             }
             break;
       }
 
-      // 4. Do the green refinement (to get conforming mesh).
+      // 3. Do the green refinement (to get conforming mesh).
       int need_refinement;
-      int refined_edge[5][3] =
-      {
-         {0, 0, 0},
-         {1, 0, 0},
-         {1, 1, 0},
-         {1, 0, 1},
-         {1, 1, 1}
-      };
-      int faces_in_group, max_faces_in_group = 0;
+      int max_faces_in_group = 0;
       // face_splittings identify how the shared faces have been split
-      int **face_splittings = new int*[GetNGroups()-1];
-      for (i = 0; i < GetNGroups()-1; i++)
+      Array<unsigned> *face_splittings = new Array<unsigned>[GetNGroups()-1];
+      for (int i = 0; i < GetNGroups()-1; i++)
       {
-         faces_in_group = GroupNFaces(i+1);
-         face_splittings[i] = new int[faces_in_group];
+         const int faces_in_group = GroupNFaces(i+1);
+         face_splittings[i].Reserve(faces_in_group);
          if (faces_in_group > max_faces_in_group)
          {
             max_faces_in_group = faces_in_group;
          }
       }
-      int neighbor, *iBuf = new int[max_faces_in_group];
+      int neighbor;
+      Array<unsigned> iBuf(max_faces_in_group);
 
-      Array<int> group_faces;
-
-      MPI_Request request;
+      MPI_Request *requests = new MPI_Request[GetNGroups()-1];
       MPI_Status  status;
 
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
@@ -2211,12 +2224,12 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       do
       {
          need_refinement = 0;
-         for (i = 0; i < NumOfElements; i++)
+         for (int i = 0; i < NumOfElements; i++)
          {
-            if (elements[i]->NeedRefinement(v_to_v, middle))
+            if (elements[i]->NeedRefinement(v_to_v))
             {
                need_refinement = 1;
-               Bisection(i, v_to_v, NULL, NULL, middle);
+               Bisection(i, v_to_v);
             }
          }
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
@@ -2239,96 +2252,85 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
             const int tag = 293;
 
             // (a) send the type of interface splitting
-            for (i = 0; i < GetNGroups()-1; i++)
+            int req_count = 0;
+            for (int i = 0; i < GetNGroups()-1; i++)
             {
-               group_sface.GetRow(i, group_faces);
-               faces_in_group = group_faces.Size();
+               const int *group_faces = group_sface.GetRow(i);
+               const int faces_in_group = group_sface.RowSize(i);
                // it is enough to communicate through the faces
                if (faces_in_group == 0) { continue; }
 
-               for (j = 0; j < faces_in_group; j++)
+               face_splittings[i].SetSize(0);
+               for (int j = 0; j < faces_in_group; j++)
                {
-                  face_splittings[i][j] =
-                     GetFaceSplittings(shared_faces[group_faces[j]], v_to_v,
-                                       middle);
+                  GetFaceSplittings(shared_faces[group_faces[j]], v_to_v,
+                                    face_splittings[i]);
                }
                const int *nbs = gtopo.GetGroup(i+1);
                neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
-               MPI_Isend(face_splittings[i], faces_in_group, MPI_INT,
-                         neighbor, tag, MyComm, &request);
+               MPI_Isend(face_splittings[i], face_splittings[i].Size(),
+                         MPI_UNSIGNED, neighbor, tag, MyComm,
+                         &requests[req_count++]);
             }
 
             // (b) receive the type of interface splitting
-            for (i = 0; i < GetNGroups()-1; i++)
+            for (int i = 0; i < GetNGroups()-1; i++)
             {
-               group_sface.GetRow(i, group_faces);
-               faces_in_group = group_faces.Size();
+               const int *group_faces = group_sface.GetRow(i);
+               const int faces_in_group = group_sface.RowSize(i);
                if (faces_in_group == 0) { continue; }
 
                const int *nbs = gtopo.GetGroup(i+1);
                neighbor = gtopo.GetNeighborRank(nbs[0] ? nbs[0] : nbs[1]);
-               MPI_Recv(iBuf, faces_in_group, MPI_INT, neighbor,
-                        tag, MyComm, &status);
+               MPI_Probe(neighbor, tag, MyComm, &status);
+               int count;
+               MPI_Get_count(&status, MPI_UNSIGNED, &count);
+               iBuf.SetSize(count);
+               MPI_Recv(iBuf, count, MPI_UNSIGNED, neighbor, tag, MyComm,
+                        MPI_STATUS_IGNORE);
 
-               for (j = 0; j < faces_in_group; j++)
+               for (int j = 0, pos = 0; j < faces_in_group; j++)
                {
-                  if (iBuf[j] == face_splittings[i][j]) { continue; }
-
-                  int *v = shared_faces[group_faces[j]]->GetVertices();
-                  for (int k = 0; k < 3; k++)
-                  {
-                     if (refined_edge[iBuf[j]][k] != 1 ||
-                         refined_edge[face_splittings[i][j]][k] != 0)
-                     { continue; }
-
-                     int ind[2] = { v[k], v[(k+1)%3] };
-                     int ii = v_to_v(ind[0], ind[1]);
-                     if (middle[ii] == -1)
-                     {
-                        need_refinement = 1;
-                        middle[ii] = NumOfVertices++;
-                        vertices.Append(Vertex());
-                        AverageVertices(ind, 2, vertices.Size()-1);
-                     }
-                  }
+                  const int *v = shared_faces[group_faces[j]]->GetVertices();
+                  need_refinement |= DecodeFaceSplittings(v_to_v, v, iBuf, pos);
                }
             }
 
-            i = need_refinement;
-            MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
+            int nr = need_refinement;
+            MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
+
+            MPI_Waitall(req_count, requests, MPI_STATUSES_IGNORE);
          }
       }
       while (need_refinement == 1);
 
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
-      i = ref_loops_all;
-      MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
-      if (MyRank == 0)
       {
-         mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
-                   << ref_loops_all << ", ref_loops_par = " << ref_loops_par
-                   << '\n' << endl;
+         int i = ref_loops_all;
+         MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
+         if (MyRank == 0)
+         {
+            mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
+                      << ref_loops_all << ", ref_loops_par = " << ref_loops_par
+                      << '\n' << endl;
+         }
       }
 #endif
 
-      delete [] iBuf;
-      for (i = 0; i < GetNGroups()-1; i++)
-      {
-         delete [] face_splittings[i];
-      }
+      delete [] requests;
+      iBuf.DeleteAll();
       delete [] face_splittings;
 
-
-      // 5. Update the boundary elements.
+      // 4. Update the boundary elements.
       do
       {
          need_refinement = 0;
-         for (i = 0; i < NumOfBdrElements; i++)
+         for (int i = 0; i < NumOfBdrElements; i++)
          {
-            if (boundary[i]->NeedRefinement(v_to_v, middle))
+            if (boundary[i]->NeedRefinement(v_to_v))
             {
                need_refinement = 1;
-               Bisection(i, v_to_v, middle);
+               BdrBisection(i, v_to_v);
             }
          }
       }
@@ -2342,31 +2344,16 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       DeleteLazyTables();
 
-      // 5a. Update the groups after refinement.
+      // 5. Update the groups after refinement.
       if (el_to_face != NULL)
       {
-         RefineGroups(v_to_v, middle);
+         RefineGroups(v_to_v);
          // GetElementToFaceTable(); // Called by RefineGroups
          GenerateFaces();
       }
+      NumOfVertices = vertices.Size();
 
-      // 6. Un-mark the Pf elements.
-      int refinement_edges[2], type, flag;
-      for (i = 0; i < NumOfElements; i++)
-      {
-         Tetrahedron* el = (Tetrahedron*) elements[i];
-         el->ParseRefinementFlag(refinement_edges, type, flag);
-
-         if (type == Tetrahedron::TYPE_PF)
-         {
-            el->CreateRefinementFlag(refinement_edges, Tetrahedron::TYPE_PU,
-                                     flag);
-         }
-      }
-
-      // 7. Free the allocated memory.
-      middle.DeleteAll();
-
+      // 6. Update element-to-edge relations.
       if (el_to_edge != NULL)
       {
          NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
@@ -2393,15 +2380,15 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       int *edge2  = new int[nedges];
       int *middle = new int[nedges];
 
-      for (i = 0; i < nedges; i++)
+      for (int i = 0; i < nedges; i++)
       {
          edge1[i] = edge2[i] = middle[i] = -1;
       }
 
-      for (i = 0; i < NumOfElements; i++)
+      for (int i = 0; i < NumOfElements; i++)
       {
          int *v = elements[i]->GetVertices();
-         for (j = 0; j < 3; j++)
+         for (int j = 0; j < 3; j++)
          {
             int ind = v_to_v(v[j], v[(j+1)%3]);
             (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
@@ -2409,7 +2396,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       }
 
       // 3. Do the red refinement.
-      for (i = 0; i < marked_el.Size(); i++)
+      for (int i = 0; i < marked_el.Size(); i++)
       {
          RedRefinement(marked_el[i], v_to_v, edge1, edge2, middle);
       }
@@ -2419,7 +2406,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       int edges_in_group, max_edges_in_group = 0;
       // edge_splittings identify how the shared edges have been split
       int **edge_splittings = new int*[GetNGroups()-1];
-      for (i = 0; i < GetNGroups()-1; i++)
+      for (int i = 0; i < GetNGroups()-1; i++)
       {
          edges_in_group = GroupNEdges(i+1);
          edge_splittings[i] = new int[edges_in_group];
@@ -2443,12 +2430,14 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       do
       {
          need_refinement = 0;
-         for (i = 0; i < nedges; i++)
+         for (int i = 0; i < nedges; i++)
+         {
             if (middle[i] != -1 && edge1[i] != -1)
             {
                need_refinement = 1;
                GreenRefinement(edge1[i], v_to_v, edge1, edge2, middle);
             }
+         }
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
          ref_loops_all++;
 #endif
@@ -2468,14 +2457,14 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
             // MPI_Barrier(MyComm);
 
             // (a) send the type of interface splitting
-            for (i = 0; i < GetNGroups()-1; i++)
+            for (int i = 0; i < GetNGroups()-1; i++)
             {
                group_sedge.GetRow(i, group_edges);
                edges_in_group = group_edges.Size();
                // it is enough to communicate through the edges
                if (edges_in_group != 0)
                {
-                  for (j = 0; j < edges_in_group; j++)
+                  for (int j = 0; j < edges_in_group; j++)
                   {
                      edge_splittings[i][j] =
                         GetEdgeSplittings(shared_edges[group_edges[j]], v_to_v,
@@ -2496,7 +2485,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
             }
 
             // (b) receive the type of interface splitting
-            for (i = 0; i < GetNGroups()-1; i++)
+            for (int i = 0; i < GetNGroups()-1; i++)
             {
                group_sedge.GetRow(i, group_edges);
                edges_in_group = group_edges.Size();
@@ -2514,15 +2503,18 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
                   MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
                            MPI_ANY_TAG, MyComm, &status);
 
-                  for (j = 0; j < edges_in_group; j++)
+                  for (int j = 0; j < edges_in_group; j++)
+                  {
                      if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
                      {
                         int *v = shared_edges[group_edges[j]]->GetVertices();
                         int ii = v_to_v(v[0], v[1]);
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
                         if (middle[ii] != -1)
+                        {
                            mfem_error("ParMesh::LocalRefinement (triangles) : "
                                       "Oops!");
+                        }
 #endif
                         need_refinement = 1;
                         middle[ii] = NumOfVertices++;
@@ -2532,27 +2524,30 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
                         }
                         vertices.Append(V);
                      }
+                  }
                }
             }
 
-            i = need_refinement;
-            MPI_Allreduce(&i, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
+            int nr = need_refinement;
+            MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR, MyComm);
          }
       }
       while (need_refinement == 1);
 
 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
-      i = ref_loops_all;
-      MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
-      if (MyRank == 0)
       {
-         mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
-                   << ref_loops_all << ", ref_loops_par = " << ref_loops_par
-                   << '\n' << endl;
+         int i = ref_loops_all;
+         MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0, MyComm);
+         if (MyRank == 0)
+         {
+            mfem::out << "\n\nParMesh::LocalRefinement : max. ref_loops_all = "
+                      << ref_loops_all << ", ref_loops_par = " << ref_loops_par
+                      << '\n' << endl;
+         }
       }
 #endif
 
-      for (i = 0; i < GetNGroups()-1; i++)
+      for (int i = 0; i < GetNGroups()-1; i++)
       {
          delete [] edge_splittings[i];
       }
@@ -2563,7 +2558,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       // 5. Update the boundary elements.
       int v1[2], v2[2], bisect, temp;
       temp = NumOfBdrElements;
-      for (i = 0; i < temp; i++)
+      for (int i = 0; i < temp; i++)
       {
          int *v = boundary[i]->GetVertices();
          bisect = v_to_v(v[0], v[1]);
@@ -2579,8 +2574,10 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
                boundary.Append(new Segment(v2, boundary[i]->GetAttribute()));
             }
             else
+            {
                mfem_error("Only bisection of segment is implemented for bdr"
                           " elem.");
+            }
          }
       }
       NumOfBdrElements = boundary.Size();
@@ -2611,9 +2608,9 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       elements.SetSize(NumOfElements);
       CoarseFineTr.embeddings.SetSize(NumOfElements);
 
-      for (j = 0; j < marked_el.Size(); j++)
+      for (int j = 0; j < marked_el.Size(); j++)
       {
-         i = marked_el[j];
+         int i = marked_el[j];
          Segment *c_seg = (Segment *)elements[i];
          int *vert = c_seg->GetVertices(), attr = c_seg->GetAttribute();
          int new_v = cnv + j, new_e = cne + j;
@@ -2962,6 +2959,226 @@ void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
    }
 }
 
+void ParMesh::RefineGroups(const HashTable<Hashed2> &v_to_v)
+{
+   int i, attr, ind, *v;
+
+   int group;
+   Array<int> group_verts, group_edges, group_faces;
+
+   // To update the groups after a refinement, we observe that:
+   // - every (new and old) vertex, edge and face belongs to exactly one group
+   // - the refinement does not create new groups
+   // - a new vertex appears only as the middle of a refined edge
+   // - a face can be refined multiple times producing new edges and faces
+
+   Array<Segment *> sedge_stack;
+   Array<Triangle *> sface_stack;
+
+   Array<int> I_group_svert, J_group_svert;
+   Array<int> I_group_sedge, J_group_sedge;
+   Array<int> I_group_sface, J_group_sface;
+
+   I_group_svert.SetSize(GetNGroups()+1);
+   I_group_sedge.SetSize(GetNGroups()+1);
+   if (Dim == 3)
+   {
+      I_group_sface.SetSize(GetNGroups()+1);
+   }
+
+   I_group_svert[0] = I_group_svert[1] = 0;
+   I_group_sedge[0] = I_group_sedge[1] = 0;
+   if (Dim == 3)
+   {
+      I_group_sface[0] = I_group_sface[1] = 0;
+   }
+
+   for (group = 0; group < GetNGroups()-1; group++)
+   {
+      // Get the group shared objects
+      group_svert.GetRow(group, group_verts);
+      group_sedge.GetRow(group, group_edges);
+      group_sface.GetRow(group, group_faces);
+
+      // Check which edges have been refined
+      for (i = 0; i < group_sedge.RowSize(group); i++)
+      {
+         v = shared_edges[group_edges[i]]->GetVertices();
+         ind = v_to_v.FindId(v[0], v[1]);
+         if (ind == -1) { continue; }
+
+         // This shared edge is refined: walk the whole refinement tree
+         attr = shared_edges[group_edges[i]]->GetAttribute();
+         do
+         {
+            ind += NumOfVertices;
+            // Add new shared vertex
+            group_verts.Append(svert_lvert.Append(ind)-1);
+            // Put the right sub-edge on top of the stack
+            sedge_stack.Append(new Segment(ind, v[1], attr));
+            // The left sub-edge replaces the original edge
+            v[1] = ind;
+            ind = v_to_v.FindId(v[0], ind);
+         }
+         while (ind != -1);
+         // Process all edges in the edge stack
+         do
+         {
+            Segment *se = sedge_stack.Last();
+            v = se->GetVertices();
+            ind = v_to_v.FindId(v[0], v[1]);
+            if (ind == -1)
+            {
+               // The edge 'se' is not refined
+               sedge_stack.DeleteLast();
+               // Add new shared edge
+               shared_edges.Append(se);
+               group_edges.Append(sedge_ledge.Append(-1)-1);
+            }
+            else
+            {
+               // The edge 'se' is refined
+               ind += NumOfVertices;
+               // Add new shared vertex
+               group_verts.Append(svert_lvert.Append(ind)-1);
+               // Put the left sub-edge on top of the stack
+               sedge_stack.Append(new Segment(v[0], ind, attr));
+               // The right sub-edge replaces the original edge
+               v[0] = ind;
+            }
+         }
+         while (sedge_stack.Size() > 0);
+      }
+
+      // Check which faces have been refined
+      for (i = 0; i < group_sface.RowSize(group); i++)
+      {
+         v = shared_faces[group_faces[i]]->GetVertices();
+         ind = v_to_v.FindId(v[0], v[1]);
+         if (ind == -1) { continue; }
+
+         // This shared face is refined: walk the whole refinement tree
+         attr = shared_faces[group_faces[i]]->GetAttribute();
+         const int edge_attr = 1;
+         do
+         {
+            ind += NumOfVertices;
+            // Add the refinement edge to the edge stack
+            sedge_stack.Append(new Segment(v[2], ind, edge_attr));
+            // Put the right sub-triangle on top of the face stack
+            sface_stack.Append(new Triangle(v[1], v[2], ind, attr));
+            // The left sub-triangle replaces the original one
+            v[1] = v[0]; v[0] = v[2]; v[2] = ind;
+            ind = v_to_v.FindId(v[0], v[1]);
+         }
+         while (ind != -1);
+         // Process all faces (triangles) in the face stack
+         do
+         {
+            Triangle *st = sface_stack.Last();
+            v = st->GetVertices();
+            ind = v_to_v.FindId(v[0], v[1]);
+            if (ind == -1)
+            {
+               // The triangle 'st' is not refined
+               sface_stack.DeleteLast();
+               // Add new shared face
+               shared_faces.Append(st);
+               group_faces.Append(sface_lface.Append(-1)-1);
+            }
+            else
+            {
+               // The triangle 'st' is refined
+               ind += NumOfVertices;
+               // Add the refinement edge to the edge stack
+               sedge_stack.Append(new Segment(v[2], ind, edge_attr));
+               // Put the left sub-triangle on top of the face stack
+               sface_stack.Append(new Triangle(v[2], v[0], ind, attr));
+               // The right sub-triangle replaces the original one
+               v[0] = v[1]; v[1] = v[2]; v[2] = ind;
+            }
+         }
+         while (sface_stack.Size() > 0);
+         // Process all edges in the edge stack (same code as above)
+         do
+         {
+            Segment *se = sedge_stack.Last();
+            v = se->GetVertices();
+            ind = v_to_v.FindId(v[0], v[1]);
+            if (ind == -1)
+            {
+               // The edge 'se' is not refined
+               sedge_stack.DeleteLast();
+               // Add new shared edge
+               shared_edges.Append(se);
+               group_edges.Append(sedge_ledge.Append(-1)-1);
+            }
+            else
+            {
+               // The edge 'se' is refined
+               ind += NumOfVertices;
+               // Add new shared vertex
+               group_verts.Append(svert_lvert.Append(ind)-1);
+               // Put the left sub-edge on top of the stack
+               sedge_stack.Append(new Segment(v[0], ind, attr));
+               // The right sub-edge replaces the original edge
+               v[0] = ind;
+            }
+         }
+         while (sedge_stack.Size() > 0);
+      }
+
+      I_group_svert[group+1] = I_group_svert[group] + group_verts.Size();
+      I_group_sedge[group+1] = I_group_sedge[group] + group_edges.Size();
+      if (Dim == 3)
+      {
+         I_group_sface[group+1] = I_group_sface[group] + group_faces.Size();
+      }
+
+      J_group_svert.Append(group_verts);
+      J_group_sedge.Append(group_edges);
+      if (Dim == 3)
+      {
+         J_group_sface.Append(group_faces);
+      }
+   }
+
+   // Fix the local numbers of shared edges and faces: sedge_ledge, sface_lface
+   {
+      NumOfVertices = vertices.Size();
+      DSTable new_v_to_v(NumOfVertices);
+      GetVertexToVertexTable(new_v_to_v);
+      for (i = 0; i < shared_edges.Size(); i++)
+      {
+         v = shared_edges[i]->GetVertices();
+         sedge_ledge[i] = new_v_to_v(v[0], v[1]);
+      }
+   }
+   if (Dim == 3)
+   {
+      STable3D *faces_tbl = GetElementToFaceTable(1);
+      for (i = 0; i < shared_faces.Size(); i++)
+      {
+         v = shared_faces[i]->GetVertices();
+         sface_lface[i] = (*faces_tbl)(v[0], v[1], v[2]);
+      }
+      delete faces_tbl;
+   }
+
+   group_svert.SetIJ(I_group_svert, J_group_svert);
+   group_sedge.SetIJ(I_group_sedge, J_group_sedge);
+   if (Dim == 3)
+   {
+      group_sface.SetIJ(I_group_sface, J_group_sface);
+   }
+   I_group_svert.LoseData(); J_group_svert.LoseData();
+   I_group_sedge.LoseData(); J_group_sedge.LoseData();
+   if (Dim == 3)
+   {
+      I_group_sface.LoseData(); J_group_sface.LoseData();
+   }
+}
+
 void ParMesh::QuadUniformRefinement()
 {
    DeleteFaceNbrData();

mesh/pmesh.hpp

@@ -63,8 +63,12 @@ protected:
 
    /// Return a number(0-1) identifying how the given edge has been split
    int GetEdgeSplittings(Element *edge, const DSTable &v_to_v, int *middle);
-   /// Return a number(0-4) identifying how the given face has been split
-   int GetFaceSplittings(Element *face, const DSTable &v_to_v, int *middle);
+   /// Append codes identifying how the given face has been split to @a codes
+   void GetFaceSplittings(Element *face, const HashTable<Hashed2> &v_to_v,
+                          Array<unsigned> &codes);
+
+   bool DecodeFaceSplittings(HashTable<Hashed2> &v_to_v, const int *v,
+                             const Array<unsigned> &codes, int &pos);
 
    void GetFaceNbrElementTransformation(
       int i, IsoparametricTransformation *ElTr);
@@ -185,8 +189,10 @@ public:
    /// Utility function: sum integers from all processors (Allreduce).
    virtual long ReduceInt(int value) const;
 
-   /// Update the groups after tet refinement
+   /// Update the groups after triangle refinement
    void RefineGroups(const DSTable &v_to_v, int *middle);
+   /// Update the groups after tetrahedron refinement
+   void RefineGroups(const HashTable<Hashed2> &v_to_v);
 
    /// Load balance the mesh. NC meshes only.
    void Rebalance();

mesh/tetrahedron.cpp

@@ -153,22 +153,14 @@ void Tetrahedron::GetMarkedFace(const int face, int *fv)
    }
 }
 
-int Tetrahedron::NeedRefinement(DSTable &v_to_v, int *middle) const
+int Tetrahedron::NeedRefinement(HashTable<Hashed2> &v_to_v) const
 {
-   int m;
-
-   if ((m = v_to_v(indices[0], indices[1])) != -1)
-      if (middle[m] != -1) { return 1; }
-   if ((m = v_to_v(indices[1], indices[2])) != -1)
-      if (middle[m] != -1) { return 1; }
-   if ((m = v_to_v(indices[2], indices[0])) != -1)
-      if (middle[m] != -1) { return 1; }
-   if ((m = v_to_v(indices[0], indices[3])) != -1)
-      if (middle[m] != -1) { return 1; }
-   if ((m = v_to_v(indices[1], indices[3])) != -1)
-      if (middle[m] != -1) { return 1; }
-   if ((m = v_to_v(indices[2], indices[3])) != -1)
-      if (middle[m] != -1) { return 1; }
+   if (v_to_v.FindId(indices[0], indices[1]) != -1) { return 1; }
+   if (v_to_v.FindId(indices[1], indices[2]) != -1) { return 1; }
+   if (v_to_v.FindId(indices[2], indices[0]) != -1) { return 1; }
+   if (v_to_v.FindId(indices[0], indices[3]) != -1) { return 1; }
+   if (v_to_v.FindId(indices[1], indices[3]) != -1) { return 1; }
+   if (v_to_v.FindId(indices[2], indices[3]) != -1) { return 1; }
    return 0;
 }
 

mesh/tetrahedron.hpp

@@ -64,7 +64,7 @@ public:
    void SetRefinementFlag(int rf) { refinement_flag = rf; }
 
    /// Return 1 if the element needs refinement in order to get conforming mesh.
-   virtual int NeedRefinement(DSTable &v_to_v, int *middle) const;
+   virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const;
 
    /// Set the vertices according to the given input.
    virtual void SetVertices(const int *ind);

mesh/triangle.cpp

@@ -34,13 +34,11 @@ Triangle::Triangle(int ind1, int ind2, int ind3, int attr)
    transform = 0;
 }
 
-int Triangle::NeedRefinement(DSTable &v_to_v, int *middle) const
+int Triangle::NeedRefinement(HashTable<Hashed2> &v_to_v) const
 {
-   int m;
-
-   if ((m = v_to_v(indices[0], indices[1])) != -1 && middle[m] != -1) { return 1; }
-   if ((m = v_to_v(indices[1], indices[2])) != -1 && middle[m] != -1) { return 1; }
-   if ((m = v_to_v(indices[2], indices[0])) != -1 && middle[m] != -1) { return 1; }
+   if (v_to_v.FindId(indices[0], indices[1]) != -1) { return 1; }
+   if (v_to_v.FindId(indices[1], indices[2]) != -1) { return 1; }
+   if (v_to_v.FindId(indices[2], indices[0]) != -1) { return 1; }
    return 0;
 }
 

mesh/triangle.hpp

@@ -42,7 +42,7 @@ public:
    virtual int GetType() const { return Element::TRIANGLE; }
 
    /// Return 1 if the element needs refinement in order to get conforming mesh.
-   virtual int NeedRefinement(DSTable &v_to_v, int *middle) const;
+   virtual int NeedRefinement(HashTable<Hashed2> &v_to_v) const;
 
    /// Set the vertices according to the given input.
    virtual void SetVertices(const int *ind);