Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17382) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17383) @@ -2617,5 +2617,5 @@ break; case P2Enum: - numnodes = 15; + numnodes = 18; penta_node_ids = xNew(numnodes); penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; @@ -2634,4 +2634,7 @@ 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; break; case P1P1Enum: case P1P1GLSEnum: @@ -2671,5 +2674,5 @@ break; case TaylorHoodEnum: - numnodes = 21; + numnodes = 24; penta_node_ids = xNew(numnodes); penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; @@ -2688,11 +2691,14 @@ 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->elements[6*index+0]; - penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+1]; - penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+2]; - penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+3]; - penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+4]; - penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+5]; + 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+iomodel->elements[6*index+0]; + penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+1]; + penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+2]; + penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+3]; + penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+4]; + penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+5]; break; case OneLayerP4zEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 17382) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 17383) @@ -23,5 +23,5 @@ #define NUMNODESP1xP3 12 #define NUMNODESP2xP1 12 -#define NUMNODESP2 15 +#define NUMNODESP2 18 #define NUMNODESP2xP4 30 @@ -259,14 +259,18 @@ basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.; /*mid-sides of quads*/ - basis[ 6]=gauss->coord1*(1.-zeta*zeta); - basis[ 7]=gauss->coord2*(1.-zeta*zeta); - basis[ 8]=gauss->coord3*(1.-zeta*zeta); + 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]=2.*gauss->coord3*gauss->coord2*zeta*(zeta-1.); - basis[10]=2.*gauss->coord3*gauss->coord1*zeta*(zeta-1.); - basis[11]=2.*gauss->coord1*gauss->coord2*zeta*(zeta-1.); - basis[12]=2.*gauss->coord3*gauss->coord2*zeta*(zeta+1.); - basis[13]=2.*gauss->coord3*gauss->coord1*zeta*(zeta+1.); - basis[14]=2.*gauss->coord1*gauss->coord2*zeta*(zeta+1.); + 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); return; case P2xP4Enum : @@ -599,15 +603,15 @@ /*Nodal function 7*/ - dbasis[NUMNODESP2*0+6 ] = -0.5*(1.-zeta*zeta); - dbasis[NUMNODESP2*1+6 ] = -SQRT3/6.*(1.-zeta*zeta); - dbasis[NUMNODESP2*2+6 ] = -2.*zeta*gauss->coord1; + dbasis[NUMNODESP2*0+6 ] = (-2.*gauss->coord1 + 0.5)*(1.-zeta*zeta); + dbasis[NUMNODESP2*1+6 ] = (-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.)*(1.-zeta*zeta); + dbasis[NUMNODESP2*2+6 ] = -2.*zeta*gauss->coord1*(2.*gauss->coord1-1.); /*Nodal function 8*/ - dbasis[NUMNODESP2*0+7 ] = 0.5*(1.-zeta*zeta); - dbasis[NUMNODESP2*1+7 ] = -SQRT3/6.*(1.-zeta*zeta); - dbasis[NUMNODESP2*2+7 ] = -2.*zeta*gauss->coord2; + dbasis[NUMNODESP2*0+7 ] = (+2.*gauss->coord2 - 0.5)*(1.-zeta*zeta); + dbasis[NUMNODESP2*1+7 ] = (-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.)*(1.-zeta*zeta); + dbasis[NUMNODESP2*2+7 ] = -2.*zeta*gauss->coord2*(2.*gauss->coord2-1.); /*Nodal function 9*/ dbasis[NUMNODESP2*0+8 ] = 0.; - dbasis[NUMNODESP2*1+8 ] = SQRT3/3.*(1.-zeta*zeta); - dbasis[NUMNODESP2*2+8 ] = -2.*zeta*gauss->coord3; + dbasis[NUMNODESP2*1+8 ] = (+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.)*(1.-zeta*zeta); + dbasis[NUMNODESP2*2+8 ] = -2.*zeta*gauss->coord3*(2.*gauss->coord3-1.); /*Nodal function 10*/ @@ -635,4 +639,17 @@ dbasis[NUMNODESP2*1+14] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2)); dbasis[NUMNODESP2*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.); + + /*Nodal function 16*/ + dbasis[NUMNODESP2*0+15] = 2.*gauss->coord3*(1.-zeta*zeta); + dbasis[NUMNODESP2*1+15] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2*2+15] = -2.*zeta*4.*gauss->coord3*gauss->coord2; + /*Nodal function 17*/ + dbasis[NUMNODESP2*0+16] = -2.*gauss->coord3*(1.-zeta*zeta); + dbasis[NUMNODESP2*1+16] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3); + dbasis[NUMNODESP2*2+16] = -2.*zeta*4.*gauss->coord3*gauss->coord1; + /*Nodal function 18*/ + dbasis[NUMNODESP2*0+17] = 2.*(gauss->coord1-gauss->coord2)*(1.-zeta*zeta); + 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 P2xP4Enum : Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 17382) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 17383) @@ -654,5 +654,9 @@ case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break; case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break; - default: _error_("node index should be in [0 14]"); + + 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; + default: _error_("node index should be in [0 17]"); } break; Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 17382) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 17383) @@ -33,5 +33,7 @@ /*intermediary: */ - int i,j,count,elementnbv; + bool isnan; + int i,j,count,elementnbv,numfacevertices; + int* faceverticesid = NULL; IssmDouble value; IssmDouble *times = NULL; @@ -70,4 +72,5 @@ case P2Enum: EdgesPartitioning(&my_edges,iomodel); + FacesPartitioning(&my_faces,iomodel); break; case P2xP4Enum: @@ -109,4 +112,24 @@ dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); count++; + } + } + } + for(i=0;inumberoffaces;i++){ + if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(my_faces[i]){ + FaceGetVertexIndices(iomodel,&numfacevertices,&faceverticesid,i); + isnan=0; + for(j=0;j(spcdata[faceverticesid[j]])) isnan=1; + } + if(isnan==0){ + value=0; + for(j=0;jAddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1, + dof,value,analysis_type)); + count++; + } } } @@ -431,4 +454,5 @@ /*Free ressources:*/ + xDelete(faceverticesid); xDelete(times); xDelete(values); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp (revision 17382) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp (revision 17383) @@ -240,2 +240,68 @@ iomodel->elementtofaceconnectivity = element_face_connectivity; }/*}}}*/ +void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber){/*{{{*/ + + int numbervertices; + if(iomodel->meshtype==Mesh3DEnum){ + if((iomodel->faces[6*facenumber+5]-1)==1){ + numbervertices=3; + } + else if((iomodel->faces[6*facenumber+5]-1)==2){ + numbervertices=4; + } + else _error_("face marker not supported yet "); + } + else _error_("mesh type not supported yet"); + + int* facevertices = xNew(numbervertices); + if(numbervertices==3){ + for(int i=0;i<3;i++) facevertices[i]=iomodel->faces[6*facenumber+i]-1; + } + else if(numbervertices==4){ + int i,j,cols,faceid; + int maxnbf,nbf,elementnbf,elementnbv,facemaxnbv; + int *elementfaces = NULL; + int *elementfaces_markers = NULL; + int elementid=iomodel->faces[6*facenumber+4]-1; + int counter=0; + + /*Recreate element properties*/ + elementnbv = 6; /*Number of vertices per element*/ + elementnbf = 5; /*Number of faces per element*/ + facemaxnbv = 4; /*Maximum number of vertices per face*/ + cols = facemaxnbv + 1; + elementfaces = xNew(elementnbf*cols); + elementfaces_markers = xNew(elementnbf); + /*2 triangles*/ + elementfaces_markers[0] = 1; + elementfaces_markers[1] = 1; + elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0; elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2; + elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3; elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5; + /*3 quads*/ + elementfaces_markers[2] = 2; + elementfaces_markers[3] = 2; + elementfaces_markers[4] = 2; + elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 1; elementfaces[cols*2+2] = 2; elementfaces[cols*2+3] = 5; elementfaces[cols*2+4] = 4; + elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 2; elementfaces[cols*3+2] = 0; elementfaces[cols*3+3] = 3; elementfaces[cols*3+4] = 5; + elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 0; elementfaces[cols*4+2] = 1; elementfaces[cols*4+3] = 4; elementfaces[cols*4+4] = 3; + + for(faceid=2;faceid<5;faceid++){ + counter=0; + for(j=0;j<3;j++){ + if(iomodel->elements[elementid*6+elementfaces[cols*faceid+j]] == iomodel->faces[6*facenumber+j]) counter++; + } + if(counter==3) break; + } + if(!counter==3) _error_("face not found in element"); + + for(j=0;j<3;j++) facevertices[i]=iomodel->elements[elementid*6+elementfaces[cols*faceid+j]]; + + /*Delete*/ + xDelete(elementfaces); + xDelete(elementfaces_markers); + } + + /*Output results*/ + *pverticesid=facevertices; + *pnumvertices=numbervertices; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 17382) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 17383) @@ -132,4 +132,5 @@ case P2Enum: EdgesPartitioning(&my_edges,iomodel); + FacesPartitioning(&my_faces,iomodel); for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ @@ -140,4 +141,19 @@ if(my_edges[i]){ nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,approximation)); + } + } + id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; + for(i=0;inumberoffaces;i++){ + if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(my_faces[i]){ + node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,approximation); + nodes->AddObject(node); + } + } + else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/ + /*Nothing*/ + } + else{ + _error_("not supported"); } } @@ -278,4 +294,5 @@ /*P2 velocity*/ EdgesPartitioning(&my_edges,iomodel); + FacesPartitioning(&my_faces,iomodel); for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ @@ -288,10 +305,25 @@ } } + id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; + for(i=0;inumberoffaces;i++){ + if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(my_faces[i]){ + node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum); + nodes->AddObject(node); + } + } + else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/ + /*Nothing*/ + } + else{ + _error_("not supported"); + } + } /*P1 pressure*/ - vnodes = id0+iomodel->numberofvertices+iomodel->numberofedges; - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,i,iomodel,analysis,FSpressureEnum)); + vnodes = id0+iomodel->numberoffaces; + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+i,lid++,i,iomodel,analysis,FSpressureEnum)); } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 17382) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 17383) @@ -33,4 +33,5 @@ void CreateFaces(IoModel* iomodel); void CreateFaces3d(IoModel* iomodel); +void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber); /*Connectivity*/