#include "./StressbalanceAnalysis.h" #include "../toolkits/toolkits.h" #include "../classes/classes.h" #include "../shared/shared.h" #include "../modules/modules.h" #include "../solutionsequences/solutionsequences.h" #include "../cores/cores.h" //#define FSANALYTICAL 12 //#define LATERALFRICTION 1 /*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->FindConstant(&g,"md.constants.g"); iomodel->FindConstant(&rho_ice,"md.materials.rho_ice"); iomodel->FindConstant(&FSreconditioning,"md.stressbalance.FSreconditioning"); iomodel->FindConstant(&isSIA,"md.flowequation.isSIA"); iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); /*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->FindConstant(&finiteelement,"md.flowequation.fe_SSA"); else if(isL1L2) finiteelement = P1Enum; else if(isHO) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_HO"); else if(isFS){ iomodel->FindConstant(&finiteelement,"md.flowequation.fe_FS"); /*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,"md.flowequation.vertex_equation"); iomodel->FetchData(&nodeonFS,NULL,NULL,"md.flowequation.borderFS"); iomodel->FetchData(&nodeonbase,NULL,NULL,"md.mesh.vertexonbase"); iomodel->FetchData(&groundedice_ls,NULL,NULL,"md.mask.groundedice_levelset"); if(iomodel->domaintype==Domain3DEnum){ iomodel->FetchData(&spcvz,&Mz,&Nz,"md.stressbalance.spcvz"); } else if (iomodel->domaintype==Domain2DverticalEnum){ iomodel->FetchData(&spcvz,&Mz,&Nz,"md.stressbalance.spcvy"); } else{ _error_("not supported yet"); } if(iomodel->domaintype==Domain3DEnum){ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,1); IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2); iomodel->DeleteData(spcvz,"md.stressbalance.spcvz"); } else if (iomodel->domaintype==Domain2DverticalEnum){ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,1); iomodel->DeleteData(spcvz,"md.stressbalance.spcvy"); } else{ _error_("not supported yet"); } iomodel->DeleteData(vertices_type,"md.flowequation.vertex_equation"); iomodel->DeleteData(nodeonFS,"md.flowequation.borderFS"); iomodel->DeleteData(nodeonbase,"md.mesh.vertexonbase"); iomodel->DeleteData(groundedice_ls,"md.mask.groundedice_levelset"); /*Pressure spc*/ count = constraints->Size(); iomodel->FetchData(&vertices_type,NULL,NULL,"md.flowequation.vertex_equation"); iomodel->FetchData(&surface,NULL,NULL,"md.geometry.surface"); iomodel->FetchData(&z,NULL,NULL,"md.mesh.z"); switch(finiteelement){ case P1Enum: for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ if(IoCodeToEnumVertexEquation(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(IoCodeToEnumVertexEquation(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(IoCodeToEnumVertexEquation(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(IoCodeToEnumVertexEquation(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,"md.flowequation.vertex_equation"); iomodel->DeleteData(surface,"md.geometry.surface"); iomodel->DeleteData(z,"md.mesh.z"); } else{ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); if(iomodel->domaintype!=Domain2DverticalEnum){ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,1); } } return; } /*Constraints: fetch data: */ iomodel->FetchData(&spcvx,&Mx,&Nx,"md.stressbalance.spcvx"); iomodel->FetchData(&spcvy,&My,&Ny,"md.stressbalance.spcvy"); iomodel->FetchData(&spcvz,&Mz,&Nz,"md.stressbalance.spcvz"); iomodel->FetchData(&nodeonSSA,NULL,NULL,"md.flowequation.borderSSA"); if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonHO,NULL,NULL,"md.flowequation.borderHO"); if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonFS,NULL,NULL,"md.flowequation.borderFS"); if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonbase,NULL,NULL,"md.mesh.vertexonbase"); if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&groundedice_ls,NULL,NULL,"md.mask.groundedice_levelset"); iomodel->FetchData(&vertices_type,NULL,NULL,"md.flowequation.vertex_equation"); iomodel->FetchData(&surface,NULL,NULL,"md.geometry.surface"); iomodel->FetchData(&z,NULL,NULL,"md.mesh.z"); /*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(IoCodeToEnumVertexEquation(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 (IoCodeToEnumVertexEquation(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 (IoCodeToEnumVertexEquation(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 (IoCodeToEnumVertexEquation(reCast(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 (IoCodeToEnumVertexEquation(reCast(vertices_type[i]))==SIAApproximationEnum){ constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum)); count++; } if (IoCodeToEnumVertexEquation(reCast(vertices_type[i]))==FSApproximationEnum || (IoCodeToEnumVertexEquation(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 (IoCodeToEnumVertexEquation(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,"md.stressbalance.spcvx"); iomodel->DeleteData(spcvy,"md.stressbalance.spcvy"); iomodel->DeleteData(spcvz,"md.stressbalance.spcvz"); iomodel->DeleteData(nodeonSSA,"md.flowequation.borderSSA"); if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonHO,"md.flowequation.borderHO"); if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonFS,"md.flowequation.borderFS"); if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonbase,"md.mesh.vertexonbase"); if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(groundedice_ls,"md.mask.groundedice_levelset"); iomodel->DeleteData(vertices_type,"md.flowequation.vertex_equation"); iomodel->DeleteData(surface,"md.geometry.surface"); iomodel->DeleteData(z,"md.mesh.z"); /*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->FindConstant(&isL1L2,"md.flowequation.isL1L2"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&numrifts,"md.rifts.numrifts"); /*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,"md.stressbalance.vertex_pairing"); 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,"md.stressbalance.vertex_pairing"); /*Create Riffront loads for rifts: */ if(numrifts){ iomodel->FetchData(&riftinfo,&numriftsegments,NULL,"md.rifts.riftstruct"); iomodel->FetchData(5,"md.rifts.riftstruct","md.geometry.thickness","md.geometry.base","md.geometry.surface","md.mask.groundedice_levelset"); 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,"md.rifts.riftstruct","md.geometry.thickness","md.geometry.base","md.geometry.surface","md.mask.groundedice_levelset"); xDelete(riftinfo); } }/*}}}*/ void StressbalanceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ /*Intermediary*/ bool isSSA,isL1L2,isHO,isFS,iscoupling; int finiteelement=-1,approximation=-1; /*Fetch parameters: */ iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); /*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->FindConstant(&finiteelement,"md.flowequation.fe_SSA"); } else if(isL1L2){ approximation = L1L2ApproximationEnum; finiteelement = P1Enum; } else if(isHO){ approximation = HOApproximationEnum; iomodel->FindConstant(&finiteelement,"md.flowequation.fe_HO"); } else if(isFS){ approximation = FSApproximationEnum; iomodel->FindConstant(&finiteelement,"md.flowequation.fe_FS"); } iomodel->FetchData(3,"md.flowequation.borderSSA","md.flowequation.vertex_equation","md.stressbalance.referential"); if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(3,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderFS"); ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation); iomodel->DeleteData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.vertex_equation", "md.stressbalance.referential","md.flowequation.borderFS"); } else{ /*Coupling: we are going to create P1 Elements only*/ Node* node = NULL; int lid=0; if(!nodes) nodes = new Nodes(); iomodel->FetchData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS", "md.flowequation.vertex_equation","md.stressbalance.referential"); if(isFS){ /*P1+ velocity*/ for(int i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ approximation=IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[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=IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[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,IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i])))); } } } iomodel->DeleteData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS", "md.flowequation.vertex_equation","md.stressbalance.referential"); } }/*}}}*/ int StressbalanceAnalysis::DofsPerNode(int** pdoftype,int domaintype,int approximation){/*{{{*/ /*output*/ int *doftype = NULL; int numdofs; switch(approximation){ case SSAApproximationEnum: switch(domaintype){ case Domain3DEnum: numdofs=2; break; case Domain2DhorizontalEnum: numdofs=2; break; case Domain2DverticalEnum: numdofs=1; break; default: _error_("mesh type not supported yet"); } break; case L1L2ApproximationEnum: numdofs =2; break; case HOApproximationEnum: switch(domaintype){ case Domain3DEnum: numdofs=2; break; case Domain2DverticalEnum: numdofs=1; break; default: _error_("mesh type not supported yet"); } break; case SIAApproximationEnum: numdofs =2; break; case FSvelocityEnum: switch(domaintype){ case Domain3DEnum: numdofs=3; break; case Domain2DverticalEnum: numdofs=2; break; default: _error_("mesh type not supported yet"); } break; case FSpressureEnum: numdofs=1; break; case NoneApproximationEnum: switch(domaintype){ case Domain3DEnum: numdofs=4; break; case Domain2DverticalEnum: numdofs=3; break; default: _error_("mesh type not supported yet"); } break; case SSAHOApproximationEnum: numdofs=4; doftype=xNew(numdofs); doftype[0]=SSAApproximationEnum; doftype[1]=SSAApproximationEnum; doftype[2]=HOApproximationEnum; doftype[3]=HOApproximationEnum; break; case HOFSApproximationEnum: numdofs=5; doftype=xNew(numdofs); doftype[0]=HOApproximationEnum; doftype[1]=HOApproximationEnum; doftype[2]=FSvelocityEnum; doftype[3]=FSvelocityEnum; doftype[4]=FSvelocityEnum; break; case SSAFSApproximationEnum: numdofs=5; doftype=xNew(numdofs); doftype[0]=SSAApproximationEnum; doftype[1]=SSAApproximationEnum; doftype[2]=FSvelocityEnum; doftype[3]=FSvelocityEnum; doftype[4]=FSvelocityEnum; break; default: _error_("Approximation " << EnumToStringx(approximation) << " not implemented yet"); } /*Assign output pointer and return*/ *pdoftype = doftype; return numdofs; }/*}}}*/ void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ /*Intermediaries*/ int materials_type,finiteelement,fe_FS; int approximation,frictionlaw; int FrictionCoupling; int* finiteelement_list=NULL; bool isSSA,isL1L2,isHO,isFS,iscoupling; bool control_analysis; bool dakota_analysis; bool ismovingfront; /*Fetch constants needed: */ iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol"); iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota"); iomodel->FindConstant(&materials_type,"md.materials.type"); iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront"); iomodel->FindConstant(&frictionlaw,"md.friction.law"); /*return if no processing required*/ if(!isSSA & !isL1L2 & !isHO & !isFS) return; /*Fetch data needed and allocate vectors: */ iomodel->FetchData(1,"md.flowequation.element_equation"); finiteelement_list=xNewZeroInit(iomodel->numberofelements); /*Do we have coupling*/ if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) iscoupling = true; else iscoupling = false; /*Get finite element type*/ if(!iscoupling){ if(isSSA) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_SSA"); else if(isL1L2) finiteelement = P1Enum; else if(isHO) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_HO"); else if(isFS) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_FS"); for(int i=0;inumberofelements;i++){ finiteelement_list[i]=finiteelement; } } else{ if(isFS){ for(int i=0;inumberofelements;i++){ approximation=IoCodeToEnumElementEquation(reCast(iomodel->Data("md.flowequation.element_equation")[i])); if(approximation==FSApproximationEnum || approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ finiteelement_list[i]=MINIcondensedEnum; } else{ finiteelement_list[i]=P1Enum; } } } else{ finiteelement = P1Enum; for(int i=0;inumberofelements;i++){ finiteelement_list[i]=finiteelement; } } } /*Update elements: */ int counter=0; for(int i=0;inumberofelements;i++){ if(iomodel->my_elements[i]){ Element* element=(Element*)elements->GetObjectByOffset(counter); element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]); counter++; } } /*Create inputs: */ iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum); iomodel->FetchDataToInput(elements,"md.geometry.surface",SurfaceEnum); iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum); iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0); iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum,0.); iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum,0.); iomodel->FetchDataToInput(elements,"md.stressbalance.loadingforcex",LoadingforceXEnum); iomodel->FetchDataToInput(elements,"md.stressbalance.loadingforcey",LoadingforceYEnum); #ifdef LATERALFRICTION iomodel->FetchDataToInput(elements,"md.mesh.vertexonboundary",MeshVertexonboundaryEnum); #endif if(iomodel->domaintype!=Domain2DhorizontalEnum){ iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); } if(iomodel->domaintype==Domain3DEnum){ iomodel->FetchDataToInput(elements,"md.flowequation.borderFS",FlowequationBorderFSEnum); iomodel->FetchDataToInput(elements,"md.stressbalance.loadingforcez",LoadingforceZEnum); iomodel->FetchDataToInput(elements,"md.initialization.vz",VzEnum,0.); } if(isFS){ iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum,0.); iomodel->FetchDataToInput(elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum,0.); } /*LATH parameters*/ iomodel->FindConstant(&fe_FS,"md.flowequation.fe_FS"); if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){ InputUpdateFromConstantx(elements,0.,SigmaNNEnum); } /*Friction law variables*/ switch(frictionlaw){ case 1: iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling"); iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); if(FrictionCoupling==1){ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); } break; case 2: iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum); iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum); break; case 3: iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling"); iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum); iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); if(FrictionCoupling==1){ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); } break; case 4: iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum); iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum); break; case 5: iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); iomodel->FetchDataToInput(elements,"md.friction.water_layer",FrictionWaterLayerEnum); break; case 6: iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum); iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum); iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum); iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum); break; case 7: iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum); iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); break; case 9: iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum); iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum); InputUpdateFromConstantx(elements,1.,FrictionPEnum); InputUpdateFromConstantx(elements,1.,FrictionQEnum); break; default: _error_("not supported"); } #ifdef _HAVE_ANDROID_ elements->InputDuplicate(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum); #endif /*Free data: */ iomodel->DeleteData(1,"md.flowequation.element_equation"); xDelete(finiteelement_list); }/*}}}*/ 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("md.flowequation.isSIA",FlowequationIsSIAEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isSSA",FlowequationIsSSAEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isL1L2",FlowequationIsL1L2Enum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isHO",FlowequationIsHOEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.fe_FS",FlowequationFeFSEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.restol",StressbalanceRestolEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.reltol",StressbalanceReltolEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.abstol",StressbalanceAbstolEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isnewton",StressbalanceIsnewtonEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.maxiter",StressbalanceMaxiterEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.penalty_factor",StressbalancePenaltyFactorEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.rift_penalty_threshold",StressbalanceRiftPenaltyThresholdEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.FSreconditioning",StressbalanceFSreconditioningEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.shelf_dampening",StressbalanceShelfDampeningEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.viscosity_overshoot",StressbalanceViscosityOvershootEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum)); /*XTH LATH parameters*/ iomodel->FindConstant(&fe_FS,"md.flowequation.fe_FS"); if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){ parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.augmented_lagrangian_r",AugmentedLagrangianREnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.augmented_lagrangian_rlambda",AugmentedLagrangianRlambdaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.XTH_theta",AugmentedLagrangianThetaEnum)); } iomodel->FindConstant(&materials_type,"md.materials.type"); if(materials_type==MatdamageiceEnum){ parameters->AddObject(iomodel->CopyConstantObject("md.damage.law",DamageLawEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.damage.kappa",DamageKappaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.damage.stress_threshold",DamageStressThresholdEnum)); } /*Requested outputs*/ iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.stressbalance.requested_outputs"); parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs)); if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs)); iomodel->DeleteData(&requestedoutputs,numoutputs,"md.stressbalance.requested_outputs"); /*Deal with friction parameters*/ int frictionlaw; iomodel->FindConstant(&frictionlaw,"md.friction.law"); if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); if(frictionlaw==3 || frictionlaw==1) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); if(frictionlaw==5) parameters->AddObject(iomodel->CopyConstantObject("md.friction.f",FrictionFEnum)); if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); }/*}}}*/ /*Finite Element Analysis*/ void StressbalanceAnalysis::Core(FemModel* femmodel){/*{{{*/ /*Intermediaries*/ bool isSSA,isL1L2,isHO,isFS; bool conserve_loads = true; int newton,domaintype,fe_FS; /* recover parameters:*/ femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum); femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum); femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); femmodel->parameters->FindParam(&fe_FS,FlowequationFeFSEnum); femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum); if(isFS && !(isSSA || isHO || isL1L2)){ if(VerboseSolution()) _printf0_(" computing velocities\n"); femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); if (fe_FS==XTaylorHoodEnum) solutionsequence_la_theta(femmodel); else if (fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum) solutionsequence_la(femmodel); else if(newton>0) solutionsequence_newton(femmodel); else solutionsequence_nonlinear(femmodel,conserve_loads); } else if(!isFS && (isSSA || isHO || isL1L2)){ if(VerboseSolution()) _printf0_(" computing velocities\n"); femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); if(newton>0) solutionsequence_newton(femmodel); else solutionsequence_nonlinear(femmodel,conserve_loads); if(domaintype==Domain2DverticalEnum && isSSA){ femmodel->parameters->SetParam(VxEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); femmodel->parameters->SetParam(VelEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); } } else if ((isSSA || isL1L2 || isHO) && isFS){ if(VerboseSolution()) _printf0_(" computing coupling between lower order models and FS\n"); solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads); } else{ _error_("not supported"); } }/*}}}*/ ElementVector* StressbalanceAnalysis::CreateDVector(Element* element){/*{{{*/ int approximation; element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ case FSApproximationEnum: return CreateDVectorFS(element); default: return NULL; //no need for doftypes outside of FS approximation } return NULL; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ int approximation; element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ case SSAApproximationEnum: return CreateJacobianMatrixSSA(element); case HOApproximationEnum: return CreateJacobianMatrixHO(element); case FSApproximationEnum: return CreateJacobianMatrixFS(element); case NoneApproximationEnum: return NULL; default: _error_("Approximation "<GetInputValue(&approximation,ApproximationEnum); switch(approximation){ case SIAApproximationEnum: return NULL; case SSAApproximationEnum: return CreateKMatrixSSA(element); case L1L2ApproximationEnum: return CreateKMatrixL1L2(element); case HOApproximationEnum: return CreateKMatrixHO(element); case FSApproximationEnum: return CreateKMatrixFS(element); case SSAHOApproximationEnum: return CreateKMatrixSSAHO(element); case HOFSApproximationEnum: return CreateKMatrixHOFS(element); case SSAFSApproximationEnum: return CreateKMatrixSSAFS(element); case NoneApproximationEnum: return NULL; default: _error_("Approximation "<GetInputValue(&approximation,ApproximationEnum); switch(approximation){ case SIAApproximationEnum: return NULL; case SSAApproximationEnum: return CreatePVectorSSA(element); case L1L2ApproximationEnum: return CreatePVectorL1L2(element); case HOApproximationEnum: return CreatePVectorHO(element); case FSApproximationEnum: return CreatePVectorFS(element); case SSAHOApproximationEnum: return CreatePVectorSSAHO(element); case HOFSApproximationEnum: return CreatePVectorHOFS(element); case SSAFSApproximationEnum: return CreatePVectorSSAFS(element); case NoneApproximationEnum: return NULL; default: _error_("Approximation "<* solution,Element* element){/*{{{*/ int approximation; element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ case FSApproximationEnum: case NoneApproximationEnum: GetSolutionFromInputsFS(solution,element); return; case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum: GetSolutionFromInputsHoriz(solution,element); return; case L1L2ApproximationEnum: GetSolutionFromInputsHoriz(solution,element); return; case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum: /*the elements around will create the solution*/ return; default: _error_("Approximation "<* solution,Element* element){/*{{{*/ IssmDouble vx,vy; int domaintype,dim,approximation,dofpernode; int* doflist = NULL; /*Get some parameters*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DhorizontalEnum: dim = 2; dofpernode = 2; break; case Domain2DverticalEnum: dim = 2; dofpernode = 1; break; case Domain3DEnum: dim = 3; dofpernode = 2; break; default: _error_("mesh "<GetNumberOfNodes(); int numdof = numnodes*dofpernode; element->GetInputValue(&approximation,ApproximationEnum); /*Fetch dof list and allocate solution vector*/ element->GetDofList(&doflist,approximation,GsetEnum); IssmDouble* values = xNew(numdof); /*Get inputs*/ Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input=NULL; if(domaintype!=Domain2DverticalEnum){vy_input=element->GetInput(VyEnum); _assert_(vy_input);} /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ Gauss* gauss=element->NewGauss(); for(int i=0;iGaussNode(element->FiniteElement(),i); /*Recover vx and vy*/ vx_input->GetInputValue(&vx,gauss); values[i*dofpernode+0]=vx; if(dofpernode==2){ vy_input->GetInputValue(&vy,gauss); values[i*dofpernode+1]=vy; } } solution->SetValues(numdof,doflist,values,INS_VAL); /*Free ressources:*/ delete gauss; xDelete(values); xDelete(doflist); }/*}}}*/ void StressbalanceAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ _error_("Not implemented yet"); }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int approximation; element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ case FSApproximationEnum: case NoneApproximationEnum: InputUpdateFromSolutionFS(solution,element); return; case SIAApproximationEnum: return; case SSAApproximationEnum: InputUpdateFromSolutionSSA(solution,element); return; case HOApproximationEnum: InputUpdateFromSolutionHO(solution,element); return; case L1L2ApproximationEnum: InputUpdateFromSolutionL1L2(solution,element); return; case SSAHOApproximationEnum: InputUpdateFromSolutionSSAHO(solution,element); return; case HOFSApproximationEnum: InputUpdateFromSolutionHOFS(solution,element); return; case SSAFSApproximationEnum: InputUpdateFromSolutionSSAFS(solution,element); return; default: _error_("Approximation "<FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DhorizontalEnum: basalelement = element; break; case Domain3DEnum: if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; default: _error_("mesh "<GetNumberOfNodes(); /*Initialize Element matrix, vectors and Gaussian points*/ ElementMatrix* Ke=this->CreateKMatrixSSA(element); //Initialize Jacobian with regular SSA (first part of the Gateau derivative) IssmDouble* dbasis = xNew(2*numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* thickness_input = basalelement->GetInput(ThicknessEnum);_assert_(thickness_input); Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); /* Start looping on the number of gaussian points: */ Gauss* gauss = basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); thickness_input->GetInputValue(&thickness, gauss); basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); basalelement->material->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; for(int i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps1dotdphii; Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps1dotdphii; Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps2dotdphii; Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps2dotdphii; } } } /*Transform Coordinate System*/ basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(dbasis); /*clean-up and return*/ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries*/ int domaintype; Element* basalelement; /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DhorizontalEnum: basalelement = element; break; case Domain3DEnum: case Domain2DverticalEnum: if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; default: _error_("mesh "<DeleteMaterials(); delete basalelement;}; delete Ke1; delete Ke2; return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/ /*Return if element is inactive*/ if(element->IsFloating() || !element->IsIceInElement()) return NULL; /*Intermediaries*/ int dim,domaintype; bool mainlyfloating; int migration_style,point1; IssmDouble alpha2,Jdet,fraction1,fraction2; IssmDouble gllevelset,phi=1.; IssmDouble *xyz_list = NULL; Gauss* gauss = NULL; /*Get problem dimension*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DverticalEnum: dim = 1;break; case Domain2DhorizontalEnum: dim = 2;break; case Domain3DEnum: dim = 2;break; default: _error_("mesh "<GetNumberOfNodes(); int numdof = numnodes*dim; /*Initialize Element matrix and vectors*/ ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum); IssmDouble* B = xNew(dim*numdof); IssmDouble* D = xNewZeroInit(dim*dim); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); element->FindParam(&migration_style,GroundinglineMigrationEnum); Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); Input* gllevelset_input = NULL; /*build friction object, used later on: */ Friction* friction=new Friction(element,dim); /*Recover portion of element that is grounded*/ if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list); 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); } else{ gauss = element->NewGauss(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->GetBSSAFriction(B,element,dim,xyz_list,gauss); element->JacobianDeterminant(&Jdet, xyz_list,gauss); for(int i=0;iweight*Jdet; TripleMultiply(B,dim,numdof,1, D,dim,dim,0, B,dim,numdof,0, &Ke->values[0],1); } /*Transform Coordinate System*/ if(dim==2) element->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ delete gauss; delete friction; xDelete(xyz_list); xDelete(B); xDelete(D); return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSALateralFriction(Element* element){/*{{{*/ /*Return if element is inactive*/ if(!element->IsIceInElement()) return NULL; /*If no boundary, return NULL*/ if(!element->IsFaceOnBoundary()) return NULL; /*Intermediaries*/ IssmDouble alpha2; IssmDouble Jdet; int domaintype; IssmDouble icefront; IssmDouble *xyz_list = NULL; IssmDouble *xyz_list_boundary = NULL; /*Get problem dimension*/ element->FindParam(&domaintype,DomainTypeEnum); if(domaintype==Domain2DverticalEnum) return NULL; //not supported yet /*Fetch number of nodes and dof for this finite element*/ int dim = 2; int numnodes = element->GetNumberOfNodes(); int numdof = numnodes*dim; /*Initialize Element matrix and vectors*/ ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum); IssmDouble* B = xNew(dim*numdof); IssmDouble* D = xNewZeroInit(dim*dim); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); element->GetLevelCoordinates(&xyz_list_boundary,xyz_list,MeshVertexonboundaryEnum,1.); Input* icelevelset_input = element->GetInput(MaskIceLevelsetEnum); _assert_(icelevelset_input); /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); this->GetBSSAFriction(B,element,dim,xyz_list,gauss); element->JacobianDeterminantSurface(&Jdet,xyz_list_boundary,gauss); icelevelset_input->GetInputValue(&icefront, gauss); if(icefront==0.) alpha2=0.; else alpha2=2.e+12; for(int i=0;iweight*Jdet; TripleMultiply(B,dim,numdof,1, D,dim,dim,0, B,dim,numdof,0, &Ke->values[0],1); } /*Transform Coordinate System*/ if(dim==2) element->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(xyz_list_boundary); xDelete(B); xDelete(D); return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries*/ int dim,domaintype,bsize; IssmDouble viscosity,newviscosity,oldviscosity; IssmDouble viscosity_overshoot,thickness,Jdet; IssmDouble D_scalar; IssmDouble *xyz_list = NULL; /*Get problem dimension*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DverticalEnum: dim = 1; bsize = 1; break; case Domain2DhorizontalEnum: dim = 2; bsize = 3; break; case Domain3DEnum: dim = 2; bsize = 3; break; default: _error_("mesh "<GetNumberOfNodes(); int numdof = numnodes*dim; /*Initialize Element matrix and vectors*/ ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum); IssmDouble* B = xNew(bsize*numdof); IssmDouble* Bprime = xNew(bsize*numdof); IssmDouble* D = xNewZeroInit(bsize*bsize); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input); Input* vy_input = NULL; Input* vyold_input = NULL; if(dim==2){ vy_input = element->GetInput(VyEnum); _assert_(vy_input); vyold_input = element->GetInput(VyPicardEnum); _assert_(vyold_input); } element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); /* Start looping on the number of gaussian points: */ Gauss* gauss = element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); this->GetBSSA(B,element,dim,xyz_list,gauss); this->GetBSSAprime(Bprime,element,dim,xyz_list,gauss); element->material->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); element->material->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); thickness_input->GetInputValue(&thickness, gauss); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet; for(int i=0;ivalues[0],1); } /*Transform Coordinate System*/ if(dim==2) element->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(D); xDelete(Bprime); xDelete(B); return Ke; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries*/ int domaintype; Element* basalelement; /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DhorizontalEnum: basalelement = element; break; case Domain3DEnum: case Domain2DverticalEnum: if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; default: _error_("mesh "<DeleteMaterials(); delete basalelement;}; delete pe1; delete pe2; return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorSSADrivingStress(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries */ int dim,domaintype; IssmDouble thickness,Jdet,slope[2]; IssmDouble* xyz_list = NULL; /*Get problem dimension*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DverticalEnum: dim = 1;break; case Domain2DhorizontalEnum: dim = 2;break; case Domain3DEnum: dim = 2;break; default: _error_("mesh "<GetNumberOfNodes(); /*Initialize Element vector and vectors*/ ElementVector* pe = element->NewElementVector(SSAApproximationEnum); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis, gauss); thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0); surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); for(int i=0;ivalues[i*dim+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; if(dim==2) pe->values[i*dim+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; } } /*Transform coordinate system*/ if(dim==2) element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ xDelete(xyz_list); xDelete(basis); delete gauss; return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*If no front, return NULL*/ if(!element->IsIcefront()) return NULL; /*Intermediaries*/ int dim,domaintype; IssmDouble Jdet,thickness,base,sealevel,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; IssmDouble *xyz_list = NULL; IssmDouble *xyz_list_front = NULL; IssmDouble normal[2]; /*Get problem dimension*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DverticalEnum: dim = 1;break; case Domain2DhorizontalEnum: dim = 2;break; case Domain3DEnum: dim = 2;break; default: _error_("mesh "<GetNumberOfNodes(); /*Initialize Element vector and other vectors*/ ElementVector* pe = element->NewElementVector(SSAApproximationEnum); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); element->GetVerticesCoordinates(&xyz_list); element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); element->NormalSection(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); thickness_input->GetInputValue(&thickness,gauss); sealevel_input->GetInputValue(&sealevel,gauss); base_input->GetInputValue(&base,gauss); element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); element->NodalFunctions(basis,gauss); surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level base_under_water = min(0.,base-sealevel); // 0 if the bottom of the glacier is above water level water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water); ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness; pressure = ice_pressure + water_pressure; for(int i=0;ivalues[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; if(dim==2) pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; } } /*Transform coordinate system*/ if(dim==2) element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ xDelete(xyz_list); xDelete(xyz_list_front); xDelete(basis); delete gauss; return pe; }/*}}}*/ 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 * by: * 2D 1D * Bi=[ dN/dx 0 ] Bi=[ dN/dx ] * [ 0 dN/dy ] * [ 1/2*dN/dy 1/2*dN/dx ] * where N is the finiteelement function for node i. * * We assume B has been allocated already, of size: 3x(NDOF2*numnodes) */ /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Get nodal functions derivatives*/ IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B: */ if(dim==2){ for(int i=0;i(dbasis); }/*}}}*/ 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 * by: * 2D 1D * 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==2){ for(int i=0;i(basis); }/*}}}*/ 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 * by: * 2D 1D * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] * [ dN/dx 2*dN/dy ] * [ dN/dy dN/dx ] * where hNis the finiteelement function for node i. * * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes) */ /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Get nodal functions derivatives*/ IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B': */ if(dim==2){ for(int i=0;i(dbasis); }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/ int i,dim,domaintype; IssmDouble rho_ice,g; int* doflist=NULL; IssmDouble* xyz_list=NULL; Element* basalelement=NULL; /*Deal with pressure first*/ int numvertices = element->GetNumberOfVertices(); IssmDouble* pressure = xNew(numvertices); IssmDouble* thickness = xNew(numvertices); IssmDouble* surface = xNew(numvertices); element->FindParam(&domaintype,DomainTypeEnum); rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); g =element->GetMaterialParameter(ConstantsGEnum); switch(domaintype){ case Domain2DhorizontalEnum: element->GetInputListOnVertices(thickness,ThicknessEnum); for(i=0;iGetVerticesCoordinates(&xyz_list); element->GetInputListOnVertices(surface,SurfaceEnum); for(i=0;iGetVerticesCoordinates(&xyz_list); element->GetInputListOnVertices(surface,SurfaceEnum); for(i=0;iAddInput(PressureEnum,pressure,P1Enum); xDelete(pressure); xDelete(thickness); xDelete(surface); /*Get basal element*/ switch(domaintype){ case Domain2DhorizontalEnum: basalelement = element; break; case Domain3DEnum: case Domain2DverticalEnum: if(!element->IsOnBase()){xDelete(xyz_list); return;} basalelement=element->SpawnBasalElement(); break; default: _error_("mesh "<GetNumberOfNodes(); int numdof = numnodes*dim; /*Fetch dof list and allocate solution vectors*/ basalelement->GetDofList(&doflist,SSAApproximationEnum,GsetEnum); IssmDouble* values = xNew(numdof); IssmDouble* vx = xNew(numnodes); IssmDouble* vy = xNew(numnodes); IssmDouble* vz = xNew(numnodes); IssmDouble* vel = xNew(numnodes); /*Use the dof list to index into the solution vector: */ for(i=0;iTransformSolutionCoord(&values[0],XYEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[i])) _error_("Inf found in solution vector"); if(dim==2){ vy[i]=values[i*dim+1]; if(xIsNan(vy[i])) _error_("NaN found in solution vector"); if(xIsInf(vy[i])) _error_("Inf found in solution vector"); } } /*Get Vz and compute vel*/ if(dim==2){ basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); for(i=0;iGetInputListOnNodes(&vy[0],VyEnum,0.); for(i=0;iInputChangeName(VxEnum,VxPicardEnum); if(dim==2)element->InputChangeName(VyEnum,VyPicardEnum); /*Add vx and vy as inputs to the tria element: */ element->AddBasalInput(VxEnum,vx,element->GetElementType()); if(dim==2)element->AddBasalInput(VyEnum,vy,element->GetElementType()); element->AddBasalInput(VelEnum,vel,element->GetElementType()); /*Free ressources:*/ xDelete(vel); xDelete(vz); xDelete(vy); xDelete(vx); xDelete(values); xDelete(xyz_list); xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ /*L1L2*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*compute all stiffness matrices for this element*/ ElementMatrix* Ke1=CreateKMatrixL1L2Viscous(element); ElementMatrix* Ke2=CreateKMatrixL1L2Friction(element); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); /*clean-up and return*/ delete Ke1; delete Ke2; return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2Friction(Element* element){/*{{{*/ if(!element->IsOnBase() || element->IsFloating()) return NULL; Element* basalelement = element->SpawnBasalElement(); ElementMatrix* Ke = CreateKMatrixSSAFriction(basalelement); /*clean-up and return*/ basalelement->DeleteMaterials(); delete basalelement; return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2Viscous(Element* element){/*{{{*/ /*Intermediaries*/ IssmDouble viscosity,Jdet; IssmDouble *xyz_list = NULL; /*Get element on base*/ Element* basalelement = element->GetBasalElement()->SpawnBasalElement(); /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); int numdof = numnodes*2; /*Initialize Element matrix and vectors*/ ElementMatrix* Ke = basalelement->NewElementMatrix(L1L2ApproximationEnum); IssmDouble* B = xNew(3*numdof); IssmDouble* Bprime = xNew(3*numdof); IssmDouble* D = xNewZeroInit(3*3); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); /* Start looping on the number of gaussian points: */ Gauss* gauss = element->NewGauss(5); Gauss* gauss_base = basalelement->NewGauss(); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); gauss->SynchronizeGaussBase(gauss_base); element->JacobianDeterminant(&Jdet,xyz_list,gauss); this->GetBSSA(B,basalelement,2,xyz_list,gauss_base); this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base); element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet; TripleMultiply(B,3,numdof,1, D,3,3,0, Bprime,3,numdof,0, &Ke->values[0],1); } /*Transform Coordinate System*/ basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ delete gauss; delete gauss_base; basalelement->DeleteMaterials(); delete basalelement; xDelete(xyz_list); xDelete(D); xDelete(Bprime); xDelete(B); return Ke; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorL1L2(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries*/ int domaintype; Element* basalelement; /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DhorizontalEnum: basalelement = element; break; case Domain3DEnum: case Domain2DverticalEnum: if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; default: _error_("mesh "<DeleteMaterials(); delete basalelement;}; delete pe1; delete pe2; return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorL1L2DrivingStress(Element* element){/*{{{*/ /*Intermediaries */ IssmDouble thickness,Jdet,slope[2]; IssmDouble* xyz_list = NULL; /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Initialize Element vector and vectors*/ ElementVector* pe = element->NewElementVector(L1L2ApproximationEnum); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis, gauss); thickness_input->GetInputValue(&thickness,gauss); surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); for(int i=0;ivalues[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; } } /*Transform coordinate system*/ element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ xDelete(xyz_list); xDelete(basis); delete gauss; return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorL1L2Front(Element* element){/*{{{*/ /*If no front, return NULL*/ if(!element->IsIcefront()) return NULL; /*Intermediaries*/ IssmDouble Jdet,thickness,bed,sealevel,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; IssmDouble *xyz_list = NULL; IssmDouble *xyz_list_front = NULL; IssmDouble normal[2]; /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Initialize Element vector and other vectors*/ ElementVector* pe = element->NewElementVector(L1L2ApproximationEnum); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); element->GetVerticesCoordinates(&xyz_list); element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); element->NormalSection(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); thickness_input->GetInputValue(&thickness,gauss); base_input->GetInputValue(&bed,gauss); sealevel_input->GetInputValue(&sealevel,gauss); element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); element->NodalFunctions(basis,gauss); surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level base_under_water = min(0.,bed-sealevel); // 0 if the bottom of the glacier is above water level water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water); ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness; pressure = ice_pressure + water_pressure; for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; } } /*Transform coordinate system*/ element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ xDelete(xyz_list); xDelete(xyz_list_front); xDelete(basis); delete gauss; return pe; }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/ int i,dim,domaintype; IssmDouble rho_ice,g; int* doflist=NULL; IssmDouble* xyz_list=NULL; Element* basalelement=NULL; /*Deal with pressure first*/ int numvertices = element->GetNumberOfVertices(); IssmDouble* pressure = xNew(numvertices); IssmDouble* thickness = xNew(numvertices); IssmDouble* surface = xNew(numvertices); element->FindParam(&dim,DomainDimensionEnum); element->FindParam(&domaintype,DomainTypeEnum); rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); g =element->GetMaterialParameter(ConstantsGEnum); if(dim==2){ element->GetInputListOnVertices(thickness,ThicknessEnum); for(i=0;iGetVerticesCoordinates(&xyz_list); element->GetInputListOnVertices(surface,SurfaceEnum); for(i=0;iAddInput(PressureEnum,pressure,P1Enum); xDelete(pressure); xDelete(thickness); xDelete(surface); /*Get basal element*/ switch(domaintype){ case Domain2DhorizontalEnum: basalelement = element; break; case Domain3DEnum: if(!element->IsOnBase()){xDelete(xyz_list); return;} basalelement=element->SpawnBasalElement(); break; default: _error_("mesh "<GetNumberOfNodes(); int numdof = numnodes*2; /*Fetch dof list and allocate solution vectors*/ basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); IssmDouble* values = xNew(numdof); IssmDouble* vx = xNew(numnodes); IssmDouble* vy = xNew(numnodes); IssmDouble* vz = xNew(numnodes); IssmDouble* vel = xNew(numnodes); /*Use the dof list to index into the solution vector: */ for(i=0;iTransformSolutionCoord(&values[0],XYEnum); basalelement->FindParam(&domaintype,DomainTypeEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[i])) _error_("Inf found in solution vector"); if(xIsNan(vy[i])) _error_("NaN found in solution vector"); if(xIsInf(vy[i])) _error_("Inf found in solution vector"); } /*Get Vz and compute vel*/ basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); for(i=0;iInputChangeName(VxEnum,VxPicardEnum); element->InputChangeName(VyEnum,VyPicardEnum); /*Add vx and vy as inputs to the tria element: */ element->AddBasalInput(VxEnum,vx,element->GetElementType()); element->AddBasalInput(VyEnum,vy,element->GetElementType()); element->AddBasalInput(VelEnum,vel,element->GetElementType()); /*Free ressources:*/ xDelete(vel); xDelete(vz); xDelete(vy); xDelete(vx); xDelete(values); xDelete(xyz_list); xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ /*HO*/ ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixHO(Element* element){/*{{{*/ /*Intermediaries */ IssmDouble Jdet; IssmDouble eps1dotdphii,eps1dotdphij; IssmDouble eps2dotdphii,eps2dotdphij; IssmDouble mu_prime; IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ IssmDouble eps1[3],eps2[3]; IssmDouble *xyz_list = NULL; /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Initialize Element matrix, vectors and Gaussian points*/ ElementMatrix* Ke=this->CreateKMatrixHO(element); //Initialize Jacobian with regular HO (first part of the Gateau derivative) IssmDouble* dbasis = xNew(3*numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = element->GetInput(VyEnum); _assert_(vy_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->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); element->material->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; for(int i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*eps1dotdphij*eps1dotdphii; Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*eps2dotdphij*eps1dotdphii; Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*eps1dotdphij*eps2dotdphii; Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*eps2dotdphij*eps2dotdphii; } } } /*Transform Coordinate System*/ element->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(dbasis); return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixHO(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*compute all stiffness matrices for this element*/ ElementMatrix* Ke1=CreateKMatrixHOViscous(element); ElementMatrix* Ke2=CreateKMatrixHOFriction(element); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); /*clean-up and return*/ delete Ke1; delete Ke2; 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,dim==3?2:1); /*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); } 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){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries*/ int dim,bsize; IssmDouble viscosity,newviscosity,oldviscosity; IssmDouble viscosity_overshoot,thickness,Jdet; IssmDouble D_scalar; IssmDouble *xyz_list = NULL; /*Get problem dimension*/ element->FindParam(&dim,DomainDimensionEnum); if(dim==2) bsize = 2; else bsize = 5; /*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(bsize*numdof); IssmDouble* Bprime = xNew(bsize*numdof); IssmDouble* D = xNewZeroInit(bsize*bsize); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); Input* vxold_input = element->GetInput(VxPicardEnum); _assert_(vxold_input); Input* vy_input = NULL; Input* vyold_input = NULL; if(dim==3){ vy_input=element->GetInput(VyEnum); _assert_(vy_input); vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); } element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); /* 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->GetBHO(B,element,dim,xyz_list,gauss); this->GetBHOprime(Bprime,element,dim,xyz_list,gauss); element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=2.*newviscosity*gauss->weight*Jdet; for(int i=0;ivalues[0],1); } /*Transform Coordinate System*/ if(dim==3) element->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(D); xDelete(Bprime); xDelete(B); return Ke; }/*}}}*/ #ifdef FSANALYTICAL ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries */ int dim; IssmDouble x_coord,y_coord,z_coord; IssmDouble Jdet,forcex,forcey,forcez; IssmDouble* xyz_list = NULL; /*Get problem dimension*/ element->FindParam(&dim,DomainDimensionEnum); /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Initialize Element vector and vectors*/ ElementVector* pe=element->NewElementVector(HOApproximationEnum); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(3); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis, gauss); x_coord=element->GetXcoord(xyz_list,gauss); y_coord=element->GetYcoord(xyz_list,gauss); if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss); else z_coord=0.; forcex=fx(x_coord,y_coord,z_coord,FSANALYTICAL); forcey=fy(x_coord,y_coord,z_coord,FSANALYTICAL); for(int i=0;ivalues[i*(dim-1)+0]+=forcex*Jdet*gauss->weight*basis[i]; pe->values[i*(dim-1)+1]+=forcey*Jdet*gauss->weight*basis[i]; } } /*Transform coordinate system*/ if(dim==3) element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ xDelete(basis); xDelete(xyz_list); delete gauss; return pe; }/*}}}*/ #else ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*compute all load vectors for this element*/ ElementVector* pe1=CreatePVectorHODrivingStress(element); ElementVector* pe2=CreatePVectorHOFront(element); ElementVector* pe =new ElementVector(pe1,pe2); /*clean-up and return*/ delete pe1; delete pe2; return pe; }/*}}}*/ #endif ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Intermediaries */ int dim; IssmDouble Jdet,slope[3]; IssmDouble* xyz_list = NULL; /*Get problem dimension*/ element->FindParam(&dim,DomainDimensionEnum); /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Initialize Element vector and vectors*/ ElementVector* pe=element->NewElementVector(HOApproximationEnum); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(3); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis, gauss); surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); for(int i=0;ivalues[i*(dim-1)+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i]; if(dim==3) pe->values[i*(dim-1)+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i]; } } /*Transform coordinate system*/ if(dim==3) element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ xDelete(basis); xDelete(xyz_list); delete gauss; return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*If no front, return NULL*/ if(!element->IsIcefront()) return NULL; /*Intermediaries*/ int dim; IssmDouble Jdet,surface,sealevel,z,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; IssmDouble* xyz_list = NULL; IssmDouble* xyz_list_front = NULL; IssmDouble normal[3]; 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 numvertices = element->GetNumberOfVertices(); /*Initialize Element vector and other vectors*/ ElementVector* pe = element->NewElementVector(HOApproximationEnum); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); element->GetVerticesCoordinates(&xyz_list); element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); element->NormalSection(&normal[0],xyz_list_front); /*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,10);//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); surface_input->GetInputValue(&surface,gauss); sealevel_input->GetInputValue(&sealevel,gauss); if(dim==3) z=element->GetZcoord(xyz_list,gauss); else z=element->GetYcoord(xyz_list,gauss); element->NodalFunctions(basis,gauss); element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); water_pressure = rho_water*gravity*min(0.,z-sealevel);//0 if the gaussian point is above water level ice_pressure = rho_ice*gravity*(surface-z); pressure = ice_pressure + water_pressure; for (int i=0;ivalues[(dim-1)*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; if(dim==3) pe->values[(dim-1)*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; } } /*Transform coordinate system*/ if(dim==3)element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ xDelete(basis); xDelete(xyz_list); xDelete(xyz_list_front); delete gauss; return pe; }/*}}}*/ 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 * by: * 3D 2D * * Bi=[ dh/dx 0 ] Bi=[ dh/dx] * [ 0 dh/dy ] [ dh/dy] * [ 1/2*dh/dy 1/2*dh/dx ] * [ 1/2*dh/dz 0 ] * [ 0 1/2*dh/dz ] * where h is the interpolation function for node i. * * We assume B has been allocated already, of size: 5x(NDOF2*numnodes) */ /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Get nodal functions derivatives*/ IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B: */ if(dim==2){ for(int i=0;i(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::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 * by: * 3D 2D * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] * [ dN/dx 2*dN/dy ] [ dN/dy ] * [ dN/dy dN/dx ] * where hNis the finiteelement function for node i. * * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes) */ /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Get nodal functions derivatives*/ IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B': */ if(dim==3){ for(int i=0;i(dbasis); }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/ int i,domaintype,dim; int* doflist=NULL; IssmDouble* xyz_list=NULL; /*Deal with pressure first*/ int numvertices = element->GetNumberOfVertices(); IssmDouble* pressure = xNew(numvertices); IssmDouble* surface = xNew(numvertices); element->FindParam(&domaintype,DomainTypeEnum); IssmDouble rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble g =element->GetMaterialParameter(ConstantsGEnum); switch(domaintype){ case Domain3DEnum: element->GetVerticesCoordinates(&xyz_list); element->GetInputListOnVertices(surface,SurfaceEnum); for(i=0;iGetVerticesCoordinates(&xyz_list); element->GetInputListOnVertices(surface,SurfaceEnum); for(i=0;iAddInput(PressureEnum,pressure,P1Enum); xDelete(pressure); xDelete(surface); /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); int numdof = numnodes*(dim-1); /*Fetch dof list and allocate solution vectors*/ element->GetDofList(&doflist,HOApproximationEnum,GsetEnum); IssmDouble* values = xNew(numdof); IssmDouble* vx = xNew(numnodes); IssmDouble* vy = xNew(numnodes); IssmDouble* vz = xNew(numnodes); IssmDouble* vel = xNew(numnodes); /*Use the dof list to index into the solution vector: */ for(i=0;iTransformSolutionCoord(&values[0],XYEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[i])) _error_("Inf found in solution vector"); if(dim==3){ vy[i]=values[i*(dim-1)+1]; if(xIsNan(vy[i])) _error_("NaN found in solution vector"); if(xIsInf(vy[i])) _error_("Inf found in solution vector"); } } /*Get Vz and compute vel*/ if(dim==3){ element->GetInputListOnNodes(&vz[0],VzEnum,0.); for(i=0;iGetInputListOnNodes(&vy[0],VyEnum,0.); for(i=0;iInputChangeName(VxEnum,VxPicardEnum); if(dim==3)element->InputChangeName(VyEnum,VyPicardEnum); /*Add vx and vy as inputs to the element: */ element->AddInput(VxEnum,vx,element->GetElementType()); if(dim==3)element->AddInput(VyEnum,vy,element->GetElementType()); element->AddInput(VelEnum,vel,element->GetElementType()); /*Free ressources:*/ xDelete(vel); xDelete(vz); xDelete(vy); xDelete(vx); xDelete(values); xDelete(xyz_list); xDelete(doflist); }/*}}}*/ /*FS*/ ElementVector* StressbalanceAnalysis::CreateDVectorFS(Element* element){/*{{{*/ int dim; /*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(); /*Initialize output vector*/ ElementVector* de = element->NewElementVector(FSvelocityEnum); for(int i=0;ivalues[i*dim+j]=VelocityEnum; } for(int i=0;ivalues[vnumnodes*dim+i]=PressureEnum; } return de; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixFS(Element* element){/*{{{*/ /*Intermediaries */ int i,j; IssmDouble Jdet; IssmDouble eps1dotdphii,eps1dotdphij; IssmDouble eps2dotdphii,eps2dotdphij; IssmDouble eps3dotdphii,eps3dotdphij; IssmDouble mu_prime; IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ IssmDouble eps1[3],eps2[3],eps3[3]; IssmDouble *xyz_list = NULL; /*Fetch number of nodes and dof for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; /*Prepare coordinate system list*/ int* cs_list = xNew(vnumnodes+pnumnodes); for(i=0;iCreateKMatrixFS(element); //Initialize Jacobian with regular FS (first part of the Gateau derivative) IssmDouble* dbasis = xNew(3*vnumnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); 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); /* 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->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss); //element->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); //eps1[0]=epsilon[0]; eps2[0]=epsilon[3]; eps3[0]=epsilon[4]; //eps1[1]=epsilon[3]; eps2[1]=epsilon[1]; eps3[1]=epsilon[5]; //eps1[2]=epsilon[4]; eps2[2]=epsilon[5]; eps3[2]=epsilon[2]; element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; element->material->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]); for(i=0;ivalues[numdof*(3*i+0)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; Ke->values[numdof*(3*i+0)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; Ke->values[numdof*(3*i+0)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; Ke->values[numdof*(3*i+1)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; Ke->values[numdof*(3*i+1)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; Ke->values[numdof*(3*i+1)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; Ke->values[numdof*(3*i+2)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; Ke->values[numdof*(3*i+2)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; Ke->values[numdof*(3*i+2)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; } } } /*Transform Coordinate System*/ element->TransformStiffnessMatrixCoord(Ke,cs_list); /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(dbasis); xDelete(cs_list); return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixFS(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; /*Get type of algorithm*/ int fe_FS; element->FindParam(&fe_FS,FlowequationFeFSEnum); /*compute all stiffness matrices for this element*/ ElementMatrix* Ke1=NULL; if(fe_FS==XTaylorHoodEnum) Ke1=CreateKMatrixFSViscousXTH(element); else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum) Ke1=CreateKMatrixFSViscousLA(element); else Ke1=CreateKMatrixFSViscous(element); ElementMatrix* Ke2=CreateKMatrixFSFriction(element); ElementMatrix* Ke3=CreateKMatrixFSShelf(element); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); /*clean-up and return*/ delete Ke1; delete Ke2; delete Ke3; 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* base_input = element->GetInput(BaseEnum); _assert_(base_input); /* Start looping on the number of gaussian points: */ gauss=element->NewGaussBase(3); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); base_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->material->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){/*{{{*/ /*Intermediaries*/ int i,dim,epssize; IssmDouble r,rl,Jdet,viscosity,DU,DUl; IssmDouble normal[3]; IssmDouble *xyz_list = NULL; IssmDouble *xyz_list_base = NULL; /*Get problem dimension*/ element->FindParam(&dim,DomainDimensionEnum); element->FindParam(&r,AugmentedLagrangianREnum); 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->GetNumberOfNodes(P1Enum); int lnumnodes = element->GetNumberOfNodes(P2Enum); int numdof = vnumnodes*dim; int pnumdof = pnumnodes; int lnumdof = lnumnodes; /*Prepare coordinate system list*/ int* cs_list = xNew(vnumnodes); if(dim==2) for(i=0;iNewElementMatrix(FSvelocityEnum); IssmDouble* B = xNew(epssize*numdof); IssmDouble* Bprime = xNew(epssize*numdof); IssmDouble* BtBUzawa = xNewZeroInit(numdof*pnumdof); IssmDouble* BU = xNew(pnumdof); IssmDouble* BprimeU = xNew(numdof); IssmDouble* D = xNewZeroInit(epssize*epssize); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); Input* vz_input = NULL; 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->GetBFSvel(B,element,dim,xyz_list,gauss); this->GetBFSprimevel(Bprime,element,dim,xyz_list,gauss); element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); for(i=0;iweight*Jdet; TripleMultiply(B,epssize,numdof,1, D,epssize,epssize,0, Bprime,epssize,numdof,0, &Ke->values[0],1); this->GetBFSUzawa(BU,element,dim,xyz_list,gauss); this->GetBFSprimeUzawa(BprimeU,element,dim,xyz_list,gauss); DU = gauss->weight*Jdet*sqrt(r); TripleMultiply(BU,1,pnumdof,1, &DU,1,1,0, BprimeU,1,numdof,0, BtBUzawa,1); } /*The pressure augmentation should not be transformed*/ MatrixMultiply(BtBUzawa,pnumdof,numdof,1, BtBUzawa,pnumdof,numdof,0, &Ke->values[0],1); if(element->IsOnBase() && 0){ element->FindParam(&rl,AugmentedLagrangianRlambdaEnum); element->GetVerticesCoordinatesBase(&xyz_list_base); element->NormalBase(&normal[0],xyz_list_base); IssmDouble* Dlambda = xNewZeroInit(dim*dim); IssmDouble* C = xNewZeroInit(dim*lnumdof); IssmDouble* Cprime = xNewZeroInit(dim*numdof); IssmDouble* CtCUzawa = xNewZeroInit(numdof*lnumdof); delete gauss; gauss = element->NewGaussBase(5); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); this->GetCFS(C,element,dim,xyz_list,gauss); this->GetCFSprime(Cprime,element,dim,xyz_list,gauss); for(i=0;iweight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl); TripleMultiply(C,dim,lnumdof,1, Dlambda,dim,dim,0, Cprime,dim,numdof,0, CtCUzawa,1); } /*The sigma naugmentation should not be transformed*/ MatrixMultiply(CtCUzawa,lnumdof,numdof,1, CtCUzawa,lnumdof,numdof,0, &Ke->values[0],1); /*Delete base part*/ xDelete(Dlambda); xDelete(C); xDelete(Cprime); xDelete(CtCUzawa); xDelete(xyz_list_base); } /*Transform Coordinate System*/ element->TransformStiffnessMatrixCoord(Ke,cs_list); /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(D); xDelete(Bprime); xDelete(B); xDelete(BprimeU); xDelete(BU); xDelete(BtBUzawa); xDelete(cs_list); return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousXTH(Element* element){/*{{{*/ /*Intermediaries*/ int i,dim,epssize; IssmDouble r,FSreconditioning,Jdet; IssmDouble *xyz_list = NULL; /*Get problem dimension*/ element->FindParam(&dim,DomainDimensionEnum); element->FindParam(&r,AugmentedLagrangianREnum); 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); 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; }/*}}}*/ #ifdef FSANALYTICAL ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/ if(element->IsFloating() || !element->IsOnBase()) return NULL; /*If on water or not FS, skip stiffness: */ int approximation; element->GetInputValue(&approximation,ApproximationEnum); if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; /*Intermediaries*/ int i,dim; IssmDouble alpha2,Jdet; IssmDouble x_coord,y_coord,z_coord; 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* B = xNew(dim*numdof); IssmDouble* D = xNewZeroInit(dim*dim); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinatesBase(&xyz_list_base); element->GetVerticesCoordinates(&xyz_list); Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); Input* vz_input = NULL; if(dim==3){ vz_input = element->GetInput(VzEnum); _assert_(vz_input);} /* Start looping on the number of gaussian points: */ gauss=element->NewGaussBase(10); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); x_coord=element->GetXcoord(xyz_list,gauss); y_coord=element->GetYcoord(xyz_list,gauss); if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss); else z_coord=0.; alpha2=alpha(x_coord,y_coord,z_coord,FSANALYTICAL); this->GetBFSFriction(B,element,dim,xyz_list_base,gauss); element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); for(int i=0;iweight*Jdet; //taub_x = -alpha2 v_x (same for y) TripleMultiply(B,dim,numdof,1, D,dim,dim,0, B,dim,numdof,0, &Ke->values[0],1); } /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ /*Clean up and return*/ delete gauss; xDelete(xyz_list); xDelete(xyz_list_base); xDelete(B); xDelete(D); return Ke; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; ElementVector* pe = NULL; ElementVector* pe1=CreatePVectorFSViscous(element); ElementVector* pe2=CreatePVectorFSFriction(element); ElementVector* pe3=CreatePVectorFSStress(element); pe =new ElementVector(pe1,pe2,pe3); delete pe1; delete pe2; delete pe3; /*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; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/ int i,dim,fe_FS; IssmDouble x_coord,y_coord,z_coord; IssmDouble Jdet,forcex,forcey,forcez; IssmDouble *xyz_list = NULL; element->FindParam(&fe_FS,FlowequationFeFSEnum); 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->GetVerticesCoordinates(&xyz_list); /* 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->NodalFunctionsVelocity(vbasis,gauss); x_coord=element->GetXcoord(xyz_list,gauss); y_coord=element->GetYcoord(xyz_list,gauss); if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss); else z_coord=0.; forcex=fx(x_coord,y_coord,z_coord,FSANALYTICAL); forcey=fy(x_coord,y_coord,z_coord,FSANALYTICAL); forcez=fz(x_coord,y_coord,z_coord,FSANALYTICAL); for(i=0;ivalues[i*dim+0] += forcex *Jdet*gauss->weight*vbasis[i]; pe->values[i*dim+1] += forcey *Jdet*gauss->weight*vbasis[i]; if(dim==3) pe->values[i*dim+2] += forcez *Jdet*gauss->weight*vbasis[i]; } } /*Transform coordinate system*/ element->TransformLoadVectorCoord(pe,cs_list); /*Clean up and return*/ delete gauss; xDelete(cs_list); xDelete(vbasis); xDelete(xyz_list); if(fe_FS==XTaylorHoodEnum){ ElementVector* pe2=CreatePVectorFSViscousXTH(element); ElementVector* pe3 = new ElementVector(pe,pe2); delete pe; delete pe2; return pe3; } else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){ ElementVector* pe2=CreatePVectorFSViscousLA(element); ElementVector* pe3 = new ElementVector(pe,pe2); delete pe; delete pe2; return pe3; } return pe; }/*}}}*/ #else ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/ if(element->IsFloating() || !element->IsOnBase()) return NULL; /*If on water or not FS, skip stiffness: */ int approximation; element->GetInputValue(&approximation,ApproximationEnum); if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; /*Intermediaries*/ bool mainlyfloating; int dim,domaintype; 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); element->FindParam(&domaintype,DomainTypeEnum); /*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* B = xNew(dim*numdof); IssmDouble* D = xNewZeroInit(dim*dim); /*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,dim==3?3:1); /*Recover portion of element that is grounded*/ if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base); if(migration_style==SubelementMigration2Enum){ if(domaintype==Domain2DverticalEnum) _error_("Subelement Migration 2 not implemented yet for Flowline"); 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(3); } else{ gauss=element->NewGaussBase(3); } /* 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->GetBFSFriction(B,element,dim,xyz_list_base,gauss); element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); for(int i=0;iweight*Jdet; //taub_x = -alpha2 v_x (same for y) TripleMultiply(B,dim,numdof,1, D,dim,dim,0, B,dim,numdof,0, &Ke->values[0],1); } /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ /*Clean up and return*/ delete gauss; delete friction; xDelete(xyz_list_base); xDelete(B); xDelete(D); return Ke; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; ElementVector* pe = NULL; int fe_FS; element->FindParam(&fe_FS,FlowequationFeFSEnum); if(fe_FS==XTaylorHoodEnum){ ElementVector* pe1=CreatePVectorFSViscous(element); ElementVector* pe2=CreatePVectorFSShelf(element); ElementVector* pe3=CreatePVectorFSFront(element); ElementVector* petemp =new ElementVector(pe1,pe2,pe3); ElementVector* pe4=CreatePVectorFSViscousXTH(element); pe = new ElementVector(petemp,pe4); delete pe1; delete pe2; delete pe3; delete petemp; delete pe4; } else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){ ElementVector* pe1=CreatePVectorFSViscous(element); ElementVector* pe2=CreatePVectorFSShelf(element); ElementVector* pe3=CreatePVectorFSFront(element); ElementVector* petemp =new ElementVector(pe1,pe2,pe3); ElementVector* pe4=CreatePVectorFSViscousLA(element); pe = new ElementVector(petemp,pe4); delete pe1; delete pe2; delete pe3; delete petemp; delete pe4; } else{ ElementVector* pe1=CreatePVectorFSViscous(element); ElementVector* pe2=CreatePVectorFSShelf(element); ElementVector* pe3=CreatePVectorFSFront(element); pe =new ElementVector(pe1,pe2,pe3); delete pe1; delete pe2; delete pe3; } /*clean-up and return*/ return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/ int i,dim; IssmDouble Jdet,forcex,forcey,forcez; IssmDouble *xyz_list = 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->GetVerticesCoordinates(&xyz_list); IssmDouble rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); Input* loadingforcex_input=element->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); Input* loadingforcey_input=element->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); Input* loadingforcez_input=NULL; if(dim==3){ loadingforcez_input=element->GetInput(LoadingforceZEnum); _assert_(loadingforcez_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->NodalFunctionsVelocity(vbasis,gauss); loadingforcex_input->GetInputValue(&forcex,gauss); loadingforcey_input->GetInputValue(&forcey,gauss); if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss); for(i=0;ivalues[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; if(dim==3){ pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i]; pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; } else{ pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; } } } /*Transform coordinate system*/ element->TransformLoadVectorCoord(pe,cs_list); /*Clean up and return*/ delete gauss; xDelete(cs_list); xDelete(vbasis); xDelete(xyz_list); return pe; }/*}}}*/ #endif ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/ /*If no front, return NULL*/ if(!element->IsIcefront()) return NULL; /*Intermediaries*/ int i,dim; IssmDouble Jdet,pressure,surface,sealevel,z; IssmDouble normal[3]; IssmDouble *xyz_list = NULL; IssmDouble *xyz_list_front = 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 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); Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_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); sealevel_input->GetInputValue(&sealevel,gauss); if(dim==3) z=element->GetZcoord(xyz_list,gauss); else z=element->GetYcoord(xyz_list,gauss); pressure = rho_water*gravity*min(0.,z-sealevel);//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,base; 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); base_input->GetInputValue(&base, gauss); water_pressure=gravity*rho_water*base; for(i=0;ivalues[i*dim+(dim-1)]+=-water_pressure*gauss->weight*Jdet*vbasis[i]; } } /* 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; 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); mb_input->GetInputValue(&mb, gauss); for(i=0;ivalues[i*dim+(dim-1)] += -dt*rho_water*gravity*mb*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(cs_list); xDelete(vbasis); xDelete(xyz_list_base); return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLA(Element* element){/*{{{*/ int i,dim; IssmDouble Jdet,pressure; IssmDouble bed_normal[3]; IssmDouble *xyz_list = NULL; 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(); /*Prepare coordinate system list*/ int* cs_list = xNew(numnodes); if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); IssmDouble* dbasis = xNew(3*numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); /*Get pressure and sigmann*/ Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input); Input* sigmann_input =element->GetInput(SigmaNNEnum); _assert_(sigmann_input); gauss=element->NewGauss(5); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); pressure_input->GetInputValue(&pressure, gauss); element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss); for(i=0;ivalues[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i]; pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i]; if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i]; } } if(element->IsOnBase() && 0){ IssmDouble sigmann; IssmDouble* vbasis = xNew(numnodes); element->GetVerticesCoordinatesBase(&xyz_list_base); element->NormalBase(&bed_normal[0],xyz_list_base); delete 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); sigmann_input->GetInputValue(&sigmann, gauss); for(i=0;ivalues[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i]; pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i]; if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i]; } } xDelete(xyz_list_base); xDelete(vbasis); } /*Transform coordinate system*/ element->TransformLoadVectorCoord(pe,cs_list); /*Clean up and return*/ delete gauss; xDelete(cs_list); xDelete(xyz_list); xDelete(dbasis); return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/ int i,tausize,dim; IssmDouble Jdet,r; IssmDouble epsxx,epsyy,epszz,epsxy,epsxz,epsyz; IssmDouble sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz; IssmDouble *xyz_list = NULL; Gauss* gauss = NULL; /*Get problem dimension*/ element->FindParam(&dim,DomainDimensionEnum); if(dim==2) tausize = 3; else tausize = 6; /*Fetch number of nodes and dof for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG /*Prepare coordinate system list*/ int* cs_list = xNew(vnumnodes+pnumnodes); if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); IssmDouble* Dstar = xNewZeroInit((dim*vnumnodes)*(tausize*tnumnodes)); IssmDouble* tau = xNew(tausize*tnumnodes); IssmDouble* d = xNew(tausize*tnumnodes); IssmDouble* vdbasis = xNew(dim*vnumnodes); IssmDouble* tbasis = xNew(tnumnodes); IssmDouble* D = xNewZeroInit(tausize*tnumnodes*tausize*tnumnodes); /*Retrieve all inputs and parameters*/ element->FindParam(&r,AugmentedLagrangianREnum); element->GetVerticesCoordinates(&xyz_list); /*Get d and tau*/ Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input); Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input); Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input); Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL; Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input); Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input); Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input); Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL; if(dim==3){ epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input); epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input); epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input); sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input); sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input); sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input); } gauss = element->NewGauss(); for(int i=0;iGaussNode(P1DGEnum,i); epsxx_input->GetInputValue(&epsxx,gauss); sigmapxx_input->GetInputValue(&sigmapxx,gauss); epsyy_input->GetInputValue(&epsyy,gauss); sigmapyy_input->GetInputValue(&sigmapyy,gauss); epsxy_input->GetInputValue(&epsxy,gauss); sigmapxy_input->GetInputValue(&sigmapxy,gauss); if(dim==2){ d[i*tausize+0]=epsxx; tau[i*tausize+0]=sigmapxx; d[i*tausize+1]=epsyy; tau[i*tausize+1]=sigmapyy; d[i*tausize+2]=epsxy; tau[i*tausize+2]=sigmapxy; } else{ epszz_input->GetInputValue(&epszz,gauss); sigmapzz_input->GetInputValue(&sigmapzz,gauss); epsxz_input->GetInputValue(&epsxz,gauss); sigmapxz_input->GetInputValue(&sigmapxz,gauss); epsyz_input->GetInputValue(&epsyz,gauss); sigmapyz_input->GetInputValue(&sigmapyz,gauss); d[i*tausize+0]=epsxx; tau[i*tausize+0]=sigmapxx; d[i*tausize+1]=epsyy; tau[i*tausize+1]=sigmapyy; d[i*tausize+2]=epszz; tau[i*tausize+2]=sigmapzz; d[i*tausize+3]=epsxy; tau[i*tausize+3]=sigmapxy; d[i*tausize+4]=epsxz; tau[i*tausize+4]=sigmapxz; d[i*tausize+5]=epsyz; tau[i*tausize+5]=sigmapyz; } } /* Start looping on the number of gaussian points: */ delete gauss; gauss=element->NewGauss(5); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); /*Create Dstar*/ /*In dim = 2 * * <----------------- tausize ---------------> x tnumnodes * | gamma_ij^x 0 gamma_ij^y | ^ * Dij = | | dim * | 0 gamma_ij^y gamma_ij^x | v * x * vnumnodes * *In dim = 3 * * | gamma_ij^x 0 0 gamma_ij^y gamma_ij^z 0 | * | | * Dij = | 0 gamma_ij^y 0 gamma_ij^x 0 gamma_ij^z | * | | * | 0 0 gamma_ij^z 0 gamma_ij^x gamma_ij^y | * * gamma_ij^x = zeta_j dphi_i/dx * * where: * - zeta_j is the nodal function for the j^th node of the tensor (P1DG) * - phi_i is the nodal function for the i^th node of the velocity (P2)*/ element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); element->NodalFunctionsTensor(tbasis,gauss); if(dim==2){ for(int i=0;iweight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i]; Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i]; Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i]; Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i]; } } } else{ for(int i=0;iweight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i]; Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+3] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i]; Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+4] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i]; Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i]; Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+3] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i]; Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+5] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i]; Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i]; Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+4] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i]; Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+5] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i]; } } } } /*contribution -Dstar tau*/ for(i=0;ivalues[0],1); /*contribution + r Dstar d*/ for(i=0;ivalues[0],1); /*Transform coordinate system*/ element->TransformLoadVectorCoord(pe,cs_list); /*Clean up and return*/ delete gauss; xDelete(cs_list); xDelete(xyz_list); xDelete(Dstar); xDelete(d); xDelete(D); xDelete(tau); xDelete(vdbasis); xDelete(tbasis); return pe; }/*}}}*/ 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 * by: Bvi=[ dphi/dx 0 ] * [ 0 dphi/dy ] * [ 1/2*dphi/dy 1/2*dphi/dx] * [ 0 0 ] * [ dphi/dx dphi/dy ] * * Bpi=[ 0 ] * [ 0 ] * [ 0 ] * [ phi_p ] * [ 0 ] * * In 3d: * Bvi=[ dh/dx 0 0 ] * [ 0 dh/dy 0 ] * [ 0 0 dh/dz ] * [ 1/2*dh/dy 1/2*dh/dx 0 ] * [ 1/2*dh/dz 0 1/2*dh/dx ] * [ 0 1/2*dh/dz 1/2*dh/dy ] * [ 0 0 0 ] * [ dh/dx dh/dy dh/dz ] * * Bpi=[ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ h ] * [ 0 ] * where phi is the finiteelement function for node i. * Same thing for Bb except the last column that does not exist. */ /*Fetch number of nodes for this finite element*/ int pnumnodes = element->NumberofNodesPressure(); int vnumnodes = element->NumberofNodesVelocity(); /*Get nodal functions derivatives*/ IssmDouble* vdbasis=xNew(dim*vnumnodes); IssmDouble* pbasis =xNew(pnumnodes); element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); element->NodalFunctionsPressure(pbasis,gauss); /*Build B: */ if(dim==2){ for(int i=0;i(vdbasis); xDelete(pbasis); }/*}}}*/ 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 * by: * Bvi' = [ dphi/dx 0 ] * [ 0 dphi/dy ] * [ dphi/dy dphi/dx ] * [ dphi/dx dphi/dy ] * [ 0 0 ] * * by: Bpi=[ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ phi ] * * 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 ] * [ dh/dx dh/dy dh/dz ] * [ 0 0 0 ] * * Bpi=[ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ 0 ] * [ h ] * where phi is the finiteelement function for node i. * In 3d: */ /*Fetch number of nodes for this finite element*/ int pnumnodes = element->NumberofNodesPressure(); int vnumnodes = element->NumberofNodesVelocity(); /*Get nodal functions derivatives*/ IssmDouble* vdbasis=xNew(dim*vnumnodes); IssmDouble* pbasis =xNew(pnumnodes); element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); element->NodalFunctionsPressure(pbasis,gauss); /*Build B_prime: */ if(dim==2){ for(int i=0;i(vdbasis); xDelete(pbasis); }/*}}}*/ 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 * by: Bvi=[ dphi/dx 0 ] * [ 0 dphi/dy ] * [ 1/2*dphi/dy 1/2*dphi/dx] * * * In 3d: * Bvi=[ dh/dx 0 0 ] * [ 0 dh/dy 0 ] * [ 0 0 dh/dz ] * [ 1/2*dh/dy 1/2*dh/dx 0 ] * [ 1/2*dh/dz 0 1/2*dh/dx ] * [ 0 1/2*dh/dz 1/2*dh/dy ] * * where phi is the finiteelement function for node i. * Same thing for Bb except the last column that does not exist. */ /*Fetch number of nodes for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); /*Get nodal functions derivatives*/ IssmDouble* vdbasis=xNew(dim*vnumnodes); element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); /*Build B: */ if(dim==2){ for(int i=0;i(vdbasis); }/*}}}*/ void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute C matrix. C=[Cp1 Cp2 ...] where: * Cpi=[phi phi]. */ /*Fetch number of nodes for this finite element*/ int lnumnodes = element->GetNumberOfNodes(P2Enum); /*Get nodal functions derivatives*/ IssmDouble* basis =xNew(lnumnodes); element->NodalFunctionsP2(basis,gauss); /*Build B: */ for(int i=0;i(basis); }/*}}}*/ void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /* Compute C' matrix. C'=[C1' C2' ...] * Ci' = [ phi 0 ] * [ 0 phi ] * * In 3d * Ci' = [ phi 0 0 ] * [ 0 phi 0 ] * [ 0 0 phi ] * where phi is the finiteelement function for node i. */ /*Fetch number of nodes for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); int vnumdof = vnumnodes*dim; IssmDouble* vbasis=xNew(vnumnodes); element->NodalFunctionsVelocity(vbasis,gauss); /*Build B: */ if(dim==3){ for(int i=0;i(vbasis); }/*}}}*/ void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector* solution,Element* element){/*{{{*/ int* vdoflist=NULL; int* pdoflist=NULL; Input* vz_input=NULL; int dim; IssmDouble vx,vy,vz,p; IssmDouble FSreconditioning; /*Get some parameters*/ element->FindParam(&dim,DomainDimensionEnum); element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); /*Fetch number of nodes and dof for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); int vnumdof = vnumnodes*dim; int pnumdof = pnumnodes*1; /*Initialize values*/ IssmDouble* vvalues = xNew(vnumdof); IssmDouble* pvalues = xNew(pnumdof); /*Get dof list: */ element->GetDofListVelocity(&vdoflist,GsetEnum); element->GetDofListPressure(&pdoflist,GsetEnum); Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} Input* p_input =element->GetInput(PressureEnum); _assert_(p_input); element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); /*Ok, we have the velocities in inputs, fill in solution */ Gauss* gauss = element->NewGauss(); for(int i=0;iGaussNode(element->VelocityInterpolation(),i); vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); vvalues[i*dim+0]=vx; vvalues[i*dim+1]=vy; if(dim==3){ vz_input->GetInputValue(&vz,gauss); vvalues[i*dim+2]=vz; } } for(int i=0;iGaussNode(element->PressureInterpolation(),i); p_input->GetInputValue(&p ,gauss); pvalues[i]=p/FSreconditioning; } /*Add value to global vector*/ solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL); if(pnumdof>0) solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL); /*Free ressources:*/ delete gauss; xDelete(pdoflist); xDelete(vdoflist); xDelete(pvalues); xDelete(vvalues); }/*}}}*/ void StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/ /*Intermediaries*/ int dim; IssmDouble dvx[3],dvy[3],dvz[3]; IssmDouble viscosity; IssmDouble *xyz_list = NULL; /*Get problem dimension*/ parameters->FindParam(&dim,DomainDimensionEnum); for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); /*Get inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); 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);} /*Allocate new inputs*/ int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG IssmDouble* epsxx = xNew(tnumnodes); IssmDouble* sigmapxx = xNew(tnumnodes); IssmDouble* epsyy = xNew(tnumnodes); IssmDouble* sigmapyy = xNew(tnumnodes); IssmDouble* epsxy = xNew(tnumnodes); IssmDouble* sigmapxy = xNew(tnumnodes); IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL; IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL; IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL; if(dim==3){ epszz = xNew(tnumnodes); sigmapzz = xNew(tnumnodes); epsxz = xNew(tnumnodes); sigmapxz = xNew(tnumnodes); epsyz = xNew(tnumnodes); sigmapyz = xNew(tnumnodes); } /*Get d and tau*/ Gauss* gauss = element->NewGauss(); for(int i=0;iGaussNode(P1DGEnum,i); vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); if(dim==3){ vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); } element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); epsxx[i] = dvx[0]; sigmapxx[i] = 2.*viscosity*epsxx[i]; epsyy[i] = dvy[1]; sigmapyy[i] = 2.*viscosity*epsyy[i]; epsxy[i] = 0.5*(dvx[1] + dvy[0]); sigmapxy[i] = 2.*viscosity*epsxy[i]; if(dim==3){ epszz[i] = dvz[2]; sigmapzz[i] = 2.*viscosity*epszz[i]; epsxz[i] = 0.5*(dvx[2] + dvz[0]); sigmapxz[i] = 2.*viscosity*epsxz[i]; epsyz[i] = 0.5*(dvy[2] + dvz[1]); sigmapyz[i] = 2.*viscosity*epsyz[i]; } } /*Add inputs*/ element->AddInput(StrainRatexxEnum,epsxx,P1DGEnum); element->AddInput(DeviatoricStressxxEnum,sigmapxx,P1DGEnum); element->AddInput(StrainRateyyEnum,epsyy,P1DGEnum); element->AddInput(DeviatoricStressyyEnum,sigmapyy,P1DGEnum); element->AddInput(StrainRatexyEnum,epsxy,P1DGEnum); element->AddInput(DeviatoricStressxyEnum,sigmapxy,P1DGEnum); if(dim==3){ element->AddInput(StrainRatezzEnum,epszz,P1DGEnum); element->AddInput(DeviatoricStresszzEnum,sigmapzz,P1DGEnum); element->AddInput(StrainRatexzEnum,epsxz,P1DGEnum); element->AddInput(DeviatoricStressxzEnum,sigmapxz,P1DGEnum); element->AddInput(StrainRateyzEnum,epsyz,P1DGEnum); element->AddInput(DeviatoricStressyzEnum,sigmapyz,P1DGEnum); } /*Clean up*/ delete gauss; xDelete(xyz_list); xDelete(epsxx); xDelete(sigmapxx); xDelete(epsyy); xDelete(sigmapyy); xDelete(epszz); xDelete(sigmapzz); xDelete(epsxy); xDelete(sigmapxy); xDelete(epsxz); xDelete(sigmapxz); xDelete(epsyz); xDelete(sigmapyz); } }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/ bool results_on_nodes; int i,dim; int* vdoflist=NULL; int* pdoflist=NULL; IssmDouble FSreconditioning; element->FindParam(&dim,DomainDimensionEnum); element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); element->FindParam(&results_on_nodes,SettingsResultsOnNodesEnum); /*Fetch number of nodes and dof for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); int vnumdof = vnumnodes*dim; int pnumdof = pnumnodes*1; /*Initialize values*/ IssmDouble* values = xNew(vnumdof+pnumdof); IssmDouble* vx = xNew(vnumnodes); IssmDouble* vy = xNew(vnumnodes); IssmDouble* vz = xNew(vnumnodes); IssmDouble* vel = xNew(vnumnodes); IssmDouble* pressure = xNew(pnumnodes); /*Prepare coordinate system list*/ int* cs_list = xNew(vnumnodes+pnumnodes); if(dim==2) for(i=0;iGetDofListVelocity(&vdoflist,GsetEnum); element->GetDofListPressure(&pdoflist,GsetEnum); /*Use the dof list to index into the solution vector: */ for(i=0;iTransformSolutionCoord(values,cs_list); /*Ok, we have vx and vy in values, fill in all arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[i])) _error_("Inf found in solution vector"); if(xIsNan(vy[i])) _error_("NaN found in solution vector"); if(xIsInf(vy[i])) _error_("Inf found in solution vector"); if(dim==3){ vz[i] = values[i*dim+2]; if(xIsNan(vz[i])) _error_("NaN found in solution vector"); if(xIsInf(vz[i])) _error_("Inf found in solution vector"); } } for(i=0;i(pressure[i])) _error_("NaN found in solution vector"); if(xIsInf(pressure[i])) _error_("Inf found in solution vector"); } /*Recondition pressure and compute vel: */ for(i=0;iInputChangeName(VxEnum,VxPicardEnum); element->InputChangeName(VyEnum,VyPicardEnum); if(pnumdof>0) element->InputChangeName(PressureEnum,PressurePicardEnum); if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum); /*Add vx and vy as inputs to the tria element: */ element->AddInput(VxEnum, vx, element->VelocityInterpolation()); element->AddInput(VyEnum, vy, element->VelocityInterpolation()); element->AddInput(VelEnum, vel, element->VelocityInterpolation()); if(pnumdof>0) element->AddInput(PressureEnum,pressure,element->PressureInterpolation()); if(dim==3) element->AddInput(VzEnum,vz, element->VelocityInterpolation()); /*Free ressources:*/ xDelete(pressure); xDelete(vel); xDelete(vz); xDelete(vy); xDelete(vx); xDelete(values); xDelete(vdoflist); xDelete(pdoflist); xDelete(cs_list); }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/ /*Intermediaries*/ int dim,tausize; IssmDouble epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar; IssmDouble epsxx_old,epsyy_old,epszz_old,epsxy_old,epsxz_old,epsyz_old; IssmDouble sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz; IssmDouble dvx[3],dvy[3],dvz[3],B,n; IssmDouble *xyz_list = NULL; IssmDouble Jdet,r; parameters->FindParam(&r,AugmentedLagrangianREnum); parameters->FindParam(&dim,DomainDimensionEnum); if(dim==2) tausize = 3; else tausize = 6; for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); /*Get inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* B_input=element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input=element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); 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);} /*Fetch number of nodes and dof for this finite element*/ int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG /*Initialize vectors*/ IssmDouble* tbasis = xNew(tnumnodes); IssmDouble* Ke = xNewZeroInit(tnumnodes*tnumnodes); IssmDouble* pe_xx = xNewZeroInit(tnumnodes); IssmDouble* pe_yy = xNewZeroInit(tnumnodes); IssmDouble* pe_xy = xNewZeroInit(tnumnodes); IssmDouble* pe_zz = NULL; IssmDouble* pe_xz = NULL; IssmDouble* pe_yz = NULL; if(dim==3){ pe_zz = xNewZeroInit(tnumnodes); pe_xz = xNewZeroInit(tnumnodes); pe_yz = xNewZeroInit(tnumnodes); } /*Get previous d*/ Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input); Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input); Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input); Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL; if(dim==3){ epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input); epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input); epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input); } /*Get tau*/ Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input); Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input); Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input); Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL; if(dim==3){ sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input); sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input); sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input); } Gauss* gauss=element->NewGauss(5); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctionsTensor(tbasis,gauss); /*Get tau from inputs*/ sigmapxx_input->GetInputValue(&sigmapxx,gauss); sigmapyy_input->GetInputValue(&sigmapyy,gauss); sigmapxy_input->GetInputValue(&sigmapxy,gauss); if(dim==3){ sigmapzz_input->GetInputValue(&sigmapzz,gauss); sigmapxz_input->GetInputValue(&sigmapxz,gauss); sigmapyz_input->GetInputValue(&sigmapyz,gauss); } /*Get previous d*/ epsxx_input->GetInputValue(&epsxx_old,gauss); epsyy_input->GetInputValue(&epsyy_old,gauss); epsxy_input->GetInputValue(&epsxy_old,gauss); if(dim==3){ epszz_input->GetInputValue(&epszz_old,gauss); epsxz_input->GetInputValue(&epsxz_old,gauss); epsyz_input->GetInputValue(&epsyz_old,gauss); } /*Calculate d from previous results*/ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); if(dim==3){ vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); } epsxx = dvx[0]; epsyy = dvy[1]; epsxy = 0.5*(dvx[1] + dvy[0]); if(dim==3){ epszz = dvz[2]; epsxz = 0.5*(dvx[2] + dvz[0]); epsyz = 0.5*(dvy[2] + dvz[1]); } /*Solve 2 eta_0 |d|^s-1 + r |d| = |rD(u) + tau|*/ IssmDouble coef1,coef2,coef3; B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2) )^(n-1)/n ) coef2 = r; if(dim==2){ coef3 = sqrt( (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx) + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy) + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy) ); } else{ coef3 = sqrt( (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx) + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy) + (r*epszz + sigmapzz)*(r*epszz + sigmapzz) + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy) + 2*(r*epsxz + sigmapxz)*(r*epsxz + sigmapxz) + 2*(r*epsyz + sigmapyz)*(r*epsyz + sigmapyz) ); } IssmDouble dnorm; if(dim==2){ dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + 2.*epsxy_old*epsxy_old ); } else{ dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old +2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old)); } /*Initial guess cannot be 0 otherwise log(0) - inf*/ if(dnorm==0.) dnorm=1.; NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n,dnorm); _assert_(dnorm>=0.); _assert_(!xIsNan(dnorm)); /*Create Ke*/ D_scalar=(coef1*pow(dnorm,(1.-n)/n)+r)*gauss->weight*Jdet; TripleMultiply(tbasis,tnumnodes,1,0, &D_scalar,1,1,0, tbasis,1,tnumnodes,0, Ke,1); /*Create Right hand sides*/ for(int ii=0;iiweight*Jdet; for(int ii=0;iiweight*Jdet; for(int ii=0;iiweight*Jdet; if(dim==3){ for(int ii=0;iiweight*Jdet; for(int ii=0;iiweight*Jdet; for(int ii=0;iiweight*Jdet; } } /*Solve the systems*/ IssmDouble* d_xx = xNew(tnumnodes); IssmDouble* d_yy = xNew(tnumnodes); IssmDouble* d_xy = xNew(tnumnodes); IssmDouble* d_zz = NULL; IssmDouble* d_xz = NULL; IssmDouble* d_yz = NULL; if(dim==2){ _assert_(tnumnodes==3); Matrix3x3Solve(&d_xx[0],Ke,pe_xx); Matrix3x3Solve(&d_yy[0],Ke,pe_yy); Matrix3x3Solve(&d_xy[0],Ke,pe_xy); for(int i=0;i<3;i++) _assert_(!xIsNan(d_xx[i])); for(int i=0;i<3;i++) _assert_(!xIsNan(d_yy[i])); for(int i=0;i<3;i++) _assert_(!xIsNan(d_xx[i])); element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum); element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum); element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum); } else{ _assert_(tnumnodes==4); d_zz = xNew(tnumnodes); d_xz = xNew(tnumnodes); d_yz = xNew(tnumnodes); Matrix4x4Solve(&d_xx[0],Ke,pe_xx); Matrix4x4Solve(&d_yy[0],Ke,pe_yy); Matrix4x4Solve(&d_xy[0],Ke,pe_xy); Matrix4x4Solve(&d_zz[0],Ke,pe_zz); Matrix4x4Solve(&d_xz[0],Ke,pe_xz); Matrix4x4Solve(&d_yz[0],Ke,pe_yz); element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum); element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum); element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum); element->AddInput(StrainRatezzEnum,d_zz,P1DGEnum); element->AddInput(StrainRatexzEnum,d_xz,P1DGEnum); element->AddInput(StrainRateyzEnum,d_yz,P1DGEnum); } /*Clean up*/ delete gauss; xDelete(xyz_list); xDelete(tbasis); xDelete(Ke); xDelete(pe_xx); xDelete(d_xx); xDelete(pe_yy); xDelete(d_yy); xDelete(pe_zz); xDelete(d_zz); xDelete(pe_xy); xDelete(d_xy); xDelete(pe_xz); xDelete(d_xz); xDelete(pe_yz); xDelete(d_yz); } }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/ /*Intermediaries*/ int dim,tausize; IssmDouble epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar; IssmDouble d_xx,d_yy,d_zz,d_xy,d_xz,d_yz; IssmDouble sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz; IssmDouble dvx[3],dvy[3],dvz[3]; IssmDouble *xyz_list = NULL; IssmDouble Jdet,r; parameters->FindParam(&r,AugmentedLagrangianREnum); parameters->FindParam(&dim,DomainDimensionEnum); if(dim==2) tausize = 3; else tausize = 6; for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); /*Get inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); Input* vz_input=NULL; if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} /*Get previous tau*/ Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input); Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input); Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input); Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL; if(dim==3){ sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input); sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input); sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input); } /*Get NEW d*/ Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input); Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input); Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input); Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL; if(dim==3){ epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input); epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input); epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input); } /*Fetch number of nodes and dof for this finite element*/ int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG /*Update tau accordingly*/ IssmDouble* tau_xx = xNew(tnumnodes); IssmDouble* tau_yy = xNew(tnumnodes); IssmDouble* tau_xy = xNew(tnumnodes); IssmDouble* tau_zz = NULL; IssmDouble* tau_xz = NULL; IssmDouble* tau_yz = NULL; if(dim==3){ tau_zz = xNew(tnumnodes); tau_xz = xNew(tnumnodes); tau_yz = xNew(tnumnodes); } Gauss* gauss = element->NewGauss(); for(int ig=0;igGaussNode(P1DGEnum,ig); /*Get D(u)*/ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); if(dim==3){ vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); } epsxx = dvx[0]; epsyy = dvy[1]; epsxy = 0.5*(dvx[1] + dvy[0]); if(dim==3){ epszz = dvz[2]; epsxz = 0.5*(dvx[2] + dvz[0]); epsyz = 0.5*(dvy[2] + dvz[1]); } /*Get tau^(n-1) from inputs*/ sigmapxx_input->GetInputValue(&sigmapxx,gauss); sigmapyy_input->GetInputValue(&sigmapyy,gauss); sigmapxy_input->GetInputValue(&sigmapxy,gauss); if(dim==3){ sigmapzz_input->GetInputValue(&sigmapzz,gauss); sigmapxz_input->GetInputValue(&sigmapxz,gauss); sigmapyz_input->GetInputValue(&sigmapyz,gauss); } /*Get new d*/ epsxx_input->GetInputValue(&d_xx,gauss); epsyy_input->GetInputValue(&d_yy,gauss); epsxy_input->GetInputValue(&d_xy,gauss); if(dim==3){ epszz_input->GetInputValue(&d_zz,gauss); epsxz_input->GetInputValue(&d_xz,gauss); epsyz_input->GetInputValue(&d_yz,gauss); } /*Get d and update tau accordingly*/ tau_xx[ig] = sigmapxx + r*(epsxx - d_xx); tau_yy[ig] = sigmapyy + r*(epsyy - d_yy); tau_xy[ig] = sigmapxy + r*(epsxy - d_xy); if(dim==3){ tau_zz[ig] = sigmapzz + r*(epszz - d_zz); tau_xz[ig] = sigmapxz + r*(epsxz - d_xz); tau_yz[ig] = sigmapyz + r*(epsyz - d_yz); } } /*Add inputs*/ element->AddInput(DeviatoricStressxxEnum,tau_xx,P1DGEnum); element->AddInput(DeviatoricStressyyEnum,tau_yy,P1DGEnum); element->AddInput(DeviatoricStressxyEnum,tau_xy,P1DGEnum); if(dim==3){ element->AddInput(DeviatoricStresszzEnum,tau_zz,P1DGEnum); element->AddInput(DeviatoricStressxzEnum,tau_xz,P1DGEnum); element->AddInput(DeviatoricStressyzEnum,tau_yz,P1DGEnum); } /*Clean up and */ delete gauss; xDelete(xyz_list); xDelete(tau_xx); xDelete(tau_yy); xDelete(tau_zz); xDelete(tau_xy); xDelete(tau_xz); xDelete(tau_yz); } }/*}}}*/ /*Coupling (Tiling)*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingHOFS(Element* element){/*{{{*/ /*Constants*/ int numnodes = 3*6+1; int numdofp = 2*6; int numdofs = 4*6 + 3; int numdoftotal = (2+4)*6+ 3; /*Intermediaries*/ int i,j,init; int* cs_list = xNew(6*3+1); int* cs_list2 = xNew(6*2+1); Node **node_list = xNew(6*3+1); /*Some parameters needed*/ init = element->FiniteElement(); /*prepare node list*/ for(i=0;i<6+1;i++){ node_list[i+6] = element->GetNode(i); cs_list[i+6] = XYZEnum; cs_list2[i] = XYZEnum; } for(i=0;i<6;i++){ node_list[i] = element->GetNode(i); node_list[i+2*6+1] = element->GetNode(i+6*1); cs_list[i] = XYEnum; cs_list[i+2*6+1] = PressureEnum; cs_list2[i+6+1] = PressureEnum; } /*compute all stiffness matrices for this element*/ ElementMatrix* Ke1=element->NewElementMatrixCoupling(6,HOApproximationEnum); ElementMatrix* Ke2=element->NewElementMatrix(FSvelocityEnum); ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); delete Ke1; delete Ke2; /*Compute HO Matrix with P1 element type\n");*/ Ke1=CreateKMatrixFS(element); element->TransformInvStiffnessMatrixCoord(Ke1,node_list,2*6+1,cs_list2); int indices[3]={18,19,20}; Ke1->StaticCondensation(3,&indices[0]); element->SetTemporaryElementType(P1Enum); // P1 needed for HO Ke2=CreateKMatrixHO(element); element->TransformInvStiffnessMatrixCoord(Ke2,XYEnum); element->SetTemporaryElementType(init); // P1 needed for HO /*Compute FS Matrix and condense it \n");*/ for(i=0;ivalues[(i+numdofp)*numdoftotal+2*j+0]+=Ke1->values[i*numdofs+3*j+0]; Ke->values[(i+numdofp)*numdoftotal+2*j+1]+=Ke1->values[i*numdofs+3*j+1]; } for(i=0;ivalues[i*numdoftotal+numdofp+3*j+0]+=Ke2->values[i*numdofp+2*j+0]; Ke->values[i*numdoftotal+numdofp+3*j+1]+=Ke2->values[i*numdofp+2*j+1]; } /*Transform Coordinate System*/ //Do not transform, already done in the matrices element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); /*clean-up and return*/ xDelete(cs_list); xDelete(cs_list2); xDelete(node_list); delete Ke1; delete Ke2; return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFS(Element* element){/*{{{*/ /*compute all stiffness matrices for this element*/ ElementMatrix* Ke1=CreateKMatrixCouplingSSAFSViscous(element); ElementMatrix* Ke2=CreateKMatrixCouplingSSAFSFriction(element); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); /*clean-up and return*/ delete Ke1; delete Ke2; return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFSFriction(Element* element){/*{{{*/ /*Constants*/ const int numdofs = (6+1)*3 + 6*1; const int numdofm = 6 *2; const int numdof2d = 3 *3; const int numdof2dm = 3 *2; const int numdoftot = 6*2 + (6+1)*3 +6; // HO + FS vel + FS Pressure /*Intermediaries */ int i,j,approximation; int dim=3; IssmDouble FSreconditioning,viscosity,alpha2_gauss,Jdet2d; IssmDouble bed_normal[3]; IssmDouble LSSAFS[8][numdof2dm]; IssmDouble LprimeSSAFS[8][numdofs]; IssmDouble DLSSAFS[8][8]={0.0}; IssmDouble LFSSSA[4][numdof2d]; IssmDouble LprimeFSSSA[4][numdof2dm]; IssmDouble DLFSSSA[4][4]={0.0}; IssmDouble Ke_drag[numdof2dm][numdofs]={0.0}; IssmDouble Ke_drag2[numdof2d][numdof2dm]={0.0}; IssmDouble *xyz_list = NULL; IssmDouble *xyz_list_tria = NULL; /*If on water or not FS, skip stiffness: */ element->GetInputValue(&approximation,ApproximationEnum); if(element->IsFloating() || !element->IsOnBase()) return NULL; int vnumnodes = element->NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); int numnodes = 2*vnumnodes-1+pnumnodes; /*Prepare node list*/ int* cs_list = xNew(2*vnumnodes-1+pnumnodes); Node **node_list = xNew(2*vnumnodes-1+pnumnodes); for(i=0;iGetNode(i); cs_list[i] = XYEnum; } for(i=0;iGetNode(i); cs_list[i+vnumnodes-1] = XYZEnum; } for(i=0;iGetNode(vnumnodes+i); cs_list[2*vnumnodes-1+i] = PressureEnum; } ElementMatrix* Ke1=element->NewElementMatrixCoupling(6,SSAApproximationEnum); ElementMatrix* Ke2=element->NewElementMatrix(FSvelocityEnum); ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); delete Ke1; delete Ke2; /*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); /*build friction object, used later on: */ Friction* friction=new Friction(element,3); /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGaussBase(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss); this->GetLSSAFS(&LSSAFS[0][0], element,gauss); this->GetLprimeSSAFS(&LprimeSSAFS[0][0], element,xyz_list, gauss); this->GetLFSSSA(&LFSSSA[0][0],element, gauss); this->GetLprimeFSSSA(&LprimeFSSSA[0][0], element,xyz_list, gauss); element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); element->NormalBase(&bed_normal[0],xyz_list_tria); friction->GetAlpha2(&alpha2_gauss,gauss); DLSSAFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d; DLSSAFS[1][1]=alpha2_gauss*gauss->weight*Jdet2d; DLSSAFS[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; DLSSAFS[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; DLSSAFS[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]; DLSSAFS[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]; DLSSAFS[6][6]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[0]; DLSSAFS[7][7]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[1]; DLFSSSA[0][0]=alpha2_gauss*gauss->weight*Jdet2d; DLFSSSA[1][1]=alpha2_gauss*gauss->weight*Jdet2d; DLFSSSA[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; DLFSSSA[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; TripleMultiply( &LSSAFS[0][0],8,numdof2dm,1, &DLSSAFS[0][0],8,8,0, &LprimeSSAFS[0][0],8,numdofs,0, &Ke_drag[0][0],1); TripleMultiply( &LFSSSA[0][0],4,numdof2d,1, &DLFSSSA[0][0],4,4,0, &LprimeFSSSA[0][0],4,numdof2dm,0, &Ke_drag2[0][0],1); } for(i=0;ivalues[i*numdoftot+j+numdofm]+=Ke_drag[i][j]; for(i=0;ivalues[(i+numdofm)*numdoftot+j]+=Ke_drag2[i][j]; /*Transform Coordinate System*/ element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); /*Clean up and return*/ xDelete(cs_list); xDelete(node_list); xDelete(xyz_list); xDelete(xyz_list_tria); delete gauss; delete friction; return Ke; }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFSViscous(Element* element){/*{{{*/ /*Constants*/ const int numdofm = 2 *3; const int numdofs = 4 *6+ 3; const int numdoftotal = 2 *numdofm+numdofs; /*Intermediaries */ int i,j; int dim=3; IssmDouble Jdet,viscosity,FSreconditioning,D_scalar; IssmDouble B[4][numdofs]; IssmDouble Bprime[4][numdofm]; IssmDouble B2[3][numdofm]; IssmDouble Bprime2[3][numdofs]; IssmDouble D[4][4]={0.0}; // material matrix, simple scalar matrix. IssmDouble D2[3][3]={0.0}; // material matrix, simple scalar matrix. IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix IssmDouble *xyz_list = NULL; /*Find penta on bed as FS must be coupled to the dofs on the bed: */ Element* pentabase=element->GetBasalElement(); Element* basaltria=pentabase->SpawnBasalElement(); int vnumnodes = element->NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); int numnodes = 2*vnumnodes-1+pnumnodes; /*Prepare node list*/ int* cs_list = xNew(2*vnumnodes-1+pnumnodes); Node **node_list = xNew(2*vnumnodes-1+pnumnodes); for(i=0;iGetNode(i); cs_list[i] = XYEnum; } for(i=0;iGetNode(i); cs_list[i+vnumnodes-1] = XYZEnum; } for(i=0;iGetNode(vnumnodes+i); cs_list[2*vnumnodes-1+i] = PressureEnum; } /*Initialize Element matrix and return if necessary*/ ElementMatrix* Ke1=pentabase->NewElementMatrixCoupling(6,SSAApproximationEnum); ElementMatrix* Ke2=element->NewElementMatrix(FSvelocityEnum); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); delete Ke1; delete Ke2; /* Get node coordinates and dof list: */ 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); /* 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->GetBSSAFS(&B[0][0],element,xyz_list, gauss); this->GetBprimeSSAFSTria(&Bprime[0][0], basaltria,xyz_list, gauss_tria); this->GetBSSAFSTria(&B2[0][0], basaltria,xyz_list, gauss_tria); this->GetBprimeSSAFS(&Bprime2[0][0], element,xyz_list, gauss); element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); D_scalar=2*viscosity*gauss->weight*Jdet; for (i=0;i<3;i++) D[i][i]=D_scalar; D[3][3]=-gauss->weight*Jdet*FSreconditioning; for (i=0;i<3;i++) D2[i][i]=D_scalar; TripleMultiply( &B[0][0],4,numdofs,1, &D[0][0],4,4,0, &Bprime[0][0],4,numdofm,0, &Ke_gg[0][0],1); TripleMultiply( &B2[0][0],3,numdofm,1, &D2[0][0],3,3,0, &Bprime2[0][0],3,numdofs,0, &Ke_gg2[0][0],1); } for(i=0;ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j]; /*Transform Coordinate System*/ element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); /*Clean-up and return*/ xDelete(cs_list); xDelete(node_list); xDelete(xyz_list); delete basaltria->material; delete basaltria; delete gauss; delete gauss_tria; 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->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); element->material->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->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); } else if (approximation==SSAFSApproximationEnum){ element->material->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=CreatePVectorCouplingHOFSViscous(element); ElementVector* pe2=CreatePVectorCouplingHOFSFriction(element); ElementVector* pe =new ElementVector(pe1,pe2); /*clean-up and return*/ 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->material->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->material->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->material->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->material->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; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorSSAFS(Element* element){/*{{{*/ /*compute all load vectors for this element*/ int init = element->FiniteElement(); element->SetTemporaryElementType(P1Enum); // P1 needed for HO ElementVector* pe1=CreatePVectorSSA(element); element->SetTemporaryElementType(init); // P1 needed for HO ElementVector* pe2=CreatePVectorFS(element); int indices[3]={18,19,20}; element->SetTemporaryElementType(MINIcondensedEnum); // P1 needed for HO ElementMatrix* Ke = CreateKMatrixFS(element); element->SetTemporaryElementType(init); // P1 needed for HO pe2->StaticCondensation(Ke,3,&indices[0]); delete Ke; ElementVector* pe3=CreatePVectorCouplingSSAFS(element); ElementVector* pe =new ElementVector(pe1,pe2,pe3); /*clean-up and return*/ delete pe1; delete pe2; delete pe3; return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/ /*compute all load vectors for this element*/ ElementVector* pe1=CreatePVectorSSA(element); ElementVector* pe2=CreatePVectorHO(element); ElementVector* pe =new ElementVector(pe1,pe2); /*clean-up and return*/ delete pe1; delete pe2; return pe; }/*}}}*/ 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::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(); /*Get nodal functions*/ IssmDouble* dbasis=xNew(2*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build Bprime: */ for(int i=0;i(dbasis); }/*}}}*/ 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 * by: * Bi=[ dh/dx 0 0 0 ] * [ 0 dh/dy 0 0 ] * [ 1/2*dh/dy 1/2*dh/dx 0 0 ] * [ 0 0 0 h ] * where h is the interpolation function for node i. * * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) */ int i; IssmDouble dbasismini[3][7]; IssmDouble basis[6]; /*Get dbasis in actual coordinate system: */ element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss); element->NodalFunctionsP1(basis,gauss); /*Build B: */ for(i=0;i<6;i++){ B[(3*7+6)*0+3*i+0] = dbasismini[0][i]; B[(3*7+6)*0+3*i+1] = 0.; B[(3*7+6)*0+3*i+2] = 0.; B[(3*7+6)*1+3*i+0] = 0.; B[(3*7+6)*1+3*i+1] = dbasismini[1][i]; B[(3*7+6)*1+3*i+2] = 0.; B[(3*7+6)*2+3*i+0] = 0.5*dbasismini[1][i]; B[(3*7+6)*2+3*i+1] = 0.5*dbasismini[0][i]; B[(3*7+6)*2+3*i+2] = 0.; B[(3*7+6)*3+3*i+0] = 0.; B[(3*7+6)*3+3*i+1] = 0.; B[(3*7+6)*3+3*i+2] = 0.; } for(i=0;i<1;i++){ B[(3*7+6)*0+3*(6+i)+0] = 0.; B[(3*7+6)*0+3*(6+i)+1] = 0.; B[(3*7+6)*0+3*(6+i)+2] = 0.; B[(3*7+6)*1+3*(6+i)+0] = 0.; B[(3*7+6)*1+3*(6+i)+1] = 0.; B[(3*7+6)*1+3*(6+i)+2] = 0.; B[(3*7+6)*2+3*(6+i)+0] = 0.; B[(3*7+6)*2+3*(6+i)+1] = 0.; B[(3*7+6)*2+3*(6+i)+2] = 0.; B[(3*7+6)*3+3*(6+i)+0] = 0.; B[(3*7+6)*3+3*(6+i)+1] = 0.; B[(3*7+6)*3+3*(6+i)+2] = 0.; } for(i=0;i<6;i++){ //last column not for the bubble function B[(3*7+6)*0+7*3+i] = 0; B[(3*7+6)*1+7*3+i] = 0; B[(3*7+6)*2+7*3+i] = 0; B[(3*7+6)*3+7*3+i] = basis[i]; } }/*}}}*/ 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 * by: * Bi=[ dN/dx 0 ] * [ 0 dN/dy ] * [ 1/2*dN/dy 1/2*dN/dx ] * where N is the finiteelement function for node i. * * We assume B has been allocated already, of size: 3x(NDOF2*numnodes) */ /*Fetch number of nodes for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Get nodal functions derivatives*/ IssmDouble* dbasis=xNew(2*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B': */ for(int i=0;i(dbasis); }/*}}}*/ 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: * 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 B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) */ int numnodes = element->GetNumberOfNodes(); IssmDouble* dbasis=xNew(3*numnodes); /*Get dbasis in actual coordinate system: */ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B: */ for(int i=0;i(dbasis); }/*}}}*/ 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: * 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]; /*Cast gauss to GaussPenta*/ _assert_(gauss_in->Enum()==GaussPentaEnum); GaussPenta* gauss = xDynamicCast(gauss_in); /*Get basis in actual coordinate system: */ basis[0]=gauss->coord1*(1-gauss->coord4)/2.0; basis[1]=gauss->coord2*(1-gauss->coord4)/2.0; basis[2]=gauss->coord3*(1-gauss->coord4)/2.0; /*Build LprimeFS: */ 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::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 * by: * Lpi=[ h 0 0 0] * [ 0 h 0 0] * [ 0 0 h 0] * [ 0 0 h 0] * [ 0 0 dh/dz 0] * [ 0 0 dh/dz 0] * [ 0 0 0 h] * [ 0 0 0 h] * where h is the interpolation function for node i. */ int num_dof=3; int num_dof_vel=3*7; int num_dof_total=3*7+1*6; IssmDouble basis[3]; IssmDouble dbasis[3][6]; /*Cast gauss to GaussPenta*/ _assert_(gauss_in->Enum()==GaussPentaEnum); GaussPenta* gauss = xDynamicCast(gauss_in); /*Get basis in actual coordinate system: */ basis[0]=gauss->coord1*(1-gauss->coord4)/2.0; basis[1]=gauss->coord2*(1-gauss->coord4)/2.0; basis[2]=gauss->coord3*(1-gauss->coord4)/2.0; element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); /*Build LprimeFS: */ for(int i=0;i<3;i++){ LprimeFS[num_dof_total*0+num_dof*i+0] = basis[i]; LprimeFS[num_dof_total*0+num_dof*i+1] = 0.; LprimeFS[num_dof_total*0+num_dof*i+2] = 0.; LprimeFS[num_dof_total*1+num_dof*i+0] = 0.; LprimeFS[num_dof_total*1+num_dof*i+1] = basis[i]; LprimeFS[num_dof_total*1+num_dof*i+2] = 0.; LprimeFS[num_dof_total*2+num_dof*i+0] = 0.; LprimeFS[num_dof_total*2+num_dof*i+1] = 0.; LprimeFS[num_dof_total*2+num_dof*i+2] = basis[i]; LprimeFS[num_dof_total*3+num_dof*i+0] = 0.; LprimeFS[num_dof_total*3+num_dof*i+1] = 0.; LprimeFS[num_dof_total*3+num_dof*i+2] = basis[i]; LprimeFS[num_dof_total*4+num_dof*i+0] = 0.; LprimeFS[num_dof_total*4+num_dof*i+1] = 0.; LprimeFS[num_dof_total*4+num_dof*i+2] = dbasis[2][i]; LprimeFS[num_dof_total*5+num_dof*i+0] = 0.; LprimeFS[num_dof_total*5+num_dof*i+1] = 0.; LprimeFS[num_dof_total*5+num_dof*i+2] = dbasis[2][i]; LprimeFS[num_dof_total*6+num_dof*i+0] = 0.; LprimeFS[num_dof_total*6+num_dof*i+1] = 0.; LprimeFS[num_dof_total*6+num_dof*i+2] = 0.; LprimeFS[num_dof_total*7+num_dof*i+0] = 0.; LprimeFS[num_dof_total*7+num_dof*i+1] = 0.; LprimeFS[num_dof_total*7+num_dof*i+2] = 0.; } for(int i=3;i<7;i++){ LprimeFS[num_dof_total*0+num_dof*i+0] = 0.; LprimeFS[num_dof_total*0+num_dof*i+1] = 0.; LprimeFS[num_dof_total*0+num_dof*i+2] = 0.; LprimeFS[num_dof_total*1+num_dof*i+0] = 0.; LprimeFS[num_dof_total*1+num_dof*i+1] = 0.; LprimeFS[num_dof_total*1+num_dof*i+2] = 0.; LprimeFS[num_dof_total*2+num_dof*i+0] = 0.; LprimeFS[num_dof_total*2+num_dof*i+1] = 0.; LprimeFS[num_dof_total*2+num_dof*i+2] = 0.; LprimeFS[num_dof_total*3+num_dof*i+0] = 0.; LprimeFS[num_dof_total*3+num_dof*i+1] = 0.; LprimeFS[num_dof_total*3+num_dof*i+2] = 0.; LprimeFS[num_dof_total*4+num_dof*i+0] = 0.; LprimeFS[num_dof_total*4+num_dof*i+1] = 0.; LprimeFS[num_dof_total*4+num_dof*i+2] = 0.; LprimeFS[num_dof_total*5+num_dof*i+0] = 0.; LprimeFS[num_dof_total*5+num_dof*i+1] = 0.; LprimeFS[num_dof_total*5+num_dof*i+2] = 0.; LprimeFS[num_dof_total*6+num_dof*i+0] = 0.; LprimeFS[num_dof_total*6+num_dof*i+1] = 0.; LprimeFS[num_dof_total*6+num_dof*i+2] = 0.; LprimeFS[num_dof_total*7+num_dof*i+0] = 0.; LprimeFS[num_dof_total*7+num_dof*i+1] = 0.; LprimeFS[num_dof_total*7+num_dof*i+2] = 0.; } for(int i=0;i<3;i++){ LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*6+num_dof_vel+i] = basis[i]; LprimeFS[num_dof_total*7+num_dof_vel+i] = basis[i]; } for(int i=3;i<6;i++){ LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*6+num_dof_vel+i] = 0.; LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.; } }/*}}}*/ 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 * by: * Li=[ h 0 0 ] * [ 0 h 0 ] * [ 0 0 h ] * [ 0 0 h ] * where h is the interpolation function for node i. */ int num_dof=3; IssmDouble basis[3]; /*Cast gauss to GaussPenta*/ _assert_(gauss_in->Enum()==GaussPentaEnum); GaussPenta* gauss = xDynamicCast(gauss_in); /*Get basis in actual coordinate system: */ basis[0]=gauss->coord1*(1-gauss->coord4)/2.0; basis[1]=gauss->coord2*(1-gauss->coord4)/2.0; basis[2]=gauss->coord3*(1-gauss->coord4)/2.0; /*Build LFS: */ 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*0+num_dof*i+2] = 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*1+num_dof*i+2] = 0.; LFS[num_dof*3*2+num_dof*i+0] = 0.; LFS[num_dof*3*2+num_dof*i+1] = 0.; LFS[num_dof*3*2+num_dof*i+2] = basis[i]; LFS[num_dof*3*3+num_dof*i+0] = 0.; LFS[num_dof*3*3+num_dof*i+1] = 0.; LFS[num_dof*3*3+num_dof*i+2] = basis[i]; } }/*}}}*/ 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: * 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]; /*Cast gauss to GaussPenta*/ _assert_(gauss_in->Enum()==GaussPentaEnum); GaussPenta* gauss = xDynamicCast(gauss_in); /*Get basis in actual coordinate system: */ basis[0]=gauss->coord1*(1-gauss->coord4)/2.0; basis[1]=gauss->coord2*(1-gauss->coord4)/2.0; basis[2]=gauss->coord3*(1-gauss->coord4)/2.0; /*Build LFS: */ 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::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/ int i; IssmDouble rho_ice,g,FSreconditioning; int* doflistHO = NULL; int* doflistFSv = NULL; int* doflistFSp = NULL; /*Only works with Penta for now*/ if(element->ObjectEnum()!=PentaEnum) _error_("Coupling not supported for "<ObjectEnum())); /*Fetch number of nodes and dof for this finite element*/ int numnodes = 6; int numdofHO = 6*2; int numdofFSv = 6*3; int numdofFSp = 6; /*Fetch dof list and allocate solution vectors*/ element->GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum); element->GetDofList(&doflistHO, HOApproximationEnum, GsetEnum); element->GetDofListPressure(&doflistFSp,GsetEnum); IssmDouble* HOvalues = xNew(numdofHO); IssmDouble* FSvalues = xNew(numdofFSv+numdofFSp); IssmDouble* vx = xNew(numnodes); IssmDouble* vy = xNew(numnodes); IssmDouble* vz = xNew(numnodes); IssmDouble* vzHO = xNew(numnodes); IssmDouble* vzFS = xNew(numnodes); IssmDouble* vel = xNew(numnodes); IssmDouble* pressure = xNew(numnodes); /*Prepare coordinate system list*/ int* cs_list = xNew(2*numnodes); for(i=0;iFindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); for(i=0;iTransformSolutionCoord(FSvalues,2*numnodes,cs_list); element->TransformSolutionCoord(HOvalues,numnodes,XYEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[i])) _error_("Inf found in solution vector"); if(xIsNan(vy[i])) _error_("NaN found in solution vector"); if(xIsInf(vy[i])) _error_("Inf found in solution vector"); if(xIsNan(vzFS[i])) _error_("NaN found in solution vector"); if(xIsInf(vzFS[i])) _error_("Inf found in solution vector"); if(xIsNan(pressure[i])) _error_("NaN found in solution vector"); if(xIsInf(pressure[i])) _error_("Inf found in solution vector"); } /*Get Vz and compute vel*/ element->GetInputListOnVertices(vzHO,VzHOEnum); for(i=0;iInputChangeName(VxEnum,VxPicardEnum); element->InputChangeName(VyEnum,VyPicardEnum); element->InputChangeName(VzEnum,VzPicardEnum); element->InputChangeName(PressureEnum,PressurePicardEnum); /*Add vx and vy as inputs to element: */ element->AddInput(VxEnum,vx,P1Enum); element->AddInput(VyEnum,vy,P1Enum); element->AddInput(VzEnum,vz,P1Enum); element->AddInput(VzFSEnum,vzFS,P1Enum); element->AddInput(VelEnum,vel,P1Enum); element->AddInput(PressureEnum,pressure,P1Enum); /*Free ressources:*/ xDelete(pressure); xDelete(vel); xDelete(vz); xDelete(vzHO); xDelete(vzFS); xDelete(vy); xDelete(vx); xDelete(FSvalues); xDelete(HOvalues); xDelete(doflistFSp); xDelete(doflistFSv); xDelete(doflistHO); xDelete(cs_list); }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/ int i; IssmDouble rho_ice,g,FSreconditioning; int* doflistSSA = NULL; int* doflistFSv = NULL; int* doflistFSp = NULL; /*we have to add results of this element for FS and results from the element * at base for SSA, so we need to have the pointer toward the basal element*/ Element* basalelement=element->GetBasalElement(); if(basalelement->ObjectEnum()!=PentaEnum){ _error_("Coupling not supported for "<ObjectEnum())); } /*Fetch number of nodes and dof for this finite element*/ int numnodes = 6; int numdof2d = numnodes; int numdofSSA = 6*2; int numdofFSv = 6*3; int numdofFSp = 6; /*Fetch dof list and allocate solution vectors*/ element->GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum); element->GetDofListPressure(&doflistFSp,GsetEnum); basalelement->GetDofList(&doflistSSA, SSAApproximationEnum, GsetEnum); IssmDouble* SSAvalues = xNew(numdofSSA); IssmDouble* FSvalues = xNew(numdofFSv+numdofFSp); IssmDouble* vx = xNew(numnodes); IssmDouble* vy = xNew(numnodes); IssmDouble* vz = xNew(numnodes); IssmDouble* vzSSA = xNew(numnodes); IssmDouble* vzFS = xNew(numnodes); IssmDouble* vel = xNew(numnodes); IssmDouble* pressure = xNew(numnodes); /*Prepare coordinate system list*/ int* cs_list = xNew(2*numnodes); for(i=0;iFindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); for(i=0;iTransformSolutionCoord(FSvalues,2*numnodes,cs_list); element->TransformSolutionCoord(SSAvalues,numnodes,XYEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[i])) _error_("Inf found in solution vector"); if(xIsNan(vy[i])) _error_("NaN found in solution vector"); if(xIsInf(vy[i])) _error_("Inf found in solution vector"); if(xIsNan(vzFS[i])) _error_("NaN found in solution vector"); if(xIsInf(vzFS[i])) _error_("Inf found in solution vector"); if(xIsNan(pressure[i])) _error_("NaN found in solution vector"); if(xIsInf(pressure[i])) _error_("Inf found in solution vector"); } /*Get Vz and compute vel*/ element->GetInputListOnVertices(vzSSA,VzSSAEnum); for(i=0;iInputChangeName(VxEnum,VxPicardEnum); element->InputChangeName(VyEnum,VyPicardEnum); element->InputChangeName(VzEnum,VzPicardEnum); element->InputChangeName(PressureEnum,PressurePicardEnum); /*Add vx and vy as inputs to element: */ element->AddInput(VxEnum,vx,P1Enum); element->AddInput(VyEnum,vy,P1Enum); element->AddInput(VzEnum,vz,P1Enum); element->AddInput(VzFSEnum,vzFS,P1Enum); element->AddInput(VelEnum,vel,P1Enum); element->AddInput(PressureEnum,pressure,P1Enum); /*Free ressources:*/ xDelete(pressure); xDelete(vel); xDelete(vz); xDelete(vzSSA); xDelete(vzFS); xDelete(vy); xDelete(vx); xDelete(FSvalues); xDelete(SSAvalues); xDelete(doflistFSp); xDelete(doflistFSv); xDelete(doflistSSA); xDelete(cs_list); }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/ int i,domaintype; IssmDouble rho_ice,g; int* SSAdoflist = NULL; int* HOdoflist = NULL; IssmDouble* xyz_list = NULL; /*we have to add results of this element for HO and results from the element * at base for SSA, so we need to have the pointer toward the basal element*/ Element* basalelement=element->GetBasalElement(); if(basalelement->ObjectEnum()!=PentaEnum){ _error_("Coupling not supported for "<ObjectEnum())); } /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); int numdof = numnodes*2; int numdof2d = numnodes; /*Fetch dof list and allocate solution vectors*/ basalelement->GetDofList(&SSAdoflist,SSAApproximationEnum,GsetEnum); element ->GetDofList(&HOdoflist, HOApproximationEnum, GsetEnum); IssmDouble* HOvalues = xNew(numdof); IssmDouble* SSAvalues = xNew(numdof); IssmDouble* vx = xNew(numnodes); IssmDouble* vy = xNew(numnodes); IssmDouble* vz = xNew(numnodes); IssmDouble* vel = xNew(numnodes); IssmDouble* pressure = xNew(numnodes); IssmDouble* surface = xNew(numnodes); /*Use the dof list to index into the solution vector: */ for(i=0;iTransformSolutionCoord(SSAvalues,XYEnum); element->TransformSolutionCoord(HOvalues,XYEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(vx[i])) _error_("NaN found in solution vector"); if(xIsInf(vx[i])) _error_("Inf found in solution vector"); if(xIsNan(vy[i])) _error_("NaN found in solution vector"); if(xIsInf(vy[i])) _error_("Inf found in solution vector"); } /*Get Vz and compute vel*/ element->GetInputListOnNodes(&vz[0],VzEnum,0.); for(i=0;iGetMaterialParameter(MaterialsRhoIceEnum); g = element->GetMaterialParameter(ConstantsGEnum); element->GetVerticesCoordinates(&xyz_list); element->GetInputListOnNodes(&surface[0],SurfaceEnum); for(i=0;iInputChangeName(VxEnum,VxPicardEnum); element->InputChangeName(VyEnum,VyPicardEnum); element->InputChangeName(PressureEnum,PressurePicardEnum); /*Add vx and vy as inputs to element: */ element->AddInput(VxEnum,vx,P1Enum); element->AddInput(VyEnum,vy,P1Enum); element->AddInput(VelEnum,vel,P1Enum); element->AddInput(PressureEnum,pressure,P1Enum); /*Free ressources:*/ xDelete(surface); xDelete(pressure); xDelete(vel); xDelete(vz); xDelete(vy); xDelete(vx); xDelete(xyz_list); xDelete(SSAvalues); xDelete(HOvalues); xDelete(SSAdoflist); xDelete(HOdoflist); }/*}}}*/