Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18283) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18284) @@ -78,6 +78,8 @@ ElementVector* CreatePVectorFSViscousXTH(Element* element); ElementVector* CreatePVectorFSShelf(Element* element); ElementVector* CreatePVectorFSFront(Element* element); + ElementVector* CreatePVectorFSFriction(Element* element); + ElementVector* CreatePVectorFSStress(Element* element); void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); void GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18283) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18284) @@ -244,6 +244,7 @@ iomodel->Constant(&fe_FS,FlowequationFeFSEnum); if(fe_FS==LATaylorHoodEnum){ iomodel->FetchDataToInput(elements,PressureEnum); + InputUpdateFromConstantx(elements,0.,SigmaNNEnum); } /*Friction law variables*/ @@ -2895,13 +2896,8 @@ /*Fetch number of nodes and dof for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); - int pnumnodes; - if(dim==2) pnumnodes=3; - else pnumnodes=6; - int lnumnodes; - if(dim==2) lnumnodes=2; - else lnumnodes=3; - //int pnumnodes = element->NumberofNodes(P1Enum); + int pnumnodes = element->GetNumberOfNodes(P1Enum); + int lnumnodes = element->GetNumberOfNodes(P2Enum); int numdof = vnumnodes*dim; int pnumdof = pnumnodes; int lnumdof = lnumnodes; @@ -2919,16 +2915,13 @@ IssmDouble* BU = xNew(pnumdof); IssmDouble* BprimeU = xNew(numdof); IssmDouble* D = xNewZeroInit(epssize*epssize); - IssmDouble* CtCUzawa = xNewZeroInit(numdof*pnumdof); - IssmDouble* C = xNew(lnumdof); - IssmDouble* Cprime = xNew(numdof); /*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; - if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} + 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); @@ -2956,7 +2949,6 @@ &DU,1,1,0, BprimeU,1,numdof,0, BtBUzawa,1); - } if(element->IsOnBase()){ @@ -2964,30 +2956,37 @@ element->GetVerticesCoordinatesBase(&xyz_list_base); element->NormalBase(&normal[0],xyz_list_base); - int lsize; - IssmDouble* Dlambda = xNewZeroInit(dim*dim); - IssmDouble* C = xNewZeroInit(dim*lnumdof); - IssmDouble* Cprime = xNewZeroInit(dim*numdof); + IssmDouble* Dlambda = xNewZeroInit(dim*dim); + IssmDouble* C = xNewZeroInit(dim*lnumdof); + IssmDouble* Cprime = xNewZeroInit(dim*numdof); + IssmDouble* CtCUzawa = xNewZeroInit(numdof*lnumdof); delete gauss; - Gauss* gauss=element->NewGaussBase(5); + 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); + 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(xyz_list_base); xDelete(Dlambda); xDelete(C); xDelete(Cprime); + xDelete(CtCUzawa); + xDelete(xyz_list_base); } /*Transform Coordinate System*/ @@ -2998,11 +2997,6 @@ BtBUzawa,pnumdof,numdof,0, &Ke->values[0],1); - /*The sigma naugmentation should not be transformed*/ - MatrixMultiply(CtCUzawa,lnumdof,numdof,1, - CtCUzawa,lnumdof,numdof,0, - &Ke->values[0],1); - /*Clean up and return*/ delete gauss; xDelete(xyz_list); @@ -3012,9 +3006,6 @@ xDelete(BprimeU); xDelete(BU); xDelete(BtBUzawa); - xDelete(Cprime); - xDelete(C); - xDelete(CtCUzawa); xDelete(cs_list); return Ke; }/*}}}*/ @@ -3297,6 +3288,21 @@ }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/ + 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::CreatePVectorFSViscous(Element* element){/*{{{*/ + int i,dim,fe_FS; IssmDouble x_coord,y_coord,z_coord; IssmDouble Jdet,forcex,forcey,forcez; @@ -3370,6 +3376,116 @@ } 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; +}/*}}}*/ #else ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/ @@ -3501,7 +3617,6 @@ /*clean-up and return*/ return pe; }/*}}}*/ -#endif ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/ int i,dim; @@ -3571,6 +3686,7 @@ xDelete(xyz_list); return pe; }/*}}}*/ +#endif ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/ int i,tausize,dim; @@ -3771,7 +3887,7 @@ element->FindParam(&r,AugmentedLagrangianREnum); element->GetVerticesCoordinates(&xyz_list); - /*Get d and tau*/ + /*Get pressure and sigmann*/ Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input); Input* sigmann_input =element->GetInput(SigmaNNEnum); _assert_(sigmann_input); @@ -3786,14 +3902,11 @@ 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(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i]; } } if(element->IsOnBase()){ - IssmDouble sigmann; IssmDouble* vbasis = xNew(numnodes); @@ -3801,7 +3914,7 @@ element->NormalBase(&bed_normal[0],xyz_list_base); delete gauss; - Gauss* gauss=element->NewGaussBase(5); + gauss=element->NewGaussBase(5); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); @@ -3810,11 +3923,9 @@ 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]; - } + pe->values[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); @@ -3822,7 +3933,7 @@ } /*Transform coordinate system*/ - //element->TransformLoadVectorCoord(pe,cs_list); Do not transform pressure augmentation + //element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation /*Clean up and return*/ delete gauss; @@ -4470,19 +4581,16 @@ */ /*Fetch number of nodes for this finite element*/ - int lnumnodes; - if(dim==2) lnumnodes=3; - else lnumnodes=6; - //int pnumnodes = element->NumberofNodes(P1Enum); + int lnumnodes = element->GetNumberOfNodes(P2Enum); /*Get nodal functions derivatives*/ IssmDouble* basis =xNew(lnumnodes); - element->NodalFunctions(basis,gauss); + element->NodalFunctionsP2(basis,gauss); /*Build B: */ for(int i=0;iGetNumberOfNodes(); + int vnumnodes = element->NumberofNodesVelocity(); int vnumdof = vnumnodes*dim; IssmDouble* vbasis=xNew(vnumnodes);