Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11447) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11448) @@ -585,6 +585,9 @@ switch(analysis_type){ #ifdef _HAVE_DIAGNOSTIC_ - case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: + case DiagnosticHorizAnalysisEnum: Ke=CreateKMatrixDiagnosticHoriz(); De=CreateDVectorDiagnosticHoriz(); + break; + case AdjointHorizAnalysisEnum: + Ke=CreateKMatrixAdjointHoriz(); break; case DiagnosticHutterAnalysisEnum: @@ -4443,4 +4446,194 @@ }/*}}}*/ +/*FUNCTION Penta::CreateKMatrixAdjointHoriz{{{1*/ +ElementMatrix* Penta::CreateKMatrixAdjointHoriz(void){ + + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + + switch(approximation){ + case MacAyealApproximationEnum: + return CreateKMatrixAdjointMacAyeal2d(); + case PattynApproximationEnum: + return CreateKMatrixAdjointPattyn(); + case StokesApproximationEnum: + return CreateKMatrixAdjointStokes(); + case NoneApproximationEnum: + return NULL; + default: + _error_("Approximation %s not supported yet",EnumToStringx(approximation)); + } +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixAdjointMacAyeal2d{{{1*/ +ElementMatrix* Penta::CreateKMatrixAdjointMacAyeal2d(void){ + + /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build + the stiffness matrix. */ + if (!IsOnBed()) return NULL; + + /*Depth Averaging B*/ + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementMatrix* Ke=tria->CreateKMatrixAdjointMacAyeal(); + delete tria->matice; delete tria; + + /*Delete B averaged*/ + this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + + /*clean up and return*/ + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixAdjointPattyn{{{1*/ +ElementMatrix* Penta::CreateKMatrixAdjointPattyn(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + bool incomplete_adjoint; + 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)*/ + parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn(); + if(incomplete_adjoint) return Ke; + + /*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); + + /* 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::CreateKMatrixAdjointStokes{{{1*/ +ElementMatrix* Penta::CreateKMatrixAdjointStokes(void){ + + /*Constants*/ + const int numdof=NDOF4*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + bool incomplete_adjoint; + double xyz_list[NUMVERTICES][3]; + double Jdet; + double eps1dotdphii,eps1dotdphij; + double eps2dotdphii,eps2dotdphij; + double eps3dotdphii,eps3dotdphij; + double mu_prime; + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double eps1[3],eps2[3],eps3[3]; + double phi[NUMVERTICES]; + double dphi[3][NUMVERTICES]; + GaussPenta *gauss=NULL; + + /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/ + parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + ElementMatrix* Ke=CreateKMatrixDiagnosticStokes(); + if(incomplete_adjoint) return Ke; + + /*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* vz_input=inputs->GetInput(VzEnum); _assert_(vz_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]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; + eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; + + 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]; + eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i]; + eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j]; + + Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; + + Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; + + Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; + Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; + Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; + } + } + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/ ElementVector* Penta::CreatePVectorAdjointHoriz(void){ Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11447) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11448) @@ -216,4 +216,5 @@ ElementMatrix* CreateKMatrixCouplingPattynStokes(void); ElementMatrix* CreateKMatrixDiagnosticHoriz(void); + ElementMatrix* CreateKMatrixAdjointHoriz(void); ElementVector* CreateDVectorDiagnosticHoriz(void); ElementVector* CreateDVectorDiagnosticStokes(void); @@ -275,4 +276,7 @@ #ifdef _HAVE_CONTROL_ ElementVector* CreatePVectorAdjointHoriz(void); + ElementMatrix* CreateKMatrixAdjointMacAyeal2d(void); + ElementMatrix* CreateKMatrixAdjointPattyn(void); + ElementMatrix* CreateKMatrixAdjointStokes(void); ElementVector* CreatePVectorAdjointMacAyeal(void); ElementVector* CreatePVectorAdjointPattyn(void); Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11447) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11448) @@ -336,6 +336,9 @@ switch(analysis_type){ #ifdef _HAVE_DIAGNOSTIC_ - case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: + case DiagnosticHorizAnalysisEnum: Ke=CreateKMatrixDiagnosticMacAyeal(); + break; + case AdjointHorizAnalysisEnum: + Ke=CreateKMatrixAdjointMacAyeal(); break; case DiagnosticHutterAnalysisEnum: @@ -3063,5 +3066,5 @@ } /*}}}*/ -/*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/ +/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{1*/ ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){ @@ -3864,5 +3867,5 @@ /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ - //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight; + Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight; } @@ -4736,5 +4739,5 @@ /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ - //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight; + Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight; } @@ -4763,4 +4766,74 @@ /*Transpose and return Ke*/ Ke->Transpose(); + return Ke; +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixAdjointMacAyeal{{{1*/ +ElementMatrix* Tria::CreateKMatrixAdjointMacAyeal(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + bool incomplete_adjoint; + 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)*/ + parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); + if(incomplete_adjoint) return Ke; + + /*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; + //Ke->Transpose(); return Ke; } Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11447) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11448) @@ -217,4 +217,5 @@ #ifdef _HAVE_CONTROL_ ElementMatrix* CreateKMatrixAdjointBalancethickness(void); + ElementMatrix* CreateKMatrixAdjointMacAyeal(void); ElementVector* CreatePVectorAdjointHoriz(void); ElementVector* CreatePVectorAdjointStokes(void);