Changeset 17383


Ignore:
Timestamp:
03/05/14 08:18:13 (11 years ago)
Author:
seroussi
Message:

BUG: fixed P2 in Penta

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r17372 r17383  
    26172617                        break;
    26182618                case P2Enum:
    2619                         numnodes         = 15;
     2619                        numnodes         = 18;
    26202620                        penta_node_ids   = xNew<int>(numnodes);
    26212621                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
    … …  
    26342634                        penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
    26352635                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
     2636                        penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
     2637                        penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
     2638                        penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
    26362639                        break;
    26372640                case P1P1Enum: case P1P1GLSEnum:
    … …  
    26712674                        break;
    26722675                case TaylorHoodEnum:
    2673                         numnodes         = 21;
     2676                        numnodes         = 24;
    26742677                        penta_node_ids   = xNew<int>(numnodes);
    26752678                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
    … …  
    26882691                        penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
    26892692                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
    2690 
    2691                         penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+0];
    2692                         penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+1];
    2693                         penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+2];
    2694                         penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+3];
    2695                         penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+4];
    2696                         penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+5];
     2693                        penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
     2694                        penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
     2695                        penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
     2696
     2697                        penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+0];
     2698                        penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+1];
     2699                        penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+2];
     2700                        penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+3];
     2701                        penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+4];
     2702                        penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+5];
    26972703                        break;
    26982704                case OneLayerP4zEnum:

  • issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp

    r17225 r17383  
    2323#define NUMNODESP1xP3 12
    2424#define NUMNODESP2xP1 12
    25 #define NUMNODESP2    15
     25#define NUMNODESP2    18
    2626#define NUMNODESP2xP4 30
    2727
    … …  
    259259                        basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.;
    260260                        /*mid-sides of quads*/
    261                         basis[ 6]=gauss->coord1*(1.-zeta*zeta);
    262                         basis[ 7]=gauss->coord2*(1.-zeta*zeta);
    263                         basis[ 8]=gauss->coord3*(1.-zeta*zeta);
     261                        basis[ 6]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta*zeta);
     262                        basis[ 7]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta*zeta);
     263                        basis[ 8]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta*zeta);
    264264                        /*mid-sides of triangles*/
    265                         basis[ 9]=2.*gauss->coord3*gauss->coord2*zeta*(zeta-1.);
    266                         basis[10]=2.*gauss->coord3*gauss->coord1*zeta*(zeta-1.);
    267                         basis[11]=2.*gauss->coord1*gauss->coord2*zeta*(zeta-1.);
    268                         basis[12]=2.*gauss->coord3*gauss->coord2*zeta*(zeta+1.);
    269                         basis[13]=2.*gauss->coord3*gauss->coord1*zeta*(zeta+1.);
    270                         basis[14]=2.*gauss->coord1*gauss->coord2*zeta*(zeta+1.);
     265                        basis[ 9]=4.*gauss->coord3*gauss->coord2*zeta*(zeta-1.)/2.;
     266                        basis[10]=4.*gauss->coord3*gauss->coord1*zeta*(zeta-1.)/2.;
     267                        basis[11]=4.*gauss->coord1*gauss->coord2*zeta*(zeta-1.)/2.;
     268                        basis[12]=4.*gauss->coord3*gauss->coord2*zeta*(zeta+1.)/2.;
     269                        basis[13]=4.*gauss->coord3*gauss->coord1*zeta*(zeta+1.)/2.;
     270                        basis[14]=4.*gauss->coord1*gauss->coord2*zeta*(zeta+1.)/2.;
     271                        /*quad faces*/
     272                        basis[15]=4.*gauss->coord3*gauss->coord2*(1.-zeta*zeta);
     273                        basis[16]=4.*gauss->coord3*gauss->coord1*(1.-zeta*zeta);
     274                        basis[17]=4.*gauss->coord1*gauss->coord2*(1.-zeta*zeta);
    271275                        return;
    272276                case P2xP4Enum :
    … …  
    599603
    600604                        /*Nodal function 7*/
    601                         dbasis[NUMNODESP2*0+6 ] = -0.5*(1.-zeta*zeta);
    602                         dbasis[NUMNODESP2*1+6 ] = -SQRT3/6.*(1.-zeta*zeta);
    603                         dbasis[NUMNODESP2*2+6 ] = -2.*zeta*gauss->coord1;
     605                        dbasis[NUMNODESP2*0+6 ] = (-2.*gauss->coord1 + 0.5)*(1.-zeta*zeta);
     606                        dbasis[NUMNODESP2*1+6 ] = (-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.)*(1.-zeta*zeta);
     607                        dbasis[NUMNODESP2*2+6 ] = -2.*zeta*gauss->coord1*(2.*gauss->coord1-1.);
    604608                        /*Nodal function 8*/
    605                         dbasis[NUMNODESP2*0+7 ] = 0.5*(1.-zeta*zeta);
    606                         dbasis[NUMNODESP2*1+7 ] = -SQRT3/6.*(1.-zeta*zeta);
    607                         dbasis[NUMNODESP2*2+7 ] = -2.*zeta*gauss->coord2;
     609                        dbasis[NUMNODESP2*0+7 ] = (+2.*gauss->coord2 - 0.5)*(1.-zeta*zeta);
     610                        dbasis[NUMNODESP2*1+7 ] = (-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.)*(1.-zeta*zeta);
     611                        dbasis[NUMNODESP2*2+7 ] = -2.*zeta*gauss->coord2*(2.*gauss->coord2-1.);
    608612                        /*Nodal function 9*/
    609613                        dbasis[NUMNODESP2*0+8 ] = 0.;
    610                         dbasis[NUMNODESP2*1+8 ] = SQRT3/3.*(1.-zeta*zeta);
    611                         dbasis[NUMNODESP2*2+8 ] = -2.*zeta*gauss->coord3;
     614                        dbasis[NUMNODESP2*1+8 ] = (+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.)*(1.-zeta*zeta);
     615                        dbasis[NUMNODESP2*2+8 ] = -2.*zeta*gauss->coord3*(2.*gauss->coord3-1.);
    612616
    613617                        /*Nodal function 10*/
    … …  
    635639                        dbasis[NUMNODESP2*1+14] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
    636640                        dbasis[NUMNODESP2*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.);
     641
     642                        /*Nodal function 16*/
     643                        dbasis[NUMNODESP2*0+15] = 2.*gauss->coord3*(1.-zeta*zeta);
     644                        dbasis[NUMNODESP2*1+15] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
     645                        dbasis[NUMNODESP2*2+15] = -2.*zeta*4.*gauss->coord3*gauss->coord2;
     646                        /*Nodal function 17*/
     647                        dbasis[NUMNODESP2*0+16] = -2.*gauss->coord3*(1.-zeta*zeta);
     648                        dbasis[NUMNODESP2*1+16] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
     649                        dbasis[NUMNODESP2*2+16] = -2.*zeta*4.*gauss->coord3*gauss->coord1;
     650                        /*Nodal function 18*/
     651                        dbasis[NUMNODESP2*0+17] = 2.*(gauss->coord1-gauss->coord2)*(1.-zeta*zeta);
     652                        dbasis[NUMNODESP2*1+17] = (1.-zeta*zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
     653                        dbasis[NUMNODESP2*2+17] = -2.*zeta*4.*gauss->coord1*gauss->coord2;
    637654                        return;
    638655                case P2xP4Enum :

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r17057 r17383  
    654654                                case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
    655655                                case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
    656                                 default: _error_("node index should be in [0 14]");
     656
     657                                case 15: coord1=0.; coord2=.5; coord3=.5; coord4=0.;break;
     658                                case 16: coord1=.5; coord2=0.; coord3=.5; coord4=0.;break;
     659                                case 17: coord1=.5; coord2=.5; coord3=0.; coord4=0.;break;
     660                                default: _error_("node index should be in [0 17]");
    657661                        }
    658662                        break;

  • issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

    r17138 r17383  
    3333
    3434        /*intermediary: */
    35         int         i,j,count,elementnbv;
     35        bool        isnan;
     36        int         i,j,count,elementnbv,numfacevertices;
     37        int*        faceverticesid   = NULL;
    3638        IssmDouble  value;
    3739        IssmDouble *times            = NULL;
    … …  
    7072                case P2Enum:
    7173                        EdgesPartitioning(&my_edges,iomodel);
     74                        FacesPartitioning(&my_faces,iomodel);
    7275                        break;
    7376                case P2xP4Enum:
    … …  
    109112                                                                                        dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
    110113                                                        count++;
     114                                                }
     115                                        }
     116                                }
     117                                for(i=0;i<iomodel->numberoffaces;i++){
     118                                        if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     119                                                if(my_faces[i]){
     120                                                        FaceGetVertexIndices(iomodel,&numfacevertices,&faceverticesid,i);
     121                                                        isnan=0;
     122                                                        for(j=0;j<numfacevertices;j++){
     123                                                                if(xIsNan<IssmDouble>(spcdata[faceverticesid[j]])) isnan=1;
     124                                                        }
     125                                                        if(isnan==0){
     126                                                                value=0;
     127                                                                for(j=0;j<numfacevertices;j++){
     128                                                                        value=value+spcdata[faceverticesid[j]]/numfacevertices;
     129                                                                }
     130                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,
     131                                                                                                dof,value,analysis_type));
     132                                                                count++;
     133                                                        }
    111134                                                }
    112135                                        }
    … …  
    431454
    432455        /*Free ressources:*/
     456        xDelete<int>(faceverticesid);
    433457        xDelete<IssmDouble>(times);
    434458        xDelete<IssmDouble>(values);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp

    r17358 r17383  
    240240        iomodel->elementtofaceconnectivity = element_face_connectivity;
    241241}/*}}}*/
     242void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber){/*{{{*/
     243
     244        int numbervertices;
     245        if(iomodel->meshtype==Mesh3DEnum){
     246                if((iomodel->faces[6*facenumber+5]-1)==1){
     247                        numbervertices=3;
     248                }
     249                else if((iomodel->faces[6*facenumber+5]-1)==2){
     250                        numbervertices=4;
     251                }
     252                else _error_("face marker not supported yet ");
     253        }
     254        else _error_("mesh type not supported yet");
     255
     256        int* facevertices = xNew<int>(numbervertices);
     257        if(numbervertices==3){
     258                for(int i=0;i<3;i++) facevertices[i]=iomodel->faces[6*facenumber+i]-1;
     259        }
     260        else if(numbervertices==4){
     261                int  i,j,cols,faceid;
     262                int  maxnbf,nbf,elementnbf,elementnbv,facemaxnbv;
     263                int *elementfaces         = NULL;
     264                int *elementfaces_markers = NULL;
     265                int elementid=iomodel->faces[6*facenumber+4]-1;
     266                int counter=0;
     267
     268                /*Recreate element properties*/
     269                elementnbv = 6; /*Number of vertices per element*/
     270                elementnbf = 5; /*Number of faces per element*/
     271                facemaxnbv = 4; /*Maximum number of vertices per face*/
     272                cols       = facemaxnbv + 1;
     273                elementfaces         = xNew<int>(elementnbf*cols);
     274                elementfaces_markers = xNew<int>(elementnbf);
     275                /*2 triangles*/
     276                elementfaces_markers[0] = 1;
     277                elementfaces_markers[1] = 1;
     278                elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0;  elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
     279                elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3;  elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5;
     280                /*3 quads*/
     281                elementfaces_markers[2] = 2;
     282                elementfaces_markers[3] = 2;
     283                elementfaces_markers[4] = 2;
     284                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;
     285                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;
     286                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;
     287
     288                for(faceid=2;faceid<5;faceid++){
     289                        counter=0;
     290                        for(j=0;j<3;j++){
     291                                if(iomodel->elements[elementid*6+elementfaces[cols*faceid+j]] == iomodel->faces[6*facenumber+j]) counter++;
     292                        }
     293                        if(counter==3) break;
     294                }
     295                if(!counter==3) _error_("face not found in element");
     296
     297                for(j=0;j<3;j++) facevertices[i]=iomodel->elements[elementid*6+elementfaces[cols*faceid+j]];
     298
     299                /*Delete*/
     300                xDelete<int>(elementfaces);
     301                xDelete<int>(elementfaces_markers);
     302        }
     303
     304        /*Output results*/
     305        *pverticesid=facevertices;
     306        *pnumvertices=numbervertices;
     307}/*}}}*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp

    r17137 r17383  
    132132                case P2Enum:
    133133                        EdgesPartitioning(&my_edges,iomodel);
     134                        FacesPartitioning(&my_faces,iomodel);
    134135                        for(i=0;i<iomodel->numberofvertices;i++){
    135136                                if(iomodel->my_vertices[i]){
    … …  
    140141                                if(my_edges[i]){
    141142                                        nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,approximation));
     143                                }
     144                        }
     145                        id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
     146                        for(i=0;i<iomodel->numberoffaces;i++){
     147                                if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     148                                        if(my_faces[i]){
     149                                                node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,approximation);
     150                                                nodes->AddObject(node);
     151                                        }
     152                                }
     153                                else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
     154                                        /*Nothing*/
     155                                }
     156                                else{
     157                                        _error_("not supported");
    142158                                }
    143159                        }
    … …  
    278294                        /*P2 velocity*/
    279295                        EdgesPartitioning(&my_edges,iomodel);
     296                        FacesPartitioning(&my_faces,iomodel);
    280297                        for(i=0;i<iomodel->numberofvertices;i++){
    281298                                if(iomodel->my_vertices[i]){
    … …  
    288305                                }
    289306                        }
     307                        id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
     308                        for(i=0;i<iomodel->numberoffaces;i++){
     309                                if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     310                                        if(my_faces[i]){
     311                                                node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
     312                                                nodes->AddObject(node);
     313                                        }
     314                                }
     315                                else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
     316                                        /*Nothing*/
     317                                }
     318                                else{
     319                                        _error_("not supported");
     320                                }
     321                        }
    290322
    291323                        /*P1 pressure*/
    292                         vnodes = id0+iomodel->numberofvertices+iomodel->numberofedges;
    293                         for(i=0;i<iomodel->numberofvertices;i++){
    294                                 if(iomodel->my_vertices[i]){
    295                                         nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,i,iomodel,analysis,FSpressureEnum));
     324                        vnodes = id0+iomodel->numberoffaces;
     325                        for(i=0;i<iomodel->numberofvertices;i++){
     326                                if(iomodel->my_vertices[i]){
     327                                        nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+i,lid++,i,iomodel,analysis,FSpressureEnum));
    296328                                }
    297329                        }

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r17122 r17383  
    3333void CreateFaces(IoModel* iomodel);
    3434void CreateFaces3d(IoModel* iomodel);
     35void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber);
    3536
    3637/*Connectivity*/
Note: See TracChangeset for help on using the changeset viewer.