Changeset 15637


Ignore:
Timestamp:
07/26/13 10:03:15 (12 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixing P1xP2

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15636 r15637  
    10961096        _assert_(basis);
    10971097
     1098        /*Get current coordinates in reference element*/
     1099        IssmDouble zeta=gauss->coord4;
     1100
    10981101        switch(this->element_type){
    10991102                case P1Enum:
    11001103                case P1DGEnum:
    1101                         basis[0]=gauss->coord1*(1.-gauss->coord4)/2.;
    1102                         basis[1]=gauss->coord2*(1.-gauss->coord4)/2.;
    1103                         basis[2]=gauss->coord3*(1.-gauss->coord4)/2.;
    1104                         basis[3]=gauss->coord1*(1.+gauss->coord4)/2.;
    1105                         basis[4]=gauss->coord2*(1.+gauss->coord4)/2.;
    1106                         basis[5]=gauss->coord3*(1.+gauss->coord4)/2.;
     1104                        basis[0]=gauss->coord1*(1.-zeta)/2.;
     1105                        basis[1]=gauss->coord2*(1.-zeta)/2.;
     1106                        basis[2]=gauss->coord3*(1.-zeta)/2.;
     1107                        basis[3]=gauss->coord1*(1.+zeta)/2.;
     1108                        basis[4]=gauss->coord2*(1.+zeta)/2.;
     1109                        basis[5]=gauss->coord3*(1.+zeta)/2.;
    11071110                        return;
    11081111                case P1bubbleEnum:
    1109                         basis[0]=gauss->coord1*(1.-gauss->coord4)/2.;
    1110                         basis[1]=gauss->coord2*(1.-gauss->coord4)/2.;
    1111                         basis[2]=gauss->coord3*(1.-gauss->coord4)/2.;
    1112                         basis[3]=gauss->coord1*(1.+gauss->coord4)/2.;
    1113                         basis[4]=gauss->coord2*(1.+gauss->coord4)/2.;
    1114                         basis[5]=gauss->coord3*(1.+gauss->coord4)/2.;
    1115                         basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+gauss->coord4)*(1.-gauss->coord4);
     1112                        basis[0]=gauss->coord1*(1.-zeta)/2.;
     1113                        basis[1]=gauss->coord2*(1.-zeta)/2.;
     1114                        basis[2]=gauss->coord3*(1.-zeta)/2.;
     1115                        basis[3]=gauss->coord1*(1.+zeta)/2.;
     1116                        basis[4]=gauss->coord2*(1.+zeta)/2.;
     1117                        basis[5]=gauss->coord3*(1.+zeta)/2.;
     1118                        basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta);
    11161119                        return;
    11171120                case P2xP1Enum:
    11181121                        /*Corner nodes*/
    1119                         basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-gauss->coord4)/2.;
    1120                         basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-gauss->coord4)/2.;
    1121                         basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-gauss->coord4)/2.;
    1122                         basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*(1.+gauss->coord4)/2.;
    1123                         basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*(1.+gauss->coord4)/2.;
    1124                         basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*(1.+gauss->coord4)/2.;
     1122                        basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta)/2.;
     1123                        basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta)/2.;
     1124                        basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta)/2.;
     1125                        basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*(1.+zeta)/2.;
     1126                        basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*(1.+zeta)/2.;
     1127                        basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*(1.+zeta)/2.;
    11251128                        /*mid-sides of triangles*/
    1126                         basis[ 6]=4.*gauss->coord3*gauss->coord2*(1.-gauss->coord4)/2.;
    1127                         basis[ 7]=4.*gauss->coord3*gauss->coord1*(1.-gauss->coord4)/2.;
    1128                         basis[ 8]=4.*gauss->coord1*gauss->coord2*(1.-gauss->coord4)/2.;
    1129                         basis[ 9]=4.*gauss->coord3*gauss->coord2*(1.+gauss->coord4)/2.;
    1130                         basis[10]=4.*gauss->coord3*gauss->coord1*(1.+gauss->coord4)/2.;
    1131                         basis[11]=4.*gauss->coord1*gauss->coord2*(1.+gauss->coord4)/2.;
     1129                        basis[ 6]=4.*gauss->coord3*gauss->coord2*(1.-zeta)/2.;
     1130                        basis[ 7]=4.*gauss->coord3*gauss->coord1*(1.-zeta)/2.;
     1131                        basis[ 8]=4.*gauss->coord1*gauss->coord2*(1.-zeta)/2.;
     1132                        basis[ 9]=4.*gauss->coord3*gauss->coord2*(1.+zeta)/2.;
     1133                        basis[10]=4.*gauss->coord3*gauss->coord1*(1.+zeta)/2.;
     1134                        basis[11]=4.*gauss->coord1*gauss->coord2*(1.+zeta)/2.;
    11321135                        return;
    11331136                case P1xP2Enum:
    11341137                        /*Corner nodes*/
    1135                         basis[ 0]=gauss->coord1*gauss->coord4*(gauss->coord4-1.)/2.;
    1136                         basis[ 1]=gauss->coord2*gauss->coord4*(gauss->coord4-1.)/2.;
    1137                         basis[ 2]=gauss->coord3*gauss->coord4*(gauss->coord4-1.)/2.;
    1138                         basis[ 3]=gauss->coord1*gauss->coord4*(gauss->coord4+1.)/2.;
    1139                         basis[ 4]=gauss->coord2*gauss->coord4*(gauss->coord4+1.)/2.;
    1140                         basis[ 5]=gauss->coord3*gauss->coord4*(gauss->coord4+1.)/2.;
     1138                        basis[ 0]=gauss->coord1*zeta*(zeta-1.)/2.;
     1139                        basis[ 1]=gauss->coord2*zeta*(zeta-1.)/2.;
     1140                        basis[ 2]=gauss->coord3*zeta*(zeta-1.)/2.;
     1141                        basis[ 3]=gauss->coord1*zeta*(zeta+1.)/2.;
     1142                        basis[ 4]=gauss->coord2*zeta*(zeta+1.)/2.;
     1143                        basis[ 5]=gauss->coord3*zeta*(zeta+1.)/2.;
    11411144                        /*mid-sides of quads*/
    1142                         basis[ 6]=gauss->coord1*(1.-gauss->coord4*gauss->coord4);
    1143                         basis[ 7]=gauss->coord2*(1.-gauss->coord4*gauss->coord4);
    1144                         basis[ 8]=gauss->coord3*(1.-gauss->coord4*gauss->coord4);
     1145                        basis[ 6]=gauss->coord1*(1.-zeta*zeta);
     1146                        basis[ 7]=gauss->coord2*(1.-zeta*zeta);
     1147                        basis[ 8]=gauss->coord3*(1.-zeta*zeta);
    11451148                        return;
    11461149                case P2Enum:
    11471150                        /*Corner nodes*/
    1148                         basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*gauss->coord4*(gauss->coord4-1.)/2.;
    1149                         basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*gauss->coord4*(gauss->coord4-1.)/2.;
    1150                         basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*gauss->coord4*(gauss->coord4-1.)/2.;
    1151                         basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*gauss->coord4*(gauss->coord4+1.)/2.;
    1152                         basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*gauss->coord4*(gauss->coord4+1.)/2.;
    1153                         basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*gauss->coord4*(gauss->coord4+1.)/2.;
     1151                        basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta-1.)/2.;
     1152                        basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta-1.)/2.;
     1153                        basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta-1.)/2.;
     1154                        basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta+1.)/2.;
     1155                        basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta+1.)/2.;
     1156                        basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.;
    11541157                        /*mid-sides of quads*/
    1155                         basis[ 6]=gauss->coord1*(1.-gauss->coord4*gauss->coord4);
    1156                         basis[ 7]=gauss->coord2*(1.-gauss->coord4*gauss->coord4);
    1157                         basis[ 8]=gauss->coord3*(1.-gauss->coord4*gauss->coord4);
     1158                        basis[ 6]=gauss->coord1*(1.-zeta*zeta);
     1159                        basis[ 7]=gauss->coord2*(1.-zeta*zeta);
     1160                        basis[ 8]=gauss->coord3*(1.-zeta*zeta);
    11581161                        /*mid-sides of triangles*/
    1159                         basis[ 9]=2.*gauss->coord3*gauss->coord2*gauss->coord4*(gauss->coord4-1.);
    1160                         basis[10]=2.*gauss->coord3*gauss->coord1*gauss->coord4*(gauss->coord4-1.);
    1161                         basis[11]=2.*gauss->coord1*gauss->coord2*gauss->coord4*(gauss->coord4-1.);
    1162                         basis[12]=2.*gauss->coord3*gauss->coord2*gauss->coord4*(gauss->coord4+1.);
    1163                         basis[13]=2.*gauss->coord3*gauss->coord1*gauss->coord4*(gauss->coord4+1.);
    1164                         basis[14]=2.*gauss->coord1*gauss->coord2*gauss->coord4*(gauss->coord4+1.);
     1162                        basis[ 9]=2.*gauss->coord3*gauss->coord2*zeta*(zeta-1.);
     1163                        basis[10]=2.*gauss->coord3*gauss->coord1*zeta*(zeta-1.);
     1164                        basis[11]=2.*gauss->coord1*gauss->coord2*zeta*(zeta-1.);
     1165                        basis[12]=2.*gauss->coord3*gauss->coord2*zeta*(zeta+1.);
     1166                        basis[13]=2.*gauss->coord3*gauss->coord1*zeta*(zeta+1.);
     1167                        basis[14]=2.*gauss->coord1*gauss->coord2*zeta*(zeta+1.);
    11651168                        return;
    11661169                default:
    … …  
    13251328                case P1xP2Enum:
    13261329                        /*Nodal function 1*/
    1327                         dbasis[NUMNODESP1xP2*0+0]   = zeta*(zeta-1.)/4.;
    1328                         dbasis[NUMNODESP1xP2*1+0]   = SQRT3/12.*zeta*(zeta-1.);
    1329                         dbasis[NUMNODESP1xP2*2+0]   = .5*(1.-2.*zeta)*gauss->coord1;
     1330                        dbasis[NUMNODESP1xP2*0+0]   = -zeta*(zeta-1.)/4.;
     1331                        dbasis[NUMNODESP1xP2*1+0]   = -SQRT3/12.*zeta*(zeta-1.);
     1332                        dbasis[NUMNODESP1xP2*2+0]   = .5*(2.*zeta-1.)*gauss->coord1;
    13301333                        /*Nodal function 2*/
    1331                         dbasis[NUMNODESP1xP2*0+1]   = zeta*(1.-zeta)/4.;
    1332                         dbasis[NUMNODESP1xP2*1+1]   = SQRT3/12.*zeta*(zeta-1);
    1333                         dbasis[NUMNODESP1xP2*2+1]   = .5*(1.-2.*zeta)*gauss->coord2;
     1334                        dbasis[NUMNODESP1xP2*0+1]   = zeta*(zeta-1.)/4.;
     1335                        dbasis[NUMNODESP1xP2*1+1]   = -SQRT3/12.*zeta*(zeta-1);
     1336                        dbasis[NUMNODESP1xP2*2+1]   = .5*(2.*zeta-1.)*gauss->coord2;
    13341337                        /*Nodal function 3*/
    13351338                        dbasis[NUMNODESP1xP2*0+2]   = 0.;
    1336                         dbasis[NUMNODESP1xP2*1+2]   = SQRT3/6.*zeta*(1.-zeta);
    1337                         dbasis[NUMNODESP1xP2*2+2]   = .5*(1.-2.*zeta)*gauss->coord3;
     1339                        dbasis[NUMNODESP1xP2*1+2]   = SQRT3/6.*zeta*(zeta-1.);
     1340                        dbasis[NUMNODESP1xP2*2+2]   = .5*(2.*zeta-1.)*gauss->coord3;
    13381341                        /*Nodal function 4*/
    1339                         dbasis[NUMNODESP1xP2*0+3]   = -zeta*(1.+zeta)/4.;
    1340                         dbasis[NUMNODESP1xP2*1+3]   = -SQRT3/12.*(1.+zeta);
    1341                         dbasis[NUMNODESP1xP2*2+3]   = .5*(1.+2.*zeta)*gauss->coord1;
     1342                        dbasis[NUMNODESP1xP2*0+3]   = -zeta*(zeta+1)/4.;
     1343                        dbasis[NUMNODESP1xP2*1+3]   = -SQRT3/12.*zeta*(zeta+1.);
     1344                        dbasis[NUMNODESP1xP2*2+3]   = .5*(2.*zeta+1.)*gauss->coord1;
    13421345                        /*Nodal function 5*/
    1343                         dbasis[NUMNODESP1xP2*0+4]   = zeta*(1.+zeta)/4.;
    1344                         dbasis[NUMNODESP1xP2*1+4]   = -SQRT3/12.*zeta*(1.+zeta);
    1345                         dbasis[NUMNODESP1xP2*2+4]   = .5*(1.+2.*zeta)*gauss->coord2;
     1346                        dbasis[NUMNODESP1xP2*0+4]   = zeta*(zeta+1.)/4.;
     1347                        dbasis[NUMNODESP1xP2*1+4]   = -SQRT3/12.*zeta*(zeta+1.);
     1348                        dbasis[NUMNODESP1xP2*2+4]   = .5*(2.*zeta+1.)*gauss->coord2;
    13461349                        /*Nodal function 6*/
    13471350                        dbasis[NUMNODESP1xP2*0+5]   = 0.;
    1348                         dbasis[NUMNODESP1xP2*1+5]   = SQRT3/6.*zeta*(1.+zeta);
    1349                         dbasis[NUMNODESP1xP2*2+5]   = .5*(1.+2.*zeta)*gauss->coord3;
     1351                        dbasis[NUMNODESP1xP2*1+5]   = SQRT3/6.*zeta*(zeta+1.);
     1352                        dbasis[NUMNODESP1xP2*2+5]   = .5*(2.*zeta+1.)*gauss->coord3;
    13501353
    13511354                        /*Nodal function 7*/
Note: See TracChangeset for help on using the changeset viewer.