Changeset 5693


Ignore:
Timestamp:
09/07/10 14:36:42 (15 years ago)
Author:
Mathieu Morlighem
Message:

Removed all use of PenlatySystemMatrices except Thermal for debug

Location:
issm/trunk/src
Files:
12 edited

Legend:

Unmodified
Added
Removed

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

    r5689 r5693  
    44
    55#include "./SystemMatricesx.h"
    6 
    76#include "../../shared/shared.h"
    87#include "../../include/include.h"
    … …  
    109#include "../../EnumDefinitions/EnumDefinitions.h"
    1110
    12 void SystemMatricesx(Mat* pKgg, Vec* ppg,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters, int kflag,int pflag){
     11void SystemMatricesx(Mat* pKgg, Vec* ppg,double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){
    1312       
    1413        /*intermediary: */
    15         int gsize;
    16         int i,j;
    17         int connectivity;
    18         int numberofdofspernode;
    19         Element* element=NULL;
    20         Load*    load=NULL;
     14        int      i,j,gsize;
     15        int      connectivity, numberofdofspernode;
     16        int      analysis_type, configuration_type;
     17        Element *element = NULL;
     18        Load    *load    = NULL;
    2119       
    2220        /*output: */
    23         Mat Kgg=NULL;
    24         Vec pg=NULL;
    25 
    26         int analysis_type;
    27         int configuration_type;
    28 
    29         /*Display message*/
    30         int verbose; parameters->FindParam(&verbose,VerboseEnum);
    31         if (verbose) printf("   Generating matrices\n");
     21        Mat    Kgg  = NULL;
     22        Vec    pg   = NULL;
     23        double kmax = 0;
    3224
    3325        /*retrive parameters: */
    3426        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3527        parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    36 
    37         /*Recover parameters: */
    3828        parameters->FindParam(&connectivity,ConnectivityEnum);
    3929
    4030        /*Get size of matrix: */
    4131        gsize=nodes->NumberOfDofs(configuration_type);
    42 
    43         /*Get numberofdofspernode: */
    4432        numberofdofspernode=nodes->MaxNumDofs(configuration_type);
    4533
    46         /*Compute stiffness matrix*/
     34        /*Checks in debugging mode {{{1*/
     35        ISSMASSERT(*pKgg==NULL);
     36        ISSMASSERT(*ppg==NULL);
     37        if(penalty_kflag)ISSMASSERT(kflag);
     38        if(penalty_pflag)ISSMASSERT(pflag);
     39        /*}}}*/
     40
     41        /*Compute penalty free mstiffness matrix and load vector*/
    4742        if(kflag){
    4843
    49                 /*Allocate Kgg*/
    5044                Kgg=NewMat(gsize,gsize,NULL,&connectivity,&numberofdofspernode);
    5145
    … …  
    6761                MatCompress(Kgg);
    6862        }
    69 
    70         /*Compute Load vector*/
    7163        if(pflag){
    7264
    73                 /*Allocate pg*/
    7465                pg=NewVec(gsize);
    7566
    … …  
    8677                }
    8778
    88                 /*Assemble right hand side: */
    8979                VecAssemblyBegin(pg);
    9080                VecAssemblyEnd(pg);
    9181        }
    9282
     83        /*Now, deal with penalties*/
     84        if(penalty_kflag){
     85
     86                /*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */
     87                MatNorm(Kgg,NORM_INFINITY,&kmax);
     88
     89                /*Fill stiffness matrix from loads: */
     90                for (i=0;i<loads->Size();i++){
     91                        load=(Load*)loads->GetObjectByOffset(i);
     92                        if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kgg,kmax);
     93                }
     94
     95                /*Assemble matrix and compress matrix to save memory: */
     96                MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
     97                MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY);
     98                MatCompress(Kgg);
     99        }
     100        if(penalty_pflag){
     101
     102                /*Fill right hand side vector, from loads: */
     103                for (i=0;i<loads->Size();i++){
     104                        load=(Load*)loads->GetObjectByOffset(i);
     105                        if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pg,kmax);
     106                }
     107
     108                VecAssemblyBegin(pg);
     109                VecAssemblyEnd(pg);
     110        }
    93111       
    94112        /*Assign output pointers: */
    95113        *pKgg=Kgg;
    96114        *ppg=pg;
    97 
     115        if(pkmax) *pkmax=kmax;
    98116}

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

    r4236 r5693  
    1010
    1111/* local prototypes: */
    12 void SystemMatricesx(Mat* pKgg, Vec* ppg,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters, int kflag,int pflag);
     12void SystemMatricesx(Mat* pKgg, Vec* ppg,double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,
     13                        bool kflag=true,bool pflag=true,bool penalty_kflag=true,bool penalty_pflag=true);
    1314
    1415#endif  /* _SYSTEMMATRICESX_H */
    15 

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

    r5689 r5693  
    2626        Vec pf=NULL;
    2727
    28         kflag=1; pflag=1;
    29         SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    30         PenaltySystemMatricesx(Kgg, pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
     28        SystemMatricesx(&Kgg, &pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    3129
    3230        Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);

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

    r5689 r5693  
    3434
    3535        /*parameters:*/
    36         int kflag,pflag;
    3736        int verbose=0;
    3837
    3938        /*Recover parameters: */
    40         kflag=1; pflag=1;
    4139        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    4240        femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
    … …  
    6765                VecFree(&old_uf);old_uf=uf;
    6866
    69                 SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    70                 PenaltySystemMatricesx(Kgg, pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters,kflag,pflag);
     67                SystemMatricesx(&Kgg, &pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
    7168               
    7269                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters);

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

    r5689 r5693  
    99
    1010void solver_linear(FemModel* femmodel){
    11 
    12         /*parameters:*/
    13         int kflag,pflag;
    14         int verbose=0;
    1511
    1612        /*output: */
    … …  
    2521        Vec pf=NULL;
    2622
    27         /*Recover parameters: */
    28         kflag=1; pflag=1;
    29         femmodel->parameters->FindParam(&verbose,VerboseEnum);
    30 
    31         SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    32         PenaltySystemMatricesx(Kgg, pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
     23        SystemMatricesx(&Kgg,&pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    3324
    3425        Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);

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

    r5689 r5693  
    5454                InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,reset_penalties,ResetPenaltiesEnum);
    5555
    56                 //*Generate system matrices
    57                 if (!lowmem){
    58 
    59                         /*Compute Kgg_nopenalty and pg_nopenalty once for all: */
    60                         if (count==1){
    61                                 SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    62                         }
    63 
    64                         /*Copy K_gg_nopenalty into Kgg, same for pg: */
    65                         MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg);
    66                         VecDuplicatePatch(&pg,pg_nopenalty);
    67 
    68                         //apply penalties each time
    69                         PenaltySystemMatricesx(Kgg, pg,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    70                 }
    71                 else{
    72                         SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    73                         //apply penalties
    74                         PenaltySystemMatricesx(Kgg, pg,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    75                 }
     56                /*===================== DEBUGING ====================*/
     57                if (count==1) SystemMatricesx(&Kgg_nopenalty,&pg_nopenalty,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     58                MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg);
     59                VecDuplicatePatch(&pg,pg_nopenalty);
     60                PenaltySystemMatricesx(Kgg, pg,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kflag,pflag);
     61                /*===================== DEBUGING ====================*/
     62                //SystemMatricesx(&Kgg,&pg,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     63                /*===================== DEBUGING ====================*/
    7664
    7765                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters);

  • issm/trunk/src/m/solvers/solver_adjoint_linear.m

    r5689 r5693  
    55%      femmodel =solver_adjoint_linear(femmodel)
    66
    7         femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=1;
    8 
    9         [K_gg, p_g]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    10         [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     7        [K_gg, p_g, kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    118       
    129        [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters);

  • issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m

    r5689 r5693  
    44%   Usage:
    55%      [femmodel]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
    6        
    76
    87        %Recover parameters
    9         femmodel.parameters.Kflag=1;
    10         femmodel.parameters.Pflag=1;
    11         dim=femmodel.parameters.Dim;
    128        min_mechanical_constraints=femmodel.parameters.MinMechanicalConstraints;
    139
    … …  
    3127                old_uf=uf;
    3228               
    33                 [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
    34                 [K_gg , p_g, kmax]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
     29                [K_gg, p_g, kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
    3530               
    3631                [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters);

  • issm/trunk/src/m/solvers/solver_linear.m

    r5689 r5693  
    55%      femmodel =solver_linear(femmodel)
    66
    7         femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=1;
    8 
    9         [K_gg, p_g]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    10         [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     7        [K_gg, p_g,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    118       
    129        [K_ff, K_fs] = Reducematrixfromgtof( K_gg,  femmodel.nodesets,femmodel.parameters);

  • issm/trunk/src/m/solvers/solver_thermal_nonlinear.m

    r5689 r5693  
    1919                [femmodel.elements femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,reset_penalties,ResetPenaltiesEnum);
    2020
    21                 %system matrices
    22                 if count==1
    23                         displaystring(femmodel.parameters.Verbose,'%s',['   system matrices']);
    24                         [K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     21                % ================= DEBUG FOR NOW ====================
     22                if count==1
     23                        [K_gg_nopenalty,p_g_nopenalty,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    2524                end
    26                 displaystring(femmodel.parameters.Verbose,'%s',['   penalty system matrices']);
    2725                [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     26                % ================= DEBUG FOR NOW ====================
     27                %[K_gg,p_g,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
     28                % ================= DEBUG FOR NOW ====================
    2829
    2930                [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters);
    … …  
    3536                t_g= Mergesolutionfromftog( t_f, femmodel.ys, femmodel.nodesets,femmodel.parameters);
    3637                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,t_g);
    37        
     38
    3839                [femmodel.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters);
    3940       

  • issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp

    r5689 r5693  
    1414        Materials  *materials  = NULL;
    1515        Parameters *parameters = NULL;
    16         int         kflag,pflag;
    1716       
    1817        /* output datasets: */
    19         Mat Kgg=NULL;
    20         Vec pg=NULL;
     18        Mat    Kgg  = NULL;
     19        Vec    pg   = NULL;
     20        double kmax;
    2121
    2222        /*Boot module: */
    … …  
    3434        FetchParams(&parameters,PARAMETERS);
    3535
    36         /*parameters: */
    37         parameters->FindParam(&kflag,KflagEnum);
    38         parameters->FindParam(&pflag,PflagEnum);
    39 
    4036        /*configure: */
    4137        elements->  Configure(elements,loads, nodes,vertices, materials,parameters);
    … …  
    4541
    4642        /*!Generate internal degree of freedom numbers: */
    47         SystemMatricesx(&Kgg, &pg,elements,nodes,vertices,loads,materials,parameters,kflag,pflag);
     43        SystemMatricesx(&Kgg,&pg,&kmax,elements,nodes,vertices,loads,materials,parameters);
    4844
    4945        /*write output datasets: */
    5046        WriteData(KGG,Kgg);
    5147        WriteData(PG,pg);
     48        WriteData(KMAX,kmax);
    5249       
    5350        /*Free ressources: */
    … …  
    6865{
    6966        _printf_("\n");
    70         _printf_("   usage: [Kgg,pg] = %s(elements,nodes,vertices,loads,materials,parameters);\n",__FUNCT__);
     67        _printf_("   usage: [Kgg,pg,kmax] = %s(elements,nodes,vertices,loads,materials,parameters);\n",__FUNCT__);
    7168        _printf_("\n");
    7269}

  • issm/trunk/src/mex/SystemMatrices/SystemMatrices.h

    r4236 r5693  
    1818
    1919/* serial input macros: */
    20 #define ELEMENTS (mxArray*)prhs[0]
    21 #define NODES (mxArray*)prhs[1]
    22 #define VERTICES (mxArray*)prhs[2]
    23 #define LOADS (mxArray*)prhs[3]
    24 #define MATERIALS (mxArray*)prhs[4]
    25 #define PARAMETERS (mxArray*)prhs[5]
     20#define ELEMENTS   (mxArray *)prhs[0]
     21#define NODES      (mxArray *)prhs[1]
     22#define VERTICES   (mxArray *)prhs[2]
     23#define LOADS      (mxArray *)prhs[3]
     24#define MATERIALS  (mxArray *)prhs[4]
     25#define PARAMETERS (mxArray *)prhs[5]
    2626
    2727/* serial output macros: */
    28 #define KGG (mxArray**)&plhs[0]
    29 #define PG (mxArray**)&plhs[1]
     28#define KGG  (mxArray**)&plhs[0]
     29#define PG   (mxArray**)&plhs[1]
     30#define KMAX (mxArray**)&plhs[2]
    3031
    3132/* serial arg counts: */
    3233#undef NLHS
    33 #define NLHS  2
     34#define NLHS  3
    3435#undef NRHS
    3536#define NRHS  6
Note: See TracChangeset for help on using the changeset viewer.