Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18927) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18928) @@ -11,4 +11,588 @@ /*Model processing*/ +void StressbalanceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + int i,j; + int count,finiteelement; + IssmDouble g; + IssmDouble rho_ice; + IssmDouble FSreconditioning; + bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling; + bool spcpresent = false; + int Mx,Nx; + int My,Ny; + int Mz,Nz; + IssmDouble *spcvx = NULL; + IssmDouble *spcvy = NULL; + IssmDouble *spcvz = NULL; + IssmDouble *nodeonSSA = NULL; + IssmDouble *nodeonHO = NULL; + IssmDouble *nodeonFS = NULL; + IssmDouble *nodeonbase = NULL; + IssmDouble *groundedice_ls = NULL; + IssmDouble *vertices_type = NULL; + IssmDouble *surface = NULL; + IssmDouble *z = NULL; + IssmDouble *timesx=NULL; + IssmDouble *timesy=NULL; + IssmDouble *timesz=NULL; + IssmDouble* values=NULL; + + /*Fetch parameters: */ + iomodel->Constant(&g,ConstantsGEnum); + iomodel->Constant(&rho_ice,MaterialsRhoIceEnum); + iomodel->Constant(&FSreconditioning,StressbalanceFSreconditioningEnum); + iomodel->Constant(&isSIA,FlowequationIsSIAEnum); + iomodel->Constant(&isSSA,FlowequationIsSSAEnum); + iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); + iomodel->Constant(&isHO,FlowequationIsHOEnum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + + /*Now, is the flag macayaealHO on? otherwise, do nothing: */ + if(!isSSA && !isHO && !isFS && !isL1L2) return; + + /*Do we have coupling*/ + if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + iscoupling = true; + else + iscoupling = false; + + /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/ + if(!iscoupling){ + + /*Get finite element type*/ + if(isSSA) iomodel->Constant(&finiteelement,FlowequationFeSSAEnum); + else if(isL1L2) finiteelement = P1Enum; + else if(isHO) iomodel->Constant(&finiteelement,FlowequationFeHOEnum); + else if(isFS){ iomodel->Constant(&finiteelement,FlowequationFeFSEnum); + /*Deduce velocity interpolation from finite element*/ + switch(finiteelement){ + case P1P1Enum : finiteelement = P1Enum; break; + case P1P1GLSEnum : finiteelement = P1Enum; break; + case MINIcondensedEnum : finiteelement = P1bubbleEnum; break; + case MINIEnum : finiteelement = P1bubbleEnum; break; + case TaylorHoodEnum : finiteelement = P2Enum; break; + case XTaylorHoodEnum : finiteelement = P2Enum; break; + case LATaylorHoodEnum : finiteelement = P2Enum; break; + case LACrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break; + case OneLayerP4zEnum : finiteelement = P2xP4Enum; break; + case CrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break; + default: _error_("finite element "<FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum); + iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum); + iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum); + if(iomodel->domaintype==Domain3DEnum){ + iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum); + } + else if (iomodel->domaintype==Domain2DverticalEnum){ + iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvyEnum); + } + else{ + _error_("not supported yet"); + } + if(iomodel->domaintype==Domain3DEnum){ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0); + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1); + IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2); + iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum); + } + else if (iomodel->domaintype==Domain2DverticalEnum){ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0); + IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,1); + iomodel->DeleteData(spcvz,StressbalanceSpcvyEnum); + } + else{ + _error_("not supported yet"); + } + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum); + iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum); + iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum); + + /*Pressure spc*/ + count = constraints->Size(); + iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); + iomodel->FetchData(&z,NULL,NULL,MeshZEnum); + switch(finiteelement){ + case P1Enum: + 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,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; + case P1bubbleEnum: + 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,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; + case P2Enum: + 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->numberofedges+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; + case P2bubbleEnum: + 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->numberofedges+iomodel->numberofelements+iomodel->numberoffaces+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; + case P2xP4Enum: + //Nothing yet + break; + default: + _error_("not implemented yet"); + } + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(surface,SurfaceEnum); + iomodel->DeleteData(z,MeshZEnum); + } + else{ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0); + if(iomodel->domaintype!=Domain2DverticalEnum){ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1); + } + } + + return; + } + + /*Constraints: fetch data: */ + iomodel->FetchData(&spcvx,&Mx,&Nx,StressbalanceSpcvxEnum); + iomodel->FetchData(&spcvy,&My,&Ny,StressbalanceSpcvyEnum); + iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum); + iomodel->FetchData(&nodeonSSA,NULL,NULL,FlowequationBorderSSAEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonHO,NULL,NULL,FlowequationBorderHOEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum); + iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); + iomodel->FetchData(&z,NULL,NULL,MeshZEnum); + + /*Initialize counter: */ + count=0; + + /*figure out times: */ + timesx=xNew(Nx); + for(j=0;j(Ny); + for(j=0;j(Nz); + for(j=0;jnumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + + /*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/ + if(reCast(vertices_type[i]==SSAHOApproximationEnum)){ + /*If grionSSA, spc HO dofs: 3 & 4*/ + if (reCast(nodeonHO[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + if (!xIsNan(spcvx[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + + } + else if (reCast(nodeonSSA[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + if (!xIsNan(spcvx[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + + } + else _error_("if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA"); + } + /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/ + else if (reCast(vertices_type[i])==HOFSApproximationEnum){ + /*If grion,HO spc FS dofs: 3 4 & 5*/ + if (reCast(nodeonHO[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + if (!xIsNan(spcvx[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + + } + else if (reCast(nodeonFS[i])){ //spc HO nodes: 1 & 2 + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + if (!xIsNan(spcvx[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvz[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + } + else _error_("if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS"); + } + /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/ + else if (reCast(vertices_type[i])==SSAFSApproximationEnum){ + /*If grion,HO spc FS dofs: 3 4 & 5*/ + if (reCast(nodeonSSA[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + if (!xIsNan(spcvx[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + + } + else if (reCast(nodeonFS[i])){ //spc SSA nodes: 1 & 2 + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + if (!xIsNan(spcvx[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvz[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + } + else _error_("if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS"); + } + /*Now add the regular spcs*/ + else{ + if (Mx==iomodel->numberofvertices && !xIsNan(spcvx[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + + } + else if (Mx==iomodel->numberofvertices+1) { + /*figure out times and values: */ + values=xNew(Nx); + spcpresent=false; + for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default + } + + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,Nx,timesx,values,StressbalanceAnalysisEnum)); + count++; + } + xDelete(values); + } + else if (vertices_type[i]==SIAApproximationEnum){ + constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,StressbalanceAnalysisEnum)); + count++; + } + + if (My==iomodel->numberofvertices && !xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy. + count++; + } + else if (My==iomodel->numberofvertices+1){ + /*figure out times and values: */ + values=xNew(Ny); + spcpresent=false; + for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default + } + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Ny,timesy,values,StressbalanceAnalysisEnum)); + count++; + } + xDelete(values); + } + else if (vertices_type[i]==SIAApproximationEnum){ + constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum)); + count++; + } + + if (reCast(vertices_type[i])==FSApproximationEnum || (reCast(vertices_type[i])==NoneApproximationEnum)){ + if (Mz==iomodel->numberofvertices && !xIsNan(spcvz[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy + count++; + } + else if (Mz==iomodel->numberofvertices+1){ + /*figure out times and values: */ + values=xNew(Nz); + spcpresent=false; + for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default + } + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Nz,timesz,values,StressbalanceAnalysisEnum)); + count++; + } + xDelete(values); + } + + } + if (reCast(vertices_type[i])==NoneApproximationEnum){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy + count++; + } + } + } + } + + /*Free data: */ + iomodel->DeleteData(spcvx,StressbalanceSpcvxEnum); + iomodel->DeleteData(spcvy,StressbalanceSpcvyEnum); + iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum); + iomodel->DeleteData(nodeonSSA,FlowequationBorderSSAEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonHO,FlowequationBorderHOEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum); + if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum); + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(surface,SurfaceEnum); + iomodel->DeleteData(z,MeshZEnum); + + /*Free resources:*/ + xDelete(timesx); + xDelete(timesy); + xDelete(timesz); + xDelete(values); + +}/*}}}*/ +void StressbalanceAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + const int RIFTINFOSIZE = 12; + int i; + int count; + int penpair_ids[2]; + bool isSSA,isL1L2,isHO,isFS; + int numpenalties,numrifts,numriftsegments; + IssmDouble *riftinfo = NULL; + IssmDouble *penalties = NULL; + int assert_int; + + /*Fetch parameters: */ + iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + iomodel->Constant(&isSSA,FlowequationIsSSAEnum); + iomodel->Constant(&isHO,FlowequationIsHOEnum); + iomodel->Constant(&numrifts,RiftsNumriftsEnum); + + /*Now, is the flag macayaealHO on? otherwise, do nothing: */ + if(!isSSA && !isHO && !isFS && !isL1L2) return; + + /*Initialize counter: */ + count=0; + + /*Create Penpair for penalties: */ + iomodel->FetchData(&penalties,&numpenalties,NULL,StressbalanceVertexPairingEnum); + + for(i=0;imy_vertices[reCast(penalties[2*i+0]-1)]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + assert_int=iomodel->my_vertices[reCast(penalties[2*i+1]-1)]; _assert_(assert_int); + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(penalties[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(penalties[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],StressbalanceAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(penalties,StressbalanceVertexPairingEnum); + + /*Create Riffront loads for rifts: */ + if(numrifts){ + iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum); + iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum); + for(i=0;imy_elements[reCast(*(riftinfo+RIFTINFOSIZE*i+2))-1]){ + loads->AddObject(new Riftfront(iomodel->loadcounter+count+1,i,iomodel,StressbalanceAnalysisEnum)); + count++; + } + } + iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum); + xDelete(riftinfo); + } +}/*}}}*/ +void StressbalanceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + bool isSSA,isL1L2,isHO,isFS,iscoupling; + int finiteelement=-1,approximation=-1; + + /*Fetch parameters: */ + iomodel->Constant(&isSSA,FlowequationIsSSAEnum); + iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); + iomodel->Constant(&isHO,FlowequationIsHOEnum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + + /*Now, check that we have non SIA elements */ + if(!isSSA & !isL1L2 & !isHO & !isFS) return; + + /*Do we have coupling*/ + if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + iscoupling = true; + else + iscoupling = false; + + /*If no coupling, call Regular CreateNodes, else, use P1 elements only*/ + if(!iscoupling){ + + /*Get finite element type*/ + if(isSSA){ + approximation=SSAApproximationEnum; + iomodel->Constant(&finiteelement,FlowequationFeSSAEnum); + } + else if(isL1L2){ + approximation = L1L2ApproximationEnum; + finiteelement = P1Enum; + } + else if(isHO){ + approximation = HOApproximationEnum; + iomodel->Constant(&finiteelement,FlowequationFeHOEnum); + } + else if(isFS){ + approximation = FSApproximationEnum; + iomodel->Constant(&finiteelement,FlowequationFeFSEnum); + } + iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(3,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum); + ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation); + iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + } + else{ + /*Coupling: we are going to create P1 Elements only*/ + + Node* node = NULL; + int lid=0; + if(!nodes) nodes = new Nodes(); + + iomodel->FetchData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + if(isFS){ + /*P1+ velocity*/ + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + approximation=reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]); + if(approximation==FSApproximationEnum) approximation=FSvelocityEnum; + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation)); + } + } + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum); + node->Deactivate(); + nodes->AddObject(node); + } + } + /*P1 pressure*/ + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + approximation=reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]); + node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum); + if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){ + node->Deactivate(); + } + nodes->AddObject(node); + } + } + } + else{ + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]))); + } + } + } + iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + } +}/*}}}*/ int StressbalanceAnalysis::DofsPerNode(int** pdoftype,int domaintype,int approximation){/*{{{*/ @@ -84,56 +668,4 @@ return numdofs; }/*}}}*/ -void StressbalanceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ - - /*Intermediaries*/ - int fe_FS; - int numoutputs; - char** requestedoutputs = NULL; - int materials_type; - - parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum)); - parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum)); - parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum)); - parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum)); - parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum)); - parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum)); - parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum)); - - /*XTH LATH parameters*/ - iomodel->Constant(&fe_FS,FlowequationFeFSEnum); - if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){ - parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum)); - parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum)); - parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum)); - } - - iomodel->Constant(&materials_type,MaterialsEnum); - if(materials_type==MatdamageiceEnum){ - parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum)); - parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum)); - } - - /*Requested outputs*/ - iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum); - parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs)); - if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs)); - iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum); - - /*Deal with friction parameters*/ - int frictionlaw; - iomodel->Constant(&frictionlaw,FrictionLawEnum); - if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum)); - -}/*}}}*/ void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -304,587 +836,55 @@ xDelete(finiteelement_list); }/*}}}*/ -void StressbalanceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - - /*Intermediary*/ - bool isSSA,isL1L2,isHO,isFS,iscoupling; - int finiteelement=-1,approximation=-1; - - /*Fetch parameters: */ - iomodel->Constant(&isSSA,FlowequationIsSSAEnum); - iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); - iomodel->Constant(&isHO,FlowequationIsHOEnum); - iomodel->Constant(&isFS,FlowequationIsFSEnum); - - /*Now, check that we have non SIA elements */ - if(!isSSA & !isL1L2 & !isHO & !isFS) return; - - /*Do we have coupling*/ - if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) - iscoupling = true; - else - iscoupling = false; - - /*If no coupling, call Regular CreateNodes, else, use P1 elements only*/ - if(!iscoupling){ - - /*Get finite element type*/ - if(isSSA){ - approximation=SSAApproximationEnum; - iomodel->Constant(&finiteelement,FlowequationFeSSAEnum); - } - else if(isL1L2){ - approximation = L1L2ApproximationEnum; - finiteelement = P1Enum; - } - else if(isHO){ - approximation = HOApproximationEnum; - iomodel->Constant(&finiteelement,FlowequationFeHOEnum); - } - else if(isFS){ - approximation = FSApproximationEnum; - iomodel->Constant(&finiteelement,FlowequationFeFSEnum); - } - iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum); - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(3,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum); - ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation); - iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, - FlowequationVertexEquationEnum,StressbalanceReferentialEnum); - } - else{ - /*Coupling: we are going to create P1 Elements only*/ - - Node* node = NULL; - int lid=0; - if(!nodes) nodes = new Nodes(); - - iomodel->FetchData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, - FlowequationVertexEquationEnum,StressbalanceReferentialEnum); - if(isFS){ - /*P1+ velocity*/ - for(int i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - approximation=reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]); - if(approximation==FSApproximationEnum) approximation=FSvelocityEnum; - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation)); - } - } - for(int i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum); - node->Deactivate(); - nodes->AddObject(node); - } - } - /*P1 pressure*/ - for(int i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - approximation=reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]); - node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum); - if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){ - node->Deactivate(); - } - nodes->AddObject(node); - } - } - } - else{ - for(int i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]))); - } - } - } - iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, - FlowequationVertexEquationEnum,StressbalanceReferentialEnum); - } -}/*}}}*/ -void StressbalanceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ - - /*Intermediary*/ - int i,j; - int count,finiteelement; - IssmDouble g; - IssmDouble rho_ice; - IssmDouble FSreconditioning; - bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling; - bool spcpresent = false; - int Mx,Nx; - int My,Ny; - int Mz,Nz; - IssmDouble *spcvx = NULL; - IssmDouble *spcvy = NULL; - IssmDouble *spcvz = NULL; - IssmDouble *nodeonSSA = NULL; - IssmDouble *nodeonHO = NULL; - IssmDouble *nodeonFS = NULL; - IssmDouble *nodeonbase = NULL; - IssmDouble *groundedice_ls = NULL; - IssmDouble *vertices_type = NULL; - IssmDouble *surface = NULL; - IssmDouble *z = NULL; - IssmDouble *timesx=NULL; - IssmDouble *timesy=NULL; - IssmDouble *timesz=NULL; - IssmDouble* values=NULL; - - /*Fetch parameters: */ - iomodel->Constant(&g,ConstantsGEnum); - iomodel->Constant(&rho_ice,MaterialsRhoIceEnum); - iomodel->Constant(&FSreconditioning,StressbalanceFSreconditioningEnum); - iomodel->Constant(&isSIA,FlowequationIsSIAEnum); - iomodel->Constant(&isSSA,FlowequationIsSSAEnum); - iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); - iomodel->Constant(&isHO,FlowequationIsHOEnum); - iomodel->Constant(&isFS,FlowequationIsFSEnum); - - /*Now, is the flag macayaealHO on? otherwise, do nothing: */ - if(!isSSA && !isHO && !isFS && !isL1L2) return; - - /*Do we have coupling*/ - if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) - iscoupling = true; - else - iscoupling = false; - - /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/ - if(!iscoupling){ - - /*Get finite element type*/ - if(isSSA) iomodel->Constant(&finiteelement,FlowequationFeSSAEnum); - else if(isL1L2) finiteelement = P1Enum; - else if(isHO) iomodel->Constant(&finiteelement,FlowequationFeHOEnum); - else if(isFS){ iomodel->Constant(&finiteelement,FlowequationFeFSEnum); - /*Deduce velocity interpolation from finite element*/ - switch(finiteelement){ - case P1P1Enum : finiteelement = P1Enum; break; - case P1P1GLSEnum : finiteelement = P1Enum; break; - case MINIcondensedEnum : finiteelement = P1bubbleEnum; break; - case MINIEnum : finiteelement = P1bubbleEnum; break; - case TaylorHoodEnum : finiteelement = P2Enum; break; - case XTaylorHoodEnum : finiteelement = P2Enum; break; - case LATaylorHoodEnum : finiteelement = P2Enum; break; - case LACrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break; - case OneLayerP4zEnum : finiteelement = P2xP4Enum; break; - case CrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break; - default: _error_("finite element "<FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); - iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum); - iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum); - iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum); - if(iomodel->domaintype==Domain3DEnum){ - iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum); - } - else if (iomodel->domaintype==Domain2DverticalEnum){ - iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvyEnum); - } - else{ - _error_("not supported yet"); - } - if(iomodel->domaintype==Domain3DEnum){ - IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0); - IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1); - IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2); - iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum); - } - else if (iomodel->domaintype==Domain2DverticalEnum){ - IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0); - IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,1); - iomodel->DeleteData(spcvz,StressbalanceSpcvyEnum); - } - else{ - _error_("not supported yet"); - } - iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); - iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum); - iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum); - iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum); - - /*Pressure spc*/ - count = constraints->Size(); - iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); - iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); - iomodel->FetchData(&z,NULL,NULL,MeshZEnum); - switch(finiteelement){ - case P1Enum: - 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,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); - count++; - } - } - } - break; - case P1bubbleEnum: - 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,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); - count++; - } - } - } - break; - case P2Enum: - 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->numberofedges+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); - count++; - } - } - } - break; - case P2bubbleEnum: - 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->numberofedges+iomodel->numberofelements+iomodel->numberoffaces+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); - count++; - } - } - } - break; - case P2xP4Enum: - //Nothing yet - break; - default: - _error_("not implemented yet"); - } - iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); - iomodel->DeleteData(surface,SurfaceEnum); - iomodel->DeleteData(z,MeshZEnum); - } - else{ - IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0); - if(iomodel->domaintype!=Domain2DverticalEnum){ - IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1); - } - } - - return; - } - - /*Constraints: fetch data: */ - iomodel->FetchData(&spcvx,&Mx,&Nx,StressbalanceSpcvxEnum); - iomodel->FetchData(&spcvy,&My,&Ny,StressbalanceSpcvyEnum); - iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum); - iomodel->FetchData(&nodeonSSA,NULL,NULL,FlowequationBorderSSAEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonHO,NULL,NULL,FlowequationBorderHOEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum); - iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); - iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); - iomodel->FetchData(&z,NULL,NULL,MeshZEnum); - - /*Initialize counter: */ - count=0; - - /*figure out times: */ - timesx=xNew(Nx); - for(j=0;j(Ny); - for(j=0;j(Nz); - for(j=0;jnumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - - /*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/ - if(reCast(vertices_type[i]==SSAHOApproximationEnum)){ - /*If grionSSA, spc HO dofs: 3 & 4*/ - if (reCast(nodeonHO[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - - } - else if (reCast(nodeonSSA[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - - } - else _error_("if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA"); - } - /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/ - else if (reCast(vertices_type[i])==HOFSApproximationEnum){ - /*If grion,HO spc FS dofs: 3 4 & 5*/ - if (reCast(nodeonHO[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - - } - else if (reCast(nodeonFS[i])){ //spc HO nodes: 1 & 2 - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvz[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - } - else _error_("if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS"); - } - /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/ - else if (reCast(vertices_type[i])==SSAFSApproximationEnum){ - /*If grion,HO spc FS dofs: 3 4 & 5*/ - if (reCast(nodeonSSA[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - - } - else if (reCast(nodeonFS[i])){ //spc SSA nodes: 1 & 2 - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvz[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - } - else _error_("if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS"); - } - /*Now add the regular spcs*/ - else{ - if (Mx==iomodel->numberofvertices && !xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - - } - else if (Mx==iomodel->numberofvertices+1) { - /*figure out times and values: */ - values=xNew(Nx); - spcpresent=false; - for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default - } - - if(spcpresent){ - constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,Nx,timesx,values,StressbalanceAnalysisEnum)); - count++; - } - xDelete(values); - } - else if (vertices_type[i]==SIAApproximationEnum){ - constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,StressbalanceAnalysisEnum)); - count++; - } - - if (My==iomodel->numberofvertices && !xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy. - count++; - } - else if (My==iomodel->numberofvertices+1){ - /*figure out times and values: */ - values=xNew(Ny); - spcpresent=false; - for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default - } - if(spcpresent){ - constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Ny,timesy,values,StressbalanceAnalysisEnum)); - count++; - } - xDelete(values); - } - else if (vertices_type[i]==SIAApproximationEnum){ - constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum)); - count++; - } - - if (reCast(vertices_type[i])==FSApproximationEnum || (reCast(vertices_type[i])==NoneApproximationEnum)){ - if (Mz==iomodel->numberofvertices && !xIsNan(spcvz[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy - count++; - } - else if (Mz==iomodel->numberofvertices+1){ - /*figure out times and values: */ - values=xNew(Nz); - spcpresent=false; - for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default - } - if(spcpresent){ - constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Nz,timesz,values,StressbalanceAnalysisEnum)); - count++; - } - xDelete(values); - } - - } - if (reCast(vertices_type[i])==NoneApproximationEnum){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy - count++; - } - } - } - } - - /*Free data: */ - iomodel->DeleteData(spcvx,StressbalanceSpcvxEnum); - iomodel->DeleteData(spcvy,StressbalanceSpcvyEnum); - iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum); - iomodel->DeleteData(nodeonSSA,FlowequationBorderSSAEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonHO,FlowequationBorderHOEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum); - if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum); - iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); - iomodel->DeleteData(surface,SurfaceEnum); - iomodel->DeleteData(z,MeshZEnum); - - /*Free resources:*/ - xDelete(timesx); - xDelete(timesy); - xDelete(timesz); - xDelete(values); - -}/*}}}*/ -void StressbalanceAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - - /*Intermediary*/ - const int RIFTINFOSIZE = 12; - int i; - int count; - int penpair_ids[2]; - bool isSSA,isL1L2,isHO,isFS; - int numpenalties,numrifts,numriftsegments; - IssmDouble *riftinfo = NULL; - IssmDouble *penalties = NULL; - int assert_int; - - /*Fetch parameters: */ - iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); - iomodel->Constant(&isFS,FlowequationIsFSEnum); - iomodel->Constant(&isSSA,FlowequationIsSSAEnum); - iomodel->Constant(&isHO,FlowequationIsHOEnum); - iomodel->Constant(&numrifts,RiftsNumriftsEnum); - - /*Now, is the flag macayaealHO on? otherwise, do nothing: */ - if(!isSSA && !isHO && !isFS && !isL1L2) return; - - /*Initialize counter: */ - count=0; - - /*Create Penpair for penalties: */ - iomodel->FetchData(&penalties,&numpenalties,NULL,StressbalanceVertexPairingEnum); - - for(i=0;imy_vertices[reCast(penalties[2*i+0]-1)]){ - - /*In debugging mode, check that the second node is in the same cpu*/ - assert_int=iomodel->my_vertices[reCast(penalties[2*i+1]-1)]; _assert_(assert_int); - - /*Get node ids*/ - penpair_ids[0]=iomodel->nodecounter+reCast(penalties[2*i+0]); - penpair_ids[1]=iomodel->nodecounter+reCast(penalties[2*i+1]); - - /*Create Load*/ - loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],StressbalanceAnalysisEnum)); - count++; - } - } - - /*free ressources: */ - iomodel->DeleteData(penalties,StressbalanceVertexPairingEnum); - - /*Create Riffront loads for rifts: */ - if(numrifts){ - iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum); - iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum); - for(i=0;imy_elements[reCast(*(riftinfo+RIFTINFOSIZE*i+2))-1]){ - loads->AddObject(new Riftfront(iomodel->loadcounter+count+1,i,iomodel,StressbalanceAnalysisEnum)); - count++; - } - } - iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum); - xDelete(riftinfo); - } +void StressbalanceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*Intermediaries*/ + int fe_FS; + int numoutputs; + char** requestedoutputs = NULL; + int materials_type; + + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum)); + + /*XTH LATH parameters*/ + iomodel->Constant(&fe_FS,FlowequationFeFSEnum); + if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){ + parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum)); + parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum)); + parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum)); + } + + iomodel->Constant(&materials_type,MaterialsEnum); + if(materials_type==MatdamageiceEnum){ + parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum)); + parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum)); + } + + /*Requested outputs*/ + iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum); + parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum); + + /*Deal with friction parameters*/ + int frictionlaw; + iomodel->Constant(&frictionlaw,FrictionLawEnum); + if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum)); + }/*}}}*/ @@ -1027,5 +1027,5 @@ } }/*}}}*/ -void StressbalanceAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ int approximation; @@ -1048,5 +1048,5 @@ } }/*}}}*/ -void StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector* solution,Element* element){/*{{{*/ IssmDouble vx,vy; @@ -1098,8 +1098,8 @@ xDelete(doflist); }/*}}}*/ -void StressbalanceAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ +void StressbalanceAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int approximation; @@ -1133,5 +1133,5 @@ } }/*}}}*/ -void StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ +void StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ bool islevelset; @@ -1660,5 +1660,5 @@ return pe; }/*}}}*/ -void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system @@ -1700,5 +1700,5 @@ xDelete(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. * For node i, Bi can be expressed in the actual coordinate system @@ -1737,5 +1737,5 @@ xDelete(basis); }/*}}}*/ -void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system @@ -1777,5 +1777,5 @@ xDelete(dbasis); }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/ int i,dim,domaintype; @@ -2111,5 +2111,5 @@ return pe; }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/ int i,dim,domaintype; @@ -2286,4 +2286,84 @@ return Ke; }/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFriction(Element* element){/*{{{*/ + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + if(element->IsFloating() || !element->IsOnBase()) return NULL; + + /*Intermediaries*/ + int dim; + bool mainlyfloating; + int migration_style,point1; + IssmDouble alpha2,Jdet,fraction1,fraction2; + IssmDouble gllevelset,phi=1.; + IssmDouble *xyz_list_base = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + int numdof = numnodes*(dim-1); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum); + IssmDouble* B = xNew((dim-1)*numdof); + IssmDouble* D = xNewZeroInit((dim-1)*(dim-1)); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->FindParam(&migration_style,GroundinglineMigrationEnum); + Input* gllevelset_input = NULL; + + /*build friction object, used later on: */ + Friction* friction=new Friction(element,2); + + /*Recover portion of element that is grounded*/ + if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base); + if(migration_style==SubelementMigration2Enum){ + gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); + element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + //gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); + gauss=element->NewGaussBase(2); + } + else{ + gauss=element->NewGaussBase(2); + } + + /* Start looping on the number of gaussian points: */ + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + friction->GetAlpha2(&alpha2,gauss); + if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2; + if(migration_style==SubelementMigration2Enum){ + gllevelset_input->GetInputValue(&gllevelset, gauss); + if(gllevelset<0.) alpha2=0.; + } + + this->GetBHOFriction(B,element,dim,xyz_list_base,gauss); + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + for(int i=0;iweight*Jdet; + + TripleMultiply(B,dim-1,numdof,1, + D,dim-1,dim-1,0, + B,dim-1,numdof,0, + &Ke->values[0],1); + } + + /*Transform Coordinate System*/ + if(dim==3) element->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + delete gauss; + delete friction; + xDelete(xyz_list_base); + xDelete(B); + xDelete(D); + return Ke; +}/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOViscous(Element* element){/*{{{*/ @@ -2358,84 +2438,4 @@ return Ke; }/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFriction(Element* element){/*{{{*/ - - /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - - if(element->IsFloating() || !element->IsOnBase()) return NULL; - - /*Intermediaries*/ - int dim; - bool mainlyfloating; - int migration_style,point1; - IssmDouble alpha2,Jdet,fraction1,fraction2; - IssmDouble gllevelset,phi=1.; - IssmDouble *xyz_list_base = NULL; - Gauss* gauss = NULL; - - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - int numdof = numnodes*(dim-1); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum); - IssmDouble* B = xNew((dim-1)*numdof); - IssmDouble* D = xNewZeroInit((dim-1)*(dim-1)); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - element->FindParam(&migration_style,GroundinglineMigrationEnum); - Input* gllevelset_input = NULL; - - /*build friction object, used later on: */ - Friction* friction=new Friction(element,2); - - /*Recover portion of element that is grounded*/ - if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base); - if(migration_style==SubelementMigration2Enum){ - gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); - element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - //gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); - gauss=element->NewGaussBase(2); - } - else{ - gauss=element->NewGaussBase(2); - } - - /* Start looping on the number of gaussian points: */ - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - friction->GetAlpha2(&alpha2,gauss); - if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2; - if(migration_style==SubelementMigration2Enum){ - gllevelset_input->GetInputValue(&gllevelset, gauss); - if(gllevelset<0.) alpha2=0.; - } - - this->GetBHOFriction(B,element,dim,xyz_list_base,gauss); - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - for(int i=0;iweight*Jdet; - - TripleMultiply(B,dim-1,numdof,1, - D,dim-1,dim-1,0, - B,dim-1,numdof,0, - &Ke->values[0],1); - } - - /*Transform Coordinate System*/ - if(dim==3) element->TransformStiffnessMatrixCoord(Ke,XYEnum); - - /*Clean up and return*/ - delete gauss; - delete friction; - xDelete(xyz_list_base); - xDelete(B); - xDelete(D); - return Ke; -}/*}}}*/ #ifdef FSANALYTICAL ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/ @@ -2633,5 +2633,5 @@ return pe; }/*}}}*/ -void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system @@ -2682,5 +2682,42 @@ xDelete(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. + * For node i, Bi can be expressed in the actual coordinate system + * by: + * 3D 2D + * Bi=[ N 0 ] Bi=N + * [ 0 N ] + * where N is the finiteelement function for node i. + * + * We assume B has been allocated already, of size: 2 x (numdof*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build L: */ + if(dim==3){ + for(int i=0;i(basis); +}/*}}}*/ +void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system @@ -2727,42 +2764,5 @@ xDelete(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. - * For node i, Bi can be expressed in the actual coordinate system - * by: - * 3D 2D - * Bi=[ N 0 ] Bi=N - * [ 0 N ] - * where N is the finiteelement function for node i. - * - * We assume B has been allocated already, of size: 2 x (numdof*numnodes) - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions derivatives*/ - IssmDouble* basis=xNew(numnodes); - element->NodalFunctions(basis,gauss); - - /*Build L: */ - if(dim==3){ - for(int i=0;i(basis); -}/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/ int i,dim; @@ -2982,4 +2982,141 @@ return Ke; }/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/ + + if(!element->IsFloating() || !element->IsOnBase()) return NULL; + + /*If on not water or not FS, skip stiffness: */ + int approximation,shelf_dampening; + element->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; + element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum); + if(shelf_dampening==0) return NULL; + + /*Intermediaries*/ + bool mainlyfloating; + int j,i,dim; + IssmDouble Jdet,slope2,scalar,dt; + IssmDouble slope[3]; + IssmDouble *xyz_list_base = NULL; + IssmDouble *xyz_list = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int numdof = vnumnodes*dim + pnumnodes; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + if(dt==0) dt=1.e+5; + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + + /* Start looping on the number of gaussian points: */ + gauss=element->NewGaussBase(3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + if(dim==2) slope2=slope[0]*slope[0]; + else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1]; + scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; + for(i=0;ivalues[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j]; + } + } + } + + /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list_base); + xDelete(xyz_list); + xDelete(vbasis); + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/ + + /*Intermediaries*/ + int i,dim,epssize; + IssmDouble viscosity,FSreconditioning,Jdet; + IssmDouble *xyz_list = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + if(dim==2) epssize = 3; + else epssize = 6; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int numdof = vnumnodes*dim + pnumnodes; + int bsize = epssize + 2; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementMatrix(FSvelocityEnum); + IssmDouble* B = xNew(bsize*numdof); + IssmDouble* Bprime = xNew(bsize*numdof); + IssmDouble* D = xNewZeroInit(bsize*bsize); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input; + if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + this->GetBFS(B,element,dim,xyz_list,gauss); + this->GetBFSprime(Bprime,element,dim,xyz_list,gauss); + + element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + for(i=0;iweight*Jdet; + for(i=epssize;iweight*Jdet; + + TripleMultiply(B,bsize,numdof,1, + D,bsize,bsize,0, + Bprime,bsize,numdof,0, + &Ke->values[0],1); + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(D); + xDelete(Bprime); + xDelete(B); + xDelete(cs_list); + return Ke; +}/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLA(Element* element){/*{{{*/ @@ -3181,141 +3318,4 @@ return Ke; }/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/ - - /*Intermediaries*/ - int i,dim,epssize; - IssmDouble viscosity,FSreconditioning,Jdet; - IssmDouble *xyz_list = NULL; - - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - if(dim==2) epssize = 3; - else epssize = 6; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - int numdof = vnumnodes*dim + pnumnodes; - int bsize = epssize + 2; - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - if(dim==2) for(i=0;iNewElementMatrix(FSvelocityEnum); - IssmDouble* B = xNew(bsize*numdof); - IssmDouble* Bprime = xNew(bsize*numdof); - IssmDouble* D = xNewZeroInit(bsize*bsize); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input; - if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} - - /* Start looping on the number of gaussian points: */ - Gauss* gauss = element->NewGauss(5); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - this->GetBFS(B,element,dim,xyz_list,gauss); - this->GetBFSprime(Bprime,element,dim,xyz_list,gauss); - - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - for(i=0;iweight*Jdet; - for(i=epssize;iweight*Jdet; - - TripleMultiply(B,bsize,numdof,1, - D,bsize,bsize,0, - Bprime,bsize,numdof,0, - &Ke->values[0],1); - } - - /*Transform Coordinate System*/ - element->TransformStiffnessMatrixCoord(Ke,cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list); - xDelete(D); - xDelete(Bprime); - xDelete(B); - xDelete(cs_list); - return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/ - - if(!element->IsFloating() || !element->IsOnBase()) return NULL; - - /*If on not water or not FS, skip stiffness: */ - int approximation,shelf_dampening; - element->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum); - if(shelf_dampening==0) return NULL; - - /*Intermediaries*/ - bool mainlyfloating; - int j,i,dim; - IssmDouble Jdet,slope2,scalar,dt; - IssmDouble slope[3]; - IssmDouble *xyz_list_base = NULL; - IssmDouble *xyz_list = NULL; - Gauss* gauss = NULL; - - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - int numdof = vnumnodes*dim + pnumnodes; - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&dt,TimesteppingTimeStepEnum); - if(dt==0) dt=1.e+5; - IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); - IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); - Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); - - /* Start looping on the number of gaussian points: */ - gauss=element->NewGaussBase(3); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); - element->NodalFunctionsVelocity(vbasis,gauss); - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - if(dim==2) slope2=slope[0]*slope[0]; - else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1]; - scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; - for(i=0;ivalues[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j]; - } - } - } - - /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list_base); - xDelete(xyz_list); - xDelete(vbasis); - return Ke; -}/*}}}*/ #ifdef FSANALYTICAL ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/ @@ -3402,4 +3402,114 @@ /*clean-up and return*/ + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/ + + if(!element->IsOnBase()) return NULL; + + /*Intermediaries*/ + int dim; + IssmDouble alpha2,Jdet; + IssmDouble bed_normal[3]; + IssmDouble *xyz_list_base = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + + /*Initialize Element matrix and vectors*/ + ElementVector* pe = element->NewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + Input* alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input); + + /* Start looping on the number of gaussian points: */ + gauss=element->NewGaussBase(3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + alpha2_input->GetInputValue(&alpha2, gauss); + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + element->NormalBase(&bed_normal[0],xyz_list_base); + + for(int i=0;ivalues[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1]; + pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0]; + if(dim==3){ + pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i]; + } + } + + } + + /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list_base); + xDelete(vbasis); + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/ + + if(!element->IsOnBase()) return NULL; + + /*Intermediaries*/ + int dim; + IssmDouble sigmann,sigmant,Jdet,bedslope,beta; + IssmDouble *xyz_list_base = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + + /*Initialize Element matrix and vectors*/ + ElementVector* pe = element->NewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + Input* sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input); + Input* sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input); + Input* bedslope_input=element->GetInput(BedSlopeXEnum); _assert_(bedslope_input); + + /* Start looping on the number of gaussian points: */ + gauss=element->NewGaussBase(3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + sigmann_input->GetInputValue(&sigmann, gauss); + sigmant_input->GetInputValue(&sigmant, gauss); + bedslope_input->GetInputValue(&bedslope, gauss); + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + + beta=sqrt(1+bedslope*bedslope); + for(int i=0;ivalues[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i]; + pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i]; + if(dim==3){ + //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i]; + _error_("3d not supported yet"); + } + } + + } + + /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list_base); + xDelete(vbasis); return pe; }/*}}}*/ @@ -3477,114 +3587,4 @@ return pe3; } - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/ - - if(!element->IsOnBase()) return NULL; - - /*Intermediaries*/ - int dim; - IssmDouble alpha2,Jdet; - IssmDouble bed_normal[3]; - IssmDouble *xyz_list_base = NULL; - Gauss* gauss = NULL; - - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - - /*Initialize Element matrix and vectors*/ - ElementVector* pe = element->NewElementVector(FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - Input* alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input); - - /* Start looping on the number of gaussian points: */ - gauss=element->NewGaussBase(3); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - alpha2_input->GetInputValue(&alpha2, gauss); - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctionsVelocity(vbasis,gauss); - element->NormalBase(&bed_normal[0],xyz_list_base); - - for(int i=0;ivalues[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1]; - pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0]; - if(dim==3){ - pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i]; - } - } - - } - - /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list_base); - xDelete(vbasis); - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/ - - if(!element->IsOnBase()) return NULL; - - /*Intermediaries*/ - int dim; - IssmDouble sigmann,sigmant,Jdet,bedslope,beta; - IssmDouble *xyz_list_base = NULL; - Gauss* gauss = NULL; - - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - - /*Initialize Element matrix and vectors*/ - ElementVector* pe = element->NewElementVector(FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - Input* sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input); - Input* sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input); - Input* bedslope_input=element->GetInput(BedSlopeXEnum); _assert_(bedslope_input); - - /* Start looping on the number of gaussian points: */ - gauss=element->NewGaussBase(3); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - sigmann_input->GetInputValue(&sigmann, gauss); - sigmant_input->GetInputValue(&sigmant, gauss); - bedslope_input->GetInputValue(&bedslope, gauss); - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctionsVelocity(vbasis,gauss); - - beta=sqrt(1+bedslope*bedslope); - for(int i=0;ivalues[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i]; - pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i]; - if(dim==3){ - //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i]; - _error_("3d not supported yet"); - } - } - - } - - /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list_base); - xDelete(vbasis); return pe; }/*}}}*/ @@ -3790,4 +3790,167 @@ }/*}}}*/ #endif +ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/ + + /*If no front, return NULL*/ + if(!element->IsIcefront()) return NULL; + + /*Intermediaries*/ + int i,dim; + IssmDouble Jdet,pressure,surface,z; + IssmDouble normal[3]; + IssmDouble *xyz_list = NULL; + IssmDouble *xyz_list_front = NULL; + Gauss *gauss = NULL; + + /*Make sure current element is floating*/ + if(!element->IsFloating()) return NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int numvertices = element->GetNumberOfVertices(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + element->NormalSection(&normal[0],xyz_list_front); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + + /*Initialize gauss points*/ + IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;izmax) zmax=xyz_list[i*3+(dim-1)]; + IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity + else gauss=element->NewGauss(xyz_list,xyz_list_front,3,3); + + /* Start looping on the number of gaussian points: */ + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + surface_input->GetInputValue(&surface,gauss); + if(dim==3) z=element->GetZcoord(xyz_list,gauss); + else z=element->GetYcoord(xyz_list,gauss); + pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level + + for (int i=0;ivalues[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i]; + pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i]; + if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + xDelete(xyz_list); + xDelete(xyz_list_front); + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/ + + int i,dim; + IssmDouble Jdet,water_pressure,bed; + IssmDouble normal[3]; + IssmDouble *xyz_list_base = NULL; + + /*Get basal element*/ + if(!element->IsOnBase() || !element->IsFloating()) return NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + Input* base_input=element->GetInput(BaseEnum); _assert_(base_input); + IssmDouble rho_water=element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + + element->NormalBase(&normal[0],xyz_list_base); + _assert_(normal[dim-1]<0.); + base_input->GetInputValue(&bed, gauss); + water_pressure=gravity*rho_water*bed; + + for(i=0;ivalues[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0]; + pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1]; + if(dim==3){ + pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2]; + } + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,cs_list); + + /* shelf dampening*/ + int shelf_dampening; + element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum); + if(shelf_dampening) { + Input* mb_input=element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input); + IssmDouble dt,mb,normal_b; + element->FindParam(&dt,TimesteppingTimeStepEnum); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + element->NormalBase(&normal[0],xyz_list_base); + if (dim==2) normal_b=normal[1]; + else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]); + mb_input->GetInputValue(&mb, gauss); + for(i=0;ivalues[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b; + } + } + } + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + xDelete(xyz_list_base); + return pe; +}/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLA(Element* element){/*{{{*/ @@ -4044,168 +4207,5 @@ return pe; }/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/ - - int i,dim; - IssmDouble Jdet,water_pressure,bed; - IssmDouble normal[3]; - IssmDouble *xyz_list_base = NULL; - - /*Get basal element*/ - if(!element->IsOnBase() || !element->IsFloating()) return NULL; - - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - Input* base_input=element->GetInput(BaseEnum); _assert_(base_input); - IssmDouble rho_water=element->GetMaterialParameter(MaterialsRhoSeawaterEnum); - IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGaussBase(5); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctionsVelocity(vbasis,gauss); - - element->NormalBase(&normal[0],xyz_list_base); - _assert_(normal[dim-1]<0.); - base_input->GetInputValue(&bed, gauss); - water_pressure=gravity*rho_water*bed; - - for(i=0;ivalues[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0]; - pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1]; - if(dim==3){ - pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2]; - } - } - } - - /*Transform coordinate system*/ - element->TransformLoadVectorCoord(pe,cs_list); - - /* shelf dampening*/ - int shelf_dampening; - element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum); - if(shelf_dampening) { - Input* mb_input=element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input); - IssmDouble dt,mb,normal_b; - element->FindParam(&dt,TimesteppingTimeStepEnum); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctionsVelocity(vbasis,gauss); - element->NormalBase(&normal[0],xyz_list_base); - if (dim==2) normal_b=normal[1]; - else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]); - mb_input->GetInputValue(&mb, gauss); - for(i=0;ivalues[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b; - } - } - } - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(vbasis); - xDelete(xyz_list_base); - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/ - - /*If no front, return NULL*/ - if(!element->IsIcefront()) return NULL; - - /*Intermediaries*/ - int i,dim; - IssmDouble Jdet,pressure,surface,z; - IssmDouble normal[3]; - IssmDouble *xyz_list = NULL; - IssmDouble *xyz_list_front = NULL; - Gauss *gauss = NULL; - - /*Make sure current element is floating*/ - if(!element->IsFloating()) return NULL; - - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - int numvertices = element->GetNumberOfVertices(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); - IssmDouble* vbasis = xNew(vnumnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); - element->NormalSection(&normal[0],xyz_list_front); - Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); - IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); - IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); - - /*Initialize gauss points*/ - IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;izmax) zmax=xyz_list[i*3+(dim-1)]; - IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity - else gauss=element->NewGauss(xyz_list,xyz_list_front,3,3); - - /* Start looping on the number of gaussian points: */ - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); - element->NodalFunctionsVelocity(vbasis,gauss); - surface_input->GetInputValue(&surface,gauss); - if(dim==3) z=element->GetZcoord(xyz_list,gauss); - else z=element->GetYcoord(xyz_list,gauss); - pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level - - for (int i=0;ivalues[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i]; - pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i]; - if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i]; - } - } - - /*Transform coordinate system*/ - element->TransformLoadVectorCoord(pe,cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(vbasis); - xDelete(xyz_list); - xDelete(xyz_list_front); - return pe; -}/*}}}*/ -void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*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 @@ -4319,5 +4319,65 @@ xDelete(pbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * For node i, Li can be expressed in the actual coordinate system + * by in 3d + * Li=[ h 0 0 0 ] + * [ 0 h 0 0 ] + * in 2d: + * Li=[ h 0 0 ] + * where h is the interpolation function for node i. + */ + + /*Fetch number of nodes for this finite element*/ + int pnumnodes = element->NumberofNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumdof = pnumnodes; + int vnumdof = vnumnodes*dim; + + /*Get nodal functions derivatives*/ + IssmDouble* vbasis=xNew(vnumnodes); + element->NodalFunctionsVelocity(vbasis,gauss); + + /*Build B: */ + if(dim==3){ + for(int i=0;i(vbasis); +}/*}}}*/ +void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system @@ -4432,5 +4492,126 @@ xDelete(pbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Bvi' = [ dphi/dx dphi/dy ] + * + * In 3d + * Bvi=[ dh/dx dh/dy dh/dz ] + * where phi is the finiteelement function for node i. + */ + + /*Fetch number of nodes for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + + /*Get nodal functions derivatives*/ + IssmDouble* vdbasis=xNew(dim*vnumnodes); + element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + + /*Build B_prime: */ + if(dim==2){ + for(int i=0;i(vdbasis); +}/*}}}*/ +void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Bvi' = [ dphi/dx 0 ] + * [ 0 dphi/dy ] + * [ dphi/dy dphi/dx ] + * + * In 3d + * 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 ] + * where phi is the finiteelement function for node i. + * In 3d: + */ + + /*Fetch number of nodes for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + + /*Get nodal functions derivatives*/ + IssmDouble* vdbasis=xNew(dim*vnumnodes); + element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + + /*Build B_prime: */ + if(dim==2){ + for(int i=0;i(vdbasis); +}/*}}}*/ +void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. + */ + + /*Fetch number of nodes for this finite element*/ + int pnumnodes; + if(dim==2) pnumnodes=3; + else pnumnodes=6; + //int pnumnodes = element->NumberofNodes(P1Enum); + + /*Get nodal functions derivatives*/ + IssmDouble* basis =xNew(pnumnodes); + element->NodalFunctionsP1(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*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 @@ -4496,186 +4677,5 @@ xDelete(vdbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. - * For node i, Bi' can be expressed in the actual coordinate system - * by: - * Bvi' = [ dphi/dx 0 ] - * [ 0 dphi/dy ] - * [ dphi/dy dphi/dx ] - * - * In 3d - * 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 ] - * where phi is the finiteelement function for node i. - * In 3d: - */ - - /*Fetch number of nodes for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - - /*Get nodal functions derivatives*/ - IssmDouble* vdbasis=xNew(dim*vnumnodes); - element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); - - /*Build B_prime: */ - if(dim==2){ - for(int i=0;i(vdbasis); -}/*}}}*/ -void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. - */ - - /*Fetch number of nodes for this finite element*/ - int pnumnodes; - if(dim==2) pnumnodes=3; - else pnumnodes=6; - //int pnumnodes = element->NumberofNodes(P1Enum); - - /*Get nodal functions derivatives*/ - IssmDouble* basis =xNew(pnumnodes); - element->NodalFunctionsP1(basis,gauss); - - /*Build B: */ - for(int i=0;i(basis); -}/*}}}*/ -void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. - * For node i, Bi' can be expressed in the actual coordinate system - * by: - * Bvi' = [ dphi/dx dphi/dy ] - * - * In 3d - * Bvi=[ dh/dx dh/dy dh/dz ] - * where phi is the finiteelement function for node i. - */ - - /*Fetch number of nodes for this finite element*/ - int vnumnodes = element->NumberofNodesVelocity(); - - /*Get nodal functions derivatives*/ - IssmDouble* vdbasis=xNew(dim*vnumnodes); - element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); - - /*Build B_prime: */ - if(dim==2){ - for(int i=0;i(vdbasis); -}/*}}}*/ -void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. - * For node i, Li can be expressed in the actual coordinate system - * by in 3d - * Li=[ h 0 0 0 ] - * [ 0 h 0 0 ] - * in 2d: - * Li=[ h 0 0 ] - * where h is the interpolation function for node i. - */ - - /*Fetch number of nodes for this finite element*/ - int pnumnodes = element->NumberofNodesPressure(); - int vnumnodes = element->NumberofNodesVelocity(); - int pnumdof = pnumnodes; - int vnumdof = vnumnodes*dim; - - /*Get nodal functions derivatives*/ - IssmDouble* vbasis=xNew(vnumnodes); - element->NodalFunctionsVelocity(vbasis,gauss); - - /*Build B: */ - if(dim==3){ - for(int i=0;i(vbasis); -}/*}}}*/ -void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute C matrix. C=[Cp1 Cp2 ...] where: * Cpi=[phi phi]. @@ -4698,5 +4698,5 @@ xDelete(basis); }/*}}}*/ -void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /* Compute C' matrix. C'=[C1' C2' ...] * Ci' = [ phi 0 ] @@ -4746,5 +4746,5 @@ xDelete(vbasis); }/*}}}*/ -void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector* solution,Element* element){/*{{{*/ int* vdoflist=NULL; @@ -4809,5 +4809,5 @@ xDelete(vvalues); }/*}}}*/ -void StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/ +void StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/ /*Intermediaries*/ @@ -4888,5 +4888,5 @@ }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/ bool results_on_nodes; @@ -4978,5 +4978,5 @@ xDelete(cs_list); }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/ /*Intermediaries*/ @@ -5191,5 +5191,5 @@ } }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/ /*Intermediaries*/ @@ -5327,163 +5327,4 @@ /*Coupling (Tiling)*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3d(Element* element){/*{{{*/ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixSSA3dViscous(element); - ElementMatrix* Ke2=CreateKMatrixSSA3dFriction(element); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dFriction(Element* element){/*{{{*/ - - /*Initialize Element matrix and return if necessary*/ - if(element->IsFloating() || !element->IsOnBase()) return NULL; - - /*Build a tria element using the 3 nodes of the base of the penta. Then use - * the tria functionality to build a friction stiffness matrix on these 3 - * nodes: */ - Element* basalelement = element->SpawnBasalElement(); - ElementMatrix* Ke=CreateKMatrixSSAFriction(basalelement); - basalelement->DeleteMaterials(); delete basalelement; - - /*clean-up and return*/ - return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/ - - /*Constants*/ - const int numdof2d=2*3; - - /*Intermediaries */ - int i,j,approximation; - int dim=3; - IssmDouble Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; - IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble epsilons[6]; //6 for FS - IssmDouble B[3][numdof2d]; - IssmDouble Bprime[3][numdof2d]; - IssmDouble D[3][3]= {0.0}; // material matrix, simple scalar matrix. - IssmDouble D_scalar; - IssmDouble Ke_gg[numdof2d][numdof2d]={0.0}; - IssmDouble *xyz_list = NULL; - - /*Find penta on bed as this is a SSA elements: */ - Element* pentabase=element->GetBasalElement(); - Element* basaltria=pentabase->SpawnBasalElement(); - - /*Initialize Element matrix*/ - ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum); - element->GetInputValue(&approximation,ApproximationEnum); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); - Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); - Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input); - Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); - Input* vz_input =element->GetInput(VzEnum); _assert_(vz_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(5); - Gauss* gauss_tria=new GaussTria(); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - gauss->SynchronizeGaussBase(gauss_tria); - - element->JacobianDeterminant(&Jdet, xyz_list,gauss); - this->GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria); - this->GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria); - - if(approximation==SSAHOApproximationEnum){ - element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - } - else if (approximation==SSAFSApproximationEnum){ - element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - } - else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet"); - - D_scalar=2*newviscosity*gauss->weight*Jdet; - for (i=0;i<3;i++) D[i][i]=D_scalar; - - TripleMultiply( &B[0][0],3,numdof2d,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof2d,0, - &Ke_gg[0][0],1); - - } - for(i=0;ivalues[i*numdof2d+j]+=Ke_gg[i][j]; - - /*Transform Coordinate System*/ - basaltria->TransformStiffnessMatrixCoord(Ke,XYEnum); - - /*Clean up and return*/ - xDelete(xyz_list); - delete basaltria->material; - delete basaltria; - delete gauss_tria; - delete gauss; - return Ke; - -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixFS(element); - int indices[3]={18,19,20}; - Ke1->StaticCondensation(3,&indices[0]); - int init = element->FiniteElement(); - element->SetTemporaryElementType(P1Enum); // P1 needed for HO - ElementMatrix* Ke2=CreateKMatrixHO(element); - element->SetTemporaryElementType(init); // P1 needed for HO - ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - delete Ke3; - return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixSSA3d(element); - ElementMatrix* Ke2=CreateKMatrixHO(element); - ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - delete Ke3; - return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixFS(element); - int indices[3]={18,19,20}; - Ke1->StaticCondensation(3,&indices[0]); - int init = element->FiniteElement(); - element->SetTemporaryElementType(P1Enum); - ElementMatrix* Ke2=CreateKMatrixSSA3d(element); - element->SetTemporaryElementType(init); - ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - delete Ke3; - return Ke; -}/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingHOFS(Element* element){/*{{{*/ @@ -5551,182 +5392,4 @@ delete Ke2; return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHO(Element* element){/*{{{*/ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(element); - ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(element); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/ - - if(element->IsFloating() || !element->IsOnBase()) return NULL; - - /*Constants*/ - int numnodes = element->GetNumberOfNodes(); - int numdof = 2*numnodes; - int numdoftotal = 4*numnodes; - - /*Intermediaries */ - int i,j; - IssmDouble Jdet2d,alpha2; - IssmDouble *xyz_list_tria = NULL; - IssmDouble* L = xNewZeroInit(2*numdof); - IssmDouble DL[2][2] = {{ 0,0 },{0,0}}; //for basal drag - IssmDouble DL_scalar; - IssmDouble* Ke_gg = xNewZeroInit(numdof*numdof); - Node **node_list = xNew(2*numnodes); - int* cs_list = xNew(2*numnodes); - - /*Initialize Element matrix and return if necessary*/ - ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum); - ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); - delete Ke1; delete Ke2; - - /*Prepare node list*/ - for(i=0;iGetNode(i); - node_list[i+1*numnodes] = element->GetNode(i); - cs_list[i+0*numnodes] = XYEnum; - cs_list[i+1*numnodes] = XYEnum; - } - - /*retrieve inputs :*/ - element->GetVerticesCoordinatesBase(&xyz_list_tria); - - /*build friction object, used later on: */ - Friction* friction=new Friction(element,2); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGaussBase(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /*Friction: */ - friction->GetAlpha2(&alpha2,gauss); - element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss); - this->GetBHOFriction(L,element,3,xyz_list_tria,gauss); - - DL_scalar=alpha2*gauss->weight*Jdet2d; - for (i=0;i<2;i++) DL[i][i]=DL_scalar; - - /* Do the triple producte tL*D*L: */ - TripleMultiply( L,2,numdof,1, - &DL[0][0],2,2,0, - L,2,numdof,0, - Ke_gg,1); - } - - for(i=0;ivalues[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j]; - for(i=0;ivalues[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j]; - - /*Transform Coordinate System*/ - element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list); - - /*Clean up and return*/ - delete gauss; - delete friction; - xDelete(cs_list); - xDelete(node_list); - xDelete(xyz_list_tria); - xDelete(Ke_gg); - xDelete(L); - return Ke; -}/*}}}*/ -ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/ - - /*Constants*/ - int numnodes = element->GetNumberOfNodes(); - int numdofm = 1 *numnodes; //*2/2 - int numdofp = 2 *numnodes; - int numdoftotal = 2 *2 *numnodes;//2 dof per nodes and 2 sets of nodes for HO and SSA - - /*Intermediaries */ - int i,j; - IssmDouble Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity - IssmDouble *xyz_list = NULL; - IssmDouble* B = xNew(3*numdofp); - IssmDouble* Bprime = xNew(3*numdofm); - IssmDouble D[3][3]={0.0}; // material matrix, simple scalar matrix. - IssmDouble D_scalar; - IssmDouble* Ke_gg = xNewZeroInit(numdofp*numdofm); - Node **node_list = xNew(2*numnodes); - int* cs_list= xNew(2*numnodes); - - /*Find penta on bed as HO must be coupled to the dofs on the bed: */ - Element* pentabase=element->GetBasalElement(); - Element* basaltria=pentabase->SpawnBasalElement(); - - /*prepare node list*/ - for(i=0;iGetNode(i); - node_list[i+1*numnodes] = element ->GetNode(i); - cs_list[i+0*numnodes] = XYEnum; - cs_list[i+1*numnodes] = XYEnum; - } - - /*Initialize Element matrix*/ - ElementMatrix* Ke1= pentabase->NewElementMatrix(SSAApproximationEnum); - ElementMatrix* Ke2= element ->NewElementMatrix(HOApproximationEnum); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - delete Ke1; delete Ke2; - - /* Get node coordinates and dof list: */ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); - Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); - Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input); - Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(5); - Gauss* gauss_tria=new GaussTria(); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - gauss->SynchronizeGaussBase(gauss_tria); - - element->JacobianDeterminant(&Jdet, xyz_list,gauss); - this->GetBSSAHO(B, element,xyz_list, gauss); - this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); - element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); - element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input); - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=2*newviscosity*gauss->weight*Jdet; - for (i=0;i<3;i++) D[i][i]=D_scalar; - - TripleMultiply( B,3,numdofp,1, - &D[0][0],3,3,0, - Bprime,3,numdofm,0, - Ke_gg,1); - } - for(i=0;ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j]; - for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i]; - - /*Transform Coordinate System*/ - element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list); - - /*Clean-up and return*/ - basaltria->DeleteMaterials(); delete basaltria; - - delete gauss; - delete gauss_tria; - xDelete(B); - xDelete(Bprime); - xDelete(Ke_gg); - xDelete(xyz_list); - xDelete(node_list); - xDelete(cs_list); - return Ke; - }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFS(Element* element){/*{{{*/ @@ -5970,9 +5633,346 @@ }/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHO(Element* element){/*{{{*/ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(element); + ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(element); + ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/ + + if(element->IsFloating() || !element->IsOnBase()) return NULL; + + /*Constants*/ + int numnodes = element->GetNumberOfNodes(); + int numdof = 2*numnodes; + int numdoftotal = 4*numnodes; + + /*Intermediaries */ + int i,j; + IssmDouble Jdet2d,alpha2; + IssmDouble *xyz_list_tria = NULL; + IssmDouble* L = xNewZeroInit(2*numdof); + IssmDouble DL[2][2] = {{ 0,0 },{0,0}}; //for basal drag + IssmDouble DL_scalar; + IssmDouble* Ke_gg = xNewZeroInit(numdof*numdof); + Node **node_list = xNew(2*numnodes); + int* cs_list = xNew(2*numnodes); + + /*Initialize Element matrix and return if necessary*/ + ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum); + ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum); + ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); + delete Ke1; delete Ke2; + + /*Prepare node list*/ + for(i=0;iGetNode(i); + node_list[i+1*numnodes] = element->GetNode(i); + cs_list[i+0*numnodes] = XYEnum; + cs_list[i+1*numnodes] = XYEnum; + } + + /*retrieve inputs :*/ + element->GetVerticesCoordinatesBase(&xyz_list_tria); + + /*build friction object, used later on: */ + Friction* friction=new Friction(element,2); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Friction: */ + friction->GetAlpha2(&alpha2,gauss); + element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss); + this->GetBHOFriction(L,element,3,xyz_list_tria,gauss); + + DL_scalar=alpha2*gauss->weight*Jdet2d; + for (i=0;i<2;i++) DL[i][i]=DL_scalar; + + /* Do the triple producte tL*D*L: */ + TripleMultiply( L,2,numdof,1, + &DL[0][0],2,2,0, + L,2,numdof,0, + Ke_gg,1); + } + + for(i=0;ivalues[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j]; + for(i=0;ivalues[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j]; + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list); + + /*Clean up and return*/ + delete gauss; + delete friction; + xDelete(cs_list); + xDelete(node_list); + xDelete(xyz_list_tria); + xDelete(Ke_gg); + xDelete(L); + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/ + + /*Constants*/ + int numnodes = element->GetNumberOfNodes(); + int numdofm = 1 *numnodes; //*2/2 + int numdofp = 2 *numnodes; + int numdoftotal = 2 *2 *numnodes;//2 dof per nodes and 2 sets of nodes for HO and SSA + + /*Intermediaries */ + int i,j; + IssmDouble Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity + IssmDouble *xyz_list = NULL; + IssmDouble* B = xNew(3*numdofp); + IssmDouble* Bprime = xNew(3*numdofm); + IssmDouble D[3][3]={0.0}; // material matrix, simple scalar matrix. + IssmDouble D_scalar; + IssmDouble* Ke_gg = xNewZeroInit(numdofp*numdofm); + Node **node_list = xNew(2*numnodes); + int* cs_list= xNew(2*numnodes); + + /*Find penta on bed as HO must be coupled to the dofs on the bed: */ + Element* pentabase=element->GetBasalElement(); + Element* basaltria=pentabase->SpawnBasalElement(); + + /*prepare node list*/ + for(i=0;iGetNode(i); + node_list[i+1*numnodes] = element ->GetNode(i); + cs_list[i+0*numnodes] = XYEnum; + cs_list[i+1*numnodes] = XYEnum; + } + + /*Initialize Element matrix*/ + ElementMatrix* Ke1= pentabase->NewElementMatrix(SSAApproximationEnum); + ElementMatrix* Ke2= element ->NewElementMatrix(HOApproximationEnum); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + delete Ke1; delete Ke2; + + /* Get node coordinates and dof list: */ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); + Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input); + Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + Gauss* gauss_tria=new GaussTria(); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + gauss->SynchronizeGaussBase(gauss_tria); + + element->JacobianDeterminant(&Jdet, xyz_list,gauss); + this->GetBSSAHO(B, element,xyz_list, gauss); + this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); + element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); + element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input); + + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=2*newviscosity*gauss->weight*Jdet; + for (i=0;i<3;i++) D[i][i]=D_scalar; + + TripleMultiply( B,3,numdofp,1, + &D[0][0],3,3,0, + Bprime,3,numdofm,0, + Ke_gg,1); + } + for(i=0;ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j]; + for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i]; + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list); + + /*Clean-up and return*/ + basaltria->DeleteMaterials(); delete basaltria; + + delete gauss; + delete gauss_tria; + xDelete(B); + xDelete(Bprime); + xDelete(Ke_gg); + xDelete(xyz_list); + xDelete(node_list); + xDelete(cs_list); + return Ke; + +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixFS(element); + int indices[3]={18,19,20}; + Ke1->StaticCondensation(3,&indices[0]); + int init = element->FiniteElement(); + element->SetTemporaryElementType(P1Enum); // P1 needed for HO + ElementMatrix* Ke2=CreateKMatrixHO(element); + element->SetTemporaryElementType(init); // P1 needed for HO + ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + delete Ke3; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixFS(element); + int indices[3]={18,19,20}; + Ke1->StaticCondensation(3,&indices[0]); + int init = element->FiniteElement(); + element->SetTemporaryElementType(P1Enum); + ElementMatrix* Ke2=CreateKMatrixSSA3d(element); + element->SetTemporaryElementType(init); + ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + delete Ke3; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixSSA3d(element); + ElementMatrix* Ke2=CreateKMatrixHO(element); + ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + delete Ke3; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3d(Element* element){/*{{{*/ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixSSA3dViscous(element); + ElementMatrix* Ke2=CreateKMatrixSSA3dFriction(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dFriction(Element* element){/*{{{*/ + + /*Initialize Element matrix and return if necessary*/ + if(element->IsFloating() || !element->IsOnBase()) return NULL; + + /*Build a tria element using the 3 nodes of the base of the penta. Then use + * the tria functionality to build a friction stiffness matrix on these 3 + * nodes: */ + Element* basalelement = element->SpawnBasalElement(); + ElementMatrix* Ke=CreateKMatrixSSAFriction(basalelement); + basalelement->DeleteMaterials(); delete basalelement; + + /*clean-up and return*/ + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/ + + /*Constants*/ + const int numdof2d=2*3; + + /*Intermediaries */ + int i,j,approximation; + int dim=3; + IssmDouble Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; + IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilons[6]; //6 for FS + IssmDouble B[3][numdof2d]; + IssmDouble Bprime[3][numdof2d]; + IssmDouble D[3][3]= {0.0}; // material matrix, simple scalar matrix. + IssmDouble D_scalar; + IssmDouble Ke_gg[numdof2d][numdof2d]={0.0}; + IssmDouble *xyz_list = NULL; + + /*Find penta on bed as this is a SSA elements: */ + Element* pentabase=element->GetBasalElement(); + Element* basaltria=pentabase->SpawnBasalElement(); + + /*Initialize Element matrix*/ + ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum); + element->GetInputValue(&approximation,ApproximationEnum); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); + Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input); + Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); + Input* vz_input =element->GetInput(VzEnum); _assert_(vz_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + Gauss* gauss_tria=new GaussTria(); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + gauss->SynchronizeGaussBase(gauss_tria); + + element->JacobianDeterminant(&Jdet, xyz_list,gauss); + this->GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria); + this->GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria); + + if(approximation==SSAHOApproximationEnum){ + element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); + element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + } + else if (approximation==SSAFSApproximationEnum){ + element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + } + else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet"); + + D_scalar=2*newviscosity*gauss->weight*Jdet; + for (i=0;i<3;i++) D[i][i]=D_scalar; + + TripleMultiply( &B[0][0],3,numdof2d,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof2d,0, + &Ke_gg[0][0],1); + + } + for(i=0;ivalues[i*numdof2d+j]+=Ke_gg[i][j]; + + /*Transform Coordinate System*/ + basaltria->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + xDelete(xyz_list); + delete basaltria->material; + delete basaltria; + delete gauss_tria; + delete gauss; + return Ke; + +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFS(Element* element){/*{{{*/ /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorSSA(element); - ElementVector* pe2=CreatePVectorHO(element); + ElementVector* pe1=CreatePVectorCouplingHOFSViscous(element); + ElementVector* pe2=CreatePVectorCouplingHOFSFriction(element); ElementVector* pe =new ElementVector(pe1,pe2); @@ -5980,4 +5980,338 @@ delete pe1; delete pe2; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSFriction(Element* element){/*{{{*/ + + /*Intermediaries*/ + int i,approximation; + int dim=3; + IssmDouble Jdet,Jdet2d,FSreconditioning; + IssmDouble bed_normal[3]; + IssmDouble viscosity, w, alpha2_gauss; + IssmDouble dw[3]; + IssmDouble *xyz_list_tria = NULL; + IssmDouble *xyz_list = NULL; + IssmDouble basis[6]; //for the six nodes of the penta + + /*Initialize Element vector and return if necessary*/ + if(!element->IsOnBase() || element->IsFloating()) return NULL; + element->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=HOFSApproximationEnum) return NULL; + + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int numnodes = vnumnodes+pnumnodes; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + Node **node_list = xNew(vnumnodes+pnumnodes); + for(i=0;iGetNode(i); + } + for(i=0;iGetNode(vnumnodes+i); + } + + ElementVector* pe=element->NewElementVector(FSvelocityEnum); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinatesBase(&xyz_list_tria); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + Input* vx_input= element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input= element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input= element->GetInput(VzEnum); _assert_(vz_input); + Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input); + + /*build friction object, used later on: */ + Friction* friction=new Friction(element,3); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + Gauss* gauss=element->NewGaussBase(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss); + element->NodalFunctionsP1(basis, gauss); + + vzHO_input->GetInputValue(&w, gauss); + vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); + + element->NormalBase(&bed_normal[0],xyz_list_tria); + element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + friction->GetAlpha2(&alpha2_gauss,gauss); + + for(i=0;i<3;i++){ + pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; + pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; + pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); + + /*Clean up and return*/ + xDelete(cs_list); + xDelete(node_list); + xDelete(xyz_list); + xDelete(xyz_list_tria); + delete gauss; + delete friction; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSViscous(Element* element){/*{{{*/ + + /*Intermediaries */ + int i,approximation; + int dim=3; + IssmDouble viscosity,Jdet,FSreconditioning; + IssmDouble dw[3]; + IssmDouble *xyz_list = NULL; + IssmDouble basis[6]; //for the six nodes of the penta + IssmDouble dbasis[3][6]; //for the six nodes of the penta + + /*Initialize Element vector and return if necessary*/ + element->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=HOFSApproximationEnum) return NULL; + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + Node **node_list = xNew(vnumnodes+pnumnodes); + for(i=0;iGetNode(i); + } + for(i=0;iGetNode(vnumnodes+i); + } + ElementVector* pe = element->NewElementVector(FSvelocityEnum); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input =element->GetInput(VzEnum); _assert_(vz_input); + Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet, xyz_list,gauss); + element->NodalFunctionsP1(&basis[0],gauss); + element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); + + element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); + + for(i=0;i<6;i++){ + pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; + pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; + pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); + pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); + + /*Clean up and return*/ + xDelete(cs_list); + xDelete(node_list); + xDelete(xyz_list); + delete gauss; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFS(Element* element){/*{{{*/ + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorCouplingSSAFSViscous(element); + ElementVector* pe2=CreatePVectorCouplingSSAFSFriction(element); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSFriction(Element* element){/*{{{*/ + + /*Intermediaries*/ + int i,j,approximation; + int dim=3; + IssmDouble Jdet,Jdet2d,FSreconditioning; + IssmDouble bed_normal[3]; + IssmDouble viscosity, w, alpha2_gauss; + IssmDouble dw[3]; + IssmDouble basis[6]; //for the six nodes of the penta + IssmDouble *xyz_list_tria = NULL; + IssmDouble *xyz_list = NULL; + + /*Initialize Element vector and return if necessary*/ + if(!element->IsOnBase() || element->IsFloating()) return NULL; + element->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=SSAFSApproximationEnum) return NULL; + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + Node **node_list = xNew(vnumnodes+pnumnodes); + for(i=0;iGetNode(i); + } + for(i=0;iGetNode(vnumnodes+i); + } + ElementVector* pe=element->NewElementVector(FSvelocityEnum); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinatesBase(&xyz_list_tria); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + Input* vx_input= element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input= element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input= element->GetInput(VzEnum); _assert_(vz_input); + Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input); + + /*build friction object, used later on: */ + Friction* friction=new Friction(element,3); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + Gauss* gauss=element->NewGaussBase(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss); + element->NodalFunctionsP1(basis, gauss); + + vzSSA_input->GetInputValue(&w, gauss); + vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); + + element->NormalBase(&bed_normal[0],xyz_list_tria); + element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + friction->GetAlpha2(&alpha2_gauss,gauss); + + for(i=0;i<3;i++){ + pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; + pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; + pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); + + /*Clean up and return*/ + xDelete(cs_list); + xDelete(xyz_list); + xDelete(xyz_list_tria); + xDelete(node_list); + delete gauss; + delete friction; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSViscous(Element* element){/*{{{*/ + + /*Intermediaries */ + int i,approximation; + IssmDouble viscosity,Jdet,FSreconditioning; + IssmDouble dw[3]; + IssmDouble *xyz_list = NULL; + IssmDouble basis[6]; //for the six nodes of the penta + IssmDouble dbasis[3][6]; //for the six nodes of the penta + + /*Initialize Element vector and return if necessary*/ + element->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=SSAFSApproximationEnum) return NULL; + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + Node **node_list = xNew(vnumnodes+pnumnodes); + for(i=0;iGetNode(i); + } + for(i=0;iGetNode(vnumnodes+i); + } + ElementVector* pe=element->NewElementVector(FSvelocityEnum); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input =element->GetInput(VzEnum); _assert_(vz_input); + Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsP1(&basis[0], gauss); + element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss); + + vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); + element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); + + for(i=0;i<6;i++){ + pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; + pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; + pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); + pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); + + /*Clean up and return*/ + xDelete(cs_list); + xDelete(node_list); + xDelete(xyz_list); + delete gauss; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/ + + /*compute all load vectors for this element*/ + int init = element->FiniteElement(); + element->SetTemporaryElementType(P1Enum); + ElementVector* pe1=CreatePVectorHO(element); + element->SetTemporaryElementType(init); + ElementVector* pe2=CreatePVectorFS(element); + int indices[3]={18,19,20}; + element->SetTemporaryElementType(MINIcondensedEnum); + ElementMatrix* Ke = CreateKMatrixFS(element); + element->SetTemporaryElementType(init); + pe2->StaticCondensation(Ke,3,&indices[0]); + delete Ke; + ElementVector* pe3=CreatePVectorCouplingHOFS(element); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; return pe; }/*}}}*/ @@ -6005,32 +6339,9 @@ return pe; }/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/ +ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/ /*compute all load vectors for this element*/ - int init = element->FiniteElement(); - element->SetTemporaryElementType(P1Enum); - ElementVector* pe1=CreatePVectorHO(element); - element->SetTemporaryElementType(init); - ElementVector* pe2=CreatePVectorFS(element); - int indices[3]={18,19,20}; - element->SetTemporaryElementType(MINIcondensedEnum); - ElementMatrix* Ke = CreateKMatrixFS(element); - element->SetTemporaryElementType(init); - pe2->StaticCondensation(Ke,3,&indices[0]); - delete Ke; - ElementVector* pe3=CreatePVectorCouplingHOFS(element); - ElementVector* pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFS(Element* element){/*{{{*/ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorCouplingHOFSViscous(element); - ElementVector* pe2=CreatePVectorCouplingHOFSFriction(element); + ElementVector* pe1=CreatePVectorSSA(element); + ElementVector* pe2=CreatePVectorHO(element); ElementVector* pe =new ElementVector(pe1,pe2); @@ -6040,341 +6351,83 @@ return pe; }/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSFriction(Element* element){/*{{{*/ - - /*Intermediaries*/ - int i,approximation; - int dim=3; - IssmDouble Jdet,Jdet2d,FSreconditioning; - IssmDouble bed_normal[3]; - IssmDouble viscosity, w, alpha2_gauss; - IssmDouble dw[3]; - IssmDouble *xyz_list_tria = NULL; - IssmDouble *xyz_list = NULL; - IssmDouble basis[6]; //for the six nodes of the penta - - /*Initialize Element vector and return if necessary*/ - if(!element->IsOnBase() || element->IsFloating()) return NULL; - element->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=HOFSApproximationEnum) return NULL; - - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - int numnodes = vnumnodes+pnumnodes; - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - Node **node_list = xNew(vnumnodes+pnumnodes); - for(i=0;iGetNode(i); - } - for(i=0;iGetNode(vnumnodes+i); - } - - ElementVector* pe=element->NewElementVector(FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->GetVerticesCoordinatesBase(&xyz_list_tria); - element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input= element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input= element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input= element->GetInput(VzEnum); _assert_(vz_input); - Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input); - - /*build friction object, used later on: */ - Friction* friction=new Friction(element,3); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - Gauss* gauss=element->NewGaussBase(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss); - element->NodalFunctionsP1(basis, gauss); - - vzHO_input->GetInputValue(&w, gauss); - vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); - - element->NormalBase(&bed_normal[0],xyz_list_tria); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - friction->GetAlpha2(&alpha2_gauss,gauss); - - for(i=0;i<3;i++){ - pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; - pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; - pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; - } - } - - /*Transform coordinate system*/ - element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - xDelete(node_list); - xDelete(xyz_list); - xDelete(xyz_list_tria); - delete gauss; - delete friction; - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSViscous(Element* element){/*{{{*/ - - /*Intermediaries */ - int i,approximation; - int dim=3; - IssmDouble viscosity,Jdet,FSreconditioning; - IssmDouble dw[3]; - IssmDouble *xyz_list = NULL; - IssmDouble basis[6]; //for the six nodes of the penta - IssmDouble dbasis[3][6]; //for the six nodes of the penta - - /*Initialize Element vector and return if necessary*/ - element->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=HOFSApproximationEnum) return NULL; - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - Node **node_list = xNew(vnumnodes+pnumnodes); - for(i=0;iGetNode(i); - } - for(i=0;iGetNode(vnumnodes+i); - } - ElementVector* pe = element->NewElementVector(FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input =element->GetInput(VzEnum); _assert_(vz_input); - Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - element->JacobianDeterminant(&Jdet, xyz_list,gauss); - element->NodalFunctionsP1(&basis[0],gauss); - element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); - - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); - - for(i=0;i<6;i++){ - pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; - pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; - pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); - pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; - } - } - - /*Transform coordinate system*/ - element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - xDelete(node_list); - xDelete(xyz_list); - delete gauss; - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFS(Element* element){/*{{{*/ - - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorCouplingSSAFSViscous(element); - ElementVector* pe2=CreatePVectorCouplingSSAFSFriction(element); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSFriction(Element* element){/*{{{*/ - - /*Intermediaries*/ - int i,j,approximation; - int dim=3; - IssmDouble Jdet,Jdet2d,FSreconditioning; - IssmDouble bed_normal[3]; - IssmDouble viscosity, w, alpha2_gauss; - IssmDouble dw[3]; - IssmDouble basis[6]; //for the six nodes of the penta - IssmDouble *xyz_list_tria = NULL; - IssmDouble *xyz_list = NULL; - - /*Initialize Element vector and return if necessary*/ - if(!element->IsOnBase() || element->IsFloating()) return NULL; - element->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=SSAFSApproximationEnum) return NULL; - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - Node **node_list = xNew(vnumnodes+pnumnodes); - for(i=0;iGetNode(i); - } - for(i=0;iGetNode(vnumnodes+i); - } - ElementVector* pe=element->NewElementVector(FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->GetVerticesCoordinatesBase(&xyz_list_tria); - element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input= element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input= element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input= element->GetInput(VzEnum); _assert_(vz_input); - Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input); - - /*build friction object, used later on: */ - Friction* friction=new Friction(element,3); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - Gauss* gauss=element->NewGaussBase(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss); - element->NodalFunctionsP1(basis, gauss); - - vzSSA_input->GetInputValue(&w, gauss); - vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); - - element->NormalBase(&bed_normal[0],xyz_list_tria); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - friction->GetAlpha2(&alpha2_gauss,gauss); - - for(i=0;i<3;i++){ - pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i]; - pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i]; - pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i]; - } - } - - /*Transform coordinate system*/ - element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - xDelete(xyz_list); - xDelete(xyz_list_tria); - xDelete(node_list); - delete gauss; - delete friction; - return pe; -}/*}}}*/ -ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSViscous(Element* element){/*{{{*/ - - /*Intermediaries */ - int i,approximation; - IssmDouble viscosity,Jdet,FSreconditioning; - IssmDouble dw[3]; - IssmDouble *xyz_list = NULL; - IssmDouble basis[6]; //for the six nodes of the penta - IssmDouble dbasis[3][6]; //for the six nodes of the penta - - /*Initialize Element vector and return if necessary*/ - element->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=SSAFSApproximationEnum) return NULL; - int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes = element->NumberofNodesPressure(); - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - Node **node_list = xNew(vnumnodes+pnumnodes); - for(i=0;iGetNode(i); - } - for(i=0;iGetNode(vnumnodes+i); - } - ElementVector* pe=element->NewElementVector(FSvelocityEnum); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input =element->GetInput(VzEnum); _assert_(vz_input); - Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctionsP1(&basis[0], gauss); - element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss); - - vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); - element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); - - for(i=0;i<6;i++){ - pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i]; - pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i]; - pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]); - pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i]; - } - } - - /*Transform coordinate system*/ - element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - xDelete(node_list); - xDelete(xyz_list); - delete gauss; - return pe; -}/*}}}*/ -void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. - * For node i, Bi can be expressed in the actual coordinate system +void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. + * For node i, Bprimei can be expressed in the actual coordinate system * by: - * Bi=[ dh/dx 0 ] - * [ 0 dh/dy ] - * [ 1/2*dh/dy 1/2*dh/dx ] + * Bprimei=[ 2*dh/dx dh/dy 0 0 ] + * [ dh/dx 2*dh/dy 0 0 ] + * [ dh/dy dh/dx 0 0 ] * where h is the interpolation function for node i. * - * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) + * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1) */ + int i; + IssmDouble dbasismini[3][7]; + + /*Get dbasis in actual coordinate system: */ + element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss); + + /*Build Bprime: */ + for(i=0;i<6;i++){ + Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i]; + Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i]; + Bprime[(3*7+6)*0+3*i+2] = 0.; + Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i]; + Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i]; + Bprime[(3*7+6)*1+3*i+2] = 0.; + Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i]; + Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i]; + Bprime[(3*7+6)*2+3*i+2] = 0.; + } + + for(i=0;i<1;i++){ //Add zeros for the bubble function + Bprime[(3*7+6)*0+3*(6+i)+0] = 0.; + Bprime[(3*7+6)*0+3*(6+i)+1] = 0.; + Bprime[(3*7+6)*0+3*(6+i)+2] = 0.; + Bprime[(3*7+6)*1+3*(6+i)+0] = 0.; + Bprime[(3*7+6)*1+3*(6+i)+1] = 0.; + Bprime[(3*7+6)*1+3*(6+i)+2] = 0.; + Bprime[(3*7+6)*2+3*(6+i)+0] = 0.; + Bprime[(3*7+6)*2+3*(6+i)+1] = 0.; + Bprime[(3*7+6)*2+3*(6+i)+2] = 0.; + } + + for(i=0;i<6;i++){ //last column not for the bubble function + Bprime[(3*7+6)*0+7*3+i] = 0.; + Bprime[(3*7+6)*1+7*3+i] = 0.; + Bprime[(3*7+6)*2+7*3+i] = 0.; + } +}/*}}}*/ +void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. + * For node i, Bprimei can be expressed in the actual coordinate system + * by: + * Bprimei=[ dN/dx 0 ] + * [ 0 dN/dy ] + * [ dN/dy dN/dx ] + N [ dN/dx dN/dy ] + * where N is the finiteelement function for node i. + * + * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); - IssmDouble* dbasis=xNew(3*numnodes); - - /*Get dbasis in actual coordinate system: */ + + /*Get nodal functions*/ + IssmDouble* dbasis=xNew(2*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - /*Build B: */ + /*Build Bprime: */ for(int i=0;i(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system @@ -6440,54 +6493,5 @@ } }/*}}}*/ -void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. - * For node i, Bprimei can be expressed in the actual coordinate system - * by: - * Bprimei=[ 2*dh/dx dh/dy 0 0 ] - * [ dh/dx 2*dh/dy 0 0 ] - * [ dh/dy dh/dx 0 0 ] - * where h is the interpolation function for node i. - * - * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1) - */ - - int i; - IssmDouble dbasismini[3][7]; - - /*Get dbasis in actual coordinate system: */ - element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss); - - /*Build Bprime: */ - for(i=0;i<6;i++){ - Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i]; - Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i]; - Bprime[(3*7+6)*0+3*i+2] = 0.; - Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i]; - Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i]; - Bprime[(3*7+6)*1+3*i+2] = 0.; - Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i]; - Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i]; - Bprime[(3*7+6)*2+3*i+2] = 0.; - } - - for(i=0;i<1;i++){ //Add zeros for the bubble function - Bprime[(3*7+6)*0+3*(6+i)+0] = 0.; - Bprime[(3*7+6)*0+3*(6+i)+1] = 0.; - Bprime[(3*7+6)*0+3*(6+i)+2] = 0.; - Bprime[(3*7+6)*1+3*(6+i)+0] = 0.; - Bprime[(3*7+6)*1+3*(6+i)+1] = 0.; - Bprime[(3*7+6)*1+3*(6+i)+2] = 0.; - Bprime[(3*7+6)*2+3*(6+i)+0] = 0.; - Bprime[(3*7+6)*2+3*(6+i)+1] = 0.; - Bprime[(3*7+6)*2+3*(6+i)+2] = 0.; - } - - for(i=0;i<6;i++){ //last column not for the bubble function - Bprime[(3*7+6)*0+7*3+i] = 0.; - Bprime[(3*7+6)*1+7*3+i] = 0.; - Bprime[(3*7+6)*2+7*3+i] = 0.; - } -}/*}}}*/ -void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system @@ -6521,34 +6525,30 @@ xDelete(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. - * For node i, Bprimei can be expressed in the actual coordinate system +void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. + * For node i, Bi can be expressed in the actual coordinate system * by: - * Bprimei=[ dN/dx 0 ] - * [ 0 dN/dy ] - * [ dN/dy dN/dx ] - N [ dN/dx dN/dy ] - * where N is the finiteelement function for node i. + * Bi=[ dh/dx 0 ] + * [ 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx ] + * where h is the interpolation function for node i. * - * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes) + * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) */ - /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions*/ - IssmDouble* dbasis=xNew(2*numnodes); + IssmDouble* dbasis=xNew(3*numnodes); + + /*Get dbasis in actual coordinate system: */ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - /*Build Bprime: */ + /*Build B: */ for(int i=0;i(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ - /* - * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. - * For node i, Li can be expressed in the actual coordinate system +void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ + /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. + * For node i, Lpi can be expressed in the actual coordinate system * by: - * Li=[ h 0 ] - * [ 0 h ] - * [ h 0 ] - * [ 0 h ] - * [ h 0 ] - * [ 0 h ] - * [ h 0 ] - * [ 0 h ] + * Lpi=[ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] * where h is the interpolation function for node i. */ - int num_dof=2; IssmDouble basis[3]; @@ -6584,25 +6578,17 @@ basis[2]=gauss->coord3*(1-gauss->coord4)/2.0; - /*Build LFS: */ + /*Build LprimeFS: */ for(int i=0;i<3;i++){ - LFS[num_dof*3*0+num_dof*i+0] = basis[i]; - LFS[num_dof*3*0+num_dof*i+1] = 0; - LFS[num_dof*3*1+num_dof*i+0] = 0; - LFS[num_dof*3*1+num_dof*i+1] = basis[i]; - LFS[num_dof*3*2+num_dof*i+0] = basis[i]; - LFS[num_dof*3*2+num_dof*i+1] = 0; - LFS[num_dof*3*3+num_dof*i+0] = 0; - LFS[num_dof*3*3+num_dof*i+1] = basis[i]; - LFS[num_dof*3*4+num_dof*i+0] = basis[i]; - LFS[num_dof*3*4+num_dof*i+1] = 0; - LFS[num_dof*3*5+num_dof*i+0] = 0; - LFS[num_dof*3*5+num_dof*i+1] = basis[i]; - LFS[num_dof*3*6+num_dof*i+0] = basis[i]; - LFS[num_dof*3*6+num_dof*i+1] = 0; - LFS[num_dof*3*7+num_dof*i+0] = 0; - LFS[num_dof*3*7+num_dof*i+1] = basis[i]; - } -}/*}}}*/ -void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ + LprimeFS[num_dof*3*0+num_dof*i+0] = basis[i]; + LprimeFS[num_dof*3*0+num_dof*i+1] = 0.; + LprimeFS[num_dof*3*1+num_dof*i+0] = 0.; + LprimeFS[num_dof*3*1+num_dof*i+1] = basis[i]; + LprimeFS[num_dof*3*2+num_dof*i+0] = basis[i]; + LprimeFS[num_dof*3*2+num_dof*i+1] = 0.; + LprimeFS[num_dof*3*3+num_dof*i+0] = 0.; + LprimeFS[num_dof*3*3+num_dof*i+1] = basis[i]; + } +}/*}}}*/ +void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. * For node i, Lpi can be expressed in the actual coordinate system @@ -6709,5 +6695,5 @@ } }/*}}}*/ -void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ +void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system @@ -6748,14 +6734,20 @@ } }/*}}}*/ -void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ - /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. - * For node i, Lpi can be expressed in the actual coordinate system +void StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ + /* + * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * For node i, Li can be expressed in the actual coordinate system * by: - * Lpi=[ h 0 ] - * [ 0 h ] - * [ h 0 ] - * [ 0 h ] + * Li=[ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] + * [ h 0 ] + * [ 0 h ] * where h is the interpolation function for node i. */ + int num_dof=2; IssmDouble basis[3]; @@ -6770,17 +6762,25 @@ basis[2]=gauss->coord3*(1-gauss->coord4)/2.0; - /*Build LprimeFS: */ + /*Build LFS: */ for(int i=0;i<3;i++){ - LprimeFS[num_dof*3*0+num_dof*i+0] = basis[i]; - LprimeFS[num_dof*3*0+num_dof*i+1] = 0.; - LprimeFS[num_dof*3*1+num_dof*i+0] = 0.; - LprimeFS[num_dof*3*1+num_dof*i+1] = basis[i]; - LprimeFS[num_dof*3*2+num_dof*i+0] = basis[i]; - LprimeFS[num_dof*3*2+num_dof*i+1] = 0.; - LprimeFS[num_dof*3*3+num_dof*i+0] = 0.; - LprimeFS[num_dof*3*3+num_dof*i+1] = basis[i]; - } -}/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/ + LFS[num_dof*3*0+num_dof*i+0] = basis[i]; + LFS[num_dof*3*0+num_dof*i+1] = 0; + LFS[num_dof*3*1+num_dof*i+0] = 0; + LFS[num_dof*3*1+num_dof*i+1] = basis[i]; + LFS[num_dof*3*2+num_dof*i+0] = basis[i]; + LFS[num_dof*3*2+num_dof*i+1] = 0; + LFS[num_dof*3*3+num_dof*i+0] = 0; + LFS[num_dof*3*3+num_dof*i+1] = basis[i]; + LFS[num_dof*3*4+num_dof*i+0] = basis[i]; + LFS[num_dof*3*4+num_dof*i+1] = 0; + LFS[num_dof*3*5+num_dof*i+0] = 0; + LFS[num_dof*3*5+num_dof*i+1] = basis[i]; + LFS[num_dof*3*6+num_dof*i+0] = basis[i]; + LFS[num_dof*3*6+num_dof*i+1] = 0; + LFS[num_dof*3*7+num_dof*i+0] = 0; + LFS[num_dof*3*7+num_dof*i+1] = basis[i]; + } +}/*}}}*/ +void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/ int i; @@ -6879,5 +6879,5 @@ xDelete(cs_list); }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/ int i; @@ -6985,5 +6985,5 @@ xDelete(cs_list); }/*}}}*/ -void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/ +void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/ int i,domaintype; Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18927) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18928) @@ -13,10 +13,10 @@ public: /*Model processing*/ - int DofsPerNode(int** doflist,int domaintype,int approximation); - void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); - void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - void CreateNodes(Nodes* nodes,IoModel* iomodel); void CreateConstraints(Constraints* constraints,IoModel* iomodel); void CreateLoads(Loads* loads, IoModel* iomodel); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); /*Finite element Analysis*/ @@ -26,90 +26,90 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetSolutionFromInputs(Vector* solution,Element* element); - void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void GetSolutionFromInputsHoriz(Vector* solution,Element* element); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GetSolutionFromInputsHoriz(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); /*SSA*/ ElementMatrix* CreateJacobianMatrixSSA(Element* element); ElementMatrix* CreateKMatrixSSA(Element* element); - ElementMatrix* CreateKMatrixSSAViscous(Element* element); ElementMatrix* CreateKMatrixSSAFriction(Element* element); ElementMatrix* CreateKMatrixSSALateralFriction(Element* element); + ElementMatrix* CreateKMatrixSSAViscous(Element* element); ElementVector* CreatePVectorSSA(Element* element); + ElementVector* CreatePVectorSSAFront(Element* element); ElementVector* CreatePVectorSSADrivingStress(Element* element); - ElementVector* CreatePVectorSSAFront(Element* element); - void GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element); + void GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element); /*L1L2*/ ElementMatrix* CreateKMatrixL1L2(Element* element); + ElementMatrix* CreateKMatrixL1L2Friction(Element* element); ElementMatrix* CreateKMatrixL1L2Viscous(Element* element); - ElementMatrix* CreateKMatrixL1L2Friction(Element* element); ElementVector* CreatePVectorL1L2(Element* element); + ElementVector* CreatePVectorL1L2Front(Element* element); ElementVector* CreatePVectorL1L2DrivingStress(Element* element); - ElementVector* CreatePVectorL1L2Front(Element* element); - void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element); + void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element); /*HO*/ ElementMatrix* CreateJacobianMatrixHO(Element* element); ElementMatrix* CreateKMatrixHO(Element* element); + ElementMatrix* CreateKMatrixHOFriction(Element* element); ElementMatrix* CreateKMatrixHOViscous(Element* element); - ElementMatrix* CreateKMatrixHOFriction(Element* element); ElementVector* CreatePVectorHO(Element* element); + ElementVector* CreatePVectorHOFront(Element* element); ElementVector* CreatePVectorHODrivingStress(Element* element); - ElementVector* CreatePVectorHOFront(Element* element); - void GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element); + void GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element); /*FS*/ ElementVector* CreateDVectorFS(Element* element); ElementMatrix* CreateJacobianMatrixFS(Element* element); ElementMatrix* CreateKMatrixFS(Element* element); + ElementMatrix* CreateKMatrixFSFriction(Element* element); + ElementMatrix* CreateKMatrixFSShelf(Element* element); + ElementMatrix* CreateKMatrixFSViscous(Element* element); ElementMatrix* CreateKMatrixFSViscousLA(Element* element); ElementMatrix* CreateKMatrixFSViscousXTH(Element* element); - ElementMatrix* CreateKMatrixFSViscous(Element* element); - ElementMatrix* CreateKMatrixFSFriction(Element* element); - ElementMatrix* CreateKMatrixFSShelf(Element* element); ElementVector* CreatePVectorFS(Element* element); + ElementVector* CreatePVectorFSFriction(Element* element); + ElementVector* CreatePVectorFSFront(Element* element); + ElementVector* CreatePVectorFSShelf(Element* element); + ElementVector* CreatePVectorFSStress(Element* element); ElementVector* CreatePVectorFSViscous(Element* element); ElementVector* CreatePVectorFSViscousLA(Element* element); ElementVector* CreatePVectorFSViscousXTH(Element* element); - ElementVector* CreatePVectorFSShelf(Element* element); - ElementVector* CreatePVectorFSFront(Element* element); - ElementVector* CreatePVectorFSFriction(Element* element); - ElementVector* CreatePVectorFSStress(Element* element); - void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetCFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetSolutionFromInputsFS(Vector* solution,Element* element); - void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element); - void InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters); - void InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters); - void InitializeXTH(Elements* elements,Parameters* parameters); + void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetCFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputsFS(Vector* solution,Element* element); + void InitializeXTH(Elements* elements,Parameters* parameters); + void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element); + void InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters); + void InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters); /*Coupling*/ + ElementMatrix* CreateKMatrixCouplingHOFS(Element* element); + ElementMatrix* CreateKMatrixCouplingSSAFS(Element* element); + ElementMatrix* CreateKMatrixCouplingSSAFSFriction(Element* element); + ElementMatrix* CreateKMatrixCouplingSSAFSViscous(Element* element); + ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element); + ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element); + ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element); + ElementMatrix* CreateKMatrixHOFS(Element* element); + ElementMatrix* CreateKMatrixSSAFS(Element* element); + ElementMatrix* CreateKMatrixSSAHO(Element* element); ElementMatrix* CreateKMatrixSSA3d(Element* element); ElementMatrix* CreateKMatrixSSA3dFriction(Element* element); ElementMatrix* CreateKMatrixSSA3dViscous(Element* element); - ElementMatrix* CreateKMatrixHOFS(Element* element); - ElementMatrix* CreateKMatrixSSAHO(Element* element); - ElementMatrix* CreateKMatrixSSAFS(Element* element); - ElementMatrix* CreateKMatrixCouplingHOFS(Element* element); - ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element); - ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element); - ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element); - ElementMatrix* CreateKMatrixCouplingSSAFS(Element* element); - ElementMatrix* CreateKMatrixCouplingSSAFSFriction(Element* element); - ElementMatrix* CreateKMatrixCouplingSSAFSViscous(Element* element); + ElementVector* CreatePVectorSSAFS(Element* element); ElementVector* CreatePVectorSSAHO(Element* element); - ElementVector* CreatePVectorSSAFS(Element* element); ElementVector* CreatePVectorCouplingSSAFS(Element* element); ElementVector* CreatePVectorCouplingSSAFSFriction(Element* element); @@ -119,16 +119,16 @@ ElementVector* CreatePVectorCouplingHOFSFriction(Element* element); ElementVector* CreatePVectorCouplingHOFSViscous(Element* element); - void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetLFSSSA(IssmDouble* L,Element* element,Gauss* gauss); - void GetLSSAFS(IssmDouble* L,Element* element,Gauss* gauss); - void GetLprimeFSSSA(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetLprimeSSAFS(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element); - void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element); - void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element); + void GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetLFSSSA(IssmDouble* L,Element* element,Gauss* gauss); + void GetLprimeFSSSA(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetLprimeSSAFS(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetLSSAFS(IssmDouble* L,Element* element,Gauss* gauss); + void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element); + void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element); + void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element); }; #endif