Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15778) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15779) @@ -6924,5 +6924,4 @@ } - /*Initialize Element matrix and return if necessary*/ ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES, this->parameters,SSAApproximationEnum); @@ -6988,11 +6987,10 @@ /*FUNCTION Penta::CreateKMatrixCouplingSSAFSFriction {{{*/ ElementMatrix* Penta::CreateKMatrixCouplingSSAFSFriction(void){ - /*Constants*/ - const int numdof = NUMVERTICES *NDOF4; + const int numdofs = (NUMVERTICES+1)*NDOF3 + NUMVERTICES*NDOF1; const int numdofm = NUMVERTICES *NDOF2; - const int numdof2d = NUMVERTICES2D *NDOF4; + const int numdof2d = NUMVERTICES2D *NDOF3; const int numdof2dm = NUMVERTICES2D *NDOF2; - const int numdoftot = numdof+numdofm; + const int numdoftot = NUMVERTICES*2 + (NUMVERTICES+1)*3 +NUMVERTICES; // HO + FS vel + FS Pressure /*Intermediaries */ @@ -7006,14 +7004,14 @@ IssmDouble xyz_list_tria[NUMVERTICES2D][3]; IssmDouble LSSAFS[8][numdof2dm]; - IssmDouble LprimeSSAFS[8][numdof2d]; + IssmDouble LprimeSSAFS[8][numdofs]; IssmDouble DLSSAFS[8][8]={0.0}; IssmDouble LFSSSA[4][numdof2d]; IssmDouble LprimeFSSSA[4][numdof2dm]; IssmDouble DLFSSSA[4][4]={0.0}; - IssmDouble Ke_drag_gaussian[numdof2dm][numdof2d]; + IssmDouble Ke_drag_gaussian[numdof2dm][numdofs]; IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm]; Friction* friction=NULL; GaussPenta *gauss=NULL; - Node *node_list[20]; + Node *node_list[20]; /*If on water or not FS, skip stiffness: */ @@ -7025,18 +7023,23 @@ int numnodes = 2*vnumnodes-1+pnumnodes; - int* cs_list = xNew(numnodes); - - ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum); - ElementMatrix* Ke2=new ElementMatrix(this->nodes,numnodes,this->parameters,FSvelocityEnum); + /*Prepare node list*/ + 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; + } + + ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES, this->parameters,SSAApproximationEnum); + ElementMatrix* Ke2=new ElementMatrix(this->nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum); ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); delete Ke1; delete Ke2; - - /*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; - } /*Retrieve all inputs and parameters*/ @@ -7086,5 +7089,5 @@ TripleMultiply( &LSSAFS[0][0],8,numdof2dm,1, &DLSSAFS[0][0],8,8,0, - &LprimeSSAFS[0][0],8,numdof2d,0, + &LprimeSSAFS[0][0],8,numdofs,0, &Ke_drag_gaussian[0][0],0); @@ -7093,6 +7096,5 @@ &LprimeFSSSA[0][0],4,numdof2dm,0, &Ke_drag_gaussian2[0][0],0); - - for(i=0;ivalues[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j]; + for(i=0;ivalues[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j]; for(i=0;ivalues[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j]; } @@ -7118,5 +7120,5 @@ /*Intermediaries*/ - int i,j,init; + int i,j,init; Node *node_list[NUMVERTICES*3+1]; int cs_list[NUMVERTICES*3+1]; @@ -8220,5 +8222,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*/ @@ -8226,8 +8235,8 @@ parameters->FindParam(&analysis_type,AnalysisTypeEnum); this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* vzSSA_input=inputs->GetInput(VzSSAEnum); _assert_(vzSSA_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* vzSSA_input=inputs->GetInput(VzSSAEnum); _assert_(vzSSA_input); for(i=0;ivalues[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; - pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; - pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; + pe->values[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; + pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; + pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; } } /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum); + TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15778) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15779) @@ -981,5 +981,7 @@ * where h is the interpolation function for node i. */ - int num_dof=4; + int num_dof=3; + int num_dof_vel=3*NUMNODESP1b; + int num_dof_total=3*NUMNODESP1b+1*NUMNODESP1; IssmDouble L1L2l3[NUMNODESP1_2d]; IssmDouble dbasis[3][NUMNODESP1]; @@ -994,36 +996,74 @@ /*Build LprimeFS: */ for(int i=0;i<3;i++){ - LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i]; - LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i]; - LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i]; - LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i]; - LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i]; - LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i]; - LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = L1L2l3[i]; - LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.; - LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = L1L2l3[i]; + LprimeFS[num_dof_total*0+num_dof*i+0] = L1L2l3[i]; + LprimeFS[num_dof_total*0+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*0+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*1+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*1+num_dof*i+1] = L1L2l3[i]; + LprimeFS[num_dof_total*1+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*2+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*2+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*2+num_dof*i+2] = L1L2l3[i]; + LprimeFS[num_dof_total*3+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*3+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*3+num_dof*i+2] = L1L2l3[i]; + LprimeFS[num_dof_total*4+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*4+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*4+num_dof*i+2] = dbasis[2][i]; + LprimeFS[num_dof_total*5+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*5+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*5+num_dof*i+2] = dbasis[2][i]; + LprimeFS[num_dof_total*6+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*6+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*6+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*7+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*7+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*7+num_dof*i+2] = 0.; + } + for(int i=3;i<7;i++){ + LprimeFS[num_dof_total*0+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*0+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*0+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*1+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*1+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*1+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*2+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*2+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*2+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*3+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*3+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*3+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*4+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*4+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*4+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*5+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*5+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*5+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*6+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*6+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*6+num_dof*i+2] = 0.; + LprimeFS[num_dof_total*7+num_dof*i+0] = 0.; + LprimeFS[num_dof_total*7+num_dof*i+1] = 0.; + LprimeFS[num_dof_total*7+num_dof*i+2] = 0.; + } + for(int i=0;i<3;i++){ + LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*6+num_dof_vel+i] = L1L2l3[i]; + LprimeFS[num_dof_total*7+num_dof_vel+i] = L1L2l3[i]; + } + for(int i=3;i<6;i++){ + LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*6+num_dof_vel+i] = 0.; + LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.; } } @@ -1034,12 +1074,12 @@ * For node i, Li can be expressed in the actual coordinate system * by: - * Li=[ h 0 0 0] - * [ 0 h 0 0] - * [ 0 0 h 0] - * [ 0 0 h 0] + * Li=[ h 0 0 ] + * [ 0 h 0 ] + * [ 0 0 h ] + * [ 0 0 h ] * where h is the interpolation function for node i. */ - int num_dof=4; + int num_dof=3; IssmDouble L1L2l3[NUMNODESP1_2d]; @@ -1054,17 +1094,13 @@ LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.; LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.; - LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.; LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.; LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i]; LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.; - LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.; LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.; LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.; LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i]; - LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.; LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.; LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.; LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i]; - LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.; } }