Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 4062) +++ /issm/trunk/src/c/Makefile.am (revision 4063) @@ -318,4 +318,5 @@ ./EnumDefinitions/StringAsEnum.cpp\ ./modules/ModelProcessorx/ModelProcessorx.h\ + ./modules/ModelProcessorx/ModelProcessorx.cpp\ ./modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp\ ./modules/ModelProcessorx/NodesPartitioning.cpp\ @@ -405,6 +406,6 @@ ./modules/MinVzx/MinVzx.h\ ./modules/MinVzx/MinVzx.cpp\ - ./modules/DuplicateInputx/DuplicateInputx.h\ - ./modules/DuplicateInputx/DuplicateInputx.cpp\ + ./modules/InputDuplicatex/InputDuplicatex.h\ + ./modules/InputDuplicatex/InputDuplicatex.cpp\ ./modules/ScaleInputx/ScaleInputx.h\ ./modules/ScaleInputx/ScaleInputx.cpp\ @@ -836,4 +837,5 @@ ./EnumDefinitions/StringAsEnum.cpp\ ./modules/ModelProcessorx/ModelProcessorx.h\ + ./modules/ModelProcessorx/ModelProcessorx.cpp\ ./modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp\ ./modules/ModelProcessorx/NodesPartitioning.cpp\ @@ -925,6 +927,6 @@ ./modules/MinVzx/MinVzx.h\ ./modules/MinVzx/MinVzx.cpp\ - ./modules/DuplicateInputx/DuplicateInputx.h\ - ./modules/DuplicateInputx/DuplicateInputx.cpp\ + ./modules/InputDuplicatex/InputDuplicatex.h\ + ./modules/InputDuplicatex/InputDuplicatex.cpp\ ./modules/ScaleInputx/ScaleInputx.h\ ./modules/ScaleInputx/ScaleInputx.cpp\ Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 4063) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 4063) @@ -0,0 +1,59 @@ +/*!\file ModelProcessorx + * \brief: create datasets using input binary file and a set of requested analyses + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../DataSet/DataSet.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../include/include.h" +#include "../../modules/modules.h" + +void ModelProcessorx(DataSet** pelements, DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads, Parameters** pparameters, ConstDataHandle IOMODEL,int solution_type,int nummodels, int* analysis_type_list){ + + int i; + int analysis_type; + + /*output: */ + DataSet* elements=NULL; + DataSet* nodes=NULL; + DataSet* vertices=NULL; + DataSet* materials=NULL; + DataSet* constraints=NULL; + DataSet* loads; + Parameters* parameters=NULL; + + + /*intermediary: */ + IoModel* iomodel=NULL; + + _printf_(" fill model with matlab workspace data\n"); + iomodel = new IoModel(IOMODEL); + + for(i=0;inummodels=nummodels; this->solution_type=in_solution_type; - analysis_counter=nummodels-1; //point to last analysis_type carried out. + this->analysis_counter=nummodels-1; //point to last analysis_type carried out. this->results=new DataSet(); //not initialized by CreateDataSets @@ -49,14 +49,12 @@ for(i=0;ielements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,IOMODEL,this->solution_type,nummodels,analyses); + /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ for(i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum); _printf_("%s%i%s%g\n"," value of misfit J after optimization #",n+1,": ",J[n]); Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4062) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4063) @@ -36,8 +36,8 @@ /*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/ if(qmu_analysis){ - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVxEnum,VxEnum); - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVyEnum,VyEnum); - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVzEnum,VzEnum); - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuPressureEnum,PressureEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVxEnum,VxEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVyEnum,VyEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVzEnum,VzEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuPressureEnum,PressureEnum); } @@ -72,5 +72,5 @@ //"recondition" pressure computed previously: - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); ScaleInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning); Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 4062) +++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 4063) @@ -54,5 +54,5 @@ /*Use ControlParameterEnum input to reinitialize our input parameter: */ - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type); /*Use search scalar to shoot parameter in the gradient direction: */ Index: /issm/trunk/src/c/solutions/steadystate_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 4062) +++ /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 4063) @@ -43,9 +43,9 @@ if(verbose)_printf_("%s\n","saving velocity, temperature and pressure to check for convergence at next step"); - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxOldEnum); - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyOldEnum); - if(dim==3)DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzOldEnum); - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureOldEnum); - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxOldEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyOldEnum); + if(dim==3)InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzOldEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureOldEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum); //increase counter Index: /issm/trunk/src/c/solutions/stokescontrolinit.cpp =================================================================== --- /issm/trunk/src/c/solutions/stokescontrolinit.cpp (revision 4062) +++ /issm/trunk/src/c/solutions/stokescontrolinit.cpp (revision 4063) @@ -44,5 +44,5 @@ //recondition" pressure computed previously: - DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); ScaleInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning); Index: sm/trunk/src/m/solutions/jpl/CreateFemModel.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/CreateFemModel.m (revision 4062) +++ (revision ) @@ -1,38 +1,0 @@ -%----------------------------------------------------------------------- -% CreateFEMModel.m: from a matlab ISSM @model class object, create a finite element -% model structure, with "c" code datasets (elements, nodes, loads, constraints, -% parameters, partitioning vectors. -% -%----------------------------------------------------------------------- - -function m=CreateFEMModel(md,solution_type,analysis_types) - - displaystring(md.verbose,'\n reading data from model %s...',md.name); - [m.elements,m.nodes,m.vertices,m.constraints,m.loads,m.materials,m.parameters]=ModelProcessor(md,solution_type,analysis_types); - - for i=1:length(analysis_types), - - displaystring(md.verbose,'%s',' generating degrees of freedom...'); - if ~isfield(m,'part') [m.vertices,m.part,m.tpart]=VerticesDof(vertices, m.parameters); - [m.nodes]=NodesDof(m.nodes,m.parameters); - - displaystring(md.verbose,'%s',' generating single point constraints...'); - [m.nodes,m.yg(i)]=SpcNodes(m.nodes,m.constraints,analysis_types[i]); - - displaystring(md.verbose,'%s',' generating rigid body constraints...'); - [m.Rmg(i),m.nodes]=MpcNodes(m.nodes,m.constraints,analysis_types[i]); - - displaystring(md.verbose,'%s',' generating node sets...'); - m.nodesets(i)=BuildNodeSets(m.nodes,analysis_types[i]); - - displaystring(md.verbose,'%s',' reducing single point constraints vector...'); - m.ys(i)=Reducevectorgtos(m.yg(i).vector,m.nodesets(i)); - - displaystring(md.verbose,'%s',' normalizing rigid body constraints matrix...'); - m.Gmn(i)= NormalizeConstraints(m.Rmg(i),m.nodesets(i)); - - displaystring(md.verbose,'%s',' configuring element and loads...'); - [m.elements,m.loads,m.nodes,m.parameters] = ConfigureObjects( m.elements, m.loads, m.nodes, m.vertices,m.materials,m.parameters); - end - -end Index: /issm/trunk/src/m/solutions/jpl/NewFemModel.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/NewFemModel.m (revision 4063) +++ /issm/trunk/src/m/solutions/jpl/NewFemModel.m (revision 4063) @@ -0,0 +1,50 @@ +function femmodel=NewFemModel(md,solution_type,analysis_types); +%NEWFEMMODEL - create a finite element model out of the matlab base \@model md. +% For each analysis_type contained in analysis_types, create a set of nodes, constraints +% and loads. All analyses rely on the same elements, vertices and parameters. See +% FemModel.cpp in src/c/objects for more information on the FemModel implementation in c++ +% +% Usage: +% femmodel=NewFemModel(md,solution_type,analysis_types) +% + + + femmodel.solution_type=solution_type; + femmodel.analysis_counter=nummodels; %point to last analysis_type carried out + + %Dynamically allocate whatever is a list of length nummodels: */ + femmodel.analysis_type_list=analysis_types; + + displaystring(md.verbose,'\n reading data from model %s...',md.name); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.constraints,femmodel.loads,femmodel.materials,femmodel.parameters]=ModelProcessor(md,solution_type,femmodel.analysis_type_list); + + %now, go through all analyses types and post-process datasets + for i=1:length(analysis_types), + + analysis_type=femmodel.analysis_type_list(i); + displaystring(md.verbose,'%s%s',' dealing with analysis type: ',EnumAsString(analysis_type); + + displaystring(md.verbose,'%s',' generating degrees of freedofemmodel...'); + if ~isfield(femmodel,'part') [femmodel.vertices,femmodel.part,femmodel.tpart]=VerticesDof(vertices, femmodel.parameters); %do not create partition vector twice! we only have one set of vertices! + + [femmodel.nodes]=NodesDof(femmodel.nodes,femmodel.parameters); + + displaystring(md.verbose,'%s',' generating single point constraints...'); + [femmodel.nodes,femmodel.yg(i)]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type); + + displaystring(md.verbose,'%s',' generating rigid body constraints...'); + [femmodel.Rmg(i),femmodel.nodes]=MpcNodes(femmodel.nodes,femmodel.constraints,analysis_types); + + displaystring(md.verbose,'%s',' generating node sets...'); + femmodel.nodesets(i)=BuildNodeSets(femmodel.nodes,analysis_type); + + displaystring(md.verbose,'%s',' reducing single point constraints vector...'); + femmodel.ys(i)=Reducevectorgtos(femmodel.yg(i).vector,femmodel.nodesets(i)); + + displaystring(md.verbose,'%s',' normalizing rigid body constraints matrix...'); + femmodel.Gmn(i)= NormalizeConstraints(femmodel.Rmg(i),femmodel.nodesets(i)); + + displaystring(md.verbose,'%s',' configuring element and loads...'); + [femmodel.elements,femmodel.loads,femmodel.nodes,femmodel.parameters] = ConfigureObjects( femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices,femmodel.materials,femmodel.parameters); + end + end Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m (revision 4062) +++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m (revision 4063) @@ -6,104 +6,71 @@ % -%recover parameters common to all solutions -verbose=model.dhu.parameters.Verbose; -dim=model.dh.parameters.Dim; -ishutter=model.dhu.parameters.IsHutter; -ismacayealpattyn=model.dh.parameters.IsMacAyealPattyn; -isstokes=model.ds.parameters.IsStokes; -numrifts=model.dhu.parameters.NumRifts; -qmu_analysis=model.dh.parameters.QmuAnalysis; + %some parameters + modify_loads=boolean(1); + conserve_loads=boolean(1); - -%for qmu analysis, be sure the velocity input we are starting from is the one in the parameters: -if qmu_analysis, - model=ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum); - model=ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum); - model=ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum); -end + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + ishutter=femmodel.parameters.IsHutter; + ismacayealpattyn=femmodel.parameters.IsMacAyealPattyn; + isstokes=femmodel.parameters.IsStokes; + stokesreconditioning=femmodel.parameters.stokesreconditioning; + qmu_analysis=femmodel.parameters.QmuAnalysis; -%Compute slopes: -[surfaceslopex,surfaceslopey]=slope_core(model.sl,SurfaceAnalysisEnum); -[bedslopex,bedslopey]=slope_core(model.sl,BedAnalysisEnum); - -%Update: -model=ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum); -model=ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum); -model=ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum); -model=ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum); - -if ishutter, - - displaystring(verbose,'\n%s',['computing hutter velocities...']); - u_g=diagnostic_core_linear(model.dhu,DiagnosticAnalysisEnum(),HutterAnalysisEnum()); - - displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']); - p_g=ComputePressure(model.dhu.elements,model.dhu.nodes,model.dhu.vertices,model.dhu.loads,model.dhu.materials,model.dhu.parameters,DiagnosticAnalysisEnum(),HutterAnalysisEnum()); - - displaystring(verbose,'\n%s',['update boundary conditions for macyeal pattyn using hutter results...']); - if ismacayealpattyn, - model.dh.y_g=u_g; - model.dh.ys=Reducevectorgtos(model.dh.y_g,model.dh.nodesets); + %for qmu analysis, be sure the velocity input we are starting from is the one in the parameters: + if qmu_analysis, + InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum); + InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum); + InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum); + InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum); end -end - -if ismacayealpattyn, + %Compute slopes: + if(ishutter)surfaceslope_core(femmodel);end + if(isstokes)bedslope_core(femmodel);end - displaystring(verbose,'\n%s',['computing horizontal velocities...']); - [u_g model.dh.loads]=diagnostic_core_nonlinear(model.dh,model.dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + if ishutter, - if dim==2, - displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']); - p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + displaystring(verbose,'\n%s',['computing hutter velocities...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticHutterAnalysisEnum); + femmodel=solver_linear(femmodel); + + if(ismacayealpattyn)femmodel=ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum); end + end -end + + if ismacayealpattyn, + + displaystring(verbose,'\n%s',['computing horizontal velocities...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum); + femmodel=solver_diagnostic_nonlinear(femmodel,modify_loads); + end -if dim==3, - displaystring(verbose,'\n%s',['extruding horizontal velocities...']); - u_g_horiz=FieldExtrude(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,u_g,'velocity',1); + if dim==3, + + displaystring(verbose,'\n%s',['computing vertical velocities...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticVertAnalysisEnum); + femmodel=solver_linear(femmodel); - [vx,vy]=SplitSolutionVector(u_g_horiz,model.dh.parameters.NumberOfNodes,model.dh.parameters.NumberOfDofsPerNode); - model=ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum); - model=ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum); - - displaystring(verbose,'\n%s',['computing vertical velocities...']); - u_g_vert=diagnostic_core_linear(model.dv,DiagnosticAnalysisEnum(),VertAnalysisEnum()); - model=ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum); + if isstokes, - displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']); - p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); - model=ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum); - u_g=zeros(model.dh.parameters.NumberOfNodes*3,1); %%%%%%%%%%%%%%%%%%%%%%% NEED TO BE CLEANED - u_g(1:3:end)=vx; - u_g(2:3:end)=vy; - u_g(3:3:end)=u_g_vert; - - if isstokes, + %"recondition" pressure computed previously: + InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); + ScaleInput(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning); - %"recondition" pressure - p_g=p_g/m_ds.parameters.stokesreconditioning; + displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); + femmodel=ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum); - displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); - model.ds.y_g=zeros(model.ds.nodesets.gsize,1); - model.ds.y_g(dofsetgen([1,2],4,model.ds.nodesets.gsize))=u_g; - model.ds.y_g(dofsetgen([3],4,model.ds.nodesets.gsize))=u_g_vert; - model.ds.ys=Reducevectorgtos(model.ds.y_g,model.ds.nodesets); + displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticStokesAnalysisEnum); + femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); + end + end - displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']); - u_g=diagnostic_core_nonlinear(model.ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); - - %"decondition" pressure - p_g=u_g(4:4:end)*model.dh.parameters.stokesreconditioning; - end -end -%load onto results -results.step=1; -results.time=0; -results.u_g=u_g; -results.p_g=p_g; - -if numrifts, - results.riftproperties=OutputRifts(model.dh.loads,model.dh.parameters); -end + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + if(dim==3) InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); Index: /issm/trunk/src/mex/InputDuplicate/InputDuplicate.cpp =================================================================== --- /issm/trunk/src/mex/InputDuplicate/InputDuplicate.cpp (revision 4063) +++ /issm/trunk/src/mex/InputDuplicate/InputDuplicate.cpp (revision 4063) @@ -0,0 +1,62 @@ +/*\file InputExtrude.c + *\brief: extrude input vertically + */ + +#include "./InputExtrude.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*diverse: */ + int noerr=1; + + /*input datasets: */ + DataSet* elements=NULL; + DataSet* nodes=NULL; + DataSet* vertices=NULL; + DataSet* loads=NULL; + DataSet* materials=NULL; + Parameters* parameters=NULL; + int NameEnum; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InputExtrudeUsage); + + /*Input datasets: */ + FetchData(&elements,ELEMENTS); + FetchData(&nodes,NODES); + FetchData(&vertices,VERTICES); + FetchData(&loads,LOADS); + FetchData(&materials,MATERIALS); + FetchParams(¶meters,PARAMETERS); + + FetchData(&reinitialized_enum,REINITIALIZEENUM); + FetchData(&original_enum,ORIGINALENUM); + + /*!Call core code: */ + InputDuplicatex( elements, nodes, vertices, loads, materials,parameters,reinitialized_enum, original_enum); + + /*write output : */ + WriteData(ELEMENTSOUT,elements); + + /*Free ressources: */ + delete elements; + delete nodes; + delete vertices; + delete loads; + delete materials; + delete parameters; + + /*end module: */ + MODULEEND(); + +} + +void InputExtrudeUsage(void) +{ + _printf_("\n"); + _printf_(" usage: [elements] = %s(elements, nodes, vertices, loads, materials, parameters, reinitialized_enum,original_enum);\n",__FUNCT__); + _printf_("\n"); +} Index: /issm/trunk/src/mex/InputDuplicate/InputDuplicate.h =================================================================== --- /issm/trunk/src/mex/InputDuplicate/InputDuplicate.h (revision 4063) +++ /issm/trunk/src/mex/InputDuplicate/InputDuplicate.h (revision 4063) @@ -0,0 +1,37 @@ +/* + InputExtrude.h +*/ + +#ifndef _INPUTEXTRUDE_H +#define _INPUTEXTRUDE_H + +/* local prototypes: */ +void InputExtrudeUsage(void); + +#include "../../c/modules/modules.h" +#include "../../c/DataSet/DataSet.h" +#include "../../c/shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "InputExtrude" + +/* serial input macros: */ +#define ELEMENTS (mxArray*)prhs[0] +#define NODES (mxArray*)prhs[1] +#define VERTICES (mxArray*)prhs[2] +#define LOADS (mxArray*)prhs[3] +#define MATERIALS (mxArray*)prhs[4] +#define PARAMETERS (mxArray*)prhs[5] +#define REINITIALIZEENUM (mxArray*)prhs[6] +#define ORIGINALENUM (mxArray*)prhs[7] + +/* serial output macros: */ +#define ELEMENTSOUT (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 8 + +#endif /* _INPUTEXTRUDE_H */ Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 4062) +++ /issm/trunk/src/mex/Makefile.am (revision 4063) @@ -30,4 +30,5 @@ Gradj\ HoleFiller \ + InputDuplicate\ InterpFromGridToMesh \ InterpFromMeshToMesh2d \ @@ -164,4 +165,7 @@ HoleFiller/HoleFiller.h +InputDuplicate_SOURCES = InputDuplicate/InputDuplicate.cpp\ + InputDuplicate/InputDuplicate.h + InterpFromGridToMesh_SOURCES = InterpFromGridToMesh/InterpFromGridToMesh.cpp\ InterpFromGridToMesh/InterpFromGridToMesh.h Index: /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp =================================================================== --- /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 4062) +++ /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 4063) @@ -11,5 +11,7 @@ /*inputs: */ - double* analyses=NULL; + int solution_type; + double* double_analyses=NULL; + int* analyses=NULL; int numanalyses; @@ -23,7 +25,4 @@ Parameters* parameters=NULL; - /*input model: */ - IoModel* iomodel=NULL; - /*Boot module: */ MODULEBOOT(); @@ -32,13 +31,12 @@ CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&ModelProcessorUsage); - /*Fill iomodel with matlab workspace data: */ - iomodel=new IoModel(MODEL); //this routine goes through the entire MODEL matlab class, and starts filling the corresponding "c" code model object. - FetchData(&analyses,&numanalyses,ANALYSES); + FetchData(&solution_type,SOLUTIONTYPE); + + FetchData(&double_analyses,&numanalyses,ANALYSES); + analyses=(int*)xmalloc(numanalyses*sizeof(int)); + for(i=0;i