Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15654) @@ -865,4 +865,53 @@ } /*}}}*/ +/*FUNCTION Penta::GetDofListVelocity{{{*/ +void Penta::GetDofListVelocity(int** pdoflist,int setenum){ + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodesVelocity(); + + /*First, figure out size of doflist and create it: */ + int numberofdofs=0; + for(int i=0;iGetNumberOfDofs(FSApproximationEnum,setenum); + + /*Allocate output*/ + int* doflist=xNew(numberofdofs); + + /*Populate: */ + int count=0; + for(int i=0;iGetDofList(doflist+count,FSApproximationEnum,setenum); + count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum); + } + + /*Assign output pointers:*/ + *pdoflist=doflist; +} +/*}}}*/ +/*FUNCTION Penta::GetDofListPressure{{{*/ +void Penta::GetDofListPressure(int** pdoflist,int setenum){ + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + + /*First, figure out size of doflist and create it: */ + int numberofdofs=0; + for(int i=vnumnodes;iGetNumberOfDofs(FSApproximationEnum,setenum); + + /*Allocate output*/ + int* doflist=xNew(numberofdofs); + + /*Populate: */ + int count=0; + for(int i=vnumnodes;iGetDofList(doflist+count,FSApproximationEnum,setenum); + count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum); + } + + /*Assign output pointers:*/ + *pdoflist=doflist; +} +/*}}}*/ /*FUNCTION Penta::GetGroundedPortion{{{*/ IssmDouble Penta::GetGroundedPortion(IssmDouble* xyz_list){ @@ -1034,7 +1083,8 @@ _assert_(nodes); - for(int i=0;iNumberofNodes(); + + for(int i=0;inodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; break; + case P1P1Enum: case P1P1GLSEnum: case MINIcondensedEnum: + numnodes = 12; + penta_node_ids = xNew(numnodes); + penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; + penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1]; + penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2]; + penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3]; + penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4]; + penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5]; + penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+0]; + penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+1]; + penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+2]; + penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+3]; + penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+4]; + penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+5]; + break; + case MINIEnum: + numnodes = 13; + penta_node_ids = xNew(numnodes); + penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; + penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1]; + penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2]; + penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3]; + penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4]; + penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5]; + penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+index+1; + penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+0]; + penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+1]; + penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+2]; + penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+3]; + penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+4]; + penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+5]; + break; default: _error_("Finite element "<NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -5045,5 +5128,5 @@ /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum); /*Clean up and return*/ @@ -5103,12 +5186,197 @@ ElementVector* Penta::CreatePVectorAdjointFS(void){ - if (!IsOnSurface()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). - ElementVector* pe=tria->CreatePVectorAdjointFS(); - delete tria->material; delete tria; - - /*clean up and return*/ + if(!IsOnSurface()) return NULL; + + /*Intermediaries */ + int i,j,resp; + int *responses= NULL; + int num_responses; + IssmDouble Jdet2d; + IssmDouble obs_velocity_mag,velocity_mag; + IssmDouble dux,duy; + IssmDouble epsvel = 2.220446049250313e-16; + IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr */ + IssmDouble scalex = 0,scaley=0,scale=0,S=0; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF3; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSApproximationEnum); + IssmDouble* vbasis = xNew(vnumdof); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); + Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vxobs_input = inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); + Input* vyobs_input = inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + + for(i=0;iGetInputValue(&S,SurfaceAreaEnum); break; + } + } + + /* Start looping on the number of gaussian points: */ + GaussPenta* gauss=new GaussPenta(3,4,5,4); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss); + + /*Get all parameters at gaussian point*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vxobs_input->GetInputValue(&vxobs,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + GetNodalFunctionsVelocity(vbasis,gauss); + + /*Loop over all requested responses*/ + for(resp=0;respGetInputValue(&weight,gauss,resp); + + switch(responses[resp]){ + + case SurfaceAbsVelMisfitEnum: + /* + * 1 [ 2 2 ] + * J = --- | (u - u ) + (v - v ) | + * 2 [ obs obs ] + * + * dJ + * DU = - -- = (u - u ) + * du obs + */ + for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; + } + break; + case SurfaceRelVelMisfitEnum: + /* + * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] + * J = --- | ------------- (u - u ) + ------------- (v - v ) | + * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] + * obs obs + * + * dJ \bar{v}^2 + * DU = - -- = ------------- (u - u ) + * du (u + eps)^2 obs + * obs + */ + for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; + } + break; + case SurfaceLogVelMisfitEnum: + /* + * [ vel + eps ] 2 + * J = 4 \bar{v}^2 | log ( ----------- ) | + * [ vel + eps ] + * obs + * + * dJ 2 * log(...) + * DU = - -- = - 4 \bar{v}^2 ------------- u + * du vel^2 + eps + * + */ + for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; + } + break; + case SurfaceAverageVelMisfitEnum: + /* + * 1 2 2 + * J = --- sqrt( (u - u ) + (v - v ) ) + * S obs obs + * + * dJ 1 1 + * DU = - -- = - --- ----------- * 2 (u - u ) + * du S 2 sqrt(...) obs + */ + for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; + } + break; + case SurfaceLogVxVyMisfitEnum: + /* + * 1 [ |u| + eps 2 |v| + eps 2 ] + * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | + * 2 [ |u |+ eps |v |+ eps ] + * obs obs + * dJ 1 u 1 + * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ + * du |u| + eps |u| u + eps + */ + for(i=0;ivalues[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i]; + } + break; + case DragCoefficientAbsGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAbsGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAcrossGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAlongGradientEnum: + /*Nothing in P vector*/ + break; + case RheologyBbarAbsGradientEnum: + /*Nothing in P vector*/ + break; + default: + _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); + } + } + } + + /*Clean up and return*/ + xDelete(responses); + xDelete(cs_list); + xDelete(vbasis); + delete gauss; return pe; } @@ -5494,33 +5762,48 @@ void Penta::InputUpdateFromSolutionAdjointFS(IssmDouble* solution){ - const int numdof=NDOF4*NUMVERTICES; - - int i; - IssmDouble values[numdof]; - IssmDouble lambdax[NUMVERTICES]; - IssmDouble lambday[NUMVERTICES]; - IssmDouble lambdaz[NUMVERTICES]; - IssmDouble lambdap[NUMVERTICES]; - int* doflist=NULL; + int i; + int* vdoflist=NULL; + int* pdoflist=NULL; + IssmDouble FSreconditioning; + GaussPenta *gauss; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocityFinal(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF3; + int pnumdof = pnumnodes*NDOF1; + + /*Initialize values*/ + IssmDouble* vvalues = xNew(vnumdof); + IssmDouble* pvalues = xNew(pnumdof); + IssmDouble* lambdax = xNew(vnumnodes); + IssmDouble* lambday = xNew(vnumnodes); + IssmDouble* lambdaz = xNew(vnumnodes); + IssmDouble* lambdap = xNew(pnumnodes); /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + GetDofListVelocity(&vdoflist,GsetEnum); + GetDofListPressure(&pdoflist,GsetEnum); /*Use the dof list to index into the solution vector: */ - for(i=0;i(lambdax[i])) _error_("NaN found in solution vector"); - if(xIsNan(lambday[i])) _error_("NaN found in solution vector"); - if(xIsNan(lambdaz[i])) _error_("NaN found in solution vector"); - if(xIsNan(lambdap[i])) _error_("NaN found in solution vector"); - } + for(i=0;i(lambdax[i])) _error_("NaN found in solution vector"); + lambday[i] = vvalues[i*NDOF3+1]; if(xIsNan(lambday[i])) _error_("NaN found in solution vector"); + lambdaz[i] = vvalues[i*NDOF3+2]; if(xIsNan(lambdaz[i])) _error_("NaN found in solution vector"); + } + for(i=0;i(lambdap[i])) _error_("NaN found in solution vector"); + } + + /*Recondition pressure and compute vel: */ + this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum); + for(i=0;i(doflist); + xDelete(vdoflist); + xDelete(pdoflist); + xDelete(lambdap); + xDelete(lambdaz); + xDelete(lambday); + xDelete(lambdax); + xDelete(vvalues); + xDelete(pvalues); } /*}}}*/ @@ -6051,19 +6341,23 @@ /*output: */ ElementVector* De=NULL; - /*intermediary: */ + + /*Initialize Element vector and return if necessary*/ int approximation; - int i; - - /*Initialize Element vector and return if necessary*/ inputs->GetInputValue(&approximation,ApproximationEnum); if(approximation!=FSApproximationEnum) return NULL; - De=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum); - - for (i=0;ivalues[i*4+0]=VelocityEnum; - De->values[i*4+1]=VelocityEnum; - De->values[i*4+2]=VelocityEnum; - De->values[i*4+3]=PressureEnum; + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + + De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum); + + for(int i=0;ivalues[i*3+0]=VelocityEnum; + De->values[i*3+1]=VelocityEnum; + De->values[i*3+2]=VelocityEnum; + } + for(int i=0;ivalues[vnumnodes*3+i]=PressureEnum; } @@ -6317,5 +6611,5 @@ node_list[i+1*NUMVERTICES] = this->nodes[i]; cs_list[i+0*NUMVERTICES] = XYEnum; - cs_list[i+1*NUMVERTICES] = XYZPEnum; + cs_list[i+1*NUMVERTICES] = XYZEnum; } @@ -6427,5 +6721,5 @@ node_list[i+1*NUMVERTICES] = this->nodes[i]; cs_list[i+0*NUMVERTICES] = XYEnum; - cs_list[i+1*NUMVERTICES] = XYZPEnum; + cs_list[i+1*NUMVERTICES] = XYZEnum; } @@ -6516,5 +6810,5 @@ node_list[i+1*NUMVERTICES] = this->nodes[i]; cs_list[i+0*NUMVERTICES] = XYEnum; - cs_list[i+1*NUMVERTICES] = XYZPEnum; + cs_list[i+1*NUMVERTICES] = XYZEnum; } @@ -6526,5 +6820,5 @@ delete Ke2; Ke1=CreateKMatrixDiagnosticHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum); - Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZPEnum); + Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZEnum); for(i=0;iFindParam(&fe_FS,FlowequationFeFSEnum); - - switch(fe_FS){ - case 0: - /*compute all stiffness matrices for this element*/ - Ke1=CreateKMatrixDiagnosticFSViscous(); - Ke2=CreateKMatrixDiagnosticFSFriction(); - Ke =new ElementMatrix(Ke1,Ke2); - break; - case 1: - /*compute all stiffness matrices for this element*/ - Ke1=CreateKMatrixDiagnosticFSGLSViscous(); - Ke2=CreateKMatrixDiagnosticFSFriction(); - Ke =new ElementMatrix(Ke1,Ke2); - break; - default: - _error_("Finite element" << fe_FS << " not supported yet"); - } + ElementMatrix* Ke1 = NULL; + ElementMatrix* Ke2 = NULL; + ElementMatrix* Ke = NULL; + + /*compute all stiffness matrices for this element*/ + Ke1=CreateKMatrixDiagnosticFSViscous(); + Ke2=CreateKMatrixDiagnosticFSFriction(); + Ke =new ElementMatrix(Ke1,Ke2); /*clean-up and return*/ @@ -7129,69 +7409,6 @@ } /*}}}*/ -/*FUNCTION Penta::CreateKMatrixDiagnosticFSViscous {{{*/ -ElementMatrix* Penta::CreateKMatrixDiagnosticFSViscous(void){ - - /*Intermediaries */ - int i,approximation; - IssmDouble Jdet,viscosity,FSreconditioning; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble B[8][27]; - IssmDouble B_prime[8][27]; - IssmDouble D_scalar; - IssmDouble D[8][8]={0.0}; - IssmDouble Ke_temp[27][27]={0.0}; //for the six nodes and the bubble - GaussPenta *gauss=NULL; - - /*If on water or not FS, skip stiffness: */ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum); - 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); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBFS(&B[0][0],&xyz_list[0][0],gauss); - GetBprimeFS(&B_prime[0][0],&xyz_list[0][0],gauss); - - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - - D_scalar=gauss->weight*Jdet; - for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity; - for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning; - - TripleMultiply( &B[0][0],8,27,1, - &D[0][0],8,8,0, - &B_prime[0][0],8,27,0, - &Ke_temp[0][0],1); - } - - /*Condensation*/ - ReduceMatrixFS(Ke->values, &Ke_temp[0][0]); - //Ke->Echo(); - //_error_("stop"); - - /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixDiagnosticFSGLSViscous {{{*/ -ElementMatrix* Penta::CreateKMatrixDiagnosticFSGLSViscous(void){ +/*FUNCTION Penta::KMatrixGLSstabilization{{{*/ +void Penta::KMatrixGLSstabilization(ElementMatrix* Ke){ int numdof = NUMVERTICES*NDOF4; @@ -7223,9 +7440,4 @@ int c=3; //index of pressure - /*If on water or not FS, skip stiffness: */ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum); - /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -7239,5 +7451,4 @@ diameter=MinEdgeLength(xyz_list); - if(stabilization){ gauss=new GaussPenta(); for(int iv=0;iv<6;iv++){ @@ -7251,5 +7462,4 @@ } delete gauss; - } /* Start looping on the number of gaussian points: */ @@ -7260,6 +7470,4 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBFSGLS(&B[0][0],&xyz_list[0][0],gauss); - GetBprimeFSGLS(&B_prime[0][0],&xyz_list[0][0],gauss); this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); @@ -7267,64 +7475,121 @@ D_scalar=gauss->weight*Jdet; - for (i=0;i<6;i++) D[i][i]=D_scalar*2.*2.*viscosity; - for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning; - - TripleMultiply( &B[0][0],8,24,1, - &D[0][0],8,8,0, - &B_prime[0][0],8,24,0, - &Ke_temp[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_temp[i][j]; + /*Add stabilization*/ - if(stabilization){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - dmu[0]=0.; dmu[1]=0.; dmu[2]=0.; - mu = 2.*viscosity; - for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){ - SU[p][i][j]=0.; - SW[p][i][j]=0.; - } - for(p=0;p<6;p++){ - for(i=0;i<3;i++){ - SU[p][i][c] += dbasis[i][p]; - SW[p][c][i] += dbasis[i][p]; - for(j=0;j<3;j++){ - SU[p][i][i] += -dmu[j]*dbasis[j][p]; - SU[p][i][j] += -dmu[i]*dbasis[i][p]; - for(ii=0;ii<6;ii++){ - SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; - SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; - } - SW[p][i][i] += -dmu[j]*dbasis[j][p]; - SW[p][j][i] += -dmu[j]*dbasis[i][p]; - for(ii=0;ii<6;ii++){ - SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; - SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; - } + GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + dmu[0]=0.; dmu[1]=0.; dmu[2]=0.; + mu = 2.*viscosity; + for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){ + SU[p][i][j]=0.; + SW[p][i][j]=0.; + } + for(p=0;p<6;p++){ + for(i=0;i<3;i++){ + SU[p][i][c] += dbasis[i][p]; + SW[p][c][i] += dbasis[i][p]; + for(j=0;j<3;j++){ + SU[p][i][i] += -dmu[j]*dbasis[j][p]; + SU[p][i][j] += -dmu[i]*dbasis[i][p]; + for(ii=0;ii<6;ii++){ + SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; + SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; + } + SW[p][i][i] += -dmu[j]*dbasis[j][p]; + SW[p][j][i] += -dmu[j]*dbasis[i][p]; + for(ii=0;ii<6;ii++){ + SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; + SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; } } } - IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity); - for(p=0;p<6;p++){ - for(q=0;q<6;q++){ - for(i=0;i<4;i++){ - for(j=0;j<4;j++){ - Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j]; - Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j]; - Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j]; - } + } + IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity); + for(p=0;p<6;p++){ + for(q=0;q<6;q++){ + for(i=0;i<4;i++){ + for(j=0;j<4;j++){ + Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j]; + Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j]; + Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j]; } } } - - } + } + } + + /*Clean up*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticFSViscous {{{*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticFSViscous(void){ + + /*Intermediaries */ + int i,j,approximation; + IssmDouble Jdet,viscosity,FSreconditioning,D_scalar; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + GaussPenta *gauss=NULL; + + /*If on water or not FS, skip stiffness: */ + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSApproximationEnum); + IssmDouble* B = xNew(8*numdof); + IssmDouble* Bprime = xNew(8*numdof); + IssmDouble* D = xNewZeroInit(8*8); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum); + 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); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetBFS(B,&xyz_list[0][0],gauss); + GetBprimeFS(Bprime,&xyz_list[0][0],gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + material->GetViscosity3dFS(&viscosity,&epsilon[0]); + + D_scalar=gauss->weight*Jdet; + for(i=0;i<6;i++) D[i*8+i] = +D_scalar*2.*viscosity; + for(i=6;i<8;i++) D[i*8+i] = -D_scalar*FSreconditioning; + + TripleMultiply(B,8,numdof,1, + D,8,8,0, + Bprime,8,numdof,0, + Ke->values,1); } /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list); /*Clean up and return*/ delete gauss; + xDelete(cs_list); + xDelete(B); + xDelete(Bprime); + xDelete(D); return Ke; } @@ -7333,26 +7598,28 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticFSFriction(void){ - /*Constants*/ - const int numdof=NUMVERTICES*NDOF4; - const int numdof2d=NUMVERTICES2D*NDOF4; - /*Intermediaries */ - int i,j; - int analysis_type,approximation; - IssmDouble alpha2,Jdet2d; - IssmDouble FSreconditioning,viscosity; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble LFS[2][numdof2d]; - IssmDouble DLFS[2][2]={0.0}; - IssmDouble Ke_drag_gaussian[numdof2d][numdof2d]; - Friction* friction=NULL; - GaussPenta *gauss=NULL; + int i,j; + int analysis_type,approximation; + IssmDouble alpha2,Jdet2d; + IssmDouble FSreconditioning,viscosity; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + Friction *friction = NULL; + GaussPenta *gauss = NULL; /*If on water or not FS, skip stiffness: */ inputs->GetInputValue(&approximation,ApproximationEnum); - if(IsFloating() || !IsOnBed() || (approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum)) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum); + if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum); + IssmDouble* BFriction = xNew(2*numdof); + IssmDouble* D = xNewZeroInit(2*2); /*Retrieve all inputs and parameters*/ @@ -7375,28 +7642,28 @@ GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - GetLFS(&LFS[0][0], gauss); + GetLFS(BFriction,gauss); this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); - friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); - - DLFS[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx - DLFS[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy - - TripleMultiply( &LFS[0][0],2,numdof2d,1, - &DLFS[0][0],2,2,0, - &LFS[0][0],2,numdof2d,0, - &Ke_drag_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_drag_gaussian[i][j]; + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + + D[0*2+0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx + D[1*2+1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy + + TripleMultiply(BFriction,2,numdof,1, + D,2,2,0, + BFriction,2,numdof,0, + Ke->values,1); } /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ - //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum); /*Clean up and return*/ delete gauss; delete friction; + xDelete(BFriction); + xDelete(D); return Ke; } @@ -7576,5 +7843,5 @@ /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); + TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum); /*Clean up and return*/ @@ -7650,5 +7917,5 @@ /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); + TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum); /*Clean up and return*/ @@ -7727,5 +7994,5 @@ /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); + TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum); /*Clean up and return*/ @@ -7801,5 +8068,5 @@ /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); + TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum); /*Clean up and return*/ @@ -8147,29 +8414,14 @@ ElementVector* Penta::CreatePVectorDiagnosticFS(void){ - int fe_FS; ElementVector* pe1; ElementVector* pe2; ElementVector* pe3; ElementVector* pe; - parameters->FindParam(&fe_FS,FlowequationFeFSEnum); - - switch(fe_FS){ - case 0: - /*compute all stiffness matrices for this element*/ - pe1=CreatePVectorDiagnosticFSViscous(); - pe2=CreatePVectorDiagnosticFSShelf(); - pe3=CreatePVectorDiagnosticFSFront(); - pe =new ElementVector(pe1,pe2,pe3); - break; - case 1: - /*compute all stiffness matrices for this element*/ - pe1=CreatePVectorDiagnosticFSGLSViscous(); - pe2=CreatePVectorDiagnosticFSShelf(); - pe3=CreatePVectorDiagnosticFSFront(); - pe =new ElementVector(pe1,pe2,pe3); - break; - default: - _error_("Finite element" << fe_FS << " not supported yet"); - } + + /*compute all stiffness matrices for this element*/ + pe1=CreatePVectorDiagnosticFSViscous(); + pe2=CreatePVectorDiagnosticFSShelf(); + pe3=CreatePVectorDiagnosticFSFront(); + pe =new ElementVector(pe1,pe2,pe3); /*clean-up and return*/ @@ -8180,98 +8432,9 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticFSViscous {{{*/ -ElementVector* Penta::CreatePVectorDiagnosticFSViscous(void){ - - /*Constants*/ - const int numdofbubble=NDOF4*NUMVERTICES+NDOF3*1; - - /*Intermediaries*/ - int i,j; - int approximation; - IssmDouble Jdet,viscosity; - IssmDouble gravity,rho_ice,FSreconditioning; - IssmDouble forcex,forcey,forcez; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble l1l7[7]; //for the six nodes and the bubble - IssmDouble B[8][numdofbubble]; - IssmDouble B_prime[8][numdofbubble]; - IssmDouble B_prime_bubble[8][3]; - IssmDouble D[8][8]={0.0}; - IssmDouble D_scalar; - IssmDouble Pe_gaussian[numdofbubble]={0.0}; //for the six nodes and the bubble - IssmDouble Ke_temp[numdofbubble][3]={0.0}; //for the six nodes and the bubble - GaussPenta *gauss=NULL; - - /*Initialize Element vector and return if necessary*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum); - - /*Retrieve all inputs and parameters*/ - this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - 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* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); - Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); - Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBFS(&B[0][0],&xyz_list[0][0],gauss); - GetBprimeFS(&B_prime[0][0],&xyz_list[0][0], gauss); - GetNodalFunctionsMINI(&l1l7[0], gauss); - - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - - loadingforcex_input->GetInputValue(&forcex, gauss); - loadingforcey_input->GetInputValue(&forcey, gauss); - loadingforcez_input->GetInputValue(&forcez, gauss); - - for(i=0;iweight*l1l7[i]; - Pe_gaussian[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l7[i]; - Pe_gaussian[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l7[i]; - Pe_gaussian[i*NDOF4+2]+=rho_ice*forcez*Jdet*gauss->weight*l1l7[i]; - } - - /*Get bubble part of Bprime */ - for(i=0;i<8;i++) for(j=0;j<3;j++) B_prime_bubble[i][j]=B_prime[i][j+24]; - - D_scalar=gauss->weight*Jdet; - for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity; - for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning; - - TripleMultiply(&B[0][0],8,numdofbubble,1, - &D[0][0],8,8,0, - &B_prime_bubble[0][0],8,3,0, - &Ke_temp[0][0],1); - } - - /*Condensation*/ - ReduceVectorFS(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]); - - /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ /*FUNCTION Penta::CreatePVectorDiagnosticFSFront{{{*/ ElementVector* Penta::CreatePVectorDiagnosticFSFront(void){ /*Intermediaries */ + int i; IssmDouble ls[NUMVERTICES]; IssmDouble xyz_list[NUMVERTICES][3]; @@ -8293,17 +8456,28 @@ /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - int numdof = numnodes*NDOF4; IssmDouble rho_ice,rho_water,gravity; IssmDouble Jdet,z_g,water_pressure,air_pressure; IssmDouble surface_under_water,base_under_water,pressure; - GaussPenta* gauss; - IssmDouble* basis = xNew(numnodes); IssmDouble xyz_list_front[4][3]; IssmDouble area_coordinates[4][3]; IssmDouble normal[3]; - + GaussPenta* gauss; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF3; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSApproximationEnum); + IssmDouble* vbasis = xNew(vnumdof); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,FSApproximationEnum); rho_water=matpar->GetRhoWater(); rho_ice =matpar->GetRhoIce(); @@ -8323,5 +8497,5 @@ gauss->GaussPoint(ig); z_g=GetZcoord(gauss); - GetNodalFunctions(basis,gauss); + GetNodalFunctionsVelocity(vbasis,gauss); GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss); @@ -8330,23 +8504,23 @@ pressure = water_pressure + air_pressure; - for (int i=0;ivalues[4*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; - pe->values[4*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; - pe->values[4*i+2]+= pressure*Jdet*gauss->weight*normal[2]*basis[i]; - pe->values[4*i+3]+= 0; + for(i=0;ivalues[3*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i]; + pe->values[3*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i]; + pe->values[3*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i]; } } /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,numnodes,XYZPEnum); + TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); /*Clean up and return*/ - xDelete(basis); + xDelete(cs_list); + xDelete(vbasis); delete gauss; return pe; } /*}}}*/ -/*FUNCTION Penta::CreatePVectorDiagnosticFSGLSViscous {{{*/ -ElementVector* Penta::CreatePVectorDiagnosticFSGLSViscous(void){ +/*FUNCTION Penta::PVectorGLSstabilization{{{*/ +void Penta::PVectorGLSstabilization(ElementVector* pe){ /*Constants*/ @@ -8374,10 +8548,4 @@ int c=3; //index of pressure - /*Initialize Element vector and return if necessary*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum); - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum); - /*Retrieve all inputs and parameters*/ rho_ice=matpar->GetRhoIce(); @@ -8385,7 +8553,4 @@ B=material->GetB(); GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); - Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); - Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); @@ -8395,5 +8560,4 @@ diameter=MinEdgeLength(xyz_list); - if(stabilization){ gauss=new GaussPenta(); for(int iv=0;iv<6;iv++){ @@ -8407,5 +8571,4 @@ } delete gauss; - } /* Start looping on the number of gaussian points: */ @@ -8420,53 +8583,102 @@ material->GetViscosity3dFS(&viscosity,&epsilon[0]); - loadingforcex_input->GetInputValue(&forcex, gauss); - loadingforcey_input->GetInputValue(&forcey, gauss); - loadingforcez_input->GetInputValue(&forcez, gauss); - - for(i=0;ivalues[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l6[i]; - pe->values[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l6[i]; - pe->values[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l6[i]; - pe->values[i*NDOF4+2]+=rho_ice*forcez*Jdet*gauss->weight*l1l6[i]; - } - - /*Add stabilization*/ - if(stabilization){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - dmu[0]=0.; dmu[1]=0.; dmu[2]=0.; - mu = 2.*viscosity; - for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){ - SW[p][i][j]=0.; - } - for(p=0;p<6;p++){ - for(i=0;i<3;i++){ - SW[p][c][i] += dbasis[i][p]; - for(j=0;j<3;j++){ - SW[p][i][i] += -dmu[j]*dbasis[j][p]; - SW[p][j][i] += -dmu[j]*dbasis[i][p]; - for(ii=0;ii<6;ii++){ - SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; - SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; - } + GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + dmu[0]=0.; dmu[1]=0.; dmu[2]=0.; + mu = 2.*viscosity; + for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){ + SW[p][i][j]=0.; + } + for(p=0;p<6;p++){ + for(i=0;i<3;i++){ + SW[p][c][i] += dbasis[i][p]; + for(j=0;j<3;j++){ + SW[p][i][i] += -dmu[j]*dbasis[j][p]; + SW[p][j][i] += -dmu[j]*dbasis[i][p]; + for(ii=0;ii<6;ii++){ + SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; + SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; } } } - IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity); - for(p=0;p<6;p++){ - for(j=0;j<4;j++){ - pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcex*SW[p][j][0]; - pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcey*SW[p][j][1]; - pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*(forcez-gravity)*SW[p][j][2]; - } + } + IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity); + for(p=0;p<6;p++){ + for(j=0;j<4;j++){ + pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcex*SW[p][j][0]; + pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcey*SW[p][j][1]; + pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*(forcez-gravity)*SW[p][j][2]; } - + } + } + + /*Clean up*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticFSViscous {{{*/ +ElementVector* Penta::CreatePVectorDiagnosticFSViscous(void){ + + /*Intermediaries*/ + int i,j; + int approximation; + IssmDouble Jdet,gravity,rho_ice; + IssmDouble forcex,forcey,forcez; + IssmDouble xyz_list[NUMVERTICES][3]; + GaussPenta *gauss=NULL; + + /*Initialize Element vector and return if necessary*/ + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF3; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSApproximationEnum); + IssmDouble* vbasis = xNew(vnumdof); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); + Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); + Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); + rho_ice = matpar->GetRhoIce(); + gravity = matpar->GetG(); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsVelocity(vbasis, gauss); + + loadingforcex_input->GetInputValue(&forcex,gauss); + loadingforcey_input->GetInputValue(&forcey,gauss); + loadingforcez_input->GetInputValue(&forcez,gauss); + + for(i=0;ivalues[i*NDOF3+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; + pe->values[i*NDOF3+2] += +rho_ice*forcez *Jdet*gauss->weight*vbasis[i]; } } /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); + TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); /*Clean up and return*/ delete gauss; + xDelete(cs_list); + xDelete(vbasis); return pe; } @@ -8484,5 +8696,4 @@ IssmDouble bed_normal[3]; IssmDouble dz[3]; - IssmDouble basis[6]; //for the six nodes of the penta IssmDouble Jdet2d; GaussPenta *gauss=NULL; @@ -8491,9 +8702,22 @@ if(!IsOnBed() || !IsFloating()) return NULL; inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF3; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iparameters,FSApproximationEnum); + IssmDouble* vbasis = xNew(vnumdof); + + /*Retrieve all inputs and parameters*/ this->parameters->FindParam(&shelf_dampening,DiagnosticShelfDampeningEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum); - - /*Retrieve all inputs and parameters*/ rho_water=matpar->GetRhoWater(); gravity=matpar->GetG(); @@ -8513,5 +8737,5 @@ GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(basis, gauss); + GetNodalFunctionsVelocity(vbasis, gauss); BedNormal(&bed_normal[0],xyz_list_tria); @@ -8522,19 +8746,25 @@ vy_input->GetInputValue(&vy, gauss); vz_input->GetInputValue(&vz, gauss); - dt=0; + dt=0.; normal_vel=bed_normal[0]*vx+bed_normal[1]*vy+bed_normal[2]*vz; damper=gravity*rho_water*pow(1+pow(dz[0],2)+pow(dz[1],2),0.5)*normal_vel*dt; } - else damper=0; + else damper=0.; water_pressure=gravity*rho_water*bed; - for(i=0;ivalues[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*basis[i]*bed_normal[j]; + for(i=0;ivalues[3*i+0]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[0]; + pe->values[3*i+1]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[1]; + pe->values[3*i+2]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[2]; + } } /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); + TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); /*Clean up and return*/ delete gauss; + xDelete(cs_list); + xDelete(vbasis); return pe; } @@ -8794,9 +9024,6 @@ ElementMatrix* Penta::CreateJacobianDiagnosticFS(void){ - /*Constants*/ - const int numdof=NDOF4*NUMVERTICES; - /*Intermediaries */ - int i,j; + int i,j,approximation; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble Jdet; @@ -8805,17 +9032,31 @@ IssmDouble eps3dotdphii,eps3dotdphij; IssmDouble mu_prime; - IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ IssmDouble eps1[3],eps2[3],eps3[3]; - IssmDouble dphi[3][NUMVERTICES]; GaussPenta *gauss=NULL; + + /*If on water or not FS, skip stiffness: */ + //inputs->GetInputValue(&approximation,ApproximationEnum); + //if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;i(3*vnumnodes); /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - 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* 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); /* Start looping on the number of gaussian points: */ @@ -8825,7 +9066,12 @@ gauss->GaussPoint(ig); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - + GetJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + GetNodalFunctionsDerivativesVelocity(dbasis,&xyz_list[0][0],gauss); + + //this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + //material->GetViscosityDerivativeEpsSquareFS(&mu_prime,&epsilon[0]); + //eps1[0]=epsilon[0]; eps2[0]=epsilon[3]; eps3[0]=epsilon[4]; + //eps1[1]=epsilon[3]; eps2[1]=epsilon[1]; eps3[1]=epsilon[5]; + //eps1[2]=epsilon[4]; eps2[2]=epsilon[5]; eps3[2]=epsilon[2]; this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); @@ -8834,24 +9080,24 @@ eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; - for(i=0;i<6;i++){ - for(j=0;j<6;j++){ - eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; - eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; - eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; - eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; - eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i]; - eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j]; - - Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; - Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; - Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; - - Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; - Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; - Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; - - Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; - Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; - Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; + for(i=0;ivalues[numdof*(3*i+0)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[numdof*(3*i+0)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[numdof*(3*i+0)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; + + Ke->values[numdof*(3*i+1)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[numdof*(3*i+1)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + Ke->values[numdof*(3*i+1)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; + + Ke->values[numdof*(3*i+2)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; + Ke->values[numdof*(3*i+2)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; + Ke->values[numdof*(3*i+2)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; } } @@ -8859,7 +9105,9 @@ /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list); /*Clean up and return*/ + xDelete(cs_list); + xDelete(dbasis); delete gauss; return Ke; @@ -8869,12 +9117,11 @@ void Penta::GetSolutionFromInputsDiagnosticHoriz(Vector* solution){ - const int numdof=NDOF2*NUMVERTICES; - - int i; int approximation; - int *doflist = NULL; + int *doflist = NULL; IssmDouble vx,vy; - IssmDouble values[numdof]; - GaussPenta *gauss; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; /*Get approximation enum and dof list: */ @@ -8883,15 +9130,14 @@ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - /*If the element is a coupling, do nothing: every node is also on an other elements - * (as coupling is between SSA and HO) so the other element will take care of it*/ + /*Fetch dof list and allocate solution vectors*/ GetDofList(&doflist,approximation,GsetEnum); + IssmDouble* values = xNew(numdof); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - gauss=new GaussPenta(); - for(i=0;iGaussNode(numnodes,i); /*Recover vx and vy*/ - gauss->GaussVertex(i); vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); @@ -8905,4 +9151,5 @@ /*Free ressources:*/ delete gauss; + xDelete(values); xDelete(doflist); } @@ -8980,15 +9227,23 @@ void Penta::GetSolutionFromInputsDiagnosticFS(Vector* solution){ - const int numdof=NDOF4*NUMVERTICES; - - int i; - int* doflist=NULL; - IssmDouble vx,vy,vz,p; - IssmDouble FSreconditioning; - IssmDouble values[numdof]; - GaussPenta *gauss; + int* vdoflist=NULL; + int* pdoflist=NULL; + IssmDouble vx,vy,vz,p; + IssmDouble FSreconditioning; + GaussPenta *gauss; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocityFinal(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF3; + int pnumdof = pnumnodes*NDOF1; + + /*Initialize values*/ + IssmDouble* vvalues = xNew(vnumdof); + IssmDouble* pvalues = xNew(pnumdof); /*Get dof list: */ - GetDofList(&doflist,FSApproximationEnum,GsetEnum); + GetDofListVelocity(&vdoflist,GsetEnum); + GetDofListPressure(&pdoflist,GsetEnum); Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); @@ -8996,28 +9251,35 @@ Input* p_input =inputs->GetInput(PressureEnum); _assert_(p_input); - /*Recondition pressure: */ this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum); - /*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */ - /*P1 element only for now*/ - gauss=new GaussPenta(); - for(i=0;iGaussVertex(i); + /*Ok, we have vx vy vz in values, fill in vx vy vz arrays: */ + gauss = new GaussPenta(); + for(int i=0;iGaussNode(vnumnodes,i); + vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); vz_input->GetInputValue(&vz,gauss); + vvalues[i*NDOF3+0]=vx; + vvalues[i*NDOF3+1]=vy; + vvalues[i*NDOF3+2]=vz; + } + for(int i=0;iGaussNode(pnumnodes,i); + p_input ->GetInputValue(&p ,gauss); - values[i*NDOF4+0]=vx; - values[i*NDOF4+1]=vy; - values[i*NDOF4+2]=vz; - values[i*NDOF4+3]=p/FSreconditioning; + pvalues[i]=p/FSreconditioning; } /*Add value to global vector*/ - solution->SetValues(numdof,doflist,values,INS_VAL); + solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL); + solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL); /*Free ressources:*/ delete gauss; - xDelete(doflist); + xDelete(pdoflist); + xDelete(vdoflist); + xDelete(pvalues); + xDelete(vvalues); } /*}}}*/ @@ -9343,5 +9605,5 @@ /*Transform solution in Cartesian Space*/ TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum); - TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum); + TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ @@ -9591,5 +9853,5 @@ /*Transform solution in Cartesian Space*/ TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum); - TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum); + TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ @@ -9814,36 +10076,47 @@ void Penta::InputUpdateFromSolutionDiagnosticFS(IssmDouble* solution){ - const int numdof=NDOF4*NUMVERTICES; - - int i; - IssmDouble values[numdof]; - IssmDouble vx[NUMVERTICES]; - IssmDouble vy[NUMVERTICES]; - IssmDouble vz[NUMVERTICES]; - IssmDouble vel[NUMVERTICES]; - IssmDouble pressure[NUMVERTICES]; - IssmDouble FSreconditioning; - int* doflist=NULL; + int i; + int* vdoflist=NULL; + int* pdoflist=NULL; + IssmDouble FSreconditioning; + GaussPenta *gauss; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocityFinal(); + int pnumnodes = this->NumberofNodesPressure(); + int vnumdof = vnumnodes*NDOF3; + int pnumdof = pnumnodes*NDOF1; + + /*Initialize values*/ + IssmDouble* vvalues = xNew(vnumdof); + IssmDouble* pvalues = xNew(pnumdof); + IssmDouble* vx = xNew(vnumnodes); + IssmDouble* vy = xNew(vnumnodes); + IssmDouble* vz = xNew(vnumnodes); + IssmDouble* vel = xNew(vnumnodes); + IssmDouble* pressure = xNew(vnumnodes); /*Get dof list: */ - GetDofList(&doflist,FSApproximationEnum,GsetEnum); + GetDofListVelocity(&vdoflist,GsetEnum); + GetDofListPressure(&pdoflist,GsetEnum); /*Use the dof list to index into the solution vector: */ - for(i=0;i(vx[i])) _error_("NaN found in solution vector"); - if(xIsNan(vy[i])) _error_("NaN found in solution vector"); - if(xIsNan(vz[i])) _error_("NaN found in solution vector"); + for(i=0;i(vx[i])) _error_("NaN found in solution vector"); + if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + if(xIsNan(vz[i])) _error_("NaN found in solution vector"); + } + for(i=0;i(pressure[i])) _error_("NaN found in solution vector"); } @@ -9851,6 +10124,6 @@ /*Recondition pressure and compute vel: */ this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum); - for(i=0;i(doflist); + xDelete(pressure); + xDelete(vel); + xDelete(vz); + xDelete(vy); + xDelete(vx); + xDelete(vvalues); + xDelete(pvalues); + xDelete(vdoflist); + xDelete(pdoflist); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15654) @@ -184,4 +184,6 @@ void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints); void GetDofList(int** pdoflist,int approximation_enum,int setenum); + void GetDofListVelocity(int** pdoflist,int setenum); + void GetDofListPressure(int** pdoflist,int setenum); void GetVertexPidList(int* doflist); void GetVertexSidList(int* sidlist); @@ -252,7 +254,6 @@ ElementMatrix* CreateKMatrixDiagnosticFS(void); ElementMatrix* CreateKMatrixDiagnosticFSViscous(void); - ElementMatrix* CreateKMatrixDiagnosticFSGLSViscous(void); + void KMatrixGLSstabilization(ElementMatrix* Ke); ElementMatrix* CreateKMatrixDiagnosticFSFriction(void); - ElementMatrix* CreateKMatrixDiagnosticFSGLSFriction(void); ElementMatrix* CreateKMatrixDiagnosticVert(void); ElementMatrix* CreateKMatrixDiagnosticVertVolume(void); @@ -295,5 +296,5 @@ ElementVector* CreatePVectorDiagnosticFSFront(void); ElementVector* CreatePVectorDiagnosticFSViscous(void); - ElementVector* CreatePVectorDiagnosticFSGLSViscous(void); + void PVectorGLSstabilization(ElementVector* pe); ElementVector* CreatePVectorDiagnosticFSShelf(void); ElementVector* CreatePVectorDiagnosticVert(void); Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15654) @@ -250,69 +250,134 @@ } /*}}}*/ +/*FUNCTION PentaRef::GetBFSstrainrate {{{*/ +void PentaRef::GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ + + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. + * For node i, Bi can be expressed in the actual coordinate system + * by: Bi=[ dh/dx 0 0 ] + * [ 0 dh/dy 0 ] + * [ 0 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx 0 ] + * [ 1/2*dh/dz 0 1/2*dh/dx ] + * [ 0 1/2*dh/dz 1/2*dh/dy ] + * where h is the interpolation function for node i. + * Same thing for Bb except the last column that does not exist. + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(3*numnodes); + GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B: */ + for(int i=0;i(dbasis); +} +/*}}}*/ /*FUNCTION PentaRef::GetBFS {{{*/ void PentaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. - * For node i, Bi can be expressed in the actual coordinate system - * by: Bi=[ dh/dx 0 0 0 ] - * [ 0 dh/dy 0 0 ] - * [ 0 0 dh/dy 0 ] - * [ 1/2*dh/dy 1/2*dh/dx 0 0 ] - * [ 1/2*dh/dz 0 1/2*dh/dx 0 ] - * [ 0 1/2*dh/dz 1/2*dh/dy 0 ] - * [ 0 0 0 h ] - * [ dh/dx dh/dy dh/dz 0 ] + /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. + * For node i, Bvi can be expressed in the actual coordinate system + * by: Bvi=[ dh/dx 0 0 ] + * [ 0 dh/dy 0 ] + * [ 0 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx 0 ] + * [ 1/2*dh/dz 0 1/2*dh/dx ] + * [ 0 1/2*dh/dz 1/2*dh/dy ] + * [ 0 0 0 ] + * [ dh/dx dh/dy dh/dz ] + * + * by: Bpi=[ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ h ] + * [ 0 ] * where h is the interpolation function for node i. * Same thing for Bb except the last column that does not exist. */ - int i; - - IssmDouble dbasismini[3][NUMNODESP1b]; - IssmDouble basis[NUMNODESP1]; - - /*Get dbasismini in actual coordinate system: */ - GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss); - GetNodalFunctionsP1(basis, gauss); + /*Fetch number of nodes for this finite element*/ + int pnumnodes = this->NumberofNodesPressure(); + int vnumnodes = this->NumberofNodesVelocity(); + + /*Get nodal functions derivatives*/ + IssmDouble* vdbasis=xNew(3*vnumnodes); + IssmDouble* pbasis =xNew(pnumnodes); + GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + GetNodalFunctionsPressure(pbasis,gauss); /*Build B: */ - for(i=0;i(vdbasis); + xDelete(pbasis); } /*}}}*/ @@ -388,64 +453,75 @@ * For node i, Bi' can be expressed in the actual coordinate system * by: - * Bi'=[ dh/dx 0 0 0] - * [ 0 dh/dy 0 0] - * [ 0 0 dh/dz 0] - * [ dh/dy dh/dx 0 0] - * [ dh/dz 0 dh/dx 0] - * [ 0 dh/dz dh/dy 0] - * [ dh/dx dh/dy dh/dz 0] - * [ 0 0 0 h] + * Bvi' = [ dh/dx 0 0 ] + * [ 0 dh/dy 0 ] + * [ 0 0 dh/dz ] + * [ dh/dy dh/dx 0 ] + * [ dh/dz 0 dh/dx ] + * [ 0 dh/dz dh/dy ] + * [ dh/dx dh/dy dh/dz ] + * [ 0 0 0 ] + * + * by: Bpi=[ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ 0 ] + * [ h ] * where h is the interpolation function for node i. - * - * Same thing for the bubble fonction except that there is no fourth column - */ - - int i; - IssmDouble dbasismini[3][NUMNODESP1b]; - IssmDouble basis[NUMNODESP1]; - - /*Get dbasismini in actual coordinate system: */ - GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss); - GetNodalFunctionsP1(basis, gauss); - - /*B_primeuild B_prime: */ - for(i=0;iNumberofNodesPressure(); + int vnumnodes = this->NumberofNodesVelocity(); + + /*Get nodal functions derivatives*/ + IssmDouble* vdbasis=xNew(3*vnumnodes); + IssmDouble* pbasis =xNew(pnumnodes); + GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + GetNodalFunctionsPressure(pbasis,gauss); + + /*Build B_prime: */ + for(int i=0;i(vdbasis); + xDelete(pbasis); } /*}}}*/ @@ -673,31 +749,37 @@ * For node i, Li can be expressed in the actual coordinate system * by: - * Li=[ h 0 ] - * [ 0 h ] - * [ 0 0 ] - * [ 0 0 ] + * Li=[ h 0 0 0 ] + * [ 0 h 0 0 ] * where h is the interpolation function for node i. */ - const int num_dof=4; - IssmDouble L1L2l3[NUMNODESP1_2d]; - - /*Get L1L2l3 in actual coordinate system: */ - L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; - L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; - L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; + /*Fetch number of nodes for this finite element*/ + int pnumnodes = this->NumberofNodesPressure(); + int vnumnodes = this->NumberofNodesVelocity(); + int pnumdof = pnumnodes; + int vnumdof = vnumnodes*NDOF3; + + /*Get nodal functions derivatives*/ + IssmDouble* vbasis=xNew(vnumnodes); + GetNodalFunctionsVelocity(vbasis,gauss); /*Build LFS: */ - for(int i=0;i<3;i++){ - LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i]; - 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.; - } + for(int i=0;i(vbasis); } /*}}}*/ @@ -1196,4 +1278,64 @@ } /*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsVelocity{{{*/ +void PentaRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussPenta* gauss){ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } +} +/*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsPressure{{{*/ +void PentaRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussPenta* gauss){ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = MINIEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } +} +/*}}}*/ /*FUNCTION PentaRef::GetNodalFunctionsDerivatives{{{*/ void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){ @@ -1229,4 +1371,60 @@ xDelete(dbasis_ref); +} +/*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsDerivativesVelocity{{{*/ +void PentaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){ + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1bubbleEnum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } +} +/*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsDerivativesPressure{{{*/ +void PentaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){ + switch(this->element_type){ + case P1P1Enum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1Enum; + return; + case P1P1GLSEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = P1P1GLSEnum; + return; + case MINIcondensedEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIcondensedEnum; + return; + case MINIEnum: + this->element_type = P1Enum; + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + this->element_type = MINIEnum; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } } /*}}}*/ @@ -1271,27 +1469,27 @@ case P1bubbleEnum: /*Nodal function 1*/ - dbasis[NUMNODESP1*0+0] = (zeta-1.)/4.; - dbasis[NUMNODESP1*1+0] = SQRT3/12.*(zeta-1.); - dbasis[NUMNODESP1*2+0] = -.5*gauss->coord1; + dbasis[NUMNODESP1b*0+0] = (zeta-1.)/4.; + dbasis[NUMNODESP1b*1+0] = SQRT3/12.*(zeta-1.); + dbasis[NUMNODESP1b*2+0] = -.5*gauss->coord1; /*Nodal function 2*/ - dbasis[NUMNODESP1*0+1] = (1.-zeta)/4.; - dbasis[NUMNODESP1*1+1] = SQRT3/12.*(zeta-1); - dbasis[NUMNODESP1*2+1] = -.5*gauss->coord2; + dbasis[NUMNODESP1b*0+1] = (1.-zeta)/4.; + dbasis[NUMNODESP1b*1+1] = SQRT3/12.*(zeta-1); + dbasis[NUMNODESP1b*2+1] = -.5*gauss->coord2; /*Nodal function 3*/ - dbasis[NUMNODESP1*0+2] = 0.; - dbasis[NUMNODESP1*1+2] = SQRT3/6.*(1.-zeta); - dbasis[NUMNODESP1*2+2] = -.5*gauss->coord3; + dbasis[NUMNODESP1b*0+2] = 0.; + dbasis[NUMNODESP1b*1+2] = SQRT3/6.*(1.-zeta); + dbasis[NUMNODESP1b*2+2] = -.5*gauss->coord3; /*Nodal function 4*/ - dbasis[NUMNODESP1*0+3] = -(1.+zeta)/4.; - dbasis[NUMNODESP1*1+3] = -SQRT3/12.*(1.+zeta); - dbasis[NUMNODESP1*2+3] = .5*gauss->coord1; + dbasis[NUMNODESP1b*0+3] = -(1.+zeta)/4.; + dbasis[NUMNODESP1b*1+3] = -SQRT3/12.*(1.+zeta); + dbasis[NUMNODESP1b*2+3] = .5*gauss->coord1; /*Nodal function 5*/ - dbasis[NUMNODESP1*0+4] = (1.+zeta)/4.; - dbasis[NUMNODESP1*1+4] = -SQRT3/12.*(1.+zeta); - dbasis[NUMNODESP1*2+4] = .5*gauss->coord2; + dbasis[NUMNODESP1b*0+4] = (1.+zeta)/4.; + dbasis[NUMNODESP1b*1+4] = -SQRT3/12.*(1.+zeta); + dbasis[NUMNODESP1b*2+4] = .5*gauss->coord2; /*Nodal function 6*/ - dbasis[NUMNODESP1*0+5] = 0.; - dbasis[NUMNODESP1*1+5] = SQRT3/6.*(1.+zeta); - dbasis[NUMNODESP1*2+5] = .5*gauss->coord3; + dbasis[NUMNODESP1b*0+5] = 0.; + dbasis[NUMNODESP1b*1+5] = SQRT3/6.*(1.+zeta); + dbasis[NUMNODESP1b*2+5] = .5*gauss->coord3; /*Nodal function 7*/ dbasis[NUMNODESP1b*0+6] = 27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3); @@ -1704,9 +1902,13 @@ switch(this->element_type){ - case P1Enum: return NUMNODESP1; - case P2Enum: return NUMNODESP2; - case P2xP1Enum: return NUMNODESP2xP1; - case P1xP2Enum: return NUMNODESP1xP2; - case MINIEnum: return NUMNODESP1b; + case P1Enum: return NUMNODESP1; + case P1bubbleEnum: return NUMNODESP1b; + case P2Enum: return NUMNODESP2; + case P2xP1Enum: return NUMNODESP2xP1; + case P1xP2Enum: return NUMNODESP1xP2; + case P1P1Enum: return NUMNODESP1*2; + case P1P1GLSEnum: return NUMNODESP1*2; + case MINIcondensedEnum: return NUMNODESP1*2; + case MINIEnum: return NUMNODESP1b+NUMNODESP1; default: _error_("Element type "<element_type)<<" not supported yet"); } @@ -1715,2 +1917,46 @@ } /*}}}*/ +/*FUNCTION PentaRef::NumberofNodesPressure{{{*/ +int PentaRef::NumberofNodesPressure(void){ + + switch(this->element_type){ + case P1P1Enum: return NUMNODESP1; + case P1P1GLSEnum: return NUMNODESP1; + case MINIcondensedEnum: return NUMNODESP1; + case MINIEnum: return NUMNODESP1; + default: _error_("Element type "<element_type)<<" not supported yet"); + } + + return -1; +} +/*}}}*/ +/*FUNCTION PentaRef::NumberofNodesVelocity{{{*/ +int PentaRef::NumberofNodesVelocity(void){ + + switch(this->element_type){ + case P1P1Enum: return NUMNODESP1; + case P1P1GLSEnum: return NUMNODESP1; + case MINIcondensedEnum: return NUMNODESP1b; + case MINIEnum: return NUMNODESP1b; + default: _error_("Element type "<element_type)<<" not supported yet"); + } + + return -1; +} +/*}}}*/ +/*FUNCTION PentaRef::NumberofNodesVelocityFinal{{{*/ +int PentaRef::NumberofNodesVelocityFinal(void){ + + /*When static condensation is applied, the final number of nodes might be different*/ + + switch(this->element_type){ + case P1P1Enum: return NUMNODESP1; + case P1P1GLSEnum: return NUMNODESP1; + case MINIcondensedEnum: return NUMNODESP1; + case MINIEnum: return NUMNODESP1b; + default: _error_("Element type "<element_type)<<" not supported yet"); + } + + return -1; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 15654) @@ -23,5 +23,9 @@ /*Numerics*/ void GetNodalFunctions(IssmDouble* basis, GaussPenta* gauss); + void GetNodalFunctionsVelocity(IssmDouble* basis, GaussPenta* gauss); + void GetNodalFunctionsPressure(IssmDouble* basis, GaussPenta* gauss); void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss); + void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss); + void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss); void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussPenta* gauss); void GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss); @@ -41,4 +45,5 @@ void GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); void GetBHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); + void GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); void GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); void GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss); @@ -65,4 +70,7 @@ int NumberofNodes(void); + int NumberofNodesVelocity(void); + int NumberofNodesVelocityFinal(void); + int NumberofNodesPressure(void); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15654) @@ -3353,8 +3353,6 @@ /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ gauss=new GaussTria(); for(int i=0;iGaussNode(numnodes,i); Index: /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp (revision 15654) @@ -110,2 +110,12 @@ } /*}}}*/ +/*FUNCTION DatasetInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){{{*/ +void DatasetInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){ + + /*Get requested input within dataset*/ + if(index<0 || index > inputs->Size()-1) _error_("index requested (" << index << ") exceeds dataset size (" << inputs->Size() << ")"); + Input* input=dynamic_cast(this->inputs->GetObjectByOffset(index)); + + input->GetInputValue(pvalue,gauss); +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h (revision 15654) @@ -50,5 +50,5 @@ void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index); - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index); void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp (revision 15654) @@ -134,143 +134,76 @@ /*FUNCTION PentaInput::GetVxStrainRate3d{{{*/ void PentaInput::GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){ - int i,j; - - const int numnodes=6; - const int DOFVELOCITY=3; - IssmDouble B[8][27]; - IssmDouble B_reduced[6][DOFVELOCITY*numnodes]; - IssmDouble velocity[numnodes][DOFVELOCITY]; - _assert_(this->NumberofNodes()==6); //Check Tria too + + /*Intermediary*/ + int numnodes=this->NumberofNodes(); + IssmDouble* B=xNew(numnodes*(NDOF3*numnodes)); + IssmDouble* velocity=xNew(numnodes*NDOF3); /*Get B matrix: */ - GetBFS(&B[0][0], xyz_list, gauss); - - /*Create a reduced matrix of B to get rid of pressure */ - for (i=0;i<6;i++){ - for (j=0;j<3;j++){ - B_reduced[i][j]=B[i][j]; - } - for (j=4;j<7;j++){ - B_reduced[i][j-1]=B[i][j]; - } - for (j=8;j<11;j++){ - B_reduced[i][j-2]=B[i][j]; - } - for (j=12;j<15;j++){ - B_reduced[i][j-3]=B[i][j]; - } - for (j=16;j<19;j++){ - B_reduced[i][j-4]=B[i][j]; - } - for (j=20;j<23;j++){ - B_reduced[i][j-5]=B[i][j]; - } - } + GetBFSstrainrate(B,xyz_list,gauss); /*Here, we are computing the strain rate of (vx,0,0)*/ - for(i=0;ivalues[i]; - velocity[i][1]=0.0; - velocity[i][2]=0.0; + for(int i=0;ivalues[i]; + velocity[NDOF3*i+1]=0.; + velocity[NDOF3*i+2]=0.; } /*Multiply B by velocity, to get strain rate: */ - MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvx,0); + MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0); + + /*Clean-up*/ + xDelete(B); + xDelete(velocity); } /*}}}*/ /*FUNCTION PentaInput::GetVyStrainRate3d{{{*/ -void PentaInput::GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){ - int i,j; - - const int numnodes=6; - const int DOFVELOCITY=3; - IssmDouble B[8][27]; - IssmDouble B_reduced[6][DOFVELOCITY*numnodes]; - IssmDouble velocity[numnodes][DOFVELOCITY]; - - _assert_(this->NumberofNodes()==6); //Check Tria too +void PentaInput::GetVyStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){ + + /*Intermediary*/ + int numnodes=this->NumberofNodes(); + IssmDouble* B=xNew(numnodes*(NDOF3*numnodes)); + IssmDouble* velocity=xNew(numnodes*NDOF3); /*Get B matrix: */ - GetBFS(&B[0][0], xyz_list, gauss); - /*Create a reduced matrix of B to get rid of pressure */ - for (i=0;i<6;i++){ - for (j=0;j<3;j++){ - B_reduced[i][j]=B[i][j]; - } - for (j=4;j<7;j++){ - B_reduced[i][j-1]=B[i][j]; - } - for (j=8;j<11;j++){ - B_reduced[i][j-2]=B[i][j]; - } - for (j=12;j<15;j++){ - B_reduced[i][j-3]=B[i][j]; - } - for (j=16;j<19;j++){ - B_reduced[i][j-4]=B[i][j]; - } - for (j=20;j<23;j++){ - B_reduced[i][j-5]=B[i][j]; - } - } + GetBFSstrainrate(B,xyz_list,gauss); /*Here, we are computing the strain rate of (0,vy,0)*/ - for(i=0;ivalues[i]; - velocity[i][2]=0.0; + for(int i=0;ivalues[i]; + velocity[NDOF3*i+2]=0.; } /*Multiply B by velocity, to get strain rate: */ - MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvy,0); - + MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0); + + /*Clean-up*/ + xDelete(B); + xDelete(velocity); } /*}}}*/ /*FUNCTION PentaInput::GetVzStrainRate3d{{{*/ -void PentaInput::GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss){ - int i,j; - - const int numnodes=6; - const int DOFVELOCITY=3; - IssmDouble B[8][27]; - IssmDouble B_reduced[6][DOFVELOCITY*numnodes]; - IssmDouble velocity[numnodes][DOFVELOCITY]; +void PentaInput::GetVzStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){ + + /*Intermediary*/ + int numnodes=this->NumberofNodes(); + IssmDouble* B=xNew(numnodes*(NDOF3*numnodes)); + IssmDouble* velocity=xNew(numnodes*NDOF3); /*Get B matrix: */ - GetBFS(&B[0][0], xyz_list, gauss); - - _assert_(this->NumberofNodes()==6); //Check Tria too - - /*Create a reduced matrix of B to get rid of pressure */ - for (i=0;i<6;i++){ - for (j=0;j<3;j++){ - B_reduced[i][j]=B[i][j]; - } - for (j=4;j<7;j++){ - B_reduced[i][j-1]=B[i][j]; - } - for (j=8;j<11;j++){ - B_reduced[i][j-2]=B[i][j]; - } - for (j=12;j<15;j++){ - B_reduced[i][j-3]=B[i][j]; - } - for (j=16;j<19;j++){ - B_reduced[i][j-4]=B[i][j]; - } - for (j=20;j<23;j++){ - B_reduced[i][j-5]=B[i][j]; - } - } + GetBFSstrainrate(B,xyz_list,gauss); /*Here, we are computing the strain rate of (0,0,vz)*/ - for(i=0;ivalues[i]; - } - + for(int i=0;ivalues[i]; + } /*Multiply B by velocity, to get strain rate: */ - MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvz,0); - + MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0); + + /*Clean-up*/ + xDelete(B); + xDelete(velocity); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 15654) @@ -475,5 +475,5 @@ /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZPEnum); + TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZEnum); /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15654) @@ -561,4 +561,39 @@ exz=epsilon[3]; eyz=epsilon[4]; + eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; + + mu_prime=(1-n)/(2*n) * mu/eff2; + } + + /*Assign output pointers:*/ + *pmu_prime=mu_prime; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosityDerivativeEpsSquareFS{{{*/ +void Matdamageice::GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* epsilon){ + + /*output: */ + IssmDouble mu_prime; + IssmDouble mu,n,eff2; + + /*input strain rate: */ + IssmDouble exx,eyy,ezz,exy,exz,eyz; + + /*Get visocisty and n*/ + GetViscosity3d(&mu,epsilon); + n=GetN(); + + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && + (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ + mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14); + } + else{ + /*Retrive strain rate components: */ + exx=epsilon[0]; + eyy=epsilon[1]; + ezz=epsilon[2]; + exy=epsilon[3]; + exz=epsilon[4]; + eyz=epsilon[5]; eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; Index: /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.h (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.h (revision 15654) @@ -56,4 +56,5 @@ void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon); void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); IssmDouble GetA(); Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 15654) @@ -32,4 +32,5 @@ virtual void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; virtual void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; + virtual void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; virtual void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; virtual IssmDouble GetA()=0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 15654) @@ -495,4 +495,39 @@ exz=epsilon[3]; eyz=epsilon[4]; + eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; + + mu_prime=(1-n)/(2*n) * mu/eff2; + } + + /*Assign output pointers:*/ + *pmu_prime=mu_prime; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosityDerivativeEpsSquareFS{{{*/ +void Matice::GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* epsilon){ + + /*output: */ + IssmDouble mu_prime; + IssmDouble mu,n,eff2; + + /*input strain rate: */ + IssmDouble exx,eyy,exy,exz,eyz,ezz; + + /*Get visocisty and n*/ + GetViscosity3d(&mu,epsilon); + n=GetN(); + + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && + (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ + mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14); + } + else{ + /*Retrive strain rate components: */ + exx=epsilon[0]; + eyy=epsilon[1]; + ezz=epsilon[2]; + exy=epsilon[3]; + exz=epsilon[4]; + eyz=epsilon[5]; eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 15654) @@ -63,4 +63,5 @@ void GetViscosityZComplement(IssmDouble*, IssmDouble*){_error_("not supported");}; void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon); void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); IssmDouble GetA(); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 15654) @@ -93,4 +93,5 @@ void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; IssmDouble GetA(){_error_("not supported");}; Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15654) @@ -67,5 +67,5 @@ /*Coordinate system provided, convert to coord_system matrix*/ _assert_(iomodel->Data(DiagnosticReferentialEnum)); - XZvectorsToCoordinateSystem(&this->coord_system[0][0],iomodel->Data(DiagnosticReferentialEnum)+io_index*6); + XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(DiagnosticReferentialEnum)[io_index*6]); if(iomodel->dim==3){ @@ -438,9 +438,4 @@ /*Put dof for this node into the s set (ie, this dof will be constrained * to a fixed value during computations. */ - - if(dof>=this->indexing.gsize){ - printf("dof spc = %i\n",dof); - this->Echo(); - } _assert_(dofindexing.gsize); @@ -1010,6 +1005,7 @@ for(i=0;iGetCoordinateSystem(&coord_system[0][0]); switch(cs_array[i]){ + case PressureEnum: + /*DO NOT change anything*/ + transform[(numdofs)*(counter) + counter] = 1.; + counter+=1; + break; case XYEnum: /*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/ @@ -1225,6 +1230,6 @@ counter+=2; break; - case XYZPEnum: - /*Only the first 3 coordinates are changed (x,y,z), leave the others (P) unchanged*/ + case XYZEnum: + /*The 3 coordinates are changed (x,y,z)*/ transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0]; transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1]; @@ -1236,6 +1241,5 @@ transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1]; transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2]; - transform[(numdofs)*(counter+3) + counter+3] = 1.0; - counter+=4; + counter+=3; break; default: Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 15654) @@ -448,4 +448,17 @@ } break; + case 7: //P1+ Lagrange element + switch(iv){ + case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break; + case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break; + case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break; + case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break; + case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break; + case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break; + + case 6: coord1=0.; coord2=0.; coord3=0.; coord4=0.;break; + default: _error_("node index should be in [0 5]"); + } + break; case 15: //P2 Lagrange element switch(iv){ @@ -470,5 +483,5 @@ } break; - default: _error_("supported number of nodes are 6 and 15"); + default: _error_("Number of nodes "<col_slocaldoflist=NULL; this->col_sglobaldoflist=NULL; - } /*}}}*/ @@ -423,4 +422,5 @@ _assert_(Ke); + _assert_(this); this->nrows =Ke->nrows; Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15654) @@ -11,22 +11,8 @@ /*intermediary: */ - int i,j,count; - IssmDouble yts; - FILE *fid = NULL; + FILE *fid = NULL; int code,vector_layout; - IssmDouble *times = NULL; - IssmDouble *values = NULL; - bool spcpresent = false; - - /*P2 finite elements*/ - int v1,v2; - bool *my_edges = NULL; - - /*variables being fetched: */ IssmDouble *spcdata = NULL; int M,N; - - /*Fetch parameters: */ - iomodel->Constant(&yts,ConstantsYtsEnum); /*First of, find the record for the enum, and get code of data type: */ @@ -38,7 +24,31 @@ iomodel->FetchData(&spcdata,&M,&N,vector_enum); + /*Call IoModelToConstraintsx*/ + IoModelToConstraintsx(constraints,iomodel,spcdata,M,N,analysis_type,finite_element,dof); + + /*Clean up*/ + xDelete(spcdata); +} + +void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,IssmDouble* spcdata,int M,int N,int analysis_type,int finite_element,int dof){ + + /*intermediary: */ + int i,j,count,elementnbv; + IssmDouble value; + IssmDouble *times = NULL; + IssmDouble *values = NULL; + bool spcpresent = false; + + /*P2 finite elements*/ + int v1,v2; + bool *my_edges = NULL; + switch(finite_element){ case P1Enum: /*Nothing else to do*/ + break; + case P1bubbleEnum: + if(iomodel->dim!=3) _error_("3d is the only supported dimension"); + elementnbv = 6; break; case P1xP2Enum: @@ -84,4 +94,26 @@ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); + count++; + } + } + } + break; + case P1bubbleEnum: + for(i=0;inumberofvertices;i++){ + if((iomodel->my_vertices[i])){ + if (!xIsNan(spcdata[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,dof,spcdata[i],analysis_type)); + count++; + } + } + } + for(i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + value = spcdata[iomodel->elements[i*elementnbv+0]-1]; + for(j=1;jelements[i*elementnbv+j]-1]; + value = value/reCast(elementnbv+0); + if(!xIsNan(value)){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, + dof,value,analysis_type)); count++; } @@ -145,5 +177,5 @@ /*figure out times: */ times=xNew(N); - for(j=0;j(spcdata); xDelete(times); xDelete(values); Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15653) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15654) @@ -9,4 +9,5 @@ /* local prototypes: */ void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element,int dof=1); +void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,IssmDouble* spcdata,int M,int N,int analysis_type,int finite_element,int dof=1); #endif /* _IOMODELTOELEMENTINPUTX_H */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 15654) @@ -97,10 +97,61 @@ break; - case MINIcondensedEnum: + /*Stokes elements*/ + case P1P1Enum: _assert_(approximation==FSApproximationEnum); /*P1 velocity*/ for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,approximation)); + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum)); + } + } + /*P1 pressure*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum)); + } + } + break; + case P1P1GLSEnum: + _assert_(approximation==FSApproximationEnum); + /*P1 velocity*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum)); + } + } + /*P1 pressure*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum)); + } + } + break; + case MINIcondensedEnum: + _assert_(approximation==FSApproximationEnum); + /*P1 velocity (bubble statically condensed)*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum)); + } + } + /*P1 pressure*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum)); + } + } + break; + case MINIEnum: + _assert_(approximation==FSApproximationEnum); + /*P1+ velocity*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum)); + } + } + for(i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis,FSvelocityEnum)); } } @@ -108,5 +159,5 @@ for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,i,i,iomodel,analysis,approximation)); + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,i,iomodel,analysis,FSpressureEnum)); } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15654) @@ -68,5 +68,5 @@ /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/ - if(!iscoupling && !isFS){ + if(!iscoupling){ /*Get finite element type*/ @@ -94,7 +94,60 @@ else if(isFS){ finiteelement = P1Enum; + iomodel->Constant(&temp,FlowequationFeFSEnum); + switch(temp){ + case 0 : finiteelement = P1Enum; break;//P1P1 + case 1 : finiteelement = P1Enum; break;//P1P1GSL + case 2 : finiteelement = P1Enum; break;//MINIcondensedEnum + case 3 : finiteelement = P1bubbleEnum; break;//MINIEnum + default: _error_("finite element "<FetchData(&spcvz,&Mz,&Nz,DiagnosticSpcvzEnum); + iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum); + iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + iomodel->FetchData(&nodeonicesheet,NULL,NULL,MaskVertexongroundediceEnum); + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if(reCast(nodeonbed[i]) && reCast(nodeonicesheet[i]) && reCast(nodeonFS[i])){ + if(vertices_type[i] == FSApproximationEnum){ + for(j=0;jDeleteData(spcvz,DiagnosticSpcvzEnum); + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum); + iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum); + iomodel->DeleteData(nodeonicesheet,MaskVertexongroundediceEnum); + + /*Pressure spc*/ + count = constraints->Size(); + iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); + iomodel->FetchData(&z,NULL,NULL,MeshZEnum); + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if (reCast(vertices_type[i])==NoneApproximationEnum){ + constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,DiagnosticHorizAnalysisEnum)); + count++; + } + } + } + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(surface,SurfaceEnum); + iomodel->DeleteData(z,MeshZEnum); + } *pconstraints=constraints; @@ -325,5 +378,5 @@ } if (reCast(vertices_type[i])==NoneApproximationEnum){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,g*rho_ice*(surface[i]-z[i])/FSreconditioning,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy count++; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15653) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15654) @@ -60,5 +60,12 @@ else if(isFS){ approximation = FSApproximationEnum; - finiteelement = P1Enum; + iomodel->Constant(&temp,FlowequationFeFSEnum); + switch(temp){ + case 0 : finiteelement = P1P1Enum; break; + case 1 : finiteelement = P1P1GLSEnum; break; + case 2 : finiteelement = MINIcondensedEnum; break; + case 3 : finiteelement = MINIEnum; break; + default: _error_("finite element "<Constant(&temp,FlowequationFeFSEnum); + switch(temp){ + case 0 : finiteelement = P1P1Enum; break; + case 1 : finiteelement = P1P1GLSEnum; break; + case 2 : finiteelement = MINIcondensedEnum; break; + case 3 : finiteelement = MINIEnum; break; + default: _error_("finite element "<