Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 26988) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 26989) @@ -226,5 +226,4 @@ /*Intermediaries */ - bool incomplete_adjoint; IssmDouble Jdet,mu_prime,n,thickness,mu,effmu; IssmDouble *xyz_list = NULL; @@ -238,5 +237,137 @@ IssmDouble vxshear, vyshear; - Element* basalelement; + /*Get basal element*/ + Element* basalelement = NULL; + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: case Domain2DverticalEnum: + _error_("mesh "<FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); + ElementMatrix* Ke=analysis->CreateKMatrix(basalelement); + delete analysis; + if(incomplete_adjoint) return Ke; + + /*Second part of the Gateau derivative if requested*/ + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = basalelement->GetNumberOfNodes(); + IssmDouble* dbasis = xNew(2*numnodes); // like SSA + IssmDouble* basis = xNew(numnodes); // like SSA + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); //vertically integrated vx + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); //vertically integrated vy + Input* vxbase_input = element->GetInput(VxBaseEnum); _assert_(vxbase_input); + Input* vybase_input = element->GetInput(VyBaseEnum); _assert_(vybase_input); + Input* vxshear_input = element->GetInput(VxShearEnum); _assert_(vxshear_input); + Input* vyshear_input = element->GetInput(VyShearEnum); _assert_(vyshear_input); + Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + Gauss* gauss_base = basalelement->NewGauss(); + while(gauss->next()){ + gauss->SynchronizeGaussBase(gauss_base); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); + basalelement->NodalFunctions(basis, gauss_base); + + thickness_input->GetInputValue(&thickness, gauss); + n_input->GetInputValue(&n,gauss); + + vxshear_input->GetInputValue(&vxshear,gauss); + vyshear_input->GetInputValue(&vyshear,gauss); + + element->material->ViscosityMLHOAdjoint(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input); + + effmu = 2.0*(1-n)/2.0/n; + + element->StrainRateHO(&epsilonbase[0],xyz_list,gauss,vxbase_input, vybase_input); + element->StrainRateHO(&epsilonshear[0],xyz_list,gauss,vxshear_input, vyshear_input); + + e1b[0] = 2*epsilonbase[0]+epsilonbase[1]; e1b[1] = epsilonbase[2]; + e2b[1] = epsilonbase[0]+2*epsilonbase[1]; e2b[0] = epsilonbase[2]; + + e1sh[0] = 2*epsilonshear[0]+epsilonshear[1]; e1sh[1] = epsilonshear[2]; + e2sh[1] = epsilonshear[0]+2*epsilonshear[1]; e2sh[0] = epsilonshear[2]; + + for(int i=0;ivalues[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb1i+viscosity[1]*(esh1j*eb1i+eb1j*esh1i)+viscosity[2]*esh1j*esh1i); + Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1j*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb2i+viscosity[1]*(esh1j*eb2i+eb1j*esh2i)+viscosity[2]*esh1j*esh2i); + Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1j*basis[i]); + + Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1i*basis[j]); + Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb1i+viscosity[4]*(esh1j*eb1i+eb1j*esh1i)+viscosity[5]*esh1j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(eb1i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(esh1i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2i*basis[j]); + Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb2i+viscosity[4]*(esh1j*eb2i+eb1j*esh2i)+viscosity[5]*esh1j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*eb2i*basis[j]+vyshear*eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2i*basis[j]+vyshear*esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vyshear*basis[j]*basis[i]); + + + Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb1i+viscosity[1]*(esh2j*eb1i+eb2j*esh1i)+viscosity[2]*esh2j*esh1i); + Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2j*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb2i+viscosity[1]*(esh2j*eb2i+eb2j*esh2i)+viscosity[2]*esh2j*esh2i); + Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2j*basis[i]); + + + Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1i*basis[j]); + Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb1i+viscosity[4]*(esh2j*eb1i+eb2j*esh1i)+viscosity[5]*esh2j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vyshear*eb1i*basis[j]+vxshear*eb2j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2j*basis[i]+vyshear*esh1i*basis[j])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2i*basis[j]); + Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb2i+viscosity[4]*(esh2j*eb2i+eb2j*esh2i)+viscosity[5]*esh2j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(eb2i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(esh2i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vyshear*vyshear*basis[j]*basis[i]); + } + } + } + + /*Transform Coordinate System*/ + /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/ + + /*Clean up and return*/ + delete gauss; + delete gauss_base; + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; + xDelete(xyz_list); + xDelete(dbasis); + xDelete(basis); + return Ke; +}/*}}}*/ +ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHOVerticalIntergrated(Element* element){/*{{{*/ + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*Intermediaries */ + IssmDouble Jdet,mu_prime,n,thickness,mu,effmu; + IssmDouble *xyz_list = NULL; + IssmDouble zeta, epsilon_eff; + int domaintype; + int dim=2; + + IssmDouble e1phi1i, e1phi1j, e2phi1i, e2phi1j, e1phi2i, e1phi2j, e2phi2i, e2phi2j; + IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble e1[3],e2[3]; + + Element* basalelement = NULL; /*Get basal element*/ @@ -251,10 +382,7 @@ default: _error_("mesh "<GetNumberOfNodes(); - IssmDouble* dbasis = xNew(2*numnodes); // like SSA - IssmDouble* basis = xNew(numnodes); // like SSA /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ + bool incomplete_adjoint; element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); @@ -263,133 +391,11 @@ if(incomplete_adjoint) return Ke; - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); //vertically integrated vx - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); //vertically integrated vy - Input* vxbase_input = element->GetInput(VxBaseEnum); _assert_(vxbase_input); - Input* vybase_input = element->GetInput(VyBaseEnum); _assert_(vybase_input); - Input* vxshear_input = element->GetInput(VxShearEnum); _assert_(vxshear_input); - Input* vyshear_input = element->GetInput(VyShearEnum); _assert_(vyshear_input); - Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss = element->NewGauss(5); - Gauss* gauss_base = basalelement->NewGauss(); - while(gauss->next()){ - gauss->SynchronizeGaussBase(gauss_base); - - element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); - basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); - basalelement->NodalFunctions(basis, gauss_base); - - thickness_input->GetInputValue(&thickness, gauss); - n_input->GetInputValue(&n,gauss); - - vxshear_input->GetInputValue(&vxshear,gauss); - vyshear_input->GetInputValue(&vyshear,gauss); - - element->material->ViscosityMLHOAdjoint(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input); - - effmu = 2.0*(1-n)/2.0/n; - - element->StrainRateHO(&epsilonbase[0],xyz_list,gauss,vxbase_input, vybase_input); - element->StrainRateHO(&epsilonshear[0],xyz_list,gauss,vxshear_input, vyshear_input); - - e1b[0] = 2*epsilonbase[0]+epsilonbase[1]; e1b[1] = epsilonbase[2]; - e2b[1] = epsilonbase[0]+2*epsilonbase[1]; e2b[0] = epsilonbase[2]; - - e1sh[0] = 2*epsilonshear[0]+epsilonshear[1]; e1sh[1] = epsilonshear[2]; - e2sh[1] = epsilonshear[0]+2*epsilonshear[1]; e2sh[0] = epsilonshear[2]; - - for(int i=0;ivalues[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb1i+viscosity[1]*(esh1j*eb1i+eb1j*esh1i)+viscosity[2]*esh1j*esh1i); - Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1j*basis[i]); - Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb2i+viscosity[1]*(esh1j*eb2i+eb1j*esh2i)+viscosity[2]*esh1j*esh2i); - Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1j*basis[i]); - - Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1i*basis[j]); - Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb1i+viscosity[4]*(esh1j*eb1i+eb1j*esh1i)+viscosity[5]*esh1j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(eb1i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(esh1i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); - Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2i*basis[j]); - Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb2i+viscosity[4]*(esh1j*eb2i+eb1j*esh2i)+viscosity[5]*esh1j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*eb2i*basis[j]+vyshear*eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2i*basis[j]+vyshear*esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vyshear*basis[j]*basis[i]); - - - Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb1i+viscosity[1]*(esh2j*eb1i+eb2j*esh1i)+viscosity[2]*esh2j*esh1i); - Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2j*basis[i]); - Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb2i+viscosity[1]*(esh2j*eb2i+eb2j*esh2i)+viscosity[2]*esh2j*esh2i); - Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2j*basis[i]); - - - Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1i*basis[j]); - Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb1i+viscosity[4]*(esh2j*eb1i+eb2j*esh1i)+viscosity[5]*esh2j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vyshear*eb1i*basis[j]+vxshear*eb2j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2j*basis[i]+vyshear*esh1i*basis[j])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); - Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2i*basis[j]); - Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb2i+viscosity[4]*(esh2j*eb2i+eb2j*esh2i)+viscosity[5]*esh2j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(eb2i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(esh2i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vyshear*vyshear*basis[j]*basis[i]); - } - } - } - - /*Transform Coordinate System*/ - /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/ - - /*Clean up and return*/ - delete gauss; - delete gauss_base; - if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; - xDelete(xyz_list); - xDelete(dbasis); - xDelete(basis); - return Ke; -}/*}}}*/ -ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHOVerticalIntergrated(Element* element){/*{{{*/ - - /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - - /*Intermediaries */ - bool incomplete_adjoint; - IssmDouble Jdet,mu_prime,n,thickness,mu,effmu; - IssmDouble *xyz_list = NULL; - IssmDouble zeta, epsilon_eff; - int domaintype; - int dim=2; - - IssmDouble e1phi1i, e1phi1j, e2phi1i, e2phi1j, e1phi2i, e1phi2j, e2phi2i, e2phi2j; - IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble e1[3],e2[3]; - - Element* basalelement; - - /*Get basal element*/ - element->FindParam(&domaintype,DomainTypeEnum); - switch(domaintype){ - case Domain2DhorizontalEnum: - basalelement = element; - break; - case Domain3DEnum: case Domain2DverticalEnum: - _error_("mesh "<GetNumberOfNodes(); IssmDouble* dbasis = xNew(2*numnodes); // like SSA IssmDouble* basis = xNew(numnodes); // like SSA - - /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ - element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); - StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); - ElementMatrix* Ke=analysis->CreateKMatrix(basalelement); - delete analysis; - if(incomplete_adjoint) return Ke; /*Retrieve all inputs and parameters*/ @@ -3428,12 +3434,12 @@ /*Fetch dof list and allocate solution vectors*/ element->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); - IssmDouble* values = xNew(numdof); - IssmDouble* lambdax = xNew(numnodes); - IssmDouble* lambday = xNew(numnodes); - IssmDouble* lambdabx = xNew(numnodes); - IssmDouble* lambdaby = xNew(numnodes); + IssmDouble* values = xNew(numdof); + IssmDouble* lambdax = xNew(numnodes); + IssmDouble* lambday = xNew(numnodes); + IssmDouble* lambdabx = xNew(numnodes); + IssmDouble* lambdaby = xNew(numnodes); IssmDouble* lambdashx = xNew(numnodes); IssmDouble* lambdashy = xNew(numnodes); - IssmDouble* n = xNew(numnodes); + IssmDouble* n = xNew(numnodes); /*Use the dof list to index into the solution vector: */ @@ -3485,4 +3491,5 @@ xDelete(lambdashx); xDelete(lambdashy); + xDelete(n); xDelete(doflist); }/*}}}*/