Index: ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22873) +++ ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22874) @@ -11,7 +11,7 @@ int M,N; int finiteelement; iomodel->FindConstant(&finiteelement,"md.thermal.fe"); - _assert_(finiteelement==P1Enum); + _assert_(finiteelement==P1Enum); /*Output*/ IssmDouble *spcvector = NULL; @@ -93,8 +93,8 @@ int finiteelement; iomodel->FindConstant(&finiteelement,"md.thermal.fe"); - _assert_(finiteelement==P1Enum); - + _assert_(finiteelement==P1Enum); + if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); ::CreateNodes(nodes,iomodel,ThermalAnalysisEnum,finiteelement); iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); @@ -112,7 +112,7 @@ /*Update elements: */ int finiteelement; iomodel->FindConstant(&finiteelement,"md.thermal.fe"); - _assert_(finiteelement==P1Enum); + _assert_(finiteelement==P1Enum); int counter=0; for(int i=0;inumberofelements;i++){ if(iomodel->my_elements[i]){ @@ -188,8 +188,10 @@ 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); + if (FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; case 2: @@ -201,8 +203,10 @@ 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); + if (FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; case 4: @@ -230,8 +234,10 @@ iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum); 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); + if (FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; case 9: @@ -452,7 +458,7 @@ K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz); K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz); for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j]; - GetBAdvecprime(Bprime,element,xyz_list,gauss); + GetBAdvecprime(Bprime,element,xyz_list,gauss); TripleMultiply(Bprime,3,numnodes,1, &K[0][0],3,3,0, @@ -487,7 +493,7 @@ return Ke; }/*}}}*/ ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/ - + /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; @@ -700,9 +706,9 @@ }/*}}}*/ void ThermalAnalysis::GetBAdvec(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*NDOF1. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Bi_advec =[ h ] * [ h ] * [ h ] @@ -729,9 +735,9 @@ xDelete(basis); }/*}}}*/ void ThermalAnalysis::GetBAdvecprime(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*NDOF1. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Biprime_advec=[ dh/dx ] * [ dh/dy ] * [ dh/dz ] @@ -758,9 +764,9 @@ xDelete(dbasis); }/*}}}*/ void ThermalAnalysis::GetBConduct(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*NDOF1. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Bi_conduct=[ dh/dx ] * [ dh/dy ] * [ dh/dz ] Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 22873) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 22874) @@ -610,7 +610,7 @@ } break; case L1L2ApproximationEnum: numdofs =2; break; - case HOApproximationEnum: + case HOApproximationEnum: switch(domaintype){ case Domain3DEnum: numdofs=2; break; case Domain2DverticalEnum: numdofs=1; break; @@ -784,8 +784,10 @@ 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); + if(FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; case 2: @@ -797,8 +799,10 @@ 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); + if(FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; case 4: @@ -826,8 +830,11 @@ iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum); 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); + if(FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); + } break; case 9: @@ -933,9 +940,9 @@ else if(newton>0) solutionsequence_newton(femmodel); else - solutionsequence_nonlinear(femmodel,conserve_loads); + solutionsequence_nonlinear(femmodel,conserve_loads); } - else if(!isFS && (isSSA || isHO || isL1L2)){ + else if(!isFS && (isSSA || isHO || isL1L2)){ if(VerboseSolution()) _printf0_(" computing velocities\n"); femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); @@ -942,7 +949,7 @@ if(newton>0) solutionsequence_newton(femmodel); else - solutionsequence_nonlinear(femmodel,conserve_loads); + solutionsequence_nonlinear(femmodel,conserve_loads); if(domaintype==Domain2DverticalEnum && isSSA){ femmodel->parameters->SetParam(VxEnum,InputToExtrudeEnum); @@ -965,7 +972,7 @@ int approximation; element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ - case FSApproximationEnum: + case FSApproximationEnum: return CreateDVectorFS(element); default: return NULL; //no need for doftypes outside of FS approximation @@ -978,11 +985,11 @@ int approximation; element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ - case SSAApproximationEnum: + case SSAApproximationEnum: return CreateJacobianMatrixSSA(element); - case HOApproximationEnum: + case HOApproximationEnum: return CreateJacobianMatrixHO(element); - case FSApproximationEnum: + case FSApproximationEnum: return CreateJacobianMatrixFS(element); case NoneApproximationEnum: return NULL; @@ -996,19 +1003,19 @@ switch(approximation){ case SIAApproximationEnum: return NULL; - case SSAApproximationEnum: + case SSAApproximationEnum: return CreateKMatrixSSA(element); - case L1L2ApproximationEnum: + case L1L2ApproximationEnum: return CreateKMatrixL1L2(element); - case HOApproximationEnum: + case HOApproximationEnum: return CreateKMatrixHO(element); - case FSApproximationEnum: + case FSApproximationEnum: return CreateKMatrixFS(element); - case SSAHOApproximationEnum: + case SSAHOApproximationEnum: return CreateKMatrixSSAHO(element); - case HOFSApproximationEnum: + case HOFSApproximationEnum: return CreateKMatrixHOFS(element); - case SSAFSApproximationEnum: + case SSAFSApproximationEnum: return CreateKMatrixSSAFS(element); case NoneApproximationEnum: return NULL; @@ -1023,19 +1030,19 @@ switch(approximation){ case SIAApproximationEnum: return NULL; - case SSAApproximationEnum: + case SSAApproximationEnum: return CreatePVectorSSA(element); - case L1L2ApproximationEnum: + case L1L2ApproximationEnum: return CreatePVectorL1L2(element); - case HOApproximationEnum: + case HOApproximationEnum: return CreatePVectorHO(element); - case FSApproximationEnum: + case FSApproximationEnum: return CreatePVectorFS(element); - case SSAHOApproximationEnum: + case SSAHOApproximationEnum: return CreatePVectorSSAHO(element); - case HOFSApproximationEnum: + case HOFSApproximationEnum: return CreatePVectorHOFS(element); - case SSAFSApproximationEnum: + case SSAFSApproximationEnum: return CreatePVectorSSAFS(element); case NoneApproximationEnum: return NULL; @@ -1125,12 +1132,12 @@ case FSApproximationEnum: case NoneApproximationEnum: InputUpdateFromSolutionFS(solution,element); return; - case SIAApproximationEnum: + case SIAApproximationEnum: return; - case SSAApproximationEnum: + case SSAApproximationEnum: InputUpdateFromSolutionSSA(solution,element); return; - case HOApproximationEnum: + case HOApproximationEnum: InputUpdateFromSolutionHO(solution,element); return; case L1L2ApproximationEnum: @@ -1580,7 +1587,7 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); + 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); @@ -1688,13 +1695,13 @@ 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. + /*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: + * by: * 2D 1D * Bi=[ dN/dx 0 ] Bi=[ dN/dx ] - * [ 0 dN/dy ] - * [ 1/2*dN/dy 1/2*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) @@ -1728,9 +1735,9 @@ xDelete(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. + /*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: + * by: * 2D 1D * Bi=[ N 0 ] Bi = N * [ 0 N ] @@ -1765,9 +1772,9 @@ xDelete(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. + /*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: + * by: * 2D 1D * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] * [ dN/dx 2*dN/dy ] @@ -2049,7 +2056,7 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); + 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); @@ -2676,16 +2683,16 @@ 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. + /*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: + * 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 ] + * [ 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) @@ -2725,9 +2732,9 @@ xDelete(dbasis); }/*}}}*/ void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. + /*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: + * by: * 3D 2D * Bi=[ N 0 ] Bi=N * [ 0 N ] @@ -2762,13 +2769,13 @@ xDelete(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. + /*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: + * 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 ] + * [ dN/dy dN/dx ] * where hNis the finiteelement function for node i. * * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes) @@ -3083,7 +3090,7 @@ 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; + 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]; @@ -3285,7 +3292,7 @@ BtBUzawa,pnumdof,numdof,0, &Ke->values[0],1); - if(element->IsOnBase() && 0){ + if(element->IsOnBase() && 0){ element->FindParam(&rl,AugmentedLagrangianRlambdaEnum); element->GetVerticesCoordinatesBase(&xyz_list_base); element->NormalBase(&normal[0],xyz_list_base); @@ -4084,7 +4091,7 @@ 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); @@ -4095,7 +4102,7 @@ } } - if(element->IsOnBase() && 0){ + if(element->IsOnBase() && 0){ IssmDouble sigmann; IssmDouble* vbasis = xNew(numnodes); @@ -4304,7 +4311,7 @@ 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. + /*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 ] @@ -4416,9 +4423,9 @@ 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. + /* 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 + * by in 3d * Li=[ h 0 0 0 ] * [ 0 h 0 0 ] * in 2d: @@ -4476,9 +4483,9 @@ xDelete(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. + /* 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: + * by: * Bvi' = [ dphi/dx 0 ] * [ 0 dphi/dy ] * [ dphi/dy dphi/dx ] @@ -4589,9 +4596,9 @@ 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. + /* 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: + * by: * Bvi' = [ dphi/dx dphi/dy ] * * In 3d @@ -4625,9 +4632,9 @@ xDelete(vdbasis); }/*}}}*/ void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. + /* 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: + * by: * Bvi' = [ dphi/dx 0 ] * [ 0 dphi/dy ] * [ dphi/dy dphi/dx ] @@ -4688,7 +4695,7 @@ xDelete(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. + /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. */ /*Fetch number of nodes for this finite element*/ @@ -4710,7 +4717,7 @@ xDelete(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. + /*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 ] @@ -4775,7 +4782,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute C matrix. C=[Cp1 Cp2 ...] where: - * Cpi=[phi phi]. + * Cpi=[phi phi]. */ /*Fetch number of nodes for this finite element*/ @@ -4795,7 +4802,7 @@ xDelete(basis); }/*}}}*/ void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute C' matrix. C'=[C1' C2' ...] + /* Compute C' matrix. C'=[C1' C2' ...] * Ci' = [ phi 0 ] * [ 0 phi ] * @@ -4929,15 +4936,15 @@ /*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; + 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); + epsxz = xNew(tnumnodes); sigmapxz = xNew(tnumnodes); + epsyz = xNew(tnumnodes); sigmapyz = xNew(tnumnodes); } /*Get d and tau*/ @@ -5171,9 +5178,9 @@ } epsxx = dvx[0]; epsyy = dvy[1]; - epsxy = 0.5*(dvx[1] + dvy[0]); + epsxy = 0.5*(dvx[1] + dvy[0]); if(dim==3){ - epszz = dvz[2]; + epszz = dvz[2]; epsxz = 0.5*(dvx[2] + dvz[0]); epsyz = 0.5*(dvy[2] + dvz[1]); } @@ -5182,7 +5189,7 @@ 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 ) + 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( @@ -5206,7 +5213,7 @@ 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 + 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*/ @@ -5359,9 +5366,9 @@ } epsxx = dvx[0]; epsyy = dvy[1]; - epsxy = 0.5*(dvx[1] + dvy[0]); + epsxy = 0.5*(dvx[1] + dvy[0]); if(dim==3){ - epszz = dvz[2]; + epszz = dvz[2]; epsxz = 0.5*(dvx[2] + dvz[0]); epsyz = 0.5*(dvy[2] + dvz[1]); } @@ -5613,7 +5620,7 @@ 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); @@ -5646,8 +5653,8 @@ 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 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: */ @@ -5718,7 +5725,7 @@ &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]; @@ -5802,7 +5809,7 @@ 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; + for (i=0;i<2;i++) DL[i][i]=DL_scalar; /* Do the triple producte tL*D*L: */ TripleMultiply( L,2,numdof,1, @@ -5883,7 +5890,7 @@ element->JacobianDeterminant(&Jdet, xyz_list,gauss); this->GetBSSAHO(B, element,xyz_list, gauss); - this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); + this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); D_scalar=2*viscosity*gauss->weight*Jdet; @@ -5893,7 +5900,7 @@ &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]; @@ -5902,7 +5909,7 @@ /*Clean-up and return*/ basaltria->DeleteMaterials(); delete basaltria; - + delete gauss; delete gauss_tria; xDelete(B); @@ -5940,9 +5947,9 @@ int indices[3]={18,19,20}; Ke1->StaticCondensation(3,&indices[0]); int init = element->FiniteElement(); - element->SetTemporaryElementType(P1Enum); + element->SetTemporaryElementType(P1Enum); ElementMatrix* Ke2=CreateKMatrixSSA3d(element); - element->SetTemporaryElementType(init); + element->SetTemporaryElementType(init); ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); @@ -5986,7 +5993,7 @@ /*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 + /*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(); @@ -6222,7 +6229,7 @@ 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); @@ -6468,9 +6475,9 @@ 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. + /*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: + * by: * Bprimei=[ 2*dh/dx dh/dy 0 0 ] * [ dh/dx 2*dh/dy 0 0 ] * [ dh/dy dh/dx 0 0 ] @@ -6517,9 +6524,9 @@ } }/*}}}*/ 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. + /*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: + * by: * Bprimei=[ dN/dx 0 ] * [ 0 dN/dy ] * [ dN/dy dN/dx ] @@ -6552,9 +6559,9 @@ xDelete(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. + /*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: + * by: * Bi=[ dh/dx 0 0 0 ] * [ 0 dh/dy 0 0 ] * [ 1/2*dh/dy 1/2*dh/dx 0 0 ] @@ -6610,9 +6617,9 @@ } }/*}}}*/ 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. + /*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: + * by: * Bi=[ dN/dx 0 ] * [ 0 dN/dy ] * [ 1/2*dN/dy 1/2*dN/dx ] @@ -6642,9 +6649,9 @@ xDelete(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. + /*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: + * by: * Bi=[ dh/dx 0 ] * [ 0 dh/dy ] * [ 1/2*dh/dy 1/2*dh/dx ] @@ -6673,9 +6680,9 @@ xDelete(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. + /* 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: + * by: * Lpi=[ h 0 ] * [ 0 h ] * [ h 0 ] @@ -6707,9 +6714,9 @@ } }/*}}}*/ 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. + /* 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: + * by: * Lpi=[ h 0 0 0] * [ 0 h 0 0] * [ 0 0 h 0] @@ -6812,9 +6819,9 @@ } }/*}}}*/ 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. + /* 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: + * by: * Li=[ h 0 0 ] * [ 0 h 0 ] * [ 0 0 h ] @@ -6852,9 +6859,9 @@ }/*}}}*/ 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. + * 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: + * by: * Li=[ h 0 ] * [ 0 h ] * [ h 0 ] @@ -7157,7 +7164,7 @@ element->GetInputListOnNodes(&vz[0],VzEnum,0.); for(i=0;iGetMaterialParameter(MaterialsRhoIceEnum); g = element->GetMaterialParameter(ConstantsGEnum);