Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4650) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4651) @@ -1034,7 +1034,9 @@ /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ if (analysis_type==DiagnosticHorizAnalysisEnum) - GetSolutionFromInputsDiagnosticHoriz(solution); + GetSolutionFromInputsDiagnosticHoriz(solution); else if (analysis_type==AdjointHorizAnalysisEnum) GetSolutionFromInputsAdjointHoriz(solution); + else if (analysis_type==DiagnosticHutterAnalysisEnum) + GetSolutionFromInputsDiagnosticHutter(solution); else ISSMERROR("analysis: %s not supported yet",EnumAsString(analysis_type)); @@ -5352,4 +5354,40 @@ } /*}}}*/ +/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/ +void Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){ + + int i; + + const int numvertices=3; + const int numdofpervertex=2; + const int numdof=numdofpervertex*numvertices; + double gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + + int doflist[numdof]; + double values[numdof]; + double vx; + double vy; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); + inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); + values[i*numdofpervertex+0]=vx; + values[i*numdofpervertex+1]=vy; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + +} +/*}}}*/ /*FUNCTION Tria::GradjDragStokes {{{1*/ void Tria::GradjDragStokes(Vec gradient){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 4650) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 4651) @@ -154,4 +154,5 @@ void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsAdjointHoriz(Vec solution); + void GetSolutionFromInputsDiagnosticHutter(Vec solution); void GradjDragStokes(Vec gradient); void InputUpdateFromSolutionAdjointHoriz( double* solution); Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4650) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4651) @@ -51,4 +51,8 @@ if(verbose)_printf_("%s\n"," computing hutter velocities..."); + + //Take the last velocity into account so that the velocity on the MacAyeal domain is not zero + if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHutterAnalysisEnum); + femmodel->SetCurrentConfiguration(DiagnosticHutterAnalysisEnum); solver_linear(NULL,femmodel); Index: /issm/trunk/src/c/solutions/stokescontrolinit.cpp =================================================================== --- /issm/trunk/src/c/solutions/stokescontrolinit.cpp (revision 4650) +++ /issm/trunk/src/c/solutions/stokescontrolinit.cpp (revision 4651) @@ -42,6 +42,5 @@ femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads); - - //vertical velocity + femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum); solver_linear(NULL,femmodel); Index: sm/trunk/src/m/solutions/ControlInitialization.m =================================================================== --- /issm/trunk/src/m/solutions/ControlInitialization.m (revision 4650) +++ (revision ) @@ -1,76 +1,0 @@ -function models=ControlInitialization(models); - -%recover models -m_dh=models.dh; -m_dv=models.dv; -m_ds=models.ds; -m_dhu=models.dhu; -m_sl=models.sl; - -%recover parameters common to all solutions -verbose=m_dhu.parameters.verbose; -dim=m_dhu.parameters.dim; -ishutter=m_dhu.parameters.ishutter; -ismacayealpattyn=m_dh.parameters.ismacayealpattyn; -isstokes=m_ds.parameters.isstokes; - -%If Pattyn or MacAyeal, just return, no initialization needed -if ~isstokes - models.active=m_dh; - return -end - -%1: Get slopes once for all - -%compute slopes -displaystring(verbose,'\n%s',['computing bed slope (x and y derivatives)...']); -slopex=diagnostic_core_linear(m_sl,SlopecomputeAnalysisEnum(),BedXAnalysisEnum()); -slopey=diagnostic_core_linear(m_sl,SlopecomputeAnalysisEnum(),BedYAnalysisEnum()); -slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.vertices,m_sl.loads,m_sl.materials,m_sl.parameters,slopex,'slopex',0); -slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.vertices,m_sl.loads,m_sl.materials,m_sl.parameters,slopey,'slopey',0); - -%Add to inputs -inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); -inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); - -%2: update spcs using MacAyeal or Pattyn - -%horizontal velocities -displaystring(verbose,'\n%s',['computing horizontal velocities...']); -u_g=diagnostic_core_nonlinear(m_dh,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); -displaystring(verbose,'\n%s',['extruding horizontal velocities...']); -u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,u_g,'velocity',1); - -%vertical velocities -displaystring(verbose,'\n%s',['computing vertical velocities...']); -inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes); -u_g_vert=diagnostic_core_linear(m_dv,DiagnosticAnalysisEnum(),VertAnalysisEnum()); - -%create 3d u_g -displaystring(verbose,'\n%s',['combining horizontal and vertical velocities...']); -u_g=zeros(m_dh.parameters.numberofnodes*3,1); -u_g(dofsetgen([1,2],3,m_dh.parameters.numberofnodes*3))=u_g_horiz; -u_g(dofsetgen([3],3,m_dh.parameters.numberofnodes*3))=u_g_vert; - -% get pressure (reconditionned) and create 4d u_g -displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']); -p_g=ComputePressure(m_dh.elements,m_dh.nodes,mdh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); -p_g=p_g/m_ds.parameters.stokesreconditioning; -u_g_stokes=zeros(m_ds.nodesets.gsize,1); -u_g_stokes(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g; -u_g_stokes(dofsetgen([4],4,m_ds.nodesets.gsize))=p_g; -inputs=add(inputs,'velocity',u_g_stokes,'doublevec',4,m_ds.parameters.numberofnodes); - -%update BC (spcs) -displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); -m_ds.y_g=zeros(m_ds.nodesets.gsize,1); -m_ds.y_g(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g; -m_ds.ys=Reducevectorgtos(m_ds.y_g,m_ds.nodesets); - -%Compute Stokes velocities once to have a reasonably good ug in input -displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']); -u_g=diagnostic_core_nonlinear(m_ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); -inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes); - -%assign output -models.active=m_ds; Index: /issm/trunk/src/m/solutions/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/diagnostic_core.m (revision 4650) +++ /issm/trunk/src/m/solutions/diagnostic_core.m (revision 4651) @@ -35,4 +35,8 @@ displaystring(verbose,'\n%s',['computing hutter velocities...']); + + %Take the last velocity into account so that the velocity on the MacAyeal domain is not zero + if(ismacayealpattyn)femmodel=ResetBoundaryConditions(femmodel,DiagnosticHutterAnalysisEnum); end + femmodel=SetCurrentConfiguration(femmodel,DiagnosticHutterAnalysisEnum); femmodel=solver_linear(femmodel); Index: sm/trunk/src/m/solutions/gradjcompute_core.m =================================================================== --- /issm/trunk/src/m/solutions/gradjcompute_core.m (revision 4650) +++ (revision ) @@ -1,52 +1,0 @@ -function results=gradjcompute_core(models); - -%recover active models -m=models.active; - -%recover parameters -verbose=m.parameters.Verbose; -dim=m.parameters.Dim; -extrude_param=m.parameters.ExtrudeParam; -ishutter=m.parameters.IsHutter; -ismacayealpattyn=m.parameters.IsMacAyealPattyn; -isstokes=m.parameters.IsStokes; -analysis_type=m.parameters.AnalysisType; -sub_analysis_type=m.parameters.SubAnalysisType; - -%Recover solution for this stiffness and right hand side: -displaystring(verbose,'%s',' computing velocities...'); -[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,analysis_type,sub_analysis_type); -inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes); - -%Buid Du, difference between observed velocity and model velocity. -displaystring(verbose,'%s',' computing Du...'); -[Du_g]=Du(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - -%Reduce adjoint load from g-set to f-set -[Du_f] = Reduceloadfromgtof( Du_g, m.Gmn, K_fs0, m.ys, m.nodesets,1);% 1 because we want ys0 - -%Solve for adjoint vector: -displaystring(verbose,'%s',' computing adjoint state...'); -lambda_f=Solver(K_ff0,Du_f,[],m.parameters); - -%Merge back to g set -lambda_g= Mergesolutionfromftog( lambda_f, m.Gmn, m.ys, m.nodesets,1); % 1 because we want ys0 -inputs=add(inputs,'adjoint',lambda_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes); - -%Compute gradJ -grad_g=Gradj(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); -if (dim==3 & extrude_param), - displaystring(verbose,'%s',' extruding gradient...'); - grad_g=FieldExtrude(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,grad_g,'gradj',0); -end - -if m.parameters.ControlSteady, - %compute initial velocity from diagnostic_core (horiz+vertical) - displaystring(verbose,'%s',' compute 3d initial velocity...'); - results_diag=diagnostic_core(models); - u_g=results_diag.u_g; -end - -%Assign output -results.u_g=u_g; -results.grad_g=grad_g; Index: /issm/trunk/src/m/solutions/stokescontrolinit.m =================================================================== --- /issm/trunk/src/m/solutions/stokescontrolinit.m (revision 4650) +++ /issm/trunk/src/m/solutions/stokescontrolinit.m (revision 4651) @@ -21,6 +21,5 @@ femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum); femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); - - %vertical velocity + femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum); femmodel=solver_linear(femmodel);