Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 3779) +++ /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 3780) @@ -161,5 +161,6 @@ if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); - xfree((void**)&dofset);dofset=dofsetgen(3,dof012,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug,dofset,3*numberofnodes); + xfree((void**)&dofset);dofset=dofsetgen(2,dof01,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug,dofset,2*numberofnodes); + xfree((void**)&dofset);dofset=dofsetgen(1,dof2,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug_vert,dofset,1*numberofnodes); VecFree(&fem_ds->ys); VecFree(&fem_ds->ys0); Reducevectorgtosx(&fem_ds->ys,&fem_ds->ys0, fem_ds->yg->vector,fem_ds->nodesets); Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m (revision 3779) +++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m (revision 3780) @@ -3,13 +3,13 @@ % % Usage: -% results=diagnostic_core(models); +% results=diagnostic_core(model); % -%recover models -m_dh=models.dh; -m_dv=models.dv; -m_ds=models.ds; -m_dhu=models.dhu; -m_sl=models.sl; +%recover fem models +m_dh=model.dh; +m_dv=model.dv; +m_ds=model.ds; +m_dhu=model.dhu; +m_sl=model.sl; %recover parameters common to all solutions @@ -25,7 +25,7 @@ %for qmu analysis, be sure the velocity input we are starting from is the one in the parameters: if qmu_analysis, - ModelUpdateInputsFromVector(models,m_dh.vx,VxEnum,VertexEnum); - ModelUpdateInputsFromVector(models,m_dh.vy,VyEnum,VertexEnum); - ModelUpdateInputsFromVector(models,m_dh.vz,VzEnum,VertexEnum); + ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum); + ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum); + ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum); end @@ -35,8 +35,8 @@ %Update: -ModelUpdateInputsFromVector(models,surfaceslopex,SurfaceSlopexEnum,VertexEnum); -ModelUpdateInputsFromVector(models,surfaceslopey,SurfaceSlopeyEnum,VertexEnum); -ModelUpdateInputsFromVector(models,bedslopex,BedSlopexEnum,VertexEnum); -ModelUpdateInputsFromVector(models,bedslopey,BedSlopeyEnum,VertexEnum); +ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum); +ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum); +ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum); +ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum); if ishutter, @@ -73,20 +73,14 @@ [vx,vy]=SplitSolutionVector(ug,numberofnodes,numberofdofspernode_dh); - model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum); + ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum); - - 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,inputs,DiagnosticAnalysisEnum(),VertAnalysisEnum()); - - 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; + u_g_vert=diagnostic_core_linear(m_dv,DiagnosticAnalysisEnum(),VertAnalysisEnum()); + ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum); displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']); - p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum); if isstokes, @@ -95,27 +89,12 @@ p_g=p_g/m_ds.parameters.stokesreconditioning; - displaystring(verbose,'\n%s',['computing bed slope (x and y derivatives)...']); - slopex=diagnostic_core_linear(m_sl,inputs,SlopecomputeAnalysisEnum(),BedXAnalysisEnum()); - slopey=diagnostic_core_linear(m_sl,inputs,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); - - 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); - - %recombine u_g and p_g: - 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); - 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.y_g(dofsetgen([1,2],4,m_ds.nodesets.gsize))=u_g; + m_ds.y_g(dofsetgen([3],4,m_ds.nodesets.gsize))=u_g_vert; [m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets); displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']); - u_g=diagnostic_core_nonlinear(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); + u_g=diagnostic_core_nonlinear(m_ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); %"decondition" pressure Index: /issm/trunk/src/m/solutions/jpl/slope_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/slope_core.m (revision 3780) +++ /issm/trunk/src/m/solutions/jpl/slope_core.m (revision 3780) @@ -0,0 +1,58 @@ +function [slopex, slopey]=slope_core(fem,AnalysisEnum) +%INPUT [slopex, slopey]=slope_core(fem,AnalysisEnum) +% +% Usage: [slopex, slopey]=slope_core(fem,AnalysisEnum) +% +% Ex: +% [bedslopex, bedslopey]=slope_core(fem,BedEnum); +% [surfaceslopex, surfaceslopey]=slope_core(fem,SurfaceEnum); +% + + +verbose=m_dhu.parameters.verbose; +dim=m_dh.parameters.dim; +isstokes=m_ds.parameters.isstokes; +ishutter=m_dhu.parameters.ishutter; + +displaystring(verbose,'\n%s',['computing slope (x and y derivatives...)']); + +%Specify type of computations: +if AnalysisEnum==SurfaceAnalysisEnum, + xanalysis=SurfaceXAnalysisEnum; + yanalysis=SurfaceYAnalysisEnum; +elseif AnalysisEnum==BedAnalysisEnum, + xanalysis=BedXAnalysisEnum; + yanalysis=BedYAnalysisEnum; +else + displaystring(verbose,'%s%s%s','analysis ',EnumAsString(AnalysisEnum),' not supported yet!'); +end + + +%Early return possible? +if ~ishutter && AnalysisEnum==SurfaceAnalysisEnum, + %no need to compute Surface Slope except for Hutter: + slopex=[]; + slopey=[]; + return; +end + +if ~isstokes && AnalysisEnum==BedAnalysisEnum, + %no need to compute Bed Slope except for full Stokes: + + slopex=[]; + slopey=[]; + return; +end + +%Call on core computations: +slopex=diagnostic_core_linear(fem,SlopecomputeAnalysisEnum,xanalysis); +slopey=diagnostic_core_linear(fem,SlopecomputeAnalysisEnum,yanalysis); + + +%extrude if we are in 3D: +if dim==3, + displaystring(verbose,'%s\n','extruding slopes in 3d...'); + slopex=FieldExtrude(fem.elements,fem.nodes,fem.vertices,fem.loads,fem.materials,fem.parameters,slopex,'slopex',1); + slopey=FieldExtrude(fem.elements,fem.nodes,fem.vertices,fem.loads,fem.materials,fem.parameters,slopex,'slopey',1); +end +