Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15756) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15757) @@ -6876,19 +6876,19 @@ /*Constants*/ - const int numnodes = 2 *NUMVERTICES; - const int numdofm = NDOF2 *NUMVERTICES2D; - const int numdofs = NDOF4 *NUMVERTICES; - const int numdoftotal = 2 *numdofm+numdofs; + const int numdofm = NDOF2 *NUMVERTICES2D; + const int numdofs = NDOF4 *NUMVERTICES; + const int numdofstotal = NDOF4 *NUMVERTICES + NDOF3; + const int numdoftotal = 2 *numdofm+numdofstotal; /*Intermediaries */ - int i,j; + int i,j; IssmDouble Jdet; IssmDouble viscosity,FSreconditioning; //viscosity IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[4][numdofs+3]; + IssmDouble B[4][numdofstotal]; IssmDouble Bprime[4][numdofm]; IssmDouble B2[3][numdofm]; - IssmDouble Bprime2[3][numdofs+3]; + IssmDouble Bprime2[3][numdofstotal]; IssmDouble D[4][4]={0.0}; // material matrix, simple scalar matrix. IssmDouble D2[3][3]={0.0}; // material matrix, simple scalar matrix. @@ -6896,10 +6896,9 @@ IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix - IssmDouble Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point. - IssmDouble Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point. + IssmDouble Ke_gg_gaussian[numdofstotal][numdofm]; //stiffness matrix evaluated at the gaussian point. + IssmDouble Ke_gg_gaussian2[numdofm][numdofstotal]; //stiffness matrix evaluated at the gaussian point. GaussPenta *gauss=NULL; GaussTria *gauss_tria=NULL; - Node *node_list[numnodes]; - int cs_list[numnodes]; + Node *node_list[20]; /*Find penta on bed as FS must be coupled to the dofs on the bed: */ @@ -6907,16 +6906,28 @@ Tria* tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1. + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int numnodes = 2*vnumnodes-1+pnumnodes; + /*Prepare node list*/ - for(i=0;inodes[i]; - node_list[i+1*NUMVERTICES] = this->nodes[i]; - cs_list[i+0*NUMVERTICES] = XYEnum; - cs_list[i+1*NUMVERTICES] = XYZEnum; - } + int* cs_list = xNew(2*vnumnodes+pnumnodes); + for(i=0;inodes[i]; + cs_list[i] = XYEnum; + } + for(i=0;inodes[i]; + cs_list[i+vnumnodes-1] = XYZEnum; + } + for(i=0;inodes[vnumnodes+i]; + cs_list[2*vnumnodes-1+i] = PressureEnum; + } + /*Initialize Element matrix and return if necessary*/ - ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum); - ElementMatrix* Ke2=new ElementMatrix(this->nodes ,NUMVERTICES,this->parameters,FSvelocityEnum); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); + ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES, this->parameters,SSAApproximationEnum); + ElementMatrix* Ke2=new ElementMatrix(this->nodes ,2*NUMVERTICES+1,this->parameters,FSvelocityEnum); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); delete Ke1; delete Ke2; @@ -6950,5 +6961,5 @@ for (i=0;i<3;i++) D2[i][i]=D_scalar; - TripleMultiply( &B[0][0],4,numdofs+3,1, + TripleMultiply( &B[0][0],4,numdofstotal,1, &D[0][0],4,4,0, &Bprime[0][0],4,numdofm,0, @@ -6957,9 +6968,9 @@ TripleMultiply( &B2[0][0],3,numdofm,1, &D2[0][0],3,3,0, - &Bprime2[0][0],3,numdofs+3,0, + &Bprime2[0][0],3,numdofstotal,0, &Ke_gg_gaussian2[0][0],0); - for( i=0; ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; @@ -7111,6 +7122,5 @@ ElementMatrix* Ke2=new ElementMatrix(this->nodes,2*NUMVERTICES+1,this->parameters,FSvelocityEnum); ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); - delete Ke1; - delete Ke2; + delete Ke1; delete Ke2; /*Compute HO Matrix with P1 element type\n");*/ @@ -7132,5 +7142,5 @@ } - /*Transform Coordinate System*/ //Do not transform, already sone in the matrices + /*Transform Coordinate System*/ //Do not transform, already done in the matrices //TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); @@ -8116,5 +8126,12 @@ inputs->GetInputValue(&approximation,ApproximationEnum); if(approximation!=SSAFSApproximationEnum) return NULL; - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSvelocityEnum); + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSvelocityEnum); /*Retrieve all inputs and parameters*/ @@ -8142,15 +8159,16 @@ for(i=0;ivalues[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; - pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; - pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); - pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; + pe->values[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; + pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; + pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); + pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; } } /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum); + TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); /*Clean up and return*/ + xDelete(cs_list); delete gauss; return pe; Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15756) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15757) @@ -107,24 +107,38 @@ /*Build B: */ - for(i=0;i