Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18283) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18284) @@ -245,4 +245,5 @@ if(fe_FS==LATaylorHoodEnum){ iomodel->FetchDataToInput(elements,PressureEnum); + InputUpdateFromConstantx(elements,0.,SigmaNNEnum); } @@ -2896,11 +2897,6 @@ /*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; @@ -2920,14 +2916,11 @@ 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: */ @@ -2957,5 +2950,4 @@ BprimeU,1,numdof,0, BtBUzawa,1); - } @@ -2965,11 +2957,11 @@ 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); @@ -2978,5 +2970,5 @@ 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, @@ -2984,9 +2976,16 @@ 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); } @@ -2997,9 +2996,4 @@ MatrixMultiply(BtBUzawa,pnumdof,numdof,1, 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); @@ -3013,7 +3007,4 @@ xDelete(BU); xDelete(BtBUzawa); - xDelete(Cprime); - xDelete(C); - xDelete(CtCUzawa); xDelete(cs_list); return Ke; @@ -3297,4 +3288,19 @@ }/*}}}*/ 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; @@ -3369,4 +3375,114 @@ return pe3; } + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/ + + if(!element->IsOnBase()) return NULL; + + /*Intermediaries*/ + int dim; + IssmDouble alpha2,Jdet; + IssmDouble bed_normal[3]; + IssmDouble *xyz_list_base = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + + /*Initialize Element matrix and vectors*/ + ElementVector* pe = element->NewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + Input* alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input); + + /* Start looping on the number of gaussian points: */ + gauss=element->NewGaussBase(3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + alpha2_input->GetInputValue(&alpha2, gauss); + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + element->NormalBase(&bed_normal[0],xyz_list_base); + + for(int i=0;ivalues[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1]; + pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0]; + if(dim==3){ + pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i]; + } + } + + } + + /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list_base); + xDelete(vbasis); + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/ + + if(!element->IsOnBase()) return NULL; + + /*Intermediaries*/ + int dim; + IssmDouble sigmann,sigmant,Jdet,bedslope,beta; + IssmDouble *xyz_list_base = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + + /*Initialize Element matrix and vectors*/ + ElementVector* pe = element->NewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + Input* sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input); + Input* sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input); + Input* bedslope_input=element->GetInput(BedSlopeXEnum); _assert_(bedslope_input); + + /* Start looping on the number of gaussian points: */ + gauss=element->NewGaussBase(3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + sigmann_input->GetInputValue(&sigmann, gauss); + sigmant_input->GetInputValue(&sigmant, gauss); + bedslope_input->GetInputValue(&bedslope, gauss); + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + + beta=sqrt(1+bedslope*bedslope); + for(int i=0;ivalues[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i]; + pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i]; + if(dim==3){ + //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i]; + _error_("3d not supported yet"); + } + } + + } + + /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list_base); + xDelete(vbasis); return pe; }/*}}}*/ @@ -3502,5 +3618,4 @@ return pe; }/*}}}*/ -#endif ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/ @@ -3572,4 +3687,5 @@ return pe; }/*}}}*/ +#endif ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/ @@ -3772,5 +3888,5 @@ 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); @@ -3787,12 +3903,9 @@ pe->values[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); @@ -3802,5 +3915,5 @@ delete gauss; - Gauss* gauss=element->NewGaussBase(5); + gauss=element->NewGaussBase(5); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); @@ -3811,9 +3924,7 @@ 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]; } } @@ -3823,5 +3934,5 @@ /*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*/ @@ -4471,17 +4582,14 @@ /*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; Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18283) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18284) @@ -79,4 +79,6 @@ 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);