Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18170) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18171) @@ -2234,4 +2234,27 @@ penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1; break; + case P2bubbleEnum: case P2bubblecondensedEnum: + numnodes = 19; + penta_node_ids = xNew(numnodes); + penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; + penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1]; + penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2]; + penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3]; + penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4]; + penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5]; + penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; + penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; + penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; + penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; + penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; + penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; + penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; + penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; + penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; + penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1; + penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1; + penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1; + penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1; + break; case P1P1Enum: case P1P1GLSEnum: numnodes = 12; Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 18170) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 18171) @@ -24,4 +24,5 @@ #define NUMNODESP2xP1 12 #define NUMNODESP2 18 +#define NUMNODESP2b 19 #define NUMNODESP2xP4 30 @@ -237,4 +238,30 @@ basis[16]=4.*gauss->coord3*gauss->coord1*(1.-zeta*zeta); basis[17]=4.*gauss->coord1*gauss->coord2*(1.-zeta*zeta); + return; + case P2bubbleEnum: case P2bubblecondensedEnum: + /*Corner nodes*/ + basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta-1.)/2.; + basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta-1.)/2.; + basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta-1.)/2.; + basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta+1.)/2.; + basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta+1.)/2.; + basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.; + /*mid-sides of quads*/ + basis[ 6]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta*zeta); + basis[ 7]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta*zeta); + basis[ 8]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta*zeta); + /*mid-sides of triangles*/ + basis[ 9]=4.*gauss->coord3*gauss->coord2*zeta*(zeta-1.)/2.; + basis[10]=4.*gauss->coord3*gauss->coord1*zeta*(zeta-1.)/2.; + basis[11]=4.*gauss->coord1*gauss->coord2*zeta*(zeta-1.)/2.; + basis[12]=4.*gauss->coord3*gauss->coord2*zeta*(zeta+1.)/2.; + basis[13]=4.*gauss->coord3*gauss->coord1*zeta*(zeta+1.)/2.; + basis[14]=4.*gauss->coord1*gauss->coord2*zeta*(zeta+1.)/2.; + /*quad faces*/ + basis[15]=4.*gauss->coord3*gauss->coord2*(1.-zeta*zeta); + basis[16]=4.*gauss->coord3*gauss->coord1*(1.-zeta*zeta); + basis[17]=4.*gauss->coord1*gauss->coord2*(1.-zeta*zeta); + /*bubble*/ + basis[18]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta); return; case P2xP4Enum : @@ -572,4 +599,86 @@ dbasis[NUMNODESP2*1+17] = (1.-zeta*zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2)); dbasis[NUMNODESP2*2+17] = -2.*zeta*4.*gauss->coord1*gauss->coord2; + return; + case P2bubbleEnum: case P2bubblecondensedEnum: + /*Nodal function 1*/ + dbasis[NUMNODESP2b*0+0 ] = .5*zeta*(zeta-1.)*(-2.*gauss->coord1 + 0.5); + dbasis[NUMNODESP2b*1+0 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.); + dbasis[NUMNODESP2b*2+0 ] = .5*(2.*zeta-1.)*gauss->coord1*(2.*gauss->coord1-1.); + /*Nodal function 2*/ + dbasis[NUMNODESP2b*0+1 ] = .5*zeta*(zeta-1.)*(+2.*gauss->coord2 - 0.5); + dbasis[NUMNODESP2b*1+1 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.); + dbasis[NUMNODESP2b*2+1 ] = .5*(2.*zeta-1.)*gauss->coord2*(2.*gauss->coord2-1.); + /*Nodal function 3*/ + dbasis[NUMNODESP2b*0+2 ] = 0.; + dbasis[NUMNODESP2b*1+2 ] = .5*zeta*(zeta-1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.); + dbasis[NUMNODESP2b*2+2 ] = .5*(2.*zeta-1.)*gauss->coord3*(2.*gauss->coord3-1.); + /*Nodal function 4*/ + dbasis[NUMNODESP2b*0+3 ] = .5*zeta*(zeta+1.)*(-2.*gauss->coord1 + 0.5); + dbasis[NUMNODESP2b*1+3 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.); + dbasis[NUMNODESP2b*2+3 ] = .5*(2.*zeta+1.)*gauss->coord1*(2.*gauss->coord1-1.); + /*Nodal function 5*/ + dbasis[NUMNODESP2b*0+4 ] = .5*zeta*(zeta+1.)*(+2.*gauss->coord2 - 0.5); + dbasis[NUMNODESP2b*1+4 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.); + dbasis[NUMNODESP2b*2+4 ] = .5*(2.*zeta+1.)*gauss->coord2*(2.*gauss->coord2-1.); + /*Nodal function 6*/ + dbasis[NUMNODESP2b*0+5 ] = 0.; + dbasis[NUMNODESP2b*1+5 ] = .5*zeta*(zeta+1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.); + dbasis[NUMNODESP2b*2+5 ] = .5*(2.*zeta+1.)*gauss->coord3*(2.*gauss->coord3-1.); + + /*Nodal function 7*/ + dbasis[NUMNODESP2b*0+6 ] = (-2.*gauss->coord1 + 0.5)*(1.-zeta*zeta); + dbasis[NUMNODESP2b*1+6 ] = (-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.)*(1.-zeta*zeta); + dbasis[NUMNODESP2b*2+6 ] = -2.*zeta*gauss->coord1*(2.*gauss->coord1-1.); + /*Nodal function 8*/ + dbasis[NUMNODESP2b*0+7 ] = (+2.*gauss->coord2 - 0.5)*(1.-zeta*zeta); + dbasis[NUMNODESP2b*1+7 ] = (-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.)*(1.-zeta*zeta); + dbasis[NUMNODESP2b*2+7 ] = -2.*zeta*gauss->coord2*(2.*gauss->coord2-1.); + /*Nodal function 9*/ + dbasis[NUMNODESP2b*0+8 ] = 0.; + dbasis[NUMNODESP2b*1+8 ] = (+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.)*(1.-zeta*zeta); + dbasis[NUMNODESP2b*2+8 ] = -2.*zeta*gauss->coord3*(2.*gauss->coord3-1.); + + /*Nodal function 10*/ + dbasis[NUMNODESP2b*0+9 ] = zeta*(zeta-1.)*gauss->coord3; + dbasis[NUMNODESP2b*1+9 ] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2b*2+9 ] = 2.*gauss->coord3*gauss->coord2*(2.*zeta-1.); + /*Nodal function 11*/ + dbasis[NUMNODESP2b*0+10] = -zeta*(zeta-1.)*gauss->coord3; + dbasis[NUMNODESP2b*1+10] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2b*2+10] = 2.*gauss->coord3*gauss->coord1*(2.*zeta-1.); + /*Nodal function 12*/ + dbasis[NUMNODESP2b*0+11] = zeta*(zeta-1.)*(gauss->coord1-gauss->coord2); + dbasis[NUMNODESP2b*1+11] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2)); + dbasis[NUMNODESP2b*2+11] = 2.*gauss->coord1*gauss->coord2*(2.*zeta-1.); + /*Nodal function 13*/ + dbasis[NUMNODESP2b*0+12] = zeta*(zeta+1.)*gauss->coord3; + dbasis[NUMNODESP2b*1+12] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2b*2+12] = 2.*gauss->coord3*gauss->coord2*(2.*zeta+1.); + /*Nodal function 14*/ + dbasis[NUMNODESP2b*0+13] = -zeta*(zeta+1.)*gauss->coord3; + dbasis[NUMNODESP2b*1+13] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2b*2+13] = 2.*gauss->coord3*gauss->coord1*(2.*zeta+1.); + /*Nodal function 15*/ + dbasis[NUMNODESP2b*0+14] = zeta*(zeta+1.)*(gauss->coord1-gauss->coord2); + dbasis[NUMNODESP2b*1+14] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2)); + dbasis[NUMNODESP2b*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.); + + /*Nodal function 16*/ + dbasis[NUMNODESP2b*0+15] = 2.*gauss->coord3*(1.-zeta*zeta); + dbasis[NUMNODESP2b*1+15] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2b*2+15] = -2.*zeta*4.*gauss->coord3*gauss->coord2; + /*Nodal function 17*/ + dbasis[NUMNODESP2b*0+16] = -2.*gauss->coord3*(1.-zeta*zeta); + dbasis[NUMNODESP2b*1+16] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2b*2+16] = -2.*zeta*4.*gauss->coord3*gauss->coord1; + /*Nodal function 18*/ + dbasis[NUMNODESP2b*0+17] = 2.*(gauss->coord1-gauss->coord2)*(1.-zeta*zeta); + dbasis[NUMNODESP2b*1+17] = (1.-zeta*zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2)); + dbasis[NUMNODESP2b*2+17] = -2.*zeta*4.*gauss->coord1*gauss->coord2; + + /*Nodal function 19*/ + dbasis[NUMNODESP2b*0+18] = 27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3); + dbasis[NUMNODESP2b*1+18] = 27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2); + dbasis[NUMNODESP2b*2+18] = -54*gauss->coord1*gauss->coord2*gauss->coord3*zeta; return; case P2xP4Enum : @@ -841,4 +950,6 @@ case P1bubblecondensedEnum: return NUMNODESP1b; case P2Enum: return NUMNODESP2; + case P2bubbleEnum: return NUMNODESP2b; + case P2bubblecondensedEnum: return NUMNODESP2b; case P2xP1Enum: return NUMNODESP2xP1; case P1xP2Enum: return NUMNODESP1xP2; Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 18170) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 18171) @@ -647,4 +647,32 @@ } break; + case P2bubbleEnum: case P2bubblecondensedEnum: + switch(iv){ + case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break; + case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break; + case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break; + case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break; + case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break; + case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break; + + case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break; + case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break; + case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break; + + case 9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break; + case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break; + case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break; + case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break; + case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break; + case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break; + + case 15: coord1=0.; coord2=.5; coord3=.5; coord4=0.;break; + case 16: coord1=.5; coord2=0.; coord3=.5; coord4=0.;break; + case 17: coord1=.5; coord2=.5; coord3=0.; coord4=0.;break; + + case 18: coord1=1./3.; coord2=1./3.; coord3=1./3.; coord4=0.; break; + default: _error_("node index should be in [0 18]"); + } + break; case P2xP4Enum: switch(iv){ Index: /issm/trunk-jpl/src/m/classes/flowequation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18170) +++ /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18171) @@ -133,5 +133,5 @@ md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'}); - md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P1xP3','P2xP4'}); + md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'}); md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'}); md = checkfield(md,'fieldname','flowequation.XTH_r','numel',[1],'>',0.); Index: /issm/trunk-jpl/src/m/classes/flowequation.py =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.py (revision 18170) +++ /issm/trunk-jpl/src/m/classes/flowequation.py (revision 18171) @@ -77,5 +77,5 @@ md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0,1]) md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2','P2bubble']) - md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P1xP3','P2xP4']) + md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4']) md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z']) md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1])