Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11332) @@ -6257,71 +6257,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/ -ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){ - - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; - - /*Intermediaries */ - int i,j,ig; - int approximation; - double xyz_list[NUMVERTICES][3]; - double Jdet; - double eps1dotdphii,eps1dotdphij; - double eps2dotdphii,eps2dotdphij; - double mu_prime; - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double eps1[3],eps2[3]; - double phi[NUMVERTICES]; - double dphi[3][NUMVERTICES]; - GaussPenta *gauss=NULL; - - /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/ - ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn(); - - /*Retrieve all inputs and parameters*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - - this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityDerivativeEpsSquare(&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(i=0;i<6;i++){ - for(j=0;j<6;j++){ - eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; - eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; - eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; - eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; - - Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; - Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; - Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; - Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; - } - } - } - - /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{1*/ ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){ @@ -7491,4 +7424,71 @@ } /*}}}*/ +/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/ +ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + int approximation; + double xyz_list[NUMVERTICES][3]; + double Jdet; + double eps1dotdphii,eps1dotdphij; + double eps2dotdphii,eps2dotdphij; + double mu_prime; + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double eps1[3],eps2[3]; + double phi[NUMVERTICES]; + double dphi[3][NUMVERTICES]; + GaussPenta *gauss=NULL; + + /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/ + ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn(); + + /*Retrieve all inputs and parameters*/ + inputs->GetInputValue(&approximation,ApproximationEnum); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); + + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosityDerivativeEpsSquare(&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(i=0;i<6;i++){ + for(j=0;j<6;j++){ + eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; + eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; + eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; + eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; + + Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/ void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11332) @@ -886,4 +886,37 @@ delete gauss; return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreateJacobianMatrix{{{1*/ +void Tria::CreateJacobianMatrix(Mat Jff){ + + /*retrieve parameters: */ + ElementMatrix* Ke=NULL; + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Checks in debugging {{{2*/ + _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); + /*}}}*/ + + /*Skip if water element*/ + if(IsOnWater()) return; + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ +#ifdef _HAVE_DIAGNOSTIC_ + case DiagnosticHorizAnalysisEnum: + Ke=CreateJacobianDiagnosticMacayeal(); + break; +#endif + default: + _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type)); + } + + /*Add to global matrix*/ + if(Ke){ + Ke->AddToGlobal(Jff); + delete Ke; + } } /*}}}*/ @@ -3028,4 +3061,70 @@ delete gauss; return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/ +ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double xyz_list[NUMVERTICES][3]; + double Jdet,thickness; + double eps1dotdphii,eps1dotdphij; + double eps2dotdphii,eps2dotdphij; + double mu_prime; + double epsilon[3];/* epsilon=[exx,eyy,exy];*/ + double eps1[2],eps2[2]; + double phi[NUMVERTICES]; + double dphi[2][NUMVERTICES]; + GaussTria *gauss=NULL; + + /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/ + ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss); + + thickness_input->GetInputValue(&thickness, gauss); + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosity2dDerivativeEpsSquare(&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(i=0;i<3;i++){ + for(j=0;j<3;j++){ + eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]; + eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]; + eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]; + eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]; + + Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii; + Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii; + Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii; + Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); + + /*Clean up and return*/ + delete gauss; + return Ke; } /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11332) @@ -84,5 +84,5 @@ void CreateKMatrix(Mat Kff, Mat Kfs,Vec df); void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");}; + void CreateJacobianMatrix(Mat Jff); int GetNodeIndex(Node* node); int Sid(); @@ -208,4 +208,5 @@ ElementVector* CreatePVectorDiagnosticMacAyeal(void); ElementVector* CreatePVectorDiagnosticHutter(void); + ElementMatrix* CreateJacobianDiagnosticMacayeal(void); void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsDiagnosticHutter(Vec solution); Index: /issm/trunk-jpl/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Icefront.cpp (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Loads/Icefront.cpp (revision 11332) @@ -361,4 +361,9 @@ } /*}}}*/ +/*FUNCTION Icefront::CreateJacobianMatrix{{{1*/ +void Icefront::CreateJacobianMatrix(Mat Jff){ + this->CreateKMatrix(Jff,NULL); +} +/*}}}1*/ /*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/ void Icefront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax){ @@ -373,4 +378,9 @@ } /*}}}*/ +/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{1*/ +void Icefront::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){ + this->PenaltyCreateKMatrix(Jff,NULL,kmax); +} +/*}}}1*/ /*FUNCTION Icefront::InAnalysis{{{1*/ bool Icefront::InAnalysis(int in_analysis_type){ Index: /issm/trunk-jpl/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Icefront.h (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Loads/Icefront.h (revision 11332) @@ -76,8 +76,8 @@ void CreateKMatrix(Mat Kff, Mat Kfs); void CreatePVector(Vec pf); - void CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");}; - void PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");}; + void CreateJacobianMatrix(Mat Jff); void PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax); void PenaltyCreatePVector(Vec pf, double kmax); + void PenaltyCreateJacobianMatrix(Mat Jff,double kmax); bool InAnalysis(int analysis_type); /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp (revision 11332) @@ -216,8 +216,5 @@ /*FUNCTION Penpair::CreateJacobianMatrix{{{1*/ void Penpair::CreateJacobianMatrix(Mat Jff){ - /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/ - /*No loads applied, do nothing: */ - return; - + this->CreateKMatrix(Jff,NULL); } /*}}}1*/ Index: /issm/trunk-jpl/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Materials/Matice.cpp (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Materials/Matice.cpp (revision 11332) @@ -549,4 +549,35 @@ eyz=epsilon[4]; eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; + + mu_prime=(1-n)/(2*n) * mu/eff2; + } + + /*Assign output pointers:*/ + *pmu_prime=mu_prime; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosity2dDerivativeEpsSquare{{{1*/ +void Matice::GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* epsilon){ + + /*output: */ + double mu_prime; + double mu,n,eff2; + + /*input strain rate: */ + double exx,eyy,exy,exz; + + /*Get visocisty and n*/ + GetViscosity2d(&mu,epsilon); + n=GetN(); + + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ + mu_prime=0.5*pow((double)10,(double)14); + } + else{ + /*Retrive strain rate components: */ + exx=epsilon[0]; + eyy=epsilon[1]; + exy=epsilon[2]; + eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ; mu_prime=(1-n)/(2*n) * mu/eff2; Index: /issm/trunk-jpl/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Materials/Matice.h (revision 11331) +++ /issm/trunk-jpl/src/c/objects/Materials/Matice.h (revision 11332) @@ -69,4 +69,5 @@ void GetViscosityComplement(double* pviscosity_complement, double* pepsilon); void GetViscosityDerivativeEpsSquare(double* pmu_prime, double* pepsilon); + void GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* pepsilon); double GetB(); double GetBbar(); Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp (revision 11331) +++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp (revision 11332) @@ -61,6 +61,5 @@ convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); - count++; - if(converged==true)break; + if(converged==true) break; if(count>=max_nonlinear_iterations){ _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); @@ -71,5 +70,5 @@ SystemMatricesx(&Kff,NULL,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax); - CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);VecScale(ys,0.); VecDuplicate(pf,&pJf); MatMultPatch(Kff,uf,pJf); MatFree(&Kff); @@ -81,4 +80,6 @@ Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys); InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + + count++; }