Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18476) @@ -45,106 +45,6 @@ }/*}}}*/ ElementVector* AdjointBalancethickness2Analysis::CreatePVector(Element* element){/*{{{*/ + _error_("Not implemented yet"); - /*Intermediaries */ - int num_responses,i; - IssmDouble hobs,hu2,weight,Jdet; - IssmDouble NUMxH2,NUMyH2,DENH2; - IssmDouble NUMxUbar,NUMyUbar,DENUbar; - IssmDouble vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs; - IssmDouble nu,phi,dphi[2]; - int *responses = NULL; - IssmDouble *xyz_list = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Initialize Element vector and vectors*/ - ElementVector* pe = element->NewElementVector(SSAApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* dbasis = xNew(2*numnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&num_responses,InversionNumCostFunctionsEnum); - element->FindParam(&responses,NULL,InversionCostFunctionsEnum); - Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); - Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* potential_input = element->GetInput(PotentialEnum); _assert_(potential_input); - Input* vxobs_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input); - Input* nu_input = element->GetInput(BalancethicknessNuEnum); _assert_(nu_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctions(basis,gauss); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - nu_input->GetInputValue(&nu,gauss); - potential_input->GetInputValue(&phi,gauss); - potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss); - thicknessobs_input->GetInputValue(&hobs,gauss); - - vxobsbar = nu*vxobs; - vyobsbar = nu*vyobs; - - vbarobs2 = (nu*nu*vxobs*vxobs + nu*nu*vyobs*vyobs); - vbarobs = sqrt(vbarobs2); - hu2 = hobs*hobs*vbarobs2; - - /*Loop over all requested responses*/ - for(int resp=0;respGetInputValue(&weight,gauss,responses[resp]); - - switch(responses[resp]){ - case IrrotationalH2MisfitEnum: - /*J = (H^2 - Hobs^2)^2*/ - for(i=0;ivalues[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight; - } - break; - case IrrotationalDirectionMisfitEnum: - /*J = 1/2 (vbar ^ gard(phi))^2*/ - for(i=0;ivalues[i]+= weight*Jdet*gauss->weight* - nu*nu*(vyobs*dphi[0] - vxobs*dphi[1])* - (vyobs*dbasis[0*numnodes+i] - vxobs*dbasis[1*numnodes+i]); - } - break; - case Balancethickness2MisfitEnum: - /*J = phi^2*/ - //for(i=0;ivalues[i]+= phi*basis[i]*weight*Jdet*gauss->weight; //OK - /*J = grad phi ^2*/ - //for(i=0;ivalues[i]+= (dphi[0]*dbasis[0*numnodes+i] + dphi[1]*dbasis[1*numnodes+i])*weight*Jdet*gauss->weight; //OK - /*J = (ubar - nu*uobs)^2*/ - //for(i=0;ivalues[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight; - //} - _error_("Not implemented yet"); - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - } - } - - /*Clean up and return*/ - xDelete(responses); - xDelete(xyz_list); - xDelete(basis); - xDelete(dbasis); - delete gauss; - return pe; }/*}}}*/ void AdjointBalancethickness2Analysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ @@ -174,25 +74,4 @@ /*Nothing, \partial J/\partial k = 0*/ break; - case IrrotationalH2MisfitEnum: - switch(control_type){ - case BalancethicknessApparentMassbalanceEnum: - /*Nothing, \partial J/\partial k = 0*/ - break; - case BalancethicknessNuEnum: - GradientJ1Nu(element,gradient,control_index); - break; - default: _error_("control type not supported yet: " << EnumToStringx(control_type)); - } - break; - case IrrotationalDirectionMisfitEnum: - switch(control_type){ - case BalancethicknessApparentMassbalanceEnum: - /*Nothing, \partial J/\partial k = 0*/ - break; - case BalancethicknessNuEnum: - GradientJ2Nu(element,gradient,control_index); break; //Might need to be renamed ? - default: _error_("control type not supported yet: " << EnumToStringx(control_type)); - } - break; default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); } @@ -200,6 +79,4 @@ /*Deal with second term*/ switch(control_type){ - case BalancethicknessApparentMassbalanceEnum: GradientJAdot(element,gradient,control_index); break; - case BalancethicknessNuEnum: /*Nothing: state equations do not depend on k*/ break; default: _error_("control type not supported yet: " << EnumToStringx(control_type)); } @@ -251,118 +128,4 @@ delete gauss; }/*}}}*/ -void AdjointBalancethickness2Analysis::GradientJ1Nu(Element* element,Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - IssmDouble vxobs,vyobs,thicknessobs,weight; - IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2; - IssmDouble *xyz_list= NULL; - - /*Fetch number of vertices for this finite element*/ - int numvertices = element->GetNumberOfVertices(); - - /*Initialize some vectors*/ - IssmDouble* basis = xNew(numvertices); - IssmDouble* ge = xNewZeroInit(numvertices); - int* vertexpidlist = xNew(numvertices); - - /*Retrieve all inputs we will be needing: */ - element->GetVerticesCoordinates(&xyz_list); - element->GradientIndexing(&vertexpidlist[0],control_index); - Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vxobs_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input); - Input* nu_input = element->GetInput(BalancethicknessNuEnum); _assert_(nu_input); - Input* potential_input = element->GetInput(PotentialEnum); _assert_(potential_input); - Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input); - - Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctionsP1(basis,gauss); - weights_input->GetInputValue(&weight,gauss,IrrotationalH2MisfitEnum); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - nu_input->GetInputValue(&nu,gauss); - potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss); - thicknessobs_input->GetInputValue(&thicknessobs,gauss); - - dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1]; - nuvobs2 = nu*nu*(vxobs*vxobs + vyobs*vyobs); - - /*Build gradient vector (actually -dJ/da): */ - for(int i=0;iweight*basis[i]*( - -2./pow(nu,3) * dphinorm2/(vxobs*vxobs + vyobs*vyobs) * (dphinorm2/nuvobs2 - thicknessobs*thicknessobs) - ); - _assert_(!xIsNan(ge[i])); - } - } - gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); - - /*Clean up and return*/ - xDelete(ge); - xDelete(xyz_list); - xDelete(basis); - xDelete(vertexpidlist); - delete gauss; -}/*}}}*/ -void AdjointBalancethickness2Analysis::GradientJ2Nu(Element* element,Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - IssmDouble vxobs,vyobs,thicknessobs,weight; - IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2; - IssmDouble *xyz_list= NULL; - - /*Fetch number of vertices for this finite element*/ - int numvertices = element->GetNumberOfVertices(); - - /*Initialize some vectors*/ - IssmDouble* basis = xNew(numvertices); - IssmDouble* ge = xNewZeroInit(numvertices); - int* vertexpidlist = xNew(numvertices); - - /*Retrieve all inputs we will be needing: */ - element->GetVerticesCoordinates(&xyz_list); - element->GradientIndexing(&vertexpidlist[0],control_index); - Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vxobs_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input); - Input* nu_input = element->GetInput(BalancethicknessNuEnum); _assert_(nu_input); - Input* potential_input = element->GetInput(PotentialEnum); _assert_(potential_input); - Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input); - - Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctionsP1(basis,gauss); - weights_input->GetInputValue(&weight,gauss,IrrotationalDirectionMisfitEnum); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - nu_input->GetInputValue(&nu,gauss); - potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss); - thicknessobs_input->GetInputValue(&thicknessobs,gauss); - - dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1]; - nuvobs2 = nu*vxobs*nu*vxobs + nu*vyobs*nu*vyobs; - - /*Build gradient vector (actually -dJ/da): */ - for(int i=0;iweight*nu*(vyobs*dphi[0] -vxobs*dphi[1])*(vyobs*dphi[0] -vxobs*dphi[1])*basis[i]; - _assert_(!xIsNan(ge[i])); - } - } - gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); - - /*Clean up and return*/ - xDelete(ge); - xDelete(xyz_list); - xDelete(basis); - xDelete(vertexpidlist); - delete gauss; -}/*}}}*/ void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ element->InputUpdateFromSolutionOneDof(solution,AdjointEnum); Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h (revision 18475) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h (revision 18476) @@ -29,6 +29,4 @@ void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); void GradientJAdot(Element* element,Vector* gradient,int control_index); - void GradientJ1Nu(Element* element,Vector* gradient,int control_index); - void GradientJ2Nu(Element* element,Vector* gradient,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18476) @@ -26,11 +26,11 @@ } - iomodel->FetchDataToInput(elements,BalancethicknessApparentMassbalanceEnum); - iomodel->FetchDataToInput(elements,BalancethicknessNuEnum); - iomodel->FetchDataToInput(elements,BalancethicknessVxObsEnum); - iomodel->FetchDataToInput(elements,BalancethicknessVyObsEnum); - iomodel->FetchDataToInput(elements,BalancethicknessThicknessObsEnum); - iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum); + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BaseEnum); iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum); }/*}}}*/ void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ @@ -43,5 +43,5 @@ int finiteelement = P1Enum; - IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcpotentialEnum,Balancethickness2AnalysisEnum,finiteelement); + IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,Balancethickness2AnalysisEnum,finiteelement); }/*}}}*/ @@ -64,5 +64,5 @@ /*Intermediaries */ - IssmDouble Jdet,D_scalar; + IssmDouble Jdet,D; IssmDouble* xyz_list = NULL; @@ -74,6 +74,10 @@ IssmDouble* dbasis = xNew(2*numnodes); + /*Create input D*/ + this->CreateDiffusionCoefficient(element); + /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); + Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); _assert_(D_input); /* Start looping on the number of gaussian points: */ @@ -83,8 +87,9 @@ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); element->JacobianDeterminant(&Jdet,xyz_list,gauss); + D_input->GetInputValue(&D,gauss); for(int i=0;ivalues[i*numnodes+j] += gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); + Ke->values[i*numnodes+j] += D*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); } } @@ -99,18 +104,6 @@ ElementVector* Balancethickness2Analysis::CreatePVector(Element* element){/*{{{*/ - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorVolume(element); - ElementVector* pe2=CreatePVectorBoundary(element); - ElementVector* pe =new ElementVector(pe1,pe2); - - /*clean-up and return*/ - delete pe1; - delete pe2; - return pe; -}/*}}}*/ -ElementVector* Balancethickness2Analysis::CreatePVectorVolume(Element* element){/*{{{*/ - /*Intermediaries */ - IssmDouble adot,Jdet; + IssmDouble dhdt,mb,ms,dD[2],db[2],Jdet; IssmDouble* xyz_list = NULL; @@ -124,5 +117,9 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* adot_input = element->GetInput(BalancethicknessApparentMassbalanceEnum); _assert_(adot_input); + Input* ms_input = element->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); + Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input); + Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); + Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); _assert_(D_input); + Input* bed_input = element->GetInput(BaseEnum); _assert_(bed_input); /* Start looping on the number of gaussian points: */ @@ -133,7 +130,15 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - adot_input->GetInputValue(&adot,gauss); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight*adot*basis[i]; + + ms_input->GetInputValue(&ms,gauss); + mb_input->GetInputValue(&mb,gauss); + dhdt_input->GetInputValue(&dhdt,gauss); + bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss); + D_input->GetInputDerivativeValue(&dD[0],xyz_list,gauss); + db[0]=0.; + db[1]=0.; + + /*Since grad(b) is constant div(D grad(b) ) = grad(D).grad(b)*/ + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(ms-mb-dhdt +dD[0]*db[0]+dD[1]*db[1])*basis[i]; } @@ -144,52 +149,4 @@ return pe; }/*}}}*/ -ElementVector* Balancethickness2Analysis::CreatePVectorBoundary(Element* element){/*{{{*/ - - /*If no front, return NULL*/ - if(!element->IsFaceOnBoundary()) return NULL; - - /*Intermediaries*/ - IssmDouble Jdet,thickness,vx,vy; - IssmDouble *xyz_list = NULL; - IssmDouble *xyz_list_front = NULL; - IssmDouble normal[2]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = element->NewElementVector(); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - Input* thickness_input = element->GetInput(BalancethicknessThicknessObsEnum); _assert_(thickness_input); - Input* vx_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input); - - element->GetVerticesCoordinates(&xyz_list); - element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); - element->NormalSection(&normal[0],xyz_list_front); - - /*Start looping on Gaussian points*/ - Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - thickness_input->GetInputValue(&thickness,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); - element->NodalFunctions(basis,gauss); - - for(int i=0;ivalues[i] += Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i]; - } - - /*Clean up and return*/ - xDelete(xyz_list); - xDelete(xyz_list_front); - xDelete(basis); - delete gauss; - return pe; -}/*}}}*/ void Balancethickness2Analysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ _error_("not implemented yet"); @@ -199,75 +156,5 @@ }/*}}}*/ void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - - /*Intermediaries */ - int Hinterpolation; - IssmDouble vx,vy,vbar,nu,normdphi,dphi[2],H; - IssmDouble* xyz_list = NULL; - int * doflist = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - int numvertices = element->GetNumberOfVertices(); - - /*Fetch dof list and allocate solution vector*/ - element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - IssmDouble* values = xNew(numnodes); - IssmDouble* thickness_list = xNew(numvertices); - IssmDouble* vx_list = xNew(numvertices); - IssmDouble* vy_list = xNew(numvertices); - - /*Use the dof list to index into the solution vector: */ - for(int i=0;i(values[i])) _error_("NaN found in solution vector"); - } - - element->AddInput(PotentialEnum,values,element->GetElementType()); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - Input* potential_input = element->GetInput(PotentialEnum); _assert_(potential_input); - Input* vx_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input); - Input* nu_input = element->GetInput(BalancethicknessNuEnum); _assert_(nu_input); - Input* thickness_input = element->GetInput(BalancethicknessThicknessObsEnum); _assert_(thickness_input); - - switch(element->GetElementType()){ - case P1Enum: Hinterpolation = P0Enum; break; - default: _error_("not implemented"); - } - - Gauss* gauss=element->NewGauss(); - for (int iv=0;iv<1;iv++){ - gauss->GaussNode(Hinterpolation,iv);//P0 Only for now - - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - nu_input->GetInputValue(&nu,gauss); - thickness_input->GetInputValue(&H,gauss); - potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss); - - vx = vx*nu; vy = vy*nu; - vbar = sqrt(vx*vx + vy*vy) + 1.e-10; - normdphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]); - - thickness_list[iv] = normdphi/vbar; - vx_list[iv] = -1./H * dphi[0]; - vy_list[iv] = -1./H * dphi[1]; - } - element->AddInput(ThicknessEnum,thickness_list,Hinterpolation); - element->AddInput(VxEnum,vx_list,Hinterpolation); - element->AddInput(VyEnum,vy_list,Hinterpolation); - - /*Clean up and return*/ - delete gauss; - xDelete(doflist); - xDelete(xyz_list); - xDelete(values); - xDelete(thickness_list); - xDelete(vx_list); - xDelete(vy_list); + element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum); }/*}}}*/ void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ @@ -275,2 +162,48 @@ return; }/*}}}*/ + +/*Specifics*/ +void Balancethickness2Analysis::CreateDiffusionCoefficient(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble omega,h,mu0,ds[2],Cmu,B; + const int n = 3.; + const IssmDouble Hstar = 500.; + const IssmDouble Lstar = 500.e+3; + IssmDouble *xyz_list = NULL; + + /*Fetch number of vertices and allocate output*/ + int numnodes = element->GetNumberOfNodes(); + IssmDouble* D = xNew(numnodes); + IssmDouble* Dgradb = xNew(numnodes); + + /*retrieve what we need: */ + element->GetVerticesCoordinates(&xyz_list); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); + + /*Calculate damage evolution source term: */ + Gauss* gauss=element->NewGauss(); + for (int i=0;iGaussNode(element->GetElementType(),i); + + B_input->GetInputValue(&B,gauss); + thickness_input->GetInputValue(&h,gauss); + surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + + mu0 = pow(2.,(1.-3*n)/(2.*n)) * B; + omega = pow(rhog,3) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n); + Cmu = 0.; + + D[i] = omega*(ds[0]*ds[0]+ds[1]*ds[1])*pow(h,4)*(1./5.*h - Cmu); + } + + /*Add input*/ + element->AddInput(BalancethicknessDiffusionCoefficientEnum,D,element->GetElementType()); + + /*Clean up and return*/ + xDelete(D); + delete gauss; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h (revision 18475) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h (revision 18476) @@ -26,10 +26,11 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - ElementVector* CreatePVectorVolume(Element* element); - ElementVector* CreatePVectorBoundary(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); + + /*Specifics*/ + void CreateDiffusionCoefficient(Element* element); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18476) @@ -1053,6 +1053,4 @@ name==BedEnum || name==BalancethicknessThickeningRateEnum || - name==BalancethicknessApparentMassbalanceEnum || - name==BalancethicknessNuEnum || name==SigmaNNEnum || name==SurfaceSlopeXEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18476) @@ -1032,5 +1032,4 @@ case ThicknessEnum: case FrictionCoefficientEnum: - case BalancethicknessNuEnum: if(iomodel->Data(control)){ for(j=0;jData(control)[tria_vertex_ids[j]-1]; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18476) @@ -593,9 +593,4 @@ case BalancethicknessMisfitEnum: BalancethicknessMisfitx(&double_result); break; case Balancethickness2MisfitEnum: Balancethickness2Misfitx(&double_result); break; - case IrrotationalH2MisfitEnum: IrrotationalH2Misfitx(&double_result); break; - case IrrotationalDirectionMisfitEnum: IrrotationalDirectionMisfitx(&double_result); break; - case IrrotationalVelMisfitEnum: IrrotationalVelMisfitx(&double_result); break; - case IrrotationalAlongGradientNuEnum: IrrotationalAlongGradientNux(&double_result); break; - case IrrotationalAcrossGradientNuEnum: IrrotationalAcrossGradientNux(&double_result); break; /*Vector */ @@ -1320,56 +1315,6 @@ /*If on water, return 0: */ if(!element->IsIceInElement()) continue; - - /* Get node coordinates*/ - element->GetVerticesCoordinates(&xyz_list); - Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); - Input* potential_input =element->GetInput(PotentialEnum); _assert_(potential_input); - Input* vxobs_input =element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input =element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input); - Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); - Input* nu_input = element->GetInput(BalancethicknessNuEnum); _assert_(nu_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum); - thickness_input->GetInputValue(&thickness,gauss); - thicknessobs_input->GetInputValue(&thicknessobs,gauss); - potential_input->GetInputValue(&potential,gauss); - potential_input->GetInputDerivativeValue(&dpotential[0],xyz_list,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - vx_input->GetInputValue(&vxbar,gauss); - vy_input->GetInputValue(&vybar,gauss); - nu_input->GetInputValue(&nu,gauss); - - vxbarobs = nu*vxobs; - vybarobs = nu*vyobs; - - /*J = (H^2 - Hobs^2)^2*/ - J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight; - /*J = phi^2*/ - //J +=.5*potential*potential*weight*Jdet*gauss->weight;// OK - /*J = grad phi^2*/ - //J +=.5*(dpotential[0]*dpotential[0] + dpotential[1]*dpotential[1])*weight*Jdet*gauss->weight; - /*J = (ubar - nux*uobs)^2*/ - //J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight; - /*J = 1/2 (vbar ^ gard(phi))^2*/ - //J += 0.5*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*weight*Jdet*gauss->weight; - } - - /*clean up and Return: */ - xDelete(xyz_list); - delete gauss; + _error_("Not implemented"); + } @@ -1382,124 +1327,4 @@ *presponse=J; -}/*}}}*/ -void FemModel::IrrotationalH2Misfitx(IssmDouble* presponse){/*{{{*/ - - /*output: */ - IssmDouble J=0.; - IssmDouble J_sum; - - IssmDouble weight,thicknessobs,thickness; - IssmDouble Jdet; - IssmDouble* xyz_list = NULL; - - /*Compute Misfit: */ - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - - /*If on water, return 0: */ - if(!element->IsIceInElement()) continue; - - /* Get node coordinates*/ - element->GetVerticesCoordinates(&xyz_list); - Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,IrrotationalH2MisfitEnum); - thickness_input->GetInputValue(&thickness,gauss); - thicknessobs_input->GetInputValue(&thicknessobs,gauss); - - /*J = (H^2 - Hobs^2)^2*/ - J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight; - } - - /*clean up and Return: */ - xDelete(xyz_list); - delete gauss; - } - - /*Sum all J from all cpus of the cluster:*/ - ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - J=J_sum; - - /*Assign output pointers: */ - *presponse=J; - -}/*}}}*/ -void FemModel::IrrotationalDirectionMisfitx(IssmDouble* presponse){/*{{{*/ - - /*output: */ - IssmDouble J=0.; - IssmDouble J_sum; - - IssmDouble weight,thicknessobs,thickness,dpotential[2]; - IssmDouble vx,vy,vxobs,vyobs,nu; - IssmDouble Jdet; - IssmDouble* xyz_list = NULL; - - /*Compute Misfit: */ - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - - /*If on water, return 0: */ - if(!element->IsIceInElement()) continue; - - /* Get node coordinates*/ - element->GetVerticesCoordinates(&xyz_list); - Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* potential_input =element->GetInput(PotentialEnum); _assert_(potential_input); - Input* vxobs_input =element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input =element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input); - Input* nu_input = element->GetInput(BalancethicknessNuEnum); _assert_(nu_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,IrrotationalDirectionMisfitEnum); - potential_input->GetInputDerivativeValue(&dpotential[0],xyz_list,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - nu_input->GetInputValue(&nu,gauss); - - /*J = 1/2 (vbar ^ gard(phi))^2*/ - J += 0.5*(nu*vyobs*dpotential[0] - nu*vxobs*dpotential[1])*(nu*vyobs*dpotential[0] - nu*vxobs*dpotential[1])*weight*Jdet*gauss->weight; - } - - /*clean up and Return: */ - xDelete(xyz_list); - delete gauss; - } - - /*Sum all J from all cpus of the cluster:*/ - ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - J=J_sum; - - /*Assign output pointers: */ - *presponse=J; - -}/*}}}*/ -void FemModel::IrrotationalVelMisfitx(IssmDouble* presponse){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void FemModel::IrrotationalAlongGradientNux(IssmDouble* presponse){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void FemModel::IrrotationalAcrossGradientNux(IssmDouble* presponse){/*{{{*/ - _error_("not implemented yet"); }/*}}}*/ void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18475) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18476) @@ -81,9 +81,4 @@ void BalancethicknessMisfitx(IssmDouble* pV); void Balancethickness2Misfitx(IssmDouble* pV); - void IrrotationalH2Misfitx(IssmDouble* pV); - void IrrotationalDirectionMisfitx(IssmDouble* pV); - void IrrotationalVelMisfitx(IssmDouble* pV); - void IrrotationalAlongGradientNux(IssmDouble* pV); - void IrrotationalAcrossGradientNux(IssmDouble* pV); #ifdef _HAVE_DAKOTA_ void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses); Index: /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18476) @@ -26,6 +26,6 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - int outputs[4] = {ThicknessEnum,PotentialEnum,VxEnum,VyEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4); + int outputs[1] = {ThicknessEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18476) @@ -41,5 +41,4 @@ case ThicknessEnum: case FrictionCoefficientEnum: - case BalancethicknessNuEnum: case BalancethicknessApparentMassbalanceEnum: iomodel->FetchData(1,control); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18475) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18476) @@ -301,13 +301,5 @@ BalancethicknessApparentMassbalanceEnum, Balancethickness2MisfitEnum, - IrrotationalH2MisfitEnum, - IrrotationalDirectionMisfitEnum, - IrrotationalVelMisfitEnum, - IrrotationalAlongGradientNuEnum, - IrrotationalAcrossGradientNuEnum, - BalancethicknessNuEnum, - BalancethicknessVxObsEnum, - BalancethicknessVyObsEnum, - BalancethicknessThicknessObsEnum, + BalancethicknessDiffusionCoefficientEnum, /*}}}*/ /*Surfaceforcings{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18476) @@ -309,13 +309,5 @@ case BalancethicknessApparentMassbalanceEnum : return "BalancethicknessApparentMassbalance"; case Balancethickness2MisfitEnum : return "Balancethickness2Misfit"; - case IrrotationalH2MisfitEnum : return "IrrotationalH2Misfit"; - case IrrotationalDirectionMisfitEnum : return "IrrotationalDirectionMisfit"; - case IrrotationalVelMisfitEnum : return "IrrotationalVelMisfit"; - case IrrotationalAlongGradientNuEnum : return "IrrotationalAlongGradientNu"; - case IrrotationalAcrossGradientNuEnum : return "IrrotationalAcrossGradientNu"; - case BalancethicknessNuEnum : return "BalancethicknessNu"; - case BalancethicknessVxObsEnum : return "BalancethicknessVxObs"; - case BalancethicknessVyObsEnum : return "BalancethicknessVyObs"; - case BalancethicknessThicknessObsEnum : return "BalancethicknessThicknessObs"; + case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient"; case SurfaceforcingsEnum : return "Surfaceforcings"; case SMBEnum : return "SMB"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18475) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18476) @@ -315,13 +315,5 @@ else if (strcmp(name,"BalancethicknessApparentMassbalance")==0) return BalancethicknessApparentMassbalanceEnum; else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum; - else if (strcmp(name,"IrrotationalH2Misfit")==0) return IrrotationalH2MisfitEnum; - else if (strcmp(name,"IrrotationalDirectionMisfit")==0) return IrrotationalDirectionMisfitEnum; - else if (strcmp(name,"IrrotationalVelMisfit")==0) return IrrotationalVelMisfitEnum; - else if (strcmp(name,"IrrotationalAlongGradientNu")==0) return IrrotationalAlongGradientNuEnum; - else if (strcmp(name,"IrrotationalAcrossGradientNu")==0) return IrrotationalAcrossGradientNuEnum; - else if (strcmp(name,"BalancethicknessNu")==0) return BalancethicknessNuEnum; - else if (strcmp(name,"BalancethicknessVxObs")==0) return BalancethicknessVxObsEnum; - else if (strcmp(name,"BalancethicknessVyObs")==0) return BalancethicknessVyObsEnum; - else if (strcmp(name,"BalancethicknessThicknessObs")==0) return BalancethicknessThicknessObsEnum; + else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum; else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum; else if (strcmp(name,"SMB")==0) return SMBEnum; @@ -383,8 +375,5 @@ else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum; else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; + else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; @@ -394,5 +383,8 @@ else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum; - else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; @@ -506,8 +498,5 @@ else if (strcmp(name,"Converged")==0) return ConvergedEnum; else if (strcmp(name,"Fill")==0) return FillEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; + else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; else if (strcmp(name,"Friction")==0) return FrictionEnum; else if (strcmp(name,"Internal")==0) return InternalEnum; @@ -517,5 +506,8 @@ else if (strcmp(name,"Pressure")==0) return PressureEnum; else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum; - else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum; else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum; else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum; @@ -629,8 +621,5 @@ else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum; else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum; + else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum; else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum; else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum; @@ -640,5 +629,8 @@ else if (strcmp(name,"MinVel")==0) return MinVelEnum; else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; - else if (strcmp(name,"MinVx")==0) return MinVxEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"MinVx")==0) return MinVxEnum; else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; Index: /issm/trunk-jpl/src/dox/issm.dox =================================================================== --- /issm/trunk-jpl/src/dox/issm.dox (revision 18475) +++ /issm/trunk-jpl/src/dox/issm.dox (revision 18476) @@ -47,23 +47,17 @@ - C++ 453158611818868849102898 + C++ 445159671614570397102509 - MATLAB 13788341162593864563245 + MATLAB 14228497166383953864673 - C/C++ Header 421350437531548922746 + C/C++ Header 411344335281524322214 - m4 915881511156513304 + m4 81036149975610941 - Python 14423302559934814237 - - - XML 41609230753327 - - - Java 1871989123213931 + Python 15124212623973314777 @@ -71,17 +65,8 @@ - Bourne Shell 36188266415 + Bourne Shell 25984262405 - XSD 23036229295 - - - Perl 169196211 - - - Ant 1167103126 - - - SUM: 24413262042335150451225406 + SUM: 24463142739469145294216190 Index: sm/trunk-jpl/src/m/classes/balancethickness2.m =================================================================== --- /issm/trunk-jpl/src/m/classes/balancethickness2.m (revision 18475) +++ (revision ) @@ -1,63 +1,0 @@ -%BALANCETHICKNESS2 class definition -% -% Usage: -% balancethickness2=balancethickness2(); - -classdef balancethickness2 - properties (SetAccess=public) - spcpotential = NaN; - apparent_massbalance = NaN; - nu = NaN; - vx_obs = NaN; - vy_obs = NaN; - thickness_obs = NaN; - end - methods - function obj = balancethickness2(varargin) % {{{ - switch nargin - case 0 - obj=setdefaultparameters(obj); - otherwise - error('constructor not supported'); - end - end % }}} - function obj = setdefaultparameters(obj) % {{{ - - end % }}} - function md = checkconsistency(obj,md,solution,analyses) % {{{ - %Early return - if solution~=Balancethickness2SolutionEnum(), return; end - - md = checkfield(md,'fieldname','balancethickness.spcpotential','size',[md.mesh.numberofvertices 1]); - md = checkfield(md,'fieldname','balancethickness.apparent_massbalance','size',[md.mesh.numberofvertices 1],'NaN',1); - md = checkfield(md,'fieldname','balancethickness.nu','size',[md.mesh.numberofvertices 1],'NaN',1,'>',0,'<=',1); - md = checkfield(md,'fieldname','balancethickness.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1); - md = checkfield(md,'fieldname','balancethickness.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1); - md = checkfield(md,'fieldname','balancethickness.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0); - end % }}} - function disp(obj) % {{{ - disp(sprintf(' balance thickness solution parameters:')); - - fielddisplay(obj,'spcpotential','potential constraints (NaN means no constraint)'); - fielddisplay(obj,'apparent_massbalance','Apparent mass balance [m/yr]'); - fielddisplay(obj,'nu','v_bar = nu v_s (in ]0 1])'); - fielddisplay(obj,'vx_obs','observed vx'); - fielddisplay(obj,'vy_obs','observed vy'); - fielddisplay(obj,'thickness_obs','observed H'); - - end % }}} - function marshall(obj,md,fid) % {{{ - WriteData(fid,'object',obj,'class','Balancethickness','fieldname','spcpotential','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','Balancethickness','fieldname','apparent_massbalance','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','Balancethickness','fieldname','nu','format','DoubleMat','mattype',1); - if(numel(obj.vx_obs)==md.mesh.numberofelements), - mattype=2; - else - mattype=1; - end - WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vx_obs','format','DoubleMat','mattype',mattype); - WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vy_obs','format','DoubleMat','mattype',mattype); - WriteData(fid,'object',obj,'class','Balancethickness','fieldname','thickness_obs','format','DoubleMat','mattype',1); - end % }}} - end -end Index: /issm/trunk-jpl/src/m/classes/inversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.m (revision 18475) +++ /issm/trunk-jpl/src/m/classes/inversion.m (revision 18476) @@ -239,9 +239,4 @@ pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum(); pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum(); - pos=find(obj.cost_functions==601); data(pos)=IrrotationalH2MisfitEnum(); - pos=find(obj.cost_functions==602); data(pos)=IrrotationalDirectionMisfitEnum(); - pos=find(obj.cost_functions==603); data(pos)=IrrotationalVelMisfitEnum(); - pos=find(obj.cost_functions==604); data(pos)=IrrotationalAlongGradientNuEnum(); - pos=find(obj.cost_functions==605); data(pos)=IrrotationalAcrossGradientNuEnum(); WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3); WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer'); Index: /issm/trunk-jpl/src/m/classes/inversionvalidation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/inversionvalidation.m (revision 18475) +++ /issm/trunk-jpl/src/m/classes/inversionvalidation.m (revision 18476) @@ -140,9 +140,4 @@ pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum(); pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum(); - pos=find(obj.cost_functions==601); data(pos)=IrrotationalH2MisfitEnum(); - pos=find(obj.cost_functions==602); data(pos)=IrrotationalDirectionMisfitEnum(); - pos=find(obj.cost_functions==603); data(pos)=IrrotationalVelMisfitEnum(); - pos=find(obj.cost_functions==604); data(pos)=IrrotationalAlongGradientNuEnum(); - pos=find(obj.cost_functions==605); data(pos)=IrrotationalAcrossGradientNuEnum(); WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3); WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer'); Index: /issm/trunk-jpl/src/m/classes/m1qn3inversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/m1qn3inversion.m (revision 18475) +++ /issm/trunk-jpl/src/m/classes/m1qn3inversion.m (revision 18476) @@ -170,9 +170,4 @@ pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum(); pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum(); - pos=find(obj.cost_functions==601); data(pos)=IrrotationalH2MisfitEnum(); - pos=find(obj.cost_functions==602); data(pos)=IrrotationalDirectionMisfitEnum(); - pos=find(obj.cost_functions==603); data(pos)=IrrotationalVelMisfitEnum(); - pos=find(obj.cost_functions==604); data(pos)=IrrotationalAlongGradientNuEnum(); - pos=find(obj.cost_functions==605); data(pos)=IrrotationalAcrossGradientNuEnum(); WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3); WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer'); Index: /issm/trunk-jpl/src/m/enum/BalancethicknessDiffusionCoefficientEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessDiffusionCoefficientEnum.m (revision 18476) +++ /issm/trunk-jpl/src/m/enum/BalancethicknessDiffusionCoefficientEnum.m (revision 18476) @@ -0,0 +1,11 @@ +function macro=BalancethicknessDiffusionCoefficientEnum() +%BALANCETHICKNESSDIFFUSIONCOEFFICIENTENUM - Enum of BalancethicknessDiffusionCoefficient +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=BalancethicknessDiffusionCoefficientEnum() + +macro=StringToEnum('BalancethicknessDiffusionCoefficient'); Index: sm/trunk-jpl/src/m/enum/BalancethicknessNuEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessNuEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=BalancethicknessNuEnum() -%BALANCETHICKNESSNUENUM - Enum of BalancethicknessNu -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=BalancethicknessNuEnum() - -macro=StringToEnum('BalancethicknessNu'); Index: sm/trunk-jpl/src/m/enum/BalancethicknessThicknessObsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessThicknessObsEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=BalancethicknessThicknessObsEnum() -%BALANCETHICKNESSTHICKNESSOBSENUM - Enum of BalancethicknessThicknessObs -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=BalancethicknessThicknessObsEnum() - -macro=StringToEnum('BalancethicknessThicknessObs'); Index: sm/trunk-jpl/src/m/enum/BalancethicknessVxObsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessVxObsEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=BalancethicknessVxObsEnum() -%BALANCETHICKNESSVXOBSENUM - Enum of BalancethicknessVxObs -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=BalancethicknessVxObsEnum() - -macro=StringToEnum('BalancethicknessVxObs'); Index: sm/trunk-jpl/src/m/enum/BalancethicknessVyObsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessVyObsEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=BalancethicknessVyObsEnum() -%BALANCETHICKNESSVYOBSENUM - Enum of BalancethicknessVyObs -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=BalancethicknessVyObsEnum() - -macro=StringToEnum('BalancethicknessVyObs'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18475) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18476) @@ -301,13 +301,5 @@ def BalancethicknessApparentMassbalanceEnum(): return StringToEnum("BalancethicknessApparentMassbalance")[0] def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0] -def IrrotationalH2MisfitEnum(): return StringToEnum("IrrotationalH2Misfit")[0] -def IrrotationalDirectionMisfitEnum(): return StringToEnum("IrrotationalDirectionMisfit")[0] -def IrrotationalVelMisfitEnum(): return StringToEnum("IrrotationalVelMisfit")[0] -def IrrotationalAlongGradientNuEnum(): return StringToEnum("IrrotationalAlongGradientNu")[0] -def IrrotationalAcrossGradientNuEnum(): return StringToEnum("IrrotationalAcrossGradientNu")[0] -def BalancethicknessNuEnum(): return StringToEnum("BalancethicknessNu")[0] -def BalancethicknessVxObsEnum(): return StringToEnum("BalancethicknessVxObs")[0] -def BalancethicknessVyObsEnum(): return StringToEnum("BalancethicknessVyObs")[0] -def BalancethicknessThicknessObsEnum(): return StringToEnum("BalancethicknessThicknessObs")[0] +def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0] def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0] def SMBEnum(): return StringToEnum("SMB")[0] @@ -590,4 +582,5 @@ def OneLayerP4zEnum(): return StringToEnum("OneLayerP4z")[0] def CrouzeixRaviartEnum(): return StringToEnum("CrouzeixRaviart")[0] +def LACrouzeixRaviartEnum(): return StringToEnum("LACrouzeixRaviart")[0] def SaveResultsEnum(): return StringToEnum("SaveResults")[0] def BoolExternalResultEnum(): return StringToEnum("BoolExternalResult")[0] Index: sm/trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=IrrotationalAcrossGradientNuEnum() -%IRROTATIONALACROSSGRADIENTNUENUM - Enum of IrrotationalAcrossGradientNu -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=IrrotationalAcrossGradientNuEnum() - -macro=StringToEnum('IrrotationalAcrossGradientNu'); Index: sm/trunk-jpl/src/m/enum/IrrotationalAlongGradientNuEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/IrrotationalAlongGradientNuEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=IrrotationalAlongGradientNuEnum() -%IRROTATIONALALONGGRADIENTNUENUM - Enum of IrrotationalAlongGradientNu -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=IrrotationalAlongGradientNuEnum() - -macro=StringToEnum('IrrotationalAlongGradientNu'); Index: sm/trunk-jpl/src/m/enum/IrrotationalDirectionMisfitEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/IrrotationalDirectionMisfitEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=IrrotationalDirectionMisfitEnum() -%IRROTATIONALDIRECTIONMISFITENUM - Enum of IrrotationalDirectionMisfit -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=IrrotationalDirectionMisfitEnum() - -macro=StringToEnum('IrrotationalDirectionMisfit'); Index: sm/trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=IrrotationalH2MisfitEnum() -%IRROTATIONALH2MISFITENUM - Enum of IrrotationalH2Misfit -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=IrrotationalH2MisfitEnum() - -macro=StringToEnum('IrrotationalH2Misfit'); Index: sm/trunk-jpl/src/m/enum/IrrotationalVelMisfitEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/IrrotationalVelMisfitEnum.m (revision 18475) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=IrrotationalVelMisfitEnum() -%IRROTATIONALVELMISFITENUM - Enum of IrrotationalVelMisfit -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh -% Please read src/c/shared/Enum/README for more information -% -% Usage: -% macro=IrrotationalVelMisfitEnum() - -macro=StringToEnum('IrrotationalVelMisfit'); Index: /issm/trunk-jpl/src/m/enum/LACrouzeixRaviartEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/LACrouzeixRaviartEnum.m (revision 18476) +++ /issm/trunk-jpl/src/m/enum/LACrouzeixRaviartEnum.m (revision 18476) @@ -0,0 +1,11 @@ +function macro=LACrouzeixRaviartEnum() +%LACROUZEIXRAVIARTENUM - Enum of LACrouzeixRaviart +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=LACrouzeixRaviartEnum() + +macro=StringToEnum('LACrouzeixRaviart');