Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17296) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17297) @@ -6,5 +6,5 @@ #include "../solutionsequences/solutionsequences.h" -//#define FSANALYTICAL 5 +#define FSANALYTICAL 5 /*Model processing*/ @@ -1176,4 +1176,5 @@ /*Intermediaries*/ + int dim,meshtype; bool mainlyfloating; int migration_style,point1; @@ -1183,4 +1184,13 @@ Gauss* gauss = NULL; + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 1;break; + case Mesh2DhorizontalEnum: dim = 2;break; + case Mesh3DEnum: dim = 2;break; + default: _error_("mesh "<GetNumberOfNodes(); @@ -1188,7 +1198,7 @@ /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum); - IssmDouble* B = xNew(2*numdof); - IssmDouble D[2][2] = {0.}; + ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum); + IssmDouble* B = xNew(dim*numdof); + IssmDouble* D = xNewZeroInit(dim*dim); /*Retrieve all inputs and parameters*/ @@ -1197,6 +1207,7 @@ 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); Input* gllevelset_input = NULL; + Input* vy_input = NULL; + if(dim==2){vy_input = element->GetInput(VyEnum); _assert_(vy_input);} /*build friction object, used later on: */ @@ -1225,11 +1236,11 @@ } - this->GetBSSAFriction(B, element, xyz_list, gauss); + this->GetBSSAFriction(B,element,dim,xyz_list,gauss); element->JacobianDeterminant(&Jdet, xyz_list,gauss); - for(int i=0;i<2;i++) D[i][i]=alpha2*gauss->weight*Jdet; - - TripleMultiply(B,2,numdof,1, - &D[0][0],2,2,0, - B,2,numdof,0, + for(int i=0;iweight*Jdet; + + TripleMultiply(B,dim,numdof,1, + D,dim,dim,0, + B,dim,numdof,0, &Ke->values[0],1); } @@ -1243,4 +1254,5 @@ xDelete(xyz_list); xDelete(B); + xDelete(D); return Ke; }/*}}}*/ @@ -1251,4 +1263,5 @@ /*Intermediaries*/ + int dim,meshtype,bsize; IssmDouble viscosity,newviscosity,oldviscosity; IssmDouble viscosity_overshoot,thickness,Jdet; @@ -1256,13 +1269,22 @@ IssmDouble *xyz_list = NULL; + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 1; bsize = 1; break; + case Mesh2DhorizontalEnum: dim = 2; bsize = 1; break; + case Mesh3DEnum: dim = 2; bsize = 3; break; + default: _error_("mesh "<GetNumberOfNodes(); - int numdof = numnodes*2; + int numdof = numnodes*dim; /*Initialize Element matrix and vectors*/ ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum); - IssmDouble* B = xNew(3*numdof); - IssmDouble* Bprime = xNew(3*numdof); - IssmDouble* D = xNewZeroInit(3*3); + IssmDouble* B = xNew(bsize*numdof); + IssmDouble* Bprime = xNew(bsize*numdof); + IssmDouble* D = xNewZeroInit(bsize*bsize); /*Retrieve all inputs and parameters*/ @@ -1270,7 +1292,11 @@ Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 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* 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); @@ -1281,9 +1307,9 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); - this->GetBSSA(B,element,xyz_list,gauss); - this->GetBSSAprime(Bprime,element,xyz_list,gauss); - - element->ViscositySSA(&viscosity,xyz_list,gauss,vx_input,vy_input); - element->ViscositySSA(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input); + this->GetBSSA(B,element,dim,xyz_list,gauss); + this->GetBSSAprime(Bprime,element,dim,xyz_list,gauss); + + element->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); + element->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input); thickness_input->GetInputValue(&thickness, gauss); @@ -1292,12 +1318,12 @@ for(int i=0;i<3;i++) D[i*3+i]=D_scalar; - TripleMultiply(B,3,numdof,1, - D,3,3,0, - Bprime,3,numdof,0, + TripleMultiply(B,bsize,numdof,1, + D,bsize,bsize,0, + Bprime,bsize,numdof,0, &Ke->values[0],1); } /*Transform Coordinate System*/ - element->TransformStiffnessMatrixCoord(Ke,XYEnum); + if(dim==2) element->TransformStiffnessMatrixCoord(Ke,XYEnum); /*Clean up and return*/ @@ -1348,6 +1374,16 @@ /*Intermediaries */ + int dim,meshtype; IssmDouble thickness,Jdet,slope[2]; IssmDouble* xyz_list = NULL; + + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 1;break; + case Mesh2DhorizontalEnum: dim = 2;break; + case Mesh3DEnum: dim = 2;break; + default: _error_("mesh "<values[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]; + pe->values[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]; } } @@ -1400,4 +1436,5 @@ /*Intermediaries*/ + int dim,meshtype; IssmDouble Jdet,thickness,bed,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; @@ -1405,4 +1442,13 @@ IssmDouble *xyz_list_front = NULL; IssmDouble normal[2]; + + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 1;break; + case Mesh2DhorizontalEnum: dim = 2;break; + case Mesh3DEnum: dim = 2;break; + default: _error_("mesh "<GetMaterialParameter(ConstantsGEnum); element->GetVerticesCoordinates(&xyz_list); -// element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + //element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); element->NormalSection(&normal[0],xyz_list_front); @@ -1440,7 +1486,7 @@ 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]; + 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]; } } @@ -1456,11 +1502,12 @@ return pe; }/*}}}*/ -void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: - * Bi=[ dN/dx 0 ] - * [ 0 dN/dy ] - * [ 1/2*dN/dy 1/2*dN/dx ] + * 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. * @@ -1472,15 +1519,22 @@ /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(2*numnodes); + IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B: */ - for(int i=0;i(dbasis); }/*}}}*/ -void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. * For node i, Bi can be expressed in the actual coordinate system * by: - * Bi=[ N 0 ] + * 2D 1D + * Bi=[ N 0 ] Bi = N * [ 0 N ] * where N is the finiteelement function for node i. @@ -1507,9 +1562,16 @@ /*Build L: */ - for(int i=0;i(basis); }/*}}}*/ -void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system * by: - * Bi_prime=[ 2*dN/dx dN/dy ] + * 2D 1D + * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] * [ dN/dx 2*dN/dy ] * [ dN/dy dN/dx ] @@ -1533,15 +1596,22 @@ /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(2*numnodes); + IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B': */ - for(int i=0;iGetInputListOnVertices(thickness,ThicknessEnum); for(i=0;iGetInputListOnVertices(surface,SurfaceEnum); for(i=0;iGetVerticesCoordinates(&xyz_list); + element->GetInputListOnVertices(surface,SurfaceEnum); + for(i=0;iIsOnBed()){xDelete(xyz_list); return;} basalelement=element->SpawnBasalElement(); @@ -1597,5 +1675,5 @@ /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); - int numdof = numnodes*2; + int numdof = numnodes*dim; /*Fetch dof list and allocate solution vectors*/ @@ -1611,29 +1689,35 @@ /*Transform solution in Cartesian Space*/ - basalelement->TransformSolutionCoord(&values[0],XYEnum); - basalelement->FindParam(&meshtype,MeshTypeEnum); + if(dim==2) basalelement->TransformSolutionCoord(&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(xIsNan(vy[i])) _error_("NaN found in solution vector"); + + if(dim==2){ + vy[i]=values[i*2+1]; + if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + } } /*Get Vz and compute vel*/ - basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); - for(i=0;iGetInputListOnNodes(&vz[0],VzEnum,0.); + for(i=0;iGetInputListOnNodes(&vy[0],VyEnum,0.); + for(i=0;iInputChangeName(VxEnum,VxPicardEnum); - element->InputChangeName(VyEnum,VyPicardEnum); + if(dim==2)element->InputChangeName(VyEnum,VyPicardEnum); /*Add vx and vy as inputs to the tria element: */ element->AddBasalInput(VxEnum,vx,P1Enum); - element->AddBasalInput(VyEnum,vy,P1Enum); + if(dim==2)element->AddBasalInput(VyEnum,vy,P1Enum); element->AddBasalInput(VelEnum,vel,P1Enum); @@ -1706,6 +1790,6 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); - this->GetBSSA(B,basalelement,xyz_list,gauss_base); - this->GetBSSAprime(Bprime,basalelement,xyz_list,gauss_base); + this->GetBSSA(B,basalelement,2,xyz_list,gauss_base); + this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base); element->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); @@ -2428,5 +2512,5 @@ Bprime[numnodes*1+i] = dbasis[1*numnodes+i]; } - } + } /*Clean-up*/ @@ -3696,6 +3780,6 @@ element->JacobianDeterminant(&Jdet, xyz_list,gauss); - this->GetBSSA(&B[0][0],basaltria, xyz_list, gauss_tria); - this->GetBSSAprime(&Bprime[0][0], basaltria,xyz_list, gauss_tria); + 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){ @@ -3998,5 +4082,5 @@ element->JacobianDeterminant(&Jdet, xyz_list,gauss); this->GetBSSAHO(B, element,xyz_list, gauss); - this->GetBSSAprime(Bprime, basaltria,xyz_list, gauss_tria); + this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17296) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17297) @@ -39,7 +39,7 @@ ElementVector* CreatePVectorSSADrivingStress(Element* element); ElementVector* CreatePVectorSSAFront(Element* element); - void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element); /*L1L2*/