Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 22177) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 22178) @@ -192,5 +192,6 @@ /*Intermediaries*/ - IssmDouble dlambda[2],ds[2],D0,omega,Jdet; + int n=3; + IssmDouble dlambda[2],ds[2],slopex,slopey,slope,omega,Jdet; IssmDouble *xyz_list= NULL; @@ -207,7 +208,9 @@ element->GradientIndexing(&vertexpidlist[0],control_index); Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input); - Input* s_input = element->GetInput(SurfaceEnum); _assert_(s_input); - Input* D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input); Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); + Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); + IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); Gauss* gauss=element->NewGauss(2); @@ -218,12 +221,15 @@ element->NodalFunctionsP1(basis,gauss); - D0_input->GetInputValue(&D0,gauss); omega_input->GetInputValue(&omega,gauss); adjoint_input->GetInputDerivativeValue(&dlambda[0],xyz_list,gauss); - s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + surfaceslopex_input->GetInputValue(&slopex,gauss); + surfaceslopey_input->GetInputValue(&slopey,gauss); + + slope = sqrt(slopex*slopex + slopey*slopey); /*Build gradient vector (actually -dJ/da): */ for(int i=0;iweight*basis[i]*exp(omega)*D0*(ds[0]*dlambda[0] + ds[1]*dlambda[1]); + ge[i]+= -Jdet*gauss->weight*basis[i]*pow(rhog,n)*exp(omega)*pow(slope,n-1)*(ds[0]*dlambda[0] + ds[1]*dlambda[1]); _assert_(!xIsNan(ge[i])); } Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 22177) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 22178) @@ -39,4 +39,6 @@ iomodel->FetchDataToInput(elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum); iomodel->FetchDataToInput(elements,"md.balancethickness.omega",BalancethicknessOmegaEnum); + iomodel->FetchDataToInput(elements,"md.balancethickness.slopex",SurfaceSlopeXEnum); + iomodel->FetchDataToInput(elements,"md.balancethickness.slopey",SurfaceSlopeYEnum); /*Update elements: */ @@ -63,68 +65,4 @@ return NULL; }/*}}}*/ -void Balancethickness2Analysis::CreateD0(Element* element){/*{{{*/ - - /*Intermediaries */ - IssmDouble Gamma,h,mu0,ds[2],Cmu,B,k,s,b,normds; - 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* D0 = xNew(numnodes); - - /*retrieve what we need: */ - element->GetVerticesCoordinates(&xyz_list); - Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); - Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); - 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); - if(surfaceslopex_input && surfaceslopey_input){ - surfaceslopex_input->GetInputValue(&ds[0],gauss); - surfaceslopey_input->GetInputValue(&ds[1],gauss); - } - else{ - surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); - } - - /*check slopes*/ - normds = sqrt(ds[0]*ds[0]+ds[1]*ds[1]); - if (normds==0.){ - _error_("surface slope is zero"); - } - if(normds<1.e-5){ - ds[0] = ds[0]/normds*1.e+5; - ds[1] = ds[1]/normds*1.e+5; - normds = 1.e-5; - } - - mu0 = pow(2.,(1-3*n)/(2.*n)) * B; - Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n); - - D0[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)/(n+2); - } - - /*Add input*/ - element->AddInput(BalancethicknessD0Enum,D0,element->GetElementType()); - //if(surfaceslopex_input && surfaceslopey_input){ - // element->DeleteInput(SurfaceSlopeXEnum); - // element->DeleteInput(SurfaceSlopeYEnum); - //} - - /*Clean up and return*/ - xDelete(D0); - xDelete(xyz_list); - delete gauss; -}/*}}}*/ ElementMatrix* Balancethickness2Analysis::CreateJacobianMatrix(Element* element){/*{{{*/ _error_("Not implemented"); @@ -133,5 +71,6 @@ /*Intermediaries */ - IssmDouble Jdet,D0,omega; + int n=3; + IssmDouble Jdet,omega,ds[2],slope; IssmDouble* xyz_list = NULL; @@ -145,10 +84,8 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input); - Input* D0_input = element->GetInput(BalancethicknessD0Enum); - if(!D0_input){ - this->CreateD0(element); - D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input); - } + Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input); + Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); + Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); + IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); /* Start looping on the number of gaussian points: */ @@ -158,10 +95,13 @@ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); element->JacobianDeterminant(&Jdet,xyz_list,gauss); - D0_input->GetInputValue(&D0,gauss); omega_input->GetInputValue(&omega,gauss); + surfaceslopex_input->GetInputValue(&ds[0],gauss); + surfaceslopey_input->GetInputValue(&ds[1],gauss); + + slope = sqrt(ds[0]*ds[0] + ds[1]*ds[1]); for(int i=0;ivalues[i*numnodes+j] += D0*exp(omega)*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); + Ke->values[i*numnodes+j] += pow(rhog,n)*exp(omega)*pow(slope,n-1)*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); } } Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h (revision 22177) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h (revision 22178) @@ -23,5 +23,4 @@ void Core(FemModel* femmodel); ElementVector* CreateDVector(Element* element); - void CreateD0(Element* element); ElementMatrix* CreateJacobianMatrix(Element* element); ElementMatrix* CreateKMatrix(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22177) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22178) @@ -1759,4 +1759,12 @@ for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case BalancethicknessOmegaEnum: + if(iomodel->Data("md.balancethickness.omega")){ + for(j=0;jData("md.balancethickness.omega")[tria_vertex_ids[j]-1]; + for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; + for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(BalancethicknessOmegaEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; Index: /issm/trunk-jpl/src/c/cores/adjointbalancethickness2_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/adjointbalancethickness2_core.cpp (revision 22177) +++ /issm/trunk-jpl/src/c/cores/adjointbalancethickness2_core.cpp (revision 22178) @@ -24,5 +24,5 @@ /*Call SurfaceAreax, because some it might be needed by PVector*/ - SurfaceAreax(NULL,femmodel); + //SurfaceAreax(NULL,femmodel); /*compute adjoint*/ Index: /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 22177) +++ /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 22178) @@ -13,11 +13,11 @@ /*parameters: */ - bool save_results; - IssmDouble l = 3.; + bool save_results; + //IssmDouble l = 3.; /*recover parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n"); + //if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n"); //femmodel->parameters->SetParam(l,SmoothThicknessMultiplierEnum); //femmodel->SetCurrentConfiguration(SmoothAnalysisEnum); @@ -26,5 +26,5 @@ //femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToSmoothEnum); //solutionsequence_linear(femmodel); - surfaceslope_core(femmodel); + //surfaceslope_core(femmodel); if(VerboseSolution()) _printf0_("call computational core:\n"); @@ -35,8 +35,6 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - const int numoutputs = 6; - int outputs[numoutputs] = {SurfaceEnum,SurfaceSlopeXEnum,SurfaceSlopeYEnum,VxEnum,VyEnum,VelEnum}; - //const int numoutputs = 4; - //int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum}; + const int numoutputs = 1; + int outputs[numoutputs] = {SurfaceEnum}; femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs); } Index: /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 22177) +++ /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 22178) @@ -104,3 +104,4 @@ xDelete(X); xDelete(X0); + xDelete(scaling_factors); }