Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 17841) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 17842) @@ -37,4 +37,5 @@ this->numberoffaces=-1; this->numberofedges=-1; + this->facescols=-1; this->elements=NULL; this->faces=NULL; @@ -83,4 +84,5 @@ FetchData(&this->numberofelements,MeshNumberofelementsEnum); FetchData(&this->elements,NULL,NULL,MeshElementsEnum); + this->facescols = -1; this->faces = NULL; this->edges = NULL; Index: /issm/trunk-jpl/src/c/classes/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.h (revision 17841) +++ /issm/trunk-jpl/src/c/classes/IoModel.h (revision 17842) @@ -38,4 +38,5 @@ int numberoffaces; int numberofedges; + int facescols; int *elements; int *faces; Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 17841) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 17842) @@ -35,5 +35,4 @@ bool isnan; int i,j,count,elementnbv,numfacevertices; - int* faceverticesid = NULL; IssmDouble value; IssmDouble *times = NULL; @@ -119,21 +118,18 @@ if(iomodel->meshelementtype==PentaEnum){ for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ if(my_faces[i]){ - FaceGetVertexIndices(iomodel,&numfacevertices,&faceverticesid,i); - isnan=0; + numfacevertices = iomodel->faces[i*iomodel->facescols+3]; + value=0.; for(j=0;j(spcdata[faceverticesid[j]-1])) isnan=1; + value += spcdata[iomodel->faces[i*iomodel->facescols+4+j] -1]; } - if(isnan==0){ - value=0; - for(j=0;jAddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1, + value = value/reCast(numfacevertices); + if(!xIsNan(value)){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1, + iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1, dof,value,analysis_type)); count++; } - xDelete(faceverticesid); } } @@ -179,25 +175,24 @@ } for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ if(my_faces[i]){ - FaceGetVertexIndices(iomodel,&numfacevertices,&faceverticesid,i); - isnan=0; + numfacevertices = iomodel->faces[i*iomodel->facescols+3]; + value=0.; for(j=0;j(spcdata[faceverticesid[j]-1])) isnan=1; - } - if(isnan==0){ - value=0; - for(j=0;jAddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1, + value += spcdata[iomodel->faces[i*iomodel->facescols+4+j] -1]; + } + value = value/reCast(numfacevertices); + if(!xIsNan(value)){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1, + iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1, dof,value,analysis_type)); - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2, + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2, + iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2, dof,value,analysis_type)); - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3, + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3, + iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3, dof,value,analysis_type)); count=count+3; } - xDelete(faceverticesid); } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp (revision 17841) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp (revision 17842) @@ -5,6 +5,7 @@ #include "../../classes/classes.h" #include "../../shared/shared.h" +#include "./ModelProcessorx.h" -void CreateEdges(IoModel* iomodel){ +void CreateEdges(IoModel* iomodel){/*{{{*/ /*If edges are already present, exit*/ @@ -142,3 +143,24 @@ iomodel->elementtoedgeconnectivity = element_edge_connectivity; iomodel->numberofedges = nbe; +}/*}}}*/ +void EdgeOnBoundaryFlags(IoModel* iomodel,bool** pflags){ + + /*Get faces and edges*/ + if(!iomodel->edges) CreateEdges(iomodel); + if(!iomodel->faces) CreateFaces(iomodel); + + /*Allocate output*/ + bool* flags = xNew(iomodel->numberofedges); + + /*Now, loop over the faces, whenever it is not connected to a second element, all edges are on boundary*/ + for(int i=0;inumberoffaces;i++){ + + if(iomodel->faces[i*iomodel->facescols+4]==-1){ + + /*Get all edges for this face*/ + } + } + + /*Clean up and return*/ + *pflags = flags; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp (revision 17841) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp (revision 17842) @@ -126,15 +126,45 @@ /*Intermediaries*/ bool exist; - int i,j,k,v0,cols; - int maxnbf,nbf,elementnbf,elementnbv,facemaxnbv; + int i,j,k,v0,cols,facescols; + int maxnbf,nbf,elementnbf,elementnbv,facenbv,facemaxnbv; int *elementfaces = NULL; int *elementfaces_markers = NULL; - /*Maximum number of faces*/ - maxnbf = 5*iomodel->numberofelements; + /*Mesh specific face indexing per element*/ + switch(iomodel->meshelementtype){ + case PentaEnum: + 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; + break; + default: + _error_("mesh "<< EnumToStringx(iomodel->meshelementtype) <<" not supported"); + } + + /*Allocate connectivity*/ + int *element_face_connectivity = xNew(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */ + + /*Maximum number of faces for initial allocation*/ + maxnbf = elementnbf*iomodel->numberofelements; + facescols = 4+facemaxnbv; _assert_(facescols>6); /*Initialize intermediaries*/ - int* facestemp = xNew(maxnbf*6); /*format: [vertex1 vertex2 vertex3 element1 element2 marker] */ - for(i=0;i(maxnbf*facescols); /*format: [element1 element2 marker nbv vertex1 vertex2 vertex3 ...] */ + for(i=0;inumberofvertices;i++) head_minv[i]=-1; - /*Mesh specific face indexing per element*/ - if(iomodel->meshelementtype==PentaEnum){ - 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; - } - else{ - _error_("mesh dimension not supported yet"); - } - int *element_face_connectivity = xNew(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */ - /*Initialize number of faces and list of vertex indices*/ nbf = 0; @@ -177,5 +181,6 @@ /*Get indices of current face*/ - for(k=0;kelements[i*elementnbv + elementfaces[cols*j+k+1]] - 1; } @@ -183,16 +188,15 @@ /*Sort list of vertices*/ HeapSort(v,elementfaces[cols*j+0]); - v0 = v[0]; + v0 = v[0]; _assert_(v0>=0 && v0numberofvertices); /*This face a priori has not been processed yet*/ exist = false; - /*Go through all processed faces connected to v1 and check whether we have seen this edge yet*/ - _assert_(v0>=0 && v0numberofvertices); - for(int e=head_minv[v0]; e!=-1; e=next_face[e]){ - if(facestemp[e*6+1]==v[1]+1 && facestemp[e*6+2]==v[2]+1){ + /*Go through all processed faces connected to v0 and check whether we have seen this face yet*/ + for(int f=head_minv[v0]; f!=-1; f=next_face[f]){ + if(facestemp[f*facescols+5]==v[1]+1 && facestemp[f*facescols+6]==v[2]+1){ exist = true; - facestemp[e*6+4]=i+1; - element_face_connectivity[i*elementnbf+j]=e; + facestemp[f*facescols+1]=i+1; + element_face_connectivity[i*elementnbf+j]=f; break; } @@ -204,9 +208,8 @@ /*Update faces*/ - facestemp[nbf*6+0] = v[0]+1; - facestemp[nbf*6+1] = v[1]+1; - facestemp[nbf*6+2] = v[2]+1; - facestemp[nbf*6+3] = i+1; - facestemp[nbf*6+5] = elementfaces_markers[j]; + facestemp[nbf*facescols+0] = i+1; + facestemp[nbf*facescols+2] = elementfaces_markers[j]; + facestemp[nbf*facescols+3] = facenbv; + for(k=0;k(nbf*6); /*vertex1 vertex2 vertex3 element1 element2 marker*/ - xMemCpy(faces,facestemp,nbf*6); + int* faces = xNew(nbf*facescols); + xMemCpy(faces,facestemp,nbf*facescols); xDelete(facestemp); @@ -238,74 +241,5 @@ iomodel->faces = faces; iomodel->numberoffaces = nbf; + iomodel->facescols = facescols; iomodel->elementtofaceconnectivity = element_face_connectivity; }/*}}}*/ -void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber){/*{{{*/ - - int numbervertices; - if(iomodel->domaintype==Domain3DEnum){ - if((iomodel->faces[6*facenumber+5])==1){ - numbervertices=3; - } - else if((iomodel->faces[6*facenumber+5])==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,k,cols,faceid; - int maxnbf,nbf,elementnbf,elementnbv,facemaxnbv; - int *elementfaces = NULL; - int *elementfaces_markers = NULL; - int elementid=iomodel->faces[6*facenumber+3]; - 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++){ - for(k=1;k<5;k++){ - if(iomodel->elements[(elementid-1)*6+elementfaces[cols*faceid+k]] == iomodel->faces[6*facenumber+j]) counter++; - } - } - if(counter==3) break; - } - if(counter!=3) _error_("face not found in element"); - - for(j=0;j<4;j++){ - facevertices[j]=iomodel->elements[(elementid-1)*6+elementfaces[cols*faceid+j+1]]; - } - - /*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 17841) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 17842) @@ -147,5 +147,5 @@ FacesPartitioning(&my_faces,iomodel); for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ if(my_faces[i]){ node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,approximation); @@ -153,5 +153,5 @@ } } - else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/ + else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ /*Nothing*/ } @@ -196,5 +196,5 @@ id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges; for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ if(my_faces[i]){ node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation); @@ -207,5 +207,5 @@ counter=counter+3; } - else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/ + else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ /*Nothing*/ } @@ -312,5 +312,5 @@ FacesPartitioning(&my_faces,iomodel); for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ if(my_faces[i]){ node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum); @@ -318,5 +318,5 @@ } } - else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/ + else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ /*Nothing*/ } @@ -377,5 +377,5 @@ id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges; for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*6+5]==2){/*Vertical quads*/ + if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ if(my_faces[i]){ node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,FSvelocityEnum); @@ -388,5 +388,5 @@ counter=counter+3; } - else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/ + else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ /*Nothing*/ } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 17841) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 17842) @@ -33,5 +33,5 @@ void CreateFaces(IoModel* iomodel); void CreateFaces3d(IoModel* iomodel); -void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber); +void EdgeOnBoundaryFlags(IoModel* iomodel,bool** pflags); /*Connectivity*/