Changeset 4055


Ignore:
Timestamp:
06/18/10 19:24:59 (15 years ago)
Author:
Eric.Larour
Message:

Keep simplifying solutions.
New convergence module at the input level, instead at the solution level, when solution vectors are not
available anymore.

Location:
issm/trunk/src/c
Files:
13 added
5 deleted
85 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/DataSet/DataSet.cpp

    r4053 r4055  
    684684        xfree((void**)&truedofs);
    685685        xfree((void**)&alltruedofs);
    686 
    687 }
    688 /*}}}*/
    689 /*FUNCTION DataSet::FieldExtrude{{{1*/
    690 void  DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){
    691 
    692         vector<Object*>::iterator object;
    693         Penta* penta=NULL;
    694         Node*  node=NULL;
    695 
    696         for ( object=objects.begin() ; object < objects.end(); object++ ){
    697 
    698                 if((*object)->Enum()==PentaEnum){
    699 
    700                         penta=(Penta*)(*object);
    701                         penta->FieldExtrude(field,field_serial,field_name,collapse);
    702 
    703                 }
    704                 if((*object)->Enum()==NodeEnum){
    705 
    706                         node=(Node*)(*object);
    707                         node->FieldExtrude(field,field_serial,field_name);
    708 
    709                 }
    710         }
    711686
    712687}

  • issm/trunk/src/c/DataSet/DataSet.h

    r4043 r4055  
    8080                void  InputExtrude(int enum_type);
    8181                int   DeleteObject(Object* object);
    82                 void  FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse);
    8382                void  OutputRifts(Vec riftproperties);
    8483                /*}}}*/

  • issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp

    r4051 r4055  
    2626                case VerticesEnum : return "Vertices";
    2727                case SolutionTypeEnum : return "SolutionType";
    28                 case DiagnosticSolutionEnum : return "DiagnosticSolution";
    29                 case SurfaceSlopeSolutionEnum : return "SurfaceSlopeSolution";
    30                 case BedSlopeSolutionEnum : return "BedSlopeSolution";
    3128                case AnalysisTypeEnum : return "AnalysisType";
    3229                case SubAnalysisTypeEnum : return "SubAnalysisType";
    … …  
    4643                case InverseAnalysisEnum : return "InverseAnalysis";
    4744                case ThermalAnalysisEnum : return "ThermalAnalysis";
    48                 case TransientAnalysisEnum : return "TransientAnalysis";
     45                case Transient2DAnalysisEnum : return "Transient2DAnalysis";
     46                case Transient3DAnalysisEnum : return "Transient3DAnalysis";
    4947                case SteadyAnalysisEnum : return "SteadyAnalysis";
    50                 case SlopeComputeAnalysisEnum : return "SlopeComputeAnalysis";
     48                case SlopeAnalysisEnum : return "SlopeAnalysis";
    5149                case BedSlopeXAnalysisEnum : return "BedSlopeXAnalysis";
    5250                case BedSlopeYAnalysisEnum : return "BedSlopeYAnalysis";
    … …  
    220218                case StringExternalResultEnum : return "StringExternalResult";
    221219                case JEnum : return "J";
     220                case RelativeEnum : return "Relative";
     221                case ResidualEnum : return "Residual";
     222                case AbsoluteEnum : return "Absolute";
    222223                case BetaEnum : return "Beta";
    223224                case CmGradientEnum : return "CmGradient";
    … …  
    225226                case CmMaxEnum : return "CmMax";
    226227                case CmMinEnum : return "CmMin";
     228                case AdjointEnum : return "Adjoint";
    227229                case GradientEnum : return "Gradient";
     230                case OldGradientEnum : return "OldGradient";
    228231                case ConnectivityEnum : return "Connectivity";
    229232                case ControlParameterEnum : return "ControlParameter";

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r4051 r4055  
    2222        /*Solution types {{{1 */
    2323        SolutionTypeEnum,
    24         DiagnosticSolutionEnum,
    25         SurfaceSlopeSolutionEnum,
    26         BedSlopeSolutionEnum,
    2724        /*}}}*/
    2825        /*Analysis types {{{1 */
    … …  
    4845        ThermalAnalysisEnum,
    4946        //transient
    50         TransientAnalysisEnum,
     47        Transient2DAnalysisEnum,
     48        Transient3DAnalysisEnum,
    5149        SteadyAnalysisEnum,
    5250        //slope
    53         SlopeComputeAnalysisEnum,
     51        SlopeAnalysisEnum,
    5452        BedSlopeXAnalysisEnum,
    5553        BedSlopeYAnalysisEnum,
    … …  
    254252        JEnum,
    255253        /*}}}*/
     254        /*Convergence{{{1*/
     255        RelativeEnum,
     256        ResidualEnum,
     257        AbsoluteEnum,
     258        /*}}}*/
    256259        /*Parameters{{{1*/
    257260        BetaEnum,
    … …  
    260263        CmMaxEnum,
    261264        CmMinEnum,
     265        AdjointEnum,
    262266        GradientEnum,
     267        OldGradientEnum,
    263268        ConnectivityEnum,
    264269        ControlParameterEnum,

  • issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp

    r4051 r4055  
    2424        else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    2525        else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
    26         else if (strcmp(name,"DiagnosticSolution")==0) return DiagnosticSolutionEnum;
    27         else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
    28         else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    2926        else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
    3027        else if (strcmp(name,"SubAnalysisType")==0) return SubAnalysisTypeEnum;
    … …  
    4441        else if (strcmp(name,"InverseAnalysis")==0) return InverseAnalysisEnum;
    4542        else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    46         else if (strcmp(name,"TransientAnalysis")==0) return TransientAnalysisEnum;
     43        else if (strcmp(name,"Transient2DAnalysis")==0) return Transient2DAnalysisEnum;
     44        else if (strcmp(name,"Transient3DAnalysis")==0) return Transient3DAnalysisEnum;
    4745        else if (strcmp(name,"SteadyAnalysis")==0) return SteadyAnalysisEnum;
    48         else if (strcmp(name,"SlopeComputeAnalysis")==0) return SlopeComputeAnalysisEnum;
     46        else if (strcmp(name,"SlopeAnalysis")==0) return SlopeAnalysisEnum;
    4947        else if (strcmp(name,"BedSlopeXAnalysis")==0) return BedSlopeXAnalysisEnum;
    5048        else if (strcmp(name,"BedSlopeYAnalysis")==0) return BedSlopeYAnalysisEnum;
    … …  
    218216        else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
    219217        else if (strcmp(name,"J")==0) return JEnum;
     218        else if (strcmp(name,"Relative")==0) return RelativeEnum;
     219        else if (strcmp(name,"Residual")==0) return ResidualEnum;
     220        else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
    220221        else if (strcmp(name,"Beta")==0) return BetaEnum;
    221222        else if (strcmp(name,"CmGradient")==0) return CmGradientEnum;
    … …  
    223224        else if (strcmp(name,"CmMax")==0) return CmMaxEnum;
    224225        else if (strcmp(name,"CmMin")==0) return CmMinEnum;
     226        else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
    225227        else if (strcmp(name,"Gradient")==0) return GradientEnum;
     228        else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
    226229        else if (strcmp(name,"Connectivity")==0) return ConnectivityEnum;
    227230        else if (strcmp(name,"ControlParameter")==0) return ControlParameterEnum;

  • issm/trunk/src/c/Makefile.am

    r4053 r4055  
    510510                                        ./modules/DepthAverageInputx/DepthAverageInputx.cpp\
    511511                                        ./modules/DepthAverageInputx/DepthAverageInputx.h\
    512                                         ./modules/FieldExtrudex/FieldExtrudex.cpp\
    513                                         ./modules/FieldExtrudex/FieldExtrudex.h\
     512                                        ./modules/VecExtrudex/VecExtrudex.cpp\
     513                                        ./modules/VecExtrudex/VecExtrudex.h\
    514514                                        ./modules/InputToResultx/InputToResultx.cpp\
    515515                                        ./modules/InputToResultx/InputToResultx.h\
     516                                        ./modules/InputConvergencex/InputConvergencex.cpp\
     517                                        ./modules/InputConvergencex/InputConvergencex.h\
    516518                                        ./modules/ExtrudeInputx/ExtrudeInputx.cpp\
    517519                                        ./modules/ExtrudeInputx/ExtrudeInputx.h\
    … …  
    10141016                                        ./modules/DepthAverageInputx/DepthAverageInputx.cpp\
    10151017                                        ./modules/DepthAverageInputx/DepthAverageInputx.h\
    1016                                         ./modules/FieldExtrudex/FieldExtrudex.cpp\
    1017                                         ./modules/FieldExtrudex/FieldExtrudex.h\
     1018                                        ./modules/VecExtrudex/VecExtrudex.cpp\
     1019                                        ./modules/VecExtrudex/VecExtrudex.h\
    10181020                                        ./modules/InputToResultx/InputToResultx.cpp\
    10191021                                        ./modules/InputToResultx/InputToResultx.h\
     1022                                        ./modules/InputConvergencex/InputConvergencex.cpp\
     1023                                        ./modules/InputConvergencex/InputConvergencex.h\
    10201024                                        ./modules/ExtrudeInputx/ExtrudeInputx.cpp\
    10211025                                        ./modules/ExtrudeInputx/ExtrudeInputx.h\
    … …  
    10301034                                        ./solutions/controlrestart.cpp\
    10311035                                        ./solutions/objectivefunctionC.cpp\
    1032                                         ./solutions/gradjcompute_core.cpp\
     1036                                        ./solutions/gradient_core.cpp\
     1037                                        ./solutions/adjoint_core.cpp\
    10331038                                        ./solutions/prognostic_core.cpp\
    10341039                                        ./solutions/prognostic2_core.cpp\
    … …  
    10401045                                        ./solutions/bedslope.cpp\
    10411046                                        ./solutions/bedslope_core.cpp\
    1042                                         ./solutions/transient_core.cpp\
    1043                                         ./solutions/transient_core_2d.cpp\
    1044                                         ./solutions/transient_core_3d.cpp\
     1047                                        ./solutions/transient2d_core.cpp\
     1048                                        ./solutions/transient3d_core.cpp\
    10451049                                        ./solutions/steadystate_core.cpp\
    10461050                                        ./solutions/ResetBoundaryConditions.cpp\
    … …  
    10581062bin_PROGRAMS =
    10591063else
    1060 bin_PROGRAMS = DiagnosticSolution.exe ThermalSolution.exe PrognosticSolution.exe Prognostic2Solution.exe BalancedthicknessSolution.exe Balancedthickness2Solution.exe BalancedvelocitiesSolution.exe TransientSolution.exe SteadystateSolution.exe SurfaceSlopeSolution.exe BedSlopeSolution.exe
     1064bin_PROGRAMS = DiagnosticSolution.exe ThermalSolution.exe PrognosticSolution.exe Prognostic2Solution.exe BalancedthicknessSolution.exe Balancedthickness2Solution.exe BalancedvelocitiesSolution.exe Transient2DSolution.exe Transient23Solution.exe SteadystateSolution.exe SurfaceSlopeSolution.exe BedSlopeSolution.exe
    10611065endif
    10621066
    … …  
    10931097BedSlopeSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_
    10941098
    1095 TransientSolution_exe_SOURCES = solutions/transient.cpp
    1096 TransientSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     1099Transient2DSolution_exe_SOURCES = solutions/transient2d.cpp
     1100Transient2DSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     1101
     1102Transient3DSolution_exe_SOURCES = solutions/transient3d.cpp
     1103Transient3DSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_

  • issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h

    r3913 r4055  
    1010
    1111/* local prototypes: */
    12 void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,
    13                         int analysis_type,int sub_analysis_type);
     12void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters);
    1413
    1514#endif  /* _MISFITX_H */

  • issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.cpp

    r3969 r4055  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type){
     12void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type,int average_enum_type){
    1313       
    1414        /*Intermediary*/
    … …  
    2222        for (i=0;i<elements->Size();i++){
    2323                element=(Element*)elements->GetObjectByOffset(i);
    24                 element->DepthAverageInputAtBase(enum_type);
     24                element->DepthAverageInputAtBase(enum_type,average_enum_type);
    2525        }
    2626
    2727        /*Then extrude vertically the new inputs*/
    28         elements->InputExtrude(enum_type);
    29 
     28        elements->InputExtrude(average_enum_type);
    3029}

  • issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.h

    r3913 r4055  
    99
    1010/* local prototypes: */
    11 void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type);
     11void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type,int average_enum_type);
    1212
    1313#endif  /* _DEPTHAVERAGEINPUTX_H */

  • issm/trunk/src/c/modules/Dux/Dux.h

    r3913 r4055  
    1010
    1111/* local prototypes: */
    12 void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,
    13                         int analysis_type,int sub_analysis_type);
    14 
     12void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters);
    1513#endif  /* _DUX_H */
    1614

  • issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp

    r3969 r4055  
    4343
    4444        /*Extrude field to all nodes: */
    45         nodes->FieldExtrude(field,field_serial,fieldname,0);
     45        for(i=0;i<nodes->Size();i++){
     46                node=(Node*)nodes->GetObjectByOffset(i);
     47                node->VecExtrude(field,field_serial);
     48        }
    4649
    4750        /*Assemble vector: */

  • issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h

    r4013 r4055  
    1010
    1111/* local prototypes: */
    12 void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters, int analysis_type, int sub_analysis_type);
     12void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters);
    1313
    1414#endif  /* _GETSOLUTIONFROMINPUTSXX_H */

  • issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp

    r4048 r4055  
    3838
    3939}
     40
     41void GetVectorFromInputsx( double* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, int NameEnum, int TypeEnum){
     42       
     43        /*output: */
     44        double* vector=NULL;
     45       
     46        /*intermediary: */
     47        Vec vec_vector=NULL;
     48
     49        GetVectorFromInputsx( &vec_vector, elements,nodes, vertices, loads, materials, parameters, NameEnum, TypeEnum);
     50        VecToMPISerial(&vector,vec_vector);
     51
     52        /*Free ressources:*/
     53        VecFree(&vec_vector);
     54
     55        /*Assign output pointers:*/
     56        *pvector=vector;
     57}

  • issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h

    r4048 r4055  
    1010/* local prototypes: */
    1111void    GetVectorFromInputsx( Vec* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters,int NameEnum,int TypeEnum);
     12void    GetVectorFromInputsx( double* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters,int NameEnum,int TypeEnum);
    1213
    1314#endif  /* _GETVECTORFROMINPUTSXX_H */

  • issm/trunk/src/c/modules/Gradjx/Gradjx.cpp

    r4053 r4055  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){
     12void Gradjx( Vec* pgradient, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){
    1313
    1414        int i;
    1515        int dim;
    16 
    17 
     16        int numberofvertices;
     17        Vec gradient=NULL;
     18       
    1819        /*First, get elements and loads configured: */
    1920        elements->  Configure(elements,loads, nodes,vertices, materials,parameters);
    … …  
    2324        /*retrieve some parameters: */
    2425        parameters->FindParam(&dim,DimEnum);
     26        parameters->FindParam(&numberofvertices,NumberOfVerticesEnum);
    2527
    2628        /*Compute gradients: */
    2729        for (i=0;i<elements->Size();i++){
    2830                Element* element=(Element*)elements->GetObjectByOffset(i);
    29                 element->Gradj(control_type);
     31                element->Gradj(gradient,control_type);
    3032        }
    3133
     34        /*Assemble vector: */
     35        VecAssemblyBegin(gradient);
     36        VecAssemblyEnd(gradient);
     37
    3238        /*Extrude if needed: */
    33         if(dim==3) ExtrudeInputx( elements,nodes, vertices, loads, materials, parameters,GradientEnum);
     39        if(dim==3) VecExtrudex(gradient, elements,nodes, vertices, loads, materials, parameters,0);
    3440
     41        /*Assign output pointers: */
     42        *pgradient=gradient;
    3543}

  • issm/trunk/src/c/modules/Gradjx/Gradjx.h

    r4047 r4055  
    1010
    1111/* local prototypes: */
    12 void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,  Parameters* parameters, int control_type);
     12void Gradjx(Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,  Parameters* parameters, int control_type);
    1313
    1414#endif  /* _GRADJX_H */

  • issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp

    r4034 r4055  
    5757                        break;
    5858
    59                 case SlopeComputeAnalysisEnum:
     59                case SlopeAnalysisEnum:
    6060                        CreateNodesSlopeCompute(pnodes, iomodel,iomodel_handle);
    6161                        CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle);

  • issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r4034 r4055  
    3737        if (iomodel->dim==2) parameters->AddObject(new IntParam(DimEnum,2));
    3838        else parameters->AddObject(new IntParam(DimEnum,3));
    39         parameters->AddObject(new StringParam(OutputFileNameEnum,iomodel->outputfilename));
    4039        parameters->AddObject(new BoolParam(IsHutterEnum,iomodel->ishutter));
    4140        parameters->AddObject(new BoolParam(IsMacAyealPattynEnum,iomodel->ismacayealpattyn));

  • issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateNodesSlopeCompute.cpp

    r4025 r4055  
    4646                       
    4747                        /*Add node to nodes dataset: */
    48                         nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,iomodel,SlopeComputeAnalysisEnum));
     48                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,iomodel,SlopeAnalysisEnum));
    4949
    5050                }

  • issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp

    r4042 r4055  
    1515#include "../../objects/objects.h"
    1616               
    17 void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters, char* filename){
     17void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters){
    1818
    1919        int i;
    2020
    2121        Patch* patch=NULL;
     22        char*  filename=NULL;
    2223
    2324        int solutiontype;
    … …  
    3738        //Recover solutiontype which will be output to disk:
    3839        parameters->FindParam(&solutiontype,SolutionTypeEnum);
     40        parameters->FindParam(&filename,OutputFileNameEnum);
    3941
    4042        //Process results to be output in the correct units
    … …  
    104106        /*Free ressources:*/
    105107        delete patch;
     108        xfree((void**)&filename);
    106109
    107110}

  • issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.h

    r4042 r4055  
    99
    1010/* local prototypes: */
    11 void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters, char* filename);
     11void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters);
    1212
    1313#endif  /* _OUTPUTRESULTS_H */

  • issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.cpp

    r3913 r4055  
    1111#include "../../DataSet/DataSet.h"
    1212
    13 void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface,
    14                 DataSet* elements, DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters,
    15                 Vec vec_newthickness,Vec vec_bed,Vec vec_surface){
     13void UpdateGeometryx( DataSet* elements, DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters){
    1614
    1715        int i;
    1816        int dof;
    1917
    20         /*serialized input: */
     18        /*vectorized inputs: */
     19        Vec     vec_newthickness=NULL;
     20        Vec     vec_bed=NULL;
     21        Vec     vec_surface=NULL:
    2122        double* newthickness=NULL;
    2223        double* bed=NULL;
    … …  
    2425
    2526        /*objects: */
    26         Node* node=NULL;
     27        Vertex* vertex=NULL;
    2728        Object* object=NULL;
    2829        Matpar* matpar=NULL;
    … …  
    3132        double h,b,s;
    3233        double rho_ice,rho_water;
    33 
    34         /*output: */
    35         Vec outthickness=NULL;
    36         Vec outbed=NULL;
    37         Vec outsurface=NULL;
    38 
    39         /*Duplicate inputs to outputs, do not copy values!: */
    40         VecDuplicate(vec_newthickness,&outthickness);
    41         VecDuplicate(vec_bed,&outbed);
    42         VecDuplicate(vec_surface,&outsurface);
    4334
    4435        /*First, get elements and loads configured: */
    … …  
    4940        parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
    5041
    51         /*Serialize inputs :*/
     42        /*retrieve inputs: */
     43        GetVectorFromInputsx( &vec_newthickness, elements, nodes,  vertices,  loads,  materials,  parameters, ThicknessEnum,VertexEnum);
     44        GetVectorFromInputsx( &vec_bed, elements, nodes,  vertices,  loads,  materials,  parameters, ThicknessEnum,VertexEnum);
     45        GetVectorFromInputsx( &vec_surface, elements, nodes,  vertices,  loads,  materials,  parameters, ThicknessEnum,VertexEnum);
     46
    5247        VecToMPISerial(&newthickness,vec_newthickness);
    5348        VecToMPISerial(&bed,vec_bed);
    … …  
    6762        rho_water=matpar->GetRhoWater();
    6863
    69         /*Go through nodes and for each node, update the thickness, bed and surface, using the
     64        /*Go through vertices and for each vertex, update the thickness, bed and surface, using the
    7065         * new thickness computed in the prognostic_core part of the transient solution sequence: */
    7166
    72         for(i=0;i<nodes->Size();i++){
     67        for(i=0;i<vertices->Size();i++){
    7368
    7469                /*get i'th node: */
    75                 node=(Node*)nodes->GetObjectByOffset(i);
     70                vertex=(Vertex*)vertices->GetObjectByOffset(i);
    7671
    77                 /*first, recover thickness, surface and bed for this node: */
    78                 dof=node->GetDofList1();
     72                /*first, recover thickness, surface and bed for this vertex: */
     73                dof=vertex->GetDofList1();
    7974                h=newthickness[dof];
    8075                s=surface[dof];
    … …  
    8580
    8681                //For grids on ice sheet, the new bed remains the same, the new surface becomes bed+new_thickness.
    87                 if (node->IsOnSheet()){
     82                if (vertex->IsOnSheet()){
    8883                        s=b+h;
    8984                }
    9085
    9186                //For grids on ice shelt, we have hydrostatic equilibrium (for now)
    92                 if (node->IsOnShelf()){
     87                if (vertex->IsOnShelf()){
    9388                        b=-rho_ice/rho_water*h;
    9489                        s=(1-rho_ice/rho_water)*h;
    … …  
    9691
    9792                /*Ok, plug values of thickness, surafce and bed into our outputs: */
    98                 VecSetValues(outthickness,1,&dof,&h,INSERT_VALUES);
    99                 VecSetValues(outbed,1,&dof,&b,INSERT_VALUES);
    100                 VecSetValues(outsurface,1,&dof,&s,INSERT_VALUES);
     93                VecSetValues(vec_newthickness,1,&dof,&h,INSERT_VALUES);
     94                VecSetValues(vec_bed,1,&dof,&b,INSERT_VALUES);
     95                VecSetValues(vec_surface,1,&dof,&s,INSERT_VALUES);
    10196        }
    10297
    10398        /*Assemble vectors: */
    104         VecAssemblyBegin(outthickness);
    105         VecAssemblyEnd(outthickness);
     99        VecAssemblyBegin(vec_newthickness);
     100        VecAssemblyEnd(vec_newthickness);
    106101
    107         VecAssemblyBegin(outbed);
    108         VecAssemblyEnd(outbed);
     102        VecAssemblyBegin(vec_bed);
     103        VecAssemblyEnd(vec_bed);
    109104
    110         VecAssemblyBegin(outsurface);
    111         VecAssemblyEnd(outsurface);
     105        VecAssemblyBegin(vec_surface);
     106        VecAssemblyEnd(vec_surface);
     107
     108        /*Put back into inputs:*/
     109        UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials,  parameters,vec_newthickness,ThicknessEnum,VertexEnum); //clobber existing thickness input with new one
     110        UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials,  parameters,vec_bed,BedEnum,VertexEnum);
     111        UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials,  parameters,vec_surface,SurfaceEnum,VertexEnum);
    112112
    113113        /*Free ressources:*/
     114        VecFree(&vec_newthickness);
     115        VecFree(&vec_bed);
     116        VecFree(&vec_surface);
    114117        xfree((void**)&newthickness);
    115118        xfree((void**)&bed);
    116119        xfree((void**)&surface);
    117120       
    118         /*Assign output pointers: */
    119         *poutthickness=outthickness;
    120         *poutbed=outbed;
    121         *poutsurface=outsurface;
    122121}

  • issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.h

    r3913 r4055  
    1010
    1111/* local prototypes: */
    12 void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface,
    13                 DataSet* elements, DataSet* nodes,DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters,
    14                 Vec newthickness,Vec bed,Vec surface);
     12void UpdateGeometryx(DataSet* elements, DataSet* nodes,DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters);
    1513
    1614#endif  /* _UPDATEGEOMETRYX_H */

  • issm/trunk/src/c/modules/modules.h

    r4053 r4055  
    5353#include "./ComputePressurex/ComputePressurex.h"
    5454#include "./ComputeStrainRatex/ComputeStrainRatex.h"
    55 #include "./FieldExtrudex/FieldExtrudex.h"
     55#include "./VecExtrudex/VecExtrudex.h"
    5656#include "./Qmux/Qmux.h"
    5757#include "./NodeConnectivityx/NodeConnectivityx.h"
    … …  
    8181#include "./ScaleInputx/ScaleInputx.h"
    8282#include "./AXPYInputx/AXPYInputx.h"
     83#include "./GetVectorFromInputsx/GetVectorFromInputsx.h"
     84#include "./InputConvergencex/InputConvergencex.h"
    8385
    8486#endif

  • issm/trunk/src/c/objects/Elements/Beam.cpp

    r4050 r4055  
    562562/*}}}*/
    563563/*FUNCTION Beam::Gradj{{{1*/
    564 void  Beam::Gradj(int control_type){
     564void  Beam::Gradj(Vec gradient,int control_type){
    565565        ISSMERROR(" not supported yet!");
    566566}
    567567/*}}}*/
    568568/*FUNCTION Beam::GradjB{{{1*/
    569 void  Beam::GradjB(void){
     569void  Beam::GradjB(Vec gradient){
    570570        ISSMERROR(" not supported yet!");
    571571}
    572572/*}}}*/
    573573/*FUNCTION Beam::GradjDrag{{{1*/
    574 void  Beam::GradjDrag(void){
     574void  Beam::GradjDrag(Vec gradient){
    575575        ISSMERROR(" not supported yet!");
    576576}
    … …  
    10551055}
    10561056/*}}}*/
     1057/*FUNCTION Beam::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
     1058void  Beam::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
     1059
     1060        Input** new_inputs=NULL;
     1061        Input** old_inputs=NULL;
     1062        int     converged=1;
     1063
     1064        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
     1065        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
     1066       
     1067        for(i=0;i<num_enums/2;i++){
     1068                new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
     1069                old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
     1070                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
     1071                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
     1072        }
     1073
     1074        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
     1075        for(i=0;i<num_criterionenums;i++){
     1076                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
     1077                if(eps[i]>criterionvalues[i]) converged=0;
     1078        }
     1079
     1080        /*Assign output pointers:*/
     1081        *pconverged=converged;
     1082
     1083}
     1084/*}}}*/

  • issm/trunk/src/c/objects/Elements/Beam.h

    r4050 r4055  
    6262                void  UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");}
    6363                void  UpdateInputsFromSolution(double* solution);
    64                 void  DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");};
     64                void  DepthAverageInputAtBase(int enum_type,int average_enum_type){ISSMERROR("not implemented yet");};
    6565                void  InputToResult(int enum_type,int step,double time);
    6666                void   ProcessResultsUnits(void);
    … …  
    9898                void  AXPYInput(int YEnum, double scalar, int XEnum);
    9999                void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
     100                void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
    100101
    101102                /*}}}*/
    … …  
    106107                void  GetThicknessList(double* thickness_list);
    107108                void  Du(Vec);
    108                 void  Gradj(int control_type);
    109                 void  GradjDrag(void);
    110                 void  GradjB(void);
     109                void  Gradj(Vec gradient,int control_type);
     110                void  GradjDrag(Vec gradient);
     111                void  GradjB(Vec gradient);
    111112                double Misfit(void);
    112113                double SurfaceArea(void);

  • issm/trunk/src/c/objects/Elements/Element.h

    r4050 r4055  
    3737                virtual void   GetBedList(double* bed_list)=0;
    3838                virtual void   Du(Vec du_g)=0;
    39                 virtual void   Gradj(int control_type)=0;
    40                 virtual void   GradjDrag(void)=0;
    41                 virtual void   GradjB(void)=0;
     39                virtual void   Gradj(Vec gradient,int control_type)=0;
     40                virtual void   GradjDrag(Vec gradient)=0;
     41                virtual void   GradjB(Vec gradient)=0;
    4242                virtual double Misfit(void)=0;
    4343                virtual double CostFunction(void)=0;
    4444                virtual double SurfaceArea(void)=0;
    45                 virtual void   DepthAverageInputAtBase(int enum_type)=0;
     45                virtual void   DepthAverageInputAtBase(int enum_type,int average_enum_type)=0;
    4646                virtual void   ComputeBasalStress(Vec sigma_b)=0;
    4747                virtual void   ComputePressure(Vec p_g)=0;
    … …  
    6969                virtual void   AXPYInput(int YEnum, double scalar, int XEnum)=0;
    7070                virtual void   ControlConstrainInput(int control_type,double cm_min, double cm_max)=0;
     71                virtual void   InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double criterionvalues,double* criterionvalues,int num_criterionenums)=0;
     72
    7173                /*Implementation: */
    7274

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4050 r4055  
    658658
    659659        }
    660         else if (analysis_type==SlopeComputeAnalysisEnum){
     660        else if (analysis_type==SlopeAnalysisEnum){
    661661
    662662                UpdateInputsFromSolutionSlopeCompute( solution);
    … …  
    12701270                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
    12711271        }
    1272         else if (analysis_type==SlopeComputeAnalysisEnum){
     1272        else if (analysis_type==SlopeAnalysisEnum){
    12731273                CreateKMatrixSlopeCompute( Kgg);
    12741274        }
    … …  
    23572357                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
    23582358        }
    2359         else if (analysis_type==SlopeComputeAnalysisEnum){
     2359        else if (analysis_type==SlopeAnalysisEnum){
    23602360                CreatePVectorSlopeCompute( pg);
    23612361        }
    … …  
    32633263}
    32643264/*}}}*/
    3265 /*FUNCTION Penta::FieldExtrude {{{1*/
    3266 void  Penta::FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed){
     3265/*FUNCTION Penta::VecExtrude {{{1*/
     3266void  Penta::VecExtrude(Vec vector,double* vector_serial,int iscollapsed){
    32673267
    32683268        /* node data: */
    3269         const int    numgrids=6;
    3270         int          numberofdofspernode;
     3269        int   i;
    32713270        Node* node=NULL;
    3272         int   i;
    32733271        int   extrude=0;
    32743272
    … …  
    32963294        if(onbed==0)extrude=0;
    32973295
    3298         /*Go on and extrude field: */
     3296        /*Go on and extrude vector: */
    32993297        if (extrude){
    33003298
    3301                 if (strcmp(field_name,"velocity")==0){
    3302 
    3303                         /* node data: */
    3304                         const int    numdof=2*numgrids;
    3305                         int          doflist[numdof];
    3306                         int          nodedofs[2];
    3307                         double       fieldel[2];
    3308 
    3309 
    3310                         GetDofList(&doflist[0],&numberofdofspernode);
    3311 
    3312                         /*this penta is a collapsed macayeal. For each node on the base of this penta,
    3313                          * we grab the field. Once we know the field, we follow the upper nodes,
    3314                          * inserting the same field value into field, until we reach the surface: */
    3315                         for(i=0;i<3;i++){
    3316 
    3317                                 node=nodes[i]; //base nodes
    3318 
    3319                                 /*get field for this base node: */
    3320                                 fieldel[0]=field_serial[doflist[numberofdofspernode*i+0]];
    3321                                 fieldel[1]=field_serial[doflist[numberofdofspernode*i+1]];
    3322 
    3323                                 //go throfieldn all nodes which sit on top of this node, until we reach the surface,
    3324                                 //and plfield  field in field
    3325                                 for(;;){
    3326 
    3327                                         node->GetDofList(&nodedofs[0],&numberofdofspernode);
    3328                                         VecSetValues(field,1,&nodedofs[0],&fieldel[0],INSERT_VALUES);
    3329                                         VecSetValues(field,1,&nodedofs[1],&fieldel[1],INSERT_VALUES);
    3330 
    3331                                         if (node->IsOnSurface())break;
    3332                                         /*get next node: */
    3333                                         node=node->GetUpperNode();
    3334                                 }
    3335                         }
    3336                 } //if (strcmp(field_name,"velocity")==0)
    3337                 else if (strcmp(field_name,"gradj")==0){
    3338 
    3339                         /* node data: */
    3340                         int          dof1;
    3341                         double       fieldel;
    3342 
    3343                         /*this penta is a collapsed macayeal. For each node on the base of this penta,
    3344                          * we grab the field. Once we know the field, we follow the upper nodes,
    3345                          * inserting the same field value into field, until we reach the surface: */
    3346                         for(i=0;i<3;i++){
    3347 
    3348                                 node=nodes[i]; //base nodes
     3299                /* node data: */
     3300                int          dof1;
     3301                double       vectorel;
     3302
     3303                /*this penta is a collapsed macayeal. For each node on the base of this penta,
     3304                 * we grab the vector. Once we know the vector, we follow the upper nodes,
     3305                 * inserting the same vector value into vector, until we reach the surface: */
     3306                for(i=0;i<3;i++){
     3307
     3308                        node=nodes[i]; //base nodes
     3309                        dof1=node->GetDofList1();
     3310
     3311                        /*get vector for this base node: */
     3312                        vectorel=vector_serial[dof1];
     3313
     3314                        //go throvectorn all nodes which sit on top of this node, until we reach the surface,
     3315                        //and plvector  vector in vector
     3316                        for(;;){
     3317
    33493318                                dof1=node->GetDofList1();
    3350 
    3351                                 /*get field for this base node: */
    3352                                 fieldel=field_serial[dof1];
    3353 
    3354                                 //go throfieldn all nodes which sit on top of this node, until we reach the surface,
    3355                                 //and plfield  field in field
    3356                                 for(;;){
    3357 
    3358                                         dof1=node->GetDofList1();
    3359                                         VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES);
    3360 
    3361                                         if (node->IsOnSurface())break;
    3362                                         /*get next node: */
    3363                                         node=node->GetUpperNode();
    3364                                 }
     3319                                VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES);
     3320
     3321                                if (node->IsOnSurface())break;
     3322                                /*get next node: */
     3323                                node=node->GetUpperNode();
    33653324                        }
    33663325                }
    3367                 else if (
    3368                                         (strcmp(field_name,"thickness")==0) ||
    3369                                         (strcmp(field_name,"surface")==0)  ||
    3370                                         (strcmp(field_name,"bed")==0)  ||
    3371                                         (strcmp(field_name,"slopex")==0)  ||
    3372                                         (strcmp(field_name,"slopey")==0)
    3373                                   ){
    3374 
    3375                         /* node data: */
    3376                         const int    numdof=1*numgrids;
    3377                         int          doflist[numdof];
    3378                         int          nodedofs;
    3379                         double       fieldel;
    3380 
    3381                         GetDofList(&doflist[0],&numberofdofspernode);
    3382 
    3383                         /*this penta is on the bed. For each node on the base of this penta,
    3384                          * we grab the thickness. Once we know the thickness, we follow the upper nodes,
    3385                          * inserting the same thickness value into tg, until we reach the surface: */
    3386                         for(i=0;i<3;i++){
    3387 
    3388                                 node=nodes[i]; //base nodes
    3389 
    3390                                 /*get velocity for this base node: */
    3391                                 fieldel=field_serial[doflist[numberofdofspernode*i+0]];
    3392 
    3393                                 //go through all nodes which sit on top of this node, until we reach the surface,
    3394                                 //and pltg  fieldel in field:
    3395                                 for(;;){
    3396 
    3397                                         node->GetDofList(&nodedofs,&numberofdofspernode);
    3398                                         VecSetValues(field,1,&nodedofs,&fieldel,INSERT_VALUES);
    3399 
    3400                                         if (node->IsOnSurface())break;
    3401                                         /*get next node: */
    3402                                         node=node->GetUpperNode();
    3403                                 }
    3404                         }
    3405 
    3406                 }
    3407                 else ISSMERROR("%s%s%s"," field ",field_name," not supported yet!");
    3408 
    3409         } //if (extrude)
     3326        }
    34103327}
    34113328/*}}}*/
    … …  
    34513368/*}}}*/
    34523369/*FUNCTION Penta::DepthAverageInputAtBase{{{1*/
    3453 void  Penta::DepthAverageInputAtBase(int enum_type){
     3370void  Penta::DepthAverageInputAtBase(int enum_type,int average_enum_type){
    34543371        ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)");
    34553372}
    … …  
    45554472/*}}}*/
    45564473/*FUNCTION Penta::Gradj {{{1*/
    4557 void  Penta::Gradj(int control_type){
     4474void  Penta::Gradj(Vec gradient,int control_type){
    45584475
    45594476        /*inputs: */
    … …  
    45674484
    45684485        if (control_type==DragCoefficientEnum){
    4569                 GradjDrag();
     4486                GradjDrag(gradient);
    45704487        }
    45714488        else if (control_type=RheologyBEnum){
    4572                 GradjB();
     4489                GradjB(gradient);
    45734490        }
    45744491        else ISSMERROR("%s%i","control type not supported yet: ",control_type);
    … …  
    45764493/*}}}*/
    45774494/*FUNCTION Penta::GradjDrag {{{1*/
    4578 void  Penta::GradjDrag(void){
     4495void  Penta::GradjDrag(Vec gradient){
    45794496
    45804497        int i;
    45814498        Tria* tria=NULL;
    45824499        TriaVertexInput* triavertexinput=NULL;
    4583         double gradient[6]={0,0,0,0,0,0};
     4500        double temp_gradient[6]={0,0,0,0,0,0};
    45844501
    45854502        /*inputs: */
    … …  
    46114528                /*MacAyeal or Pattyn*/
    46124529                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    4613                 tria->GradjDrag();
     4530                tria->GradjDrag(gradient);
     4531                delete tria;
    46144532        }
    46154533        else if (sub_analysis_type==StokesAnalysisEnum){
    … …  
    46174535                /*Stokes*/
    46184536                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    4619                 tria->GradjDragStokes();
     4537                tria->GradjDragStokes(gradient);
     4538                delete tria;
    46204539        }
    46214540        else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
    46224541
    4623         /*Now, the tria  has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */
    4624         triavertexinput=(TriaVertexInput*)tria->inputs->GetInput(GradientEnum);
    4625         for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i];
    4626         this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0]));
    4627        
    4628         delete tria;
    46294542
    46304543}
    46314544/*}}}*/
    46324545/*FUNCTION Penta::GradjB {{{1*/
    4633 void  Penta::GradjB(void){
     4546void  Penta::GradjB(Vec gradient){
    46344547
    46354548        int i;
    46364549        Tria* tria=NULL;
    46374550        TriaVertexInput* triavertexinput=NULL;
    4638         double gradient[6]={0,0,0,0,0,0};
    46394551
    46404552        /*inputs: */
    … …  
    46584570                 * and compute gardj*/
    46594571                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    4660                 tria->GradjB();
     4572                tria->GradjB(gradient);
     4573                delete tria;
    46614574        }
    46624575        else{
    46634576                /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
    46644577                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    4665                 tria->GradjB();
     4578                tria->GradjB(gradient);
     4579                delete tria;
    46664580        }
    46674581       
    4668         /*Now, the tria  has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */
    4669         triavertexinput=(TriaVertexInput*)tria->inputs->GetInput(GradientEnum);
    4670         for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i];
    4671         this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0]));
    4672        
    4673         delete tria;
    46744582
    46754583}
    … …  
    54075315}
    54085316/*}}}*/
     5317/*FUNCTION Penta::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
     5318void  Penta::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
     5319
     5320        Input** new_inputs=NULL;
     5321        Input** old_inputs=NULL;
     5322        int     converged=1;
     5323
     5324        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
     5325        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
     5326       
     5327        for(i=0;i<num_enums/2;i++){
     5328                new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
     5329                old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
     5330                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
     5331                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
     5332        }
     5333
     5334        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
     5335        for(i=0;i<num_criterionenums;i++){
     5336                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
     5337                if(eps[i]>criterionvalues[i]) converged=0;
     5338        }
     5339
     5340        /*Assign output pointers:*/
     5341        *pconverged=converged;
     5342
     5343}
     5344/*}}}*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r4050 r4055  
    8686                bool   GetOnBed();
    8787                void  Du(Vec du_gg);
    88                 void  Gradj(int control_type);
    89                 void  GradjDrag(void);
    90                 void  GradjB(void);
     88                void  Gradj(Vec gradient,int control_type);
     89                void  GradjDrag(Vec gradient);
     90                void  GradjB(Vec gradient);
    9191                double Misfit(void);
    9292                double SurfaceArea(void);
    … …  
    110110                void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord);
    111111                void  GetNodalFunctions(double* l1l6, double* gauss_coord);
    112                 void  FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed);
     112                void  VecExtrude(Vec vector,double* vector_serial,int iscollapsed);
    113113                void  InputExtrude(int enum_type);
    114                 void  DepthAverageInputAtBase(int enum_type);
     114                void  DepthAverageInputAtBase(int enum_type,int average_enum_type);
    115115                void  ComputeBasalStress(Vec sigma_bg);
    116116                void  ComputePressure(Vec p_gg);
    … …  
    165165                void  AXPYInput(int YEnum, double scalar, int XEnum);
    166166                void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
     167                void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
    167168                void  GetVectorFromInputs(Vec vector,int NameEnum);
    168169

  • issm/trunk/src/c/objects/Elements/Sing.cpp

    r4050 r4055  
    383383/*}}}*/
    384384/*FUNCTION Sing::Gradj {{{1*/
    385 void  Sing::Gradj(int control_type){
     385void  Sing::Gradj(Vec gradient,int control_type){
    386386        ISSMERROR(" not supported yet!");
    387387}
    388388/*}}}*/
    389389/*FUNCTION Sing::GradB {{{1*/
    390 void  Sing::GradjB(void){
     390void  Sing::GradjB(Vec gradient){
    391391        ISSMERROR(" not supported yet!");
    392392}
    393393/*}}}*/
    394394/*FUNCTION Sing::GradjDrag {{{1*/
    395 void  Sing::GradjDrag(void){
     395void  Sing::GradjDrag(Vec gradient){
    396396        ISSMERROR(" not supported yet!");
    397397}
    … …  
    761761}
    762762/*}}}*/
     763/*FUNCTION Sing::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
     764void  Sing::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
     765
     766        Input** new_inputs=NULL;
     767        Input** old_inputs=NULL;
     768        int     converged=1;
     769
     770        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
     771        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
     772       
     773        for(i=0;i<num_enums/2;i++){
     774                new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
     775                old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
     776                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
     777                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
     778        }
     779
     780        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
     781        for(i=0;i<num_criterionenums;i++){
     782                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
     783                if(eps[i]>criterionvalues[i]) converged=0;
     784        }
     785
     786        /*Assign output pointers:*/
     787        *pconverged=converged;
     788
     789}
     790/*}}}*/

  • issm/trunk/src/c/objects/Elements/Sing.h

    r4050 r4055  
    6262                void  UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");}
    6363                void  UpdateInputsFromSolution(double* solutiong);
    64                 void  DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");};
     64                void  DepthAverageInputAtBase(int enum_type,int average_enum_type){ISSMERROR("not implemented yet");};
    6565                void  InputToResult(int enum_type,int step,double time);
    6666                void   ProcessResultsUnits(void);
    … …  
    9696                void  ScaleInput(int enum_type,double scale_factor);
    9797                void  AXPYInput(int YEnum, double scalar, int XEnum);
     98                void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
    9899                void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
    99100
    … …  
    106107                void  GetThicknessList(double* thickness_list);
    107108                void  Du(Vec);
    108                 void  Gradj(int control_type);
    109                 void  GradjDrag(void);
    110                 void  GradjB(void);
     109                void  Gradj(Vec gradient,int control_type);
     110                void  GradjDrag(Vec gradient);
     111                void  GradjB(Vec gradient);
    111112                double Misfit(void);
    112113                double SurfaceArea(void);

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r4050 r4055  
    569569                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
    570570        }
    571         else if (analysis_type==SlopeComputeAnalysisEnum){
     571        else if (analysis_type==SlopeAnalysisEnum){
    572572                UpdateInputsFromSolutionSlopeCompute( solution);
    573573        }
    … …  
    979979                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
    980980        }
    981         else if (analysis_type==SlopeComputeAnalysisEnum){
     981        else if (analysis_type==SlopeAnalysisEnum){
    982982                CreateKMatrixSlopeCompute( Kgg);
    983983        }
    … …  
    23942394                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
    23952395        }
    2396         else if (analysis_type==SlopeComputeAnalysisEnum){
     2396        else if (analysis_type==SlopeAnalysisEnum){
    23972397                CreatePVectorSlopeCompute( pg);
    23982398        }
    … …  
    41254125/*}}}*/
    41264126/*FUNCTION Tria::Gradj {{{1*/
    4127 void  Tria::Gradj(int control_type){
     4127void  Tria::Gradj(Vec gradient,int control_type){
    41284128
    41294129        /*inputs: */
    … …  
    41374137
    41384138        if (control_type==DragCoefficientEnum){
    4139                 GradjDrag();
     4139                GradjDrag(gradient);
    41404140        }
    41414141        else if (control_type==RheologyBEnum){
    4142                 GradjB();
     4142                GradjB(gradient);
    41434143        }
    41444144        else ISSMERROR("%s%i","control type not supported yet: ",control_type);
    … …  
    41464146/*}}}*/
    41474147/*FUNCTION Tria::GradjDrag {{{1*/
    4148 void  Tria::GradjDrag(void){
     4148void  Tria::GradjDrag(Vec gradient){
    41494149
    41504150
    … …  
    42954295
    42964296
    4297         /*Add grade_g to the inputs of this element: */
    4298         this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
     4297        /*Add grade_g to global vector gradient: */
     4298        VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    42994299
    43004300        cleanup_and_return:
    … …  
    43084308/*}}}*/
    43094309/*FUNCTION Tria::GradjDragStokes {{{1*/
    4310 void  Tria::GradjDragStokes(void){
     4310void  Tria::GradjDragStokes(Vec gradient){
    43114311
    43124312        int i;
    … …  
    44654465        }
    44664466
     4467        /*Add grade_g to global vector gradient: */
     4468        VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
     4469
    44674470        /*Add grade_g to the inputs of this element: */
    44684471        this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
    … …  
    44784481/*}}}*/
    44794482/*FUNCTION Tria::GradjB{{{1*/
    4480 void  Tria::GradjB(void){
     4483void  Tria::GradjB(Vec gradient){
    44814484
    44824485        int i;
    … …  
    46134616        }
    46144617
    4615         /*Add grade_g to the inputs of this element: */
    4616         this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
     4618        /*Add grade_g to global vector gradient: */
     4619        VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    46174620
    46184621        cleanup_and_return:
    … …  
    50145017/*}}}*/
    50155018/*FUNCTION Tria::DepthAverageInputAtBase {{{1*/
    5016 void  Tria::DepthAverageInputAtBase(int enum_type){
     5019void  Tria::DepthAverageInputAtBase(int enum_type,int average_enum_type){
    50175020
    50185021        /*New input*/
    … …  
    50265029
    50275030        /*Assign new name (average)*/
    5028         switch(enum_type){
    5029 
    5030                 case VxEnum:
    5031                         newinput->ChangeEnum(VxAverageEnum);
    5032                         break;
    5033 
    5034                 case VyEnum:
    5035                         newinput->ChangeEnum(VyAverageEnum);
    5036                         break;
    5037 
    5038                 default:
    5039                         ISSMERROR("enum_type: %i (%s) not suppoerted yet",enum_type,EnumAsString(enum_type));
    5040         }
     5031        newinput->ChangeEnum(average_enum_type);
    50415032
    50425033        /*Add new input to current element*/
    … …  
    55035494}
    55045495/*}}}*/
     5496/*FUNCTION Tria::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
     5497void  Tria::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
     5498
     5499        Input** new_inputs=NULL;
     5500        Input** old_inputs=NULL;
     5501        int     converged=1;
     5502
     5503        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
     5504        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
     5505       
     5506        for(i=0;i<num_enums/2;i++){
     5507                new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
     5508                old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
     5509                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
     5510                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
     5511        }
     5512
     5513        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
     5514        for(i=0;i<num_criterionenums;i++){
     5515                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
     5516                if(eps[i]>criterionvalues[i]) converged=0;
     5517        }
     5518
     5519        /*Assign output pointers:*/
     5520        *pconverged=converged;
     5521
     5522}
     5523/*}}}*/

  • issm/trunk/src/c/objects/Elements/Tria.h

    r4050 r4055  
    5555                int   MyRank(void);
    5656                void  SetClone(int* minranks);
    57                 void  DepthAverageInputAtBase(int enum_type);
     57                void  DepthAverageInputAtBase(int enum_type,int average_enum_type);
    5858                void*  SpawnSing(int g0);
    5959                void*  SpawnBeam(int g0, int g1);
    … …  
    8888                void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
    8989                void  Du(Vec du_g);
    90                 void  Gradj(int control_type);
    91                 void  GradjDrag(void);
    92                 void  GradjDragStokes(void);
    93                 void  GradjB(void);
     90                void  Gradj(Vec gradient,int control_type);
     91                void  GradjDrag(Vec gradient);
     92                void  GradjDragStokes(Vec gradient);
     93                void  GradjB(Vec gradient);
    9494                void  SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
    9595                double Misfit(void);
    … …  
    142142                void  ScaleInput(int enum_type,double scale_factor);
    143143                void  AXPYInput(int YEnum, double scalar, int XEnum);
     144                void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
    144145                void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
    145146                void  GetVectorFromInputs(Vec vector,int NameEnum);

  • issm/trunk/src/c/objects/FemModel.cpp

    r4050 r4055  
    3535        /*Dynamically allocate whatever is a list of length nummodels: */
    3636        analysis_type_list=(int*)xmalloc(nummodels*sizeof(int));
    37         Rmg=(Mat*)xmalloc(nummodels*sizeof(Mat));
    38         Gmn=(Mat*)xmalloc(nummodels*sizeof(Mat));
    39         nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*));
    40         yg=(Vec*)xmalloc(nummodels*sizeof(Vec));
    41         ys=(Vec*)xmalloc(nummodels*sizeof(Vec));
     37        m_Rmg=(Mat*)xmalloc(nummodels*sizeof(Mat));
     38        m_Gmn=(Mat*)xmalloc(nummodels*sizeof(Mat));
     39        m_nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*));
     40        m_yg=(Vec*)xmalloc(nummodels*sizeof(Vec));
     41        m_ys=(Vec*)xmalloc(nummodels*sizeof(Vec));
    4242
    4343        /*Initialize: */
    4444        for(i=0;i<nummodels;i++)analysis_type_list[i]=analyses[i];
    45         for(i=0;i<nummodels;i++)Rmg[i]=NULL;
    46         for(i=0;i<nummodels;i++)Gmn[i]=NULL;
    47         for(i=0;i<nummodels;i++)nodesets[i]=NULL;
    48         for(i=0;i<nummodels;i++)yg[i]=NULL;
    49         for(i=0;i<nummodels;i++)ys[i]=NULL;
     45        for(i=0;i<nummodels;i++)m_Rmg[i]=NULL;
     46        for(i=0;i<nummodels;i++)m_Gmn[i]=NULL;
     47        for(i=0;i<nummodels;i++)m_nodesets[i]=NULL;
     48        for(i=0;i<nummodels;i++)m_yg[i]=NULL;
     49        for(i=0;i<nummodels;i++)m_ys[i]=NULL;
    5050
    5151        _printf_("   fill model with matlab workspace data\n");
    … …  
    6464
    6565                _printf_("   create single point constraints: \n");
    66                 SpcNodesx( &yg[i], nodes,constraints,analysis_type);
     66                SpcNodesx( &m_yg[i], nodes,constraints,analysis_type);
    6767
    6868                _printf_("   create rigid body constraints:\n");
    69                 MpcNodesx( &Rmg[i], nodes,constraints,analysis_type);
     69                MpcNodesx( &m_Rmg[i], nodes,constraints,analysis_type);
    7070
    7171                _printf_("   create node sets:\n");
    72                 BuildNodeSetsx(&nodesets[i], nodes,analysis_type);
     72                BuildNodeSetsx(&m_nodesets[i], nodes,analysis_type);
    7373
    7474                _printf_("   reducing single point constraints vector:\n");
    75                 Reducevectorgtosx(&ys[i], yg[i],nodesets[i]);
     75                Reducevectorgtosx(&m_ys[i], m_yg[i],m_nodesets[i]);
    7676
    7777                _printf_("   normalizing rigid body constraints matrix:\n");
    78                 NormalizeConstraintsx(&Gmn[i], Rmg[i],nodesets[i]);
     78                NormalizeConstraintsx(&m_Gmn[i], m_Rmg[i],m_nodesets[i]);
    7979
    8080                _printf_("   configuring element and loads:\n");
    … …  
    106106
    107107        for(i=0;i<nummodels;i++){
    108                 Mat temp_Rmg=Rmg[i];
     108                Mat temp_Rmg=m_Rmg[i];
    109109                MatFree(&temp_Rmg);
    110                 Mat temp_Gmn=Gmn[i];
     110                Mat temp_Gmn=m_Gmn[i];
    111111                MatFree(&temp_Gmn);
    112                 NodeSets* temp_nodesets=nodesets[i];
    113                 delete nodesets;
    114                 Vec temp_yg=yg[i];
     112                NodeSets* temp_nodesets=m_nodesets[i];
     113                delete temp_nodesets;
     114                Vec temp_yg=m_yg[i];
    115115                VecFree(&temp_yg);
    116                 Vec temp_ys=ys[i];
     116                Vec temp_ys=m_ys[i];
    117117                VecFree(&temp_ys);
    118118        }
    119119
    120120        /*Delete dynamically allocated arrays: */
    121         delete Rmg;
    122         delete Gmn;
    123         delete nodesets;
    124         delete yg;
    125         delete ys;
     121        delete m_Rmg;
     122        delete m_Gmn;
     123        delete m_nodesets;
     124        delete m_yg;
     125        delete m_ys;
    126126
    127127}
    … …  
    151151/*FUNCTION FemModel::SetCurrentAnalysis {{{1*/
    152152void FemModel::SetCurrentAnalysis(int analysis_type){
     153
    153154        int found=-1;
    154155        for(int i=0;i<nummodels;i++){
    … …  
    160161        if(found!=-1) analysis_counter=found;
    161162        else ISSMERROR("Could not find analysis_type %s in list of FemModel analyses",EnumAsString(analysis_type));
     163
     164        /*activate matrices/vectors: */
     165        Rmg=m_Rmg[analysis_counter];
     166        Gmn=m_Gmn[analysis_counter];
     167        nodesets=m_nodesets[analysis_counter];
     168        yg=m_yg[analysis_counter];
     169        ys=m_ys[analysis_counter];
     170
     171}
     172/*}}}1*/
     173/*FUNCTION FemModel::SetCurrentAnalysisAlias {{{1*/
     174void FemModel::SetCurrentAnalysisAlias(int base_analysis_type,int real_analysis_type){
     175        int found=-1;
     176        for(int i=0;i<nummodels;i++){
     177                if (analysis_type_list[i]==base_analysis_type){
     178                        found=i;
     179                        break;
     180                }
     181        }
     182        if(found!=-1) analysis_counter=found;
     183        else ISSMERROR("Could not find alias for analysis_type %s in list of FemModel analyses",EnumAsString(base_analysis_type));
    162184}
    163185/*}}}1*/

  • issm/trunk/src/c/objects/FemModel.h

    r4050 r4055  
    4040                DofVec*             tpartition;
    4141               
    42                 //multiple  sets of matrices/vectors for each analysis_type
    43                 Mat*                 Rmg; //rigid body motions matrices
    44                 Mat*                 Gmn;
    45                 NodeSets**           nodesets; //boundary conditions dof sets
    46                 Vec*                 yg; //boundary conditions in global g-set
    47                 Vec*                 ys; //boundary conditions, in reduced s-set
     42                //multiple  sets of matrices/vectors for each analysis_type. m stands for multiple
     43                Mat*                 m_Rmg; //rigid body motions matrices
     44                Mat*                 m_Gmn;
     45                NodeSets**           m_nodesets; //boundary conditions dof sets
     46                Vec*                 m_yg; //boundary conditions in global g-set
     47                Vec*                 m_ys; //boundary conditions, in reduced s-set
     48
     49                //pointers to point to sets of matrices/vectors, for a certain analysis type. activated in SetCurrentAnalysis
     50                Mat                  Rmg;
     51                Mat                  Gmn;
     52                NodeSets*            nodesets;
     53                Vec                  yg;
     54                Vec                  ys;
    4855
    4956                /*constructors, destructors: */
    … …  
    5663                /*Fem: */
    5764                void  SetCurrentAnalysis(int analysis_type);
     65                void  SetCurrentAnalysisAlias(int base_analysis_type,int real_analysis_type);
    5866                int   GetCurrentAnalysis(void);
    5967

  • issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp

    r4050 r4055  
    294294}
    295295/*}}}*/
     296/*FUNCTION BeamVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
     297void BeamVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
     298
     299        *pvalues=this->values;
     300        *pnum_values=2;
     301
     302}
     303/*}}}*/

  • issm/trunk/src/c/objects/Inputs/BeamVertexInput.h

    r4050 r4055  
    7878                void Constrain(double cm_min, double cm_max);
    7979                void GetVectorFromInputs(Vec vector,int* doflist);
     80                void GetValuesPtr(double* pvalues,int* pnum_values);
    8081                /*}}}*/
    8182

  • issm/trunk/src/c/objects/Inputs/BoolInput.cpp

    r4050 r4055  
    256256}
    257257/*}}}*/
     258/*FUNCTION BoolInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
     259void BoolInput::GetValuesPtr(double* pvalues,int* pnum_values){
     260
     261        ISSMERROR(" not supported yet!");
     262
     263}
     264/*}}}*/

  • issm/trunk/src/c/objects/Inputs/BoolInput.h

    r4050 r4055  
    7878                void Constrain(double cm_min, double cm_max);
    7979                void GetVectorFromInputs(Vec vector,int* doflist);
     80                void GetValuesPtr(double* pvalues,int* pnum_values);
    8081                /*}}}*/
    8182

  • issm/trunk/src/c/objects/Inputs/DoubleInput.cpp

    r4050 r4055  
    267267}
    268268/*}}}*/
     269/*FUNCTION DoubleInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
     270void DoubleInput::GetValuesPtr(double* pvalues,int* pnum_values){
     271
     272        ISSMERROR(" not supported yet!");
     273
     274}
     275/*}}}*/

  • issm/trunk/src/c/objects/Inputs/DoubleInput.h

    r4050 r4055  
    7979                void Constrain(double cm_min, double cm_max);
    8080                void GetVectorFromInputs(Vec vector,int* doflist);
     81                void GetValuesPtr(double* pvalues,int* pnum_values);
    8182                /*}}}*/
    8283

  • issm/trunk/src/c/objects/Inputs/Input.h

    r4050 r4055  
    5555                virtual void Constrain(double cm_min, double cm_max)=0;
    5656                virtual void GetVectorFromInputs(Vec vector,int* doflist)=0;
     57                virtual void GetValuesPtr(double* pvalues,int* pnum_values)=0;
    5758                /*}}}*/
    5859

  • issm/trunk/src/c/objects/Inputs/IntInput.cpp

    r4050 r4055  
    257257}
    258258/*}}}*/
     259/*FUNCTION IntInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
     260void IntInput::GetValuesPtr(double* pvalues,int* pnum_values){
     261
     262        ISSMERROR(" not supported yet!");
     263
     264}
     265/*}}}*/

  • issm/trunk/src/c/objects/Inputs/IntInput.h

    r4050 r4055  
    7878                void Constrain(double cm_min, double cm_max);
    7979                void GetVectorFromInputs(Vec vector,int* doflist);
     80                void GetValuesPtr(double* pvalues,int* pnum_values);
    8081                /*}}}*/
    8182

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp

    r4050 r4055  
    943943}
    944944/*}}}*/
     945/*FUNCTION PentaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
     946void PentaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
     947
     948        *pvalues=this->values;
     949        *pnum_values=6;
     950
     951}
     952/*}}}*/

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.h

    r4050 r4055  
    8787                void Constrain(double cm_min, double cm_max);
    8888                void GetVectorFromInputs(Vec vector,int* doflist);
     89                void GetValuesPtr(double* pvalues,int* pnum_values);
    8990                /*}}}*/
    9091

  • issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp

    r4050 r4055  
    258258
    259259}
    260 /*}}}*/
     260/2*}}}*/
     261/*FUNCTION SingVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
     262void SingVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
     263
     264        *pvalues=&this->value;
     265        *pnum_values=1;
     266
     267}
     268/*}}}*/

  • issm/trunk/src/c/objects/Inputs/SingVertexInput.h

    r4050 r4055  
    7777                void Constrain(double cm_min, double cm_max);
    7878                void GetVectorFromInputs(Vec vector,int* doflist);
     79                void GetValuesPtr(double* pvalues,int* pnum_values);
    7980                /*}}}*/
    8081

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

    r4050 r4055  
    517517}
    518518/*}}}*/
     519/*FUNCTION TriaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
     520void TriaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
     521
     522        *pvalues=this->values;
     523        *pnum_values=3;
     524
     525}
     526/*}}}*/

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.h

    r4050 r4055  
    8484                void Constrain(double cm_min, double cm_max);
    8585                void GetVectorFromInputs(Vec vector,int* doflist);
     86                void GetValuesPtr(double* pvalues,int* pnum_values);
    8687                /*}}}*/
    8788

  • issm/trunk/src/c/objects/Node.cpp

    r4021 r4055  
    144144                                analysis_type==PrognosticAnalysisEnum ||
    145145                                analysis_type==MeltingAnalysisEnum ||
    146                                 analysis_type==SlopeComputeAnalysisEnum ||
     146                                analysis_type==SlopeAnalysisEnum ||
    147147                                analysis_type==BalancedvelocitiesAnalysisEnum ||
    148148                                analysis_type==BalancedthicknessAnalysisEnum
    … …  
    668668}
    669669/*}}}*/
    670 /*FUNCTION Node::FieldExtrude {{{2*/
    671 void  Node::FieldExtrude(Vec field,double* field_serial,char* field_name){
     670/*FUNCTION Node::VecExtrude {{{2*/
     671void  Node::VecExtrude(Vec vector,double* vector_serial){
    672672               
    673673        /* node data: */
    … …  
    684684        if (onbed){
    685685
    686                 if (strcmp(field_name,"velocity")==0){
    687 
    688                         /* node data: */
    689                         int          vertexdof;
    690                         int          doflist[2];
    691                         double       fieldnode[2];
    692 
    693                         vertexdof=GetVertexDof();
    694 
    695                         /*get dofs for this base node velocity: we know there are two dofs in field_serial */
    696                         doflist[0]=2*vertexdof;
    697                         doflist[1]=2*vertexdof+1;
    698 
    699                         //initilaize node
    700                         node=this;
    701                
    702                         /*get field for this base node: */
    703                         fieldnode[0]=field_serial[doflist[0]];
    704                         fieldnode[1]=field_serial[doflist[1]];
    705 
    706                         //go through all nodes which sit on top of this node, until we reach the surface,
    707                         //and plfield  field in field
    708                         for(;;){
    709 
    710                                 vertexdof=node->GetVertexDof();
    711                                 doflist[0]=2*vertexdof;
    712                                 doflist[1]=2*vertexdof+1;
    713                                 VecSetValues(field,1,&doflist[0],&fieldnode[0],INSERT_VALUES);
    714                                 VecSetValues(field,1,&doflist[1],&fieldnode[1],INSERT_VALUES);
    715 
    716                                 if (node->IsOnSurface())break;
    717                                 /*get next node: */
    718                                 node=node->GetUpperNode();
    719                         }
    720                 } //if (strcmp(field_name,"velocity")==0)
    721                 else if (strcmp(field_name,"gradj")==0){
    722 
    723                         /* node data: */
    724                         int          dof1;
    725                         double       fieldel;
    726 
    727                         //initilaize node and get dof1
    728                         node=this;
     686                /* node data: */
     687                int          dof1;
     688                double       vectorel;
     689
     690                //initilaize node and get dof1
     691                node=this;
     692                dof1=node->GetVertexDof();
     693
     694                /*get vector for this base node: */
     695                vectorel=vector_serial[dof1];
     696
     697                //go throvectorn all nodes which sit on top of this node, until we reach the surface,
     698                //and plvector  vector in vector
     699                for(;;){
     700
    729701                        dof1=node->GetVertexDof();
    730 
    731                         /*get field for this base node: */
    732                         fieldel=field_serial[dof1];
    733 
    734                         //go throfieldn all nodes which sit on top of this node, until we reach the surface,
    735                         //and plfield  field in field
    736                         for(;;){
    737 
    738                                 dof1=node->GetVertexDof();
    739                                 VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES);
    740 
    741                                 if (node->IsOnSurface())break;
    742                                 /*get next node: */
    743                                 node=node->GetUpperNode();
    744                         }
    745                 }
    746                 else if (
    747                                 (strcmp(field_name,"vx")==0) ||
    748                                 (strcmp(field_name,"vy")==0) ||
    749                                 (strcmp(field_name,"thickness")==0) ||
    750                                 (strcmp(field_name,"temperature")==0) ||
    751                                 (strcmp(field_name,"surface")==0)  ||
    752                                 (strcmp(field_name,"bed")==0)  ||
    753                                 (strcmp(field_name,"slopex")==0)  ||
    754                                 (strcmp(field_name,"slopey")==0)
    755                                 ){
    756 
    757                         /* node data: */
    758                         int          doflist;
    759                         int          nodedofs;
    760                         double       fieldnode;
    761                                
    762                         GetDofList(&doflist,&numberofdofspernode);
    763 
    764                         //initialize node
    765                         node=this;
    766 
    767                         /*get field for this bed node: */
    768                         fieldnode=field_serial[doflist];
    769 
    770                         //go through all nodes which sit on top of this node, until we reach the surface,
    771                         //and plug  fieldnode in field:
    772                         for(;;){
    773 
    774                                 node->GetDofList(&nodedofs,&numberofdofspernode);
    775                                 VecSetValues(field,1,&nodedofs,&fieldnode,INSERT_VALUES);
    776 
    777                                 if (node->IsOnSurface())break;
    778                                 /*get next node: */
    779                                 node=node->GetUpperNode();
    780                         }
    781                 }
    782                 else ISSMERROR("%s%s%s"," field ",field_name," not supported yet!");
    783 
    784         } //if (extrude)
     702                        VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES);
     703
     704                        if (node->IsOnSurface())break;
     705                        /*get next node: */
     706                        node=node->GetUpperNode();
     707                }
     708        }
    785709}
    786710/*}}}*/

  • issm/trunk/src/c/objects/Node.h

    r4043 r4055  
    8989                int   IsOnShelf();
    9090                int   IsOnSheet();
    91                 void  FieldExtrude(Vec field,double* field_serial,char* field_name);
     91                void  VecExtrude(Vec vector,double* vector_serial);
    9292                /*}}}*/
    9393                /*FUNCTION DofObject routines {{{1*/

  • issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp

    r4021 r4055  
    2525                numdofs=2;
    2626        }
    27         else if (analysis_type==SlopeComputeAnalysisEnum){
     27        else if (analysis_type==SlopeAnalysisEnum){
    2828                numdofs=1;
    2929        }

  • issm/trunk/src/c/shared/Numerics/numerics.h

    r3775 r4055  
    88#include "./GaussPoints.h"
    99#include "./isnan.h"
     10class Input;
    1011
    1112struct OptArgs;
    … …  
    1819void   cross(double* result,double* vector1,double* vector2);
    1920double norm(double* vector);
     21void IsInputConverged(double* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum);
    2022
    2123#endif //ifndef _NUMERICS_H_

  • issm/trunk/src/c/solutions/adjoint_core.cpp

    r4047 r4055  
    2020        char* solverstring=NULL;
    2121        bool  isstokes=false;
     22        bool conserve_loads=true;
    2223       
    2324        /*intermediary: */
    … …  
    3839
    3940        /*set analysis type to compute velocity: */
    40         if(isstokes)femmodel->SetAnalysis(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum);
    41         else femmodel->SetAnalysis(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum);
     41        if(isstokes)femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum);
     42        else femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
    4243       
    4344        _printf_("%s\n","      recover solution for this stiffness and right hand side:");
    44         solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0,NULL, femmodel,DiagnosticAnalysisEnum,sub_analysis_type);
     45        solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0, femmodel,conserve_loads);
    4546
    4647        _printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
    47         Dux(du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     48        Dux(&du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    4849
    4950        _printf_("%s\n","      reduce adjoint load from g-set to f-set:");
    … …  
    6667
    6768        /*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */
    68         if(isstokes)femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum);
    69         else femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum);
     69        if(isstokes)femmodel->SetCurrentAnalysisAlias(DiagnosticStokesAnalysisEnum,AdjointAnalysisEnum);
     70        else femmodel->SetCurrentAnalysisAlias(DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum);
    7071       
    71         UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g,AdjointAnalysisEnum, NoneAnalysisEnum);
     72        UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g);
    7273
    7374        /*Free ressources:*/
    74         char* solverstring=NULL;
    75         VecFree(&ug);
     75        xfree((void**)&solverstring);
     76        VecFree(&u_g);
    7677        MatFree(&K_ff0);
    77         MatFree(&Mat K_fs0);
     78        MatFree(&K_fs0);
    7879        VecFree(&du_g);
    7980        VecFree(&du_f);

  • issm/trunk/src/c/solutions/balancedthickness.cpp

    r3938 r4055  
    2222        char* inputfilename=NULL;
    2323        char* outputfilename=NULL;
     24        bool  qmu_analysis=false;
    2425        char* lockname=NULL;
    25         int   numberofnodes;
    2626        bool  waitonlock=false;
    2727
    28         Model* model=NULL;
    29 
    30         double  dt;
    31         double  yts;
    32         bool    qmu_analysis;
    33 
    34         /*Results: */
    35         Results* results=NULL;
    36         Param*   param=NULL;
     28        /*FemModel: */
     29        FemModel* femmodel=NULL;
    3730
    3831        /*time*/
    … …  
    4033        double   start_core, finish_core;
    4134        double   start_init, finish_init;
     35
     36        int analyses[1]={BalancedthicknessAnalysisEnum};
     37        int solution_type=BalancedthicknessAnalysisEnum;
    4238
    4339        MODULEBOOT();
    … …  
    6662        fid=pfopen(inputfilename,"rb");
    6763
    68         /*Initialize model structure: */
    69         model=new Model();
     64        _printf_("create finite element model:\n");
     65        femmodel=new FemModel(fid,solution_type,analyses,1);
    7066
    71         _printf_("read and create finite element model:\n");
    72         model->AddFormulation(fid,BalancedthicknessAnalysisEnum);
    73 
    74         /*recover parameters: */
    75         model->FindParam(&waitonlock,WaitOnLockEnum);
    76         model->FindParam(&qmu_analysis,QmuAnalysisEnum);
     67        /*add outputfilename in parameters: */
     68        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename);
     69       
     70        /*get parameters: */
     71        femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
     72        femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
     73        femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
    7774
    7875        MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
    … …  
    8481                _printf_("call computational core:\n");
    8582                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    86                 results=balancedthickness_core(model);
     83                balancedthickness_core(model);
    8784                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
     85
     86                _printf_("write results to disk:\n");
     87                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
     88        }
    8889
    8990        }
    … …  
    9596                #ifdef _HAVE_DAKOTA_
    9697                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    97                 Qmux(model,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
     98                Qmux(femmodel,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
    9899                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    99100                #else
    … …  
    101102                #endif
    102103        }
    103 
    104         _printf_("write results to disk:\n");
    105         OutputResults(results,outputfilename);
    106104
    107105        if (waitonlock>0){
    … …  
    111109
    112110        /*Free ressources:*/
    113         delete results;
    114         delete model;
     111        delete femmodel;
    115112
    116113        /*Get finish time and close*/

  • issm/trunk/src/c/solutions/balancedthickness2.cpp

    r3938 r4055  
    2323        char* inputfilename=NULL;
    2424        char* outputfilename=NULL;
     25        bool  qmu_analysis=false;
    2526        char* lockname=NULL;
    2627        bool  waitonlock=false;
    2728
    28         Model* model=NULL;
    29 
    30         double  dt;
    31         double  yts;
    32         bool    qmu_analysis;
    33 
    34         /*Results: */
    35         Results* results=NULL;
    36         Result*  result=NULL;
    37 
    38         Param*   param=NULL;
     29        /*FemModel: */
     30        FemModel* femmodel=NULL;
    3931
    4032        /*time*/
    … …  
    4234        double   start_core, finish_core;
    4335        double   start_init, finish_init;
     36
     37        int analyses[1]={Balancedthickness2AnalysisEnum};
     38        int solution_type=Balancedthickness2AnalysisEnum;
     39
    4440
    4541        MODULEBOOT();
    … …  
    6864        fid=pfopen(inputfilename,"rb");
    6965
    70         /*Initialize model structure: */
    71         model=new Model();
     66        _printf_("create finite element model:\n");
     67        femmodel=new FemModel(fid,solution_type,analyses,1);
    7268
    73         _printf_("read and create finite element model:\n");
    74         model->AddFormulation(fid,Balancedthickness2AnalysisEnum);
    75 
    76         /*recover parameters: */
    77         model->FindParam(&waitonlock,WaitOnLockEnum);
    78         model->FindParam(&qmu_analysis,QmuAnalysisEnum);
     69        /*add outputfilename in parameters: */
     70        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename);
    7971
    8072        MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
    … …  
    8678                _printf_("call computational core:\n");
    8779                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    88                 results=balancedthickness2_core(model);
     80                balancedthickness2_core(femmodel);
    8981                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
     82
     83                _printf_("write results to disk:\n");
     84                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
    9085
    9186        }
    … …  
    9792                #ifdef _HAVE_DAKOTA_
    9893                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    99                 Qmux(model,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
     94                Qmux(femmodel,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
    10095                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    10196                #else
    … …  
    10398                #endif
    10499        }
    105 
    106         _printf_("write results to disk:\n");
    107         OutputResults(results,outputfilename);
    108100
    109101        if (waitonlock>0){
    … …  
    113105
    114106        /*Free ressources:*/
    115         delete results;
    116         delete model;
     107        delete femmodel;
    117108
    118109        /*Get finish time and close*/

  • issm/trunk/src/c/solutions/balancedthickness2_core.cpp

    r4043 r4055  
    1212#include "../solvers/solvers.h"
    1313
    14 Results* balancedthickness2_core(Model* model){
    15 
    16         extern int my_rank;
    17 
    18         /*output: */
    19         Results* results=NULL;
    20         Result* result=NULL;
    21 
    22         /*intermediary: */
    23         Vec vx_g=NULL;
    24         Vec vy_g=NULL;
    25 
    26         /*solutions: */
    27         Vec h_g=NULL;
     14void balancedthickness2_core(FemModel* femmodel){
    2815
    2916        /*flags: */
    3017        int verbose=0;
    31         int numberofdofspernode;
    32         int numberofnodes;
    33         int numberofvertices;
    34         int dofs[1]={1};
     18        int dim;
    3519
    36         /*fem balancedthickness2 model: */
    37         FemModel* fem_p=NULL;
    38 
    39         //initialize results:
    40         results=new Results();
    41 
    42         fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum);
    43 
    44         //first recover parameters common to all solutions
    45         model->FindParam(&verbose,VerboseEnum);
    46         model->FindParam(&numberofnodes,NumberOfNodesEnum);
    47         model->FindParam(&numberofvertices,NumberOfVerticesEnum);
    48         model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
     20        /*activate formulation: */
     21        femmodel->SetCurrentAnalysis(BalancedthicknessAnalysisEnum);
     22       
     23        /*recover parameters: */
     24        femmodel->parameters->FindParam(&verbose,VerboseEnum);
     25        femmodel->parameters->FindParam(&dim,DimEnum);
    4926
    5027        _printf_("depth averaging velocity...\n");
    51         //vx_g=inputs->Get("vx",&dofs[0],1);
    52         //vy_g=inputs->Get("vy",&dofs[0],1);
    53         /* NOT WORKING YET....
    54         FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx");
    55         FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy");
    56         */
    57         //inputs->Add("vx_average",vx_g,1,numberofvertices);
    58         //inputs->Add("vy_average",vy_g,1,numberofvertices);
    59         ISSMERROR("not supported yet!");
    60        
     28        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
     29        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
     30        if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum);
     31
     32
    6133        _printf_("call computational core:\n");
    62         diagnostic_core_linear(&h_g,fem_p,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
     34        solver_linear(NULL,femmodel);
    6335
    6436        _printf_("Averaging over vertices:\n");
    65         FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
     37        ISSMERROR(" not supported yet!");
     38//FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
    6639
    6740//      _printf_("extrude computed thickness on all layers:\n");
    6841//      FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
    6942
    70         /*Plug results into output dataset: */
    71         InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
    72         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(Balancedthickness2AnalysisEnum)));
    73 
    74         /*Free ressources:*/
    75         VecFree(&vx_g);
    76         VecFree(&vy_g);
    77         VecFree(&h_g);
    78        
    79         /*return: */
    80         return results;
     43        if(verbose)_printf_("saving results:\n");
     44        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    8145
    8246}

  • issm/trunk/src/c/solutions/balancedthickness_core.cpp

    r4043 r4055  
    1212#include "../solvers/solvers.h"
    1313
    14 Results* balancedthickness_core(Model* model){
     14void balancedthickness_core(FemModel* femmodel){
    1515
    16         extern int my_rank;
     16        /*parameters: */
     17        int verbose=0;
     18        int dim;
    1719
    18         /*output: */
    19         Results* results=NULL;
    20         Result* result=NULL;
    21 
    22         /*intermediary: */
    23         Vec u_g=NULL;
    24 
    25         /*solutions: */
    26         Vec h_g=NULL;
    27 
    28         /*flags: */
    29         int verbose=0;
    30         int numberofdofspernode;
    31         int numberofnodes;
    32         int dofs[2]={1,1};
    33 
    34         /*fem balancedthickness model: */
    35         FemModel* fem_p=NULL;
    36 
    37         //initialize results:
    38         results=new Results();
    39 
    40         /*recover fem model: */
    41         fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
    42 
    43         //first recover parameters common to all solutions
    44         model->FindParam(&verbose,VerboseEnum);
    45         model->FindParam(&numberofnodes,NumberOfNodesEnum);
    46         model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
     20        /*activate formulation: */
     21        femmodel->SetCurrentAnalysis(BalancedthicknessAnalysisEnum);
     22       
     23        /*recover parameters: */
     24        femmodel->parameters->FindParam(&verbose,VerboseEnum);
     25        femmodel->parameters->FindParam(&dim,DimEnum);
    4726
    4827        _printf_("depth averaging velocity...\n");
    49         //u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
    50         //FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity");
    51         //inputs->Add("velocity_average",u_g,2,numberofnodes);
    52         ISSMERROR("not supported yet!");
    53        
     28        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
     29        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
     30        if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum);
     31
    5432        _printf_("call computational core:\n");
    55         diagnostic_core_linear(&h_g,fem_p,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
     33        solver_linear(NULL,femmodel);
    5634
    5735        _printf_("extrude computed thickness on all layers:\n");
    58         FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
     36        ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,ThicknessEnum);
    5937
    60         /*Plug results into output dataset: */
    61         InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
    62         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(BalancedthicknessAnalysisEnum)));
    63 
    64         /*Free ressources:*/
    65         VecFree(&u_g);
    66         VecFree(&h_g);
    67        
    68         /*return: */
    69         return results;
     38        if(verbose)_printf_("saving results:\n");
     39        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    7040
    7141}

  • issm/trunk/src/c/solutions/balancedvelocities.cpp

    r3938 r4055  
    2323        char* inputfilename=NULL;
    2424        char* outputfilename=NULL;
     25        bool   qmu_analysis;
    2526        char* lockname=NULL;
    26         int   numberofnodes;
    2727        bool  waitonlock=false;
    2828
    29         Model* model=NULL;
    30 
    31         bool   qmu_analysis;
    32 
    33         /*Results: */
    34         Results* results=NULL;
    35         Result*  result=NULL;
    36         Param*   param=NULL;
     29        /*FemModel: */
     30        FemModel* femmodel=NULL;
    3731
    3832        /*time*/
    … …  
    4034        double   start_core, finish_core;
    4135        double   start_init, finish_init;
     36
     37        int analyses[1]={BalancedvelocitiesAnalysisEnum};
     38        int solution_type=BalancedvelocitiesAnalysisEnum;
    4239
    4340        MODULEBOOT();
    … …  
    6663        fid=pfopen(inputfilename,"rb");
    6764
    68         /*Initialize model structure: */
    69         model=new Model();
     65        _printf_("create finite element model:\n");
     66        femmodel=new FemModel(fid,solution_type,analyses,1);
    7067
    71         _printf_("read and create finite element model:\n");
    72         model->AddFormulation(fid,BalancedvelocitiesAnalysisEnum);
    73 
    74         /*recover parameters: */
    75         model->FindParam(&waitonlock,WaitOnLockEnum);
    76         model->FindParam(&qmu_analysis,QmuAnalysisEnum);
     68        /*add outputfilename in parameters: */
     69        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename);
     70       
     71        /*get parameters: */
     72        femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
     73        femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
     74        femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
    7775
    7876        MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
    … …  
    8482                _printf_("call computational core:\n");
    8583                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    86                 results=balancedvelocities_core(model);
     84                balancedvelocities_core(femmodel);
    8785                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
     86
     87                _printf_("write results to disk:\n");
     88                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
    8889
    8990        }
    … …  
    9596                #ifdef _HAVE_DAKOTA_
    9697                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    97                 Qmux(model,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum);
     98                Qmux(femmodel,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum);
    9899                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    99100                #else
    … …  
    101102                #endif
    102103        }
    103 
    104         _printf_("write results to disk:\n");
    105         OutputResults(results,outputfilename);
    106104
    107105        if (waitonlock>0){
    … …  
    111109
    112110        /*Free ressources:*/
    113         delete results;
    114         delete model;
     111        delete femmodel;
    115112
    116113        /*Get finish time and close*/

  • issm/trunk/src/c/solutions/balancedvelocities_core.cpp

    r4043 r4055  
    1111#include "../solvers/solvers.h"
    1212
    13 Results* balancedvelocities_core(Model* model){
    14 
    15         extern int my_rank;
    16 
    17         /*output: */
    18         Results* results=NULL;
    19         Result* result=NULL;
    20 
    21         /*intermediary: */
    22         Vec u_g=NULL;
    23 
    24         /*solutions: */
    25         Vec v_g=NULL;
     13void balancedvelocities_core(FemModel* femmodel){
    2614
    2715        /*flags: */
    2816        int verbose=0;
    29         int numberofdofspernode;
    30         int numberofnodes;
    31         int dofs[2]={1,1};
     17        int dim;
    3218
    33         /*fem balancedvelocities model: */
    34         FemModel* fem_p=NULL;
    35 
    36         //initialize results:
    37         results=new Results();
    38 
    39         /*recover fem model: */
    40         fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum);
    41 
    42         //first recover parameters common to all solutions
    43         model->FindParam(&verbose,VerboseEnum);
    44         model->FindParam(&numberofnodes,NumberOfNodesEnum);
    45         model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
     19        /*activate formulation: */
     20        femmodel->SetCurrentAnalysis(BalancedvelocitiesAnalysisEnum);
     21       
     22        /*recover parameters: */
     23        femmodel->parameters->FindParam(&verbose,VerboseEnum);
     24        femmodel->parameters->FindParam(&dim,DimEnum);
    4625
    4726        _printf_("depth averaging velocity...\n");
    48         //u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
    49         //FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity");
    50         //inputs->Add("velocity_average",u_g,2,numberofnodes);
    51         ISSMERROR(" not supported yet!");
    52        
     27        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
     28        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
     29        if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum);
     30
    5331        _printf_("call computational core:\n");
    54         diagnostic_core_linear(&v_g,fem_p,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum);
     32        solver_linear(NULL,femmodel);
    5533
    56         _printf_("extrude computed thickness on all layers:\n");
    57         FieldExtrudex( v_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
     34        _printf_("extrude computed velocity on all layers:\n");
     35        ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
     36        ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
    5837
    59         /*Plug results into output dataset: */
    60         InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result);
    61         InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result);
    62         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(BalancedvelocitiesAnalysisEnum)));
     38        if(verbose)_printf_("saving results:\n");
     39        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
     40        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
    6341
    64         /*Free ressources:*/
    65         VecFree(&u_g);
    66         VecFree(&v_g);
    67 
    68         /*return: */
    69         return results;
    7042}

  • issm/trunk/src/c/solutions/bedslope.cpp

    r4030 r4055  
    3737        double   start_init, finish_init;
    3838
    39         int analyses[1]={BedSlopeComputeAnalysisEnum};
    40         int solution_type=BedSlopeComputeSolutionEnum;
     39        int analyses[1]={SlopeAnalysisEnum};
     40        int solution_type=SlopeAnalysisEnum;
    4141
    4242        MODULEBOOT();
    … …  
    6767        _printf_("create finite element model, using analyses types statically defined above:\n");
    6868        femmodel=new FemModel(fid,solution_type,analyses,1);
     69
     70        /*add outputfilename in parameters: */
     71        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    6972       
    7073        /*get parameters: */
    … …  
    8386
    8487                _printf_("write results to disk:\n");
    85                 OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum);
     88                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
    8689        }
    8790        else{
    … …  
    9295                #ifdef _HAVE_DAKOTA_
    9396                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    94                 Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum);
     97                Qmux(femmodel,BedSlopeAnalysisEnum,NoneAnalysisEnum);
    9598                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    9699                #else

  • issm/trunk/src/c/solutions/bedslope_core.cpp

    r4043 r4055  
    2525
    2626        /*Call on core computations: */
    27         solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,BedSlopeXAnalysisEnum);
    28         solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,BedSlopeYAnalysisEnum);
     27        femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,BedSlopeXAnalysisEnum);
     28        solver_linear(NULL,femmodel);
     29        femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,BedSlopeYAnalysisEnum);
     30        solver_linear(NULL,femmodel);
    2931       
    3032        /*extrude inputs if we are in 3D: */
    3133        if (dim==3){
    32                 if(verbose)_printf_("%s\n","extruding slope in 3d...");
    33                 InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum,0);
    34                 InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum,0);
     34                if(verbose)_printf_("%s\n","extruding bed in 3d...");
     35                ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,BedSlopeXEnum);
     36                ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,BedSlopeYEnum);
    3537        }
     38       
     39        if(verbose)_printf_("saving results:\n");
     40        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum);
     41        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum);
     42
    3643}

  • issm/trunk/src/c/solutions/control_core.cpp

    r4052 r4055  
    104104                if (((n+1)%5)==0){
    105105                        _printf_("%s","      saving temporary results...");
    106                         controlrestart(femmodel);
     106                        controlrestart(femmodel,J);
    107107                }
    108108        }

  • issm/trunk/src/c/solutions/controlrestart.cpp

    r4052 r4055  
    77#include "../EnumDefinitions/EnumDefinitions.h"
    88
    9 void controlrestart(FemModel* femmodel){
     9void controlrestart(FemModel* femmodel,double* J){
    1010
    11         char*    outputfilename=NULL;
    12        
     11        int      control_type;
     12        int      nsteps;
     13
    1314        /*retrieve output file name: */
    14         model->FindParam(&outputfilename,OutputFileNameEnum);
     15        femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
     16        femmodel->parameters->FindParam(&nsteps,NStepsEnum);
    1517
    1618        /*we essentially want J and the parameter: */
    1719        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
    18         femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
    19         femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
     20        femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
     21        femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumAsString(control_type),1,0));
    2022
    2123        /*write to disk: */
    22         OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
    23        
    24         /*Free ressources:*/
    25         xfree((void**)&outputfilename);
     24        OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
     25
    2626}

  • issm/trunk/src/c/solutions/diagnostic.cpp

    r4042 r4055  
    3636        double   start_init, finish_init;
    3737
    38         int analyses[5]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeComputeAnalysisEnum};
     38        int analyses[5]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum};
    3939        int solution_type=DiagnosticAnalysisEnum;
    4040
    … …  
    5858        lockname=argv[4];
    5959
    60         /*Initialize femmodel structure: */
    6160        MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
    6261
    … …  
    6665        _printf_("create finite element model:\n");
    6766        femmodel=new FemModel(fid,solution_type,analyses,5);
     67
     68        /*add outputfilename in parameters: */
     69        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    6870       
    6971        /*get parameters: */
    … …  
    9496
    9597                _printf_("write results to disk:\n");
    96                 OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
     98                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
    9799        }
    98100        else{

  • issm/trunk/src/c/solutions/diagnostic_core.cpp

    r4052 r4055  
    2222        bool isstokes=false;
    2323        double stokesreconditioning;
     24        bool conserve_loads=true;
     25        bool modify_loads=true;
    2426
    25         /* recover parameters: {{{1*/
     27        /* recover parameters:*/
    2628        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    2729        femmodel->parameters->FindParam(&dim,DimEnum);
    … …  
    3133        femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    3234        femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
    33         /*}}}*/
    3435
    3536        /*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.*/
    … …  
    4849                       
    4950                if(verbose)_printf_("%s\n"," computing hutter velocities...");
    50                 solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum);
     51                femmodel->SetCurrentAnalysis(DiagnosticHutterAnalysisEnum);
     52                solver_linear(NULL,femmodel);
    5153               
    5254                if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum);
    … …  
    5658               
    5759                if(verbose)_printf_("%s\n"," computing horizontal velocities...");
    58                 solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution
     60                femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
     61                solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,modify_loads);
    5962        }
    6063       
    … …  
    6366
    6467                if(verbose)_printf_("%s\n"," computing vertical velocities...");
    65                 solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum);
     68                femmodel->SetCurrentAnalysis(DiagnosticVertAnalysisEnum);
     69                solver_linear(NULL,femmodel);
    6670
    6771                if (isstokes){
    … …  
    7579
    7680                        if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
    77                         solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum);
     81                        femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum);
     82                        solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
    7883                }
    7984        }

  • issm/trunk/src/c/solutions/gradient_core.cpp

    r4048 r4055  
    3232       
    3333        if(verbose)_printf_("%s\n","      compute gradient:");
    34         Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type,AdjointEnum);
     34        Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type);
    3535       
    3636        if(verbose)_printf_("%s\n","      retrieve old gradient:");
    3737        GetVectorFromInputsx(&old_gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, OldGradientEnum,VertexEnum);
    3838
    39         if(control_steady)diagnostic_core(model);
     39        if(control_steady)diagnostic_core(femmodel);
    4040       
    4141        if (step>0 && search_scalar==0){
    … …  
    5555
    5656        /*plug back into inputs: */
    57         UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmoel->parameters,gradient,GradientEnum,VertexEnum);
    58         UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmoel->parameters,old_gradient,OldGradientEnum,VertexEnum);
     57        UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,gradient,GradientEnum,VertexEnum);
     58        UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,old_gradient,OldGradientEnum,VertexEnum);
    5959
    6060        /*Free ressources and return:*/

  • issm/trunk/src/c/solutions/objectivefunctionC.cpp

    r4052 r4055  
    33 */
    44
    5 #include "../modules/modules.h"
    6 #include "./solutions.h"
    7 #include "../EnumDefinitions/EnumDefinitions.h"
    8 
     5/*include files: {{{1*/
    96#ifdef HAVE_CONFIG_H
    107        #include "config.h"
    … …  
    129#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    1310#endif
     11
     12#include "../toolkits/toolkits.h"
     13#include "../objects/objects.h"
     14#include "../shared/shared.h"
     15#include "../EnumDefinitions/EnumDefinitions.h"
     16#include "../solvers/solvers.h"
     17#include "./solutions.h"
     18#include "../modules/modules.h"
     19#include "../include/include.h"
     20/*}}}*/
    1421
    1522double objectivefunctionC(double search_scalar,OptArgs* optargs){
    … …  
    2229        /*parameters: */
    2330        FemModel* femmodel=NULL;
    24         int n;
     31        int     n;
    2532        double* optscal=NULL;
    2633        double* fit=NULL;
    … …  
    2936        int     control_type;
    3037        int     analysis_type;
    31         int     sub_analysis_type;
    3238        int     control_steady;
    33                
     39        bool    conserve_loads=true;
     40       
    3441        /*Recover finite element model: */
    3542        femmodel=optargs->femmodel;
    … …  
    4552        femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
    4653        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    47         femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum);
    4854
    4955        /*Use ControlParameterEnum input to  reinitialize our input parameter: */
    … …  
    5763
    5864        /*Run diagnostic with updated inputs: */
    59         if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,true,femmodel,DiagnosticAnalysisEnum,sub_analysis_type);  //true means we conserve loads at each diagnostic run
     65        if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);  //true means we conserve loads at each diagnostic run
    6066        else                diagnostic_core(femmodel);  //We need a 3D velocity!! (vz is required for the next thermal run)
    6167
    6268        /*Compute misfit for this velocity field.*/
    6369        UpdateInputsFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,fit[n],FitEnum);
    64         CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,analysis_type,sub_analysis_type);
     70        CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
    6571
    6672        /*Free ressources:*/

  • issm/trunk/src/c/solutions/prognostic.cpp

    r3938 r4055  
    2424        char* outputfilename=NULL;
    2525        char* lockname=NULL;
    26         int   numberofnodes;
     26        bool  qmu_analysis=false;
    2727        bool  waitonlock=false;
    2828
    29         Model* model=NULL;
    30 
    31         bool qmu_analysis;
    32 
    33         /*Results: */
    34         Results* results=NULL;
    35 
    36         Param*   param=NULL;
     29        /*FemModel: */
     30        FemModel* femmodel=NULL;
    3731
    3832        /*time*/
    … …  
    4034        double   start_core, finish_core;
    4135        double   start_init, finish_init;
     36
     37        int analyses[1]={PrognosticAnalysisEnum};
     38        int solution_type=PrognosticAnalysisEnum;
    4239
    4340        MODULEBOOT();
    … …  
    6259        /*Initialize model structure: */
    6360        MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
    64         model=new Model();
    6561
    6662        /*Open handle to data on disk: */
    6763        fid=pfopen(inputfilename,"rb");
    6864
    69         _printf_("read and create finite element model:\n");
    70         model->AddFormulation(fid,PrognosticAnalysisEnum);
     65        _printf_("create finite element model:\n");
     66        femmodel=new FemModel(fid,solution_type,analyses,1);
    7167
    72         /*recover parameters: */
    73         model->FindParam(&waitonlock,WaitOnLockEnum);
    74         model->FindParam(&qmu_analysis,QmuAnalysisEnum);
     68        /*add outputfilename in parameters: */
     69        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
     70
     71        /*get parameters: */
     72        femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
     73        femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
     74
    7575
    7676        MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
    … …  
    8282                _printf_("call computational core:\n");
    8383                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    84                 results=prognostic_core(model);
     84                prognostic_core(femmodel);
    8585                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
     86               
     87                _printf_("write results to disk:\n");
     88                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
    8689
    8790        }
    … …  
    9396                #ifdef _HAVE_DAKOTA_
    9497                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    95                 Qmux(model,PrognosticAnalysisEnum,NoneAnalysisEnum);
     98                Qmux(femmodel,PrognosticAnalysisEnum,NoneAnalysisEnum);
    9699                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    97100                #else
    … …  
    101104
    102105               
    103         _printf_("write results to disk:\n");
    104         OutputResults(results,outputfilename);
    105 
    106106        if (waitonlock>0){
    107107                _printf_("write lock file:\n");
    … …  
    110110
    111111        /*Free ressources:*/
    112         delete results;
    113         delete model;
     112        delete femmodel;
    114113
    115114        /*Get finish time and close*/

  • issm/trunk/src/c/solutions/prognostic2.cpp

    r3938 r4055  
    2424        char* outputfilename=NULL;
    2525        char* lockname=NULL;
    26         int   numberofnodes;
     26        bool  qmu_analysis=false;
    2727        bool  waitonlock=false;
    2828
    29         Model* model=NULL;
    30 
    31         bool   qmu_analysis;
    32 
    33         /*Results: */
    34         Results* results=NULL;
    35 
    36         Param*   param=NULL;
     29        /*FemModel: */
     30        FemModel* femmodel=NULL;
    3731
    3832        /*time*/
    … …  
    4034        double   start_core, finish_core;
    4135        double   start_init, finish_init;
     36
     37        int analyses[1]={Prognostic2AnalysisEnum};
     38        int solution_type=Prognosti2cAnalysisEnum;
    4239
    4340        MODULEBOOT();
    … …  
    6259        /*Initialize model structure: */
    6360        MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
    64         model=new Model();
    6561
    6662        /*Open handle to data on disk: */
    6763        fid=pfopen(inputfilename,"rb");
    6864
    69         _printf_("read and create finite element model:\n");
    70         model->AddFormulation(fid,Prognostic2AnalysisEnum);
     65        _printf_("create finite element model:\n");
     66        femmodel=new FemModel(fid,solution_type,analyses,1);
     67
     68        /*add outputfilename in parameters: */
     69        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    7170
    7271        /*recover parameters: */
    7372        model->FindParam(&waitonlock,WaitOnLockEnum);
    7473        model->FindParam(&qmu_analysis,QmuAnalysisEnum);
    75         model->FindParam(&numberofnodes,NumberOfNodesEnum);
     74
    7675
    7776        MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
    … …  
    8382                _printf_("call computational core:\n");
    8483                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    85                 results=prognostic2_core(model);
     84                prognostic2_core(femmodel);
    8685                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
     86
     87                _printf_("write results to disk:\n");
     88                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
    8789
    8890        }
    … …  
    9496                #ifdef _HAVE_DAKOTA_
    9597                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    96                 Qmux(model,Prognostic2AnalysisEnum,NoneAnalysisEnum);
     98                Qmux(femmodel,Prognostic2AnalysisEnum,NoneAnalysisEnum);
    9799                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    98100                #else
    … …  
    100102                #endif
    101103        }
    102 
    103         _printf_("write results to disk:\n");
    104         OutputResults(results,outputfilename);
    105104
    106105        if (waitonlock>0){
    … …  
    110109
    111110        /*Free ressources:*/
    112         delete results;
    113         delete model;
     111        delete femmodel;
    114112
    115113        /*Get finish time and close*/

  • issm/trunk/src/c/solutions/prognostic2_core.cpp

    r4043 r4055  
    1111#include "../solvers/solvers.h"
    1212
    13 Results* prognostic2_core(Model* model){
    14 
    15         extern int my_rank;
    16 
    17         /*output: */
    18         Results* results=NULL;
    19         Result* result=NULL;
    20 
    21         /*intermediary: */
    22         Vec vx_g=NULL;
    23         Vec vy_g=NULL;
    24 
    25         /*solutions: */
    26         Vec h_g=NULL;
     13void prognostic2_core(FemModel* femmodel){
    2714
    2815        /*flags: */
    2916        int verbose=0;
    30         int numberofdofspernode;
    31         int numberofnodes;
    32         int numberofvertices;
    33         int dofs[1]={1};
    3417
    35         /*fem prognostic model: */
    36         FemModel* fem_p=NULL;
     18        /*activate formulation: */
     19        femmodel->SetCurrentAnalysis(Prognostic2AnalysisEnum);
    3720
    38         //initialize results:
    39         results=new Results();
    40 
    41         /*recover fem model: */
    42         fem_p=model->GetFormulation(Prognostic2AnalysisEnum);
    43 
    44         //first recover parameters common to all solutions
    45         model->FindParam(&verbose,VerboseEnum);
    46         model->FindParam(&numberofnodes,NumberOfNodesEnum);
    47         model->FindParam(&numberofvertices,NumberOfVerticesEnum);
    48         model->FindParam(&numberofdofspernode,NumberOfNodesEnum);
     21        /*recover parameters: */
     22        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    4923
    5024        _printf_("depth averaging velocity...\n");
    5125        /*Where is it done?*/
     26
     27        _printf_("call computational core:\n");
     28        solver_linear(NULL,femmodel);
    5229       
    53         _printf_("call computational core:\n");
    54         diagnostic_core_linear(&h_g,fem_p,Prognostic2AnalysisEnum,NoneAnalysisEnum);
    55 
    5630        _printf_("Averaging over vertices:\n");
    57         FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
     31        //FieldAverageOntoVerticesx(&h_g,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    5832
    5933        //_printf_("extrude computed thickness on all layers:\n");
    6034        //FieldExtrudex(h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
    6135
    62         /*Plug results into output dataset: */
    63         InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
    64         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(Prognostic2AnalysisEnum)));
    65 
    66         /*Free ressources:*/
    67         VecFree(&vx_g);
    68         VecFree(&vy_g);
    69         VecFree(&h_g);
    70        
    71         /*return: */
    72         return results;
     36        if(verbose)_printf_("saving results:\n");
     37        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    7338
    7439}

  • issm/trunk/src/c/solutions/prognostic_core.cpp

    r4043 r4055  
    1111#include "../solvers/solvers.h"
    1212
    13 Results* prognostic_core(Model* model){
     13void prognostic_core(FemModel* femmodel){
    1414
    15         extern int my_rank;
    16 
    17         /*output: */
    18         Results* results=NULL;
    19         Result*  result=NULL;
    20 
    21         /*solutions: */
    22         Vec h_g=NULL;
    23 
    24         /*flags: */
     15        /*parameters: */
    2516        int verbose=0;
    2617
    27         /*fem prognostic model: */
    28         FemModel* fem_p=NULL;
    29 
    30         //initialize results
    31         results=new Results();
    32 
    33         /*recover fem model: */
    34         fem_p=model->GetFormulation(PrognosticAnalysisEnum); ISSMASSERT(fem_p);
    35 
    36         //first recover parameters common to all solutions
    37         model->FindParam(&verbose,VerboseEnum);
     18        /*activate formulation: */
     19        femmodel->SetCurrentAnalysis(PrognosticAnalysisEnum);
     20       
     21        /*recover parameters: */
     22        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    3823
    3924        _printf_("depth averaging velocity...\n");
    40         DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VxEnum);
    41         DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VyEnum);
     25        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum,VxAverageEnum);
     26        DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum,VyAverageEnum);
    4227       
    4328        _printf_("call computational core:\n");
    44         diagnostic_core_linear(&h_g,fem_p,PrognosticAnalysisEnum,NoneAnalysisEnum);
     29        solver_linear(NULL,femmodel);
    4530               
    46         _printf_("update inputs:\n");
    47         UpdateInputsFromSolutionx( fem_p->elements,fem_p->nodes, fem_p->vertices, fem_p->loads, fem_p->materials, fem_p->parameters,h_g,PrognosticAnalysisEnum,NoneAnalysisEnum);
     31        _printf_("extrude computed thickness on all layers:\n");
     32        ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,ThicknessEnum);
    4833
    49         _printf_("extrude computed thickness on all layers:\n");
    50         ExtrudeInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum);
     34        if(verbose)_printf_("saving results:\n");
     35        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    5136       
    52         _printf_("extract result from extruded inputs: \n");
    53         InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
    54         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(PrognosticAnalysisEnum)));
     37}
    5538
    56         /*Free ressources:*/
    57         VecFree(&h_g);
    58        
    59         //return:
    60         return results;
    61 
    62 }

  • issm/trunk/src/c/solutions/solutions.h

    r4052 r4055  
    1313
    1414/*cores: */
    15 Results* prognostic_core(FemModel* femmodel);
    16 Results* prognostic2_core(FemModel* femmodel);
    17 Results* balancedthickness_core(FemModel* femmodel);
    18 Results* balancedthickness2_core(FemModel* femmodel);
    19 Results* balancedvelocities_core(FemModel* femmodel);
    20 Results* slopecompute_core(FemModel* femmodel);
    21 Results* steadystate_core(FemModel* femmodel);
    22 Results* transient_core(FemModel* femmodel);
    23 Results* transient_core_2d(FemModel* femmodel);
    24 Results* transient_core_3d(FemModel* femmodel);
    2515void adjoint_core(FemModel* femmodel);
    2616void gradient_core(FemModel* femmodel,int step=0, double search_scalar=0);
    2717void diagnostic_core(FemModel* femmodel);
    2818void thermal_core(FemModel* femmodel);
    29 void surfaceslope_core(FemModel* femfemmodel);
    30 void bedslope_core(FemModel* femfemmodel);
     19void thermal_core_step(FemModel* femmodel,int step, double time);
     20void surfaceslope_core(FemModel* femmodel);
     21void bedslope_core(FemModel* femmodel);
    3122void control_core(FemModel* femmodel);
     23void prognostic_core(FemModel* femmodel);
     24void prognostic2_core(FemModel* femmodel);
     25void balancedthickness_core(FemModel* femmodel);
     26void balancedthickness2_core(FemModel* femmodel);
     27void balancedvelocities_core(FemModel* femmodel);
     28void slopecompute_core(FemModel* femmodel);
     29void steadystate_core(FemModel* femmodel);
     30void transient2d_core(FemModel* femmodel);
     31void transient3d_core(FemModel* femmodel);
     32double objectivefunctionC(double search_scalar,OptArgs* optargs);
    3233
    33 //int GradJOrth(WorkspaceParams* workspaceparams);
     34//convergence:
    3435void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters);
    3536int controlconvergence(double* J, double* fit, double eps_cm, int n);
    3637
    37 int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femfemmodel);
     38//optimization
     39int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);
     40int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);
     41int GradJSearch(double* search_vector,FemModel* femmodel,int step);
    3842
    39 int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femfemmodel);
    40        
    41 double objectivefunctionC(double search_scalar,OptArgs* optargs);
    42 
    43 int GradJSearch(double* search_vector,FemModel* femfemmodel,int step);
    44 //int GradJCheck(WorkspaceParams* workspaceparams,int step,int status);
    45 
    46 //int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
     43//diverse
    4744void WriteLockFile(char* filename);
    48 
    49 void stokescontrolinit(FemModel* femfemmodel);
    50 void controlrestart(FemModel* femfemmodel);
    51 
    52 void CreateFemModel(FemModel* femfemmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);
    53 //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femfemmodel,char* filename);
    54 void ResetBoundaryConditions(FemModel* femfemmodel, int analysis_type, int sub_analysis_type);
     45void stokescontrolinit(FemModel* femmodel);
     46void controlrestart(FemModel* femmodel,double* J);
     47void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type);
    5548
    5649#endif

  • issm/trunk/src/c/solutions/steadystate.cpp

    r3938 r4055  
    2424        char* outputfilename=NULL;
    2525        char* lockname=NULL;
    26         int   numberofnodes;
    2726        bool  qmu_analysis=false;
    2827        bool  control_analysis=false;
    29         char* control_type=NULL;
     28        bool waitonlock=false;
    3029
    31         /*Model: */
    32         Model* model=NULL;
    33         FemModel* fem_dh=NULL;
    34         FemModel* fem_ds=NULL;
    35 
    36         /*Results: */
    37         Results* results=NULL;
    38        
    39         bool waitonlock=false;
    40        
    41         double* u_g_initial=NULL;
    42         double* p_g_initial=NULL;
    43         double* u_g_obs=NULL;
    44         double* weights=NULL;
    45         double  dt;
    46         BoolParam*  param=NULL;
     30        /*FemModel: */
     31        FemModel* femmodel=NULL;
    4732
    4833        /*time*/
    … …  
    5035        double   start_core, finish_core;
    5136        double   start_init, finish_init;
     37
     38        /*intermediary: */
     39        BoolParam* param=NULL;
     40
     41        int analyses[7]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum};
     42        int solution_type=SteadystateAnalysisEnum;
    5243
    5344        MODULEBOOT();
    … …  
    7061        lockname=argv[4];
    7162
     63        /*Initialize model structure: */
     64        MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
     65
    7266        /*Open handle to data on disk: */
    7367        fid=pfopen(inputfilename,"rb");
    7468
    75         /*Initialize model structure: */
    76         MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
    77         model=new Model();
     69        _printf_("create finite element model:\n");
     70        femmodel=new FemModel(fid,solution_type,analyses,5);
    7871
    79         _printf_("read and create finite element model:\n");
    80        
    81         _printf_("\n   reading diagnostic horiz model data:\n");
    82         model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum);
    83 
    84         _printf_("\n   reading diagnostic vert model data:\n");
    85         model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum);
    86        
    87         _printf_("\n   reading diagnostic stokes model data:\n");
    88         model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum);
    89        
    90         _printf_("\n   reading diagnostic hutter model data:\n");
    91         model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum);
    92        
    93         _printf_("\n   reading surface and bed slope computation model data:\n");
    94         model->AddFormulation(fid,SlopecomputeAnalysisEnum);
    95 
    96         _printf_("\n   read and create thermal finite element model:\n");
    97         model->AddFormulation(fid,ThermalAnalysisEnum);
    98         _printf_("\n   read and create melting finite element model:\n");
    99         model->AddFormulation(fid,MeltingAnalysisEnum);
     72        /*add outputfilename in parameters: */
     73        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    10074
    10175        /*recover parameters: */
    … …  
    11387                        _printf_("call computational core:\n");
    11488                        MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    115                         results=steadystate_core(model);
     89                        steadystate_core(femmodel);
    11690                        MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    11791
    … …  
    11993                else{
    12094                        /*change control_steady to 1 to know we are doing steadystate*/
    121                         fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
    122                         fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
    123 
    124                         param=(BoolParam*)fem_dh->parameters->FindParamObject(ControlSteadyEnum); param->value=true;
    125                         param=(BoolParam*)fem_ds->parameters->FindParamObject(ControlSteadyEnum); param->value=true;
     95                        param=(BoolParam*)femmodel->parameters->FindParamObject(ControlSteadyEnum); param->value=true;
    12696
    12797                        /*run control analysis: */
    12898                        _printf_("call computational core:\n");
    12999                        MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    130                         results=control_core(model);
     100                        control_core(femmodel);
    131101                        MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    132 
    133102                }
    134103
    135104                _printf_("write results to disk:\n");
    136                 OutputResults(results,outputfilename);
     105                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
     106
    137107        }
    138108        else{
    … …  
    142112                #ifdef _HAVE_DAKOTA_
    143113                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    144                 Qmux(model,SteadystateAnalysisEnum,NoneAnalysisEnum);
     114                Qmux(femmodel,SteadystateAnalysisEnum,NoneAnalysisEnum);
    145115                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    146116                #else
    … …  
    155125       
    156126        /*Free ressources */
    157         xfree((void**)&u_g_initial);
    158         xfree((void**)&u_g_obs);
    159         xfree((void**)&p_g_initial);
    160         delete model;
    161         delete results;
     127        delete femmodel;
    162128
    163129        /*Get finish time and close*/

  • issm/trunk/src/c/solutions/steadystate_core.cpp

    r4043 r4055  
    99#include "./solutions.h"
    1010#include "../modules/modules.h"
     11#include "../include/include.h"
    1112#include "../solvers/solvers.h"
    1213
    13 Results* steadystate_core(Model* model){
     14void steadystate_core(FemModel* femmodel){
    1415
    15         extern int my_rank;
     16        /*intermediary: */
     17        int step;
    1618
    17         /*fem models: */
    18         FemModel* fem_dh=NULL;
    19         FemModel* fem_dv=NULL;
    20         FemModel* fem_dhu=NULL;
    21         FemModel* fem_ds=NULL;
    22         FemModel* fem_sl=NULL;
    23         FemModel* fem_t=NULL;
    24         FemModel* fem_m=NULL;
     19        /*parameters: */
     20        int verbose;
     21        int dim;
     22       
     23        /* recover parameters:*/
     24        femmodel->parameters->FindParam(&verbose,VerboseEnum);
     25        femmodel->parameters->FindParam(&dim,DimEnum);
    2526
    26         /*output: */
    27         Results* results=NULL;
    28         Results* results_thermal=NULL;
    29         Results* results_diagnostic=NULL;
    30 
    31         /*solutions: */
    32         Vec u_g=NULL;
    33         Vec old_u_g=NULL;
    34         Vec t_g=NULL;
    35         Vec t_g_average=NULL;
    36         Vec old_t_g=NULL;
    37         Vec p_g=NULL;
    38         Vec m_g=NULL;
    39         Vec du_g=NULL;
    40         Vec dt_g=NULL;
    41         double ndu,nu;
    42         double normdt,normt;
    43         double eps_rel;
    44 
    45         /*flags: */
    46         int verbose=0;
    47         int isstokes=0;
    48         int numberofnodes;
    49         int ndof;
    50         int converged;
    51         int step;
    52 
    53         /*recover fem models: */
    54         fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
    55         fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum);
    56         fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
    57         fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum);
    58         fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum);
    59         fem_t=model->GetFormulation(ThermalAnalysisEnum);
    60         fem_m=model->GetFormulation(MeltingAnalysisEnum);
    61 
    62 
    63         //first recover parameters common to all solutions
    64         model->FindParam(&verbose,VerboseEnum);
    65         model->FindParam(&numberofnodes,NumberOfNodesEnum);
    66         model->FindParam(&eps_rel,EpsResEnum);
    67         model->FindParam(&isstokes,IsStokesEnum);
    68 
    69         //initialize:
    70         converged=0;
     27        /*intialize counters: */
    7128        step=1;
    72 
    73         if (isstokes)ndof=4;
    74         else ndof=3;
    7529
    7630        for(;;){
    7731       
    7832                if(verbose)_printf_("%s%i\n","   computing temperature and velocity for step: ",step);
     33                thermal_core(femmodel);
    7934
    80                 //first compute temperature at steady state.
    81                 results_thermal=thermal_core(model);
     35                if(verbose)_printf_("%s\n","computing depth average temperature");
     36                DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum);
    8237       
    83                 //get t_g and m_g;
    84                 VecFree(&t_g);results_thermal->FindResult(&t_g,"t_g");
    85                 VecFree(&m_g);results_thermal->FindResult(&m_g,"m_g");
    86                 delete results_thermal;
     38                if(verbose)_printf_("%s\n","computing new velocity");
     39                diagnostic_core(femmodel);
    8740
    88                 //Add temperature to inputs.
    89                 //compute depth averaged temperature and add to inputs
    90                 VecDuplicatePatch(&t_g_average,t_g);
    91                 FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->vertices,fem_t->loads, fem_t->materials,fem_t->parameters,"temperature");
    92                 model->UpdateInputsFromVector(t_g_average,TemperatureAverageEnum,VertexEnum);
    93                 model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum);
    94                 VecFree(&t_g_average); //not needed anymore
    95 
    96                 //now compute diagnostic velocity using the steady state temperature.
    97                 results_diagnostic=diagnostic_core(model);
    98 
    99                 //get p_g and u_g
    100                 VecFree(&u_g);results_diagnostic->FindResult(&u_g,"u_g");
    101                 VecFree(&p_g);results_diagnostic->FindResult(&p_g,"p_g");
    102                 delete results_diagnostic;
    103 
    104                 //convergence?
    105                 if(step>1){
    106                         VecDuplicatePatch(&du_g,old_u_g);VecAYPX(du_g,-1.0,u_g);
    107                         VecNorm(du_g,NORM_2,&ndu); VecNorm(old_u_g,NORM_2,&nu); VecFree(&du_g);
    108 
    109                         VecDuplicatePatch(&dt_g,old_t_g); VecAYPX(dt_g,-1.0,t_g);
    110                         VecNorm(dt_g,NORM_2,&normdt); VecNorm(old_t_g,NORM_2,&normt);VecFree(&dt_g);
    111                                        
    112                         if (verbose) _printf_("%-60s%g\n                                     %s%g\n                                     %s%g%s\n",
    113                                           "      relative convergence criterion: velocity -> norm(du)/norm(u)=   ",ndu/nu*100," temperature -> norm(dt)/norm(t)=",normdt/normt*100," eps_rel:                        ",eps_rel*100," %");
     41                if(verbose)_printf_("%s\n","checking velocity, temperature and pressure convergence");
     42                if (step>1) if(steadystateconvergence(femmodel))goto cleanup_and_return;
    11443               
    115                         if ((ndu/nu<=eps_rel)  && (normdt/normt<=eps_rel)) converged=1;
    116                         else converged=0;
    117                 }
    118                 else{
    119                         converged=0;
    120                 }
    121 
    122                 VecFree(&old_u_g);VecDuplicatePatch(&old_u_g,u_g);
    123                 VecFree(&old_t_g);VecDuplicatePatch(&old_t_g,t_g);
    124 
     44                if(verbose)_printf_("%s\n","saving velocity, temperature and pressure to check for convergence at next step");
     45                DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,OldVxEnum);
     46                DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,OldVyEnum);
     47                if(dim==3)DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,OldVzEnum);
     48                DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,OldPressureEnum);
     49                DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,OldTemperatureEnum);
     50               
     51                //increase counter
    12552                step++;
    126                 if (converged)break;
    12753        }
    12854
    129         /*Plug results into output dataset: */
    130         results->AddObject(new Result(results->Size()+1,0,1,"u_g",u_g));
    131         results->AddObject(new Result(results->Size()+1,0,1,"p_g",p_g));
    132         results->AddObject(new Result(results->Size()+1,0,1,"t_g",t_g));
    133         results->AddObject(new Result(results->Size()+1,0,1,"m_g",m_g));
    134         results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(SteadyAnalysisEnum)));
    135 
    136 
    137         /*Free ressource*/
    138         VecFree(&old_u_g);
    139         VecFree(&old_t_g);
    140         VecFree(&u_g);
    141         VecFree(&p_g);
    142         VecFree(&t_g);
    143         VecFree(&m_g);
     55        cleanup_and_return:
     56       
     57        if(verbose)_printf_("saving results:\n");
     58        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
     59        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
     60        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum);
     61        if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
    14462}

  • issm/trunk/src/c/solutions/stokescontrolinit.cpp

    r4052 r4055  
    1818        int isstokes=0;
    1919        double stokesreconditioning;
    20 
     20        bool conserve_loads=true;
     21       
    2122        /*first recover parameters common to all solutions:*/
    22         model->FindParam(&verbose,VerboseEnum);
    23         model->FindParam(&isstokes,IsStokesEnum);
    24         model->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     23        femmodel->parameters->FindParam(&verbose,VerboseEnum);
     24        femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
     25        femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    2526
    2627        /*if no Stokes analysis carried out, assign output and return*/
    27         if (!isstokes)femmodel->SetAnalysisType(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum);
     28        if (!isstokes)femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
    2829
    2930        /* For Stokes inverse control method, we are going to carry out the inversion only on the Stokes part. So we need
    30          * to solve the Hutter or MacAyeal/Pattyn model here, and constrain the Stokes model using the Hutter
     31         * to solve the Hutter or MacAyeal/Pattyn femmodel here, and constrain the Stokes femmodel using the Hutter
    3132         * or MacAyeal/Pattyn at the boundary. We don't want to have to do that at every inversion step, as
    3233         * it needs be done only once: */
    … …  
    3637       
    3738        /*Run a complete diagnostic to update the Stokes spcs: */
    38         solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution
     39        femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
     40        solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
    3941
    4042        //vertical velocity
    41         solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum);
     43        solver_linear(NULL,femmodel);
    4244
    4345        //recondition" pressure computed previously:
    … …  
    4951
    5052        if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
    51         solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum);
    52 
    53         /*Assign output*/
    54         model->SetActiveFormulation(fem_ds);
     53        femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum);
     54        solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
    5555}

  • issm/trunk/src/c/solutions/surfaceslope.cpp

    r4030 r4055  
    3737        double   start_init, finish_init;
    3838
    39         int analyses[1]={SurfaceSlopeComputeAnalysisEnum};
    40         int solution_type=SurfaceSlopeComputeSolutionEnum;
     39        int analyses[1]={SlopeAnalysisEnum};
     40        int solution_type=SlopeAnalysisEnum;
    4141
    4242        MODULEBOOT();
    … …  
    6767        _printf_("create finite element model, using analyses types statically defined above:\n");
    6868        femmodel=new FemModel(fid,solution_type,analyses,1);
     69
     70        /*add outputfilename in parameters: */
     71        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    6972       
    7073        /*get parameters: */
    … …  
    8386
    8487                _printf_("write results to disk:\n");
    85                 OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum);
     88                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
    8689        }
    8790        else{
    … …  
    9295                #ifdef _HAVE_DAKOTA_
    9396                MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
    94                 Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum);
     97                Qmux(femmodel,SurfaceSlopeAnalysisEnum,NoneAnalysisEnum);
    9598                MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
    9699                #else

  • issm/trunk/src/c/solutions/surfaceslope_core.cpp

    r4043 r4055  
    2525
    2626        /*Call on core computations: */
    27         solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeXAnalysisEnum);
    28         solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeYAnalysisEnum);
     27        femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum);
     28        solver_linear(NULL,femmodel);
     29        femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum);
     30        solver_linear(NULL,femmodel);
    2931       
    3032        /*extrude inputs if we are in 3D: */
    3133        if (dim==3){
    3234                if(verbose)_printf_("%s\n","extruding slope in 3d...");
    33                 InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum,0);
    34                 InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum,0);
     35                ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum);
     36                ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum);
    3537        }
     38       
     39        if(verbose)_printf_("saving results:\n");
     40        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum);
     41        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum);
     42
    3643}

  • issm/trunk/src/c/solutions/thermal.cpp

    r4042 r4055  
    6868        femmodel=new FemModel(fid,solution_type,analyses,2);
    6969
     70        /*add outputfilename in parameters: */
     71        femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
     72
    7073        /*get parameters: */
    7174        femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
    … …  
    8588               
    8689                _printf_("write results to disk:\n");
    87                 OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
     90                OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
    8891
    8992        }

  • issm/trunk/src/c/solutions/thermal_core.cpp

    r4043 r4055  
    4242                time=(i+1)*dt;
    4343
    44                 if(verbose)_printf_("computing temperatures:\n");
    45                 solver_thermal_nonlinear(NULL,NULL,femmodel,ThermalAnalysisEnum,NoneAnalysisEnum);
    46 
    47                 if(verbose)_printf_("computing melting:\n");
    48                 solver_linear(NULL,femmodel,MeltingAnalysisEnum,NoneAnalysisEnum);
     44                /*call thermal_core_step: */
     45                thermal_core_step(femmodel,i,time);
    4946
    5047                if(verbose)_printf_("saving results:\n");

  • issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp

    r4029 r4055  
    99#include "./solutions.h"
    1010
    11 void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads, int analysis_type,int sub_analysis_type){
     11void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads){
    1212
    1313
    … …  
    6161
    6262        /*Start non-linear iteration using input velocity: */
    63         GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters, analysis_type, sub_analysis_type);
     63        GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters);
    6464        Reducevectorgtofx(&uf, ug, fem->nodesets);
    6565
    … …  
    7272                if (verbose) _printf_("   Generating matrices\n");
    7373                //*Generate system matrices
    74                 SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     74                SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag);
    7575
    7676                if (verbose) _printf_("   Generating penalty matrices\n");
    7777                //*Generate penalty system matrices
    78                 PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     78                PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag);
    7979
    8080                if (verbose) _printf_("   reducing matrix from g to f set\n");
    … …  
    102102
    103103                //Update inputs using new solution:
    104                 UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type);
     104                UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug);
    105105
    106106                //Deal with penalty loads
    107107                if (verbose) _printf_("   penalty constraints\n");
    108                 PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type);
     108                PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters);
    109109
    110110                if(verbose)_printf_("   number of unstable constraints: %i\n",num_unstable_constraints);
    … …  
    146146                kflag=1; pflag=0; //stiffness generation only
    147147       
    148                 SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     148                SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag);
    149149                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
    150150                MatFree(&Kgg);VecFree(&pg);

  • issm/trunk/src/c/solvers/solver_linear.cpp

    r4029 r4055  
    88#include "../modules/modules.h"
    99
    10 void solver_linear(Vec* pug, FemModel* fem,int analysis_type,int sub_analysis_type){
     10void solver_linear(Vec* pug, FemModel* fem){
    1111
    1212        /*parameters:*/
    … …  
    3333        //*Generate system matrices
    3434        if (verbose) _printf_("   Generating matrices\n");
    35         SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     35        SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
    3636
    3737        if (verbose) _printf_("   Generating penalty matrices\n");
    3838        //*Generate penalty system matrices
    39         PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     39        PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
    4040
    4141        /*!Reduce matrix from g to f size:*/
    … …  
    5656
    5757        //Update inputs using new solution:
    58         UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type);
     58        UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug);
    5959
    6060        /*free ressources: */

  • issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp

    r4029 r4055  
    88#include "../modules/modules.h"
    99
    10 void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type){
     10void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem){
    1111
    1212        /*solution : */
    … …  
    6363                        /*Compute Kgg_nopenalty and pg_nopenalty once for all: */
    6464                        if (count==1){
    65                                 SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     65                                SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
    6666                        }
    6767
    … …  
    7171
    7272                        //apply penalties each time
    73                         PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     73                        PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
    7474                }
    7575                else{
    76                         SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     76                        SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
    7777                        //apply penalties
    78                         PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type);
     78                        PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
    7979                }
    8080
    … …  
    108108                //Deal with penalty loads
    109109                if (verbose) _printf_("   penalty constraints\n");
    110                 PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type);
     110                PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters);
    111111               
    112112                //Update inputs using new solution:
    113113                UpdateInputsFromVectorx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,TemperatureEnum,VertexEnum);
    114                 UpdateInputsFromSolutionx(fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,analysis_type, sub_analysis_type);
     114                UpdateInputsFromSolutionx(fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg);
    115115
    116116                if (!converged){

  • issm/trunk/src/c/solvers/solvers.h

    r4047 r4055  
    1212class FemModel;
    1313
    14 void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type);
    15 void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type);
    16 void solver_linear(Vec* pug, FemModel* femmodel,int  analysis_type,int sub_analysis_type);
    17 void solver_adjoint(Vec* pug, FemModel* femmodel,int  analysis_type,int sub_analysis_type);
     14void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* femmodel);
     15void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* femmodel,bool conserve_loads);
     16void solver_linear(Vec* pug, FemModel* femmodel);
     17void solver_adjoint(Vec* pug, FemModel* femmodel);
    1818
    1919
Note: See TracChangeset for help on using the changeset viewer.